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