Adaptation robuste de la méthode du propagateur à l’identification récursive des sous-espaces

30/09/2017
Publication e-STA e-STA 2005-1
OAI : oai:www.see.asso.fr:545:2005-1:20027
DOI :

Résumé

Adaptation robuste de la méthode du propagateur à l’identification récursive des sous-espaces

Métriques

18
6
204.86 Ko
 application/pdf
bitcache://9d3cf77cbc176360316b8c2f2332190c9aee7fe4

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:2005-1/20027</identifier><creators><creator><creatorName>Christian Vasseur</creatorName></creator><creator><creatorName>Guillaume Mercere</creatorName></creator><creator><creatorName>Stéphane Lecoeuche</creatorName></creator></creators><titles>
            <title>Adaptation robuste de la méthode du propagateur à l’identification récursive des sous-espaces</title></titles>
        <publisher>SEE</publisher>
        <publicationYear>2017</publicationYear>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><dates>
	    <date dateType="Created">Sat 30 Sep 2017</date>
	    <date dateType="Updated">Sat 30 Sep 2017</date>
            <date dateType="Submitted">Fri 20 Jul 2018</date>
	</dates>
        <alternateIdentifiers>
	    <alternateIdentifier alternateIdentifierType="bitstream">9d3cf77cbc176360316b8c2f2332190c9aee7fe4</alternateIdentifier>
	</alternateIdentifiers>
        <formats>
	    <format>application/pdf</format>
	</formats>
	<version>34043</version>
        <descriptions>
            <description descriptionType="Abstract"></description>
        </descriptions>
    </resource>
.

