The Bio-economic Model:

We consider a fishery in which two fishing fleets (a fleet of trawlers and a fleet of gill-netters) exploit a stock of soles (Solea solea). The simplified multi-fleet bio-economic model is described in discrete time with an annual step.

Code Modèle

///// BIO-ECONOMIC MODEL
///// Fishery with 1 species (sole) et 2 fishing fleets (trawlers and gill-netters)
//The bio-economic model is in discrete time with an annual time step
N_species=1; // sole
N_fleets=2; // f=1 : trawlers ; f=2 : gill-netters
t_0=2008; // Initial time for the simulation

Biological Parameters: 

BioParameters

// BIOLOGICAL PARAMETERS - SOLE ----------------------------------
// Number of age classes for the sole in the model (model structured in age)
Age=7;
// Abundances at t0 (in 2008) at different age
N_0=[23191;17416;10707;4864;3425;2627;2590]; // in thousands of individuals
// Natural mortality
M=ones(Age,1)*0.3;
// Mean annual recruitment
Rbar=23414;
// Average weight by ages (kg/individual)
W_stock= [0.189;0.241;0.297;0.352;0.423;0.449;0.599];
// Proportion of mature individuals by age
O=[0.32;0.83;0.97;1;1;1;1];

Economical Parameters: 

EcoParameters

TECHNICAl-ECONOMICAL PARAMETERS
// Mean market price by age and per kilo
price=[8.6;10.2;12.3;14.2;14.2;14.2;13.9]; // euro per kilo
///////////// Fishing mortality (related to the 2008 effort levels)
// Fishing mortality by age due to trawlers
F_trawlers=[0.0486085 ; 0.0614180; 0.0463817; 0.0404101; 0.0392817; 0.0393315; 0.0393021];
// Fishing mortality by age due to gill-netters
F_gillnetters=[0.0431997; 0.1422847; 0.1957702; 0.2115014; 0.3135258; 0.3084042; 0.3037340];
F0_a_f=[F_trawlers,F_gillnetters];
// Initial effort per fleet
K_0=[241,133]; // number of vessels // first: trawlers // second: gill-netters
e_0=[175,171] // mean annual number of days at sea per vessel
///Economic data
V_cost=[481,421]; // variable cost per effort unit (in euros per day at sea)
F_cost=[113419,124627]; // annual fixed cost per vessel
// Proportion of income coming from catches of other species (not explicitly modelled here)

alpha_Inc_autres_sp=[5.22, 0.64];

The Spawning Biomass Function (SSB):

SSB

function SSB=Spawning_biomass(N)
// N and W_stocks are 2 "column vectors"
//SSB: spawning stock biomass in tonnes (if N is in thousands of individuals)
SSB=sum(O.*W_stock.*N);
endfunction

The Stock Dynamic Function aims  to return the abundances by age at time t+1 based on abundance at time t and according to a vector of fishing effort multiplier: u.

DynStock

function y=dyn_stock(N,u)
y(1,1)=Rbar;
u_a_f=ones(Age,1)*u; // matrix of fishing effort multiplier. Dimensions: number of ages x number of fleets
F_a_total=sum(F0_a_f.*u_a_f,'c'); // fishing mortality per age
Z=M+F_a_total; //Total mortality per age
y(2:Age-1,1)=int(N(1:Age-2).*exp(-Z(1:Age-2))); // int() to have integer.
y(Age,1)=int(N(Age-1)*exp(-Z(Age-1))+N(Age)*exp(-Z(Age)));
y=max(0,y); // "y" corresponds to a vector of abundances by age

endfunction

The Production Function that calculates the annual catches per fleet

Production

// Production / Catches
function C=production(N,u);
F_a_f=F0_a_f.*(ones(Age,1)*u); // fishing mortalities per age and per fleet.
//Dimensions: number of ages * number of fleets
F_a=sum(F_a_f,'c'); // total fishing mortality by age
Z=M+F_a; // total mortality by age
C_trawlers=int(N.*(F_a_f(:,1).*((1-exp(-Z))./Z)));
C_gillnetters=int(N.*(F_a_f(:,2).*((1-exp(-Z))./Z)));
C=[C_trawlers,C_gillnetters];
endfunction

