MINLP instead of indicator constraints?

In this post, I want to discuss indicator constraints and how to replace them with simple multiplications. As we shall see, this is a somewhat harebrained but still interesting idea. 

Indicator constraints are implications of the form: \[\delta=0 \implies \text{linear constraint}\] or \[\delta=1 \implies \text{linear constraint}\] where \(\delta \in \{0,1\}\) is a binary decision variable.

There are two aspects of indicator constraints: 

• Indicator constraints help with MIP models where otherwise we would use big-M constraints. This will help address the numerical issues that result from using big-M constraints. These include small (and sometimes not-so-small) solution values where you expect zero, and poor solver performance. With big-M constraints, we need to pay much attention to the size of the big-M constants. This is the algorithmic aspect.
• Indicator constraints provide a convenient modeling construct. They form a useful abstraction that makes MIP modeling easier and more straightforward. This is the modeling aspect.

Another example where a solver feature has such a dual benefit (algorithmic and modeling) are SOS2 sets for piecewise linear functions. 

Convex hull models

Convex hull as an optimization problem

In the previous posts, the construction of a convex hull played a significant role: it was an easy way to reduce the size of the data sets, often by a large amount. Somehow, it escaped me that we can try to formulate this as what turns out to be a conceptually rather simple optimization problem. This has probably little or no practical value. But it remains an interesting exercise.

MINLP vs MIQCP

In [1] a small non-convex MINLP is stated. The elements are:

• We have binary variables \(x_i \in \{0,1\}\) for \(i=1,\dots,52\)
• The data is comprised of three arrays: \(a_i\), \(b_i\) and \(c_i\)
• We want to minimize the following objective: \[\begin{align} - &\left(\frac{1166}{2000} \sum_i a_i \cdot x_i\right) \times \\ & \left(\frac{1}{2100} \sum_i b_i \cdot x_i + 0.05\right) \times \\ & \left(\frac{1}{1500} \sum_i c_i \cdot x_i + 1.5\right)\end{align}\] 
• Of the first four elements \(x_i\), only one has the value 1, and the other three are 0. Same thing for the second, third, etc. groups of four. 

Data

----     66 SET i  items

i1 ,    i2 ,    i3 ,    i4 ,    i5 ,    i6 ,    i7 ,    i8 ,    i9 ,    i10,    i11,    i12,    i13,    i14,    i15
i16,    i17,    i18,    i19,    i20,    i21,    i22,    i23,    i24,    i25,    i26,    i27,    i28,    i29,    i30
i31,    i32,    i33,    i34,    i35,    i36,    i37,    i38,    i39,    i40,    i41,    i42,    i43,    i44,    i45
i46,    i47,    i48,    i49,    i50,    i51,    i52


----     66 SET g  groups

group1 ,    group2 ,    group3 ,    group4 ,    group5 ,    group6 ,    group7 ,    group8 ,    group9 ,    group10
group11,    group12,    group13


----     66 SET group  grouping of items

                 i1          i2          i3          i4          i5          i6          i7          i8          i9

group1          YES         YES         YES         YES
group2                                                          YES         YES         YES         YES
group3                                                                                                          YES

      +         i10         i11         i12         i13         i14         i15         i16         i17         i18

group3          YES         YES         YES
group4                                              YES         YES         YES         YES
group5                                                                                              YES         YES

      +         i19         i20         i21         i22         i23         i24         i25         i26         i27

group5          YES         YES
group6                                  YES         YES         YES         YES
group7                                                                                  YES         YES         YES

      +         i28         i29         i30         i31         i32         i33         i34         i35         i36

group7          YES
group8                      YES         YES         YES         YES
group9                                                                      YES         YES         YES         YES

      +         i37         i38         i39         i40         i41         i42         i43         i44         i45

group10         YES         YES         YES         YES
group11                                                         YES         YES         YES         YES
group12                                                                                                         YES

      +         i46         i47         i48         i49         i50         i51         i52

group12         YES         YES         YES
group13                                             YES         YES         YES         YES


----     66 PARAMETER data  

              a           b           c

