Correction Devoir maison 2 - Exposition au musée

Présentation

L’énoncé complet du problème est ici.

Essentiellement on nous donne une liste de n points dans le plan. Le but est de les couvrir avec des disques . Deux rayons sont possibles 4 et 8, qui coûtent respectivement 1 et 2 Euros. Le but est de produire une couverture avec un coût minimal.

Modélisation

Le problème est un problème de set covering avec un aspect géométrique. L’univers est l’ensemble des points données, et chaque disque i correspond à un ensemble de points \(S_i\) et a un coût \(w_i\). Le but est de trouver une collection des ensembles de moindre coût. Le programme linéaire modélisant ce problème est fort simple. Soit U l’ensemble des points et C la collection des ensembles potentiels.

\[ \min \sum_{i\in C} w_i \cdot x_i \] tel que \[ \forall p\in U: \sum_{i\in C: p\in S_i} x_i \geq 1 \] et \[ \forall i\in C: x_i \in \{ 0, 1 \}. \]

Discrétisation

Il existe un nombre infini de disques de rayon fixé dans le plan. Pour avoir un nombre fini et assez petit de disques dans C on restreint le choix. On dit qu’un disque est bon si soit il contient un seul point de U et que celui-ci est le centre, soit il contient au moins deux points de U et ces points sont sur son bord.

Sans perte de généralité une solution optimale n’est composée que de bons disques. Pour cela il suffit de translater les disques pour qu’ils deviennent bons. Par exemple soit un disque qui contient plusieurs points, et aucun n’est sur son bord. On peut le translater dans une direction arbitraire jusqu’à ce qu’un point soit sur son bord. À ce moment on peut faire une rotation du disque autour de ce point, jusqu’à ce qu’un deuxième point soit sur son bord. Ces opérations préservent la solution.

Nous générons C alors de la manière suivante. D’abord chaque point p on génère un disque avec centre p et rayon 4 (le rayon le moins cher). Puis pour chaque rayon r et chaque paire de points p1, p2 à distance au plus 2r, on génère tous les deux bons disques définis par p1,p2,r. Dans le cas dégénéré où la distance est exactement 2r les deux disques sont identiques et on pourrait en générer qu’un seul. Cette optimisation n’a pas été implémentée.

Voici comment on détermine le centre d’un disque défini par deux points et un rayon. Les deux ingrédients géométriques sont le théorème de Pythagore et le fait qu’une rotation de 90 dégrées d’un vecteur (a,b) soit (-b,a).

def dist2(p, q):
    """distance L2 au carré entre deux points
    """
    dx = p[0] - q[0]
    dy = p[1] - q[1]
    return dx * dx + dy * dy


def circle(p1, p2, r):
    """retourne le centre d'un des deux cercles de rayon r définis par les points p1, p2 sur le cercle.
    Quand p1 et p2 sont échangés, la fonction devrait retourner l'autre cercle.
    """
    x1, y1 = p1
    x2, y2 = p2
    bx = (x1 + x2) / 2.
    by = (y1 + y2) / 2.
    q2 = dist2( p1, p2 )
    q = sqrt(q2)
    x = bx + sqrt( r * r - q2 / 4. ) * (y1 - y2) / q
    y = by + sqrt( r * r - q2 / 4. ) * (x2 - x1) / q
    return x, y

Filter pour déterminer les points proches

Étant donné une distance d et un point p1 on cherche tous les points p2 de distance au plus d de p1. Cette étape est importante pour la génération des bons disques ainsi que pour déterminer la liste de points couverts par une caméra. Pour ne pas boucler sur tous les $n^2$ paires plusieurs techniques sont possibles. Le plus simple peut-être de définir une grille booléenne, qui indique pour chaque point de coordonnées entière la présence d’une œuvre d’art. Il suffit alors de chercher p2 dans le carré de côté 2d centré autour de p1. Cependant nous proposons ici une technique alternative, qui fonctionnerait également si les points d’entrée n’avaient pas des coordonnées entières.

L’idée consiste à placer une grille au dessus du plan. Chaque point se trouve dans une case de la grille. Il suffit de chercher les points p2 (en bleu clair ci-dessous, et artificiellement agrandi pour la visibilité) dans les 9 cases autour de la case contenant le point p (en bleu foncé). Ce filtrage améliore les performances de la génération de C. Pour un rayon constant, et une distribution uniforme des points dans un rectangle, cette technique permet de répondre à la requête des points proches en temps linéaire en n plutôt que quadratique pour la technique naïve.