Function Income that calculates the annual income per fleet

Income

//Income per fleet
function Inc=Calcul_Income(N,Price,u)
y=production(N,u);
Income_trawlers=y(:,1)'*(W_stock.*Price); // in thousands of euros
Income_gillnetters=y(:,2)'*(W_stock.*Price); // in thousands of euros
Inc=[Income_trawlers,Income_gillnetters];
endfunction

Function Profit which calculates the annual profit per fleet, the effort multiplier u is applied to the number of vessels

Profit

//PROFITS per fleet
function Pi=Calcul_Profits(N,Price,u)
Inc=Calcul_Income(N,Price,u);
effort_total_f=K_0.*e_0.*u;
Inc_tot=Inc.*(ones(1,N_fleets)+alpha_Inc_autres_sp)*10^3; // en euros
Pi=Inc_tot-(V_cost.*effort_total_f)-(F_cost.*(K_0.*u));
Pi=Pi/10^3; // en milliers d'euros
endfunction

Using the function Trajectories , you can plot the trajectories of the spawning stock biomass and profits for a period of 10 years, with a status quo scenario 

Trajectories

//// Function  trajectories
function [N,SSB,Inc,Profit]=trajectories(N_0,u,Price,T)
N=[N_0];
SSB=[Spawning_biomass(N_0)];
Inc=[Calcul_Income(N_0,Price,ones(1,2))];
Profit=[Calcul_Profits(N_0,Price,ones(1,2))];
for i=1:T
Nt=dyn_stock(N(:,i),u);
N=[N, Nt];
SSB=[SSB,Spawning_biomass(Nt)];
Inc=[Inc;Calcul_Income(Nt,Price,u)]; //in thousands of euros
Profit=[Profit;Calcul_Profits(Nt,Price,u)]; //in thousands of euros
end
endfunction

