(!******************************************************
   Xpress-SP Example
   ======================
   ``````````````````
   FILE: ForestRelax.mos
   
   DESCRIPTION: Stochastic 2-stage linear model
   
   TYPE: Forest Planning
         
   MODEL:    	   	
   	The goal is to maximize the value of timber, both cut 
   and remaining, after the last period of planning. The constraints 
   in this model pertain to the amount of forest cut in a single period, 
   the age of trees, and the likelihood that trees remaining uncut are 
   destroyed by fire. Each tree belongs to one of K age classes of equal 
   length. The length that defines the age classes is also used to define 
   the length of stages of the planning period. Thus, trees belonging to 
   age class A in a given stage will belong to age class A+1 in the next 
   stage if they are not cut or destroyed by fire. Whenever any of these 
   events take place, new trees are planted and, the number of trees in 
   the first age class equals the sum of these two quantities. 
   
   The variables of the model are the total area of forest for each age class, 
   for each period, the amount of forest harvested for each age class in each 
   period, and penalty factors for not satisfying constraints on the amount 
   of harvested area in each period. 
    
   We show how scenario tree can be reduced by aggregation or deletion. Scnerios
   that have small probabilities are aggregated. Similar scenarios that do not yield
   objective value close to expected value of objective function may be deleted.
   
   FURTHER INFO: see Xpress-SP guide
   see K. A. Ariyawansa and A. J. Felt, “Test – Problem Collection for Stochastic 
   Linear Programming" (http://www.uwsp.edu/math/afelt/slptestset/doc.pdf).
   
   DATE: Oct 2006			
   
   (c) 2006 Dash Associates
       Author: Nitin Verma
       Dash Optimization Inc, NJ  
       
*******************************************************!)

model forest_relax_doc
uses "mmsp"


parameters
	EPS = 1e-5
	DIR = 'Data/'
end-parameters

declarations
	T           = 7
	K           = 8
	MAXDISCRET  = 3
	
	index                : integer
	alpha, beta,
	delta, gama          : real
	ScenariosAgg,
	ScenariosDel         : set of integer
	yield, v             : array (1..K) of real
	ProbDiscret,
	ValuesDiscret        : array (1..MAXDISCRET, 1..MAXDISCRET) of real
	s1                   : array (1..K) of real
end-declarations

declarations
	P                    : array (1..K, 2..T+1) of sprand
	x                    : array (1..K, 1..T) of spvar
	s                    : array (1..K, 1..(T+1)) of spvar
	p1, p2               : array (2..T) of spvar
	
	HarvestLimit         : array (1..K, 1..T) of splinctr
	Fix                  : array (1..K) of splinctr
	TotalArea            : array (1..K, 2..(T+1)) of splinctr
	YLowerLimit          : array (2..T) of splinctr
	YUpperLimit          : array (2..T) of splinctr
end-declarations

! Read data
initializations from DIR+'forest.dat'
	yield  v  alpha  beta  delta  gama
	ProbDiscret  ValuesDiscret  s1
end-initializations

!---------------------------------------------------------
! Create scenario tree
declarations
	Stages = 1..(T+1)
	Branches : array (1..T) of integer
	Probabilities : array (1..MAXDISCRET) of real
	PValues : array (1..MAXDISCRET) of real
end-declarations

spsetstages(Stages)
forall (t in 2..T+1, i in 1..K) spsetstage(P(i,t), t)

Branches := [1, 3,3,3, 3,3,3]
spcreatetree(Branches)

forall (t in 2..T+1) do
	forall (i in 1..MAXDISCRET) PValues(i) := ValuesDiscret(Branches(t-1),i)
	forall (i in 1..K)
		spsetrand(P(i,t), PValues)
end-do

forall (t in 1..T) do
	forall (i in 1..MAXDISCRET)
		Probabilities(i) := ProbDiscret(Branches(t),i)
	spsetprobcond(t+1,Probabilities)
end-do

!setparam('xsp_implicit_stage', true)
setparam('xsp_scen_based', false)
spgentree
!---------------------------------------------------------

! Objective function
ObjFnc := sum(t in 1..T, i in 1..K) (delta^(t-1))*yield(i)*x(i,t) +
          sum(i in 1..K)            (delta^T)*v(i)*s(i,T+1) -
          gama*sum(t in 2..T) (delta^(t-1))*(p1(t) + p2(t))
    ! spvars s(T+1) should be associated to stage T

