Fermer X


1 Optimal Allocation in Random Environment (Theory)

1.1 Optimisation model in random environment

We extend the model in constant environment introduced at Optimal Resources Allocation in Ecology: to Grow or to Reproduce? (I) to a random environment as follows. We suppose that

  • the total biomass generated during one period is now a function of environmental state R; thus, x`→f(x) is now replaced by a function x`→f(x,R) with the same properties;
  • the environmental state Rt occurs during period [t,t+1[; we assume that it is a random variable having values in a finite set {r1,...,rm};
  • the random variables sequence R0,...,RT1 forms a homogeneous Markov chain; we denote by π the matrix of general term π(i,j) = (Rt+1 = rj|R t = ri);
  • at the beginning of period [t,t + 1[, the realization of the random variable Rt is experienced by the plant.

Stochastic dynamics

The stochastic dynamics of the plant is the following. Consider a plant starting at the beginning of period [t,t + 1[ with vegetative biomass x + and subject to environmental state ri. Applying the control v [0, 1], at the end of period [t,t + 1[, the plant

  • either dies, i.e. transits towards vegetative biomass 0, with probability 1 β ;
  • transits towards vegetative biomass vf(x,rj) with probability βπ(i,j).

In parallel, the environmental state ri transits towards state rj with probability π(i,j).


The control v is here the fraction of total biomass allocated to growth.

Now, a (Markovian) strategy v() is a mapping {0,...,T 1+ ×{r1,...,rm}→ [0, 1]: v0,...,vT1 are in feedback on both vegetative biomass and on environmental state.

Feedback strategies

An optimal strategy for the plant is a strategy which maximizes the mathematical expectation of the fitness, here the cumulated offspring (reproductive biomass) at the end of the season:

J(v(.)) = 𝔼(   (1 − vt)f(xt,Rt))  with    v(.) = (v0,...,vT−1) ∈ [0,1]T .

Since the environmental states have values in the finite set {r1,...,rm}, and since mortality may be modeled by a random variable with values in {0, 1}, the expectation is here a finite sum.

1.2 Resolution by stochastic dynamic programming

The stochastic dynamic programming equation (sdpe) is

V (x,r,t) = max  𝔼{(1 − v)f(x, R   )+ (1 − β)V (0,t+ 1) + βV (vf(x, R   ),t+ 1) | R =  r}.
           0≤v≤1                t+1                                  t+1           t

Question 1 Show by induction that V (0,t) = 0, so that

      i          ∑                      j               j
V(x, r,t) = m0a≤xv≤1    π(i,j){(1 − v)f(x,r ) + βV (vf(x, r ),t + 1)}.

2 Computation of optimal strategies in random environment (PW Scilab)

We take

f(x,r) = rxγ   with  0 < γ ≤  1.

The reader has less informations than in the constant environment case. However, the Scilab code includes many comments.

2.1 Discretization of the problem

In this model, contrarily to the constant environment case, the state is no longer the single vegetative biomass x, but is the couple (x,R) + ×{r1,...,rm}.

First, we discretize the space + in {x1,...,xn}, where x1 corresponds to x = 0 (dead plant) and xn is a maximal vegetative biomass for the plant. We also discretize the function (x,r,v)`→vf(x,r) in a function

F : {x1, ...,xn } × {r1,...,rm} × [0,1] → {x1, ...,xn }.

Then, we express the stochastic dynamics above by a family, indexed by the controls v, of transition probabilities Mv on the finite product set {x1,...,xn}×{r1,...,rm}.

2.2 Using hypermatrices

Since the state is two dimensional, after discretization it becomes a couple in a finite product set. Thus, transition probabilities are not naturally matrices, but are rather hypermatrices.

The following Scilab code expresses these transition probabilities Mv by means of four indexes hypermatrices. An hypermatrix is an extension with more than two indexes of Scilabmatrices.

The vector taille has for elements the sequence (x1,...,xm) of values taken by discretized biomass: taille(i) = xi, i = 1,...,m.

2.2.1 Environmental states transitions

The vector env has for elements the sequence (r1,...,rm) of values taken by the environment: env(i) = ri, i = 1,...,m.

Open a new file plantII.sci and copy the following code. It provides different examples of transition matrices π for the environment Markov chain.

  function pi=tr_env_auc_cor_unif(cardinal_alea)
    // no correlation + uniformity
    // environment = i.i.d. with uniform common law
  function pi=tr_env_cor_proche_vois(cardinal_alea)
    // correlation with closest neighbours
    pi=diag([0.5,1/3*ones(1:(dd-2)),0.5])+diag([0.5,1/3*ones(1:(dd-2))],1)+ ...
  function pi=tr_env_cor_crois(cardinal_alea,rho)
    // correlation with increasing trend
    // rho is the correlation intensity 0.8
  function pi=tr_env_cor_decrois(cardinal_alea,rho)
    // correlation with decreasing trend
    // rho is the correlation intensity 0.8
  function pi=tr_env_cor_et_chute(cardinal_alea,rho)
    // correlation with probability of falling
    // rho is the correlation intensity 0.9
  function pi=tr_env_cor_et_ascension(cardinal_alea,rho)
    // correlation with a probability of rise
    // rho is the correlation intensity 0.9

2.2.2 State transitions

In the file plantII.sci, copy the following code which encapsulates both state and environment transitions in the so called Scilab object hypermatrices.

  function Hypermatrices=constr_HyperM(taille,env,control,pi,controlled_dynamics)
    // Construction of a family of matrices of transition probabilities on
    //  \{1,...,n\} X \{1,...,m\} from the problem data:
    // Hypermatrices is a list of hypermatrices indexed by the control
    // control is a vector
    // controlled_dynamics(x,r,u) is a function with real values
    for l=1:cardinal_control do
      // allows a natural expression for transitions from a product space
      // towards itself
      for y=1:dims(2) do
        // vector of the images of the vector state for r fixed
        for i=1:cardinal_taille do
        for z=1:dims(2) do
          // Alive (index>1), the plant faces a death probability
          // and the environment transitions are those of the matrix pi
      // and not Hypermat(x,dyn,y,:)=(1-ateb)*pi(y,:) end, end
      // because, at the previous step, one may have attributed
      // Hypermat(x,dyn,y,:) with Hypermat(x,1,y,:)=ateb*pi.

2.3 Conversion operations to fit with the Scilab function Bell_stoch

The Scilab function Bell_stoch allows numerical resolution of the stochastic dynamic programming equation. It is written for a Markov chain on {1,2...,taille_state}. Indeed, the algorithm is powerful because state is identified with a vector index.

This is why, the four indexes transition hypermatrices will be transformed into ordinary matrices with two indexes by “unrolling” the indexes. For instance, the following matrix (particular case of hypermatrix) (  1  3 )
(       )
   7  5
   2  8 will be unrolled in (1, 7, 2, 3, 5, 8).

This is how to a couple (x,r) ∈{x1,...,xn}×{r1,...,rm} we associate a state e ∈{1,...,n×m}. We numerically solve the optimization problem with this state space (intermediary and one-dimensional). Then, we go back to the original two-dimensional state space.

In the file plantII.sci, copy the following code in which conversion Scilab macros are defined.

  function Matrices=conv_HyperM(Hypermatrices)
    //Matrices is a list of transition matrices indexed by
    //the control
    //dimensions of the hypermatrix
    for i=1:prod(size(control)) do
      for j=1:dims(2) do
        if dims(2) <> 1 then
        // Transformation into transition matrix (prod(dims),prod(dims))
        // This latter matrix M is adapted to the resolution of the SDPE
  function n=convert1(i,j,dims)
    // i and j are coordinates in a hypermatrix of dimension dims
    // n is the position of (i,j) in a vector of dimension prod(dims)
    if i > dims(1) | j > dims(2) then
      n='i ou j &gt; dims'
  function [i,j]=convert2(n,dims)
    // n is the position of (i,j) in a vector of dimension prod(dims)
    // i and j are scoordinates in a hypermatrix of dimension dims
  function Hmat=convert2_mat(matrix,dims)
    // matrix is a matrix (ex:feed) whose elements come from
    // a vector space with dimension greater than 2
    // Hmat is the hypermatrix corresponding to this matrix (ex:feed_hyp)
    for t=1:T-1 do
      for index=1:prod(dims) do
  function instant_cost=conv_cost(taille,env,control,cost,horizon)
    //instant_cost is a list of instantaneous costs indexed for
    // each value of the control
    //taille is a vector
    //env is a vector
    //control is a vector
    //cost is a function(x,u,t)
    for j=1:prod(size(control)) do
      for t=1:horizon do
        // cost_it = hypermat(dims,sparse([],[],[dims(1),dims(2)]));
        for state1=1:dims(1) do
          for state2=1:dims(2) do
        // Transformation into vector : size(cost_i)=[1,prod(dims)]
        // adapted to the resolution of the SDPE

2.4 Transitions and costs adapted to the Scilab function Bell_stoch

Open a file plantII.sce and copy the following code.

  // env=[10:12]
  // control=0:0.5:1;
  pi=tr_env_auc_cor_unif(cardinal_env);// no correlation + uniformity
  //[pi]=tr_env_cor_proche_vois(cardinal_env); //cor with close neighbours
  //cor (rho)  with trend to grow
  //cor (rho) with trend to decrease
  //cor (rho) with a probability of fall
  //cor (rho) with a probability of rise
  //[pi]=tr_env_cor_crois(cardinal_env,1); //constant env: total correlation
  function b=dyn_plant_control(x,r,v)
    //b = total biomasse
    //x = vegetative biomass
    //r = environment
    //v = control
  function ci=offspring2(taille,env,control,time)
    // fitness of the plant
  // the final cost is zero

Question 2 Check that the list of matrices matrix_transition indeed consists of transition matrices, that is with nonnegative coefficients summing to 1 on each row.

2.5 Simulation of optimal trajectories

Copy the following code in the file plantII.sce.

  // computation of optimal trajectories (indexes)
  // optimal state trajectories in naturel units
  xset("window",2);xbasc();plot2d2(1:prod(size(e)),e,rect = [0,0,T+1,max(e)+1]);

Question 3 Execute different simulations. Change the matrix pi in the file plantII.sce.

3 Scilab Code

L'École des Ponts ParisTech est membre fondateur de

L'École des Ponts ParisTech est certifiée 

ParisTech ParisEst ParisTech


Copyright 2014
École des Ponts ParisTech
All rights reserved