## Monday, March 9, 2020

### Tiling

This problem is original from [1]. I made it a bit more general.

#### Problem statement

Given a collection of square tiles, what is the maximum square space we can fill with them?

#### Example

Let's consider our inventory:

----     29 PARAMETER input  inventory data

count       width        area

tile1           6           1           1
tile2           5           2           4
tile3           4           3           9
tile4           3           4          16
tile5           2           5          25
tile6           1           6          36
total          21                     196


So we have a total of 21 tiles of different sizes. Their total area is 196. This means that the best we can hope for is to fill square space of $$14 \times 14 = 196$$.

The optimization model shows we can fill this whole $$14 \times 14$$ area, using up all the tiles. No leftovers. As we will see later, this is somewhat of a special case. The optimal solution can look like:

----    110 PARAMETER tiles  solution

x           y       width        area

tile1-1                      12           1           1
tile1-2                      13           1           1
tile1-3          13                       1           1
tile1-4          13           1           1           1
tile1-5          13           2           1           1
tile1-6          13           3           1           1
tile2-1           3          10           2           4
tile2-2           7           6           2           4
tile2-3           1          12           2           4
tile2-4           3          12           2           4
tile2-5           7           8           2           4
tile3-1                                   3           9
tile3-2                       3           3           9
tile3-3                       6           3           9
tile3-4                       9           3           9
tile4-1           3           6           4          16
tile4-2           9                       4          16
tile4-3           5          10           4          16
tile5-1           9           4           5          25
tile5-2           9           9           5          25
tile6-1           3                       6          36
space                                    14         196


The solution presented in a numeric table is a bit difficult to grasp, so let's make a picture:

 An optimal solution

It is noted that this solution is not unique.

A proper tiling has no overlapping tiles and no uncovered spaces, i.e., no holes. The picture confirms this is the case here.

Below I will discuss two different formulations. Warning: they are a bit tricky.

#### Formulation A: no-overlap constraints

One way to model this is to consider no-overlap constraints between two squares. The idea is as follows.  Assume we have two squares $$i,j$$ with size $$s_k$$ and location $$(x_k,y_k)$$, then we can say: \begin{align} &x_j \ge x_i + s_i \>\mathbf{or} \\ & y_j \ge y_i+ s_i \>\mathbf{or} \\ & x_i \ge x_j + s_j \>\mathbf{or} \\ & y_i \ge y_j + s_j\end{align} A standard way to shoehorn this into a MIP model is to use big-M constraints: \begin{align} &x_j \ge x_i + s_i - M (1-\delta^{(1)}_{i,j}) \\ & y_j \ge y_i + s_i - M (1-\delta^{(2)}_{i,j}) \\ & x_i \ge x_j + s_j - M (1-\delta^{(3)}_{i,j}) \\ & y_i \ge y_j + s_j - M (1-\delta^{(4)}_{i,j})\\ & \sum_d \delta^{(d)}_{i,j} \ge 1\end{align} Alternatively, using the complement: \begin{align} &x_j \ge x_i + s_i - M \delta^{(1)}_{i,j} \\ & y_j \ge y_i + s_i - M\delta^{(2)}_{i,j} \\ & x_i \ge x_j + s_j - M \delta^{(3)}_{i,j} \\ & y_i \ge y_j + s_j - M \delta^{(4)}_{i,j}\\ & \sum_d \delta^{(d)}_{i,j} \le 3\end{align} The disadvantage of this formulation is that we need good, tight values for $$M$$.

Using this framework we can formulate the model:

Model A: No-overlap Constraints
Sets\begin{align}&i,j&&\text{tiles}\\ &a&&\text{possible widths for space to fill: 1,2,\dots} \\ & c =\{x,y\} && \text{coordinates}\\ & k = \{i\lt j,j\lt i\} && \text{no-overlap cases} \end{align}
Parameters\begin{align} & \color{darkblue}{\mathit{size}}_{i} && \text{size of tile i}\\ &\color{darkblue} M &&\text{big-M constants}\\ &\color{darkblue}{\mathit{sizeSpace}}_a = 1,2,\dots,\color{darkblue}{\mathit{maxSpace}} &&\text{possible sizes of total space we can fill}\end{align}
Variables\begin{align} &\color{darkred}t_i \in \{0,1\}&&\text{select tile i}\\ &\color{darkred}{p}_{i,c} \ge 0 &&\text{position of tile i, (x,y) coordinates} \\ &\color{darkred}W && \text{width of space we can fill} \\ &\color{darkred}A && \text{area of space to fill} \\ & \color{darkred}\delta_{i,j,k,c} \in \{0,1\} && \text{no-overlap cases} \\ & \color{darkred}\theta_a \in \{0,1\} && \text{select width of space to fill}\end{align}
Objective$\max \color{darkred}W$Maximize width of space we can fill
Constraints$\color{darkred}p_{i,c} + \color{darkblue}{\mathit{size}}_{i}\le\color{darkred}{\mathit{W}}\>\>\forall i,c$Each tile should be located inside the space we want to fill. Really only important for selected tiles, but it is easier to apply this constraint to all tiles.
\begin{align} & \color{darkred}p_{j,c} \ge \color{darkred}p_{i,c} + \color{darkblue}{\mathit{size}}_{i} - \color{darkblue}M (1-\color{darkred}\delta_{i,j,i\lt j,c}) - \color{darkblue}M(1-\color{darkred}t_i) - \color{darkblue}M(1-\color{darkred}t_j)&& \forall i\lt j,c \\ & \color{darkred}p_{j,c} + \color{darkblue}{\mathit{size}}_{j} \le \color{darkred}p_{i,c} + \color{darkblue}M (1-\color{darkred}\delta_{i,j,j\lt i,c}) + \color{darkblue}M(1-\color{darkred}t_i) + \color{darkblue}M(1-\color{darkred}t_j)&& \forall i\lt j,c \\ & \sum_{k,c} \color{darkred}\delta_{i,j,k,c} \ge 1 && \forall i\lt j\end{align}No-overlap constraints. Only compare when both tiles $$i$$ and $$j$$ are selected.
\begin{align} &\color{darkred}W = \sum_a \color{darkblue}{\mathit{sizeSpace}}_a \cdot\color{darkred}\theta_a\\ &\color{darkred}A = \sum_a \color{darkblue}{\mathit{sizeSpace}}_a^2 \cdot\color{darkred}\theta_a \\& \sum_a \color{darkred}\theta_a = 1 \end{align}Select width of space to fill. Both size and area are calculated. We use binary variables to prevent the model to become quadratic.
$\sum_i \color{darkblue}{\mathit{size}}_{i}\cdot \color{darkred}t_i = \color{darkred}A$Area of the selected tiles should be equal to the area of the outer space. This will ensure the whole space is covered with tiles.

In the model, we need both the area and the length of each side (called size in the model). The easiest would be to add $A = W^2$ However this would make the model non-linear. We use here that the outer area can only have integer-valued sizes.

