### Rämö Janne ### # Bollandsås et al. (2008) growth model # ### Species, 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 := sp..sp; set SIs ordered; 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 y1 {t in T, i in 1..sc, j in pl} >= 0; #Stand before harvests 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; #Harvest interval param r; #Interest rate param R := 1/(1+r); #Discount factor param ps {j in 1..4}; #Saw log price param pp {j in 1..4}; #Pulp price # Multipliers for sensitivity analyses param pkerroin {j in 1..4}; #Price param Ikerroin {j in 1..4}; #Growth param dkerroin {j in 1..4}; #Ingrowth ### Volume-functions, 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}; ### 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; #Portion of trees that grow 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]; #Portion of trees that stay 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]; ### 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]; ### Revenues ### var c {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]); ### 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]+a[t,1,j]*y[t,1,j]-Hbool[t]*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] = b[t,i,j]*y[t,i,j]+a[t,i+1,j]*y[t,i+1,j]-Hbool[t]*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] = b[t,sc-1,j]*y[t,sc-1,j]+y[t,sc,j]*(1-m[t,sc,j])-Hbool[t]*H[t,sc,j]; # subject to yrajoite {t in T, i in 1..sc, j in pl}: # y[t,i,j] <= 400; ### Setting the boolean-matrix's effect on harvests ### subject to harvests {t in T, i in 1..sc, j in pl}: H[t,i,j] = Hbool[t]*H[t,i,j]; ### Maximizing the profits ### maximize objective: sum {t in 0..maxt-5 by 5} c[t]*R^t/1000; param i;