*****************************************************************
*** Chapter 8
*** Section: Batch Plant Design Under Uncertainty
*** Test Problem 1
***
*** Single-Product Campaign (2 products, 3 stages, 1 scenario)
*****************************************************************


SET
     i             products                                  /1*2/
     j             equipment                                 /1*3/
     q1            quadrature index                          /1*5/
     q2            quadrature index                          /1*5/
     p             periods                                   /1*1/;

****************************************************************************
* c = penalty parameter
* d, a(j), beta(j) = equipment cost parameters
* omega(p) = weighting parameter for scenario p
* pr(i) = price of product i
* v#(q#), W#(q#) = Gaussian quadrature parameters
* k(q1,q2) = joint probability for quadrature point
* S(i,j,p) = size parameter
* t(i,j,p) = processing time parameter
* H = planning time horizon
* theta(i,q1,q2) = demand for product i at quadrature point
* tlow(i,p) = time limiting step for product i in scenario p
* mean#, sigma# = statistics for demand distribution

PARAMETERS
d, c, a(j), beta(j), omega(p), pr(i),
v1(q1), v2(q2), W1(q1), W2(q2), k(q1,q2),
S(i,j,p), theta(i,q1,q2), H, t(i,j,p), tlow(i,p),
sigma1, sigma2, mean1, mean2,
t1max, t1min, t2max, t2min,
PI;

PI = 3.1415927;

mean1 = 200;   mean2 = 100;
sigma1 = 10;   sigma2 = 10;
t1max=mean1+4*sigma1;
t1min=mean1-4*sigma1;
t2max=mean2+4*sigma2;
t2min=mean2-4*sigma2;

**** Quadrature Points
v1('1')=-0.9061798459;   W1('1') = 0.23692885;
v1('2')=-0.5384693101;   W1('2') = 0.4786286705;
v1('3')=0;               W1('3') = 0.56888889;
v1('4')=0.5384693101;    W1('4') = 0.4786286705;
v1('5')=0.9061798459;    W1('5') = 0.23692885;
v2('1')=-0.9061798459;   W2('1') = 0.23692885;
v2('2')=-0.5384693101;   W2('2') = 0.4786286705;
v2('3')=0;               W2('3') = 0.56888889;
v2('4')=0.5384693101;    W2('4') = 0.4786286705;
v2('5')=0.9061798459;    W2('5') = 0.23692885;

theta('1',q1,q2) = 0.5*(t1max*(1+v1(q1))+t1min*(1-v1(q1)));
theta('2',q1,q2) = 0.5*(t2max*(1+v2(q2))+t2min*(1-v2(q2)));

k(q1,q2) = EXP(-0.5*(
		     ((theta('1',q1,q2)-mean1)*(theta('1',q1,q2)-mean1))
		     /(sigma1*sigma1)
		    +((theta('2',q1,q2)-mean2)*(theta('2',q1,q2)-mean2))
		     /(sigma2*sigma2)
		     )
		)/(2*PI*sigma1*sigma2);

d = 0.3;
c = 0;
a(j) = 10;
beta(j) = 0.6;
omega(p) = 1;
H = 8;

pr('1') = 5.5;   pr('2') = 7;

S('1','1','1') = 2; S('1','2','1') = 3; S('1','3','1') = 4;
S('2','1','1') = 4; S('2','2','1') = 6; S('2','3','1') = 3;

tlow('1','1') = LOG(20);  tlow('2','1') = LOG(16);

****************************************************************************
**** Variable bounds
PARAMETERS
vubd(j), vlbd(j),
blbd(i), bubd(i),
Qlbd(i,q1,q2,p), Qubd(i,q1,q2,p);

Qlbd('1',q1,q2,p) = t1min;
Qlbd('2',q1,q2,p) = t2min;
Qubd(i,q1,q2,p) = theta(i,q1,q2);

vubd(j) = LOG(4500);
vlbd(j) = LOG(500);

bubd('1') = LOG(4500/4);
bubd('2') = LOG(4500/6);
blbd('1') = LOG(500/4);
blbd('2') = LOG(500/6);

******************************************************************
* z = objective value
* v(j) = capacity of equipment j (volume)
* b(i) = batch size of produce i
* Q(i,q1,q2,p) = production of i at quadrature point (q1,q2)
*                for scenario p

VARIABLES z, v(j), Q(i,q1,q2,p), b(i);

******************************************************************
* obj = objective function
* con1 = equipment size constraint
* con2 = planning horizon constraint

EQUATIONS
         obj                        Convex objective function
	 con1(i,j,p)                Convex Inequality
	 con2(q1,q2,p)     Nonconvex Inequality;

*--------------------------------------------------------------------

obj..    z =e= d*SUM(j,a(j)*EXP(beta(j)*v(j)))
	        -SUM(p, SUM(q1, SUM(q2, SUM(i,
		 (c+1)*((t1max-t1min)/2)*((t2max-t2min)/2)
		      *omega(p)*W1(q1)*W2(q2)
		      *k(q1,q2)*pr(i)*Q(i,q1,q2,p)
		 -c*((t1max-t1min)/2)*((t2max-t2min)/2)
		      *omega(p)*W1(q1)*W2(q2)
		      *k(q1,q2)*pr(i)*theta(i,q1,q2)
		     )      )       )       );

con1(i,j,p)..   v(j) - b(i) =g= LOG(S(i,j,p));                    
con2(q1,q2,p)..  SUM(i, Q(i,q1,q2,p)*EXP(tlow(i,p)-b(i))) =l= H;


*-------------------------------------------------------

MODEL spc / all /;

*------------------------------------------------------