Thursday, April 2, 2020

On the importance of well-written error messages

This is funny:
"I also tried to get Fortran programs running on Imperial's IBM 7010, since I sort of knew Fortran, certainly better than I knew Cobol, and Fortran would have been better suited for data analysis. It was only after weeks of fighting JCL, IBM's Job Control Language, that I deduced that there was no Fortran compiler on the 7010, but the JCL error messages were so inscrutable that no one else had figured that out either." [1]

Indeed, even today error messages are often poorly formulated and leaving users in the dark.

IBM 7010


References

  1. Brian Kernighan, UNIX, A History and a Memoir, 2020. 

Wednesday, April 1, 2020

Mondriaan Tilings


Introduction



Piet Mondriaan (later Mondrian) (1872-1944) was a Dutch painter, famous for his abstract works [1]. 



A tiling problem is related to this [2]:



The problem can be stated as follows:
  • Given is a square \((n \times n)\) area,
  • Partition this square into multiple, integer-sided, non-congruent rectangles,
  • Minimize the difference in size between the largest and the smallest rectangle.
Non-congruent means: only one \((p \times q)\) or \((q\times p)\) rectangle can be used.

It is interesting to make some predictions.

  • We try not to use a combination of very small and very large rectangles. That is obvious.
  • Only using small rectangles will always lead to a larger difference, as we need many small rectangles to cover the complete square. This will lead automatically to a larger range of sizes.
  • Very large rectangles are probably no good either: high probability we need some very small rectangles to fill up any remaining holes.
  • So I guess, we will use medium-sized rectangles of similar size. 


Example


Consider a small problem with \(n=5\). The possible sides for each rectangle are \(1,\dots,n\). From this, we can form a set of possible rectangles (I automated this step):



----     65 SET s  sizes of rectangles (height or width)

s1,    s2,    s3,    s4,    s5


----     65 SET k  rectangle id

k1 ,    k2 ,    k3 ,    k4 ,    k5 ,    k6 ,    k7 ,    k8 ,    k9 ,    k10,    k11,    k12,    k13,    k14


----     65 SET rec  rectangles

                s1          s2          s3          s4

k1 .s1         YES
k2 .s2         YES
k3 .s2                     YES
k4 .s3         YES
k5 .s3                     YES
k6 .s3                                 YES
k7 .s4         YES
k8 .s4                     YES
k9 .s4                                 YES
k10.s4                                             YES
k11.s5         YES
k12.s5                     YES
k13.s5                                 YES
k14.s5                                             YES

We see there are 14 different rectangles we can use. Note that the rectangle with size \((5 \times 5)\) is excluded. It would completely fill the area, while we need to place multiple rectangles. Also note: we only have \((p \times q)\) but not \((q \times p)\).  We cannot use both such rectangles as they are congruent. We can rotate a rectangle when placing it.

Optimal solution


An optimal solution for this problem is:


These are rectangles \(k5, k6, k12\) with \(k5\) being rotated. Obviously, this solution is not unique.

Optimization model


To model this, I use a binary variable: \[ x_{k,r,i,j} = \begin{cases} 1 & \text{if rectangle $k$ (rotated if $r=\mathrm{y}$) is placed at location $(i,j)$} \\ 0 & \text{otherwise}\end{cases}\] Here \(r = \{\mathrm{n},\mathrm{y}\}\) indicates if a rectangle is rotated. Note that not all rectangles can be rotated (it does not make sense to rotate a square, or there may not be space left to rotate for a given position).

For convenience, I also introduce \[y_{k} = \begin{cases} 1 & \text{if rectangle $k$ is used} \\  0 & \text{otherwise}\end{cases}\] This is essentially an aggration of \(x\).

In addition, I use a few more complex data structures: \[\mathit{ok}_{k,r,i,j} = \begin{cases} \mathit{Yes} & \text{if we can place rectangle $k$ at location $(i,j)$} \\ \mathit{No} & \text{otherwise}\end{cases}\] Again, \(r\) indicates if rectangle \(k\) is rotated. I implemented this as a GAMS set. You can consider it as a Boolean data structure. Furthermore: \[\mathit{cover}_{k,r,i,j,i',j'} = \begin{cases} \mathit{Yes} & \text{if cell $(i',j')$ is covered by placing rectangle $k$ at location $(i,j)$} \\ \mathit{No} & \text{otherwise}\end{cases}\]