i1      251.000     179.000     179.000
i2      179.000     251.000     179.000
i3      215.000     215.000     118.000
i4      251.000                 179.000
i5       63.000      45.000      45.000
i6       45.000      63.000      45.000
i7       54.000      54.000      30.000
i8       63.000                  45.000
i9       47.000      34.000      34.000
i10      34.000      47.000      34.000
i11      40.000      40.000      22.000
i12      47.000                  34.000
i13     141.000     101.000     101.000
i14     101.000     141.000     101.000
i15     121.000     121.000      67.000
i16     141.000                 101.000
i17      47.000      34.000      34.000
i18      34.000      47.000      34.000
i19      40.000      40.000      22.000
i20      47.000                  34.000
i21      94.000      67.000      67.000
i22      67.000      94.000      67.000
i23      81.000      81.000      44.000
i24      94.000                  67.000
i25      47.000      34.000      34.000
i26      34.000      47.000      34.000
i27      40.000      40.000      22.000
i28      47.000                  34.000
i29     157.000     108.000     108.000
i30     108.000     157.000     108.000
i31     133.000     133.000      71.000
i32     157.000                 108.000
i33     126.000      85.000      85.000
i34      85.000     126.000      85.000
i35     106.000     106.000      56.000
i36     126.000                  85.000
i37     126.000      85.000      85.000
i38      85.000     126.000      85.000
i39     106.000     106.000      56.000
i40     126.000                  85.000
i41     110.000      74.000      74.000
i42      74.000     110.000      74.000
i43      92.000      92.000      49.000
i44     110.000                  74.000
i45     110.000      74.000      74.000
i46      74.000     110.000      74.000
i47      92.000      92.000      49.000
i48     110.000                  74.000
i49      63.000      40.000      40.000
i50      40.000      63.000      40.000
i51      52.000      52.000      27.000
i52      63.000                  40.000

The set group was populated using the GAMS statement:

set group(g,i) 'grouping of items'; group(g,i) = 1+floor((ord(i)-0.5)/4)=ord(g);

The term 0.5 is to make sure we don't take the floor of integer values (the floor function can then deliver different results depending on the very last bit - possibly dangerous when we deal with floating-point computations).  

Direct MINLP formulation

A direct MINLP model formulation can look like:

MINLP model A

\[\begin{align}\min\> - &\left(\frac{1166}{2000} \sum_i \color{darkblue}a_i \cdot \color{darkred}x_i\right) \times \\ & \left(\frac{1}{2100} \sum_i \color{darkblue}b_i \cdot \color{darkred}x_i + 0.05\right) \times \\ & \left(\frac{1}{1500} \sum_i \color{darkblue}c_i \cdot \color{darkred}x_i + 1.5\right)\\ & \sum_{i|\color{darkblue}{\mathit{group}}(g,i)} \color{darkred}x_i = 1 &&\forall g\\ & \color{darkred}x_i \in \{0,1\}\end{align}\]

The solution can look like: 

 

----     89 VARIABLE x.L  

i1  1.000,    i5  1.000,    i9  1.000,    i14 1.000,    i17 1.000,    i21 1.000,    i25 1.000,    i29 1.000
i34 1.000,    i38 1.000,    i41 1.000,    i46 1.000,    i49 1.000


----     89 VARIABLE z.L                   =     -889.346  obj

Of course, we see 13 items selected (one for each group of 4).

This is a small model, so we can try a few different global and local solvers.

 

Solver	Type	Obj	Time	Nodes
Couenne	global	-889.3463	18.25	9972
Baron	global	-889.3463	7.18	389
Antigone	global	Stuck in preprocessing
Bonmin	local	-889.3463	89	2280

Interesting:

• Antigone has problems during preprocessing
• The local solver Bonmin finds the global optimum but is somewhat slow 

Alternative formulation

There is an alternative formulation. We can substitute out the linear parts of the objective function. This will make the model larger (extra variables and equations), but less non-linear (e.g. in terms of non-linear variables). 

MINLP model B

