Fermer X

Viable Harvesting of a Renewable Resource

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

Contents

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

B (t + 1) = Biol(B (t) − h (t))
(1)

with the admissibility constraint

0 ≤ h(t) ≤ B(t).
(2)

The policy goal is to guarantee at each time a minimal harvesting

hlim ≤  h(t),
(3)

together with a non extinction level for the resource

Blim ≤  B (t).
(4)

The combination of constraints (2) and (3) yields the state constraint

hlim ≤  B(t).
(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 defined by
    h = Sust (B ) ⇐ ⇒  B = Biol (B −  h) and 0 ≤ h ≤ B.
    (6)

  • The maximum sustainable biomass BMSE and maximum sustainable yield hMSE are defined by
    hMSE  = Sust (BMSE ) = max Sust (B ).
                       B ≥0
    (7)

  • When hlim hMSE, the viability barrier BV is the solution of
    BV  =     min     B.
      B,hlim=Sust(B)
    (8)

Result 1 Assume that Biol is an increasing continuous function on +. The viability kernel is characterized by

       {
          [BV ,+ ∞ [  if hlim ≤ hMSE
𝕍iab =
          ∅           otherwise.
(9)

Viable controls are those h (depending on B) such that

Biol (B −  h) ∈ [BV ,+ ∞ [ and 0 ≤ h ≤  B.
(10)

3 The Beverton-Holt case

3.1 Dynamics and equilibrium

We consider the natural evolution governed by

            -RB----
Biol (B ) = 1 + bB .
(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 = 250 000 metric tons. Since the carrying capacity solves Biol(K) = K, we obtain

R  = 2.25 metric tons per year and  b = R-−-1-=  5 10−6 (metric tons per year)− 1.
                                         K
(12)

Question 1

  1. Define a Scilab function for the dynamics Biol in (11), and a Scilab function for the equilibrium harvesting h = Sust(B) defined in (6).
  2. Plot the sustainable yield function B↦→Sust(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 different 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;
  yearss=[years,years($)+1];
  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,
    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: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,
    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')

PIC
(a) Population viability constraint violated (stationary harvesting)
PIC
(b) Population viability constraint violated (random harvesting)

Figure 1: Biomass trajectories violating the viability state constraint (5) (horizontal line at hlim)

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
    h♯(B ) = B − ---BV----.
             R  − bBV
    (13)

  4. For different initial conditions B0, compute different trajectories for a harvesting policy of the form
                                   ♯
h (t,B ) = α(t)hlim + (1 − α (t))h (B ).
    (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;
  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: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,
    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

PIC
(a) Population viability constraint satisfied or violated
PIC
(b) Minimal guaranteed harvesting viability constraint satisfied

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

           √ ----
Biol(B ) =   KB.
(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 

ParisTech ParisEst ParisTech

 

Copyright 2014
École des Ponts ParisTech
All rights reserved