#AMPL mod file for ECMOD 13-251. Optimal harvesting of an age-structured, two sex herbivore-plant system
#August 14. 2013
#For Knitro 7.0.0


#parameters

param maxt > 0;         			#time horizon
param r;                			#interest rate
param b=1/(1+r);				#discount factor
param p;      	         			#producer price of meat
param fn;               			#number of female age classes
param mn;               			#number of male age classes
param wfs{s in 0..fn};  			#female autumn weights
param wms{s in 0..mn};  			#male autumn weights
param K;                			#scaling parameter
param g;                			#("gamma" the fraction of carcass weight on autumn weight
param xf0 {s in 1..fn}; 			#initial number of females                              	
param xm0 {s in 1..mn}; 			#initial number of males                              	
param pxf0 {s in 1..fn};			#initial number of females at the end of period t=-1                              	
param pxm0 {s in 1..mn};			#initial number of males at the end of period t-1
param z0;               			#initial number of lichen biomass
param A;                			#the area of winter lichen pastures
param f{s in 0..fn};    			#baseline number of calves per female
param fcalvw{s in 0..fn}; 			#baseline weight of calves
param u;                			#the fraction of female calves
param alpha;					#concavity parameter in the objective function	 
param G3;					#lichen wastage parameter
param mm;      					#parameter in the male mortality function
param fm{s in 0..mn};   			#the harem size in different age classes
param U;					#auxliary parameter in the objective function for producing initial trial solutions
param ck;					#unit cost per ha
param cep;					#unit cost per adult animal
param ctp;					#unit cost per harvested animal
param EDf {s in 1..fn}:=0.683*(wfs[s-1]^0.75);	#overwinter daily energy requirement for females
param EDm {s in 1..mn}:=0.683*(wms[s-1]^0.75);  #overwinter daily energy requirement for males
param Yck:=ck*(A+2000)/K;			#unit land area cost times total land area					


#variables


#harvested females	
	var hf {s in 0..fn, t in 0..maxt-1}>=0;

#harvested males
	var hm {s in 0..mn, t in 0..maxt-1}>=0;

#number of females per age class
	var xf {s in 0..fn, t in -1..maxt}>= 0;

#number of males per age class
	var xm {s in 0..mn, t in -1..maxt}>=0;

#total number of females
	var Xf {t in 0..maxt}=sum{s in 1..fn}xf[s,t];

#total number of males
	var Xm {t in 0..maxt}=sum{s in 1..mn}xm[s,t];

#total number of individuals
	var X {t in 0..maxt}=sum{s in 1..fn}xf[s,t]+sum{s in 1..mn}xm[s,t];

#lichen biomass ha^-1
	var z {t in 0..maxt}>=0;

#the energy intake from lichen (fraction of requirement, cf. equation 20)
	var Lie {t in 0..maxt}=1.3242-4.0292*(1+exp(1*(z[t]+1000)/495.3806))^-0.522;

#consumption of lichen per winter day, females (cf. eq. 27)
	var LCf{s in 1..fn, t in 0..maxt}=G3*EDf[s]*Lie[t]/10.8;	

#consumption of lichen per winter day, males (cf. eq. 27)
	var LCm{s in 1..mn, t in 0..maxt}=G3*EDm[s]*Lie[t]/10.8;

#total lichen consumption per winter ha^-1
	var TC{t in 0..maxt}=(sum{s in 1..fn}LCf[s,t]*xf[s,t]*181+sum{s in 1..mn}LCm[s,t]*xm[s,t]*181)/A;

#fraction of energy from dwarf shrubs etc (cf. eq. 21)
	var Sie {t in 0..maxt}=0.4914-0.397*(1-exp(-0.0011*z[t]))^0.8646;

#total fraction of energy from requirement
	var IE{t in 0..maxt}=Lie[t]+Sie[t]; 

#the overwinter weight decrease (cf. eq. 22)
	var WS{t in 0..maxt}=0.488*(1+exp((IE[t]-0.6493)/0.133))^-1; 

#mortality adult females (cf. eq. 24)
	var MF{t in 0..maxt}=(1+exp(((0.36-WS[t])/0.011)))^-0.25;

#mortality, males, first winter females
	var MM{t in 0..maxt}=(1+exp(((0.36-mm*WS[t])/0.011)))^-0.25; 

#the decrease in the number of calves (fraction, cf. eq. 23)
	var Edf{t in 0..maxt}=1.2272-1.2272*(1+exp((0.2715-WS[t])/0.0239))^-0.1488;

#the dependence of calf weight on female weight decrease (cf. eqs. 28, 29) 
	var DFCAL{t in 0..maxt}=1.0275*(1+exp((WS[t]-0.3146)/0.0876))^-1; 

#the harmonic mean mating variable
	var Mdf{t in -1..maxt}>=0;

#auxliary variable for eq. (1)
	var df{t in 0..maxt}=Mdf[t-1]*Edf[t];

#number of calves born (cf. eq 1)
	var V {t in 0..maxt-1}=sum{s in 2..fn} f[s]*(1-MF[t])*xf[s,t]*df[t];

#the weight of female
	var FCALW{s in 1..fn, t in 0..maxt}=fcalvw[s]*DFCAL[t];

#the weight of male calves
	var MCALW{s in 1..fn, t in 0..maxt}=1.08*fcalvw[s]*DFCAL[t];

