Identification Des Systemes PDF

Identification des systèmes
Gonzalo Cabodevila
gonzalo.cabodevila@femto-st.fr
ème
3 année
Option Mécatronique
École Nationale Supérieure de

Mécanique et des Microtechniques
26, chemin de l’Épitaphe
25030 Besançon cedex – FRANCE
http://intranet-tice.ens2m.fr
3
" Je distingue deux moyens de cultiver les sciences : l'un

d'augmenter la masse des connaissances par des découvertes ;
et c'est ainsi qu'on mérite le nom d'inventeur ; l'autre de
rapprocher les découvertes et de les ordonner entre elles, an
que plus d'hommes soient éclairés, et que chacun participe,
selon sa portée, à la lumière de son siècle ... "
Diderot
4
Table des matières
I Identication des systèmes 9
1 Modèles de connaissance 11
1.1 Modèles simples linéaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.1 Exemple : le moteur à courant continu . . . . . . . . . . . . . . . . . . . . . . . 11
1.2 Méthodes systématiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.1 Un exemple simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Détermination des constantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4 Linéarisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Identication de modèles non paramétriques 15

2.1 Quelle entrée ? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.1 Analyse harmonique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.2 Réponse impulsionnelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.3 Réponse indicielle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Etude des réponses indicielles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.1 Systèmes du premier ordre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.2 Systèmes du second ordre résonnant . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.3 Systèmes du premier ordre retardés ou du second ordre non résonnant . . . . . 16
2.2.4 Systèmes d'ordre supérieur à 2 non résonnants . . . . . . . . . . . . . . . . . . 17
2.2.5 Système avec intégrateur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3 Méthode de Ziegler-Nichols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Autres modèles ? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Algorithme général d'identication 21

3.1 Choix du signal d'excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Les séquences binaires pseudo aléatoires (SBPA) . . . . . . . . . . . . . . . . . . . . . 22
3.2.1 Choix des paramètres d'une SBPA . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Identication basée sur l'erreur de sortie . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 Identication basée sur l'erreur de prédiction . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.1 Méthode des moindres carrés simples . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.2 Calcul du biais de l'estimateur . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4.3 Méthode des moindres carrés généralisés . . . . . . . . . . . . . . . . . . . . . . 29
3.4.4 Méthode de la Matrice Instrumentale . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5 Identication en boucle fermée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.6 Et Matlab
r dans tout ça ? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4 Estimations récursives 35
4.1 Moindres carrés récursifs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Variable instrumentale récursive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5
6 TABLE DES MATIÈRES
5 Algorithmes d'optimisation paramétrique. 39

5.1 Présentation générale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2 Exemple introductif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.3 Problème posé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.4 Revue bibliographique succincte. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.5 Eléments de comparaison des diérentes méthodes. . . . . . . . . . . . . . . . . . . . . 41
5.5.1 Critères. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.5.2 Comparaisons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.5.3 Récapitulatif. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.5.4 Synthèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6 Description des algorithmes utilisés. 45

6.1 Gradient et quasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.1.1 Méthode du gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.1.2 Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.1.3 Méthode de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.1.4 Méthode de Levenberg-Markardt . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.1.5 l'algorithme BFGS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.2 Simplex de Nelder & Mead. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.2.1 Extensions - Améliorations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.3 Algorithmes génétiques. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.3.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.3.2 Principe des algorithmes génétiques binaires. . . . . . . . . . . . . . . . . . . . 52
6.3.4 Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.4 Algorithmes génétiques codés réels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.4.1 Principe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.4.2 Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.5 Méthode du Recuit simulé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.5.2 Principe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.5.4 Mise en ÷uvre du recuit simulé. . . . . . . . . . . . . . . . . . . . . . . . . . . 56
Bibliographie 56
7 Travaux dirigés 59
7.1 Strecj - Broida - Programmation non linéaire . . . . . . . . . . . . . . . . . . . . . . . 59
7.1.1 Strecj - Broida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.1.2 Programmation non linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.1.3 Limites de la Programmation non linéaire . . . . . . . . . . . . . . . . . . . . . 59
7.2 Identication des systèmes instables ou intégrateurs . . . . . . . . . . . . . . . . . . . . 59
7.3 Choix de la source d'excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.4 Les modèles paramétriques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.5 Cas pratique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
II Annexes 63
Examen nal janvier 2004 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
Correction examen nal janvier 2004 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
TABLE DES MATIÈRES 7

8 TABLE DES MATIÈRES
Première partie
Identication des systèmes
9
Chapitre 1
Modèles de connaissance
1.1 Modèles simples linéaires
Ce sont les modèles issus de la physique. Lorsque le système est peu complexe il est possible d'écrire
les relations entre les diérentes grandeurs physiques décrivant les diérents sous-systèmes.
1.1.1 Exemple : le moteur à courant continu
R L
i
-
6 6'$
J, f
u E
&%
Figure 1.1 Schéma électromécanique d'un moteur à courant continu
avec : u = E + Ri + Ldi/dt
C = ki
E = kω
J dω
dt = C − fω
De l'ensemble des équations on peut tirer diérents graphes :

