Showing posts with label Set Covering Model. Show all posts
Showing posts with label Set Covering Model. Show all posts

Monday, September 14, 2020

Cover points with tiles

In [1] a problem is posted:

Given \(n\) points and a collection of tiles (rectangles) of given size, try to cover all all points with the minimum number of tiles.

There are a number of things to consider:

  • Are the tiles of equal size or have they different sizes?
  • Can tiles overlap?
  • After minimizing the number of tiles, should we minimize the total size of the selected tiles, so we select the smaller tiles when possible.

Data


Let's generate some data, assuming tiles have different sizes:


----     19 PARAMETER p  location of points

              x           y

p1       17.175      84.327
p2       55.038      30.114
p3       29.221      22.405
p4       34.983      85.627
p5        6.711      50.021
p6       99.812      57.873
p7       99.113      76.225
p8       13.069      63.972
p9       15.952      25.008
p10      66.893      43.536
p11      35.970      35.144
p12      13.149      15.010
p13      58.911      83.089
p14      23.082      66.573
p15      77.586      30.366
p16      11.049      50.238
p17      16.017      87.246
p18      26.511      28.581
p19      59.396      72.272
p20      62.825      46.380
p21      41.331      11.770
p22      31.421       4.655
p23      33.855      18.210
p24      64.573      56.075
p25      76.996      29.781
p26      66.111      75.582
p27      62.745      28.386
p28       8.642      10.251
p29      64.125      54.531
p30       3.152      79.236


----     19 PARAMETER r  size of rectangles

              x           y

r1       12.183      15.270
r2       25.769      32.506
r3       15.344      11.024
r4       27.554      28.637
r5       21.681      20.761
r6       17.291      17.393
r7       13.915      38.003
r8       21.398      33.502
r9       19.001      13.764
r10      32.466      12.077


Random points

Model


To model this problem, we introduce binary variables \(\mathit{select}_j \in \{0,1\}\): \[\mathit{select}_j = \begin{cases} 1 & \text{if tile $j$ is selected}\\ 0 & \text{otherwise}\end{cases}\] and a continuous variable \(\mathit{place}_{j,c}\) indicating the \(x\) and \(y\) coordinates of (selected) tile \(j\). To be precise: it is the left-lower corner of the rectangle. The set \(c\) is simply \(c=\{x,y\}\). In addition we need a variable indicating if a point is covered by a tile:\[\mathit{cov}_{i,j}=\begin{cases} 1 & \text{if point $i$ is covered by tile $j$} \\ 0 &\text{otherwise}\end{cases} \]

The model can look like:

MIP Model
\[\begin{align}\min &\sum_j \color{darkred}{\mathit{select}}_j \\ & \color{darkred}{\mathit{place}}_{j,c} \le \color{darkblue}p_{i,c}+\color{darkblue}M\cdot(1-\color{darkred}{\mathit{cov}}_{i,j}) && \forall i,j,c \\ &\color{darkred}{\mathit{place}}_{j,c}+\color{darkblue}r_{j,c}\ge\color{darkblue}p_{i,c}-\color{darkblue}M\cdot(1-\color{darkred}{\mathit{cov}}_{i,j}) && \forall i,j,c \\ & \color{darkred}{\mathit{cov}}_{i,j} \le \color{darkred}{\mathit{select}}_j && \forall i,j \\ & \sum_j \color{darkred}{\mathit{cov}}_{i,j} \ge 1 && \forall i \\ & \color{darkred}{\mathit{select}}_j \in \{0,1\} \\ & \color{darkred}{\mathit{cov}}_{i,j} \in \{0,1\} \end{align}\]


Our points are drawn from \[p_{i,c}\sim U(0,100)\] So it make sense to restrict the placement of tiles to: \[0 \le  \mathit{place}_{j,c} \le 100-r_{j,c}\] Using this, a good value for \(M\) is \(M=100\).

Results



----     55 VARIABLE z.L                   =        7.000  objective

----     55 PARAMETER tile  placement of tiles

              x           y           w           h

r1       53.928      67.819      12.183      15.270
r2        5.652                  25.769      32.506
r4       27.484       6.507      27.554      28.637
r5       78.131      55.464      21.681      20.761
r6        5.791      49.181      17.291      17.393
r8       56.188      22.573      21.398      33.502
r10       2.517      75.169      32.466      12.077


Points are covered by tiles


No overlap


We can add no-overlap constraints [2] to the model soo that no pair of tiles can overlap. Just to be sure, I'll allow tiles to go outside the area \([0,100]\times [0,100]\). Of course, that also has an impact on the value for \(M\). The new solution looks like:


With no-overlap constraints


Same-sized tiles