These sets are large. The set ok looks like:


----    103 SET ok  usable position

INDEX 1 = k1

              j1          j2          j3          j4          j5

n.i1         YES         YES         YES         YES         YES
n.i2         YES         YES         YES         YES         YES
n.i3         YES         YES         YES         YES         YES
n.i4         YES         YES         YES         YES         YES
n.i5         YES         YES         YES         YES         YES

INDEX 1 = k2

              j1          j2          j3          j4          j5

n.i1         YES         YES         YES         YES         YES
n.i2         YES         YES         YES         YES         YES
n.i3         YES         YES         YES         YES         YES
n.i4         YES         YES         YES         YES         YES
y.i1         YES         YES         YES         YES
y.i2         YES         YES         YES         YES
y.i3         YES         YES         YES         YES
y.i4         YES         YES         YES         YES
y.i5         YES         YES         YES         YES

...

INDEX 1 = k13

              j1          j2          j3

n.i1         YES         YES         YES
y.i1         YES
y.i2         YES
y.i3         YES

INDEX 1 = k14

              j1          j2

n.i1         YES         YES
y.i1         YES
y.i2         YES


The set cover is even larger:


----    103 SET cover  coverage of cells

INDEX 1 = k1  INDEX 2 = n  INDEX 3 = i1  INDEX 4 = j1

            j1

i1         YES

INDEX 1 = k1  INDEX 2 = n  INDEX 3 = i1  INDEX 4 = j2

            j2

i1         YES

...

INDEX 1 = k14  INDEX 2 = y  INDEX 3 = i1  INDEX 4 = j1

            j1          j2          j3          j4          j5

i1         YES         YES         YES         YES         YES
i2         YES         YES         YES         YES         YES
i3         YES         YES         YES         YES         YES
i4         YES         YES         YES         YES         YES

INDEX 1 = k14  INDEX 2 = y  INDEX 3 = i2  INDEX 4 = j1

            j1          j2          j3          j4          j5

i2         YES         YES         YES         YES         YES
i3         YES         YES         YES         YES         YES
i4         YES         YES         YES         YES         YES
i5         YES         YES         YES         YES         YES


The reason to develop these sets is to make the equations simpler. In addition, we can develop, test and debug the calculation of these sets in isolation of the test of the model and before the equations are written.

With this, we can formulate:

MIP Model
\[\begin{align}\min\> & \color{darkred}u - \color{darkred}\ell\\ & \sum_{k,r,i,j|\color{darkblue}{\mathit{cover}}(k,r,i,j,i',j')}\color{darkred}x_{k,r,i',j'}=1 && \forall i',j'&&\text{cover cell $(i',j')$}\\ &\color{darkred}y_k = \sum_{r,i,j|\color{darkblue}{\mathit{ok}}(k,r,i,j)}\color{darkred}x_{k,r,i,j}&& \forall k && \text{select max one type of rectangle} \\ & \color{darkred}\ell \le \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k + \color{darkblue}M (1-\color{darkred}y_k) && \forall k && \text{bound smallest area}\\ & \color{darkred}u \ge \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k  && \forall k && \text{bound largest area}\\ & \sum_k \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k = \color{darkblue} n^2 && && \text{extra constraint}\\ & \color{darkred}x_{k,r,i,j} \in \{0,1\} \\&  \color{darkred}y_k \in \{0,1\} \\ & \color{darkred}u, \color{darkred}\ell \ge 0 \\& \color{darkblue}M = \max_k \color{darkblue}{\mathit{area}}_k\end{align}\]


When solving this model, we see:


----    146 VARIABLE x.L  place tile

                  j1          j4

k5 .y.i4           1
k6 .n.i1           1
k12.n.i1                       1


----    149 VARIABLE y.L  tile is used

k5  1,    k6  1,    k12 1


----    157 PARAMETER v  tile numbers 

            j1          j2          j3          j4          j5