\[\begin{align}\max & \prod_k\color{darkred}v_k \\ & \color{darkred}v_1 = \frac{1166}{2000} \sum_i \color{darkblue}a_i \cdot \color{darkred}x_i \\ & \color{darkred}v_2 = \frac{1}{2100} \sum_i \color{darkblue}b_i \cdot \color{darkred}x_i + 0.05 \\ & \color{darkred}v_3 = \frac{1}{1500} \sum_i \color{darkblue}c_i \cdot \color{darkred}x_i + 1.5\\ & \sum_{i|\color{darkblue}{\mathit{group}}(g,i)} \color{darkred}x_i = 1 &&\forall g\\ & \color{darkred}x_i \in \{0,1\}\end{align}\]

This formulation has more linear constraints, but in return, we get a much simpler non-linear objective. Note that we also got rid of the minus sign in the objective. This was probably put in there to convert the problem to a minimization problem.

The solution is obviously the same:

----    122 VARIABLE x.L  

i1  1.000,    i6  1.000,    i9  1.000,    i13 1.000,    i17 1.000,    i22 1.000,    i25 1.000,    i29 1.000
i33 1.000,    i38 1.000,    i42 1.000,    i46 1.000,    i49 1.000


----    122 VARIABLE v.L  obj factors

obj1 713.592,    obj2   0.582,    obj3   2.140


----    122 VARIABLE z.L                   =      889.346  obj

Some of the solvers perform better on this model:

 

Solver	Type	Obj	Time	Nodes
Couenne	global	889.3463	5.1	4398
Baron	global	889.3463	9.5	309
Antigone	global	889.3463	2	5
Bonmin	local	889.3463	105	3336

Couenne and especially Antigone do much better on this formulation.

Non-convex quadratic formulation

Gurobi can solve non-convex quadratic models, so it is interesting to reformulate the problem into a form that Gurobi can digest.

 

MIQCP model

\[\begin{align}\max\> & \color{darkred}v_1 \cdot \color{darkred}v_{23} \\& \color{darkred}v_{23} = \color{darkred}v_{2} \cdot \color{darkred}v_{3}\\ & \color{darkred}v_1 = \frac{1166}{2000} \sum_i \color{darkblue}a_i \cdot \color{darkred}x_i \\ & \color{darkred}v_2 = \frac{1}{2100} \sum_i \color{darkblue}b_i \cdot \color{darkred}x_i + 0.05 \\ & \color{darkred}v_3 = \frac{1}{1500} \sum_i \color{darkblue}c_i \cdot \color{darkred}x_i + 1.5\\ & \sum_{i|\color{darkblue}{\mathit{group}}(g,i)} \color{darkred}x_i = 1 &&\forall g\\ & \color{darkred}x_i \in \{0,1\}\end{align}\]

Note that this requires the Gurobi option nonConvex=2.

 

Solver	Type	Obj	Time	Nodes
Gurobi	global	889.3463	0.26	6985

Gurobi solves this model very quickly, although the node count is not that low.

Conclusion

Even this small MINLP model has some interesting wrinkles and reformulations.

References

1. R Binary Integer Optimization with Groups, https://stackoverflow.com/questions/63584707/r-binary-integer-optimization-with-groups

Non-convex quadratic models II: rectangular covering

Rectangular Covering

This problem is from [1].

Problem statement: given \(N\) data points, arrange \(k\) boxes such that every point is inside a box and the size of the boxes is minimized. We assume the boxes have sides parallel to the axes. I.e. the box is not angled.

A 2d data set with \(N=50\) points can look like:

N=50 data points 

The best configuration of \(K=5\) boxes for this data set is:

Optimal solution for K=5 rectangles

The problem is like the opposite of the "largest empty box" [2], where we try to find the largest box without any data points inside.

2D Model

Let's define some symbols:

• We will denote the data points by \(i\) and the rectangles by \(k\).
• The coordinates of the data points are stored in \(p_{i,c} \in [L,U]\) with \(c \in \{x,y\}\).
• The location and size of the rectangles are \(r_{k,s}\in [L,U]\) . Here \(s\in \{x,y,w,h\}\).
• The assignment of (data) points to rectangles is modeled with a binary variable \(x_{i,k} \in \{0,1\}\).  