If all tiles are identical, we just have a special case of our original problem. It may make sense to add an ordering constraint: \[\mathit{select}_j \le \mathit{select}_{j-1}\] We can also consider additional symmetry-breaking constraints, such as \[\mathit{place}_{j,x} \ge \mathit{place}_{j,x-1}\] It requires some experimentation if these constraints pay off in performance.


All tiles have the same size

In this example all the times are \((20 \times 20)\). Below is a larger example with \(n=50\) points and \((10 \times 10)\) tiles.

Larger example


Set-covering formulation for equal-sized tiles


If the problem has all equal-sized tiles, we can employ a very different algorithm.  For this we need to generate candidate locations in advance.  One way to generate those is as follows:

  1. \(k:=0\)
  2. For each point with coordinates \((x_1,y_1)\):
    1. For each point with coordinates \((x_2,y_2)\) such that \(x_1 \le x_2 \le x_1+w\) and \(y_1-h\le y_2\le y_1+h\):
      1. Form the candidate rectangle \((x_1,y_2)-(x_1+w,y_2+h)\)
      2. \(k := k+1\)
      3. Store as \(\mathit{loc}_{k,c}=(x_1,y_2)\)
      4. Find all points \(i\) inside the candidate rectangle, and store this inside a mapping \(\mathit{contain}_{k,i}\)
Now we can solve the set-covering problem [3]:

Set-covering Model
\[\begin{align}\min & \sum_j \color{darkred}{\mathit{select}}_j \\ & \sum_j \color{darkblue}{\mathit{contain}}_{j,i} \cdot   \color{darkred}{\mathit{select}}_j  \ge 1 && \forall i \\ & \color{darkred}{\mathit{select}}_j \in \{0,1\}\end{align}\]

where \[\mathit{contain}_{j,i} = \begin{cases} 1 & \text{if candidate rectangle $j$ contains point $i$}\\ 0 & \text{otherwise}\end{cases}\]


This model solves very fast.


A multi-criteria version


Instead of just minimizing the number of tiles, we can subsequently minimize the total area covered by the tiles. I.e. use smaller tiles if possible. This is essentially a multi-objective problem. We used here a lexicographic approach (i.e. do two solves in sequence). In our experiment, we reduced the total size of the selected tiles from 3672.5809 to 3004.0812.


Multi-objective model results


Obviously, this additional step makes only sense if the tiles have different sizes. 


GAMS Models


As requested, here are some of the source files.

Model for tiles with different sizes.

set
  i
'points to cover'  /p1*p30/
  j
'rectangles' /r1*r10/
  c
'coordinates' /x,y/
;

parameter
   A
'side of area' / 100 /
   M
'big-M'
   p(i,c)
'location of points'
   r(j,c)
'size of rectangles'
;
p(i,c) = uniform(0,A);
r(j,c) = uniform(10,40);
M = A;
display p,r;

binary variable select(j)  'select tile';
positive variable place(j,c) 'placement of tile';
place.up(j,c) = A-r(j,c);
variable z 'objective';
binary variable cov(i,j);

equations
   obj
   cover1
   cover2
   cover3
   cover4
;

obj..  z =e=
sum(j, select(j));
cover1(i,j,c)..  place(j,c) =l= p(i,c)+M*(1-cov(i,j));
cover2(i,j,c)..  place(j,c)+r(j,c) =g= p(i,c)-M*(1-cov(i,j));
cover3(i,j)..  cov(i,j) =l= select(j);
cover4(i)..
sum(j,cov(i,j)) =g= 1;

model m1 /all/;
option optcr=0, threads=8;
solve m1 minimizing z using mip;

select.l(j) = round(select.l(j));

parameter tile(j,*) 'placement of tiles';
tile(j,c) = place.l(j,c)*select.l(j);
tile(j,
'w')$select.l(j) =  r(j,'x');
tile(j,
'h')$select.l(j) =  r(j,'y');
display z.l,tile;



Set-covering model for equal-sized tiles.

set
  i
'points to cover'  /p1*p50/
  j
'candidate rectangles' /r1*r1000/
  contain(j,i)
'points contained in rectangle j'
  c
'coordinates' /x,y/
;

parameter
   A
'side of area' / 100 /
   p(i,c)
'location of points'
   r(c)
'size of rectangles'  /x 10, y 10/
   loc(j,c) 
'locations of candidate rectangles'
;
p(i,c) = uniform(0,A);
display p,r;

*------------------------------------------------------------
* generate candidate rectangle locations
*------------------------------------------------------------

alias (i,ii,iii);
scalar k / 0 /;

