bilinear.gms : Convexification of bilinear term binary times x

Description

The model demonstrates various forumlations 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 roduction 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:


Reference

  • Kallrath, J, Combined Strategic Design and Operative Planning in the Process Industry, 2009. Submitted to Computers & Chemical Engineering

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 forumlations 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 roduction 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
$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 //

Sets 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;

Variables
     z        objective variable
     x(i)     production
     y(i)     waste material produced on special machine
     delta(i) indicator for production on special machine
Positive variables x, y;
Binary variable delta;

Equations
     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 /
model 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, Couenne, 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;

Equations
     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 /;

model 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);

Equations
     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
Equations
     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 Variables S1(i,j), S2(i,j);

Equations
     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);