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

Non-convex quadratic models

Some non-convex quadratic models I

The new version of Gurobi v 9 introduced support for non-convex quadratic models. Cplex already supported some classes of such models (such as non-convex QPs). Here are just some models we can now try to solve with these solvers.

It is noted that you need to specify a solver option for Cplex and Gurobi to accept most non-convex models, Some models can be reformulated by the solver, but "real" con-convex models require the option OptimalityTarget 3 (Cplex) or NonConvex 2 (Gurobi).

1a. Largest Empty Rectangle

Given is an area with \(n\) points:

n=50 random points

The problem is to find the largest empty rectangle, i.e., without any points inside. The solution looks like:

Largest empty rectangle

The idea behind the model below, it that a data point \(p_i\) must be either left, right, below or above the rectangle. So it the coordinates of the rectangle are denoted by \(r_{x,lo},r_{x,up},r_{y_lo},r_{y_up}\), then at least one of the following constraints must hold: \[\begin{align} & p_{i,x} \le r_{x,lo} \\& p_{i,y} \le r_{y,lo} \\ & p_{i,x} \ge r_{x,up} \\& p_{i,y} \ge r_{y,up}\end{align}\] The complete mathematical model to find this rectangle can be stated as:

Smallest empty rectangle
\[\begin{align} \max\> & \color{darkred}{\mathit{area}} = \prod_c (\color{darkred}r_{c,up}-\color{darkred}r_{c,lo}) \\ & \color{darkblue}p_{i,c} \le \color{darkred}r_{c,lo}+\color{darkblue}M_c(1-\color{darkred}\delta_{i,c,lo}) && \forall i,c\\ & \color{darkblue}p_{i,c} \ge \color{darkred}r_{c,up}-\color{darkblue}M_c(1-\color{darkred}\delta_{i,c,up}) && \forall i,c\\ & \sum_{c,s} \color{darkred}\delta_{i,c,s} \ge 1 && \forall i \\ & \color{darkblue} L_c \le \color{darkred}r_{c,lo} \le \color{darkred}r_{c,up}\le \color{darkblue} U_c \\ & \color{darkred}\delta_{i,c,s} \in \{0,1\} \end{align} \]
\[\begin{align} & i = \{1,\dots,n\}\\ & c = \{x,y\}\\ & s = \{lo,up\} \\ & \color{darkblue} L_c = 0 \\ & \color{darkblue} U_c = 1 \\ & \color{darkblue} M_c = \color{darkblue}U_c-\color{darkblue}L_c \end{align}\]

For the 2d case, the objective can be written as \[\min\> \mathit{area} = (r_{x,up}-r_{x,lo})\cdot (r_{y,up}-r_{y,lo})\] The outer area is a given rectangle with limits \(L_x, U_x, L_y, U_x\). I used here the unit box: \(L_x=L_y = 0\) and \(U_x=U_y = 1\). The big-M constants \(M_x, M_y\) need tight values. We can use the area dimensions for that. For a formulation without big-M constraints, we can implement implications using indicator constraints: \[\begin{align} & \delta_{i,c,lo}=1 \implies p_{i,c} \le r_{i,c,lo} \\ &\delta_{i,c,up}=1 \implies p_{i,c} \ge r_{i,c,up}\end{align}\]

The problem is not convex because of the objective. It is noted that we can simplify the objective a bit by writing \[\begin{align}\max\>&\mathit{area} =  \mathit{side}_x \cdot \mathit{side}_y\\ & \mathit{side}_c = r_{c,up} - r_{c,lo}\end{align}\] This is a small problem with just 4 continuous quadratic variables and \(n\times 4\) binary variables. Gurobi or Cplex solve this model very quickly.

This problem can also be solved with available global solvers. Both Antigone and Baron have no problems with it. The optimal solution looks like:

----     50 VARIABLE area.L                =        0.113  

----     50 VARIABLE r.L  coordinates of rectangle

           lo          up

x       0.403       0.792
y       0.285       0.575

I also tried the open-source solver Couenne. That was a surprise:

               S O L V E      S U M M A R Y

     MODEL   m                   OBJECTIVE  area
     TYPE    MIQCP               DIRECTION  MAXIMIZE
     SOLVER  COUENNE             FROM LINE  49

**** SOLVER STATUS     1 Normal Completion         
**** MODEL STATUS      1 Optimal                   
**** OBJECTIVE VALUE                0.0121


----     50 VARIABLE area.L                =        0.012  

----     50 VARIABLE r.L  coordinates of rectangle

           lo          up

x       0.939       0.982
y       0.720       1.000

This is obviously a bad and non-optimal solution. I don't know what is going on here, but this looks like a bug.

