### 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";