Estimation d'état des systèmes non-linéaires à temps continu par une approche ensembliste

01/10/2017
Auteurs : Luc Jaulin
Publication e-STA e-STA 2003-1
OAI : oai:www.see.asso.fr:545:2003-1:20079
DOI : You do not have permission to access embedded form.

Résumé

Estimation d'état des systèmes non-linéaires à temps continu par une approche ensembliste

Métriques

12
5
294.67 Ko
 application/pdf
bitcache://48e5a984e10fcc83fa3270c381ae01448bdf4d83

Licence

Creative Commons Aucune (Tous droits réservés)
<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xmlns="http://datacite.org/schema/kernel-4"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/545:2003-1/20079</identifier><creators><creator><creatorName>Luc Jaulin</creatorName></creator></creators><titles>
            <title>Estimation d'état des systèmes non-linéaires à temps continu par une approche ensembliste</title></titles>
        <publisher>SEE</publisher>
        <publicationYear>2017</publicationYear>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><dates>
	    <date dateType="Created">Sun 1 Oct 2017</date>
	    <date dateType="Updated">Sun 1 Oct 2017</date>
            <date dateType="Submitted">Sat 24 Feb 2018</date>
	</dates>
        <alternateIdentifiers>
	    <alternateIdentifier alternateIdentifierType="bitstream">48e5a984e10fcc83fa3270c381ae01448bdf4d83</alternateIdentifier>
	</alternateIdentifiers>
        <formats>
	    <format>application/pdf</format>
	</formats>
	<version>34096</version>
        <descriptions>
            <description descriptionType="Abstract"></description>
        </descriptions>
    </resource>
.

