(!******************************************************
   Xpress-SP Example
   ======================
   ``````````````````
   FILE: OptionContract.mos
   
   DESCRIPTION: Stochastic multi-stage linear model
   
   TYPE: Supply Chain: Option contracts

   	We consider the problem of a buyer who is engaged with a supplier 
   in a contract having periodical commitment with options, for a given 
   number of stages. At the beginning of the horizon, both the parties 
   agree to a fixed amount of product that would be delivered at each 
   stage to the buyer at a specified price. The supplier has a limited 
   amount of options available for each stage, with each unit of an option 
   giving the buyer a right to purchase a unit of product. Whenever this 
   option is exercised, the product is delivered at the next stage. 
   Initially, the buyer may also buy these options for each of the stages 
   at a specified price. At each stage, the buyer creates finished goods 
   from the total amount of product available, and sells them in the market. 
   The demand of finished goods at each stage is uncertain. Hence, at each 
   stage and at a specified exercise price, the buyer may buy additional 
   product (by exercising one or more of the options bought for that stage), 
   in order to meet the demands in the following stages. The unmet demand 
   in each stage is penalized and carried over to the next stage, whereas 
   excess quantity is stored in an inventory. If there is excess inventory 
   of products left at the end of the horizon, they are sold at a salvage 
   value, otherwise the unmet demand is lost.
            
   MODEL:
	Decision variables: 
		Q(t): Fixed order of products decided at stage 1, and delivered at stage t in {2,..,T} 
		M(t): Number of options bought at stage 1, which may be exercised at stage t in {2,..,T-1} 
		m(t): Number of options exercised at stage t, and delivered at stage t+1, for all t in {2,..,T-1} 
	State variables: 
		I(t): Finished goods inventory at stage t in {2,..,T} 
		I(t)+: Physical finished goods inventory at stage t in {2,..,T} 
		I(t)-: Backorder of finished goods inventory at stage t in{2,..,T} 
	Uncertainty: 
		D(t): Demand of finished goods at stage t in {2,..,T} 
		It is assumed that the demand at a stage is correlated with the demand in its previous stage, therefore they: 
		• form a conditionally heteroskedastic Gaussian process 
		• are normally distributed
		E[D(t+1) | D(t)=d(t)]= mu +sigma.eps(t+1), t=1
						= mu +rho.(d(t)-mu)+sigma.sqrt(1-rho^2).eps(t+1), 1<t<=T
		where mu and sigma are the unconditional mean and the variance of D(t) respectively, 
		and rho is the correlation coefficient between D(t) and D(t+1).
		eps(t) is assumed to be distributed Normally with mean 0 and varinace 1.
		
		In order to create scenario tree eps(t) is discretized into Neps(t) sample points.
		The n-th point will have a 
			value: x = -3+(n-0.5)*6/Neps(t) and 
			probability: p(x)= (1/sqrt(2*PI))*exp(-1/2.x^2)*6/Neps(t)
	Data: 
		v: Unit salvage value of finished goods (in $) 
		M: Bound on number of options that can be bought at each stage 
		o: Unit price for buying an option (in $) 
		e: Unit price of exercising an already bought option (in $) 
		r: Unit selling price of finished goods (in $) 
		p: Unit purchasing cost of a product from the supplier (in $) 
		s: Unit shortage cost for finished goods (in $) 
		h: Unit holding cost for finished goods (in $)
	Optimization:
		R(I+,I-): revenue= r[D(2)-(I(2)-)]+r.sum(t in 3..T) [D(t)+(I(t-1)-)-(I(t)-)]+v.(I(T)+)
		E(I+,I-,m,Q,M) : expenses=sum(t in 1..T-1) [e.m(t)+o.M(t)]+sum(t in 1..T) [h.(I(t)+)+s.(I(t)-)+p.Q(t)]
		max R(I+,I-) - E(I+,I-,m,Q,M) 
		s.t
		I(t)=(I(t)+) - (I(t)-), 1<t<=T
		I(t)=Q(t)-D(t) t=2
		I(t)=I(t-1)+Q(t)+m(t-1)-D(t)} 2<t<=T
		0<=m(t)<=M(t)<=M,  1<t<T
		Q(t), I(t)+,I(t)->=0, 1<t<=T

	The scenario tree may be symmetrix and exhaustive based on all Neps(t) realizations of eps(t), or
	assymetric, that is if the probability of visiting a node in a given stage is much smaller as 
	compared to other nodes in that stage, one may curtail the number of branches emerging from 
	that node because the detailed future information collected on that node by having many branches 
	will neither be fruitful nor affect the optimality of the solution.	Therefore one may SWITCH from
	'fine' branching to 'coarse' branching.
		
		#branches from a node (t,n) = Nfine(t) 		if Pr{visiting node(t,n)>=SwitchLevel(t)/#Nodes(t)
									= Ncoarse(t) 	otherwise
		
   
   	We may also model the chance constraint: Pr(net inventory in last stage >=0} > 1-beta by modeling
   	z<=beta as a 'global' constraint, where z=1	if I(T)- < 0
   											=0	otherwise 
   	which is expanded as sum(s in Sceanrios) Pr(s).z(s)<=beta internally in SP
   
   
   	Finally, we may calculate CVar(alpha) for a given alpha by solving:
   		CVar(alpha)=min {S + E[max(L-S,0)]/alpha : -inf<S<+inf}
   			where L is Loss = E(I+,I-,m,Q,M) - R(I+,I-)
   			
   	Hence, in the stochastic programming context, CVaR can be calculated by introducing variables CVaR, 
   	and S belonging to stage 1, and a non-negative variable z belonging to stage T. 
   	z represents max(L-S,0) therefore a constraint z>=L-S must also be introduced. 
   	Next, a bound c is introduced on CVaR. During minimization the term thata.CVaR
   	is added to the objective function where theta is sufficiently small, so that it does not perturb 
   	the optimal solution to the original problem. The constraint S+z/apha<=CVar is added and set as ‘global’ 
   	implying S+1/aplha sum(s in Sceanrios) Pr(s).z(s)<=CVar
   
   
   FURTHER INFO: 
   see Xpress-SP guide
   
   see A. Dormer, A. Vazacopoulos, N. Verma, H. Tipi – “Modeling & Solving 
   Stochastic Programming Problems in Supply Chain Management using Xpress-SP”, 
   “Supply Chain Optimization”, Editors: Joseph Geunes and Panos Pardalos
   
   see Ch. van Delft and J.-Ph.Vial, “A practical implementation of stochastic 
   programming: an application to the evaluation of option contracts in supply 
   chains”(http://www.unige.ch/hec/logilab/templeet/template/papiers/papier1ScVvDRev3.pdf)
   
   DATE: Oct 2006			
   
   (c) 2006 Dash Associates
       Author: Nitin Verma
       Dash Optimization Inc, NJ  
       
*******************************************************!)
model 'SCM Option Contract' !model name
uses 'mmsp' !mosel library for stochastic programming

