Description
The model demonstrates various formulations to represent bilinear
product terms of one continuous and one binary variable.
A set of 60 products i is produced on a set of machine with a given
total capacity. Some machine are special in the sense that if a
product is produced on one of them, cleaning treatment costs apply
caused by a set of cleaning treatment machines t.
A binary variable, delta(i), indicates that product i is produced on
one of the special machines. The model is simplified regarding the
machine-product relations.
Here we mimic a larger production problem, and just require that
E1.. sum(iE, delta(iE)) =e= 2;
E2.. sum(iO, delta(iO)) =e= 5;
which represents the fact that it cannot be avoided to use the special
machine and, thus, to have some cleaning treatment.
If product i is produced on a special machine, then the amount, y(i),
of the by-product is given by the recipe constraint y(i)=0.164*p(i),
where the non-negative variable p(i) is the amount produced on special
machines. For each product there is a specific yield of YS(i) $/ton.
The by-product is burnt and leads to an energy yield of YB(i) $/ton,
where YB(i)<YS(i). The by-product also passes the treatment plant.
The production is limited by the production capacity C, where x(i),
100+i <= x(i) <= XUB, is the amount of product i produced.
The amount produced on special machines is p(i)=x(i)*delta(i).
We compare the non-convex MINLP formulation to equivalent linear
forms of p(i)=x(i)*delta(i) using big-M, convex hull, and indicator
forumlations. Moreover, a new SOS-1 formulation is presented which is
described in:
Large Model of Type : MINLP
Category : GAMS Model library
Main file : bilinear.gms
$title Convexification of bilinear term binary times x (BILINEAR,SEQ=346)
$onText
The model demonstrates various formulations to represent bilinear
product terms of one continuous and one binary variable.
A set of 60 products i is produced on a set of machine with a given
total capacity. Some machine are special in the sense that if a
product is produced on one of them, cleaning treatment costs apply
caused by a set of cleaning treatment machines t.
A binary variable, delta(i), indicates that product i is produced on
one of the special machines. The model is simplified regarding the
machine-product relations.
Here we mimic a larger production problem, and just require that
E1.. sum(iE, delta(iE)) =e= 2;
E2.. sum(iO, delta(iO)) =e= 5;
which represents the fact that it cannot be avoided to use the special
machine and, thus, to have some cleaning treatment.
If product i is produced on a special machine, then the amount, y(i),
of the by-product is given by the recipe constraint y(i)=0.164*p(i),
where the non-negative variable p(i) is the amount produced on special
machines. For each product there is a specific yield of YS(i) $/ton.
The by-product is burnt and leads to an energy yield of YB(i) $/ton,
where YB(i)<YS(i). The by-product also passes the treatment plant.
The production is limited by the production capacity C, where x(i),
100+i <= x(i) <= XUB, is the amount of product i produced.
The amount produced on special machines is p(i)=x(i)*delta(i).
We compare the non-convex MINLP formulation to equivalent linear
forms of p(i)=x(i)*delta(i) using big-M, convex hull, and indicator
forumlations. Moreover, a new SOS-1 formulation is presented which is
described in:
Kallrath, J, Combined Strategic Design and Operative Planning in the
Process Industry, 2009. Submitted to Computers & Chemical Engineering
Keywords: mixed integer nonlinear programming, mixed integer quadratic constraint
programming, extended mathematical programming, special ordered sets,
mathematics, production planning, modeling techniques, indicator constraints
$offText
$if not set solveNC $set solveNC 1
$if not set solvebigM1 $set solvebigM1 1
$if not set solvebigM2 $set solvebigM2 0
$if not set solveIndic $set solveIndic 0
$if not set solveEMPCH $set solveEMPCH 0
$if not set solveEMPI $set solveEMPI 0
$if not set solveEMPBM1 $set solveEMPBM1 0
$if not set solveEMPBM2 $set solveEMPBM2 0
$if not set solveSOS1 $set solveSOS1 0
* Modell dimensions
$if not set MaxI $set MaxI 60
$if not set MaxT $set MaxT 10
$eolCom //
Set
i 'products to be produced and sold' / i1*i%MaxI% /
iE(i) 'products with even ordinal number'
iO(i) 'products with odd ordinal number'
t 'cleaning treatment facilities' / t1*t%MaxT% /;
iE(i) = mod(ord(i),2) = 0;
iO(i) = not iE(i);
Parameter
Capacity 'total machine capacity' / 20000 /
C(i,t) 'cleaning treatment costs'
XUB(i) 'upper bound on production'
XLB(i) 'lower bound on production'
YS(i) 'yield from selling product i'
YB(i) 'yield from burning extra waste';
C(i,t) = sqrt(ord(i))*ord(t);
C(iE,t) = -C(iE,t) + 5;
XUB(i) = 10000;
XLB(i) = 100 + ord(i);
YS(i) = 0.04 + 0.001*sqrt(ord(i));
YB(i) = 0.007;
Variable
z 'objective variable'
x(i) 'production'
y(i) 'waste material produced on special machine'
delta(i) 'indicator for production on special machine';
Positive Variable x, y;
Binary Variable delta;
Equation
E1, E2 'force use some of the special machine'
ByProductNC(i) 'by-product produced on special machine'
ProdCap 'production capacity'
ObjFuncNC 'objective function: yield minus cleaning treatment costs';
ObjFuncNC.. z =e= sum(i, YS(i)*x(i) + YB(i)*y(i))
- sum(t, sqr(sum(i, C(i,t)*x(i)*delta(i) + y(i))));
ProdCap.. sum(i, x(i)) =l= Capacity;
ByProductNC(i).. y(i) =e= 0.164*x(i)*delta(i);
E1.. sum(iE, delta(iE)) =e= 2;
E2.. sum(iO, delta(iO)) =e= 5;
Model
core / ProdCap, E1, E2 /
NC 'non-convex model' / core, ByProductNC, ObjFuncNC /;
x.lo(i) = XLB(i);
x.up(i) = XUB(i);
* We need a global solver to find optimum of non-convex model
* Solver alternatives: Baron, LindoGlobal, SCIP
option miqcp = cplex, optCr = 0;
NC.workFactor = 10;
if(%solveNC%, solve NC max z using minlp;);
* First bigM Convexification
Positive Variable
p(i) 'product x times delta';
Equation
ByProduct(i) 'by-product produced on special machine'
ObjFunc 'objective function: yield minus cleaning treatment costs'
bigM1_1, bigM1_2, bigM1_3 'bigM convexification of binary times bounded continuous';
ByProduct(i).. y(i) =e= 0.164*p(i);
ObjFunc.. z =e= sum(i, YS(i)*x(i) + YB(i)*y(i))
- sum(t, sqr(sum(i, C(i,t)*p(i) + y(i))));
bigM1_1(i).. p(i) =l= x(i); // this is not needed because of the sign of p in the objective
bigM1_2(i).. p(i) =l= XUB(i)*delta(i);
bigM1_3(i).. p(i) =g= x(i) - XUB(i)*(1 - delta(i));
Model
coreConv / core, ByProduct, ObjFunc /
bigM1 / coreConv, bigM1_1, bigM1_2, bigM1_3 /;
p.up(i) = XUB(i);
$onEcho > cplex.opt
mipEmphasis 3
$offEcho
if(%solvebigM1%, bigM1.optFile = 1; solve bigM1 max z using miqcp;);
* Alternative bigM forumulation
Positive Variable slack(i);
Equation bigM2_1, bigM2_2, bigM2_3 'bigM convexification of binary times bounded continuous';
bigM2_1(i).. p(i) =e= x(i) - slack(i);
bigM2_2(i).. p(i) =l= XUB(i)*delta(i); // this is not needed because of the sign of p in the objective
bigM2_3(i).. slack(i) =l= XUB(i)*(1 - delta(i));
Model bigM2 / coreConv, bigM2_1, bigM2_2, bigM2_3 /;
slack.up(i) = XUB(i);
if(%solvebigM2%, bigM2.optFile = 1; solve bigM2 max z using miqcp;);
* Cplex Indicator Formulation
Equation disj1, disj2 'indicator convexification of binary times bounded continuous';
disj1(i).. p(i) =e= x(i);
disj2(i).. p(i) =e= 0; // this is not needed because of the sign of p in the objective
Model indic / coreConv, disj1, disj2 /;
$onEcho > cplex.op2
indic disj1(i)$delta(i) 1
indic disj2(i)$delta(i) 0
cuts 3
$offEcho
if(%solveIndic%, indic.optFile = 2; solve indic max z using miqcp;);
* The EMP (Extended Math Programming) framework explores modeling
* extensions that result in non-traditional math programs (like
* disjunctions) and automate the reformulation into traditional math
* programs (like MIPs). The manually generated big-M and indicator
* formulations above are automatically produced by EMP from a model
* with disjunctions. Moreover, EMP provides a convex hull formulation
* (which is independent of a bigM) for disjunctions.
* EMP Formulations
File femp / "%emp.info%" /;
put femp;
$onEcho > jams.opt
SubSolver cplex
SubSolverOpt 1
$offEcho
* Convex Hull Convexification
putClose 'modeltype miqcp disjunction delta disj1 else disj2';
if(%solveEMPCH%, indic.optFile = 1; solve indic max z using emp;);
* Cplex Indicator Convexification
putClose 'modeltype miqcp disjunction indic delta disj1 else disj2';
if(%solveEMPI%, indic.optFile = 1; solve indic max z using emp;);
* Big-M Convexification type 1 (similar to bigM1 formulation)
put 'modeltype miqcp';
loop(i, put / 'disjunction bigM' XUB(i) delta(i) disj1(i) 'else' disj2(i));
putClose;
if(%solveEMPBM1%, indic.optFile = 1; solve indic max z using emp;);
* Big-M Convexification type 2 (similar to bigM2 forumlation)
put 'modeltype miqcp';
loop(i, put / 'disjunction bigM' XUB(i) 1e-4 1 delta(i) disj1(i) 'else' disj2(i));
putClose;
if(%solveEMPBM2%, indic.optFile = 1; solve indic max z using emp;);
* SOS1 Formulation
delta.prior(i) = inf; // relax binary requirement of delta
Set j 'binary choice' / 0, 1 /;
SOS1 Variable S1(i,j), S2(i,j);
Equation defS1_0, defS1_1, defS2_0, defS2_1 'selection constraints';
defS1_0(i).. S1(i,'0') =e= delta(i);
defS1_1(i).. S1(i,'1') =e= x(i) - p(i);
defS2_0(i).. S2(i,'0') =e= 1 - delta(i);
defS2_1(i).. S2(i,'1') =e= p(i);
Model sos1conv / coreConv, defS1_0, defS1_1, defS2_0, defS2_1 /;
if(%solveSOS1%, solve sos1conv max z using miqcp;);