Estimation d’état des systèmes non-linéaires à temps continu par une approche ensembliste Luc Jaulin1 1Laboratoire d’Ingénierie des Systèmes Automatiséss, 62 avenue Notre Dame du Lac 49 000 Angers, France Luc.Jaulin@univ-angers.fr http://www.istia.univ-angers.fr/~jaulin/ Résumé— Cet article illustre l’utilisation de l’analyse par in- tervalles et des méthodes de propagation de contraintes pour l’estimation d’état de systèmes à temps continu décrit par des équations différentielles ordinaires. Cette approche est présentée dans un contexte ensembliste, c’est-à-dire que les incertitudes sur les quantités manipulées sont représentées par un ensemble et non plus par une loi de probabilité. Mots-clés— estimation à erreurs bornées, propagation de contraintes, CSP, identification non-linéaire, analyse par in- tervalles, estimation d’état, estimation ensembliste. I. Introduction Considérons le problème de l’estimation de l’état d’un système non-linéaire décrit par ses équations d’état ½ ẋ(t) = f(x(t)), y(t) = g(x(t)) (1) où – t est le temps, supposé dans l’intervalle [t, t̄] ⊂ R – x(t) ∈ Rn est le vecteur d’état à l’instant t, – et y(t) ∈ Rm est le vecteur de sortie à l’instant t. Ce problème a déjà été abondamment étudié dans le cas où les fonctions f et g sont toutes les deux linéaires [1], [2], ou lorsque le temps t est discret [3], [4] ou encore dans un contexte probabiliste. Pourtant, à ma connaissance, dans un contexte à erreurs bornées et pour le cas des systèmes non-linéaires à temps continus, l’estimation d’état n’a ja- mais été traitée, sauf bien sur dans l’article [5] dont cet article reprend les principaux résultats. Les outils que nous allons utiliser sont l’analyse par intervalles [6] qui nous per- mettra d’obtenir des résultats garantis et les techniques de propagation de contraintes [7] , [8] qui nous permettront d’améliorer l’efficacité des algorithmes par intervalles. Nous allons supposer que la fonction f est différentiable et que pour tout t, des domaines [xi] (t) et [yj] (t) contenant les variables xi(t) et yj(t) sont disponibles a priori. Ces domaines peuvent être arbitrairement grands, par exemple [xi] (t) =] − ∞, ∞[ si aucune information n’est disponible sur xi(t). Ils peuvent aussi être arbitrairement petits. Ainsi, si une mesure ŷj(t0) de la variable de sortie yj(t0) a été recueillie à l’instant t = t0 et si une borne supérieure ē sur la valeur absolue de l’erreur ŷj(t0) − yj(t0) est connue, alors, le domaine pour yj(t0) sera [yj] (t0) = [ŷj(t0) − ē, ŷj(t0) + ē]. (2) Le paragraphe II rappelle les notions de base sur l’ana- lyse par intervalles qui seront utilisées dans cet article. Le paragraphe III introduit une méthode élémentaire permet- tant de traiter simplement avec des équations différentielles ordinaires. Un algorithme permettant de faire l’estimation d’état du système (1) sera proposé dans le paragraphe IV. Enfin, un exemple sera donné au paragraphe V. II. Calcul sur les intervalles Un intervalle [x] = [x, x] est un ensemble fermé et connexe de R. Les opérations sur les nombres réels peuvent être étendues aux intervalles [6]. En effet, tout opérateur binaire ¦ tel que +, −, ∗, / et toute fonction élémentaire f telles que exp, sin, sqr, sqrt, sur les nombres réels se généralisent aux intervalles de la façon suivante : [x] ¦ [y] = [infx∈[x],y∈[y] x ¦ y, supx∈[x],y∈[y] x ¦ y] f([x]) = [infx∈[x] f(x), supx∈[x] f(x)] (3) Par exemple, [x] + [y] = [x + y, x + y] [x] ∗ [y] = [min(xy, xy, xy, xy) , max(xy, xy, xy, xy)] exp ([x]) = [exp (x) , exp (x)], p [x] = [ p max(x, 0), √ x] si x ≥ 0 et ∅ si x < 0 (4) Un pavé [x] de Rn est le produit cartésien de n intervalles. Il peut s’écrire sous la forme [x] = [x1, x̄1] × · · · × [xn, x̄n] = [x, x̄], (5) où x = (x1, · · · , xn)T et x̄ = (x1, · · · , x̄n)T . L’ensemble de tous les pavés de Rn sera noté IRn . Remarquons que l’ensemble Rn =] − ∞, ∞[× · · · ×] − ∞, ∞[ (6) est un élément de IRn . La longueur w([x]) du pavé [x] est la longueur de son plus long côté. Couper [x] signifie : ”couper en deux suivant l’axe de symétrie perpendiculaire au côté le plus long”. Soit f : Rn → Rp une fonction vectorielle. La fonction ensembliste [f] : IRn → IRp est une fonction d’inclusion de f si f([x]) ⊂ [f]([x]) (7) pour tout [x] ∈ IRn . Une fonction d’inclusion est conver- gente si, pour toute suite de pavés [x] de IRn , on a w([x]) → 0 ⇒ w([f]([x])) → 0, (8) Différentes méthodes existent pour obtenir des fonctions d’inclusion convergentes, [9]. La plus simple d’entre elles et qui permet de traiter une très grande classe de fonc- tion consiste à remplacer toutes les variables d’entrées de la fonction par des intervalles et tous les opérateurs et les fonc- tions élémentaires définissant la fonction par les opérateurs et fonctions intervalles associés. Illustrons cela par deux exemples. Exemple 1: Soit la fonction f(x1, x2) = 1 √ 2πx2 e− 1 2 (t−x1)2 x2 , (9) où t est un nombre réel donné. Une fonction d’inclusion [f] de f est donnée par [f]([x1], [x2]) = 1 p 2π[x2] e − 1 2 (t−[x1])2 [x2] , (10) où toutes les opérations utilisées sont celles de l’arithmétique des intervalles. Exemple 2: Considérons la fonction vectorielle f(x1, x2) , µ (1 − 0.01x2) x1 (−1 + 0.02x1) x2 ¶ . (11) Le calcul suivant montre comment se calcule un pavé conte- nant l’ensemble f ([x]), où [x] = ([10, 20], [40, 50]) : f ([x]) ⊂ µ (1 − 0.01 ∗ [40, 50]) ∗ [10, 20] (−1 + 0.02 ∗ [10, 20]) ∗ [40, 50] ¶ = µ (1 − [0.4, 0.5]) ∗ [10, 20] (−1 + [0.2, 0.4]) ∗ [40, 50] ¶ = µ [0.5, 0.6] ∗ [10, 20] [−0.8, −0.6] ∗ [40, 50] ¶ = µ [5, 12] [−40, −24] ¶ . Bien sûr, ce calcul peut être appliqué à tout pavé [x] = [x1]×[x2]. L’algorithme résultant correspond à une fonction d’inclusion pour f. ¥ Soient [x] et [y] deux pavés de Rn . Nous noterons [x]t[y] le plus petit pavé contenant l’union [x]∪[y] des deux pavés. Gonfler un pavé [x] = [x1, x1] × · · · × [xn, xn] de ε > 0 consiste à le remplacer par le pavé gonfle ([x], ε) , [x1 − ε, x1 + ε] × · · · × [xn − ε, xn + ε]. III. Théorème de Picard L’utilisation de l’analyse par intervalles pour le trai- tement d’équations différentielles a été proposée pour la première fois par Moore [6] (voir [10] pour un survol du su- jet). Elle permet de développer des méthodes numériques qui fournissent des ensembles contenant à coup sûr les so- lutions exactes des équations différentielles aux instants d’échantillonnage t0, t1, . . . lorsqu’un pavé [x](0) contenant l’état initial x(t0) est donné. L’approche est fondée sur le théorème de Picard que nous allons maintenant énoncer. Théorème 1: Soient t1 et t2 deux nombres réels. Suppo- sons que x(t1) appartienne au pavé [x](t1). Soit [w] un pavé Fig. 1. Illustration du principe de l’algorithme [φ] ; les flèches représentent le champs de vecteur associé à la fonction d’évolution f du système de IRn . Dans notre contexte, ce pavé sera choisi pas trop gros et de telle façon qu’il ait de fortes chances de contenir la trajectoire x(τ) pour tout τ ∈ [t1, t2]. Si [x](t1) + [0, t2 − t1] ∗ [f]([w]) ⊂ [w], où [f](.) est une fonction d’inclusion pour f et où [0, t2 −t1] est le plus petit intervalle contenant 0 et t2 −t1, alors pour tout τ ∈ [t1, t2], on a x(τ) ∈ [x](t1) + [0, t2 − t1] ∗ [f]([w]), (12) De plus x(t2) ∈ [x](t1) + (t2 − t1) ∗ [f]([w]). (13) Nous pouvons maintenant construire un algorithme [φ] capable de calculer un pavé [x](t2) contenant x(t2) à partir d’un pavé [x](t1) dont on sait qu’il contient x(t1). La figure 1 illustre le principe de cet algorithme. Fonction [φ] (in : t1, t2, [x](t1), out : [x](t2)) 1 [x̂](t2) := [x](t1) + (t2 − t1) ∗ [f]([x](t1)); 2 [v] := [x](t1) t [x̂](t2); 3 [w] := gonfle([v], η); //où η > 0 est petit 4 si [x](t1) + [0, t2 − t1] ∗ [f]([w])) * [w] 5 {[x](t2) := Rn ; fin} ; 6 [x](t2) := [x](t1) + (t2 − t1) ∗ [f]([w]). Le pas 1 calcule un pavé [x̂](t2) qui pourrait contenir l’ensemble de tous les x(t2) tels que x(t1) ∈ [x](t1), mais aucune garantie n’est offerte. Le pas 2 calcule le plus petit pavé [v] contenant [x](t1) et [x̂](t2). Le pavé [v] est alors légèrement gonflé (de η > 0) au pas 3 pour obtenir le pavé [w]. Si, au pas 4, la condition du théorème 1 n’est pas satis- faite, aucun pavé fini contenant x(t2) ne peut être calculé et Rn se trouve alors renvoyé au pas 5. Sinon, le pas 6 calcule un pavé contenant x(t2) grâce à la propriété (13). Remarque 1: Le coefficient de gonflement η de la fonc- tion [φ] doit être choisi suffisamment faible de façon à avoir [w] à peine plus gros que [v], mais pas nul car sinon, la condition du pas 4 a peu de chance d’être satisfaite. Dans l’exemple traité dans la suite de cet article, nous avons choisi η = 0.1.w([v]) + 0.0001 mais beaucoup d’autres stratégies pour l’obtention de ce coefficient donnent des Fig. 2. Superposition de pavés contenant toutes les trajectoires pos- sibles du système sachant que le vecteur initial appartient au pavé [x] (0) ; cette superposition a été générée par l’algorithme Picard résultats tout aussi satisfaisants. Une recherche dichoto- mique pourrait aussi être envisagée pour obtenir un η convenable de façon automatique. L’algorithme suivant utilise l’opérateur [φ] pour calculer des pavés contenant à coup sûr les quantités x(δ), x(2δ), . . . , x(k̄δ), où δ est la période d’échantillonnage et où k est le plus grand entier inférieur à t/δ. Ce calcul est fait à partir de la connaissance d’un pavé [x](0) contenant x(0). Pour une question de concision, la notation [x](k) sera utilisée à la place [x](kδ). Picard (in : [x](0), out : [x](k), k ≥ 1) 1 pour k := 0 jusqu’à k̄ − 1 2 [x](k + 1) := [φ] (kδ, (k + 1)δ, [x](k)) Exemple 3: Considérons le modèle de Lotka-Volterra donné par : ½ ẋ1 = (1 − 0.01x2) x1 ẋ2 = (−1 + 0.02x1) x2 (14) Pour [x](0) = [49.5; 50.5] × [49.5; 50.5], δ = 0.005 et k̄ = 400, l’algorithme Picard génère un ensemble de pavés [x](k), dont la superposition est représentée sur la figure 2. Pour tout k ∈ {1, . . . , k̄}, nous avons l’implication x(0) ∈ [x](0) ⇒ x(k) ∈ [x](k). (15) Les x(k) en représentés en blanc sur la figure ont été obte- nus pour x(0) = (50, 50). L’exemple 3 illustre le principal inconvénient des méthodes par intervalles pour le traitement des équations différentielles, à savoir, l’explosion rapide de la taille des pavés [x](k) lorsque le temps augmente. Cette explo- sion peut être réduite considérablement en utilisant les méthodes de Lohner [11]. Dans le contexte de l’estimation d’état, cette explosion peut aussi être réduite grâce aux mesures prélevées sur le système à différents instants. Ce dernier point fera l’objet du prochain paragraphe. IV. Algorithme pour l’estimation d’état Considérons à nouveau le système non-linéaire (1), où des domaines (qui sont en fait des intervalles) [xi] (t) et [yj] (t) pour les variables xi(t) et yj(t) sont supposés connus. Rap- pelons que ces intervalles sont égaux à ] − ∞, ∞[ si au- cune information a priori n’est disponible sur la variable as- sociée. Aux instants d’échantillonnage, t = kδ, les variables xi(t) et yj(t) seront notées xi(k) et yj(k). Les techniques de propagation permettent de contracter les domaines pour les variables en supprimant des parties des domaines qui sont incompatibles avec (1) et les domaines pour les autres variables. L’utilisation de ces méthodes pour le traitement des équations différentielles a été proposée pour la première fois dans [12]. Les contractions considérées utilisent les im- plications suivantes. (i) x(k) ∈ [x](k) =⇒ x(k + 1) ∈ [φ](kδ, (k + 1)δ, [x](k)), (ii) x(k) ∈ [x](k) =⇒ y(k) ∈ g([x](k)), (iii) y(k) ∈ [y](k) =⇒ x(k) ∈ g−1 ([y](k)), (iv) x(k + 1) ∈ [x](k + 1) =⇒ x(k) ∈ [φ]((k + 1)δ, kδ, [x](k + 1)). Chacune de ces implications engendre une procédure de contraction de l’algorithme présenté ci-dessous. Contract(inout : [x](k), [y](k), k ≥ 0) 1 pour k := 0 jusqu’à k̄, 2 [x](k + 1) := [x](k + 1) ∩[φ](kδ, (k + 1)δ, [x](k)); 3 [y](k) := [y](k) ∩ [g]([x](k)); 4 pour k := k̄ jusqu’à 0, 5 [x](k) := [x](k) ∩ [g−1 ]([y](k)); 6 [x](k) := [x](k) ∩[φ]((k + 1)δ, kδ, [x](k + 1)); Au pas 3, [g](.) représente une fonction d’inclusion pour g(.). La contraction pour [x](k) au pas 5 peut se faire grâce à l’algorithme de propagation-rétropropagation (voir par exemple [4]). Cette procédure doit être appelée plu- sieurs fois, tant que des contractions significatives peuvent être effectuées. Lorsque l’équilibre est atteint, des do- maines contenant encore de grosses parties invraisem- blables peuvent encore subsister. Un partitionnement doit alors être fait pour les éliminer. La stratégie que nous al- lons suivre est le partitionnement de [x](k0), pour un k0 donné, en petits pavés. L’algorithme résultant est décrit ci-dessous. Strangle(in : k0, ε, inout : [x](k), [y](k), k ≥ 0) 1 partitionner [x](k0) en ` pavés [x1](k0), . . . , [x`](k0) de longueur plus petite que ε ; 2 pour i := 1 jusqu’à `, 3 pour tout k 6= k0 4 {[xi](k) := [x](k); [yi](k) := [y](k)} ; 5 Contract([xi](k), [yi](k), k ≥ 0) ; 6 pour k := 0 jusqu’à k̄, 7 [x](k) := [x1](k) t · · · t [x`](k); [y](k) := [y1](k) t · · · t [y`](k); Même pour des valeurs infiniment petites de ε, les do- maines pour [x](k) et [y](k) peuvent encore contenir de grosses parties inconsistantes. Le principe de l’algorithme StateBound ci-dessous est d’appeler Strangle pour différents k0 et de faire décroı̂tre ε dans le but de contracter les domaines le plus possible. StateBound(in : ε, inout : [x](k), [y](k), k ≥ 0) 1 faire 2 faire 3 choisir aléatoirement k0 ∈ {0, 1, . . . , k̄}; 4 Strangle(k0, ε, [x](k), [y](k), k ≥ 0) ; 5 tant que la contraction est significative ; 6 ε := ε/2 ; 7 tant que la contraction est significative. Remarque 2: Par simplicité, dans l’exemple du para- graphe V, la boucle externe a été remplacée par la boucle ”pour n1 := 1 jusqu’à n̄1”, où n̄1 est un entier préalablement défini. La boucle interne a été remplacée par la boucle ”pour n2 := 1 jusqu’à 30”. ¥ V. Exemple Considérons à nouveau le modèle de Lotka-Volterra. L’équation d’observation est donnée par y = x1. Les données ŷ(t) ont été obtenues en simulant le modèle avec x(0) = (50, 50) et en y ajoutant un bruit blanc avec une dis- tribution uniforme dont le support est l’intervalle [−1, −1]. Les domaines [y](t) pour y(t) ont été pris égaux à [y](t) = [ŷ(t) − 1.5, ŷ(t) + 1.5], (16) pour prendre en compte le fait que les bornes exactes pour les erreurs de mesures sont rarement disponibles. Les domaines a priori disponibles pour x(t) ont été pris égaux à [0, 100] × [0, 200]. Pour δ = 0.03 et k = 200, Strangle(k0 = 0, ε = 4) génère les [x](k) représentés sur la figure 3 (a), en 1.33 sec. sur un Pentium 300. Notons que pour k faible, les domaines [x](k) sont plus précis que pour les grands k. Lançons Strangle ¡ k0 = k̄, ε = 4 ¢ . Pour des grands k les domaines ont maintenant été considérablement contractés (voir la figure 3 (b)). Pour n̄1 = 1 (c’est-à-dire que la boucle externe n’est parcourue qu’une seule fois), StateBound(ε = 4) génère des domaines de la Figure 3 (c) en 11.4 sec. VI. Conclusion Cette article étudie, pour la première fois, les applica- tions de l’analyse par intervalles et des méthodes de propa- gation pour l’estimation d’état de systèmes non-linéaires décrits par des équations différentielles dans un contexte ensembliste. Un algorithme efficace, capable d’enfermer tous les vecteurs d’état vraisemblables à l’intérieur d’en- sembles constitués d’unions de pavés, a été proposé. Son efficacité pourrait être améliorée grâce à l’utilisation de méthodes de consistances plus fortes [12] et les techniques de simulation ensemblistes plus précises [11]. Notons que même si la méthode d’enfermement pro- posée semble efficace pour notre exemple relativement académique, nous n’avons, à ce stade, aucune possi- bilité d’évaluer la qualité de l’enfermement généré, ni même d’imaginer les performances de l’approche pour des systèmes plus complexes. De plus, l’extension de l’approche pour l’estimation d’état de systèmes incertains n’est a Fig. 3. Domaines [x](k) contractés par Strangle priori immédiate que pour des incertitudes de nature pa- ramétrique. Qu’en est-il pour des perturbations d’une autre nature ? Beaucoup de questions restent donc à résoudre avant de pouvoir espérer traiter des applications réelles. Les codes C++ (en Builder 3) complets correspondant à l’exemple du paragraphe V sont disponibles sur http :// www.istia.univ-angers.fr /~jaulin/voltera cpp.zip. Références [1] F. C. Schweppe. Recursive state estimation : unknown but boun- ded errors and system inputs. IEEE Transactions on Automatic Control, 13(1) :22–28, 1968. [2] F. L. Chernousko. State Estimation for Dynamic Systems. CRC Press, Boca Raton, FL, 1994. [3] G. Chen, J. Wang, et L. S. Shieh. Interval Kalman filtering. IEEE Trans. on Aerospace and Electronic Systems, 33(1) :250– 258, 1997. [4] L. Jaulin, M. Kieffer, O. Didrit, et E. Walter. Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics. Springer-Verlag, London, 2001. [5] L. Jaulin. Nonlinear bounded-error state estimation of continuous-time systems. to appear in Automatica, 2002. [6] R. E. Moore. Interval Analysis. Prentice-Hall, Englewood Cliffs, NJ, 1966. [7] E. Davis. Constraint propagation with interval labels. Artificial Intelligence, 32(3) :281–331, 1987. [8] D.J. Sam-Haroud et B. Faltings. Consistency techniques for continuous constraints. Constraints, 1(1-2) :85–118, 1996. [9] R. E. Moore. Methods and Applications of Interval Analysis. SIAM, Philadelphia, PA, 1979. [10] M. Berz, C. Bischof, G. Corliss, et Griewank A., editors. Com- putational Differentiation : Techniques, Applications and Tools. SIAM, Philadelphia, Penn., 1996. [11] R. Lohner. Enclosing the solutions of ordinary initial and boun- dary value problems. E. Kaucher, U. Kulisch, et Ch. Ullrich, editors, Computer Arithmetic : Scientific Computation and Pro- gramming Languages, pages 255–286. BG Teubner, Stuttgart, Germany, 1987. [12] Y. Deville, M. Janssen, et P. Van Hentenryck. Consistency techniques in ordinary differential equations. Proceedings of the Fourth International Conference on Principles and Practice of Constraint Programming, Lecture Notes in Computer Science. Springer Verlag, 1998.