(!******************************************************
   Xpress-SP Example
   ======================
   ``````````````````
   FILE: HT_simple.mos
   
   DESCRIPTION:   multi-stage stochastic linear model
   
   TYPE:  Energy: Power management in a hydro-thermal 
   						system under uncertainty
   
   	There are T time periods. I and J are the numbers of thermal 
   and hydro units, respectively. For a thermal unit i in period t
   u(i,t) in {0,1} is its commitment (1 if on, 0 if off) and p(i,t)
   its production, with p(i,t)=0 if u(i,t)=0, p(i,t) in 
   [P_min(i,t),P_max(i,t)] if u(i,t)=1. Additionally there are minimum 
   up/down-time requirements : when unit i is switched on (off) it must 
   remain on (off) for at least tau_up(i) (tau_dn(i)) periods. For a
   hydro plant j, s(j,t) and w(j,t) are its generation and pumping 
   levels in period t with upper bounds S_max(j,t) and W_max(j,t) resp,
   and l(j,t) is the storage volume in the upper dam at the end of period t
   with upper bound L_max(j,t). The water balance relates l(j,t) with l(j,t-1),
   s(j,t), w(j,t) and water inflow gamma(j,t) using the pumping efficiency eta(j).
   The initial and final volumes are specified by l_in(j) and l_end(j).
   
   The basic system requirement is to meet the electric load. Another
   important requirement is the spinning reserve constraint. To maintain
   reliability (compensate sudden load peaks or unforeseen outages of units)
   the total commited capacity should exceed the load in every period by a
   certain amount (e.g., a fraction of the demand). The load and the spinning
   reserve during period t are denoted by d(t) and r(t), respectively.
   
   Efficient operation of pumped storage hydro plants exploits daily cycles of 
   load curves by generating during peak load periods and pumping during off-peak 
   periods. Since the operating costs of hydro plants are usually negligible, the
   total system cost is given by the sum of startup and operating costs of
   all thermal units over the whole scheduling horizon. The fuel cost C(i,t) for
   operating thermal unit i during period t has the form
   			C(i,t)=max(l in 1..L) {a(i,l,t).p(i,t) + b(i,l,t).u(i,t)}
   	
   	The startup cost of unit i depends on its downtime; it may vary from a maximum cold
   	start value to a much smaller value when the unit is still relatively close to its 
   	operating temperature. This is modeled by the startup cost:
   		S(i,t)=max(tau in 0..tau_c(i)) c(i,tau).{u(i,t)-sum(k in 1..tau) u(i,k-tau)}
   	where 0= c(i,0) < ... < c(i,tau_c(i)) are fixed cost coefficients. tau_c(i) is the
   	cool-down time of unit i, c(i,tau_c(i)) is its maximum cold-start cost,
   	u(i,tau) for tau in tau_ini..0 are given initial values.
   		where tau_ini=1-max(i)(tau_c(i),1-tau_up(i),tau_dn(i)-1)
   	
   MODEL:
   			min E{sum(t,i) [C(i,t)+S(i,t)]}
   			s.t
   			P_min(i,t).u(i,t)<=p(i,t)<=P_max(i,t).u(i,t)
   			u(i,tau)-u(i,tau-1)<=u(i,t), 				t-tau_up(i)<tau<t
   			u(i,tau-1)-u(i,tau)<=1-u(i,t), 				t-tau_dn(i)<tau<t
   			0<=s(j,t)<=S_max(j,t), 0<=w(j,t)<=W_max(j,t), 0<=l(j,t)<=L_max(j,t) 
   			l(j,t)=l(j,t-1)-s(j,t)+eta(j).w(j,t)+gamma(j,t)
   			l(j,0)=l_in(j), l(j,T)=l_end(j)
   			sum(i) p(i,t)+ sum(j) [s(j,t)-w(j,t)] >= d(t)
   			sum(i) [u(i,t).P_max(i,t) - p(i,t)] >= r(t)
   
   	Uncertainty:
   	
   	1. Aggregate stages: assign 
   			set stage of 1..t1		to 1
   			set stage of t1+1..t2	to 2
   			.....
   			set stage of tK+1..T	to NumStages
   	2. Assume demand process (may be identified based on a time series or regression model)
   		d(t)=alpha(t)+si(t)*beta*d(t-1)+lambda*eps(t)		: alpha(t),beta,lambda are given
   		where	si(t) ~Unif(min_si, max_si)		: min_si, max_si are given
				eps(t) ~Normal(0,1)
	3. Based on the model determine an initial structure of the load tree. Compute scenario 
	   	values, using the sample means and standard deviations of the simulated scenarios. Run 
	   	NumSimRun simulations and estimate mean demand: mean_dem(t) and variance in demand: var_dem(t)
	4. Re-build a stable and smaller scenario tree using the estimates and aggregated stages.
		a) set si(t), eps(t) to stage 0
		b) set distribution of delta(1),delta(2),...,delta(NumStages) = {1 w.p 0.5, -1 w.p 0.5}
		c) ensure that means and variances still hold by re-modeling
			d(t)= mean_dem(t)+ sum(t_ in 2..t) delta(stage(t_))*sqrt(var_dem(t_)/2^(t_-1))
		d) generate exhaustive tree based on distributions of delta(1),delta(2),...,delta(NumStages)
	5. Reduce scenario tree
		a) define distance between vectors v1 and v2: c(h,v1,v2)=max{1,||v1||^(h-1),||v2||^(h-1)}.||v1-v2||
			where h is a parameter, and ||.|| is a Eucleadian norm
		b) Let realizations of delta=[delta(1),delta(2),...,delta(NumStages)] in scenario k be vector del(k)
		   Then delete scenario k = arg min (l) { Pr(l). min (s<>l) c(h,del(l),del(s)) }
		   Then, roughly speaking, deletion occurs where scenarios are close as measured by the
		   distance c(h,.,.) or where probabilities are small.
		   The reduced scenario tree has one less scenario; scenario k is deleted, and Pr(s)+=Pr(k).
		   This reduction procedure may be repeated until a prescribed # scenarios in the reduced 
		   measure is attained.
		c) Create a scenario tree with final prescribed # scenarios
			1. create tree with 2 branches from node in 1st stage, and 1 branch per node from each node in other stages
			2. set the revised realized values of delta, and probabilities of scenarios
	6. Restructure scenario tree
		a)	For each node(t,n) let PossScen(t,n) be set of original scenarios that branches off from node(t,n) that have same realized values of del(t+1)
		b)  Find all such possible sets of the PossScen(t,n)'s corresponding to the realizations of del(t+1)
		c)	Each of the set PossScen(t,n) forms a new child of node(t,n) with the corresponding realized value of delta(t+1) and conditional probability = 1/ unconditional probabilty of node(t,n) * sum(s in PossScen(t,n)) Pr(s)
		
		In above steps keep track of #children of node(t,n), #node in stage t, and un conditional probaility of node(t,n)
		
	7. Solve optimization problem with reduced scenario tree
			
   
   FURTHER INFO:  
   see http://www.dashoptimization.com/home/services/publications/casestudies/casestudy_22.html
   
   see N. Gröwe-Kuska, K.C. Kiwiel, M.P. Nowak, W. Römisch and I. Wegner: Power management 
   in a hydro-thermal system under uncertainty by Lagrangian relaxation, in: Decision Making 
   under Uncertainty: Energy and Power (C. Greengard, A. Ruszczynski eds.), 
   IMA Volumes in Mathematics and its Applications Vol. 128, Springer-Verlag, New York 2002, 39-70
   (http://www-iam.mathematik.hu-berlin.de/~romisch/papers/pmhsfin.ps)
   
   DATE: Oct 2006			
   
   (c) 2006 Dash Associates
       Author: Nitin Verma
       Dash Optimization Inc, NJ  
       
*******************************************************!)


