Equation de Hamilton-Jacobi-Bellman et fonction valeur

Méthode directe : optimisation non-linéaire sous contraintes 

La méthode directe est la méthode la plus évidente lorsque l’on doit résoudre un problème de commande optimale. En discrétisant l’état et la commande dans le problème (3.1), on se ramène à un problème d’optimisation non-linéaire en dimension finie (N variables) :    min uk∈Uk J(u) := N−1 ∑ k=0 Lk(uk)∆t +φ(xN) avec : xk+1 = fk(xk,uk), x(0) = x0. xmin ≤ x(t) ≤ xmax umin(t) ≤ u(t) ≤ umax(t), (3.2) où xk+1 = fk(xk,uk) représente la version discrète de l’équation d’état x˙ = f(u(t),t), discrétisée en utilisant par exemple la méthode d’Euler explicite. Pour résoudre ce problème, on suppose que le temps est subdivisé de manière égale, tel que 0 = t0 < t1 < … < tN = T, le pas de discrétisation étant noté ∆t. On suppose que la commande reste constante par morceau durant le pas de temps ∆t. Les contraintes sur la commande ou sur l’état sont appliquées sur les valeurs discrétisées. 3.1.1 Résolution numérique La résolution de ce problème peut s’effectuer par exemple via un algorithme de type SQP (Sequential Quadratic Programming), [Nocedal and Wright, 2000]. Cette méthode a été implémentée, mais la mise en oeuvre du problème et la prise en compte des contraintes est assez complexe et contraignante, et le temps nécessaire à la convergence est important. Pour ces raisons, la programmation non-linéaire n’a pas été retenue pour l’étude paramétrique présentée dans la section 4. 

Paramétrisation simplifiée du contrôle (découpage par zones)

 Dans la précédente méthode, les paramètres du contrôle u sont les uk, k = 0,…,N −1. Au lieu de rechercher les N commandes optimales, problème d’optimisation complexe si N est élevé, on peut imposer des règles pour la définition de u, règles définies par des paramètres dont le nombre serait Nr ¿ N. La difficulté est alors de trouver empiriquement des règles intuitives permettant de définir la commande u d’une façon simple, à l’aide d’un nombre restreint de paramètres. Dans l’application étudiée, cela pourrait revenir à définir des seuils de couples et/ou de régime audessus desquels certaines valeurs de u seraient figées. Plutôt que d’optimiser une commande en fonction du temps, on peut imposer des zones de fonctionnement (voir [Voise et al., 2005]) pour le moteur électrique, ainsi que pour le fonctionnement simultané des deux moteurs, etc. Dans un second temps, une loi de gestion d’énergie propre à chaque zone pourra être définie, permettant une utilisation de cette stratégie en temps-réel. Cette paramétrisation permet de réduire grandement la taille du problème d’optimisation, puisque seules quelques variables devront être utilisées dans la définition du problème d’optimisation. Cette technique ne peut se faire sans connaissance du système : on définit, par exemple, une commande u˜ telle que : u˜ =    u1 si Trq < Ta, u2 si Ta ≤ Trq < Tb, u3 si Tb ≤ Trq. (3.3) La stratégie obtenue, après optimisation, reste évidemment sous-optimale, puisque la paramétrisation enlève des degrés de liberté à la commande. 3.2 Méthodes de tir Les méthodes de tir consistent à résoudre directement les conditions d’optimalité du problème (PO). Elles sont basées sur le Principe de Pontryagin. Ces méthodes sont notamment très utilisées dans le contrôle des véhicules spatiaux, où la précision est primordiale. Elles ont pour avantage d’être numériquement très précises, et relativement rapides comparées aux méthodes indirectes. 3.2.1 Méthode de tir simple La méthode de tir consiste à itérer sur la valeur de p pour atteindre une valeur objectif. On définit une fonction de tir, notée Γ, telle que Γ(p0) = p(T)−φ 0 (x(T),T). (3.4) La résolution des conditions d’optimalité revient alors à ajuster la valeur de p0 telle que Γ(p0) = 0. (3.5) En pratique, pour résoudre ce type de problème, on utilisera une méthode numérique itérative, qui consiste à itérer sur la valeur de p0 lorsque l’on possède une expression analytique pour les équations (2.3b) et (2.3d). Cette méthode convient donc bien à la résolution de problèmes dont la forme est analytique, permettant de calculer une forme elle aussi analytique pour l’expression de u ∗ . Lorsque la valeur initiale du paramètre p0 est bien choisie, sa convergence est assez rapide. Néanmoins, cette méthode n’est plus valide en présence de contraintes actives sur l’état. En effet, l’état adjoint p ne suit plus la loi d’évolution définie par (2.3b) lorsque l’état x touche ponctuellement ou longe une contrainte, puisqu’un nouveau multiplicateur associé à la contrainte sur l’état doit être introduit.