i1           6           6           6          12          12
i2           6           6           6          12          12
i3           6           6           6          12          12
i4           5           5           5          12          12
i5           5           5           5          12          12


----    159 VARIABLE smallest.L            =        6.000  smallest used tile
            VARIABLE largest.L             =       10.000  largest used tile
            VARIABLE z.L                   =        4.000  objective


It is noted that these types of formulations lead to very large models [5].

One remaining question is whether we can use some symmetry-breaking constraint to speed up solution times.

A larger problem


It looks like \(n=19\) is interesting [3]. This is about the limit that we can solve comfortably (a few hours to global optimality) with a MIP solver. This model had about 36,000 binary variables.

An optimal 19x19 problem

The colors are not part of the model: they were added by hand to make the picture more Mondrianic.

A list of optimal objective values can be found in [4].

References


Thursday, March 26, 2020

Special Sudoku (2)

In [1], a variant of Sudoku was discussed, where a side constraint is:

Orthogonally neigboring cells differ in value by at most 5

This max 5 difference is a convex constraint, so it can be expressed quite easily: \[-5 \le v_{i,j} - v_{i',j'} \le 5\] for neighboring cells \((i,j)\) and \((i',j')\).

A more difficult restriction is: neighboring cells must differ by at least 2. Here is such a problem, again designed by Aad van de Wetering:


Problem data

The rules are:
  • Standard Sudoku-rules apply. cells in each row, column and sub-block have unique values from 1 through 9.
  • In addition, neighboring (horizontal and vertical) cells have values that differ by 2 or more.
  • Also, values 1 and 9 are considered to have a difference of 1. So neighboring cells, cannot have values 1 and 9.

I will try to solve this problem in two different ways. The first is an extension of what we did in [1]. The second one is a bit different.

Approach I


The constraint: neighbors must be different by at least two, can be formalized as \[\left| v_{i,j} - v_{i',j'} \right| \ge 2\] for neighboring cells \((i,j)\) and \((i',j')\). The non-convexity of this constraint  \[\begin{align} & v_{i,j} - v_{i',j'} \ge 2 \\ & \mathbf{or} \\&v_{i,'j'} - v_{i,j} \ge 2 \end{align} \] requires extra binary variables and big-M constraints: \[\begin{align} & v_{i,j} - v_{i',j'} \ge 2 - M\cdot \delta_{i,j,i',j'} \\ & v_{i,'j'} - v_{i,j} \ge 2 - M\cdot (1-\delta_{i,j,i',j'})\\ &  \delta_{i,j,i',j'} \in \{0,1\}\end{align}\] So, how large should \(M\) be? The maximum difference between two different values is 8. So we need \(2-M\le -8\) or \(M=10\).

We ignored that values 1 and 9 also have a difference of 1. As all other combinations are allowed, we can handle this by \[\left|v_{i,j}-v_{i',j'}\right| \le 7\] This is a constraint we already know how to handle: \[-7 \le v_{i,j}-v_{i',j'} \le 7\] Interestingly, this addition restriction also can help us in tightening our big-M value. We can now use \(M = 9\). Note that with this smaller big-M we still need this constraint: the big-M constraint only covers one of the needed two sides of this constraint.

The complete model can look like:


Mixed-Integer Programming Model I
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ &\color{darkred}v_{i,j} =\sum_k \color{darkblue}k\cdot\color{darkred}x_{i,j,k} && \forall i,j &&\text{Value of cell}\\ & \color{darkred}v_{i,j} - \color{darkred}v_{i',j'} \ge 2 - \color{darkblue}M\cdot \color{darkred}\delta_{i,j,i',j'}&&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} &&\text{Min difference of 2} \\ & \color{darkred}v_{i,'j'} - \color{darkred}v_{i,j} \ge 2 - \color{darkblue}M\cdot (1-\color{darkred}\delta_{i,j,i',j'}) &&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'}\\ & -7 \le \color{darkred}v_{i,j} - \color{darkred}v_{i',j'} \le 7 &&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} && \text{Neighbors differ by max 7} \\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\&\color{darkred}x_{i,j,k} \in \{0,1\} \\ &\color{darkred}v_{i,j} \in \{1,\dots,9\} \\ & \color{darkred}\delta_{i,j,i',j'} \in \{0,1\} \\ & \color{darkblue}M = 9\end{align}\]