model "Power management in a hydro-thermal system under uncertainty"

uses 'mmsp','mmrng'

parameters
	T=24
	NumStages=5
	NumSim=50000
	NumSimRun=3
	RedSize=3
	
	I=3
	J=2
	L=5
	
	RelStop=0.02
	use_sos=false
	solve_EV_with_rec=false
	
	DIR="Data/"
	RandomDatFile=DIR+"HTdist.dat"
	DatFile=DIR+"HT.dat"	
end-parameters

forward procedure SetupProblem(scenario_count: integer)
forward procedure SolveProblem
forward procedure runActualSimulation(simRunNum:integer,numSimRun:integer)
forward procedure buildScenTreeFromSim
forward procedure ReduceScenarios(h:integer,S_:integer)
forward function EucledianNorm(v:array(range) of integer):real
forward function getMinInd(a:array(range) of real,r:set of integer,minval:real):integer
forward procedure ReStructureTree(w:array(range,range) of integer,prob_scen:array(range) of real)
forward procedure UnHideFirstStage


!##################################### Stochastic Data ##############################################

declarations
	TimeBlocks=1..T
	Stages=1..NumStages
	AggStage:array(TimeBlocks) of integer

	alpha:array(TimeBlocks) of real
	beta,lambda,min_si,max_si:real
	si:array(1..T) of sprand !~Unif(min_si, max_si)
	eps:array(1..T) of sprand !~Normal(0,1)
	delta:array(1..NumStages) of sprand
	d: array (TimeBlocks) of sprandexp			
	mu_dem,sigma_sq_dem:array(TimeBlocks,1..NumSimRun) of real	
	mean_dem,var_dem:array(TimeBlocks) of real
