Correction Devoir maison 2

Présentation

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

Essentiellement on nous donne un ensemble de points V avec des coordonnées GPS et un sous-ensemble \(P\subseteq V\). Le but est de trouver un sous-ensemble \(S\subseteq P\) et une affectation \(f: V\rightarrow S\) qui minimise \(500|S|\) plus la somme des distances entre v et f(v) pour tout \(v\in V\). Deux points particuliers sont forcés à être dans S.

Extraire l’ensemble P et calculer les distances

Pour cela on lit le graphe, pour pour chaque sommet \(v \in V\) et chaque sommet \(p\in V\) — dont le degré sortant est suffisamment grand — on affiche dans une ligne le couple de sommets avec leur distance.

from math import sin, cos, asin, sqrt, radians
from sys import argv

# --- données globales

graph = {}         # graph[u] est la liste des sommets v pour tous les arcs (u,v)
GPS = {}           # associe identifiant à un point (latitude, longitude)
selected = []      # liste de centres


def dist(p, q):
    """distance à vol d'oiseau entre deux sommets, en kilomètres (haversine)
    """
    dlon = p[1] - q[1]
    dlat = p[0] - q[0]
    a = sin(radians(dlat / 2))**2 + cos(radians(p[0])) * cos(radians(q[0])) * sin(radians(dlon / 2))**2
    return asin( sqrt(a) ) * 2 * 6373  # le rayon pourrait être ajusté en fonction de la latitude


def read_gpsfile(gpsfile):
    """lit un fichier texte dont les lignes sont dans le format
    v [identificateur] [latitude] [longitude]
    """
    for line in open(gpsfile,'r'):
        if line[0] == 'v':
            tab = line.split()
            identficateur = int(tab[1])
            latitude  = float(tab[2])
            longitude = float(tab[3])
            GPS[identficateur] = (latitude, longitude)
            graph[identficateur] = []
        if line[0] == 'a':
            tab = line.split()
            u = int(tab[1])
            v = int(tab[2])
            graph[u].append(v)
        if line[0] == 's':
            print("Warning: input file should not have lines starting with 's'")


if __name__ == "__main__":
    if len(argv) != 2:
        print("usage: haversine.py <gpsfilename>")
        exit(1)

    read_gpsfile(argv[1])
    for v in GPS:
        for p in GPS:
             if len(graph[p]) >= 4:
               d = dist(GPS[v], GPS[p])
               print("%i %i %.10g" % (v, p, d)

Le programme linéaire

On introduit des variables de décision open[p] pour chaque point \(p\in P\) qui défissent l’ensemble S avec \(p \in S\) ssi open[p]=1. Puis on introduit aussi des variables de décision client[v,p] pour chaque couple de points \(v\in V, p\in P\), avec f(v)=p ssi client[v,p]=1. Il n’y a que deux types de contraintes, une qui dit que tout client \(v\in V\) est bien affecté à un point de service \(p\in P\) et une autre qui dit que l’affectation ne peut se faire que vers les points de service ouverts.

Comme indiqué dans l’énoncé, seuls les variables open sont imposées à être entiers. Notez que ceci est important pour les performances. Si toutes les variables sont imposées à être entières, le temps de résolution serait trop grand. Si aucune variable n’est imposée à être entière, alors la solution au programme linéaire ne serait pas forcément entière.

\( \min \sum_{p\in P} 500 \cdot \textrm{open}[p] + \sum_{v\in V, p\in P} \textrm{dist}[v,p] \cdot \textrm{client}[v,p]) \)

\( \forall v\in V: \sum_{p\in P} \textrm{client}[v,p] = 1 \tag{ matching} \)

\( \forall v\in V, \forall p\in P: \textrm{client}[v,p] \leq \textrm{open}[p] \tag{client_poste_ouverte} \)

\( \forall p\in\{283492455, 283494345\}: \textrm{open}[p]=1 \tag{existing} \)

\( \forall p\in P: \textrm{open}[p]\in\{0,1\}, \forall v\in V, p\in P: \textrm{client}[v,p] \geq 0 \)

Notez: ce n’est pas la peine d’imposer la borne supérieure \(\textrm{client}[v,p] \leq 1\), car elle impliquée par la contrainte matching.

L’implémentation en ZIMPL

# Ecole Centrale Supélec, C. Dürr, 2016
# DM2 - la poste
# facility location
# ZIMPL format

set already := {283492455, 283494345};

set V := { read "dist.txt" as "<1n>"};
set P := { read "dist.txt" as "<2n>"};
param dist[C cross P] := read "dist.txt" as "<1n, 2n> 3n";

var open[P] binary;          # poste ouverte
var client[V cross P] >= 0;  # client associé à une poste

