Commande optimale des systèmes dynamiques
Commande optimale des systèmes dynamiques
J. FRÉDÉRIC BONNANS
DIRECTEUR DE RECHERCHE INRIA
INRIA ET CENTRE DE MATHÉMATIQUES APPLIQUÉES,
ECOLE POLYTECHNIQUE, PALAISEAU
Résumé. L’objet de la commande optimale est l’optimisation de systèmes dynamiques suivant différents
objectifs : atteinte d’une cible en temps ou énergie minimale, maximisation du rendement d’un processus
industriel par exemple. Pour cela on joue à la fois sur des paramètres indépendants du temps et sur
les commandes qui, elles, dépendent du temps. L’article analyse les conditions d’optimalité du premier
et second ordre, et leur résolution par discrétisation temporelle, algorithme de tir, ou programmation
dynamique.
Abstract. The optimal control theory analyzes how to optimize dynamical systems with various crite-
ria : reach a target in minimal time or minimal energy, maximize the efficiency of an industrial process
for instance. This involves the optimization of both time independent parameters, and the control va-
riables that are function of time. The article analyzes the first and second order optimality conditions,
and the ways to solve them, by time discretization, the shooting algorithm, or dynamic programming.
Éclairage historique. La commande optimale est née dans les années 1950, avec les travaux
de l’école de L.S. Pontriaguine [19] et de R. Bellman [5]. Pontriaguine a obtenu des conditions
d’optimalité d’une trajectoire, qui ont permis la résolution explicite ou numérique de nombreux
problèmes, en particulier de trajectoires aériennes ou spatiales. Bellman a énoncé le principe de
programmation dynamique, qui s’applique également dans le cadre discret ou stochastique. Les
approfondissements théoriques et algorithmiques de la théorie, et l’amélioration des moyens de
calcul permettent maintenant de résoudre facilement un grand nombre de problèmes appliqués.
Introduction
On dit qu’un système dynamique est commandé s’il est possible d’agir sur lui par des variables
dépendant du temps, appelées commandes. Illustrons ce concept dans le cas d’un engin spatial, décrit
par des variables de position et vitesse (dans R3 ) h et V , et une masse m > 0, soit 7 variables d’état.
La dynamique est, omettant l’argument temps, ḣ = V , mV̇ = F (h, V ) + u, et ṁ = −c|u|. Ici c est une
constante positive et F (h, V ) correspond aux forces de gravité et (le cas échéant) aérodynamiques. La
commande est la force appliquée, dont on a noté la norme euclidienne par |u|, soumise à une contrainte
du type |u| ≤ U . Étant donné un point initial fixé, on cherche à atteindre une cible (partie de l’espace
d’état) en minimisant un compromis entre temps de parcours et énergie dépensée.
Pour l’implémentation en temps réel d’une commande il est nécessaire de prendre en compte les moyens
d’observations et la reconstitution de l’état, prenant en compte des aspects de traitement du signal et
le choix de l’électronique de commande. Au contraire, dans cet article, nous ne considérons que l’étude
amont, dans laquelle on se donne un cadre déterministe et on calcule donc hors ligne une commande
optimale. L’allure de cette dernière pourra guider la conception de la commande en temps réel.
L’exposé suivra d’abord l’approche de Lagrange et Pontriaguine qui consiste à étudier les variations
d’une trajectoire optimale pour déterminer les propriétés de cette dernière. On analysera les conditions
d’optimalité du premier et second ordre, en lien avec l’algorithme de tir, et avec une attention toute par-
ticulière pour les problèmes avec contraintes sur l’état. Dans une seconde partie, on abordera l’approche
par programmation dynamique et équation de Hamilton-Jacobi-Bellman (HJB). Après avoir introduit
2
la notion de solution de viscosité de l’équation HJB, on présentera les problèmes de temps d’arrêt et de
commande impulsionnelle, et on finira par l’étude des schémas de résolution numérique.
1. Notations
n
PnR l’espace euclidien de dimension n, dont les éléments sont des vecteurs verticaux, de norme
On note
|x| := ( i=1 x2i )1/2 ; son dual Rn∗ est identifié à l’ensemble des vecteurs horizontaux. On note A† la
Pn
transposée d’une matrice A, x · y le produit scalaire de x et y, ainsi que kxks := ( i=1 |xi |s )1/s , pour
s ∈ [1, ∞[, et kxk∞ := maxni=1 |xi |. L’ensemble réduit au zéro de Rn est noté {0}Rn ; l’ensemble des
vecteurs de coordonnées toutes positives (resp. négatives) est appelé cône positif (resp. négatif) de Rn et
noté Rn+ (resp. Rn− ). Ainsi, une partie de Rn définie par un ensemble fini d’égalités et d’inégalités, soit
x ∈ Rn ; gi (x) = 0, i = 1 . . . , n0 ;
(1)
gi (x) ≤ 0, i = n0 + 1 . . . , n0 + n00
peut se réécrire comme g(x) ∈ K[n0 , n00 ], où
00
(2) K[n0 , n00 ] := {0}Rn0 × Rn− .
L’ensemble K[n0 , n00 ] est convexe et fermé. On considèrera plus généralement des contraintes du type
g(x) ∈ K, où K est une partie convexe fermée de Rn . Soit x̄ ∈ K. On appelle (i) cône normal à K en x̄
l’ensemble
(3) NK (x̄) := {y ∈ Rn ; y · (x − x̄) ≤ 0, pour tout x ∈ K},
(ii) cône des directions radiales à K en x̄ l’ensemble
(4) RK (x̄) := {z ∈ Rn ; x̄ + εz ∈ K, pour ε > 0 assez petit},
et (iii) cône tangent à K en x̄ la fermeture notée TK (x̄) de RK (x̄). Quand le bord de K est lisse, NK (x̄)
se réduit à la demi droite engendrée par la normale extérieure en x̄ à K, et TK (x̄) est le demi espace
fermé des directions faisant un produit scalaire négatif avec cette normale. Si K = K[n0 , n00 ], et x̄ ∈ K,
on vérifie facilement que
0 00
(5) NK (x̄) = {y ∈ Rn × Rn+ ; yj x̄j = 0, j = n0 + 1, . . . , n0 + n00 },
00
(6) RK (x̄) = TK (x̄) = {0}Rn0 × {y ∈ Rn ; yj ≤ 0 si x̄j = 0, j = n0 + 1, . . . , n0 + n00 }.
Soient A et B des matrices symétriques de taille n. On note A B si A − B est semi définie positive, et
A B si A − B est définie positive. Si f : Rn → Rp est dérivable on identifie sa jacobienne, notée Df (x)
ou f 0 (x), en x ∈ Rn , à une matrice de taille p × n. Quand p = 1, on note par ∇f (x) := f 0 (x)† le gradient
de f en x ; c’est l’élément de Rn de coordonnées égales aux dérivées partielles de f en x. Si f (x, y) est
fonction de x ∈ Rn et y ∈ Rq on note par Dx f (x, y) ou fx (x, y) sa dérivée partielle, et par ∇x f (x, y) son
gradient partiel. Si f est à valeurs scalaires, fxy (x, y) désigne la matrice de taille n × q de terme général
∂f (x, y)/∂xi ∂yj .
On note la dérivée temporelle d’une fonction f du temps par f˙(t). Notons D(0, T ), où T > 0, l’ensemble
des fonctions de classe C ∞ à support compact sur ]0, T [. Soit f ∈ L1 (0, T ). S’il existe g ∈ L1 (0, T ) tel
que
Z T Z T
g(t)ϕ(t)dt = − f (t)ϕ̇(t)dt, pour tout ϕ ∈ D(0, T ),
0 0
on dit que g est la dérivée faible de f . Si f est de classe C 1 , alors g = f˙ p.p. (presque partout). Soit
s ∈ [0, ∞]. On note W 1,s (0, T ) l’ensemble des fonction de Ls (0, T ) de dérivée faible dans Ls (0, T ), de
norme
kyk1,s := kykLs (0,T ) + kẏkLs (0,T ) .
On sait que W 1,∞ (0, T ) est l’ensemble des fonctions lipschitziennes sur [0, T ]. On définit les espaces de
commande et d’état
(7) Us := Ls (0, T )m , Ys := W 1,s (0, T )n ,
de norme associées notées kuks , kyk1,s (les fonctions du temps telles que commandes et état sont notées
en gras), et on pose U := U∞ , Y := Y∞ , kuk := kuk∞ , kyk := kyk1,∞ .
3
2. Cadre mathématique
2.1. Système dynamique commandé et critère associé. Considérons des systèmes dynamiques
commandés du type
(8) (i) ẏ(t) = f (u(t), y(t)), t ∈ [0, T ]; (ii) y(0) = y 0 .
Les données sont l’horizon ou temps final T > 0, la dynamique f : Rm ×Rn → Rn , supposée lipschitzienne
et de classe C ∞ , et la condition initiale y 0 . On appelle u(t) et y(t) la commande et l’état à l’instant t.
On choisit U et Y comme espaces pour la commande et pour l’état. Alors (8)(i) est une égalité dans
L∞ (0, T )n , donc aussi p.p. Une variante du théorème de Cauchy-Lipschitz permet de vérifier que, pour
tout (u, y 0 ) ∈ U1 × Rn , l’équation d’état (8) a dans Y une solution unique notée y[u, y 0 ] (ou y[u] si y 0
est fixé), et d’après le lemme de Gronwall, on a pour un certain Cf ne dépendant que de la constante de
Lipschitz de f :
ky[u0 , (y 0 )0 ] − y[u, y 0 ]k∞ ≤ Cf ku0 − uk1 + |(y 0 )0 − y 0 | .
(9)
De plus (u, y 0 ) 7→ y[u, y 0 ] est de classe C ∞ : U × Rn → Y (appliquer le théorème des fonctions implicites
(TFI) à l’équation d’état).
On dira que (u, y) ∈ U × Y est une trajectoire si y = y[u, y(0)]. Dans la suite on effectue l’analyse
autour de la trajectoire nominale (ū, ȳ), et on note
(10) f¯(t) := f (ū(t), ȳ(t)), Df¯(t) := f 0 (ū(t), ȳ(t))
la dynamique nominale et sa dérivée, et de même pour les dérivées de tous ordres et les autres fonctions
de (ū(t), ȳ(t)), les dérivées partielles étant notées f¯u (t) par exemple. La dérivée directionnelle de y[u, y 0 ]
au point (ū, ȳ 0 ) dans la direction (v, z 0 ) ∈ U × Rn , notée
(11) z[v, z 0 ] := lim (y[ū + sv, ȳ 0 + sz 0 ] − y[ū, ȳ 0 ])/s
s→0
est solution unique de l’équation d’état linéarisée
(12) (i) ż(t) = Df¯(t)(v̄(t), z(t)), t ∈ [0, T ]; (ii) z(0) = z 0 .
Si z 0 = 0, on notera z[v] := z[v, 0]. On associe au système commandé (8) le critère
Z T
(13) J(u, y) := `(u(t), y(t))dt + ϕ(y(T )),
0
somme d’un coût intégral, d’intégrande `, et d’un coût final ϕ ; on suppose ` et ϕ de classe C ∞ . Alors
J : U × Y → R est aussi de classe C ∞ , de même que le critère réduit
(14) F (u, y 0 ) := J(u, y[u, y 0 ]).
De plus, si (v, z 0 ) ∈ U × Rn , notant z = z[v, z 0 ], la dérivée du critère réduit a pour expression
Z T
(15) 0 0 0
F (ū, ȳ )(v, z ) = `¯0 (t)(v(t), z(t))dt + ϕ0 (ȳ(T ))z(T ).
0
2.2. Différents types de contraintes. Les contraintes mixtes (commande et état) sont du type
(16) g(u(t), y(t)) ∈ Kg , p.p. t ∈ [0, T ],
avec g : R × R → R et Kg partie convexe fermée de Rng . Si g ne dépend que de la première (resp.
m n ng
seconde) variable on parle de contraintes sur la commande (resp. de contraintes (pures) sur l’état). On
peut aussi considérer des contraintes abstraites sur la commande du type
(17) u(t) ∈ Uad , p.p. t ∈ [0, T ],
m
où Uad est un fermé de R , traduisant par exemple l’intervalle de variation admissible de chaque compo-
sante : alors Uad est de la forme Πmi=1 [ai , bi ]. Un autre cas fréquent est celui où Uad est après normalisation
P3
égal à la boule euclidienne fermée B̄ := {u ∈ R3 ; 2
i=1 ui ≤ 1}.
Les contraintes sur l’état final, étant donné ΦF : Rn → RnΦ , sont de la forme
(18) ΦF (y(T )) ∈ KF ,
nΦ
avec KF partie convexe fermée de R . Soit une subdivision de [0, T ] :
(19) 0 = t0 < t1 < · · · < tnpi < tnpi+1 = T.
On considère les contraintes de points intérieurs
(20) Φi (y(ti )) ∈ Ki , i = 1, . . . , npi ,
4
avec Ki partie convexe fermée de Rni . Elles peuvent traduire une contrainte d’interpolation, le passage
d’un véhicule par des points imposés, etc. Dans certains cas les ti sont des variables de décision.
2.3. Un format général. À une trajectoire (u, y) on peut associer un critère incluant un coût aux deux
bouts et des contraintes aux deux bouts :
Z T
IF
(21) J (u, y) := `(u(t), y(t))dt + ϕ(y(0), y(T )),
0
Le problème réduit correspondant est obtenu en exprimant l’état comme fonction de la commande et
de l’état initial, soit y = y[u, y 0 ]. Montrons comment ramener au format (23) plusieurs autres types de
problèmes.
(a) Paramètres de décision. Soit le problème de minimisation par rapport à (u, y, π), où π ∈ Rnπ ,
l’état initial y 0 étant donné, du critère
Z T
(24) J(u, y, π) := `(u(t), y(t), π)dt + ϕ(y(T ), π),
0
m n nπ
avec ici ` : R × R × R → R et ϕ : Rn × Rnπ → R, sous la contrainte
(b) Horizon variable. Si T est une variable de décision, notons τ ∈ [0, 1], le temps normalisé, yN (τ ) :=
y(τ T ) l’état en temps normalisé, et de même pour les autres fonctions du temps. Alors l’équation
d’état (8) se traduit par
Dans cette écriture en temps normalisé, l’horizon apparaît comme un paramètre de décision. On
a vu au point (a) comment se ramener au format (23) en augmentant la taille de l’état.
(c) Points intérieurs (voir (20)). Nous ne donnons que le principe. Reprenant l’analyse qui précède,
on peut ramener chaque intervalle ]ti , ti+1 [ à une durée 1, en introduisant comme paramètre la
durée réelle Ti := ti+1 − ti . On intègre alors en temps fictif τ ∈ [0, 1], en parallèle, les états dans
chaque intervalle, sous contrainte de recollement : la valeur finale de l’état sur un intervalle doit
coïncider avec la valeur initiale de l’état sur l’intervalle suivant.
5
3. Contraintes sur la commande
3.1. Dérivée du critère. Définissons le préhamiltonien H : Rm × Rn × Rn → R par
(30) H(u, y, p) := `(u, y) + p · f (u, y).
On adopte la notation H̄(t) dans l’esprit de (10). À la trajectoire (ū, ȳ) on associe l’état adjoint p̄ ∈ Y,
solution unique de l’équation rétrograde dite équation adjointe
˙
(i) −p̄(t) = ∇y H̄(t), p.p. t ∈ [0, T ];
(31)
(ii) p̄(T ) = ∇ϕ(ȳ(T )).
¯ + f¯y (t)† p̄(t), et donc
Noter que ∇y H̄(t) = ∇y `(t)
Lemme 3.1. Le critère réduit F , défini en (14), a en (ū, ȳ 0 ) une dérivée caractérisée par
Z T
(33) F 0 (ū, ȳ 0 )(v, z 0 ) = p̄(0) · z 0 + ∇u H̄(t) · v(t)dt.
0
6
3.2. Conditions du premier ordre avec contraintes sur la commande. On suppose la condition
initiale fixée, et on ne considère que des contraintes sur la commande du type (17). Le problème s’écrit
donc (supprimant l’état initial comme argument du critère réduit)
(39) Min F (u) soumis à (17).
u
Si la relation précédente est satisfaite on dit aussi que (ū, ȳ, p̄) (ou (ȳ, p̄) s’il n’y a pas d’ambiguité) est
une biextrémale de Pontriaguine. Notons que, si Uad est convexe, une extrémale de Pontriaguine satisfait
la condition du premier ordre (41). Le PMP s’énonce ainsi :
7
Théorème 3.7. Soit ū solution de (39). Posons ȳ := y[ū]. Alors (ū, ȳ) est une extrémale de Pontria-
guine.
La démonstration nécessite le lemme suivant. Dans la suite, M ≥ kūk ; posons
(48) U M := {u ∈ U; kuk ≤ M }.
Lemme 3.8. Il existe cM > 0 tel que, pour tout u ∈ U M , il existe r(u) ∈ R satisfaisant :
Z T
(49) F (u) = F (ū) + (H(u(t), ȳ(t), p̄(t)) − H̄(t))dt + r(u), et |r(u)| ≤ cM ku − ūk21 .
0
on a |∆3 | = O ((ku − ūk1 + kδyk∞ ) kδyk∞ ). On conclut en notant que kδyk∞ = O(ku − ūk1 ).
Démonstration du théorème 3.7. (a) On montre d’abord que, si la trajectoire (u, y) est admissible,
alors
(53) H(ū(t), ȳ(t), p̄(t)) ≤ H(u(t), ȳ(t), p̄(t)) p.p. sur [0, T ].
En effet, dans le cas contraire il existe ε > 0 et une partie mesurable I de [0, T ] de mesure strictement
positive tels que
H(u(t), ȳ(t), p̄(t)) ≤ H(ū(t), ȳ(t), p̄(t)) − ε p.p. sur I.
Soient M ≥ max(kuk, kūk), I 0 une partie mesurable de I, et u0 ∈ U tel que u0 (t) = u(t) si t ∈ I 0 , et
u0 (t) = ū(t) sinon. Comme ku0 − ūk1 ≤ 2M mes(I 0 ), le lemme 3.8 implique que
(54) F (u0 ) ≤ F (ū) − ε mes(I 0 ) + 4cM M 2 mes(I 0 )2 ,
ce qui donne une contradiction en choisissant la mesure de I 0 strictement positive et assez petite.
(b) Soient uk une suite dense dans Uad , et uk la suite de U définie par u0 := ū, et
(55) uk (t) = uk si H(uk , ȳ(t), p̄(t)) < H(uk−1 (t), ȳ(t), p̄(t)), et uk (t) = uk−1 (t) sinon.
La densité de uk implique que H(uk (t), ȳ(t), p̄(t)) → inf u∈Uad H(u, ȳ(t), p̄(t)) p.p. Prenant u = uk dans
(53), on conclut en appliquant le théorème de convergence monotone.
On vérifie facilement que si Uad = Rm , tout extrémale de Pontriaguine (en particulier quand ū est
solution de (39)) satisfait la condition de Legendre-Clebsch suivante :
(56) H̄uu (t) 0 p.p. sur [0, T ].
Introduisons une nouvelle notion de minimum.
8
Définition 3.9. On dit qu’une suite uk de U converge vers ū au sens de Pontriaguine s’il existe M > 0 tel
que kuk k ≤ M , et que uk → ū p.p. (par convergence dominée, uk → ū aussi dans Ls (0, T )m pour tout
s ∈ [1, ∞[). On dit que ū est un minimum de Pontriaguine si uk → ū au sens de Pontriaguine implique
F (ū) ≤ F (uk ) pour k assez grand.
Remarque 3.10. On vérifiera aisément que la démonstration du théorème 3.7. est encore valide quand ū
est un minimum de Pontriaguine. À tout minimum de Pontriaguine correspond donc une extrémale de
Pontriaguine.
Le lagrangien est défini comme une combinaison linéaire du critère et des contraintes par des coefficients
appelés multiplicateurs de Lagrange. L’approche lagrangienne des conditions d’optimalité consiste à for-
mer et à exprimer la stationnarité du lagrangien par rapport aux variables à optimiser. Redéfinissons le
préhamiltonien, sous la forme dite non qualifiée, H : R+ × Rm × Rn × Rn par
(58) H(β, u, y, p) := β`(u, y) + p · f (u, y).
nΨ
Notons Ψ ∈ R le multiplicateur associé aux contraintes aux deux bouts. On définit le lagrangien aux
deux bouts
(59) LIF (β, y 0 , yT , Ψ) := βϕ(y 0 , yT ) + Ψ · Φ(y 0 , yT ).
Le lagrangien du problème (57) est par définition
Z T
L(β, u, y, p, Ψ) := βJ IF (u, y) + p(t) · (f (u(t), y(t)) − ẏ(t))dt + Ψ · Φ(y(0), y(T ))
(60) Z T 0
= (H(β, u(t), y(t), p(t)) − p(t) · ẏ(t)) dt + LIF (β, y(0), y(T ), Ψ).
0
Cherchons les multiplicateurs (β̄, p̄, Ψ̄) rendant le lagrangien stationnaire (c’est-à-dire de dérivée par-
tielle nulle) par rapport à l’état, sur la trajectoire nominale. Ceci revient après intégration par parties
(supposant p̄ assez régulier pour cela) à vérifier, pour tout z ∈ Y :
Z T
˙
∇y H(β̄, ū(t), ȳ(t), p̄(t)) + p̄(t) · z(t)dt
(61) 0
+(∇y0 LIF (β̄, ȳ(0), ȳ(T ), Ψ̄) + p̄(0)) · z(0)
+(∇yT LIF (β̄, ȳ(0), ȳ(T ), Ψ̄) − p̄(T )) · z(T ) = 0.
On vérifie que ceci équivaut à la nullité des coefficients ci-dessus de z(t), z(0) et z(T ). Ceci conduit à
définir l’équation de l’état adjoint comme
(i) −p̄(t) ˙ = ∇y H(β̄, ū(t), ȳ(t), p̄(t)), p.p. t ∈ [0, T ],
(62) (ii) −p̄(0) = ∇y0 LIF (β̄, ȳ(0), ȳ(T ), Ψ̄),
(iii) p̄(T ) = ∇yT LIF (β̄, ȳ(0), ȳ(T ), Ψ̄).
Définition 4.1. On appelle multiplicateur de Pontriaguine, associé à la trajectoire nominale (ū, ȳ), tout
triplet λ̄ := (β̄, Ψ̄, p̄) vérifiant (62), β̄ ≥ 0, Ψ̄ ∈ NKΦ (Φ(ȳ(0), ȳ(T ))), la relation de non nullité
(63) β̄ + |Ψ̄| > 0,
ainsi que l’inégalité hamiltonienne, qui généralise (47) :
(64) H(β̄, ū(t), ȳ(t), p̄(t)) = inf H(β̄, u, ȳ(t), p̄(t)) p.p. sur [0, T ].
u∈Uad
On dit que (ū, ȳ) est une extrémale de Pontriaguine s’il existe au moins un multiplicateur de Pontriaguine
associé.
Définition 4.2. Si Uad est convexe, on définit de manière similaire un multiplicateur de Lagrange, en
changeant l’inégalité hamiltonienne (64) par la condition du premier ordre correspondante (comparer à
(41)), soit
(65) Hu (β̄, ū(t), ȳ(t), p̄(t))(u − ū(t)) ≥ 0, pour tout u ∈ Uad , p.p. sur [0, T ].
9
L’ensemble ΛP (ū, ȳ) des multiplicateurs de Pontriaguine est donc un cône (il est stable par multipli-
cation par un nombre strictement positif). Noter que la relation (63) empêche un multiplicateur d’être
nul. Si ΛP (ū, ȳ) n’est pas vide, il contient donc un élément avec, soit β̄ = 0 (multiplicateur singulier),
soit β̄ = 1 (multiplicateur régulier).
On trouvera une démonstration courte de ce résultat dans [13] (avec état initial fixé, mais l’extension
est facile), basée sur le “principe variationnel d’Ekeland”.
Exemple 4.4. Le problème avec état initial fixé et sans contraintes sur l’état final correspond au cas où
ϕ n’est fonction que de l’état final, Φ(y 0 , yT ) = y 0 , et KΦ est réduit à un singleton {ȳ 0 }. Dans l’équation
adjointe, les conditions initiale et finale se réduisent resp. à p̄(T ) = β̄∇ϕ(ȳ(T )) et p̄(0) = −Ψ̄ ; cette
dernière relation exprime Ψ̄ en fonction de p̄. Si β̄ = 0, la condition finale sur l’adjoint est nulle et la
dynamique de l’état adjoint étant linéaire, p̄ est nul donc aussi Ψ̄, ce qui contredit (63). Il existe donc
un multiplicateur avec β̄ = 1 ; on retrouve ainsi la conclusion du théorème 3.7.
Exemple 4.5. Plus généralement, si l’état initial est fixé, on peut écrire ϕ comme fonction de yT seul, et
supposer les contraintes aux deux bouts de la forme (y 0 , ΦF (y(T ))) ∈ {ŷ 0 }Rn × KF , avec KF convexe
fermé. On vérifie facilement qu’on obtient un système équivalent au PMP en remplaçant les conditions
aux deux bouts de l’équation adjointe par
Définition 4.6. On appelle condition de qualification toute hypothèse sur les contraintes, assurant la
positivité stricte de β̄.
Donnons-en un exemple simple (rappelons que z[v] est la solution de (12) avec z 0 = 0).
Définition 4.7. On dira que la contrainte finale ΦF (y(T )) ∈ KF est fortement qualifiée sur la trajectoire
nominale si
En l’absence de contraintes sur la commande et avec un état initial fixé, cette hypothèse assure bien
que β̄ 6= 0, car dans le cas contraire on aurait Ψ̄F 6= 0 d’après (67), donc notant z = z[v], pour un choix
convenable de v :
Z T Z T
† 0 ˙
(69) 0 6= Ψ̄F ΦF (ȳ(T ))z[v](T ) = p̄(T ) · z(T ) = (p̄(t) · z(t) + p̄(t) · ż(t))dt = H̄u (t)v(t)dt.
0 0
Or le dernier terme est nul car H̄u (t) = 0 p.p. en raison du PMP. Par le même type de calcul on vérifie
que la forte qualification entraîne l’unicité du multiplicateur pour lequel β̄ = 1. Si ΦF est l’identité, la
qualification forte s’interprète comme la contrôlabilité de l’état final du système linéarisé.
Exemple 4.8. Par un raisonnement analogue on vérifie plus généralement que, si les seules contraintes
sont initiales-finales, et si
on a nécessairement β > 0 et on peut donc se ramener à la forme qualifiée (β̄ = 1), pour laquelle on a
l’unicité du multiplicateur.
Remarque 4.9. Dans le cas où les contraintes finales modélisent un ensemble d’égalités et d’inégalités,
comme l’analyse est locale, on peut éliminer de la formulation du problèmes les contraintes d’inégalités
non actives sur la trajectoire de référence, ce qui conduit à une condition de qualification plus faible.
Voir d’autres conditions de qualification dans [10, Ch. 2].
10
4.2. Conservation du préhamiltonien. Soient (ū, ȳ) une trajectoire et p̄ ∈ Y satisfaisant (62)(i). On
dit que (ū, ȳ, p̄) est une extrémale de Pontriaguine sur ]0, T [ si l’inégalité hamiltonienne (64) est vérifiée.
Posons h(t) := H(β̄, ū(t), ȳ(t), p̄(t)).
Lemme 4.10. Soit (ū, ȳ, p̄) une extrémale de Pontriaguine sur ]0, T [. Alors h(t) est (à un ensemble de
mesure nulle près) constante sur ]0, T [.
Démonstration. a) Supposons d’abord ū(t) dérivable en tout point, et l’absence de contraintes sur la
commande. Alors Hu (t) = 0 p.p. d’après le PMP, et h(t) a pour dérivée, utilisant les équation d’état et
d’état adjoint :
˙
ḣ(t) = Hy (t)ȳ(t) ˙
+ Hp (t)p̄(t) ˙
= −p̄(t) ˙
· ȳ(t) ˙
+ ȳ(t) ˙
· p̄(t) = 0,
d’où le résultat.
b) Pour le cas général, adapter les arguments de [9, Partie II, Section 1.4].
4.3. Paramètres de décision. Reprenons le problème de minimisation de (24) sous la contrainte (25),
et sa reformulation (26)-(27). À l’état augmenté ya := (y, π) associons l’état adjoint pa := (p, pπ ). On
écrit le système d’optimalité sous forme qualifiée (β̄ = 1). Posons H(u, y, π, p) := `(u, y, π) + p · f (u, y, π).
Comme π̄(t) = π̄ pour tout t, l’équation adjointe (62) implique
˙
−p̄(t) = ∇y H(ū(t), ȳ(t), π̄, p̄(t)), p.p. t ∈ [0, T ],
−p̄˙ π (t) = ∇π H(ū(t), ȳ(t), π̄, p̄(t)), p.p. t ∈ [0, T ],
(71) p̄π (0) = 0,
p̄(T ) = ∇yT ϕ(ȳ(T ), π̄),
p̄π (T ) = ∇π ϕ(ȳ(T ), π̄).
RT
Comme p̄π (0) = − 0 p̄˙ π (t)dt + p̄π (T ), on en déduit que
Z T
0 = p̄π (0) = ∇π H(ū(t), ȳ(t), π̄, p̄(t))dt + ∇π ϕ(ȳ(T ), π̄),
0
qui n’est rien d’autre que la condition de stationnarité du lagrangien (du problème (24)-(25)) par rapport
à π, et on a d’autre part l’inégalité hamiltonienne comme conséquence du PMP : p.p. en t on a
H(ū(t), ȳ(t), π̄, p̄(t)) = inf H(u, ȳ(t), π̄, p̄(t)).
u∈Uad
4.4. Horizon variable. On reprend le problème (57), mais avec ici un horizon libre. Nous avons vu
comment se ramener au format standard en passant en temps normalisé, voir (28)-(29), dans lequel
l’horizon apparaît comme un paramètres de décision. D’après ce qui précède on sait que la condition
d’optimalité est identique à celle du cas à horizon fixé (théorème 4.3), avec en plus la condition de
stationnarité du lagrangien (du problème en temps normalisé) par rapport à l’horizon, ce qui équivaut
RT
à 0 H(β̄, ū(t), ȳ(t), p̄(t))dt = 0. Or le préhamiltonien est constant d’après le lemme 4.10 ; on retiendra
donc que le préhamiltonien est nul si l’horizon est libre. Ceci ne vaut, bien entendu, que si l’horizon
n’intervient pas dans les dynamique, coûts et contraintes du problème.
Exemple 4.11. Supposons que l’horizon intervienne dans la contrainte initiale-finale (rendez-vous avec un
engin spatial) qui s’écrit alors Φ(y(0), y(T ), T ) = 0. Notant par h̄ la valeur constante du préhamiltonien,
on voit que la condition de stationnarité du lagrangien par rapport à T donne alors
h̄ + Ψ̄ · DT Φ(ȳ(0), ȳ(T ), T ) = 0.
4.5. Temps minimal, géodésiques. C’est le cas particulier du problème précédent dans lequel le coût
final est nul, et l’intégrande du coût intégral est `(u, y) = 1, avec état initial fixé et contrainte finale
ΦF (ȳ(T )) ∈ KF . Le préhamiltonien β + p · f (u, y) est nul le long d’une extrémale de Pontriaguine car
l’horizon est libre, et la dynamique adjointe est (voir (66)-(67)) :
˙ = f¯y (t)† p̄(t) p.p.;
−p̄(t)
p̄(T ) = DΦF (ȳ(T ))† Ψ̄F ,
avec Ψ̄F ∈ NKF (ȳ(T )), de sorte que p̄ s’annule en un instant ssi il est identiquement nul. Le résultat
essentiel est :
Lemme 4.12. Si (ū, ȳ) est solution du problème en temps minimal, alors Ψ̄F 6= 0, et si Φ0F (ȳ(T )) est
surjective, alors p̄ ne s’annule en aucun instant.
11
Démonstration. On a Ψ̄F 6= 0, sinon p̄(t) = 0 pour tout t d’après l’équation adjointe, et comme
le temps final est libre, le préhamiltonien est nul donc β̄ = 0, ce qui contredit (67). Si Φ0F (ȳ(T )) est
surjective, sa transposée est injective, donc p̄(T ) = Φ0F (ȳ(T ))† Ψ̄F n’est pas nul, d’où le résultat.
L’inégalité hamiltonienne (64) se ramène à
p̄(t) · f (ū(t), ȳ(t)) = inf p̄(t) · f (u, ȳ(t)),
u∈Uad
le membre de gauche ci-dessus étant négatif, et nul ssi β̄ = 0 (hamiltonien nul et β̄ ≥ 0).
Remarque 4.13. Une alternative à la formulation avec coût intégral ci-dessus est la formulation avec
coût final T et coût intégral nul. Nous laissons le lecteur vérifier qu’elle aboutit aux mêmes conditions
d’optimalité. Par préhamiltonien d’un problème de temps minimal, on entend souvent celui de cette
seconde formulation, qui est p · f (u, y).
Exemple 4.14. Cas du système dynamique linéaire : f (u, y) = Ay+Bu avec A et B matrices de dimension
†
convenable. Alors −p̄(t) ˙ = A† p̄(t) donc p̄(t) = e(T −t)A p̄(T ) est une fonction analytique. Notons p̄(t)(k)
la dérivée kième de p̄(t). La commande doit minimiser sur Uad la fonction p̄(t)† Bu. Si celle-ci est
constamment nulle, alors pour k = 0 à n − 1, on a 0 = B † p̄(t)(k) = (−1)k B † (A† )k p̄(t) et donc p̄(t) est
orthogonal à A(k) B. Si le système est commandable, ce qui signifie que la famille des A(k) B, k = 0 à
n − 1, est de rang n, p̄ est identiquement nul : on sait que ceci est impossible si Φ0F (ȳ(T )) est surjective.
Si m = 1 et Uad = [−1, 1], supposant le système commandable et Φ0F (ȳ(T )) surjective, utilisant le fait
qu’une fonction analytique non constante a un nombre fini de changements de signes sur un intervalle
borné, on en déduit que la commande est bang-bang (elle vaut ±1 p.p.), plus précisément égale à −1
quand B † p̄(t) > 0, et à 1 quand B † p̄(t) < 0, avec un nombre fini d’instants de commutations.
Exemple 4.15 (Optique géométrique). La dynamique nominale est ȳ(t) ˙ = F (ȳ(t))ū(t), avec F : Rn → R
1
de classe C à valeurs strictement positives (indice variable du milieu) et ū(t) ∈ B̄ (boule unité euclidienne
fermée) p.p. On suppose Φ0F (ȳ(T )) surjective, donc l’état adjoint ne s’annule pas. Le préhamiltonien
β̄ + F (ȳ(t))p̄(t) · ū(t) atteint son minimum sur B̄ en −p̄(t)/|p̄(t)| ; comme il est nul ceci impose β̄ > 0
(on peut donc supposer que β̄ = 1). Le couple état-état adjoint a donc pour dynamique, supprimant le
temps en argument :
(72) ȳ˙ = −F (ȳ)p̄/|p̄|; p̄˙ = |p̄|∇F (ȳ).
Si sur une partie ouverte F est constant, p̄ et ū le sont aussi et la trajectoire est localement une ligne
droite.
RT
Remarque 4.16. Le problème de transfert en énergie minimale consiste à minimiser 21 0 u(t)2 dt avec les
contraintes finales et dynamique précédentes. Le préhamiltonien est 21 βu2 + F (y)p · u. On peut vérifier
comme précédemment que β̄ > 0 et donc supposer β̄ = 1. La commande minimisant le hamiltonien est
ū(t) = −F (ȳ(t))p̄(t), et le couple état-état adjoint a pour dynamique
ȳ˙ = −F (ȳ)2 p̄; p̄˙ = F (ȳ)|p̄|2 ∇F (ȳ).
Comme cette dynamique est proportionnelle à celle de (72), les lieux des trajectoires associées de (ȳ, p̄)
sont identiques, donc les extrémales de Pontriaguine de l’optique géométrique coincident avec celles du
transfert en énergie minimale. Par ailleurs, dans le problème de transfert en énergie minimale, comme
le préhamiltonien est constant, on vérifie que |ū(t)| = F (ȳ(t))|p̄(t)| l’est aussi : le transfert à énergie
minimale se fait à vitesse constante.
Exemple 4.17 (Géométrie riemannienne). Soit le problème de transfert en temps minimal avec la dyna-
mique nominale ȳ˙ = F (ȳ)ū, où la matrice symétrique F (y), de taille n, est inversible pour tout y, et
fonction de classe C ∞ de y, et la contrainte ū(t) ∈ B̄ p.p., ce qui revient à imposer ȳ˙ † F (ȳ)−2 ȳ˙ ≤ 1. On
peut encore se ramener à β̄ = 1 et montrer que p̄ ne peut s’annuler. Le préhamiltonien 1+p† F (y)u atteint
sur une extrémale de Pontriaguine son minimum en −F (ȳ)p̄/|F (ȳ)p̄|, donc ȳ˙ = −F 2 (ȳ)p̄/|F (ȳ)p̄|. De
plus, pour k = 1 à n, notant Fk (y) := ∂F (y)/∂yk , on a p̄˙ k = p̄† Fk (ȳ)F (ȳ)p̄/|F (ȳ)p̄|.
RT
Le problème correspondant en énergie minimale consiste à minimiser 21 0 u(t)2 dt sous les dynamique
et contrainte précédentes. Le préhamiltonien 21 u2 + p† F (y)u atteint son minimum en u = −F (y)p, et
donc ȳ˙ = −F 2 (ȳ)p̄ et p̄˙ k = p̄† Fk (ȳ)F (ȳ)p̄. Ici encore, cette dynamique est proportionnelle à celle du
temps minimal : les lieux de (ȳ, p̄) pour ces deux problèmes coïncident donc. On appelle géodésiques les
courbes donnant le lieu de l’état pour toute paire d’état initial et final possible. Le préhamiltonien du
problème en énergie minimale étant constant, la vitesse |ū(t)| = |F (ȳ(t))p̄(t)| l’est aussi.
Nous reverrons des exemples de problème en temps minimal en section 8.1.
12
5. Algorithme de tir et discrétisation du problème
5.1. Problèmes sans contraintes. On suppose ici que u 7→ F (u) a un minimum faible ū, fonction
continue du temps, satisfaisant la condition de Legendre-Clebsch forte (uniforme en temps) suivante, où
Id est la matrice identité (comparer à la condition nécessaire (56)) : Il existe α > 0 tel que
(73) H̄uu (t) αId , p.p. sur [0, T ].
Appliquant à la relation Hu (u, y, p) = 0 le théorème des fonctions implicites, on déduit que, pour τ ∈ [0, T ]
et (u, y, p) voisin de (ū(τ ), ȳ(τ ), p̄(τ )), il existe Υ : Rn × Rn → Rm de classe C ∞ tel que la relation
Hu (u, y, p) = 0 équivaut à u = Υ(y, p). Ajoutant si nécessaire le temps comme état supplémentaire, on
peut “recoller” les fonctions Υ pour différents instants et on obtient ainsi que pour un certain ε > 0, si
|y − ȳ(t)| + |p − p̄(t)| < ε, alors
Hu (u, y, p) = 0 ssi u = Υ(y, p).
On notera Υ(t) := Υ(ȳ(t), p̄(t)).
Remarque 5.1. Quand dans la suite on fera référence à (73), il faudra comprendre H̄uu (t) comme la
dérivée seconde par rapport à la commande du préhamiltonien spécifique de chaque problème.
Remarque 5.2. Différenciant en temps la relation ∇u H(ū(t), ȳ(t), p̄(t)) = 0, il vient
H̄uu (t)Υ̇(t) + H̄uy (t)f¯(t) − f¯u (t)† ∇y H̄(t) = 0,
˙
ce qui permet d’exprimer ū(t) = Υ̇(t) en fonction de (ū(t), ȳ(t), p̄(t)).
Éliminant la commande dans l’équation d’état et d’état adjoint on voit que (ȳ(t), p̄(t)) est solution
de l’équation différentielle autonome
˙
ȳ(t) = f (Υ(t), ȳ(t)),
(74) ˙ p.p. t ∈ [0, T ],
−p̄(t) = ∇y H(Υ(t), ȳ(t), p̄(t)),
avec les conditions aux deux bouts
(75) ȳ(0) = ȳ 0 ; p̄(T ) = ∇ϕ(ȳ(T )).
Dans le cas particulier d’un problème quadratique, on retrouve les relations (38). Introduisons le para-
mètre p0 ∈ Rn , et notons (si elle existe) par (y[p0 ], p[p0 ]) la solution de (74) avec la condition initiale
(ȳ 0 , p0 ). Alors le système au deux bouts (74)-(75) se ramène à trouver un zéro de la fonction de tir
Rn → Rn définie par
(76) T (p0 ) := ∇ϕ(y[p0 ](T )) − p[p0 ](T ).
Une telle équation est habituellement résolue par une variante de la méthode de Newton. Cette dernière
nécessite l’évaluation du jacobien de T . Donnons son expression quand p0 = p̄(0) et donc (y[p0 ], p[p0 ]) =
(ȳ, p̄). On vérifie que
(77) DT (p0 )q 0 = D2 ϕ(ȳ(T ))z(T ) − q(T ),
où (z(T ), q(T )) est solution de la linéarisation de (74) avec la condition initiale (z(0), q(0)) = (0, q 0 ).
Explicitons les relations satisfaites par cette linéarisation. On a, pour un certain v ∈ U :
ż(t) = f¯0 (t)(v(t), z(t)),
(78) −q̇(t) = f¯y (t)† q(t) + H̄yu (t)v(t) + H̄yy (t)z(t) p.p. t ∈ [0, T ],
0 = H̄uu (t)v(t) + H̄uy (t)z(t) + f¯u (t)† q(t).
5.2. Problèmes avec contraintes aux deux bouts. On se restreint au cas de la minimisation de J IF
avec contraintes d’égalités aux deux bouts, et donc KΦ = {0}RnΦ , sous la condition de qualification (70)
qui permet de fixer β̄ à 1. La fonction de tir a pour paramètres (y 0 , p0 , Ψ), les valeurs initiales de l’état
et de l’état adjoint, ainsi que le multiplicateur associé aux contraintes. Elle exprime les conditions aux
deux bouts pour l’état adjoint ainsi que les contraintes :
p0 + ∇y0 LIF (1, y 0 , y(T ), Ψ)
5.3. Discrétisation du problème. Discrétisons l’équation d’état (8) par la méthode d’Euler avec un
pas de temps constant, supposant la commande constante sur chaque pas de temps, soit
(80) yk+1 = yk + hf (uk , yk ), k = 0, . . . , N − 1; y0 = y 0 ,
avec h = T /N , N entier strictement positif, uk ∈ Rm et yk ∈ Rn . L’espace des commandes discrètes est
U N := RN m , celui des états discrets est Y N := R(N +1)n . À chaque u ∈ U N est associé un unique état
y[u, N ] ∈ Y N solution de (80). La discrétisation correspondante du critère (13) est
N
X −1
N
(81) J (u, y) := h `(uk , yk ) + ϕ(yN ).
k=0
on obtient l’équation adjointe discrétisée en rendant le lagrangien stationnaire par rapport à l’état, soit
pk = pk+1 + h ∇y `(uk , yk ) + fy (uk , yk )† pk+1 , k = 0, . . . , N − 1; pN = ∇ϕ(yN ).
(84)
Le critère réduit est F N (u) := J N (u, y[u, N ]). Opérant de manière analogue au cas du temps continu,
et notant y = y[u, N ], on peut vérifier qu’il satisfait
N
X −1
(85) DF N (u)v = Hu (uk , yk , pk+1 )vk .
k=0
14
6. Conditions d’optimalité du second ordre
6.1. Développements du critère réduit. Analysons un problème sans contrainte autre que l’état
initial fixé. Soit ū ∈ U. Nous avons noté en section 2.1 que l’application ū 7→ F (ū) est de classe C ∞ sur
U. Notons par Q0 (v) := F 00 (ū)(v, v) la form quadratique associée au hessien (dérivée seconde) de F . Le
développement de F au second ordre, en v ∈ U, s’écrit :
F (ū + v) = F (ū) + F 0 (ū)v + 21 Q0 (v) + R2 (v),
où le terme de reste
Z 1
R2 (v) = (1 − σ)(F 00 (ū + σv) − F 00 (ū))(v, v)dσ
0
vérifie |R2 (v)| = O(kvk3∞ ). Une condition nécessaire d’optimalité locale dans U est donc
F 0 (ū) = 0 et Q0 (v) ≥ 0, pour tout v ∈ U.
Or, du lemme 6.4 ci-dessous on déduit aisément que |R2 (v)| = O(kvk∞ kvk22 ), d’où la condition suffisante
d’optimalité locale suivante :
Théorème 6.1. Si F 0 (ū) = 0, et si
(92) Il existe α > 0 tel que Q0 (v) ≥ 2αkvk22 pour tout v ∈ U,
alors ū est un minimum faible de F , et vérifie la condition de croissance quadratique suivante : si kvk∞
est assez petit, alors F (ū + v) ≥ F (ū) + 21 αkvk22 .
On dira que le préhamiltonien satisfait une hypothèse de croissance quadratique si, pour tout M > 0,
il existe αM > 0 tel que si |u| < M , alors
H(u, ȳ(t), p̄(t)) ≥ H̄(t) + 2αM |u − ū(t)|2 .
On peut déduire du lemme précédent le résultat suivant :
Théorème 6.3. Supposons (92) satisfaite ainsi que l’hypothèse de croissance quadratique du préhamil-
tonien. Soit uk → ū au sens de Pontriaguine. Alors il existe α > 0, fonction de M := supk kuk k∞ , telle
que pour k assez grand :
F (uk ) ≥ F (ū) + αkuk − ūk22 .
6.2. Expression du hessien du critère réduit. Considérons le problème de minimisation, sans
contrainte autre que l’état inital fixé et l’équation d’état (8), du critère (13). Le lagrangien associé
est, cf. (60) :
Z T
J(u, y) + p(t) · (f (u(t), y(t)) − ẏ(t))dt
Z T 0 Z T
= H(u(t), y(t), p(t))dt + ϕ(y(T )) − p(t) · ẏ(t))dt.
0 0
Notons par D2 H̄(t) la dérivée seconde du hamiltonien par rapport à (u, y), au point (ū(t), ȳ(t), p̄(t)).
La dérivée seconde du lagrangien dans la direction (v, z) ∈ U × Y est
Z T
Ω0 (v, z) := D2 H̄(t)(v(t), z(t))2 dt + ϕ00 (ȳ(T ))(z(T ))2 .
0
où J IF est définie en (21), avec KΦ = K[n0Φ , n00Φ ] (nombre fini d’égalités et d’inégalités). Le lagrangien a
pour expression (60). Sa dérivée seconde (hessien) en (β, ū, ȳ, p, Ψ), dans la direction (v, z) ∈ U × Y est,
le préhamiltonien étant sous forme non qualifiée (58), et posant λ := (β, p, Ψ) :
Z T
Ω0 [λ](v, z) := D2 H̄(t)(v(t), z(t))2 dt + D2 LIF (β, ȳ(0), ȳ(T ), Ψ)(z(0), z(T ))2 ,
0
2 IF
où par D L on entend les dérivées secondes par rapport à (y 0 , yT ). Posons
Q[λ](v, z 0 ) := Ω[λ](v, z[v, z 0 ]).
Comme cette forme quadratique est continue dans la norme de U2 × Rn , on peut l’étendre à cet espace.
Les conditions du second ordre font intervenir le cône critique, dont les éléments sont appelés directions
critiques, défini comme suit, où z = z[v, z 0 ] (on a défini le cône tangent TKΦ (·) en section 1) :
(v, z 0 ) ∈ U2 × Rn ; ϕ0 (ȳ(0), ȳ(T ))(z(0), z(T )) ≤ 0
C2 (ū, ȳ) := .
Φ0 (ȳ(0), ȳ(T ))(z(0), z(T )) ∈ TKΦ (Φ(ȳ(0), ȳ(T )))
Théorème 6.5. Si (ū, ȳ) est solution faible (resp. de Pontriaguine) de (93), à toute direction critique
(v, z 0 ) est associée un multiplicateur de Lagrange (resp. de Pontriaguine) tel que Q[λ](v, z 0 ) ≥ 0.
Donnons maintenant des conditions suffisantes. Rappelons la notation ΛP (ū, ȳ) de l’ensemble des
multiplicateurs de Pontriaguine. Nous utiliserons les hypothèses suivantes :
Pour tout M > kūk, il existe λ̂M ∈ ΛP (ū, ȳ), λ̂M = (β̂, p̂, Ψ̂), et αM > 0 tels que
(94)
H(β̂, ū(t), ȳ(t), p̂(t)) + αM |u − ū(t)|2 ≤ H(β̂, u, ȳ(t), p̂(t)), si |u| ≤ M , p.p. t.
On choisit C([0, T ])ng comme espace de Banach pour la contrainte sur l’état. On sait que toute forme
RT
linéaire continue sur C([0, T ]) est du type h 7→ 0 h(t)dµ(t), où dµ ∈ M (0, T ), l’espace des mesures de
Borel sur [0, T ]. On peut identifier l’action d’une telle mesure dµ à l’intégrale de Stieltjes d’une fonction
à variation bornée nulle au temps T , notée µ, qu’on suppose nulle en T . Rappelons que la définition de
la variation d’une fonction µ sur [0, T ] est
( npi )
X
(99) var(µ) := sup |µ(ti+1 ) − µ(ti )|; où (ti ) est une subdivision de [0, T ] ,
i=0
(sur la notion de subdivision, voir (19)), et que l’intégrale de Stieltjes de h ∈ C([0, T ]) associé à la
fonction à variation bornée µ est
Z T npi
X
(100) h(t)dµ(t) := lim h(τi )(µ(ti+1 ) − µ(ti ))
0 i=0
où τi ∈ [ti , ti+1 ], la limite étant prise sur les subdivisions d’incrément maximum tendant vers 0.
Le lagrangien du problème (98) est identique à celui du problème (57) (défini en (60)), à l’ajout près
de la contribution ∆g : Y × M (0, T )ng → R de la contrainte sur l’état, soit
ng Z T
X
∆g (y, µ) := gj (y(t))dµj (t).
j=1 0
Notons la relation, où z ∈ Y :
ng Z T
X
∆gy (y, µ)z = gj0 (y(t))z(t)dµj (t).
j=1 0
On obtient encore l’équation adjointe en rendant le lagrangien stationnaire par rapport à l’état. On
cherche maintenant l’état adjoint dans l’espace P := BV (0, T )n des fonctions à variation bornée, ce qui
permet d’utiliser la formule d’intégration par parties
Z T n Z T
X
(101) p(t) · ẏ(t)dt = p(T ) · y(T ) − p(0) · y(0) − yi (t)dpi (t).
0 i=1 0
L’équation de l’état adjoint a alors pour expression, LIF (·) étant défini en (59) :
ng
X
−dp̄(t) = ∇y H(β̄, ū(t), ȳ(t), p̄(t))dt + ∇gj (ȳ(t))dµ̄j (t), p.p. t ∈ [0, T ],
(102) j=1
−p̄(0) = ∇y0 LIF (β̄, ȳ(0), ȳ(T ), Ψ̄),
p̄(T ) = ∇yT LIF (β̄, ȳ(0), ȳ(T ), Ψ̄).
On appelle multiplicateur de Pontriaguine, associé à la trajectoire nominale, tout λ = (β̄, Ψ̄, µ̄, p̄) vérifiant
(102), Ψ̄ ∈ NKΦ (Φ(ȳ(0), ȳ(T ))), β̄ ≥ 0, l’inégalité hamiltonienne (64), et les relations de complémentarité
Z T
(103) dµ̄ ≥ 0; gj0 (ȳ(t))dµ̄j (t) = 0, j = 1, . . . , ng ,
0
et de non nullité
(104) β̄ + |Ψ̄| + kµ̄k > 0.
Théorème 7.1. Si (ū, ȳ) est solution de (39), c’est une extrémale de Pontriaguine.
On peut prouver ce résultat en étendant la technique de [13].
17
7.2. Ordre des contraintes sur l’état et conditions de jonction. Pour simplifier l’exposé on supp-
pose l’absence de contraintes sur la commande et l’extrémale qualifiée (β̄ = 1). Notons les sauts des
fonctions du temps ayant des limites à droite et à gauche (c’est le cas des fonctions à variation bornée)
à l’instant τ ∈]0, T [, par un crochet, ainsi [p̄(τ )] := p̄(τ+ ) − p̄(τ− ), et posons ν̄j (τ ) := [µ̄j (τ )]. En raison
de (102), on a pour τ ∈ [0, T ] :
ng
X
(105) − [p̄(τ )] = ν̄j (τ )∇gj (ȳ(τ )).
j=1
On appelle dérivée totale de g (suivant la dynamique f ) la fonction g (1) (u, y) := g 0 (y)f (u, y). Le long
(1)
d’une trajectoire (u, y), g (1) (u, y) est la dérivée de g(y). Si g (1) (u, y) dépend de u, au sens où gu (u, y)
n’est pas identiquement nulle, on dit que g est une contrainte du premier ordre sur l’état. Dans le cas
contraire on peut écrire g (1) comme fonction de y seul, et former la dérivée totale à l’ordre deux, omettant
les arguments de f = f (u, y) :
g (2) (u, y) := Dg (1) (y)f (u, y) = g 00 (y)(f, f ) + g 0 (y)fy f.
Si g (2) (u, y) dépend de u, on dit que g est d’ordre 2 et ainsi de suite pour un ordre quelconque. Si g est
vectorielle, on note qi l’ordre supposé bien défini de sa ième composante, Ij l’ensemble des contraintes
d’ordre j, et Ij (t) celles actives à l’instant t. Considérons l’hypothèse
(1)
ū est continu et {∇u gj (ū(t), ȳ(t)), j ∈ I1 (t)}
(106)
est une famille libre, pour tout t ∈ [0, T ].
Lemme 7.2. Soient (ū, ȳ) une extrémale de Pontriaguine et (β̄, p̄, µ̄) un multiplicateur de Pontriaguine
associé. Alors (i) Si H(·, ȳ(t), p̄(t)) est (uniformément en temps) fortement convexe, ū est continue. (ii)
Si (106) est satisfait, les composantes de µ̄ associées aux contraintes du premier ordre sont continues
sur ]0, T [.
Le point (ii) est conséquence de la relation suivante, qui utilise (105) :
ng
(1)
X
0 = [H̄u (t)] = [p̄(t)] · f¯u (t) = ν̄j (t)∇u gj (ū(t), ȳ(t)).
j=1
Montrons maintenant le lien entre le saut de la dérivée de la commande et les sauts de µ ou de sa dérivée.
Lemme 7.3. Si (106) est satisfait, alors (i) H̄uu (t) est continue, et (ii) si les sauts ci-dessous sont bien
définis, on a pour tout τ ∈]0, T [ :
˙ j (τ )]∇u gj(1) (τ ) − (2)
X X
(107) ˙ )] =
H̄uu (τ )[ū(τ [µ̄ ν̄j (τ )∇u gj (τ ).
j∈I1 (τ ) j∈I2 (τ )
(1)
Démonstration. (i) Si j ∈ I1 (τ ), alors ν̄j (τ ) = 0 d’après le lemme 7.2 (ii), et sinon, gj (u, y) ne
dépend pas de u, de sorte que, pour 1 ≤ i, k ≤ m :
ng ng
(1)
X X
−[H̄ui uk (τ )] = −[p̄(τ )] · f¯ui uk (τ ) = ν̄j (τ )ḡj0 (τ )f¯ui uk (τ ) = ν̄j (τ )Du2 i uk gj (ū(τ ), ȳ(τ )) = 0.
j=1 j=1
(ii) Soit i ∈ {1, . . . , m}. On a, supposant pour simplifier ` = 0 et omettant le temps en argument :
d d
p̄ · f¯ui = H̄ui u ū˙ + p̄† f¯ui y f¯ + p̄˙ · f¯ui .
Hui =
dt dt
Utilisant l’équation adjointe et le fait que ḡj0 f¯u = 0 si j 6∈ I1 (τ ), il vient :
d X
(108) Hui = H̄ui u ū˙ + p̄ · (f¯ui y f¯ − f¯y f¯ui ) − ˙ j ḡj0 f¯ui .
µ̄
dt
j∈I1
Comme H̄ui = H̄˙ ui = 0 p.p., on a [H̄˙ ui ] = 0. La commande et H̄uu étant continus, ceci implique
X
(109) ˙ = −[p̄] · (f¯ui y f¯ − f¯y f¯ui ) +
H̄ui u [ū] ˙ j ]ḡj0 f¯ui .
[µ̄
j∈I1
18
0 ¯
Or, si j 6∈ I1 (τ ), on a 0 = ∂
∂y (ḡj fui ) f¯ = ḡj00 (f¯, f¯ui ) + ḡj0 f¯ui y f¯, et donc
(2)
gj,ui = Dui (ḡj00 (f¯, f¯) + ḡj0 f¯y f¯)
(110) = 2ḡj00 (f¯, f¯ui ) + ḡj0 f¯ui y f¯ + ḡj0 f¯y f¯ui
= ḡj0 (f¯y f¯ui − f¯ui y f¯).
ng
X
On conclut en remplaçant [p̄] = − ν̄j ∇ḡj dans (109), et en notant que ν̄j = 0 quand qj = 1, et
j=1
(2)
ḡj,ui = 0 si qj > 2.
Exemple 7.4. Supposons (106) satisfait, H̄uu (τ ) inversible, et la contrainte sur l’état scalaire. Si elle est
(1)
du premier ordre et gu (τ ) 6= 0, on obtient H̄uu (τ )[ū(τ ˙ )]∇u g (1) (τ ) et donc ū(τ
˙ )] = [µ̄(τ ˙ ) est discontinue
(2)
ssi µ̄ l’est. Si la contrainte est du second ordre et gu (τ ) 6= 0, on obtient H̄uu (τ )[ū(τ )] = −ν̄(τ )∇u g (2) (τ )
˙ ˙
˙ ) est discontinue ssi µ̄ l’est. Si la contrainte n’est pas d’ordre un ou deux, alors ū(τ
et ū(τ ˙ ) est continue.
On voit que le comportement aux jonctions dépend d’une façon essentielle de l’ordre de la contrainte.
Exemple 7.5. Soit le problème de calcul de l’équilibre d’une corde élastique pesante de longueur T ,
soumise à un obstacle et fixée à ses extrémités. Plus spécifiquement, on minimise la somme des énergies
RT
élastique et potentielle : E(y) := 0 ( 12 ȳ(t)
˙ 2 + ȳ(t))dt, où ȳ(t) représente la déformation verticale, sous
les contraintes ȳ(0) = ȳ(T ) = 1 et ȳ(t) ≥ 0 pour tout t. On se ramène au cadre de la commande optimale
avec l’équation d’état ȳ(t)˙ = ū(t), un critère intégral d’intégrande `(u, y) = 21 u2 + y, et la contrainte sur
l’état g(y) = −y. Quand T est assez grand, la solution optimale est, pour un certain t0 ∈]0, 12 T [, nulle sur
[t0 , t1 ] avec t1 = T −t0 , et strictement positive en dehors. Le préhamiltonien est H(u, y, p) = 12 u2 +y +pu,
et la dynamique adjointe est −dp̄(t) = dt − dµ̄(t), donc (comme µ̄(T ) = 0), p̄(t) = T − t + µ̄(t) + p̄(T ).
L’inégalité hamiltonienne donne ū(t) = −p̄(t). Sur ]t0 , t1 [, on a donc 0 = ū(t) = −p̄(t), ce qui implique
˙ = 1. Sur chaque arc (intervalle de temps) non contraint, on a p̄(t) = c − t pour un certain c ∈ R et
µ̄
donc ū(t) = t − c. D’après le lemme 7.2, la commande ainsi que µ̄ sont continus en t0 et t1 , donc la
commande vaut t − t0 sur [0, t0 [ et t − t1 sur [t1 , T [. Le saut de ū˙ vaut -1 en t0 et 1 en t1 , et (comme
˙ = 0 hors de [t0 , t1 ]) le saut de µ̄
µ̄ ˙ vaut 1 en t0 et -1 en t1 , de sorte que (107) est bien satisfait.
Exemple 7.6. Soit le problème de calcul de l’équilibre d’une poutre élastique pesante de longueur T ,
soumise à un obstable et fixée à ses extrémités. On minimise la somme des énergies élastique et po-
RT
tentielle : E(x) := 0 ( 21 x̄ ¨ (t)2 + x̄(t))dt, où x̄(t) représente la déformation verticale, sous les contraintes
x̄(0) = x̄(T ) = 1 et x̄(t) ≥ 0 pour tout t. On se ramène au cadre de la commande optimale avec l’équation
d’état ȳ˙ 1 (t) = ȳ2 (t), ȳ˙ 2 (t) = ū(t), un critère intégral d’intégrande `(u, y) = 21 u2 + y1 , et la contrainte
sur l’état g(y) = −y1 , et les contraintes aux extrémités, par exemple ȳ1 (0) = ȳ1 (T ) = 1. Quand T est
assez grand, la solution optimale est, pour un certain t0 ∈]0, 21 T [, nulle sur [t0 , t1 ] avec t1 = T − t0 , et
strictement positive en dehors. Le préhamiltonien est H(u, y, p) = 21 u2 + y1 + p1 y2 + p2 u, et la dynamique
adjointe est −dp̄1 (t) = dt−dµ̄(t), −p̄˙ 2 (t) = p̄1 (t), donc (comme µ̄(T ) = 0), p̄1 (t) = T −t+ µ̄(t)+ p̄1 (T ).
L’inégalité hamiltonienne donne ū(t) = −p̄2 (t). Sur ]t0 , t1 [, on a donc
0 = ū(t) = −p̄2 (t) = −p̄˙ 2 (t) = p̄1 (t),
ce qui implique µ̄ ˙ = 1. Sur chaque arc non contraint, on a p̄1 (t) = c1 − t pour un certain c1 , donc
p̄2 (t) = 2 t − c1 t − c2 pour un certain c2 , et finalement ū(t) = c1 t + c2 − 21 t2 . D’après le lemme 7.2,
1 2
On suppose toutes les fonctions de classe C ∞ . Le préhamiltonien est p · (f0 (y) + uf1 (y)). D’après le
théorème 4.3 et l’exemple 4.5, une solution (ū, ȳ) est une extrémale de Pontriaguine : il existe un multi-
plicateur (β̄, Ψ̄F ) ∈ R+ × NKF (ΦF (ȳ(T ))) non nul tel que l’état adjoint solution de
(
˙
−p̄(t) = (f00 (ȳ(t)) + ū(t)f10 (ȳ(t))† p̄(t) p.p.,
p̄(T ) = β̄∇ϕ(ȳ(T )) + DΦF (ȳ(T ))† Ψ̄F
vérifie l’inégalité hamiltonienne (64), qui se ramène ici à
(115) (ū(t) − u)p̄(t) · f1 (ȳ(t)) ≤ 0, pout tout u ∈ [−1, 1].
Transfert en temps minimal. La théorie qui suit s’applique aussi si on remplace le coût final par la
minimisation de l’horizon ; dans la dynamique adjointe, il faut alors remplacer la condition finale par
p̄(T ) = DΦ(ȳ(T ))† Ψ̄F .
En raison du lemme 4.12, Ψ̄F 6= 0.
Non nullité de l’état adjoint. La dynamique de l’état adjoint fait que ce dernier est nul à un instant donné
ssi il est identiquement nul, et alors l’inégalité hamiltonienne ne donne pas d’information. On peut exclure
ce cas : (i) dans le problème avec coût final, comme (β̄, Ψ̄) 6= 0, si la famille {∇ϕ(ȳ(T )), ∇Φi (ȳ(T )), i =
20
1, . . . , nΦ } est linéairement indépendante, et (ii) pour le problème de transfert en temps minimal si
DΦ(ȳ(T )) est surjective (lemme 4.12). Nous ferons donc dans certains cas l’hypothèse d’un état adjoint
non nul.
8.2. Arcs bang et singuliers. Notant Ξ(t) := p̄(t) · f1 (ȳ(t)) la fonction de commutation, on peut
exprimer l’inégalité hamiltonienne (115), pour t ∈ [0, T ], comme
ū(t) = 1 si Ξ(t) < 0,
(116) ū(t) = −1 si Ξ(t) > 0,
ū(t) ∈ [−1, 1] si Ξ(t) = 0.
Soient 0 ≤ a < b ≤ T . On dira que (a, b) est un arc singulier (arc S) si ū(t) ∈] − 1, 1[ p.p. sur ]a, b[, et
un arc bang si ū est constant, égal à ±1 sur ]a, b[ (on parlera d’arc B− s’il vaut -1, d’arc B+ sinon).
Soient X et Y deux champs de vecteurs. Leur crochet de Lie est le champ de vecteurs défini par
[X, Y ](y) := X 0 (y)Y (y) − Y 0 (y)X(y),
0
avec Xij (y) = ∂Xi (y)/∂yj , soit
n
X
0
[X, Y ]i (y) := Xij (y)Yj (y) − Yij0 (y)Xj (y).
j=1
L’application (X, Y ) 7→ [X, Y ] est bilinéaire et antisymétrique et, si X est un champ de vecteurs, alors
d
(p̄(t) · X(ȳ(t))) = p̄(t) · [X, f0 ](ȳ(t)) + ū(t)p̄(t) · [X, f1 ](ȳ(t)).
dt
Comme [f1 , f1 ] = 0, prenant X = f1 , il vient successivement
Ξ̇(t) = p̄(t) · [f1 , f0 ](ȳ(t)),
Ξ̈(t) = p̄(t) · [[f1 , f0 ], f0 ](ȳ(t)) + ū(t)p̄(t) · [[f1 , f0 ], f1 ](ȳ(t)).
Sur un arc singulier, Ξ̈(t) = 0, donc si p̄(t) · [[f1 , f0 ], f1 ](ȳ(t)) 6= 0, on peut exprimer la commande en
fonction de l’état et de l’état adjoint, soit, utilisant l’antisymétrie du crochet : ū(t) = Γ(ȳ(t), p̄(t)), où
(117) Γ(y, p) := p · [[f1 , f0 ], f0 ](y)/p · [[f0 , f1 ], f1 ](y).
Il faut quelquefois dériver la fonction de commutation plus de deux fois pour obtenir une expression de
la commande en fonction de l’état et de l’état adjoint, et cette expression est toujours obtenue après un
nombre pair de dérivations, comme l’illustre l’exemple suivant.
RT
Exemple 8.1. Soit le problème de Fuller [15] : minimisation de 21 0 y1 (t)2 dt, avec équation d’état ẏ1 (t) =
y2 (t), ẏ2 (t) = u(t), état initial fixé, et u(t) ∈ [−1, 1]. Pour se ramener à un coût final, on introduit l’état
y3 de dynamique ẏ3 (t) = 21 y1 (t)2 , avec y3 (0) = 0, et le coût devient y3 (T ). Le préhamiltonien est
p1 y2 + p2 u + 21 p3 y12 , et l’équation adjointe est −p̄˙ 1 (t) = p̄3 (t)ȳ1 (t), −p̄˙ 2 (t) = p̄1 (t), −p̄˙ 3 (t) = 0, avec
la condition finale p̄(T ) = (0 0 1)† , donc p̄3 (t) = 1 pour tout t. Il faut dériver quatre fois la fonction
de commutation Ξ(t) = p̄2 (t) pour obtenir Ξ(4) (t) = ū(t). L’état et la commande sont nuls sur un
arc singulier. On généralise facilement cet exemple à celui d’un intégrateur n fois, pour lequel il faut
dériver 2n fois la fonction de commutation pour retrouver une expression de la commande. Notons que
la commande optimale comporte en général un nombre infini de commutations entre ±1.
8.3. Cas particuliers - Problèmes de dimension 2 et 3. On dira que t ∈ [0, T ] est critique si
Ξ(t) = 0, et bicritique si de plus Ξ̇(t) = 0. Notons T l’ensemble des instants critiques. Comme Ξ̇(·) est
continue, un instant critique mais non bicritique est un point isolé du fermé T . En particulier, s’il n’existe
aucun instant bicritique, alors [0, T ] \ T est une union finie d’arcs bang. Si t est bicritique, alors p̄(t) est
orthogonal à f1 (t) et [f1 , f0 ](t).
Cas de deux états (n = 2). On suppose satisfaite l’hypothèse d’état adjoint non nul. Si t ∈ [0, T ] est tel
que f1 (t) et [f1 , f0 ](t) sont linéairement indépendants, t ne peut être bicritique, et donc :
— Si f1 (t) et [f1 , f0 ](t) sont linéairement indépendants pour tout t ∈ [0, T ], alors [0, T ] \ T est une
union finie d’arcs bang.
— Sinon, les états aux instants bicritiques (et en particulier les arcs singuliers) sont contenus dans
la variété singulière définie par (notant par ∧ le produit vectoriel)
{y ∈ R2 ; f1 (y) ∧ [f1 , f0 ](y) = 0}.
21
RT
Exemple 8.2. Soit le problème de minimiser 21 0 x̄(t)2 dt sous la contrainte x̄(t) ˙ = ū(t) ∈ [−1, 1] (qui se
résout facilement !). On se ramène à un coût final en posant ȳ(t) ∈ R2 , ȳ˙ 1 (t) = u(t), ȳ˙ 2 (t) = 12 ȳ1 (t)2 , de
dynamique f0 (ȳ) = (0 12 ȳ12 )† , f1 (ȳ) = (1 0)† . Il vient
[f1 , f0 ](y) = −f00 (y)f1 (y) = (0 − y1 )† ,
et finalement f1 (y) ∧ [f1 , f0 ](y) = −y1 . La variété singulière est donc l’axe vertical y1 = 0 (sur lequel la
commande est nulle).
Cas de trois états (n = 3). Posons G(y) := f1 (y) ∧ [f1 , f0 ](y). Si t est un instant bicritique tel que f1 (t)
et [f1 , f0 ](t) sont linéairement indépendants, alors l’état adjoint est colinéaire à G(ȳ(t)). En particulier,
quand p̄(t) · [[f0 , f1 ], f1 ](ȳ(t)) ne s’annule pas sur un arc singulier, comme seule la direction de l’adjoint
intervient dans l’expression (117) de la commande, on obtient l’expression de la commande fonction de
l’état : ū(t) = Γ(ȳ(t), G(ȳ(t))).
Si de plus l’horizon est libre, le préhamiltonien s’annule le long de la trajectoire optimale. Sur un arc
singulier, comme p̄(y) · f1 (ȳ(t)) = 0, on en déduit que p̄(y) · f0 (ȳ(t)) = 0. Si l’état adjoint ne s’annule
pas, les états aux instants bicritiques (et en particulier les arcs singuliers) sont contenus dans la variété
singulière {y ∈ Rn ; G(y) ∧ f0 (y) = 0}.
Exemple 8.3 (Problème de Goddard). Il s’agit d’une variante « monodimensionnelle » de l’exemple de
l’introduction. On désire maximiser l’altitude finale d’un vol vertical de fusée. Les variables d’état sont
l’altitude h, la vitesse V et la masse m, de dynamique ḣ = V , mV̇ = cu − D(h, V ), ṁ = −u pour un
certain c > 0, D(h, V ) ≥ 0 représentant la traînée aérodynamique, avec contrainte u(t) ∈ [0, 1] pour
tout t, et m(T ) ≥ mT . On maximise donc h(T ) avec horizon libre. On peut appliquer à cet exemple la
théorie précédente (nous ne détaillerons pas les calculs), ce qui permet d’exprimer la commande fonction
de l’état sur l’arc singulier, et d’expliciter l’équation de la variété singulière.
8.4. Contraintes sur l’état. En présence de contraintes sur l’état de type (97), le théorème 7.1 s’ap-
plique et l’expression (105) des sauts de l’état adjoint reste valable. Nous allons reprendre l’étude des
conditions de jonction, en tenant compte de la contrainte sur la commande ainsi que de la possibilité
d’un saut de cette dernière. Si τ ∈]0, T [, notons d’abord que (116) implique [Hu (τ )][ū(τ ] ≤ 0. On note,
pour t ∈ [0, T ], l’ensemble de contraintes du premier ordre régulières en τ par :
I1+ (t) := {j ∈ I1 (t); gj0 (ȳ(t))f1 (ȳ(t)) 6= 0}.
Lemme 8.4. Si la commande est discontinue en τ ∈]0, T [, alors ν̄j := [µ̄j (τ )] est nul, pour tout j ∈
I1+ (τ ).
Démonstration. Comme [H̄u (τ )] = [p̄(τ )] · f¯1 (τ ) = − j νj ḡj0 (τ )f¯1 (τ ), et que ḡj0 (τ )f¯1 (τ ) = 0 si
P
Exemple 8.5. Soit τ un point de jonction entre un arc singulier et un arc où la contrainte sur l’état (et
non celle sur la commande) est active. Alors [Ξ̇(τ )] = 0. Supposons g scalaire, et envisageons deux cas.
(i) Si ḡ 0 (τ )f¯1 (τ ) 6= 0, la commande étant discontinue à l’instant τ , alors ν = 0 d’après le lemme 8.4, et
˙ )]ḡ 0 (τ )f¯1 (τ ),
0 = [Ξ̇(τ )] = −[µ̄(τ
donc µ̄˙ a pour limite zéro en τ .
(ii) Si ḡ 0 (τ )f¯1 (τ ) = 0, alors
0 = [Ξ̇(τ )] = −νḡ 0 (τ ) · [f¯1 , f¯0 ](τ ).
22
Si ḡ 0 (τ ) · [f¯1 , f¯0 ](τ ) 6= 0, on en déduit que µ̄ est continue à l’instant τ .
8.4.2. Illustration : un problème proie-prédateur. Dans ce modèle inspiré de [17] (on a ajouté la contrainte
sur l’état), y1 est le nombre de proies et y2 est le nombre de prédateurs, de dynamique ẏ1 = (1 − y2 )y1 −
RT
0.9y1 u, ẏ2 = (y1 −1)y2 . On minimise 0 (1+u(t))dt, sous les contraintes u ≥ 0, y1 ≤ 2, et y(T ) = (1, 1)†
(point invariant de la dynamique non commandée). L’horizon est libre et la contrainte sur l’état est du
premier ordre. La trajectoire calculée avec le logiciel Bocop (développé dans l’équipe de l’auteur) est
constituée de 4 arcs : les premier et troisième avec commande nulle, le second avec contrainte d’état
active, et le dernier est un arc singulier. La commande ne semble discontinue qu’à la première jonction.
9. Programmation dynamique
9.1. Principe de programmation dynamique. On considère une famille de systèmes dynamiques
commandés paramétrée par (x, t) ∈ Rn × R, l’état et l’instant initial : la dynamique est
(119) ẏ(s) = f (u(s), y(s)), p.p. s ∈ [t, T ]; y(t) = x,
avec T > 0 fixé, t ∈ [0, T ], et les contraintes sur la commande
(120) u(s) ∈ Uad , p.p. s ∈ [t, T ].
On supposera dans le reste de l’article Uad compact, et f lipschitzienne. Considérons la valeur de Bellman
du problème de minimisation suivant, paramétré par l’état et le temps initiaux :
Z T
(121) V (x, t) := inf `(u(s), y(s))ds + ϕ(y(T )) soumis à (119)-(120),
u,y t
∞
où ` et ϕ sont continues. Ici u ∈ L (t, T )m et y ∈ W 1,∞ (t, T )n . En particulier, V (x, T ) = ϕ(x). On peut
vérifier que V est localement lipschitzienne. Posons
Uad [t, τ ] := {u ∈ L∞ (t, τ )m ; u(s) ∈ Uad p.p.}.
Le principe de programmation dynamique s’énonce ainsi, pour tout τ ∈]t, T ] :
Théorème 9.1. Soient 0 ≤ t ≤ τ ≤ T . Notant par y[u] la solution de l’équation d’état sur [t, τ ], avec
condition initiale x, on a :
Z τ
(122) V (x, t) = inf `(u(s), y[u](s))ds + V (y[u](τ ), τ ) .
u∈Uad [t,τ ] t
0
Démonstration. Notons V (x, t) le membre de droite de (122). Soit u ∈ Uad [t, τ ] ; posons y = y[u] et
RT Rτ
J(u, y) := t `(u(s), y(s))ds + ϕ(y(T )). Alors J(u, y) = t `(u(s), y(s))ds + ∆0 , avec
Z T
0
∆ := `(u(s), y(s))ds + ϕ(y(T )) ≥ V (y[u](τ ), τ ),
τ
23
et donc
Z τ
J(u, y) ≥ `(u(s), y(s))ds + V (y[u](τ ), τ ).
t
Minimisant en u ∈ Uad [t, τ ], on obtient V (x, t) ≥ V 0 (x, t). Réciproquement, soient ε > 0 et u ∈ Uad [t, τ ]
tel que
Z τ
0
V (x, t) + ε ≥ `(u(s), y[u](s))ds + V (y[u](τ ), τ ).
t
Par définition de V (y[u](τ ), τ ), il existe u0 ∈ Uad [τ, T ] tel que, notant par y0 la solution associée de
l’équation d’état, partant du point y[u](τ ) à l’instant τ :
Z T
V (y[u](τ ), τ ) + ε ≥ `(u0 (s), y0 (s))ds + ϕ(y0 (T )).
τ
Notant par u le recollement de u et u (fonction égale à u sur [t, τ [ et à u0 sur [τ, T ]), commande
00 0
admissible dont l’état associé y00 est le recollement de y et y0 , on vérifie en additionnant les inégalités
ci-dessus que
V 0 (x, t) + 2ε ≥ J(u00 , y00 ) ≥ V (x, t).
On conclut en faisant tendre ε vers 0.
9.2. Equation HJB. Rappelons que le préhamiltonien a pour expression H(u, y, p) := `(u, y)+p·f (u, y).
et
Z τ
(124) y[u](τ ) = x + f (u(t), x)dt + o(τ − t),
t
avec ici et dans la suite, o(τ − t)/(τ − t) → 0 quand τ ↓ t, uniformément par rapport à la commande.
Comme V est différentiable en (x, t), on a donc
Z τ
(125) V (y[u](τ ), τ ) = V (x, t) + Vt (x, t)(τ − t) + ∇x V (x, t) · f (u(t), x)ds + o(τ − t).
t
Comme l’infimum ci-dessus est atteint en minimisant séparément pour chaque t, il vaut
δ inf u∈Uad H(u, x, ∇x V (x, t)), d’où le résultat.
On appellera équation de Hamilton-Jacobi-Bellman (HJB), l’équation aux dérivées partielles (EDP)
du premier ordre suivante :
(i) vt (x, t) + inf u∈Uad H(u, x, ∇x v(x, t)) = 0, pour tout t ∈ [0, T ] et x ∈ Rn ,
(127)
(ii) v(x, T ) = ϕ(x), pour tout x ∈ Rn .
Exemple 9.3. Soit la dynamique ẏ(t) = u(t), avec |u(t)| ≤ 1, et un coût intégral nul. L’équation HJB
a pour expression vt (x, t) − |∇x v(x, t)| = 0. Si le coût final est ϕ(x) := e−|x| , la stratégie optimale est
u(t) = y(t)/|y(t)| quand y(t) 6= 0, et donc V (x, t) = et−T −|x| , qui vérifie bien l’équation HJB aux point
où elle est dérivable (hors de 0).
24
9.3. Solution de viscosité. On sait que V (x, t) est localement lipschitzienne mais en général non
différentiable partout. Une fonction lipschitzienne étant différentiable p.p. (théorème de Rademacher),
d’après le lemme 9.2, V (x, t) vérifie la relation (127)(i) p.p. et bien entendu (127)(ii), mais ce n’est pas la
seule fonction qui vérifie cette propriété. Il est donc nécessaire de donner un sens plus restreint qui assure
l’unicité de la solution ; c’est l’objet de la théorie des EDP au sens de viscosité. Pour nous conformer aux
conventions de signe de cette théorie, réécrivons l’équation HJB sous la forme
(i) −vt (x, t) + H(x, ∇x v(x, t)) = 0 pour tout t ∈ [0, T ] et x ∈ Rn ,
(128)
(ii) v(x, T ) = ϕ(x), pour tout x ∈ Rn ,
le hamiltonien H(·) étant défini par
H(x, p) := sup (−`(u, x) − p · f (u, x)) .
u∈Uad
Définition 9.4. Soit v une fonction continue de Rn × R vers R. On dit que v est sous solution (resp.
sur solution) de (128)(i) au sens de viscosité si, pour tout (x, t) ∈ Rn ×]0, T [, et toute fonction Θ :
Rn × R → R de classe C 1 , telle que v − Θ a un point de maximum (resp. minimum) en (x, t), on a
−Θt (x, t) + H(x, ∇x Θ(x, t)) ≤ 0 (resp. −Θt (x, t) + H(x, ∇x Θ(x, t)) ≥ 0). On dit que v est solution de
(128)(i) au sens de viscosité si c’est une sur et sous solution de (128).
Remarque 9.5. Si v(x, t) est de classe C 1 au voisinage de (x0 , t0 ), et vérifie l’équation HJB en ce point,
alors elle vérifie aussi en ce point les inégalités de la définition ci-dessus.
Théorème 9.6. La fonction valeur V est solution de viscosité de (128)(i).
Démonstration. (a) Montrons que V est une sous solution. Soient (x̄, t̄) et Θ comme dans la définition
de sous solution ; on peut supposer que V (x, t) = Θ(x, t), et donc V (x0 , t0 ) ≤ Θ(x0 , t0 ) si (x0 , t0 ) est voisin
de (x, t). Soit τ ∈]t, T ]. Si δ = τ − t est assez petit, comme Uad est compact, on a V (y[u](τ ), τ ) ≤
Θ(y[u](τ ), τ ) pour toute commande u ∈ Uad [t, τ ], donc d’après le principe de programmation dynamique
Z τ
(129) Θ(x, t) = V (x, t) ≤ inf `(u(s), y[u](s))ds + Θ(y[u](τ ), τ ) .
u∈Uad [t,τ ] t
1
Comme Θ est de classe C , reprenant le principe des calculs de la preuve du lemme 9.2, mais avec Θ à
la place de V et des inégalités au lieu d’égalités, on déduit que −Θt (x, t) + H(x, ∇x Θ(x, t)) ≤ 0 comme
il fallait le montrer. La preuve de sur solution est analogue.
Soient F1 et F2 deux ensembles de fonctions sur Rn × [0, T ]. On dit que l’équation HJB vérifie un
principe de comparaison relatif à (F1 , F2 ) si, pour toute sous solution u ∈ F1 et toute sur solution v ∈ F2
de (128)(i), avec u(x, T ) ≤ v(x, T ) pour tout x ∈ Rn , on a u(x, t) ≤ v(x, t), pour tout (x, t) ∈ Rn × [0, T [.
Comparant les éléments de F1 et F2 à V , on obtient :
Lemme 9.7. On suppose que l’équation HJB vérifie un principe de comparaison relatif à (F1 , F2 ), avec
V ∈ F1 ∩ F2 . Alors V est la seule solution au sens de viscosité de (128)(i), égale à ϕ au temps T , dans
l’ensemble F1 ∩ F2 .
Reste à établir des principes de comparaison. Donnons l’idée motivant une telle propriété en consi-
dérant le cas d’un sous solution u et d’une sur solution v, toutes deux de classe C 1 . Posons Φ(x, t) :=
et (u(x, t) − v(x, t)). Si u n’est pas toujours inférieure à v, le supremum de Φ est strictement positif ;
supposons-le atteint en (x, t). Alors t < T et Φ a une dérivée nulle en (x, t), soit ux (x, t) = vx (x, t) et
(130) u(x, t) − v(x, t) + ut (x, t) − vt (x, t) = 0.
D’après la remarque 9.5, on a
vt (x, t) ≤ H(x, ∇x v(x, t)) = H(x, ∇x u(x, t)) ≤ ut (x, t),
donc avec (130), u(x, t) ≤ v(x, t), ce qui donne la contradiction recherchée.
Dans les faits les semi solutions ne sont pas C 1 et Φ(x, t) n’atteint pas nécessairement son maximum.
Nous admettrons donc le résultat suivant, renvoyons à [3] pour les détails. On note BU C(Rn ) l’espace
des fonctions bornées et uniformément continues sur Rn , et on supposera que
(131) f , `, ϕ sont bornées et lipschitziennes.
Théorème 9.8. Si (131) est satisfait, l’équation HJB vérifie un principe de comparaison avec F1 =
F2 = BU C(Rn ), et donc, la fonction valeur V est l’unique élément de BU C(Rn ), solution de viscosité
de (128)(i), égal à ϕ quand t = T .
25
9.4. Problèmes en horizon infini. Soit le problème en horizon infini de fonction valeur
Z ∞
(132) V (x) := inf `(u(t), y[u, x](t))e−λt dt.
u∈Uad [0,∞] 0
Ici λ > 0 est le coefficient d’actualisation, et y[u, x] est solution de l’équation différentielle (bien définie
quand (131) est satisfait)
ẏ(t) = f (u(t), y(t)), p.p. t ∈]0, ∞[; y(0) = x.
Théorème 9.9. Si (131) est satisfait, la fonction valeur du problème en horizon infini est l’unique
élément de BU C(Rn ), solution de viscosité de l’équation
(133) λv(x) + H(x, ∇x v(x)) = 0, pour tout x ∈ Rn .
9.5. Problèmes de temps d’arrêt. On considère le problème de temps d’arrêt de fonction valeur
"Z #
θ
−λt −λθ
(134) V (x) := inf `(u(t), y[u, x](t))e dt + ϕ(y[u, x](θ))e .
(u,θ) 0
Ici λ > 0 est le coefficient d’actualisation, θ ∈]0, ∞[ est le temps d’arrêt, ϕ(·) est le coût d’arrêt,
u ∈ Uad [0, θ] est la commande, et y[u, x] est solution de
ẏ(t) = f (u(t), y(t)), p.p. t ∈]0, θ[; y(0) = x.
Théorème 9.10. Si (131) est satisfait, la fonction valeur du problème de temps d’arrêt est l’unique
élément de BU C(Rn ), solution de viscosité du problème d’obstacle
(135) max (λv(x) + H(x, ∇x v(x)), v(x) − ϕ(x)) = 0, pour tout x ∈ Rn .
Le terme d’obstacle est dû au fait que certains problèmes d’obstacle en mécanique se modélisent par
des équations similaires. Faisant θ ↓ 0 dans (134), on vérifie que v(x) ≤ ϕ(x) ; donc (135) (pris au sens
ponctuel, et non de viscosité) revient à dire que
λv(x) + H(x, ∇x v(x)) ≤ 0, avec égalité si v(x) < ϕ(x).
Intuitivement, la première (seconde) relation exprime la sous optimalité (l’optimalité) de ne pas faire
d’arrêt au point x.
9.6. Commande impulsionnelle. Dans cette modélisation la décision combine à une commande, la
possibilité de sauts de l’état, de coût positif. Plus précisément, on se donne une suite strictement crois-
sante, et tendant vers l’infini, de temps d’arrêt θi > 0, i ≥ 1 ; on pose θ0 := 0. La dynamique est
(136) ẏ(t) = f (u(t), y(t)), θi < t < θi+1 , y(θi+ ) = y(θi− ) + ξi , i ≥ 1, y(0) = x.
La décision est donc (u, θ, ξ) avec u ∈ Uad [0, ∞], et la valeur associé est
Z ∞ X
(137) V (x) = inf `(u(t), y(t))e−λt dt + (c0 + c(ξ))e−λθi .
(u,θ,ξ) 0 i≥1
Mutipliant par h0 , ajoutant vj de chaque côté et posant β := (1 + h0 λ)−1 , on obtient la forme de point
fixe
(146)
vj = β min h0 `(u, xj ) + (1 − h0 h−1 −1 −1
1 |f (u, xj )|)vj + h0 h1 f (u, xj )+ vj+1 − h0 h1 f (u, xj )− vj−1 .
u∈Uad
27
On suppose (141) satisfaite. Notons B(Z) l’espace des fonctions bornées sur Z, muni de la norme de
la convergence uniforme. On peut alors vérifier que le membre de droite ci-dessus définit un opérateur
contractant (de facteur β < 1) de B(Z) sur lui même.
Le schéma a donc une solution unique (qui ne dépend pas de h0 ) qu’on peut calculer en itérant cet
opérateur de point fixe : c’est l’approche par itération sur les valeurs. On peut également résoudre (146)
en alternant le calcul du minimum en u, pour v fixé, et la résolution du système linéaire à (uj ), j ∈ Z
fixé, soit
vj = β h0 `(uj , xj ) + (1 − h0 h−1 −1 −1
1 |f (uj , xj )|)vj + h0 h1 f (uj , xj )+ vj+1 − h0 h1 f (uj , xj )− vj−1 ,
pour tout j ∈ Z. C’est l’approche par itération sur les stratégies, ou algorithme de Howard. Cette seconde
approche nécessite (après troncature du domaine de calcul) la résolution de grands systèmes linéaires.
10.4. Problèmes de temps d’arrêt. On reprend le problème (134), quand n = 1. Soient encore h0 > 0
et β := (1 + h0 λ)−1 . L’équation HJB (135) peut se réécrire comme
v(x) = min(β(v(x) − h0 H(x, ∇x v(x, t))), ϕ(x)),
= min(β(v(x) + h0 inf H(u, x, ∇x v(x, t))), ϕ(x)).
u∈Uad
Discrétisant cette relation dans l’esprit du cas de l’horizon infini, on obtient le schéma
(147)
h0 h0 h0
vj = min β min h0 `(u, xj ) + (1 − f (u, xj )|)vj + f (u, xj )+ vj+1 − f (u, xj )− vj−1 , ϕ(xj ) .
u∈Uad h1 h1 h1
Le membre de droite correspond encore à un point fixe contractant quand (141) est satisfait, et le schéma
est donc bien posé.
où I k+1 est une formule d’interpolation des valeurs de v k+1 aux points de la grille. Si
X
I k+1 (x) = αq (x)vqk+1
q∈Zn
avec les αq (x) positifs et de somme un, le schéma est monotone et (142) est encore satisfait. On peut
raffiner le schéma en introduisant des schémas d’intégration temporelle et d’interpolation spatiale d’ordre
plus élevés, au risque de perdre la monotonie.
10.6. Bilan sur les méthodes numériques pour l’équation HBJ. La minimisation par rapport
à u s’effectue habituellement en discrétisant les valeurs des commandes. Ces algorithmes convergent
uniformément sur tout compact vers la fonction valeur. Cependant ils deviennent coûteux quand la
dimension (nombre de variables d’état et de commandes) augmente. On les appliquera donc si possible
à une version simplifiée du problème de commande optimale, dont la dimension est réduite.
Références
[1] M.S. Aronna, J.F. Bonnans, and P. Martinon. A shooting algorithm for optimal control problems with singular arcs.
J. Optim. Theory Appl., 158(2) :419–459, 2013.
[2] M. Bardi and I. Capuzzo-Dolcetta. Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations.
Birkhäuser, Boston, MA, 1997.
[3] G. Barles. Solutions de viscosité des équations de Hamilton-Jacobi. Springer-Verlag, Paris, 1994.
[4] T. Bayen, A. Rapaport, and M. Sebbah. Minimal time control of the two tanks gradostat model under a cascade inputs
constraint. SIAM J. Control Optim, 2014. In print.
[5] R. Bellman. Dynamic programming. Princeton University Press, Princeton, N. J., 1957.
[6] J.T. Betts. Practical methods for optimal control and estimation using nonlinear programming. SIAM, Philadelphia,
PA, second edition, 2010.
[7] L.T. Biegler. Nonlinear programming. SIAM, Philadelphia, PA, 2010.
[8] J.F. Bonnans, X. Dupuis, and L. Pfeiffer. Second-order sufficient conditions for strong solutions to optimal control
problems. ESAIM-COCV, 20(3) :704–724, 2013.
[9] J.F. Bonnans and P. Rouchon. Commande et optimisation de systèmes dynamiques. Editions de l’Ecole Polytechnique,
Palaiseau, 2005.
[10] J.F. Bonnans and A. Shapiro. Perturbation analysis of optimization problems. Springer-Verlag, New York, 2000.
[11] B. Bonnard, L. Faubourg, and E. Trélat. Mécanique céleste et contrôle des véhicules spatiaux. Springer-Verlag, Berlin,
2006.
[12] B. Bonnard and D. Sugny. Optimal control with applications in space and quantum dynamics. American Institute of
Mathematical Sciences (AIMS), Springfield, MO, 2012.
[13] I. Ekeland. Nonconvex minimization problems. Bull. Amer. Math. Soc. (N.S.), 1(3) :443–474, 1979.
[14] M. Falcone and R. Ferretti. Semi-Lagrangian Approximation Schemes for Linear and Hamilton–Jacobi Equations.
SIAM Publications, 2013.
[15] A. T. Fuller. Study of an optimum non-linear control system. J. Electronics Control (1), 15 :63–71, 1963.
[16] B.S. Goh. Management and analysis of biological populations. Elsevier, 1980.
[17] B.S. Goh, G. Leitmann, and T.L. Vincent. Optimal control of a prey-predator system. Math. Biosci., 19 :263–286,
1974.
[18] U. Ledzewicz, H. Schättler, A. Friedman, and E. Kashdan, editors. Mathematical methods and models in biomedicine.
Springer, New York, 2013.
[19] L.S. Pontryagin, V.G. Boltyanskiı̆, R.V. Gamkrelidze, and E.F. Mishchenko. The mathematical theory of optimal
processes. Gordon & Breach Science Publishers, New York, 1986. Reprint of the 1962 English translation.
[20] H. Schättler and U. Ledzewicz. Geometric optimal control. Springer, New York, 2012.
[21] Baratchart (L.) et Bernhard (P.). - Automatique et systèmes [A 1370]. Base documentaire "Mathématiques
pour l’ingénieur” (1992).
[22] Borne (P.) et Rotella (F.). - Commande optimale [R 7427]. Base documentaire ”Régulation et commande
des systèmes asservis” (1996).
[23] Brezinski (C.). - Aspects numériques du contrôle linéaire [AF 1400]. Base documentaire "Mathématiques pour
l’ingénieur” (2007).
[24] Dacorogna (B.). - Calcul des Variations [AF 111]. Base documentaire "Mathématiques pour l’ingénieur” (2007).
[25] Lévine (J.) et Rouchon (P.). Systèmes dynamique et commande [S 7430]. Base documentaire "Mathématiques
pour l’ingénieur” (2008).
[26] Yvon (J.P.). - Contrôle des systèmes à paramètres distribués [AF 1372]. Base documentaire "Mathématiques
pour l’ingénieur” (2009).