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