////////////////////////////////
//
// Un modele structuré en âge
// Stratégie optimale d'abattage
//
// (taper resume pour passer l'étape suivante)
////////////////////////////////
// ETAPE 1 : evolution d'un troupeau sans abattage
//////////////////////////////
n=12;
F0=[0,0.045,0.391,0.472, 0.484, 0.546, 0.511, 0.412, 0.352, 0.122, 0.088, 0.0];
//
T0=[0.845,0.975 0.965 0.950 0.926 0.895 0.850 0.786 0.691 0.561 0.370];
//
x0init=10*[10;12;4;3;5;4;3;4;5;2;1;1]  // troupeau initial
//    
profit=[1,1,1,1];         // evolution du bénéfice en fonction de l'age du mouton
////////////////////////////////////
figure(1)
clf()
plot(1:n,F0,'b-*','thickness',3)
plot(1:(n-1),T0,'g-*','thickness',5)
xtitle('taux de fertilité et de survie')
////////////////////////////////////
T=diag(T0,-1)
F=zeros(n,n)
F(1,:)=F0;
A=F+T
x0=x0init;
Xt=x0/norm(x0,1);Nt=[norm(x0,1)]
Nit=30;
for i=1:Nit
    x0=A*x0
    N=norm(x0,1)
    Xt=[Xt,x0/norm(x0,1)]
    Nt=[Nt,norm(x0,1)]
end
///////////////////////////////////////
disp('plus grande valeur propre valeur de A');disp(max(abs(spec(A))));
disp('population finale après 30 ans (sans abattage)'); disp(norm(x0,1))
figure(2)
clf()
plot(0:Nit,Xt,'thickness',5)
xtitle('repartition par classe (sans abattage)')
pause
////////////////////////////////////////////
// ETAPE 2 : stratégie d'abattage optimale
/////////////////////////////////////////////
//
Ftot=[];Phitot=[];
//
for age=1:4        // abattage possible entre 1 et 4 ans
//////////////
   deff('R=abattage(phi)','Tbis=T0,Tbis(age)=Tbis(age)-phi*Tbis(age),R=F0(n),for i=n-1         :-1:1,R=F0(i)+Tbis(i)*R,end,R=R-1')
   phiaba=fsolve(0.1,abattage) // taux d'abattage pour garder troupeau constant
   Phitot=[Phitot,phiaba];  
   //////////////
   T0bis=T0;T0bis(age)=T0bis(age)-phiaba*T0bis(age);
   T=diag(T0bis,-1)
   F=zeros(n,n)
   F(1,:)=F0;
   A=F+T;
   x0=x0init;Xt=x0/norm(x0,1);Nt=[norm(x0,1)]
   Nit=50;
   for i=1:Nit
    x0=A*x0;N=norm(x0,1)
    Xt=[Xt,x0]
    Nt=[Nt,norm(x0,1)]
   end
   //
   frac=x0(age)/norm(x0,1)
   Ftot=[Ftot,frac];  /// nouvelle fraction limite de la classe d'age abattue
//
end
/////////////////
[s1,s2]=max(Ftot.*Phitot(1:4).*profit(1:4))  
ageopt=s2; disp('age optimal pour l''abattage:');disp(ageopt)
///////////////
figure(3)
clf()
plot(1:4,norm(x0init,1)*Ftot.*Phitot(1:4),'r-*','thickness',4)
xtitle('nombre de tetes abattues par an en fonction de l''age d''abattage')
pause
//////////////////////////////////
// ETAPE 3 : troupeau initial avec abattage optimal
// 
T0(ageopt)=T0(ageopt)-Phitot(ageopt)*T0(ageopt);
T=diag(T0,-1)
F=zeros(n,n)
F(1,:)=F0;
A=F+T
x0=x0init;
Xt=x0/norm(x0,1);Nt=[norm(x0,1)]
Nit=50;
for i=1:Nit
    x0=A*x0
    N=norm(x0,1)
    Xt=[Xt,x0/norm(x0,1)]
    Nt=[Nt,norm(x0,1)]
end
//
disp('plus grande valeur propre valeur de A');disp(max(abs(spec(A))))
disp('population finale après 30 ans (avec abattage)'); disp(norm(x0,1))
figure(4)
clf()
plot(0:Nit,Xt,'thickness',5)
xtitle('repartition par classe')
