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:
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).
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:
 |
| y matrix |
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. | |
---- 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
| Data set 1 | Data set 2 | |
|---|---|---|
| Data | ||
| Tiles | 21 | 45 |
| Area | 196 | 824 |
| Solution | ||
| Width | 14 | 28 |
| Area | 196 | 784 |
| Tiles used | 21 | 40 |
| Model A | ||
| Variables | 919 | 4,125 |
| Binary variables | 875 | 4,033 |
| Constraints | 1,126 | 5,116 |
| Solution time | 39 | NA |
| Nodes | 73,975 | NA |
| Model B | ||
| Variables | 4,378 | 36,196 |
| Binary variables | 4,151 | 35,353 |
| Constraints | 950 | 3,473 |
| Solution time | 1.3 | 256 |
| Nodes | 0 | 10 |
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:
- In the first example, we could use all the tiles of the inventory.
- 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.
References
- Max square side length of smaller squares, https://stackoverflow.com/questions/60526058/max-square-side-lenght-of-smaller-squares
No comments:
Post a Comment