//// Function to plot trajectories
function plot_trajectories(N_0,u,Price,T,num_window)
[Stock_a_t,SSB_t,Inc_t_f,Profit_t_f]=trajectories(N_0,u,Price,T);
clf(num_window);
xset("window",num_window);
//// SSB evolution
subplot(1,2,1)
aa=get("current_axes");
aa.font_style = 4 ; aa.font_size = 3;
xl = aa.x_label ; xl.font_style = 4 ; xl.font_size = 3 ;
yl = aa.y_label ; yl.font_style = 4 ; yl.font_size = 3 ;
tl=aa.title; tl.font_style=4; tl.font_size=3;
plot2d([t_0:t_0+T],[SSB_t'], style=[2],rect=[t_0,min(5000,min(SSB_t)),t_0+T,max(SSB_t)]);
bb= gce();
bb.children.thickness=2;
xtitle('','time','sole SSB (in tons)')

//// Profits evolution
subplot(1,2,2)
aa=get("current_axes");
aa.font_style = 4 ; aa.font_size = 3;
xl = aa.x_label ; xl.font_style = 4 ; xl.font_size = 3 ;
yl = aa.y_label ; yl.font_style = 4 ; yl.font_size = 3 ;
tl=aa.title; tl.font_style=4; tl.font_size=3;
plot2d([t_0:t_0+T],[Profit_t_f, zeros(T+1,1)], style=[13,27,1]);
bb= gce();
bb.children.thickness=2;
bb.children(1).line_style=3; // dashed line
xtitle('u= '+string(u)+'','time','Profits (in thousands of euros)')
legend('trawlers','gill-netters')
endfunction

u_SQ=[1,1]; // effort equal to the one in 2008
T=10; // temporal horizon
plot_trajectories(N_0,u_SQ,price,T,1)

Stochasticity

Here we add stochasticity to the recruitment. We create a recruitment matrix  to compare different management strategies.

Matrix

R_historical=[25137 26434 23797 29612 23792 22647 24494 25167 16670 24615 24281 19723 18937 22490];
N_simu_recruitment=50; // number of simulations
max_r=length(R_historical);
recruitment_list=grand(T,N_simu_recruitment,'uin',1,max_r);
matrix_R_random=matrix(R_historical([recruitment_list]),T,N_simu_recruitment);

Random Dynamic stock function and Random trajectories function.

RandomStockTrajectories

// Dynamic function where the recruitment value can change and is not fixed to Rbar anymore
function y=dyn_stock_random(N,u,R_random)
y(1,1)=R_random;
u_a_f=ones(Age,1)*u;
F_total=sum(F0_a_f.*u_a_f,'c');
Z=M+F_total; // Total mortality
y(2:Age-1,1)=int(N(1:Age-2).*exp(-Z(1:Age-2)));
y(Age,1)=int(N(Age-1)*exp(-Z(Age-1))+N(Age)*exp(-Z(Age)));
y=max(0,y);
endfunction

function [N,SSB,Inc,Profit]=trajectories_random(N_0,u,Price,T,vec_R_random)
N=[N_0];
SSB=[Spawning_biomass(N_0)];
Inc=[Calcul_Income(N_0,Price,ones(1,2))];
Profit=[Calcul_Profits(N_0,Price,ones(1,2))];
for i=1:T
Nt=dyn_stock_random(N(:,i),u,vec_R_random(i,1));
N=[N, Nt];
SSB=[SSB,Spawning_biomass(Nt)];
Inc=[Inc;Calcul_Income(Nt,Price,u)]; //in thousands of euros

Profit=[Profit;Calcul_Profits(Nt,Price,u)]; //in thousands of euros
end
endfunction

Function Plot_traj_random, which plots the different SSB and profits trajectories with stochastic recruitment

Trajectory plot random

function plot_traj_random(N_0,u,Price,T,matrix_R,nb_simu,num_window)
[Stock_a_t,SSB_t,Inc_t_f,Profit_t_f]=trajectories_random(N_0,u,Price,T,matrix_R(:,1));
clf(num_window);
xset("window",num_window);
subplot(1,2,2)
aa=get("current_axes");
aa.font_style = 4 ; aa.font_size = 3;
xl = aa.x_label ; xl.font_style = 4 ; xl.font_size = 3 ;
yl = aa.y_label ; yl.font_style = 4 ; yl.font_size = 3 ;
tl=aa.title; tl.font_style=4; tl.font_size=3;
plot2d([t_0:t_0+T],[Profit_t_f, zeros(T+1,1)], style=[13,27,1]);
bb= gce();

bb.children.thickness=2;
bb.children(1).line_style=3; // dashed line
xtitle('u= '+string(u)+'','time','Profits (in thousands of euros)')
legend('trawlers','gillnetters')
for i=1:nb_simu
[Stock_a_t,SSB_t,Inc_t_f,Profit_t_f]=trajectories_random(N_0,u,Price,T,matrix_R(:,i));
// SSB
subplot(1,2,1)
plot2d([t_0:t_0+T],[SSB_t'], style=[2],rect=[t_0,min(5000,min(SSB_t)),t_0+T,max(SSB_t)]);
bb= gce();
bb.children.thickness=2;
// profits
subplot(1,2,2)
plot2d([t_0:t_0+T],[Profit_t_f], style=[13,27]);
end

subplot(1,2,1)
aa=get("current_axes");
aa.font_style = 4 ; aa.font_size = 3;
xl = aa.x_label ; xl.font_style = 4 ; xl.font_size = 3 ;
yl = aa.y_label ; yl.font_style = 4 ; yl.font_size = 3 ;
tl=aa.title; tl.font_style=4; tl.font_size=3;
xtitle('','time','sole SSB (in tons)');
endfunction

Viability approach

Biological constraint: SSB(t) is maintained above a threshold value B

pa

Economic constraint:  maintaining positive profits for each fleet over time ∏f(t)>0

The co-viability probability  of the fishery under a statu quo scenario with random recruitment:

Status Quo

Blim=9706; // corresponds to Bloss which is the smallest observed value of SSB (here it is the SSB of 1998)
Bpa=13000; // precautionary biomass, corresponding to Blim*1.39
///// CO-VIABILITY PROBABILITY CALCUL
u_SQ=[1,1];
PVA_tot=0;
EVA_tot=0;
CVA_tot=0;
for simu=1:N_simu_recruitment // loop for the different simulations of recruitment
[Stock_a_t,SSB_t,Inc_t_f,Profit_t_f]=trajectories_random(N_0,u_SQ,price,T,matrix_R_random(:,simu));
// we allow the first time step to be below the constraints (it is why we have "2:T+1")
// check if SSB is above the threshold or not for each time
PVAt=bool2s((SSB_t(:,2:T+1) -Bpa*ones(1,T) >=0)) ;
// check if profits of each fleet at each time is strictly superior to 0 or not
EVAt_i=bool2s((Profit_t_f(2:T+1,:) >0));
// we do the product, because we want that the profit o each fleet is strictly positive and not only for one.
EVAt=prod(EVAt_i,'c');
PVA_simu=prod(PVAt);
EVA_simu=prod(EVAt);
PVA_tot=PVA_tot+PVA_simu;
EVA_tot=EVA_tot+EVA_simu;
CVA_tot=CVA_tot+prod([PVA_simu;EVA_simu]);
end
PVA=PVA_tot/N_simu_recruitment // average - we have a value between 0 and 1.
EVA=EVA_tot/N_simu_recruitment
CVA=CVA_tot/N_simu_recruitment

Modified "plot trajecotries random" function to display the viability thresholds.

Viability thresholds

function plot_traj_random(N_0,u,Price,T,matrix_R,nb_simu,num_window)
[Stock_a_t,SSB_t,Inc_t_f,Profit_t_f]=trajectories_random(N_0,u,Price,T,matrix_R(:,1));

clf(num_window);
xset("window",num_window);
subplot(1,2,2)
aa=get("current_axes");
aa.font_style = 4 ; aa.font_size = 3;
xl = aa.x_label ; xl.font_style = 4 ; xl.font_size = 3 ;
yl = aa.y_label ; yl.font_style = 4 ; yl.font_size = 3 ;
tl=aa.title; tl.font_style=4; tl.font_size=3;
plot2d([t_0:t_0+T],[Profit_t_f, zeros(T+1,1)], style=[13,27,1]);
bb= gce();
bb.children.thickness=2;
bb.children(1).line_style=3; // dashed line
xtitle('u= '+string(u)+'','time','Profits (in thousands of euros)')
legend('trawlers','gill-netters')
for i=1:nb_simu
[Stock_a_t,SSB_t,Inc_t_f,Profit_t_f]=trajectories_random(N_0,u,Price,T,matrix_R(:,i));
// SSB
subplot(1,2,1)
plot2d([t_0:t_0+T],[ones(T+1,1)*Bpa, SSB_t'], style=[1,2],rect=[t_0,min(5000,min(SSB_t)),t_0+T,max(SSB_t)]);
bb= gce();
bb.children.thickness=2;
// profits
subplot(1,2,2)
plot2d([t_0:t_0+T],[Profit_t_f], style=[13,27]);
end

subplot(1,2,1)
aa=get("current_axes");
aa.font_style = 4 ; aa.font_size = 3;
xl = aa.x_label ; xl.font_style = 4 ; xl.font_size = 3 ;
yl = aa.y_label ; yl.font_style = 4 ; yl.font_size = 3 ;
tl=aa.title; tl.font_style=4; tl.font_size=3;
xtitle('','time','sole SSB (in tons)');
legend('Bpa','SSB',3)
endfunction

/// SOME EXAMPLE OF MANAGEMENT STRATEGIES

/// trajectories with effort multipliers associated to a probability of co-viability equal to 1
plot_traj_random(N_0,[1.2, 0.3],price,T,matrix_R_random,N_simu_recruitment,22)
plot_traj_random(N_0,[0.4, 0.5],price,T,matrix_R_random,N_simu_recruitment,23)
/// trajectories with effort multipliers associated to a probability of co-viability not equal to 0.28
plot_traj_random(N_0,[1.2, 0.6],price,T,matrix_R_random,N_simu_recruitment,24)
/// trajectories with effort multipliers associated to a probability of co-viability not equal to 0.5
plot_traj_random(N_0,[0.2, 0.8],price,T,matrix_R_random,N_simu_recruitment,24)

Updated on 24/08/2016