Méthode de tir multiple

Dans certains problèmes d’optimisation, le système différentiel régissant l’évolution des variables d’état n’est pas le même au cours de la trajectoire. L’évolution des variables d’état pour ce type de problème s’écrit par exemple : (x˙(t), p˙(t)) =    f0(t,x(t), p(t)) si 0 ≤ t ≤ t1 f1(t,x(t), p(t)) si t1 ≤ t ≤ t2 . . . fN(t,x(t), p(t)) si tN ≤ t ≤ T, (3.6) x étant l’état, p étant l’état adjoint, et t1,t2,…tN ∈ [0,T] étant des temps de commutation. Lorsqu’il est possible de dériver une expression analytique de la commande optimale, la méthode de tir multiple ([Bryson and Ho, 1975], [Evans, 2000]) se trouve être bien adaptée. Elle consiste à résoudre 45 Chapitre 3. Méthodes numériques le problème d’optimisation comme dans la méthode de tir simple, avec des conditions de continuité (ou conditions de jonction) sur les variables d’état et d’état adjoint, lors d’un changement du système différentiel. Ce type d’algorithme se trouve donc bien adapté à la résolution de problèmes d’optimisation avec contraintes d’état. En effet, dans ce cas les temps t1,t2,…tN ∈ [0,T] peuvent être des temps de jonction avec un arc frontière (contrainte sur l’état active pendant un temps donné), ou bien des temps de contact avec la frontière (l’état touche une contrainte, mais la quitte aussitôt). Un très bon exemple se trouve dans [Bonnans and Hermant, 2008], exemple qui est repris dans [Rousseau et al., 2008c], voir section (3.4.6) : il s’agit de connaître la trajectoire optimale d’une corde élastique attachée à ses deux extrémités, et qui repose sur un support plan. Lorsque la corde ne touche pas le support, l’équation représentant la trajectoire de la corde est parfaitement connue. En revanche, lorsque la corde touche le support, il devient difficile de connaître sa trajectoire, notamment à quels endroits elle touche le support puis le quitte. On peut alors « découper » la trajectoire en 3 tronçons, chaque partie étant décrite par une équation différentielle différente. Les temps de commutation représentent alors les abscisses auxquels la corde touche puis quitte le support. Dans certains problèmes très simples, on peut résoudre le problème, en écrivant les conditions d’optimalité associées. Différents cas sont énumérés dans [Bryson and Ho, 1975] : système continu avec état final fixé à un temps donné, état final fixé à un temps final libre (incluant les problèmes à temps minimal), etc. L’auteur résout les problèmes posés en adjoignant l’équation d’état et les contraintes associées à l’état, à la fonction à minimiser, puis définit les conditions d’optimalité associées au problème via une méthode de perturbation (telle que celle qui est utilisée dans (2.8)). Le lecteur pourra aussi se référer à [Sethi and Thompson, 2006] pour des applications relatives aux problèmes économiques. Dans le cas de contraintes du premier ordre sur l’état, [Bonnans and Hermant, 2008] a prouvé que l’état adjoint était continu sous certaines hypothèses. Pour des contraintes d’ordre supérieur cependant, des « conditions de saut » sur l’état adjoint doivent être considérées dans le problème d’optimisation, la continuité de celui-ci n’étant pas assurée. L’inconvénient majeur de ces méthodes provient du fait qu’il soit nécessaire de connaître la forme de la trajectoire optimale pour pouvoir intégrer une équation différentielle valide, ce qui revient à connaître le nombre de contraintes actives.

Approche par l’équation de Hamilton-Jacobi-Bellman 

Equation de Hamilton-Jacobi-Bellman et fonction valeur

