#AMPL mod file for ECOMOD-14-746. #Reindeer management and winter pastures in the presence of supplementary feeding and government subsidies #May 22, 2015 #For Knitro 7.0.0 and 9.1.0 #Antti-Juhani Pekkarinen, Olli Tahvonen, Jouko Kumpula #parameters param maxt > 0; #time horizon param r; #interest rate 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 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 lichen biomass param f{s in 0..fn}; #baseline number of calves per female param wc{s in 0..fn}; #baseline weight of calves param u; #the fraction of female calves param alpha; #concavity parameter in the objective function param sigmam; #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 CL; #unit cost per ha param Cx; #unit cost per adult animal param Cs; #unit cost per harvested animal param CV; #price of supplementary food param EDf {s in 0..fn}:=0.683*(wfs[s]^0.75); #daily energy requirement for females from spring to autumn param EDm {s in 0..mn}:=0.683*(wms[s]^0.75); #daily energy requirement for males from spring to autumn param q; #availability (kg/ha) of arboreal lichens param IQpure:=q*0.03*10.8/17.6; #Intake rate of pure arboreal lichen diet in winter period wII param IV:=0.57; #Energy intake rate for supplementary food param supplementaryfeeding; #1=supplementary feeding possible, 0=supplementary feeding not possible # Parameters for type of the pasture rotation param dwi; # number of the winter days when reindeer are in lichen pastures param dsp; # number of the spring days when reindeer are in lichen pastures (0=with pasture rotation, 31=without pasture rotation) param dsu; # number of the summer days when reindeer are in lichen pastures (0=with pasture rotation, 92=without pasture rotation) param dau; # number of the autumn days when reindeer are in lichen pastures (0=with pasture rotation, 61=without pasture rotation) param wwi; # total consumption of lichen in winter relative to that consumed for energy param wsp; # total consumption of lichen in spring relative to that consumed for energy param wsu; # total consumption of lichen in summer relative to that consumed for energy param wau; # total consumption of lichen in autumn relative to that consumed for energy param tausp; # the use of lichen for energy in spring relative to the use in winter with same biomass param tausu; # the use of lichen for energy in summer relative to the use in winter with same biomass param tauau; # the use of lichen for energy in autumn relative to the use in winter with same biomass # Parameters for pasture classes and land areas param Aomf; #area of old or mature coniferous forest param Aybf; #area of young pine forest, logging area or mountain birch forest param Amoh; #area of mountain heath param A:=Aomf+Aybf+Amoh; #total area of winter lichen pastures param Atotal:=A*3; #total land area of the reindeer herding district param YCL:=CL*(Atotal)/K; #unit land area cost times total land area param AQ; #the land area (ha) of arboreal lichen pastures param g:=(Aomf+Aybf*0.6+Amoh*0.4)/A; #Parameter g in eq.(23) ##### variables ##### #harvested females var hf {s in 0..fn, t in 0..maxt}>=0; #harvested males var hm {s in 0..mn, t in 0..maxt}>=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; # Length of the winter period b, eq.(1) see also var SQ[t], restriction28 and restriction29 var db{t in 0..maxt}>=0; # Supplementary feeding per day in winter period a (kg/day) var va {t in 0..maxt}>=0; # Supplementary feeding per day in winter period b (kg/day) var vb {t in 0..maxt}>=0; #### ENERGY INTAKE IN LATE WINTER WITH ARBOREAL LICHENS (b), energy intake model #### # Intake rate of pure cratered food diet in winter period b var IZpureb {t in 0..maxt}=(0.055+0.00005*z[t])*2.4/3; # The share of foraging time of natural food used for cratered food in winter period b var TZb{t in 0..maxt}=((IZpureb[t])^5)/(((IZpureb[t])^5)+((IQpure)^5)); # Intake rate of natural food in winter period wII var IZQb {t in 0..maxt}=TZb[t]*IZpureb[t]*(TZb[t])^-0.2+(1-TZb[t])*IQpure*(1-TZb[t])^-0.2; # The maximum share of foraging time used for supplementary food in winter period wII var TVbmax{t in 0..maxt}=0.3*(0.2+0.9/(1+exp((IZQb[t]-0.1)/0.05))); # The share of foraging time used for supplementary food in winter period wII var TVb{t in 0..maxt}=0.1*vb[t]*17.6/(sum{s in 1..fn}xf[s,t]*EDf[s]+sum{s in 1..mn}xm[s,t]*EDm[s]); # Total average intake rate in winter period wII var Ib{t in 0..maxt}=IZQb[t]*(1-TVb[t])+IV*TVb[t]; # Total daily foraging time in winter period b var Fb {t in 0..maxt}=1.8508+8.1492/((1+exp((Ib[t]-0.0953)/0.0013))^0.0066); # Total energy intake relative to energy requirement in winter period wII var ETb{t in 0..maxt}=Ib[t]*Fb[t]; # Energy intake from lichen relative to energy requirement in winter period wII var ELb {t in 0..maxt}=(0.3621*(1-exp(-0.0048985*z[t]))+0.5603*(1-exp(-0.0015299*z[t])))*Fb[t]*IZpureb[t]*(TZb[t]^-0.2)*(1-TVb[t])*TZb[t]; #### LENGTHS OF WINTER PERIODS #### # Suffisiency of arboreal lichens (for days per winter) var SQ {t in 0..maxt}=AQ/(X[t]*0.03*Fb[t]); # Number of winter days without arboreal lichens, eq(2) var da {t in 0..maxt}=181-db[t]; #### ENERGY INTAKE IN EARLY WINTER(wI), energy intake model #### # Intake rate of cratered food in winter period wI var IZa {t in 0..maxt}=(0.055+0.00005*z[t])*((3.3*121+2.4*(60-db[t]))/(121+(60-db[t])))/3; # The maximum share of foraging time used for supplementary food in winter period wI var TVamax{t in 0..maxt}=0.3*(0.2+0.9/(1+exp((IZa[t]-0.1)/0.05))); # The share of foraging time used for supplementary food in winter period wI var TVa{t in 0..maxt}=0.1*va[t]*17.6/(sum{s in 1..fn}xf[s,t]*EDf[s]+sum{s in 1..mn}xm[s,t]*EDm[s]); # Total average intake rate in winter period a var Ia{t in 0..maxt}=IZa[t]*(1-TVa[t])+IV*TVa[t]; # Total daily eating time in winter period a var Fa {t in 0..maxt}=1.8508+8.1492/((1+exp((Ia[t]-0.0953)/0.0013))^0.0066); # Total energy intake relative to energy requirement in winter period a var ETa{t in 0..maxt}=Ia[t]*Fa[t]; # Energy intake from lichen relative to energy requirement in winter period wI var ELa {t in 0..maxt}=(0.3621*(1-exp(-0.0048985*z[t]))+0.5603*(1-exp(-0.0015299*z[t])))*Fa[t]*IZa[t]*(1-TVa[t]); #### AVERAGE ENERGY INTAKE IN WINTER AND THE ITS CONSEQUENCES, population model #### # Average total energy intake in winter var ET {t in 0..maxt}=(da[t]*ETa[t]+db[t]*ETb[t])/181; # Average energy intake from lichen in winter var ELwi {t in 0..maxt}=(da[t]*ELa[t]+db[t]*ELb[t])/181; #the overwinter weight decrease, eq.(5) var wd{t in 0..maxt}=0.5*exp(-exp((ET[t]-0.72)/0.22)); #mortality adult females, eq.(T6) var mf{t in 0..maxt}=(1+exp(((0.36-wd[t])/0.011)))^-0.25; #mortality, males, first winter females, eq.(T6) var mm{t in 0..maxt}=(1+exp(((0.36-sigmam*wd[t])/0.011)))^-0.25; #the decrease in the number of calves, eq.(T5) var fwd{t in 0..maxt}=1.2272-1.2272*(1+exp((0.2715-wd[t])/0.0239))^-0.1488; #the dependence of calf weight on female weight decrease, see eq.(T8) var wcwd{t in 0..maxt}=1.0275*(1+exp((wd[t]-0.3146)/0.0876))^-1; #the harmonic mean mating variable, see restriction8, restriction9 and restriction10 var beta{t in -1..maxt}>=0; #auxliary variable for eq.(T1) var fwdbeta{t in 0..maxt}=beta[t-1]*fwd[t]; #number of calves born, eq.(T1) var X0 {t in 0..maxt}=sum{s in 2..fn} f[s]*(1-mf[t])*xf[s,t]*fwdbeta[t]; #the weight of female calves, eq.(T8) var wcf{s in 1..fn, t in 0..maxt}=wc[s]*wcwd[t]; #the weight of male calves, eq.(T8) var wcm{s in 1..fn, t in 0..maxt}=1.08*wc[s]*wcwd[t]; #average weight of female calves var w0f{t in 0..maxt}=(sum{s in 2..fn} f[s]*(1-mf[t])*xf[s,t]*fwdbeta[t]*wcf[s,t])/(X0[t]); #average weight of male calves var w0m{t in 0..maxt}=(sum{s in 2..fn} f[s]*(1-mf[t])*xf[s,t]*fwdbeta[t]*wcm[s,t])/(X0[t]); #### REVENUES AND COSTS, economic model #### #annual meat gross revenues, eq.(36) var R {t in 0..maxt}=(gamma*8*p*w0f[t]*hf[0,t]+gamma*8*p*w0m[t]*hm[0,t]+gamma*(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 YCs {t in 0..maxt}=Cs*((sum{s in 0..fn}hf[s,t])+(sum{s in 0..mn}hm[s,t]))/K; #annual management cost var YCx {t in 0..maxt}=Cx*X[t]/K; # Annual feeding cost var YCVa {t in 0..maxt}=CV*va[t]*da[t]/K; var YCVb {t in 0..maxt}=CV*vb[t]*db[t]/K; var YCV {t in 0..maxt}=YCVa[t]+YCVb[t]; #annual total costs, eq.(37) var YC {t in 0..maxt}=YCs[t]+YCx[t]+YCL+YCV[t]; #annual net revenues (scaled) var Y {t in 0..maxt}=R[t]-YC[t]; #annual net revenues in euros per A var Y2 {t in 0..maxt}=Y[t]*K; #### PASTURE ROTATION, lichen model##### # Lichen consumption in winter # # Lichen consumption per female per winter day var lfwi{s in 1..fn, t in 0..maxt}=wwi*EDf[s-1]*ELwi[t]/10.8; # Lichen consumption per male per winter day var lmwi{s in 1..mn, t in 0..maxt}=wwi*EDm[s-1]*ELwi[t]/10.8; # Total lichen consumption per winter per hectare, eq.(31) var lwi{t in 0..maxt}=dwi*(sum{s in 1..fn}lfwi[s,t]*xf[s,t]+sum{s in 1..mn}lmwi[s,t]*xm[s,t])/A; # Lichen consumption in spring # # Lichen biomass in the beging of the spring period, eq.(19) var zsp {t in 0..maxt}=z[t]-lwi[t]; # Energy intake from lichen in spring relative to energy requirement, eq.(28) var ELsp {t in 0..maxt}=tausp*(1.3242-4.0292/(1+exp(-1*(zsp[t]+1000)/-495.3806))^0.522); # Lichen consumption per female per spring day, eq.(24) var lfsp{s in 1..fn, t in 0..maxt}=wsp*EDf[s]*ELsp[t]/10.8; # Lichen consumption per male per spring day, eq.(24) var lmsp{s in 1..mn, t in 0..maxt}=wsp*EDm[s]*ELsp[t]/10.8; # Total lichen consumption per spring per hectare, eq.(32) var lsp{t in 0..maxt}=dsp*(lfsp[1,t]*(1-mm[t])*xf[1,t]+sum{s in 2..fn}lfsp[s,t]*(1-mf[t])*xf[s,t]+sum{s in 1..mn}lmsp[s,t]*(1-mm[t])*xm[s,t])/A; # Lichen consumption in summer # # Lichen biomass in the beging of the summer period, eq.(20) var zsu {t in 0..maxt}=zsp[t]-lsp[t]; # Energy intake from lichen in summer relative to energy requirement, eq.(28) var ELsu {t in 0..maxt}=tausu*(1.3242-4.0292/(1+exp(-1*(zsu[t]+1000)/-495.3806))^0.522); # Lichen consumption per female per summer day, eq.(24) var lfsu{s in 1..fn, t in 0..maxt}=wsu*EDf[s]*ELsu[t]/10.8; # Lichen consumption per male per summer day, eq.(24) var lmsu{s in 1..mn, t in 0..maxt}=wsu*EDm[s]*ELsu[t]/10.8; # Total lichen consumption per summer per hectare, eq.(33) var lsu{t in 0..maxt}=dsu*(lfsu[1,t]*(1-mm[t])*xf[1,t]+sum{s in 2..fn}lfsu[s,t]*(1-mf[t])*xf[s,t]+sum{s in 1..mn}lmsu[s,t]*(1-mm[t])*xm[s,t])/A; # Growht of lichen in summer # # Lichen growth, eq.(23) var G {t in 0..maxt}=g*zsu[t]*(-0.7008+(1+zsu[t]/100.5832)^-0.0853); # Lichen consumption in autumn # # Lichen biomass in the beging of the autumn period, eq.(21) var zau {t in 0..maxt}=zsu[t]-lsu[t]+G[t]; # Energy intake from lichen in autumn relative to energy requirement, eq.(28) var ELau {t in 0..maxt}=tauau*(1.3242-4.0292/(1+exp(-1*(zau[t]+1000)/-495.3806))^0.522); # Lichen consumption per female per autumn day, eq.(24) var lfau{s in 0..fn, t in 0..maxt}=wau*EDf[s]*ELau[t]/10.8; # Lichen consumption per male per autumn day, eq.(24) var lmau{s in 0..mn, t in 0..maxt}=wau*EDm[s]*ELau[t]/10.8; # Total lichen consumption per autumn per hectare, eq.(34) var lau{t in 0..maxt}=dau*(lfau[0,t]*X0[t]*u*0.98+lmau[0,t]*X0[t]*(1-u)*0.98+lfau[1,t]*(1-mm[t])*xf[1,t]+sum{s in 2..fn}lfau[s,t]*(1-mf[t])*xf[s,t]+sum{s in 1..mn}lmau[s,t]*(1-mm[t])*xm[s,t])/A; #### VARIABLES FOR DISPLAYING THE RESULTS FOR FIGURE 9a #### var LichenBiomass {t in 0..maxt}=z[t]; var NumberOfReindeer {t in 0..maxt}=X[t]; var SlaugtheredReindeer {t in 0..maxt}=sum{s in 0..fn}hf[s,t]+sum{s in 0..mn}hm[s,t]; var SlaugtheredCalves {t in 0..maxt}=hf[0,t]+hm[0,t]; var NetRevenues {t in 0..maxt}=Y2[t]; var CalvesPer100Females {t in 0..maxt}=100*X0[t]/Xf[t]; var AverageBirthWeight {t in 0..maxt}=(xf[0,t]*w0f[t]+xm[0,t]*w0m[t])/X0[t]; var SupplemFoodPerReindeer {t in 0..maxt}=(vb[t]*db[t]+va[t]*da[t])/X[t]; ########################################################################################################################## #RESTRICTIONS ########################################################################################################################## #### Objective funtion, eq.(35) maximize objective: sum{t in 0..maxt} ((1/(1+r))^t*(Y[t]+U)^alpha); #use U=6 to obtain an inital trial solution if: 0=0; subject to restriction12 {t in 0..maxt}: (1-mm[t])*xf[1,t]-hf[1,t]>=0; subject to restriction13 {s in 2..fn, t in 0..maxt}: (1-mf[t])*xf[s,t]-hf[s,t]>=0; subject to restriction14 {t in 0..maxt}: 0.98*X0[t]*(1-u)-hm[0,t]>=0; subject to restriction15 {s in 1..mn, t in 0..maxt}: (1-mm[t])*xm[s,t]-hm[s,t]>=0; #### Restriction for lichen model subject to restriction16 {t in 0..maxt-1}: # eq.(22) z[t+1]=zau[t]-lau[t]; #### Initial levels for the state variables subject to restriction17 {s in 1..fn}: #eq.(43) xf[s,-1] = pxf0[s]; subject to restriction18{ s in 1..mn}: #eq.(43) xm[s,-1] = pxm0[s]; subject to restriction19 {s in 1..fn}: #eq.(42) xf[s,0] = xf0[s]; subject to restriction20{ s in 1..mn}: #eq.(42) xm[s,0] = xm0[s]; subject to restriction21: #eq.(44) z[0]=z0; ##### auxliary restrictions for controlling the inital trial solutions subject to restriction22 {t in 0..maxt}: z[t]>=200; subject to restrriction23 {t in 0..maxt}: z[t]<=5000; subject to restriction24 {s in 1..fn, t in 0..maxt}: xf[s,t]<=500; subject to restriction25 {s in 1..mn, t in 0..maxt}: xm[s,t]<=500; subject to restriction26 {s in 1..fn, t in 0..maxt}: hf[s,t]<=500; subject to restriction27 {s in 1..mn, t in 0..maxt}: hm[s,t]<=500; subject to restriction27b {t in 0..maxt}: Y[t]>=0; ##### Length of winter period b subject to restriction28 {t in 0..maxt}: db[t]<=60; subject to restriction29 {t in 0..maxt}: db[t]<=SQ[t]; ##### Restrictions for the share of foraging time used for supplementary food subject to restriction30 {t in 0..maxt}: TVa[t]<=TVamax[t]*supplementaryfeeding; subject to restriction31 {t in 0..maxt}: TVb[t]<=TVbmax[t]*supplementaryfeeding;