Let us consider some regulating agency aiming at the sustainable use and harvesting of a
renewable resource. The biomass of this resource at time t is denoted by B(t) while the
harvesting level is h(t). We assume that the natural dynamics is described by some growth
function Biol. Under exploitation, the following controlled dynamics is obtained for
t ∈ ℕ
(1)
with the admissibility constraint
(2)
The policy goal is to guarantee at each time a minimal harvesting
(3)
together with a non extinction level for the resource
(4)
The combination of constraints (2) and (3) yields the state constraint
(5)
For the sake of simplicity, we assume that hlim> Blim.
2 The viability result
We use the viability approach to handle the problem. We aim at “showing” and recovering
numerically the characterization of the viability kernel described in Result 1. We introduce the
following notations.
The sustainable yield functionSust is defined by
(6)
The maximum sustainable biomass BMSE and maximum sustainable yield hMSE are
defined by
(7)
When hlim≤ hMSE, the viability barrier BV is the solution of
(8)
Result 1Assume thatBiolis an increasing continuous function on ℝ+. The viability kernel ischaracterized by
(9)
Viable controls are those h (depending on B) such that
(10)
3 The Beverton-Holt case
3.1 Dynamics and equilibrium
We consider the natural evolution governed by
(11)
For numerical computations and simulations, we consider the case of the Pacific yellowfin tuna,
with intrinsic growth R = 2.25 metric tons per year and carrying capacity K = 250000 metric
tons. Since the carrying capacity solves Biol(K) = K, we obtain
(12)
Question 1
Define a Scilab function for the dynamicsBiolin(11), and a Scilab function for theequilibrium harvesting h = Sust(B) defined in (6).
Plot the sustainable yield function BSust(B).
Compute the maximum sustainable population equilibrium BMSEand the maximumsustainable yield hMSE.
// Population dynamics parameters R_tuna=2.25; R=R_tuna; K_tuna=250000;// carrying capacity K=K_tuna; // BEVERTON-HOLT DYNAMICS b=(R-1)/K; function y=dynamics(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; disp('The maximum sustainable population equilibrium is '+string(B_MSE)+ ... ' metric tons (MT)'); h_MSE=SY(B_MSE); disp('The maximum sustainable yield is '+string(h_MSE)+' metric tons (MT)'); xset("window",10);xbasc(); // xbasc() is no longer valid with scilab, but works with scicos abcisse_B=linspace(0,K,20); plot2d(abcisse_B,SY(abcisse_B)) xtitle('Sustainable yield associated to Beverton-Holt dynamics','biomass (MT)', ... 'yield (MT)');
3.2 The unsustainable case: 𝕍iab = ∅
Question 2
Choose a guaranteed harvesting hlimstrictly larger thanSust(BMSE).
For several initial conditions B0, compute different trajectories for the smallestadmissible harvesting, namely h(t) = hlim.
What happens with respect to viability constraint (5)?
Show that the situation is more catastrophic with sequences of harvesting h(t) larger than theguaranteed one hlim. Try for instance admissible policies:
h(t,B) = B;
h(t,B) = αhlim + (1 − α)B with α ∈ [0, 1] (see Figure1(a));
h(t,B) = α(t)hlim+(1−α(t))B, where (α(t))t∈ℕis an i.i.d. sequence of uniformrandom variables in [0, 1] (see Figure1(b)).
// number of random initial conditions nbsimul=20; // The unsustainable case: h_lim=(1+rand()/5)*h_MSE Horizon=30; years=1:Horizon; yearss=[years,years($)+1]; traj_B=[]; xset("window",11);xbasc(); for i=1:nbsimuldo traj_B(1)=K*rand(1); for t=yearsdo traj_B(t+1)=max(0,dynamics(traj_B(t)-h_lim)); // max is no longer valid with scilab (use maxi), but works with scicos end, plot2d(yearss,traj_B); end, plot2d(yearss,h_lim*ones(yearss)) xtitle('Population trajectories violating the constraint','years (t)', ... 'biomass B (metric tons)') legends('guaranteed harvesting',1,'ur') ///////////////////////////////////////////////// function h=policy(t,B,alpha) h=max(h_lim,alpha*h_lim+(1-alpha)*B); endfunction Horizon=8; years=1:Horizon; yearss=[years,years($)+1]; traj_B=[]; traj_h=[]; alpha=rand(years); xset("window",12);xbasc(); for i=1:nbsimuldo traj_B(1)=K*rand(1); for t=yearsdo traj_h(t)=policy(t,traj_B(t),alpha(t)); traj_B(t+1)=max(0,dynamics(traj_B(t)-traj_h(t))); end, plot2d(yearss,traj_B); end, plot2d(yearss,h_lim*ones(yearss)) xtitle('Population trajectories violating the constraint','years (t)', ... 'biomass B (metric tons)') legends('guaranteed harvesting',1,'ur')
(a) Population viability constraint violated(stationary harvesting)
(b) Population viability constraint violated (randomharvesting)
Figure 1: Biomass trajectories violating the viability state constraint (5) (horizontal line
at hlim)
3.3 The sustainable case: 𝕍iab≠∅
Question 3
Now choose a guaranteed harvesting hlimstrictly smaller than hMSE =
Sust(BMSE).
Compute the viability barrier BV.
Prove that the viable policies are those h(t,B) which lie within the set [hlim,h♯(B)]
where
(13)
For different initial conditions B0, compute different trajectories for a harvesting policy of theform
(14)
What happens with respect to viability constraint (5)? See Figure2.
// The sustainable case h_lim=(1-rand()/5)*h_MSE ///////////////////////////////////////////////// // numerical estimation of the viability barrier function residual=viabbarrier(B) residual=SY(B)-h_lim 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_lim,alpha*h_lim+(1-alpha)*h_max(B)); endfunction Horizon=10; years=1:Horizon; yearss=[years,years($)+1]; xset("window",21);xbasc(); xtitle('Population trajectories ','years (t)','biomass B (metric tons)') traj_B=[]; xset("window",22);xbasc(); xtitle("Harvesting trajectories satisfying the constraint",'years (t)', ... 'harvesting h (metric tons)') traj_h=[]; alpha=rand(years); for i=1:nbsimuldo traj_B(1)=K*rand(1); for t=yearsdo traj_h(t)=viab(t,traj_B(t),alpha(t)); traj_B(t+1)=max(0,dynamics(traj_B(t)-traj_h(t))); end, xset("window",21);plot2d(yearss,traj_B); plot2d(yearss,B_V*ones(yearss)); legends('viability barrier',1,'ur') // xset("window",22);plot2d(years,traj_h); plot2d(years,h_lim*ones(years)); legends('guaranteed harvesting',1,'ur') plot2d(years,0*ones(years)); // for the scale end
(a) Population viability constraint satisfied orviolated
Figure 2: Population trajectories violating or satisfying the viability state constraint (5),
depending whether the original biomass is lower or greater than the viability barrier
(horizontal line at BV). Harvest trajectories satisfying the minimal harvesting viability
constraint (3) (horizontal line at hlim).
4 Another population dynamics
The natural evolution is now governed by
(15)
sqrt(max(0,K*B))
Question 4Do the same as in the previous Beverton-Holt case.