Monday, March 2, 2020

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

8 comments:

  1. Alternatively, you can model the problem linearly as a cardinality-constrained set covering problem, as follows. Let R be the set of all possible rectangles with a grid point on each boundary. Here, |R| = (50 choose 2)^2 = 1,500,625. For r in R, let binary variable U(r) indicate whether rectangle r is used. The objective is to minimize the total area sum {r in R} area(r) * U(r). The cardinality constraint is sum {r in R} U(r) = 5. The cover constraint for each point i is sum {r in R: i in R} U(r) >= 1. This formulation avoids the symmetry issue, and the resulting optimal solution matches yours.

    ReplyDelete
    Replies
    1. Great idea. I think the number of rectangle is not 100% correct. We miss rectangles with a single point. Also some rectangles may not be of value, because one of the two "y" points (or both) may be outside the "x" window. And vice versa.

      Delete
    2. This comment has been removed by the author.

      Delete
    3. Yes, you can omit any rectangles that do not contain any points. In that case, I get 1,342,573 nonempty rectangles; I don't understand how to get 50,658. You can also allow side length 0, in which case I get 1,386,723 nonempty but possibly 0 area rectangles.

      Delete
    4. Yes, I confirm 50,658 rectangles now, including those with 0 area. I wrote "with a grid point on each boundary" but had not checked this condition.

      Delete
    5. This comment has been removed by the author.

      Delete
    6. For the 3D instance, I get 6039 useful boxes.

      Delete
  2. This comment has been removed by the author.

    ReplyDelete