parameters	
	Chance=false
	SolveCVaR=false
	Asymm=false	!if true, create a sceanrio tree based on discrete values of 
end-parameters

!procedure to generate symmetric scenario tree
forward procedure GenSymScenTree(N:array(range) of integer)
!procedure to generate asymmetric scenario tree
forward procedure GenAsymScenTree(Nfine:array(range) of integer,
				Ncoarse:array(range) of integer,SwitchLevel:array(range) of real)
!procedure to get discretized values and probabilities on Normal distribution
!for a given cardinality of state space
forward procedure GetNormDist(N:integer,x:array(range) of real,p:array(range) of real) 

declarations
	T=3 !number of stages
	Stages=1..T !set of stages

	r=12 !selling price of the finished goods
	h=0.5 !inventory holding cost
	s=6 !shortage cost
	v=2 !salvage value of finished goods
	o=1.5 !unit cost of an option
	e=8 !unit cost of exercising the option
	p=8 !unit purchase cost of product
	Mmax=10000 !upper bounds on the number of option
	
	mu=1500 !mean demand
	sigma=330 !deviation in demand
	rho=0.5 !correlation coefficient
	eps:array(2..T) of sprand !eps~Normal(0,1)
	D:array(2..T) of sprandexp !demand in stage t
	
	Neps:array(2..T) of integer !size of discrete state space of eps(t)
	Nfine,Ncoarse:array(1..T-1) of integer !# coarse and fine grid from each node in each stage
	GridSwitchLevel:array(1..T-1) of real !swithcing level from fine to coarse level
	
	M:array(2..T-1) of spvar !# of options bought in 1,to exercise in t
	m: array(2..T-1) of spvar !option exercised in t, effective in t+1
	Q:array(2..T) of spvar !fixed quantities ordered in 1 and delivered in t
	I:array(2..T) of spvar !net inventory at stage t
	I_p:array(2..T) of spvar !physical inventory at stage t
	I_m:array(2..T) of spvar !backorder at stage t
			
	R:splinctr !total revenue
	E:splinctr !total expenses
	P:splinctr !total profit
	NetInv:array(2..T) of splinctr !net inventory balance 
	InvBal:array(2..T) of splinctr !inventory balance across stages
	MaxNumOpt:array(2..T-1) of splinctr !bound on # of options bought
	MaxNumExOpt:array(2..T-1) of splinctr !bound on # of options exercised