Notes:

  • This model solves quickly: the presolve can remove all rows and columns.
  • The solution is unique (see [1]).
  • Using the above data set we actually can drop the "max 7" constraint. I think that is just a property of this data set. In general, we should include this constraint.


Approach II


A simpler approach is to consider the pairs of values that cannot be neighbors: \[(1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,1)\] We also need to add their symmetrical counterparts \[(2,1),(3,2),\dots\] Detail: this is due to the way we organized the neighbors. The neighbor set did not have both: i.e. cells \(\text{$(1,2)$-$(1,3)$}\) are neighbors, so this is stored. But, I did not also store the combination \(\text{$(1,3)$-$(1,2)$}\).

After automatically generating the combinations \((k,k')\) that we need to prevent, the forbid data can look like:



----     55 SET forbid  combinations of values we want to exclude

            v1          v2          v3          v4          v5          v6          v7          v8          v9

v1                     YES                                                                                 YES
v2         YES                     YES
v3                     YES                     YES
v4                                 YES                     YES
v5                                             YES                     YES
v6                                                         YES                     YES
v7                                                                     YES                     YES
v8                                                                                 YES                     YES
v9         YES                                                                                 YES


This table contains elements \(\mathit{forbid}_{k,k'}\) we want to forbid. Note that we can exclude the diagonal from this: neighboring cells already cannot have the same value. For this set of value-pairs, we want to add the constraint:\[x_{i,j,k} + x_{i',j',k'} \le 1 \] for all neighboring cells \((i,j)\) and \((i',j')\) and for all forbidden value combinations \((k,k')\). We have 18 of these combinations and 144 neighbors to check, yielding 2592 of these constraints.

The full model can look like:

Mixed-Integer Programming Model II
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ & \color{darkred}x_{i,j,k} + \color{darkred}x_{i',j',k'} \le 1 && \begin{aligned}& \forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} \\ & \forall k,k'|\color{darkblue}{\mathit{forbid}}_{k,k'} \end{aligned}&&\text{Forbid combinations} \\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\&\color{darkred}x_{i,j,k} \in \{0,1\} \end{align}\]


When we run this model, we will see the solution:


----     84 PARAMETER v  values (postprocessing)

            c1          c2          c3          c4          c5          c6          c7          c8          c9

r1           9           2           6           4           1           7           3           5           8
r2           3           5           8           2           6           9           7           1           4
r3           1           7           4           8           3           5           2           6           9
r4           7           4           1           6           9           2           5           8           3
r5           5           8           3           1           7           4           9           2           6
r6           2           6           9           3           5           8           4           7           1
r7           6           9           2           5           8           3           1           4           7
r8           8           3           5           7           4           1           6           9           2
r9           4           1           7           9           2           6           8           3           5

This solution is unique. The presolver completely solves the model, so no actual iterations are performed.

References



  1. Special Sudoku, http://yetanothermathprogrammingconsultant.blogspot.com/2020/03/special-sudoku.html

Sunday, March 22, 2020

Special Sudoku

Here is a difficult Sudoku variant designed by Aad van de Wetering. Aad is a former colleague of mine.

Problem
The standard Sudoku rules apply:

  • Each cell has an integer value between 1 and 9 
  • In each row, column and sub-block cells must be unique
  • Some cells have a predefined value

In addition, we have the following restriction:

  • The maximum difference between the value of two neighboring cells (orthogonal, that is horizontal or vertical) is 5

To solve this you can follow this video [1]:



Of course, we want to try to solve this as a Mixed-Integer Programming model.

Review: a linear integer programming model for standard Sudoku


When we want to solve Sudokus, the easiest approach is to define the following binary decision variables [2]: \[x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has value $k$} \\ 0 & \text{otherwise}\end{cases}\] Here \(k \in \{1,\dots,9\}\). We have 27 areas we need to check for unique values: rows, columns and sub-blocks. We can organize this as a set:  \[u_{a,i,j}\>\text{exists if and only if area $a$ contains cell $(i,j)$}\] This is data. We also have a set of given cells, which we can fix. The resultimg MIP model can look like [2]:

Mixed-Integer Programming Model for standard Sudoku
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\ &\color{darkred}x_{i,j,k} \in \{0,1\} \end{align}\]

During post-processing, we can calculate the completed grid using the optimal values of \(x\): \[v_{i,j} = \sum_k k \cdot x^*_{i,j,k}\]

Linear integer programming model for Sudoku variant


We need to introduce the concept of "neighbors". We can define the set: \[nb_{i,j,i',j'} \> \text{exists iff cells $(i,j)$ and $(i',j')$ are neighbors}\] I excluded symmetric cases. To be precise: in GAMS, I populated this set with:

alias (i,ii),(j,jj);
set
  nb(i,j,ii,jj)
'neighbors'
;

nb(i,j,i-1,j) =
yes;
nb(i,j,i,j-1) =
yes;


We can introduce the post-processing step directly into the MIP and then enforce our difference constraint on the variables \(v_{i,j}\). This can look like \[\begin{align} & v_{i,j} = \sum_k k \cdot x_{i,j,k} && \text{cell values as constraint} \\ & -5 \le v_{i,j} - v_{i',j'} \le 5 && \text{for neigboring cells $(i,j)$ and $(i',j')$}\end{align}\]
With this, we can formulate:


Mixed-Integer Programming Model for Sudoku Variant
\[\begin{align}\min\>& 0 && && \text{Dummy objective}\\ & \sum_k \color{darkred}x_{i,j,k} = 1 &&\forall i,j && \text{One $k$ in each cell}\\ & \sum_{i,j|\color{darkblue}u_{a,i,j}} \color{darkred}x_{i,j,k} = 1 && \forall a,k && \text{Unique values in each area}\\ &\color{darkred}v_{i,j} =\sum_k \color{darkblue}k\cdot\color{darkred}x_{i,j,k} && \forall i,j &&\text{Value of cell}\\ & -5 \le \color{darkred}v_{i,j} - \color{darkred}v_{i',j'} \le 5 &&\forall i,j,i',j'|\color{darkblue}{nb}_{i,j,i',j'} && \text{Neighbors differ by max 5} \\ & \color{darkred}x_{i,j,k} = 1 && \text{where } k=\color{darkblue}{\mathit{Given}}_{i,j} &&\text{Fix given values}\\&\color{darkred}x_{i,j,k} \in \{0,1\} \\ &\color{darkred}v_{i,j} \in \{1,\dots,9\} \end{align}\]

This model solves very easily. The presolver can solve this model completely, so we don't even have to start the branch-and-bound phase. As a result, we can make a much shorter movie:




The solution of this model looks like:


----     80 VARIABLE v.L  value of each cell

            c1          c2          c3          c4          c5          c6          c7          c8          c9

r1           2           4           9           6           7           5           8           3           1
r2           3           8           5           1           4           9           6           7           2
r3           7           6           1           2           3           8           9           5           4
r4           8           9           6           3           1           4           5           2           7
r5           5           7           3           8           6           2           1           4           9
r6           1           2           4           5           9           7           3           8           6
r7           6           1           2           4           8           3           7           9           5
r8           4           3           7           9           5           1           2           6           8
r9           9           5           8           7           2           6           4           1           3


Unique solution


A proper Sudoku puzzle has one unique solution. We can verify this by adding a cut to the problem after solving it. The cut should forbid the current solution (and allow all other solutions). In this case, the cut can look like \[\sum_{i,j,k} x^*_{i,j,k} x_{i,j,k}\le 9^2-1\] where \(x^*\) is the previously found solution. Indeed, this puzzle has a single, unique solution: after adding the cut, the problem becomes (integer) infeasible.

References


  1. A Sudoku With Just 11 Digits = Dutch Genius!, https://www.youtube.com/watch?v=ZU5fSDHJq8k
  2. MIP Modeling: from Sudoku to KenKen via Logarithms, https://yetanothermathprogrammingconsultant.blogspot.com/2016/10/mip-modeling-from-sudoku-to-kenken.html

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