Schéma bloc
Graphe de uence (voir polycopié de cours de Bernard Lang)
Graphes à Liens ou Bond graph (voir gure 1.2) (Hors du cadre de ce cours : cf Techniques de
l'ingénieur et autres bons bouquins sur la question)
De ces graphes, par réduction, on obtient la fonction de transfert, préliminaire à la majorité des
méthodes d'analyse et de synthèse de correcteurs que vous connaissez.
La gure 1.3 rappelle les analogies entre les diérents domaines de la physique.
11
12 CHAPITRE 1. MODÈLES DE CONNAISSANCE
Figure 1.2 Système et son modèle en Bond Graph. (tiré de : http ://www.univ-
lille1.fr/l2ep/commande/ee/presentations/p-cs.htm)
Figure 1.3 Analogie des équations diérentielles linéaires.
1.2 Méthodes systématiques
Lorsque le système s'avère plus complexe, il est souvent utile de mettre en ÷uvre une méthode sys-
tématique. Les équations de Lagrange sont l'une de ses méthodes, particulièrement adaptée dans le
cadre de la robotique mais presque tous les systèmes sont modélisables par cette méthode.
Les équations de Lagrange s'écrivent sous la forme :

d ∂T ∂T ∂D ∂V
− + + = Qi
dt ˙
∂ qi ˙
∂qi ∂ qi ∂qi
avec : T = Σ énergies cinétiques

1
D =
2 Σ pertes dépendantes de la vitesse
V = Σ énergies potentielles
Q = Forces extérieures
1.2.1 Un exemple simple

L0
On suppose L(x) = x
1.2. MÉTHODES SYSTÉMATIQUES 13
i(t)
R - L0
e(t) 6 L(x) = x(t)

?x(t)

F = Mg
?
Figure 1.4 Sustentation magnétique.
Mécanique Electrique
Coordonnées q2 = x(t) q1 = q(t)
Vitesse q̇2 = ẋ(t) q̇1 = q̇(t) = i(t)
Travail des forces extérieures τ2 = M gx(t) τ1 = e(t).q(t) = e.q1
∂τ2 ∂τ1
Forces extérieures Q2 = ∂q 2
= M g Q 1 = ∂q1 = e
Energie cinétique T2 = 21 M q̇22 T1 = 12 L(x)q̇12 = 21 Lq20 q̇12
Energie potentielle V2 = 0 V1 = 0 (pas de capacité)
Fonction de dissipation D2 = 0 D1 = 12 Rq̇12
T = T2 + T1 = 12 M q̇22 + 21 Lq20 q̇12
V = V2 + V1 = 0
D = D2 + D1 = 12 Rq̇12
∂T L0
∂ q˙i
M q̇2 q2 q̇1
d ∂T q̈1 q2 −q̇1 q̇2
(
dt ∂ q˙i ) M q̈ 2 L 0 q2 2
∂D
∂ q̇2 0 Rq̇1
∂V
∂qi 0 0
∂T L0 q̇12
∂qi 2 (− q22 ) 0
Qi Mg e
d'où on en déduit les formules suivantes :
q̇12

L0
M q̈2 + = Mg
2 q22

q̈1 q2 − q̇1 q̇2
L0 + Rq̇1 = e
q22
q̇12

L0
q̈2 = g −
2M q22

q2 L0 q̇1 q̇2
q̈1 = e − Rq̇1 +
L0 q22
L'utilisation de logiciels de calcul analytique (Maple, Mathematica, ...) permettent une modélisation
rapide et sans erreurs.
14 CHAPITRE 1. MODÈLES DE CONNAISSANCE
1.3 Détermination des constantes
Evidemment, avant de pouvoir simuler le système, il faut déterminer la valeur numérique des diérents
paramètres entrant en jeu. La bonne méthode consiste alors à soumettre le système à une série de
mesures permettant d'"isoler" le paramètre que l'on désire mesurer. La littérature regorge de méthodes
permettant des mesures ables.
ATTENTION : l'utilisation de méthodes de type "identication paramétrique" (voir chapitre 3) ne
permet en aucun cas de mesurer ces paramètres !
1.4 Linéarisation
Ce système étant non linéaire et notre but étant d'en déterminer un modèle linéaire (fonction de
transfert, représentation d'état), la première étape consiste donc à le linéariser. Il existe d'innombrables
méthodes de linéarisation, la littérature vous fournira une méthode appropriée à votre problème. Je
n'en présenterai ici qu'une :
la linéarisation autour d'un point de fonctionnement.
Cela consiste en un développement limité au premier ordre (a = a0 + â) autour d'un point.
En reprenant l'exemple du paragraphe 1.2.1
q̈1 = F1 (q1 , q2 , q̇1 , q̇2 , e)
La linéarisation autour du point E est donnée par :

¨q̂ = ∂F1 q̂1 + ∂F1 q̂2 + ∂F1 q̂˙ + ∂F1 q̂˙ + ∂F1 ê

1 1 2
∂q1 E ∂q2 E ∂ q̇1 E ∂ q̇2 E ∂e E
Notez bien que cette étape de linéarisation n'est nécessaire que si vous comptez utiliser une méthode
de synthèse fondée sur une représentation linéaire !
Concernant la simulation : implantez les deux modèles (linéarisé et non-linéarisé), le premier vous
permettra de vérier votre synthèse de correcteur, le second (un peu plus proche de la réalité) vous
permettra de simuler la robustesse de votre correcteur vis-à-vis des non-linéarités du système.
1.5 Conclusions
Il faut toujours avoir un modèle de connaissances, si mauvais soit-il pour pouvoir "sentir" le procédé et
interpréter vos résultats. Ce modèle de connaissances est le préalable à toute identication de système
quelle que soit la méthode employée.
Chapitre 2
Identication de modèles non

paramétriques
Si le modèle connaissance s'avère trop peu précis ou impossible à obtenir, on peut alors se sur un
modèle de comportement entrée - sortie. On parle alors d' identication du système. L'objectif est
d'obtenir un modèle qui "se comporte comme" le système. La plupart du temps, les paramètres de ce
modèles n'ont qu'un rapport lointain avec les paramètres d'un modèle de connaissances.
En premier lieu, il faut donc répondre aux deux questions suivantes :
Quelle entrée ?
Comment traiter la sortie ?
2.1 Quelle entrée ?
2.1.1 Analyse harmonique
Entrée sinusoïdale de type u = A sin(ωt), ω balaye l'espace des pulsations susceptibles de contenir
une pulsation de coupure du système. En notant l'amplitude et le déphasage de la sortie vis-à-vis de
l'entrée on trace un diagramme de Bode. De l'analyse de ce diagramme on détermine le modèle
Les résultats sont diciles à exploiter si les constantes de temps sont proches. Bonne excitation sur
tout le spectre de fréquences. Ce n'est pas une commande industrielle classique et par conséquent elle
est dicile à mettre en oeuvre.
2.1.2 Réponse impulsionnelle
Idéalement la meilleure méthode car le spectre est constant. Mais il est impossible de réaliser un Dirac
correct, sauf peut-être en électronique...
2.1.3 Réponse indicielle
Le spectre est correct, la commande est facile à implanter car c'est une commande classique. C'est la
méthode la plus utilisée.
2.2 Etude des réponses indicielles
2.2.1 Systèmes du premier ordre
Pas de retard pur, démarrage franc (à comparer aux systèmes d'ordre 2 ou supérieur).
En cas de bruit, la méthode de mesure de τ via l'intersection de la pente de la tangente au démarrage
avec la valeur nale est meilleure que le relevé du temps à 63% de la valeur nale.
15
16 CHAPITRE 2. IDENTIFICATION DE MODÈLES NON PARAMÉTRIQUES
Table 2.1 Réponse d'un système du premier ordre à un échelon.
6s(t)
sf inal
63% sf inal
A
H(p) = 1+τ p
-
τ t
2.2.2 Systèmes du second ordre résonnant
Table 2.2 Réponse d'un système du deuxième ordre à un échelon.
D%
sf inal
A
- H(p) = p2
Tpseudo−periode 1+ 2m
w0 p+ w 2
0
-
tpic t
Le facteur d'amortissement m se détermine à l'aide de la mesure du dépassement :
√−πm
D% = 100e 1−m2
La pulsation propre w0 s'obtient via l'une des deux formules :
π
tpic = √
w0 1 − m2
ou
2π
Tpseudo-période = √
w0 1 − m2
2.2.3 Systèmes du premier ordre retardés ou du second ordre non résonnant

Ae−T p
La méthode de Broïda consiste à "faire coller" un modèle de la forme H(p) = 1+τ p sur la réponse
du système.
Les valeurs de T et de τ sont calculées à partir des relations suivantes :
2.2. ETUDE DES RÉPONSES INDICIELLES 17
Table 2.3 Réponse d'un système du deuxième ordre non résonnant à un échelon.
6 s(t)
sf inal
Ae−T p
40% sf inal H(p) = 1+τ p
28% sf inal
t1 t2 -
t
avec : τ = 5, 5 × (t2 − t1 )
T = 2, 8 × t1 − 1, 8 × t2
En pratique, la prise en compte du retard pur est dicile dans une synthèse de correcteur. On utilise
donc l'approximation suivante :
Ae−T p A(1 − T p) A
H(p) = ' '
1 + τp 1 + τp (1 + τ p)(1 + T p)
Cette approximation reste valable jusqu'à la pulsation ω = 1/T
2.2.4 Systèmes d'ordre supérieur à 2 non résonnants
La méthode la plus connue est la méthode de Strejc. Le modèle est :
Ae−T p
H(p) =
(1 + τ p)n
Tu
1. Mesurer
Ta.
Tu
2. Dans la colonne
T a du tableau 2.4 trouver la valeur immédiatement inférieure à ce ratio.
3. Sur la ligne de ce ratio déterminer n.
Ta
4. Toujours à l'aide des valeurs numériques de cette ligne, calculer τ avec
τ .
5. Calculer la nouvelle valeur de T u0 avec
Tu
τ .
6. En déduire T avec T = T u − T u0 .
Note :
Le tableau 2.4 vient du calcul de la réponse indicielle des fonctions de transfert de la forme
Ae−T p

−1
s(t) = L
p(1 + τ p)n
Pour n=4 t t t
te− τ 1 t2 e− τ 1 t3 e− τ

1 t
s(t) = L−1 = 1 − e− τ − − −
p(1 + τ p)4 τ 2 τ2 6 τ3
la dérivée seconde s'écrit :
t t
d2 s(t) 1 t2 e− τ 1 t3 e− τ
= − −
dt2 2 τ4 6 τ5
Table 2.4 Méthode de Strejc : Réponse indicielle.

0.18
0.16
n Tu/Ta Tu/τ Ta/τ

0.14
1 0 0 1
0.12 2 0.104 0.282 2.718
3 0.218 0.805 3.695
0.1
4 0.319 1.425 4.465
0.08 5 0.410 2.100 5.119
6 0.493 2.811 5.699
0.06
7 0.570 3.549 6.226
0.04
8 0.642 4.307 6.711
0.02 9 0.709 5.081 7.164

10 0.773 5.869 7.590
0
0 1 2 3 4 5 6 7 8
Tu Ta
qui s'annule pour t = 3τ .

la tangente de plus grande pente est donc la tangente à la courbe en t = 3τ . Le calcul de la droite
tangente donne :
e−3 53 −3
∆ = 9/2 t− e +1
τ 2
53
Cette droite ∆ intercepte la droite d'équation s = 1 en t = 9 τ et l'axe des abscisses en t =
τ (53 e−3 −2) τ (53 e−3 −2)
1/9 e−3
T a = 53
donc
9 τ − 1/9 e−3
et T u/T a = .3193573120
Note : si T u/T a < 0.1, donc n = 1, On applique plutôt la méthode de Broïda.
2.2.5 Système avec intégrateur
Table 2.5 Méthode de Strejc (système avec intégrateur)

1.4
1.2
n x/y
0.8 1 0.368
∆y
2 0.271
0.6 ∆t 3 0.224
4 0.195
0.4
5 0.175
Y
0.2
X
0
0 1 2 3 4 5 6 7 8
Tu
2.3. MÉTHODE DE ZIEGLER-NICHOLS 19
Identié sous la forme
ku 1
H(p) = ×
p (1 + τ p)n
1. Avec x/y : déterminer n

Tu
2. Calcul de τ avec τ= n
∆y 1
3. Calcul de ku = ∆u × ∆t
En pratique il est assez rare d'identier un système intégrateur en boucle ouverte car le système dérive
rapidement en dehors de son domaine de linéarité (saturation, ...). En pratique, on procède à une
identication en boucle fermée.
2.3 Méthode de Ziegler-Nichols
Tous comptes fait, le but principal de l'identication en automatique est la mise en place d'un correc-
teur ou d'un régulateur an d'améliorer les performances du système en boucle fermée. Aussi l'on est
en droit de se demander si cette étape d'identication est bien nécessaire...
C'est exactement le principe de la méthode de Ziegler-Nichols, qui n'est pas une méthode d'identica-
tion à proprement parler mais une méthode de réglage fondée sur des mesures eectuées directement
sur le système. Le tableau 2.6 résume la méthode.
2.4 Autres modèles ?
Les méthodes précédentes sont très performantes si toutefois le comportement du système est proche
des réponses indicielles attendues.
Qu'en est-il pour les systèmes décrits par les fonctions de transfert suivantes :
Ae−T p 1 − Tp
H1 (p) = p2
et H2 (p) =
1+ 2z (1 + τ1 p)(1 + τ2 p)(1 + τ − 2p)
w0 p + w0 2
Dans la plupart des cas, l'un des modèles précédents sera susant. N'oubliez pas que nous cherchons
un modèle de comportement uniquement. Dans les quelques cas résistant aux méthodes précédentes,
un peu d'astuce permet de trouver une méthode pour déterminer les constantes du système. Et si rien
de ce qui précède ne vous aide, Il reste les méthodes de modélisation paramétrique.
Table 2.6 Coecients d'un PID réglé par les méthodes de Ziegler-Nichols et Chien-Hrones-Reswick :
essai indiciel et méthode du pompage.
Méthode de pompage Méthodes apériodiques

Boucle fermée Boucle ouverte
Systèmes stables ou instables Systèmes stables, instables

en boucle ouverte ou intégrateurs
Ziegler-Nichols Ziegler-Nichols Chien-Hrones-Reswick
Régulation ou Régulation ou
Poursuite Poursuite
Régulation Poursuite
P K = 0.5Kosc K = Tτ K = 0.3 Tτ K = 0.3 Tτ
P.I K = 0.45Kosc K = 0.9 Tτ K = 0.6 Tτ K = 0.35 Tτ
Ti = 0.83Tosc Ti = 3.3τ Ti = 4τ Ti = 1.2T
P.I.D K = 0.6Kosc K = 1.2 Tτ K = 0.95 Tτ K = 0.6 Tτ
Ti = 0.5Tosc Ti = 2τ Ti = 2.4τ Ti = T
Td = 0.125Tosc Td = 0.5τ Td = 0.4τ Td = 0.5τ
Chapitre 3
Algorithme général d'identication

Il faut répondre aux questions suivantes :
Quel genre de modèle ?

Paramétrique ou non paramétrique ?
Quelles mesures eectuer ?
Quel signal d'excitation ?
Quelle structure adopter ?
Doit-on y inclure des connaissances a priori ?
Comment choisir le bon modèle parmi tous ceux calculés ?
Ce modèle est-il adapté à ce que l'on va lui demander ?

Connaissances a priori

? ? ? 6
Choix de la
Mesures Critère de
structure
sélection
?
Données
? ? ?
Calcul d'un
modèle
Validation
?
Modèle nal
Figure 3.1 Algorithme général d'identication des systèmes. Notez que le processus d'identication
est un processus itératif, chaque passage apportant un peu plus de connaissances.
21
22 CHAPITRE 3. ALGORITHME GÉNÉRAL D'IDENTIFICATION
b0 b1 b2 b3 b4 b5 b6 b7
- -
Figure 3.2 Principe de génération de nombre aléatoires fondé sur un registre à décalage.
Table 3.1 Tableau des bits à utiliser pour obtenir une séquence de longueur maximale
n bits
L
3 b2 L b1
4 b3 L b2
5 b4 L b2
6 b5 L b4
7
L b6 L b5 L
8 b7 b3 L b2 b1
9 b8 L b4
10 b9 b6
3.1 Choix du signal d'excitation
Pour bien identier, il faut bien exciter dans tout le spectre de fréquences susceptible de contenir des
constantes de temps du système.
sin(wt) : parfait d'un point de vue spectre (balayage en fréquence) mais peu de systèmes acceptent
ce genre d'entrées.
δ(t) : parfait du point de vue théorique, mais, hormis en électronique, il est impossible de réaliser
un tel signal.
u(t) : moins bon d'un point de vue spectral (u(f ) = sinc(f )), mais facile à implanter
b(t) : bruit blanc idéal d'un point de vue spectral mais comment le réaliser ?
Il arrive que vous n'ayez aucune possibilité d'exciter le système (ex : machine en production), il faudra
alors proter des commandes "naturelles" du système comme signal d'entrée du système. Dans ce cas,
le premier travail consiste à calculer le spectre du signal d'entrée (FFT par exemple). Il faudra vérier
a posteriori que les constantes de temps identiées sont bien dans des domaines de fréquences dans
lesquels le système a été excité.
3.2 Les séquences binaires pseudo aléatoires (SBPA)
L'un des moyens de réaliser un signal "aléatoire" est la mise en ÷uvre de Séquences Binaires Pseudo
Aléatoires (en anglais PRBS : Pseudo Random Binary Sequence)
Principe :
Pour que le système fonctionne, on doit initialiser le registre à n'importe quelle valeur binaire sauf
zéro. le nombre codé en binaire A = b7 b6 b5 b4 b3 b2 b1 b0 est "aléatoire. de la même façon, le bit b7 semble
1
indépendant de ses valeurs précédentes.
Une séquence sur N bits à une longueur de 2N −1 = L dont 2N −1 1et 2N −1 −1 0, sa valeur moyenne
est donc non nulle.
a
E(s(t)) =
L
1. Cela ne semble pas très sérieux, mais sachez que les générateurs de bruit en électronique utilisent ce principe mais
sur 40 bits et que le cryptage de textes utilise un principe tout à fait similaire sur 128 bits !
3.2. LES SÉQUENCES BINAIRES PSEUDO ALÉATOIRES (SBPA) 23
amplitude
6
"1"
a
-
t
"0"
Figure 3.3 diagramme temporel d'une séquence binaire pseudo aléatoire.
a2
−Te Te LTe t-
2
− aL -
Figure 3.4 Fonction d'auto-corrélation d'une SBPA.
Pour déterminer la forme du spectre il faut en premier lieu calculer sa fonction d'auto-corrélation, le
spectre étant la transformée de Fourier de la fonction d'auto-corrélation
Fonction d'auto-corrélation
Z LT
1
Cxx (τ ) = x(τ )x(t − τ )dt
LT 0
si |τ | ≤ T
L + 1 |τ |
Cxx (τ ) = a2 (1 − . )
L T
si |τ | ≥ T
a2
Cxx (τ ) = −
L
On en déduit le spectre qui est la transformée de Fourier de la fonction d'auto-corrélation.
+∞
2L
+1X n sin π Ln 2 a2
P (f ) = a δ(v − )( ) − δ(v)
L2 −∞ LT π Ln L
Le spectre de la SBPA représenté en gure 3.5 est donc un peigne de Dirac modulé par un sinus
cardinal. On considère que le bruit peut être considéré comme un bruit blanc dans le premier tiers du
premier lobe (±1/3Te ). On remarque donc que :
plus L grand ⇒ plus la valeur moyenne est faible
plus L grand ⇒ plus de raies, plus proche du bruit blanc
SBP A(f )
6
6
6
6
6
6 6
6 -f
1
LTe 2 1/Te 2/Te
LTe
Figure 3.5 Spectre d'une SBPA.
3.2.1 Choix des paramètres d'une SBPA
Il faut donc déterminer a, L etT e "optimum" vis-à-vis du système à identier.

Le premier compromis concerne L et Te .
Plus la séquence est longue :

plus il y a d'informations
plus le transitoire est négligeable
plus les organes sourent
plus les dérives s'accentuent
Le deuxième compromis est : plus a est petit,

plus le signal est noyé dans le bruit (le vrai celui là !)
moins les organes sourent
moins les dérives s'accentuent
Le choix de a est donc avant tout lié au bruit présent dans le système. On choisira a juste assez grand
pour avoir du signal en sortie du système.
Choix de L et Te Une SBPA de longueur L = 2N − 1 envoyée à la fréquence 1/Te présente un

plateau de longueur maximale de longueur N × Te . Le plus petit plateau étant de longueur Te
Le choix de N (donc L) et Te est donc, encore une fois, un compromis entre :
une bonne identication du gain statique
une bonne excitation sur la bande de fréquence du système
Soit τmax la plus grande constante de temps du système et τmin la plus petite constante de temps du
système.
Pour une bonne identication du gain statique on choisira :
NTe = 3 à 5 τmax
Pour une bonne excitation sur le spectre entre 0 et la plus haute fréquence de coupure on choisira :
1 1
0, 3 =
Te 2πτmin
La résolution de ses deux équations donne N et Te
Note 1 : si τmax τmin N devient vite très grand donc LTe devient prohibitif !
3.3. IDENTIFICATION BASÉE SUR L'ERREUR DE SORTIE 25
SBP A(f )
9 zone où la SBPA est un bruit blanc
6
6
6 b
f = 0, 3/T e 6
6
plateau maximum=N Te
6 z }| {
a
6
6
6 6
6 -f
t
1
-
LTe 2 1/Te 2/Te t= 3 à 5 τmax
LTe
Figure 3.6 Paramètres d'une SBPA.
Note 2 : Si fe = 1/Te est bloquée et trop grande (beaucoup de systèmes sont sur-échantillonnés), il
est parfaitement possible choisir Te = n/fe en répétant chaque bit n fois
Note 3 : Si fe est bloquée et trop petite, abandonnez l'idée d'identier la petite constante de temps :
de toutes façons vous ne pourrez pas la contrôler !
3.3 Identication basée sur l'erreur de sortie
Le principe de cette méthode d'identication est illustré en gure 3.7. Le modèle est une fonction de
n paramètres θi , i variant de 1 à n. Il s'agit alors de déterminer les paramètres θi tels que le critère
soit minimum. Le critère est en général choisi de la forme J = Σε2 .
yi
Système
ui
+
εi
Critère

−

Modèle

yî

Figure 3.7 Principe d'identication fondé sur l'erreur de sortie.
La recherche des paramètres optimaux θbi se fait par programmation non linéaire. IL s'agit d'utiliser
un algorithme qui à partir de paramètres non optimaux θi et un critère J donne les paramètres θbi .
Ces algorithmes sont nombreux, à titre de d'exemple voici les plus utilisés.
gradient
quasi Newton
Nelder et Mead
algorithmes génétiques
recuit simulé
Avantages de cette méthode :

pas d'hypothèse sur la forme du modèle : non linéaire
adapté à la recherche de paramètres physiques si modèle de connaissances (modèle continu !)

L'inconvénient majeur est qu'aucun de ces algorithmes n'est capable de garantir que le résultat est
réellement l'optimum.
3.4 Identication basée sur l'erreur de prédiction
Dans cette partie, nous supposerons que le modèle obtenu est un prédicteur, c'est à dire qu'il permet
de calculer la sortie à l'instant i en fonction des entrées et des sorties réelles aux instants précédents
ui−k et yi−k .
yi
Système
ui
+
εi
Critère

−

Modèle

yî

Figure 3.8 Principe d'identication foncée sur l'erreur de prédiction.
3.4.1 Méthode des moindres carrés simples
Les calculs suivants seront fondés sur le modèle présenté en gure 3.9.
βi
vi yi
ui B(z −1 )
A(z −1 )
Figure 3.9 Modèle du système étudié.
La mise en équation du système donne :
Avi = Bui
yi = vi + βi
A(yi − βi ) = Bui
Ayi = Bui + Aβi
posons :
ei = Aβi
ei sont les résidus de l'estimation. On obtient nalement :
Ayi = Bui + ei (3.1)
2
Si nous possédons N mesures consécutives, on peut écrire N −n fois l'équation 3.1.
2. n est l'ordre, supposé connu du polynôme A, on appellera p l'ordre du polynôme B.

3.4. IDENTIFICATION BASÉE SUR L'ERREUR DE PRÉDICTION 27
Sous forme matricielle, on obtient :
 
a1  

yN
 
−yN −1 −yN −2 . . . −yN −n uN . . . uN −p

 a2  eN
yN −1   −yN −2 −yN −3 −y u uN −p−1   
. 
eN −1 
N −n−1 N −1

 . 
 
 .   eN −2
  
yN −2 . .
 
. .
  
  −yN −3 . . uN −2 uN −p−2

.

  an 
 
yN −3 .
 
= + .
  
. .
  
. .   b0 
 
. . −yn+1 . uN −p−3
 
. .
  
. .
  
  b1   .
   
.
 
.
  
yn+2   −yn+1 −yn . uN −p  .   
en+2

 .. 
 
yn+1 −yn −yn−1 . . . −y1 un+1 en+1
bm
soit
y = Xθ + e
le critère J est : X
J= e2 = eT e
donc
J = (y − Xθ)T (y − Xθ) = y T y − θT X T y − y T Xθ + θT X T Xθ
Nous cherchons la valeur θb de θ qui minimise J . Ainsi

∂J
= 0 = −2X T y + 2X T Xθθ=θb

∂θ θ=θb

D'où on en déduit :
θb = (X T X)−1 X T y
Il reste à vérier que la valeur obtenue est bien un minimum
∂ 2 J

= 2X T X
∂θ2
matrice dénie positive : c'est bien un minimum !
3.4.2 Calcul du biais de l'estimateur
b = E[(X T X)−1 X T y]
E[θ] et y = Xθ + e
= E[θ + (X T X)−1 X T e]
= θ + E[(X T X)−1 X T e]
L'estimateur est non biaisé si

E[θ]
b =θ
donc si
E[(X T X)−1 X T e] = 0
il faut donc que : X et e soient non corrélés
e soit centré
Mais ce n'est pas le cas ! Le calcul de ce biais sur l'exemple présenté en gure 3.10 permet de le
montrer.
βi
vi yi
ui bz −1
1+az −1
Figure 3.10 Système du premier ordre bruité.
Les équations du système sont :

vi + avi−1 = bui−1
yi = vi + βi
où βi est un bruit répondant au deux équations suivantes :
E[βi ] = 0
E[βi .βi−k ] = σ 2 δ(k)
En réécrivant les équations du système on obtient :
yi + ayi−1 = bui−1 + βi + aβi−1

= bui−1 + ei
E[ei ] = E[βi + aβi−1 ] = E[βi ] + aE[βi−1 ] = 0
E[ei .ei−k ] = E[(βi + aβi−1 )(βi−k + aβi−k−1 )]

= E[βi .βi−k ] + aE[βi .βi−k−1 ] + aE[βi−1 .βi−k ] + a2 E[βi−1 .βi−k−1 ]
= σ 2 [(1 + a2 )δ(k) + aδ(k − 1) + aδ(k + 1)]
6σ 2 (1 + a2 )
aσ 2 6 6
aσ 2
-
-1 0 1 k
Figure 3.11 Fonction d'auto-corrélation des résidus d'estimation.
La fonction d'auto-corrélation (voir g. 3.11) montre bien que les résidus ei ne sont pas un bruit
blanc, notre estimateur est bien biaisé ! Il reste donc à trouvé une méthode qui permette de rendre ce
biais nul. Deux méthodes sont proposées ci-après : La méthode des moindres carrés généralisés et la
méthode de la matrice instrumentale.
3.4. IDENTIFICATION BASÉE SUR L'ERREUR DE PRÉDICTION 29
3.4.3 Méthode des moindres carrés généralisés
Le système étudié est décrit sur la gure 3.12
ξi
P (z −1 )
Q(z −1 )
βi
B(z −1 ) vi yi
m
A(z −1 )
ui ei
n
yî
modèle
Figure 3.12 Système étudié.
Avi = Bui
yi = v i + β i
A(yi − βi ) = Bui
Ayi = Bui + Aβi (3.2)
les résidus sont donc :
ei = Aβi
écrits en fonction du bruit blanc ξi

P
ei = A ξi
Q
équation que nous écrirons dans un soucis de simplication d'écriture :
F ei = ξi
En reprenant l'équation 3.2 et en multipliant tous les termes par F on obtient :
AF yi = BF ui + AF βi
Q
= BF ui + A βi
AP
= BF ui + ξi
En réécrivant cette dernière équation,
Ayi? = Bu?i + ξi
où yi? et u?i sont, respectivement les sorties yi et les commandes ui ltrées par le ltre F.
En appliquant la méthode des moindres carrés simples sur ces entrées-sorties ltrées le biais de l'esti-
mateur sera nul car le bruit ξi est bien un bruit blanc. L'estimation des paramètres sera donc parfaite.
Toutefois il reste à déterminer le ltre F.
Détermination de F Reprenons 3.4.3 sous forme matricielle (avec F de la forme : F = 1 + f1 z −1 +

f2 z −2 + . . . + fn z −n )


  f 1
eN 
−eN −1 −eN −2 . . . −eN −n

 f2 

ξN

eN −1
 
  −eN −2 −eN −3 −eN −n−1   .   ξN −1
  
  ..  
 
eN −2

. .
  ..   ξN −2
  
. .
   
  −eN −3 . .
  
eN −3


=  .  +  .
  .

.
  ..   .

. .
 

 .
.
 
  . −en+1   .   ..


−en+1 −en   .
  
en+2

 .. 
  
en+1 −en −en−1 . . . −e1  .  ξn
fn
soit
e = Xe f + ξ
qui est une forme connue et dont on sais que la solution est :
fb = (XeT Xe )−1 XeT e
Les paramètres f de F sont donc déterminés à condition de disposer des résidus e. Mais c'est justement
pour les avoir que l'on déroule toutes ces équations...
La solution proposée consiste alors à procéder par itérations en utilisant une première estimation des
résidus e obtenue par un premier modèle utilisant
1 θ une première approximation des paramètres du
système obtenue par la méthode des moindres carrés simples.
L'algorithme est alors, à partir des entrées-sorties

1u = ui et
1y = yi
i i
1. Estimation des paramètres

kθ = (k X T k X)−1k X T k y
2. Estimation des résidus
ke = ky − kXkθ
3. Construction du ltre F1
k fb = (k X T k X )−1k X T k e
e e e
4. Filtrage les données
k y? = k yi et
k u? = k F k ui
i i
5. Rebouclage en utilisant ce nouvel ensemble d'entrées sorties
Quelques remarques :
k F sont du premier ordre (n

Si les ltres = 1), la matrice Xe se réduit à un réel ce qui simplie la
k k Tk
formule fb = ( Xe Xe )
−1k X T k e = k e .
e kX
e
Critère d'arrêt de l'algorithme : Le premier réexe consiste à calculer la "blancheur" des résidus, mais
ce critère est assez dicile à mettre en ÷uvre. Un méthode plus simple consiste à calculer un critère
d'arrêt de la forme :
∆J k J − k−1 J
Ca = = kJ
avec J = eT e
J
En pratique on utilisera Ca ≤ 5%
A chaque itération le bruit "devient plus blanc" et pas conséquent le critère décroît.
3.5. IDENTIFICATION EN BOUCLE FERMÉE 31
3.4.4 Méthode de la Matrice Instrumentale
Puisque le biais de l'estimateur des moindres carrés est biaisé à cause de la corrélation entre X et
e, on se propose de déterminer une autre matrice qui permette le calcul de θ tout en évitant cette
corrélation et donc le biais de l'estimateur.
Posons :
θb = (Z T X)−1 Z T y
où Z est la matrice instrumentale.

Quelles sont les conditions sur Z pour que l'équation précédente ait un sens ?
En poursuivant le calcul :
θb = (Z T X)−1 Z T (Xθ + e)
θb = θ + (Z T X)−1 Z T e
Pour que
E[θ]
b =θ
il faut :
E[Z T X] 6= 0 et E[Z T e] = 0
Plusieurs choix de Z sont possibles, en voici deux :

le premier choix consiste à faire deux campagnes de mesures avec la même entrée ui , le plus souvent
consécutives (décalées dans le temps). Vous avez alors la possibilité de déterminer deux matrices X
soit X1 et X2 . En posant
θb = (X2T X1 )−1 X2T y1

il n'y a plus de corrélation entre X2 et y1 (le bruit est stochastique et ergodique) le biais est bien nul.
La deuxième proposition consiste à créer les données nécessaires à la création de la matrice X2 à l'aide
d'un modèle de type moindre carrés simples issu de la première campagne de mesure.
3.5 Identication en boucle fermée
Ce domaine de l'identication est actuellement très actif. En eet ces modèles sont nécessaires dans le
cadre de la commande prédictive. Plus pragmatiquement, rappelons que le plus souvent l'identication
en milieu industriel s'eectue sur une machine en production. Il est par conséquent dicile voire
impossible d'ouvrir la boucle pour procéder à une identication. C'est tout particulièrement le cas des
systèmes instables.
Il reste néanmoins possible d'obtenir un modèle. Le système de la gure 3.13 représente le système à
C(z −1 )
identier en boucle fermée. le correcteur du système est représenté par la fraction rationnelle .
D(z −1 )
Trois approches sont alors possibles.
La méthode directe Dans ce cas on utilise les entrées vi et les sorties yi . L'identication se fait
alors comme en boucle ouverte. L'avantage de cette approche est de ne pas nécessiter la connaissance
du régulateur. Les méthodes ARX, ARMAX et les modèles d'état donnent de bons résultats. Les
meilleurs résultats seront obtenus si :
Le modèle de bruit est bon
la boucle de retour n'aecte pas ou peu le spectre du signal d'entrée (vi ))
le rapport signal à bruit est important (peu de bruit)
Dans le cas contraire, l'estimateur est biaisé.
ξi
P (z −1 )
Q(z −1 )
ui vi + yi
+ m C(z −1 ) B(z −1 ) m
H(z −1 ) = D(z −1 )
G(z −1 ) = A(z −1 )
− +
Figure 3.13 Identication des systèmes en boucle fermée.
La méthode indirecte Dans ce cas on utilise les entrées ui et les sorties yi . Bien entendu on identie
alors le système en boucle fermée soit :
B 0 (z −1 ) H(z −1 )G(z −1 ) C(z −1 )B(z −1 )

M (z −1 ) = = =
A0 (z −1 ) 1 + H(z −1 )G(z −1 ) C(z −1 )B(z −1 ) + D(z −1 )A(z −1 )
La connaissance du régulateur du régulateur nous permet de déterminer le modèle du système :
M (z −1 )
G(z −1 ) =
H(z −1 ) − M (z −1 )H(z −1 )
L'avantage de cette méthode est que n'importe quelle méthode d'identication est applicable et don-
nera un modèle du système en boucle fermée. Par contre, la moindre erreur sur le régulateur (pas
forcement bien connu, paramètres, saturations ...) se retrouve dans le modèle.
La méthode "Joint Input-Output" Avec cette méthode on se propose de déterminer simultané-

ment les deux modèles :
le modèle du système
le modèle du correcteur.
Ceci s'eectue facilement en utilisant les modèles paramétriques d'état sous Matlab.

y
= Au + Bξ
v
La matrice A contient alors les dynamiques du correcteur Guv et celle du système en boucle fermée
Guy . Comme pour la méthode précédente on remonte au modèle du système par calcul.
r
3.6 Et Matlab dans tout ça ?
Matlab est avant tout un outil de calcul ! Scilab, gratuit, présente les mêmes fonctionnalités, il est
vrai, sans les boites à outils du moins pas aussi complètes.
Matlab propose donc toute une série d'outils d'aide à l'identication. En particulier, quatre modèles
linéaires, les plus couramment utilisés sont des fonctions de Matlab.
Dans le tableau 3.2 ces modèles sont présentés avec leurs utilisations "classique" et la méthode de
détermination des paramètres. Notez qu'une lecture de la documentation concernant la boite à outils
"identication"est une excellente activité !
Face à la facilité de calcul de modèles, le problème n'est plus vraiment de choisir une structure adaptée
mais faire le tri parmi les quantités de modèles calculés.
3.6. ET MATLAB
r DANS TOUT ÇA ? 33
Table 3.2 Modèles paramétriques implantés sous Matlab.

ξi Modèle ARX
Le modèle le plus simple, donne souvent de bons résultats.
1 le seul hic est le traitement du bruit qui est soumis à
A(z −1 )
la même dynamique que l'entrée. A utiliser en première
approximation ou lorsque le bruit est surtout à l'entrée.
ui +
yi
B(z −1 ) m
A(z −1 ) +
Méthode de détermination des paramètres :
Moindres carrés ou Matrice instrumentale
ξi
Modèle ARMAX
C(z −1 ) Proche du modèle ARX, il s'utilise dans les mêmes cas. Il
A(z −1 )
permet en outre de créer un modèle de bruit un peu plus
réaliste. C'est le modèle le plus utilisé.
ui +
yi
B(z −1 ) m
A(z −1 ) +
Maximum de vraisemblance
ξi
Modèle OE
Bien que semblant plus simple que les précédents, le calcul
des paramètres s'avère plus dicile. Parfait lorsque le bruit
est surtout un bruit de capteur, donc proche de la sortie.
ui +
yi
B(z −1 ) m
A(z −1 ) +
ξi
P (z −1 ) Modèle de Box-Jenkins
Q(z −1 )
Le modèle complet par excellence, la dynamique diérente
pour l'entrée et le bruit en font un bon modèle.
ui +
yi
B(z −1 ) m
A(z −1 ) +
Chapitre 4
Estimations récursives
4.1 Moindres carrés récursifs
L'estimation de paramètres par la méthode des moindres carrés simples présente un inconvénient
majeur, la nécessité de calculer l'inverse d'une matrice, ce qui est long et parfois impossible sur un
microcontrôleur. On se propose à travers cet exercice de déterminer une forme récursive de cette
estimation.
Les avantages d'une formulation récursive tiennent essentiellement en deux points.
La possibilité de traiter un plus grand nombre de données que dans le cas de la formulation directe
(pas de pseudo-inverse à calculer), notamment dans le cas de l'implantation sur un microcontrôleur.
Dans le cas des systèmes variants dans le temps, la forme récursive permet de "suivre" les paramètres
du système. Dans ce cas, l'identication en ligne est généralement suivie d'une commande qui elle
aussi "s'adapte" aux paramètres courants, on entre alors dans le domaine des commandes auto-
adaptatives qui dépassent le cadre de ce cours.
Comme dans tout problème récursif on s'intéresse d'abord à la boucle, ensuite comment on en sort et
enn comment on y entre. Supposons que nous possédions une estimation des paramètres à l'instant
N : θN
T
θN = (XN XN )−1 XN
T
YN
à l'instant suivant la nouvelle estimation est :
T −1 T
θN +1 = (XN +1 XN +1 ) XN +1 YN +1
avec :

xN +1 yN +1 eN +1
XN +1 = , YN +1 = et EN +1 =
XN YN EN
où xN +1 = [−yN − yN −1 . . . − yN −n+1 uN +1 . . . uN −p+1 ]

θN +1 s'écrit alors :
T −1 T
θbN +1 = (XN +1 XN +1 ) XN +1 YN +1
T
= (XN XN + xTN +1 xN +1 )−1 (XN
T
YN + xTN +1 yN +1 )
En utilisant le lemme d'inversion matricielle :
(A + BCD)−1 = A−1 − A−1 B(C −1 + DA−1 B)−1 DA−1
on obtient :
35
36 CHAPITRE 4. ESTIMATIONS RÉCURSIVES
−1
T
XN )−1 − (XN
T
XN )−1 xTN +1 1 + xN +1 (XN
T
XN )−1 xTN +1 T
XN )−1

