(!******************************************************
   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. 
    
   
   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
uses "mmsp"
(!-----------------------------------------------------------------------------------
 Based on the forest planning model described in "On a new collection of
 stochastic linear programming problems", K.A. Ariyawansa and A.J. Felt.

 Dash Optimization
-----------------------------------------------------------------------------------!)
declarations
	DIR="Data/"
end-declarations


declarations
	T           = 7
	K           = 8
	MAXDISCRET  = 3
	
	index                : integer
	alpha, beta,
	delta, gama          : real
	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

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

! 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)
(!
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