end-declarations


setparam('xsp_implicit_stage',true) !assume last index set of array as stage
setparam('xsp_disp_warnings',false) !don't display warning messages
spsetstages(Stages) !set the stage set

!---------------Define demand process & generate tree ------------------------
forall(t in 2..T)
	if t>2 then D(t):=mu+rho*(D(t-1)-mu)+eps(t)*sigma*(1-rho^2)^.5
	else D(t):=mu+eps(t)*sigma
	end-if				  	

if Asymm then
	Nfine:=[3,3]
	Ncoarse:=[3,1]
	GridSwitchLevel:=[.5,1.5]
	GenAsymScenTree(Nfine,Ncoarse,GridSwitchLevel)
else
	Neps:=[41,31] !define discretized grid size
	GenSymScenTree(Neps) !generate scenario tree
end-if

!----------------------------Model formulation -------------------------------
forall(t in 2..T-1) spsetstage(M(t),1) !explicitly set to 1-st stage
forall(t in 2..T) spsetstage(Q(t),1) !explicitly set to 1-st stage
forall(t in 2..T) I(t) is_free !set as free variable

R:=r*(D(2)-I_m(2))+sum(t in 3..T) r*(D(t)+I_m(t-1)-I_m(t))+v*I_p(T)
E:=sum(t in 2..T-1) (e*m(t)+o*M(t))+sum(t in 2..T) (h*I_p(t)+s*I_m(t)+p*Q(t))
P:=R-E
forall(t in 2..T) NetInv(t):=I(t)=I_p(t)-I_m(t)
forall(t in 2..T)
	if (t>2) then InvBal(t):=I(t)=I(t-1)+Q(t)+m(t-1)-D(t);
	else InvBal(t):=I(t)=Q(t)-D(t)
	end-if
forall(t in 2..T-1) MaxNumExOpt(t):=m(t)<=M(t)
forall(t in 2..T-1) MaxNumOpt(t):=M(t)<=Mmax

!-----------------declarations for Chance and CVaR constraints ---------------
declarations
	z:spvar
	theta=5*10^(-6) !penalty
	
	ShortageProb,ShortageBound:splinctr
	beta=0.25 !max shortage prob
	
	S,CVaR:spvar
	c=6000
	alpha=.01
	CondLoss,CvarLoss,GlbBnd:splinctr		
end-declarations
!-----------------------------Chance constraints -----------------------------
if Chance and not SolveCVaR then
	spsetstage(z,T)
	z is_binary
	ShortageBound:=I_m(T)<=(sum(t in 2..T) D(t))*z
		!max shortage=D(2)+..+D(T)
	ShortageProb:=z<=beta
	spsettype(ShortageProb,"global")
	P-=theta*z
	
	!set xprs control parameters for better performance
	setparam("xprs_cutstrategy",0)
	setparam("xprs_treegomcuts",10000)
	setparam("xprs_miprelstop",0.01) 
end-if

