Polyetat
Polyetat
Luc Jaulin
5 décembre 2008
2 Commande par espace d’état
Table des matières
1 Introduction 13
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.1 Intégrateur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 Méthodologie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3
4 Commande par espace d’état
3 Simulation 41
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4.1 Motif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.5 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4 Systèmes linéaires 61
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3.1 Dé…nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.4.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.5.6 Composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.5.7 Exemple 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.5.8 Exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.2.2 Commandabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.2.3 Observabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.4.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
8 Identi…cation 137
Bibliographie 160
Table des matières 11
Glossaire 163
Index 165
12 Commande par espace d’état
Chapitre 1
Introduction
Une très grande classe de systèmes physiques, biologiques, économiques, etc. qui nous entourent peuvent
être décrits par une équation di¤érentielle du type
(
x(t)
_ = f (x(t); u(t))
(1.1)
y(t) = g(x(t); u(t));
sous l’hypothèse que le temps t dans lequel évolue le système est continu (c’est-à-dire appartient à R). Le
vecteur u(t) est l’entrée (ou commande) du système. Sa valeur peut être choisie arbitrairement pour tout
t. Le vecteur y(t) est la sortie du système et peut être mesuré avec une certaine précision. Le vecteur x(t)
est appelé état du système. Comme nous l’illustrerons dans la suite, il représente la mémoire du système,
c’est-à-dire, l’ensemble des informations dont le système a besoin pour prédire son propre avenir, pour une
entrée u(t) connue. La première des deux équations de (1.1) s’appelle équation d’évolution. Il s’agit d’une
équation di¤érentielle qui permet de savoir vers où va se diriger l’état x(t) sachant sa valeur à l’instant
présent t et la commande u(t) que nous sommes en train d’exercer. La deuxième équation s’appelle équation
d’observation. Elle permet de calculer le vecteur de sortie y(t), connaissant l’état et la commande à l’instant
t. Attention, contrairement à l’équation d’évolution, cette équation n’est pas une équation di¤érentielle car
elle ne fait pas intervenir les dérivées des signaux. Les équations (1.1) forment la représentation d’état du
système.
Il est parfois utile de considérer aussi un temps k discret, où k sera un élément de Z, l’ensemble des entiers
relatifs. En e¤et, si l’univers que nous considérons est un ordinateur, il est concevable de considérer que le
temps k qui le régit est discret, synchronisé sur l’horloge du microprocesseur. De tels systèmes, dits à temps
discret, obéissent généralement à une équation de récurrence du type
(
x(k + 1) = f (x(k); u(k))
(1.2)
y(k) = g(x(k); u(k)):
Dans ce cours (qui est fortement inspiré de mon livre [8], paru en 2005), sauf mention particulière, nous
allons uniquement nous intéresser aux systèmes à temps continu (c’est-à-dire de type (1.1)).
Le premier objectif de ce cours est de bien comprendre cette représentation d’état. Pour cela, nous allons
considérer dans le chapitre 2 un grand nombre d’exemples variés et montrer comment naturellement nous
13
14 Commande par espace d’état
arrivons à une représentation d’état pour modéliser le système qui nous intéresse. Ensuite, nous montrerons
au chapitre 3 comment simuler sur ordinateur un système pour lequel une représentation d’état est connue.
Le second objectif de ce cours est de proposer des méthodes pour commander les systèmes décrits par des
équations d’état. C’est-à-dire que nous allons tenter de fabriquer des machines automatiques (où l’homme
n’intervient quasiment pas, sauf pour donner ses ordres, où consignes), appelées régulateurs capables de
domestiquer (changer le comportement dans le sens que l’on souhaite) les systèmes considérés. Pour cela, le
régulateur devra calculer les entrées u(t) à appliquer au système à partir de la connaissance (plus ou moins
bruitées) des sorties y(t) et de la consigne w(t) donnée (voir …gure).
Vu de l’utilisateur, le système, dit système bouclé, d’entrée w(t) et de sortie y(t) aura un comportement
convenable. On dira que l’on a asservi le système.
Dans cet objectif de régulation, nous allons, dans une première phase, nous intéresser uniquement aux
systèmes linéaires, c’est-à-dire que les fonctions f et g sont supposées linéaires. Ainsi, dans le cas du temps
continu, (1.1) s’écrit (
x(t)
_ = Ax(t) + Bu(t);
(1.3)
y(t) = Cx(t) + Du(t);
et dans le cas du temps discret, (1.2) devient
(
x(k + 1) = Ax(k) + Bu(k);
(1.4)
y(k) = Cx(k) + Du(k):
Les matrices A; B; C; D sont appelées matrices d’évolution, de commande, d’observation et directe. Une
analyse détaillée de ces systèmes sera faite au chapitre 4. Nous expliquerons ensuite, dans le chapitre 5
comment stabiliser ces systèmes.
Dans une seconde phase, nous montrerons, au chapitre 6, qu’autour de certains points, dits de fonction-
nement, les systèmes non linéaires se comportent comme des systèmes linéaires et qu’il alors possible de
les stabiliser par des méthodes développées pour le cas linéaire. Nous donnerons aussi d’autres types d’ap-
proches, moins générales mais plus adaptées, pour stabiliser des systèmes non linéaires.
Les programmes informatiques utilisés pour tester nos régulateurs et e¤ectuer nos simulations seront tous
donnés pour un environnement Scilab[2]. Ces derniers sont disponibles sur
http ://[Link]/jaulin/[Link].
Rappelons que Scilab est un logiciel de calcul numérique qui intègre les opérations courantes du calcul
matriciel (tout comme Matlab), et o¤re l’avantage d’être gratuit (contrairement à Matlab). Il peut être
Introduction 15
téléchargé à l’adresse
http ://[Link] /scilab/.
Le calcul formel sera quant à lui e¤ectué sous Mupad[1] qui est un logiciel proche de Maple. Notons que
contrairement à Maple, Mupad est gratuit et peut être téléchargé sur
http ://[Link].
Dans le …chier compressé [Link], vous trouverez les programmes Scilab suivants :
[Link], [Link],
BateauHybride, [Link],
BrasManipulateur, CharLibre,
[Link], [Link],
[Link], [Link],
[Link], SegWayAsservi,
SegwayLibre, VerinLibre,
VoitureAsservieCircuit, VoitureAsservieCycle,
[Link], VoitureAsservieNord,
[Link], [Link],
[Link], [Link],
[Link].
Dans les noms des …chiers, on retrouve le nom du système traité (voiture, pendule, etc), un quali…catif (libre
signi…e non asservi) et l’extension sce rattachée aux programmes Scilab.
16 Commande par espace d’état
Chapitre 2
2.1 Introduction
Rappelons que de nombreux systèmes peuvent se mettre sous une représentation d’état,
(
x(t)
_ = f (x(t); u(t))
y(t) = g(x(t); u(t));
pour les systèmes à temps continus et
(
x(k + 1) = f (x(k); u(k))
y(k) = g(x(k); u(k));
pour les systèmes à temps discret. Nous appellerons modélisation l’étape qui consiste à trouver une repré-
sentation d’état, plus ou moins précise, du système qui nous intéresse. En général, apparaissent dans les
équations d’états des paramètres constants (comme la masse ou le moment d’inertie d’un corps, le coe¢ -
cient du frottement visqueux, la capacité d’un condensateur, etc). Dans de tels cas, une étape d’identi…cation
peut s’avérer nécessaire. Dans ce cours, nous supposerons que tous les paramètres sont connus et dans le cas
contraire, nous invitons le lecteur à consulter le livre d’Eric Walter et Luc Pronzato [7] pour une panoplie
des méthodes d’identi…cation.
Bien sûr, il n’existe pas de méthodologie systématique à appliquer pour modéliser un système. Le but de ce
chapitre est de présenter sur quelques exemples variés comment arriver à une représentation d’état a…n de
permettre au lecteur d’acquérir une certaine expérience qui l’aidera à modéliser ses propres systèmes.
Dans ce paragraphe, nous présentons quelques exemples de systèmes linéaires à temps continus, qui, rappelons-
le, peuvent être décrits par les équations d’état
(
x(t)
_ = Ax(t) + Bu(t);
y(t) = Cx(t) + Du(t):
17
18 Commande par espace d’état
2.2.1 Intégrateur
L’intégrateur est un système linéaire décrit par l’équation di¤érentielle y_ = u. Une représentation d’état
possible pour ce système est (
x(t)
_ = u(t)
y(t) = x(t):
Les matrices associées à ce système sont A = (0), B = (1), C = (1) et D = (0) : Elles sont toutes de
dimension 1 1.
y• + a1 y_ + a0 y = bu:
ou encore, sous une forme standard faisant apparaître les matrices d’état A, B, C et D :
8 ! !
>
> 0 1 0
< x_ = x+ u
a0 a1 b
>
>
: y = 1 0 x:
2.3.1 Méthodologie
Le principe fondamental de la dynamique permet de trouver simplement les équations d’état de systèmes
mécaniques (ou robots). Les calculs engendrés sont relativement lourds pour les systèmes complexes et
l’utilisation du logiciel de calcul formel (Mupad par exemple) peut s’avérer utile. La méthode sera illustrée
à travers plusieurs exemples.
Pour obtenir les équations d’état d’un système mécanique composé de plusieurs sous-systèmes S1 ; S2 ; : : : ; Sm ;
supposés rigides, on procède en trois étapes.
1. Obtention des équations di¤ érentielles. Il faut appliquer pour chaque sous-système Sk , de masse m et
de matrice d’inertie J, les relations
P
f = ma
P i i
i Mfi = J! _
Modélisation des systèmes 19
où les fi sont les forces agissant sur le sous-système Sk , Mfi représente le moment engendré par la force
fi sur Sk , par rapport à son centre de gravité. Le vecteur a représente l’accélération tangentielle de Sk
et le vecteur ! _ représente l’accélération angulaire de Sk . Après décomposition de ces 2 m équations
vectorielles suivant leurs composantes, on obtient 6 m équations di¤érentielles scalaires dont certaines
peuvent être dégénérées (c’est-à-dire du type 0 = 0).
2. Suppression des composantes des forces intérieures. Dans les équations di¤érentielles interviennent
des forces, dites de liaisons, qui sont intérieures au système mécanique global, bien qu’extérieures à
chaque sous-système qui le compose. Elles représentent l’action d’un sous-système Sk sur un autre
sous-système S` . D’après le principe d’action-réaction, l’existence d’une telle force, notons la f k;` ,
implique l’existence d’une autre force f `;k , représentant l’action de S` sur Sk ; telle que f `;k = f k;` .
Par une manipulation formelle des équations di¤érentielles et en prenant en compte les équations dues
au principe d’action-réaction, il est possible de se débarrasser des composantes fx ; fy et fz de ces forces
intérieures. Le nombre des équations di¤érentielles résultantes doit se réduire au nombre n de degrés
de liberté q1 ; : : : ; qn du système.
3. Obtention des équations d’état. Il nous faut alors isoler les dérivées secondes q•1 ; : : : ; q•n de l’ensemble
des n équations di¤érentielles de façon à obtenir une relation vectorielle du type
• = f (q; q;
q _ u) ;
où u est le vecteur des forces extérieures et ne dérivant pas d’un potentiel (c’est-à-dire, celles que nous
appliquons au système). Les équations d’état s’écrivent alors
! !
d q q_
= :
dt q_ f (q; q;
_ u)
Le lecteur intéressé par plus de détails sur la modélisation des systèmes mécaniques pourra consulter
[5] et [13]. Nous allons maintenant illustrer la méthodologie présentée sur divers exemples.
Considérons le satellite, représenté sur la …gure ci dessous, et repéré par ses coordonnées (x1 ; x2 ). Ce satellite
est en orbite autour de la terre.
20 Commande par espace d’état
La loi de gravitation universelle de Newton a¢ rme que la terre exerce une force d’attraction f sur le satellite
donnée par
!
x1
f = G xM m
2 +x2 u = G xM m 1
2 +x2 : p 2
1 2 1 2 x1 +x22 x2
!
x1
= G p M2 m 2 3 ;
x1 +x2 x2
Puisque
(
x_ 1 = x3
x_ 2 = x4 ;
On considère le système d’entrée f et de sortie z1 de la …gure ci-dessus (f : force appliquée sur le deuxième
chariot, zi écart du iième chariot par rapport à sa position d’équilibre, ki raideur du iième ressort,
coe¢ cient de frottement visqueux). Prenons pour vecteur d’état
(
f z_2 k2 (z2 z1 ) = m2 z•2 ;
k1 z1 z_1 + k2 (z2 z1 ) = m1 z•1 :
8 0 1 0 1
>
> 0 1 0 0 0
>
> B k1 +k2 k2 C B C
>
> B 0 C B 0 C
< x_ = B m1 m1 m1
C x +B Cf
@ 0 0 0 1 A @ 0 A
>
> k2 k2 1
>
> m2 0 m2 m2 m2
>
>
: z1 = 1 0 0 0 x
On considère le pendule de la …gure ci dessous. L’entrée de ce système est le couple c exercé sur le pendule.
La sortie est y(t); la distance entre la masse m et l’axe vertical. Cherchons les équations d’état de ce système.
22 Commande par espace d’état
Pendule simple
Comme indiqué au point 3 du paragraphe 2.3.1 de la page 18, nous devons prendre pour vecteur d’état
x = (q; q).
_ Les équations d’état du système s’écrivent alors :
! !
d q q_
dt = `mg sin q+c ;
q_ m`2
(2.1)
y = ` sin q:
1 2 2
Em = m` q_ + mg` (1 cos q):
|2 {z } | {z }
énergie potentielle
énergie cinétique
dEm 1 2
dt = 2 m` (2q_q •) + mg`q_ sin q
`mg sin q
= m`2 q_ m`2 + mg`q_ sin q = 0:
L’énergie mécanique du pendule reste donc bien constante, ce qui est cohérent avec le fait que le pendule
sans frottement est un système conservatif.
Modélisation des systèmes 23
Pendule inversé
On considère le système, appelé pendule inversé, formé d’un pendule posé en équilibre instable sur un chariot
roulant, comme représenté sur la …gure. La quantité u est la force exercée sur le chariot de masse M , x
indique la position du chariot, est l’angle entre le pendule et la verticale et R ~ est la force exercée par le
chariot sur le pendule. A l’extrémité B du pendule est …xée une masse m. On négligera la masse de la tige
du pendule. En…n, A est le point d’articulation entre la tige et le chariot et ~ = _ ~k est le vecteur de rotation
associé à la tige.
Pour obtenir les équations d’état de ce système, nous allons à nouveau suivre la procédure présentée au
paragraphe 2.3.1 de la page 18. Le principe fondamental de la dynamique appliqué sur le chariot et le
pendule nous donne
où vB est le vecteur vitesse du point B. Pour la troisième équation, le moment d’inertie du pendule a été
posé nul. Puisque
!
OB = (x ` sin )~i + ` cos ~j;
nous avons
vB = x_ ` _ cos ~i ` _ sin ~j:
Donc, l’accélération du point B est donnée par
2 2
v_ B = x
• `• cos + ` _ sin ~i `• sin + ` _ cos ~j:
Ces quatre équations décrivent respectivement (i) le chariot en translation, (ii) le pendule en translation
suivant ~i; (iii) le pendule en translation suivant ~j et (iv) le pendule en rotation. On véri…e bien que le nombre
de degré de liberté (ici x et ) additionné au nombre de composantes des forces intérieures (ici Rx et Ry )
est égal au nombre d’équations. Sous forme matricielle, ces équations s’écrivent
0 10 1 0 1
M 0 1 0 x• u
B m m` cos 1 0 CB • C B 2
m` _ sin C
B CB C B C
B CB C=B C:
@ 0 m` sin 0 1 A @ Rx A @ mg m` _ 2 cos A
0 0 cos sin Ry 0
Donc
0 1 1 0 1
! ! M 0 1 0 u
x
• 1 0 0 0 B m m` cos 1 0 C B 2
m` _ sin C
B C B C
= :B C :B C:
• 0 1 0 0 @ 0 m` sin 0 1 A @ mg m` _ 2 cos A
0 0 cos sin 0
B :=matrix([u,m*l*w^2*sin(a),m*g-m*l*w^2*cos(a),0])
P :=matrix([[1,0,0,0],[0,1,0,0]])
A :=matrix([[M,0,1,0],[-m,m*l*cos(a),1,0],
[0,m*l*sin(a),0,1],[0,0,cos(a),sin(a)]])
normal(P/A*B)
1
L’instruction normal(P/A*B) renvoie la matrice PA B avec des coe¢ cients réduits au même dénominateur.
Mupad nous donne
! 0 2 1 !
m sin (` _ g cos ) `
x
• M +m sin2 `(M +m sin2 )
• =@ 2 A+
cos u:
sin ((M +m)g m` _ cos )
2
`(M +m sin ) `(M +m sin2 )
0 1 0 1 0 1
x x_ 0
B C B _ C B C
B C B C B 0 C
d
dt B C = B
B
2
m sin (` _ g cos ) C+B
C B ` C u:
C (2.3)
@ x_ A @ M +m sin2 A @ `(M +m sin2 ) A
2
_ sin ((M +m)g m` _ cos ) cos
`(M +m sin2 ) `(M +m sin2 )
Modélisation des systèmes 25
Le Segway représenté sur la …gure est un véhicule à deux roues pour un seul essieu. Il est stable car il est
régulé. Dans une phase de modélisation, nous allons bien sûr supposer que l’engin n’est pas régulé. Son
comportement en boucle ouverte est très proche du monocycle planaire représenté sur la …gure ci-dessous.
Sur cette …gure, u représente le couple exercé entre le corps et la roue.
La liaison entre ces deux éléments est une liaison pivot. Nous noterons B le centre de gravité du corps et
!
A celui de la roue. C est un point …xé sur le disque. Notons l’angle entre le vecteur AC et l’horizontale
et l’angle entre le corps du monocycle et la verticale. Ce système possède deux degrés de liberté et .
26 Commande par espace d’état
D’après les explications données au paragraphe 2.3.1 de la page 18, l’état de notre monocycle est donné par
le vecteur x = ; ; _ ; _ . Les paramètres de notre système sont
A…n de trouver les équations d’état, nous appliquons le principe fondamental de la dynamique sur chaque
sous-système, à savoir la roue et le corps. Nous avons :
8
>
> R x + Fx = M a• (roue en translation)
>
< Fx a + u = JM • (roue en rotation)
~ ~ ~ (2.4)
>
> Rx i + Ry j mg j = mv_ B (corps en translation)
>
:
Rx ` cos + Ry ` sin u = Jp • (corps en rotation),
On véri…e bien que le nombre de degré de liberté (ici et ) additionné au nombre de composantes des
forces intérieures (ici Rx , Ry et Fx ) est égal au nombre d’équations. Sous forme matricielle, ces équations
s’écrivent 0 10 1 0 1
Ma 0 1 0 1 • 0
B J a C B C B C
B M 0 0 0 CB • C B u C
B CB C B B _ 2 C
C:
B ma m` cos 1 0 0 C B Rx C = B m` sin C
B CB C B
@ 0 m` sin 0 1 0 A @ Ry A @ mg m` _ 2 cos C A
0 Jp ` cos ` sin 0 Fx u
Donc
0 1 1
0 1
Ma 0 1 0 1 0
! ! B C B C
B JM 0 0 0 a C B u C
x
• 1 0 0 0 0 B C B 2 C
• = :B ma m` cos 1 0 0 C :B
B m` _ sin C:
C
0 1 0 0 0 B C B 2 C
@ 0 m` sin 0 1 0 A @ mg m` _ cos A
0 Jp ` cos ` sin 0 u
Modélisation des systèmes 27
B :=matrix([0,u,m*l*w^2*sin(T),m*g-m*l*w^2*cos(T),u])
P :=matrix([[1,0,0,0],[0,1,0,0]])
A :=matrix([[M*a,0,-1,0,1],[JM,0,0,0,-a],
[m*a,m*l*cos(T),1,0,0],[0,m*l*sin(T),0,1,0],
[0,-Jp,l*cos(T),l*sin(T),0]])
N0 :=normal(P/A*B),
N1 :=subs(N0,sin(T)^2=1-cos(T)^2)
N2 :=normal(N1)
avec
1 = JM + a2 (m + M ) ; 2 = Jp + m`2 ;
3 = aml; g = glm:
Le char (ou le fauteuil roulant) de la …gure comprend deux chenilles (ou roues) parallèles motrices dont les
accélérations (qui forment les entrées u1 et u2 du système) sont commandées par deux moteurs indépendants.
28 Commande par espace d’état
Dans le cas où le fauteuil roulant est considéré, la stabilité du système est assurée par une ou deux roues
folles, non représentées sur la …gure. Les variables d’état sont constituées par les coordonnées x; y du centre
de l’essieu, de l’orientation du char, et des vitesses v1 et v2 du centre de chacune des roues motrices.
Remarquons que le vecteur d’état ne peut être choisi égal à x; y; ; x; _ _ , ce qui pourtant nous aurait
_ y;
semblé naturel au regard de la théorie lagrangienne brièvement rappelée au paragraphe 2.3.1. En, e¤et, si
tel était notre choix, certains états n’auraient aucune signi…cation physique. Par exemple l’état
x = 0; y = 0; = 0; x_ = 1; y_ = 1; _ = 0
n’a pas de sens car le char n’est pas autorisé à déraper. Ce phénomène est dû à l’existence des roues qui
entraîne des contraintes entre les variables d’état naturelles. Ici, nous avons nécessairement la contrainte de
dite de non holonomie :
y_ = x_ tan :
Les systèmes mécaniques pour lesquels il existe de telles contraintes d’égalités sur les variables d’état na-
turelles (par variable d’état naturelles, nous entendons le vecteur (q; q) _ où q est le vecteur des degrés de
liberté de notre système) sont dit non holonomes. Lorsqu’une telle situation arrive, il convient d’utiliser ces
contraintes a…n de réduire le nombre de variables d’état et ceci, jusqu’à ce qu’aucune contrainte de non ho-
lonomie ne subsiste. Ici, nous choisissons pour vecteur d’état (x; y; ; v1 ; v2 ). Ce choix se comprend aisément
dans le sens où ces variables nous permettent de dessiner le char (x; y; ) et la connaissance de v1 ,v2 nous
permet de calculer les variables x; _ _ : De plus, tout choix arbitraire du vecteur (x; y; ; v1 ; v2 ) correspond à
_ y;
une situation physiquement possible. Les équations d’état du système s’obtiennent directement sans aucun
calcul. Elles s’écrivent 0 1 0 v1 +v2 1
x_ 2 cos
B y_ C B v1 +v2 sin C
B C B C
B _ C B 2v2 v1 C
B C=B ` C:
B C B C
@ v_ 1 A @ Ru1 A
v_ 2 Ru2
Le conducteur de la voiture représentée sur la …gure à gauche possède deux commandes : l’accélération
des roues avant (supposées motrices) et la vitesse de rotation du volant. Le frein apparaît ici comme une
Modélisation des systèmes 29
accélération négative. Nous noterons l’angle entre les roues avant et l’axe de la voiture, l’angle que fait
la voiture par rapport à l’horizontale et (x; y) les coordonnées du milieu de l’essieu arrière. En raisonnant
comme dans le paragraphe 2.4.2, il est clair que les variables d’état de notre système sont constituées par
– les coordonnées de position, c’est-à-dire, toute la connaissance nécessaire pour dessiner la voiture, à savoir
les coordonnées x; y du centre de l’essieu arrière, l’orientation de la voiture, et l’angle des roues avant ;
– et la coordonnée cinématique v représentant la vitesse du milieu de l’essieu avant (en e¤et, la seule
connaissance de cette quantité et des coordonnées de position permet de retrouver toutes les vitesses de
tous les autres éléments de la voiture).
Obtention des équations d’état [12]. Supposons pour simpli…er que les deux roues avant tournent à la
même vitesse v. Comme illustré sur la …gure à droite, tout de passe comme s’il existait seulement deux roues
virtuelles situées au centre des essieux. Considérons un observateur …xe par rapport au sol. D’après la règle
de composition des vitesses, nous avons
!
vA = vM + AM ^ ! !;
où !
! est le vecteur de rotation instantanée de la voiture. Notons que cette relation est une relation vectorielle
qui dépend de l’observateur mais qui est indépendante du repère. Exprimons la dans le repère de la voiture,
représenté sur la …gure. Attention, il ne faut pas confondre l’observateur …xe sur le sol et le repère dans
lequel s’exprime la relation. Cette équation s’écrit :
0 1 0 1 0 1 0 1
v cos vM L 0
B C B C B C B C
@ v sin A = @ 0 A + @ 0 A ^ @ 0 A ;
0 0 0 _
Ainsi,
_ (2.6) v sin
= ;
L
et 8
< x_ = v cos (2.6)
= v cos cos
M
: y_ = v sin (2.6)
= v cos sin :
M
Rajout d’une remorque. Rajoutons une remorque à cette voiture dont le point d’attache est le milieu de
l’essieu arrière de la voiture. Les équations régissant la remorque peuvent être obtenues directement à partir
de celles de la voiture.
30 Commande par espace d’état
0 1 0 1
x_ r vr cos r cos r
B C B C
_ =
@ r A @ vr cos r sin
y r A
_r vr sin r
Lr
avec
p q
vr = x_ 2 + y_ 2 = (v cos cos )2 + (v cos sin )2 = v cos ;
r = r:
0 1 0 1
x_ v cos cos
B C B C
B y_ C B v cos sin C
B _ C B v sin C
B C B C
B C=B L C:
B _rC B v cos sin( r) C
B C B
B Lr C
C
@ v_ A @ u1 A
_ u2
Modélisation des systèmes 31
La …gure représente le voilier que nous voulons modéliser. Le vecteur d’état x, de dimension 7, est composé
– des coordonnées de position, c’est-à-dire les coordonnées x; y du centre de gravité G du bateau (la dérive
se trouve en G), l’orientation , et les angles v et g de la voile et du gouvernail,
– et des coordonnées cinématiques v et ! représentant respectivement la vitesse du centre de gravité G et
la vitesse angulaire du bateau autour de G.
Les entrées u1 et u2 du système sont les dérivées des angles v et g. Les paramètres (supposés connus et
constants) sont
– V la vitesse du vent,
– rg la distance du gouvernail à G,
– rv la distance du mât à G,
– g la portance du gouvernail (si le gouvernail se trouve perpendiculaire à la marche du bateau, l’eau
exerce une force de g v Newton sur le gouvernail),
– v la portance de la voile (si la voile se trouve immobile, perpendiculaire au vent, ce dernier exerce une
force de v V Newton),
– f le coe¢ cient de frottement du bateau sur l’eau dans le sens de la marche (l’eau exerce sur le bateau
une force opposée au sens de la marche égale à f v),
– le coe¢ cient angulaire de frottement (l’eau exerce sur le bateau un couple de frottement égal à !;
étant donné la forme du bateau, plutôt pro…lé pour garder un cap, sera grand devant f ),
– J le moment d’inertie du bateau,
– ` la distance entre le centre de poussée de la voile et le mât,
– le coe¢ cient de dérive (lorsque la voile du bateau est relâchée, le bateau tend à dériver, dans le sens du
vent, à une vitesse égale à V ).
Cherchons à obtenir des équations d’état pour notre système, c’est-à-dire, une équation di¤érentielle de la
forme
x_ = f (x; u);
où x = (x; y; ; v ; g ; v; !)T et u = (u1 ; u2 )T : Pour cela, nous allons utiliser le principe fondamental de la
dynamique [14] en translation (a…n d’obtenir une expression de l’accélération tangentielle v) _ puis en rotation
32 Commande par espace d’état
Accélération tangentielle v_ : Le vent exerce sur la voile une force orthogonale à celle-ci et d’intensité
égale à
fv = v (V cos ( + v) v sin v) :
fg = g v sin g ;
orthogonale à celui-ci. La force de frottement qu’elle exerce sur le bateau est supposée visqueuse, c’est-à-dire
d’intensité proportionnelle à la vitesse. L’équation fondamentale de la dynamique, projetée suivant l’axe du
bateau donne
mv_ = sin v fv sin g fg f v:
L’accélération radiale peut être considérée comme nulle si nous supposons que la dérive est parfaite.
L’accélération angulaire !_ : Parmi les forces qui agissent sur la rotation du bateau, on trouve les forces fv
et fg exercées par la voile et le gouvernail, mais aussi une force de frottement angulaire que nous supposons
visqueuse. L’équation fondamentale de la dynamique nous donne
J !_ = dv fv dg fg !;
où (
dv = ` rv cos v
dg = rg cos g :
8
>
> x_ = v cos ; (i)
>
>
>
> y_ = v sin V; (ii)
>
>
>
> _ = !; (iii)
>
>
>
> _v
>
< = u1 ; (iv)
_g = u2 ; (v) (2.8)
>
> fv sin v fg sin g fv
>
> v_ = ; (vi)
>
> m
>
> !_ =
(` rv cos v )fv rg cos g fg !
; (vii)
>
> J
>
>
>
> fv = v (V cos ( + v ) v sin v ) ; (viii)
>
: fg = g v sin g : (ix)
Notons que les deux dernières équations ne sont pas di¤érentielles mais algébriques. A…n de coller complè-
tement avec une équation d’état comme dé…nie par l’équation (1.1) de la page 13, il conviendrait de faire
disparaître ces deux équations ainsi que les deux forces internes fv et fg apparaissant dans les équations (vi)
et (vii): Nous allons toutefois les garder a…n de conserver des équations relativement simples et de pouvoir
visualiser ces forces lors de nos simulations. Nous reviendrons plus tard à cet exemple du bateau à voile
pour sa régulation.
Modélisation des systèmes 33
Dans le bac de gauche, l’eau s’écoule sans frottement en direction du bac de droite. Dans le bac de gauche,
l’eau s’écoule de façon ‡uide, contrairement au bac de droite, où des turbulences existent. Ce sont ces
turbulences qui absorbent l’énergie cinétique de l’eau et la transforme en chaleur. Sans ces turbulences, nous
aurions un éternel mouvement de va et vient de l’eau entre les deux bacs. Dans le but d’appliquer la relation
de Bernouilli dans le bac de gauche, considérons un tube de courant, c’est à dire, une tube virtuel (voir
…gure) à l’intérieur duquel l’eau possède un mouvement ‡uide et ne traverse pas les parois. La relation de
Bernouilli nous dit que dans ce tube, en tout point,
v2
P+ + gz = constante;
2
où P est la pression au point considéré, z son altitude, v la vitesse de l’eau en ce point. Le coe¢ cient est la
masse volumique de l’eau et g est la constante de gravitation. D’après la relation de Bernouilli, nous avons
2
vD v2
PD + + gzD = PA + A + gzA ;
2 2
c’est-à-dire
2
vD
PD = PA + g (zA zD ) : (2.9)
2
De plus, nous pouvons supposer que C est loin de la zone de turbulence et que l’eau ne bouge pas. Ainsi,
nous avons, d’après Bernouilli
PC + gzC = PB + gzB ;
34 Commande par espace d’état
c’est-à-dire
PC = PB + g (zB zC ) : (2.10)
Or, dans cette zone de turbulence, l’eau se trouve ralentie, mais nous pouvons supposer que la pression de
varie pas, c’est-à-dire, PC = PD . Ainsi, d’après (2.9) et (2.10), nous avons
2
vD
PB + g (zB zC ) = PA + g (zA zD ) :
|{z} |{z} 2
Pa tm Pa tm
Dans le cas ou le niveau du bac droit est supérieur à celui du bac gauche, une étude similaire nous donne
p
vD = 2g (zB zA ):
Le signe moins de l’expression indique que le courant circule maintenant du bac droit vers le bac gauche.
Donc, la relation générale pour la vitesse de l’eau dans le canal est
p
vD = sign (zA zB ) 2gjzA zB j: (2.11)
Si a est la section du canal, le débit d’eau du bac droit vers le bac gauche est
p
QD = a:sign (zA zB ) 2gjzA zB j: (2.12)
Remarque. Initialement, cette loi a été démontrée dans un contexte plus simple où l’eau s’écoule dans
le vide. L’énergie totale d’un élement de ‡uide de masse m se conserve si on considére ce dernier comme
chutant librement dans le tube de courant. Ainsi, pour les deux points A et B, où A est à la surface et B
dans le tube, on a
1 2 1 2
mghA + mvA = mghB + mvB
2
| {z } 2
=0
et donc
p
vB = 2g (hA hB ):
On peut alors en déduire la relation (2.11). Mais ce raisonnement est plus délicat si on prend en considération
les forces de pressions tangencielles au tube se produisant dans le ‡uide.
L’eau des bacs 1 et 3 peut se déverser vers le bac 2, mais aussi vers l’extérieur se trouvant à pression
atmosphérique. Les débits associés sont, d’après la relation (2.12), donnés par
p
Q1ext = a: 2gh1 ;
p
Q3ext = a: 2gh3 :
Les variables d’état de ce système qui peuvent être considérées sont les hauteurs dans les bacs. Pour simpli…er,
nous supposerons que la surface des bacs sont toutes égales à 1 m2 , ainsi, le volume d’eau dans un bac se
confond avec la hauteur. Les équations d’état sont obtenues en écrivant que la variation du volume d’eau
dans un bac est égale à la somme des débits entrants moins la somme des débits sortants, c’est-à-dire,
h_ 1 = Q1ext Q12 + u1 ;
_h2 = Q12 Q23 ;
h_ 3 = Q3ext + Q23 + u2 ;
ou encore
p p
h_ 1 = a: 2gh1 a:sign (h1 h2 ) 2gjh1 h2 j + u1
p p
h_ 2 = a:sign (h1 h2 ) 2gjh1 h2 j a:sign (h2 h3 ) 2gjh2 h3 j
p p
h_ 3 = a: 2gh3 + a:sign (h2 h3 ) 2gjh2 h3 j + u2 :
On considère le vérin pneumatique avec ressort de rappel représenté sur la …gure. Un tel vérin est souvent
quali…é de simple e¤et car l’air sous pression n’existe que dans une des deux chambres.
36 Commande par espace d’état
Ce type de vérin est d’usage courant en robotique pour bouger les éléments mécaniques. Les paramètres de
ce système sont la raideur du ressort k, la surface du piston a et la masse m en bout de piston (les masses
de tous les autres objets sont négligées). On suppose que tout se passe à température constante T0 . Nous
prendrons pour vecteur d’état x = (z; z; _ p) où z est la position du vérin, z_ sa vitesse et p la pression dans
la chambre. Pour simpli…er, nous supposerons que le vide règne dans la chambre du ressort et que lorsque
pour z = 0 (le vérin est en butée gauche) le ressort se trouve en position d’équilibre. L’entrée du système
est le débit volumique u d’air vers la chambre du vérin. Nous avons donc
V
u= n;
_ (2.13)
n
où n est le nombre de moles de gaz dans la chambre et V est le volume de la chambre. Pour bien comprendre
la relation (2.13), il su¢ t de considérer le petit volume d’air V qui s’apprête à entrer dans la chambre entre
l’instant t et t + t. la loi des gaz parfait appliquée au petit volume et à la chambre, tous les deux à la même
pression, nous donne
p: V = n:R:T0 et p:V = n:R:T0 .
n n V
D’où V = V; c’est-à-dire V = n n. Donc
V V n V
u, = = n;
_
t n t n
qui se trouve être la relation (2.13). Attention, ici, nous avons approximé nt par n_ car n correspond bien
à un accroissement de n. En revanche, nous n’aurions pas pu écrire V_ = Vt car V ne correspond pas
nécessairement à un accroissement du volume V . On peut par exemple concevoir un V positif alors que
le volume V diminue. Le principe fondamental de la dynamique nous donne pa kz = m• z . Donc, les deux
premières équations d’état sont (
z_ = z;
_
ap kz
z• = m :
La loi des gaz parfaits (pV = nRT ) est donnée par pza = nRT . En di¤érenciant, on obtient
a (pz
_ + pz) _ + nT_ :
_ = R nT
En isolant p,
_ on obtient la troisième équation d’état pour notre système, à savoir
p u
p_ = z_ :
z a
Les équations d’état du système sont donc
8
>
< z_ = z_
ap kz
z• = m
>
: p_ = p u
z a z_
Le système de la …gure ci-dessous a pour entrée la tension u(t) et pour sortie la tension y(t). Cherchons à
modéliser ce système.
Soit i1 le courant électrique dans la résistance R1 (de haut en bas). D’après les lois des mailles et des nœuds,
on a 8
>
< u(t) v(t) R1 i1 (t) = 0 (loi des mailles)
di
L dt + R2 i(t) R1 i1 (t) = 0 (loi des mailles)
>
: i(t) + i1 (t) C dv
dt = 0 (loi des nœuds).
données par les quantités i(t) (proportionnelle au ‡ux) et v(t) (proportionnelle à la charge). On obtient les
di
équations d’état en cherchant à se débarrasser de i1 dans les équations précédentes et en isolant dt et dv
dt .
Bien sûr, une équation doit disparaître. On obtient
(
dv 1 1 1
dt = CR1 v(t) + C i(t) + CR1 u(t)
di 1 R2 1
dt = L v(t) L i(t) + L u(t):
Or, la sortie est donnée par y(t) = R2 i(t). Finalement, nous arrivons à la représentation d’état d’un système
linéaire donnée par
! ! ! !
1 1 1
d v(t) CR C v(t) CR
dt = 1
1
R2 + 1
1 u(t)
i(t) i(t)
L L ! L
v(t)
y(t) = 0 R2 :
i(t)
Un moteur à courant continu peut être décrit par la …gure suivante, où u est la tension d’alimentation du
moteur, i est le courant absorbé par le moteur, R est la résistance de l’induit, L est l’inductance de l’induit,
e est la force électromotrice, est le coe¢ cient de frottement dans le moteur, ! est la vitesse angulaire du
moteur et Tr est le couple exercé par le moteur sur la charge.
Dans le cas d’un moteur à excitation indépendante, ou à aimants permanents, est constant. On alors les
équations d’état linéaires suivantes
(
di R K u
dt = Li L !+ L
!_ = KJ i J ! TJr
Modélisation des systèmes 39
où les entrées sont u et Tr et les variables d’état sont i et !. Ces équations peuvent s’écrire sous forme
matricielle ! ! ! ! !
R K 1
d i L L i L 0 u
= K + 1 :
dt ! J J ! 0 J Tr
Dans le cas d’un moteur à excitation série, le ‡ux est proportionnel à i, c’est-à-dire = ki. On obtient
les équations d’état suivantes
(
di R Kk u
dt = Li L i! + L
Tr
!_ = Kk J i
2
J! J
qui sont devenues non-linéaires. De tels moteurs sont utilisés lorsque l’on souhaite un fort couple au démar-
rage, comme par exemple pour démarrer les moteurs à explosion.
Lorsque le moteur fonctionne, le couple Tr devient imposé. La table suivante donne quelques caractéristiques
mécaniques, en régime permanent, du couple (Tr ; !).
En fonctionnement, un moteur est donc un système qui ne possède plus qu’une seule entrée qui est la tension
d’induit u(t).
Bien que dans ce cours, nous ayons décidé de nous limiter aux systèmes à temps continu, nous allons, dans ce
paragraphe, modéliser un exemple de système à temps discret. Ceci permettra d’illustrer que la méthodologie
utilisée pour les systèmes à temps continus peut s’étendre aisément au cas discret.
Il s’agit d’étudier l’évolution du nombre y(k) de couples de lapins dans un élevage en fonction de l’année k.
L’année 0, il y a un seulement un couple de lapins nouveau-nés dans l’élevage (donc y(0) = 1). Les lapins ne
deviennent fertiles que un an après leur naissance. Donc à l’année 1, il y a toujours un seul couple de lapins,
mais ce couple est fertile (donc y(1) = 1). Un couple fertile donne naissance chaque année à un autre couple
de lapins. Donc à l’année 2, il y a un couple de lapin fertile et un couple de nouveaux-nés. Cette évolution
peut être décrite par le tableau suivante, où N signi…e nouveau né et A signi…e adulte.
Appelons x1 (k) le nombre de couples nouveaux-nés, x2 (k) le nombre de couples fertiles et y(k) le nombre
total de couples. Les équations d’état sont données par
8
>
< x1 (k + 1) = x2 (k)
x2 (k + 1) = x1 (k) + x2 (k)
>
: y(k) = x1 (k) + x2 (k);
avec x1 (0) = 1 et x2 (0) = 0 pour conditions initiales. Ce système est appelé système de Fibonacci. Cherchons
maintenant à trouver la relation de récurrence associée à ce système. Il nous faut pour cela exprimer y(k);
y(k + 1) et y(k + 2) en fonction de x1 (k) et x2 (k): Le calcul résultant est le suivant
C’est-à-dire 0 1 0 1
y(k) 1 1 !
B C B C x1 (k)
@ y(k + 1) A = @ 1 2 A :
x2 (k)
y(k + 2) 2 3
En éliminant x1 (k) et x2 (k) de ce système de trois d’équations linéaires, on obtient plus qu’une seule équation
donnée par
y(k + 2) y(k + 1) y(k) = 0:
Les conditions initiales sont y(0) = y(1) = 1: C’est généralement sous cette forme qu’est décrit le système
de Fibonacci.
Chapitre 3
Simulation
3.1 Introduction
Dans ce chapitre, nous allons montrer comment e¤ectuer une simulation sur ordinateur d’un système non
linéaire décrit par ses équations d’état
(
x(t)
_ = f (x(t); u(t));
y(t) = g(x(t); u(t)):
Cette étape est importante pour tester le comportement d’un système (régulé ou non). De plus, les tech-
niques que nous présenterons aux chapitres 5 et 6 nous permettront d’obtenir des équations d’état pour les
régulateurs que nous concevrons. L’ordinateur pourra nous être utile pour réaliser nos régulateurs et devra
pour cela simuler ces derniers.
Avant de présenter la méthode de simulation, nous introduirons dans les paragraphes 3.2 et 3.3 la notion
de champ de vecteur. Cette notion nous permettra de mieux comprendre la méthode de simulation qui sera
présentée au paragraphe 3.5 ainsi que certains comportements pouvant apparaître dans les systèmes non
linéaires. Au paragraphe 3.4, nous donnerons quelques notions de graphisme nécessaires à la représentation
graphique de nos systèmes. L’intérêt de l’interactivité dans les programmes de simulation ainsi que la façon
de l’implémenter seront présentés au paragraphe 3.6.
Dans ce paragraphe, nous allons présenter la notion de champ de vecteur et montrer en quoi elle est utile
pour mieux comprendre les di¤érents comportements des systèmes. Nous invitons le lecteur à consulter
l’ouvrage de Lamnabhi-Lagarrigue [10] pour plus de détails à ce sujet.
Nous appellerons champ de vecteur une fonction continue f de Rn vers Rn : Lorsque n = 2, une représentation
graphique de la fonction f peut être imaginée. Nous allons maintenant traiter deux exemples a…n d’illustrer la
notion de champ de vecteur. Le premier traite le cas d’une fonction linéaire et nous permettra de donner une
41
42 Commande par espace d’état
interprétation graphique de la notion de vecteur propre. Pour le deuxième exemple, la fonction f représente
le gradient d’une fonction à valeur dans R.
Considérons la fonction
R2 ! ! R2 !
f: x1 x1 + x2
! :
x2 x1 x2
Une illustration graphique de cette fonction est donnée par la …gure ci-dessous. Pour obtenir cette …gure,
nous avons pris un ensemble de vecteurs de l’ensemble de départ, suivant un maillage. Puis, pour chacun des
vecteurs du maillage x, nous avons dessiné son vecteur image f (x) en lui donnant pour origine le vecteur x.
Le programme Scilab qui nous a permis de faire ce dessin est donné par le tableau 3.1. Au pas 2, on dé…nit
le maillage pour x1 et x2 , a représente la précision du maillage. Aux pas 3 et 4, on dé…nit les fonctions
coordonnées f1 (u; v) = u + v et f2 (u; v) = u v de notre fonction f . Ici, u et v sont des variables muettes.
Le pas 5 e¤ectue l’évaluation de f sur tout le maillage. En…n le pas 6 trace le champ de vecteur.
Essayons maintenant d’interpréter le champ de vecteurs représenté sur la …gure. La matrice de notre appli-
cation linéaire est donnée par !
1 1
A= :
1 1
p p
Ses valeurs propres sont 2 et 2 et les vecteurs propres associés sont
! !
0:9239 0:3827
v1 = et v2 = :
0:3827 0:9239
Simulation 43
Programme [Link]
1 a=0.5
2 x1=-2 :a :2 ; x2=-2 :a :2 ;
3 deff(’z=f1(u,v)’,’z=u+v’) ;
4 deff(’z=f2(u,v)’,’z=u-v’) ;
5 y1=feval(x1,x2,f1) ; y2=feval(x1,x2,f2) ;
6 champ(x1,x2,y1,y2) ;
Pour obtenir ce résultat sous Scilab, il faut d’abord saisir la matrice en faisant A=[1 1 ;1 -1]. Puis on
exécute l’instruction [D,V]=bdiag(A). La matrice D contient la matrice diagonale (et donc les valeurs propres
de A) alors que V est la matrice de passage dont les colonnes sont les vecteurs propres. Notons que le vecteur
x représenté sur la …gure n’est pas un vecteur propre car x et f (x) ne sont pas colinéaires. En revanche, tous
les vecteurs appartenant aux sous-espaces propres (représentés en pointillé sur la …gure) sont des vecteurs
propres. Le long du sous-espace propre associé à la valeur propre négative, les vecteurs du champ tendent à
pointer vers 0 alors que ces vecteurs pointent vers l’in…ni le long du sous-espace propre associé à la valeur
propre positive.
Considérons la fonction
R2 ! ! R
f: x1
! x1 exp x21 x22 :
x2
Le champ de vecteur associé est donné par la …gure ci-dessous. Notons que le gradient s’annule en deux
points, représentés par les petits disques noirs de la …gure. Celui de droite correspond au maximum de la
fonction f (le maximum se trouve là où convergent les ‡èches) alors que celui de gauche correspond au
minimum (là où les ‡èches divergent). Notons que si nous simulons le système
! !
x_ 1 (1 2x21 ) exp( x21 x22 )
= ;
x_ 2 2x1 x2 exp( x21 x22 )
!
nous avons une forte chance de converger vers le maximum de la fonction f . Simuler le système x_ = grad
f (x), dans le but de trouver le maximum de f correspond à la méthode du gradient.
44 Commande par espace d’état
Considérons un système autonome (c’est-à-dire sans entrée) dont l’évolution est donnée par l’équation x(t)
_ =
2 2
f (x(t)). Lorsque f est une fonction de R vers R , nous pouvons obtenir une représentation graphique de
f en traçant le champ de vecteur associé à f . Le graphique nous permettra alors de mieux comprendre le
comportement de notre système. A titre d’illustration nous allons considérer deux exemples : le système
proies-prédateurs et le pendule simple.
où = 0:01 et = 0:02. Les variables d’état x1 (t) et x2 (t) représentent la taille des populations des proies
et des prédateurs. Le terme quadratique représente des interactions entre les deux espèces. Notons que la
population des proies augmente de façon exponentielle lorsque qu’il n’y a pas de prédateur. De même, la
population des prédateurs décroît lorsqu’il n’y a pas de proie. Pour tracer le champ de vecteur associé à la
fonction d’évolution !
(1 x2 (t)) x1 (t)
f (x) = ;
( 1 + x1 (t)) x2 (t)
Simulation 45
Programme [Link]
a=10 ; x1=0 :a :100 ; x2=0 :a :200 ;
deff(’z=f1(u,v)’,’z=(1-0.01*v)*u’) ;
deff(’z=f2(u,v)’,’z=(-1+0.02*u)*v’) ;
y1=feval(x1,x2,f1) ; y2=feval(x1,x2,f2) ;
champ(x1,x2,y1,y2) ;
Tableau 3.2 –Programme Scilab pour tracer le champ de vecteur associé au système proies-prédateurs
sur le pavé [0; 100] [0; 200]; on tape le programme du tableau 3.2 et on obtient la …gure ci-dessous.
Cette …gure illustre le comportement du système. En e¤et, l’évolution du vecteur d’état du système se fait
dans le sens des ‡èches (car x_ = f (x)). Ainsi, nous pouvons comprendre que l’évolution de système est
périodique et que l’état parcourt une courbe quasi-circulaire (dont le centre est (50; 100)) dans le sens direct
trigonométrique. Le point (50; 100) semble être un point d’équilibre pour notre système. Véri…ons-le par
calcul. Les points d’équilibre se trouvent là où le champ de vecteur s’annule. Il y en a deux, représentés par
les petits disques noirs sur la …gure. Ils satisfont l’équation f (x) = 0, c’est-à-dire
(
(1 x2 )x1 = 0;
( 1 + x1 )x2 = 0:
! !
x1 0
=
x2 0
qui correspond à une situation où aucune des deux espèces n’existe. Le deuxième est donné par
! ! !
1
x1 50
= 1 = ;
x2 100
Programme [Link]
a=0.4 ; x1=-7 :a :7 ; x2=-4 :a :4 ;
L=2 ;g=9.81 ;m=1 ;
deff(’y1=f1(x1,x2)’,’y1=x2’) ;
deff(’y2=f2(x1,x2)’,’y2=-g*sin(x1)’) ;
y1=feval(x1,x2,f1) ; y2=feval(x1,x2,f2) ;
champ(x1,x2,y1,y2) ;
La notion de point d’équilibre sera revue de façon plus approfondie dans le paragraphe 6.2, page 97.
Considérons le pendule simple modélisé au chapitre précédent au paragraphe 2.3.4 page 21. On suppose ici
que le couple d’entrée est nul. Avec Scilab, traçons le champ de vecteur associé à la fonction d’évolution
f (x), pour u = 0, m = 1, g = 9:81 et L = 1: Le programme [Link], qui e¤ectue cette tâche, est
donné sur le tableau 3.3.
Les petits disques noirs représentent les points d’équilibre. Le premier, le troisième et le dernier correspondent
à la situation où le pendule est en bas, en équilibre stable. Le deuxième et le quatrième correspondent à la
situation où le pendule est en haut, en équilibre instable. Autour d’un point d’équilibre stable, le vecteur
d’état tend à tourner autour de ce point, formant ainsi un cycle. Si on initialise le pendule en haut de
la …gure (par exemple en ( 7; 3)), ce qui correspond à une situation où _ est élevé, la variable tendra
à croître indé…niment. Cela correspond à une situation où le pendule tourne toujours dans le sens direct
trigonométrique autour de son axe et sans jamais s’arrêter.
Dans ce paragraphe, nous donnons quelques notions nécessaires à la représentation graphique des systèmes
lors des simulations.
3.4.1 Motif
Un motif est une matrice à deux ou trois lignes (suivant que l’objet est dans le plan ou l’espace) et n colonnes
qui représentent les n sommets d’un polygone indéformable, censé représenter l’objet. Il est important que
l’union de tous les segments formés par deux points consécutifs du motif forme toutes les arêtes du polygone
que l’on souhaite représenter.
Exemple 3.4.1 Le motif M du châssis (voir …gure ci-dessous) de la voiture (avec les roues arrières) est
donné par
!
1 4 5 5 4 1 1 0 0 1 1 0 0 1 1 0 0 3 3 3
(3.1)
2 2 1 1 2 2 2 2 3 3 3 3 3 3 3 3 2 2 3 3
Il est clair que sur le dessin de la voiture en mouvement, les roues avant peuvent bouger par rapport au
châssis, mais aussi l’une par rapport à l’autre. Elles ne peuvent donc être incorporées au motif du châssis.
48 Commande par espace d’état
Pour le dessin de la voiture, nous devrons donc utiliser 3 motifs : celui du châssis, celui de la roue avant
gauche et celui de la roue avant droit. Sous Scilab, le motif M (ici en deux dimensions) peut être dessiné
par l’instruction
xpoly(M(1, :),M(2, :)).
Sous Scilab, M(i, :) renvoie la i ième ligne de la matrice M.
est associé au cube unité [0; 1]3 de R3 . Notons que ce motif dispose de 16 colonnes, alors qu’on aurait pu
s’attendre à 13 colonnes (en e¤ et, le cube possède 12 arêtes). Cela vient du fait que pour dessiner dans
l’espace toutes les arêtes d’un cube, sans lever le stylo, nous devons nécessairement passer par un minimum
de 16 sommets du cube.
Il est parfois utile de concevoir un modeleur a…n de faciliter l’élaboration de nos motifs. Voici un exemple
simple de modeleur
M=[ ] ;isoview(-10,10,-10,10) ;
for k=1 :10,
M=[M,round(locate(1))] ;
xpoly(M(1, :),M(2, :)) ;
endfunction
La fonction locate(1) retourne les coordonnées du point de la fenêtre graphique cliqué à la souris. Après
10 clics, la programme s’arrête. La matrice M contient alors les coordonnées des 10 points cliqués.
Rappelons que la jième colonne de la matrice d’une application linéaire de Rn ! Rn représente l’image du
jième vecteur ej de la base canonique. Ainsi, l’expression d’une matrice de rotation d’angle dans le plan
Simulation 49
En ce qui concerne les rotations dans l’espace R3 , il est important de préciser l’axe de rotation. Nous pouvons
distinguer 3 rotations principales : la rotation d’angle autour de l’axe Ox, celle autour de Oy et celle autour
de Oz. Les matrices associées sont respectivement données par
Rotation d’angle
0 1 0 1 0 1
1 0 0 cos 0 sin cos sin 0
B C B C B C
Rx = @ 0 cos sin A ; Ry = @ 0 1 0 A et Rz = @ sin cos 0 A: (3.4)
0 sin cos sin 0 cos 0 0 1
Le dessin d’objets bidimensionnels ou tridimensionnels sur un écran nécessite une suite de transformations
a¢ nes (rotations, translations, homothéties) de la forme
Rn ! Rn
fi :
x 7! Ai x + bi
avec n = 2 ou 3. Or, la manipulation de compositions de fonctions a¢ nes est moins aisée que celle d’applica-
tions linéaires. L’idée la transformation en coordonnées homogènes est de transformer un système d’équations
a¢ nes en système d’équations linéaires. Remarquons tout d’abord qu’une équation a¢ ne du type y = Ax+b
peut se récrire ! ! !
y A b x
= :
1 0 1 1
Nous dé…nirons donc la transformation homogène d’un vecteur comme suit
!
x
x 7! xh = :
1
50 Commande par espace d’état
En nous aidant des matrices de la formule (3.4), nous pouvons écrire une fonction Scilab pour générer une
matrice homogène de rotation de R3 autour du iième axe de la base canonique.
function R=Rot(a,i)
R=eye(4,4) ;
J=modulo([i,i+1],3)+[1,1] ;
R(J,J)=[cos(a) -sin(a) ;sin(a) cos(a)] ;
endfunction
La fonction eye génère une matrice identité de dimension 4 4. La fonction modulo(a,b) renvoie le reste de
la division euclidienne de a par b (puisque ici a est un vecteur, ce calcul se fait composante par composante).
Le vecteur J correspond aux indices de la matrice à modi…er. Par exemple, si i = 1; alors J= [2; 3], si i = 2;
alors J= [3; 1], et si i = 3; alors J= [1; 2].
Une fonction pour générer une matrice homogène de translation d’un vecteur v dans R3 est donnée ci-dessous
function T=Trans(v)
T=eye(4,4) ;
T(1 :3,4)=v ;
endfunction
Pour dessiner le châssis de la voiture, il nous faut tout d’abord prendre le motif de la voiture (voir équation
(3.1)) et le rendre homogène, en lui rajoutant une ligne de 1. Le motif M ainsi obtenu est donné par
0 1
1 4 5 5 4 1 1 0 0 1 1 0 0 1 1 0 0 3 3 3
B C
@ 2 2 1 1 2 2 2 2 3 3 3 3 3 3 3 3 2 2 3 3 A
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Ensuite, il nous faut lui faire subir une rotation d’angle et une translation de vecteur (x; y). Pour il su¢ t
de multiplier M à gauche par la matrice
0 1
cos sin x
B C
R = @ sin cos y A:
0 0 1
Simulation 51
1 function draw_voiture(x)
2 M=[-1,4,5,5,4,-1,-1,-1,0,0,-1,
3 1,0,0,-1,1,0,0,3,3,3 ;
4 -2,-2,-1,1,2,2,-2,-2,-2,-3,-3,
5 -3,-3,3,3,3,3,2,2,3,-3 ;
6 ones(1 :21)] ;
7 Rav=[-1 1 ;0 0 ;1 1] ;
8 R=[cos(x(3)),-sin(x(3)),x(1) ;
9 sin(x(3)),cos(x(3)),x(2) ;0,0,1] ;
10 M=R*M ;
11 Ravd=R*[cos(x(5)),-sin(x(5)),3 ;
12 sin(x(5)),cos(x(5)),3 ;0 0 1]*Rav ;
13 Ravg=R*[cos(x(5)),-sin(x(5)) 3 ;
14 sin(x(5)),cos(x(5)) -3 ;0 0 1]*Rav ;
15 xbasc() ;xset(’thickness’,2) ;
16 isoview(-20,40,-10,30) ;
17 xpoly(M(1, :),M(2, :)) ;
18 xpoly(Ravd(1, :),Ravd(2, :)) ;
19 xpoly(Ravg(1, :),Ravg(2, :)) ;
20 endfunction ;
Pour dessiner la roue avant gauche, on dé…nit son motif (en coordonnées homogènes)
0 1
1 1
B C
M=@ 0 0 A:
1 1
On lui a¤ecte une rotation d’angle suivie d’une translation de (3; 3), puis à nouveau une rotation de
suivi d’une translation de (x; y). La matrice de transformation résultante est
0 10 1
cos sin x cos sin 3
B CB C
@ sin cos y A @ sin cos 3 A:
0 0 1 0 0 1
Une matrice similaire peut être obtenue pour la roue avant droite. Sur ces principes, une fonction Scilab
qui dessine la voiture dans un état donné peut être conçue. La fonction draw_voiture du tableau 3.4 en
donne un exemple. Les lignes 2 à 6 dé…nissent le motif du châssis de la voiture et des roues arrières. La ligne
7 dé…nit le motif d’une roue avant.
52 Commande par espace d’état
Le robot manipulateur représenté sur la …gure est composé de trois bras en série. Le premier, de longueur
3, peut pivoter en l’origine autour de l’axe Oz. Le second, de longueur 2, placé au bout du premier peut
lui aussi pivoter autour de l’axe Oz. Quant au troisième, de longueur 1, placé au bout du second, il peut
pivoter autour de l’axe formé par le second bras. Ce robot admet 3 degrés de liberté x = ( 1 ; 2 ; 3 ), où
les i représentent les angles formés par chacun des bras. Le motif de base choisi pour la représentation de
chacun des bras est le cube unité. Chacun des bras est supposé être un parallélépipède d’épaisseur 0.3. A…n
de prendre la forme du bras, le motif doit subir une a¢ nité, représentée par une matrice diagonale. Ensuite,
il doit subir une suite de rotations et de translations a…n de le positionner correctement.
Le tableau 3.5 représente la suite des transformations à faire subir au motif a…n de représenter chaque bras.
Comme indiqué sur ce tableau, le bras 2 doit subir, au …nal, toutes les transformations faites au bras 1 et
le bras 3 doit subir toutes celles faites au bras 2. Le tableau 3.6 donne le programme Scilab qui simule un
déplacement de ce bras.
Tableau 3.5 –Tableau des transformations à e¤ectuer pour chacun des bras
Simulation 53
Programme [Link]
function draw(x)
M=[0 1 1 0 0 0 1 1 0 0 0 0 1 1 1 1 ;
0 0 0 0 0 1 1 1 1 1 1 0 0 1 1 0;
0 0 1 1 0 0 0 1 1 0 1 1 1 1 0 0;
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1] ;
M=[1 0 0 0 ;0 1 0 -0.5 ; 0 0 1 -0.5 ;0 0 0 1]*M ;
M1=diag([3,0.3,0.3,1])*M ; R1=Rot(0.5,1)*Rot(x(1),3) ;
M2=diag([2,0.3,0.3,1])*M ; R2=Rot(x(2),3) ;
M3=diag([1,0.3,0.3,1])*M ; R3=Rot(x(3),1) ;
T2=Trans([3 ;0 ;0]) ; T3=Trans([2 ;0 ;0]) ;
M1=R1*M1 ; M2=R1*T2*R2*M2 ; M3=R1*T2*R2*T3*R3*M3 ;
xbasc() ;xset(’thickness’,2) ;
isoview(-4,4,-3,3) ;
xpoly(M1(1, :),M1(3, :)) ;
xpoly(M2(1, :),M2(3, :)) ;
xpoly(M3(1, :),M3(3, :))
endfunction ;
x=[0 ;1 ;0] ;w=[3 ;-1 ;7] ; dt=0.001
for k=1 :10000, x=x+dt*w ; draw(x) ;end ;
3.5 Simulation
Dans ce paragraphe, nous allons présenter la méthode d’Euler pour e¤ectuer une simulation sur ordinateur
d’un système non linéaire décrit par ses équations d’état
(
x(t)
_ = f (x(t); u(t));
(3.5)
y(t) = g(x(t); u(t)):
Cette méthode est assez approximative, mais reste simple à comprendre et su¢ t pour décrire les com-
portements de la plupart des systèmes. Après avoir brièvement présenté cette méthode, nous traiterons la
simulation d’un pendule simple puis d’une voiture.
Soit dt, un nombre très petit devant les constantes de temps du système et qui correspond à la période
d’échantillonnage de la méthode (par exemple dt = 0:01). L’équation d’évolution de l’équation (3.5) s’ap-
proxime par
x(t+dt) x(t)
dt ' f (x(t); u(t));
c’est-à-dire
x (t + dt) ' x (t) + f (x(t); u(t)):dt:
54 Commande par espace d’état
Cette équation peut être interprétée comme une formule de Taylor à l’ordre 1. On en déduit l’algorithme
de simulation (dit méthode d’Euler) du tableau 3.7. L’échantillonneur produit une interruption périodique
Algorithme Euler(in : x0 )
1 x := x0 ; t := 0; dt = 0:01;
2 répéter
3 saisir u à l’entrée de l’ordinateur ;
4 y := g(x; u);
5 sortir y;
6 x := x + f (x; u):dt;
7 attendre une interruption de l’échantillonneur ;
8 t = t + dt;
9 éternellement.
toutes les dt secondes. Ainsi, si l’ordinateur est su¢ samment rapide, la simulation se déroule à la même
vitesse que notre système physique. On parlera alors de simulation temps réel . Dans certaines circonstances,
ce qui nous intéresse est d’obtenir le résultat de la simulation le plus rapidement possible (par exemple pour
prédire comment va se comporter un système dans le futur). Dans ce cas, il n’est pas nécessaire de ralentir
l’ordinateur pour le synchroniser avec notre temps physique.
L’erreur locale et est ici d’ordre 2 et la méthode d’intégration est par conséquent beaucoup plus précise.
On considère à nouveau le pendule de longueur ` de masse négligeable et portant une masse m à son
extrémité (voir paragraphe 2.3.4 page 21): Rappelons, que, en son axe, on peut exercer un couple u et que
sa représentation d’état est donnée par
! !
x_ 1 x2
= g u :
x_ 2 ` sin x1 + m`2
Le programme Scilab [Link] du tableau 3.8 simule ce pendule par la méthode d’Euler.
Ce programme, qui débute en ligne 8, est composé des fonctions draw_pendule(x), qui dessine le pendule,
Simulation 55
Programme [Link]
1 function draw_pendule(x)
2 xbasc() ;isoview(-4,4,-3,3) ;
3 xpoly([0,L*sin(x(1))],[0 -L*cos(x(1))])
4 endfunction ;
5 function v=f(x,u)
6 v=[x(2) ; -g/L*sin(x(1))+u/(m*L^2)] ;
7 endfunction ;
8 L=2 ;g=9.81 ;m=1 ;u=0 ;dt=0.006 ;x=[1 ;0] ;
9 for k=1 :1500,
10 draw_pendule(x) ;
11 x=x+f(x,u)*dt ;
12 end ;
et de la fonction d’évolution f(x,u). Le pas 8 initialise les constantes et le vecteur initial. La boucle 9-12
lance la méthode d’Euler.
La …gure ci-dessous représente la trajectoire du pendule dans l’espace d’état pour la condition initiale = 1
et _ = 0. Sur cet exemple, le pendule, qui ne possède aucun frottement, devrait normalement e¤ectuer
un cycle. Or, la méthode d’Euler donne à chaque pas un peu d’énergie au pendule, ce qui explique sa
trajectoire qui tend à diverger. Cela vient du fait que le champ, bien que tangent à la trajectoire du système,
tend à sortir vers l’extérieur celle-ci. Le programme [Link], qui nous a permis
de dessiner la …gure ci-dessous, illustre ce phénomène. Il apparaît clairement que la méthode d’Euler ne
doit pas être utilisée pour simuler des systèmes conservatifs (comme le système planétaire), qui n’ont pas de
frottements et tendent à e¤ectuer des cycles. En revanche, pour les systèmes dissipatifs (avec frottements),
une méthode d’Euler s’avère souvent su¢ sante pour décrire correctement le comportement du système. Pour
notre pendule, la méthode d’Euler tend à donner de l’énergie au pendule, qui, après une dizaine d’oscillations,
se met à e¤ectuer des révolutions complètes autour de son axe, dans le sens indirect trigonométrique.
56 Commande par espace d’état
Considérons à nouveau la voiture modélisée au paragraphe 2.4.3 de la page 28. Nous allons ici tenter de
la simuler et visualiser graphiquement son mouvement. Rappelons que les équations d’état de cette voiture
sont données par 0 1 0 1
x_ v cos cos
B y_ C B v cos sin C
B C B C
B _ C B v sin C
B C=B L C:
B C B C
@ v_ A @ u1 A
_ u2
A un vecteur d’état donné x = (x; y; ; v; )T correspond une con…guration possible pour la voiture. Le
programme complet de la simulation de la voiture est donné par le tableau 3.9. En plus de la fonction
draw_voiture qui dessine la voiture, il appelle la fonction d’évolution du système. L’état initial de la
voiture est (0 0 0 7 0)T ; ce qui signi…e qu’à l’instant 0, la voiture est centrée à l’origine, avec un angle nul,
sa vitesse est de 7 ms 1 et les roues avant sont parallèles à l’axe de la voiture. La commande vectorielle u(t)
reste constante et égale à (0 0:2)T . Ce qui signi…e que la voiture n’accélère pas (car u1 = 0) et que le volant
tourne à une vitesse constante de 0:2 rad:s 1 (car u2 = 0:2).
Au moment de tester le programme de simulation, il est fondamental de pouvoir agir sur le programme
de façon à déceler des réactions anormales du système, à comprendre mieux son comportement, à détecter
Simulation 57
Programme [Link]
1 function v=f(x,u)
2 v=[x(4)*cos(x(5))*cos(x(3)) ;
3 x(4)*cos(x(5))*sin(x(3)) ;
4 x(4)*sin(x(5))/3 ; u(1) ;u(2) ] ;
5 endfunction ;
6 x=[0 ;0 ;0 ;7 ;0] ;u=[0 ;0.2] ;dt=0.02 ;
7 for k=1 :3000,
8 draw_voiture(x) ; x1=x+f(x,u)*dt ; x=x1 ;
9 end ;
d’éventuelles erreurs de programmation. Dans ce but, nous pouvons autoriser l’utilisateur à créer des événe-
ments (par exemple en appuyant sur une touche du clavier, ou en bougeant la souris) qui seront interprétées
par le programme grâce à un gestionnaire d’événement (eventhandler ). Le gestionnaire d’événement tourne
en parallèle avec le programme principal. Il attend qu’un événement se produise et s’exécute dès que cela
arrive. Le programme ([Link]) listé et commenté ci-dessous simule une voiture avec gestion
d’événement. Ce programme contient trois fonctions et un programme principal. Lorsque l’utilisateur ap-
puie sur une touche, la fonction my_eventhandler s’exécute immédiatement. Cette dernière modi…e une
variable globale (u, sortie ou ech) en permettant ainsi à l’utilisateur d’agir sur le programme. A…n de se
rapprocher d’un comportement temps réel de notre voiture, nous chronométrons, à l’aide d’un timer qui
compte en microsecondes, le temps de calcul d’un parcours de la boucle while. La valeur de ce temps de
calcul est alors a¤ectée à la variable dt.
//---------------------------------------------
function my_eventhandler(win,x,y,a)
//variable partagees avec le programme principal
global u sortie ech;
xinfo(’Volant:w,x; Acceleration:p,l; Zoom:z,Z; Fin:F’)
if ascii(a)==’l’ then u=[-500,0]; end //freinage
if ascii(a)==’p’ then u=[ 500,0]; end //acceleration
if ascii(a)==’w’ then u=[0,20]; end //a gauche
if ascii(a)==’x’ then u=[0,-20]; end //a droite
if ascii(a)==’z’ then ech=ech*0.95; end //zoom in
if ascii(a)==’Z’ then ech=ech/0.95; end //zoom out
if ascii(a)==’F’ then sortie=%T; end //fin
endfunction
//---------------------------------------------
function draw_voiture(x)
//Chassis
58 Commande par espace d’état
M=[-1 4 5 5 4 -1 -1 -1 0 0 -1 1 0 0 -1 1 0 0 3 3 3;
-2 -2 -1 1 2 2 -2 -2 -2 -3 -3 -3 -3 3 3 3 3 2 2 3 -3;
ones(1:21)] ;
Rav=[-1 1;0 0;1 1]; //motif d’une roue avant
R=[cos(x(3)),-sin(x(3)),x(1);
sin(x(3)), cos(x(3)),x(2);
0 0 1] ;
M=R*M;
Ravd=R*[cos(x(5)),-sin(x(5)) 3;
sin(x(5)),cos(x(5)) 3 ;
0 0 1]*Rav; // roue avant droite
Ravg=R*[cos(x(5)),-sin(x(5)) 3;
sin(x(5)),cos(x(5)) -3;
0 0 1]*Rav; // roue avant gauche
xbasc(); // efface l’ecran
xset(’thickness’,2); // regle l’epaisseur du stylo//
isoview(-ech,ech,-ech,ech); // regle les echelles
xpoly(M(1,:),M(2,:));
xpoly(Ravd(1,:),Ravd(2,:))
xpoly(Ravg(1,:),Ravg(2,:))
endfunction;
//---------------------------------------------
//---------------------------------------------
//debut du programme principal
//---------------------------------------------
while sortie==%F,
pvm_set_timer(); //initialise le chronometre
draw_voiture(x);
x=x+f(x,u)*dt; //methode d’Euler
u=[0;0];
//converti en seconde la valeur du chronometre
dt=pvm_get_timer()/1000000;
end;
seteventhandler(’ ’) //supprime le gestionnaire d’evenement
xdel(); //ferme la fenetre graphique
//---------------------------------------------
60 Commande par espace d’état
Chapitre 4
Systèmes linéaires
4.1 Introduction
Dans ce chapitre, nous allons étudier les systèmes linéaires (voir page 14, formule (1.3) et (1.4)). Les notions
présentées seront fondamentales pour la bonne compréhension des chapitres qui suivent et plus particu-
lièrement pour la conception de régulateurs linéaires présentée au chapitre 5. Les systèmes linéaires sont,
rappelons-le, de la forme
(
x(t)
_ = Ax(t) + Bu(t)
(4.1)
y(t) = Cx(t) + Du(t);
(
x(k + 1) = Ax(k) + Bu(k)
(4.2)
y(k) = Cx(k) + Du(k);
pour les systèmes à temps discret. Nous allons tout d’abord donner une solution analytique pour ces équations
d’état au paragraphe 4.2. Ces solutions seront utilisées au paragraphe 4.2 pour établir des critères de stabilité.
Ensuite, au paragraphe 4.4, nous montrerons que la représentation d’état d’un système linéaire n’est pas
unique et comment il est possible de passer d’une représentation d’état à une autre. En…n aux paragraphes
4.6 et 4.7, nous expliciterons les liens forts qui existent entre la représentation entrées-sorties (utilisant la
notion de matrice de transfert) et la représentation d’état.
Dans ce paragraphe, nous allons chercher à calculer sous forme analytique la solution des équations d’état.
Rappelons que ces dernières sont composées de l’équation di¤érentielle (4.1), dans le cas continu, ou de
l’équation de récurrence (4.2) dans le cas discret.
61
62 Commande par espace d’état
Rappel sur les exponentielles de matrices. A…n de résoudre notre équation di¤érentielle, nous allons
utiliser la notion d’exponentielle de matrice, que nous allons brièvement rappeler. L’exponentielle d’une
matrice carrée M de dimension n se dé…nit par son développement en séries entière :
1
X
1 2 1 1 i
eM = In + M + M + M3 + = M;
2! 3! i!
i=0
où In est la matrice identité de dimension n. Il est clair que eM est de la même dimension que M. Voici
quelques-unes des propriétés importantes concernant les exponentielles de matrices. Si 0n est la matrice
nulle de dimension n n et si M et N sont deux matrices n n, alors
e0n = In ;
eM :eN = eM+N ; (si les matrices commutent)
d Mt = MeMt :
dt e
Théorème 4.2.1 Notons x(0), l’état à l’instant initial t = 0. Pour le système linéaire à temps continu
(4.1), l’état à l’instant t est donné par
Z t
At
x(t) = e x(0) + eA(t ) Bu( )d ; (4.3)
0
Rt ) Bu(
La fonction CeAt x(0) est appelée solution homogène, libre ou transitoire. La fonction 0 CeA(t )d +
Du(t) est appelée solution forcée.
Preuve. Posons z(t) = e At x(t). On a donc x(t) = eAt z(t) et donc, par dérivation x(t)
_ = AeAt z(t)+eAt z_ (t).
L’équation d’évolution x(t)
_ = Ax(t) + Bu(t); se transforme en
Donc Rt
x(t) = eAt z (0) + 0 e A Bu( )d
Rt
= eAt z (0) + 0 eAt e A Bu( )d
Rt
= eAt x(0) + 0 eA(t ) Bu( )d :
Pour obtenir la sortie, il su¢ t de noter que y(t) = Cx(t) + Du(t):
Systèmes linéaires 63
Théorème 4.2.2 Notons x(0) l’état à l’instant initial k = 0. Pour le système linéaire à temps discret (4.2),
l’état à l’instant k est donné par
k 1
X
k
x(k) = A x(0) + Ak ` 1
Bu(`); (4.4)
`=0
Pk
De même que dans le cas continu, la fonction CAk x(0) est la solution homogène `=0 CA
k `
Bu(`)+Du(k)
est la solution forcée.
Preuve. La preuve peut se faire aisément par récurrence. Tout d’abord, si k = 0; la relation (4.4) est
véri…ée. Véri…ons que si elle est véri…é pour k elle l’est aussi pour k + 1: On a
La relation est donc aussi valide pour k+1. Pour obtenir la sortie, il su¢ t de noter que y(k) = Cx(k)+Du(k):
4.3 Stabilité
4.3.1 Dé…nition
Un système linéaire est stable (appelé aussi asymptotiquement stable dans la littérature) si, au bout d’un
temps su¢ samment long, l’état ne dépend plus des conditions initiales et ceci, quelles que soient ces dernières.
Cela revient à dire que
limt!1 eAt = 0n si le système à temps continu,
(4.5)
limk!1 Ak = 0n si le système à temps discret.
Nous allons maintenant chercher à proposer un critère de stabilité ne dépendant que de la matrice A. Pour
cela, on rappelle tout d’abord le théorème de correspondance des valeurs propres. Ensuite, nous donnerons
le critère de stabilité.
64 Commande par espace d’état
Théorème 4.3.1 (de correspondance des valeurs propres). Si f est un polynôme (ou plus généralement une
série entière) et si A est une matrice de Rn n . Les vecteurs propres de A sont aussi des vecteurs propres pour
f (A). De plus si les valeurs propres de A sont f 1 ; : : : ; n g alors celles de f (A) sont ff ( 1 ) ; : : : ; f ( n )g :
Preuve. Soit x un vecteur propre de A de valeur propre . Nous allons montrer par récurrence que x est
aussi un vecteur propre de f (A) de valeur propre f ( ). Tout d’abord, cette propriété est vraie si f (A) = Ai .
En e¤et, puisque
Ai :x = A:A: : : : ::A:x = i
:x:
nous avons bien la propriété f (A):x = f ( ):x. Supposons que cette propriété soit vraie pour deux polynômes
f1 et f2 ; nous allons maintenant montrer qu’elle l’est aussi pour f1 + f2 et f1 . Puisque qu’elle est vraie
pour f1 et f2 , on a
f1 (A):x = f1 ( ):x
f2 (A):x = f2 ( ):x
et donc
Par récurrence, nous pouvons en déduire que la propriété est vraie pour toutes les fonctions f (A) qui peuvent
être générées à partir des Ai par compositions d’additions et de multiplications par un scalaire, c’est-à-dire
pour les fonctions f qui sont des polynômes.
Théorème 4.3.2 Un système linéaire à temps continu est stable si et seulement si toutes les valeurs propres
de sa matrice d’évolution A sont à parties réelles strictement négatives. Un système linéaire à temps discret
est stable si et seulement si toutes les valeurs propres de A sont strictement dans le cercle unité.
P
Preuve. Commençons par le cas du temps continu. Si on prend f (A) = eAt = 1 1 i
i=0 i! (At) ; le théorème
4.3.1 de correspondance des valeurs propres (que nous supposerons s’appliquer même pour des polynômes de
degré in…ni, c’est-à-dire, des séries entières), nous dit que les valeurs propres de eAt sont de la forme e i t , où
At = 0
i denote la iième valeur propre de A. Or, la stabilité du système se traduit par la condition limt!1 e n
(voir (4.5)). Nous avons
t!1 it
t!1
eAt ! 0n , 8i 2 f1; : : : ; ng ; e ! 0 (continuité de l’exponentielle et des valeurs propres)
t!1
, 8i 2 f1; : : : ; ng ; jje(Re i +i Im i )t jj ! 0
, 8i 2 f1; : : : ; ng ; Re ( i ) < 0
Dans le cas discret, il nous faut prendre f (A) = Ak . On aura stabilité si limk!1 Ak = 0n , c’est-à-dire
si toutes les valeurs propres ki de la matrice Ak convergent vers zéro, c’est-à-dire si toutes les i sont de
module strictement inférieur à un.
Systèmes linéaires 65
Que ce soit pour les systèmes à temps continu (4.1) ou à temps discret (4.2), la position des valeurs propres de
A est d’une importance capitale pour l’étude de la stabilité d’un système linéaire. Le polynôme caractéristique
d’un système linéaire est dé…ni comme étant le polynôme caractéristique de la matrice A qui est donné par
la formule
P (s) = det (sIn A) :
Ses racines sont les valeurs propres de A. En e¤et, si s est une racine de P (s), alors det (sIn A) = 0;
c’est-à-dire qu’il existe un vecteur v non nul tel que (sIn A) v = 0. Ce qui signi…e que sv Av = 0 ou
encore Av = sv. Donc s est une valeur propre de A. Un corollaire du théorème 4.3.2 est donc le suivant.
Corollaire 4.3.1 Un système linéaire à temps continu est stable si et seulement si toutes les racines de son
polynôme caractéristique sont à parties réelles négatives. Un système linéaire à temps discret est stable si et
seulement si toutes les racines de son polynôme caractéristique sont dans le cercle unité.
4.4.1 Principe
Considérons le système linéaire à temps continu donné par ses équations d’état
(
x_ = Ax + Bu
(4.6)
y = Cx + Du:
Tous les développements qui vont suivre s’appliquent de la même façon au cas des systèmes linéaires à temps
discret. Posons v = P 1 x, où P est une matrice de passage (c’est-à-dire carrée et inversible). Remplaçons
x par Pv. On obtient (
Pv_ = APv + Bu
y = CPv + Du;
c’est-à-dire (
v_ = P 1 APv + P 1 Bu
y = CPv + Du:
Si on pose
1 1
A=P AP; B = P B; C = CP et D = D; (4.7)
le système s’écrit (
v_ = Av + Bu
(4.8)
y = Cv + Du
qui est bien une représentation d’état. Ainsi, un système linéaire possède autant de représentations d’état
qu’il existe de matrices de passage. Bien sûr, certaines représentations sont préférables suivant l’application
visée.
66 Commande par espace d’état
Remarque 4.4.1 Si nous disposons de deux représentations d’état données par (4.6) et (4.8), censées
représenter un même système, la matrice de passage qui relie ces deux représentations peut être obtenue en
résolvant le système linéaire
PA = AP;
PB = B;
C = CP;
où les n2 inconnues sont les coe¢ cients pij de la matrice de passage P. Puisque ce système est linéaire
en les coe¢ cients pij , sa résolution est aisée. Si aucune solution n’est obtenue, il est clair que les deux
représentations ne sont pas équivalentes.
Considérons un système avec une entrée. Ainsi, la matrice B deviendra un vecteur b. Prenons pour matrice
de changement de base
P = b j Ab j A2 b j : : : j An 1 b ; (4.9)
qui sera supposée inversible. Notons que si ei désigne le vecteur ne contenant que des zéros, sauf un 1 à la
position i, on a Pei = Ai 1 b; soit
ei = P 1 :Ai 1 b: (4.10)
Ainsi,
(4.7) 1 (4.10)
b = P b = e1 :
Pour l’obtention de A, nous écrivons,
A = a1 j a2 j a3 j : : : j an 1 j an
(4.7) 1
= P AP
(4.9) 1 1 2 1 3 1 n 1 1 n
= P Ab j P A bjP A b j ::: j P A bjP A b
(4.10) 1 n
= e2 j e3 j e4 j : : : j en j P A b ;
où les ai désignent les iièmes colonnes de A. Or, le théorème de Cayley-Hamilton nous dit que
An + an 1A
n 1
+ + a1 A + a0 I = 0:
La matrice d’évolution, est appelée matrice compagne . Une telle matrice possède une sous-diagonale com-
posée de un, et, sur sa dernière colonne, fait apparaître les coe¢ cients du polynôme caractéristique.
La transformée de Laplace est un outil, qui pour l’automaticien, est très utile pour la manipulation de
systèmes décrits par des équations di¤érentielles. Nous allons tout d’abord donner un exemple introductif
pour montrer l’intuition sous-jacente.
Première approche. On considère deux systèmes en parallèle comme illustré sur la …gure ci-dessous
Nous cherchons à obtenir l’équation di¤érentielle qui relie u à y. Ce système peut être décrit les équations
suivantes 8
>
< y_ 1 + y1 = u
y_ 2 + 2y2 = u (4.11)
>
: y +y
1 2 = y:
Il nous faut donc nous débarrasser des y1 et y2 . Par dérivation de équations (4.11), nous pouvons parvenir
aux équations suivantes 8
>
> y•1 + y_ 1 = u_
>
< y• + 2y_
2 2 = u_
>
> y_ 1 + y
_ 2 = y
_
>
:
y• y•1 y•2 = 0:
Nous avons en tout 7 équations avec 6 variables en trop : y1 ; y2 ; y_ 1 ; y_ 2 ; y•1 ; y•2 . Nous pouvons les éliminer par
une méthode de substitution et ainsi obtenir l’équation di¤érentielle recherchée :
y• + 3y_ + 2y = u:
Deuxième approche. Une autre approche serait de réécrire notre système (4.11) sous la forme :
8
> d
< dt + 1 (y1 ) = u
d
> dt + 2 (y2 ) = u
: y +y = y
1 2
68 Commande par espace d’état
d’où
1 1
y = y1 + y2 = d
:u d
:u
dt +1 + 2
! dt
1 1
= d d
:u
dt +1 dt + 2
Donc
d2 d
+3 +2 y =u
dt2 dt
ou encore
y• + 3y_ + 2y = u:
Ce raisonnement permet d’arriver au résultat escompté en peu de calcul et de façon élégante. Toutefois,
ce raisonnement demande à être placé dans un cadre mathématique. C’est ce qu’apporte la transformée de
Laplace.
Remarque. Le raisonnement précédent aurait en e¤et pu aboutir à une incohérence que l’on retrouve avec
aussi avec l’utilisation de la transformée de Laplace. Par exemple, le raisonnement suivant est clairement
faux
d
d d
y_ = u_ , y = u , y = dt d
u , y = u:
dt dt dt
d d
On a clairement pas le droit de simpli…er par dt car l’ensemble des opérateurs di¤érentiels générés par dt
est un anneau et non un corps. De même, dans le cadre des équations di¤érentielles non linéaires (comme
par exemple y• + y y_ = u) on arrive vite à des raisonnement complètement absurdes.
d
L’espace opérateurs di¤érentiels en dt est un anneau et admet de bonnes propriétés comme l’associativité
d4 d3 d d4 d3 d4 d d7 d5
+ = + = + ;
dt4 dt3 dt dt4 dt3 dt4 dt dt7 dt5
ou la commutativité,
d4 d3 d d3 d d4
+ = + :
dt4 dt3 dt dt 3 dt dt4
d d4 d3 d
On peut associer à l’opérateur dt le symbole s appelée variable de Laplace. Ainsi l’opérateur dt4 dt3
+ dt
sera représenté par le polynôme s4 s3 + s .
Systèmes linéaires 69
Soit un système linéaire d’entrée u et de sortie y reliés par une relation di¤érentielle du type
d
y(t) = H :u(t):
dt
La fonction H(s) est appelée fonction de transfert du système. Prenons par exemple le système décrit par
l’équation di¤érentielle
y• + 2y_ + 3y = 4u_ 5u:
On a !
d
4 dt 5
y(t) = d2 d
:u(t):
dt2
+ 2 dt +3
Sa fonction de transfert est donc
4s 5
H(s) = :
s2 + 2s + 3
Si la fonction de transfert H(s) d’un système est une fonction rationnelle, son dénominateur P (s) est appelé
polynôme caractéristique.
Soit une fonction y(t) analytique. Développons cette fonction autour du point t. On a
1 1
!
X 1 (i) i
X 1 d i d
y(t )= y (t): ( ) = (y(t)) = e dt (y(t)) :
i! i! dt
i=0 i=0
d
Ainsi, l’opérateur di¤érentiel e dt correspond à un retard de du signal y(t). La fonction de transfert du
système associé est e s .
On appelle transformée de Laplace y^(s) du signal y(t), la fonction de transfert du système générant y(t)
à partir de l’impulsion de Dirac (t). On notera y^(s) = L (y(t)). Ainsi, la transformée de (t ) est
s
P1 P1 s
e . De même, la transformée de Laplace d’une somme de Dirac i=0 i (t i ) est i=0 i e
I . En…n,
comme toute fonction intégrable f (t) peut s’approximer comme la somme d’une in…nité de Dirac in…niment
R1
rapprochés, c’est-à-dire, f (t) = 1 f ( ) (t )d , la transformée de Laplace de f (t) s’écrit
Z 1
f^(s) = f ( )e s
d :
1
4.5.6 Composition
On a
d d d
y(t) = H( ):u(t) = H( ):^
u : (t)
dt dt dt
Donc la transformée de Laplace de y(t) est
y^(s) = H(s):^
u(s):
Systèmes en série. Considérons deux systèmes de fonction de transfert H1 (s) et H2 (s) mis en série comme
sur la …gure.
La fonction de transfert du système résultant est H(s) = H2 (s):H1 (s). En e¤et, y^(s) = H2 (s):H1 (s)^
u(s).
Systèmes en parallèle. Mettons les deux systèmes de fonction de transfert H1 (s) et H2 (s) en parallèle
comme sur la …gure.
Système bouclé
Systèmes linéaires 71
4.5.7 Exemple 1
y^(s)
Calculons la fonction de transfert G(s) = u^(s) de ce système. Pour cela, notons, xi la sortie du ième
intégrateur. On a.
x
^1 x^1
y^ = 2^
x3 + 7^
x2 + 3^
x1 = 2
2+ 7 + 3^
x1 (4.13)
s s
^1 x
x ^1
s^
x1 = u
^(s) + 7^
x1 4^
x2 + x
^3 = u
^ + 7^
x1 4 + 2 :
s s
En isolant dans chacune de ces deux équations la variable x1 , il vient
y^(s) u
^(s)
x
^1 = = : (4.14)
2 s12 + 7 1s + 3 s 7 + 4 1s 1
s2
2 s12 + 7 1s + 3 3s2 + 7s 2
H(s) = = : (4.15)
s 7+ 4 1s 1
s2
s3 7s2 + 4s 1
4.5.8 Exemple 2
y^(s)
La fonction de transfert H(s) = u^(s) de ce système peut être calculée sous Scilab en tapant les lignes
suivantes.
s = poly(0,’s’)
H1 = 1/(s-1)
H2 = (s^2/(s+1))/(1+(s^2/(s+1))*(1/s)-s^2)
H3 = (s+1)/(1+(s+1)*(1/s))
G = G1*G2*G3
La première de ces lignes a pour but de créer le polynôme P (s) = s qui a pour unique racine 0. Les fonctions
H1 , H2 et H3 correspondent aux fonctions de transfert de chacun des trois blocs en série composant le
système. On obtient
s3 + s4
H(s) = :
2s5 + s4 6s3 s2 + 3s + 1
Le but de ce paragraphe est de montrer comment obtenir les équations di¤érentielles entrée-sortie (c’est-à-
dire que l’état n’intervient pas) d’un système linéaire donné par ses équations d’état :
(
x_ = Ax + Bu
y = Cx + Du:
Donc
1
^ = C(sI
y A) B+D u
^: (4.17)
La matrice
1
G(s) = C(sI A) B+D (4.18)
est appelée matrice de transfert. Il s’agit d’une matrice de fonctions de transfert (c’est-à-dire de fonctions
rationnelles en s) dont tous les dénominateurs sont des diviseurs du polynôme caractéristique PA (s) de
d
A. En multipliant de chaque côté par PA (s) et en remplaçant s par dt on obtient un système d’équations
di¤érentielles entrées-sorties pour notre système. L’état x n’y apparaît plus.
Exemple 4.6.1 On considère le système linéaire à temps continu décrit par sa représentation d’état :
8 ! !
>
> 1 3 1
>
> = x(t) + u(t);
< x(t)
_
2 0 1
! !
>
> 1 2 2
>
>
: y(t) = 1 0
x(t) +
0
u(t):
(4.18) 1B
G(s) = C(sI A) +D
! ! 1 ! !
1 2 s 1 3 1 2
= +
1 0 2 s 1 0 (4.19)
!
2s2 +s 7
= s2 s 6 :
s+3
s2 s 6
Ce calcul peut être e¤ ectué directement sous Scilab en tapant les lignes suivantes :
s=poly(0,’s’)
A=[13,20] ; B=[1,1] ; C=[12 ;10] ;
D=[2 ;0] ; I=eye(2,2)
G=C*(s*I-A)^(-1)*B+D
La première de ces lignes a pour but de créer le polynôme P (s) = s qui a pour unique racine 0. La fonction
eye(2; 2) génère une matrice identité de dimension 2 2. Remarquons que Scilab est ici capable d’e¤ ectuer
un peu de calcul symbolique (la variable s est ici un symbole). Mais les capacités de Scilab dans ce domaine
ne vont pas beaucoup plus loin. Ainsi, d’après (4.19), la relation y
^ = G(s)^ u s’écrit
!
2s2 + s 7
s2 s 6 y
^= u
^:
s+3
d
Soit, en substituant s par dt ; (
y•1 y_ 1 6y1 = 2•
u + u_ 7
y•2 y_ 2 6y2 = u_ + 3u:
74 Commande par espace d’état
Exemple 4.6.2 Considérons à nouveau le système linéaire formé deux masses reliées par un ressort intro-
duit au paragraphe 2.3.3 de la page 21. A partir de sa représentation d’état, la relation (4.17) nous dit que
sa matrice de transfert est
00 1 0 11 1 0 1
s 0 0 0 0 1 0 0 0
BB 0 s 0 0 C B k1 +k2 k2
0 C C B C
BB C B m1 m1 m1 CC B 0 C
1 0 0 0 BB C B CC B C:
@@ 0 0 s 0 A @ 0 0 0 1 AA @ 0 A
k2 k2 1
0 0 0 s m2 0 m2 m2 m2
La matrice de transfert est donc une fonction de transfert. Pour ce calcul, Scilab ne peut rien faire. En
revanche, il peut être e¤ ectué aisément par Mupad grâce aux instructions suivantes
M :=Dom : :Matrix() ;
A :=matrix([[0,1,0,0],[-(k1+k2)/m1,-a/m1,k2/m1,0],
[0,0,0,1],[k2/m2,0,-k2/m2,-a/m2]]) ;
B :=matrix([0,0,0,1/m2]) ;
C :=matrix([1,0,0,0]) ;C :=M : :transpose(C) ;
Id :=matrix([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]) ;
G :=C*(s*Id-A)^(-1)*B ;
Nous obtenons
1
2
:
m1 m2 4 m1 +m2 3 k1 +k2 2k2 +k1
k2 s + k2 s + m1 + k2 + k2 m2 s2 + k2 s + k1
L’équation di¤ érentielle entrée sortie est donc
m1 m2 d4 y m1 +m2 d3 y 2 k1 +k2
k2 dt4 + k2 dt3
+ m1 + k2 + k 2 m2 y•
+ 2k2k+k2
1
y_ + k1 y = u:
Dans ce paragraphe, nous montrons, dans le cas des systèmes linéaires monovariables (c’est-à-dire avec
une seule entrée et une seule sortie), comment obtenir une représentation d’état à partir d’une équation
di¤érentielle. Plusieurs méthodes sont proposées aboutissant chacune à une forme di¤érente.
On considère le système linéaire d’ordre 3 (bien que tout le raisonnement puisse être e¤ectué pour un ordre
quelconque) avec une seule entrée et une seule sortie, décrit par l’équation di¤érentielle suivante :
d3 y d2 y d2 u du
3
+ a2 2
+ a1 y
_ + a0 y = b 2 2
+ b1 + b0 u: (4.20)
dt dt dt dt
La fonction de transfert de notre système est
b2 s2 + b1 s + b0
G(s) = :
s3 + a2 s2 + a1 s + a0
Systèmes linéaires 75
C’est-à-dire, (
s3 x ^ a2 s2 x
^1 = u ^1 a1 s^ x1 a0 x ^1
2 (4.21)
y^ = b2 s x^1 + b1 s^
x1 + b0 x
^1 :
Cherchons à dessiner le câblage associé à ces deux équations. Le seul opérateur di¤érentiel que l’on s’autorise
est l’intégrateur, de fonction de transfert 1s : On construit d’abord une chaîne de 3 intégrateurs a…n de
fabriquer s3 x
^ 1 ; s2 x
^1 et s^
x1 . Puis, on câble les deux équations (4.21). On obtient le câblage de la …gure ci
dessous.
Les variables d’état de ce câblage sont les valeurs x1 ; x2 ; x3 mémorisées par chacun des intégrateurs (les
additionneurs et les ampli…cateurs ne mémorisent rien). En lisant le schéma, on peut directement écrire les
équations d’état de ce système
8 0 1 0 10 1 0 1
>
> x_ 1 0 1 0 x1 0
> B
> C B CB C B C
>
> @ x_ 2 A = @ 0 0 1 A @ x2 A + @ 0 A u
>
>
< x_ 3 a0 a1 a2 x 1
0 3 1 (4.22)
>
> x1
>
> B C
>
> y = b0 b1 b2 @ x2 A :
>
>
: x3
Rappelons que ce raisonnement peut s’appliquer a…n de trouver la représentation d’état pour tout système
linéaire monovariable d’ordre n quelconque. Cette forme particulière pour la représentation d’état, qui fait
76 Commande par espace d’état
apparaître dans les matrices les coe¢ cients de la fonction de transfert est appelée forme canonique de
commande.
Remarque 4.7.1 Pour obtenir la forme canonique de commande équivalente à un système linéaire mono-
variable donné, il nous faut calculer sa fonction de transfert (voir paragraphe 4.6) sous sa forme développée.
Ensuite, nous écrivons directement sa forme canonique de commande en nous in‡uençant par exemple de
la formule (4.22). La matrice de passage peut, elle-aussi, être obtenue aisément en résolvant un système
d’équations linéaires comme expliqué en remarque 4.4.1 page 66.
s3 y^ + a2 s2 y^ + a1 s^
y + a0 y^ = b2 s2 u
^ + b1 s^
u + b0 u
^:
8 0 1 0 1 0 1 0 1
>
> x_ 1 0 0 a 0 x 1 b 0
>
> B C B C B C B C
>
> @ x_ 2 A = @ 1 0 a1 A @ x2 A + @ b1 A u
>
>
< x_ 3 0 1 a2 x b2
0 31
>
> x1
>
> B C
>
> y = 0 0 1 @ x2 A :
>
>
: x3
Cette forme particulière pour la représentation d’état s’appelle forme canonique d’observation. Notons que
la transformation A ! AT ; B ! CT ; C ! BT , nous donne la forme canonique de commande (4.22).
Cherchons à expliquer maintenant pourquoi cette transformation, qui nous fait passer de la forme cano-
nique de commande à celle d’observation et vice-versa, ne change pas le système. Pour cela, considérons
le système dont les matrices d’état sont (A; B; C; D) et le système 0 dont les matrices d’état sont
A0 = AT ; B0 = CT ; C0 = BT ; D0 = DT : La matrice de transfert associée à est
(4.18) 1
G(s) = C(sI A) B + D:
(4.18) 1
G0 (s) = C0 (sI A0 ) 1 B0 + D0 = BT sI AT CT + DT
T T
= BT (sI A) 1 CT + DT = C(sI A) 1 B + DT
T
= C(sI A) 1 B + D = GT (s):
Un système linéaire monovariable est sous forme modale s’il s’écrit sous la forme :
8 0 1 0 1
>
> 0 0 1
>
> 1
B C
>
> B 0 0 C B 1 C
>
< x_ = B 2 C
B Cx + B C
B .. C u
@ A @ . A (4.23)
>
>
>
> 0 0 n 1
>
>
>
: y = c1 c2 cn x + d u:
Son câblage est donné par la …gure ci-dessus. La fonction de transfert du système est donnée par
1
G(s) = C(sI A) B+d
0
1
0 1 1
1
s 1 0 ::: B C
B .. C B 1 C
= c1 cn @ 0 . 0 A B C+d
@ A
::: 0 s n 1
0 10 1
1
s 1
0 0 1
B CB
B 0 s
1
0 CB 1 C C
= c1 cn B 2 CB C+d
B .. C@ A
@ . A
0 0 1 1
s n
c1 c2 cn
= + + + + d:
s 1 s 2 s n
Systèmes linéaires 79
0 1
s 1 0 0
B 0 s 0 C
B 2 C
det(sI A) = det B C = (s 1 ) (s 2 ) : : : (s n)
@ A
0 0 s n
Remarque 4.7.2 Si nous disposons d’un système linéaire monovariable donné par sa représentation d’état
et si nous souhaitons obtenir sa forme modale, il nous faudra calculer sa fonction de transfert (voir pa-
ragraphe 4.6) et la décomposer en éléments simples. Si tous ces éléments simples sont de première espèce
(c’est-à-dire que la fonction de transfert a la forme s c1 1 + s c2 2 + + s cn n + d) il est possible de mettre
le système sous forme modale (attention, fréquemment, cela n’est pas possible). On peut écrire directement
cette forme modale en s’in‡uençant par exemple de la formule (4.23). La matrice de passage P peut être
obtenue aisément en résolvant un système linéaire comme expliqué en remarque 4.4.1. Notons en…n que cette
matrice P a pour colonnes les vecteurs propres de la matrice d’évolution A.
0 1 0 1
2 1 0 0 0 0
B 0 2 1 0 0 C B 0 C
B C B C
B C B C
x_ = B 0 0 2 0 0 Cx + B 1 Cu
B C B C
@ 0 0 0 3 1 A @ 0 A
0 0 0 0 3 1
y = 2 1 3 4 7 x + 2u
est sous forme de Jordan car sa matrice d’évolution A est une matrice de Jordan. C’est-à-dire qu’elle est
diagonale par bloc et que chaque bloc (ici il y en a deux) possède des zéros partout, sauf sur sa diagonale,
qui contient des éléments tous égaux et sur sa sur-diagonale qui ne contient que des uns. De plus, la matrice
de commande ne possède, pour élément non nuls, que des uns (autant que de blocs), positionnés au niveau
de la dernière ligne de chaque bloc. Un câblage pour ce système est donné sur la …gure ci-dessous.
80 Commande par espace d’état
On cherche à obtenir une représentation plus simple, c’est-à-dire, avec plus de zéros et plus de uns, a…n de
limiter de nombre d’ampli-op nécéssaires à la réalisation du circuit. La décomposition sous forme de jordan
de la matrice d’évolution est
0 1 10 10 1 0 1
1 1
1 1 1 4 2 2 1 1 1 2 1 0
B C B CB C B C
@ 2 1 2 A @ 4 1 1 A@ 2 1 2 A = @ 0 2 0 A:
2 1 0 4 2 2 2 1 0 0 0 3
| {z }| {z }| {z } | {z }
P 1 A P A
5.1 Introduction
Dans ce chapitre, nous allons étudier la conception de régulateurs pour les systèmes donnés par des équations
d’état linéaires. Nous montrerons dans le chapitre suivant qu’autour de points bien particuliers de l’espace
d’état, dit points de fonctionnement, de nombreux systèmes non linéaires se comportent comme de véritables
systèmes linéaires. Les techniques développées dans ce chapitre seront alors utilisées pour la commande de
systèmes non linéaires.
Notons m; n; p les dimensions respectives des vecteurs u, x et y. Rappelons que A est appelée matrice
d’évolution, B est la matrice de commande et C est la matrice d’observation.
Remarque 5.1.1 Nous avons ici supposé, dans un but de simpli…cation que la matrice directe D intervenant
dans l’équation d’observation (voir 1.3 page 14) était nulle. C’est souvent le cas en pratique. Si une telle
matrice directe existe, on peut facilement s’en débarrasser. En e¤ et, considérons le système
(
x_ = Ax + Bu
y = Cx + Du;
et réalisons la nouvelle sortie z = y Du comme indiqué sur la …gure. On obtient un nouveau système dont
les équations d’état sont
(
x_ = Ax + Bu
z = Cx:
81
82 Commande par espace d’état
Après avoir dé…ni les notions fondamentales de commandabilité et d’observabilité dans les paragraphes
5.2.2 et 5.2.3, nous proposerons deux approches pour la conception de régulateurs. Tout d’abord, dans le
paragraphe 5.3, nous allons supposer que l’état x est accessible à la mesure. Bien que cette hypothèse ne soit
généralement pas véri…ée, elle nous permettra d’établir les principes de la méthode par placement de pôles.
Dans une deuxième phase, nous ne supposerons plus que l’état est accessible. Il nous faudra alors développer
des estimateurs d’état, au paragraphe 5.4, capables d’estimer le vecteur état a…n de pouvoir utiliser les outils
développés lors de la première phase.
Le lecteur pourra consulter l’ouvrage de Philippe de Larminat [4] pour avoir une vision large de l’ensemble des
méthodes utilisées pour la commande des systèmes linéaires. Un cours complet et pédagogique accompagné
de nombreux exercices peut être trouvé dans les livres de Maurice Rivoire et Jean-Louis Ferrier [6].
Ce petit rappel nous sera utile pour comprendre certains éléments de ce paragraphe comme par exemple les
preuves des critères de commandabilité et d’observabilité. Le système linéaire
Ax = b
5.2.2 Commandabilité
Il existe de multiples dé…nitions équivalentes pour la commandabilité des systèmes linéaires. Une dé…nition
simple pour la commandabilité est la suivante.
Dé…nition 5.2.1 Le système linéaire (5.1) est dit commandable si pour tout couple de vecteur d’état
Commande des systèmes linéaires 83
(x0 ; x1 ), on peut trouver un temps t1 et une commande u(t), t 2 [0; t1 ], tel que le système, initialisé en
x0 ; atteigne l’état x1 ; à l’instant t1 .
Théorème 5.2.1 (critère de commandabilité) Le système linéaire 5.1 est commandable si et seulement si
rang B j AB j A2 B j : : : j An 1
B = n;
c’est-à-dire que la matrice com , dite de commandabilité, obtenue en juxtaposant les unes à côté des autres
les n matrices B, AB, : : : ; An 1 B; est de rang plein.
La dé…nition et le théorème sont aussi valables pour les systèmes linéaires à temps discret.
Le rang de com est égal à 2 (les deux dernières lignes sont identiques) : le système est donc non comman-
dable. En revanche, pour une matrice de commande donnée par
0 1
0 0
B C
B = @ 1 0 A;
1 1
Le rang de com est égal à trois car on peut engendrer les trois vecteurs e1 = (1 0 0)T ; e2 = (0 1 0)T ; e3 =
(0 0 1)T par combinaisons linéaires des vecteurs colonnes vi de C1 . En e¤ et, e1 = v3 v2 , e2 = v1 v2 et
e3 = v2 . Ainsi, le système est commandable.
Preuve du critère de commandabilité. Nous allons nous limiter à donner une preuve pour les systèmes
linéaires à temps discret x(k + 1) = Ax(k) + Bu(k), la démonstration dans le cas continu n’étant pas aussi
simple. Pour cela, nous allons montrer que l’on peut imposer arbitrairement x(k + n) en jouant sur les
entrées u(k); : : : ; u(k + n 1) et ceci quel que soit x(k). Si x(k + n) peut être choisi comme on le désire, il
84 Commande par espace d’état
semble clair qu’il est possible de se diriger dans l’espace d’état, et donc que le système est commandable.
Nous avons
Ainsi,
0 1
u(k + n 1)
B u(k + n 2) C
B C
B .. C
x(k + n) = An x(k) + B j AB j : : : j An 1
B B . C
| {z }B
B
C
C
com @ u(k + 1) A
u(k)
= An x(k) + com v;
5.2.3 Observabilité
Une dé…nition simple pour l’observabilité d’un système linéaire est la suivante.
Dé…nition 5.2.2 Le système linéaire 5.1 est dit observable si la connaissance de y(t) et de u(t) pour t 2 R
nous permet de déterminer de façon unique l’état x(t), pour tout t.
c’est-à-dire que la matrice, dite d’observabilité obs , obtenue en mettant les une en dessous des autres les n
matrices C, CA, : : : , CAn 1 ; est de rang plein.
Commande des systèmes linéaires 85
Ce théorème est aussi valable pour les systèmes linéaires à temps discret.
Preuve du critère d’observabilité. Pour donner une idée de la preuve de ce critère, nous allons montrer
que la connaissance des n 1 premières dérivées des sorties et des n 2 dérivées des entrées nous permet
de retrouver le vecteur d’état: Pour cela, dérivons n 1 fois l’équation d’observation. Nous obtenons
y = Cx
y_ = CAx + CBu
• = CA2 x + CABu + CBu_
y
...
y(n 1)
= CAn 1
x + CAn 2
Bu + CAn 3
Bu_ + + CABu(n 3)
+ CBu(n 2)
:
z= obs x + v;
où z est le vecteur de toutes les sorties et de leurs dérivées, v le vecteur de toutes les entrées et leurs dérivées,
obs la matrice d’observabilité et la matrice restante. Le système qu’il nous faut résoudre pour retrouver
l’état x est donné par
obs x = z v:
Cette équation admet au plus une solution si obs est de rang plein. L’absence de solution signi…erait que
v et z sont incompatibles avec les équations de notre système, ce qui est incompatible avec nos hypothèses.
Cette solution est donnée par
1 T
x= T obs obs obs : (z v) :
1
La matrice T obs obs
T
obs s’appelle l’inverse généralisée de la matrice obs . Elle n’existe que si obs est
de rang plein.
La décomposition de Kalman va nous permettre de mieux comprendre comment interviennent les problèmes
de commandabilité et d’observabilité dans les systèmes linéaires.
Notons que sur la …gure, il n’existe aucun chemin (en respectant le sens des ‡èches) qui mène de l’entrée u
vers un système non commandable. De même, il n’existe aucun chemin qui mène d’un système non observable
vers y. Le schéma de la …gure ci-dessous montre, sous une forme développée la décomposition en 4 sous-
systèmes Si ; chacun étant associé au sous-vecteur d’état xi . D’après ce schéma, nous pouvons en déduire
qu’un système linéaire à temps continu peut toujours se mettre (après un changement de base adéquat) sous
la forme suivante :
8 0 1 0 1
>
> A11 A12 0 0 B1
>
> B C B C
>
> B 0 A22 0 0 C B 0 C
< x(t)
_ = B C x(t) + B C u(t)
@ A31 A32 A33 A34 A @ B3 A
>
>
>
> 0 A42 0 A44 0
>
>
: y(t) = C1 C2 0 0 x(t) + (D) u(t)
Il est légitime de vouloir choisir la matrice de régulation K de façon à imposer les pôles du système bouclé.
Ce problème est équivalent à imposer le polynôme caractéristique du système. Soit Pcom (s) le polynôme
désiré, que l’on supposera bien sûr de degré n. Il nous faut résoudre l’équation polynomiale
dite de placement de pôles. Cette équation peut se traduire en n équations scalaires. Rappelons en e¤et
que deux polynômes de degré n et unitaires sn + an 1 s + + a0 et sn + bn 1 s + + b0 sont égaux si et
seulement si leurs coe¢ cients sont tous égaux, c’est-à-dire si an 1 = bn 1 ; : : : ; a0 = b0 : Notre système de
n équations possède m:n inconnues qui sont les coe¢ cients kij , i 2 f1; : : : ; mg; j 2 f1; : : : ; ng. En fait, une
seule matrice solution K nous su¢ t. On peut donc …xer (m 1) éléments de K a…n qu’il ne nous reste plus
que n inconnues. Mais le système obtenu n’est pas toujours linéaire. L’instruction ppol de Scilab permet
de résoudre l’équation (5.2).
Tout est simple lorsque le système possède une seule entrée. En e¤et, l’équation polynomiale (5.2) se traduit
forcément par un système de n équations linéaires à n inconnues qui admet une et une seule solution (car le
système est commandable).
Nous allons ici illustrer la résolution de l’équation polynomiale (5.2), lorsque le système n’admet qu’une
seule entrée. Les méthodes proposées ici nécessitent des calculs assez fastidieux. Un logiciel de calcul formel
comme Mupad peut s’avérer bien utile. Considérons par exemple le système
0 1 0 1
1 4 1 2
B C B C
x_ = @ 6 1 3 A x + @ 3 A u;
2 2 5 1
que l’on cherche à stabiliser par un retour d’état de la forme
u=w Kx;
avec
K= k1 k2 k3 :
Cherchons K de façon à ce que ce polynôme caractéristique Pcom (s) du système en boucle fermée ait pour
racines 1; 1 2i; 1 + 2i, c’est-à-dire,
c’est-à-dire 0 1
s 1 + 2k1 4 + 2k2 1 + 2k3
B C
det @ 6 + 3k1 s + 1 + 3k2 3 + 3k3 A = s3 + 3s2 + 7s + 5
2 k1 2 k2 s + 5 k3
ou encore
s3 + (2k1 + 3k2 k3 + 5) s2 + (25k1 + 21k2 + 10k3 29) s
+41k1 + 72k2 + 71k3 129 = s3 + 3s2 + 7s + 5:
On obtient le système linéaire suivant
0 10 1 0 1 0 1
2 3 1 k1 5 3
B CB C B C B C
@ 25 21 10 A @ k2 A + @ 29 A = @ 7 A :
41 72 71 k3 129 5
Ainsi,
0 1 0 1 1 00 1 0 11 0 1
k1 2 3 1 3 5 1:4227
B C B C BB C B CC B C
k =
@ 2 A @ 25 21 10 A @@ 7 A @ 29 AA = @ 0:94158 A :
k3 41 72 71 5 129 2:0206
Soit
K= 1:4227 0:94158 2:0206 :
Nous aurions obtenir directement ce résultat en utilisant l’instruction ppol (pole placement) de Scilab en
tapant
A = [1; 4; 1; 6; 1; 3; 2; 2; 5];
B = [2; 3; 1];
K = ppol(A; B; [ 1; 1 2 %i; 1 + 2 %i]);
5.4.1 Principe
mais cette fois-ci, l’état x du système n’est plus supposé accessible à la mesure. Seules la consigne w et la
sortie du système y peuvent être utilisées par le régulateur. Nous allons chercher à estimer l’état a…n de
pouvoir calculer la commande u. Pour cela, on intègre à notre régulateur un simulateur de notre système.
L’erreur "y entre la sortie du simulateur y ^ et la sortie du système y nous permet de venir corriger, par
l’intermédiaire d’une matrice de correction L, l’évolution de l’état estimé x
^. Le simulateur corrigé s’appelle
observateur . Son seul rôle est de nous donner une bonne estimation du vecteur d’état x(t) a…n que l’on puisse
appliquer une technique de retour d’état. Le schéma de principe du régulateur par retour de sortie est donné
sur la …gure. Les inconnues du régulateur sont les matrices K et L. Pour le calcul K nous pouvons utiliser la
méthode par placement de pôles décrite précédemment, qui consiste à résoudre det(sI A+BK) = Pcom (s),
où Pcom (s) est le polynôme caractéristique de degré n choisi pour la dynamique de commande. Il nous reste
maintenant à trouver la matrice de correction L.
A…n de calculer L, extrayons du système régulé de la …gure, le sous-système d’entrée u formé du système à
réguler et de son observateur (voir …gure ci-dessous).
Principe de l’observateur
d
(^
x x) = A^
x + Bu L(C^
x Cx) Ax Bu = A (^
x x) LC(^
x x):
dt
Ainsi, les nouvelles équations d’état s’écrivent
(
x_ = Ax + Bu
"_ x = A"x LC"x
Les pôles de ce système correspondent aux valeurs propres de A et à celles de A LC. La structure de
l’équation (5.4) nous montre que "x est un sous-vecteur d’état non commandable. Il obéit à l’équation
di¤érentielle
"_ x = (A LC) "x (5.5)
sur laquelle la commande u n’intervient pas. L’erreur d’estimation "x sur l’état tends vers zéro si toutes les
valeurs propres de A LC sont à parties réelles négatives. Imposer la dynamique de l’erreur (c’est-à-dire,
sa rapidité de convergence) revient à résoudre
où Pobs (s) est choisi comme on le désire, de façon à avoir les pôles requis. Puisque le déterminant d’une
matrice est égal à celui de sa transposée, cette équation est équivalente à
On obtient une équation de type placement de pôles (voir équation (5.2) page 87). La même méthode que
celle exposée à la section 5.3 page 87 peut donc s’appliquer pour trouver LT et donc L: On pourra donc
utiliser la fonction ppol de Scilab.
Ce sont les équations qu’il nous faut câbler ou programmer pour réguler notre système. La méthode pour
fabriquer un régulateur par retour de sortie est récapitulée sur le tableau 5.1, où pcom et pobs sont les
vecteurs des pôles désirés pour le régulateur et pour l’observateur.
Commande des systèmes linéaires 93
et cherchons un régulateur par retour de sortie qui place tous les pôles en 1. Pour trouver K et L il nous
faut résoudre
det (s 3 + 2K) = s + 1
det (s 3 + L:4) = s + 1
On obtient donc K = 2 et L = 1. Le régulateur recherché est donc donné par
(
d
R := dt x
^ = 5^
x + 2w + y
u = 2^
x+w
On peut alors véri…er que les pôles de ce système sont bien ceux qui nous avions placés (à savoir 1 et 1).
Ce phénomène est expliqué dans le paragraphe suivant.
Les pôles pcom placés pour la conception du régulateur par retour d’état et ceux pobs placés pour l’observateur
sont-ils bien les pôles du système en boucle fermée ? Le but de ce paragraphe est de répondre (positivement)
à cette question.
Les équations d’évolution associées au système bouclé d’entrée w et de sortie y sont données par
(
x_ = Ax + B (w K^ x)
d^
x (5.6)
dt = (A BK LC) x ^ + LCx + Bw:
94 Commande par espace d’état
Insistons sur le fait que le système bouclé possède pour vecteur d’état ceux du système, rangés dans x, et
ceux du régulateur rangés dans x ^. Son vecteur d’état est donc donné par (xT x
^T )T . Sous forme matricielle,
(5.6) s’écrit ! ! ! !
d x A BK x B
= + w
dt x
^ LC A BK LC x
^ B
Posons "x = x
^ x. Puisque ! ! !
x I 0 x
=
"x I I x
^
ou de façon équivalente ! ! !
x I 0 x
= ;
x
^ I I "x
un autre vecteur d’état possible pour le système bouclé est (x; "x ). Les équations d’état deviennent après
changement de base
! ! ! ! ! ! !
x_ I 0 A BK I 0 x I 0 B
= + w
"_ x I I LC A BK LC I I "x I I B
ou encore ! ! ! !
x_ A BK BK x B
= + w (5.7)
"_ x 0 A LC "x 0
Notons que l’entrée w ne peut pas agir sur "x , ce qui est compatible avec le fait que "x n’est toujours pas
(voir équation 5.5) un sous-vecteur commandable (le fait de boucler le système ne peut, bien sûr, pas rendre
commandable une variable d’état non commandable).
Donc les pôles du système bouclé sont constitués des pôles placés pour la commande et des pôles placés
pour l’observation. C’est le principe de séparation.
L’algorithme présenté sur le tableau 5.1 page 93 pour synthétiser un régulateur nous assure que si le vecteur
de consigne w est nul, l’état du système x va converger vers 0 avec une dynamique déterminée par les
Commande des systèmes linéaires 95
pôles placés. Lorsque w n’est plus nul, l’état converge vers une valeur qui n’est plus forcément nulle. Un
précompensateur est une matrice carrée H; que l’on place juste après le vecteur de consigne, comme sur la
…gure.
Ce précompensateur ne change pas les pôles du système bouclé. Il permet de mettre en correspondance
certaines composantes de la consigne avec certaines variables d’état préalablement choisies. Les équations
d’évolution du système bouclé (voir équation (5.7) où w a été remplacé par Hw) sont
! ! ! !
d x A BK BK x BH
dt = + w:
"x 0 A LC "x 0
Si w(t) est une constante w, une fois atteint le point d’équilibre, nous avons
! ! !
A BK BK x BH
0 = + w:
0 A LC "x 0
Puisque (A LC) est inversible (car tous les pôles placés pour l’observateur sont strictement stables), "x
est nécessairement nul. L’équation précédente devient
(A BK) x + BHw = 0:
Puisque A BK est aussi inversible (car tous les pôles placés pour la commande sont strictement stables),
nous avons
1
x= (A BK) BHw: (5.8)
On appelle variables consignées, xc un ensemble de m variables d’état (où m = dim(w) = dim (u)) pour
lesquelles on souhaiterait que, pour w = w constant, xc converge vers xc = w. Supposons que ces variables
puissent être obtenues par une combinaison linéaire des composantes de x par une relation du type
def
xc = Ex; (5.9)
96 Commande par espace d’état
Ainsi
1
1 1
xc = w , E (A BK) BH = I , H = E (A BK) B :
L’insertion d’un précompensateur permet donc d’assigner à chacune des consignes, constituant w, une
variable d’état particulière. Les variables ainsi consignées peuvent être alors commandées indépendamment
les unes des autres. L’algorithme du tableau 5.2 récapitule la méthode permettant de calculer un régulateur
par retour de sortie avec précompensateur. La fonction Scilab associée est donnée par le tableau 5.3.
Tableau 5.2 –Algorithme de calcul d’un régulateur par retour de sortie avec précompensateur
function [Ar,Br,Cr,Dr]=RegulKLH(A,B,C,E,pcom,pobs) ;
K=ppol(A,B,pcom)
L=ppol(A’,C’,pobs)’
H=-inv(E*inv(A-B*K)*B)
Ar=A-B*K-L*C
Br=[B*H L]
Cr=-K
Dr=[H,zeros(C*B)]
endfunction ;
Tableau 5.3 –Fonction Scilab permettant de construire un régulateur pour un système linéaire
Chapitre 6
6.1 Introduction
Au chapitre 5, nous avons montré comment concevoir des régulateurs pour les systèmes linéaires. Or, en
pratique, les systèmes sont rarement linéaires. En revanche si leur vecteur d’état reste localisé dans une
zone de petite taille de l’espace d’état, le système peut être assimilé à un système linéaire et les techniques
développées au chapitre 5 peuvent alors être utilisées. Dans le paragraphe 6.2, nous allons tout d’abord
montrer comment linéariser un système non linéaire autour d’un point donné de l’espace d’état. Ensuite, au
paragraphe 6.3, nous montrerons comment stabiliser les systèmes non linéaires par les régulateurs proposés
au chapitre précédent. Quelques exemples illustratifs seront alors traités en détail aux paragraphes 6.4, 6.5
et 6.6.
df
f (x) ' f (x) + (x) (x x) ; (6.1)
dx
avec 0 1
@f1 @f1 @f1
@x1 (x) @x2 (x) : : : @xn (x)
B @f2 @f2 @f2 C
df B @x1 (x) @x2 (x) : : : @xn (x) C
(x) = B
B .. .. .. C:
C
dx @ . . . A
@fp @fp @fp
@x1 (x) @x2 (x) : : : @xn (x)
97
98 Commande par espace d’état
Cette matrice est appelée matrice jacobienne. Considérons par exemple la fonction
! !
x1 x21 x2
f = :
x2 x21 + x22
Le calcul de matrice jacobienne est facilité par le calcul formel. Mupad peut ainsi d’obtenir cette matrice
par les instructions suivantes :
export(linalg)
J :=jacobian([x1^2*x2,x1^2+x2^2],[x1,x2])
J0 :=subs(J(x1,x2),x1=1,x2=2)
La fonction h (x; u) sera appelée la fonction d’évolution/observation. Autour du point z = (x; u), nous avons
dh
h(z) ' h(z) + (z) (z z) ;
dz
avec 0 1
@h1 @h1
@z1 (z) @z2 (z) : : :
dh B @h2
(z) @h2
(z) : : : C
(z) = @ @z1 @z2 A:
dz .. ..
. .
Commande linéaire des systèmes non linéaires 99
Autour du point (x; u); le comportement de S s’approxime donc par les équations d’état suivantes :
(
x_ = f (x; u) + A (x x) + B (u u)
(6.2)
y = g(x; u) + C (x x) + D (u u) :
Il s’agit d’un système a¢ ne (car sa fonction d’évolution/observation est a¢ ne) qui est appelé système tangent
à S au point (x; u):
Un point (x; u) est un point de fonctionnement (aussi appelé point de polarisation) si f (x; u) = 0. Si u = 0,
on parle de point d’équilibre. Remarquons tout d’abord que si x = x et si u = u, alors x_ = 0, c’est-à-dire
que le système n’évolue plus si on maintient la commande u = u et s’il est dans l’état x. Dans ce cas, la
sortie y a pour valeur y = y = g(x; u). Autour du point de fonctionnement (x; u), d’après (6.2), le système
S admet pour système tangent :
(
x_ = A (x x) + B (u u)
y = y + C (x x) + D (u u) :
Posons ue = u u, x e = x x et y e=y y. Ces vecteurs sont appelés les variations de u, x et y. Pour des
e, x
petites variations u e, y
e ; on a (
d
e
dt x = Ae
x + Be
u
e
y = Ce
x + De
u:
Le système ainsi formé est appelé système linéarisé de S autour du point de fonctionnement (x; u).
100 Commande par espace d’état
est un point d’équilibre. Autour de ce dernier, f (x) peut s’approximer par son système tangent comme suit
df
f (x) ' f (x) + (x) (x x)
! dx ! !
0 1 x2 x1 x1 x1
= +
0 x2 1 + x1 x2 x2
! ! ! !
0 x1 1 x2 + 1 1
2 x2 + 50
= = = :
0 x2 1 x1 1 2x1 100
Comme le montre la …gure ci-dessous, l’approximation reste valable seulement autour du point autour duquel
nous avons linéarisé.
e1 = x1
Le système linéarisé s’obtient en posant x e2 = x2
x1 et x x2 . Il est donné par
!
d 0
e=
x e:
x
dt 0
Commande linéaire des systèmes non linéaires 101
Les valeurs propres s’obtiennent en calculant les racines du polynôme caractéristique. Puisque
!
s
det = s2 + 1;
s
les valeurs propres sont i. Elles correspondent à un système oscillant. Sous Mupad le calcul précédent se
traduit par :
export(linalg)
J :=jacobian([(1-a*x2)*x1,(-1+b*x1)*x2],[x1,x2])
J0 :=subs(J(x1,x2),x1=1/b,x2=1/a)
charpoly(J0,s)
où charpoly(J0,s) calcule le polynôme caractéristique de la matrice J0.
On rappelle (voir page 21) que les équations d’état d’un pendule simple sont données par
8 ! !
>
< x_ 1 x2
= f (x; u) = `mg sin x1 +u
x_ 2 m`2
>
:
y = g(x; u) = ` sin x1 :
Puisque la matrice jacobienne de h(x; u), la fonction évolution/observation, est
0 1
0 1 0
dh B 1 C
(x; u) = @ g cos `
x1
0 m` 2 A;
d(x;u)
` cos x1 0 0
au point (x = (0; 0);u = 0), le système linéarisé est décrit par les équations d’état suivantes :
8 ! !
>
> 0 1 0
< x_ = x+ u
g 1
` 0 m`2
>
>
: y = ` 0 x:
Sous Mupad, ce calcul peut être e¤ectué grâce aux instructions suivantes :
export(linalg)
J :=jacobian([x2,(-l*m*g*sin(x1)+u)/(m*l^2),l*sin(x1)]
,[x1,x2,u])
J0 :=subs(J(x1,x2,u),x1=0,x2=0,u=0)
simplify(J0)
On considère à nouveau le vérin modélisé au paragraphe 2.5.3 de la page 35. Ses équations d’état sont
données par 8
>
< x_ 1 = x2
ax3 kx1
S: x_ 2 = m
>
: x_ = x3 u
3 x1 x 2 a :
102 Commande par espace d’état
Supposons que x1 > 0 et calculons l’ensemble des points de fonctionnement (x; u) possibles pour S. La
condition x_ = 0 donne 8
>
< x2 = 0
ax3 kx1 = 0
>
: x x u
1 2 a = 0
soit 8
>
< x2 = 0
ax3 kx1 = 0
>
: u = 0
Les points de fonctionnement sont donc de la forme
k
(x; u) = x1 ; 0; x1 ; 0 :
a
Di¤érencions maintenant S autour du point de fonctionnement (x; u) a…n d’en obtenir une approximation
a¢ ne. Sous Mupad , on tape
export(linalg)
J :=jacobian([x2,(a*x3-k*x1)/m,-(x3/x1)*(x2-u/a)]
,[x1,x2,x3,u])
J0 :=subs(J(x1,x2,x3,u),x1=x1bar,x2=0,x3=k*x1bar/a,u=0)
simplify(J0)
on obtient 0 1
0 1 0 0
df B k a C
(x; u) = @ m 0 m 0 A:
d(x;u) k k
0 a 0 a2
Le système linéarisé s’exprime donc par
0 1 01
0 1 0 0
B k a C B C
x_ = @ m 0 m A (x x) + @ 0 A u:
k k
0 a 0 a2
On remarque que les coe¢ cients des matrices A et B ne dépendent pas de x: Ceci est assez rare, et signi…e
que le comportement d’une commande linéaire ne dépendra pas du point de fonctionnement choisi. Notons
que l’absence de matrice C et D dans le système linéarisé est une conséquence du fait que le système non
linéaire considéré soit autonome (c’est-à-dire sans sortie).
6.3.1 Principe
dont (x; u) constitue un point de fonctionnement. A…n que notre système se comporte comme un système
e = u u et y
linéaire autour de (x; u), construisons les variables u e = y y, comme représenté sur la …gure
ci-dessous.
e et de sortie y
Le système d’entrée u e ainsi réalisé est appelé système polarisé. Les équations d’état du système
polarisé peuvent s’approximer par
(
x_ = f (x; u) ' A: (x x) + B: (u u) = A: (x x) + B:e u
e =
y y + y = y + g(x) ' y + C (x x)
qui est en fait le système linéarisé du système non linéaire (6.3). Soit xc = Ex le sous-vecteur des variables
consignées. Nous pouvons construire un régulateur RL pour ce système linéaire grâce à l’algorithme Re-
gulKLH(A; B; C; E; pcom ; pobs ) de la page 96. Nous savons que, lorsque la consigne w e entrée dans RL
est constante, nous avons Ee x = w.e Or, nous voudrions que l’entrée du régulateur w que nous fabriquons
satisfasse w = Ex: Nous devons donc construire w e à partir de w de telle façon que w = Ex, à l’équilibre.
On a
we = Ee
x = E (x x) = w w;
où w = Ex. Le régulateur ainsi obtenu est représenté sur la …gure ci-dessous, dans le cadre épais.
Ce régulateur stabilise et découple notre système non linéaire autour de son point de fonctionnement. Un
récapitulatif de la méthode pour calculer un régulateur pour un système non linéaire est donné par le tableau
6.1. Cet algorithme nous renvoie notre régulateur d’entrées y; w et de sortie u sous la forme de ses équations
d’état.
Tableau 6.1 – Méthodologie à suivre pour la conception d’un régulateur linéaire pour un système non-
linéaire
Exemple 6.3.1 Considérons le système non linéaire donné par les équations d’état
(
x_ = 2x2 + u
S:
y = 3x;
que l’on cherche à stabiliser autour de l’état x = 2. L’intérêt de cet exemple est que les vecteurs d’état,
d’entrée et de sortie sont tous des scalaires. Ainsi, toutes les matrices sont aussi scalaires et les calculs
peuvent se faire à la main. A…n d’avoir f (x; u) = 0, il nous faut 2x2 + u = 0, soit u = 8. Si nous voulons
qu’à l’équilibre y = w, il faudra prendre E = 3. Si de plus nous voulons que tous les pôles du système bouclé
soient égaux à 1; il faudra pcom = pobs = 1. La linéarisation autour de (x; u) nous donne A = 8; B = 1
et C = 3: Pour K et L; il faut résoudre les deux équations polynomiales
(
det (sI A + BK) = s+1
det sI AT + C T LT = s + 1
Rappelons (voir paragraphe 5.4 de la page 90) que dans le régulateur renvoyé par l’algorithme RegulNL
du tableau 6.1, se cache un simulateur du système linéarisé de notre système non linéaire S à réguler. Il
Commande linéaire des systèmes non linéaires 105
est possible de remplacer ce simulateur par un simulateur non linéaire ayant les mêmes équations que S
lui-même. Le régulateur généré par l’algorithme RegulNL, de la page 104, peut s’écrire :
(
d
^ = A^
dt x x + B ( K^ x + H (w w)) L (C^ x + y) + Ly
R:
u = u K^ x + H (w w) :
Or,
A^x + B ( K^ x + H (w w)) // voir équation d’évolution de R
= A^ x + B (u u) // équation d’observation de R
= f (x; u) + A^x + B (u u) // puisque f (x; u) = 0
' f (^
x + x; u) // Taylor d’ordre 1 en (x; u).
= f (^
x + x; u K^ x + H (w w)) // équation d’observation de R
et
C^x+y // voir équation d’évolution de R
= g(x) + C^x // puisque y = g (x)
' g(^
x + x) // Taylor à l’ordre 1 autour du point x.
Nous pouvons donc remplacer le régulateur R par le suivant
(
d
^ = f (^
dt x x + x; u K^ x + H (w w)) L (g(^
x + x) + y)
Re :=
u = u K^ x + H (w w) ;
qui intègre non pas un simulateur du système linéarisé, mais un simulateur du système non linéaire (les
matrices A; B; C n’apparaissent plus dans ce régulateur). L’observateur du régulateur ainsi obtenu est plus
proche de la réalité et se trompe donc moins. On obtient ainsi un régulateur plus performant.
On considère à nouveau le système modélisé au paragraphe 2.3.5 de la page 23. Rappelons que l’entrée u est
la force exercée sur le chariot de masse M , x est la position du chariot et est l’angle entre le pendule et la
verticale. En prenant pour vecteur d’état x = x; ; x; _ _ et en supposant que seule la position du chariot x
est mesurée, les équations d’état sont données par :
8 0 1 0 1
>
> x_ 1 x3
>
> B C
> B x_ C x4
< B 2 C = B
>
B
C
C
B C B m sin x2 (g cos x2 `x24 )+u C
@ x_ 3 A @ M +m sin 2x
A
>
> 2
2
>
> x_ 4
sin x2 ((M +m)g m`x4 cos x2 )+cos x2 u
>
> `(M +m sin2 x2 )
:
y = x1
(voir (2.3) page 24). A…n de linéariser ce système, il nous faut calculer la matrice jacobienne de la fonction
évolution/observation h(x; u) puis l’évaluer au point de fonctionnement x = (0; 0; 0; 0) et u = 0. Nous
obtenons la matrice 0 1
0 0 1 0 0
B 0 0 0 1 0 C
dh B C
B mg 1 C
(x;u) = B 0 M 0 0 M C
d (x;u) B C
@ 0 (MM +m)g
` 0 0 1 A
M`
1 0 0 0 0
106 Commande par espace d’état
Ici, comme x = 0, les quantités xe et x se confondent et c’est pour cela que nous avons gardé la notation x
dans l’équation d’état ci-dessus. Tous ces calculs peuvent être obtenus par Mupad grâce aux instructions
suivantes
export(linalg)
J :=jacobian(
[x3, x4,
(m*sin(x2)*(g*cos(x2)-l*x4^2)+u)
/(M+m*sin(x2)^2),
(sin(x2)*((M+m)*g-m*l*x4^2*cos(x2))+cos(x2)*u)
/(l*(M+m*sin(x2)^2)),
x1],[x1,x2,x3,x4,u])
J0 :=subs(J(x1,x2,x3,x4,u),x1=0,x2=0,x3=0,x4=0,u=0)
J1 :=simplify(J0)
Il est aisé de véri…er que le système linéarisé de notre pendule inversé est observable et commandable pour
les valeurs nominales des paramètres. L’algorithme RegulKLH de la page 96 peut donc être poursuivi
a…n d’obtenir notre régulateur. Un programme de la simulation du pendule inversé régulé est donné dans
le …chier [Link]. Une version complète et commenté de ce programme est donnée
ci-dessous. Le programme est composé de trois fonctions et du programme principal.
//------------------------------------------------------
// fonction qui dessine le pendule inverse
function draw(x,couleur);
xbasc(); //efface l’ecran
isoview(-2,3,-1,2); //regle les echelles
//dessin du pendule
xpoly([ x(1), x(1)-l*sin(x(2))], [0, l*cos(x(2))]);
//dessin du chariot
xpoly([x(1)-.5,x(1)+.5,x(1)+.5,x(1)-.5,x(1)-.5],
[0,0,-.25,-.25,0]);
endfunction
//------------------------------------------------------
//fonction pour le calcul du regulateur
function [Ar,Br,Cr,Dr]=RegulKLH(A,B,C,E,pcom,pobs);
K=ppol(A,B,pcom); //matrice pour la regulation
L=ppol(A’,C’,pobs)’; //matrice pour l’observation
H=-inv(E*inv(A-B*K)*B); //matrice de pre-compensation
Commande linéaire des systèmes non linéaires 107
On rappelle (voir équation (2.5) page 27) que les équations d’état du monocycle sont de la forme
8 0 1
> 0 1 x3
>
> x_ 1
>
> B x4 C
>
> B x_ C B C
< B 2 C = B 3( 2 cos x2 ) sin x2 +( 2 +
C
B C B 2 x4 g 3 cos x2 )u C
@ x_ 3 A B 2 cos2 x C
>
> @ 1 2
2 x2 cos x
3 2
A
>
> x_ 4 ( 1 g 3 4 2 ) sin x2 ( 1 + 3 cos x2 )u
>
> 2 cos2 x
>
: 1 2 3 2
y = x1
où 1 ; 2 ; 3 ; g sont des constantes. Nous avons rajouté l’équation d’observation y = x1 a…n de supposer
une situation où seul l’angle de la roue est mesuré. Linéarisons avec Mupad ce système autour du point
de fonctionnement u = 0, x = 0 :
export(linalg)
J :=jacobian( [
x3, x4,
(m3*(m2*x4^2-mg*cos(x2))*sin(x2)+(m2+m3*cos(x2))*u)
/(m1*m2-m3^2*cos(x2)^2),
((m1*mg-m3^2*x4^2*cos(x2))*sin(x2)-(m1+m3*cos(x2))*u)
/(m1*m2-m3^2*cos(x2)^2),
x1],
[x1,x2,x3,x4,u])
J0 :=subs(J(x1,x2,x3,x4,u),x1=0,x2=0,x3=0,x4=0,u=0)
J1 :=simplify(J0)
On obtient 0 1
0 0 1 0 0
B 0 0 0 1 0 C
B C
dh B 2+
C
(x;u) = B C
3 g
0 2 0 0 3
2
d (x;u) B 1 2 3 1 2 3 C
B 0 1 g
0 0 1+ 3 C
@ 1 2
2
3 1 2
2
3
A
1 0 0 0 0
au point (x = (0; 0); u = 0). Le système linéarisé est décrit par les équations d’état suivantes :
8 0 1 0 1
>
> 0 0 1 0 0
>
> B C B C
>
> B 0 0 0 1 C B 0 C
>
< x_ = B C B Cu
B 0 3 g
0 0 C x + B 2 + 3
C
@ 1 2
2
3 A @ 1 2
2
3 A
>
> 1 g +
>
> 0 2 0 0 1 3
2
>
> 1 2 3 1 2 3
>
: y = 1 0 0 0 x
L’exécution complète de l’algorithme RegulNL du tableau 6.1 page 104 nous génère ainsi notre régulateur.
Un programme complet de la simulation du Segway stabilisé par ce régulateur est donné par le …chier
[Link].
Commande linéaire des systèmes non linéaires 109
La commande des systèmes non holonomes (comme les robots à roues) est loin d’être facile car souvent il
arrive que le système considéré soit commandable alors que son linéarisé ne l’est pas. Prenons le cas d’une
voiture à l’arrêt, il est possible de bouger le véhicule latéralement par une suite de manœuvres semblables à
celles e¤ectuées lorsqu’on essaie de garer sa voiture alors la place laissée libre est à peine su¢ sante. Pourtant,
cette même suite de manœuvres e¤ectuées sur le linéarisé de la voiture n’engendrera aucun mouvement
latéral. Une autre di¢ culté que nous allons rencontrer dans ce paragraphe est que parfois, on ne souhaite
pas stabiliser le système autour d’un point de fonctionnement mais lui donner un mouvement précis. C’est
le cas pour la conduite automatique : aller plein nord avec une vitesse de 50 km.h 1 ne correspond pas à un
point de fonctionnement pour la voiture car sa variable d’état y n’est pas constante (elle augmente à une
vitesse 50 km.h 1 et tend donc vers l’in…ni).
A travers ce paragraphe, nous allons montrer comment utiliser les techniques linéaires pour arriver à conce-
voir un régulateur pour la conduite automatique. Nous aurons besoin pour cela d’imaginer un monde idéal
pour le régulateur où le modèle associé admet un point de fonctionnement qui correspond au comportement
voulu pour notre système réel.
Dans ce paragraphe, nous allons traiter le cas d’une voiture (celle modélisée au paragraphe 2.4.3 page 28)
qui roule sur une route inconnue. La voiture est équipée
– d’un télémètre mesurant la distance latérale d entre le milieu de l’essieu arrière de la voiture et le bord
de la route (voir …gure),
– d’un capteur de vitesse qui mesure la vitesse v des roues avant,
– et d’un capteur d’angle qui mesure l’angle du volant (pour simpli…er, nous allons supposer que cor-
respond aussi à l’angle que forment les roues avant avec l’axe de la voiture).
On souhaite que la voiture se déplace à vitesse constante le long de cette route. Bien sûr, il est hors de
question de mettre dans la connaissance du régulateur la forme de la route car elle est inconnue. De plus, les
110 Commande par espace d’état
variables de position et d’orientation de la voiture ne sont pas mesurées (c’est-à-dire qu’on ne dispose ni de
boussole, ni de capteur GPS). Or ces quantités sont souvent inutiles pour atteindre l’objectif …xé, à savoir
suivre la route. En e¤et, ne sommes nous pas, nous-mêmes, capables de conduire une voiture sur une route
sans avoir la carte des environs, sans savoir où nous sommes et où est le nord ? Ce qui nous intéresse lorsque
nous conduisons une voiture, c’est la position relative de la voiture par rapport à la route. Le monde idéal
pour le régulateur est évoqué par la …gure.
Il doit être tel que le modèle utilisé par le régulateur possède des variables d’état décrivant cette position
relative de la voiture par rapport à la route et que la constance des variables d’état de ce modèle corresponde
à une situation réelle où la voiture suit la route à une vitesse constante.
Tout d’abord, dans le paragraphe 6.6.1, nous allons donner la fonction d’observation de notre voiture roulant
sur un circuit en forme de polygone (voir …gure).
Ensuite, nous imaginerons dans le paragraphe 6.6.2 un monde idéal pour le régulateur. Au paragraphe 6.6.3
nous implémenterons le régulateur ainsi conçu sur une voiture tournant autour d’un circuit.
Commande linéaire des systèmes non linéaires 111
Rappelons que la fonction d’évolution de notre voiture est déjà connue (voir formule (2.7) page 29). Pour
l’expression de la fonction d’observation, nous avons besoin de calculer la distance renvoyée par le télémètre
et qui correspond à la distance entre le milieu de l’essieu arrière et le polygone formé par le circuit. Pour
cela, le théorème suivant nous sera utile.
Remarque 6.6.1 Pour comprendre la preuve qui suit, il est important de se rappeler la signi…cation du
signe du déterminant de deux vecteurs ! u et ! v de R2 . On a (i) det(! u;!v ) > 0 si !
v est sur la gauche de
u , (ii) det(!
! u;!
v ) < 0 si !
v est sur la droite de !
u et (iii) det(!u;!v ) = 0 si !
u et !
v son colinéaires. Ainsi,
par exemple, sur la …gure ci-dessous à gauche, det(a m; ! u ) > 0 et det(b m; ! u ) < 0: Rappelons aussi
que le déterminant est une forme multilinéaire, c’est-à-dire que
det(au + bv; cx + dy) = a det(u; cx + dy) + b det(v; cx + dy)
= ac det(u; x) + bc det(v; x)
+ad det(u; y) + bd det(v; y):
Preuve : La droite D m;! u qui passe par le point m et de vecteur directeur ! u coupe le plan en
!
2 demi-plans : ceux qui satisfont det(z m; u ) 0 et ceux qui satisfont det(z m; ! u) 0. Elle
coupe donc le segment [ab] si a et b sont dans des demi-plans di¤érents (voir …gure), c’est-à-dire si
det a m; ! u : det b m; ! u 0:
La droite D (m; u
~ ) coupe le segment [ab]
112 Commande par espace d’état
det(m a; b a) + det(!
u;b a) = 0:
En isolant , on obtient
det(a m; b a)
= :
det(!
u ; b a)
Si 0, alors représente la distance d parcourue par le faisceau partant de m dans la direction ! u avant
de rencontrer le segment. Si < 0, cela signi…e que le faisceau ne rencontrera jamais le segment car ce
dernier est du mauvais côté. La condition 0 correspond à la deuxième inégalité de la condition (6.4).
La fonction d’observation pour notre système, donnée par le programme Scilab du tableau 6.2, est une
conséquence directe de ce théorème. Rappelons que cette fonction nous renvoie la distance d mesurée par
le télémètre, la vitesse v des roues avant et l’angle du volant . La ligne 2 calcule le vecteur u qui indique
la direction du télémètre laser et le point m, d’où part le laser. La distance d renvoyée doit être la plus
petite (voir ligne 7) parmi toutes les distances susceptibles (voir lignes 5 et 6) d’être renvoyées par chacun
des segments. Dans ce programme, Aj et Bj représentent des matrices contenant les points extrêmes aj et
bj des sommets du jième segment. Notons en…n que cette fonction ne sera pas utilisée par le régulateur,
mais uniquement a…n de simuler notre système réel. En e¤et, notre régulateur ne connaît pas la forme du
polygone autour duquel la voiture est censée tourner : il pense que la voiture est en train de longer un mur
droit et in…ni. Sa fonction d’observation est donnée au paragraphe suivant.
Commande linéaire des systèmes non linéaires 113
1 function y=g(x)
2 u=[sin(x(3)) ;-cos(x(3))] ; m=[x(1) ;x(2)] ; d=%inf ;
3 for j=1 :size(Aj,2)
4 a=Aj( :,j) ; b=Bj( :,j) ;
5 if (det([a-m u])*det([b-m u]) <= 0)&
6 (det([a-m b-a])*det([u b-a])>=0)
7 d=min(det([a-m b-a])/det([u b-a]),d) ;
8 end
9 end
10 y=[d ;x(4) ;x(5)] ;
11 endfunction
A…n de concevoir le modèle pour notre régulateur, intéressons nous à la position relative de la voiture
par rapport au bord de la route. Nous allons pour cela supposer que notre voiture roule à une distance x
du bord. Ce modèle doit posséder les mêmes entrées et sorties que le système réel, à savoir deux entrées
(l’accélération des roues avant v_ et la vitesse angulaire du volant _ ) et trois sorties (la distance d du milieu
de l’essieu arrière au bord de la route, la vitesse v des roues avant et l’angle du volant). Du fait que seule
la position relative de la voiture nous intéresse, notre modèle ne doit posséder que quatre variables d’états :
x; ; v; , comme représenté sur la …gure. Attention, la signi…cation de la variable x a légèrement changé par
rapport au contexte du paragraphe 2.4.3 page 28. En s’aidant de la modélisation faite au paragraphe 2.4.3,
on obtient pour notre modèle, les équations d’état suivantes
8 0 1 0 1
>
> x_ v cos cos
>
> B _ C B C
>
> B C B
v sin
C
>
> B C = B L C
>
> @ v_ A @ u1 A
<
_ u2 1
> 0
>
> x
>
> sin
>
> B C
>
> y = @ v A;
>
:
où L est la distance entre les essieux avant et arrière. Choisissons pour point de fonctionnement
qui correspond à une vitesse de 7 ms 1 et une distance de 5 m entre le milieu de l’essieu arrière et le bord
de la route. La matrice jacobienne de la fonction d’évolution-observation
8
>
> R6 ! 0 R7 1
>
>
>
> 0 1 v cos cos
>
> x B C
>
> B v sin
C
>
> B C B L C
< B C B u1 C
h: B v C B C
B C
>
> B C ! B
B u2 C
C
>
> B C B C
>
> B C B
x
C
>
> @ u A B sin
C
>
> 1
@ v A
>
> u2
:
Cette matrice a été obtenue par Mupad grâce aux instructions suivantes.
export(linalg)
J :=jacobian([
-v*cos(delta)*cos(theta),
v*sin(delta)/L,
u1, u2, x/sin(theta), v, delta],
[x,theta,v,delta,u1,u2])
J0 :=subs(J(x,theta,v,delta,u1,u2),
x=5,theta=PI/2,v=7,delta=0,u1=0,u2=0)
J1 :=subs(J0(L),L=3)
J2 :=simplify(J1)
Mupad obtient 0 1
0 7 0 0 0 0
B 7 C
B 0 0 0 3 0 0 C
B C !
B 0 0 0 0 1 0 C
dh B C A B
(x; u) = B
B 0 0 0 0 0 1 C=
C :
d (x; u) B C C D
B 1 0 0 0 0 0 C
B C
@ 0 0 1 0 0 0 A
0 0 0 1 0 0
Commande linéaire des systèmes non linéaires 115
Prenons pour variables consignées x et v, ce qui signi…e que nous voulons nos consignes w1 et w2 corres-
pondent à la distance au bord de la route et la vitesse de la voiture. La matrice de consigne est donc
!
1 0 0 0
E= :
0 0 1 0
l’algorithme RegulKLH(in : A; B; C; E; pcom ; pobs ) du tableau 5.2 page 96 nous génère le régulateur sui-
vant :
8 0 1 0 1 0 1
>
> 4 7 0 0 0 0 4 0 0 ! 20
>
> B C B
>
>
> d B 0:6 0 0 0 C B 0 0 0:6 0 2:3 C
C w B
B 2:9 C
C
>
< dt x
^ =B Cx^+B C B C
@ 0 0 4 0 A @ 0 2 0 2 0 A y @ 28 A
>
> 0:5 5:1 0 8 ! 0:5 0 !0 0 2 !! 2:4
>
>
>
> 0 0 2 0 0 2 5
>
>
: u= 0:49 5:14 0 6
^+
x
0:48 0
w
7
:
Une simulation du comportement de ce régulateur sur ce modèle de voiture simpli…é (où un seul mur droit
existe) est donnée par le …chier [Link].
6.6.3 Bouclage
Le régulateur que nous venons juste de proposer a été conçu pour un modèle idéalisé de notre voiture, mais
qui ne correspond pas exactement à notre système réel. Rappelons que notre voiture est un système d’ordre
5 alors que sa version idéalisée utilisée par le régulateur est d’ordre 4. Or, ces deux systèmes admettent les
mêmes entrées (u1 et u2 ) et les mêmes sorties (d; v; ). Il est donc possible de brancher le régulateur conçu
pour notre système idéal sur le système réel. Il nous reste à espérer que le comportement engendré pour
notre système sera convenable.
Le programme de simulation de notre commande est proche de celui du tableau 3.9, page 57. La même
fonction d’évolution f est utilisée. En ce qui concerne la fonction draw, il nous faut rajouter le dessin du
116 Commande par espace d’état
Programme [Link]
A=[0 7 0 0 ;0 0 0 7/3 ;0 0 0 0 ;0 0 0 0] ;
B=[0 0 ;0 0 ;1 0 ;0 1] ;
C=[1 0 0 0 ;0 0 1 0 ;0 0 0 1] ;
E=[1 0 0 0 ;0 0 1 0] ;
K=ppol(A,B,-2*[1,1,1,1]) ;
L=(ppol(A’,C’,-2*[1,1,1,1]))’ ;
H=-inv(E*inv(A-B*K)*B) ;
Ar=A-B*K-L*C ;Br=[B*H L] ;Cr=-K ;
ubar=[0 ;0] ;xbar=[5 ;1.57 ;7 ;0] ;wbar=E*xbar ;
ybar=[xbar(1)/sin(xbar(2)) ;xbar(3) ;xbar(4)] ;
xr=[0 ;0 ;0 ;0] ;w=[5 ;7] ;
x=[-15 ;0 ;%pi/2 ;7 ;0.1] ;dt=0.07 ;
for k=1 :500,
u=ubar+Cr*xr+H*(w-wbar) ;y=g(x) ;
x1=x+f(x,u)*dt ;
xr1=xr+(Ar*xr+Br*[w-wbar ;y-ybar])*dt ;
x=x1 ; xr=xr1 ;
end ;
faisceau lancé par le télémètre ainsi que le dessin du circuit. Cela correspond aux instructions suivantes :
Aj=[-10 -10 0 10 20 32 35 30 20 0 ;
-5 5 15 20 20 15 10 0 -3 -6 ] ;
Bj=[Aj( :,2 :$),Aj( :,1)] ;
xfpoly(Aj(1, :),Aj(2, :)) ;
xpoly([x(1),x(1)+sin(x(3))*d],
[x(2),x(2)-cos(x(3))*d]) ;
Le programme de simulation de la voiture sur circuit est donné par le tableau 6.3. En lançant ce programme,
on s’aperçoit que la voiture tourne correctement autour du circuit, mais que lorsqu’elle franchit un sommet
anguleux du polygone, elle est déstabilisée pendant une petite seconde. Ceci est dû au fait que le régulateur
(pensant que la route est droite à l’in…ni) est surpris par ce changement brutal de direction. Ces changements
de direction dans la route sont considérés par le régulateur comme des perturbations.
Chapitre 7
7.1 Introduction
Au chapitre 6, nous avons montré comment concevoir des régulateurs pour les systèmes non linéaires, lorsque
le vecteur d’état reste localisé dans une zone de petite taille de l’espace d’état. L’idée était de linéariser le
système non linéaire à réguler et d’utiliser les techniques développées au chapitre 5 pour la régulation des
systèmes linéaires. Dans ce chapitre, nous allons chercher à proposer des régulateurs non linéaires a…n de
contraindre le vecteur d’état à rester dans une petite zone de l’espace d’état. Contrairement à l’approche
linéaire qui o¤re une méthodologie générale mais limitée au voisinage d’un point de l’espace d’état, les
approches non linéaires ne s’appliquent qu’à des classes limitées de systèmes, mais elles permettent d’étendre
la plage de bon fonctionnement du système.
Il n’existe en e¤et pas de méthode générale pour stabiliser globalement les systèmes non linéaires. En
revanche, il existe une multitude de méthodes qui s’appliquent à des cas particuliers. Le but de ce chapitre
est d’en présenter une assez représentative des méthodes existantes. Il s’agit de la méthode de linéarisation
par bouclage. Après l’avoir introduite par un petit exemple au paragraphe 7.2 et présenté son principe au
paragraphe 7.3, nous appliquerons cette méthode à la régulation (par retour d’état) du problème des 3 bacs
modélisé au paragraphe 2.5.2 de la page 34, à la commande de la voiture modélisée au paragraphe 2.4.3 de
la page 28, puis à la commande du bateau à voile modélisé au paragraphe 2.4.4 page 31.
117
118 Commande par espace d’état
On souhaiterait que sa sortie y soit égale à une consigne supposée constante v. En dérivant la sortie y nous
obtenons
y_ = x_ exp (x) = x3 exp (x) + x exp (x) u:
Posons
1
u= x3 exp (x) + v ; (7.2)
x exp (x)
où v correspond à la nouvelle entrée (ou consigne). Le bouclage (7.2) transforme le système (7.1) en un
simple intégrateur décrit par l’équation di¤érentielle
SL : y_ = v:
Le bouclage (7.2) s’appelle bouclage linéarisant car il transforme le système non linéaire en un système
linéaire. Le système ainsi obtenu peut être stabilisé par des techniques linéaires classiques. Nous choisirons
ici une commande PI (proportionnelle et intégrale) de la forme
Rt
v = P (w y) + I 0 (w ( ) y ( )) d ;
où w est la consigne. Les équations d’état de cette commande PI sont données par
(
z_ = w y
v = I z+ P (w y) ;
où z est la variable d’état associée à l’intégrateur. Donc les équations d’état d’un régulateur par retour d’état
RNL pour notre système non linéaire (7.1) sont données par
(
z_ = w exp(x)
RNL : 1
u = x exp(x) x3 exp (x) + Iz + P (w exp(x)) :
7.3.1 Principe
Nous cherchons à généraliser ici la méthode décrite au paragraphe précédent. Une approche plus détaillée
pourra être trouvée dans le livre d’Andréa-Novel [3]. Considérons le système non linéaire décrit par
(
x_ = f (x) + g(x)u
y = h(x)
dont le nombre d’entrées et de le nombre de sorties sont tous les deux égaux à m. L’idée de la linéarisation
par bouclage, est de boucler le système par une commande du type u = r(x; v), où v est la nouvelle entrée,
aussi de dimension m. Cette opération nécessite que l’état soit complètement accessible. Si cela n’est pas le
cas, on se doit de fabriquer un observateur, mais dans un contexte non linéaire, cela est une opération très
di¢ cile. Puisque ici l’état est supposé accessible, y ne doit pas vraiment être considérée comme une sortie,
mais plutôt comme le vecteur des variables consignées.
Pour e¤ectuer ce bouclage, il nous faut exprimer les dérivées successives de chacun des yi en fonction de
l’état et de l’entrée. On s’arrête de dériver yi ; dès que les entrées commencent à intervenir dans l’expression
de la dérivée. Nous disposons ainsi d’une équation du type
0 1
(k )
y1 1
B . C
B .. C = A (x) u + b (x) ; (7.3)
@ A
(k )
ym m
où ki désigne le nombre de fois qu’il nous faut dériver yi pour y voir apparaître une entrée (voir les exemples
dans les paragraphes qui suivent pour une meilleure compréhension). Sous l’hypothèse que la matrice A(x)
soit inversible, le bouclage suivant
u = A 1 (x) (v b (x)) ;
où v est notre nouvelle entrée (voir …gure ci-dessous), forme un système linéaire SL de m entrées à m sorties
décrit par les équations di¤érentielles 8
> (k1 )
< y1
> = v1
SL : .
.. ..
> = .
>
: y (km ) = v :
m m
120 Commande par espace d’état
Ce système est linéaire et complètement découplé (c’est-à-dire que chaque entrée vi agit sur une et une seule
sortie yi ). Il est donc très facile à commander par les techniques classiques de commande par retour de sortie
(voir paragraphe 5.4 page 90). Ici, comme le système à commander est constitué de chaînes d’intégrateurs
découplés, nous utiliserons m régulateurs de types PID (proportionnel intégral et dérivées) dont nous rap-
pellerons les principes au paragraphe 7.3.2. Notons que pour utiliser de tels régulateurs, il est nécessaire de
disposer des dérivées des sorties. La remarque suivante explique comment les obtenir dans notre contexte
de linéarisation par bouclage.
Remarque 7.3.1 En analysant bien le cheminement pour obtenir l’équation (7.3), on se rend compte que
(j)
la j ieme dérivée de la iieme sortie yi s’exprime sous la forme
(j)
yi = ^bij (x) si j < ki
(j) ^
yi = ^aT
ij (x):u + bij (x) si j = ki
(j)
yi = ^
aij (x; u; u;
_ u• ; : : : ) si j > ki
Le coe¢ cient ki s’appelle le degré relatif de la iieme sortie. Si on mesure l’état du système x et son entrée u,
(j)
on peut donc disposer de toutes les dérivées successives des sorties yi , tant que j reste inférieur ou égal à ki .
En e¤ et, en raison des bruits hautes fréquences apparaissant dans les signaux, on ne peut pas disposer de la
dérivée des signaux par l’utilisation de dérivateurs de façon …able. Nous avons donc une fonction analytique
qui nous permet de disposer de toutes les dérivées des sorties (jusqu’à leur degré relatif ), et ceci sans avoir
à utiliser de dérivateurs de signaux. Une autre approche pour obtenir les dérivées des sorties aurait été de
construire un observateur d’état (voir paragraphe 5.4 page 90).
Nous avons
y = 2x
y_ = 2x_ = 2xu + 2x3
y• = 2xu _ 2 = 2(xu + x3 )u + 2xu_ + 6(xu + x3 )x2 :
_ + 2xu_ + 6xx
Nous avons donc un degré relatif k = 1 pour la sortie y: On peut donc disposer de y_ sans avoir à utiliser de
dérivateurs de signaux. Cela n’est pas le cas pour y• car disposer de u avec une bonne précision ne signi…e
pas que l’on dispose de u.
_ Ici,
(x; u) = (2x; 2xu + 2x3 ):
Commande non linéaire des systèmes non linéaires 121
y (n) = u:
u= 0 (w y) + 1 (w
_ y)
_ + n 1 w(n 1)
y (n 1)
+ w(n) :
où w est la consigne souhaitée pour y. Notons que w n’est pas ici supposée constante, mais peut dépendre
du temps. Le fait que ce régulateur nécessite les dérivées de y n’est pas un problème dans le cadre dé…ni par
la linéarisation par bouclage. En e¤et, toutes ces dérivées peuvent être écrite comme fonction analytique de
l’état x du système et de l’entrée u (voir remarque 7.3.1). En ce qui concerne la consigne w(t), elle est choisie
par l’utilisateur et une expression analytique de w(t) peut être supposée connue (par exemple w(t) = sin(t)).
Ainsi, le calcul des dérivées de w se fait de façon formelle et aucune sensibilité de l’opérateur dérivation par
rapport au bruit n’est à craindre.
y (n) = u = 0 (w y) + 1 (w
_ y)
_ + n 1 w(n 1)
y (n 1)
+ w(n) :
e(n) + n 1e
(n 1)
+ + 1 e_ + 0e = 0:
Cette équation di¤érentielle est appelée dynamique de l’erreur. Son polynôme caractéristique
P (s) = sn + n 1s
n 1
+ + 1s + 0; (7.4)
peut donc être choisi arbitrairement parmi les polynômes de degré n: Bien sûr, on choisira un polynôme
dont les racines sont toutes à parties réelles négatives, a…n d’assurer la stabilité du système. Par exemple,
si n = 3 et si on souhaite que tous les pôles soient égaux à 1, on posera
s3 + 2s
2
+ 1s + 0 = (s + 1)3 = s3 + 3s2 + 3s + 1:
D’où
2 = 3; 1 = 3; 0 = 1:
Le régulateur PI alors obtenu est donné par
...
u = (w y) + 3 (w_ y)
_ + 3 (w
• y•) + w: (7.5)
Régulateur proportionnel intégral et dérivées : A…n de compenser les perturbations constantes, nous
pouvons décider de rajouter un terme intégral. Nous obtenons une commande de type PID de la forme
Z t
u = 1 (w( ) y( )) d (7.6)
=0
+ 0 (w y) + 1 (w
_ y)
_ + n 1 w(n 1)
y (n 1)
+ w(n) :
122 Commande par espace d’état
+ 0 (w y) + 1 (w
_ y)
_ + n 1 w(n 1)
y (n 1)
+ w(n) :
e(n+1) + n 1e
(n)
+ + 1e
•+ 0 e_ + 1e = 0:
P (s) = sn+1 + n 1s
n
+ + 1s
2
+ 0s + 1
Pour une étude plus complète des régulateurs de type PID, le lecteur est encouragé à consulter le livre de
Landau [11].
A…n d’illustrer le principe de la linéarisation par bouclage, reprenons à nouveau le système composé de trois
bacs modélisé au paragraphe 2.5.2 page 34. Les équations d’état obtenues pour ce système sont
8
>
> h_ 1 = (h1 ) (h1 h2 ) + u1
>
> _
>
< 2h = (h 1 h 2 ) (h2 h3 )
_h3 = (h3 ) + (h2 h3 ) + u2 : (7.7)
>
>
>
> y 1 = h1
>
:
y 2 = h3
p
où (h) = a:sign (h) 2gjhj: Nous avons ici choisi pour sorties les hauteurs dans le premier et le troisième
bacs. Les dérivées des sorties y1 et y2 s’expriment par
y_ 1 = h_ 1 = (h1 ) (h1 h2 ) + u 1
y_ 2 = h_ 3 = (h3 ) + (h2 h3 ) + u 2
Le bouclage suivant
!
1 (h1 ) (h1 h2 )
u=A (x) (v b (x)) = v (7.8)
(h3 ) + (h2 h3 )
Commande non linéaire des systèmes non linéaires 123
où v est notre nouvelle entrée, rend notre système linéaire. Plus précisément ce dernier a la forme
(
y_ 1 = v1
SL : (7.9)
y_ 2 = v2 :
Cherchons maintenant à réguler le système linéaire (7.9) par un régulateur composé de deux régulateurs PI
(proportionnel et intégrale) de la forme
( Rt
v1 (t) = 0 (w1 (t) y1 (t)) + 1 0 (w1 ( ) y1 ( )) d + w_ 1
RL : Rt
v2 (t) = 0 (w2 (t) y2 (t)) + 1 0 (w2 ( ) y2 ( )) d + w_ 2
où w1 et w2 sont les nouvelles consignes pour y1 et y2 . Si on souhaite avoir pour pôles uniquement des 1,
il faut (voir équation (??)) que
(
s2 + 0 s + 1 = (s + 1)2 = s2 + 2s + 1
s2 + 0 s + 1 = (s + 1)2 = s2 + 2s + 1;
soit
1 = 1 = 1; 0 = 0 = 2:
En prenant en compte l’équation (7.8), nous obtenons que les équations d’état d’un régulateur par retour
d’état RNL pour notre système non linéaire (7.7) sont données par
8
>
> z_1 = w1 h1
>
< z_ = w
2 2 h3
RNL :
>
> u1 = z1 + 2 (w1 h1 ) + w_ 1 + (h1 ) + (h1 h2 )
>
:
u2 = z2 + 2 (w2 h3 ) + w_ 2 + (h3 ) (h2 h3 ) :
On rappelle que les équations d’état du char (voir paragraphe 2.4.2 page 27) sont données par
0 1 0 vg +vd 1
x_ 2 cos
B y_ C B vg +vd
sin C
B C B 2 C
B _ C B v d vg C
B C=B ` C:
B C B C
@ v_ g A @ Ru1 A
v_ d Ru2
124 Commande par espace d’état
y_ 1 = _ = vd vg
(on dérive encore une fois)
`
R
y•1 = • = (u2 u1 )
`
R
y_ 2 = (u2 + u1 ) ;
2
soit ! ! !
R R
y•1 ` ` u1
= R R :
y_ 2 2 2 u2
| {z }
A(x)
Si on e¤ectue le bouclage
! ! ! !
u1 1 v1 1 ` 2 v1
=A (x) = ;
u2 v2 2R ` 2 v2
qui peut lui aussi être commandé par un deuxième bouclage. Choisissons une commande proportionnelle et
dérivée de façon à avoir uniquement des pôles égaux à 1 :
(
v d vg
v1 = (w1 y1 ) + 2 (w_ 1 y_ 1 ) + w
•1 = w1 + 2 w_ 1 ` +w
•1
vg +vd
v2 = w2 y2 + w_ 2 = w2 2 + w_ 2 :
On rappelle que l’équation d’évolution de la voiture (voir page 29) est donnée par
0 1 0 1
x_ v cos cos
B y_ C B v cos sin C
B C B C
B _ C B v sin C
B C=B L C:
B C B C
@ v_ A @ u1 A
_ u2
Le but de ce paragraphe est de faire décrire à notre voiture diverses trajectoires en utilisant le principe de
la linéarisation par bouclage.
Choisissons pour sortie le vecteur y = (v; )T : La dérivée des sorties y1 et y2 s’exprime par
y_ 1 = v_ = u1 ;
v sin
y_ 2 = _ = :
L
Comme la dérivée y_ 2 de y2 ne fait pas apparaître l’entrée, on la redérive une nouvelle fois :
En imposant le bouclage
1
u=A (x) v;
et devient donc linéaire et découplé. Nous avons donc a¤aire à deux systèmes monovariables découplés. Le
premier, d’ordre 1, peut être stabilisé par une commande proportionnelle. Pour le second, du deuxième ordre,
une commande de type proportionnelle et dérivée) est adaptée. Si w = (w1 ; w2 )T représente la consigne pour
y, cette commande s’exprime par
(
v1 = 0 (w1 y1 ) + w_ 1
RL :
v2 = 0 (w2 y2 ) + 1 (w_ 2 y_ 2 ) + w
•2 :
126 Commande par espace d’état
Si on souhaite avoir pour pôles que des 1, il faut (voir équation (7.4))
s+ 0 = s+1
s2 + 1s + 0 = (s + 1)2 = s2 + 2s + 1;
soit
0 = 0 = 1; 1 = 2:
Donc les équations d’un régulateur par retour d’état RNL pour notre système non linéaire sont données par
! !
1 0 0 (w1 v) + w_ 1
u= tan L (7.10)
v v cos 0 (w2 )+ 1 w_ 2 v sin
L +w
•2
Notons que ce régulateur n’admet aucune variable d’état. Il s’agit donc d’un régulateur statique.
A…n de faire parcourir à notre voiture un hexagone en marche arrière, à une vitesse de 2ms 1, nous
prendrons pour consigne
w1 (t) = 2
w2 (t) = E (t=10) :
4
Nous voulons maintenant faire tourner la voiture autour d’un cercle de rayon r0 et à la vitesse v0 . Le nouveau
vecteur de consigne devient c = (r0 ; v0 ). Pour cela, choisissons dans le plan (x; y) un champ de vecteur h
donné par
R2 ! ! R2 !
p
h: x vr x tan 1 ( (r0 p x2 + y 2 )) + v y
0
7! p 1 :
y x2 +y 2 vr y tan 1 ( (r0 x2 + y 2 )) v0 x
Le programme qui trace ce champ a été appelé [Link]. Le système associé à ce champ
(c’est-à-dire qu’il est décrit par l’équation di¤érentielle z_ = h(z)) admet un cycle limite qui est un cercle
de rayon r0 . Sur ce cercle, nous avons jj_zjj = v0 . Les quantités vr et indiquent respectivement avec quelle
vitesse l’état z se rapproche du cercle limite lorsque ce dernier est loin ou proche du disque. A…n de faire
suivre aux coordonnées (x; y) de la voiture, cette dynamique, il semble naturel d’imposer
! !
x_ x
=h = h (c; x; y) :
y_ y
La notation h (c; x; y) nous rappelle la dépendance du champ en c = (r0 ; v0 ) : Or,
! ! !
x_ v cos cos w1 cos w2
= '
y_ v cos sin w1 sin w2
car w1 est une consigne pour la vitesse de la voiture et w2 pour son cap. Pour choisir w nous devons donc
tenter de résoudre !
w1 cos w2
= h (c; x; y) :
w1 sin w2
Ainsi, nous obtenons
w1 = kh (c; x; y)k et w2 = arg (h (c; x; y)) ;
où arg retourne l’angle entre + et + du vecteur h(c; x; y). Ainsi, d’après la formule (7.10), les
équations de notre régulateur sont
! !
1 0 0 (kh (c; x; y)k v)
u= tan L v sin :
v v cos 0 (arg (h (c; x; y)) ) 1 L
128 Commande par espace d’état
Ici, nous avons pris w_ 2 = 0 alors que w2 est loin d’être constante. Pour plus de précision, il faudrait la calculer
et mettre son expression dans celle de la commande u. Le régulateur obtenu est à nouveau un régulateur
statique (c’est-à-dire, sans état). Si nous appliquons cette commande, sur la voiture, nous nous apercevons
que la voiture tourne e¤ectivement autour d’un cercle, mais qu’il existe un large biais. Plus précisément,
le rayon du cercle parcouru et la vitesse de la voiture sur ce cercle sont di¤érents de ceux rentrés dans
la consigne c et que cette di¤érence n’est pas négligeable. Ce phénomène est principalement dû au fait la
consigne a été supposée constante alors qu’elle ne l’est pas. A…n de supprimer simplement ce biais, nous
pouvons rajouter un e¤et intégrateur. En choisissant un coe¢ cient de 15 pour l’intégrateur, les équations
d’état obtenues pour le régulateur sont
8 p !
>
> x2 + y 2
>
> = c
< z_ v
! !
>
> 1 0 h c + z5 ; x; y v
>
>
: u = tan L
arg h c + z5 ; x; y 2 v sin
:
v v cos L
Considérons à nouveau le bateau à voile modélisé au paragraphe 2.4.4 de la page 31. Rappelons ses équations
d’état : 8
>
> x_ = v cos ;
>
>
>
> y_ = v sin V;
>
>
>
> _ = !;
>
>
>
> _
< v = u1 ;
_g = u2 ; (7.11)
>
>
>
> v_ = 1
m (fv sin v fg sin g f v) ;
>
>
>
> 1
!_ = J ((` rv cos v ) fv rg cos g fg !) ;
>
>
>
> fv =
>
> v (V cos ( + v ) v sin v ) ;
:
fg = g v sin g :
Nous allons ici chercher à stabiliser les valeurs de x et y autour de certaines valeurs désirées xd et yd : Notons
tout d’abord qu’il n’est pas possible d’obtenir une stabilisation exacte si V 6= 0, même dans le cas idéal,
car on ne peut annuler simultanément x, _ y;
_ x•; y•. On peut le montrer formellement, bien que cette opération
soit loin d’être triviale, même avec l’aide d’un logiciel de calcul formel comme Mupad[1]. On peut aussi
prouver l’absence de point d’équilibre grâce à un solveur numérique utilisant le calcul par intervalles [9].
Puisqu’une stabilisation exacte ne peut être envisagée, nous allons tenter une stabilisation approximative de
x et y autour des valeurs de consigne.
A…n de bien comprendre l’approche choisie, considérons la situation d’un moniteur de bateau à voile, dans
son Zodiac, dictant à distance à un apprenti (qui apprend à manœuvrer le bateau) comment amener le
bateau sur une cible (xd ; yd ). Le moniteur donnera à l’apprenti des consignes du type « prend la direction
de la pointe» , ou bien « relâche un peu la voile» et l’apprenti cherchera à suivre ces consignes. On fait alors
apparaître deux types de régulation.
– la régulation de bas niveau, faite par l’apprenti qui consiste à tirer sur l’écoute de la voile et manœuvrer
le gouvernail dans le but de respecter les consignes du moniteur ;
– la régulation de haut niveau faite par le moniteur, qui consiste à générer une suite de consignes (sur le
cap et sur l’ouverture de la voile) dans le but d’arriver à l’objectif choisi, à savoir atteindre la cible.
Dans une première étape, pour réaliser la commande de bas niveau, nous allons utiliser une méthode de
linéarisation par bouclage statique qui nous permettra de stabiliser le cap du bateau ainsi que l’ouverture de
la voile. Dans une deuxième étape, pour la régulation de haut niveau, nous proposerons un second régulateur
à états discrets qui nous permettra de maintenir le bateau dans une zone autour du point (xd ; yd ).
Rappelons que pour appliquer une méthode de linéarisation par bouclage, il nous faut tout d’abord dériver
les variables d’état, une ou plusieurs fois par rapport au temps t, et ceci jusqu’à faire apparaître les entrées
u1 ou u2 .
Dans les quantités exprimées par les équations (7.11), seules _ v et _ g sont liées algébriquement à u. Les autres
130 Commande par espace d’état
quantités sont reliées à u, mais indirectement (c’est-à-dire di¤érentiellement). Rangeons-les dans la liste L
des variables dépendant algébriquement de u. On a donc L = f _ v ; _ g g. Il nous faut dériver uniquement celles
qui ne dépendent pas algébriquement de u, c’est-à-dire, x; _ _ ; v_ et !.
_ y; _ On obtient
8
>
> x• = v_ cos v _ sin ;
>
>
>
> y• = v_ sin + v _ cos ;
>
>
>
< • = !;
_
f_v sin v +fv u1 cos vf_g sin g fg u2 cos g f v_ (7.12)
>
> v• = m ;
>
> u1 rv sin v fv +(` rv cos v )f_v
>
> !
• =
>
> J
>
: rg (u2 sin g fg cos g f_g ) !_
+ J :
avec
8
< f_v =
> vV (! + u1 ) sin ( + v )
v v
_ sin v v vu1 cos v
>
: f_ =
g g (v_ sin g + vu2 cos g ) :
Notons que, puisque les grandeurs v;_ _ ; !_ sont des fonctions analytiques de l’état (voir (7.11)), on peut consi-
dérer que nous avons ici une expression analytique de x •; y•; •; v•; !
• en fonction de l’état et des entrées. Rangeons
dans L les quantités qui dépendent algébriquement de u; c’est-à-dire v• et ! • , d’où L = f _ v ; _ g ; v•; !
• ; f_v ; f_g g.
Et dérivons à nouveau celles qui ne dépendent pas algébriquement de u; c’est-à-dire x •, y• et •: On obtient
8 ... 2
>
< x = v• cos 2v_ _ sin v • sin v _ cos
... 2
y = v• sin + 2v_ _ cos + v • cos v _ sin (7.13)
>
: ...
= !
•:
Puisque toutes ces quantités dépendent algébriquement de u, on s’arrête ici de dériver. Notons à nouveau,
que les équations (7.13) peuvent s’interpréter (via les équations (7.12) et (7.11)) comme des expressions
... ... ...
analytiques des grandeurs x ; y ; en fonction de l’état et des entrées.
Rappelons que les sorties (ce sont en fait des variables consignées) choisies sont l’ouverture de la voile y1 = v
et le cap y2 = : Nous avons
! ! ! !
y_ 1 _ 1 0 u1
... = ...v = rv rg
y2 J fv sin v J fg sin g u2
| {z }
A1 (x)
! ! !
0 0 f_v 0
+ rg + :
`
J
rv
J cos v J cos g f_g J !_
| {z } | {z }
A2 (x) b1 (x)
Commande non linéaire des systèmes non linéaires 131
Or
! ! !
f_v v (V sin ( + v) + v cos v) 0 u1
=
f_g 0 g v cos g u2
| {z }
A3 (x)
!
v (V ! sin ( + v ) + v_ sin v)
+ :
g v_ sin g
| {z }
b2 (x)
qui sont linéaires et découplées. Le système linéarisé est d’ordre 4 au lieu de 7. Nous avons donc perdu le
contrôle sur 3 variables qui se trouvent être x; y et v. La perte de contrôle sur x et y était prévisible (on
veut que le bateau avance et il est naturel que cela corresponde à une instabilité pour ces deux variables
x et y). Quant à la perte de contrôle sur v, elle est sans conséquence car la dynamique associée est stable.
Comment en e¤et concevoir que le bateau puisse maintenir un cap et une ouverture de voile …xes, sans que
sa vitesse ne converge vers une valeur …nie ?
Déterminons maintenant les singularités de notre bouclage linéarisant. On montre aisément que
rg
det (A(x)) = fg sin g v g cos2 g ;
J
est nul si
v 2 sin2 g 1 = 0;
c’est-à-dire si
v = 0 ou bien+k : g = (7.15)
4 2
Une telle con…guration correspond à une singularité qu’il faudra tenter d’éviter.
Nous avons a¤aire à deux systèmes monovariables découplés. Notons w = (w1 ; w2 )T la consigne pour y.
Nous noterons parfois w = (^v ; ^) pour rappeler que w1 et w2 sont les consignes correspondant à l’angle
d’ouverture de la voile et au cap. Choisissons le régulateur donné de type PID (en fait, il faudrait dire PID2,
car la dérivée deuxième est utilisée) donné par
8
>
< v1 = 0 (w1 y1 ) + w_ 1 ;
RL : v2 = 0 (w2 y 2 ) + 1 (w _ 2 y_ 2 )
>
: ...
+ 2 (w •2 y•2 ) + w 2 :
132 Commande par espace d’état
s+ 0 = s + 1;
s3 + 2s
2
+ 1s + 0 = (s + 1)3 = s3 + 3s2 + 3s + 1;
soit
0 = 0 = 1; 1 = 3; 2 = 3:
En supposant la consigne w constante, les équations d’état du régulateur par retour d’état pour notre
système non linéaire sont données par
! !
1 w1 v
u=A (x) b(x) : (7.16)
w2 3_ 3•
Or _ et • sont des fonctions analytiques de l’état x. En e¤et, d’après les équations (7.11) et (7.12), on a
_ = !;
• = (` rv cos v )fv rg cos g fg !
J :
u = r (x; w) = r x;^v ; ^ :
Le régulateur r x;^v ; ^ développé dans la section précédente permet de stabiliser notre bateau pour un
cap ^ et une ouverture de voile ^v désirés. Or, ce qui nous intéresse est de stabiliser le bateau autour d’une
zone centrée en un point de coordonnées (xd ; yd ). En manœuvrant le gouvernail et l’angle d’ouverture de la
voile, il est possible de rester autour du point désiré.
Nous allons ici proposer un automate capable de générer les consignes w = ^v ; ^ pour notre régulateur
a…n que notre bateau se trouve comme attiré par le point (xd ; yd ) : Contrairement aux autres approches de
régulation classiques, l’élaboration de cet automate sera guidée par une connaissance empirique du pilotage
d’un voilier, plus que par une prise en compte des équations d’états du système.
Ainsi, nous allons considérer quatre caps privilégiés (à 2k près) donnés par
^1 = 7 ; ^2 = 5 ; ^3 = et ^4 = 5 ;
4 4 6 6
comme illustré par la …gure. A chacun de ces caps désirés, nous associons un état discret q 2 f1; 2; 3; 4g.
Commande non linéaire des systèmes non linéaires 133
Les caps ^1 et ^2 correspondent à une situation de vent arrière alors que ^3 et ^4 correspondent à une
situation de remonté au vent. Le passage d’un cap désiré à l’autre se fait suivant la position relative du
bateau par rapport au point désiré. Par exemple, si q = 3 et que x xd > a (où a est une constante de
tolérance positive), alors, q prend la valeur 4 (qui fera décroître x). L’automate qui réalise ces changements
d’état apparaît dans le bloc du bas de la …gure ci-dessous.
134 Commande par espace d’état
A…n d’éviter le passage brutal d’un cap désiré à l’autre, ce qui risque de déstabiliser notre commande
r x;^v ; ^ (qui suppose une consigne constante), il convient d’imposer un changement en douceur, c’est le
rôle des systèmes du premier ordre situés dans les cases de l’automate qui e¤ectuent un …ltrage des hautes
fréquences. La fonction
!
^ ^i
arg ^; ^i = ^i + 2 :round ;
2
renvoie l’angle égal à ^i (à 2k près), qui est le plus proche de ^. Son rôle est d’éviter les sauts brutaux
dans la génération de ^ et le phénomène de déroulement, qui tend à faire tourner inutilement les éléments
(voile et bateau) a…n de respecter scrupuleusement les égalités entre les angles alors qu’une égalité à 2 -près
su¢ t. La sortie de notre automate est le cap désiré ^, il nous reste alors à trouver l’ouverture de la voile
^v à demander à notre régulateur. L’angle d’ouverture de la voile ne servant qu’à la propulsion du bateau,
il est concevable de la rendre fonction du cap. Une fonction assez naturelle est celle donnée par la …gure
ci-dessous.
Commande non linéaire des systèmes non linéaires 135
^v = ' ( ) = :‡oor 1
+ +
2 4 4 2
où ‡oor est la fonction qui renvoie la partie entière d’un réel. Notons que lorsque nous demandons au bateau
d’aller face au vent (^ ' =2 + k ), l’angle d’ouverture demandé ^v sera nul ('(^) ' 0) alors qu’en vent
arrière (^ ' =2 + 2k ), la voile sera grand ouverte ('(^) ' =2 + k ). Le schéma complet du régulateur
ainsi conçu est représenté en pointillé sur la …gure ci-dessous. Le système bouclé est admet pour entrées xd
et yd et le régulateur se charge d’amener le bateau autour du point désiré et de l’y maintenir.
La …gure ci-dessous illustre le fait que pour toute condition initiale, le bateau régulé semble converger vers
un cycle limite unique, centré sur le point de consigne.
Identi…cation
Le but de l’identi…cation est d’estimer, avec une certaine précision, des quantités non-mesurées à partir
d’autres quantités qui elles, sont mesurées. Dans le cas particulier, où la quantité à estimer est le vecteur
d’état d’un système linéaire invariant, les observateurs d’état vus au chapitre 5 peuvent être considérés
comme un outil e¢ cace pour l’identi…cation. Dans ce chapitre, nous allons présenter quelques notions de
base sur l’estimation, dans le but d’introduire le …ltrage de Kalman. Rapidement, le …ltrage de Kalman
peut être vu comme un observateur d’état pour des systèmes linéaires à coe¢ cients variant dans le temps.
Mais, contrairement aux observateurs plus classiques utilisant une approche par placement de pôles (voir le
chapitre 5), le …ltrage de Kalman utilise les propriétés probabilistes des signaux.
8.1.1 Dé…nition
f (x) = xT Q x + Lx + b;
où Q est une matrice symétrique. Cette dé…nition est équivalente à dire que f (x) est une combinaison
linéaire d’une constante b, des xi , de leur carré x2i et des produits croisés xi xj où i 6= j. Par exemple, la
fonction f (x1 ; x2 ) = 2x21 6x1 x2 + x22 2x1 + x2 + 1 est une fonction quadratique. On a
! ! !
2 3 x1 x1
f (x) = (x1 x2 ) + ( 2 1) + 1: (8.1)
3 1 x2 x2
Nous allons montrer dans la suite que la fonction dérivée de f au point x est une fonction a¢ ne. Pour notre
exemple, la dérivée de f au point x est donnée par
df @f @f
(x) = (x) (x)
dx @x1 @x2
137
138 Commande par espace d’état
@f @f
avec @x1 (x) = 4x1 6x2 2 et @x2 (x) = 6x1 + 2x2 + 1; c’est-à-dire
df
(x) = (4x1 6x2 2 ; 6x1 + 2x2 + 1) :
dx
Il s’agit d’une fonction a¢ ne en x. La fonction, x 7! xT Qx qui compose f (x) ne comporte que des termes
en xi xj et en x2i . Une telle fonction est appelée forme quadratique.
f (x + x) = (x + x)T :Q: (x + x)
= xT :Q:x + xT :Q: x + xT :Q:x + xT :Q: x (8.3)
= xT :Q:x + 2xT :Q: x + o (jj xjj) ;
car Q est symétrique et xT :Q: x = o (jj xjj). Par unicité du développement de Taylor, et d’après les
formules (8.2) et (8.3), nous avons
df @f @f
(x) = (x); : : : ; (x) = 2xT Q:
dx @x1 @xn
Par exemple, la dérivée de la fonction quadratique (8.1) est donnée par
!
2 3
2 (x1 x2 ) + ( 2 1) = (4x1 6x2 2 ; 6x1 + 2x2 + 1) :
3 1
Ce sont les valeurs propres de Q: Les valeurs propres sont toutes réelles et les vecteurs propres sont tous
deux à deux orthogonaux. Les courbes de niveaux d’une fonction quadratique f (x) = sont de la forme
xT Qx + Lx = b
et sont appelées quadriques. Ce sont des ellipsoïdes si toutes les valeurs propres sont de même signe ou des
hyperboloïdes si elles ont des signes distincts. Si toutes les valeurs de Q sont positive, on dit que la forme
quadratique xT Qx est positive. Si elle sont toutes non-nulles, on dit que la forme quadratique est dé…nie.
Si elles sont toutes strictement positives, on dira que la forme quadratique est dé…nie positive. La fonction
quadratique f admet un et un seul minimum si et seulement si sa forme quadratique associée est dé…nie
positive.
Commande non linéaire des systèmes non linéaires 139
Estimer, c’est obtenir un ordre de grandeur sur certaines quantités d’un système à partir de mesures d’autres
quantités de ce même système. Le problème d’estimation que nous allons considérer dans ce chapitre est le
suivant. Soit un système sur lequel on a e¤ectué des mesures y = (y1 ; : : : ; yp ) et un modèle M(p) dépendant
d’un vecteur de paramètres p. Il nous faut estimer p tel que les sorties ym (p) générées par M(p) ressemblent
le plus possible à y.
ym (p) = Mp:
140 Commande par espace d’état
Le modèle est alors quali…é de linéaire par rapport aux paramètres. Nous voudrions avoir
ym (p) = y
mais cela n’est généralement pas possible à cause de la présence du bruit et du fait que le nombre de
mesures est généralement supérieur au nombre de paramètres (c’est-à-dire dim(y) > dim(p)). Nous allons
donc chercher le meilleur p, c’est-à-dire, celui qui minimise le critère, dit des moindres-carrés,
Nous avons
On cherche à trouver une parabole p1 t2 + p2 t + p3 qui passe par n points donnés dans le tableau suivant
t 3 1 0 2 3 6
y 17 3 1 5 11 46
p
Pour obtenir ces mesures, nous avons pris p1 = 2; p2 = 1; p3 = 1, pour en déduire les mesures non
bruitées
y = (16:72; 3:41; 1; 4:65; 10:73; 45:91)T :
Ensuite, nous avons tronqué à l’entier le plus proche. La ligne d’instruction sous Scilab est donnée par
Le vecteur des mesures est donc y = (17 3 1 5 11 46)T et le vecteur des paramètres est p = (p1 p2 p3 )T : La
sortie modèle est
0 1 0 1
ym,1 (p) = 9p1 3p2 + p3 9 3 1
B ym,2 (p) = p1 p2 + p3 C B 1 1 1 C 0 1
B C B C p
B y (p) = 0p 0p2 + p3 C B 1 C
1
B m,3 1 C B 0 0 CB C
ym (p) = B C=B C @ p2 A
B ym,4 (p) = 4p1 + 2p2 + p3 C B 4 2 1 C
B C B C p3
@ ym,5 (p) = 9p1 + 3p2 + p3 A @ 9 3 1 A
ym,6 (p) = 36p1 + 6p2 + p3 36 6 1
Commande non linéaire des systèmes non linéaires 141
^ = ym (^
y p = (16:76 ; 3:46 ; 1:06 ; 4:76 ; 10:84 ; 46:11)T :
p) = M^
La vitesse angulaire d’un moteur à courant continu en régime permanent dépend linéairement de la tension
d’alimentation U et du couple résistant Tr :
= p1 U + p2 T r :
On a
ym (p) = M:p
avec
0 1 0 1
40 0 490
B 100 0 C B 1210 C
B C ! B C
B 100 20 C p1 B 1090 C
B C B C
M=B C; p = et y = B C:
B 130 20 C p2 B 1410 C
B C B C
@ 150 10 A @ 1730 A
150 30 1630
Donc
!
1 11:98
^ = MT M
p MT y = :
6:02
Nous pouvons alors en déduire la vitesse angulaire du moteur pour d’autres valeurs de U et Tr . Par exemple
pour U = 200 V et Tr = 10 Nm, on a
!
^= 11:98
U Tr ^=
p 200 10 = 2335tr= min
6:02
142 Commande par espace d’état
Si y est le vecteur des mesures et si ym (p) est la sortie générée par le modèle, alors l’estimée au sens des
moindres carrés est dé…ni par
^ = arg minn jjym (p)
p yjj2 :
p2R
Lorsque ym (p) est linéaire par rapport à p, c’est-à-dire ym (p) = Mp alors le vecteur des paramètres
1
p ^ = MT M
^ estimé au sens des moindres-carrés est p My et le vecteur des mesures …ltrées est y ^ =
T 1
M M M My. En général, et même lorsque ym (p) est non-linéaire, on peut faire l’interprétation géo-
métrique suivante :
Lorsque ym (p) est non-linéaire, on peut utiliser un algorithme d’optimisation local pour espérer obtenir p
^.
La table ci-dessous en propose une version simple, qui est censé converger vers un optimum local du critère
j(p) = jjym (p) yjj2 . La quantité p représente un petit vecteur tiré aléatoirement dans Rn .
input p
1 j + = j(p);
2 q = p + p;
3 if j (q) < j + then {p = q; j + = j(q)} ;
4 goto 2.
On considère le robot de la …gure ci-dessous. Ce robot possède 8 capteurs capables de lui donner les distances
des premiers obstacles suivant les directions d’angles k4 , k 2 f0; : : : ; 7g. Les obstacles sont supposés être n
segments de droites [ai bi ] ; i = 1; : : : ; n, où les ai et bi , exprimés dans le repère de la pièce sont connus du
robot. Les 8 distances sont rangées dans le vecteur des mesures y. Proposons nous de retrouver la position
et l’orientation du robot à partir du vecteur des distances mesurées y.
144 Commande par espace d’état
Pour cela, il nous faut, tout d’abord, développer un simulateur ym (p) qui nous génère nos 8 distances à
partir de la connaissance de p. Ce simulateur peut être représenté par un algorithme ayant pour entrée le
vecteur p = (x; y; )T et pour sortie le vecteur formé des 8 distances fournies par les capteurs. Dans la suite,
les coordonnées du centre m du robot seront notées (x; y) et ! u sera un vecteur unitaire représentant la
direction du faisceau. Pour le kième capteur, l’expression du vecteur directeur ! u du faisceau lancé par le
capteur est
!
! cos k4 +
u = ; k 2 f0; : : : ; 7g:
sin k4 +
En s’inspirant de l’algorithme de calcul de distance élaboré au paragraphe 6.6.1 de la page 111, nous obtenons
le simulateur ym (p) donné par la table ci-dessous
input : (x; y; )
for i = 1 to 8
!
u := cos (i 1)
+ ; sin (i 1)
+ ;
4 4
T
m := (x y) ; `i := 1;
for j = 1 to n
det(a m;bj aj ) (8.5)
:= det(j! u ;bj aj )
;
!
if det aj m; u : det bj m; !
u 0 and ( 0)
then `i := min (`i ; ) ;
next j
next i;
return (`1 ; : : : ; `8 ) ;
8.4.1 Principe
Le point p
^ k qui minimise le critère
k 1
X 2
j(p) , jjMk p 2
yk jj = mT
i p y(i)
i=0
La suite du problème consiste à trouver des équations récursives pour p ^ k , c’est-à-dire des équations d’état
de la forme
(
x (k + 1) = f (x (k) ; mk ; y (k))
p
^k = g(x (k))
Posons k = MT T
k Mk et vk = Mk yk . Notons que, contrairement aux quantités Mk et yk , les dimensions de
k et vk ne dépendent pas de k. Nous avons
8 !
>
> Mk
>
> k+1 = MT mk = MT T
k Mk + mk mk
>
> k
mT
< k !
> yk
>
> vk+1 = MT
k mk = MT
k yk + mk y(k)
>
> y(k)
>
: 1
p
^k = k vk :
On peut donc exprimer k+1 , vk+1 en fonction de k et vk et mk et y(k): Ainsi, les équations d’état pour
générer les p
^ k sont données par
8
>
< (i) k+1 = k + mk mk
T
Le vecteur d’état x (k) contient les coe¢ cients de la matrice k et ceux du vecteur vk . Les entrées de ce
système sont mk et y(k).
146 Commande par espace d’état
^ k+1 = p
p ^ k + kk y(k) mT
kp^k
où apparaissent
la prédiction y^(k) = mTkp
^k
l’erreur de prédiction e(k) = y(k) mT kp
^k
T 1
le gain de correction ou gain de Kalman kk = k + mk mk mk
Pour cela, on procède comme suit :
(8.7,iii) 1
p
^ k+1 = k+1 vk+1
(8.7,ii,iii) 1
= k+1 ( kp
^k + mk y(k))
1
= ^k +
p k+1 (( k k+1 ) p
^ k + mk y(k))
(8.7,i) 1
= ^k +
p k+1 : m k mT
kp^ k + mk y(k)
(8.7,i) 1
= ^k +
p k + mk mT
k mk y(k) mT
kp^k :
Les équations d’état pour notre estimateur s’écrivent alors
(
T
k+1 = k + mk m k
T 1
^ k+1 = p
p ^k + k + mk mk mk y(k) mT
kp^k
k , k
1
Lemme d’inversion matricielle : Si A est une matrice carrée n n et inversible, B est une matrice
colonne n 1 et C une matrice ligne 1 n,
1 1 1 1 1 1
(A + BC) =A A B 1 + CA B CA :
1 1 1 1 1 1 1
k+1 = k + mk mT
k = k k mk 1 + mT
k k mk mT
k k
1 1
1 mk mT
k k k
= k 1
1 + mT
k k m k
soit
8 T
< k mk mk k
k+1 = k 1+mT k k mk
T
: p k mk mk k
^ k+1 = p
^k + k T
1+mk k mk
mk y(k) mT
kp^k :
En remarquant que
T
k mk mk k k mk
k mk = ;
1 + mTk k mk 1+ mTk k mk
on en déduit que
8 T
< k+1 = k
k mk mk k
T
1+mk k mk
: p k mk
^ k+1 = p
^k + 1+mT
y(k) mT
kp^k
k k mk
Ce sont des équations d’état où les entrées est y(k) et mk et la sortie est p
^ k . Remarquons que, tout comme
k la matrice k est symétrique.
(p ^ k )T
p k
1
(p ^k )
p 1:
où " est très proche de zéro (par exemple 0:001). Ce qui revient à supposer que p est dans une sphère
centrée en zéro et de rayon très grand ( 1" ).
Remarque (facteur d’oubli) : Pour donner plus d’importance aux dernières mesures et ainsi prendre en
compte une légère évolution de p dans le temps, on multiplie chacune des équations mi p = yi par k i ,
où est un coe¢ cient d’oubli légèrement inférieur à 1 (par exemple 0.99). Ainsi, les équations de récurence
148 Commande par espace d’état
(8.7) deviennent
0 1
k
mT
0
B . C k
B .. C X 2
= k
m0 mk mk B C= k i
mi m T
k+1 1 B C i
@ mT k 1 A i=0
mT k
k 1
!
X 2
= 2 (k 1) i
mi mT
i + mk m T
k =
2
k + mk mT
k
i=0
0 1
k
y(0)
B . C k
B .. C X 2
vk+1 = k
m0 mk mk B C= k i
mi y(i)
1 B C
@ y(k 1) A i=0
y(k)
k 1
X
2 (k 1) i 2
= mi y(i) + mk y(k) = vk + mk y(k):
i=1
8 T
< 2 k mk mk k
k+1 = k 2
+mTk k mk
: p
^ k+1 = p
^k + 2 k mk
y(k) mT
1+mT kp^k
k k mk
Soit x un vecteur aléatoire, on dé…nit la matrice de covariance de x, la matrice dont les éléments sont
ij = E ((xi xi ) (xj xj )) :
Par exemple, traçons sous Scilab un nuage de point gaussien de matrice de covariance égale à l’identité.
Cela se fait par la suite d’instruction suivante :
On obtient e¤ectivement une matrice proche de l’identité. Fabriquons maintenant une nouveau vecteur
aléatoire y à l’aide de x à l’aide de l’instruction suivante :
La distribution de probabilité d’un vecteur aléatoire gaussien x est entièrement caractérisée par son espérance
x et sa matrice de covariance x :
n 1 1
px (x) = (2 ) 2 det( x) 2 exp (x x)T x
1
(x x) . (8.8)
2
150 Commande par espace d’état
Le caractère gaussien d’une variable aléatoire se conserve par toute transformation linéaire. Considérons, par
exemple, le vecteur y = Ax, où A est une matrice carrée inversible. Un élément de volume dx se transforme,
par l’application linéaire associée à la matrice A, en un élément de volume dy = j det(A)jdx. On a
Donc
px (A 1 y) (8.8) 1 n 1 1 1 1 T 1 1 1
py (y) = = (2 ) 2 det( x) 2 exp A y A y x A y A y
j det(A)j j det(A)j 2
Or
T 1
A 1
y A 1
y x
1
A 1
y A 1
y = (y y)T A xA
T
(y y)
et
1
det( x ) 2 1 1 1 1
= j det(A)j 2 : det( x) 2 :j det(AT )j 2 = det A xA
T 2
j det(A)j
donc
n 1 1 1
py (y) = (2 ) 2 det(A xA
T
) 2 exp (y y)T A xA
T
(y y)
2
Ainsi, y est aussi une variable aléatoire gaussienne. Son espérance est y = Ax et sa matrice de covariance
est y = A: x :AT :
Cette propriété se généralise même lorsque A est rectangulaire, de rang quelconque. Soient x1 et x2 sont
deux vecteurs aléatoires gaussiens indépendants, d’espérances x1 et x2 et de matrices de covariance x1 et
T T T est aussi gaussien d’espérance x = xT ; xT T et de matrice de covariance
x2 . Le vecteur x = x1 ; x2 1 2
!
x1 0
x = :
0 x2
Le vecteur aléatoire !
x1
y = A1 x1 + A2 x2 = A1 A2
x2
satisfait !
x1
y= A1 A2 = A1 x1 + A2 x2 :
x2
et ! !
x1 0 AT
1 T T
y = A1 A2 = A1 x1 A 1 + A2 x2 A 2 .
0 x2 AT
2
Les opérations linéaires sur les vecteurs aléatoires gaussiens sont donc relativement aisées, du fait que la
gaussianité se conserve par transformation linéaire. Elles se font par calcul direct sur les espérances et sur
les matrices de covariance, qui su¢ sent à caractériser la distribution gaussienne considérée. Ces opérations
se retrouvent dans le …ltre de Kalman auquel nous allons nous intéressés.
Commande non linéaire des systèmes non linéaires 151
avec
T T 1
K(k) = x (k):C (k): z (k) + C(k) x (k)C (k) :
Le vecteur x
^(k) représente une estimée pour le vecteur x(k). Les matrices x (k), v (k); z (k) sont respec-
tivement les matrices de covariance sur l’erreur d’estimation x
^(k) x(k), sur le bruit d’état v(k) et sur le
bruit de mesure z(k). La matrice K(k) est appelée gain de Kalman.
Initialisation du …ltre : En général, nous n’avons aucune idée de la valeur de x(0). Ainsi, il semble assez
logique de poser
0 1
1
0 0
B " 1 C
B 0 " C
^(0) = (0; 0; : : : ; 0)T et x (0) = B
x B ..
C
C
@ 0 . 0 A
1
0 "
où " est très proche de zéro (par exemple 0:001). Ce qui revient à supposer que x(0) est dans une sphère
centrée en zéro et de rayon très grand ( 1" ).
T T 1
x =A x xC z +C xC C x AT + v
L’équation ci-dessus est une équation de Ricatti que l’on peut résoudre une fois pour toute en lançant
l’itération
T T 1
x (k + 1) = A x (k) x (k)C z +C x (k)C C x (k) AT + v
jusqu’à atteindre l’équilibre. Ainsi, dans le cas stationnaire, on peut faire l’économie du calcul en temps réel
de x (k).
Lien avec les moindres-carrés sans oubli : Remarquons que si A(k) = I, B = 0, x ^ , C(k) = mT
^=p k,
v (k) = 0 et z (k) = 1 on retrouve les moindres-carrés récursifs sans oublis :
8
< p 1
^ (k + 1) = p
^ (k) + x (k)mk 1 + mT
k x (k)mk y(k) mT
kp^ (k)
: 1
x (k + 1) = x (k) x (k)mk 1 + mT
k x (k)mk mT
k x (k) + v (k)
152 Commande par espace d’état
8.6.1 Principe
Soit un système d’entrée u(k) et de sortie y(k) d’ordre 3 décrit par l’équation de récurrence
On cherche à estimer en temps réel les coe¢ cients de la fonction de transfert du système, sachant que ceux-ci
évoluent doucement dans le temps (voir …gure ci-dessous).
où v(k) est un bruit blanc de matrice de covariance v (souvent, on la prend diagonale). On peut construire
un observateur de Kalman basé sur le système
(
p(k + 1) = p(k) + v(k)
y(k) = C(k)p(k) + z(k)
avec
C(k) = ( y(k 1); y(k 2); y(k 3); u(k 1); u(k 2); u(k 3))
et où z(k) est le bruit de mesure de variance z . Ce …ltre est donné par
(
^ (k + 1) = p
p ^ (k) + K(k) (y(k) C(k)^ p(k))
x (k + 1) = x (k) K(k)C(k) x (k) + v
avec
T
T T 1 x (k)C (k)
K(k) = x (k):C (k): z + C(k) x (k)C (k) = T
:
z + C(k) x (k)C (k)
puisque z + C(k) x (k)CT (k) est un scalaire. Il ne nous reste qu’à trouver une relation récursive qui nous
génère les C(k). Puisque
C(k + 1) = ( y(k); y(k 1); y(k 2); u(k); u(k 1); u(k 2))
Commande non linéaire des systèmes non linéaires 153
des équations récursives pour générer les C(k) sont donc données par
0 1
0 1 0 0 0 0
B 0 0 1 0 0 0 C
B C !
B 0 0 0 0 0 0 C 0 0 0 1 0 0
B C
C(k + 1) = C(k) B C + (u(k) y(k))
B 0 0 0 0 0 0 C 1 0 0 0 0 0
B C
@ 0 0 0 0 1 0 A
0 0 0 0 0 1
que nous noterons
C(k + 1) = C(k)N + (u(k) y(k)) M
Les équations d’état d’un identi…cateur basé sur le …ltre de Kalman sont donc :
8 T
>
> ^ (k + 1) = p x (k)C (k)
^ (k) + z +C(k) (y(k) C(k)^
p(k))
< p T
x (k)C (k)
T
x (k)C (k)
> x (k + 1) = x (k) +C(k) (k)CT (k) C(k) x (k) + v
>
: C(k + 1)
z x
= C(k)N + (u(k) y(k)) M
Nous proposons dans cette section deux applications possibles de l’estimateur de fonction de transfert.
Commande adaptative : Considérons un système d’ordre connu mais de fonction de tranfert inconnu.
Le schéma décrit sur la …gure ci-dessous propose un régulateur permettant de donner au système le com-
portement désiré. L’estimateur de fonction de transfert utilise les moindres-carrés récursifs pour estimer la
fonction de transfert du processus. Ensuite, le module utilise cette la fonction de transfert estimée pour cal-
culer un régulateur (par exemple par la méthode du placement de pôle). Il vient alors modi…er les équations
d’état du système du régulateur classique.
Surveillance : A…n de surveiller qu’un procesus se comporte bien, on identi…e en temps réel sa fonction
de transfert, par un …ltre de Kalman. Si une panne survient, si un élément casse,: : : alors le système ne
comportera di¤érement. Cela se traduira par une nouvelle fonction de transfert pour notre système et donc
un saut dans le vecteur des paramètres estimé. Un détecteur de saut pourra alors nous signaler un tel
changement de comportement du système.
154 Commande par espace d’état
On considère un robot sous marin (dont les équations d’état ne sont pas nécessairement connues) se déplaçant
dans une piscine rectangulaire de demie longueur Ry et de demie largeur Ry . Un sonar positionné juste en
dessous du robot tourne avec une vitesse angulaire constante. Le but de ce paragraphe est de proposer
l’utilisation d’un …ltre de Kalman pour localiser le robot. La profondeur z s’obtient facilement avec un
capteur de pression et nous supposerons donc cette grandeur connue. Le robot est supposé être lesté de
telle façon que les angles de roulis et de tangage puissent être considérés comme nuls. L’angle de cap est
mesuré par une boussole. Dans notre contexte, se localiser signi…e donc estimer les coordonnées (x; y) du
robot. L’origine du repère sera pris au centre de la piscine. Pour cette localisation, nous supposons que nous
mesurons l’angle du sonar relativement au corps du robot, l’angle par une boussole et les accélérations
tangentielles aT et normales aN par l’intermédiaire d’accéléromètres. Le sonar nous renvoie toutes les 0:1s
la longueur ` du faiseau sonar. La …gure qui suit représente la longueur `(t) du faiseau sonar, obtenue par
simulation, lorsque le sonar e¤ectue sept tours sur lui-même, alors que le robot se déplace.
Commande non linéaire des systèmes non linéaires 155
A partir du signal récupéré par le sonar, et dont la …gure ci-dessus donne une allure dans une situation sans
bruit, nous allons devoir détecter les minimums locaux, qui correspondent à la situation où le sonar pointe
perpendiculairement à un des quatre murs de la piscine. Lorsque le sonar pointe sur un mur, la distance `
renvoyée satisfait
a = `: cos ; (8.9)
où est l’angle entre la normale au mur et la faiseau sonar et a est la distance entre le sonar et le mur.
Nous allons supposer que le robot est immobile et que seul le sonar tourne (cela revient à supposer que les
vitesses tangentielle v et angulaires _ du robot sont négligeables devant la vitesse de rotation _ du sonar).
Considérons tout d’abord, la situation où à l’instant t, le sonar est normal au mur. Si est un réel positif
su¢ samment petit, c’est-à-dire tel que à l’instant t le sonar pointe toujours sur le même mur, nous
devrions avoir, d’après la relation (8.9),
a = `(t ): cos ( _ )
Rappelons que la vitesse de rotation du sonar _ est supposée connue et constante. Prenons = k ; k 2
f0; 1; 2; : : : ; N 1g; où est la durée séparant deux pings du sonar et N est un entier tel que, à l’instant
t N , le sonar pointe nécessairement sur le mur qui se trouve orthogonal au faiseau sonar à l’instant t.
Posons
a
~k = `(t k ): cos ( k _ ) :
La quantité a
~k devrait correspondre à la distance a entre le robot et le mur pointé. Mais du fait de la présence
d’un bruit de mesure, il est préférable d’obtenir une estimation a ^ de la distance a par une moyenne :
N 1
1 X
a
^= a
~k :
N
k=0
156 Commande par espace d’état
Nous pouvons véri…er le fait que le sonar pointe bien perpendiculairement à un mur en véri…ant que la
variance des a
~k est faible, c’est-à-dire que
N 1
1 X
(~
ak ^)2 < ";
a
N
k=1
Or, en pratique il y a beaucoup de données aberrantes. Il nous faut donc modi…er notre méthode. Une ap-
proche plus robuste que celle présentée précédemment, consiste à calculer la médiane plutot que la moyenne.
Pour cela, il faut trier les a
~k par ordre croissant. Ensuite, on prend le milieu a de la liste et on supprime de
la liste les éléments de cette liste qui sont les plus éloignés de a. On en supprime par exemple une moitié.
Il sont facile à trouver car ils se trouvent soit en début, soit en …n de liste. On fait alors la moyenne des
éléments restant pour obtenir a ^. On fait alors de test de la variance sur les éléments restant pour véri…er
que le sonar est normal à un mur. Notons que si le robot dispose d’une boussole …able, le test de la variance
devient inutile. En e¤et, la boussole nous donne et l’angle est connu, ce qui nous permet de savoir si le
sonar pointe ou non suivant la normale à un mur et aussi de quel mur il s’agit.
8.7.3 Localisation
Nous supposons maintenant que nous ayons recueillies (soit par simulation, soit par une expérience réelle)
pour chaque t 2 [0; tmax ], les accélérations (aT ; aN ), l’angle de cap et l’angle du sonar . Nous allons
maintenant chercher à nous localiser et à estimer notre vitesse à l’aide d’un …ltre de Kalman. Nous allons
modéliser les relations di¤érentielles entre les mesures par les équations d’état suivantes
8
>
> x_ = vx
>
< y_ = v
y
(8.10)
>
> v_ x = aT cos aN sin
>
:
v_ y = aT sin + aN cos
où vx = x,
_ vy = y._ Pour les deux dernières équations il su¢ t de remarquer que l’accélération absolue est
obtenue par la relation ! ! !
x• cos sin aT
= :
y• sin cos aN
A…n de pouvoir utiliser un …ltre de Kalman il nous faut tout d’abord discrétiser le temps. Remarquons que
lorsque + est un multiple de 2 , le faiseau sonar est orienté face à un des quatres murs de la piscine (car
cette dernière est supposée rectangulaire). Dans un tel cas, la longueur mesurée peut nous permettre de
calculer soit x soit y. Pour notre problème, le temps k discret s’incrémente chaque fois que le faiseau sonar
pointe face à un mur de la piscine, c’est-à-dire que k = E( + =2 ), où E désigne la partie entière d’un nombre
réel. Une discrétisation par Euler de (8.10) se traduit par
8
>
> x(k + 1) = x(k) + vx (k):T (k)
>
< y(k + 1) = y(k) + v (k):T (k)
y
>
> vx (k + 1) = vx (k) + (aT (k) cos (k) aN (k) sin (k)) :T (k)
>
:
vy (k + 1) = vy (k) + (aT (k) sin (k) + aN (k) cos (k)) :T (k)
Commande non linéaire des systèmes non linéaires 157
où T (k) est le temps qui s’est écoulé entre deux incrémentations successives de k. Sous forme matricielle,
ces équations d’état deviennent
0 1 0 1
1 0 T (k) 0 0 0
B 0 1 0 T (k) C B 0 0 C
B C B C
x(k + 1) = B C x(k) + B C u(k)
@ 0 0 1 0 A @ T (k) sin (k) T (k) cos (k) A (8.11)
0 0 0 1 T (k) cos (k) T (k) sin (k)
r(k) = C(k):x(k);
avec
– Cas 0, mur droit ( (k)+ (k) = 2k ). Dans ce cas, nous avons une mesure sur x : r(k) = x(k) ' Rx d(k),
où d(k) est la distance retournée par le sonar. Ainsi, la matrice d’observation sera
C(k) = 1 0 0 0 :
– Cas 1, mur fond ( (k) + (k) = 2k + 2 ). Dans ce cas, nous avons une mesure sur y : r(k) = y(k) '
Ry d(k) et donc
C(k) = 0 1 0 0 :
– Cas 2, mur gauche ( (k) + (k) = 2k + ). Nous avons à nouveau une mesure sur x : r(k) = x(k) '
Rx + d(k). La matrice d’observation sera
C(k) = 1 0 0 0 :
– Cas 3, mur devant ( (k) + (k) = 2k + 32 ). Nous avons une mesure sur y : r(k) = y(k) ' Ry + d(k).
Ainsi
C(k) = 0 1 0 0 :
Si nous connaissons + , pour savoir dans quel cas i nous sommes, il faut résoudre
i
9k 2 N; (k) + (k) = 2k +
2
c’est-à-dire
2
9k 2 N; ( (k) + (k)) = i + 4k
Donc, pour avoir le cas i le plus probable, on calcule l’entier le plus proche de la quantité 2 ( (k) + (k)),
et on regarde le reste de la division euclidienne de cet entier par 4. Sous Scilab, cela se fait pas le calcul
suivant
i=modulo(round((thetak+alphak)*2/%pi),4)
Notons que le système (8.11) n’a pas pour but de reproduire le comportement dynamique du système au
niveau de la commande, mais plutôt de permettre l’utilisation d’un …ltre de Kalman, dans le but d’estimer
158 Commande par espace d’état
la position et la vitesse du robot. Pour cela, nous nous sommes arrangés pour avoir un système discret décrit
par des équations d’état linéaires. On a donc bien un système de la forme
(
x(k + 1) = A(k)x(k) + B(k)u(k) + v(k)
r(k) = C(k)x(k) + z(k):
Nous venons de rajouter deux signaux de bruit v et z qui seront supposés blancs, de matrice de covariance
v et z (notons que ici, z est un scalaire). Ces deux matrices ont pour but de modéliser les incertitudes
de modèle et les bruits de mesure. Le …ltre de Kalman est donc décrit par :
(
^(k + 1) = A(k): (^
x x(k) + K(k) (r(k) C(k)^ x(k))) + B(k)u(k)
x (k + 1) = A(k) ( x (k) K(k)C(k) x (k)) AT (k) + v (k);
avec
T T 1
K(k) = x (k):C (k): z (k) + C(k): x (k):C (k) :
Pour chaque k, le …ltre de Kalman nous donne un couple (^ x; ), qui peut s’interpréter comme un ellipsoïde
4
de R contenant l’ensemble des x = (x; y; vx ; vy ) compatibles avec nos mesures. Cet ellipsoïde est donné par
(x ^)T
x 1
(x ^)
x 1:
En pratique, on ne trace qu’une représentation bi-dimensionnelle de cet ellipsoïde. Dans notre exemple, on
pourra se limiter aux composantes de x correspondant au sous-vecteur w = (x; y) qui a une signi…cation
claire dans l’espace de travail du robot. Ainsi, nous considérerons la matrice de covariance Q, sous matrice
de , donnée par !
11 12
Q= :
21 22
ou encore
1
jjQ 2 (w w) jj2 1:
1
En posant s = Q 2 (w w), cette inégalité devient jjsjj2 1. L’ellipsoïde en w peut être considérée comme
1
l’image par l’application a¢ ne w(s) = w + Q 2 s du disque unité. Pour tracer sa frontière, il su¢ t de rentrer
les instructions Scilab suivantes
La …gure qui suit représente le robot à un intant t; le faiseau sonar associé et une suite d’ellipsoïdes de
con…ance engendrés par le …ltre de Kalman. Le grand cercle représente l’ellipsoïde de con…ance initial.
Remarque. Soit b une distribution gaussienne centrée en zéro. On s’intéresse à la probabilité que b appar-
tienne à la sphère de rayon , c’est-à-dire,
P bT b 2
La variable aléatoire z = bT b suit une loi du 2 . Sa densité de probabilité est donnée par
( n
p n1 z 2 1 exp z
2 si z 0
fn (z) = 2 (n=2)
0 sinon
où Z 1
(a) = e t ta 1
dt = (x 1)!
0
est la fonction d’Euler. Par exemple, si on veut que cette probabilité soit égale à 0:95, et pour n = 2, à
l’aide d’une table du 2 , nous pouvons voir qu’il nous faut prendre 2 = 6. Si maintenant, on considère un
vecteur gaussien de matrice de covariance et d’espérance x. On s’intéresse à l’ellipsoïde de con…ance
(x x)T : 1
: (x x) 2
c’est-à-dire
1
jj 2 (x x) jj2 2
:
1
Le vecteur b = 2 (x x) est un vecteur blanc gaussien, donc
P (x x)T : 1
: (x x) 2
= P (jjbjj2 2
):
Pour n = 2, on obtient on obtient que l’ellipsoïde de con…ance qui contient x avec une probabilité de 0.95
est donnée par
(x x)T : 1 : (x x) 6:
160 Commande par espace d’état
Bibliographie
161
162 Commande par espace d’état
Glossaire
bdiag(A) : renvoie une matrice diagonale contenant les valeurs propres de A et une matrice dont les colonnes
sont les vecteurs propres de A.
champ(x1,x2,y1,y2) : trace le champ de vecteur (y1,y2) sur la grille formée par les vecteurs x1 et x2.
deff(’y=f(x1,x2)’,’y=x1+sin(x2)’) : dé…nit une fonction f à l’aide d’une chaîne de caractère.
eye(n,n) : génère une matrice identité de dimension n n.
feval(x1,x2,f) : évalue la fonction f sur la grille formée par les vecteurs x1 et x2. Le résultat renvoyé est
une matrice.
isoview(xmin,xmax,ymin,ymax) : dé…nit les échelles sur la fenêtre graphique.
modulo(a,b) : renvoie le reste de la division euclidienne de l’entier a par l’entier b.
ones(v) : renvoie une matrice de même dimension que v qui ne contient que des 1.
pvm_set_timer() : initialise le chronomètre.
pvm_get_timer() : renvoie la valeur du chronomètre en microsecondes.
poly([a1,a2,...],’s’) : renvoie le polynôme (s a1 )(s a2 ) : : :
round(a) : renvoie l’entier le plus proche du nombre réel a.
seteventhandler() : active ou désactive le gestionnaire d’événements.
xbasc() : e¤ace et initialise la fenêtre graphique.
xdel() : supprime la fenêtre graphique.
xfpoly(x,y) : trace un polygone, tout comme xpoly, mais cette fois, le polygone est plein.
xpoly(x,y) : trace, sur la fenêtre graphique, un polygone dont les abscisses et ordonnées sont rangées dans
les vecteur x et y.
xset(’thickness’,a) : con…gure l’épaisseur a du stylo, lors du dessin graphique.
zeros(v) : renvoie une matrice de même dimension que v qui ne contient que des 0.
163
164 Commande par espace d’état
165
166 Commande par espace d’état
régulateur, 14
régulateur statique, 115, 121
représentation d’état, 13, 17
robot manipulateur, 49
robots à roues, 24, 98
satellite, 19
Scilab, 14
simulation, 13, 39, 50
solution des équations d’état, 57
solution forcée, 58, 59
solution homogène, 58, 59
solution libre, 58
solution transitoire, 58
sortie, 13
stabilité, 59