end-declarations

!HT
initializations from RandomDatFile
	AggStage alpha beta lambda min_si max_si
end-initializations

setparam('xsp_disp_warnings',false)
spsetstages(Stages)

!##################################### Real Data ##############################################

declarations
	ThermalUnits=1..I
	HydroUnits=1..J
	
	tau_up,tau_dn,tau_c: array (ThermalUnits) of integer
	P_max,P_min: array (ThermalUnits) of real
	S_max,W_max,SL_max,l_in,l_end,eta: array (HydroUnits) of real
	r: array (TimeBlocks) of real
	a,b: array (ThermalUnits,TimeBlocks,1..L) of real
	c: array (ThermalUnits,range) of real
	gamma: array(HydroUnits,TimeBlocks) of real
	price: array(TimeBlocks) of real
	prodt,cost:array(ThermalUnits,TimeBlocks,1..L) of real
	tau_ini:integer	
end-declarations

initializations from DatFile
	tau_up tau_dn tau_c P_max P_min S_max W_max SL_max l_in l_end eta gamma r a b c price 
end-initializations

tau_ini:=1-max(i in ThermalUnits)(maxlist(tau_c(i),1-tau_up(i),tau_dn(i)-1))	
forall(i in ThermalUnits,t in TimeBlocks,l_ in 1..L) prodt(i,t,l_):=P_min(i)+(P_max(i)-P_min(i))*(l_-1)/(L-1)
forall(i in ThermalUnits,t in TimeBlocks,l_ in 1..L) cost(i,t,l_):=a(i,t,l_)*prodt(i,t,l_)+b(i,t,l_)


!###################################### Stochastic Variables and Constraints #################################
	
declarations
	p,S,C:array(ThermalUnits,TimeBlocks) of spvar
	u:array(ThermalUnits,(tau_ini..0)+TimeBlocks) of spvar
	s,w,l:array(HydroUnits,TimeBlocks) of spvar
	lam:array (ThermalUnits,TimeBlocks,1..L) of spvar
	
	TotalCost:splinctr
	UBprod:array(ThermalUnits,TimeBlocks) of splinctr
	LBprod:array(ThermalUnits,TimeBlocks) of splinctr
	UBstoreageBal:array(HydroUnits,TimeBlocks) of splinctr
	Minuptime,Mindowntime: dynamic array(ThermalUnits,range,TimeBlocks) of splinctr
	Load:array(TimeBlocks) of splinctr
	SpinningReserve:array(TimeBlocks) of splinctr
	ProdLevel,FuelCost,SumToOne:array(ThermalUnits,TimeBlocks)  of splinctr
	Prodtn:dynamic array(ThermalUnits,TimeBlocks)  of splinctr
	StartupCost:array(ThermalUnits,range,TimeBlocks) of splinctr
end-declarations

!##############################################  Main model ##################################################

SetupProblem(10) ! set up problem to reduce # scenarios to 10
SolveProblem ! set up model and solve
exit(0)

