sysid blog

Square Tiling

May 02, 2020 ☕️ 7 min read

5_5

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

The picture above consists of the following set of tiles:

length123
count432

8 out of the 9 given tiles can be used to form the resulting (5x5) square. One (3x3) tile cannot be used. With the given set this is the larges square area we can fill.

Yet Another Mathprogramming Consultant provides two different formulations of the problem. Here I will focus on the better performing ‘grid approach’. Also this approach has similarity with Mondriaan Puzzle so we can reuse and extend some of already learned techniques here.

Grid Approach

We split the potential area which we want to fill into a grid of element cells. This allows us to impose constraints per cell and formulate packing and coverage constraints elegantly.

Main Idea

We introduce a binary data structure to encode the information which grid cells are covered by a tile ii at position (p,q)(p,q):

coveri,p,q,p,q={1,  if (p,q) is covered when tile i is placed at (p,q)0,  else cover_{i,p,q,p',q'} = \begin{cases} 1, \ \text{ if $(p',q')$ is covered when tile i is placed at $(p,q)$}\\ 0, \ \text{ else }\\ \end{cases}

The main idea is now the introduction of a boolean ‘grid’ variable which marks the area to fill:

yp,q={1,  if (p,q) is within our current area to fill 0,  else y_{p',q'} = \begin{cases} 1, \ \text{ if $(p',q')$ is within our current area to fill }\\ 0, \ \text{ else }\\ \end{cases}

This is a variable, because the area is shrinking during the solution process.

main_idea

Here we have a tile at (6,3)(6,3) which conforms to the covering constraint (CC), i.e. lies within the viable squared area and is therefor a solution candidate.

To keep track of yy and make sure there are not holes we need another (ordered) binary variable:

δpδp1, ppδp=W\delta_p \le \delta_{p-1}, \ \forall p\\ \sum_p \delta_p = W \\

To form a ‘square constrain’ for our grid variable we formulated:

yp,q=δpδqy_{p,q} = \delta_p \delta_q \\

This non-linearity can be linearized by standard means (see below).

Upper Bound for solution

The maximal possible area which we can fill if we use all tiles from inventory is just the sum of all individual tile areas. So an upper bound for the side length of the enclosing square is:

P=floor(isizei2)P = floor(\sqrt{\sum_i size_i^2})

However, the resulting square can be significantly smaller than this upper bound. Example for this given set:

length123456789
count111111111

We could expect to fill a square of size:

P=floor(285)=16P = floor(\sqrt{285}) = 16

Here we can only use the (9x9) tile:

16_9

All other tile combinations of that set do not result in a square area.

Tooling

Model

Sets

i: tilesp,q:  gridpoints p,q(1..P)p,q:  gridpoints p,q(1..P+1)i: \ \text{tiles}\\ p,q: \ \text{ gridpoints } p,q \in (1..P)\\ p',q': \ \text{ gridpoints } p',q' \in (1..P+1)\\

Parameters/Data

Maximal possible length of area to be filled if all tiles from inventory are used:

P=floor(isizei2)P = floor(\sqrt{\sum_i size_i^2})sizei:  size of tile icoveri,p,q,p,q:  1 if (p’,q’) is covered by tile i at (p,q) size_i: \ \text{ size of tile i}\\ cover_{i,p,q,p',q'}: \ \text{ 1 if (p',q') is covered by tile i at (p,q) }\\

Parameter to indicate coverage of field (p,q)(p',q'):

coveri,p,q,p,q={1,  if (p,q) is covered when tile i is placed at (p,q)0,  else cover_{i,p,q,p',q'} = \begin{cases} 1, \ \text{ if $(p',q')$ is covered when tile i is placed at $(p,q)$}\\ 0, \ \text{ else }\\ \end{cases}

Variables

xi,p,q={1,  if tile i is placed at location (p,q) 0,  else x_{i,p,q} = \begin{cases} 1, \ \text{ if tile i is placed at location (p,q) }\\ 0, \ \text{ else }\\ \end{cases}

Then for each location (p,q)(p,q) we check that the number tiles covering this location is equal to one within the area to fill and zero outside. This results in two additional variables:

Variable to indicate whether the covered field is within current area to fill:

yp,q={1,  if (p,q) is within our current area to fill 0,  else y_{p',q'} = \begin{cases} 1, \ \text{ if $(p',q')$ is within our current area to fill }\\ 0, \ \text{ else }\\ \end{cases}

This matrix yy starts as large as possible, i.e. the number of rows and columns is determined by the max area provided by the inventory or squares. During the optimization this area can shrink. Due to this dynamic yy must be a variable.

Starting value for matrix y is (P+1 x P+1).

The way to keep track of yy is to introduce a binary varialbe:

δ{0,1}δpδp1pδp=W\delta \in \{0, 1\}\\ \delta_p \le \delta_{p-1}\\ \sum_p \delta_p = W\\W:  width of space that can be filled W: \ \text{ width of space that can be filled }\\

Objective

maxWmax W

Constraints

Each tile can be placed max once:

q,pxi,p,q1 i\sum_{q,p} x_{i,p,q} \le 1 \ \forall i\\

Cover Constraint (CC):

Each cell inside the fill area is covered exactly once, outside of the area there must be no coverage.

i,p,qcoveri,p,q,p,qxi,p,q=yp,q p,q\sum_{i,p,q} cover_{i,p,q,p',q'} x_{i,p,q} = y_{p',q'} \ \forall p', q' \\

Pattern constraint for matrix yy:

The way to keep track of yy is to introduce a binary varialbe:

δ{0,1}δpδp1pδp=W\delta \in \{0, 1\}\\ \delta_p \le \delta_{p-1}\\ \sum_p \delta_p = W\\

Constraint for matrix yy:

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.

yp,q=δpδqy_{p,q} = \delta_p \delta_q\\

This non-linear constraint needs to be linearized:

yp,qδpyp,qδqyp,qδp+δqqy_{p',q'} \le \delta_{p'}\\ y_{p',q'} \le \delta_{q'}\\ y_{p',q'} \ge \delta_{p'} + \delta_{q'} -q\\

Pyomo Implementation

The Pyomo implementation is close to the mathematical model. Here an example for the constraint definition:

    # tile can be placed max once:
    model.tile_c = Constraint(
        model.K,
        rule=lambda model, k: sum(model.x[k, i, j] for i in model.I for j in model.J) <= 1
    )

    # cover constraint:
    def all_covered_c(model, ii, jj):
        covers = [key for key, _ in cover.items() if (key[3] == ii and key[4] == jj)]
        return sum(model.cover[k, i, j, ii, jj] * model.x[k, i, j]
                   for (k, i, j, ii, jj) in covers) == model.y[ii, jj]
    model.all_covered_c = Constraint(model.I, model.J, rule=all_covered_c)

Index mapping: cover[k, i, j, ii, jj] corresponds to coveri,p,q,p,qcover_{i,p,q,p',q'}.

Result

Given set:

length1234567
count7654321

Resulting square area: 18_18

It is noted in 1 that the binary variable yp,qy_{p,q} can be relaxed to yp,q[0,1]y_{p,q} \in [0,1]. However, for me this constraint relaxation of yy increased the solution time from 20s to 2:17min. An unexpected observation which I have no explanation for (ideas welcome).

Here another result after 26min running on 6 cores:

length123459
count10108541

Total area is 383, which results in a theoretical maximum square width of 19.

19

Here 2 (1x1) tiles and one (2x2) tile could not be used for the solution, however this is very close to perfection. The side length of the filled area actually hits the theoretical maximum of 19.

Summary

This challenge helped to deepen the understanding of tiling problems. Extending and practicing skills and techniques from Mondriaan Puzzle allowed for a much quicker progress in sovling this puzzle. And again many thanks to Yet Another Math Programming Consultant for publishing his work.

An interesting extension of the tiling principle to a scheduling problem: Patient Scheduling.

If you are interested in the Pyomo model or the Python code contact me via mail.


always be building, by sysid.

© 2022