θbN +1 = (XN xN +1 (XN
T
YN + xTN +1 yN +1

XN
En posant
T X )−1 xT
a = 1 + xN +1 (XN N N +1
−1
T
XN )−1 − (XN
T
XN )−1 xTN +1 1 + xN +1 (XN
T
XN )−1 xTN +1 T
XN )−1

T
YN + xTN +1 yN +1

XN
T
XN )−1 − (XN
T
XN )−1 xTN +1 a−1 xN +1 (XN
T
XN )−1 XN
T
YN + xTN +1 yN +1

θbN +1 = (XN
T
= (XN XN )−1 XN
T T
YN +(XN XN )−1 xTN +1 yN +1 − (XN
T
XN )−1 xTN +1 a−1 xN +1 (XN
T
XN )−1 XN
T
YN
| {z } | {z }
θbN θbN
−(XN XN )−1 xTN +1 a−1 xN +1 (XN

T T
XN )−1 xTN +1 yN +1
= θbN + (XN T
XN )−1 xTN +1 yN +1 − (XNT
XN )−1 xTN +1 a−1 xN +1 θbN
−(XN T
XN )−1 xTN +1 a−1 xN +1 (XNT
XN )−1 xTN +1 yN +1

T
= θbN − (XN XN )−1 xTN +1 a−1 xN +1 θbN
T
XN )−1 xTN +1 a−1 −xN +1 (XN T
XN )−1 xTN +1 yN +1 + (XN T
XN )−1 xTN +1 yN +1

