### Rämö Janne ###

# Rämö, J. & Tahvonen, O. #
# OPTIMIZING THE HARVEST TIMING IN CONTINUOUS COVER FORESTRY #
# DOI: 10.1007/s10640-016-0008-4 #

reset;
model norwaysp-hc.mod.txt; 
data norwaysp-hc.dat.txt;

option solver knitroampl;
option display_1col 10;
option randseed 0;
option knitro_options "maxit=8000 opttol=1.0e-5 ms_enable=1 ms_maxsolves=5 outlev=0 xtol=1.0e-20 par_numthreads=1";

let totalRuns := 5;
display totalRuns;
for runControl {q in 1..totalRuns} {
	#Randomly generate the starting point, Harvest 1: [0,50], harvest 2: (H1,H1+50], harvest 3: (H2,H2+50] etc.
	let ch[1]:=ceil(Uniform(0,10))*5;
	let ch[2]:=ch[1]+ceil(Uniform(0,10))*5;
	let ch[3]:=ch[2]+ceil(Uniform(0,10))*5;
	let ch[4]:=ch[3]+ceil(Uniform(0,10))*5;
	let ch[5]:=ch[4]+ceil(Uniform(0,10))*5;
	#Steady state interval
	let s:=ceil(Uniform(1,6));
	
	display ch,s,C;
	
	#Set the times of harvests for Knitro, H=g*H (in model file)
	let {k in 1..numberOfH} g[ch[k]] := 1;

	let {i in ch[numberOfH]+(5*s)..maxt by 5*s} g[i] := 1;

	solve;
	if solve_result_num <> 0 then {
		print "Run number ", q, "of", totalRuns > (filename);
		print "Locally optimal solution not found, continuing with next step" > (filename);
		continue runControl;
	}

	let objCandidate := objective;
	
	#Checking for optimality, if better, save in parameters, hs and ss include the best solution of this run.
	if objCandidate >= BESTobj then {
		let BESTobj := objCandidate;
		let {k in 1..numberOfH} hs[k]:=ch[k];
		
		let ss:=s;
	}

	#Set the optimality to zero (this is needed here for subsequent starts)
	let {k in 1..numberOfH} ho[k]:=0;
	

	let finished := 0;	#When Finished=1 we have found local optimum
	let dir:=1;	#Determines whether the harvests should adjust upwards or downwards
	
	repeat while finished < 1 {	#Repeat until local optimum is found
		for {k in 1..numberOfH} {	#Optimize the times for all harvests, numberOfH is set in the model file.
			let o:=0;
			repeat while ho[k] < 1 {	#When ho[k]=1 we have the best result for harvest number k, i.e. repeat this until objective cannot be improved anymore
				
				let {t in T} g[t] := 0;	#Reset harvesting strategy
				if ch[k]+dir*5 >= 0 then { #Times of harvests cannot be below 0
					let {j in 1..numberOfH} g[ch[j]] := 1;
					
					let g[ch[k]] := 0;	#Set the change in harvest k (in this and next row)
					let g[ch[k]+dir*5] := 1;

					let {i in ch[numberOfH]+(5*s)..maxt by 5*s} g[i] := 1;


					solve;
					if solve_result_num <> 0 then {
						print "Run number ", q, "of", totalRuns > (filename);
						print "Locally optimal solution not found, continuing with next step" > (filename);
						continue runControl;
					}
					let objCandidate := objective;
					
					for {j in 1..numberOfH} {	#Display so we can track where optimization is going
						if j = k then {
							display ch[j]+dir*5;
						}
						else display ch[j];
					}
					display s;
					display objCandidate-BESTobj;	
				} else {
					print "Harvest no.", k, " is zero:", ch[k]+dir*5, "< 0";
				}
				
				if objCandidate > BESTobj then {			#Check for optimality
					let BESTobj := objCandidate;			#Save for best result
					let ch[k]:=ch[k]+dir*5;					#If optimal, set the time of harvest k to new result
					let {j in 1..numberOfH} ho[j] := 0;		#Reset optimality of other harvests
					let {j in 1..numberOfH} hs[j] := ch[j]; #Save the times of harvests of the best result
					
					let o:=1;
					let hso:=0;
					let ss:=s;
				}
				else {
					if dir = 1 and o <> 1 then {	#If no steps have been taken, change direction
						let dir:=-1;
					}
					else {	#If already have taken 1 or more steps, or dir=-1, declare the result as optimal
						let dir:=1;
						let ho[k]:=1;
						print "Harvest", k ,"is optimal:";
						display ch[k], BESTobj;
					}
				}
				display ho, hso;
			}
		}
		let dir:=1;
		repeat while hso < 1 {	#The steady state harvesting interval optimization
			if s+dir > 0 then {
				let {t in T} g[t] := 0;	#Reset	
				let {j in 1..numberOfH} g[ch[j]] := 1;		#Set the times of harvests
				
				let {i in ch[numberOfH]+(5*(s+dir))..maxt by 5*(s+dir)} g[i] := 1;

				solve;
				if solve_result_num <> 0 then {
					print "Run number ", q, "of", totalRuns > (filename);
					print "Locally optimal solution not found, continuing with next step" > (filename);
					continue runControl;
				}
				
				let objCandidate := objective;
				
				display ch, s+dir, ho, hso;
				display objCandidate-BESTobj;
				if objCandidate > BESTobj then {	#Check and save if better
					let BESTobj := objCandidate;
					let {j in 1..numberOfH} ho[j] := 0;
					let {j in 1..numberOfH} hs[j] := ch[j];

					
					let s:=s+dir;
					let ss:=s;
					let o:=1;
				}
				else {
					if dir = 1 and o <> 1 then {	#Change direction
						let dir:=-1;
					}
					else {	#Optimal solution found if dir=-1 and BESTobj > objCandidate
						let dir:=1;
						let hso:=1;
						print "s is optimal";
						display s, BESTobj;
					}
				}
			} else {
				print "S is 1, direction -1:", s+dir, "< 0 -> S is optimal";
				let dir:=1;
				let hso:=1;
				display s, BESTobj;
			}
		}
		if ho[1] = 1 and ho[2] = 1 and ho[3] = 1 and ho[4] = 1 and ho[5] = 1 and hso = 1 then {	#Everything optimal
			print "Optimal harvesting strategy obtained, results pending";
			print "";
			print "Run number", q, "of", totalRuns;
			display hs, ss;
			let finished := 1;
			let {t in T} g[t] := 0;
			let {j in 1..numberOfH} g[hs[j]] := 1;
			
			let {i in hs[5]+(5*ss)..maxt by 5*ss} g[i] := 1;

			solve;
			if solve_result_num <> 0 then {
				print "Run number ", q, "of", totalRuns > (filename);
				print "Locally optimal solution not found, continuing with next step" > (filename);
				continue runControl;
			}
			let objCandidate := objective;
			display objCandidate-BESTobj;
			if objCandidate > BESTobj then let BESTobj := objCandidate;
			if BESTobj > GBestObj then { #Check if current best result is the overall best of all runs.
				let GBestObj := BESTobj;
				let {j in 1..numberOfH} Ghs[j] := hs[j];#If so, save the times of harvests
				let Gss := ss;
			}
			#Save the different runs into files
			print "Run number ", q, "of", totalRuns > (filename);
			print "Candidate:" > (filename);
			display hs, ss, BESTobj > (filename);
			print "Current best:" > (filename);
			display Ghs, Gss, GBestObj > (filename);
			print "----------" > (filename);
			print "" > (filename);
			let BESTobj := 0;
		}
	}

} #End of for-loop

#Calculate once more the solution based on best optimal strategy so we can get all the details out
let {t in T} g[t] := 0;

let {j in 1..numberOfH} g[Ghs[j]] := 1;

let {i in Ghs[numberOfH]+5*Gss..maxt by 5*Gss} g[i] := 1;
	
solve;
display C, Ghs, Gss, objective, GBestObj;

close;