Linearize a nonlinear function using piecewise linear approximation

Frequently asked questions about GAMS

Moderator: aileen

Forum rules
Please ask questions in the other sub-forums
Locked
aileen
User
User
Posts: 136
Joined: 3 years ago

Linearize a nonlinear function using piecewise linear approximation

Post by aileen »

Find below some approaches to linearize the cubic function

Code: Select all

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.

Code: Select all

$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:

Code: Select all

----     13 VARIABLE y.L                   =       -2.113  
            VARIABLE x.L                   =        3.215 
Using binary variables to model a piecewise linear function

Code: Select all

$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

Code: Select all

$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

Code: Select all

$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

Code: Select all

$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:

Code: Select all

----     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
Locked