!##############################################  All Procedures and Functions ################################
procedure SetupProblem(scenario_count: integer)

	!1. set stages
	!2. set demand model
	!3. sim simualtion
	!4. build another scenario tree using simulation
	!5. reduce size of scenario tree

	forall(t in 1..T) do
		spsetstage(si(t),AggStage(t))
		spsetstage(eps(t),AggStage(t))
	end-do
	
	forall(t in TimeBlocks) do
		if t=1 then 
			d(t):=alpha(t)
		else
			d(t):=alpha(t)+si(t)*beta*d(t-1)+lambda*eps(t)
		end-if
	end-do
	
	forall(sim in 1..NumSimRun) do
		RNGsetseed(sim)
		runActualSimulation(sim,NumSim)
	end-do
	
	forall(t in TimeBlocks) mean_dem(t):=(1/NumSimRun)*sum(sim in 1..NumSimRun) mu_dem(t,sim)
	forall(t in TimeBlocks) var_dem(t):=(1/NumSimRun)*sum(sim in 1..NumSimRun) sigma_sq_dem(t,sim)
	buildScenTreeFromSim
	
	if scenario_count<spgetscencount and scenario_count>0 then
		ReduceScenarios(2,scenario_count)
	end-if
end-procedure

procedure SolveProblem
	
	!1. set stages of variables
	!2. build constraints
	!3. optimize
	
	forall(t in TimeBlocks) do
		forall(i in ThermalUnits) do
			spsetstage(p(i,t),AggStage(t))
			spsetstage(u(i,t),AggStage(t))
			spsetstage(S(i,t),AggStage(t))
			spsetstage(C(i,t),AggStage(t))
			forall(l_ in 1..L) spsetstage(lam(i,t,l_),AggStage(t))
		end-do
		forall(j in HydroUnits) do
			spsetstage(s(j,t),AggStage(t))
			spsetstage(w(j,t),AggStage(t))
			spsetstage(l(j,t),AggStage(t))
		end-do		
	end-do
	
	forall(i in ThermalUnits, t in (tau_ini..0)+TimeBlocks) u(i,t) is_binary
	
	forall(t in tau_ini..0,i in ThermalUnits) do
		spsetstage(u(i,t),1)
		u(i,t)=0
	end-do
	
	
	forall(i in ThermalUnits,t in TimeBlocks) UBprod(i,t):=p(i,t)<=P_max(i)*u(i,t)	
	forall(i in ThermalUnits,t in TimeBlocks) LBprod(i,t):=p(i,t)>=P_min(i)*u(i,t)
	forall(j in HydroUnits,t in TimeBlocks) s(j,t)<=S_max(j)
	
	forall(j in HydroUnits,t in TimeBlocks) w(j,t)<=W_max(j)
	forall(j in HydroUnits,t in TimeBlocks) l(j,t)<=SL_max(j)
	forall(j in HydroUnits,t in TimeBlocks)
		if t=1 then
			UBstoreageBal(j,t):=l(j,t)=l_in(j)-s(j,t)+eta(j)*w(j,t)+gamma(j,t)
		elif t=T then
			UBstoreageBal(j,t):=l_end(j)=l(j,t-1)-s(j,t)+eta(j)*w(j,t)+gamma(j,t)
		else
			UBstoreageBal(j,t):=l(j,t)=l(j,t-1)-s(j,t)+eta(j)*w(j,t)+gamma(j,t)
		end-if
	
	forall(i in ThermalUnits, t in TimeBlocks, 
		tau in (t-tau_up(i)+1)..(t-1) | t-tau_up(i)+1<=t-1 and tau-1>=tau_ini)
			Minuptime(i,tau,t):=u(i,tau)-u(i,tau-1)<=u(i,t)
	forall(i in ThermalUnits, t in TimeBlocks, 
		tau in (t-tau_dn(i)+1)..(t-1) | t-tau_dn(i)+1<=t-1 and tau-1>=tau_ini)
			Mindowntime(i,tau,t):=u(i,tau-1)-u(i,tau)<=1-u(i,t)
	
	forall(t in TimeBlocks) Load(t):=sum(i in ThermalUnits) 
								p(i,t)+sum(j in HydroUnits) (s(j,t)-w(j,t))>=d(t)
	forall(t in TimeBlocks) SpinningReserve(t):=sum(i in ThermalUnits) (P_max(i)*u(i,t)-p(i,t))>=r(t)

	forall(i in ThermalUnits, t in TimeBlocks) do
		ProdLevel(i,t):=p(i,t)=sum(l_ in 1..L) prodt(i,t,l_)*lam(i,t,l_)
		FuelCost(i,t):=C(i,t)=sum(l_ in 1..L) cost(i,t,l_)*lam(i,t,l_)
		SumToOne(i,t):=sum(l_ in 1..L) lam(i,t,l_)=u(i,t)
		if use_sos then
			Prodtn(i,t):=sum(l_ in 1..L) prodt(i,t,l_)*lam(i,t,l_)
			Prodtn(i,t) is_sos2
		end-if
	end-do

	forall(i in ThermalUnits,t in TimeBlocks,tau in 0..tau_c(i))
		StartupCost(i,tau,t):=S(i,t)>=c(i,tau)*(u(i,t)-sum(k in 1..tau) u(i,t-k))
	
	TotalCost:=sum(i in ThermalUnits,t in TimeBlocks) (C(i,t)+S(i,t))
	
	setparam('xsp_scen_based',false)
	setparam('xprs_miprelstop',RelStop)
		
	if solve_EV_with_rec then
		minimize(TotalCost,XSP_EV+XSP_REC)
		UnHideFirstStage
	else
		minimize(TotalCost)
	end-if
	writeln("Total Cost= ",getobjval)
