Find below some approaches to linearize the cubic function
y =e= x*sqr(x) - 7*sqr(x) + 14*x - 8
using piecewise linear approximation. In general, piecewise linear functions can be modeled using binary variables, SOS1 type formulations or SOS2 type formulations. Check textbooks like Nemhauser and Wolsey, Integer and Combinatorial Optimization (p. 10) and H.P. Williams, Model Building in Mathematical Programming (Sections 7.3, 9.3) for more information.
$offsymxref offsymlist
option limrow = 0,
limcol = 0,
solprint = off ;
variable y,x ;
equations eq ;
eq.. y =e= x*sqr(x) - 7*sqr(x) + 14*x - 8 ;
x.lo = 1 ; x.up = 10.5 ;
model cubic /all/ ;
solve cubic using nlp minimizing y ;
display y.l, x.l ;
The global solution to that problem is:
---- 13 VARIABLE y.L = -2.113
VARIABLE x.L = 3.215
Using binary variables to model a piecewise linear function
$ontext
Using Binary variables to model a piecewise linear function in GAMS.
This formulation uses points on the curve and is particularly good for developing approximations of
nonlinear functions
Ref: "Integer and Combinatorial Optimization" by
Nemhauser, G.L., and Wolsey, L.A., John Wiley
and Sons. Page 11.
$offtext
option limrow = 0,
limcol = 0
;
* this contains the set of points
set i /
1*20
/ ;
set notlast(i) ;
notlast(i)$(ord(i) lt card(i)) = yes ;
* this contains the values of x-coordinates
parameter xval(i) /
1 1
2 1.5
3 2
4 2.5
5 3
6 3.5
7 4
8 4.5
9 5
10 5.5
11 6
12 6.5
13 7
14 7.5
15 8
16 8.5
17 9
18 9.5
19 10
20 10.5
/ ;
* this contains the nonlinear function of x
parameter yval(i) ;
yval(i) = xval(i)*sqr(xval(i)) - 7*sqr(xval(i)) + 14*xval(i) - 8 ;
variables x,y,lam(i),bin(i) ;
positive variables lam(i) ;
binary variables bin(i) ;
equations xdef, ydef, norm, lamdef(i), soslam ;
xdef.. x =e= sum(i,lam(i)*xval(i)) ;
ydef.. y =e= sum(i,lam(i)*yval(i)) ;
norm.. sum(i,lam(i)) =e= 1 ;
lamdef(i).. lam(i) =l= bin(i-1) + bin(i)$notlast(i) ;
soslam.. sum(i$notlast(i),bin(i)) =e= 1 ;
x.lo = smin(i,xval(i)) ;
x.up = smax(i,xval(i)) ;
model piece /all/ ;
piece.optcr = 0.0 ;
solve piece minimizing y using mip ;
display x.l, y.l ;
display lam.l, bin.l ;
A SOS1 type formulation of a piecewise linear function
$ontext
The SOS1 type formulation of a piecewise linear function in GAMS.
This formulation uses points on the curve and is particularly good for developing approximations of nonlinear functions
Ref: "Integer and Combinatorial Optimization" by
Nemhauser, G.L., and Wolsey, L.A., John Wiley
and Sons. Page 11.
$offtext
option limrow = 0,
limcol = 0
;
* this contains the set of points
set i /
1*20
/ ;
set notlast(i) ;
notlast(i)$(ord(i) lt card(i)) = yes ;
* this contains the values of x-coordinates
parameter xval(i) /
1 1
2 1.5
3 2
4 2.5
5 3
6 3.5
7 4
8 4.5
9 5
10 5.5
11 6
12 6.5
13 7
14 7.5
15 8
16 8.5
17 9
18 9.5
19 10
20 10.5
/ ;
* this contains the nonlinear function of x
parameter yval(i) ;
yval(i) = xval(i)*sqr(xval(i)) - 7*sqr(xval(i)) + 14*xval(i) - 8 ;
variables x,y,lam(i),bin(i) ;
positive variables lam(i) ;
sos1 variables bin(i) ;
equations xdef, ydef, norm, lamdef(i), soslam ;
xdef.. x =e= sum(i,lam(i)*xval(i)) ;
ydef.. y =e= sum(i,lam(i)*yval(i)) ;
norm.. sum(i,lam(i)) =e= 1 ;
lamdef(i).. lam(i) =l= bin(i-1) + bin(i)$notlast(i) ;
soslam.. sum(i$notlast(i),bin(i)) =e= 1 ;
x.lo = smin(i,xval(i)) ;
x.up = smax(i,xval(i)) ;
model piece /all/ ;
piece.optcr = 0.0 ;
solve piece minimizing y using mip ;
display x.l, y.l ;
display lam.l, bin.l ;
A SOS2 type formulation of a piecewise linear function
$ontext
The SOS2 type formulation of a piecewise linear function in GAMS.
This formulation uses points on the curve and is particularly good for developing approximations of nonlinear functions
Ref: "Integer and Combinatorial Optimization" by
Nemhauser, G.L., and Wolsey, L.A., John Wiley
and Sons. Page 11.
$offtext
option limrow = 0
limcol = 0
;
* this contains the set of points
set i /
1*20
/ ;
set notlast(i) ;
notlast(i)$(ord(i) lt card(i)) = yes ;
* this contains the values of x-coordinates
parameter xval(i) /
1 1
2 1.5
3 2
4 2.5
5 3
6 3.5
7 4
8 4.5
9 5
10 5.5
11 6
12 6.5
13 7
14 7.5
15 8
16 8.5
17 9
18 9.5
19 10
20 10.5
/ ;
* this contains the nonlinear function of x
parameter yval(i) ;
yval(i) = xval(i)*sqr(xval(i)) - 7*sqr(xval(i)) + 14*xval(i) - 8 ;
variables x,y,lam(i) ;
sos2 variables lam(i) ;
equations xdef, ydef, norm ;
xdef.. x =e= sum(i,lam(i)*xval(i)) ;
ydef.. y =e= sum(i,lam(i)*yval(i)) ;
norm.. sum(i,lam(i)) =e= 1 ;
x.lo = smin(i,xval(i)) ;
x.up = smax(i,xval(i)) ;
model piece /all/ ;
piece.optcr = 0.0 ;
solve piece minimizing y using mip ;
display x.l, y.l ;
display lam.l ;
Another example of modeling piecewise linear cost functions using binary variables
$title piecewise linear cost functions
$ontext
Piecewise liner functions are used to model production cost.
We use simple binary variables. Production cost and volume data are presented using the
tulple cd with the following definition:
output.s0 minimum production
cost.s0 fixed cost
output.sN maximum total output for the Nth segment
cost.sN unit cost for this segment
Note product p04. Because of the minimum production we have to pay the fixed
charge and then try to minimize the loss.
$offtext
sets p products / p01*p05 /
s segments / s0*s4 /
i items / output, cost /
ps(p,s) product segment mapping
table cd(p,i,s) cost data
s0 s1 s2 s3 s4
p01.output 0 200 500 1000 4000
p01.cost 0 20 15 10 8
p02.output 0 50 100 500
p02.cost 0 10 15 20
p03.output 0 100 200
p03.cost 100 2 1
p04.output 0 100
p04.cost 100 2
p05.output 200 400 800
p05.cost 900 3 1
table market(p,*) market data
price min max
p01 16 3000
p02 12 inf
p03 5 inf
p04 3 50 80
p05 10 300 500
parameter fixed(p,s) fixed segment cost;
ps(p,s+1) = cd(p,'output',s) < cd(p,'output',s+1);
fixed(p,'s0') = cd(p,'cost','s0');
loop(ps(p,s), fixed(p,s) = fixed(p,s-1) + cd(p,'cost',s)*cd(p,'output',s-1) );
variables out(p) total output
seg(p,s) segment output
b(p,s) active segment
cost(p) tpotal product cost
profit
binary variable b; positive variables seg,out;
equations defout(p),defcost(p),defprofit,defseg(p,s),muex(p);
defprofit.. profit =e= sum(p, market(p,'price')*out(p) - cost(p));
defout(p).. out(p) =e= sum(ps(p,s), cd(p,'output',s-1)*b(ps) + seg(ps));
defseg(ps(p,s)).. seg(ps) =l= (cd(p,'output',s)-cd(p,'output',s-1))*b(ps);
defcost(p).. cost(p) =e= sum(ps(p,s), fixed(ps)*b(ps) + cd(p,'cost',s)*seg(ps));
muex(p).. sum(ps(p,s), b(ps)) =l= 1;
model sku / all /;
out.lo(p) = market(p,'min');
out.up(p) = market(p,'max');
sku.optcr = 0.001;
solve sku maximizing profit using mip;
parameter rep quick and dirty report;
rep(p,'sales') = out.l(p);
rep(p,'profit') = out.l(p)*market(p,'price') - cost.l(p);
display ps, fixed, rep;
Here are the results:
---- 85 SET ps product segment mapping
s1 s2 s3 s4
p01 YES YES YES YES
p02 YES YES YES
p03 YES YES
p04 YES
p05 YES YES
---- 85 PARAMETER fixed fixed segment cost
s0 s1 s2 s3 s4
p01 3000.000 8000.000 16000.000
p02 750.000 2750.000
p03 100.000 100.000 200.000
p04 100.000 100.000
p05 900.000 1500.000 1900.000
---- 85 PARAMETER rep quick and dirty report
sales profit
p01 3000.000 16000.000
p02 50.000 100.000
p03 300.000 1100.000
p04 80.000 -20.000
p05 500.000 3000.000