Adaptation robuste de la méthode du propagateur à l’identification récursive des sous-espaces Guillaume MERCERE1,2 , Stéphane LECŒUCHE1,3 , Christian VASSEUR1 1 Laboratoire d’Automatique, Génie Informatique et Signal Bâtiment P2, USTL, 59655 Villeneuve d’Ascq Cedex 2 Ecole d’Ingénieurs du Pas de Calais Campus de la Malassise, BP39, 62967 Longuenesse Cedex 3 Ecole des Mines de Douai 941, rue Charles Bourseul, BP 838, 59508 Douai Cedex gmercere@eipc.fr, lecoeuche@ensm-douai.fr, Christian.Vasseur@univ-lille1.fr Résumé— L’objectif de cet article est de présenter un nouvel algorithme d’identification récursive des sous-espaces dont la qualité première est sa stabilité numérique. Cette amé- lioration concerne plus précisément la phase de mise à jour d’un opérateur linéaire appelé propagateur, opérateur à par- tir duquel il est facile de construire une base de l’espace d’observabilité. Nous montrons que, par application de rota- tions hyperboliques et de Givens, il est possible d’imposer la définie-positivité de certaines matrices, condition nécessaire à la robustesse numérique de l’algorithme. Par ailleurs, la description d’une technique particulière d’estimation récur- sive d’un vecteur concentrant l’ensemble des informations des données d’entrée-sortie est également développée. Son principal avantage est de posséder un coût calculatoire plus faible que la mise à jour d’une factorisation QR utilisée jus- qu’ici. La partie application met en évidence les bénéfices de ces deux techniques sur un exemple de simulation. Mots-clés— Identification récursive, méthodes des sous- espaces, lemme d’inversion matricielle, rotations hyperbo- liques, rotations de Givens, stabilité numérique. I. Introduction Le problème d’identification que nous souhaitons ré- soudre consiste à estimer récursivement une réalisation d’état du système étudié à chaque mise à jour des don- nées d’entré-sortie perturbées. Les solutions proposées au sein de cet article sont fondées sur une approche de type sous-espaces. Elles s’inspirent plus particulièrement de la démarche suivie par les techniques d’identification MOESP [1]. Leur objectif premier est d’obtenir une estimée consis- tante de la matrice d’observabilité du système, directement à partir des données d’entrée-sortie, sans calcul explicite du vecteur d’état. Il est en effet aisé d’en extraire ensuite des estimées fiables des matrices d’état du système [2]. Puisque la représentation d’état est plurielle, nous nous intéressons plus précisément à l’estimation récursive d’une base de l’es- pace d’observabilité. Plusieurs solutions ont été proposées dans la littérature [3], [4], [5], [6]. L’une d’entre elles [7] repose sur l’adaptation de la méthode du propagateur [8] à notre problème d’identification. Cette technique, nommée EIVPM, permet d’obtenir une estimée consistante d’une base du sous-espace d’observabilité à partir d’un problème d’optimisation de type moindres carrés qu’il est facile de minimiser par des algorithmes récursifs simples. Elle s’ap- puie plus précisément sur deux étapes : – la mise à jour, en ligne et à partir des données d’entrée- sortie perturbées, d’un vecteur particulier nommé vec- teur d’observation, – l’obtention d’une estimée de la matrice d’observabilité par estimation récursive de l’opérateur propagateur. Jusqu’ici, la première phase reposait soit sur une approxi- mation [9], soit sur la mise à jour d’une factorisation QR [10]. Afin d’éviter l’imprécision de la première solution tout en atténuant le coût calculatoire de la seconde, nous pro- posons d’actualiser, à chaque acquisition, une projection orthogonale particulière (cf. §III) à l’aide du lemme d’in- version matricielle [11]. Enfin, puisque nous avons remar- qué qu’une des matrices de covariance exploitées par l’algo- rithme EIVPM pouvait être mal conditionnée et non semi- définie positive, nous exposons une amélioration de l’algo- rithme de minimisation en imposant des conditions de sta- bilité numérique (cf. §IV-B). Cette robustesse est atteinte en contraignant la définie positivité de cette matrice de co- variance via l’utilisation de rotations hyperboliques et de Givens [11]. II. Problématique et notations Considérons la représentation d’état d’un modèle discret linéaire d’ordre n défini par : ½ x(t + 1) = Ax(t) + Bũ(t) + w(t) ỹ(t) = Cx(t) + Dũ(t) (1) où ũ(t) ∈ Rm×1 et ỹ(t) ∈ Rl×1 sont respectivement les vec- teurs d’entrée et de sortie non perturbées, x(t) ∈ Rn×1 le vecteur d’état et w(t) ∈ Rn×1 le bruit d’état. Les vecteurs d’entrée et de sortie sont caractérisés de la façon suivante : u(t) = ũ(t) + f(t) y(t) = ỹ(t) + v(t) (2) avec f(t) ∈ Rm×1 et v(t) ∈ Rl×1 des bruits de mesure. Les bruits sont supposés être trois vecteurs de bruit blanc statistiquement indépendants. Nous faisons de plus l’hypo- thèse que la réalisation d’état du modèle est minimale. Comme précisé dans l’introduction, les méthodes déve- loppées au sein de cet article s’inspirent d’une technique particulière du traitement du signal. Afin de mettre en évi- dence l’analogie entre la problématique d’identification et celle rencontrée en traitement d’antennes [12], introduisons le vecteur de sorties décalées yi(t̄) ∈ Rli×1 suivant : yi(t̄) = £ yT (t) · · · yT (t + i − 1) ¤T (3) où i est un entier fixé par l’utilisateur tel que i > n et t̄ = t + i − 1. Il est alors facile de montrer récursivement que ce vecteur vérifie1 : yi(t̄) = Γix(t̄ − i + 1) + Hiui(t̄) + Hw i wi(t̄) − Hifi(t̄) + vi(t̄) | {z } bi(t̄) (4) avec Γi ∈ Rli×n la matrice d’observabilité du système : Γi = £ CT (CA)T · · · (CAi−1 )T ¤ , (5) Hi ∈ Rli×mi la matrice de Toeplitz contenant les para- mètres de Markov du procédé et Hw i ∈ Rli×ni la matrice de Toeplitz associée au bruit d’état. La similitude entre le problème d’identification par approche des sous-espaces et celui de la détection des directions d’arrivée est alors natu- relle [10] si l’équation (4) est réécrite de la façon suivante : zi(t̄) = yi(t̄) − Hiui(t̄) = Γix(t̄ − i + 1) + bi(t̄). (6) Cette analogie permet alors de décomposer notre problème d’identification en deux étapes : 1. la mise à jour du vecteur d’observation zi à partir des nouvelles données d’entrée-sortie (cf. §III) : zi(t̄) = yi(t̄) − Hiui(t̄), (7) 2. l’estimation d’une base de la matrice d’observabilité (cf. §IV) à partir de ce nouveau vecteur en adaptant la méthode du propagateur [7], [8]. Sur ces deux étapes vient se greffer une troisième phase im- pérative en identification des systèmes : la suppression des effets des perturbations. En effet, en traitement d’antennes, la variance du bruit est, la plupart du temps, supposée pro- portionnelle à l’identité, hypothèse à laquelle ne déroge pas la méthode du propagateur [8], [13]. Cette condition étant trop restrictive pour être appliquée à des processus réels, il est essentiel d’inclure une étape de traitement des pertur- bations afin de fournir des estimations consistantes quelque soit le type de bruit agissant sur le procédé. Cette opération est réalisée en parallèle à l’estimation d’une base de l’espace d’observabilité en introduisant une variable instrumentale [14] particulière au sein du critère minimisé (cf. §IV-A). III. Estimation du vecteur d’observation par mise à jour d’une projection orthogonale Comme précisé dans la section précédente, la première étape de notre technique d’identification récursive consiste 1ui(t), wi(t), fi(t) et vi(t) sont définis de manière similaire à yi(t). en la mise à jour, à chaque nouvelle acquisition, de l’es- timée du vecteur d’observation zi. La solution proposée dans la suite de cette section s’appuie sur l’adaptation, au cas récursif, d’un problème similaire rencontré en identi- fication hors ligne des sous-espaces : éliminer de l’espace Imcol{Yt,i,j} le terme relatif à Imcol{Hi} : Zt,i,j = Yt,i,j − HiUt,i,j (8) où Yt,i,j ∈ Rli×j et Ut,i,j ∈ Rmi×j sont les matrices de Hankel des données d’entrée-sortie : Yt,i,j =      y(t) · · · y(t + j − 1) y(t + 1) · · · y(t + j) . . . ... . . . y(t + i − 1) · · · y(t + i + j − 2)      . (9) En effet, l’analogie entre ces deux situations est facilement mise en évidence en considérant les égalités suivantes2 : Zt,i,j+1 = Yt,i,j+1 − HiUt,i,j+1 = £ Yt,i,j yi(t̄ + 1) ¤ − Hi £ Ut,i,j ui(t̄ + 1) ¤ = £ Yt,i,j − HiUt,i,j yi(t̄ + 1) − Hiui(t̄ + 1) ¤ = £ Zt,i,j zi(t̄ + 1) ¤ . (10) La mise à jour de la matrice de Hankel Zt,i,j permet donc d’avoir accès au vecteur d’observation à l’instant t̄ + 1. Il serait donc intéressant d’ajuster la technique d’estimation hors ligne de Zt,i,j+1 à notre problème récursif afin d’ob- tenir zi(t̄ + 1). En identification des sous-espaces, dans le cas déterministe, le calcul de Zt,i,j+1 est réalisé à l’aide de la projection orthogonale sur le noyau de Ut,i,j+1 [2] : Zt,i,j+1 = Yt,i,j+1ΠU⊥ t,i,j+1 . (11) Or, dans la plupart des situations rencontrées en identifica- tion, ces matrices sont composées de données bruitées. Pour pallier le manque de consistance de l’estimateur précédent dans le cas stochastique, la stratégie suivie au sein de cet ar- ticle consiste à traiter l’effet des perturbations au cours de la seconde phase d’estimation. Nous pouvons donc considé- rer, au cours de cette première étape, que le problème envi- sagé est déterministe et ainsi analyser la projection (11), le biais d’estimation étant annulé lors du calcul de la matrice Γi. La méthode proposée jusqu’ici [4], [10], [16] consistait à mettre à jour la factorisation QR suivante3 : · Ut,i,j+1 Yt,i,j+1 ¸ = · R11(t̄ + 1) 0mi×li R21(t̄ + 1) R22(t̄ + 1) ¸ · Q1(t̄ + 1) Q2(t̄ + 1) ¸ (12) sachant que [2] : Yt,i,j+1ΠU⊥ t,i,j+1 = R22(t̄ + 1)Q2(t̄ + 1). (13) Cette technique, robuste numériquement, présente l’incon- vénient de nécessiter un nombre d’opérations assez consé- quent (cf. §III-B). La solution développée dans cet article repose sur la mise à jour directe de la projection orthogo- nale Yt,i,j+1ΠU⊥ t,i,j+1 . 2t̄ = t + i + j − 2. 3Il est démontré dans [10] que cette mise à jour fournit, au signe près, une estimation fiable du vecteur d’observation. A. Mise à jour directe de la projection orthogonale L’objectif est d’ajuster, à notre problème de mise à jour, une technique proposée par H. Oku [15]. Pour cela, consi- dérons qu’un nouveau couple de données est accessible à l’instant t̄ + 1 = t + i + j − 1. Les matrices de Hankel d’en- trée et de sortie sont alors modifiées de la façon suivante4 : Yj+1 z }| { Yt,i,j+1 = [λ Yj z }| { Yt,i,j yi z }| { yi(t̄ + 1)] (14) Ut,i,j+1 | {z } Uj+1 = [λ Ut,i,j | {z } Uj ui(t̄ + 1) | {z } ui ]. (15) L’idée est de mettre à jour directement la relation suivante : Yj+1ΠU⊥ j+1 = Yj+1 n I − UT j+1 ¡ Uj+1UT j+1 ¢−1 Uj+1 o (16) afin d’estimer le vecteur d’observation zi. Pour cela, calcu- lons, dans un premier temps, la projection orthogonale sur le noyau de Uj+1 étape par étape : ΠU⊥ j+1 = I − UT j+1 ¡ Uj+1UT j+1 ¢−1 Uj+1. (17) * A partir de (15) : Uj+1UT j+1 = λ2 UjUT j + uiuT i . (18) * En appliquant le lemme d’inversion matricielle [11] à (18), nous obtenons : ¡ Uj+1UT j+1 ¢−1 = ³¡ UjUT j ¢−1 − γiγT i λ2+δi ´ /λ2 (19) en posant : γi = ¡ UjUT j ¢−1 ui (20) δi = uT i γi. (21) * En multipliant (19) à droite par Uj+1 : ¡ Uj+1UT j+1 ¢−1 Uj+1 = ¡ Uj+1UT j+1 ¢−1 £ λUj ui ¤ h 1 λ n¡ UjUT j ¢−1 Uj − γiγT i λ2+δi Uj o γi λ2+δi i . (22) * En multipliant (22) à gauche par UT j+1 : UT j+1 ¡ Uj+1UT j+1 ¢−1 Uj+1 = " UT j ¡ UjUT j ¢−1 Uj − UT j γiγT i λ2+δi Uj λUT j γi λ2+δi λ γT i λ2+δi Uj δi λ2+δi # . (23) * La projection orthogonale sur le noyau de Uj+1 vaut : ΠU⊥ j+1 = Ij+1 − UT j+1(Uj+1UT j+1)−1 Uj+1 =     Π⊥ Uj z }| { Ij − UT j (UjUT j )−1 Uj +UT j γiγT i λ2+δi Uj −λUT j γi λ2+δi −λ γT i λ2+δi Uj λ2 λ2+δi     = · Π⊥ Uj 0 0 0 ¸ + 1 λ2 + δi · UT j γi −λ ¸ £ γT i Uj −λ ¤ . (24) 4λ est un facteur d’oubli facultatif permettant de pondérer l’effet des données passées. t, i et t̄ + 1 sont supprimées afin d’alléger les équations. * Finalement, nous obtenons : Yj+1ΠU⊥ j+1 = £ λYjΠ⊥ Uj 0 ¤ + 1 λ2 + δi £ λYjUT j γi −λyi ¤ £ γT i Uj −λ ¤ = h λYjΠ⊥ Uj − 1 √ λ2+δi žiγT i Uj λ √ λ2+δi ži i (25) avec : ži = λ √ λ2 + δi ¡ yi − YjUT j γi ¢ . (26) La relation entre ži et zi est démontrée en considérant le produit Zj+1ZT j+1. En effet, à partir de l’équation (11), nous avons : Zj+1ZT j+1 = Yj+1Π⊥ Uj+1 ³ Yj+1Π⊥ Uj+1 ´T = Yj+1Π⊥ Uj+1 YT j+1. (27) Or, en associant les équations (14) et (25), il est facile de montrer que : Yj+1Π⊥ Uj+1 YT j+1 = λ2 YjΠ⊥ Uj YT j + žižT i . (28) Ainsi : Zj+1ZT j+1 = λ2 YjΠ⊥ Uj YT j + žižT i . (29) De plus, nous savons que : Zj+1ZT j+1 = £ λZj zi ¤ · λZT j zT i ¸ = λ2 ZjZT j + zizT i . (30) Puisque YjΠ⊥ Uj YT j = ZjZT j , la relation suivante est alors évidente : ži = ±zi. (31) La mise à jour du vecteur zi est ainsi réalisable en appli- quant l’algorithme récursif (nommé IM) suivant : γi = ¡ UjUT j ¢−1 ui (32a) δi = uT i γi (32b) ži = λ √ λ2 + δi ¡ yi − YjUT j γi ¢ (32c) Yj+1UT j+1 = λ2 YjUT j + yiuT i (32d) ¡ Uj+1UT j+1 ¢−1 = ³¡ UjUT j ¢−1 − γiγT i λ2+δi ´ /λ2 . (32e) Il est important de noter que les dimensions des matrices composant l’algorithme (32) n’évoluent pas au cours de la récursion. B. Complexité Il est intéressant de remarquer que la technique proposée précédemment présente certains avantages calculatoires en comparaison avec le coût numérique de la méthode basée sur la mise à jour de la factorisation QR (12) [10], [16]. Le nombre de multiplications nécessaires est réduit ainsi que la capacité de stockage mémoire utile (cf. tab. I). C’est l’intérêt majeur de l’algorithme IM. Complexité Mémoire QR O(4mi2 (l + m 2 )) (m + l)2 i2 IM O(2mi2 (l + m)) mi2 (m + l) TABLE I Compléxité et espace mémoire de deux techniques d’estimation du vecteur d’observation. IV. Algorithme robuste pour l’estimation récursive de la matrice d’observabilité Nous venons de voir qu’il est possible d’estimer le vecteur d’observation à chaque instant à l’aide d’un algorithme de faible coût calculatoire. La seconde étape consiste donc à estimer récursivement une base de l’espace d’observabilité, une estimation consistante des matrices d’état du système étant réalisable à partir de cette dernière [2]. Plus précisé- ment, il s’agit de développer un algorithme de minimisation plus robuste que celui proposé jusqu’ici [10] afin d’assu- rer la fiabilité des estimées. Il est malgré tout nécessaire de faire quelques rappels sur la technique du propagateur appliquée en identification afin de mettre en évidence les lacunes possibles de l’algorithme EIVPM [9], [10]. A. Algorithme EIVPM Supposons que {A, C} soit observable. On définit le pro- pagateur P ∈ Rli×n comme l’unique opérateur tel que : Γi = · Γi1 Γi2 ¸ = · Γi1 PT Γi1 ¸ = · In PT ¸ Γi1 (33) avec Γi1 ∈ Rn×n une sous-matrice de Γi contenant n lignes linéairement indépendantes de Γi et Γi2 ∈ Rli−n×n la sous- matrice complémentaire. Il est facile de montrer que : spancol{Γi} = spancol ½ In PT ¾ . (34) Ainsi, en supposant que l’ordre est a priori connu, une es- timée du sous-espace engendré par les colonnes de Γi est accessible via l’estimation de P. Deux algorithmes d’esti- mation du propagateur ont été proposés dans [7]. Le plus efficace est l’algorithme EIVPM dont les formules de mise à jour sont les suivantes : g(t̄ + 1) = £ Rzi2 ξ(t̄)ξ(t̄ + 1) zi2 (t̄ + 1) ¤ (35a) Λ(t̄ + 1) = · −ξT (t̄ + 1)ξ(t̄ + 1) λ λ 0 ¸ (35b) Ψ(t̄ + 1) = £ Rzi1 ξ(t̄)ξ(t̄ + 1) zi1 (t̄ + 1) ¤ (35c) K(t̄ + 1) = ¡ Λ(t̄ + 1) + ΨT (t̄ + 1) M(t̄)Ψ(t̄ + 1)) −1 ΨT (t̄ + 1)M(t̄) (35d) PT (t̄ + 1) = PT (t̄) + ¡ g(t̄ + 1) − PT (t̄)Ψ(t̄ + 1) ¢ K(t̄ + 1) (35e) Rzi1 ξ(t̄ + 1) = λRzi1 ξ(t̄) + zi1 (t̄ + 1)ξT (t̄ + 1) (35f) Rzi2 ξ(t̄ + 1) = λRzi2 ξ(t̄) + zi2 (t̄ + 1)ξT (t̄ + 1) (35g) M(t̄ + 1) = 1 λ2 (M(t̄) − M(t̄)Ψ(t̄ + 1)K(t̄ + 1)) (35h) avec M(t̄) = ³ Rzi1 ξ(t̄)RT zi1 ξ(t̄) ´−1 . (35i) Il découle de la minimisation du critère suivant : J(P̂) = t X k=1 λt−k kzi2 (k)ξT (k) − P̂T zi1 (k)ξT (k)k 2 F (36) pour lequel λ est un facteur d’oubli permettant de pondé- rer les données passées, ξ ∈ Rγ×1 (γ > n) une variable instrumentale supposée non corrélée avec le bruit mais suf- fisamment corrélée avec l’état et {zi1 , zi2 } des sous-vecteurs de zi correspondant au découpage {Γi1 , Γi2 } de Γi [7]. B. Amélioration de EIVPM : algorithme EIVsqrtPM Comme tout algorithme de type moindres carrés récur- sifs, une condition nécessaire de stabilité numérique est d’imposer, à chaque itération, que la matrice M soit semi- définie positive. De par sa formule de mise à jour (cf. (35h)), cette condition ne peut être vérifiée. En effet, même si la matrice M est semi-définie positive à l’instant t̄, l’expres- sion : M(t̄ + 1) = 1 λ2 ³ M(t̄) − M(t̄)Ψ(t̄ + 1) (Λ(t̄ + 1) +ΨT (t̄ + 1)M(t̄)Ψ(t̄ + 1) ¢−1 ΨT (t̄ + 1)M(t̄) ´ (37) ne respecte pas cette propriété à l’instant t̄ + 1 puisque : Λ(t̄ + 1) + ΨT (t̄ + 1)M(t̄)Ψ(t̄ + 1) = " −ξT (t̄ + 1) ³ I − RT zi1 ξ(t̄)M(t̄)Rzi1 ξ(t̄) ´ ξ(t̄ + 1) λ + zT i1 (t̄ + 1)M(t̄)Rzi1 ξ(t̄)ξ(t̄ + 1) λ + ξT (t̄ + 1)RT zi1 ξ(t̄)M(t̄)zi1 (t̄ + 1) zT i1 (t̄ + 1)M(t̄)zi1 (t̄ + 1) ¸ (38) possède deux valeurs propres de signe opposé, I − RT zi1 ξ(t̄)M(t̄)Rzi1 ξ(t̄) étant une matrice semi-définie po- sitive. Il est alors impossible de garantir une bonne conver- gence de EIVPM (ou tout autre algorithme utilisant la technique EIV [14]). La solution proposée jusqu’ici reposait sur l’évaluation de M à partir de l’équation (35h), de l’ex- traction de la partie triangulaire supérieure puis de l’utili- sation de sa propriété de symétrie pour compléter sa mise à jour [3]. La technique développée au sein de cette section est plus fiable puisqu’elle utilise des outils mathématiques ro- bustes. Elle consiste en l’adaptation, au cas multivariable, d’un algorithme de mise à jour de la méthode de la variable instrumentale étendue proposée par B. Porat et B. Fried- lander [17]. L’idée sous-jacente est de calculer, à chaque itération, la racine carrée de M, condition suffisante à la vérification de la semi-définie positivité de M. Supposons donc que la matrice M soit semi-définie posi- tive à l’instant t̄. Si Λ était également semi-définie positive, il serait possible d’écrire M(t̄ + 1) de la façon suivante : M(t̄ + 1) = 1 λ2 ³ M1/2 (t̄)MT/2 (t̄)− M1/2 (t̄) ³ ΨT (t̄ + 1)M1/2 (t̄) ´T µ Λ1/2 (t̄ + 1)ΛT/2 (t̄ + 1) +ΨT (t̄ + 1)M1/2 (t̄) ³ ΨT (t̄ + 1)M1/2 (t̄) ´T ¶−1 ΨT (t̄ + 1)M1/2 (t̄)MT/2 (t̄) ´ . (39) La mise à jour des différents éléments composant cette ex- pression serait alors accessible en considérant la matrice suivante : W(t̄ + 1) = · Λ1/2 (t̄ + 1) ΨT (t̄ + 1)M1/2 (t̄) 0 M1/2 (t̄) ¸ (40) puisque : W(t̄ + 1)WT (t̄ + 1) = · Λ(t̄ + 1) + ΨT (t̄ + 1)M(t̄)Ψ(t̄ + 1) ΨT (t̄ + 1)M(t̄) M(t̄)Ψ(t̄ + 1) M(t̄) ¸ (41) contient les composantes de M(t̄ + 1). Malheureusement, Λ (cf. equ. (35b)) n’est pas semi- définie positive. Elle peut cependant être réécrite de la fa- çon suivante : Λ(t̄ + 1) = "p ξT (t̄ + 1)ξ(t̄ + 1) 0 − λ √ ξT (t̄+1)ξ(t̄+1) λ √ ξT (t̄+1)ξ(t̄+1) # · −1 0 0 1 ¸   p ξT (t̄ + 1)ξ(t̄ + 1) − λ √ ξT (t̄+1)ξ(t̄+1) 0 λ √ ξT (t̄+1)ξ(t̄+1)   = ∆(t̄ + 1)J∆T (t̄ + 1) (42) avec J matrice de signature de Λ [11]. Ainsi, en posant : W̄(t̄ + 1) = · ∆(t̄ + 1) ΨT (t̄ + 1)M1/2 (t̄) 0 M1/2 (t̄) ¸ , (43) l’expression (41) est modifiable comme suit : W̄(t̄ + 1) · J 0 0 In ¸ W̄T (t̄ + 1) = · Λ(t̄ + 1) + ΨT (t̄ + 1)M(t̄)Ψ(t̄ + 1) ΨT (t̄ + 1)M(t̄) M(t̄)Ψ(t̄ + 1) M(t̄) ¸ . (44) L’étape suivante consiste donc à mettre à jour la matrice W̄(t̄ + 1). En se référant à [11], nous savons qu’il existe une matrice RotH ∈ R2×2 composée de rotations hyper- boliques telle que : (RotH)J(RotH)T = J et ∆RotH triangulaire inférieure. (45) De même, pour toute rotation de Givens RotG [11] : (RotG)I(RotG)T = I. (46) La méthode adaptée de [17] consiste à associer ces deux types de rotations afin de les appliquer à la matrice W̄(t̄+1) de telle sorte que leurs actions cumulées fournissent une matrice triangulaire : W̄(t̄ + 1) · RotH 0 0 RotG ¸ | {z } Rot = · L11 0 L21 L22 ¸ (47) avec L11 et L22 deux matrices triangulaires inférieures de dimension respective 2 × 2 et n × n. Ainsi : W̄(t̄ + 1) · J 0 0 In ¸ W̄T (t̄ + 1) = W̄(t̄ + 1)Rot · J 0 0 In ¸ RotT W̄T (t̄ + 1) = · L11 0 L21 L22 ¸ · J 0 0 In ¸ · LT 11 LT 21 0 LT 22 ¸ . (48) Les égalités suivantes sont alors vérifiées : ∆(t̄ + 1)J∆T (t̄ + 1) + ΨT (t̄ + 1)M(t̄)Ψ(t̄ + 1) = L11JLT 11 (49) M(t̄)Ψ(t̄ + 1) = L21JLT 11 (50) M(t̄) = L21JLT 21 + L22LT 22. (51) A partir, de l’expression (50) : L21 = M(t̄)Ψ(t̄ + 1)L−T 11 J−1 . (52) Ainsi, avec (49), (52) et (35d), l’équation (51) devient : L22LT 22 = M(t̄) − L21JLT 21 = M(t̄)− M(t̄)Ψ(t̄ + 1) ¡ L11JLT 11 ¢−1 ΨT (t̄ + 1)MT (t̄) = M(t̄) − M(t̄)Ψ(t̄ + 1)K(t̄ + 1) = λ2 M(t̄ + 1). (53) Donc : M1/2 (t̄ + 1) = 1 λ L22. (54) De plus : K(t̄ + 1) = ¡ L11JLT 11 ¢−1 ΨT (t̄ + 1)MT (t̄) = ¡ L11JLT 11 ¢−1 L11JLT 21 = L−T 11 LT 21 (55) donc : K(t̄ + 1) = ¡ L21L−1 11 ¢T . (56) L’algorithme complet devient alors : g(t̄ + 1) = £ Rzi2 ξ(t̄)ξ(t̄ + 1) zi2 (t̄ + 1) ¤ (57a) ∆(t̄ + 1) = "p ξT (t̄ + 1)ξ(t̄ + 1) 0 − λ √ ξT (t̄+1)ξ(t̄+1) λ √ ξT (t̄+1)ξ(t̄+1) # (57b) Ψ(t̄ + 1) = £ Rzi1 ξ(t̄)ξ(t̄ + 1) zi1 (t̄ + 1) ¤ (57c) Trouver la matrice de rotations Rot telle que : · ∆(t̄ + 1) ΨT (t̄ + 1)M1/2 (t̄) 0 M1/2 (t̄) ¸ Rot = · L11 0 L21 L22 ¸ (57d) M1/2 (t̄ + 1) = 1 λ L22 (57e) K(t̄ + 1) = ¡ L21L−1 11 ¢T (57f) PT (t̄ + 1) = PT (t̄) + ¡ g(t̄ + 1) − PT (t̄)Ψ(t̄ + 1) ¢ K(t̄ + 1) (57g) Rzi1 ξ(t̄ + 1) = λRzi1 ξ(t̄) + zi1 (t̄ + 1)ξT (t̄ + 1) (57h) Rzi2 ξ(t̄ + 1) = λRzi2 ξ(t̄) + zi2 (t̄ + 1)ξT (t̄ + 1). (57i) Il sera désormais nommé EIVsqrtPM. V. Application L’objectif de cette simulation est de mettre en évidence la robustesse numérique de l’algorithme EIVsqrtPM (57) et de comparer le coût calculatoire de l’algorithme IM (32) avec celui de la mise à jour de la factorisation QR [10]. Pour cela, considérons le système suivant : ½ x(t + 1) = ¡ 0.9 0.5 0 0.3 ¢ x(t) + ¡ 0.5 −0.5 ¢ ũ(t) y(t) = ( −0.3 0.6 ) x(t) + v(t) (58) pour lequel ũ et v sont deux bruits blancs gaussiens de moyenne nulle et de variance respective 1 et 0.1. Dans un premier temps, les deux algorithmes étudiés utilisent l’algo- rithme (32) pour estimer le vecteur d’observation. La figure 0 100 200 300 400 500 600 700 800 900 1000 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Valeurs propres théoriques Valeurs propres avec EIVPM Valeurs propres avec EIVsqrtPM Echantillons Valeurs propres Fig. 1. Comparaison des algorithmes EIVPM et EIVsqrtPM. 1 illustre la fiabilité de l’algorithme EIVsqrtPM en compa- raison avec EIVPM. Sur 1000 appels de fonctions, le condi- tionnement de la matrice M est en moyenne de l’ordre de 150. On peut remarquer que la technique EIVPM présente des difficultés à estimer les valeurs propres du système alors que EIVsqrtPM fournit des estimées non biaisées. De plus, le tracé des valeurs propres obtenues avec EIVPM est beau- coup plus chaotique que celui produit avec EIVsqrtPM. En revanche, EIVsqrtPM présente l’inconvénient d’augmenter le temps de calcul par comparaison avec EIVPM. En ef- fet, sur un Pentium IV 2.4GHz, l’algorithme (35) est, en moyenne, 3 fois plus rapide que l’algorithme (57) pour réa- liser le même calcul. Ainsi, une étude du conditionnement de la matrice M est justifiée pour choisir parmi ces deux algorithmes, EIVPM ayant montré de bonnes dispositions sur des systèmes à matrice M bien conditionnée [10]. Méthodes Sub IM QR Durée adimensionnée 1 0.7 1.1 TABLE II Comparaison du temps de calcul des trois méthodes. La seconde partie de simulation concerne le temps de calcul accordé à l’estimation du vecteur d’observation. Sur cet exemple à entrée bien conditionnée (bruit blanc), toutes les techniques fournissent des estimées consistantes. Le ta- bleau II présente la comparaison des trois algorithmes sur 1000 appels de fonction : la méthode de soustraction (Sub) [7], la méthode QR [10] et l’algorithme IM. Ces résultats confirment les avantages numériques de l’algorithme (32). VI. Conclusion Le problème de l’identification récursive de procédés mo- délisés sous forme d’état est considéré au sein de cet article. La solution proposée se compose de deux étapes : – La mise à jour en ligne d’une projection orthogonale à l’aide du lemme d’inversion matricielle aboutissant à l’estimation fiable d’un vecteur concentrant les infor- mations contenues dans les données d’entrée-sortie. – L’estimation d’une base de l’espace d’observabilité à partir d’une mise à jour récursive d’un opérateur li- néaire nommé propagateur, estimation réalisée à l’aide d’un algorithme robuste de minimisation imposant la semi-définie positivité de certaines matrices. Les performances de ces algorithmes sont mises en évidence à l’aide d’une simulation. Références [1] M. Verhaegen et P. Dewilde. Subspace model identification : part 1 and part 2. International Journal of Control, 56 :1187–1241, 1992. [2] M. Verhaegen. Identification of the deterministic part of MIMO state space models given in innovations form from input output data. Automatica, 30 :61–74, 1994. [3] T. Gustafsson. Instrumental variable subspace tracking using projection approximation. IEEE Transactions on Signal Pro- cessing, 46 :669–681, 1998. [4] M. Lovera, T. Gustafsson, et M. Verhaegen. Recursive subspace identification of linear and non linear Wiener state-space models. Automatica, 36 :1639–1650, 2000. [5] H. Oku et H. Kimura. Recursive 4SID algorithms using gradient type subspace tracking. Automatica, 38 :1035–1043, 2002. [6] M. Lovera. Recursive subspace identification based on projector tracking. The 13th IFAC Symposium on System Identification, August 2003. [7] G. Mercère, S. Lecœuche, et C. Vasseur. Output noise insensitive identification of systems : a new recursive state-space method. The 3rd workshop on Physics in Signal and Image Processing (PSIP), Grenoble, France, January 2003. [8] J. Munier et G. Y. Delisle. Spatial analysis using new proper- ties of the cross spectral matrix. IEEE Transactions on Signal Processing, 39 :746–749, 1991. [9] G. Mercère, S. Lecœuche, et C. Vasseur. A new recursive method for subspace identification of noisy systems : EIVPM. The 13th IFAC Symposium on System Identification, Rotterdam, The Ne- therlands, August 2003. [10] G. Mercère et S. Lecœuche. Identification en ligne de sys- tèmes bruités : une nouvelle méthode récursive des sous-espaces. Journées Doctorales d’Automatique, Valenciennes, France, June 2003. [11] G. H. Golub et C. F. Van Loan. Matrix computations. John Hopkins University Press, Baltimore MD, 3rd edition edition, 1996. [12] H. Krim et M. Viberg. Two decades of array signal processing research : the parametric approach. IEEE Signal Processing Magazine, 13 :67–94, 1996. [13] S. Marcos, A. Marsal, et M. Benidir. The propagator method for source bearing estimation. Signal Processing, 42 :121–138, 1995. [14] T. Söderström et P. Stoica. System identification. Prentice Hall International Series in Systems and Control Engineering, New York, 1989. [15] H. Oku et H. Kimura. A recursive 4SID from the input-output point of view. Asian Journal of control, 1 :258–269, 1999. [16] M. Verhaegen et E. Deprettere. A fast, recursive MIMO state space model identification algorithm. The 30th IEEE Conference on Decision and Control, Brighton,England, December 1991. [17] B. Porat et B. Friedlander. The square-root overdetermined re- cursive instrumental variable algorithm. IEEE Transactions on Automatic Control, 34 :656–658, 1989.