Scilab Code

Numerical implementations of the model

Using the scilab code below and the Excell datasheet Data_Guiana.xls, plot the historical catches 2006-2010 of the four fleets.

code 1

clear();
N_species=13; // number of fished species
N_fleet=4; // number of fleets
fleet=["Canots creoles";"Canots creoles ameliores ";"Pirogues";"Tapouilles"]; // French name of fleets

t_0=2006; // first data year
T=48;  // number of data months
nb_year_data=T/12;

sheets=readxls("Data_Guiana.xls"); data_fishing=sheets(1); // load excell sheet of data

//Catches in kg of species by four fleets
H_hist=[];
for i=1:N_species
    H_hist=[H_hist,data_fishing(3+(i-1)*(T+2):i*(T+2),4:4+N_fleet-1)];
end // size  H_hist= T*N_species*N_Fleet

function H_f=Catch_by_fleet(H) //Total catches in kg  by four fleets
    for fl=1:N_fleet
    H_fl=0;
    for i=1:N_species
        H_fl=H_fl+ H(:,fl+4*(i-1));
    end
    H_f(:,fl)=H_fl
end// size T*N_Fleet
endfunction

H_fleet_hist=Catch_by_fleet(H_hist)

clf(1);scf(1); // open and clear window 1
max_graphe=max(H_fleet_hist)./1000;
min_graphe=min(H_fleet_hist)./1000;
t_data_final=t_0+nb_year_data-1/12;
t_data=nb_year_data*12;
t_data_plot=[t_0:1/12:t_data_final];

for fl=1:N_fleet
  subplot(N_fleet/2,N_fleet/2,fl)
  xtitle(string(fleet(fl)),"","Catches  (tons)");
  plot2d(t_data_plot,H_fleet_hist(:,fl)./1000,style=[1],rect=[t_0,min_graphe-10,t_data_final+1,max_graphe+10]);
end

Plot and compare on the same figure the estimated catches 2006-2010 of the four fleets.

code 2

// data of fishing efforts  (hours at sea) by fleet
e_hist=data_fishing(3:2+T,9:9+N_fleet-1); // size e = 48 *4

Calib_r=sheets(2); r=Calib_r(2:N_species+2,2);
Calib_S=sheets(3); S=Calib_S(3:N_species+3,2:N_species+2)/10^(15);
Calib_Q=sheets(4); Q=Calib_Q(3:N_species+3,2:N_fleet+1)/10^(8);
Calib_B0=sheets(5); B_0=Calib_B0(2:N_species+2,2);

function B=dyn(B,e)
// Dynamics of the ecosystem
    B=(ones(B)+r+S*B-Q*e’).*B ;    ///  Lokta Volterra dynamics
    B=max(zeros(B),B); // Ensure positivity of biomasses
endfunction

function Y=Catch(B,e)
// Harvesting function (Schaefer type)
    Y=[];
    for (i=1:N_species)
     Y_i=Q(i,:)’.*e’*B(i); Y=[Y; Y_i];
    end
   Y=Y’; // transpose of the vector
endfunction

// Simulations provided by the model for 2006-2010
B(:,1)=B_0;//Initial biomass
for t = 1:t_data
    e(t,:)=e_hist(t,:);
    H(t,:)=Catch(B(:,t),e(t,:));
    B(:,t+1)= dyn(B(:,t),e(t,:));
end

// Plot by fleet of catches  estimated by the model
scf(1); // open and clear window 1
H_fleet=Catch_by_fleet(H)
for fl=1:N_fleet
  subplot(N_fleet/2,N_fleet/2,fl)
  plot2d(t_data_plot,H_fleet(:,fl)./1000,style=[2],rect=[t_0,min_graphe-10,t_data_final+1,max_graphe+10]);
end

Plot the biomass of the 14 species over the period 2006-2050 for the total closure scenario. What do you note regarding species viability ? What are the name of species at stake ?

code 3

t_horizon=2050;
month=(t_horizon-t_0)*12;
time=[t_0:1/12:t_horizon];

function e=effort_CL(t,B)
if t< t_data+1 then
e=e_hist(t,:);
else
e=zeros(1,N_fleet);
end;
endfunction

// Simulations Closure for 2010-2050
B_CL(:,1)=B_0;
for t = 1:month
e=effort_CL(t,B_CL(:,t))
e_hist(t,:)=e;
H_CL(t,:)=Catch(B_CL(:,t),e_hist(t,:));
B_CL(:,t+1)= dyn(B_CL(:,t),e_hist(t,:));
end

// Plot the biomass by species
clf(2);scf(2); // open and clear window 2

min_graph=0;
for s=1:N_species+1
max_graph=max(B_CL(s,:))*10^(-6);
subplot(int(N_species/3),int(N_species/3),s)
plot2d(time,B_CL(s,:)*10^(-6),style=[3],rect=[t_0,min_graph,max(time)+1,max_graph]);
xtitle("Species " +string(s),"","Biomass (Ktons)");
end