+(XN

= θbN − (XNT
T
XN )−1 xTN +1 −a−1 xN +1 (XN T
XN )−1 xTN +1 yN +1 + yN +1

+(XN

= θbN − (XNT
T
XN )−1 xTN +1 a−1 −xN +1 (XN T
XN )−1 xTN +1 + a yN +1

+(XN

= θbN − (XNT
+(XN T
XN )−1 xTN +1 a−1 yN +1

= θbN + (XNT
XN )−1 xTN +1 a−1 yN +1 − xN +1 θbN
On obtient nalement une formule de récurrence :

T
θbN +1 = θbN + (XN XN )−1 xTN +1 a−1 yN +1 − xN +1 θbN

Le terme yN +1 − xN +1 θbN représente l'erreur d'estimation à l'aide des paramètres précédents. L'équa-
tion est alors : la nouvelle estimation est l'ancienne estimation corrigée par un terme proportionnel à
l'erreur d'estimation précédente, que l'on peut récrire sous la forme.

θbN +1 = θbN + KN +1 yN +1 − xN +1 θbN
L'algorithme général est donc :

θbN +1 = θbN + KN +1 yN +1 − xN +1 θbN (4.1)
−1
KN +1 = PN xTN +1 I + xN +1 PN xTN +1 (4.2)
PN +1 = PN − KN +1 xN +1 PN (4.3)
4.2. VARIABLE INSTRUMENTALE RÉCURSIVE 37
On remarquera une analogie avec le ltre de Kalman, en fait c'est bien l'équivalent d'un ltre de
Kalman appliqué sur le système suivant :
θN +1 = θN
yN +1 = xN +1 θN + eN +1
La condition de convergence est la même : eN doit être un bruit blanc.

L'initialisation de l'algorithme se fait de deux façons :
Si on connais une première estimation de θ (en provenance d'une méthode non récursive par
exemple), on l'utilise. Dans ce cas on prend P0 = λI avec λ petit (petite variance du bruit)
Si on ne connais pas de première approximation, on prend θ0 quelconque et P0 = λI avec λ grand
(grande variance du bruit).
Cet algorithme ne peut être utilisé que si les paramètres du système sont constants, en eet lorsque N
tend vers l'inni, PN tend vers 0, alors une variation même importante des paramètres n'inue plus
sur l'estimation de θ . Notons encore que sur le plan numérique la forme PN +1 = PN − KN +1 xN +1 PN
est très mal conditionnée, on lui préfère la forme :
PN +1 = (I − KN +1 xN +1 )PN (I − KN +1 xN +1 )T + KN +1 KN
T
+1
Si les paramètres évoluent brusquement, une solution consiste à réinitialiser PN = λI avec λ grand Si
les paramètres évoluent lentement on peut utiliser :

−1
KN +1 = PN xTN +1 λI + xN +1 PN xTN +1 (4.5)
1
PN +1 = (PN − KN +1 xN +1 PN ) (4.6)
λ
Cependant cet algorithme présente l'inconvénient de faire croître PN de façon exponentielle s'il n'y a
plus d'excitation.
Voici d'autres choix pertinents
à gain constant : on force alors PN +1 = PN ; cela évite la diminution du gain en cours de recherche
des paramètres, et on donne ainsi plus de poids aux acquisitions les plus récentes ; cette option
convient bien à un système dont les paramètres varient.
à gain décroissant : PN = cst, cas classique, qui convient à un système à paramètres constants.
à trace constante : on garde tr(PN ) = cst (multiplication par un facteur correctif à chaque itération) ;
on a les mêmes avantages que dans la recherche à gain constant mais on module en cours de recherche
le poids relatif de chaque paramètre.
Cette méthode présente les mêmes avantages que la méthode des moindre carrées non récursive,
mais aussi le même inconvénient : l'estimateur est biaisé ! Aussi on lui préfère généralement d'autres
méthodes telle que la méthode de la variable instrumentale récursive.
4.2 Variable instrumentale récursive
L'algorithme, très proche de celui énoncé en

T 2 T
−1
KN +1 = PN ZN +1 σ + xN +1 PN ZN +1 (4.8)
PN +1 = PN − KN +1 xN +1 PN (4.9)
Dans ce cas ZN +1 représente un vecteur ligne instrumental composé, un peu comme la matrice ins-
trumentale, soit d'observations retardées de la sortie soit de la sortie d'un modèle auxiliaire.
38 CHAPITRE 4. ESTIMATIONS RÉCURSIVES
Beaucoup d'autres algorithmes existent, le lecteur intéressé en trouvera pléthore dans la littérature.
Néanmoins ils sont tous a peu près fondés sur le même principe. En fait, ils se résument à l'expression
d'un asservissement à gain réglable : on asservit les paramètres pour annuler une erreur, d'équation
ou de sortie selon la méthode. La valeur du gain inue, comme sur n'importe quel système, sur la
stabilité, soit la convergence de la méthode, et sur la rapidité.
Moduler le gain revient donc à moduler la rapidité de convergence : si l'on a aaire à un procédé dont il
faut poursuivre les paramètres, qui varient eectivement dans le temps, on aura intérêt à avoir un gain
fort ; si par contre le système est fortement bruité, on a intérêt à avoir un gain faible, sinon la sortie
du modèle va poursuivre le bruit en faisant varier les paramètres à chaque période d'échantillonnage,
ce qui n'a pas de sens physique. Un bon réglage du gain demande plusieurs essais et pas mal de bon
sens !
Chapitre 5
Algorithmes d'optimisation paramétrique.
5.1 Présentation générale.
Ce chapitre est consacré au choix d'une méthode d'optimisation des paramètres. Une étude biblio-
graphique succincte montre que le domaine de l'optimisation paramétrique est en pleine évolution
notamment dans le cadre des méthodes stochastiques.
Nous présenterons les méthodes les plus connues an d'en extraire les principes d'optimisation les plus
aptes à résoudre notre problème.
Les méthodes présentées sont en fait celles qui reviennent le plus souvent dans la littérature. Parfois
anciennes, ces méthodes, en raison de leur capacité à optimiser des fonctions très diérentes, sont
encore souvent utilisées aujourd'hui. Par ailleurs, ces méthodes sont souvent à la base de nouvelles
méthodes d'optimisation.
Le domaine d'utilisation de chaque méthode est un point délicat à aborder. Une façon simple de
déterminer ce domaine serait de le limiter au domaine de convergence garantie, en général les fonctions
convexes. Mais cette approche est quelque peu réductrice et ne tient pas compte de la robustesse de la
méthode, c'est pourquoi les spécialistes de l'optimisation paramétrique préfèrent se référer à des tests
sur des fonctions connues et présentant des caractéristiques des plus diverses.
5.2 Exemple introductif
Imaginons que la fonction à minimiser soit l'altitude et que l'espace sur lequel nous cherchons à
optimiser l'altitude soit la surface de notre planète. Il existe plusieurs façons de procéder.
La première façon consiste à parcourir l'ensemble des points de la surface et de mesurer l'altitude en
chaque point. En supposant un maillage de la planète de 10m sur 10m et à raison d'une mesure par
seconde : cela représente environs 160 000 ans de travail (énumération).
La deuxième méthode consiste à partir d'un point et de descendre dans le sens de la plus grande pente.
Nous évitons ainsi tous les sommets locaux mais un observateur qui partirai d'un point des Alpes a
peu de chances de trouver le fond de la mer (gradient). Nous pouvons améliorer cette méthode en
prenant trois observateurs qui se communiquent leur altitude. Celui qui se situe au point le plus haut
se déplace en direction d'un point situé entre les deux autres observateurs et parcoure deux fois la
distance qui le sépare de ce point. Ces observateurs s'épargnent le calcul de la pente mais ont toujours
peu de chances de trouver la mer (Nelder et Mead).
Une troisième méthode consiste à donner à notre observateur initial une soucoupe volante (et amphibie)
qui lui permet de se déplacer très vite. Pendant un temps il se pose au hasard sur la planète et note
l'altitude et la position de chaque point. puis au cours du temps il concentre sa recherche autour des
points les plus bas en eectuant des déplacements aléatoires mais de plus en plus cours en moyenne.
Cet observateur trouverai très rapidement un océan et déterminerai son point le plus bas. Par contre
si il se concentre trop vite sur l'un des océans il risque de rater le point le plus bas de la planète (recuit
simulé).
39
40 CHAPITRE 5. ALGORITHMES D'OPTIMISATION PARAMÉTRIQUE.
La quatrième méthode consiste naturellement à prendre plusieurs observateurs qui feraient le même
travail que le précédent. Rapidement il se constitueraient en groupes qui déterminent le relief de chaque
mers et océans. Lorsqu'un groupe atteint le minimum de son étendue d'eau et remarque qu'un autre
groupe obtient un minimum plus faible, le premier groupe se déplace pour aider un autre groupe
n'ayant toujours pas atteint le fond (algorithmes génétiques).
5.3 Problème posé
L'idée générale de l'optimisation paramétrique consiste à déterminer un ensemble de paramètres qui

optimise un critère soumis à des contraintes. Dans la suite de ce mémoire, nous parlerons exclusivement
de minimisation, la maximisation étant la minimisation de l'opposé de la fonction.
minimiser f (x),
avec x = {x1 , x2 , x3 , · · · , xn },
(5.1)
soumis à Ck (x) ≥ 0 k = 1 à m,
et Ck (x) = 0 k = m à p.
Dans la première partie de ce chapitre, nous nous consacrerons à la détermination de trajectoires
optimales sans contrainte. Ce cas, bien qu'inapplicable dans le cas des robots marcheurs, peut très
bien l'être pour des robots manipulateurs.
5.4 Revue bibliographique succincte.
Les principales méthodes d'optimisation peuvent être classées en trois catégories principales :
1. Les méthodes analytiques, qui déterminent le point de calcul suivant en fonction des caractéris-
tiques de la fonction au point considéré.
2. Les méthodes stochastiques, qui intègrent une part de hasard dans le choix du point de calcul
suivant.
3. L'énumération, qui consiste à discrétiser l'espace engendré par le vecteur à optimiser et à calculer
l'ensemble des possibilités.
La gure 5.1 présente les diérents groupes d'algorithmes d'optimisation paramétrique et leur classi-
cation.
L'énumération : Elle consiste à prendre un nombre ni de valeurs pour chaque paramètre et à
calculer le critère associé pour l'ensemble des possibilités. Dans notre cas elle est inutilisable compte
tenu du nombre de paramètres. En utilisant 50 valeurs par paramètre et 13 paramètres par axe, il faut
examiner 5013×5 possibilités !
Les algorithmes d'ordre n : Fondés sur une décomposition à l'ordre n en séries de Taylor de la
fonction à minimiser, ces algorithmes utilisent la dérivée d'ordre n pour la détermination du point de
calcul suivant. Ils sont très rapides lorsque la fonction est continue, dérivable et convexe. Par contre,
ils sont peu robustes si la fonction présente des minima locaux.
Les algorithmes d'ordre 0 : N'utilisant que la valeur de la fonction en certains points, ces algo-
rithmes sont lents mais robustes en cas de discontinuités de la fonction.
Les algorithmes génétiques : Transposition à l'informatique de la sélection naturelle selon Charles

Darwin, ils utilisent plusieurs vecteurs de paramètres (une population) qui évoluent vers le minimum
par la mise en ÷uvre de trois opérations stochastiques, la sélection, le croisement et la mutation.
5.5. ELÉMENTS DE COMPARAISON DES DIFFÉRENTES MÉTHODES. 41
Figure 5.1 Classes d'algorithmes d'optimisation paramétrique
Les algorithmes de recuit simulé : Algorithmes stochastiques par excellence, ils sont statisti-
quement capables de déterminer le minimum global d'une fonction quelles que soient ses caractéris-
tiques. Par contre, le temps de calcul nécessaire à l'obtention du minimum est long.
Les stratégies d'évolution : Fondées sur une évaluation itérative de la fonction, le déplacement
est aléatoire mais il s'adapte en fonction de la topologie locale de la fonction à minimiser. La diérence
fondamentale entre stratégies d'évolution et le recuit simulé est que la première méthode explore la
fonction de façon locale alors que la deuxième explore l'ensemble de l'espace.
5.5 Eléments de comparaison des diérentes méthodes.
5.5.1 Critères.
La robustesse est l'aptitude de la méthode à converger vers un minimum global sans être perturbée
par des minima locaux. Cette notion couvre les notions de robustesse à la rugosité de la fonction et la
notion d'aptitude à trouver le minimum global d'une fonction.
La vitesse est mesurée en nombre d'itérations sur la fonction de coût, le temps de calcul nécessaire
pour l'algorithme de la méthode elle-même est considéré comme négligeable.
La globalité du minimum est donnée pour une minimisation d'une fonction multimodale sans aucune
connaissance a priori de la position du minimum global, donc d'une initialisation aléatoire.
5.5.2 Comparaisons.
Nous reprendrons ici quelques comparaisons trouvées dans la littérature an d'établir un classement
des diérentes méthodes par rapport aux trois critères dénis précédemment, à savoir : la robustesse,
la vitesse et la globalité du minimum.
Les algorithmes d'ordre n sont très rapides, d'autant plus rapides que les fonctions optimisées ré-
pondent aux hypothèses formulées, par contre ils sont peu robustes en cas de discontinuités. Dans le
cas de discontinuités, parmi les méthodes analytiques, ce sont les algorithmes d'ordre 0 qui sont les
meilleurs.
0% épistasie 100%
| {z } | {z }
gradient aléatoire
| {z }
Quasi Newton
| {z }
Simplex
| {z }
algorithmes génétiques
| {z }
recuit simulé
Figure 5.2 Récapitulatif des diérentes méthodes en fonction de l'épistasie.
5.5.3 Récapitulatif.
Le tableau 5.1 est une synthèse des données recueillies dans la littérature. Les notations sous forme
de + et − permettent seulement une comparaison qualitative entre les méthodes sans augurer de la
qualité intrinsèque d'une méthode.
Principe robustesse vitesse minimum

ordre 1 : Gradient - - - ++ local
ordre 2 : Quasi Newton - - +++ local
ordre 0 : Simplex - + local
Algorithmes d'évolution + - - global -
Algorithmes génétiques + - - - global ++
Recuit simulé ++ - - - global ++
Table 5.1 Comparaison des méthodes d'optimisation.

La gure 5.2 donne le domaine d'utilisation de chaque méthode en fonction de l'épistasie. Cette no-
tion directement tirée du vocabulaire de la génétique indique la propension d'un gène à inhiber ou
augmenter l'eet d'un autre gène sur la valeur de la fonction. Ainsi une fonction non linéaire, multi-
modale présente une forte épistasie, par contre près d'un minimum local de cette fonction, l'épistasie
est minimale voire nulle.
An de xer les idées, voici trois exemples de fonctions :
f1 (x) = x21 + x22 + x23 , épistasie nulle (5.2)
f2 (x) = 100(x21 − x2 )2 + (1 − x3 )2 , épistasie moyenne (5.3)
f3 (x) = (x22 − x3 )2 si sin(x1) > 0

(5.4)
= (x21 − x2 )2 si sin(x1) < 0, épistasie forte
Ainsi f1 (x) présente une épistasie nulle puisque la valeur d'un gène ne change pas l'inuence d'un
autre gène sur le critère. Dans f2 (x) l'inuence existe entre x1 et x2 . Dans f3 (x) l'inuence est forte
puisque la valeur de x1 décide de quels gènes interviennent dans le critère.
5.5.4 Synthèse
La meilleure méthode serait une méthode qui évolue en fonction de la connaissance de la fonction
à optimiser. Ainsi la recherche dans ce domaine s'oriente vers le couplage de plusieurs méthodes
5.5. ELÉMENTS DE COMPARAISON DES DIFFÉRENTES MÉTHODES. 43
d'optimisation. Le problème d'optimisation revient alors à eectuer successivement :
1. Dénition du problème et en particulier des limites de l'espace à étudier notamment par la prise
en compte de limites sur la valeur des diérents paramètres à optimiser.
2. Génération aléatoire d'une population respectant les limites précédentes et évaluation de chaque
individu vis-à-vis du critère.
3. Optimisation lente permettant :

une exploration de l'espace total,
la détermination de sous-espaces convexes présentant un minimum local.
4. Optimisation rapide dans les sous-espaces convexes.
Le passage d'une méthode à l'autre peut être fondé sur la détermination de l'épistasie.
Chapitre 6
Description des algorithmes utilisés.
6.1 Gradient et quasi-Newton
Ces méthodes consistent à écrire un développement limité de la fonction autour du point courant, puis
d'en déduire par la minimisation de ce développement limité, le point suivant. Le lecteur intéressé par
ce chapitre lira avec bonheur Numerical recipes in C [5] où les "détails" des algorithmes (initialisation,
choix de décroissance, ...) seront parfaitement explicités.
6.1.1 Méthode du gradient
La dérivée de f : Rn → R s'écrit :
∂f
Df (x)i =
∂xi
Df (x) est un vecteur ligne.
Le gradient de la fonction est :
T
∂f
∇f (x) = Df (x)T =
∂x
∇f (x) est un vecteur colonne.

La fonction f à minimiser peut être approchée autour d'un point x0 par :
f (x0 + ∆x) ' f (x0 ) + ∇f (x0 )T ∆x
le point de calcul suivant sera calculé par :
xi+1 = xi + ∆x
le pas étant déterminé par :
∆x = −k∇f (x)
k étant un réel positif. Le choix de k donne lieu à de nombreuses variantes de l'algorithme.

si k est trop petit, l'algorithme converge lentement,
si k est trop grand, l'algorithme oscille autour du minimum,
En pratique on choisi k variable au cours de l'optimisation, la plus connue est le line search. Une
méthode simple consiste à faire varier k suivant :
Si f (xi+1 ) < f (xi ) alors augmenter k (on est dans la bonne direction, on va voir plus loin)
Si f (xi+1 ) = f (xi ) alors augmenter k (on est sur un plateau, on va voir plus loin)
Si f (xi+1 ) > f (xi ) alors diminuer k et abandonner f (xi+1 ) (on remonte, le minimum est plus près
de xi )
45
46 CHAPITRE 6. DESCRIPTION DES ALGORITHMES UTILISÉS.
6.1.2 Méthode de Newton
C'est en fait, l'extension à l'ordre 2 de la méthode du gradient. Ainsi la fonction f à minimiser est
approchée par :
1
f (x0 + ∆x) ' f (x0 ) + ∇f (x0 )T ∆x + ∆xT ∇2 f (x0 )∆x
2
ou ∇2 f (x0 ) est le Hessien de la fonction (la dérivée seconde)
∂2f ∂2f ∂2f

 
∂x21 ∂x1 ∂x2 ··· ∂x1 ∂xn
∂2f ∂2f ∂2f
 
∂ 2f 
∂x2 ∂x1 ∂x2 ∂x2 ··· ∂x2 ∂xn

D∇f (x) = ∇2 f (x)ij = =
 
. . .. ∂2f
∂xi ∂xj

. . .
. .
 
 ∂xi ∂xj 
∂2f ∂2f ∂2f
∂xn ∂x1 ∂xn ∂x2 ··· ∂xn ∂xn
Le point où la fonction approchée est minimum est tel que la dérivée de la fonction est nulle soit
l'équation :

∂ T 1 T 2
f (x0 ) + ∇f (x0 ) ∆x + ∆x ∇ f (x0 )∆x = 0
∂∆x 2
∇f (x0 )T + ∇2 f (x0 )∆x = 0
le pas est alors :

∆x = −∇2 f (x)−1 ∇f (x)
On vois que par rapport à la méthode du gradient, la direction de recherche est pondérée par l'in-
verse du hessien. En pratique cette méthode est très peu utilisée car la détermination du hessien est
longue et coûteuse en temps de calcul, son inverse aussi. On lui préfèrera des méthodes fondées sur
l'approximation du hessien.
6.1.3 Méthode de Newton-Raphson
Comme pour la méthode de Newton le pas est calculé par :
∆x = −∇2 f (x)−1 ∇f (x)

1
Mais au lieu de calculer le hessien puis son inverse on se contente de calculer ∆x tel que :
∇2 f (x)∆x = ∇f (x)
Notez bien que l'on a transformé le calcul de l'inverse de ∇2 f (x)−1 en une résolution d'un système
linéaire à chaque itération.
6.1.4 Méthode de Levenberg-Markardt
Les deux méthodes précédentes présentent l'inconvénient (que quelques variantes contournent) de ne
pas converger lorsque le point initial est loin du minimum recherché.
L'algorithme Levenberg-Markardt lève cet inconvénient on proposant un compromis entre la méthode
du gradient (robuste mais lente à l'approche du minimum) et celle de Newton (peu robuste loin du
minimum mais très ecace près). Le pas ∆xi s'obtient en résolvant le système d'équations linéaires :
(∇f (xi )∇f (xi )T + λi I)∆xi = −∇f (xi )T f (xi )
ainsi, lorsque λi est grand devant le hessien, l'algorithme tend vers l'algorithme du gradient, lorsque
λi est petit, il converge aussi vite que celui de Newton.
1. Pour une fonction de 100 paramètres le hessien est une matrice 100 × 100 !
6.1. GRADIENT ET QUASI-NEWTON 47
6.1.5 l'algorithme BFGS
BFGS sont les initiales de Broyden-Fletcher-Goldfarb-Shanno, c'est sans doute l'algorithme le plus
utilisé.
A partir d'un point initial x0 et une première approximation du hessien H0
1. Calculer le pas par :
Hi ∆xi = −∇f (xi )
2. Déterminer αi optimal par une méthode de type "line search" et calculer
xi+1 = xi + αi ∆xi
3. calculer la variation du gradient :
∇f (xi+1 ) − ∇f (xi )
yi =
αi
4. puis la nouvelle estimation du hessien
yi yiT Hi ∆xi ∆xi T HiT

Hi+1 = Hi + +
yiT ∆xi ∆xi T Hi ∆xi
L'ensemble de ses méthodes convergent bien vers un minimum à condition que la fonction soit convexe !
Le principal inconvénient de ses méthodes est de devoir calculer le gradient de la fonction à chaque
itération aussi les méthodes suivantes se passent de cette étape de calcul et sont nettement plus
robustes si la fonction présente des minima locaux.
6.2 Simplex de Nelder & Mead.
Parmi les algorithmes d'ordre 0, la méthode d'optimisation Simplex de Nelder & Mead [5], est l'une
des plus anciennes connues et est encore aujourd'hui souvent utilisée ( fminsearch de matlab). La
démonstration de la convergence de cet algorithme dans le cas de fonctions convexes, n'a été réalisée
qu'en 1995.
Un Simplex est un tableau de n+1 vecteurs (
i x) de n paramètres (αi ). La méthode procède à une
série de :
réexions (le volume du Simplex reste constant)
pour déplacer le barycentre du Simplex en direction du point de critère le plus faible,
expansions (le volume augmente)
qui étendent le Simplex en direction du point de critère le plus faible,
contractions (le volume diminue)
qui lui permettent de passer dans des goulets.
Algorithme
L'algorithme est présenté par l'ordinogramme de la gure 6.5 Le vecteur donnant le critère le plus
grand à la k ème itération (xkmax ) est projeté à travers le barycentre des autres vecteurs du Simplex
(xk ).
xr = (1 + α)xk − αxkmax , avec α = 1, réexion, (6.1)
La gure 6.1 représente la position des vecteurs du Simplex dans l'espace des paramètres (ici représenté
dans un espace de dimension 3).
X Simplex après réexion

Simplex original
B xr
B XXXXX
1 B

B

,,
B B

HH B HH B ,,

H H

HB HB,
Figure 6.1 Image d'un Simplex de 4 vecteurs avant et après réexion.

Si la valeur du critère f (xr ), associée au vecteur xr est la plus faible obtenue jusqu'à présent, on tente
une expansion du Simplex dans la direction de la projection. Sinon on conserve xr et on recommence
une nouvelle itération. La procédure d'expansion (voir gure 6.2 est donnée par :
xe = γxr + (1 − γ)xk , avec γ = 2, expansion, (6.2)
xe
1
((
Simplex original
B (( (((
Simplex
après
expansion
B (
B
B B

HH B

HHB B

HH B
HHB
Figure 6.2 Image d'un Simplex de 4 vecteurs avant et après expansion.

Si la valeur du critère f (xe ) est plus faible que celle associée aux autres vecteurs du Simplex on la
conserve, sinon on conserve xr et l'on passe à l'itération suivante.
Si f (xr ) est plus grande que la valeur maximale des vecteurs du Simplex, on tente une compression
(voir gure 6.3) suivant la direction de la projection.
6.2. SIMPLEX DE NELDER & MEAD. 49
xc = βxmax + (1 + β)xk , avec β = 1/2, compression, (6.3)
Simplex après compression

Simplex original
B B
B B
B B
HH B HH B
xc

1

` ```
HHB H`B
H
Figure 6.3 Image d'un Simplex de 4 vecteurs avant et après compression.

Si f (xc ) est plus grande que la plus grande valeur associée des autres vecteurs du Simplex, on procède
alors à une contraction (voir gure 6.4) de tous les vecteurs du Simplex en direction du vecteur donnant
le critère le plus faible (xkmin ).
ix + xkmin
i
xk+1 =
k
, contraction. (6.4)
2
Simplex après contraction

Simplex original
B
B B
B H B
HH B
HB

1
}
Z
H H
ZB
Figure 6.4 Image d'un Simplex de 4 vecteurs avant et après contraction.
Initialisation. L'initialisation du Simplex peut être quelconque à condition que les vecteurs soient
linéairement indépendants. L'initialisation classique est de la forme :
xi = x0 + δui , i = 1..n, (6.5)
où : 

 x0 = : vecteur initial quelconque,
xi = : vecteurs du Simplex,


 δ= : paramètre réel,
ui = :

vecteurs linéairement indépendants (en général normés).
Les vecteurs ui forment donc une base de l'espace des paramètres.
Critère d'arrêt . Le critère d'arrêt de l'optimisation est :
λmax − λmin
s= < (6.6)
λmax + λmin
étant une valeur réelle positive choisie au départ.
Le choix de est limité par :
La borne supérieure est liée à la valeur approximative du critère optimal. L'application de 6.6 pour
une précision relative voulue de 1% sur un critère optimal environ égal à 10 conduit au choix de
= 5 × 10−3 .
La borne inférieure est liée à la précision du calculateur. On montre que doit être choisi comme
égal à la racine carrée du plus petit nombre qui peut être représenté, ce qui donne, pour des calculs
p
eectués en virgule ottante et en double précision, une limite inférieure de = 1, 7 × 10−308 =
1, 3 × 10−154 .
Figure 6.5 Ordinogramme de l'algorithme du Simplex de Nelder et Mead.
6.2.1 Extensions - Améliorations.
Les améliorations sont souvent liées à la vitesse de convergence de l'algorithme. Notons parmi les plus
connues :
La méthode du barycentre pondéré (Weighted Centroïd) :

Le point rééchi est calculé par réexion non pas par rapport à l'isobarycentre mais par rapport au
barycentre pondéré par l'inverse de la valeur du critère. Ainsi le point rééchi tend à se rapprocher
des points ayant le plus faible critère.
La méthode du Simplex super modiée (Super Modied Simplex) : : Après calcul du critère au point
rééchi et du critère du barycentre, nous faisons passer une équation du second degré par les trois
points (original, le barycentre et le point rééchi) dont on détermine l'optimum. C'est ce dernier
point qui sera choisi comme point rééchi.
L'introduction d'une variable aléatoire dans la recherche de nouveaux points rend l'optimisation plus
globale. Cette amélioration sort déjà du cadre des méthodes d'ordre 0 pour entrer dans le cadre des
Algorithmes Génétiques (voir 6.3) ou des méthodes du recuit simulé (voir 6.5).
6.2. SIMPLEX DE NELDER & MEAD. 51
Conclusions.
Cette méthode est la plus robuste des méthodes de type analytique. Les capacités d'adaptation du
volume du Simplex ainsi que l'absence de calcul du gradient en font une méthode rapide, peu sensible
aux discontinuités de la fonction de coût et facile à mettre en ÷uvre.
6.3 Algorithmes génétiques.
6.3.1 Introduction.
Développés par Holland en 1970, ces algorithmes sont la transposition algorithmique de la théorie de
Ch. Darwin sur la survivance du plus fort ou la sélection naturelle des espèces. Selon Ch. Darwin,
les individus les plus faibles disparaissent, seuls les individus les plus forts se reproduisent et parfois
ceux-ci mutent.
6.3.2 Principe des algorithmes génétiques binaires.
Les algorithmes génétiques binaires travaillent sur une entité fondamentale appelée chromosome qui
est en fait une chaîne de bits, l'unité d'information la plus simple. L'utilisation de cette structure
fondamentale liée au calculateur et non pas au processus à optimiser permet une grande adaptabilité
aux diérents problèmes (combinatoires, optimisation de critères fonction de variables réelles). Ainsi
le vecteur de paramètres initial est transformé en une chaîne de bits.
x = (α1 , α2 , α3 , · · · , αn ) −→ x = (b1 , b2 , b3 , b4 , b5 , · · · , bm ) (6.7)
où :
x : est un individu ou chromosome, (6.8)
bi : sont les gènes. (6.9)
Le passage de l'un à l'autre fera l'objet du codage.

L'algorithme travaille sur une population de taille constante (30 à 300 individus typiquement) qui
évolue à chaque itération, pour créer une nouvelle génération.
Initialisation. La population initiale est créée aléatoirement.
Sélection. La sélection des individus qui vont se reproduire est aléatoire mais fonction de la valeur
de leur critère associé. L'individu j est sélectionné suivant une probabilité p sj
f ( j x)
psj = Pn i
. (6.10)
i=1 f ( x)
Il existe une autre possibilité : la sélection des parents à croiser suivant leur rang (après classement
de l'ensemble de la population suivant la valeur de leur critère associé) et non la valeur de leur critère
associé. Cette deuxième méthode présente l'avantage d'empêcher une convergence trop rapide vers un
individu nettement meilleur que les autres.
Croisement. Le croisement est l'étape principale de l'optimisation. Cette opération consiste à en-
gendrer 2 enfants ayant un chromosome composé de parties des chromosomes des 2 parents.
Exemple avec 2 points de croisement (k) :

parent1 = 001k00110010k1011 
enf ant1 = 001k10011000k1011
⇑ ⇓ −→ (6.11)
enf ant2 = 101k00110110k1101
parent2 = 101k10011100k1101

Si la position des points de croisement était constante, combiné à la sélection, le croisement aurait pour
eet d'optimiser chaque partie de la chaîne. Mais la position des points de croisement est aléatoire, ce
qui participe à l'exploration de l'espace total. La probabilité de croisement, c'est-à-dire la probabilité
que deux individus sélectionnés pour le croisement subissent eectivement le croisement est en général
de pc = 0, 5 à 0, 7.
6.4. ALGORITHMES GÉNÉTIQUES CODÉS RÉELS. 53
Mutation. La mutation d'un individu permet de réaliser l'aristogénèse et sert ainsi à empêcher la
perte de diversité. En eet si au cours de l'optimisation un gène de rang donné avait la même valeur
pour toute une génération, tous les descendants conserveraient ce gène. La probabilité de mutation
est en général de pm = 0, 01 à 0, 05.
individu = 001100k1k10001011 −→ individu = 001100k0k10001011 (6.12)
Codage. Les fonctions que l'on veut minimiser étant dans la grande majorité des cas des fonctions
de variables réelles, un codage des chaînes binaires est nécessaire. Le codage le plus utilisé est le codage
binaire à virgule xe, par exemple :
1001, 1010 −→ 1 × 23 + 0 × 22 + 0 × 21 + 1 × 20 + 1 × 2−1 + 0 × 2−2 + 1 × 2−3 + 0 × 2−4 = 9, 625 (6.13)
La longueur de la chaîne dépend donc de la précision souhaitée sur les valeurs réelles. Certains auteurs
préfèrent le code de Gray qui présente l'avantage d'avoir des codes consécutifs toujours distants de 1,
Mais ce codage diminue l'aptitude à l'exploration de l'espace total par le croisement. Ces méthodes
de codage donnent une grande facilité pour imposer des contraintes sur les valeurs d'un paramètre.
Les algorithmes génétiques sont actuellement en pleine expansion. Le principal problème actuel des
algorithmes génétiques étant la faible propension à l'exploration de l'espace complet, les recherches
portent en particulier sur ce point.
Longueur de la chaîne de gènes dynamique.
Notion de niche écologique.
Sur les croisements :
croisement multipoint (3,4,..) au lieu de 2 comme décrit précédemment,
croisement uniforme qui est un croisement bit à bit avec une probabilité d'échanger les bits
qui dépend du nombre d'échanges déjà eectués. Le croisement uniforme augmente la capacité
d'exploration de l'espace total, il est surtout utilisé dans le cas de populations faibles (dans le cas
de fortes populations, le nombre élevé d'individus assure déjà un bonne exploration de l'espace
total).
Conservation d'une génération sur l'autre du ou des meilleurs individus.
6.3.4 Conclusion.
Cette méthode est bien adaptée à la minimisation dans le cas de fonctions fortement non linéaires.
La mutation est surtout utile dans le cas d'un critère à forte épistasie, le croisement étant plutôt
utile dans la phase de détermination des minima locaux. Les algorithmes génétiques sont une bonne
méthode intermédiaire entre les méthodes analytiques et les méthodes stochastiques.
6.4 Algorithmes génétiques codés réels.
An de résoudre le problème du codage et de la précision des paramètres dans les algorithmes géné-
tiques, certains auteurs ont transposé les idées des algorithmes génétiques à des gènes de type réel.
6.4.1 Principe.
La séquence des opérations reste identique à celle des algorithmes génétiques binaires, seules les opé-
rations et le codage dièrent.
Un chromosome est composé d'un ensemble de gènes réels de la forme :
chromosome = (α1 , α2 , α3 , α4 , α5 , · · · , αn ), gènes = αi , i = 1..n. (6.14)

Croisement. Le croisement, qui consiste à engendrer des descendants par échange de gènes, de-
vient :

parent1 = (α1 , kα2 , α3 , kα4 , α5 ) 
enf ant1 = (α1 , kβ2 , β3 , kα4 , α5 )
⇑ ⇓ −→ (6.15)
enf ant2 = (β1 , kα2 , α3 , kβ4 , β5 )
parent2 = (β1 , kβ2 , β3 , kβ4 , β5 )

ou

parent1 = (α1 , kα2 , α3 , kα4 , α5 ) 
enf ant1 = (α1 , k α2 +β2 α3 +β3
2 , 2 , kα4 , α5 )
⇑ ⇓ −→ α2 +β2 α3 +β3 (6.16)
enf ant2 = (β1 , k 2 , 2 , kβ4 , β5 )
parent2 = (β1 , kβ2 , β3 , kβ4 , β5 )

Mutation. Là encore, la mutation consiste à empêcher la perte de matériel génétique. Transposé

au codage réel, cela revient à ajouter un bruit sur un ou plusieurs chromosomes choisis aléatoirement.
αi0 = αi ± p × αi . (6.17)
où p est un nombre aléatoire de moyenne nulle.

Renders propose une mutation dont l'amplitude décroît avec le temps, de la forme suivante :
5
1− kk
αi0 = αi + (αimax − αi ) 1 − r 0 si b = 1,
5 (6.18)
0 1− kk
αi = αi + (αimin − αi ) 1 − r 0 si b = 0,
où :
αimax : valeur maximum du paramètre αi ,
αimin : valeur minimum du paramètre αi ,
r : nombre réel aléatoire donnant l'amplitude de la mutation,
k : indice de la génération courante,
k0 : indice de la dernière génération subissant une mutation,
b : nombre binaire aléatoire donnant le sens de la mutation.
Par analogie avec les algorithmes de type recuit simulé (voir 6.5), l'amplitude des mutations décroît
au fur et à mesure que l'algorithme converge.
Améliorations. An d'améliorer la vitesse de convergence, il est courant de mettre en ÷uvre une
combinaison de plusieurs méthodes d'optimisation : algorithmes génétiques et méthodes analytiques.
L'algorithme génétique recherche les minima locaux et la méthode analytique détermine rapidement
la valeur du minimum de la fonction dans chaque minimum local.
6.4.2 Conclusion.
Les algorithmes génétiques réels présentent l'avantage de ne pas être limités par le codage des variables
réelles et la forme générale de l'algorithme permet un passage aisé des algorithmes génétiques vers des
méthodes de type Simplex.
6.5. MÉTHODE DU RECUIT SIMULÉ 55
6.5 Méthode du Recuit simulé
6.5.1 Introduction
La méthode du recuit simulé ( Simulated Annealing) repose sur un principe de la thermodynamique :

la loi de refroidissement des métaux. Un atome donné est animé par un mouvement brownien dont
l'amplitude dépend de la température. Lorsque la température décroît ce mouvement s'estompe et
l'atome tend à se positionner dans un état d'énergie minimale.
6.5.2 Principe.
Le problème posé est la détermination d'un vecteur x qui minimise le critère λ(x). La méthode du
recuit simulé la plus simple résout le problème par un processus itératif et stochastique comme suit.
Après une initialisation aléatoire du premier point de calcul x0 , l'algorithme génère, à chaque itération,
un pas aléatoire δx suivant une loi de probabilité de type gaussienne g(δx).
1 −(δx)2
g(δx) = e 2T (k) (6.19)
(2πT (k))−n/2
Le nouveau point de calcul est :
xk+1 = xk + δx. (6.20)
Deux cas peuvent se présenter :
1. Le nouveau point est meilleur : λ(xk+1 ) < λ(xk ),

le point est accepté et le processus recommence.
2. Le nouveau point est plus mauvais : λ(xk+1 ) > λ(xk ),

le nouveau point peut être accepté suivant une probabilité :

λ(xk+1 ) − λ(xk )
pk = exp − (6.21)
T (k)
A chaque itération la température T (k) diminue suivant la loi :
T0
T (k) = . (6.22)
ln(1 + k)
Une valeur minimale de T0 est proposée par Hajcek qui assure une convergence statistique vers un
minimum global.
Ainsi lorsque la température est grande, l'algorithme peut quitter un minimum local, soit par un saut
direct, soit par un point intermédiaire de critère plus grand. Au fur et à mesure que la température
décroît, la longueur du pas diminue et la probabilité d'accepter un pas δx qui fait augmenter le
critère décroît. L'algorithme nit par générer des pas de faible amplitude et n'accepter pratiquement
plus de pas qui font augmenter le critère, les pas qui le font diminuer étant toujours acceptés. Cette
méthode permet statistiquement de trouver un minimum global mais demande un très grand nombre
d'itérations.
Par extension, de nombreuses méthodes se rapportent à celle-ci du moment que la probabilité d'ac-
ceptation du pas est de la forme de celle présentée par l'équation (6.21).
Parmi les améliorations proposées dans la littérature, notons celle appelée : Fast Annealing où g(δx)
et T (k) sont choisis tels que :
T (k)
g(δx) = (n+1)/2
(6.23)
((δx)2 + T (k)2 )
et
T0
T (k) = (6.24)
k
L'auteur propose donc une loi de décroissance plus rapide, tout en conservant la mobilité du point de
calcul.
D'autre part, notons aussi la modication de [5] qui introduisent dans la méthode du Simplex une
probabilité d'acceptation du pas. Cet algorithme présente l'avantage de converger vers l'algorithme
classique du Simplex (voir [5]) lorsque la température décroît.
6.5.4 Mise en ÷uvre du recuit simulé.
Pour la mise en ÷uvre de la méthode du recuit simulé, nous avons choisi l'algorithme proposé par
[5].
L'algorithme présente trois phases distinctes mais imbriquées :
1. lorsque la température est grande vis-à-vis du critère, pratiquement tous les déplacements sont
acceptés, le Simplex s'agrandit et explore l'espace,
2. lorsque la température est du même ordre de grandeur que le critère, le Simplex converge vers
les grandes vallées tout en acceptant des déplacements qui augmentent le critère, ce qui lui
permet de sortir des rugosités de celui-ci,
3. lorsque la température est faible vis-à-vis du critère, l'algorithme détermine le minimum local
qui, selon toute probabilité, est le minimum global.
Choix de la température. Il est nécessaire de choisir une température initiale T0 et une loi de
décroissance au cours du temps. B. Hajcek a démontré qu'une décroissance de la forme :
T0
T (k) = (6.25)
ln(1 + k)
associée à une probabilité d'acceptation du pas de la forme :

λ(xk+1 ) − λ(xk )
pk = exp − (6.26)
T (k)
assure statistiquement la détermination d'un minimum global à condition de choisir T0 supérieur au

maximum de profondeur de tous les minima locaux, minimum global exclu.
Le maximum de profondeur de tous les minima locaux est impossible à déterminer, c'est pourquoi
certains auteurs [5] préconisent l'emploi d'une valeur largement supérieure à une approximation de ce
second membre.
Nous avons borné cette profondeur maximale par :
T0 = 10 × la valeur maximale du critère à l'initialisation, (6.27)
la valeur minimale du critère étant minimisée par 0.

Bibliographie
[1] C. Foulard, S.Gentil, and J.P. Sandraz. Commande et régulation par calculateur numérique : de la
théorie aux applications. Eyrolles, Paris, 1987.
[2] Lennart Ljung. System identication (2nd ed.) : theory for the user. Prentice Hall PTR, Upper
Saddle River, NJ, USA, 1999.
[3] Lennart Ljung and Torkel Glad. Modeling of dynamic systems. Prentice-Hall, Inc., Upper Saddle
River, NJ, USA, 1994.
[4] Lennart Ljungs. System identication toolbox for use with matlab, 1997.
[5] William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical
Recipes : The Art of Scientic Computing. Cambridge University Press, Cambridge (UK) and New
York, 2nd edition, 1992. Disponible gratuitement sur Internet.
Liens
http://www.mathworks.com/products/fuzzylogic.html
http://www.mathtools.net/MATLAB/Fuzzy_Logic/index.html
http://www.dbai.tuwien.ac.at/marchives/fuzzy-mail/
http://europa.eu.int/comm/research/rtdinfo/fr/24/02.html
http://www.gala.univ-perp.fr/~polit/chap0.html
http://www.eru.ulaval.ca/ptt15225
http://vcampus.u-strasbg.fr/uticeweb/mhiri/projet_V2
http://www.mathworks.fr/products/controldesign/modanal.shtml
57
58 BIBLIOGRAPHIE
Chapitre 7
Travaux dirigés
TDAO : Illustration de plusieurs aspect de l'identication des systèmes
7.1 Strecj - Broida - Programmation non linéaire
7.1.1 Strecj - Broida
Ouvrez le chier simulink intitulé "TDAOident.mdl", lancez la simulation,

Appliquez la méthode de Broïda puis celle de Strecj et implantez-les,
relancez la simulation, génial non ?
7.1.2 Programmation non linéaire
Ouvrez le chier "pnl.m" et comprenez le contenu.

Dans la fenêtre de commande matlab, tapez "save ES Uin tout Yout" puis,
"X=fminsearch('pnl',[-1 -1 -1 1])" (tapez "help fminsearch" et lisez la doc.)
la valeur de X présente le résultat sous la forme "3 pôles, gain"
implantez ce nouveau modèle et relancez la simulation.
7.1.3 Limites de la Programmation non linéaire
lancez "X=fminsearch('pnl',[0 0 0 0])", comparez les solutions.

Dans le chier "pnl.m" changez l'utilisation de x0 en introduisant des zéros et plus de pôles.
7.2 Identication des systèmes instables ou intégrateurs
Ouvrez le chier simulink intitulé "instable.mdl". Ce système est instable mais stabilisé par le gain de
0.1 et le bouclage. Les entrées - sorties sont dans le "workspace" de Matlab et s'importent facilement
dans l'interface ident (Tapez "ident" dans la fenêtre de commande de Matlab), voir gure 7.1.
Note : les modèles calculés dans "ident" s'exportent vers le workspace par glisser déplacer dans la case
"to workspace" puis du workspace vers simulink par une commande comme :
[A1,B1,C1,D1]=ssdata(amx2122)
Observez les caractéristiques spectrales de la SBPA et comparez à celles du signal carré.
Après avoir annulé le bruit, importez les données puis calculez un modèle de type ARX752.
observez la qualité de la simulation.
observez les pôles et les zéros de ce modèle.
déduisez-en un modèle plus simple et tout aussi correct
vérier que vous obtenez bien les coecients du modèle en boucle fermée
Observez la diérence d'estimation des paramètres du système en fonction de :
la puissance de bruit (0 0,01 et 0,0001)
59
60 CHAPITRE 7. TRAVAUX DIRIGÉS
Figure 7.1 Comment importer des variables dans ident de Matlab.
les diérents modèles (en particulier, ARX, ARMAX, BOX-JENKINS)

les signaux utilisés pour l'estimation (carré ou SBPA)
la position de la source d'excitation (entrée de la boucle ou juste devant le système à identier)
Reprenez ce qui précède en prenant comme entrée la sortie du correcteur.
7.3 Choix de la source d'excitation
L'objectif est de déterminer la bonne séquence binaire pseudo-aléatoire pour exciter correctement le
système.
Ouvrez le chier "sbpa_ident.mdl"
Le réglage de la période de la SBPA s'eectue en déclarant la variable TeSBPA=1 (ici 1 seconde).
Déterminer un premier modèle (ARX441 par exemple).
Déduisez-en les constantes de temps du système.
Déterminer un SBPA susceptible de bien exciter le système dans toute la plage de fréquences que
vous voulez identier.
Identier à nouveau le système.
7.4 Les modèles paramétriques
Il s'agit de prédire le cours d'une action du CAC40 à partir des cours des autres actions.
Exécuter le chier "trait_don".
Importez les variables 'entree" comme entrée et "sortie" comme sortie
Calculez et observez les réponses des diérents modèles (ARX, ARMAX, BOX-JENKINS...)
En particulier, notez l'erreur de prédiction 5 pas en avance
Peut-on devenir riche ? Si oui comment
Valider vos modèles à l'aide des données 'entreeV" comme entrée et "sortieV"
7.5 Cas pratique
Le système modélisé dans le chier "MEMSswitch.mdl" contient un modèle "réaliste" d'un interrupteur
fabriqué en technologies MEMS que l'on désire utiliser comme un bras de levier commandé en position.
L'objectif est de déterminer un modèle paramétrique du système.
L'ordre de grandeur de l'entrée est de 50V
Le comportement est du type second ordre.
Les caractéristiques dynamiques varient un peu en fonction de la tension de commande continue.
7.5. CAS PRATIQUE 61
La mise en place des données physiques du systèmes est faite par le chier "donneesMEMS.m".
Débrouillez-vous !
62 CHAPITRE 7. TRAVAUX DIRIGÉS
Deuxième partie
Annexes
63
JANVIER 2004 65
3ème année
Option automatique
Examen d'Identication et Commande Floue
Cet examen consiste à modéliser et à commander le système de commande d'altitude d'un avion. Il
s'agit dans un premier temps de déterminer un modèle simple du système de façon graphique, puis de
commander le système par un régulateur ou. Enn, le système étant corrigé par ce régulateur ou
on déterminera un modèle paramétrique plus précis.
1 Méthodes graphiques
Step Response Step Response

From: U(1) From: U(1)
5 0.35
0.3
4
0.25
3
0.2
Amplitude
Amplitude
To: Y(1)
To: Y(1)
2 0.15
0.1
1
0.05
0
0
−1 −0.05
0 1 2 3 4 5 6 0 0.5 1 1.5 2 2.5
Time (sec.) Time (sec.)
Figure 2 Réponses du système à un échelon : à gauche en boucle ouverte, à droite en boucle fermée.
c + lε correcteur
u système
s
Figure 3 Identication en boucle fermée

1. La réponse à un échelon du système en boucle ouverte est donnée sur la gure 2 (gauche). En
faisant l'hypothèse que le système est un système du premier ordre avec un retard pur, donnez
les valeurs numériques d'un tel modèle.
2. Toujours sur cette même gure appliquez la méthode de Strejc an d'obtenir un modèle plus
précis.
3. Ce même système est bouclé par un retour unitaire et un gain de 0,1 en guise de correcteur
comme indiqué sur la gure 3. Sa réponse à l'échelon unité est donnée sur la gure 2(droite).
(a) Déterminer un modèle du système en boucle fermée.
(b) En négligeant le retard pur, déduisez-en un modèle en boucle ouverte.
4. Finalement, on insère un gain variable dans la boucle et il s'avère que le système oscille pour un
gain k = 1.
(a) Ce résultat est-il concordant avec vos modèles (raisonnement qualitatif !) ?

66 ANNEXE .ANNALES
(b) Au vu des résultats précédents, proposez un modèle à trois constantes de temps
5. Ce système présente de façon évidente un zéro dit instable !
(a) Proposez alors un modèle du système incluant ce zéro instable.
(b) Donnez une valeur numérique très approchée (pratiquement sans calcul) de ce zéro agré-
menté d'une ligne de commentaire sur ce choix.
2 Commande par contrôleur ou
Finalement, le système précédent est commandé par un contrôleur ou.
µ(ε) µ(∆ε) µ(u)

N Z
6 P N Z
6 P NG N Z
6 P PG
ε
- ∆ε
- u -
-20 0 20 -10 0 10 -10 -5 0 5 10
Figure 4 Fuzzication des variables

ε \ ∆ε N Z P
opérateurs logiques : ET = min OU = max
N NG N Z
avec méthode d'inférence : somme-produit
Z N Z P
défuzzication : moyenne des maximums.
P Z P PG
1. En supposant ε = 10 et ∆ε = 0, calculez la sortie u du contrôleur.
2. Toujours en supposant ∆ε = 0, calculez la sortie u du contrôleur en fonction de la valeur de ε.

Tracez la courbe u = f (ε) en précisant les points importants et la forme de la courbe entre ses
points (droites, paraboles, ...)
3. Même question mais en utilisant les hypothèses suivantes :

méthode d'inférence : max-min
défuzzication : centre de gravité
3 Modèles paramétriques
Le système étant corrigé par le régulateur ou précédent, on décide de déterminer un modèle paramé-
trique du système corrigé. Dans ce qui suit, et quels que soient vos résultats précédents, on supposera
que le système présente :
3 constantes de temps de l'ordre de : 0,1, 0,05 et 0,02
1. Calculer la fréquence d'échantillonnage nécessaire pour pouvoir identier correctement ces trois
constantes de temps.
2. Calculer la longueur de la séquence binaire pseudo aléatoire en fonction de cette fréquence

d'échantillonnage pour bien identier le gain du système.
3. Quelle est alors le temps d'excitation du système ?

3. MODÈLES PARAMÉTRIQUES 67
4. On opte pour une méthode d'estimation des paramètres dite de la matrice instrumentale, dénie
par :
 ∧



 θ, Paramètres estimés ;
∧  Z,
 Matrice instrumentale ;
T −1 T
θ = (Z X) Z y avec y, Vecteur des sorties mesurées ;

X, Construite à partir du jeu d'E/S (y = Xθ + e) ;




e, Résidus.

Montrez que pour que cet estimateur soit sans biais, il faut que
E[Z T X] 6= 0
E[Z T e] = 0
68 ANNEXE .ANNALES
3ème année
Option automatique
Corrigé de l'Examen d'Identication et Commande Floue
Step Response Step Response

From: U(1) From: U(1)
5 0.35
0.3
4
0.25
3
0.2
Amplitude
Amplitude
To: Y(1)
To: Y(1)
2 0.15
0.1
1
0.05
τ 0
Ta
T
−1 −0.05
Tu
0 1 2 3 4 5 6 0 0.5 1 1.5 2 2.5
Time (sec.) Time (sec.)
Figure 5 Réponses du système à un échelon : à gauche en boucle ouverte, à droite en boucle fermée.
1. Modèle du premier ordre
On mesure T = 0, 3s et τ = 1, 9 − 0, 3 = 1, 6s le modèle est donc :
Ae−T p 5e−0,3p
G(p) = =
1 + τp 1 + 1, 6p
2. Modèle Strejc en boucle ouverte

On mesure T u = 0, 3s et T a = 1, 6s. En appliquant la méthode :
T u/T a = 0, 1875 ⇒ n = 2
T a/τ = 2, 718 ⇒ τ = T a/2, 718 = 0, 5887
T u/τ = 0, 282 ⇒ T u0 = 0.1660
T r = T u − T u0 ⇒ T r = 0.134
d'où le modèle :
Ae−T rp 5e−0,166p
Gs (p) = =
(1 + τ p)n (1 + 0, 5887p)2
3. En boucle fermée
(a) Modèle du système en boucle fermée,

On mesure T u = 0, 33s et T a = 1s. En appliquant la méthode :
T u/T a = 0, 330 ⇒ n = 4
T a/τ = 4, 463 ⇒ τ = T a/4, 463 = 0, 2241
T u/τ = 1, 425 ⇒ T u0 = 0, 319
T r = T u − T u0 ⇒ T r = 0.01
2. COMMANDE PAR CONTRÔLEUR FLOU 69
d'où le modèle :
Ae−T rp 0, 34e−0,01p
Gs (p) = =
(1 + τ p)n (1 + 0, 224p)4
(b) Modèle en boucle ouverte déduit :
A/k
Gs (p) =
τ 4 p4 + 4τ 3 p3 + 6τ 2 p2 + 4τ p + 1 − A
4. Finalement, on insère un gain variable dans la boucle et il s'avère que le système oscille pour un
gain k = 1.
(a) A priori oui, les deux modèles obtenus sont susceptibles d'être instables dans une boucle
fermée si le gain est susamment grand.
(b) le plus simple est de reprendre le modèle obtenu par la méthode de Strejc sur la boucle
ouverte :
Ae−T rp A
Gs (p) = n
' 2
(1 + τ p) (1 + τ p) (1 + T rp)
5. Introduction du zéro instable
(a) Encore une fois, en reprenant le modèle précédent et en développant au premier ordre
l'exponentielle, on obtient :
Ae−T rp A(1 − T rp)

Gs (p) = n
'
(1 + τ p) (1 + τ p)2
(b) T r = 0.166. En fait pour tenir compte des termes que l'on néglige on augmente "un peu"
la valeur de Tr. Ainsi T r0 = 1, 5T r donne de bien meilleurs résultats.
2 Commande par contrôleur ou
1. Si ε = 10 et ∆ε = 0, la sortie vaut 2,5.
2. En supposant ∆ε = 0, calculez la sortie u du correcteur en fonction de ε

u
6
5
2,5
-20 -10 10 20
-ε
-2,5
-5
Figure 6 u = f (ε)(courbe exagérée !)
3. avec :
méthode d'inférence : max-min
défuzzication : centre de gravité
70 ANNEXE .ANNALES
3 Modèles paramétriques
Le système présente 3 constantes de temps de l'ordre de : 0,1, 0,05 et 0,02.
1. Une SBPA est un bruit blanc sur 0, 3fe donc :
1 1 1
0, 3fe ≥ ⇒ fe ≥ = = 26, 5Hz
2πτmin 0, 3 × 2πτmin 0, 3 × 2π × 0.02
2. Le plateau maximum est de N Te avec Te = 1/fe et 2N − 1 la longueur de la SBPA.

Il faut donc que :
N Te ≥ 3 à 5τmax ⇒ N ≥ 3 à 5fe τmax = 3 à 5 × 26, 5 × 0, 1 = 8 à 13
N = 10 semble un bon compromis.
3. Le temps d'excitation est :
Texcitation = (2N − 1)Te = 1023Te = 38, 6s
mais en général on envoie deux séquences consécutives pour éviter de prendre en compte le
transitoire au début de la séquence, donc Texcitation ' 77, 2s
4. Méthode de la matrice instrumentale :
 ∧



 θ, Paramètres estimés ;
∧  Z,
 Matrice instrumentale ;
T −1 T
θ = (Z X) Z y avec y, Vecteur des sorties mesurées ;

X, Construite à partir du jeu d'E/S (y = Xθ + e) ;




e, Résidus.

∧
T −1 T
θ = (Z X) Z (Xθ + e)
= θ + (Z T X)−1 Z T e
∧
pour que E[θ ] = θ, il faut que :
E[(Z T X)−1 Z T e] = 0
donc que :
E[Z T X] 6= 0
E[Z T e] = 0
1. MÉTHODES GRAPHIQUES 71
3ème année
Option mécatronique
Examen d'Identication
Cet examen consiste à modéliser un bras de robot présentant deux types de dérives lentes. En premier
lieu, la constante de temps principale varie dans le temps puis le système présente une dérive lente de
la valeur moyenne de sortie.
40 10
9
35
8
30
7
25
6
20 5
4
15
3
10
2
5
1
0 0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
Figure 7 Réponse du système à un échelon au début de l'identication (gauche) et après une heure
(droite).
14
12
10
−2
−4
0 10 20 30 40 50 60 70 80 90 100
Figure 8 Réponse du système à une excitation de type SBPA

72 ANNEXE .ANNALES
2 Modèles de Strecj
1. La réponse à un échelon du système en boucle ouverte est donnée sur la gure 7(gauche). En
faisant l'hypothèse que le système est un système du premier ordre avec un retard pur, donnez
les valeurs numériques d'un tel modèle.
2. Toujours sur cette même gure appliquez la méthode de Strejc an d'obtenir un modèle plus
précis.
3. Appliquez Strecj sur la gure 7(droite) et déduisez-en un nouveau modèle du système.
4. En faisant l'hypothèse d'une variation linéaire de la constante de temps dominante et à l'aide

des deux modèles précédents, proposez un modèle du système H(p, t)
3 Modèles ARX
An d'étudier plus précisément les autres constantes de temps du système, vous décidez de soumettre
le système à une séquence binaire pseudo-aléatoire. Le résultat est donné g. 8
1. Rappeler le principe de calcul des modèles ARX et la formule appliquée pour obtenir une esti-
mation des paramètres du système.
2. La première tentative de modélisation ARX est très mauvaise ! Quels sont les deux phénomènes
qui "empêchent" d'avoir un modèle ARX correct.
3. An de compenser la dérive lente de la valeur moyenne, vous envisagez deux méthodes
(a) Boucler le système, et faire une identication en boucle fermée,
(b) Enlever cette dérive par un traitement numérique hors ligne (un peu comme on enlève la
valeur moyenne).
4. Proposez un correcteur simple pour le bouclage du système
5. Donnez la relation entre y(t) valeurs de la gure 8 et y1 (t) qui représente y(t) sans la dérive. Les
valeurs numériques approchées seront prises sur la gure 8.
6. Après calcul des modèles ARX, les modèles restent passablement mauvais avec un léger avantage
pour le modèle calculé à partir des valeurs en boucle fermée, pourquoi ?
7. Les méthodes précédentes ne donnant pas de modèles susants pour l'application envisagée,
proposez une méthode qui permet d'identier ce système malgré toutes ses imperfections.
4. FUZZIFICATION DES VARIABLES 73
3ème année
Examen de commande oue
Cet examen est destiné à déterminer un estimateur ou de pluie en fonction de la pression atmosphé-
rique, du taux d'hygrométrie et de la température. Les valeurs numériques seront estimées grossière-
ment à partir des graphiques donnés où les abscisses sont en jours.
4 Fuzzication des variables
1. Après avoir observé les relations entre quantité de pluie et les variables d'entrée en gure
??(pression atmosphérique, taux d'hygrométrie, température et leurs dérivées ) peu de temps
avant le début des précipitations, proposez un choix d'ensembles ous pertinents. justiez briè-
vement le choix du nombre d'ensembles ous par variable et la forme de ses ensembles.
2. La forme de la variable de sortie est donnée gure 9. Etablissez les règles d'inférence pour chaque
précipitation.
L'estimateur s'avère être trop cher (trop de capteurs) et votre cahier des charges devient donc "esti-
mateur de pluie uniquement fonction de la pression atmosphérique et de ses variations".
1. proposez, dans ce cas, les règles d'inférence sous la forme d'un tableau :
p \ ∆p N Z P
N
Z
P
2. Proposez une méthode pour chaque point suivant :

opérateurs logiques : ET = ? OU = ?,
méthode d'inférence,
défuzzication.
3. Après un premier test en simulation il s'avère que votre estimateur prédit correctement dans
70% des cas mais le service marketing vous explique qu'il serait bon de donner outre le fait qu'il
pleuve ou pas une estimation de la probabilité d'orage. L'orage est probable s'il y a une baisse
de 1 à 2 hPa de la pression atmosphérique en une heure. Proposez un deuxième estimateur ou
donnant cette indication.
4. Question subsidiaire pour ceux qui ont suivit le cours d'identication. D'après vous, un modèle
de type ARX serait-il meilleur ou moins bon pour ce type de prévisions ?
fonction d'appartenance
6
beau variable pluie
pluviométrie
-
0 30 60 100
Figure 9 Fuzzication de la variable de sortie

74 ANNEXE .ANNALES
Figure 10 Données météorologiques à St-Philibert (abscisses en jours)

1. MODÈLE DE CONNAISSANCES 75
3ème année
Janvier 2006
Examen d'Identication des systèmes

Tous documents autorisés
Cet examen consiste à modéliser un système thermique particulièrement pénible à modéliser. Les
questions suivantes sont liées mais peuvent être traitées dans le désordre.
1 Modèle de connaissances
La mise en équation du comportement dynamique du système donne le schéma linéaire donné en

gure 11 avec les valeurs suivantes : R2 = 0, 2; R1 = 2, 1; C = 1; L = 1. 1 P est la puissance thermique
introduite dans le système et T est la température que l'on veut asservir.
R2 L
R1
+
T
C
P −
Figure 11 Modèle thermique du système.

T (p)
1. Déterminer la fonction de transfert linéaire continue (en p) du système
P (p) .
2. Calculer le facteur d'amortissement m du système.
3. Ce système s'avère trop peu amorti, on cherche alors à augmenter cet amortissement en agissant
sur la seule variable R1 (résistance thermique entre l'intérieur et l'extérieur de la chambre).
Comment doit varier R1 ?
2 Modèle non-paramétrique
Le modèle précédent étant particulièrement grossier, vous décider de soumettre le système à un échelon
unité démarrant à l'instant t = 0.
1. La réponse à un échelon du système en boucle ouverte est donnée sur la gure 12. En faisant
l'hypothèse que le système est un système du second ordre avec un retard pur, donnez la forme
générale du modèle.
2. Donnez les valeurs numériques des diérents coecients du modèle.
1. Les inductances thermiques n'existent pas ! Il s'agit en fait d'un système actif modélisable en première approxi-
mation par une inductance.

76 ANNEXE .ANNALES
1.4
1.2
1
0.91
0.8
0.6
0.4
0.2
0
0 2 4 6 8 10 12 14 16 18 20
Figure 12 Réponse indicielle du système.

bi
y
H(z) m i
vi
ui
Intercorrélateur
Figure 13 Identication par intercorrélation.
3 Modèle par intercorrélation
Les deux modèles précédents ayant donné quelques nouvelles connaissances sur le système, vous décidez
de mettre en ÷uvre une identication paramétrique fondée sur une identication par intercorrélation.
Dans cette partie les réponses demandent 1 ligne de calcul et une ou plusieurs interprétations du
P
résultat. Pour les pinailleurs mathématiques, et ils ont raison de l'être, on intervertira les signes et
lim sans démonstration.
∞
X
y(i) = (h(j) u(i − j)) + b(i)
j=0
avec :
y(i) = la sortie à l'instant i

h(j) = la séquence de pondération
u(i) = l'entrée à l'instant i

b(i) = le bruit à l'instant i
3. MODÈLE PAR INTERCORRÉLATION 77
Par dénition la fonction d'intercorrélation entre deux signaux y et u est :
N
1 X
Cyu (k) = lim y(i) u(i − k)
N →∞ 2N + 1
i=−N
De la même façon, la fonction d'autocorrélation d'un signal u est :
N
1 X
Cuu (k) = lim u(i) u(i − k)
N →∞ 2N + 1
i=−N
2
1. Quel est le lien entre H(z) et h(i) (cours des années précédentes)
2. Calculez Cyu (k) en fonction de l'entrée, du système et du bruit.
3. Introduisez dans ce résultat la fonction d'autocorrélation de l'entrée.
4. Montrez que si l'entrée et le bruit ne sont pas corrélés on obtient :
∞
X
Cyu (k) = h(j) Cuu (j − k)
j=0
5. Montrez que si l'entrée u est un bruit blanc parfait alors :
Cyu (k) = A h(k)
ou A est un coecient lié à la puissance du bruit.
6. Que proposez vous pour réaliser ce bruit blanc ?
7. La réponse du système est donnée en gure 14. Est-elle compatible avec celle présentée en gure
12.
0.7
0.6
0.5
0.4
0.3
0.2
0.1
−0.1
−0.2
−0.3
0 2 4 6 8 10 12 14 16 18 20
Figure 14 Réponse de l'intercorrélation entrée-sortie du système.

1
2. Un peu d'aide ? Calculer la sortie yk pour k variant de 0 à 5 de G(z) = z−0,5
soumis à δ(z)
78 ANNEXE .ANNALES
4 Filtrage des mesures
Tout comptes fait, vous décider de procéder à une identication paramétrique classique, et vous lancez
une campagne de mesures entrée-sortie du système. Malheureusement, le stagiaire qui devait faire cela,
trouvant que le signal de sortie était trop bruité n'a rien trouvé de plus malin que de ltrer les données
de la sortie avec un ltre passe bas.
1. Que pouvez-vous faire pour utilisez au mieux ces données ?
2. Essayer de le démontrer.
3. le premier résultat de calcul d'un modèle ARX 10 10 3 est donné en gure 15. Quel modèle
tenteriez-vous ensuite ?
Poles (x) and Zeros (o)
0.8
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
−1
−1 −0.5 0 0.5 1 1.5
Figure 15 Pôles et zéros du modèle ARX 10 10 3.

1. MODÈLE DE CONNAISSANCES 79
3ème année
Janvier 2006
Corrigé de l'Examen d'Identication
1 Modèle de connaissances
1. Fonction de transfert linéaire continue
T (p) R1
=
P (p) R1 LCp2 + (L + R1 R2 C)p + R1 + R2
2. Facteur d'amortissement
1 L + R1 R2 C
m= p
2 (R1 + R2 )R1 LC
3. R1 doit diminuer, il faut donc augmenter les pertes.
1. Forme générale du modèle.

Ae−τ p
H(p) = p2
2m
1+ w0 p + w0 2
2. Coecients du modèle.
A = 0, 91
ω0 = 1, 04
m = 0, 32
τ = 2
τ est "visible" mais s'obtient par calcul en appliquant :
√−πm
D% = 100e 1−m2 donne m
Tpseudo-période = w √2π1−z 2 donne w0
0
tpic − tpic = tpic − √π donne τ
mesuré calculé mesuré w0 1−z 2
3 Modèle par intercorrélation
1. h(i) est la séquence de pondération de la fonction de transfert H(z).

En temporel, tout système est un convoluteur : sa sortie est le produit de convolution entre
l'entrée et sa séquence de pondération.
En z, cette convolution est transformée en un produit (c'est pour cela qu'on utilise les trans-
formée de Laplace et en z !).
La séquence de pondération s'obtient facilement en calculant la sortie du système soumis à une
impulsion de Dirac.
2. Cyu (k)
80 ANNEXE .ANNALES
N
1 X
Cyu (k) = lim y(i) u(i − k)
N →∞ 2N + 1
i=−N
N ∞
1 X X
= lim u(i − k) h(j) u(i − j) + b(i)
N →∞ 2N + 1
i=−N j=0
∞ N N
X 1 X 1 X
= h(j) lim u(i − k)u(i − j) + lim b(i) u(i − k)
N →∞ 2N + 1 N →∞ 2N + 1
j=0 i=−N i=−N
∞ N N
X 1 X 1 X
= h(j) lim u(i)u(i − (j − k)) + lim b(i) u(i − k)
N →∞ 2N + 1 N →∞ 2N + 1
j=0 i=−N i=−N
3. Introduisez dans ce résultat la fonction d'autocorrélation de l'entrée
∞
X
Cyu (k) = h(j) Cuu (j − k) + Cbu (k)
j=0
4. Pour que :
∞
X
Cyu (k) = h(j) Cuu (j − k)
j=0
il faut donc que

Cbu (k) = 0
Ce qui est le cas si l'entrée est le bruit sont décorrélés.
5. Si l'entrée u est un bruit blanc parfait alors
N
1 X
Cuu (k) = lim u(i) u(i − k) = Aδ(k)
N →∞ 2N + 1
i=−N
Par dénition d'un bruit blanc, sa fonction d'autocorrélation est une impulsion de Dirac et donc :
Cyu (k) = A h(k)
ou A est la "puissance" du bruit, en fait c'est le carré de sa valeur ecace.
6. Que proposez vous pour réaliser ce bruit blanc ? Une séquence binaire pseudo aléatoire !
7. Oui, l'une est la réponse indicielle (entrée = échelon), l'autre la réponse impulsionnelle (entrée
= impulsion). Observez la pseudo-période, pour le gain, il faut un peu plus de calcul...
4 Filtrage des mesures
1. Que pouvez-vous faire pour utilisez au mieux ces données ?

Filtrer les entrées avec le même ltre !
2. Essayer de le démontrer.
Bonne chance ! (une démonstration est donnée dans le Ljung).
3. On observe facilement que 6 ou 8 zéros compensent pratiquement 6 ou 8 pôles. La prochaine

tentative serai donc ARX 4 4 3 ou ARX 2 2 3
1. MODÈLE NON-PARAMÉTRIQUE (5 PTS) 81
3ème année
Janvier 2007

1 Modèle non-paramétrique (5 pts)
Figure 16 Réponse du système à un créneau unité de durée t = 0, 5s.

1. Identier le système avec le méthode de Strejc, comme si la réponse était la réponse à un échelon.
2. Déduisez-en la fonction de transfert du système.
2 Modèle paramétrique (5 pts)
Un système sans bruit soumis à une rampe unité donne la sortie suivante.
k 0 1 2 3 4 5
uk 0 1 2 3 4 5
yk 0 0 7 7 42 −14
On se propose d'identier un modèle paramétrique de la forme
yk = −a1 yk−1 + b1 u2k−1

1. Ecrivez la relation précédente pour k=4 et k = 3.
2. Montrez que ces deux équations peuvent se mettre sous une forme matricielle de la forme :
y = Xθ + e
3. Déterminez les coecients a1 et b1 .

82 ANNEXE .ANNALES
3 Moindres carrés récursifs (5 pts)
estimation.
A l'instant N on peut écrire :
 
a1  

yN
 
−yN −1 −yN −2 . . . −yN −n uN . . . uN −p

 a2  eN
yN −1   −yN −2 −yN −3 −y u uN −p−1   
. 
eN −1 
N −n−1 N −1

 . 
 
 .   eN −2
  
yN −2 . .
 
. .
  
  −yN −3 . . uN −2 uN −p−2

.

  an 
 
yN −3 .
 
= + . (1)
  
. .
  
. .   b0 
 
. . −yn+1 . uN −p−3
 
. .
  
. .
  
  b1   .
   
.
 
.
  
yn+2   −yn+1 −yn . uN −p  .   
en+2

 .. 
 
yn+1 −yn −yn−1 . . . −y1 un+1 en+1
bm
sous forme condensée :
YN = XN θN + EN
dont l'estimation optimale est :

T
θbN = (XN XN )−1 XN
T
YN
1. Ecrivez la forme (1) à l'instant N +1 en utilisant

xN +1 yN +1 eN +1
XN +1 = , YN +1 = et EN +1 =
XN YN EN
2. Donnez le vecteur ligne xN +1

N + 1 s'écrit T
θbN +1 = (XN −1 T
3. L'estimation à l'instant : +1 XN +1 ) XN +1 YN +1 , écrivez θbN +1 en
fonction de XN , xN +1 , YN , yN +1 .
Note :
T
a a
= AT A + aT a
A A
4. Sachant que :
(A + BCD)−1 = A−1 − A−1 B(C −1 + DA−1 B)−1 DA−1
Montrez que :
−1
T
XN )−1 − (XN
T
XN )−1 xTN +1 1 + xN +1 (XN
T
XN )−1 xTN +1 T
XN )−1

T
YN + xTN +1 yN +1

XN
5. En posant
T X )−1 xT
a = 1 + xN +1 (XN N N +1 , montrez que :

T
θbN +1 = θbN + (XN XN )−1 xTN +1 a−1 yN +1 − xN +1 θbN (2)

6. Que représente le terme yN +1 − xN +1 θbN , interprétez alors l'équation (2)
4. MODÉLISATION PARAMÉTRIQUE (5 PTS) 83
Figure 17 Modèle issus d'un modèle de connaissances
4 Modélisation paramétrique (5 pts)
Un premier modèle issu d'un modèle de connaissance, est donné en gure 17.
1. Voyant que le système est intégrateur, comment procédez vous pour obtenir un modèle paramé-
trique du système ?
2. Déterminer les paramètres d'une séquence binaire pseudo aléatoire pour procéder à une identi-
cation paramétrique de ce système.
3. Quel modèle entre ARX et ARMAX donnera les meilleurs résultats d'après vous.
84 ANNEXE .ANNALES
3ème année
Janvier 2006
Corrigé de l'Examen d'Identication
1 Modèle non-paramétrique (5 pts)
1. Modèle Strejc
On mesure T u = 0, 3s et T a = 1, 3s. En appliquant la méthode :
T u/T a = 0, 2308 ⇒ n = 2
T a/τ = 2, 718 ⇒ τ = T a/2, 718 = 0, 4783
T u/τ = 0, 282 ⇒ T u0 = 0, 1349
T r = T u − T u0 ⇒ T r = 0, 1651
Gain : 0.83
d'où le modèle :
Ae−T rp 0.83e−0,134p
G(p) = =
(1 + τ p)n (1 + 0, 5887p)2
Figure 18 comparaison sortie réelle et sortie du modèle.

2. MODÈLE PARAMÉTRIQUE (5 PTS) 85
2. On écrit l'égalité : réponse indicielle du modèle = réponse à un créneau du système.
1 1 − e−0.5p Ae−T rp 1 2Ae−T rp

G(p) × = G2 (p) × =⇒ G2 (p) = '
p p (1 + τ p)n 1 − e−0.5p p(1 + τ p)n
2 Modèle paramétrique (5 pts)
On se propose d'identier un modèle paramétrique de la forme

42 −7 9 a1
= = y = Xθ + |{z}
e
7 −7 4 b1
=0
la solution est :

T −1 T 3
θb = (X X) X y=
7
3 Moindres carrés récursifs (5 pts)
estimation.
A l'instant N on peut écrire :
 
a1  

yN
 
−yN −1 −yN −2 . . . −yN −n uN . . . uN −p

 a2  eN
yN −1   −yN −2 −yN −3 −y u uN −p−1   
. 
eN −1 
N −n−1 N −1

 . 
 
 .   eN −2
  
yN −2 . .
 
. .
  
  −yN −3 . . uN −2 uN −p−2

.

  an 
 
yN −3 .
 
= + . (3)
  
. .
  
. .   b0 
 
. . −yn+1 . uN −p−3
 
. .
  
. .
  
  b1   .
   
.
 
.
  
yn+2   −yn+1 −yn . uN −p  .   
en+2

 .. 
 
yn+1 −yn −yn−1 . . . −y1 un+1 en+1
bm
sous forme condensée :

YN = XN θN + EN
dont l'estimation optimale est :
T
θbN = (XN XN )−1 XN
T
YN
1.
    a1
 
eN +1

yN +1 −yN −yN −1 . . . −yN −n+1 uN +1 . . . uN −p+1
 yN   −yN −1 −yN −2 . . . −yN −n uN . . . uN −p  a2   eN 
  
. 

eN −1
  
yN −1 −yN −2 −yN −3 −yN −n−1 uN −1 uN −p−1  . 
 
.
  
  
eN −2
  
yN −2 . .    
. . an  
  
  −yN −3 . . uN −2 uN −p−2  
.

= +

yN −3  .

b0  .

. .   
. .
  
 . . −yn+1 . uN −p−3    
 .. b1   .
 
 .
 
.
 
.   
. . 
  
 yn+2   −yn+1 −yn . uN −p  . 
 
 .  en+2 
yn+1 −yn −yn−1 ... −y1 un+1 bm en+1
(4)
2. xN +1 = −yN − yN −1 . . . − yN −n+1 uN +1 . . . uN −p+1

86 ANNEXE .ANNALES
3.
T −1 T
θbN +1 = (XN +1 XN +1 ) XN +1 YN +1
T
= (XN XN + xTN +1 xN +1 )−1 (XN
T
YN + xTN +1 yN +1 )
4. Il sut de développer
−1
T
XN )−1 − (XN
T
XN )−1 xTN +1 1 + xN +1 (XN
T
XN )−1 xTN +1 T
XN )−1

T
YN + xTN +1 yN +1

XN
5.
−1
T
XN )−1 − (XN
T
XN )−1 xTN +1 1 + xN +1 (XN
T
XN )−1 xTN +1 T
XN )−1

T
YN + xTN +1 yN +1

XN
T
XN )−1 − (XN
T
XN )−1 xTN +1 a−1 xN +1 (XN
T
XN )−1 XN
T
YN + xTN +1 yN +1

θbN +1 = (XN
T
= (XN XN )−1 XN
T T
YN +(XN XN )−1 xTN +1 yN +1 − (XN
T
XN )−1 xTN +1 a−1 xN +1 (XN
T
XN )−1 XN
T
YN
| {z } | {z }
θbN θbN
−(XN XN )−1 xTN +1 a−1 xN +1 (XN

T T
XN )−1 xTN +1 yN +1
= θbN + (XN T
XN )−1 xTN +1 yN +1 − (XNT
−(XN T
XN )−1 xTN +1 a−1 xN +1 (XNT
XN )−1 xTN +1 yN +1

T
= θbN − (XN XN )−1 xTN +1 a−1 xN +1 θbN
T
XN )−1 xTN +1 a−1 −xN +1 (XN T
XN )−1 xTN +1 yN +1 + (XN T
XN )−1 xTN +1 yN +1

+(XN

= θbN − (XNT
T
XN )−1 xTN +1 −a−1 xN +1 (XN T
XN )−1 xTN +1 yN +1 + yN +1

+(XN

= θbN − (XNT
T
XN )−1 xTN +1 a−1 −xN +1 (XN T
XN )−1 xTN +1 + a yN +1

+(XN

= θbN − (XNT
+(XN T
XN )−1 xTN +1 a−1 yN +1

= θbN + (XNT
XN )−1 xTN +1 a−1 yN +1 − xN +1 θbN

6. Le terme yN +1 − xN +1 θbN représente l'erreur d'estimation à l'aide des paramètres précédents.
L'equation est alors : la nouvelle estimation est l'ancienne estimation corrigée par un terme
proportionnel à l'erreur d'estimation précédente.
4 Modélisation paramétrique (5 pts)
1. Par un bouclage simple, suivi d'une identication.
2. Les paramètres dépendent nettement du choix fait pour le bouclage.

Par exemple, bouclage avec un gain A=0,25 =⇒ deux pôles simples en -0,5 =⇒ τ ' 2s
Pour une bonne identication du gain statique on choisira :
NTe = 3 à 5 τ
JANVIER 2007 87
Pour une bonne excitation sur le spectre entre 0 et la plus haute fréquence de coupure on choisira :
1 1
0, 3 =
Te 2πτ
5τ
Bon échantillonnage =⇒ Te ' 3, 7s Bonne identication du gain statique =⇒ n ' Te = 3.
3. ARMAX donnera les meilleurs résultats.

88 ANNEXE .ANNALES
3ème année
Rattrapage 2012 Option mécatronique

Cet examen consiste à déterminer un modèle de la peau humaine soumise à un échelon de compression.
La peau est un matériaux complexe (viscohyperélastique) très variable d'un individu à l'autre et même
d'une région à l'autre sur un même individu. L'objectif et de déterminer un modèle de peau "local".
Les valeurs numériques des diérents modèles sont volontairement diérentes. Les questions suivantes
sont liées mais peuvent être traitées dans le désordre.
Figure 19 Système étudié.
1 Modèle de connaissances (5 pts)
Un modèle simplié de la peau est donné par la gure 20
Figure 20 Modèle de peau viscoélastique.

σ = contrainte
dε(t) ε = déformation = ε = ∆x/x
σ(t) = Eε(t) + η avec :
dt E = module d'Young équivalent
η = coecient de frottement visqueux
le modèle de la force exercée par le poinçon sur la peau est donné par :
F = force
F (p) Kf
= avec : i = courant dans la bobine
i(p) 1 + τp
τ = L/R = constante de temps =inductance/résistance
Applications numériques : E =300 Mpa, η/E = 100 ms, Kf =10, τ =1ms, épaisseur totale de la peau :
1mm, diamètre du poinçon : Φ=1mm.
x(p)
1. Déterminer la fonction de transfer du système H(p) = i(p) .
2. Déterminer les termes T ,D ,V et Q utilisés dans le formalisme de Lagrange.

d ∂T ∂T ∂D ∂V
− + + = Qi
dt ˙
∂ qi ˙
∂qi ∂ qi ∂qi
2. MODÈLE NON-PARAMÉTRIQUE 89
Figure 21 Réponse indicielle de la peau.

1. Identier le système avec la méthode de Broïda.
2. Identier le système avec la méthode de Strejc.
3. Quel modèle vous semble le plus pertinent ?
3 Signal d'excitation
Notre système étant très sensible à la valeur moyenne du signal d'excitation, on décide d'utiliser une
séquence ternaire pseudo aléatoire (STPA). Une méthode de synthèse de séquences ternaires pseudo
aléatoires consiste à les déduire d'une SBPA par l'algorithme suivant :
1
y(i) = (x(i) − x(i − 1))
2
z = {y(0), y(2), y(4), · · · , y(1), y(3), · · ·}
avec : x = séquence binaire pseudo aléatoire (valeurs : -1 ou 1)

y = séquence intermédiaire
z = séquence ternaire pseudo aléatoire
1. Représentez le diagramme temporel d'une séquence binaire pseudo aléatoire utilisant n=3 bits
L = 2n − 1
L
soit de longueur en prenant b0 = b2 b1 .
2. En prenant une période d'échantillonnage Te =100 ms, quelles sont les constantes de temps τmin
et τmax que l'on peux correctement identier en utilisant la SBPA précédente ?
3. Avec la méthode proposée, déduisez-en une séquence ternaire pseudo aléatoire.
4. Quelle est la constante de temps τmax que l'on est susceptible de bien identier ?
5. Calculer la valeur moyenne de la STPA calculée précédemment.
6. Calculez sa fonction d'autocorrélation ( ! long !).

90 ANNEXE .ANNALES
4 Modèle paramétrique
On décide de procéder à une identication paramétrique récursive. Le modèle contient 5 parametres :

a1 , a2 , a3 , b0 et b1 .
b0 + b1 z −1
H(z) =
1 + a1 z −1 + a2 z −2 + a3 z −3
1. Donner la taille des vecteurs et matrices suivants (notation du cours) : θbN , KN , PN . Attention
3
il y à des erreurs dans le poly !
2. Le terme b0 tends vers 0, que peut-on en déduire concernant le retard pur du système ?
5 Modèle d'ordre non entier
En fait il s'agit d'un modèle incluant une partie approchant un modèle d'ordre fractionnaire
(p
α, α ∈ R).
le modèle est de la forme :
n
Y z − zz i
H(z) = K
z − zp i
i=0
avec :
zp i+1 = zp i ∗ λ
zz i = z p i ∗ ν
zp 0 = z 0
ν < λ
Figure 22 Caractéristiques d'un modèle d'odre fractionnaire : réponse indicielle et pôles et zéros.
1. Proposez une méthode pour déterminer les paramètres λ , ν , z0 et K optimums vis-à-vis d'un
critère que vous dénirez.
3. désolé !
1. MODÉLISATION DE LA GÉNÉRATRICE 91
Figure 23 Réponse indicielle de la peau.
3ème année
Janvier 2009

Les parties suivantes sont totalement indépendantes et au sein des parties certaines questions sont
indépendantes.
Ce sujet d'examen traite de l'identication et de la commande d'une éolienne utilisant l'eet Magnus.
Sur cette éolienne, les pales sont en fait des cylindres qui tournent sur leur axe de symétrie. La force
engendrée par cette rotation combinée à la vitesse du vent est illustrée en gure 24. L'avantage de ce
système est une plus grande plage d'utilisation en fonction de la vitesse du vent, la vitesse de rotation
de l'éolienne étant contrôlée directement par la vitesse de rotation des pales.
1 Modélisation de la génératrice
L'éolienne se compose d'une génératrice multipolaire de vitesse de rotation nominale de 1500 tr/mn.
La génératrice est modélisée comme indique sur la gure . On procède à deux essais de lâcher, c'est à
dire que l'on amène la génératrice à sa vitesse nominale puis on coupe l'alimentation.
un essai en induit ouvert (u = E, la force électromotrice)
un essai induit en court-circuit (u = 0)
92 ANNEXE .ANNALES
Figure 24 Eolienne utilisant l'eet Magnus.

R L
i
-
6 6'$
J, f
u E
&%
Figure 25 Schéma de la génératrice
1.1 Essai de lâcher : induit ouvert
Figure 26 Essai de lâcher en circuit ouvert.

2. IDENTIFICATION DE LA FORCE DES PALES 93
1. A partir du schéma de la gure 25, déterminer le modèle analytique de la génératrice en tenant

compte des conditions initiales : ω = ω0 et i = i0 .
2. la génératrice à été amenée à sa vitesse nominale (régime permanent) ω = ω0 = 1500tr.mn−1
en appliquant un échelon de tension U = 120V Par une méthode de votre choix, déterminer un
modèle dynamique de ce système. On justiera la méthode utilisée !
3. Montrer que pour un système du premier ordre de constante de temps τ la vitesse pendant un
essai de lâcher ω(t) = L−1 [ω(p)] suit la loi :
ω(t + T ) = e−T /τ ω(t) ∀t > 0 et T >0
4. Tracer pour T =100s et t = {0, 100, 200, 300} les points d'abscisse ω(t) et d'ordonnée ω(t + T ).
4
5. Déterminer graphiquement la droite qui passe au mieux par les points .
6. Outre le coecient de corrélation du nuage de points, que représente la pente de la droite ?
7. Déduisez-en la valeur de τ.
8. Comparer cette dernière valeur de τ à la constante de temps obtenue à la deuxième question.
1.2 Essai de lâcher : induit en court-circuit
En faisant l'hypothèse que la génératrice peut être modélisée par une fonction de transfert de la forme :
Go ω(p)
H(p) = =
(1 + τ1 p)(1 + τ2 p) u(p)
1. Déterminer à partir de la gure 27 les valeurs de τ1 , τ2 et G0 par une méthode de votre choix
mais justiée.
2. En notant a1 = e−T /τ1 et a2 = e−T /τ2 , montrer que l'équation diérentielle
ω(t + 2T ) − (a1 + a2 )ω(t + T ) + a1 a2 ω(t) = 0
ω(t + 2T ) ω(t + T )
− (a1 + a2 ) + a1 a2 = 0
ω(t) ω(t)
ω(t+T )
3. En prenant les points sur la gure 28, tracer dans une repère, les points d'abscisse et
ω(t)
ω(t+2T )
d'ordonnée
ω(t)
4. Déterminer les relations entre la pente de la droite qui passe au mieux par les points, l'ordonnée
à l'origine et les valeurs a1 et a2 .
5. Déterminer les constantes de temps τ1 et τ2 .
2 Identication de la force des pales
La gure 29 représente le diagramme de Bode obtenu en faisant varier la vitesse de rotation des pales
sur elles-mêmes autour de 1000 tr.mn
−1 avec un signal sinusoïdal d'amplitude 100 tr.mn−1 (soit entre
900 et 1100 tr.mn

−1 ).
1. Déterminer les caractéristiques (nombre de bits, période d'échantillonnage, amplitude) d'une

séquence binaire pseudo aléatoire convenable pour espérer avoir une bonne identication de ce
système.
2. Combien de temps prendra l'identication du système.
3. Quel modèle ARX tenterez vous en premier ?
4. Si c'est le fruit d'un calcul propre, cette droite est la droite de corrélation.
94 ANNEXE .ANNALES
Figure 27 Essai de lâcher en cout-circuit.
Figure 28 Essai de lâcher en cout-circuit (zoom).

3 Contrôle de la vitesse de l'éolienne
Le contôle de la vitesse d'un éolienne est très important à deux titres :

4. RIEN À VOIR 95
Figure 29 Diagramme de Bode : couple en fonction de la vitesse de rotation des pales.
Le rendement de la génératrice chute rapidement en cas de survitesse

La sécurité n'est plus assurée en cas de survitesse.
Dans le cas d'un éolienne à eet Magnus, il y a deux variables qui permettent de contrôler la vitesse
de l'éolienne.
Le courant consommé (donc le couple résistant)

La vitesse de rotation des pâles, donc le couple exercé par les pales sur la génératrice.
Les ensembles ous des variables d'entrée V : vitesse du vent et ω vitesse de rotation de l'éolienne
sont dénis en gure 30.
1. Pour le réglage de la puissance, proposer un ensemble minimal de règles, une fuzzication des
variables de sortie et une méthode de défuzzication qui respecte le cahier des charges suivant :
la vitesse de rotation des pâles doit être nulle lorsque la vitesse du vent est nulle (nul) ou trop
forte (arret)
la vitesse de rotation des pâles est comprise entre 0 et 1000 tr/mn.
le courant est compris entre 0 et 10 A.
le courant doit être maximal.
4 Rien à voir
1. Déterminer la parabole y = ax2 + bx + c qui passe "au plus près" au sens du moindre carré des
erreurs par les point :
x −3 −1 0 2 3 6
y 17 3 1 5 11 46
2. Quelle est la valeur de cette parabole en x = 2.

96 ANNEXE .ANNALES
Figure 30 Ensembles ous de la variable d'entrée vitesse du vent.
Liens
http://www.mathworks.com/products/fuzzylogic.html
http://www.mathtools.net/MATLAB/Fuzzy_Logic/index.html
http://www.dbai.tuwien.ac.at/marchives/fuzzy-mail/
http://europa.eu.int/comm/research/rtdinfo/fr/24/02.html
http://www.gala.univ-perp.fr/~polit/chap0.html
http://www.eru.ulaval.ca/ptt15225
http://vcampus.u-strasbg.fr/uticeweb/mhiri/projet_V2
http://www.mathworks.fr/products/controldesign/modanal.shtml
Bibliographie
[1] C. Foulard, S.Gentil, and J.P. Sandraz. Commande et régulation par calculateur numérique : de la
théorie aux applications. Eyrolles, Paris, 1987.
[2] Lennart Ljung. System identication (2nd ed.) : theory for the user. Prentice Hall PTR, Upper
Saddle River, NJ, USA, 1999.
[3] Lennart Ljung and Torkel Glad. Modeling of dynamic systems. Prentice-Hall, Inc., Upper Saddle
River, NJ, USA, 1994.
[4] Lennart Ljungs. System identication toolbox for use with matlab, 1997.
[5] William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical
Recipes : The Art of Scientic Computing. Cambridge University Press, Cambridge (UK) and New
York, 2nd edition, 1992. Disponible gratuitement sur Internet.
97

Identification Des Systemes PDF

Transféré par

Droits d'auteur :

Formats disponibles

Identification Des Systemes PDF

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Identification Des Systemes PDF

Transféré par

Droits d'auteur :

Formats disponibles

Identification des systèmes

École Nationale Supérieure de

" Je distingue deux moyens de cultiver les sciences : l'un

2 Identication de modèles non paramétriques 15

3 Algorithme général d'identication 21

5 Algorithmes d'optimisation paramétrique. 39

6 Description des algorithmes utilisés. 45

Examen nal janvier 2007 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Identication des systèmes

1.1 Modèles simples linéaires

1.1.1 Exemple : le moteur à courant continu

Figure 1.1  Schéma électromécanique d'un moteur à courant continu

De l'ensemble des équations on peut tirer diérents graphes :

Figure 1.3  Analogie des équations diérentielles linéaires.

1.2 Méthodes systématiques

avec : T = Σ énergies cinétiques

1.2.1 Un exemple simple

Figure 1.4  Sustentation magnétique.

d'où on en déduit les formules suivantes :

1.3 Détermination des constantes

q̈1 = F1 (q1 , q2 , q̇1 , q̇2 , e)

La linéarisation autour du point E est donnée par :

Identication de modèles non

2.1 Quelle entrée ?

2.1.1 Analyse harmonique

2.1.2 Réponse impulsionnelle

2.1.3 Réponse indicielle

2.2 Etude des réponses indicielles

2.2.1 Systèmes du premier ordre

Table 2.1  Réponse d'un système du premier ordre à un échelon.

2.2.2 Systèmes du second ordre résonnant

Table 2.2  Réponse d'un système du deuxième ordre à un échelon.

Le facteur d'amortissement m se détermine à l'aide de la mesure du dépassement :

La pulsation propre w0 s'obtient via l'une des deux formules :

2.2.3 Systèmes du premier ordre retardés ou du second ordre non résonnant

Cette approximation reste valable jusqu'à la pulsation ω = 1/T

2.2.4 Systèmes d'ordre supérieur à 2 non résonnants

La méthode la plus connue est la méthode de Strejc. Le modèle est :

Table 2.4  Méthode de Strejc : Réponse indicielle.

n Tu/Ta Tu/τ Ta/τ

0.02 9 0.709 5.081 7.164

qui s'annule pour t = 3τ .

2.2.5 Système avec intégrateur

Table 2.5  Méthode de Strejc (système avec intégrateur)

Identié sous la forme

1. Avec x/y : déterminer n

2.3 Méthode de Ziegler-Nichols

2.4 Autres modèles ?

Méthode de pompage Méthodes apériodiques

Systèmes stables ou instables Systèmes stables, instables

Ziegler-Nichols Ziegler-Nichols Chien-Hrones-Reswick

Algorithme général d'identication

 Quel genre de modèle ?

3.1 Choix du signal d'excitation

3.2 Les séquences binaires pseudo aléatoires (SBPA)

Figure 3.3  diagramme temporel d'une séquence binaire pseudo aléatoire.

Figure 3.4  Fonction d'auto-corrélation d'une SBPA.

Figure 3.5  Spectre d'une SBPA.

3.2.1 Choix des paramètres d'une SBPA

2 Identication de modèles non paramétriques 15

3 Algorithme général d'identication 21

Examen nal janvier 2007 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Identication des systèmes

Figure 1.1 Schéma électromécanique d'un moteur à courant continu

De l'ensemble des équations on peut tirer diérents graphes :

Figure 1.3 Analogie des équations diérentielles linéaires.

Figure 1.4 Sustentation magnétique.

Identication de modèles non

Table 2.1 Réponse d'un système du premier ordre à un échelon.

Table 2.2 Réponse d'un système du deuxième ordre à un échelon.

Table 2.4 Méthode de Strejc : Réponse indicielle.

Table 2.5 Méthode de Strejc (système avec intégrateur)

Identié sous la forme

Algorithme général d'identication

Quel genre de modèle ?

Figure 3.3 diagramme temporel d'une séquence binaire pseudo aléatoire.

Figure 3.4 Fonction d'auto-corrélation d'une SBPA.

Figure 3.5 Spectre d'une SBPA.

Il faut donc déterminer a, L etT e "optimum" vis-à-vis du système à identier.

Figure 3.6 Paramètres d'une SBPA.

3.3 Identication basée sur l'erreur de sortie

Figure 3.7 Principe d'identication fondé sur l'erreur de sortie.

adapté à la recherche de paramètres physiques si modèle de connaissances (modèle continu !)

3.4 Identication basée sur l'erreur de prédiction

Figure 3.8 Principe d'identication foncée sur l'erreur de prédiction.

Figure 3.9 Modèle du système étudié.

matrice dénie positive : c'est bien un minimum !

Figure 3.10 Système du premier ordre bruité.

Figure 3.11 Fonction d'auto-corrélation des résidus d'estimation.

Le système étudié est décrit sur la gure 3.12

Figure 3.12 Système étudié.

3.5 Identication en boucle fermée

Figure 3.13 Identication des systèmes en boucle fermée.

Table 3.2 Modèles paramétriques implantés sous Matlab.

On obtient nalement une formule de récurrence :

Figure 5.1 Classes d'algorithmes d'optimisation paramétrique

5.5 Eléments de comparaison des diérentes méthodes.