Description
This problem computes minimal cost solutions satisfying the
demand of pre-given product portfolios. It determines the number
and size of reactors and gives a schedule of how may batches of
each product run on each reactor. There are two scenarios (s1 20
products. s2 40 products), add --scenario s2 as a GAMS
parameter to specify the second scenario.
The global optimal reactor volumes are:
data set s1
vr.fx('r1') = 132.5; vr.fx('r2') = 250;
data set s2
vr.fx('r1') = 20; vr.fx('r2') = 100; vr.fx('r3') = 250;
Two formulations are presented, a compact MINLP formulations
and a linearized MIP formulation using special ordered sets.
Problem sizes for data set s1
MINLP MIP
variables 102 548
equations 28 418
Non-zeros 190 1574
discrete variables 40 186
Large Model of Type : MINLP
Category : GAMS Model library
Main file : kport.gms
$title Product Portfolio Optimization (KPORT,SEQ=270)
$onText
This problem computes minimal cost solutions satisfying the
demand of pre-given product portfolios. It determines the number
and size of reactors and gives a schedule of how may batches of
each product run on each reactor. There are two scenarios (s1 20
products. s2 40 products), add --scenario s2 as a GAMS
parameter to specify the second scenario.
The global optimal reactor volumes are:
data set s1
vr.fx('r1') = 132.5; vr.fx('r2') = 250;
data set s2
vr.fx('r1') = 20; vr.fx('r2') = 100; vr.fx('r3') = 250;
Two formulations are presented, a compact MINLP formulations
and a linearized MIP formulation using special ordered sets.
Problem sizes for data set s1
MINLP MIP
variables 102 548
equations 28 418
Non-zeros 190 1574
discrete variables 40 186
Kallrath, J. Exact Computation of Global Minima of a Nonconvex
Portfolio Optimization Problem. In Frontiers in Global Optimization.
Eds Floudas C A and Pardalos P M. Kluwer Academic Publishers,
Dortrecht, 2003.
Keywords: mixed integer nonlinear programming, mixed integer linear programming,
portfolio optimization, special ordered sets, global optimization,
concave objective functions, chemical engineering
$offText
$eolCom //
Set
s 'scenario' / s1,s2 /
rR 'reactors' / R1*R3 /
pP 'products' / L1*L37 /
r(rR) 'reactors considered in scenario'
p(pP) 'products considered in scenario';
Table RData(rR,s,*) 'reactor data'
s1.VMIN s1.VMAX s2.VMIN s2.VMAX
R1 102.14 250 20 50
R2 176.07 250 52.5 250
R3 151.25 250;
Table PData(pP,s,*) 'product data'
s1.Dem s1.PTime s2.Dem s2.Ptime
L1 2600 6 2600 6
L2 2300 6 2300 6
L3 1700 6 450 6
L4 530 6 1200 6
L5 530 6 560 6
L6 280 6 530 6
L7 250 6 530 6
L8 230 6 140 6
L9 160 6 110 6
L10 90 6 110 6
L11 70 6 10 6
L12 390 6 110 6
L13 250 6 90 6
L14 160 6 90 6
L15 100 6 90 6
L16 70 6 70 6
L17 50 6 50 6
L18 50 6 30 6
L19 50 6 10 6
L20 10 6
L21 10 6
L22 190 6
L23 180 6
L24 70 6
L25 70 6
L26 40 6
L27 40 6
L28 40 6
L29 30 6
L30 20 6
L31 20 6
L32 20 6
L33 10 6
L34 10 6
L35 10 6
L36 10 6
L37 10 6;
Parameter
VMIN(rR) 'volume flow of products in m^3 per week'
VMAX(rR) 'volume flow of products in m^3 per week'
DEMAND(pP) 'volume flow of products in m^3 per week'
PRODTIME(pP) 'production time in hours per batch';
Scalar
WHRS 'hours in a week' / 168 /
CSTI 'in kEuro depreciation per m^3 reactor and week' / 0.97 /
CSTF 'in kEuro per week and reactor' / 2.45 /
ESF 'economies of scale factor' / 0.5 /;
$if not set scenario $set scenario s1
VMIN(rR) = RDATA(rR,'%scenario%','VMIN');
VMAX(rR) = RDATA(rR,'%scenario%','VMAX');
DEMAND(pP) = PDATA(pP,'%scenario%','Dem');
PRODTIME(pP) = PDATA(pP,'%scenario%','PTime');
* determine scenario sets
r(rR) = VMAX(rR) > 0;
p(pP) = DEMAND(pP) > 0;
* definition of compact MINLP model
Variable
cTotal 'total costs'
cInvest 'invest cost'
cFixed 'fix costs'
f(rR,pP) 'utilization rate'
vR(rR) 'reactor volume in m^3'
pS(pP) 'surplus production'
bvr(rR) 'indicating whether reactor r is active'
nB(rR,pP) 'number of batches of product p in reactor r';
Positive Variable f, vR, pS;
Integer Variable nB;
Binary Variable bvr;
Equation
DEFcT 'total costs'
DEFcF 'fix costs'
DEFcI 'invest cost'
TR(rR) 'production time of reactor r'
SPP(pP) 'compute surplus production p'
RVUB(rR) 'maximal volume of reactor r'
RVLB(rR) 'minimal volume reactor r';
* define the total cost
DEFcT.. cTotal =e= cFixed + cInvest;
DEFcF.. cFixed =e= sum(r, CSTF*bvr(r));
DEFcI.. cInvest =e= sum(r, CSTI**ESF*vR(r)**ESF);
* limit the total production time of reactor r
TR(r).. sum(p, PRODTIME(p)*nB(r,p)) =l= WHRS*bvr(r);
* compute the surplus production
SPP(p).. pS(p) =e= sum(r, nB(r,p)*f(r,p)*vR(r))/DEMAND(p) - 1;
* lower and upper bounds on reactor volume
RVLB(r).. vR(r) =g= VMIN(r)*bvr(r);
RVUB(r).. vR(r) =l= VMAX(r)*bvr(r);
Model portfolioMINLP / DEFcT, DEFcF, DEFcI, TR, SPP, RVLB, RVUB /;
f.lo(r,p) = 0.4; f.up(r,p) = 1; // bounds on the utilization rates
pS.lo(p) = 0; pS.up(p) = 1; // bounds on the surplus production
* bounds on the number of batches
nB.lo(r,p) = 0;
nB.up(r,p) = min(WHRS/PRODTIME(p),floor(2*DEMAND(p)/(VMIN(r)*f.lo(r,p))));
nB.up(r,p)$(2*DEMAND(p) < f.lo(r,p)*VMIN(r)) = 0;
vR.l(rR) = 99;
vR.lo(r) = VMIN(r);
solve portfolioMINLP using minlp minimizing cTotal;
* additional variables and equations to define the MIP formulation
* first we need to linearize the product terms:
Set
i 'dyadic represenation set' / 0*10 /
j 'discretization points for SOS2' / 0*10 /
rpi(rR,pP,i) 'i required for representing np';
Parameter
vRj(rR,j) 'x part of SOS2 construct'
ESFvRj(rR,j) 'y part of SOS2 construct';
Positive Variable
pT(rR,pP) 'number of batches x reactor volume in m^3'
pT2(rR,pP,i) 'same for in dyadic representation'
ESFvR(rR) 'economies of scale for vR';
Binary Variable
nBx(rR,pP,i) 'dyadic represenation of nB';
SOS2 Variable
lambda(rR,j) 'approximation of economies of scale function';
Equation
CNP(rR,pP) 'compute the nonlinear products nB(rp)*f(rp)*vR(r)'
SPPx(pP) 'compute surplus production p'
CNPl0(rR,pP) 'linearized version of CNP'
CNPl1(rR,pP) 'linearized version of CNP'
CNPl2(rR,pP,i) 'linearized version of CNP'
CNPl3(rR,pP,i) 'linearized version of CNP'
CNPl4(rR,pP,i) 'linearized version of CNP'
DEFSOSx(rR) 'SOS2 x construct'
DEFSOSy(rR) 'SOS2 y construct'
DEFSOSone(rR) 'SOS2 sum construct'
DEFcIlp 'linearized version of DEFcI';
* new surplus production equation
SPPx(p).. pS(p) =e= sum(r, pT(r,p))/DEMAND(p) - 1;
* computes batches x volume
CNP(r,p).. pT(r,p) =e= nB(r,p)*f(r,p)*vR(r);
* Linearized version of CNP
CNPl0(r,p).. pT(r,p) =e= sum(rpi(r,p,i),2**(ord(i) - 1)*pT2(rpi));
CNPl1(r,p).. nB(r,p) =e= sum(rpi(r,p,i),2**(ord(i) - 1)*nBx(rpi));
CNPl2(rpi(r,p,i)).. pT2(rpi) =l= VMAX(r)*nBx(rpi);
CNPl3(rpi(r,p,i)).. pT2(rpi) =l= vR(r);
CNPl4(rpi(r,p,i)).. pT2(rpi) =g= f.lo(r,p)*(vR(r) - VMAX(r)*(1 - nbx(rpi)));
* SOS2 approximation of economies of scale function
DEFSOSx(r).. vR(r) =e= sum(j, vRj(r,j)*lambda(r,j));
DEFSOSy(r).. ESFvR(r) =e= sum(j, ESFvRj(r,j)*lambda(r,j));
DEFSOSone(r).. sum(j, lambda(r,j)) =e= 1;
DEFcIlp.. cInvest =e= sum(r, CSTI**ESF*ESFvR(r));
rpi(r,p,i) = ord(i) <= ceil(log(max(1,nB.up(r,p)))/log(2)) + 1$nB.up(r,p);
vRj(r,j) = (VMAX(r) - VMIN(r))*(ord(j) - 1)/(card(j) - 1) + VMIN(r);
ESFvRj(r,j) = vRj(r,j)**ESF;
Model PortfolioMIP / TR, SPPx, RVLB, RVUB, DEFcF
DEFcIlp, DEFcT, CNPl0, CNPl1, CNPl2
CNPl3, CNPl4, DEFSOSx, DEFSOSy, DEFSOSone /;
portfolioMIP.optCr = .05;
solve portfolioMIP using mip minimizing cTotal;