minimize objective: sum <p> in P: 500 * open[p]
                  + sum <v> in V: sum <p> in P: dist[v,p] * client[v,p] ;

subto matching:  # chaque homme est client d'exactement un bureau de poste
    forall <v> in V:
        sum <p> in P: client[v,p] == 1;

subto client_poste_ouverte:
    forall <v> in V:
        forall <p> in P:
            client[v,p] <= open[p];

subto existing:
    forall <p> in already:
        open[p] == 1;

Une solution

SCIP calcule une solution en moins d’une minute. Elle consiste en l’ouverture de 13 bureaux et atteint une valeur objectif d’à peu près 13 915. La visualisation de la solution montre de nombreux bureaux ouverts dans la partie la plus dense de l’île, comme attendu.

Preuve d’intégralité

Soit open, client une solution point extrême au programme linéaire. Les variables open ont été imposés entiers. Nous allons montrer que les variables client sont également entiers. Pour cela il suffit de montrer que le programme linéaire restreints aux variables client est entier pour les variables open entiers fixées.

Soit S l’ensemble décrit par open, c’est-à-dire S est l’ensemble des points \(p\in P\) tel que open[p]=1. Alors par la contrainte client_poste_ouverte dans le programme linéaire les variables client[v,p] sont forcées à 0 pour tout \(p\in P \setminus S\).

Finalement on a le programme linéaire suivant

\( \min \sum_{v\in V, p\in S} \textrm{dist}[v,p] \cdot \textrm{client}[v,p]) \)

\( \forall v\in V: \sum_{p\in S} \textrm{client}[v,p] = 1 \tag{matching} \)

\( \forall v\in V, p\in S: \textrm{client}[v,p]\geq 0 \)

Notez: aucune contrainte ne relie des variables client[v,p] et client[v',p'] pour deux sommets distincts \(v,v’\in V\). Les variables de la forme client[v,*] forment donc des sous problèmes indépendants. La solution optimale d’un de ces sous problème consiste à affecter v à un des bureaux de postes p le plus proche et ouvert.

Preuve directe

Supposont qu’il existe une variable client[v,p] qui ne soit pas entière pour des sommets \(v\in V, p\in P\). Nous avons 0 < client[v,p] < 1. Alors par la contrainte matching il existe un autre sommet \(p’\in P\) avec 0 < client[v,p'] < 1. Soit \(\epsilon>0\) une valeur suffisamment petite, disons \(\min\{\textrm{client}[v,p], \textrm{client}[v,p’], 1-\textrm{client}[v,p],1-\textrm{client}[v,p’]\}\). Alors ajouter \(\epsilon\) à client[v,p] et enlever \(\epsilon\) de client[v,p'] préserve toutes les contraintes. Similairement ajouter \(\epsilon\) à client[v,p'] et enlever \(\epsilon\) de client[v,p] préserve toutes les contraintes. Ceci montre que open, client n’est pas une solution point extrême au programme linéaire.

Preuve alternative

D’abord reformulons le programme linéaire afin de n’avoir que des inégalités.

\( \min \sum_{v\in V, p\in S} \textrm{dist}[v,p] \cdot \textrm{client}[v,p]) \)

\( \forall v\in V: \sum_{p\in S} \textrm{client}[v,p] \geq 1 \tag{matching-lower-bound} \)

\( \forall v\in V: \sum_{p\in S} -\textrm{client}[v,p] \geq -1 \tag{matching-upper-bound} \)

\( \forall v\in V, p\in S: \textrm{client}[v,p]\geq 0 \)

Alors la matrice contraintes-variables a dans la colonne associée à une variable client[v,p] seulement deux entrées non-nulles, une étant +1 et une étant -1. Il s’agit donc de la transposée d’une matrice d’incidence d’un graphe orienté. Et il est connu qu’une telle matrice est totalement unimodulaire (TUM). Et on sait qu’un polytope est entier s’il est décrit par un programme linéaire de la forme \(Ax \geq b \) avec b entier et A TUM.

Autres preuves

On pourrait aussi observer que les contraintes matching-upper-bound ne sont pas utiles pour les solutions optimales, car si pour un sommet v on a \(\sum_{p\in S} \textrm{client}[v,p] > 1\), alors il est possible de diminuer une des variables, préserver les contraintes et en même diminuer la valeur objectif. Le programme linéaire sans la contrainte matching-upper-bound a une matrice contraintes-variables qui contient des valeurs non-nulles (+1 en fait) pour chaque ligne v, seulement dans dans des colonnes consécutives si les variables sont ordonnées dans l’ordre lexicographique (v,p). Ces matrices sont aussi connues pour être totalement unimodulaires.