A high-level model for 2 dimensions can look like:

High-level 2D MIQP model

\[\begin{align}\min & \sum_k \color{darkred} r_{k,w} \cdot \color{darkred} r_{k,h} \\ & \sum_k \color{darkred}x_{i,k} = 1 && \forall i  \\ & \color{darkred}x_{i,k}=1 \implies \begin{cases} \color{darkblue}p_{i,x} \ge \color{darkred}r_{k,x} \\  \color{darkblue}p_{i,y} \ge \color{darkred}r_{k,y} \\ \color{darkblue}p_{i,x} \le \color{darkred}r_{k,x}+\color{darkred}r_{k,w} \\ \color{darkblue}p_{i,y} \le \color{darkred}r_{k,y}+\color{darkred}r_{k,h}  \end{cases} && \forall i,k\\ & \color{darkred} x_{i,k}\in \{0,1\} \\ & \color{darkred}r_{k,h} \in [\color{darkblue}L,\color{darkblue}U] \\ &\color{darkred}r_{k,x}+\color{darkred}r_{k,w} \le \color{darkblue}U \\ &\color{darkred}r_{k,y}+\color{darkred}r_{k,h} \le \color{darkblue}U   \end{align}\]

The implication says "if \(x_{i,k}=1\) then point \(i\) should be in box \(k\)". This can be implemented using indicator constraints or with a big-M construct: \[\begin{align} & p_{i,x} \ge r_{k,x} -M(1-x_{i,k})\\ & p_{i,y} \ge r_{k,y}-M(1-x_{i,k}) \\ & p_{i,x} \le r_{k,x}+r_{k,w} +M(1-x_{i,k})\\ & p_{i,y} \le r_{k,y}+r_{k,h}+M(1-x_{i,k}) \end{align}\] A value for \(M\) can be \(M=U-L\). We can of course scale the data \(p_{i,c}\) to \(p_{i,c} \in [0,1]\)  beforehand, in which case we can set \(M=1\). This is what I did in my experiments.

I also add the constraint: \[r_{k,x} \ge r_{k-1,x}\] to make the solution more predictable. The solution will be ordered by the \(x\)-coordinate of the box. It also reduces symmetry in the model.

The objective has non-convex quadratic terms. So we need a non-convex quadratic solver.

The detailed data and the solution look like:

----     56 PARAMETER p  data points

              x           y

i1      0.80601     0.17283
i2      0.52952     0.14887
i3      0.64826     0.69240
i4      0.35222     0.02016
i5      0.43130     0.55356
i6      0.64087     0.77456
i7      0.23545     0.78074
i8      0.26800     0.08158
i9      0.97328     0.11440
i10     0.87431     0.66711
i11     0.75621     0.96771
i12     0.19863     0.24001
i13     0.21990     0.26101
i14     0.98863     0.17156
i15     0.06625     0.93008
i16     0.80603     0.83235
i17     0.10543     0.02902
i18     0.22927     0.09355
i19     0.12983     0.90322
i20     0.43682     0.72756
i21     0.24754     0.57512
i22     0.36012     0.51628
i23     0.71015     0.74601
i24     0.70367     0.74621
i25     0.18517     0.93583
i26     0.81681     0.67338
i27     0.46350     0.57754
i28     0.08932     0.65688
i29     0.97335     0.69059
i30     0.89449     0.07780
i31     0.14124     0.39550
i32     0.99382     0.49870
i33     0.05406     0.56047
i34     0.32927     0.79803
i35     0.22975     0.85134
i36     0.78532     0.27930
i37     0.18014     0.28874
i38     0.38628     0.82680
i39     0.48970     0.57654
i40     0.46709     0.14913
i41     0.19949     0.83893
i42     0.25438     0.57046
i43     0.55116     0.95672
i44     0.02634     0.93260
i45     0.49950     0.08040
i46     0.44529     0.43742
i47     0.55100     0.41666
i48     0.29626     0.69398
i49     0.93426     0.58556
i50     0.28083     0.92941


