Piecewise linear functions and formulations for interpolation (part 2)

Pyomo supports a number of methods to perform piecewise linear interpolation. In part 1 [1], the following methods for piecewise linear interpolation were discussed:

1. SOS2: use SOS2 variables to model the piecewise linear functions. This is an easy modeling exercise.
2. BIGM_BIN: use binary variables and big-M constraints to enable and disable constraints. This is a more complex undertaking, especially if we want to use the smallest possible values for the big-M constants. Pyomo has bug in this version.
3. BIGM_SOS1: a slight variation of the BIGM_BIN model. Pyomo has the same problem with this model.

Here I show a few more formulations:

4. DCC: Disaggregated Convex Combination formulation is a simulation of the SOS2 model by binary variables. 
5. CC:  Convex Combination formulation is a similar SOS2 like approach using binary variables only.
6. MC: Multiple Choice model uses a semi-continuous approach.
7. INCR: Incremental formulation is using all previous segments.

I'll state the formulation in mathematical notation, hopefully a little bit more accessible than in some papers I looked at, and illustrate the formulation with GAMS code to make it more concrete.

Note that part 3 [2] will discuss some newer formulations: LOG and DLOG that only need a logarithmic number of binary variables.

For demonstration purposes we use the same small example with 4 breakpoints and 3 segments.

Method 4: DCC - Disaggregated Convex Combination Model

This is an impressive name for a model that is actually not very complicated. 

First we define a binary variable for each segment: $$\delta_s = \begin{cases} 1 & \text{if segment $s$ is selected}\\ 0 & \text{otherwise}\end{cases}$$ Then we perform interpolation on the selected segment. For this we define a weight variable: $$\lambda_{s,k} \ge 0$$ for the two end points belonging to segment $s$.  This way we achieve the "only two neighboring weights $\lambda$ can be nonzero" constraint from our SOS2 model. The mathematical model can look like:

DCC - Disaggregated Convex Combination Formulation

$$\begin{align} & \color{DarkRed}x = \sum_{s,k|\color{DarkBlue}SK(s,k)} \color{DarkBlue}{\bar{x}}_k \color{DarkRed}\lambda_{s,k} \\ & \color{DarkRed}y = \sum_{s,k|\color{DarkBlue}SK(s,k)} \color{DarkBlue}{\bar{y}}_k \color{DarkRed}\lambda_{s,k}\\& \color{DarkRed}\delta_s = \sum_{k|\color{DarkBlue}SK(s,k)} \color{DarkRed} \lambda_{s,k} \\& \sum_s \color{DarkRed}\delta_s=1 \\ & \color{DarkRed}\lambda_{s,k} \ge 0 \\ & \color{DarkRed}\delta_s \in \{0,1\} \end{align}$$

Here the mapping set $SK(s,k)$ is TRUE for the two breakpoints $k$ belonging to segment $s$.

We can also write the model more directly, without the mapping $SK(s,k)$:

DCC - Disaggregated Convex Combination Formulation 2

$$\begin{align} & \color{DarkRed}x = \sum_s \left( \color{DarkBlue}{\bar{x}}_s \color{DarkRed}\lambda_{s,s} + \color{DarkBlue}{\bar{x}}_{s+1} \color{DarkRed}\lambda_{s,s+1} \right) \\ & \color{DarkRed}y = \sum_s \left(  \color{DarkBlue}{\bar{y}}_s \color{DarkRed}\lambda_{s,s}+\color{DarkBlue}{\bar{y}}_{s+1} \color{DarkRed}\lambda_{s,s+1}\right)\\& \color{DarkRed}\delta_s =  \color{DarkRed} \lambda_{s,s}+ \color{DarkRed}\lambda_{s,s+1} \\& \sum_s \color{DarkRed}\delta_s=1 \\ & \color{DarkRed}\lambda_{s,k} \ge 0 \\ & \color{DarkRed}\delta_s \in \{0,1\} \end{align}$$

Here is the GAMS version:

set   k 'breakpoints' /point1*point4/   s 'segments' /segment1*segment3/   sk(s,k) 'mapping' ; sk(s,k) = ord(s)=ord(k) or ord(s)=ord(k)-1; display sk; table data(k,*)            x   y   point1   1   6   point2   3   2   point3   6   8   point4  10   7 ; positive variable lambda(s,k) 'interpolation'; binary variable delta(s) 'select segment'; variable x,y; equations    xdef       'x'    ydef       'y'    link(s)    'link delta-lambda'    sumdelta   'select segment' ; xdef.. x =e= sum(sk(s,k), lambda(s,k)*data(k,'x')); ydef.. y =e= sum(sk(s,k), lambda(s,k)*data(k,'y')); link(s).. delta(s) =e= sum(sk(s,k),lambda(s,k)); sumdelta.. sum(s, delta(s)) =e= 1; x.fx = 5; model m /all/; option optcr=0; solve m maximizing y using mip; display x.l,y.l,delta.l,lambda.l;

In this model I made a strong distinction between segments $s$ and breakpoints $k$. This will help the GAMS model to perform domain checking (type checking), so we get a bit better protection against errors in the model. Of course, if you prefer you just can indicate segments by $1,\dots,K-1$.

Note that the linking constraints link, perform two functions. First: they make sure that for unselected segments (with $\delta_s=0$), we have $\lambda_{s,k}=0$.  Second, for the selected segment, we automatically sum the $\lambda$'s to 1.

The output of the model looks like:

----     13 SET sk  mapping

              point1      point2      point3      point4

segment1         YES         YES
segment2                     YES         YES
segment3                                 YES         YES

----     45 VARIABLE x.L                   =        5.000  
            VARIABLE y.L                   =        6.000  

----     45 VARIABLE delta.L  

segment2 1.000


----     45 VARIABLE lambda.L  

              point2      point3

segment2       0.333       0.667

The display of the set SK confirms we setup the topography correctly: segment $k$ uses points $k, k+1$.

We fixed $x=5$ which yields $y=6$. We see $\delta_2=1$ so the second segment is selected, and this is also visible from the variables $\lambda_{s,k}$ indicating we interpolate between points 2 and 3.

We can have a look at the generated equations:

---- xdef  =E=  

xdef..  - lambda(segment1,point1) - 3*lambda(segment1,point2) - 3*lambda(segment2,point2) - 6*lambda(segment2,point3)
     
      - 6*lambda(segment3,point3) - 10*lambda(segment3,point4) + x =E= 0 ; (LHS = 5, INFES = 5 ****)
     

