function[b]=dyn_plant(x)
// b = total biomass generated by the plant in one period
//x = vegetative biomass
     b=r*x^{power}
endfunction


function[b]=dyn_plant_control(x,v,s)
// Growth function, with control v and random term s
// b = future vegetative biomass
// x = vegetative biomass
// v = fractiwe havellocated to growth (control)
// s = survival factor (s=1 survival, s=0 death)
// if s=0, x transits towards 0 
// if s=1, x transits towards v*dyn\_plant(x)
  if v < 0 | v > 1, 
     b = 'ERROR : control beyond bounds'
  else
           if s==0,
             b=0 
           //transition towards 0 if the survival factor is zero
           else  
           b=v*dyn_plant(x) 
            end 
  end 
endfunction


function [cost]=offspring(state,control,time)
cost=(1-control)*ateb^{time}*dyn_plant(state);
endfunction

function[cost]=final_cost_zero(state,time)
cost=0;
endfunction

function[v]=strategy_D(t,x)
v=0
endfunction

function[v]=strategy_R(t,x)
v=rand()
endfunction

function[v]=strategy_G(t,x)
  if t<T then
// t=1:T corresponds to t=0,...,T-1
     v=1
  else
     v=0
  end
endfunction


function[x,v]=traj_feedback...
(initial_state,controlled_dynamics,feedback,alea)
// initial_state
// controlled_dynamics(x,v,s) : returns a state
// feedback(t,x) : returns a control
// alea : sequence of 0 or 1, with length horizon
// x : state trajectory
// v : control trajectory
  horizon=size(alea,2);
  x=initial_state
  v=[]
  for t=1:horizon
  // t=0,...,T-1
    v=[v,feedback(x($),t)]
    // x($) is the last value of vector x
    // v has dimension horizon=T
    x=[x,controlled_dynamics(x($),v($),alea(t))]
    // x has dimension 1+horizon=1+T
  end
endfunction


function[cost]=total_cost(trajectory,control,inst_cost,final_cost)
// trajectory is a vector of dimension 1+horizon, 
// control is a vector of control of dimension horizon
// inst_cost(state,control,time) is a function 
// final_cost(state,time) is a function 
// cost is a vector of dimension 1+horizon, comprising the sequence of
// cumulated sums of instantaneous cost (and of final cost)

  horizon=prod(size(control))

if 1+horizon ~= prod(size(trajectory)),
     cost = 'ERROR : dim(trajectory) different from 1+dim(control)'
else
  cost=[]
  for t=1:horizon
// i.e. t=0,...,T-1
    cost=[cost,cost($)+ inst_cost(trajectory(t),control(t),t)] 
// cumulated sums of instantaneous cost
  end
  cost=[cost,cost($)+final_cost(trajectory(1+horizon),1+horizon)]
end
endfunction


function[v]=strategy_O(t,x)
  if x<xm then
     v=1
  else
     v=xp/dyn_plant(x)
  end
endfunction