Le pas de la grille (côté d’une cellule) est choisi juste un peu plus grand que deux fois le plus grand des rayons.

from collections import defaultdict

grid = defaultdict(set)  # associe à une coordonnée de cellule de grille, la liste de ses points


def make_disks():
    print( "2. generer les emplacements potentiels" )
    global grid, camera
    # ranger les points dans une grille
    for (i, (x, y)) in enumerate( art ):
        cell = (x // size, y // size)
        grid[cell].add( i )
    # boucler sur les paires de points proches, donc issues de cellules voisines
    for i1, (x1, y1) in enumerate( art ):
        camera.append( art[i1] + (0,) )           # bon disque avec i1 au centre
        cx1 = x1 // size
        cy1 = y1 // size
        for cx2 in [cx1 - 1, cx1, cx1 + 1]:       # chercher point i2 dans les cellules voisines
            for cy2 in [cy1 - 1, cy1, cy1 + 1]:
                for i2 in grid[cx2, cy2]:
                    if i1 != i2:
                        for t in [0, 1]:          # tous les rayons
                            if dist2( art[i1], art[i2] ) <= (2 * radius[t]) ** 2:
                                x, y = circle( art[i1], art[i2], radius[t] )
                                camera.append( (x, y, t) )

Enfin la grille est encore utilisée pour déterminer les ensembles, donc pour chaque disque de C l’ensemble des points qu’il contient.

def detect_inclusions():
    print( "3. detecter les inclusions" )
    global camera_over_art
    camera_over_art = [[] for _ in art]
    for i, (x, y, t) in enumerate( camera ):
        cx = int( x ) // size
        cy = int( y ) // size
        for cx2 in [cx - 1, cx, cx + 1]:
            for cy2 in [cy - 1, cy, cy + 1]:
                for j in grid[cx2, cy2]:
                    # toujours ajouter une tolérance quand on compare des flottants
                    if dist2( (x, y), art[j] ) <= radius[t] * radius[t] + 1e-10:
                        camera_over_art[j].append( i )

Appeler le solveur de programme linéaire à variables entières

La création d’un modèle pour SCIP et l’extraction d’une solution se fait comme suit.

def solve_MIP():
    print( "4. creer le modele MIP" )
    global sol
    m = Model()
    var = []
    for k, (x, y, t) in enumerate( camera ):
        var.append( m.addVar( "v%i" % k, vtype='B', obj=price[t] ) )
    for i, can_cover in enumerate( camera_over_art ):
        m.addCons( quicksum( var[k] for k in can_cover ) >= 1 )
    m.setMinimize()
    print( "5. appeler le solveur SCIP" )
    m.optimize()
    print( m.getObjVal() )
    print( m.getStatus() )
    print( "6. ecrire la solution" )
    f = open( "sol.txt", "w" )
    sol = set()
    for k, v in enumerate( var ):
        if m.getVal( v ) >= 0.5:  # devrait être 0 ou 1, mais il peut avoir des erreurs d'arrondis
            print( "{0},{1},{2}".format(camera[k][2] + 1, camera[k][0], camera[k][1]), file=f )
            sol.add( k )

Performances

Le ficher d’entrée contient certains points en double. Après élimination des doublons il reste 4982 points. Ceux ci forment 43036 bons disques et autant de variables de décision. L’exécution du programme, la résolution par SCIP comprise, prend 5 secondes sur un MacBook 1,6 GHz Intel Core i5. Le score obtenu est 2586.

La recherche locale

Une possibilité que j’avais en tête consiste à choisir une fenêtre de petite taille, peut-être 30 fois 30 et de la place aléatoirement sur la grille. Puis on enlève tous les disques inclus dans la fenêtre et on calcule la solution optimale pour couvrir les points ainsi découverts. On espère que l’espace de solution est petit, par exemple \(2^{20}\) s’il y a 20 disques pouvant entrer dans la solution, et ainsi avoir un temps de ré-optimisation locale de l’ordre de la seconde. Cette direction n’a pas été expérimentée.

Références