## Friday, February 22, 2019

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

1. DCC: Disaggregated Convex Combination formulation is a simulation of the SOS2 model by binary variables.
2. CC:  Convex Combination formulation is a similar SOS2 like approach using binary variables only.
3. MC: Multiple Choice model uses a semi-continuous approach.
4. 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(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(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