! Constraints
forall (j in 1..K, t in 1..T) do
	HarvestLimit(j,t) := x(j,t) <= s(j,t)
end-do
forall (j in 1..K) Fix(j) := s(j,1) = s1(j)

forall (i in 1..K, t in 2..(T+1)) do
	if (i = 1) then
		TotalArea(i,t) := s(i,t) = sum(j in 1..K) (1 - P(j,t))*x(j,t-1) +
	    	                sum(j in 1..K)      P(j,t) *s(j,t-1)
	else
		if (i = K) then
			TotalArea(i,t) := s(i,t) = -(1 - P(K-1,t))*x(K-1,t-1) - (1 - P(K,t))*x(K,t-1) +
	    		                 (1 - P(K-1,t))*s(K-1,t-1) + (1 - P(K,t))*s(K,t-1)
		else
			TotalArea(i,t) := s(i,t) = -(1 - P(i-1,t))*x(i-1,t-1) +
	    		                 (1 - P(i-1,t))*s(i-1,t-1)
	    end-if
	end-if
end-do

forall (t in 2..T) do
	YUpperLimit(t) := - sum(j in 1..K) yield(j)*x(j,t) + alpha*sum(j in 1..K) yield(j)*x(j,t-1) <= p1(t)
	YLowerLimit(t) :=   sum(j in 1..K) yield(j)*x(j,t) -  beta*sum(j in 1..K) yield(j)*x(j,t-1) <= p2(t)
end-do

! Associate Smpvars to stages
forall (t in 1..T) do
	forall (i in 1..K) spsetstage(x(i,t),t)
	forall (i in 1..K) spsetstage(s(i,t),t)
end-do
forall (i in 1..K) spsetstage(s(i,T+1),T)
forall (t in 2..T) do
	spsetstage(p1(t),t)
	spsetstage(p2(t),t)
end-do

setparam('xsp_verbose', true)
setparam('xsp_ive_enable',true)

! Optimize
maximize(ObjFnc)

! Report
writeln("---->",getobjval)
!forall (i in 1..K)
!	writeln(getsol(x(i,1)))


writeln('------------------------------------------------')
Scenarios := 1
forall (t in 1..T) Scenarios := Scenarios * Branches(t)

ScenariosAgg := {}
forall (scen in 1..Scenarios) do
	if (spgetprobscen(scen)<0.001) then
		ScenariosAgg := ScenariosAgg + {scen}
	end-if
end-do


 

scen := getsize(ScenariosAgg)
AggOcurred := false
while (scen >= 1) do
	siblings := {ScenariosAgg(scen)}
	forall (j in 1..(MAXDISCRET-1) | scen-j>0 ) do
		if (ceil(ScenariosAgg(scen-j)/Branches(T))=ceil(ScenariosAgg(scen)/Branches(T))) then
			siblings := siblings + {ScenariosAgg(scen-j)}
		else
			break
		end-if
	end-do
	if (getsize(siblings)>1) then
		writeln('aggregating ', siblings)
		spaggregate(siblings)
		AggOcurred := true
	end-if
	Scenarios := Scenarios - getsize(siblings) + 1
	scen := scen - getsize(siblings)
end-do
	
maximize(ObjFnc)
writeln("---->",getobjval)

!(!
!============================================================
while (true) do
	ScenariosDel := {}
	forall (sce in 1..Scenarios)
		if (speval(ObjFnc,sce)<0.95*getobjval) then
			ScenariosDel := ScenariosDel + {sce}
		end-if
	if (getsize(ScenariosDel)=0) then
		break
	end-if
	spdelscen(ScenariosDel)
	Scenarios := Scenarios - getsize(ScenariosDel)
	maximize(ObjFnc)
	writeln(getobjval)
end-do
!)

!(!
forall (i in 1..Scenarios, t in 2..T) do
	if (getsolinscen(p1(t),i) <> 0) then
		writeln('penalty1(stage=',t,',scen=',i,') ',getsolinscen(p1(t),i))
	end-if
	if (getsolinscen(p2(t),i) <> 0) then
		writeln('penalty2(stage=',t,',scen=',i,') ',getsolinscen(p2(t),i))
	end-if
end-do!)

end-model