Considérons à nouveau le problème d’optimisation général (3.1), qu’on englobe dans la classe des problèmes avec état et instantsinitiaux (x,t) quelconques, où x suit la loi d’évolution x˙(t) = f(x(t),u(t),t). On définit la fonction valeur (ou fonction coût) V(x(t),t) V(x(t),t) := min u∈U ·Z T t L(x(τ),u(τ), τ)dτ+φ(x(T),T) ¸ . (3.7) On suppose que V(x(t),t) existe et est C 1 . Soit ∆t l’intervalle de temps, d’après l’équation d’évolution on a alors : x(t +∆t) ≈ x(t)+ f(x(t),u(t),t)∆t. (3.8) 46 3.3 Approche par l’équation de Hamilton-Jacobi-Bellman Considérons qu’une commande quelconque soit appliquée durant le laps de temps ∆t. Alors on aura au premier ordre : V(x(t),t) ≤ V(x(t)+ f(x(t),u(t),t)∆t,t +∆t)+L(x(t),u(t),t)∆t. (3.9) L’égalité dans (3.9) n’est obtenue que si une commande optimale est appliquée durant l’intervalle de t à t +∆t. On peut donc écrire : V(x(t),t) = min u∈U {V(x(t)+ f(x(t),u(t),t)∆t,t +∆t)+L(x(t),u(t),t)∆t} (3.10) Par un développement de Taylor au premier ordre, on obtient : V(x(t),t) = min u∈U {V(x(t),t)+ ∂V ∂x (x(t),t)f(x(t),u(t),t)∆t + ∂V ∂t (x(t),t)∆t +L(x(t),u(t),t)∆t}. (3.11) Comme V et ∂V ∂t ne dépendent pas de u, on peut écrire, lorsque ∆t → 0 − ∂V ∂t (x(t),t) = min u∈U h L(x(t),u(t),t)+ ∂V ∂x (x(t),t)f(x(t),u(t),t) i . (3.12) L’équation (3.12) est appelée équation de Hamilton-Jacobi-Bellman (ou équation HJB). Elle correspond à l’équation représentant le problème de minimisation (PO) (défini dans la section 2). La fonction valeur V est ainsi solution d’une équation aux dérivées partielles où apparaissent l’équation du système x˙(t) = f(x,u,t) et la fonction à minimiser L(x,u,t). Trouver V(x,t) dans (3.7) revient donc à résoudre l’équation (3.12). 