!-------------------------------CVaR constraints -----------------------------
if SolveCVaR and not Chance then
	spsetstage(z,T) !by default z>=0
	S is_free !S belongs to stage 1 when implicit_stg is true
	CVaR is_free
	L:=E-R
	CvarLoss:=z>=L-S
	GlbBnd:=S+z/alpha<=CVaR 
	spsettype(GlbBnd,XSP_GLOBAL) !specify as global constraint
	CVaR<=c !bound
	L+=theta*CVaR !add this term to objective function
end-if

!--------------------------------Solve the problem ---------------------------
!generate extensive form with entities for each node in the event tree
setparam('xsp_scen_based',false)
if SolveCVaR and not Chance then
	minimize(L) !minimize the expected loss
else
	maximize(P) !maximize the expected profit
end-if

!---------------------------------- Procedures -------------------------------
procedure GenSymScenTree(N:array(range) of integer)
	Nmax:=max(t in 2..T) N(t) !max size of discretized state space
	declarations
		!store the discretized values and probabilities
		val,prob:array(1..Nmax) of real
	end-declarations
	forall(t in 2..T) do
		forall(n in 1..Nmax) do val(n):=0;prob(n):=0;end-do
		GetNormDist(N(t),val,prob) !get discretized val & probs
		spsetdist(eps(t),val,prob) !set them as distribution of eps(t)
	end-do	
	spgenexhtree !generate exhaustive tree based on the distribution
end-procedure

procedure GetNormDist(N:integer,x:array(range) of real,p:array(range) of real)
	(!
	if(not isodd(N) or N<3) then
		writeln("N must be an odd number and >=3")
		exit(1) !the cardinality of state space must be odd and >=3
	end-if
	!)
	delta:=6/N !delta element for distribution in (-3,3)
	p_:=0.0
	forall(n in 1..N) do
		x(n):=-3+(n-0.5)*delta !n-th discretized value
		!n-th discretized probability
		p(n):=(1/(2*M_PI)^0.5)*exp(-(x(n))^2/2)*delta
		p_+=p(n)
	end-do
	forall(n in 1..N) p(n):=p(n)/p_
end-procedure

procedure GenAsymScenTree(Nfine:array(range) of integer,
			Ncoarse:array(range) of integer,SwitchLevel:array(range) of real)
	Brmax:=max (t in 1..T-1) Nfine(t) !max # branches from a node
	Nmax:=integer(Brmax^(T-1)) !max # of nodes in a stage
	declarations
		N:array(1..T) of integer !# of nodes in each stage
		val,prob:array(1..Brmax) of real !assumed vals and probs by eps
		fine:dynamic array(1..T,1..Nmax) of boolean !is this a fine node
		UnconProb:dynamic array(1..T,1..Nmax) of real !unconditional prob
	end-declarations
	!initialize
	N(1):=1;UnconProb(1,1):=1;
	if(SwitchLevel(1)<=1) then fine(1,1):=true
	else fine(1,1):=false
	end-if
	!begin creating tree
	forall(t in 1..T-1) do
		forall(n in 1..N(t)) do
			if(fine(t,n)) then Br:=Nfine(t);else Br:=Ncoarse(t);end-if
			spaddchildren(t,n,Br) !create Br # branches from node(t,n)
			GetNormDist(Br,val,prob) !get realized vals and probs
			forall(b in 1..Br) do
				n_:=spgetchild(t,n,b)!node number in next stage
				spsetrandatnode(eps(t+1),n_,val(b)) !set val at node
				spsetprobcondatnode(t+1,n_,prob(b)) !set prob at node
				UnconProb(t+1,n_):=UnconProb(t,n)*prob(b) !update
			end-do
		end-do
		N(t+1):=spgetnodecount(t+1) !update # nodes in next stage			
		!update whether next stage nodes are fine or coarse
		if(t+1<T) then forall(n in 1..N(t+1)) 
			if(N(t+1)*UnconProb(t+1,n)>=SwitchLevel(t+1))
			then fine(t+1,n):=true
			else fine(t+1,n):=false
			end-if
		end-if
	end-do
	spgentree !generate the tree	
end-procedure

end-model