----     56 VARIABLE x.L  assignment variables

             k1          k2          k3          k4          k5

i1                                                            1
i2                                    1
i3            1
i4                                    1
i5                                                1
i6            1
i7            1
i8                                    1
i9                                                            1
i10                                                           1
i11           1
i12                       1
i13                       1
i14                                                           1
i15           1
i16           1
i17                                   1
i18                                   1
i19           1
i20           1
i21                       1
i22                                               1
i23           1
i24           1
i25           1
i26                                                           1
i27                                               1
i28                       1
i29                                                           1
i30                                                           1
i31                       1
i32                                                           1
i33                       1
i34           1
i35           1
i36                                                           1
i37                       1
i38           1
i39                                               1
i40                                   1
i41           1
i42                       1
i43           1
i44           1
i45                                   1
i46                                               1
i47                                               1
i48           1
i49                                                           1
i50           1


----     56 VARIABLE r.L  rectangles

             x           y           w           h

k1     0.02634     0.69240     0.77969     0.27530
k2     0.05406     0.24001     0.20032     0.41687
k3     0.10543     0.02016     0.42409     0.12897
k4     0.36012     0.41666     0.19087     0.16088
k5     0.78532     0.07780     0.20850     0.61279


----     56 VARIABLE z.L                   =      0.51132  objective

This model was solved to proven global optimality with Gurobi in about 670 seconds. Cplex also can solve models of this type.

To higher dimensions

The above model dealt with the 2D problem. It is not very difficult to extend this to more dimensions. Let's try to make this general, for any dimension.

First of all, I reorganized my sets a bit:

• We have \(d \in \{1,\dots,D\}\) indicating the dimension we are looking at. Think about this as \(1=x,2=y,3=z,...\).
• To characterize a box, for each dimension, we have a coordinate and a side-length. So, \(cs \in \{\mathrm{coord},\mathrm{side}\}\).

With this, we can write:

• \(p_{i,d}\) are the coordinates of the points,
• \(r_{k,d,cs}\) indicate the location and size of the boxes.

The model can now look like:

High-level n-dimensional MINLP model

\[\begin{align}\min & \sum_k \prod_d \color{darkred} r_{k,d,side} \\ & \sum_k \color{darkred}x_{i,k} = 1 && \forall i  \\ & \color{darkred}x_{i,k}=1 \implies \begin{cases} \color{darkblue}p_{i,d} \ge \color{darkred}r_{k,d,coord} \\   \color{darkblue}p_{i,d} \le \color{darkred}r_{k,d,coord}+\color{darkred}r_{k,d,side}   \end{cases} && \forall i,k,d\\ & \color{darkred} x_{i,k}\in \{0,1\} \\ & \color{darkred}r_{k,d,cs} \in [\color{darkblue}L,\color{darkblue}U] \\ &\color{darkred}r_{k,d,coord}+\color{darkred}r_{k,d,side} \le \color{darkblue}U   \end{align}\]

The objective is not quadratic. We can make it quadratic by chaining[2]. We introduce a variable \(v\) which is defined by:\[v_{k,d} = \begin{cases} v_{k,d-1}\cdot r_{k,d,side} & \text{for $d\ge 3$}\\ r_{k,2,side} \cdot r_{k,1,side} & \text{for $d=2$}\end{cases}\] Now we can just minimize \(\sum_k v_{k,D}\). For the special case of a 3d problem, we can simplify this to: \[\begin{align} \min & \sum_k \mathit{area}_k \cdot r_{k,3,side} \\ & \mathit{area}_k = r_{k,2,side}\cdot r_{k,1,side}\end{align}\]

A small data set (\(N=20\) points, and \(K=5\)) leads to:

----     62 PARAMETER p  data points

           dim1        dim2        dim3

