(!******************************************************
   Xpress-SP Example
   ======================
   ``````````````````
   FILE: Airlift.mos
   
   DESCRIPTION: Stochastic 2-stage linear model
   
   TYPE: Transportation: Airlift operations scheduling
         
   MODEL:    	   	
   	The problem consists of scheduling airlift operations along 
   a set of routes during a period of one month. The demands of routes 
   can be predicted but are subject to changes. A random variable is 
   assigned for each demand. If the predicted demands do not agree with 
   the actual demands, a recourse action is taken. To meet the route 
   demands, several types of aircrafts are available. Each aircraft has 
   its own characteristics, namely availability of number of flight hours 
   per month and carrying capacity when flying a specific route. The 
   recourse actions that can be taken are to switch aircrafts from one 
   route to another, to leave some flights unused, and to contract commercial 
   flights to meet the demands of routes. 
   
   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 airlift
   uses "mmsp"
   
   parameters
   	DIR="Data/"
   	solve_mip=false !solve MIP or pure LP version of the problem
   	mipgap=0.01 ! if not solve_relax
   end-parameters

   
   
   declarations
       AIRCRAFTS = 2                                                    
       ROUTES    = 2                                                    
       SCEN      = 25                                                   

       A         : array(1..AIRCRAFTS,1..ROUTES) of integer             
       Aswitch   : array(1..AIRCRAFTS,1..ROUTES,1..ROUTES) of
                   integer 
       F         : array (1..AIRCRAFTS) of integer                      
       b         : array (1..AIRCRAFTS,1..ROUTES) of integer            
       Cost      : array (1..AIRCRAFTS,1..ROUTES) of real               
       CostSwitch: array (1..AIRCRAFTS,1..ROUTES,1..ROUTES) of real     
       cplus     : array (1..ROUTES) of real                            
       cminus    : array (1..ROUTES) of real                            
       Probabilities : array (1..SCEN) of real                             
       DValues   : array (1..ROUTES,1..SCEN) of real                    
   end-declarations

   declarations       
       d         : array (1..ROUTES) of sprand

       MaxHours  : array (1..AIRCRAFTS) of splinctr
       MaxSwitch : array (1..AIRCRAFTS, 1..ROUTES) of splinctr
       DemandCtr : array (1..ROUTES) of splinctr

       x         : array (1..AIRCRAFTS,1..ROUTES) of spvar            
       xswitch   : dynamic array (1..AIRCRAFTS,1..ROUTES,1..ROUTES) of spvar
       yplus     : array (1..ROUTES) of spvar                          
       yminus    : array (1..ROUTES) of spvar                          
   end-declarations
   
   ! Read data
   initializations from DIR+'airlift.dat'
      A Aswitch F b Cost CostSwitch
      cplus cminus DValues Probabilities
   end-initializations
   !---------------------------------------------------------
   ! Create scenario tree
   declarations
      Stages   = {"First","Recourse"}
      Branches : integer
      Values   : array (1..SCEN) of real
   end-declarations
   
   !setparam('xsp_loadnames',false) !uncomment to prevent exported files from having real names
   spsetstages(Stages)
   forall (i in 1..ROUTES) spsetstage(d(i),"Recourse")
   spcreatetree(SCEN)
   spsetprobcond(2,Probabilities)
   forall (i in 1..ROUTES, j in 1..SCEN)
      spsetrandatnode(d(i),j,DValues(i,j))
   spgentree
   !---------------------------------------------------------
   
   forall (i in 1..AIRCRAFTS, j in 1..ROUTES, k in 1..ROUTES | k<>j)
      create(xswitch(i,j,k))

   ExpTotalCost := 
      sum(i in 1..AIRCRAFTS, j in 1..ROUTES) Cost(i,j)*x(i,j) +
      sum(i in 1..AIRCRAFTS, j in 1..ROUTES, k in 1..ROUTES | k<>j)
         (CostSwitch(i,j,k) - Cost(i,j)*Aswitch(i,j,k)/A(i,j))*
           xswitch(i,j,k)+
      sum(j in 1..ROUTES) cplus(j)*yplus(j) +
      sum(j in 1..ROUTES) cminus(j)*yminus(j)
   
   ! First stage constraints:
   forall (i in 1..AIRCRAFTS)
      MaxHours(i) := sum(j in 1..ROUTES) A(i,j)*x(i,j) <= F(i)

   ! Second stage constraints:
   forall (i in 1..AIRCRAFTS, j in 1..ROUTES)
      MaxSwitch(i,j) := 
         sum(k in 1..ROUTES | k<>j) 
            Aswitch(i,j,k)*xswitch(i,j,k) <= A(i,j)*x(i,j)
   ! Demand constraint for recourse
   forall (j in 1..ROUTES)
      DemandCtr(j) := -sum(i in 1..AIRCRAFTS,k in 1..ROUTES | k<>j)
                   (b(i,j)*Aswitch(i,j,k)/A(i,j))*xswitch(i,j,k) +
                   sum(i in 1..AIRCRAFTS,k in 1..ROUTES | k<>j)
                      b(i,j)*xswitch(i,k,j) +
                   yplus(j) - yminus(j) =
                   d(j)- sum(i in 1..AIRCRAFTS) b(i,j)*x(i,j)
   ! Associate spvars to stages
   forall(i in 1..AIRCRAFTS,j in 1..ROUTES)
      spsetstage(x(i,j),'First')
   forall(i in 1..AIRCRAFTS,j in 1..ROUTES,k in 1..ROUTES | k<>j)
      spsetstage(xswitch(i,j,k),'Recourse')
   forall(j in 1..ROUTES) spsetstage(yplus(j),'Recourse')
   forall(j in 1..ROUTES) spsetstage(yminus(j),'Recourse')
   
   if solve_mip then   
	   ! Set variables as integer
	   forall(i in 1..AIRCRAFTS,j in 1..ROUTES) do
	      x(i,j) is_integer
	      forall(k in 1..ROUTES | exists(xswitch(i,j,k)))
	         xswitch(i,j,k) is_integer
	   end-do
	   setparam('xprs_miprelstop',mipgap)
   end-if
   
   spexportprob(XSP_MIN,ExpTotalCost,XSP_X_SMPS,"Airlift")
   minimize(ExpTotalCost)
   
   if solve_mip then
		if getparam('xprs_mipstatus') not in {XSP_MIP_SOLUTION,XSP_MIP_OPTIMAL} then
			exit(getparam('xprs_mipstatus'))
		end-if
   else
		if getparam('xprs_lpstatus') not in {XSP_LP_UNFINISHED,XSP_LP_OPTIMAL} then
			exit(getparam('xprs_lpstatus'))
		end-if
   end-if
   
   writeln(getobjval)

end-model