---- ydef  =E=  

ydef..  - 6*lambda(segment1,point1) - 2*lambda(segment1,point2) - 2*lambda(segment2,point2) - 8*lambda(segment2,point3)
     
      - 8*lambda(segment3,point3) - 7*lambda(segment3,point4) + y =E= 0 ; (LHS = 0)
     

---- link  =E=  

link(segment1)..  - lambda(segment1,point1) - lambda(segment1,point2) + delta(segment1) =E= 0 ; (LHS = 0)
     
link(segment2)..  - lambda(segment2,point2) - lambda(segment2,point3) + delta(segment2) =E= 0 ; (LHS = 0)
     
link(segment3)..  - lambda(segment3,point3) - lambda(segment3,point4) + delta(segment3) =E= 0 ; (LHS = 0)
     

---- sumdelta  =E=  

sumdelta..  delta(segment1) + delta(segment2) + delta(segment3) =E= 1 ; (LHS = 0, INFES = 1 ****)

In a sense we simulated the SOS2 model using binary variables. More visible here is the two-step approach:

1. Select the segment $s$ using the binary variables $\delta_s$.
2. Interpolate between the breakpoints of this segment using the positive variables $\lambda_{s,k}$.

Of course in a MIP model we have simultaneous equations, so these things happen actually at the same time. The 2-step paradigm is more of a useful mental model.

For completeness, here is the GAMS model without set SK(s,k):

* DCC2: Disaggregated Convex Combination formulation 2 * without mapping SK(s,k) set   k 'breakpoints' /point1*point4/   s(k) 'segments' /point1*point3/ ; table data(k,*)            x   y   point1   1   6   point2   3   2   point3   6   8   point4  10   7 ; positive variable lambda(s,k) 'interpolation'; binary variable delta(s) 'select segment'; variable x,y; equations    xdef       'x'    ydef       'y'    link(s)    'link delta-lambda'    sumdelta   'select segment' ; xdef.. x =e= sum(s(k), lambda(s,k)*data(k,'x')+lambda(s,k+1)*data(k+1,'x')); ydef.. y =e= sum(s(k), lambda(s,k)*data(k,'y')+lambda(s,k+1)*data(k+1,'y')); link(s(k)).. delta(s) =e= lambda(s,k)+lambda(s,k+1); sumdelta.. sum(s, delta(s)) =e= 1; x.fx = 5; model m /all/; option optcr=0; solve m maximizing y using mip; display x.l,y.l,delta.l,lambda.l;

Method 5: CC - Convex Combination Model

This formulation is very similar to the previous one. The main difference is the structure of the $\lambda$ variables. We index them by $k$ only. This is more like the SOS2 model.  The manner in which we link $\lambda_k$ to $\delta_s$ becomes somewhat different. Instead of an equality we use an inequality $$\lambda_k \le \sum_{s|SK(s,k)}\delta_s$$ This make sure that only for a selected segment with $\delta_s=1$, the corresponding $\lambda_k$'s can be nonzero. We need to add explicitly that $\sum_k \lambda_k=1$ as this is no longer implied by the linking constraints. The model looks like:

CC - Convex Combination Formulation

$$\begin{align} & \color{DarkRed}x = \sum_k \color{DarkBlue}{\bar{x}}_k \color{DarkRed}\lambda_{k} \\ & \color{DarkRed}y = \sum_k \color{DarkBlue}{\bar{y}}_k \color{DarkRed}\lambda_{k} \\ & \sum_k \color{DarkRed}\lambda_k = 1 \\& \color{DarkRed} \lambda_k \le \sum_{s|\color{DarkBlue}SK(s,k)} \color{DarkRed}\delta_s \\& \sum_s \color{DarkRed}\delta_s=1 \\ & \color{DarkRed}\lambda_k \ge 0 \\ & \color{DarkRed}\delta_s \in \{0,1\} \end{align}$$

This model has fewer continuous variables than the DCC model. The GAMS version looks like:

set   k 'breakpoints' /point1*point4/   s 'segments' /segment1*segment3/   sk(s,k) 'mapping' ; sk(s,k) = ord(s)=ord(k) or ord(s)=ord(k)-1; display sk; table data(k,*)            x   y   point1   1   6   point2   3   2   point3   6   8   point4  10   7 ; positive variable lambda(k) 'interpolation'; binary variable delta(s) 'select segment'; variable x,y; equations    xdef       'x'    ydef       'y'    link(k)    'link delta-lambda'    sumlambda  'interpolation'    sumdelta   'select segment' ; xdef.. x =e= sum(k, lambda(k)*data(k,'x')); ydef.. y =e= sum(k, lambda(k)*data(k,'y')); link(k).. lambda(k) =l= sum(sk(s,k),delta(s)); sumlambda.. sum(k, lambda(k)) =e= 1; sumdelta.. sum(s, delta(s)) =e= 1; x.fx = 5; model m /all/; option optcr=0; solve m maximizing y using mip; display x.l,y.l,delta.l,lambda.l;

The output (including the generated equations) is:

Generated equations

---- xdef  =E=  x

xdef..  - lambda(point1) - 3*lambda(point2) - 6*lambda(point3) - 10*lambda(point4) + x =E= 0 ; (LHS = 5, INFES = 5 ****)
     

---- ydef  =E=  y

ydef..  - 6*lambda(point1) - 2*lambda(point2) - 8*lambda(point3) - 7*lambda(point4) + y =E= 0 ; (LHS = 0)
     

---- link  =L=  link delta-lambda

link(point1)..  lambda(point1) - delta(segment1) =L= 0 ; (LHS = 0)
     
link(point2)..  lambda(point2) - delta(segment1) - delta(segment2) =L= 0 ; (LHS = 0)
     
link(point3)..  lambda(point3) - delta(segment2) - delta(segment3) =L= 0 ; (LHS = 0)
     
link(point4)..  lambda(point4) - delta(segment3) =L= 0 ; (LHS = 0)
     

Solution

---- sumlambda  =E=  interpolation

sumlambda..  lambda(point1) + lambda(point2) + lambda(point3) + lambda(point4) =E= 1 ; (LHS = 0, INFES = 1 ****)
     

---- sumdelta  =E=  select segment

sumdelta..  delta(segment1) + delta(segment2) + delta(segment3) =E= 1 ; (LHS = 0, INFES = 1 ****)