i1      0.80601     0.17283     0.52952
i2      0.14887     0.64826     0.69240
i3      0.35222     0.02016     0.43130
i4      0.55356     0.64087     0.77456
i5      0.23545     0.78074     0.26800
i6      0.08158     0.97328     0.11440
i7      0.87431     0.66711     0.75621
i8      0.96771     0.19863     0.24001
i9      0.21990     0.26101     0.98863
i10     0.17156     0.06625     0.93008
i11     0.80603     0.83235     0.10543
i12     0.02902     0.22927     0.09355
i13     0.12983     0.90322     0.43682
i14     0.72756     0.24754     0.57512
i15     0.36012     0.51628     0.71015
i16     0.74601     0.70367     0.74621
i17     0.18517     0.93583     0.81681
i18     0.67338     0.46350     0.57754
i19     0.08932     0.65688     0.97335
i20     0.69059     0.89449     0.07780


----     62 VARIABLE x.L  assignment variables

             k1          k2          k3          k4          k5

i1                                                            1
i2                        1
i3                                    1
i4                                                1
i5                        1
i6            1
i7                                                1
i8                                                            1
i9                                    1
i10                                   1
i11           1
i12           1
i13                       1
i14                                                           1
i15                                               1
i16                                               1
i17                       1
i18                                               1
i19                       1
i20           1


----     62 VARIABLE r.L  rectangles

    dim1.coord   dim1.side  dim2.coord   dim2.side  dim3.coord   dim3.side

k1     0.02902     0.77701     0.22927     0.74401     0.07780     0.03660
k2     0.08932     0.14613     0.64826     0.28757     0.26800     0.70534
k3     0.17156     0.18066     0.02016     0.24086     0.43130     0.55733
k4     0.36012     0.51418     0.46350     0.24017     0.57754     0.19701
k5     0.72756     0.24015     0.17283     0.07471     0.24001     0.33511


----     62 VARIABLE z.L                   =      0.10539  objective

This model can be solved with Gurobi. Cplex does not allow non-convex quadratic constraints.

3d solution

Linear formulation by set covering

This formulation was inspired by the comments by Rob Pratt below. The idea is to generate all possible rectangles in advance. We pick two arbitrary points to determine the horizontal part of the rectangle and two arbitrary points to determine the vertical part of the rectangle. This would mean \[{50 \choose 2} \times {50 \choose 2}= 1,500,625\] rectangles. In addition, we need to keep track of which points are covered by each rectangle \(k\). Let's introduce: \[s_{i,k} = \begin{cases} 1 & \text{if point $i$ is covered by rectangle $k$}\\ 0 & \text{otherwise}\end{cases}\] The resulting model is \[\begin{align} \min& \sum_k x_k \mathit{area}_k\\ & \sum_k x_k=5\\ & \sum_k s_{i,k} x_k \ge 1&& \forall i \\ & x_k \in \{0,1\} \end{align}\] This is now a very large, but very easy linear MIP model.

Notes:

• As we generate all rectangles in advance, we can also calculate the area (or volume) in advance. That explains why the model becomes linear.
• Note that the horizontal pair of points \(x_1,x_2\) can be identical to the vertical pair \(y_1,y_2\). They are drawn independently.
• We did not generate rectangles with just a single point. We can easily add them. 
• When generating two "horizontal" points \(x_1, x_2\) and two "vertical" points \(y_1,y_2\), one or more of the points may be outside the resulting box. We can just drop that candidate rectangle. See the picture below. Note that elsewhere in the generation process we will form the rectangle: \(y_1,x_2\) in the \(x\) direction and \(x_1,y_1\) in the \(y\) direction.  
• So: we generated more rectangles than needed if we did not drop cases like in the previous bullet. A little bit more logic is needed to prevent these unneeded rectangles. Note this will not change the solution. We just will generate a smaller model.  
• I think if you do this right, then just 50,658 rectangles will be generated. This means we could discard 96.6% of the candidate rectangles!
• This can be extended to larger dimensions. I expect the total number of possible boxes can become very large, very quickly.

This candidate rectangle is dropped from consideration.

References 

1. How can I find k minimum bounding rectangles to enclose all the given points?, https://stackoverflow.com/questions/60401171/how-can-i-find-k-minimum-bounding-rectangles-to-enclose-all-the-given-points
2. Non-convex Quadratic Models, http://yetanothermathprogrammingconsultant.blogspot.com/2020/02/non-convex-quadratic-models.html

