IDENTIFICATION DES SYST ÈMES

IDENTIFICATIONDES SYSTÈMES Akram Ghorayeb

y1 u1 Entrées ou actions Système Sorties ou réponses up yq w Vitesse de rotation v Vitesse de translation q Débits I. INTRODUCTION 1. DÉFINITIONS • Système u = (u1, …, up)’ y = (y1, …, yq)’ Exemple(missile) • Modèle entrée-sortie

i C v R G j F a mg A q Exemple(accéléromètre) y k Entrée: u = accélération du véhicule m c u Sortie: y = déplacement relatif de m avec h fonction linéaire (1) est un modèle linéaire Exemple(missile dans un plan vertical) Entrée: u = (F, q) = force de propulsion Résistance de l’air: R =  (a.|vx|vx, b.|vy|vy)’ + e (e = effet du vent = perturbation) Sortie: y = (vx, vy, a) = vitesse et inclinaison Mouvement selon i (2) Modèle de la forme: Mouvement selon j (3) Mouvement de rotation Non-linéaire perturbé (4)

Equations d’état d’un modèle Elles sont de la forme x: vecteur d’état, y: sortie, u: entrée, w et v: perturbations Exemple (missile dans un plan vertical) Supposons que y = a (on s’intéresse seulement à l’inclinaison) En posant le modèle (2), (3), (4) s’écrit sous la forme Le modèle est linéaire ssi f et g sont des applications linéaires Dans ce cas,

Exemple (accéléromètre) Equation entrée-sortie: En posant (1) sera équivalente aux équations d’état: 2. LINÉARISATION Soit Si x est voisin de x0 Taylor d’ordre 1  (5) ou avec Si chaque variable d’un modèle reste voisine d’une valeur fixe, on peut linéariser ce modèle en appliquant (5) à chaque terme de ses équations.

Exemple (missile dans un plan vertical) On désire un vol stationnaire d’inclinaison a0 >0 et de vitesse selon i vx0 > 0 Déterminons d’abord la propulsion F0 et son angle q0 correspondants (2) Modèle: (3) (4) En régime stationnaire, les dérivées sont nulles. D’où, dans l’absence des perturbations (pas de vent), (2), (3) et (4) impliquent (vx0, vy0, a0, F0, q0) = point de fonctionnement

A cause des perturbations, la réponse y s’écarte de y0 = (vx0, vy0, a0)’ Un régulateur écarte l’entrée de (F0, q0) de sorte à ramener y vers y0 Ainsi, (vx, vy, a, F, q) fluctue autour de (vx0, vy0, a0, F0, q0) et les dérivées des variables fluctuent autour de 0. vx et vy étant voisins de vx0 et vy0, sont positives d’où |vx|vx = vx2 et |vy|vy = vy2 En appliquant (5) à chaque terme des équations (2), (3) et (4), on obtient: Ce sont des équations linéaires reliant les fluctuations des entrées aux fluctuations des sorties • Propriété essentielle des systèmes linéaires La réponse y à une combinaison linéaire d’entrées u1, …, uk est égale à la même combinaison linéaire des réponses y1, …, yk à ces entrées prises séparément. Cette propriété n’est pas vraies pour des systèmes non linéaires.

3. DISCRETISATION Pour une période T suffisamment petite, (9) (10) (11) La discrétisation d’un modèle entrée-sortie ou d’état s’obtient en remplaçant les dérivées par les approximations précédentes. Dans la suite on choisit T comme l’unité du temps (T = 1). Exemple Modèle continu : Pour T = 1, (9) et (10) impliquent D’où le modèle discrétisé (12)

y(t) (ut, yt-1) v(t) j g[, q] yt-1 ut Espace de la sortie Espace des données Espace du régresseur ( : espérance mathématique) On admettra que [v(t)] = 0 => Prédiction de y(t) [y(t)] = g[j(t), q] = En général, pour un système 1/1 (à une entrée et une sortie), posons ut = [u(0), …, u(t)]’, yt = [y(0), …, y(t)]’ et q = [q1, …, qm]’ un vecteur de paramètres • Forme générale d’un modèle discrétisé 1/1 y(t) = g[j(t), q] + v(t) avec j(t) = j(ut, yt-1)  k et g(, q) : k   v(t) : effet aléatoire des perturbations sur la sortie car si [v(t)] = m  0, on fait : v(t)  v(t) – m et g(,q)  g(, q) + m

Cas d’un système p/q En posant ut= [u1(0), …, up(0), …, u1(t), …, up(t)]’ yt= [y1(0), …, yq(0), …, y1(t), …, yq(t)]’ qj = [qj1, …, qjmj]’ et jj(t) = jj(ut,yt-1)  k yj(t) = gj[jj(t), qj] + vj(t) on a : j = 1, …, q q équations de même forme que 1/1 4. OBJECTIF DE L’IDENTIFICATION • Critique de la modélisation analytique • Avantage : • - Explique physiquement le fonctionnement du système. • Inconvénients : • Nécessite des hypothèses simplificatrices. • Le système peut être physiquement complexe • Certaines lois physiques régissant le fonctionnement du système peuvent être mal connues. Dans tous les cas l’expérience est nécessaire pour valider le modèle et déterminer ses paramètres.

