Fermer X

Viable Harvesting of a Renewable Resource

Luc Doyen, Michel De Lara
(last modiﬁcation date: October 10, 2017)
Version pdf de ce document
Version sans bandeaux

1 The modelling

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)

 (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 function Sust is deﬁned by  (6)

• The maximum sustainable biomass BMSE and maximum sustainable yield hMSE are deﬁned by  (7)

• When hlim hMSE, the viability barrier BV is the solution of  (8)

Result 1 Assume that Biol is an increasing continuous function on +. The viability kernel is characterized 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 Paciﬁc yellowﬁn tuna, with intrinsic growth R = 2.25 metric tons per year and carrying capacity K = 250 000 metric tons. Since the carrying capacity solves Biol(K) = K, we obtain

 (12)

Question 1

1. Deﬁne a Scilab function for the dynamics Biol in (11), and a Scilab function for the equilibrium harvesting h = Sust(B) deﬁned in (6).
2. Plot the sustainable yield function BSust(B).
3. Compute the maximum sustainable population equilibrium BMSE and the maximum sustainable 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

1. Choose a guaranteed harvesting hlim strictly larger than Sust(BMSE).
2. For several initial conditions B0, compute diﬀerent trajectories for the smallest admissible harvesting, namely h(t) = hlim.
3. What happens with respect to viability constraint (5)?
4. Show that the situation is more catastrophic with sequences of harvesting h(t) larger than the guaranteed one hlim. Try for instance admissible policies:
• h(t,B) = B;
• h(t,B) = αhlim + (1 α)B with α [0, 1] (see Figure 1(a));
• h(t,B) = α(t)hlim+(1α(t))B, where (α(t))t is an i.i.d. sequence of uniform random variables in [0, 1] (see Figure 1(b)).

// number of random initial conditions
nbsimul=20;
// The unsustainable case:
h_lim=(1+rand()/5)*h_MSE

Horizon=30;
years=1:Horizon;
traj_B=[];

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

for i=1:nbsimul do
traj_B(1)=K*rand(1);
for t=years do
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,
end,
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;
traj_B=[];
traj_h=[];
alpha=rand(years);

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

for i=1:nbsimul do
traj_B(1)=K*rand(1);
for t=years do
traj_h(t)=policy(t,traj_B(t),alpha(t));
traj_B(t+1)=max(0,dynamics(traj_B(t)-traj_h(t)));
end,
end,
xtitle('Population trajectories violating the constraint','years (t)', ...
'biomass B (metric tons)')
legends('guaranteed harvesting',1,'ur')

3.3 The sustainable case: 𝕍iab≠∅

Question 3

1. Now choose a guaranteed harvesting hlim strictly smaller than hMSE = Sust(BMSE).
2. Compute the viability barrier BV .
3. Prove that the viable policies are those h(t,B) which lie within the set [hlim,h(B)] where
 (13)

4. For diﬀerent initial conditions B0, compute diﬀerent trajectories for a harvesting policy of the form
 (14)

5. What happens with respect to viability constraint (5)? See Figure 2.

// 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;

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:nbsimul do
traj_B(1)=K*rand(1);
for t=years do
traj_h(t)=viab(t,traj_B(t),alpha(t));
traj_B(t+1)=max(0,dynamics(traj_B(t)-traj_h(t)));
end,
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

4 Another population dynamics

The natural evolution is now governed by

 (15)

sqrt(max(0,K*B))

Question 4 Do the same as in the previous Beverton-Holt case.

 L'École des Ponts ParisTech est membre fondateur de L'École des Ponts ParisTech est certifiée