How to model y=min(x1,x2)

Introduction

The problem of modeling \(y=\min(x_1,x_2)\) in optimization models is not completely trivial. Here are some notes.

Small example model

To experiment, here is a small test model:

MINLP Model
\[\begin{align} \max\>& \color{darkred} z= \color{darkred}x_1 + 2\color{darkred}x_2\\ & 2\color{darkred}x_1 + \color{darkred}x_2 = 5 + \min(\color{darkred}x_1,\color{darkred}x_2)\\ & \color{darkred}x_1,\color{darkred}x_2 \in [0,4] \end{align}\]

Global MINLP solvers

Interestingly the global MINLP solvers Baron, Couenne, and Antigone (under GAMS) have no direct facilities for this function:

Solver	Error message
Baron	*** Cannot handle function 'min'
Couenne	Gams function min not supported.
Antigone	ERROR: Unsupported function min in equation e0

This is unfortunate. I think it would be better if this function was supported directly. Of course, we are not defeated that easily. As the solvers do support the \(\mathrm{abs}()\) function, lets see if we can use that instead. A possible reformulation is to replace \(\min(x_1,x_2)\) by \[\min(x_1,x_2)=\frac{x_1+x_2-|x_1-x_2|}{2}\] The derivation is as follows: \[\frac{x_1+x_2-|x_1-x_2|}{2} = \begin{cases}\displaystyle\frac{x_1+x_2-x_1+x_2}{2}=x_2 & \text{if $x_1\ge x_2$}\\ \displaystyle\frac{x_1+x_2+x_1-x_2}{2}=x_1 & \text{if $x_1\lt x_2$}\end{cases}\] Using this we see:

                           LOWER          LEVEL          UPPER

---- VAR z                 -INF            9.0000        +INF         
---- VAR x1                  .             1.0000         4.0000      
---- VAR x2                  .             4.0000         4.0000      

MIP Formulation

The construct \(y=\min(x_1,x_2)\) can be interpreted as \[\begin{align} & y\le x_1 \> \mathbf{and}\> y\le x_2\\ & y\ge x_1 \>\mathbf{or}\> y\ge x_2 \end{align}\] With this, we can form the set of constraints: 

MIP big-M reformulation of min function
\[\begin{align} & \color{darkred} y \le \color{darkred}x_1 \\ & \color{darkred} y \le \color{darkred}x_2\\ & \color{darkred} y \ge \color{darkred}x_1 - \color{darkblue}M \color{darkred}\delta \\ & \color{darkred} y \ge \color{darkred}x_2 - \color{darkblue}M (1-\color{darkred}\delta)\\ & \color{darkred}\delta \in \{0,1\} \end{align}\]

The main problem here is that we need a good value for \(M\). In this case, \(M\) should be the largest difference between \(x_1\) and \(x_2\) we can observe. In our case \(x_1,x_2\in[0,4]\) so a good value would be \(M=4\). The solution can look like:

                           LOWER          LEVEL          UPPER         MARGINAL

---- VAR z                 -INF            9.0000        +INF             .          
---- VAR x1                  .             1.0000         4.0000          .          
---- VAR x2                  .             4.0000         4.0000         1.0000      
---- VAR y                   .             1.0000        +INF             .          
---- VAR delta               .              .             1.0000         EPS         

It may not always be so easy to derive good values. In that case, an alternative is to use SOS1 variables (SOS1 stands for Special Ordered Set of type 1). This type of variable is supported by most MIP solvers (but not all).

MIP SOS1 reformulation of the min function
\[\begin{align} & \color{darkred} y \le \color{darkred}x_1 \\ & \color{darkred} y \le \color{darkred}x_2\\ & \color{darkred} y \ge \color{darkred}x_1 - \color{darkred}s_1 \\ & \color{darkred} y \ge \color{darkred}x_2 - \color{darkred}s_2\\ & \color{darkred}s_1,\color{darkred}s_2 \ge 0 \\ & (\color{darkred}s_1,\color{darkred}s_1) \in SOS1\end{align}\]