L’identification d’un système est la détermination d’un modèle directement à partir de l’expérience sans le recours obligatoire à une analyse physique. Il s’agit de déterminer des fonctions j(t) = j(ut, yt-1) et g(, q) telles que où y(t) = la sortie mesurée à l’instant t • Plan du cours • Systèmes linéaires • - Fonction de transfert • - Réponse à un processus aléatoire • - Identification non paramétrique • - Modèles linéaires • - Observation de l’état (filtre de Kalman) • II. Identification paramétrique des systèmes linéaires • - Moindres carrés récursifs • - Méthodes récursives d’identification : • de la prédiction de l’erreur • du maximum de vraisemblance • de la variable instrumentale • - Validation du modèle • Identification des systèmes non linéaires • - Par gradient ajusté • - Par ondelettes • - Par logique flou

f(t) h(t) n t g(t) d(t) Syst. 0 0 II. SYSTEMES LINEAIRES 1. FONCTION DE TRANSFERT • Transformée-z d’une fonction • Pour f(t), t = 0, 1, 2,… Soit • Réponse impulsionnelle • = la réponse g(t) à d(t) d’un syst. • initialement inerte • Stabilité Système stable si

Fonction de transfert (FT) Car, pour t < 0, g(t) = 0. Pôle de G : est un point p de  tel que G(p) =  • Réponse d’un syst. linéaire à une entrée quelconque Pour u(t), t = 0, 1, 2, …., (u(t) = 0 si t < 0), on a par superposition (cf. fig) g(1)u(t-1) g(2)u(t-2) … t-2 t-3 t-1 t (* = produit de convolution) Remarque: comme pour n > t, u(t – n) = 0, Donc Zy = G(z).Zu

A|G(ejw)| = Réponse fréquentielle. = fonction de Bode A y wt -j u • Réponse fréquentielle Pour u(t) = ejwt • Equivalence entre équation et FT d’un système linéaire Forme de l’équation : Par transformée z : Si a0 = 0 et b0 0, y(t – 1) serait affecté par u(t). Impossible (principe de causalité)

Exemple sur Matlab >> n = [1 0 2]; >> d = [1 0.5 0.6]; >> sys = tf(n,d,1,'variable','z^-1') Transfer function: 1 + 2 z^-2 ----------------------- 1 + 0.5 z^-1 + 0.6 z^-2 Sampling time: 1 >> t = 1:20; >> u = 1 - exp(-2*t); >> y = lsim(sys,u); >> plot(t,[u', y]); >> impulse (sys) >> step(sys) >> bode(sys)

2. REPONSE A UN PROCESSUS ALEATOIRE • Rappels Déf. 1 : [u() processus aléatoire]  [ t, u(t) est une variable aléatoire]. Soit ft1,…tn() la densité de probabilité de [u(t1),…,u(tn)]. Déf. 2 : [u() est stationnaire]  [ft1,…tn() indép. de l’origine du temps]  Déf. 3 F[R(t)] : Tf. de Fourier Densité spectrale :  F1[F(w)] Pour u() échantillonné de période T,

Déf. 4 [u() est blanc]  [u() stationnaire et ft1,…,tn () = [f()]n]  et R(t) = s2.d(t) Déf. 5 [u() est ergodique]  [u() stat et ]  Pour u() échantillonné, les intégrales deviennent des sommes. Pour bruit blanc échantillonne (T=1 et t entier relatif)

u(t) y(t) G(z) • Propriétés de la réponse à une entrée stochastique 1) u(t) stationnaire => y(t) stationnaire car Comme la distribution de probabilité de u(t + k – n) est égale à celle de u(t – n), les distributions de y(t + k) et de y(t) sont égales. 2) u(t) ergodique => y(t) ergodique car a) b) Par un calcul semblable, on démontre que

3) my = G(1).mu 4) Autocorrélation (on supposera que mu = 0 (=> my = 0)) Dans le cas où u est blanc,

C1 z 1 5) Spectre de y En posant t – n2 + n1 = n, on obtient: (1) Si u est blanc, Fu(w) = s2. Mais

= s2 Si G(z) est stable, les pôles de G(z-1) sont extérieurs à C1 6)Corrélation entre u et y Si u est blanc, 7) Spectre couplé entre y et u

Donc et si u est blanc, • Modélisation d’une perturbation Si b(t) est un signal ergodique (par exple une perturbation), on peut le voir comme la réponse à un bruit blanc d’un système linéaire (1) vérifiant Noter que les coefficients c0 et d0 de zn et zm dans Gb(z) sont choisis égaux à 1 puisque le choix de s2 dans (1) est arbitraire. En désignant par li et mj les pôles et les zeros de Gb(z), on a:

g(n) 0 1 c 1 c 0 t = 1 t = -1 t = 0 • Exemples 1) Rb(t) et Fb(w) ? a) R’b(t) = Rb(t)/Rb(0) b)

2) a) b) ,t  0 et Vérification:

Analyse u(t) & G(z) y(t), R(t), F(w),.. 3. IDENTIFICATION NON PARAMETRIQUE Exemples précédents La suite du cours Identification G(z) u(t) & y(t) On sait que pour un système linéaire, G(z) est de la forme Si b0 = b1 = … = br-1 = 0, r < k, on peut écrire: Dans ce cas, le système a un retard pur de r périodes. Identification paramétrique : détermination des ai, bj, n, m et r (Chapitre suivant) Identification non paramétrique : détermination de g(t) et de G(ejw)

v + Syst. u + y • Réponse impulsionnelle : g(t) Soit l’équation générale d’un système linéaire d’entrée u et de sortie y perturbée par le signal v. On a : Or [u non corrélé avec v]  E[v(t).u(t – t)] = 0  t et si [u blanc]  Donc N grand Numériquement :

v u y G(z) F(z-1) F(z-1) yf uf Si l’entrée u(t) n’est pas un bruit blanc, [u coloré], • on détermine un filtre, dit de blanchiment, ayant la forme tel que e = bruit blanc c.à.d vérifiant Les paramètres fi de ce processus peuvent être déterminés par la méthode des moindres carrés qui sera donnée au chapitre suivant. • En multipliant les deux membres de l’équation du système, Zy = G(z).Zu + Zv, par F(z-1), on obtient: F(z-1).Zy = G(z). F(z-1).Zu + F(z-1).Zv ou uf, yf et vf sont les filtrages de u, y, et v par F avec uf blanc. 

Réponse fréquentielle L’équation du système peut s’écrire sous la forme (y0 ind de v) avec v(t) une perturbation et Spectre de y. Or, (1) Couplage entre u et y [v et u non corrélés]  [Rvu(t) = 0  t] (2)  (1) et (2) 

Val de y t=t N t=N-t Val de u wg 1 t -g +g Calcul numérique deFyu(w), Fy(w) et Fu(w). Quand Meilleure résolution Mais Si t grand  estimation de Ryu moins précise Remède : on donne moins d’importance (un poids plus faible) aux grandes valeurs de t. Fenêtre de Hamming (formules semblables pour Fu et Fy)

Calcul de la variance de G(ejw) Remarque : généralement, Fv(w)/Fu(w)  avec w Var[G(ejw)]  avec w • Exemple de MATLAB Taper « iddemo » pour apprendre à utiliser le Toolbox « system identification ». dryer2 = tableau de y2 et u2 >> load dryer2 >> dry = iddata(y2, u2, 0.08); >> ze = dry(100 : 300); >> ze = dtrend(ze); >> plot(ze) Toutes les données, 0.08 = T Données pour identification Pour supprimer les moyennes

>> impulse(ze, 'sd', 3, 'fill') Estimation de la rép. impul. Réponse à un échelon unité >> step(ze) >> gs = spa(ze); Analyse spectrale de ze Amplitude et phase de G(ejw) >> bode(gs)

e Modèle de la perturbation H(z) Système v u y0 y + G(z) + 4. MODELES LINEAIRES • Modèles entrée-sortie y0 = réponse propre y =réponse mesurée v = perturbation e = bruit blanc (bb) (1) En posant Forme générale d’une FT Les racines de C et de D sont de modules < 1 (stables) Modèle Box-Jenkens (B-J)

Autorégressif (AR) Moyenne mobile (MA) Exogène(X) bb Si F et D ont des racines communes, on pose : F(z) = A(z).F(z) et D(z) = A(z).D(z) La forme la plus générale et B-J  (2) • Modèles particuliers 1) Si F(z) = D(z) = 1, on a : (2)  (3) ou, en forme temporelle, (3)  Modèle ARMAX

Autorégressif (AR) Exogène(X) bb 2) Si F(z) = D(z) = C(z) = 1, on a : (4) (2)  ou, en temporelle, (4)  Modèle ARX 3) Si F(z) = D(z) = 1 et C(z) = A(z), on a : (2)  Sortie erronée ou, en temporelle, Modèle OE (output error) (5) 

Modèle d’état discrétisé En remplaçant dans 5. OBSERVATION DE L’ETAT FILTRE DE KALMAN On admet que R1: matrice de corrélation (n, n) entre les composantes de w. R2: matrice de corrélation (q, q) entre les composantes de v. R12 : matrice de corrélation (n, q) entre w et v.

J(t+1) = E[T(t+1).(t+1)] = E{[xT(t+1)  T(t+1)]. [x(t+1)  (t+1)]} = E[ ||x(t+1)||2  2.xT(t+1). (t+1) + T(t+1). (t+1)]. Forme de l’estimateur de x: F, G, K sont des matrices à déterminer. Pour cela on pose et on minimise par rapport aux éléments de F, G et K. Or, Mais car dF/dfij = (1ij)

(t+1) = A. (t) + B.u(t) + K(t).[y(t) - C. ] Equations d’orthogonalité E[(t+1). T(t)] = 0. E[(t+1).u T(t)] = 0, E[(t+1).y T(t)] = 0 3 équations matricielles à 3 inconnues F, G et K Maintenant, en posant: Matrice de corrélation de l’erreur (t) On démontre que Equations de Ricatti avec

IDENTIFICATION DES SYST ÈMES