Le vendeur de journaux sur un horizon T

Jean-Philippe Chancelier
(last modification date: October 10, 2017)
Version pdf de ce document
Version sans bandeaux

1 Le problème du vendeur de journaux à un pas de temps

Contents

1 Le problème du vendeur de journaux à un pas de temps
 1.1 On cherche maintenant à résoudre le problème du vendeur de journaux sur un horizon T

1 Le problème du vendeur de journaux à un pas de temps

On reprend ici le code du T.P précédent et on va utiliser la loi de demande discrète sur trois points que l’on avait considéré sur le problème à un pas de temps.

// -*- Mode:scilab -*-  
 
// une loi discrète  
wi=[30,50,80];  
pi=[1/2,1/4,1/4];  
titre = sprintf("loi discrète sur %d valeurs",size(wi,'*'));  
mu=pi*wi'; // la moyenne  
// un echantillong de taille N  
N=100000;  
P=ones(size(wi,'*'),1)*pi;  
Y=grand(N,'markov',P,1);  
W=wi(Y);  
 
if %f then  
  function graphique_loi(titre,valeurs,probas)  
  // dessin de la loi  
    xset('window',1);  
    xbasc();  
    plot2d3(wi,pi,rect=[0,0,max(valeurs)+10,min(max(probas)+0.2,1)],style=2);  
    xtitle(titre);  
    // dessin de la fonction de repartition  
    xset('window',2);  
    xbasc();  
    plot2d2([0,wi],[0,cumsum(pi)],rect=[0,0,max(valeurs)+10,1.1],style=2)  
    plot2d([0,wi],[0,cumsum(pi)],rect=[0,0,max(valeurs)+10,1.1],style=-3)  
    xtitle(sprintf('Fonction de repartition de la %s',titre));  
  endfunction  
 
  graphique_loi(titre,wi,pi);  
  xclick();  
end  
 
// paramètres de la fonction cout  
 
c=10,cm=20,cs=5;cf=200;  
 
// la fonction coût j(u,w)  
 
function y=j(u,w)  
  y= c*u + cs*max(u - w,0) + cm*max(w - u,0)  
endfunction  
 
// La fonction J(u) pour la loi binomiale  
 
function y=J(u)  
  y= c*u + cs*pi*max((u- wi'),0)+ cm*pi*max(( wi'-u),0);  
endfunction  
 
// graphique de la fonction J(u) et recherche du minimum  
 
if %t then  
  U=-10:100;  
  Ju=[];  
  for u=U, Ju=[Ju,J(u)];end  
  xset('window',1);  
  xbasc();  
  plot2d(U,Ju,style=2);  
  xtitle("La fonction J");  
  // recherche du minimum  
  [m,kmin]=min(Ju);  
  uopt=U(kmin);  
  S=uopt;  
  printf("Nombre optimal de journaux a commander %f\n",U(kmin));  
  printf("Moyenne de la demande %f\n",mu);  
  plot2d(U(kmin),m,style=-4);// dessin  
  plot2d3(U(kmin),m);  
  xclick();  
end  
 
// On regarde maintenant un cas ou le vendeur à déjà des journaux  
// et ou il paye un coup fixe si il commande des journaux  
// on voit une stratégie [s,S]  
 
function y=Jtilde(u,x)  
  y= cf*(u>0) + J(u+x) -c*x  
endfunction  
 
// On calcule pour x variant de 0 à 2*max(wi)  
// la commande optimale de journaux correspondante  
 
xuopt=[];  
xv=0:1:2*max(wi);  
for x0=xv  
  U=0:2*max(wi);  
  Ju=[];  for u=U, Ju=[Ju,Jtilde(u,x0)];end  
  g=[];  
  [m,kmin]=min(Ju);  
  xuopt=[xuopt,U(kmin)];  
end  
 
I=find(xuopt==0);  
s=xv(I(1)-1);  
 
// vérifier que la stratégie optimale est [s,S];  
 
xset('window',1);  
xbasc();  
plot2d2(xv,xuopt,style=2);  
plot2d2(xv(1:s),xv(1:s)+xuopt(1:s),style=1);  
xtitle(sprintf("Valeur de s=%d et de S=%d",xv(s),xv(s-1)+xuopt(s-1)));  
 
//  
 
