Optimal Resources Allocation in Ecology:
to Grow or to Reproduce? (II)
The Stochastic Environment Case
Les questions

Michel de lara1
(last modification date: October 10, 2017)
Version pdf de ce document
Version sans bandeaux

1 Optimal Allocation in Random Environment (Theory)
2 Computation of optimal strategies in random environment (PW Scilab)
3 Scilab Code

Contents

1 Optimal Allocation in Random Environment (Theory)
 1.1 Optimisation model in random environment
 1.2 Resolution by stochastic dynamic programming
2 Computation of optimal strategies in random environment (PW Scilab)
 2.1 Discretization of the problem
 2.2 Using hypermatrices
  2.2.1 Environmental states transitions
  2.2.2 State transitions
 2.3 Conversion operations to fit with the Scilab function Bell_stoch
 2.4 Transitions and costs adapted to the Scilab function Bell_stoch
 2.5 Simulation of optimal trajectories
3 Scilab Code

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

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

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

Control

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,...,vT−1 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:

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

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
(2)

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

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

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

We take

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

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 }.
(5)

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
    dd=cardinal_alea
    pi=1/dd*ones(dd,dd)
  endfunction
  
  function pi=tr_env_cor_proche_vois(cardinal_alea)
    // correlation with closest neighbours
    dd=cardinal_alea
    pi=diag([0.5,1/3*ones(1:(dd-2)),0.5])+diag([0.5,1/3*ones(1:(dd-2))],1)+ ...
       diag([1/3*ones(1:(dd-2)),0.5],-1)
  endfunction
  
  function pi=tr_env_cor_crois(cardinal_alea,rho)
    // correlation with increasing trend
    // rho is the correlation intensity 0.8
    dd=cardinal_alea
    pi=diag([rho*ones(1:(dd-1)),1])+diag([(1-rho)*ones(1:(dd-1))],1)
  endfunction
  
  function pi=tr_env_cor_decrois(cardinal_alea,rho)
    // correlation with decreasing trend
    // rho is the correlation intensity 0.8
    dd=cardinal_alea
    pi=diag([1,rho*ones(1:(dd-1))])+diag([(1-rho)*ones(1:(dd-1))],-1)
  endfunction
  
  function pi=tr_env_cor_et_chute(cardinal_alea,rho)
    // correlation with probability of falling
    // rho is the correlation intensity 0.9
    dd=cardinal_alea
    pi=diag([rho*ones(1:dd)])
    pi(2:$,1)=1-rho
    pi(1,2)=1-rho
  endfunction
  
  function pi=tr_env_cor_et_ascension(cardinal_alea,rho)
    // correlation with a probability of rise
    // rho is the correlation intensity 0.9
    dd=cardinal_alea
    pi=diag([rho*ones(1:dd)])
    pi(1:$,$)=1-rho
    pi(dd,dd-1)=1-rho
  endfunction

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
    dims=[prod(size(taille)),prod(size(env))];
    ddims=[dims(1),dims(1),dims(2),dims(2)];
    Hypermatrices=list();
  
    cardinal_control=prod(size(control));
    cardinal_taille=prod(size(taille));
    for l=1:cardinal_control do
      //Hypermat=hypermat([ddims],sparse([],[],[prod(dims),prod(dims)]));
      Hypermat=hypermat([ddims]);
      // allows a natural expression for transitions from a product space
      // towards itself
  
      for y=1:dims(2) do
        image=controlled_dynamics(taille,y,control(l));
        // vector of the images of the vector state for r fixed
  
        indexes_image_discretisee=predec_succes(taille,image);
        indexes1=indexes_image_discretisee(1);
        indexes2=indexes_image_discretisee(2);
        probabilites=indexes_image_discretisee(3);
  
        M1=zeros(cardinal_taille,cardinal_taille);
        M2=zeros(M1);
        for i=1:cardinal_taille do
          M1(i,indexes1(i))=probabilites(i);
          M1(i,indexes2(i))=1-probabilites(i);
        end
  
        for z=1:dims(2) do
          Hypermat(:,1,y,z)=(1-ateb)*pi(y,z);
          // Alive (index>1), the plant faces a death probability
          // and the environment transitions are those of the matrix pi
          Hypermat(:,:,y,z)=Hypermat(:,:,y,z)+(M1+M2)*ateb*pi(y,z);
        end
      end
      // 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.
      Hypermatrices(l)=Hypermat
    end
  endfunction

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
    dd=size(Hypermatrices(1))
    dims=[dd(1),dd(3)]
    //dimensions of the hypermatrix
    Matrices=list()
    for i=1:prod(size(control)) do
      Mat=[]
      for j=1:dims(2) do
        if dims(2) <> 1 then
          a=Hypermatrices(i)(:,:,j,:).entries
          b=matrix(a,dims(1),prod(dims))
        else
          b=matrix(Hypermatrices(i)(:,:,j,:),dims(1),prod(dims))
        end
        Mat=[Mat;b]
        // Transformation into transition matrix (prod(dims),prod(dims))
        // This latter matrix M is adapted to the resolution of the SDPE
      end
      Matrices(i)=full(Mat)
    end
  endfunction
  
  
  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'
    else
      n=(j-1)*dims(1)+i;
    end
  endfunction
  
  
  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
    i=modulo(n,dims(1))
    ind=find(i==0)
    i(ind)=dims(1)
    j=int(n/dims(1))+1
    ind2=find((dims(1))\n==int((dims(1))\n))
    jj=(dims(1))\n
    j(ind2)=jj(ind2)
  endfunction
  
  
  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)
    Hmat=hypermat([dims,T-1])
    n=1:prod(dims)
    [i,j]=convert2(n,dims)
    for t=1:T-1 do
      for index=1:prod(dims) do
        Hmat(i(index),j(index),t)=matrix(n(index),t)
      end
    end
  endfunction
  
  
  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)
    dims=[prod(size(taille)),prod(size(env))]
    instant_cost=list()
    for j=1:prod(size(control)) do
      cost_i=[]
      for t=1:horizon do
        // cost_it = hypermat(dims,sparse([],[],[dims(1),dims(2)]));
        cost_it=hypermat(dims);
        for state1=1:dims(1) do
          for state2=1:dims(2) do
            cost_it(state1,state2)=cost(taille(state1),env(state2),control(j),t)
          end
        end
        cost_it=matrix(cost_it,1,prod(dims))
        // Transformation into vector : size(cost_i)=[1,prod(dims)]
        // adapted to the resolution of the SDPE
        cost_i=[cost_i(cost_it)']
      end
      instant_cost(j)=full(cost_i)
    end
  endfunction

2.4 Transitions and costs adapted to the Scilab function Bell_stoch

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

  exec('plantI2.sci');
  exec('plantII.sci');
  
  T=10;
  ateb=0.9;
  rr=1.2;
  power=0.5;
  xp=(ateb*rr*power)^{1/(1-power)};
  xm=(xp/rr)^{1/power};
  taille=[0:(xm/3):(1.2*xp)];
  //taille=0:2;
  env=[0.6,0.8,1,1.2,1.4];
  // env=[10:12]
  control=0:0.05:1;
  // control=0:0.5:1;
  
  cardinal_env=prod(size(env));
  cardinal_taille=prod(size(taille));
  
  pi=tr_env_auc_cor_unif(cardinal_env);// no correlation + uniformity
  //[pi]=tr_env_cor_proche_vois(cardinal_env); //cor with close neighbours
  //[pi]=tr_env_cor_crois(cardinal_env,rho);
  //cor (rho)  with trend to grow
  //[pi]=tr_env_cor_decrois(cardinal_env,rho);
  //cor (rho) with trend to decrease
  //[pi]=tr_env_cor_et_chute(cardinal_env,rho);
  //cor (rho) with a probability of fall
  //[pi]=tr_env_cor_et_ascension(cardinal_env,rho);
  //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
    b=v*r*x^{power}
  endfunction
  
  controlled_dynamics=dyn_plant_control;
  
  Hypermatrices=constr_HyperM(taille,env,control,pi,controlled_dynamics);
  
  transition_matrix=conv_HyperM(Hypermatrices);
  
  
  function ci=offspring2(taille,env,control,time)
    // fitness of the plant
    ci=(1-control)'*env*taille^(power)*ateb^{time-1}
    ind=find(taille==1),ci(ind)=0
  endfunction
  
  offspring=conv_cost(taille,env,control,offspring2,T);
  
  final_cost_nul=zeros(cardinal_env*cardinal_taille,1);
  // 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.

  dims=[cardinal_taille,cardinal_env];
  
  stacksize(3000000);
  
  instant_cost=offspring;
  final_cost=final_cost_nul;
  
  [value,feedback]=Bell_stoch(transition_matrix,instant_cost,final_cost,1);
  
  index_taille_init=grand(1,1,'uin',2,cardinal_taille);
  index_env_init=grand(1,1,'uin',1,cardinal_env);
  
  initial_state=convert1(index_taille_init,index_env_init,dims);
  
  z=trajopt(transition_matrix,feedback,instant_cost,final_cost,initial_state);
  
  // computation of optimal trajectories (indexes)
  
  [index_taille,index_env]=convert2(z(1),dims)
  
  x=taille(index_taille);
  e=env(index_env);
  // optimal state trajectories in naturel units
  
  xset("window",1);xbasc();plot2d2(1:prod(size(x)),x);
  xtitle("taille")
  
  xset("window",2);xbasc();plot2d2(1:prod(size(e)),e,rect = [0,0,T+1,max(e)+1]);
  xtitle("environment")
  
  xset("window",3);xbasc();plot2d2(1:prod(size(control(z(2)))),control(z(2)));
  xtitle("control")
  
  xset("window",4);xbasc();plot2d2(1:prod(size(z(3))),z(3));
  xtitle("fitness")

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

3 Scilab Code

plantII.sci plantII.sce