----     49 VARIABLE x.L                   =        5.000  
            VARIABLE y.L                   =        6.000  

----     49 VARIABLE delta.L  select segment

segment2 1.000


----     49 VARIABLE lambda.L  interpolation

point2 0.333,    point3 0.667

The methods DCC and CC simulate our SOS2 model from [1] using binary variables. This means the model can handle step functions (we don't need to form a slope). The models are not too difficult to setup, as is illustrated by the GAMS models implementing them.

An early reference to this method is [4].

Method 6: MC - Multiple Choice formulation

This is a well-known method. We assume we can calculate slopes and intercepts for each segment. I.e. we have $$y = a_s + b_s x\>\> \text{ if $\bar{x}_s \le x \le \bar{x}_{s+1}$}$$   with $$\begin{align} & a_s = \frac{\bar{y}_{s+1}-\bar{y}_s}{\bar{x}_{s+1}-\bar{x}_s} \\ & b_s = \bar{y}_s -  a_s \bar{x}_s \end{align}$$

For the model we introduce binary variables $\delta_s$ to indicate which segment is selected, and so-called semi-continuous variables $v_s \in {0} \cup [\bar{x}_s,\bar{x}_{s+1}]$. The variable $v_s$ can be modeled as $$\bar{x}_s \delta_s \le v_s \le \bar{x}_{s+1} \delta_s$$

With this we can formulate the complete model:

MC - Multiple Choice Formulation

$$\begin{align} & \color{DarkRed}x = \sum_s \color{DarkRed}v_s \\ & \color{DarkRed}y = \sum_s \left( \color{DarkBlue}a_s \color{DarkRed} v_s + \color{DarkBlue}b_s \color{DarkRed} \delta_s \right) \\ &\color{DarkBlue}{\bar{x}}_s \color{DarkRed}\delta_s \le \color{DarkRed}v_s \le \color{DarkBlue}{\bar{x}}_{s+1} \color{DarkRed}\delta_s \\& \sum_s \color{DarkRed}\delta_s=1 \\ & \color{DarkRed}\delta_s \in \{0,1\}\\ & \color{DarkBlue}a_s = \frac{\color{DarkBlue}{\bar{y}}_{s+1}-\color{DarkBlue}{\bar{y}}_s}{\color{DarkBlue}{\bar{x}}_{s+1}-\color{DarkBlue}{\bar{x}}_s} \\ &\color{DarkBlue}b_s = \color{DarkBlue}{\bar{y}}_s - \color{DarkBlue} a_s \color{DarkBlue}{\bar{x}}_s \end{align}$$

Notes:

• The sandwich equation  $\bar{x}_s \delta_s \le v_s \le \bar{x}_{s+1} \delta_s$ must likely be implemented as two inequalities.
• This construct says: $v_s=0$ or $v_s \in [\bar{x}_s, \bar{x}_{s+1}]$. $v_s$ is sometimes called semi-continuous.
• The variables $\delta_s$ and $v_s$ are connected: $$\begin{align} &\delta_s = 0 \Rightarrow v_s = 0\\ & \delta_s = 1 \Rightarrow v_s = x\end{align}$$ I.e. they operate in parallel.
• I have seen cases where this model outperformed the SOS2 formulation.

The GAMS model is simple:

set   k 'breakpoints' /k1*k4/   s(k) 'segments' /k1*k3/ ; table data(k,*)        x   y   k1   1   6   k2   3   2   k3   6   8   k4  10   7 ; data(s(k),'dx') = data(k+1,'x')-data(k,'x'); data(s(k),'dy') = data(k+1,'y')-data(k,'y'); data(s(k),'slope') = data(k,'dy')/data(k,'dx'); data(s(k),'intercept') = data(k,'y')-data(k,'slope')*data(k,'x'); display data; variable v(s) 'equal to x or 0'; binary variable delta(s) 'select segment'; variable x,y; equations    xdef       'x'    ydef       'y'    semicont1(k)    semicont2(k)    sumdelta   'select segment' ; xdef.. x =e= sum(s, v(s)); ydef.. y =e= sum(s, data(s,'slope')*v(s)+data(s,'intercept')*delta(s)); semicont1(s(k)).. v(s) =l= data(k+1,'x')*delta(s); semicont2(s)..    v(s) =g= data(s,'y')*delta(s); sumdelta.. sum(s, delta(s)) =e= 1; x.fx = 5; model m /all/; option optcr=0; solve m maximizing y using mip; display data,x.l,y.l,delta.l,v.l;

The results look like:

----     43 PARAMETER data  

             x           y          dx          dy       slope   intercept

k1       1.000       6.000       2.000      -4.000      -2.000       8.000
k2       3.000       2.000       3.000       6.000       2.000      -4.000
k3       6.000       8.000       4.000      -1.000      -0.250       9.500
k4      10.000       7.000


----     43 VARIABLE x.L                   =        5.000  
            VARIABLE y.L                   =        6.000  

----     43 VARIABLE delta.L  select segment

k2 1.000


----     43 VARIABLE v.L  equal to x or 0

k2 5.000

Segment 2 is selected. We can see that directly from the binary variable $\delta_s$, but also from the semi-continuous variable $v_s$.

Method 7: INC - Incremental Formulation

In the incremental or delta method we add up all contributions of "earlier" segments, to find our $x$ and $y$.

Let $s'$ be the segment that contains our current $x$. We define a binary variable $\delta_{s}$ as $$\delta_{s} = \begin{cases} 1 & \text{for $s\lt s'$}\\ 0  & \text{for $s \ge s'$}\end{cases}$$ In addition we use a continuous variable $\lambda_s \in [0,1]$ indicating how much of each segment we "use up". The contribution for earlier segments is 1 and for the current segment we have a fractional value. So we have: $$\begin{cases}  \lambda_{s} = 1 & \text{for $s\lt s'$} \\\lambda_{s} \in [0,1] & \text{for $s= s'$}  \\ \lambda_{s} = 0 & \text{for $s\gt s'$}\end{cases}$$ With these definitions we can write: $$\begin{align} x = \bar{x}_1 + \sum_s \lambda_s (\bar{x}_{s+1} - \bar{x}_s) \\ y = \bar{y}_1 + \sum_s \lambda_s (\bar{y}_{s+1} - \bar{y}_s) \end{align}$$

