Gestion optimale d’une ressource renouvelable
Les questions

Luc Doyen
(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 Présentation du problème
2 Cas déterministe: croissance certaine
3 Cas aléatoire: croissance incertaine

1 Présentation du problème
2 Cas déterministe: croissance certaine
2.1 Résolution analytique par programmation dynamique
2.2 Programmation sous Scilab
2.3 Résolution par dynoptim, routine d’optimisation dynamique
3 Cas aléatoire: croissance incertaine
3.1 Résolution analytique par programmation dynamique
3.2 Programmation sous Scilab

1 Présentation du problème

On considère un planiﬁcateur souhaitant utiliser de manière optimale une ressource naturelle renouvelable (forêt, poissons...) sur T périodes. Le niveau de la ressource à l’instant t est noté P(t) et son prélèvement sur la période [t,t + 1[ est noté h(t). On suppose que la ressource non extraite à l’instant t se renouvelle à un taux de croissance R(t), ce qui donne la dynamique suivante :

P (t + 1) = R (t)(P (t) − h (t)) pour t = 0,..,T − 1.

(1)

Le bien-être que le planiﬁcateur souhaite optimiser est représenté par la somme actualisée des utilités de ses prélèvements successifs h(t) (supposés liés à la consommation par exemple) et du niveau de la ressource ﬁnale P(t) :

T∑− 1 ρtL(h(t)) + ρTL (P (T )), t=0

où ρ ∈ [0, 1] est un facteur d’actualisation et L(.) une fonction d’utilité.

2 Cas déterministe: croissance certaine

On suppose tout d’abord que la ressource naturelle est constituée d’une seule ressource de croissance certaine et constante R:

R(t) = R pour t = 0,..,T − 1.

(2)

Le planiﬁcateur souhaite optimiser le bien-être soit

{ } T∑−1 V(0,P0 ) = max ρtL (h (t)) + ρTL (P (T )) . h(0),h(1),..,h(T −1) t=0

(3)

On suppose maintenant que l’utilité est de type exponentiel

L(h ) = 1 − e− h.

(4)

2.1 Résolution analytique par programmation dynamique

Question 1

Que vaut la fonction valeur V (T,P) à l’instant ﬁnal ?
Quelle est la relation de récurrence liant les fonctions valeurs V (t − 1,.) et V (t,.)?
Montrer que, pour certaines valeurs de t et de P, le prélèvement optimal en boucle fermée h^♯(t,P) et la fonction valeur V (t,P) satisfont
$( | h♯(t,P ) = b(t)P + f(t), { ( −h♯(t,P )) | V(t,P ) = ρt a(t) − e------- , ( b(t)$ (5)

où les paramètres a(t),b(t),f(t) sont déﬁnis par les récurrences

$( || --Rb-(t)-- ||| b(t − 1) = 1 + Rb (t),b(T) = 1, |{ a(t − 1 ) = 1 + ρa (t),a(T ) = 1, ||| |||( f(t − 1) = f(t) −-log-(ρR-),f(T ) = 0. 1 + Rb (t)$ (6)

En déduire que les prélèvements optimaux en boucle fermée h(t) := h^♯(t,P^♯(t)) où P^♯(t + 1) = R(P^♯(t) − h^♯(t,P^♯(t))) satisfont la relation :

$-- --♯ h(t + 1) − h (t) = log(ρR ).$ (7)

Donner alors des conditions sur ρ et R pour que les prélèvements optimaux soient croissants dans le temps. Que se passe-t-il en particulier si ρR = 1 ?

2.2 Programmation sous Scilab

Question 2 Ouvrir un ﬁchier nom_de_fichier.sce (par exemple, TP_renouv.sce) et y recopier les paramètres suivants.

Y construire les vecteurs b et f donnés par l’équation (6), puis créer les vecteurs de stocks et de prélèvements optimaux Popt et hopt en utilisant la dynamique (1) et l’équation (5) du feedback.

Tracer trajectoire et décisions optimales en fonction du temps, à l’aide de la commande Scilab plot2d2. On notera que le vecteur hopt est de taille (1,T), alors que Popt est de taille (1,T + 1) ; on rajoutera donc la dernière valeur Popt($) de Popt à hopt pour signiﬁer que le dernier prélèvement consiste à épuiser totalement la ressource.

  //paramètres

  Horizon=4;
  R=1.1;
  P0=10;
  rho=0.8;
  //rho=0.1

  // Construction de b et f

  b=[];
  b(Horizon+1)=1;
  for t=Horizon:-1:1 do
    b(t)=R*b(t+1)/(1+R*b(t+1));
  end;

  f=[];
  f(Horizon+1)=0;
  for t=Horizon:-1:1 do
    f(t)=(f(t+1)-log(rho*R))/(1+R*b(t+1));
  end;


  // Prélèvements  et ressource optimales

  Popt=[];
  hopt=[];
  Popt(1)=P0;
  for t=1:Horizon do
    hopt(t)=b(t)*Popt(t)+f(t);
    Popt(t+1)=R*(Popt(t)-hopt(t));
  end,


  // Affichage
  xbasc();
  plot2d2([0:(Horizon+1);0:(Horizon+1)]',[[hopt;Popt($);0]';[Popt;0]']', ...
          rect = [0,0,Horizon+2,Popt(1)+1])

Question 3 Répéter l’opération précédente pour diﬀérentes valeurs du taux d’actualisation ρ entre 0.5 et 1. Quel résultat analytique retrouve-t-on ?

2.3 Résolution par dynoptim, routine d’optimisation dynamique

On introduit, pour des facilités de programmation, une nouvelle commande u ∈ [0, 1] telle que h = uP. On considère alors le problème

( || T∑−1 |{ max ρtL (h(t)) + ρT L(P (T)) t=0 ||| P (t + 1) = R(t)(1 − u(t))P (t) pour t = 0, ..,T − 1 ( 0 ≤ u(t) ≤ 1 pour t = 0,..,T − 1.

(8)

comme un problème d’optimisation en la variable (P(0),...,P(T),u(0),...,u(T − 1)) ∈ ℝ^T+1 × [0, 1]^T sous les contraintes d’égalité P(t + 1) = R(t)(1 − u(t))P(t) et d’inégalité 0 ≤ u(t) ≤ 1.

La macro dynoptim permet de traiter ce type de problèmes.

Question 4 Charger dynoptim (menu toolboxes sous Scilab) et consulter le help. Déﬁnir la dynamique dyn et ses dérivées partielles dyn_P par rapport à l’état et dyn_u par rapport à la commande. Faire de même avec le coût instantané et le coût ﬁnal. Écrire les contraintes sur les commandes.

  help dynoptim

  ////////////////////////
  //dynamique
  ////////////////////////

  function z=dyn(u,P,t) z=R*P*(1-u),endfunction;
  function D=dyn_P(u,P,t) D=R*(1-u),endfunction;
  function D=dyn_u(u,P,t) D=-R*P,endfunction;

  ////////////////////////
  // Fonction d'utilité
  ////////////////////////

  function z=Utilite(h) z=1-exp(-h),endfunction;
  function D=DUtilite(h) D=exp(-h),endfunction;

  ////////////////////////
  // cout instantané (à minimiser !)
  ////////////////////////

  function z=L(u,P,t) z=-(rho^t)*Utilite(u*P),endfunction;
  function D=L_P(u,P,t) D=-(rho^t)*DUtilite(u*P)*u,endfunction;
  function D=L_u(u,P,t) D=-(rho^t)*DUtilite(u*P)*P,endfunction;

  ////////////////////////
  // cout final (à minimiser !)
  ////////////////////////

  function z=g(P,t) z=-(rho^t)*Utilite(P),endfunction;
  function D=g_P(P,t) D=-(rho^t)*DUtilite(P),endfunction;

  ////////////////////////
  // contrainte controles
  ////////////////////////

  u_min=zeros(1,Horizon);
  u_max=ones(1,Horizon);

Question 5 Choisir un niveau initial de ressource S0 et un horizon T. Initialiser l’algorithme d’optimisation dynamique dynoptim avec un vecteur u_init. Faire un appel à dynoptim et tracer les solutions. Vériﬁer dans quels cas les solutions coïncident avec celles des questions précédentes.

  Horizon=4;
  //Horizon=40;

  //conditions initiales
  P0=[10];
  u_init=0.5*ones(1,Horizon);

  ////////////////////
  // Appel dynoptimsc
  //////////////////////
  L_noms=["L","L_u","L_P"];
  dyn_noms=["dyn","dyn_u","dyn_P"];
  g_noms=["g","g_P"];

  [u_opt,P_opt]=dynoptim(L_noms,dyn_noms,g_noms,u_min,u_init,u_max,P0),

  h_opt=u_opt .*P_opt(1:($-1));

  // Affichage
  xbasc();
  plot2d2([0:(Horizon+1);0:(Horizon+1)]',[[h_opt';P_opt($);0]';[P_opt';0]']', ...
          rect = [0,0,Horizon+2,P_opt(1)+1])

Question 6 Que constatez-vous en faisant varier le facteur d’actualisation ρ ?

3 Cas aléatoire: croissance incertaine

On suppose maintenant que la ressource totale P(t) est constituée de deux ressources P₁(t) et P₂(t). Le taux de croissance de la première ressource, supposé aléatoire, est notée M(t), tandis que la seconde ressource a une croissance certaine, de taux noté N. On appelle α(t) la proportion de la ressource P₂ dans le total des ressources P(t) à la période t. Ainsi, le taux de croissance total, aussi aléatoire, s’écrit:

R (t) = α (t)N + (1 − α(t))M (t).

(9)

Pour simpliﬁer, on suppose que les variables aléatoires (M(t))_{t=0,...,T−1} sont indépendantes et de même loi.

On considère aussi le niveau de prélèvement agrégé h(t) comme décision, ce qui donne la dynamique

P (t + 1) = (α(t)N + (1 − α (t))M (t))(P (t) − h(t)).

(10)

On suppose que le planiﬁcateur optimise l’espérance de la somme actualisée des utilités de ses prélèvements successifs et de la ressource ﬁnale

{T −1 } ∑ t V (0,P0) = ({ max 𝔼 ρ L(h (t)) h (0),h(1),..,h(T − 1) t=0 ( α (0),α(1),..,α (T − 1)

(11)

où h(t) et α(t) correspondent à une stratégie de prélèvement et de répartition.

On se restreint ici au cas dit isoélastique correspondant à une utilité

L (h) = hγ avec 0 < γ < 1.

(12)

3.1 Résolution analytique par programmation dynamique

Question 7

Montrer que l’équation de Bellman s’écrit, pour t = 0,...,T − 1,
$( ) V (P,t) = max ρtL(h ) + 𝔼 (V(t + 1,(αN + (1 − α)M )(P − h)) , α∈[0,1],h≥0$ (13)

où M est une variable aléatoire de même loi que la loi commune des (M(t)).
Montrer, par récurrence, que le prélèvement optimal h^♯(t,P), la proportion optimale α^♯(t) et la fonction valeur V (t,P) satisfont
$({ t γ−1 γ V(t,P ) = ρ b(t) P , ( ♯ , h (t,P ) = b(t)P , t = 0,...,T − 1,$ (14)

et

$( ♯ ♯ )γ−1 𝔼((N − M ) α (t)N + (1 − α (t))M ) = 0.$ (15)

Montrer que l’équation de récurrence satisfaite par b(t) est

$b(t) = --ab-(t-+-1)--, b(T ) = 1 1 + ab(t + 1)$ (16)

où on utilise , l’équivalent certain de croissance optimale, c’est-à-dire

$L(^R ) = 𝔼[L(R ♯)] avec R ♯ = α♯N + (1 − α♯)M$ (17)

et le terme

$a = (ρ^Rγ) 1γ−1.$ (18)

3.2 Programmation sous Scilab

On suppose que

( || P0 = 10, ||{ T = 10, γ = 0.5, || ||( N = 1.1, ρ = 1∕N,

Question 8 Ouvrir un ﬁchier nom_de_fichier.sce et y recopier les paramètres précédents. Écrire une fonction scilab util(x) pour la fonction d’utilité puissance.

  //paramètres
  puis=0.5;
  Horizon=10;
  N=1.1;rho=1/N;
  P0=10;

  // Fonction d'utilité puissance
  function u=util(x) u=x^puis,endfunction;

On suppose aussi que le taux de croissance risqué M suit une loi du type

{ M = M ♭ avec une probabilit´e p M ♯ avec une probabilit´e q = 1 − p

(19)

Question 9 Montrer que la solution de (15) est donnée par

♯ -----−-p2(N--−-M-♭)2M-♯ +-q2(N-−-M-♯)2M-♭------ α = p2(N − M )2(N − M ) − q2(N − M )2(N − M ). ♭ ♯ ♯ ♭

(20)

Prendre les paramètres suivants pour p,q,M_♭,M_♯. Calculer la proportion optimale α^♯ de la ressource 1.

  // Données variable aléatoire

  p=0.35;q=1-p;
  Mb=1.3;M#=1;

  //Construction de  proportion optimale

  alpopt=(-p^2*(N-Mb)^2*M#+q^2*(N-M#)^2*Mb)/.(p^2*(N-Mb)^2*(N-M#)-q^2*(N-M#)^2*(N-Mb));

Question 10 Construire le vecteur b donné par l’équation (14) en utilisant l’équivalent certain déﬁni en (17) et le paramètre a déﬁni en (18).

  //Construction de Rstar, Rchap et a

  function R=R_OPT(M) R=alpopt*N+(1-alpopt)*M,endfunction;

  Rchap=(p*util(R_OPT(Mb))+q*util(R_OPT(M#)))^(1/puis);

  a=(rho*util(Rchap))^(1/(puis-1));

  //Construction de b

  b=[];
  b(Horizon+1)=1;
  for t=Horizon:-1:1 do
    b(t)=a*b(t+1)/(1+a*b(t+1));
  end;

Question 11 Écrire une fonction Scilab qui, à un niveau global de ressource ressource P et à une valeur M de taux de croissance, associe le niveau total de la ressource au pas de temps suivant résultant de l’application de la commande optimale (h^♯(t,P),α^♯).

  // Prélèvement  et ressource optimaux

  function h=PREL_OPT(P,t) h=b(t)*P,endfunction;
  function x=DYN_OPT(P,M,t) x=R_OPT(M)*(P-PREL_OPT(P,t)),endfunction;

Question 12 Écrire une procédure de simulation de la variable aléatoire M.

  // Simuler variable aléatoire M
  z=rand(1,Horizon,'uniform');
  for t=1:Horizon do
    if (z(t) <= p) then
      Msimu(t)=Mb;
    else
      Msimu(t)=M#;
    end;
  end;

Question 13 À partir de P₀, générer une trajectoire optimale stochastique P^♯(t) de ressources. Tracer les trajectoires t(P^♯(t),h^♯(t)) pour une réalisation de M.

  Popt(1)=P0;J=0;
  for t=1:Horizon do
    hopt(t)=PREL_OPT(Popt(t),t);
    Popt(t+1)=DYN_OPT(Popt(t),Msimu(t),t);
    J=J+rho^(t-1)*util(hopt(t));
  end,
  J=J+rho^(Horizon)*util(Popt(Horizon+1));

  xbasc();plot(hopt);plot(Popt);

Question 14 Écrire une procédure Scilab fournissant S réalisations des variables aléatoires J(h^♯(.)), P^♯(t) et h^♯(t). En déduire des valeurs approchées pour leur moyenne et variance. Donner un histogramme de J(h^♯(.)).

Formations

Départements

Pratique

Recherche

Management of Energy

Management of Renewable Resources and Biodiversity

Management of Exhaustible Resources

Models in Life Science

Miscellaneous

Gestion optimale d’une ressource renouvelable
Les questions

Contents

1 Présentation du problème

2 Cas déterministe: croissance certaine

2.1 Résolution analytique par programmation dynamique

2.2 Programmation sous Scilab

2.3 Résolution par dynoptim, routine d’optimisation dynamique

3 Cas aléatoire: croissance incertaine

3.1 Résolution analytique par programmation dynamique

3.2 Programmation sous Scilab

L'École des Ponts ParisTech est membre fondateur de	L'École des Ponts ParisTech est certifiée

Formations

Départements

Pratique

Recherche

Management of Energy

Management of Renewable Resources and Biodiversity

Management of Exhaustible Resources

Models in Life Science

Miscellaneous

Gestion optimale d’une ressource renouvelable Les questions

Contents

1 Présentation du problème

2 Cas déterministe: croissance certaine

2.1 Résolution analytique par programmation dynamique

2.2 Programmation sous Scilab

2.3 Résolution par dynoptim, routine d’optimisation dynamique

3 Cas aléatoire: croissance incertaine

3.1 Résolution analytique par programmation dynamique

3.2 Programmation sous Scilab

Gestion optimale d’une ressource renouvelable
Les questions