Programmation Dynamique

 La programmation dynamique est fréquemment utilisée pour résoudre les problèmes du type (3.1) (voir [Wu et al., 2002], [Scordia, 2004], [Sundström et al., 2008b] pour des applications sur véhicules hybrides) : cette méthode repose sur le principe d’optimalité de Bellman, énoncé par R. Bellman : « Une suite de commandes optimales est telle que, quels que soient l’état et l’instant considérés sur une trajectoire optimale, les commandes ultérieures constituent pour le problème ayant cet état et cet instant comme éléments initiaux une suite de commandes optimales. » Le principe d’optimalité de Bellman indique qu’une commande optimale peut être construite séquentiellement, tout d’abord en calculant la commande optimale pour le dernier pas de temps, ensuite pour les deux derniers, puis pas à pas jusqu’à considérer le problème dans sa globalité. Principe La programmation dynamique est donc une méthode d’optimisation en temps rétrograde, qui permet de calculer une approximation de la fonction V solution de l’équation HJB, discrétisée en espace et en temps. Pour cela, on calcule de façon séquentielle la trajectoire optimale entre deux instants, pour un état initial donné. On définit le Principe de Programmation Dynamique, valide ∀x ∈ R et ∀tk,tk+1 tels que tk < tk+1 ≤ T, avec tk+1 = tk +∆t par V(x,tk) = min u∈U {∆tL(x,u,tk)+V(x+∆t f(x,u,tk),tk+1)}. (3.13) et la fonction V(x,tk) s’écrit, au temps tk = T V(x,T) = φ(xN). (3.14) 47 Chapitre 3. Méthodes numériques Justification de l’équation (3.13) Pour tout (x0,t0) fixé et tout u0 ∈ U donné, on considère la trajectoire t 7→ x˜(t;x0,t0,u0) définie par le système différentiel ( d dt x˜(t; x0,t0,u0) = f(x˜(t; x0,t0,u0),u0,t0) x˜(t = t0; x0,t0,u0) = x0 (3.15) dont la solution existe au voisinage de t0. Soit la fonction composée V˜(t; x0,t0,u0) = V(x˜(t;x0,t0,u0),t) qui représente la valeur de V sur la trajectoire x˜. Alors, pour tout t, on a d dt V˜(t; x0,t0,u0) = ∂ ∂t V(x˜(t;x0,t0,u0),t)+ ∂ ∂x V(x˜(t; x0,t0,u0),t) d dt x˜(t;x0,t0,u0) = ∂ ∂t V(x˜(t;x0,t0,u0),t)+ ∂ ∂x V(x˜(t; x0,t0,u0),t)f(x˜(t;x0,t0,u0),t0,u0). En spécifiant t = t0, on obtient d dt V˜(t = t0;x0,t0,u0) = ∂ ∂t V(x0,t0)+ ∂ ∂x V(x0,t0)f(x0,t0,u0). Cette relation est valide pour tout (x0,t0,u0). La fonction composée V˜ nous aide à discrétiser l’équation HJB en temps. En effet, entre l’instant tk et tk+1 = tk +∆t, on a, pour tout u ∈ U ∂ ∂t V(x,tk)+ f(x,u,tk) ∂V ∂x (x,tk) = d dt V˜(t = tk;x,tk,u), qu’on peut approcher par la différence finie V˜(t = tk+1;x,tk,u)−V˜(t = tk; x,tk,u) ∆t . Or, il est immédiat que V˜(t = tk;x,tk,u) = V(x,tk,u). D’autre part V˜(t = tk+1; x,tk,u) = V(x˜(t = tk+1;x,tk,u),tk+1), et on peut approcher la position x˜(t = tk+1,x,tk,u) par sa version discrète x˜(t = tk+1; x,tk,u)−x ∆t = f(x,u,tk) du système différentiel (3.15) pour x0 = x, t0 = tk. Finalement, on a l’approximation ∂V ∂t (x,tk)+ f(x,u,tk) ∂V ∂x (x,tk) ≈ V(x+∆t f(x,u,tk),tk+1)−V(x+∆t f(x,u,tk),tk+1) ∆t . 48 3.3 Approche par l’équation de Hamilton-Jacobi-Bellman En reportant cette approximation dans l’équation HJB, et en multipliant par ∆t, on aboutit à min u∈U {V(x+∆t f(x,u,tk),tk+1)−V(x,tk)+∆tL(x,u,tk)} = 0. Comme V(x,tk) ne dépend pas de u, on peut réécrire la relation précédente sous la forme d’un schéma rétrograde V(x,tk) = min u∈U {V(x+∆t f(x,u,tk),tk+1)+∆tL(x,u,tk)}. (3.16) On retrouve ainsi le principe de Programmation Dynamique.

Méthode de résolution 

Pour l’application de la programmation dynamique, on doit résoudre l’équation (3.13) à chaque pas de temps. On se place maintenant dans l’ensemble [xmin,xmax] ⊂ R. Il n’est cependant pas possible de résoudre cette équation ∀x ∈ [xmin,xmax]. Aussi, on a recours à une discrétisation de l’espace d’état, l’équation (3.13) étant alors évaluée à chaque pas de temps, et à chaque pas d’espace. Cette discrétisation dans le cas d’un problème à un seul état impose de parcourir une grille, l’axe horizontal représentant le temps, l’axe vertical l’espace. Pour chaque point de la grille, un certain nombre de commandes admissibles sont testées, la commande optimale retenue étant celle qui minimise le terme de droite dans (3.13). L’algorithme de programmation dynamique s’écrit : Algorithme 1: Programmation Dynamique Données: xmin, xmax, ∆x, ∆t, x0, T Définir V(x,T) = β(x0 −x) 2 début pour tk = (T −∆t) : −∆t : 0 faire pour x = xmin : ∆x : xmax faire Calculer u ∗ (x,tk) = argmin u∈U {V(x+∆t f(x,u,tk),tk+1)+∆tL(x,u,tk)} fin fin x(0) = x0 ; pour tk = 0 : ∆t : (T −∆t) faire x(tk+1) = f(x(tk),u ∗ (x(tk),tk),tk) fin fin Cet algorithme comporte donc une partie backward, qui résout le problème d’optimisation défini par (3.13) un certain nombre de fois (voir section 3.3.4) du tempst = T −1 au tempst = 0. Cette étape terminée, l’ensemble des commandes optimales sont connues pour chaque noeud de la grille de discrétisation. Il suffit ensuite de parcourir la grille dans le sens forward, de l’état au temps initial x(0) = x0 jusqu’au temps final x(T), en sélectionnant la commande optimale du noeud près duquel on se trouve. Le calcul deV(x+∆t f(x,u,tk),tk+1) nécessite l’utilisation d’un schéma numérique. Nous choisissons un schéma explicite, tel qu’il est décrit dans [Guilbaud, 2002] : le schéma décentré ou schéma upwind.