Now we need to formulate a structure that enforce these rules on $\delta$ and $\lambda$. The following will do that: $$\lambda_{s+1} \le \delta_s \le \lambda_s$$ Typically you will need to implement this as two separate constraints.

The complete model looks like:

INC - Incremental Formulation

$$\begin{align} & \color{DarkRed}x = \color{DarkBlue}{\bar{x}}_1 + \sum_s \color{DarkRed}\lambda_s (\color{DarkBlue}{\bar{x}}_{s+1} - \color{DarkBlue}{\bar{x}}_s) \\ & \color{DarkRed}y = \color{DarkBlue}{\bar{y}}_1 + \sum_s \color{DarkRed}\lambda_s (\color{DarkBlue}{\bar{y}}_{s+1} - \color{DarkBlue}{\bar{y}}_s) \\ & \color{DarkRed} \lambda_{s+1} \le \color{DarkRed}\delta_s \le \color{DarkRed}\lambda_s \\ & \color{DarkRed} \delta_s \in \{0,1\} \\ & \color{DarkRed} \lambda_s \in [0,1] \end{align}$$

We can see that we can fix the last $\delta_s=0$. The model will behave correctly without this, but we may help the presolver a bit with this.


The GAMS model can look like:

set   k 'breakpoints' /k1*k4/   s(k) 'segments' /k1*k3/ ; table data(k,*)        x   y   k1   1   6   k2   3   2   k3   6   8   k4  10   7 ; positive variable lambda(s) 'contribution factor of segment'; lambda.up(s) = 1; binary variable delta(s) 'previously contributing segments'; variable x,y; equations    xdef       'x'    ydef       'y'    link1(s)   'links lambda delta'    link2(s)   'links lambda delta' ; * fix last delta(s) to 0. * this is not strictly needed delta.fx(s)$(ord(s)=card(s)) = 0; xdef.. x =e= data('k1','x')+ sum(s(k), lambda(s)*(data(k+1,'x')-data(k,'x'))); ydef.. y =e= data('k1','y')+ sum(s(k), lambda(s)*(data(k+1,'y')-data(k,'y'))); link1(s)..  lambda(s+1) =l= delta(s); link2(s)..  delta(s) =l= lambda(s); x.fx = 5; model m /all/; option optcr=0; solve m maximizing y using mip; display x.l,y.l,delta.l,lambda.l;

Note that in equation link1, we go one position too far when addressing $\lambda_{s+1}$ for the last $s$. GAMS will make that reference zero, and that is correct for this case. We see in the equation listing:

---- link1  =L=  links lambda delta

link1(k1)..  lambda(k2) - delta(k1) =L= 0 ; (LHS = 0)
     
link1(k2)..  lambda(k3) - delta(k2) =L= 0 ; (LHS = 0)
     
link1(k3)..  - delta(k3) =L= 0 ; (LHS = 0)

The solution is

----     36 VARIABLE x.L                   =        5.000  
            VARIABLE y.L                   =        6.000  

----     36 VARIABLE delta.L  previously contributing segments

k1 1.000


----     36 VARIABLE lambda.L  contribution factor of segment

k1 1.000,    k2 0.667

Notes:

• A slightly different formulation is the following. Use the definition: $$\delta_{s} = \begin{cases} 1 & \text{for $s\le s'$}\\ 0  & \text{for $s \gt s'$}\end{cases}$$ Now we use the sandwich equation: $$\delta_{s+1} \le \lambda_s \le \delta_s$$  With this formulation we can fix: $\delta_1=1$.
• This approach goes back to [5].

References

1. Piecewise linear functions and formulations for interpolation (part 1). This post discusses formulations SOS2, BIGM_BIN and BIGM_SOS1. The SOS2 formulation is the most important. Also illustrated how a Pyomo implementation is not implementing BIGM_BIN and BIGM_SOS1 correctly.  http://yetanothermathprogrammingconsultant.blogspot.com/2019/02/piecewise-linear-functions-and.html
2. Piecewise linear functions and formulations for interpolation (part 3). This post describes the LOG and DLOG formulations, which have a logarithmic number of binary variables. We see also some issues with the Pyomo implementation. http://yetanothermathprogrammingconsultant.blogspot.com/2019/03/piecewise-linear-functions-and.html
3. Keely L. Croxton, Bernard Gendron, Thomas L. Magnanti, A Comparison of Mixed-Integer Programming Models for Non-Convex Piecewise Linear Cost Minimization Problems, Management Science, Volume 49, Issue 9, 2003, pages 1121-1273
4. G. Dantzig, On the significance of solving linear programming problems with some integer variables, Econometrica 28, 30-44, 1960
5. H. Markowitz, A. Manne, On the solution of discrete programming problems, Econometrica 25, 84-110, 1957

Piecewise linear functions and formulations for interpolation (part 1)

Pyomo has a tool to use different formulations for piecewise linear functions. From the help:

>>> from pyomo.core import * >>> help(Piecewise) Help on class Piecewise in module pyomo.core.base.piecewise: class Piecewise(pyomo.core.base.block.Block)  |  Piecewise(*args, **kwds)  |  |      Adds piecewise constraints to a Pyomo model for functions of the  |      form, y = f(x).  |  |      Usage:  |              model.const = Piecewise(index_1,...,index_n,yvar,xvar,**Keywords)  |              model.const = Piecewise(yvar,xvar,**Keywords)  |  |      Keywords:  |  |  -pw_pts={},[],()  |            A dictionary of lists (keys are index set) or a single list  |            (for the non-indexed case or when an identical set of  |            breakpoints is used across all indices) defining the set of  |            domain breakpoints for the piecewise linear  |            function. **ALWAYS REQUIRED**  |  |  -pw_repn=''  |            Indicates the type of piecewise representation to use. This  |            can have a major impact on solver performance.  |            Choices: (Default 'SOS2')  |  |               ~ + 'SOS2'      - Standard representation using sos2 constraints  |               ~   'BIGM_BIN'  - BigM constraints with binary variables.  |                                 Theoretically tightest M values are automatically  |                                 determined.  |               ~   'BIGM_SOS1' - BigM constraints with sos1 variables.  |                                 Theoretically tightest M values are automatically  |                                 determined.  |               ~*+ 'DCC'       - Disaggregated convex combination model  |               ~*+ 'DLOG'      - Logarithmic disaggregated convex combination model  |               ~*+ 'CC'        - Convex combination model  |               ~*+ 'LOG'       - Logarithmic branching convex combination  |               ~*  'MC'        - Multiple choice model  |               ~*+ 'INC'       - Incremental (delta) method  |  |             + Supports step functions  |             * Source: "Mixed-Integer Models for Non-separable Piecewise Linear  |                        Optimization: Unifying framework and Extensions" (Vielma,  |                        Nemhauser 2008)  |             ~ Refer to the optional 'force_pw' keyword.  |  |  -pw_constr_type=''  |            Indicates the bound type of the piecewise function.  |            Choices: -- More  --