end-procedure



procedure runActualSimulation(simRunNum:integer,numSimRun:integer)
	
	declarations
		demVal:array(TimeBlocks,1..numSimRun) of real
	end-declarations

	forall(sim in 1..numSimRun) do
		forall(t in TimeBlocks) do			
			if t=1 then
				demVal(t,sim):=alpha(t)
			else
				!d(t):=alpha(t)+si(t)*beta*d(t-1)+lambda*eps(t)
				demVal(t,sim):=alpha(t)+RNGuniform(min_si,max_si)*beta*demVal(t-1,sim)+lambda*RNGnormal(0,1)				
			end-if
		end-do
	end-do
		
	forall(t in TimeBlocks) mu_dem(t,simRunNum):=1/numSimRun*sum(sim in 1..numSimRun) demVal(t,sim)
	forall(t in TimeBlocks) sigma_sq_dem(t,simRunNum):=(1/(numSimRun-1))*sum(sim in 1..numSimRun) (demVal(t,sim)-mu_dem(t,simRunNum))^2
end-procedure


procedure buildScenTreeFromSim
	declarations		
		bin_val,bin_prob:array(1..2) of real
		val1,prob1:array(1..1) of real
	end-declarations
	
	forall(t in 1..T) do
		spsetstage(eps(t),XSP_STAGE_ZERO)
		spsetstage(si(t),XSP_STAGE_ZERO)	
	end-do
	forall(stg in 1..NumStages) spsetstage(delta(stg),stg)	
	forall( t in TimeBlocks) do
		if t=1 then
			d(t):=mean_dem(t)
		else
			d(t):=mean_dem(t)+sum(t_ in 2..t) delta(AggStage(t_))*sqrt(var_dem(t_))/2^((t_-1)/2)
		end-if
	end-do
	
	val1:=[0];prob1:=[1];bin_val:=[1,-1];bin_prob:=[0.5,0.5]
	forall(stg in 1..NumStages)
		if stg=1 then
			spsetdist(delta(stg),val1,prob1)
		else
			spsetdist(delta(stg),bin_val,bin_prob)
		end-if
	setparam('xsp_implicit_stage',false)
	spgenexhtree
end-procedure