Using the functions Richness _index(B) and Trophic _index(B) below , plot both the species richness and the trophic index of the ecosystem over the period 2006-2050 for the total closure scenario. Use subplot function.Compute and plot also the Simpson index 4 by creating a function Simpson _index(B) and using subplot.

code 4

weight=[7 4 2.9 5 15 5 3.5 1.5 8 0.7 60 5 1] ;
//average weight of species (except 14)

function y =Richness_index(B)
  N=int(B(1:N_species)’./weight);
  S=bool2s(N>0)// test the positivity of abundance
y=sum(S);
endfunction

trophic=[4.05,4.35,4.03,4.2,4.5,3.3,4.5,4.04,4.11,3.5,4.09,2.13,2];
// Trophic level of species

function y =Trophic_index(B)
  N=int(B(1:N_species)’./weight);
  N_tot=sum(N);
  y=sum(N.*trophic)/N_tot;
endfunction

Plot the biomass of the different 14 species over the period 2006-2050 for the status quo (2009) scenario. Regarding biodiversity scores, plot also the trophic index and species richness

code 5

function e=effort_SQ(t,B)
  if t< t_data+1 then
    e=e_hist(t,:);
else
    e=mean(e_hist(t_data-11:t_data,:),"r");// mean effort 2009
    end;
endfunction

// Simulations  Status Quo for 2010-2050
B_SQ(:,1)=B_0;
for t = 1:month
    e_SQ(t,:)=effort_SQ(t,B_SQ(:,t))
    H_SQ(t,:)=Catch(B_SQ(:,t),e_SQ(t,:));
    B_SQ(:,t+1)= dyn(B_SQ(:,t),e_SQ(t,:));
end

// Plot the biomass by species  on window 2
scf(2);
min_graph=0;
for s=1:N_species
  subplot(int(N_species/3),int(N_species/3),s)
  max_graph=max(B_SQ(s,:))*10^(-6);
  plot2d(time,B_SQ(s,:)*10^(-6),style=[1],rect=[t_0,min_graph,max(time)+1,max_graph]);
   xtitle("Species " +string(s),"","Biomass  (Ktons)");
end

Using following scilab code relying on function profit(B,e), plot the four fleets’ rents. Use the unit variable costs, fixed costs, crew share earnings and prices as below.

code 6

Calib_cf=sheets(6); cf=Calib_cf(2:N_fleet+1,2); //monthly fixed cost
Calib_cv=sheets(7);
cv=Calib_cv(2:N_fleet+1,2); //monthly variable cost
sal_share=Calib_cv(2:N_fleet+1,3); //salary share
Calib_price=sheets(8); price=Calib_price(3:N_species+3,2:N_fleet+1);

function y=profit(B,e)
  H=Catch(B,e)  // catches
  H=matrix(H,N_fleet,N_species)’
  Income=H.*price(1:N_species,:);
  Income=sum(Income,"r")’
  Cost_var=e’.*cv;
  y=(1-sal_share).*(Income - Cost_var)-cf;  // including salaries
endfunction

// Simulations Status Quo for 2010-2050
B_SQ(:,1)=B_0;
for t = 1:month
    e_SQ(t,:)=effort_SQ(t,B_SQ(:,t))
    H_SQ(t,:)=Catch(B_SQ(:,t),e_SQ(t,:));
    rent_SQ(t,:)=profit(B_SQ(:,t),e_SQ(t,:))’
    B_SQ(:,t+1)= dyn(B_SQ(:,t),e_SQ(t,:));
end

// Plot the rent by fleet
clf(6);scf(6); // open and clear window 6

for f=1:N_fleet
  subplot(N_fleet/2,N_fleet/2,f)
  max_graph=max(rent_SQ(:,f))/1000;
  min_graph=min(rent_SQ(:,f))/1000;;
  plot2d(time(1:month),rent_SQ(:,f)’/1000,1,rect=[t_0,min_graph,max(time),max_graph]);
    plot2d(time(1:month),zeros(rent_SQ(:,f)’),5); // rent viability thresholds equals 0
   xtitle("fleet " +string(f),"","Profits (Keuros)");
end

Plot the global catch (proxy of food supply) by defining a function supply(B,e) as suggested below (Hint: use catch function Catch(B,e) as in profit function together with the function sum). Plot in red (color 5) the demand defined by H(2009)(1 + d)t where d ≈ 3.5 12 % stands for the monthly growth rate of the population. What do you observe for the SQ scenario in terms of food security ?

Implement the CVA strategy as suggested below. Compare through biodiversity scores (species richness, trophic index) and economic or productive scores with your own strategy as well as status quo strategy. Use colors green (3) for CL, black (1) for SQ, and blue (2) for CVA scenarios. What are the main changes underlying the CVA strategy in terms of fishing effort ?

Updated on 02/02/2016

Scilab Code

Marine Ecosystem in French Guiana

code 1

code 2

code 3

code 4

code 5

code 6

code 7

code 8

Our partners