### Rämö Janne ### # Rämö, J. & Tahvonen, O. # # OPTIMIZING THE HARVEST TIMING IN CONTINUOUS COVER FORESTRY # # DOI: 10.1007/s10640-016-0008-4 # ### Amount of species, 1=Spruce, 2=Spruce+Birch; 3=Spruce+Birch+Pine; 4=Spruce+Birch+Pine+Others ### ### Set in .dat ### param sp; ## Amount of size classes, default 12, max 15 ### param sc; ### 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 := 1..sp; set SIs ordered; param g {t in T} binary default 0; param C; param BESTobj default 0; param SI; #Site index param pi = 4 * atan(1); #Set the value of pi to parameter param LAT; #Latitude param w; #Width of the size classes param y0 {i in 1..15, j in 1..4}; #Initial stand, number of trees param DBH {i in 1..15}; #Diameter at breast height, mm var y {t in T, i in 1..sc, j in pl} >= 0; #Stand development parameter var BAsc {t in T, i in 1..sc, j in pl} = pi*(DBH[i]/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..sc} BAsc[t,e,i]; #Basal area of species var BA {t in T} = sum {e in 1..sc, u in pl} BAsc[t,e,u]; #Basal area, m^2/ha var BAL {t in T, i in 1..sc} = sum {e in (i+1)..sc, 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 param s; #Hakkuusyklin pituus, 1=5v, 2=10v jne. param r; #Korkokanta param R := 1/(1+r); #Diskonttotekijä param ps {j in 1..4}; #Tukin hinta, tienvarressa #param ps {i in 1..sc, j in 1..4}; param pp {j in 1..4}; #Kuidun hinta, tienvarressa param pkerroin {j in 1..4}; #Koivun hintakerroin koodin debuggausta varten, vektori muun koodin rakenteen takia param Ikerroin {j in 1..4}; #Growth multiplier for sensitivity analysis. param dkerroin {j in 1..4}; #Ingrowth multiplier for sensitivity analysis. ### Volumes, Vs=Sawn timber, Vp=Pulp ### param Vs {e in SIs, i in 1..15, j in 1..4}; param Vp {e in SIs,i in 1..15, j in 1..4}; ### Stand volume ### var standVol {t in T, j in pl} = sum {i in 1..sc} (y[t,i,j]*Vp[SI,i,j]+y[t,i,j]*Vs[SI,i,j]); ### Mortality-function ### var m {t in T, i in 1..sc, j in pl} = (1+exp(-(mp[1,j]+mp[2,j]*DBH[i]+mp[3,j]*10^-5*DBH[i]^2+mp[4,j]*BA[t])))^-1; ### Diameter increment -function ### var I {t in T, i in 1..sc, j in pl} = Ikerroin[j]*(di[1,j]+di[2,j]*DBH[i]+di[3,j]*10^-5*DBH[i]^2+di[4,j]*10^-8*DBH[i]^3+di[5,j]*BAL[t,i]+di[6,j]*SI+di[7,j]*BA[t]+di[8,j]*LAT); ### Size class variation functions ### var b {t in T, i in 1..sc, j in pl} = I[t,i,j]/w; #Possibility that species grows to next size class var a {t in T, i in 1..sc, j in pl} = 1-b[t,i,j]-m[t,i,j]; #Possibility that tree stays in the same size class ### Recruitment-functions ### 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]*dkerroin[i]; #Amount recruited ### Harvests, amount of trees ### var H {t in T, i in 1..sc, j in pl} >= 0; param Hbool {t in T} default 0; ### Total cubic meters of wood harvested per size class ### var Vtotal {t in T, j in pl} = sum {i in 1..sc} H[t,i,j]*(Vs[SI,i,j]+Vp[SI,i,j]); ### Total number of trees harvested per size class ### var Htotal {t in T, j in pl} = sum {i in 1..sc} H[t,i,j]; ### Total amount of trees per size class ### var ytotal {t in T, j in pl} = sum {i in 1..sc} y[t,i,j]; #var d {t in T} = exp(6.154-1.683*log(BA[t])+0.0642*sqrt(ytotal[t,3]))-1; ### Pukkala et al. 2009 ### Hakkuumäärät, m3, Hs = sawn timer, Hp = Pulp ### var Hs {t in T, i in 1..sc, j in pl} = H[t,i,j]*Vs[SI,i,j]; var Hp {t in T, i in 1..sc, j in pl} = H[t,i,j]*Vp[SI,i,j]; # Cubic meters per size class # param Vsc {i in 1..sc, j in pl, k in SIs} = (Vs[SI,i,j]+Vp[SI,i,j]); ### Harvesting costs ### var cuttingCosts {t in T} = sum {i in 1..sc, j in pl} H[t,i,j]*(0.412+0.758*Vsc[i,j,SI]-0.180*Vsc[i,j,SI]^2)*1.15*2.1; var haulingCosts {t in T} = sum {j in pl} Vtotal[t,j]*(2.272+500/(33.7*Vtotal[t,j]+0.0001)+0.211/(0.39459*(Vtotal[t,j]+0.0001)^0.3))*1; var varCostsTotal {t in T} = cuttingCosts[t]+haulingCosts[t]; # var costsperscNL {t in T, i in 1..sc} = ((0.44472+0.94*sum {u in pl} Vsc[i,u,SI]-(146.17/(sum {u in pl} Vsc[i,u,SI]*1000+862.05)))*2.1001366*H[t,i,j])*Hbool[t]; param fixedCosts {t in T} = g[t]*C; var Revenues {t in T} = sum {i in 1..sc, j in pl} (Hs[t,i,j]*ps[j]+Hp[t,i,j]*pp[j])*pkerroin[j]; ### Assisting transition-matrix, number of trees that move to next size class ### var trans {t in T, i in 1..sc, j in pl} = b[t,i,j]*y[t,i,j]; ### Assisting stand-state matrix, trees that say in the same size class ### var stay {t in T, i in 1..sc, j in pl} = a[t,i,j]*y[t,i,j]; ### Revenues ### var c {t in T} = Revenues[t]-fixedCosts[t]-varCostsTotal[t]; ### Set the initial state ### subject to initial_state {i in 1..sc, j in pl}: y[0,i,j] = y0[i,j]; ### Stand matrices ### subject to standstate_1sc {t in 0..maxt-5 by 5, j in pl}: y[t+5,1,j] = d[t,j]+stay[t,1,j]-H[t,1,j]; subject to standstate {t in 0..maxt-5 by 5, i in 1..sc-2, j in pl}: y[t+5,i+1,j] = trans[t,i,j]+stay[t,i+1,j]-H[t,i+1,j]; ### Trees from biggest size class cannot grow to next size class, so it gets its own stand-constraint to prevent trees from disappearing ### subject to standstate_12sc {t in 0..maxt-5 by 5, j in pl}: y[t+5,sc,j] = trans[t,sc-1,j]+y[t,sc,j]*(1-m[t,sc,j])-H[t,sc,j]; ### Setting the boolean-matrix's effect on harvests ### subject to harvests3 {t in T, i in 1..sc, j in pl}: H[t,i,j] = g[t]*H[t,i,j]; subject to yrajoite {t in T, i in 1..sc, j in pl}: y[t,i,j] <= 1750; ### Objective ### maximize objective: sum {t in 0..maxt-5 by 5} c[t]*R^t/1000; # param i; # param k; param numberOfH := 5; param ch {1..numberOfH} default 0; param finished; param ho {1..numberOfH} default 0; param hso default 0; param hs {1..numberOfH} default 0; param ss default 0; param dir; param o default 0; param objCandidate default 0; param GBestObj default 0; param Ghs {1..numberOfH} default 0; param Gss default 0; param totalRuns; # param i; param filename symbolic := "UEA-300-3%-optimalRuns-2.txt";