Pavage 
Introduction
On souhaite paver (sans chevauchement) une grille rectangulaire en la remplissant au maximum à l'aide de tuiles (chacune disponible en quantité infinie).
Voici un exemple de pavage optimal (pourquoi est-il optimal?) d'un carré 5x5 à l'aide de tuiles en L (dans n'importe quel sens) :
Dans cet exemple, seule la case en haut à gauche n'est pas couverte.
Modélisation par PLNE
Considérons l'ensemble $\mathcal{P}$ des positions possibles d'une tuile sur la grille (exercice: que vaut $\vert \mathcal{P} \vert$ sur l'exemple ci-dessus?).
Définissons une variable $x_k$ par élément $p_k$ de $\mathcal{P}$, valant 1 si la tuile correspondante est posée, 0 sinon.
Définissons $w_k$ comme le nombre de cases occupées par la tuile $p_k$.
On souhaite alors maximiser le nombre de cases couvertes, c'est-à-dire:
$$\sum_k x_k \times w_k$$
La contrainte de non-chevauchement s'exprime en disant que chaque case de coordonnée $(i, j)$ est recouverte par au plus une tuile, c'est-à-dire:
$$\displaystyle\sum_{p_k \in \mathcal{P},\\(i, j) \in p_k} x_k \leq 1$$
Il s'agit d'un programme linéaire en nombres entiers (PLNE).
Résolution avec MIP
Malheureusement, la résolution d'un PLNE est plus difficile que celle d'un programme linéaire où les variables sont réelles.
En effet, dans ce deuxième cas, la méthode de l'ellipsoïde/du simplexe donne un algorithme efficace en théorie/en pratique.
Heureusement, il existe quand même des méthodes pour résoudre un PLNE comme celle de Branch & Bound (et ses variations comme Branch & Cut). Le package MIP en Python permet d'utiliser un solveur de PLNE s'appuyant sur ce genre de méthodes.
Par défaut, MIP utilise le solveur COIN-OR.
from mip import *
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({
"axes.grid": True,
"image.aspect": "equal"
})
Définissons un PLNE avec MIP.
Pour cela, on commence par définir un modèle m.
Les variables (binaires: 0 ou 1) sont ajoutées avec m.add_var.
On ajoute une contrainte en écrivant m += ... <= ....
constraints[i][j] contient la liste des variables (i.e tuiles) couvrant la case $(i, j)$.
def pl(tiles, n, p, var_type): # renvoie un modèle de PL pour paver une grille nxp avec des tuiles de tiles
m, x, w = Model(), [], []
constraints = [[[] for j in range(p)] for i in range(n)]
for t in range(len(tiles)):
a, b = tiles[t].shape
for i in range(n - a + 1):
for j in range(p - b + 1):
x.append(m.add_var(var_type=var_type, name=str(t)+","+str(i)+","+str(j)))
w.append(tiles[t].sum())
for k, l in np.argwhere(tiles[t]):
constraints[i + k][j + l].append(x[-1])
for i in range(len(constraints)):
for j in range(len(constraints[i])):
m += xsum(constraints[i][j]) <= 1
m.objective = maximize(xsum(x[i]*w[i] for i in range(len(x))))
return m
Une fois ce modèle obtenu, on peut le résoudre:
def solve(tiles, n, p, var_type=BINARY): # affiche un pavage optimal d'une grille nxp avec des tuiles de tiles
m = pl(tiles, n, p, var_type)
img = np.zeros((n, p))
m.optimize(max_seconds=60)
ntile = 0
for v in m.vars:
if abs(v.x) > 1e-6:
t, i, j = map(int, v.name.split(",")) # tuile correspondante à la variable v
ntile += v.x
for k, l in np.argwhere(tiles[t]):
img[i + k, j + l] = ntile
print(f"Solution optimale avec {ntile} tuiles")
if var_type == BINARY:
print(f"{(img == 0).sum()} cases ne sont pas couvertes")
plt.imshow(img, cmap = "gist_ncar", extent=[0,p,0,n])
Essayons:
L = np.full((2, 2), True)
L[1, 1] = False # L[i, j] vaut True si (i, j) est occupé par la tuile, False sinon
Ls = [L, L[::-1, :], L[:, ::-1], L[::-1, ::-1]] # les 4 rotations possibles de L
solve(Ls, 5, 5)
On retrouve bien l'exemple initial !
On peut constater que la relaxation continue (où on peut sélectionner une tuile un nombre fractionnaire de fois, par exemple $0.5$ fois) a le même optimum :
solve(Ls, 5, 5, var_type=CONTINUOUS)
Essayons avec d'autres tuiles, par exemple avec des carrés 2x2:
C = np.full((2, 2), True)
solve([C], 5, 5)
Avec des bâtons 4x1:
I = np.full((4, 1), True)
Is = [I, I.T]
solve(Is, 5, 5)
Avec un peu n'importe quoi:
solve([C, I, L], 11, 21)
solve([C, I, L], 11, 21, var_type=CONTINUOUS)
On voit sur l'exemple ci-dessus que l'optimum du PLNE n'est pas égal à l'optimum de sa relaxation.
Dual
Écrivons le dual du PL des pavages. Pour cela, on définit une variable duale $y_{i, j}$ pour chaque contrainte de $P$ et on obtient les contraintes suivantes : $$\min \sum_{(i, j)} y_{i, j}$$ $$\forall \ell \in \mathcal{P}, ~\sum_{(i, j) \in \ell} y_{i, j} \geq 1$$ $$\forall (i, j), ~y_{i, j} \in \{0, 1\}$$ Interprétation : le dual demande de choisir un nombre minimum de cases (des << pièges >>) pour empêcher tout placement de L sur la grille.
def pl_dual(tiles, n, p, var_type):
m = Model()
y = np.array([[m.add_var(var_type=var_type, name=str(i)+","+str(j)) for j in range(p)] for i in range(n)])
for t in tiles:
a, b = t.shape
for i in range(n - a + 1):
for j in range(p - b + 1):
ind = np.array([i, j]) + np.argwhere(t)
pos = y[ind[:, 0], ind[:, 1]].flatten()
m += xsum(pos) >= 1
m.objective = minimize(xsum(y.flatten()))
return m
def solve_dual(tiles, n, p, var_type=BINARY):
m = pl_dual(tiles, n, p, var_type)
img = np.ones((n, p))
m.optimize(max_seconds=60)
n = 0
for v in m.vars:
if abs(v.x) > 1e-6:
i, j = map(int, v.name.split(",")) # tuile correspondante à la variable v
n += v.x
img[j, i] = 0
plt.text(i + .35, j + .35, str(v.x), color="white")
print("Solution optimale avec", n, "pièges")
plt.imshow(img, cmap = "gray", extent=[0,5,0,5])
On trouve que le nombre minimum de pièges à poser sur une grille $5\times 5$ pour empêchement le placement d'un L est $10$ :
solve_dual(Ls, 5, 5, BINARY)
Il n'y a donc pas égalité de l'optimum du PLNE primal et du PLNE dual.
Par contre, en accord avec le théorème de dualité forte, les relaxations du primal et dual ont le même optimum, à savoir 8 :
solve_dual(Ls, 5, 5, CONTINUOUS)
solve(Ls, 2, 10)
solve(Ls, 2, 11)
solve(Ls, 2, 12)
solve(Ls, 2, 13)
solve(Ls, 16, 16)
solve(Is, 10, 10)
