### Rämö Janne ### # Bollandsås et al. (2008) growth model, age structured # ### Species, set in .dat ### param sp; ## Amount of age classes, default 12 ### param ac; ### Recruit-function, probability ### param rp {1..4, 1..4}; ### Recruit-function parameters, conditional ### param rc {1..4, 1..4}; ### Diameter increment -function parameters ### param di {1..8, 1..4}; ### Mortality-function parameters ### param mp {1..4, 1..4}; param maxt; #Max time set T := 0..maxt by 5; set pl := sp..sp; param q := 0.5; param r; param R := 1/(1+r); param s; param SI; #Site index param pi = 4 * atan(1); #Set the value of pi to parameter param LAT; #Latitude param y0 {i in 1..ac, j in 1..4} default 0; #Initial stand, number of trees var DBH {t in T, i in 1..ac, j in pl} >= 0; param DBH0 {i in 1..ac} default 0; param ps {j in 1..4}; #Tukin hinta, kantohinta param pp {j in 1..4}; #Kuidun hinta, kantohinta var y {t in T, i in 1..ac, j in pl} >= 0; #Stand development parameter var H {t in T, i in 1..ac, j in pl} >= 0; param Hbool {t in T} default 0; var BAsc {t in T, i in 1..ac, j in pl} = pi*(DBH[t,i,j]/2/1000)^2*y[t,i,j]; #Basal area of individual size classes, array, calculated from y & DBH var BAsp {t in T, i in pl} = sum {e in 1..ac} BAsc[t,e,i]; #Basal area of species var BA {t in T} = sum {u in pl} BAsp[t,u]; #Basal area, m^2/ha var BAL {t in T, i in 1..ac} = sum {e in (i+1)..ac, u in pl} BAsc[t,e,u]; #Basal area of larger trees, m^2/ha var PBA {t in T, i in pl} = BAsp[t,i]/BA[t]*100; #Percentage of Basal area for the subject species #Uncomment the species-sitetype combination that is optimized, estimated from the data used in size-structured model. # Spruce volume functions, SI=15 var Vs {t in T, i in 1..ac, j in pl} = (116.0906-31.1854*(DBH[t,i,j]/10)+1.9407*(DBH[t,i,j]/10)^2-0.0121*(DBH[t,i,j]/10)^3)/1000; var Vp {t in T, i in 1..ac, j in pl} = (0.0068176*(DBH[t,i,j]/10)^3-0.660699*(DBH[t,i,j]/10)^2+18.2853*(DBH[t,i,j]/10)-72.8905)/1000; # Pine volume functions, SI=15 # var Vs {t in T, i in 1..ac} = ((-32.777+3623.353/(1+(DBH[t,i,3]/10/48.547)^-3.256))/1000; # var Vp {t in T, i in 1..ac} = (24.954+110.575/(1+((DBH[t,i,3]/10-15.797)/12.562)^2))/1000; # Pine volume functions, SI=11 # var Vs {t in T, i in 1..ac} = (-31.764+3557.296/(1+(DBH[t,i,j]/10/49.073)^-3.233))/1000; # var Vp {t in T, i in 1..ac} = (24.6+107.658/(1+((DBH[t,i,j]/10-15.771)/2.470)^2))/1000; var standVol {t in T, j in pl} = sum {i in 1..ac} (y[t,i,j]*Vp[t,i,j]+y[t,i,j]*Vs[t,i,j]); var Hs {t in T, i in 1..ac, j in pl} = H[t,i,j]*Vs[t+5,i+1,j]; #Harvests occur at the end of the period var Hp {t in T, i in 1..ac, j in pl} = H[t,i,j]*Vp[t+5,i+1,j]; var m {t in T, i in 1..ac, j in pl} = (1+exp(-(mp[1,j]+mp[2,j]*DBH[t,i,j]+mp[3,j]*10^-5*DBH[t,i,j]^2+mp[4,j]*BA[t])))^-1; var I {t in T, i in 1..ac, j in pl} = di[1,j]+di[2,j]*DBH[t,i,j]+di[3,j]*10^-5*DBH[t,i,j]^2+di[4,j]*10^-8*DBH[t,i,j]^3+di[5,j]*BAL[t,i]+di[6,j]*SI+di[7,j]*BA[t]+di[8,j]*LAT; # ) <= 0 then 0 else (di[1,j]+di[2,j]*DBH[t,i,j]+di[3,j]*10^-5*DBH[t,i,j]^2+di[4,j]*10^-8*DBH[t,i,j]^3+di[5,j]*BAL[t,i]+di[6,j]*SI+di[7,j]*BA[t]+di[8,j]*LAT); var CR {t in T, i in pl} = rc[1,i]*BA[t]^rc[2,i]*SI^rc[3,i]*PBA[t,i]^rc[4,i]; #Recruitment, combability var Rprob {t in T, i in pl} = (1+exp(-(rp[1,i]+rp[2,i]*BA[t]+rp[3,i]*SI+rp[4,i]*PBA[t,i])))^-1; #Recruitment, probability var d {t in T, i in pl} = CR[t,i]*Rprob[t,i]; #Amount recruited ### Total cubic meters of wood harvested per size class ### var Vtotal {t in T, j in pl} = sum {i in 1..ac} H[t,i,j]*(Vs[t+5,i+1,j]+Vp[t+5,i+1,j]); #Harvests occur at the end of the period ### Total number of trees harvested per size class ### var Htotal {t in T, j in pl} = sum {i in 1..ac} H[t,i,j]; ### Total amount of trees per size class ### var ytotal {t in T, j in pl} = sum {i in 1..ac} y[t,i,j]; ### Revenues ### var c {t in T} = sum {i in 1..ac, j in pl} (Hs[t,i,j]*ps[j]+Hp[t,i,j]*pp[j]); ### Set the initial state ### subject to initial_state {i in 1..ac, j in pl}: y[0,i,j] = y0[i,j]; subject to initialDBH {i in 2..ac, j in pl}: DBH[0,i,j] = DBH0[i]; subject to diameter1 {t in T, j in pl}: DBH[t,1,j] = 50; subject to diameter {t in 0..maxt-5 by 5, i in 1..ac-1, j in pl}: DBH[t+5,i+1,j] = DBH[t,i,j]+I[t,i,j]; subject to harvests3 {t in T, i in 1..ac, j in pl}: H[t,i,j] = Hbool[t]*H[t,i,j]; subject to standstate1 {t in 0..maxt-5 by 5, j in pl}: y[t+5,1,j] = d[t,j]-H[t,1,j]; subject to standstate {t in 0..maxt-5 by 5, i in 1..ac-2, j in pl}: y[t+5,i+1,j] = (1-m[t,i,j])*y[t,i,j]-H[t,i+1,j]; subject to standstateF {t in 0..maxt-5 by 5, j in pl}: y[t+5,ac,j] = (1-m[t,ac-1,j])*y[t,ac-1,j]+(1-m[t,ac,j])*y[t,ac,j]-H[t,ac,j]; subject to ylimiter {t in T, i in 1..ac, j in pl}: y[t,i,j] <= 300; maximize objective: (sum {t in 0..maxt-5 by 5} (c[t]*R^t/1000)); param i;