1 Programmation dynamique stochastique sur des entiers (cas des matrices de
transition)
Le problème d’optimisation considéré ici est
(1)
avec
coût instantané, L : {0,...,T − 1}×{1,...,n}×{1,...,p}→ ℝ, où
n ∈ ℕ∗ est le cardinal de l’espace d’état 𝕏 := {1,...,n}⊂ ℕ∗,
p ∈ ℕ∗ est le cardinal de l’espace de commande 𝕌 := {1,...,p}⊂ ℕ∗,
T ∈ ℕ∗ est l’horizon,
coût final, Φ : {1,...,n}→ ℝ,
matrices de transition, M : {1,...,p}→ℳ(n,n), application à valeurs dans les matrices
(n,n) stochastiques, c’est-à-dire ∀u = 1,...,p,
le terme Mu(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′ :
fonction valeur
(2)
la fonction valeur évaluée sur l’état x à l’instant t,
l’équation de la programmation dynamique
(3)
feedback optimal
(4)
si l’inf est atteint en un unique argument,
trajectoires d’état, de commande, de coût
(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 final
dynamique bruitée, F : {0,...,T − 1}×{1,...,n}×{1,...,p}× 𝕎 →{1,...,n}, application
telle que
(6)
où w(0),…,w(T − 1) est une suite de variables aléatoires i.i.d.
fonction valeur
(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
(8)
feedback optimal
(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) :
(10)
où on notera le décalage dans les indices temporels, et où
le cardinal n de l’espace d’état est un entier cardinal_etat
(11)
le cardinal p de l’espace de commande est un entier
(12)
l’horizon T est un entier horizon
(13)
le coût final est un vecteur cout_final de dimension (n, 1) :
(14)
Transitions
les matrices de transition sont une liste mat_transition à p éléments composée de
matrices (n,n) stochastiques :
(15)
la dynamique bruitée est une fonction dynamique(t,i,l,w)
la fonction valeur est une matrice valeur de dimension (n,T + 1) :
(16)
où on notera le décalage dans les indices temporels,
le feedback optimal est une matrice feedback de dimension (n,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) :
(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:2do //é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_commandedo 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-1do // 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 fini
Le problème d’optimisation est toujours (1), mais avec
coût instantané, L : {x1,...,xn}×{u1,...,up}×{0,...,T}→ ℝ, où
𝕏 := {x1,...,xn} est l’espace d’état,
𝕌 := {u1,...,up} est l’espace de commande,
T ∈ ℕ∗ est l’horizon,
coût final, Φ : {x1,...,xn}→ ℝ,
matrices de transition
(19)
On se ramène au cas des entiers en introduisant
nouveau coût instantané, : {1,...,n}×{1,...,p}×{0,...,T − 1}→ ℝ, avec
(20)
nouveau coût final, : {1,...,n}→ ℝ, avec
(21)
transition
nouvelles matrices de transition définies par
(22)
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 suffit alors d’introduire un vecteur etat de
dimension (1,n) et tel que
(23)
et un vecteur commande de dimension (1,p) et tel que
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
(25)
Une dynamique stochastique commandée est une application
(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ù 𝕏 = {x1,...,xn} est un ensemble fini et où
(27)
il est immédiat de définir une famille de matrices de transition sur 𝕏 = {x1,...,xn}
par
(28)
7.2 Discrétisation d’un problème d’optimisation stochastique sur ℝ
Nous supposons donnée une dynamique stochastique commandée
(29)
Nous discrétisons l’espace d’état ℝ en {x1,...,xn}, avec x1<< xn. Nous discrétisons
également l’espace de commande 𝕌 en {u1,...,ul}.
Nous définissons les applications prédécesseur Π : ℝ →{x1,...,xn} et successeur
Σ : ℝ →{x1,...,xn} par
(30)
avec la convention que sup ∅ = x1 et inf ∅ = xn.
Nous allons distribuer aléatoirement tout x vers son prédécesseur et vers son successeur par la
fonction ψ : ℝ × [0, 1] →{x1,...,xn} définie par
(31)
On discrétise alors la dynamique stochastique F par
(32)
où
(33)
7.3 Discrétisation d’un problème d’optimisation stochastique sur ℝn
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 xi et appliquant la commande ul, on transite vers les deux voisins
encadrant dynamique_commandee(xi,ul) (ou un seul si on est en deçà ou au delà de
{x1,...,xn}) 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_etatdo 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
finie
(34)
et en définissant 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 fini 𝕏, 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...