1b. 3d largest box 

It is easy to make a 3d version of the model above. Just replace \(c = \{x,y\}\) by \(c = \{x,y,z\}\). Unfortunately, there is one complication: now the model has an objective that is no longer quadratic. However, this can be fixed by chaining: \[\begin{align} \max\>&\mathit{volume} = \mathit{area}\cdot \mathit{side}_{z}\\ & \mathit{area} = \mathit{side}_{x} \cdot \mathit{side}_{y} \\ & \mathit{side}_c = r_{c,up}-r_{c,lo} && \forall c \end{align}\]

A solution can look like:

----     52 VARIABLE area.L                =        0.340  
            VARIABLE volume.L              =        0.146  

----     52 VARIABLE r.L  coordinates of 3D box

           lo          up

x       0.608       0.998
y       0.091       0.963
z       0.140       0.571

Note that Cplex can not handle this problem: Cplex allows non-convexities in the objective only (apart for some cases where it can reformulate things).

3d version: largest box not clashing with given points

2. Hostile Brothers Problem

A father is facing the following problem. He has has a plot of land. And he has \(n\) sons, who don't get along at all. In order to minimize friction, he wants to build a house for each of his sons, such that the minimum distance between them is maximized. The mathematical model can look like:

Hostile Brothers Problem
\[\begin{align} \max\> & \color{darkred}z \\ & \color{darkred}z \le \sum_{c} (\color{darkred}x_{i,c}-\color{darkred}x_{j,c})^2&& \forall i\lt j\\ & \color{darkred}x_{i,c} \in [\color{darkblue}L,\color{darkblue}U]\end{align} \]

Here \(i\) indicates the point, \(i \in \{1,\dots,n\}\) and \(c\) is the coordinate. For a 2D problem, we have \(c \in \{x,y\}\).

When we use \(n=50\) random points, Gurobi is not able to even come close to proven optimality, I used a time limit of 2 hours, and the final gap is enormous:

 1560307 1354785    0.40487   42 1653    0.01935    0.67173  3371%   332 7171s
 1562085 1356368     cutoff  283         0.01935    0.67173  3371%   332 7177s
 1563023 1356980     cutoff  285         0.01935    0.67173  3371%   332 7180s
 1564565 1358071    0.01937  248  930    0.01935    0.67173  3371%   332 7190s
 1565720 1359205    0.60071   48 1510    0.01935    0.67173  3371%   332 7196s
 1567335 1359826     cutoff  272         0.01935    0.67173  3371%   332 7200s


Cutting planes:
  RLT: 1312

Explored 1567341 nodes (520860935 simplex iterations) in 7200.00 seconds
Thread count was 8 (of 16 available processors)

Summary:

• Best found solution: 0.01935
• Best possible objective value: 0.6173
• Relative gap: 3371%

The big gap is partly due to the small best-found solution (0.01935). But mainly: this is a very difficult job for Gurobi. Even when the problem is relatively small: just 100 continuous variables,

By inspecting the solution, we see this solution is not very good at all:

Solution found by Gurobi after 2 hours (obj=0.01935)

I expect a much more regular pattern.

We can try to improve performance by noticing there is much symmetry in the model. The same solution can be obtained by just renumbering the points. To break symmetry, I added the constraint: \[\sum_c x_{i,c} \ge \sum_c x_{i-1,c}\] This was not very successful:

 2978972 2407171    0.01708  130  813    0.00885    0.64271  7159%   222 7180s
 2981109 2409232    0.52720   64 2378    0.00885    0.64271  7159%   222 7185s
 2984648 2411500    0.20399  107 2156    0.00885    0.64271  7159%   222 7191s
 2986302 2412811    0.12220  127 1721    0.00885    0.64271  7159%   222 7195s
 2987852 2413879    0.06942  175 1728    0.00885    0.64271  7159%   222 7200s

Explored 2988184 nodes (662244035 simplex iterations) in 7200.01 seconds
Thread count was 8 (of 16 available processors)

So we see:

• The best-found solution (0.00885) is worse than before (0.01935)
• The best possible solution (0.64271) is worse than before (0.6173)
• The relative gap is worse (7159% vs 3371%)

As an alternative, let's try Baron on this model. I dropped the symmetry breaking constraint here (adding this constraint was not beneficial, at least not during the initial period of the search). We see in the first few minutes:

  Iteration    Open nodes         Time (s)    Lower bound      Upper bound