The SOS1 condition says: only one of the variables \((s_1,s_2)\) can be non-zero (note that they can be both zero). This means: if one of the variables is non-zero, the other one will be zero. As you can see, there are no bounds needed on \(s_i\), so this formulation is useful when we can not find a good value for \(M\) in the previous formulation.

I want to mention that some solvers have a built-in function for \(\min()\). That makes life easier, as there is no reformulation required (the solver will do this behind the scenes).

How not to handle the min function

I see often the following suggestion to handle the above model:

Don't do this
\[\begin{align} \max\>& \color{darkred} z= \color{darkred}x_1 + 2\color{darkred}x_2 + \color{darkblue} \alpha \cdot\color{darkred}y \\ & 2\color{darkred}x_1 + \color{darkred}x_2 = 5 + \color{darkred}y \\ & \color{darkred}y \le \color{darkred}x_1 \\ & \color{darkred}y \le \color{darkred}x_2\\ & \color{darkred}x_1,\color{darkred}x_2 \in [0,4] \end{align}\]

Here \(\alpha\gt 0\) is a constant. The idea is that the objective will now push \(y\) upwards, so we don't need to worry anymore about \[y \ge x_1 \> \mathbf{or} \> y \ge x_2\] The problem is that we are changing the problem. Indeed the results (with \(\alpha=10\)) look like:

                           LOWER          LEVEL          UPPER         MARGINAL

---- VAR z                 -INF           32.5000        +INF             .          
---- VAR x1                  .             2.5000         4.0000          .          
---- VAR x2                  .             2.5000         4.0000          .          
---- VAR y                   .             2.5000        +INF             .          

This means the objective \(z_0=x_1+2x_2\) is equal to 7.5 which is less than the objective we saw earlier (we had before \(z=9\)). So this trick is just not the right thing to do.

There are values for \(\alpha\) that actually work. Using \(\alpha=0.1\) we observe:

                           LOWER          LEVEL          UPPER         MARGINAL

---- VAR z                 -INF            9.1000        +INF             .          
---- VAR x1                  .             1.0000         4.0000          .          
---- VAR x2                  .             4.0000         4.0000         0.9000      
---- VAR y                   .             1.0000        +INF             .          

This corresponds to the solution we found earlier. With some experimentation it looks like we have: \[\begin{cases} \alpha \le 1 \Rightarrow x_1=1,x_2=4\\ \alpha \gt 1 \Rightarrow x_1,x_2=2.5\end{cases}\] I am not sure if in general there is an easy way to find the threshold or even if there is always some value of \(\alpha\) that works. I don't particularly like this approach as it really changes the problem. As a result, I never use this method.

Octeract

I don't know what or who this is.

This seems to be a parallel deterministic solver for non-convex MINLPs.  Some things I noticed:

• !np.easy : cute (but nonsense of course: some problems just remain difficult).
• "The first massively parallel Deterministic Global Optimization solver for general non-convex MINLPs."
• Symbolic manipulation: like some other global solvers they need the symbolic form of the problem so they can reason about this. I.e. no black-box problems.
• Support for AMPL and Pyomo
•  "Octeract Engine has implementations of algorithms that guarantee convergence to a global optimum in finite time. The correctness of the majority of the calculations are ensured through a combination of interval arithmetic and infinite precision arithmetic."
• It looks like the benchmarks [3] are against itself (so always a winner). 
• "Massively parallel" makes me think about thousands of cores. The benchmarks do not seem to report anything like this (mostly 1,4 and 8 cores).
• I don't see any names on the web site. The About Company section is unusually vague. I received some more links with additional background info [4,5].

Some of the competing solvers are Baron, Couenne, and Antigone.

References

1. https://octeract.com/ 
2. Manual: https://octeract.com/wp-content/uploads/2019/08/user_manual.pdf
3. Benchmarks: https://octeract.com/benchmarks/
4. Presentation: https://www.youtube.com/watch?v=ayZ3iHmQUak
5. PhD thesis: https://spiral.imperial.ac.uk/bitstream/10044/1/45359/1/Kazazakis-N-2017-PhD-Thesis.pdf