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

/////////////////////////////////////////////
// Population dynamics 
/////////////////////////////////////////////

// Population dynamics parameters
R=1.5; 	
b=0.01;

function [y]=g(B)  
	y=R*B./(1+b*B) 
endfunction

function h=SY(B)  
	h= B-(B./(R-b*B)) 				 
endfunction

B_MSE=(R-sqrt(R))/b
h_MSE=SY(B_MSE)
K=(R-1)/b

xset("window",0); xbasc();
B=[0:1:K]; plot2d(B,SY(B))
xtitle('Sustainable yield curve associated to ...
	   Beverton-Holt dynamics','biomass B')


// number of random initial conditions
nbsimul=30;


/////////////////////////////////////////////
// The unsustainable case
/////////////////////////////////////////////

h_min= h_MSE +0.5	// for instance 

Horizon=60;
B=zeros(1,Horizon+1);


// // Stationary harvesting
xset("window",11); xbasc();

for i=1:nbsimul
	B(1)=K*rand(1);
	for t=1:Horizon
		B(t+1)=max(0,g(B(t)-h_min));
	end,
	plot2d(1:(Horizon+1),B);	
end,
plot2d(1:(Horizon+1),h_min * ones(1:(Horizon+1)),style=-4)
legends("Viability threshold",-4,'ur')

xtitle('Biomass trajectories violating the constraint',...
	   'time t','biomass B(t)')


// // Random harvesting

function h=feedback(t,B,alpha)  
	h=max(h_min,alpha*h_min+(1-alpha)*B); 				 
endfunction

Horizon=10;

alpha=rand(1,Horizon);
B=zeros(1,Horizon+1);

xset("window",12); xbasc();

for i=1:nbsimul
	B(1)=K*rand(1);
	for t=1:Horizon
		h(t)=feedback(t,B(t),alpha(t));
		B(t+1)=max(0,g(B(t)-h(t)));
	end,
	plot2d(1:(Horizon+1),B);	
end,
plot2d(1:0.2:(Horizon+1),...
	   h_min*ones(1:0.2:(Horizon+1)),style=-4)
legends("Viability threshold",-4,'ur')

xtitle('Biomass trajectories violating the constraint',...
	   'time t','biomass B(t)')



/////////////////////////////////////////////
// The sustainable case
/////////////////////////////////////////////

h_min= h_MSE *0.5	// for instance 

Horizon=10;

// numerical estimation of the viability barrier

function error=viabbarrier(B)
error=SY(B)-h_min
endfunction

B_V=fsolve(0,viabbarrier)



function h=h_max(B)
    h= B- ( B_V / (R - b*B_V) )
endfunction
    
function h=viab(t,B,alpha)  
	h=max(h_min,alpha*h_min+(1-alpha)*h_max(B))
endfunction

Horizon=10;

xset("window",21); xbasc();
xtitle("Biomass trajectories",'time t','biomass B');

xset("window",22);xbasc();
xtitle("Harvesting trajectories ...
	   satisfying the constraint",'time t','harvest h');

alpha=rand(1,Horizon);
B=zeros(1,Horizon+1);
h=zeros(1,Horizon);

for i=1:nbsimul
	B(1)=K*rand(1)/2;
	for t=1:Horizon
		h(t)=viab(t,B(t),alpha(t));
		B(t+1)=max(0,g(B(t)-h(t)));
	end,
	xset("window",21);plot2d(1:(Horizon+1),B);
		plot2d(1:0.2:(Horizon+1),...
			   B_V*ones(1:0.2:(Horizon+1)),style=-4)
		legends("Viability barrier",-4,'ur')
		//
		xset("window",22);plot2d(1:Horizon,h);	
		plot2d(1:0.2:Horizon,h_min*ones(1:0.2:Horizon),style=-4)
		legends("Minimal harvest",-4,'ur')
end

// \end{verbatim}