## 3 ILP Model Formulation
### Assembly of the Linear Systems  
For the reader's convenience, we briefly review how to assemble the linear system of equations that makes up the ILP formulation of the polyomino tiling problem (see [46] for further details) and then discuss the modified procedure required when tiling with parioninoes. Assume we wish to tile a target region \(R\) with area \(c\) using \(F\) free polyominoes \(\{P_{i}\}_{i=1}^{F}\) with areas \(\{c_{i}\}_{i=1}^{F}\) and the numbers of copies \(\{n_{i}\}_{i=1}^{F}\) (\(n_{i}\geq 1\)). The total number of polyominoes used is \(N:=\sum_{i=1}^{F}n_{i}\). Clearly, the total area of \(R\) must equal the sum of the areas of the polyominoes, i.e., \(\sum_{i=1}^{F}c_{i}n_{i}=c\). Using a combination of rotations, reflections and translations, assume each polyomino \(P_{i}\) fits \(s_{i}\) ways into \(R\). Recall the definitions of the binary matrices \(A^{i,j}\) and \(B\) (see Section 3.1). We define _Series_\(i\) to be the set of all binary matrices \(A^{i,j}\) associated with fitting the polyominoes \(P_{i}\) into \(R\), given by  
\[\text{Series }i:=\{A^{i,j}\in\{0,1\}^{r\times c}\mid j=1,\dots,s_{i}\},\quad i =1,\dots,F.\]  
The total number of binary matrices \(A^{i,j}\) is given by \(n:=\sum_{i=1}^{F}s_{i}\), which also corresponds to the total number of possible placements of all polyominoes in the region \(R\). We introduce a variable \(\alpha_{i,j}\in\{0,1\}\), associated with \(A^{i,j}\), for whether the \(j\)-th placement of \(P_{i}\) is used to tile \(R\). The mathematical description of tiling the region \(R\) with copies of the \(F\) free polyominoes in \(P_{i}\) is a natural one: we would like to know if there is a linear combination of the binary matrices \(A^{i,j}\), weighted by the respective \(\alpha_{i,j}\) that is equal to \(B\). In this context, it is helpful to think of the binary matrix \(B\) as representing the fully tiled region \(R\). The problem formulation is as follows. Seek parameters \(\alpha_{i,j}\in\{0,1\}\), \(1\leq i\leq F\), \(1\leq j\leq s_{i}\) such that  
\[\sum_{i=1}^{F}\sum_{j=1}^{s_{i}}\alpha_{i,j}A^{i,j}=B, \tag{3}\] \[\text{subject to }\sum_{j=1}^{s_{i}}\alpha_{i,j}=n_{i},\quad 1 \leq n_{i}\leq s_{i},\quad i=1,\dots,F. \tag{4}\]  
The constraints (4) enforce the conditions that we must use exactly \(n_{i}\) polyominoes from each series \(i\) to tile \(R\). It is easy to verify that the problem formulation automatically enforces the condition that the polyomino areas sum to the area of the target region. Figure 11: Constructing the binary matrices for two placements of a parionino in a checkerboard coloured region \(R\): (**a**) first placement of \(P_{1}^{+}\); (**b**) binary matrix \(\tilde{A}^{1,1}\); (**c**) second placement of \(P_{1}^{+}\); (**d**) binary matrix \(\tilde{A}^{1,2}\). The parioninoes are outlined in red. It is a standard first-year university linear algebra exercise to convert matrix equations of the type (3) into a system of linear equations. After re-arranging the left-hand side of Equation (3) into a single matrix, we equate the non-zero entries of \(B\) with the corresponding entries in this matrix yielding \(c\) linear equations. The ordering of equations corresponds to reading entries in \(B\) left-to-right and top-to-bottom. Together with the constraints (4), this yields the following binary linear system of \(m\) equations in \(n\) unknowns:  
Seek \(\boldsymbol{\alpha}\in\{0,1\}^{n}\) such that  
\[M\boldsymbol{\alpha} =\widetilde{\boldsymbol{b}},\quad M\in\{0,1\}^{m\times n},\quad m =c+F,\quad n=\sum_{i=1}^{F}s_{i},\] (5) and \[\quad\{\widetilde{\boldsymbol{b}}\}_{i} = \begin{cases}1&\text{for }i=1,\ldots,c\\ n_{i}&\text{for }i=c+1,\ldots,m^{\prime},\end{cases} \] \[\boldsymbol{\alpha} =\bigl{(}\alpha_{1,1}\ldots\alpha_{1,s_{1}}\mid\alpha_{2,1} \ldots\alpha_{2,s_{2}}\mid\ldots\mid\alpha_{F,1}\ldots\alpha_{F,s_{F}}\bigr{)} ^{T}.\]  
The vertical lines in the solution vector separate the unknowns according to their corresponding series.