procedure ReduceScenarios(h:integer,FinalNumScen:integer)

	InitialNumScen:=spgetscencount
	
	declarations
		NotExistsScen:array(1..InitialNumScen) of boolean
		closest_scen:array(1..InitialNumScen) of integer
		initial_realized_vals:array(2..NumStages,1..InitialNumScen) of integer
		final_realized_vals:array(2..NumStages,range) of integer
		InitScenProbs:array(1..InitialNumScen) of real
		FinalScenProbs:array(range) of real
		c_h:dynamic array(1..InitialNumScen,1..InitialNumScen) of real
		vec_k,vec_sk,vec_diff:array(2..NumStages) of integer
		vec_c_sk,vec_mindist:dynamic array(1..InitialNumScen) of real
	end-declarations
	
	forall(k in 1..InitialNumScen) InitScenProbs(k):=spgetprobscen(k)
	forall(stg in 2..NumStages,k in 1..InitialNumScen) initial_realized_vals(stg,k):=integer(speval(delta(stg),k))

	forall(k in 1..InitialNumScen,sk in k+1..InitialNumScen) do
		forall(stg in 2..NumStages) do
			vec_k(stg):=initial_realized_vals(stg,k) !all realizations in scenario k
			vec_sk(stg):=initial_realized_vals(stg,sk)!all realizations in scenario sk
			vec_diff(stg):=vec_k(stg)-vec_sk(stg)
		end-do
		c_h(k,sk):=maxlist(1,EucledianNorm(vec_k)^(h-1),EucledianNorm(vec_sk)^(h-1))*EucledianNorm(vec_diff)
	end-do
	
	ExistingScens:=union(k in 1..InitialNumScen) {k}
	repeat		
		forall(k in ExistingScens) do
			set_sk:=union(sk in k+1..InitialNumScen | not NotExistsScen(sk)) {sk}
			forall(sk in set_sk) vec_c_sk(sk):=c_h(k,sk)
			minval:=min(sk in set_sk) vec_c_sk(sk)
			closest_scen(k):=getMinInd(vec_c_sk,set_sk,minval)
			vec_mindist(k):=InitScenProbs(k)*minval
		end-do
		
		minval_delscen:=min(k in ExistingScens) vec_mindist(k)
		delscen:=getMinInd(vec_mindist,ExistingScens,minval_delscen)
		InitScenProbs(closest_scen(delscen))+=InitScenProbs(delscen)
		InitScenProbs(delscen):=0
		NotExistsScen(delscen):=true
		ExistingScens-={delscen}
		
	until getsize(ExistingScens)<=FinalNumScen
	
	declarations
		Branches:array(2..NumStages) of integer
	end-declarations
	
	Branches(2):=FinalNumScen; forall(stg in 3..NumStages) Branches(stg):=1
	spcreatetree(Branches)
	node:=0	
	forall(k in ExistingScens) do
		node+=1
		forall(stg in 2..NumStages) do
			spsetrandatnode(delta(stg),node,initial_realized_vals(stg,k))
			spsetprobcondatnode(stg,node,if(stg=2,InitScenProbs(k),1))				
		end-do
	end-do
	
	spgentree	

	forall(stg in 2..NumStages,k in 1..spgetscencount) final_realized_vals(stg,k):=integer(speval(delta(stg),k))
	forall(k in 1..spgetscencount) FinalScenProbs(k):=spgetprobscen(k)
	ReStructureTree(final_realized_vals,FinalScenProbs)	
	
	
end-procedure

function EucledianNorm(v:array(range) of integer):real
	returned:=sqrt(sum(t in 2..NumStages) v(t)^2)	
end-function

function getMinInd(arr:array(range) of real,ind_set:set of integer,minval:real):integer
	returned:=0
	forall(i in ind_set)
		if arr(i)=minval then
			returned:=i
			break
		end-if
end-function