loop(i,
  
loop(ii$(p(ii,'x') >= p(i,'x') and p(ii,'x') <= p(i,'x')+r('x') and
               abs(p(ii,
'y')-p(i,'y')) <= r('y')),
      k = k + 1;
     
loop(j$(ord(j)=k),
         loc(j,
'x') = p(i,'x');
         loc(j,
'y') = p(ii,'y')-r('y');
         contain(j,iii) = p(iii,
'x') >= loc(j,'x') and p(iii,'x') <= loc(j,'x')+r('x') and
                          p(iii,
'y') >= loc(j,'y') and p(iii,'y') <= loc(j,'y')+r('y')
      );
   );
);

*------------------------------------------------------------
* set covering model
*------------------------------------------------------------

binary variable select(j)  'select tile';
variable z 'objective';

equations
   obj
   cover(i)
;

obj..  z =e=
sum(j, select(j));
cover(i).. 
sum(contain(j,i),select(j)) =g= 1;

model m1 /all/;
option optcr=0;
solve m1 minimizing z using mip;

*------------------------------------------------------------
* reporting
*------------------------------------------------------------

parameter tile(j,*) 'placement of tiles';
tile(j,c) = loc(j,c)*select.l(j);
tile(j,
'w')$select.l(j) =  r('x');
tile(j,
'h')$select.l(j) =  r('y');
display z.l,tile;


The R code to plot the results looks like:

library(ggplot2)

