Programmation dynamique
Macros générales

Michel de lara
(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 Programmation dynamique stochastique sur des entiers (cas des matrices de transition)
2 Programmation dynamique stochastique sur des entiers (cas dynamique bruitée)
3 Transcription en objets Scilab
4 Macros Scilab pour la résolution numérique (cas des matrices de transition)
4.1 Résolution de l’équation de la programmation dynamique stochastique
4.2 Simulation de trajectoires optimales
5 Macros Scilab pour la résolution numérique (cas dynamique bruitée)
6 Programmation dynamique stochastique sur un ensemble ﬁni
7 Problèmes de discrétisation
7.1 Dynamique stochastique commandée
7.2 Discrétisation d’un problème d’optimisation stochastique sur ℝ
7.3 Discrétisation d’un problème d’optimisation stochastique sur ℝⁿ
7.4 Discrétisation d’un problème d’optimisation déterministe sur ℝ et mise sous forme stochastique
7.5 Discrétisation d’un problème d’optimisation stochastique sur ℝ^N, N ≥ 2

1 Programmation dynamique stochastique sur des entiers (cas des matrices de transition)

Le problème d’optimisation considéré ici est

(T −1 ) ∑ u0∈𝕌,i..n.,fuT−1∈𝕌𝔼 L(t,x(t),u(t)) + Φ (T, x(T )) , t=0

(1)

avec

coût instantané, L : {0,...,T − 1}×{1,...,n}×{1,...,p}→ ℝ, où
1. n ∈ ℕ^∗ est le cardinal de l’espace d’état 𝕏 := {1,...,n}⊂ ℕ^∗,
2. p ∈ ℕ^∗ est le cardinal de l’espace de commande 𝕌 := {1,...,p}⊂ ℕ^∗,
3. T ∈ ℕ^∗ est l’horizon,
coût ﬁnal, Φ : {1,...,n}→ ℝ,
matrices de transition, M : {1,...,p}→ℳ(n,n), application à valeurs dans les matrices (n,n) stochastiques, c’est-à-dire ∀u = 1,...,p, $′ 2 u ′ ∑ u ′ ∀(x, x) ∈ {1,...,n } ,M (x,x ) ≥ 0 et ∀x ∈ {1,...,n}, M (x,x ) = 1; x′∈{1,...,n}$
le terme M^u(x,x′) représente la probabilité que, partant de l’état x et appliquant la commande u à l’instant t, le futur état au temps t + 1 soit x′ :
$M u(x,x′) = ℙ(x(t + 1) = x′ | x (t) = x,u (t) = u ),$
fonction valeur
$(T − 1 ) ∑ V (t,x) := ut∈𝕌,.in..f,uT−1∈𝕌𝔼 L(s,x (s),u (s )) + Φ (T, x(T )) | x (t) = x , s=t$ (2)

la fonction valeur évaluée sur l’état x à l’instant t,
l’équation de la programmation dynamique
$( ) ( ∑ u ′ ′) V (t,x) = mui∈n𝕌 L(t,x,u ) + M (x,x )V (t + 1,x ) x′∈{1,...,n}$ (3)
feedback optimal
$( ) ∑ u ♯(t,x) := argmin (L (t,x,u) + M u(x,x′)V(t + 1,x′)) , u∈𝕌 x′∈{1,...,n}$ (4)

si l’inf est atteint en un unique argument,
trajectoires d’état, de commande, de coût
$( || (x(0),...,x(T − 1),x(T)) ||| |{ (u(0),...,u(T − 1)) ||| T∑−2 T∑− 1 |||( (L (0,x(0),u(0)),..., L (t,x(t),u (t)), L(t,x(t),u(t)) + Φ (T,x(T ))) t=0 t=0$ (5)

2 Programmation dynamique stochastique sur des entiers (cas dynamique bruitée)

Le problème d’optimisation considéré ici identique avec

coût instantané,
coût ﬁnal
dynamique bruitée, F : {0,...,T − 1}×{1,...,n}×{1,...,p}× 𝕎 →{1,...,n}, application telle que
$x (t + 1) = F (t,x (t),u (t),w (t))$ (6)

où w(0),…,w(T − 1) est une suite de variables aléatoires i.i.d.
fonction valeur
$( ) T∑− 1 V (t,x) := inf 𝔼 L(s,x (s),u (s )) + Φ (T, x(T )) | x (t) = x , ut∈𝕌,...,uT−1∈𝕌 s=t$ (7)

la fonction valeur évaluée sur l’état x à l’instant t, prise sur une loi produit ℙ sur 𝕎^T;
l’équation de la programmation dynamique
$( ) V (t,x) = mui∈n𝕌 L (t,x, u) + 𝔼w(t) (V (t + 1,F (t,x, u,w(t))))$ (8)
feedback optimal
$♯ ( ) u (t,x) := argmui∈n𝕌 L (t,x,u) + 𝔼w(t)(V(t + 1,F (t,x, u,w (t))))$ (9)

si l’inf est atteint en un unique argument,
trajectoires d’état, de commande, de coût

3 Transcription en objets Scilab

Le coût instantané est une liste cout_instantane à p éléments composée de matrices de dimension (n,T) :
$coutxinstantane (l)(i,j) = L(i,l,j − 1), l = 1,...,p, i = 1,...,n, j = 1,...,T ,$ (10)

où on notera le décalage dans les indices temporels, et où
1. le cardinal n de l’espace d’état est un entier cardinal_etat
  $cardinalxetat = size (coutxinstantane (1)) (1) = n,$ (11)
2. le cardinal p de l’espace de commande est un entier
  $cardinalxcommande = size(coutxinstantane ) = p,$ (12)
3. l’horizon T est un entier horizon
  $horizon = 1 + size (coutxinstantane (1) )(2) = T$ (13)
le coût ﬁnal est un vecteur cout_final de dimension (n, 1) :
$coutxfinal (i) = Φ(i), l = 1,...,p,$ (14)
Transitions
1. les matrices de transition sont une liste mat_transition à p éléments composée de matrices (n,n) stochastiques :
  $matxtransition (l)(i,i′) = M l(i,i′), l = 1,...,p, i = 1,...,n, i′ = 1,...,n,$ (15)
2. la dynamique bruitée est une fonction dynamique(t,i,l,w)
la fonction valeur est une matrice valeur de dimension (n,T + 1) :
$valeur (i,j) = V (i,j − 1), i = 1,...,n, j = 1,...,T + 1,$ (16)

où on notera le décalage dans les indices temporels,
le feedback optimal est une matrice feedback de dimension (n,T) :
$feedback (i,j) = u ♯(i,j − 1), i = 1,...,n, j = 1, ...,T ,$ (17)

où on notera le décalage dans les indices temporels,
les trajectoires d’état, de commande, de coût sont des vecteurs trajectoire_etat, trajectoire_commande et trajectoire_cout, de dimensions respectives (1,T + 1), (1,T) et (1,T + 1) :
$( || trajectoirexetat (j) = x(j − 1), j = 1,...,T + 1 ||| ||| trajectoirexcommande (j) = u(j − 1), j = 1,...,T |||{ j∑−1 | trajectoirexcout (j) = L (x(t),u(t),t), j = 1,...,T |||| t=0 ||| ∑T ||| trajectoirexcout (T + 1) = L (x(t),u(t),t) + Φ (x (T)) ( t=0$ (18)

4 Macros Scilab pour la résolution numérique (cas des matrices de transition)

4.1 Résolution de l’équation de la programmation dynamique stochastique

L’astuce employée ici consiste à utiliser un indice dans une matrice pour représenter l’état.

  function [valeur,feedback]=Bell_stoch(matrice_transition,cout_instantane,cout_final, ...
                                        miniOUmaxi)
    // algorithme de programmation dynamique stochastique
    //     pour une chaîne de Markov sur \{1,2...,cardinal_etat\}
    // matrice_transition = liste à cardinal_commande éléments
    //     composée de matrices de dimension (cardinal_etat,cardinal_etat)
    // cout_instantane = liste à cardinal_commande éléments
    //     composée de matrices de dimension (cardinal_etat,horizon-1)
    // cout_final = vecteur de dimension (cardinal_etat,1)
    // valeur = fonction valeur de Bellman,
    //     matrice de dimension (cardinal_etat,horizon)
    // feedback = commande en feedback sur l'état et le temps
    //     matrice de dimension (cardinal_etat,horizon-1)
    // MINIMISATION si miniOUmaxi=-1
    // MAXIMISATION si miniOUmaxi=1

    ee=size(cout_instantane(1));cardinal_etat=ee(1);
    cardinal_commande=size(cout_instantane);
    hh=size(cout_instantane(1));horizon=1+hh(2);

    valeur=0*ones(cardinal_etat,horizon);
    // tableau ou l'on stocke la fonction de Bellman
    // valeur(:,t) est un vecteur
    // : représente ici le vecteur d'état

    valeur(:,horizon)=cout_final;
    // La fonction valeur au temps T est cout_final

    feedback=0*ones(cardinal_etat,horizon-1);
    // tableau ou l'on stocke le feedback u(x,t)

    // Calcul rétrograde de la fonction de Bellman et
    // de la commande optimale à la date 'temp' par l'équation de Bellman
    // On calcule le coût pour la commande choisie connaissant la valeur de la
    // fonction coût à la date suivante pour l'état obtenu

    for temp=horizon:-1:2 do
      // équation rétrograde
      // (attention, pour des questions d'indices, bien finir en :2)
      loc=zeros(cardinal_etat,cardinal_commande);
      // variable locale contenant les valeurs de la fonction à minimiser
      // loc est une matrice
      for j=1:cardinal_commande do
        loc(:,j)=matrice_transition(j)*valeur(:,temp)+cout_instantane(j)(:,temp-1);
      end;

      [mmn,jjn]=mini(loc,'c');
      [mmx,jjx]=maxi(loc,'c');
      // mm est l'extremum atteint
      // jj est la commande qui réalise l'extremum

      valeur(:,temp-1)=(1-miniOUmaxi)/2*mmn+(1+miniOUmaxi)/2*mmx;
      // mm;
      // coût optimal
      feedback(:,temp-1)=(1-miniOUmaxi)/2*jjn+(1+miniOUmaxi)/2*jjx;
      //jj;
      // feedback optimal
    end
  endfunction

4.2 Simulation de trajectoires optimales

  function z=trajopt(matrice_transition,feedback,cout_instantane,cout_final,etat_initial)
    // z est une liste :
    // z(1) est la trajectoire optimale obtenue partant de etat_initial
    // z(2) est la trajectoire correspondante des commandes optimales
    // z(3) est la trajectoire correspondante des coûts
    // etat_initial = un entier
    // pour le reste, voir la macro Bell_stoch

    ee=size(cout_instantane(1));cardinal_etat=ee(1);
    cardinal_commande=size(cout_instantane);
    hh=size(cout_instantane(1));horizon=1+hh(2);

    z=list();
    x=etat_initial;
    u=[];
    c=[];

    for j=1:horizon-1 do
      // remplissage des trajectoires : commande, coût, état
      u=[u,feedback(x($),j)];
      c=[c,c($)+cout_instantane(u($))(x($),j)];
      mat_trans=full(matrice_transition(u($)));
      x=[x,grand(1,'markov',mat_trans,x($))];
    end
    c=[c,c($)+cout_final(x($))];

    z(1)=x;z(2)=u;z(3)=c;
  endfunction

5 Macros Scilab pour la résolution numérique (cas dynamique bruitée)

6 Programmation dynamique stochastique sur un ensemble ﬁni

Le problème d’optimisation est toujours (1), mais avec

coût instantané, L : {x¹,...,xⁿ}×{u¹,...,u^p}×{0,...,T}→ ℝ, où
1. 𝕏 := {x¹,...,xⁿ} est l’espace d’état,
2. 𝕌 := {u¹,...,u^p} est l’espace de commande,
3. T ∈ ℕ^∗ est l’horizon,
coût ﬁnal, Φ : {x¹,...,xⁿ}→ ℝ,
matrices de transition
$M ul(xi,xi′) = ℙ(x(t + 1) = xi′ | x(t) = xi,u(t) = ul).$ (19)

On se ramène au cas des entiers en introduisant

nouveau coût instantané, : {1,...,n}×{1,...,p}×{0,...,T − 1}→ ℝ, avec
$i l ^L (i,l,t) := L (x ,u ,t), l = 1,...,p, i = 1,...,n, t = 0,...,T − 1,$ (20)
nouveau coût ﬁnal, : {1,...,n}→ ℝ, avec
$^Φ(i) := Φ (xi), i = 1,...,n,$ (21)
transition
1. nouvelles matrices de transition déﬁnies par
  $^M l(i,i′) := M ul(xi,xi′).$ (22)
2. dynamique bruitée...

On transcrit ces données , et en objets Scilab cout_instantane, cout_final, matrice_trans. L’exécution de la macro Bell_stoch retourne en particulier feedback.

Pour la simulation de trajectoires optimales, il suﬃt alors d’introduire un vecteur etat de dimension (1,n) et tel que

etat (i) = xi, i = 1,...,n,

(23)

et un vecteur commande de dimension (1,p) et tel que

commande (l) = ul, l = 1,...,p,

(24)

et d’exécuter une instruction du type

  etat_initial=i;
  // correspond à l'état original x^i
  z=trajopt(matrice_trans,feedback,cout_instantane,cout_final,etat_initial)

  zz=list();
  zz(1)=etat(z(1));
  zz(2)=commande(z(2));
  zz(3)=z(3);

7 Problèmes de discrétisation

Les problèmes de discrétisation de l’espace d’état se posent assez naturellement quand le problème d’optimisation n’est pas formulé directement avec des matrices de transition, mais plutôt par le biais d’une fonction de dynamique. C’est ainsi qu’un problème d’optimisation déterministe est généralement formulé, mais aussi de nombreux problèmes d’optimisation stochastique.

7.1 Dynamique stochastique commandée

Une dynamique commandée est une application

F : 𝕏 × 𝕌 → 𝕏.

(25)

Une dynamique stochastique commandée est une application

F : 𝕏 × 𝕌 × Ω → 𝕏

(26)

où 𝕏 est un espace mesurable, (Ω,μ) est un espace de probabilité, et où F(x,u,⋅) est supposée mesurable pour tous (x,u) ∈ 𝕏 × 𝕌.

Dans le cas où 𝕏 = {x¹,...,xⁿ} est un ensemble ﬁni et où

F : {x1,...,xn } × 𝕌 × Ω → {x1,...,xn },

(27)

il est immédiat de déﬁnir une famille de matrices de transition sur 𝕏 = {x¹,...,xⁿ} par

M u(xi,xi′) = ℙ (x (t + 1) = xi′ | x (t) = xi,u(t) = u) = μ (ω ∈ Ω | F (xi,u,ω ) = xi′}.

(28)

7.2 Discrétisation d’un problème d’optimisation stochastique sur ℝ

Nous supposons donnée une dynamique stochastique commandée

F : ℝ × 𝕌 × Ω → ℝ.

(29)

Nous discrétisons l’espace d’état ℝ en {x¹,...,xⁿ}, avec x¹ < ⋅⋅⋅ < xⁿ. Nous discrétisons également l’espace de commande 𝕌 en {u¹,...,u^l}.

Nous déﬁnissons les applications prédécesseur Π : ℝ →{x¹,...,xⁿ} et successeur Σ : ℝ →{x¹,...,xⁿ} par

{ ′ 1 n ′ ∀x ∈ ℝ Π (x ) = sup {x ′∈ {x1,...,xn} | x ′≤ x } Σ (x ) = inf {x ∈ {x ,...,x } | x > x }

(30)

avec la convention que sup ∅ = x¹ et inf ∅ = xⁿ.

Nous allons distribuer aléatoirement tout x vers son prédécesseur et vers son successeur par la fonction ψ : ℝ × [0, 1] →{x¹,...,xⁿ} déﬁnie par

( || Π (x) = Σ (x ) ⇒ ψ (x,ω′) = Π (x) = Σ (x ) ||| |{ (| ′ x − Π(x) |{ Π (x) si ω > Σ-(x-) −-Π-(x) ||| Π (x) < Σ (x ) ⇒ ψ (x,ω′) = ||| ||( ′ --x −-Π(x)--- ( Σ (x) si ω ≤ Σ (x ) − Π (x).

(31)

On discrétise alors la dynamique stochastique F par

F^ : {x1,...,xn} × 𝕌 × (Ω × [0,1],μ ⊗ dω′) → {x1, ...,xn }

(32)

où

^F (xi,u,ω,ω ′) = ψ (F(xi,u,ω ),ω′).

(33)

7.3 Discrétisation d’un problème d’optimisation stochastique sur ℝⁿ

7.4 Discrétisation d’un problème d’optimisation déterministe sur ℝ et mise sous forme stochastique

Nous reprenons ce qui est fait ci-dessus, mais avec F déterministe.

La dynamique commandée dynamique_commandee est transformée en une famille de matrices de transition discretisation(etat,commande,dynamique_commandee) : partant d’un état xⁱ et appliquant la commande u^l, on transite vers les deux voisins encadrant dynamique_commandee(xⁱ,u^l) (ou un seul si on est en deçà ou au delà de {x¹,...,xⁿ}) avec une probabilité égale à la distance normalisée à ces points (ou égale à 1).

La macro Scilab suivantepredec_succes utilise la fonction dsearch, pas forcément disponible sur les versions les plus anciennes de Scilab. En cas de problème, la remplacer par le code à télécharger ici.

  function indices_image_discretisee=predec_succes(etat,image)
    //etat = vecteur ordonné
    //image = vecteur
    // indices_image_discretisee = liste
    //indices_image_discretisee(1)(i)=j
    // SSI etat(j) est le prédécesseur de image(i)
    //indices_image_discretisee(2)(i)= j
    // SSI etat(j) est le successeur de image(i)
    //indices_image_discretisee(3)(i)=
    //probabilité d'aller vers le prédécesseur de image(i)

    cardinal_etat=prod(size(etat));
    ind=dsearch(image,etat);
    // image(i) \in [etat(ind(i)), etat(ind(i)+1)[ sauf si
    // ind(i)=0, auquel cas image(i)<etat(1) ou image(i)>etat($)
    //ind_nd=ind(image_ind);
    image_ind_h=find(image > etat($));
    ind(image_ind_h)=cardinal_etat;
    // ind(i)=0 ssi image(i)<etat(1)
    // ind(i)=n ssi image(i)>etat($)
    image_ind_m=find(ind > 0 & ind < cardinal_etat);
    // indices du milieu : image(image_ind_m) \subset [etat(1),etat($)]
    // ind(image_ind_m) \subset ]0,n[
    image_ind_b=find(ind==0);
    // indices du bas : image(image_ind_b) <etat(1)
    image_ind_h=find(ind==cardinal_etat);
    // indices du haut : image(image_ind_h) >etat($)

    ind1=zeros(image);ind2=ind1;
    proba=zeros(image);

    ind1(image_ind_h)=ind(image_ind_h);
    // envoie les indices pour lesquels l'image de l'etat correspond est
    // >etat($) vers cardinal_etat, dernier indice de etat
    ind2(image_ind_h)=ind1(image_ind_h);
    // prédécesseur = successeur = etat($)
    proba(image_ind_h)=ones(image_ind_h);
    // probabilité 1 d'aller vers prédécesseur = successeur

    ind1(image_ind_m)=ind(image_ind_m);
    ind2(image_ind_m)=1+ind(image_ind_m);
    proba(image_ind_m)=(etat(1+ind(image_ind_m))-image(image_ind_m)) ./ ...
                       (etat(1+ind(image_ind_m))-etat(ind(image_ind_m)));

    ind1(image_ind_b)=ind(image_ind_b)+1;
    // +1 à cause de ind(image_ind_b) composé de 0
    // envoie les indices pour lesquels l'image de l'etat correspond est
    // <etat(1) vers 1, premier indice de etat
    ind2(image_ind_b)=ind1(image_ind_b);
    // prédécesseur = successeur = etat(1)
    proba(image_ind_b)=ones(image_ind_b);
    // probabilité 1 d'aller vers prédécesseur = successeur

    indices_image_discretisee=list();
    indices_image_discretisee(1)=ind1;
    indices_image_discretisee(2)=ind2;
    indices_image_discretisee(3)=proba;
  endfunction


  function z=egal(x,y)
    z=bool2s(x==y)
  endfunction


  function matrice_transition=discretisation(etat,commande,dynamique_commandee)
    matrice_transition=list()
    //etat = vecteur ordonné
    //commande = vecteur
    //dynamique_commandee = fonction(x,u)
    //matrice_transition = liste de matrices de transition
    //  sur les indices de etat
    cardinal_etat=prod(size(etat));
    for l=1:prod(size(commande)) do
      //la liste matrice_transition est indicée par les indices du vecteur commande
      image=dynamique_commandee(etat,commande(l));
      // vecteur des images du vecteur etat
      indices_image_discretisee=predec_succes(etat,image)
      indices1=indices_image_discretisee(1);
      indices2=indices_image_discretisee(2);
      probabilites=indices_image_discretisee(3);

      M1=zeros(cardinal_etat,cardinal_etat);
      M2=zeros(M1);
      for i=1:cardinal_etat do
        M1(i,indices1(i))=probabilites(i);
        M2(i,indices2(i))=1-probabilites(i);
      end
      matrice_transition(l)=M1+M2

      //coincidmoins=feval(zmoins(1,:),etat,egal)
      // coincidmoins(i,j)=1 ssi zmoins(1,i) == etat(j)
      // ssi partant de etat(i), etat(j) est visité
      //coincidplus=feval(zplus(1,:),etat,egal)
      //coincidmoins=1-sign ( abs ( log( ( exp(-zmoins(1,:)') ) * exp(etat) )));
      //coincidplus=1-sign ( abs ( log( ( exp(-zplus(1,:)') ) * exp(etat) ) ) );
      //matrice_transition(l)= coincidmoins .* ( ones(etat)' * zmoins(2,:) )' +...
      //coincidplus .* ( ones(etat)' * zplus(2,:) )'
      // affectation des probabilités
    end
  endfunction

7.5 Discrétisation d’un problème d’optimisation stochastique sur ℝ^N, N ≥ 2

On peut généraliser ce qui a été fait pour N = 1, en discrétisant ℝ^N sur une grille ﬁnie

𝕏 = {x1,...,xn1} × ⋅⋅⋅ × {x1 ,...,xnN } 1 1 N N

(34)

et en déﬁnissant des voisins 𝒱(x) ⊂ 𝕏 pour tout point x ∈ ℝ^N, avec des probabilités de transiter vers ces voisins.

On est alors ramené à une chaîne de Markov sur un espace ﬁni 𝕏, qui a la particularité d’être un espace produit. Il reste alors à numéroter les éléments de 𝕏 pour se ramener à un problème d’optimisation stochastique sur des entiers.

Bien sûr, en terme de programmation informatique, ceci n’est pas si simple...

Formations

Départements

Pratique

Recherche

Management of Energy

Management of Renewable Resources and Biodiversity

Management of Exhaustible Resources

Models in Life Science

Miscellaneous

Programmation dynamique
Macros générales

Contents

1 Programmation dynamique stochastique sur des entiers (cas des matrices de transition)