// trouver s, S=uopt (calculé avant)  
x=0:2*max(wi);  
Jv=[];  
for i=1:size(xv,'*');  Jv=[Jv,J(x(i))];end  
costs= Jv - (cf+ J(uopt));  
I=find(costs <= 0);  
if isempty(I) then  
  s=0;  
else  
  s=x(I(1)-1)  
end  
 
printf("On trouve s=%d et S=%d\n",s,S);  
 
 

Question 1 Faites éxecuter le code précédent dans scicoslab et vérifiez que la stratégie optimale est bien de la forme [s,S].

1.1 On cherche maintenant à résoudre le problème du vendeur de journaux sur un horizon T

On propose ici un squelette de programme qui implémente l’équation de Programmation dynamique. Il faut compléter ce programme et résoudre le problème sur un horizon T. À l’instant t = 1 le vendeur de journaux fait sa première commande et à l’instant T il ne commande plus. Le cas T = 2 correspond donc au problème du vendeur sur une journée.

// -*- Mode:scilab -*-  
// Vendeur de journaux à T-pas de temps  
 
exec('vendeur-T-un.sci');  
 
alpha=1.0;  
 
function y=ct(u,x)  
// coût instantané pour t autre que la date de départ  
  y= cf*(u>0) + c*u + alpha*cs*max(x,0) + alpha*cm*max(-x,0)  
endfunction  
 
function y=c0(u)  
// coût instantané pour t=1  
  y= cf*(u>0) + c*u  
endfunction  
 
function y=K(x)  
// coût final  
  y= cs*max(x,0) + cm*max(-x,0)  
endfunction  
 
// on transforme la dynamique pour rester dans un domaine  
// borné [-100,100]  
 
xmin=- 100; xmax = 100;  
 
function y=f(x,u,w)  
  y = min(max(x+u-w,xmin),xmax)  
endfunction  
 
// VT  
// On commencera avec T=2 qui doit donner le problème à un pas de temps  
 
 
T=2; // l'Horizon  
 
etats= xmin:xmax; // les etats possibles dans un domaine borné  
allu=  0:xmax;    // les commandes de journaux possible  
 
// La fonction v_T(x) = K(x)  
VT = alpha^T*K(etats);  
 
// On va calculer V_t pour T-1,...,1  
// et stocker les résultats dans une liste  
// On stocke aussi la commande optimale  
 
L=list(VT); // liste qui permet de stocker (V_T,V_{T-1},.....,V_1)  
U=list(0);  // liste qui permet de stocker (-- ,U_{T-1},.....,U_1)  
 
// boucle rétrograde en temps  
for t=T-1:-1:1  
 
  printf("t=%d\n",t);  
  Vt = zeros(1,size(etats,'*'));  
  Ut = %nan*ones(1,size(etats,'*'));  
  Vtp1= L($);  
  // On cherche a calculer ici Vt et Ut  
  // Vtp1 est la fonction de Bellman au temps t+1  
  for i=1:size(etats,'*')  
    x = etats(i);  
    mcost= %inf;  
    // boucle sur les commandes possibles  
    for k = 1:size(allu,'*')  
      u = allu(k);  
      // on calcule dans cost le cout associé a la commande u  
      // initialisation:  
      cost = 0 ;  
      // faire une boucle sur les aléas pour calculer  
      // E[ ct(u,x) + V_{t+1}(f(x,u,W))]  
      for ...  
 // xtp1: utiliser la dynamique de la chaine  
 xtp1 = ...  
 ix = find(xtp1 == etats) // donne l'indice de xtp1 dans le  
                                  // vecteurs des etats  
         // actualiser le cout  
 cost = cost + ...  
      end  
      // Il faut maintenant comparer  
      // cost au meilleur trouvé jusque la (on minimise)  
      // mcost et mettre à jour mcost  
      // et kumin (l'indice de la commande qui a donné le meilleur coût)  
      if ...  
 ...  
      end  
    end  
    // On stocke pour l'état x d'indice i le coût optimal  
    Vt(i) = mcost;  
    // et la commande optimale  
    Ut(i) = allu(kumin);  
  end  
  // On range Vt et Ut dans la liste  
  L($+1)=Vt;  
  U($+1)=Ut;  
end  
 
// dessin de la fonction de Bellman au temps t=1;  
xset('window',1);  
xbasc();  
plot2d2(etats,L($))  
 
// dessin de la commande optimale au temps t=1;  
xset('window',1);  
xbasc();  
plot2d2(etats,U($))  
 

Question 2