pt <- read.table(text="p x y
p1       17.175      84.327
p2       55.038      30.114
p3       29.221      22.405
p4       34.983      85.627
p5        6.711      50.021
p6       99.812      57.873
p7       99.113      76.225
p8       13.069      63.972
p9       15.952      25.008
p10      66.893      43.536
p11      35.970      35.144
p12      13.149      15.010
p13      58.911      83.089
p14      23.082      66.573
p15      77.586      30.366
p16      11.049      50.238
p17      16.017      87.246
p18      26.511      28.581
p19      59.396      72.272
p20      62.825      46.380
p21      41.331      11.770
p22      31.421       4.655
p23      33.855      18.210
p24      64.573      56.075
p25      76.996      29.781
p26      66.111      75.582
p27      62.745      28.386
p28       8.642      10.251
p29      64.125      54.531
p30       3.152      79.236
",header=T)


r <- read.table(text="r x y w h
r1       53.928      67.819      12.183      15.270
r2       15.562       2.638      25.769      32.506
r4       50.032      27.438      27.554      28.637
r5        7.540       4.247      21.681      20.761
r6        5.791      49.181      17.291      17.393
r7       85.897      38.222      13.915      38.003
r10       2.517      75.169      32.466      12.077
",header=T)


ggplot(pt, aes(x=x , y=y)) + geom_point() + 
  geom_rect(data=r, aes(xmin=x,xmax=x+w,ymin=y,ymax=y+h,color=r,fill=r,alpha=0.1))+
  guides(alpha=FALSE)


Conclusion


The small MIP model presented above can be used as a basis for alternative, related problems. For equal-sized tiles, a faster set-covering model is shown.

References




Sunday, June 7, 2020

A facility location problem: debunking a myth

In [1] the following problem is proposed:

Consider a grid with 40,000 cells. We have customers located on this grid (this is data). We want to place stores on the grid such that each customer has at least one store within reach. Finally, we want to minimize the number of stores we need to use.

Furthermore:

so in a situation like mine I could easily end up with a linear system with millions of rows and columns, which with integer/binary variables will take about the age of the universe to solve.

Well, let's see...


Data


The first thing we do is generating some data. We introduce a grid of  \(200 \times 200\) with \(40,000\) cells. Then we randomly place \(M=500\) customers on this grid.



----     22 PARAMETER cloc  customer locations

                  x           y

cust1            35          75
cust2           169          84
cust3           111          18
cust4            61         163
cust5            59         102
cust6            45         147
cust7            70         165
cust8           172          83
cust9            14         185
cust10          101          79
. . .
cust495          83         186
cust496         174         159
cust497         196         148
cust498         115         136
cust499          63         101
cust500          92          87


Reach


Secondly, we create a large, sparse boolean parameter \(\mathit{reach}_{c,i,j}\) indicating if cell \((i,j)\) is within reach of customer \(c\). I used the rule: if the Manhattan distance between \(c\) and \((i,j)\) is 10 or less, we are within reach: \[ |i-\mathit{cloc}_{c,x}|+|j-\mathit{cloc}_{c,y}| \le 10\] This generates a large parameter. A fragment can look like:



----     29 PARAMETER reach  location (i,j) is within reach of customer

                      j1          j2          j3          j4          j5          j6          j7          j8          j9

cust3  .i110                                                                                                           1
cust3  .i111                                                                                               1           1
cust3  .i112                                                                                                           1
cust39 .i116                                                                                   1
cust39 .i117                                                                       1           1           1
cust39 .i118                                                           1           1           1           1           1
cust39 .i119                                               1           1           1           1           1           1
cust39 .i120                                   1           1           1           1           1           1           1
cust39 .i121                       1           1           1           1           1           1           1           1

There are 106,376 elements in this data structure.


MIP model


With this, we can formulate the set covering model:


Set covering model
\[\begin{align}\min\>& \color{darkred}{\mathit{numStores}} \\ & \color{darkred}{\mathit{numStores}} = \sum_{i,j} \color{darkred}{\mathit{placeStore}}_{i,j}\\ & \sum_{i,j} \color{darkblue}{\mathit{reach}}_{c,i,j} \cdot \color{darkred}{\mathit{placeStore}}_{i,j} \ge 1 & \forall c \\ &\color{darkred}{\mathit{placeStore}}_{i,j} \in \{0,1\} \end{align}\]


This is a large, but easy MIP model. The size is large:


MODEL STATISTICS

BLOCKS OF EQUATIONS           2     SINGLE EQUATIONS          501
BLOCKS OF VARIABLES           2     SINGLE VARIABLES       40,001
NON ZERO ELEMENTS       146,377     DISCRETE VARIABLES     40,000

But it solves very fast:

Found incumbent of value 495.000000 after 0.00 sec. (1.89 ticks)
Tried aggregator 8 times.
MIP Presolve eliminated 306 rows and 39764 columns.
MIP Presolve modified 6 coefficients.
Aggregator did 29 substitutions.
Reduced MIP has 166 rows, 208 columns, and 658 nonzeros.
Reduced MIP has 208 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.36 sec. (72.59 ticks)
Probing time = 0.00 sec. (0.05 ticks)
Tried aggregator 1 time.
Detecting symmetries...
MIP Presolve eliminated 1 rows and 5 columns.
Reduced MIP has 165 rows, 203 columns, and 644 nonzeros.
Reduced MIP has 203 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.42 ticks)
Probing time = 0.00 sec. (0.05 ticks)
Clique table members: 35.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 8 threads.
Parallel mode: deterministic, using up to 2 threads for parallel tasks at root LP.
Tried aggregator 1 time.
No LP presolve or aggregator reductions.
Presolve time = 0.00 sec. (0.09 ticks)

Iteration log . . .
Iteration:     1   Dual objective     =            57.000000
Iteration:    62   Dual objective     =            93.000000
Iteration:   124   Dual objective     =           109.000000
Iteration:   186   Dual objective     =           111.500000
Root relaxation solution time = 0.06 sec. (0.94 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                          495.0000       56.0000            88.69%
Found incumbent of value 495.000000 after 0.44 sec. (77.77 ticks)
*     0+    0                          118.0000       56.0000            52.54%
Found incumbent of value 118.000000 after 0.44 sec. (78.15 ticks)
      0     0      111.7778    78      118.0000      111.7778      203    5.27%
*     0+    0                          116.0000      111.7778             3.64%
Found incumbent of value 116.000000 after 0.45 sec. (79.32 ticks)
      0     0      112.5000    78      116.0000      Cuts: 10      209    3.02%
*     0+    0                          113.0000      112.5000             0.44%
Found incumbent of value 113.000000 after 0.52 sec. (85.02 ticks)
      0     0        cutoff            113.0000      112.5000      209    0.44%
Elapsed time = 0.55 sec. (85.03 ticks, tree = 0.01 MB, solutions = 3)

Zero-half cuts applied:  5
Lift and project cuts applied:  2
Gomory fractional cuts applied:  3

Root node processing (before b&c):
  Real time             =    0.55 sec. (86.90 ticks)
Parallel b&c, 8 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.55 sec. (86.90 ticks)
MIP status(101): integer optimal solution


We can see that the presolver is especially successful here: the model size is reduced by an extremely large amount. This will be less pronounced when the number of customers increases.

Solution


The solution looks like:



----     51 VARIABLE numStores.L           =          113  number of stores

----     51 VARIABLE placeStore.L  store locations

              j1          j6          j7          j8          j9         j15         j16         j17         j18

i4                                                                                                             1
i18                                                            1
i40                                                                                    1
i70            1
i79                                                                                                1
i80                        1
i107                                                                                   1
i118                                   1
i136                                                                       1
i157                                               1
i167                                                                       1
i193                                               1

   +         j21         j23         j26         j28         j29         j31         j32         j36         j38

i10                                    1
i28                                                                        1
i54                                                            1
i72                                                                                    1
i96                                                1
i113                       1
i147                                                           1
i158                                                                                                           1
i179                                                                                               1
i184                       1
i198           1

   +         j39         j44         j45         j46         j49         j50         j56         j58         j59

i5                                                                                                             1
i18                                                1
i39            1
i62                                                            1
i85                        1
i102                                                                       1
i104                                   1
i133                                               1
i166                                                                                               1
i195                                                                                   1

   +         j62         j66         j67         j68         j69         j73         j74         j76         j80

i11                                                                                                1
i16                                                1
i36                                                            1
i61                                                                                    1
i76                                                                        1
i105                                                                       1
i112           1
i117                                                                                                           1
i128                                   1
i146                                               1
i190                       1

   +         j82         j84         j85         j88         j90         j92         j95         j96         j97

i17                                    1
i26            1
i35                                                                                                1
i48                                    1
i68                                                                        1
i79                                                                                                            1
i97            1
i136                       1
i156                                               1
i170                                   1
i183                                                                                   1
i191                                                           1

   +         j98        j102        j107        j111        j112        j114        j115        j116        j118

i4                                                                         1
i22                                                                                    1
i36                                                            1
i56                        1
i63                                                1
i68                                                                                                            1
i88                                                                                                1
i100                                                                                   1
i101           1
i111                                                                                                           1
i129                                                           1
i140                                   1

   +        j119        j121        j126        j127        j132        j133        j134        j136        j139

i11                                                            1
i30                                    1
i53                                                1
i72                                                                                                1
i111                                                                                                           1
i129                                                                                   1
i144                                                                       1
i159           1
i183                                                                                               1
i191                       1

   +        j140        j147        j149        j150        j152        j153        j154        j156        j158

i14                                                            1
i35                        1
i48                                                                        1
i83                                                                                                1
i98                                                1
i117                                                                                                           1
i158                                                                                   1
i174                                   1
i194           1

   +        j161        j162        j163        j164        j166        j170        j172        j174        j175

i5                                                                         1
i32            1
i42                                                                                                1
i61                                                                                    1
i69                                                                        1
i103                                                           1
i143                                   1
i145                                                                                                           1
i158                                                                                               1
i192                       1
i198                                               1

   +        j176        j178        j179        j180        j182        j183        j184        j188        j191

i6                                                                         1
i13                                                                                                            1
i23                                                1
i47                                                                                                            1
i61                                                                                    1
i81                        1
i93            1
i103                                                                                                           1
i125                                   1
i182                                                           1
i193                                                                                               1

   +        j192        j193        j196

i73                                    1
i120           1
i138           1
i167                       1



We have really debunked that this is an intractable problem for a good MIP solver.

Note that the age of the universe is 13.7 billion years or 4.3e17 seconds. So we have achieved a speed-up of about 1e17. This is a record for me.

It is noted that this model does not really answer the question: what are the best locations for my stores? It only answers: how many stores do we need? In particular, it does not try to place stores as close as possible to customers. A subsequent optimization model can find good (or optimal) locations for the stores. Reference [2] discusses a similar issue for the continuous case, where facilities can be placed in any \((x,y)\) location instead of on a grid \((i,j)\).

Conclusion


The question: what is the minimum number of stores needed to service all customers can be answered by formulating and solving a set covering problem. This will give you locations for the stores. However, this model only looks at this objective, so no further attempts are undertaken to minimize the customer-store distances.

We can solve this problem for this particular problem quite efficiently using a standard MIP solver. Quite often I see that these solution approaches are dismissed, due to a misunderstanding of complexity-theoretic arguments. Complexity theory says very little about the performance of a particular MIP model using a particular solver on a given data set. It is better not to reject a MIP approach even before trying.

References


  1. Minimize number of shops while reaching all customers, https://stackoverflow.com/questions/62224042/minimize-number-of-shops-while-reaching-all-customers
  2. Solving a facility location problem as an MIQP, https://yetanothermathprogrammingconsultant.blogspot.com/2018/01/solving-facility-location-problem-as.html (this is essentially a continuous version of the current problem)

Thursday, March 12, 2020

CellBlock or Shikaku puzzle

This is a simple puzzle: given a grid with some cells containing a positive integer, form rectangles such that

  • Each rectangle contains exactly one cell with a number
  • The size of the rectangle is equal to this number
  • Every cell in the grid is covered by exactly one rectangle 
Example:

Shikaku Problem and Solution

Approach


One possible approach is a two-step algorithm:

  1. Generate all feasible rectangles,
  2. Use a simple MIP model to select a subset of these rectangles such that each cell is covered exactly once.

For the problem above, we just have 31 feasible rectangles. Feasible means here: (a) rectangle contains exactly one cell with a number, (b) the size of the rectangle is equal to that number.

All feasible rectangles


Designing an algorithm that generates all feasible rectangles is easier said than done, especially if you want to do this efficiently for larger instances. This is mainly a programming exercise, rather than a modeling problem. 

The second part is surprisingly simple. Assume all \(k\in K\) rectangles are stored in a single large boolean data structure: \[\mathit{allRectangles}_{k,i,j} = \begin{cases} 1 & \text{if cell \((i,j)\) is part of rectangle \(k\)}\\  0 & \text{otherwise}\end{cases}\] then the optimization model is a simple set covering model:

Set Covering Problem
\[\begin{align}\min\>& 0\\ &\sum_k \color{darkblue}{\mathit{allRectangles}}_{k,i,j} \color{darkred}x_{k} = 1&& \forall i,j\\ &\color{darkred}x_{k} \in \{0,1\} \end{align}\]

In this case, we don't have an objective. We simply state that every cell \((i,j)\) must be covered by exactly one rectangle.

This is a tiny model and solves very easily. Typically it can be solved completely in the presolve phase:


Tried aggregator 1 time.
MIP Presolve eliminated 37 rows and 32 columns.
MIP Presolve modified 3 coefficients.
All rows and columns eliminated.

The solution can look like:

----     89 VARIABLE x.L  selection of rectangles

k1  1.000,    k3  1.000,    k5  1.000,    k8  1.000,    k9  1.000,    k12 1.000,    k20 1.000,    k22 1.000
k26 1.000,    k29 1.000,    k31 1.000


----     97 PARAMETER result  solution: selected rectangles

            c1          c2          c3          c4          c5          c6

r1           1           3           5           5           8           9
r2           1           3           5           5           8           9
r3          12           3           5           5           8           9
r4          12           3          20          20          22          22
r5          12           3          26          26          26          26
r6          29          29          29          29          31          31

The model selects 11 rectangles. We could have predicted this: there are 11 numbered cells in the input data. The parameter result displays the selected rectangle numbers on the grid. You should be able to match this up with the picture of the feasible rectangles shown earlier.

GAMS Model


$ontext

  
Cell Block or Shikaku puzzle
  
Set Covering model for (small) square grids
  
erwin@amsterdamoptimization.com

$offtext

sets
  i
'rows' /r1*r6/
  j
'columns' /c1*c6/
;

scalars
   m
'rows'
   n
'columns'
;
m =
card(i);
n =
card(j);

alias(i,p,h);
alias(j,q,w);

table grid(i,j)
    
c1 c2 c3 c4 c5 c6
r1    2              3
r2             6  3
r3       5
r4    3     2     2
r5                4
r6          4        2
;

scalar nn 'sum of numbers';
nn =
sum((i,j),grid(i,j));
display nn;

abort$(nn <> card(i)*card(j)) "error in grid specification";


*------------------------------------------
* generate all possible rectangles
* requirements:
*  1. block contains exactly one number from grid
*  2. size of block is equal to this number
*------------------------------------------

set
  k
'rectangle number' /k1*k10000/
  allRectangles(k,i,j)
  r(i,j)
'current rectangle'
;
singleton set
  kk(k)
'current index' /k1/
;
loop((p,q),
  
loop((w,h)$(ord(p)+ord(h)-1<=m and ord(q)+ord(w)-1<=n),
* candidate rectangle
      r(i,j) = (
ord(i)>=ord(p) and ord(i)<=ord(p)+ord(h)-1 and
               
ord(j)>=ord(q) and ord(j)<=ord(q)+ord(w)-1);
* check: only one grid number and size = gridnumber
     
if(sum(r$grid(r),1)=1 and card(r)=sum(r,grid(r)),
         allRectangles(kk,r) =
yes;
         kk(k)=kk(k-1);
      );
   );
);

*------------------------------------------
* model: select rectangles such that each
* cell is covered exactly once
*------------------------------------------

binary variable x(k) 'selection of rectangles';
variable z 'dummy objective';

equations
   dummy        
'dummy objective'
   covering(i,j)
'each cell is covered by exactly one rectangle'
;

covering(i,j).. 
sum(allRectangles(k,i,j),x(k)) =e= 1;
dummy.. z =e= 0;

model cover /all/;
* just need feasibility
solve cover minimizing z using mip;

display x.l;

* reporting
parameter result(i,j) 'solution: selected rectangles';
loop(k$(x.l(k)>0.5),
   result(i,j)$allRectangles(k,i,j) =
ord(k);
);
option result:0;
display result;


Notes:
  • The data check makes sure the sum of the numbers is equal to the size (area) of the grid.
  • The enumeration of the feasible rectangles is using a simplistic and crude algorithm. It may need to be refined for larger instances. 
  • allRectangles(k,i,j) is implemented as a set.
  • As we are looking for a feasible solution only, there is no need to specify a gap with the optcr option.

A larger problem


A larger problem is taken from [1].

17 x 15 problem

Some statistics:
  • This problem has 17 rows and 15 columns. 
  • The number of feasible rectangles is 253.
  • The number of numbered cells is 66. This is also the number of rectangles we need to select.
  • If one wants to make a big impression, mention that the number of ways we can choose 66 out of 253 is: \[ \begin{align} {253 \choose 66} &= 66 281 207 211 442 081 043 167 193 447 008 690 678 844 358 496 354 619 442 275 550 \\ &\approx 6.6 \times 10^{61}\end{align}\]
The results look like:


----    105 VARIABLE x.L  selection of rectangles

k1   1.000,    k3   1.000,    k6   1.000,    k11  1.000,    k13  1.000,    k16  1.000,    k30  1.000,    k34  1.000
k36  1.000,    k37  1.000,    k41  1.000,    k44  1.000,    k45  1.000,    k51  1.000,    k54  1.000,    k63  1.000
k65  1.000,    k69  1.000,    k73  1.000,    k82  1.000,    k84  1.000,    k87  1.000,    k90  1.000,    k91  1.000
k94  1.000,    k103 1.000,    k107 1.000,    k110 1.000,    k112 1.000,    k125 1.000,    k128 1.000,    k130 1.000
k132 1.000,    k134 1.000,    k136 1.000,    k140 1.000,    k143 1.000,    k150 1.000,    k152 1.000,    k154 1.000
k159 1.000,    k162 1.000,    k171 1.000,    k174 1.000,    k190 1.000,    k192 1.000,    k195 1.000,    k199 1.000
k201 1.000,    k202 1.000,    k203 1.000,    k204 1.000,    k209 1.000,    k216 1.000,    k219 1.000,    k222 1.000
k226 1.000,    k228 1.000,    k235 1.000,    k238 1.000,    k244 1.000,    k245 1.000,    k247 1.000,    k249 1.000
k251 1.000,    k253 1.000


----    113 PARAMETER result  solution: selected rectangles

             c1          c2          c3          c4          c5          c6          c7          c8          c9

r1            1           3           3           6           6          11          11          11          11
r2            1           3           3           6           6          11          11          11          11
r3           30          30          30          30          34          34          36          37          37
r4           30          30          30          30          34          34          36          51          51
r5           30          30          30          30          34          34          63          65          65
r6           73          73          73          73          34          34          63          65          65
r7           87          87          90          90          91          91          94          65          65
r8          103         103         103         103         103         103          94         107         107
r9          103         103         103         103         103         103          94         125         125
r10         134         136         136         140         140         143         143         125         125
r11         134         154         154         159         159         159         159         159         159
r12         171         154         154         174         174         174         174         174         174
r13         171         190         190         192         195         195         195         199         199
r14         204         190         190         192         195         195         195         209         209
r15         204         190         190         216         216         219         219         222         222
r16         204         190         190         216         216         235         235         235         238
r17         204         244         244         245         245         245         247         247         249

  +         c10         c11         c12         c13         c14         c15

r1           11          11          13          13          16          16
r2           11          11          13          13          16          16
r3           41          41          13          13          44          45
r4           54          54          13          13          44          45
r5           54          54          69          69          44          45
r6           82          84          84          84          44          45
r7           82          84          84          84          44          45
r8          107         110         110         112          44          45
r9          128         128         130         112         132         132
r10         128         128         130         150         150         152
r11         162         162         130         150         150         152
r12         162         162         130         150         150         152
r13         199         201         201         202         203         203
r14         209         209         209         202         203         203
r15         222         222         222         226         226         228
r16         238         238         238         226         226         228
r17         249         249         251         251         253         253


Again the presolver could eliminate all rows and columns.

Questions


I still have some lingering questions about this problem.

  1. How would a better algorithm that enumerates all feasible rectangles look like? I guess it would start with the cells containing numbers.
  2. Can we use the solution pool to automate this step?
  3. Is the solution unique?
  4. Are there cases where the presolver is not solving the complete problem?
  5. Nerdy: when disabling the presolver. it looks like we can solve the problem as an LP. I.e. the binary variables are integer-valued automatically. Is this true?

Q2: Solution Pool to generate feasible rectangles


In the comments, a MIP model is proposed that will deliver feasible rectangles. It is worthwhile to reproduce it here. A high-level version of this model can look like:

High-level Feasible Rectangle Model
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ &\sum_{i,j|\color{darkblue}{\mathit{Grid}}_{i,j}} \color{darkred}{\mathit{useCell}}_{i,j} = 1 && && \text{One numbered cell in rectangle}\\ & \sum_{i,j} \color{darkblue}{\mathit{Grid}}_{i,j} \color{darkred}{\mathit{useCell}}_{i,j} = \sum_{i,j} \color{darkred}{\mathit{useCell}}_{i,j} &&&&\text{Number equal to area of rectangle} \\ & \color{darkred}{\mathit{useCell}}_{i,j} = \color{darkred}{\mathit{useRow}}_i \cdot \color{darkred}{\mathit{useColumn}}_j &&\forall i,j && \text{Interior of rectangle} \\ & \color{darkred}{\mathit{useRow}}_i \ge \color{darkred}{\mathit{useRow}}_{i1} \cdot \color{darkred}{\mathit{useRow}}_{i2} && \forall i1 \lt i \lt i2 && \text{Contiguous rows} \\ & \color{darkred}{\mathit{useColumn}}_j  \ge \color{darkred}{\mathit{useColumn}}_{j1} \cdot \color{darkred}{\mathit{useColumn}}_{j2} && \forall j1 \lt j \lt j2 && \text{Contiguous columns}\\ & \color{darkred}{\mathit{useCell}}_{i,j},\color{darkred}{\mathit{useRow}}_i,\color{darkred}{\mathit{useColumn}}_j \in \{0,1\} \end{align}\]


This model is nonlinear: there are quadratic terms. We are lucky, however. The quadratic terms involve only binary variables. This means that we can use a standard linearization:\[z=x\cdot y \iff \begin{cases} z \le x \\ z \le y \\ z \ge x+y-1 \end{cases}\] We assume here that \(x,y,x\) are binary variables. When we apply this, the linearized model looks like:

Linearized Feasible Rectangle Model
\[\begin{align}\min\>& 0 \\ &\sum_{i,j|\color{darkblue}{\mathit{Grid}}_{i,j}} \color{darkred}{\mathit{useCell}}_{i,j} = 1 \\ & \sum_{i,j} \color{darkblue}{\mathit{Grid}}_{i,j} \color{darkred}{\mathit{useCell}}_{i,j} = \sum_{i,j} \color{darkred}{\mathit{useCell}}_{i,j} \\ & \color{darkred}{\mathit{useCell}}_{i,j} \le \color{darkred}{\mathit{useRow}}_i &&\forall i,j \\ & \color{darkred}{\mathit{useCell}}_{i,j} \le \color{darkred}{\mathit{useColumn}}_j &&\forall i,j \\ & \color{darkred}{\mathit{useCell}}_{i,j} \ge \color{darkred}{\mathit{useRow}}_i + \color{darkred}{\mathit{useColumn}}_j - 1 &&\forall i,j \\ & \color{darkred}{\mathit{useRow}}_i \ge \color{darkred}{\mathit{useRow}}_{i1} + \color{darkred}{\mathit{useRow}}_{i2} - 1 && \forall i1 \lt i \lt i2  \\ & \color{darkred}{\mathit{useColumn}}_j \ge \color{darkred}{\mathit{useColumn}}_{j1} + \color{darkred}{\mathit{useColumn}}_{j2} - 1 && \forall j1 \lt j \lt j2 \\ & \color{darkred}{\mathit{useCell}}_{i,j},\color{darkred}{\mathit{useRow}}_i,\color{darkred}{\mathit{useColumn}}_j \in \{0,1\} \end{align}\]

After using a solution pool algorithm (this is part of solvers like Cplex and Gurobi) on this model, we get an enumeration of all feasible rectangles:

Feasible rectangles
For data set 1, we get 31 rectangles and for the second data set, we get 253 rectangles.

Q3. Uniqueness


To establish whether solutions are unique, we can use a standard approach:


  1. Solve the optimization problem. Record the optimal values as \(x_k^*\).
  2. Add the cut: \[  \sum_k x_k^* x_k - \sum_k (1-x_k^*)x_k \le \sum_k x_k^* - 1\] to the problem. This cut forbids the solution \(x_k^*\) but allows any other solution.
  3. Resolve. If infeasible then the solution \(x_k^*\) is unique.


The cut can be specialized. Let \(M\) be the number of numbered cells. We know that \[\sum_k x_k = M\] With this our cut can look like: \[\sum_k x_k^* x_k \le M-1\]

An alternative, somewhat heavy-handed approach would be to use the solution pool to prove uniqueness.

When we try this, we see that both data sets yield unique solutions.

References