Some of these representations are not immediately obvious. Let's try to get a intuitive understanding of them. For a more formal discussion see [3].

I just consider the case I am using most often: a one-dimensional piecewise function $y = f(x)$. Note that in the context of an optimization model we need not necessarily to interpret this as: given an $x$, calculate an $y$,  but rather as part of a system of equations. I.e. $x$ and $y$ are determined at same time (if $y$ is known, we find $x$, or even: we find $x,y$ simultaneously).

The piecewise linear function is defined by its breakpoints:

Some observations:

1. We have $K$ breakpoints $(\bar{x}_k,\bar{y}_k)$ and $K-1$ segments,
2. we assume $\bar{x}_k$ is ordered: $\bar{x}_{k+1}\ge \bar{x}_k$,  
3. we assume we have lower and upper bounds on $x$ and $y$,
4. we assume the curve is connected (i.e. no forbidden regions for $x$, these have to be handled separately,
5. we allow step functions by adding a vertical segment (not all formulation allow this).

Finally, we note that interpolating between two points $(\bar{x}_1,\bar{y}_1)$, $(\bar{x}_2,\bar{y}_2)$ can be stated as choosing two weights $\lambda_1, \lambda_2$ with: $$\begin{align} & x = \lambda_1 \bar{x}_1 + \lambda_2 \bar{x}_2\\ & y = \lambda_1 \bar{y}_1 + \lambda_2 \bar{y}_2 \\ &  \lambda_1+\lambda_2 = 1 \\ & \lambda_1, \lambda_2 \ge 0 \end{align}$$

Method 1: SOS2 formulation

This one is the easiest to remember. Many more advanced MIP solvers have facilities for SOS2 variables [4]. The definition is: 

Let $x_1,\dots,x_n$ be members of a Special Ordered Set of Type 2 (SOS2), then only two variables $x_i$ can be nonzero and they have to be neighbors. All other variables are zero.

This looks like a strange animal, but in fact this construct is especially designed for cases like our piecewise linear interpolation scheme. Using SOS2 constraints, we can easily formulate a model: 

SOS2 Formulation
$$\begin{align} & \color{DarkRed}x = \sum_k \color{DarkBlue}{\bar{x}}_k \color{DarkRed}\lambda_k \\ & \color{DarkRed}y = \sum_k \color{DarkBlue}{\bar{y}}_k \color{DarkRed}\lambda_k\\& \sum_k \color{DarkRed} \lambda_k = 1 \\ & \color{DarkRed}\lambda_k \ge 0 \\ & \mathit{SOS2}(\color{DarkRed}\lambda_1,\dots,\color{DarkRed}\lambda_K) \end{align}$$

Note that $\color{DarkRed}x$ and $\color{DarkRed}y$ are decision variables, while $\color{DarkBlue}{\bar{x}}_k$ and $\color{DarkBlue}{\bar{y}}_k$ are constants. The color coding helps to distinguish variables from data.

In the solution, two adjacent $\lambda_s,\lambda_{s+1}$ will be between 0 and 1, and for these we have $\lambda_s+\lambda_{s+1}=1$. All the others will be zero. I.e. we will be interpolating between the two corresponding breakpoints $(\bar{x}_s, \bar{y}_s)$ and $(\bar{x}_{s+1}, \bar{y}_{s+1})$.  This SOS2 formulation does two things in one swoop: it selects the segment $s$ the variables $x$ and $y$ belong to, and it interpolates between the two selected breakpoints.

Notes:

• This method is used quite a lot in practice, not in the least because it makes life easy for the modeler. It is useful to have this formulation in your toolbox.
• This approach handles step functions without a problem.
• There are no issues with big-M constants: there are none.
• However, binary variables may offer performance benefits. 

The Pyomo model is as follows:

#
# expected solution X=5, Y=6
#

xdata = [1., 3., 6., 10.]
ydata  = [6.,2.,8.,7.]

from pyomo.core import *

model = ConcreteModel()

model.X = Var(bounds=(1,10))
model.Y = Var(bounds=(0,100))

model.con = Piecewise(model.Y,model.X,
                      pw_pts=xdata,
                      pw_constr_type='EQ',
                      f_rule=ydata,
                      pw_repn='SOS2')

# see what we get for Y when X=5
def con2_rule(model):
    return model.X==5

model.con2 = Constraint(rule=con2_rule)

model.obj = Objective(expr=model.Y, sense=maximize)

It is also instructive to see how this can be modeled in a more traditional modeling language like GAMS:

set   k 'breakpoints' /point1*point4/ ; table data(k,*)            x   y   point1   1   6   point2   3   2   point3   6   8   point4  10   7 ; sos2 variable lambda(k) 'sos2 variable'; variable x,y; equations    refrow    'reference row'    funrow    'function row'    convexity ; refrow.. x =e= sum(k, lambda(k)*data(k,'x')); funrow.. y =e= sum(k, lambda(k)*data(k,'y')); convexity.. sum(k, lambda(k)) =e= 1; x.fx = 5; model m /all/; option optcr=0; solve m maximizing y using mip; display x.l,y.l,lambda.l;

Results will look like:

----     29 VARIABLE x.L                   =        5.000  
            VARIABLE y.L                   =        6.000  

----     29 VARIABLE lambda.L  sos2 variable

point2 0.333,    point3 0.667

This is the correct solution when $x=5$.

Method 2: Big-M Binary Formulation

In this formulation we turn on or off each linear equality depending if we select a segment. Looking at our simple example, we can state the problem as: $$y = \begin{cases} -2x+8 & x \in [1,3]\\ 2x-4 & x \in [3,6]\\ -0.25x+9.5 & x \in [6,10]\end{cases}$$ We can model this "case statement" with binary variables as follows: define $$\delta_s = \begin{cases} 1 & \text{if segment $s$ is selected}\\ 0 & \text{otherwise}\end{cases}$$ When $\delta_s=0$ for a given segment $s$, we need to make the corresponding equation non-binding, and also the corresponding limits on $x$. This can be accomplished using inequalities and some big-M's. As an example, consider the second segment. We can write: $$\begin{align} &   2x-4-(1-\delta_2)M\le y \le  2x-4 +(1-\delta_2)M \\ & 3 - (1-\delta_2)M \le x \le 6+ (1-\delta_2)M\end{align}$$ Finding good values for each $M$ is a bit tedious. Basically we have to choose the smallest $M$ such that each breakpoint is reachable (i.e. not cut off). For the constraint $$y \le  2x-4 +(1-\delta_2)M$$ we can write: $$\begin{align} M & = \max_k \{\bar{y}_k - 2\bar{x}_k+4\}\\ & = \max \{ 6-2\cdot 1+4,2-2\cdot 3 + 4, 8-2\cdot 6 +4 , 7-2 \cdot 10+4 \} \\ & = \max \{8,0,0,-9\}\\&=8\end{align}$$ The zero values correspond to breakpoints on the current segment.

A high level description can look like: $$\begin{align} &\delta_s = 1 \Rightarrow \begin{cases} y = a_s x + b_s \\  \bar{x}_s \le x \le \bar{x}_{s+1} \end{cases}\\ &\sum_s \delta_s = 1\\ &\delta_s \in \{0,1\}\end{align}$$ where $a_s, b_s$ are the slope and intercept for the linear function in segment $s$. These coefficients can be calculated easily from the breakpoints $(\bar{x}_s,\bar{y}_s)$ and $(\bar{x}_{s+1},\bar{y}_{s+1})$. This implication can be implemented as follows: 

Big-M Binary Formulation

$$\begin{align} & \color{DarkBlue}a_s \color{DarkRed}x + \color{DarkBlue}b_s - (1-\color{DarkRed}\delta_s) \color{DarkBlue} M \le \color{DarkRed}y \le \color{DarkBlue}a_s \color{DarkRed}x + \color{DarkBlue}b_s + (1-\color{DarkRed}\delta_s) \color{DarkBlue} M \\ & \color{DarkBlue}{\bar{x}}_s - (1-\color{DarkRed}\delta_s)\color{DarkBlue}M \le \color{DarkRed}x \le \color{DarkBlue}{\bar{x}}_{s+1} + (1-\color{DarkRed}\delta_s)\color{DarkBlue}M \\ & \sum_s \color{DarkRed} \delta_s=1\\ & \color{DarkRed}\delta_s \in \{0,1\}\\ & \color{DarkBlue}a_s = \frac{\color{DarkBlue}{\bar{y}}_{s+1}-\color{DarkBlue}{\bar{y}}_s}{\color{DarkBlue}{\bar{x}}_{s+1}-\color{DarkBlue}{\bar{x}}_s} \\ &\color{DarkBlue}b_s = \color{DarkBlue}{\bar{y}}_s - \color{DarkBlue} a_s \color{DarkBlue}{\bar{x}}_s \end{align}$$

This looks a bit more intimidating, but it is conceptually simple. Unfortunately we see some problems with the Pyomo implementation of this formulation.

Notes:

• We need 4 inequalities per segment, and one binary variable.
• This method assumes we can calculate a slope $a_s$. For step functions we have an infinite slope when there is a vertical segment. I.e. we cannot do step functions using this method.
• It is imperative to choose reasonable values for the big-M's.
• In the model above I just used one $M$, but actually each of them can have a different value.
• This is not a formulation I typically use: it is somewhat cumbersome to setup. 
• Pyomo implements this method incorrectly, and we will get wrong results.

Method 3: Big-M SOS1 formulation

This is just a minor variant on the previous method. Instead of using binary variables, we use continuous variables that are members of a Special Ordered Set of Type 1. SOS1 sets are defined by:

Let $x_1,\dots,x_n$ be members of a Special Ordered Set of Type 1 (SOS1), then only one variable $x_i$ can be nonzero. All other variables are zero.

Basically we replace the binary variables: $$\begin{align} &\sum_s \delta_s = 1\\ & \delta_s \in \{0,1\}\end{align}$$ by $$\begin{align} &\sum_s \lambda_s = 1 &\\ & \lambda_s \ge 0 \\ & \mathit{SOS1}(\lambda_1,\dots,\lambda_{K-1})\end{align}$$ This does not change the model a lot:  

Big-M SOS1 Formulation

$$\begin{align} & \color{DarkBlue}a_s \color{DarkRed}x + \color{DarkBlue}b_s - (1-\color{DarkRed}\lambda_s) \color{DarkBlue} M \le \color{DarkRed}y \le \color{DarkBlue}a_s \color{DarkRed}x + \color{DarkBlue}b_s + (1-\color{DarkRed}\lambda_s) \color{DarkBlue} M \\ & \color{DarkBlue}{\bar{x}}_s - (1-\color{DarkRed}\lambda_s)\color{DarkBlue}M \le \color{DarkRed}x \le \color{DarkBlue}{\bar{x}}_{s+1} + (1-\color{DarkRed}\delta_s)\color{DarkBlue}M \\ & \sum_s \color{DarkRed} \lambda_s=1\\ & \color{DarkRed}\lambda_s\ge 0\\& \mathit{SOS1}(\color{DarkRed}\lambda_1,\dots\,\color{DarkRed}\lambda_{K-1})\\ & \color{DarkBlue}a_s = \frac{\color{DarkBlue}{\bar{y}}_{s+1}-\color{DarkBlue}{\bar{y}}_s}{\color{DarkBlue}{\bar{x}}_{s+1}-\color{DarkBlue}{\bar{x}}_s} \\ &\color{DarkBlue}b_s = \color{DarkBlue}{\bar{y}}_s - \color{DarkBlue} a_s \color{DarkBlue}{\bar{x}}_s \end{align}$$

Notes:

• This method shows the same Pyomo bug as the big-M binary variable model: it is incorrectly implemented. See the analysis below for what is wrong.
• If a solver support SOS1 variables, it also will typically support SOS2 variables. The SOS2 approach may be easier to use.
• This method is supposedly marked for removal in future versions of Pyomo. 
• We can use the SOS1 variables to prevent the need for big-M's. We can write something like $$\begin{align} & a_s x+b_s -\lambda_s \le y \le a_s x+b_s +\lambda_s\\ & \bar{x}_s - \lambda_s \le x \le \bar{x}_{s+1} + \lambda_s \\ & \sum_s \delta_s = 1\\& \lambda_s \ge 0\\& \delta_s \in \{0,1\} \\ & \mathit{SOS1}(\lambda_s,\delta_s) \end{align}$$ If you look at this carefully, you can see this implements the implications: $$\begin{align} & \delta_s = 1 \Rightarrow \begin{cases} y=a_s x + b_s \\ \bar{x}_s \le x \le \bar{x}_{s+1}\end{cases}\\ &\sum_s \delta_s = 1\\ & \delta_s \in \{0,1\} \end{align}$$ This is similar to how a MIP solver may rewrite implications. This particular way of dealing with an implication is a useful concept to know.
• I never have used the BIGM_SOS1 formulation in my models. It just does not look very appetizing.

A nasty Pyomo bug

To test the bigM-bin method on our simple data set I formulated the model: $$\begin{align}\max \> & y \\ & \text{piecewise on our data} \\ & x=5 \\ & x\in [1,10]\end{align}$$ We expect to see as solution: $x=5, y=6$ and an objective of 6. Instead we see:

D:\Python\Python37\Scripts>pyomo solve --solver=glpk bigm.py [    0.00] Setting up Pyomo environment [    0.00] Applying Pyomo preprocessing actions WARNING: DEPRECATED: The 'BIGM_BIN' and 'BIGM_SOS1' piecewise representations     will be removed in a future version of Pyomo. They produce incorrect     results in certain cases [    0.60] Creating model [    0.60] Applying solver [    0.83] Processing results     Number of solutions: 1     Solution Information       Gap: 0.0       Status: optimal       Function Value: 8.25     Solver results file: results.json [    0.86] Applying Pyomo postprocessing actions [    0.86] Pyomo Finished D:\Python\Python37\Scripts>

Well, this is disappointing. We get a wrong solution!

We also get a warning that the Pyomo people are aware of problems with this method.

This is more wrong than what can be explained away with some numerical issues related to large big-M values. Large big-M values can lead to wrong results. I don't think that is the case here, but let's investigate.

The details of the Pyomo error

This section may become a bit dreary: we will dive into the details of this bug. Not for the faint of heart.

The model looks like the Pyomo model shown in the SOS2 section. The only difference is that we used a different method: BIGM_BIN instead of SOS2.

The solution file is:

{
    "Problem": [
        {
            "Lower bound": 8.25,
            "Name": "unknown",
            "Number of constraints": 9,
            "Number of nonzeros": 21,
            "Number of objectives": 1,
            "Number of variables": 6,
            "Sense": "maximize",
            "Upper bound": 8.25
        }
    ],
    "Solution": [
        {
            "number of solutions": 1,
            "number of solutions displayed": 1
        },
        {
            "Constraint": "No values",
            "Gap": 0.0,
            "Message": null,
            "Objective": {
                "obj": {
                    "Value": 8.25
                }
            },
            "Problem": {},
            "Status": "optimal",
            "Variable": {
                "X": {
                    "Value": 5.0
                },
                "Y": {
                    "Value": 8.25
                },
                "con.bin_y[3]": {
                    "Value": 1.0
                }
            }
        }
    ],
    "Solver": [
        {
            "Error rc": 0,
            "Statistics": {
                "Branch and bound": {
                    "Number of bounded subproblems": "1",
                    "Number of created subproblems": "1"
                }
            },
            "Status": "ok",
            "Termination condition": "optimal",
            "Time": 0.12992453575134277
        }
    ]
}

We see two problems:

• Clearly the Y value is wrong: it should be 6 instead of 8.25. 
• We see X=5 and con.bin_y[3]=1. The last value indicates we are in the third segment. But the third segment is $[6,10]$ and does not include $x=5$.

To debug this, we can generate and inspect the LP file. This is a bit hard to read:

\* Source Pyomo model name=unknown *\

max 
obj:
+1 Y

s.t.

c_e_con2_:
+1 X
= 5

c_u_con_BIGM_constraint1(1)_:
+2 X
+1 Y
+19 con_bin_y(1)
<= 27

c_u_con_BIGM_constraint1(2)_:
-2 X
+1 Y
+8 con_bin_y(2)
<= 4

c_u_con_BIGM_constraint1(3)_:
+0.25 X
+1 Y
<= 9.5

c_e_con_BIGM_constraint2_:
+1 con_bin_y(1)
+1 con_bin_y(2)
+1 con_bin_y(3)
= 1

c_l_con_BIGM_constraint3(1)_:
+2 X
+1 Y
>= 8

c_u_con_BIGM_constraint3(2)_:
+2 X
-1 Y
+9 con_bin_y(2)
<= 13

c_u_con_BIGM_constraint3(3)_:
-0.25 X
-1 Y
+6.75 con_bin_y(3)
<= -2.75

c_e_ONE_VAR_CONSTANT: 
ONE_VAR_CONSTANT = 1.0

bounds
   1 <= X <= 10
    -inf <= Y <= +inf
   0 <= con_bin_y(1) <= 1
   0 <= con_bin_y(2) <= 1
   0 <= con_bin_y(3) <= 1
binary
  con_bin_y(1)
  con_bin_y(2)
  con_bin_y(3)
end

Let's compare carefully what we expect and what we see:

Constraint	Expected	Generated by Pyomo	Notes
Objective	$$\max\>y$$	max obj: +1 Y
Segment 1	$$y \le -2x+8 +(1-\delta_1)M$$	+2X+1Y+19con_bin_y(1)<=27	$M=19$
	$$y \ge -2x+8 -(1-\delta_1)M$$	+2X+1Y>=8	$M=0$
	$$x \ge 1 -(1-\delta_1)M$$		This one is not needed.
	$$x \le 3 + (1-\delta_1)M$$		$M=7$
Segment 2	$$y \le 2x-4 +(1-\delta_2)M$$	-2X+1Y+8con_bin_y(2)<=4	$M=8$
	$$y \ge 2x-4 -(1-\delta_2)M$$	+2X-1Y+9con_bin_y(2)<=13	$M=9$
	$$x \ge 3 -(1-\delta_2)M$$		$M=2$
	$$x \le 6 + (1-\delta_2)M$$		$M=4$
Segment 3	$$y \le -0.25x+9.5 +(1-\delta_3)M$$	+0.25X+1Y<=9.5	$M=0$
	$$y \ge -0.25x+9.5 -(1-\delta_3)M$$	0.25X-1Y+6.75con_bin_y(3)<=-2.75	$M=6.75$
	$$x \ge 6 -(1-\delta_3)M$$		$M=5$
	$$x \le 10 + (1-\delta_3)M$$		This one is not needed.
Sum	$$\sum_s \delta_s = 1$$	+1con_bin_y(1) +1con_bin_y(2) +1con_bin_y(3)=1
Fix	$$x=5$$	+1X=5

We see we have missing inequalities relating $x$ to the segments. The cells marked in yellow indicate the missing constraints. No wonder we got the wrong results. The omission to include the condition $$\delta_s = 1 \Rightarrow \bar{x}_s \le x \le \bar{x}_{s+1}$$ explains our erroneous solution: we are looking at the wrong segment.

We cannot blame the big-M's for this problem. The Pyomo piecewise function is just forgetting to include a number of inequalities.

This was an interesting exercise: debugging without looking at the code.

References

1. Piecewise linear functions and formulations for interpolation (part 2). This post discusses formulations DCC, CC, MC and INCR. http://yetanothermathprogrammingconsultant.blogspot.com/2019/02/piecewise-linear-functions-and_22.html
2. Piecewise linear functions and formulations for interpolation (part 3). This post describes the LOG and DLOG formulations, which have a logarithmic number of binary variables. We see also some issues with the Pyomo implementation. http://yetanothermathprogrammingconsultant.blogspot.com/2019/03/piecewise-linear-functions-and.html
3. Juan Pablo Vielma, Shabbir Ahmed and George Nemhauser, Mixed-Integer Models for Nonseparable Piecewise Linear Optimization: Unifying Framework and Extensions, Operations Research 58(2):303-315, 2010
4. E.M.L. Beale, J.J.H. Forrest, Global Optimization Using Special Ordered Sets, Math Prog 10, 52-69, 1976

Appendix: script to calculate big-M's in BIGM_INT and BIGM_SOS1 models

It is easy to make mistakes in the calculate of the big M values for the BIGM_INT and BIGM_SOS1 models. I used the following model to calculate them:

set   k 'breakpoints' /k1*k4/   s(k) 'segments' /k1*k3/ ; table data(k,*)        x   y   k1   1   6   k2   3   2   k3   6   8   k4  10   7 ; data(s(k),'dx') = data(k+1,'x')-data(k,'x'); data(s(k),'dy') = data(k+1,'y')-data(k,'y'); data(s(k),'slope') = data(k,'dy')/data(k,'dx'); data(s(k),'intercept') = data(k,'y')-data(k,'slope')*data(k,'x'); display data; parameter bound(*,*) 'implied bounds on x and y'; bound('x','lo') = smin(k,data(k,'x')); bound('x','up') = smax(k,data(k,'x')); bound('y','lo') = smin(k,data(k,'y')); bound('y','up') = smax(k,data(k,'y')); display bound; parameter M(*,*,s) 'tight big-M values'; M('y','>=',s) = -smin(k, data(k,'y') - data(s,'slope')*data(k,'x') - data(s,'intercept')); M('y','<=',s) = smax(k, data(k,'y') - data(s,'slope')*data(k,'x') - data(s,'intercept')); M('x','>=',s(k)) = max(0,data(k,'x')-bound('x','lo')); M('x','<=',s(k)) = max(0,bound('x','up')-data(k+1,'x')); display M; parameter num(s) 'ordinal number of each s'; num(s) = ord(s); variable x,y; binary variable delta(s) 'selection of segment'; equations    ylo(s)    yup(s)    xlo(s)    xup(k)    select ; ylo(s)..    y =g= data(s,'slope')*x + data(s,'intercept') - M('y','>=',s)*(1-delta(s)); yup(s)..    y =l= data(s,'slope')*x + data(s,'intercept') + M('y','<=',s)*(1-delta(s)); xlo(s)..    x =g= data(s,'x') - M('x','>=',s)*(1-delta(s)); xup(s(k)).. x =l= data(k+1,'x') + M('x','<=',s)*(1-delta(s)); select ..   sum(s, delta(s)) =e= 1; x.fx = 5; model bigmbin /all/; option optcr=0; solve bigmbin using mip maximizing y;

The output is:

----     18 PARAMETER data  

             x           y          dx          dy       slope   intercept

k1       1.000       6.000       2.000      -4.000      -2.000       8.000
k2       3.000       2.000       3.000       6.000       2.000      -4.000
k3       6.000       8.000       4.000      -1.000      -0.250       9.500
k4      10.000       7.000


----     25 PARAMETER bound  implied bounds on x and y

           lo          up

x       1.000      10.000
y       2.000       8.000


----     32 PARAMETER M  tight big-M values

              k1          k2          k3

x.>=                   2.000       5.000
x.<=       7.000       4.000
y.>=                   9.000       6.750
y.<=      19.000       8.000

The generated model looks like:

---- ylo  =G=  

ylo(k1)..  2*x + y =G= 8 ; (LHS = 10)
     
ylo(k2)..  - 2*x + y - 9*delta(k2) =G= -13 ; (LHS = -10)
     
ylo(k3)..  0.25*x + y - 6.75*delta(k3) =G= 2.75 ; (LHS = 1.25, INFES = 1.5 ****)
     

---- yup  =L=  

yup(k1)..  2*x + y + 19*delta(k1) =L= 27 ; (LHS = 10)
     
yup(k2)..  - 2*x + y + 8*delta(k2) =L= 4 ; (LHS = -10)
     
yup(k3)..  0.25*x + y =L= 9.5 ; (LHS = 1.25)
     

---- xlo  =G=  

xlo(k1)..  x =G= 1 ; (LHS = 5)
     
xlo(k2)..  x - 2*delta(k2) =G= 1 ; (LHS = 5)
     
xlo(k3)..  x - 5*delta(k3) =G= 1 ; (LHS = 5)
     

---- xup  =L=  

xup(k1)..  x + 7*delta(k1) =L= 10 ; (LHS = 5)
     
xup(k2)..  x + 4*delta(k2) =L= 10 ; (LHS = 5)
     
xup(k3)..  x =L= 10 ; (LHS = 5)
     

---- select  =E=  

select..  delta(k1) + delta(k2) + delta(k3) =E= 1 ; (LHS = 0, INFES = 1 ****)