*         1             1             5.51     0.257746E-01     2.00000    
*         1             1             7.22     0.274782E-01     2.00000    
          1             1             7.41     0.274782E-01     2.00000    
         12             7            49.51     0.274782E-01     1.72727    
         31            16            80.03     0.274782E-01     1.72727    
         46            24           111.27     0.274782E-01     1.72727    

I.e. Baron finds a much better solution very quickly (but is far from getting the gap close to zero).

Another approach is to use a multi-start approach with a local NLP solver. Generate a random starting point, solve, and repeat. When I do this 20, I see when I use the CONOPT NLP solver:

----     57 PARAMETER objs  objective values of NLP solver

cycle1  0.02751,    cycle2  0.02669,    cycle3  0.02705,    cycle4  0.02041,    cycle5  0.02620,    cycle6  0.01235
cycle7  0.02675,    cycle8  0.02590,    cycle9  0.02660,    cycle10 0.02679,    cycle11 0.02558,    cycle12 0.02644
cycle13 0.02622,    cycle14 0.02625,    cycle15 0.02620,    cycle16 0.02749,    cycle17 0.02742,    cycle18 0.02692
cycle19 0.02675,    cycle20 0.02688,    best    0.02751

This is the best we have seen so far. By accident, the best results were found in the first cycle. The results look like:

Best solution with Multi-start Conopt (obj=0.02751)

This solution looks much better than the one Gurobi found. In addition, CONOPT needs only 30 seconds to solve these 20 problems.

Conclusion:

• Cplex cannot solve this problem
• Gurobi finds a solution with objective 0.01935 (time limit: 2 hours)
• Baron finds a solution of  0.02748 in 100 seconds
• Conopt finds a solution of 0.02751 in a loop of 20 solves that takes 30 seconds
• Increasing this to 200 solves gives an objective of 0.02771
• Remember that we are maximizing

Even though this is a non-convex quadratic problem, it makes sense to look beyond quadratic solvers.

Nonlinear variant of a knapsack problem

In [1] a problem is posed:

Original Problem
\[\begin{align}\max & \sum_i \log_{100}(\color{darkblue}w_i) \cdot \color{darkred}x_i \\ & \frac{\sum_i \color{darkblue}w_i \color{darkred}x_i}{\sqrt{\color{darkred}k}} \le 10,000\\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \color{darkred}x_i \in \{0,1\} \end{align}\]

The data vector \(w_i\) is assumed to integer valued with \(w_i\ge 1\). Hence the logarithm can be evaluated without a problem. Also we can assume that the number of items can be around 1,000.

I don't think I have ever seen \(\log_{100}()\) being used in a model. Most programming (and modeling) languages only support natural logarithms \(\ln()\) and may be \(\log_{10}()\). We can convert things by: \[\log_{100}(x) = \frac{\ln(x)}{\ln(100)}\] This means the objective can be written as \[\max \frac{1}{\ln(100)} \sum_i \ln(w_i) x_i\] Essentially, the \(\log_{100}()\) function just adds a scaling factor. We can simplify the objective to \[\max \sum_i \ln(w_i)x_i\] (The objective value will be different, but the optimal solution will be the same). Note that this is a linear objective.

log100(x) vs ln(x)

In this post I'll discuss a few different formulations:

• As as MINLP (Mixed Integer Nonlinear Programming) model. This is an obvious choice considering the square root in the model.
• An MIQCP (Mixed Integer Quadratically Constrained Programming) model. This (convex) quadratic reformulation allows us to use alternative solvers.
• A linear MIP model. This formulation may not be immediately obvious. Of course with this model we can use widely available MIP solvers.

MINLP Formulation

The constraint can be rewritten as: \[\sum_i w_i x_i \le 10,000 \sqrt{k}\] If we ignore the all zero solution, we can assume \(k\ge 1\). This bound will make sure the square root function is always differentiable. With this, we have a standard MINLP model.

MINLP Model
\[\begin{align}\max & \sum_i \ln(\color{darkblue}w_i) \color{darkred}x_i \\ & \sum_i \color{darkblue}w_i \color{darkred}x_i \le 10,000\sqrt{\color{darkred}k}\\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \color{darkred}x_i \in \{0,1\} \\ & \color{darkred} k \ge 1 \end{align}\]

This model has only a single, well behaved, non-linearity. Literally, there is only one nonlinear nonzero element.  In addition, the model is convex. So we don't expect any problems. For \(w_i\),  I generated 1000 random integer values from the interval \([1,10000]\). All timings are in seconds.

MINLP Results
Solver	Obj	Time	Notes
Dicopt	1057.1355	0.6	3 NLP, 2 MIP subproblems
SBB	1056.9361	95	Node limit exceeded
Bonmin	1056.9472	600	Time limit exceeded
Bonmin	1057.1355	16	Option: bonmin.algorithm B-OA