A way to get good values for the big-$$M$$ constants is the following. From our inventory of tiles, we can calculate the total area we can cover. The maximum size of the outer box is the square root of this number, rounded down to the nearest integer: $\mathit{maxSpace} = \left\lfloor \sqrt{\sum_i \mathit{size}_i^2} \right\rfloor$ We can simply use $M = \mathit{maxSpace}$ We can refine this further by taking the size of the individual tiles into account (basically replacing $$M$$ by $$M_{i,j}$$.

I also added some ordering constraints: $t_i \ge t_{i+1} \>\>\text{if \mathit{size}_i = \mathit{size}_{i+1}}$ This will select tiles according to their ordinal number. Furthermore, $\sum_c p_{i,c} \ge \sum_c p_{i+1,c} \>\>\text{if \mathit{size}_i = \mathit{size}_{i+1}}$ This will order the locations, and hopefully breaks some symmetry in the model.

This model works fine for our example data. However, for larger instances, things become very slow.

#### Formulation B: grid approach

In this formulation, we work with binary variables: $x_{i,p,q}=\begin{cases} 1 & \text{if tile i is placed at location (p,q)} \\ 0 & \text{otherwise}\end{cases}$ Then for each location $$(p,q)$$ we check that the number tiles covering this location is equal to one within the area to fill and zero outside. To model this we develop a binary parameter $\mathit{cover}_{p1,q1,i,p,q}=\begin{cases} 1 & \text{if (p1,q1) is covered when tile i is placed at (p,q)}\\ 0 & \text{otherwise}\end{cases}$ Furthermore, we maintain a boolean matrix $$y_{p1,q1}$$ with $y_{p1,q1} = \begin{cases} 1 & \text{if (p1,q1) is within our current area to fill} \\ 0 & \text{otherwise}\end{cases}$ Notice that this area is shrinking during the solution process, so this matrix is dynamic (i.e. a variable).

 y matrix
Detail: $$y$$ has $$P+1$$ rows and columns. Initially, we just have one final row and column of zeroes. The size of $$P$$  is determined from the inventory: sum all areas of the tiles, take the square root and round down. The way we keep track of $$y$$ is to introduce a binary variable $$\delta_p\in \{0,1\}$$ with \begin{align} &\delta_p \le \delta_{p-1} \\ & \sum_p \delta_p = W\end{align} Now we just can do $$y_{p,q} = \delta_p \cdot \delta_q$$, for which a standard linearization is available.

With this, the covering constraint can look like: $\sum_{i,p,q} \mathit{cover}_{p1,q1,i,p,q} x_{i,p,q} = y_{p1,q1}\>\>\forall p1,q1$ This says: each cell inside the area to fill is covered by exactly one tile, and cells outside can not be covered. The complete model can look like:

Model B: Grid Approach
Sets\begin{align}&i&&\text{tiles}\\ &p,q && \text{gridpoints, p,q\in \{1,..,P\} }\\ & p1,q1 && \text{gridpoints plus extra row and column, p1,q1\in\{1,\dots,P+1\}} \end{align}
Parameters\begin{align} & \color{darkblue}{\mathit{size}}_{i} && \text{size of tile i}\\ & \color{darkblue}{\mathit{cover}}_{p1,q1,i,p,q}\in \{0,1\} && \text{1 if (p1,q1) is covered by tile i at (p,q)} \end{align}
Variables\begin{align} &\color{darkred}x_{i,p,q} \in \{0,1\}&&\text{place tile i at location (p,q)}\\ &\color{darkred}{y}_{p1,q1} \in \{0,1\} &&\text{1 inside area to be filled, 0 outside} \\ & \color{darkred}\delta_{p1} \in \{0,1\} && \text{1 if $$p1\le W$$}\\ &\color{darkred}W && \text{width of space we can fill} \end{align}The $$y$$ variables can be relaxed to $${y}_{p1,q1} \in [0,1]$$.
Objective$\max \color{darkred}W$Maximize width of space we can fill
Constraints$\sum_{p,q} \color{darkred}x_{i,p,q} \le 1 \>\>\forall i$Each tile can be placed once
$\sum_{i,p,q} \color{darkblue}{\mathit{cover}}_{p1,q1,i,p,q} \cdot \color{darkred} x_{i,p,q} = \color{darkred}y_{p1,q1}\>\>\forall p1,q1$Cover constraint. Inside the fill-area, each cell must be covered exactly once. Outside: zero coverage.
\begin{align} &\color{darkred} \delta_{p1} \le \color{darkred}\delta_{p1-1}\\ & \sum_{p1} \color{darkred}\delta_{p1}=\color{darkred} W\end{align}Force a pattern like 111000. Number of ones is $$W$$. The last element is a 0. This is used to form the $$y$$ matrix.
\begin{align} &\color{darkred}y_{p1,q1} \le \color{darkred}\delta_{p1}\\ & \color{darkred}y_{p1,q1} \le \color{darkred}\delta_{q1} \\ & \color{darkred}y_{p1,q1} \ge \color{darkred} \delta_{p1} +\color{darkred}\delta_{q1}-1\end{align}Maintain the $$y$$ matrix. The first W rows and columns are ones. The other entries are zeros. One extra row and column of zeros for safeguarding against covering outside the area. These constraints are a linearization of a binary multiplication.

This model can handle larger problems. I tried the inventory:

----     56 PARAMETER input  inventory data

count       width        area

tile1           9           1           1
tile2           8           2           4
tile3           7           3           9
tile4           6           4          16
tile5           5           5          25
tile6           4           6          36
tile7           3           7          49
tile8           2           8          64
tile9           1           9          81
total          45                     825


This is a larger data set with 45 tiles. The total area of the tiles is 825. The square root is $\sqrt{825} = 28.723$ So the maximum square we can hope to fill is $$28 \times 28$$ with an area of 784. The generated MIP model is very large:

MODEL STATISTICS

BLOCKS OF EQUATIONS           8     SINGLE EQUATIONS        3,473
BLOCKS OF VARIABLES           5     SINGLE VARIABLES       36,196  1 projected
NON ZERO ELEMENTS       614,345     DISCRETE VARIABLES     35,353


Although it is large, it solves rather easily. We just needed 10 nodes. It has the following optimal solution:

----    161 PARAMETER tiles  solution

x           y       width        area

item1-1          27          23           1           1
item1-2                      10           1           1
item1-3          10          13           1           1
item1-4                       9           1           1
item1-5          14           5           1           1
item1-6          14          18           1           1
item1-7          20           6           1           1
item1-8          27          24           1           1
item2-1          25          23           2           4
item2-2           1           9           2           4
item2-3                      14           2           4
item2-4                      24           2           4
item2-5                      16           2           4
item2-6                      26           2           4
item3-1                                   3           9
item3-2          11          16           3           9
item3-3                       6           3           9
item3-4                      11           3           9
item3-5                       3           3           9
item3-6          25          25           3           9
item3-7          11          13           3           9
item4-1           2          14           4          16
item4-2          10           5           4          16
item4-3          10           9           4          16
item4-4           2          24           4          16
item5-1          10                       5          25
item5-2          15          18           5          25
item5-3          15          23           5          25
item5-4           6          14           5          25
item5-5          20          23           5          25
item6-1          15                       6          36
item6-2          14           6           6          36
item6-3                      18           6          36
item6-4          14          12           6          36
item7-1           3                       7          49
item7-2           3           7           7          49
item7-3          21                       7          49
item8-1          20          15           8          64
item8-2          20           7           8          64
item9-1           6          19           9          81
space                                    28         784


Indeed we find the $$28 \times 28$$ tiling.

 An optimal solution

#### Summary of performance

These results are with Cplex, 8 threads, default settings.

Data set 1Data set 2
Data
Tiles2145
Area196824
Solution
Width1428
Area196784
Tiles used2140
Model A
Variables9194,125
Binary variables8754,033
Constraints1,1265,116
Solution time39NA
Nodes73,975NA
Model B
Variables4,37836,196
Binary variables4,15135,353
Constraints9503,473
Solution time1.3256
Nodes010

Obviously, approach B is performing much better, even though it generates much larger MIP models.

#### An example with smaller covered area

Above we discussed two examples:
1. In the first example, we could use all the tiles of the inventory.
2. In the second example, the square area we could cover was equal to the maximum possible area calculated from the total area of the tiles.

Here is a small example where we are much below this maximum estimate. The data looks like:

----     44 PARAMETER input  inventory data

count       width        area

tile1           1           1           1
tile2           1           2           4
tile3           1           3           9
tile4           1           4          16
tile5           1           5          25
tile6           1           6          36
tile7           1           7          49
tile8           1           8          64
tile9           1           9          81
total           9                     285


So we have one of each size from 1 to 9. Looking at the total area, the maximum size we can expect is:$\left\lfloor \sqrt{285} \right\rfloor = 16$ When we solve this, we see:

----    125 PARAMETER tiles  solution

width        area

tile9-1           9          81
space             9          81


I.e. just the $$9 \times 9$$ tile is used.

#### Conclusion

This problem is a bit more complex to model as a MIP.  The first "no-overlap" formulation did not perform very well. The second "grid cover" formulation did much better, and I could solve some larger problems in a few minutes.