#average weight of female calves
	var MFCALVW{t in 0..maxt-1}=DFCAL[t]*(sum{s in 2..fn} f[s]*(1-MF[t])*xf[s,t]*df[t]*fcalvw[s])/V[t];

#average weight of male calves
	var MMCALVW{t in 0..maxt-1}=DFCAL[t]*(sum{s in 2..fn} f[s]*(1-MF[t])*xf[s,t]*df[t]*1.08*fcalvw[s])/V[t];

#annual meat gross revenues
	var YT {t in 0..maxt-1}=(g*8*p*MFCALVW[t]*hf[0,t]+g*8*p*MMCALVW[t]*hm[0,t]+g*(p*sum{s in 1..fn}wfs[s]*hf[s,t]+p*sum{s in 1..mn}wms[s]*hm[s,t]))/K;

#annual harvesting cost
	var Yctp {t in 0..maxt-1}=ctp*((sum{s in 0..fn}hf[s,t])+(sum{s in 0..mn}hm[s,t]))/K;

#annual management cost
	var Ycep {t in 0..maxt-1}=cep*X[t]/K;

#annual net revenues
	var Y {t in 0..maxt-1}=YT[t]-(Yctp[t])-Ycep[t]-Yck;




maximize objective:

sum{t in 0..maxt-1} ((1/(1+r))^t*(Y[t]+U)^alpha);    #use U=6 to obtain an inital trial solution



subject to restriction1 {t in 0..maxt-1}:

      xf[0,t]=V[t]*u;

subject to restriction2 {t in 0..maxt-1}:

      xf[1,t+1]=0.98*V[t]*u-hf[0,t];        

subject to restriction3 {t in 0..maxt-1}:

      xf[2,t+1]=(1-MM[t])*xf[1,t]-hf[1,t];  

subject to restriction4 {s in 2..fn-1, t in 0..maxt-1}:

      xf[s+1,t+1]=(1-MF[t])*xf[s,t]-hf[s,t];                                       				

subject to restriction5 {t in 0..maxt-1}:
      
      xm[0,t]=V[t]*(1-u);

subject to restriction6 {t in 0..maxt-1}:

      xm[1,t+1]=0.98*V[t]*(1-u)-hm[0,t];

subject to restriction7 {t in 0..maxt-1}:
      
      xm[2,t+1]=(1-MM[t])*xm[1,t]-hm[1,t]; 

subject to restriction8 {s in 2..mn-1, t in 0..maxt-1}:
      
      xm[s+1,t+1]=(1-MM[t])*xm[s,t]-hm[s,t];           

subject to restriction9:
        
        Mdf[-1]<=(2*sum{s in 1..mn}pxm0[s]*fm[s])/((sum{s in 1..mn}pxm0[s]*fm[s])+(sum{s in 2..fn}pxf0[s]));

subject to restriction10 {t in 0..maxt}:
        
         Mdf[t]<=(2*sum{s in 1..mn}(1-MM[t])*xm[s,t]*fm[s])/((sum{s in 1..mn}(1-MM[t])*xm[s,t]*fm[s])+(sum{s in 2..fn}(1-MF[t])*xf[s,t]));

subject to restriction11 {t in -1..maxt}:

      Mdf[t]<=1;

subject to restriction12 {t in 0..maxt-1}:

	z[t+1]=z[t]-TC[t]+((z[t]-TC[t])*(-0.7008)+(z[t]-TC[t])*(1-(z[t]-TC[t])/(-100.5832))^-0.0853); 
      
subject to restriction13 {s in 1..fn}: 

      xf[s,-1] = pxf0[s];

subject to restriction14{ s in 1..mn}: 

      xm[s,-1] = pxm0[s];

subject to restriction15 {t in 0..maxt-1}:
	
      0.98*V[t]*u-hf[0,t]>=0;

subject to restriction16 {t in 0..maxt-1}:
 
      (1-MM[t])*xf[1,t]-hf[1,t]>=0;

subject to restriction17 {s in 2..fn, t in 0..maxt-1}:
	
      (1-MF[t])*xf[s,t]-hf[s,t]>=0;   

subject to restriction18 {t in 0..maxt-1}: 

      0.98*V[t]*(1-u)-hm[0,t]>=0; 

subject to restriction19 {s in 1..mn, t in 0..maxt-1}:
	
      (1-MM[t])*xm[s,t]-hm[s,t]>=0;

subject to restriction20 {s in 1..fn}: 

      xf[s,0] = xf0[s];

subject to restriction21{ s in 1..mn}: 

      xm[s,0] = xm0[s];

subject to restriction22:

      z[0]=z0;

#auxliary restrictions for controlling the inital trial solutions

subject to restriction23 {t in 0..maxt-1}:

      z[t]>=200;

subject to restrriction24 {t in 0..maxt-1}:

      z[t]<=7000;

subject to restriction25 {s in 0..fn, t in 0..maxt-1}:

      xf[s,t]<=300;

subject to restriction26 {s in 0..mn, t in 0..maxt-1}:

      xm[s,t]<=300;

subject to restriction27 {s in 0..fn, t in 0..maxt-1}:

      hf[s,t]<=300;

subject to restriction28 {s in 0..mn, t in 0..maxt-1}:

      hm[s,t]<=300;