//\begin{verbatim}
// exec precaution.sce 

///////////////////////////////////////////////////////
//               Cost-effectiveness 
// for renewable resource management with learning
///////////////////////////////////////////////////////
//

function [y]=f(B,w)
// linear population dynamics
  y=w*B; 
endfunction

function z=dyna(B,e,w)
// controlled uncertain dynamics
  z=f(B*(1-e),w);
endfunction

function e=cost_effective(t,y)
// cost-effective effort
  w_r=w_hat;
  if t==1 then e=1-sqrt(rho*B_min./((w_r^2).*y(1)));end
  if t==2 then e=1-B_min./(y(2).*y(1));end
  e=max(0,e);
endfunction

function y=Obser(t,B,w)
// observation with learning
  y=zeros(1,2);
  if t==1 then y=[B,0];end
  if t==2 then y=[B,w];end
endfunction

xset("window",1);xbasc();
xtitle("Population state",'time t','biomass B(t)');

xset("window",2);xbasc();
xtitle("Catch effort",'time t','effort e(t)');

xset("window",3);xbasc();
xtitle("Uncertainty",'time t','uncertainty w(t)');

w_min=0.9; w_max=2;
// parameters dynamics 
w_hat=(w_max-w_min)/(log(w_max)-log(w_min));
// Certainty equivalent  of 1/w
Horizon=3; 
// two period problem
B_min=1; B_max=B_min*w_max/(1.1*w_min);
// Safety Constraint 
B_prec=B_min/(w_min^(Horizon-1));
// Precautionnary state 
rho=1/1.05;
//discount factor
B_simu=10;
// Simulation number

for i=1:B_simu
// Simulations
	B=B_prec+rand(1,Horizon)*(B_max-B_min);
	//  precautionnary initial conditions 
	// B(0)>= B_min/(w_min^(Horizon-1))	
  w=w_min+(w_max-w_min)*rand(1,Horizon-1);
  // Random or growth along time
	for t=1:Horizon-1
	// precautionary Trajectory B(.) e(.)
    y=Obser(t,B(t),w(t));
    e(t)=cost_effective(t,y);
    B(t+1)=dyna(B(t),e(t),w(t));
	end,

  rect1=[0,0,Horizon-1,(Horizon-1)*B_prec];
  rect2=[0,0,Horizon-2,1];
  rect3=[0,w_min,Horizon-2,w_max];
  xB=[0:1:Horizon-1];xe=[0:1:Horizon-2];
  //
  xset("window",1);
  plot2d(xB,[B' B_min+zeros(1,Horizon)' ...
			 B_max+zeros(1,Horizon)'],rect=rect1);
  // drawing diamonds, crosses, etc. to identify the curves
  // plot2d(xB,B',style=1);
  abcisse=linspace(0,Horizon-1,20);
  plot2d(abcisse,[B_min*ones(abcisse)' B_max*ones(abcisse)'],... 
			 style=-[4,5]);
 legends(["lower threshold";"upper threshold"],-[4,5],'lr');
 //  
  xset("window",2);plot2d(xe,e,rect=rect2);
  xset("window",3);plot2d(xe,w,rect=rect3);
end

// \end{verbatim}