procedure ReStructureTree(delta_val:array(range,range) of integer,prob_scen:array(range) of real)
	
	declarations
		N:array(1..NumStages) of integer
		
		Agg,AggInto:array(range) of boolean
		
		K:array(range,range) of integer !K(stg,n) = # possible scens from node(stg,n)
		C:array(range,range) of integer !C(stg,n) = # children of node(stg,n)
		
		CumProb:array(range,range) of real
		
		PossScen:array(range,range,range) of integer !{PossScen(stg,n,1),..,PossScen(stg,n,K(stg,n))}=set of possible scenarios from node(stg,n)
		Scen:array(range,range,range) of integer  !{Scen(stg,n,1),..,Scen(stg,n,C(stg,n))}=set of actual scenarios from node(stg,n)
		
		! Scen within PossScen
		! forall n_ is child of (stg,n) do
			! find (stg,n)-> K(stg+1,n_):
			! find (stg,n)-> {PossScen(stg+1,n_,1),..,PossScen(stg+1,n_,K(stg+1,n_))}
		
	end-declarations
	
	numScen:=spgetscencount
	N(1):=1;K(1,1):=numScen;forall(k in 1..numScen) PossScen(1,1,k):=k;CumProb(1,1):=1.0
	forall(stg in 1..NumStages-1) do
		forall(n in 1..N(stg)) do
			
			C(stg,n):=0
			
			forall(k in 1..numScen) Agg(k):=false
			forall(k in 1..numScen) AggInto(k):=false
			forall(sn in 1..K(stg,n)) do
				if not (Agg(PossScen(stg,n,sn)) or AggInto(PossScen(stg,n,sn))) then
					C(stg,n)+=1
					n_:=N(stg+1)+C(stg,n)
					
					Scen(stg,n,C(stg,n)):=PossScen(stg,n,sn)
					
					K(stg+1,n_)+=1					
					PossScen(stg+1,n_,K(stg+1,n_)):=Scen(stg,n,C(stg,n))
					
					Agg(PossScen(stg,n,sn)):=true
				end-if				
				forall(sk in 1..K(stg,n) | sk<>sn) do
					if not (Agg(PossScen(stg,n,sk)) or AggInto(PossScen(stg,n,sk))) then
						if delta_val(stg+1,PossScen(stg,n,sk))=delta_val(stg+1,PossScen(stg,n,sn)) then
							AggInto(PossScen(stg,n,sk)):=true
							K(stg+1,n_)+=1
							PossScen(stg+1,n_,K(stg+1,n_)):=PossScen(stg,n,sk)
														
						end-if
					end-if
				end-do				
			end-do
			
			spaddchildren(stg,n,C(stg,n))
			forall(ch in 1..C(stg,n)) do
				spsetrandatnode(delta(stg+1),N(stg+1)+ch,delta_val(stg+1,Scen(stg,n,ch)))
				CumProb(stg+1,N(stg+1)+ch):=sum(sk in 1..K(stg+1,N(stg+1)+ch)) prob_scen(PossScen(stg+1,N(stg+1)+ch,sk))
				spsetprobcondatnode(stg+1,N(stg+1)+ch,CumProb(stg+1,N(stg+1)+ch)/CumProb(stg,n))
			end-do
			N(stg+1)+=C(stg,n)
		end-do
	end-do

	spgentree

end-procedure

procedure UnHideFirstStage
	
	forall(t in TimeBlocks | AggStage(t)=1) do
		forall(i in ThermalUnits) do
			spunfix(p(i,t))	
			spunfix(S(i,t))
			spunfix(C(i,t))
			spunfix(u(i,t))
			forall(l_ in 1..L) spunfix(lam(i,t,l_))
			
			spsethidden(Load(t),false)
			spsethidden(SpinningReserve(t),false)
			spsethidden(UBprod(i,t),false)
			spsethidden(LBprod(i,t),false)
			
			spsethidden(ProdLevel(i,t),false)
			spsethidden(FuelCost(i,t),false)
			spsethidden(SumToOne(i,t),false)
			
			if exists(Prodtn(i,t)) then
				spsethidden(Prodtn(i,t),false)
			end-if
			forall(tau in 0..tau_c(i)) spsethidden(StartupCost(i,tau,t),false)
			forall(tau in (t-tau_up(i)+1)..(t-1) | t-tau_up(i)+1<=t-1 and tau-1>=tau_ini and exists(Minuptime(i,tau,t))) spsethidden(Minuptime(i,tau,t),false)
			forall(tau in (t-tau_up(i)+1)..(t-1) | t-tau_up(i)+1<=t-1 and tau-1>=tau_ini and exists(Mindowntime(i,tau,t))) spsethidden(Mindowntime(i,tau,t),false)
			
		end-do		
		forall(j in HydroUnits) do
			spunfix(s(j,t))	
			spunfix(w(j,t))
			spunfix(l(j,t))
			
			spsethidden(UBstoreageBal(j,t),false)
		end-do		
	end-do
	
end-procedure

end-model