The outer-approximation based algorithms (Dicopt, Bonmin with B-OA option) do much better than the branch & bound algorithms (SBB, default Bonmin). Even the global solvers do better:

Global Solver Results
Solver	Obj	Time
Baron	1057.1355	2
Antigone	1057.1355	1
Couenne	1057.1355	9

MIQCP Formulation

The problem can also be formulated as a convex MIQCP model:

MIQCP Model
\[\begin{align}\max & \sum_i \ln(\color{darkblue}w_i) \color{darkred}x_i \\ & \color{darkred}y^2 \le \color{darkred}k\\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \color{darkred}y = \frac{\sum_i \color{darkblue}w_i \color{darkred}x_i}{10,000} \\ & \color{darkred}x_i \in \{0,1\} \\ & \color{darkred} k \ge 1, \color{darkred}y \ge 0 \end{align}\]

Cplex has two ways to form relaxations:

• MIQCPstrat 1: SOCP-based branch & bound
• MIQCPstrat 2: OA-based branch & cut

The results are shown in the next section. We see again that the Outer Approximation approach works better.

MIP Formulation

Finally, we can also linearize this model by observing that \(k\in \{1,\dots,1000\}\). So we don't really have a continuous function \(f(k)=\sqrt{k}\), but rather only need function values at the integer points. We can exploit this by making this explicit. We can write: \[\begin{align} & k = \sum_i i\cdot \delta_i \\ & \sqrt{k} = \sum_i \sqrt{i}\cdot \delta_i \\ & \sum_i \delta_i = 1 \\ & \delta_i \in \{0,1\}\end{align}\] This is essentially a SOS1 (Special Ordered Set of Type 1) structure implementing a table lookup. The MIP model looks like:

MIP Model

\[\begin{align}\max & \sum_i \ln(\color{darkblue}w_i) \color{darkred}x_i \\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \sum_i \color{darkblue}w_i \color{darkred} x_i \le 10,000 \color{darkred}q \\ & \color{darkred} k = \sum_i i \cdot \color{darkred}\delta_i \\ & \color{darkred} q = \sum_i \sqrt{i} \cdot \color{darkred}\delta_i \\ & \sum_i \color{darkred} \delta_i = 1 \\ & \color{darkred}x_i, \color{darkred}\delta_i \in \{0,1\} \\ & \color{darkred} k, \color{darkred} q \ge 1 \end{align}\]

When we solve this problem we see:

MIQCP and MIP Results
Model	Solver	Obj	Time	Notes
MIQCP	Cplex	1057.1355	600	SOCP-based branch & bound, time limit
MIQCP	Cplex	1057.1355	0.3	OA-based branch & cut
MIP	Cplex	1057.1355	0.6

Conclusion: the question in the original post was: how to solve this problem? Here we proposed three different models: an MINLP, MIQCP and MIP model. All these models can be solved quickly. It is noted that the quadratic and linear models are not approximations: they give the same solution as the original MINLP model. Pure non-linear branch & bound methods are having a bit of a problem with the MINLP model, but Outer-Approximation works very well.

References

1. How do we solve a variant of the knapsack problem in which the capacity of the knapsack keeps increasing as we add more items into the knapsack?, https://stackoverflow.com/questions/59242370/how-do-we-solve-a-variant-of-the-knapsack-problem-in-which-the-capacity-of-the-k
2. Pierre Bonami, Andrea Tramontani, Recent improvement to MISOCP in CPLEX, Informs 2015, http://www.optimizationdirect.com/data/informs2015_bonami.pdf

Appendix: automatic reformulation with CVXPY

CVXPY accepts the MINLP model, and can automatically reformulate this into a MISOCP model:

x = cp.Variable(n,integer=True)
obj = cp.Maximize(np.log(w)*x)
con = w*x <= 10000 * cp.sqrt(cp.sum(x)) 
prob = cp.Problem(obj,[con, x >= 0, x <= 1])
prob.solve(solver=cp.ECOS_BB, max_iters=10000)
print("status:",prob.status)
print("Obj:",prob.value)

Unfortunately, the MISOCP solver that comes with CVXPY is not the most reliable one. It shows:

status: optimal_inaccurate
Obj: 1057.1085712222387
X: [-1.52640229e-10 -4.44776117e-11 -6.08671797e-11 -9.71628737e-11
 -9.95050649e-11 -1.23087763e-10 -8.63183870e-11 -4.39955343e-11
  1.00000000e+00 -6.54372650e-11 -3.95051292e-11 -5.86051991e-11
...

The reported objective is close to what we have found earlier.