Découverte du solveur GLPK

GNU/Linux Magazine n° 135 | février 2011 | Damien Balima
Creative Commons
  • Actuellement 0 sur 5 étoiles
0
Merci d'avoir participé !
Vous avez déjà noté cette page, vous ne pouvez la noter qu'une fois !
Votre note a été changée, merci de votre participation !
GLPK, the GNU Linear Programming Kit, le solveur open source qui simplifie la résolution de problèmes complexes et réconcilie les utilisateurs avec la programmation linéaire.

1. Introduction

Il arrive souvent qu'un problème difficile à résoudre se pose dans le développement d'un programme, d'un système ainsi que dans la vie courante. Par exemple, dans un lycée, il est difficile de réaliser des emplois du temps pour une cinquantaine de professeurs et une vingtaine de classes en respectant les contraintes humaines et matérielles.

De même, il est difficile d'ordonnancer un nombre important d'éléments lorsqu'on dispose de ressources limitées, de façon à maximiser ou minimiser le coût. Par exemple, lorsque l'on souhaite acquérir à moindre coût et le plus rapidement possible des fournitures scolaires, sachant que l'on dispose de deux magasins, l'un offrant un grand choix de produits de marques et ayant un temps d'attente en caisse assez long, et l'autre, moins cher, mais ne disposant pas de tous les produits nécessaires et étant plus rapide en caisse.

En optimisation combinatoire, la théorie de la complexité des problèmes de décision définit plusieurs classes de problèmes. Les problèmes de décision les plus simples, pouvant être résolus en temps polynomial par rapport à ses données d'entrée, sont des problèmes de la classe P(P comme temps Polynômiaux). Par exemple, savoir s'il existe ou non assez de place pour tous les invités.

Les problèmes de décision plus complexes, ne pouvant pas être résolus en temps polynomial, font partie de la classe NP(comme Polynomiaux Non déterministes). Ce sont des problèmes dont la solution est facile à vérifier par une machine de Turing non déterministe. C'est toute la différence entre résoudre un problème et vérifier sa solution, il est souvent plus facile de vérifier sa solution que de le résoudre. Mais est-ce qu'un problème facile à vérifier n'est pas, dans l'absolu, un problème facile à résoudre ? Est-ce que P=NP ? C'est une question à un million de dollars, un des problèmes du prix du millénaire, qui n'est toujours pas résolu aujourd'hui. Aussi, n'ayant pas prouvé que tout problème NP équivaut à un problème P, la conjecture actuelle est que l'ensemble P est inclus dans NP.

Les problèmes de la classe NP-Complet, en plus d'être NP, sont réductibles à tout autre problème NP. Par exemple, les problèmes des emplois du temps et des fournitures scolaires cités plus haut, qui deviennent difficiles dès lors que l'on commence à augmenter le nombre de leurs variables, sont NP-Complet ; ils peuvent se réduire à un problème NP bien connu de sac à dos (ou KP, Knapsack Problem) ou decouverture d'ensemble (Set Covering)[1] [2].

Pour les résoudre, le choix de la méthode dépend du contexte. Une première méthode qui va rapidement trouver une solution approchée du problème consiste à développer un algorithme, ou métaheuristique, basé, par exemple, sur un algorithme glouton, sur le recuit simulé, sur la méthode Tabou ou sur un algorithme génétique [3] [4].

Une seconde méthode est d'utiliser un solveur, ou sa bibliothèque, qui utilisera une batterie d'outils et de méthodes de calcules (algorithme du simplexe, branch and bound) pour trouver la solution exacte, optimale, du problème après un temps de calcul plus ou moins long. Comme il s'agit de problèmes avec une solution optimale, des variables et des contraintes spécifiques, la façon la plus simple de procéder consiste à modéliser un programme linéaire(PL) qui sera résolu par un solveur ad hoc. On formalise les programmes linéaires en indiquant la formule linéaire du problème, suivie des contraintes linéaires (de forme a.x+b=0, à la différence de la programmation non linéaire). Dans le cas où les données et variables sont entières, il s'agit d'un programme linéaire en nombres entiers (PLNE).

En pratique, s'il est impératif de trouver le résultat optimal, sans contrainte de temps de calcul, on se tournera vers un solveur, tandis que s'il importe juste de trouver un résultat approché, mais le plus rapidement possible, on choisira les métaheuristiques. Concernant les solveurs, GLPK, the GNU Linear Programming Kit, est un solveur de programmes linéaires, libre et open source, qui contient aussi une bibliothèque pouvant être utilisée par un programme C ou C++. On trouve également des bibliothèques GLPK pour Java, Python et R.

Dans la suite de cet article, nous nous intéresserons à l'application GLPK et tenterons de résoudre quelques problèmes NP simples.

2. GLPK et son langage

GLPK peut résoudre des problèmes linéaires dans plusieurs langages de programmation linéaire. Nous verrons en particulier le langage GMPL, le GNU Mathematical Programming Language.

2.1 Utilisation de GLPK

Les programmes linéaires (PL) sont structurés en plusieurs parties : les données, les variables, la fonction d'optimisation et les contraintes. La commande glpsol lance le calcul après vérification de l'intégrité du PL. Les langages MPS, Free MPS, CPLEX LP, GLPK et GMPL sont pris en charge respectivement par les options mps, freemps, lp glp et math(man glpsol).

En langage CPLEX LP, Minimize indique la fonction a minimiser (Maximize pour maximiser), Suject To indique la suite des contraintes (dans l'exemple ci-dessous : c1 pour la constante d, et c2 pour f=0), General pour les variables entières (Binary pour les variables binaires) et End pour la fin du programme linéaire. On peut également ajouter une section Bounds de variables bornées.

Par exemple, au format CPLEX LP, pour résoudre l'équation 3x+5y-6z+3=0 avec x, y et z des nombres entiers, on écrit le fichier ex1.lp suivant :

Minimize                                                            

        f: 3x + 5y – 6z + d                                             

                                                                    

Subject To                                                          

        c1: d = 3                                               

        c2: 3x + 5y – 6z + d = 0                                        

General                                                             

        x                                                           

        y                                                           

        z                                                           

End             

Puis on lance la commande :

glpsol --lp ex1.lp --output ex1.sol

La sortie du terminal nous renseigne sur le fait qu'une solution optimale a rapidement été trouvée :

[...]

INTEGER OPTIMAL SOLUTION FOUND

Time used:   0.0 secs

Memory used: 0.1 Mb (64425 bytes)

Writing MIP solution to `ex1.sol'...

On obtient les résultats suivants (ex1.sol) :

    No. Column name       Activity     Lower bound   Upper bound

------ ------------    ------------- ------------- -------------

     1 x            *              0             0               

     2 y            *              3             0               

     3 z            *              3             0               

     4 d                           3             0   

Les colonnes Activity et Lower bound indiquent la solution S={0;3;3}.

2.2 Le langage GMPL

Les données peuvent être incorporées au début du programme linéaire ou être placées dans un fichier à part (option data) appelé lors de la commande glpsol. Le langage requis avec l'option data est le GMPL. L'option model indique le programme linéaire à utiliser, également en GMPL. Par exemple :

$ glpsol --model program.mod --data donnees.dat --output resultats.txt

Fournis avec GLPK, la documentation GMPL explique la syntaxe du langage. En voici les principaux points :

2.2.1 Les données

Les données peuvent être placées dans un fichier à part ou directement dans le PL. Un fichier de données externes commence par data, se termine par end et ne contient que des données (i.e.pas d'instructions). Les points-virgules sont à placer en fin de déclaration. Pour affecter une valeur à un paramètre, on utilise l'opérateur d'affectation:=.

Pour les colonnes d'un tableau, d'une matrice, on place deux-points après le nom du paramètre, et les indices des colonnes avant l'opérateur :=. La dernière ligne d'une matrice se termine également par un point-virgule. Les commentaires sont à placer après un croisillon #.

Les données peuvent être déclarées sous forme de paramètres (param) ou d'ensembles (set).

Exemple de fichier de données que l'on nomme donnees.dat :

data;

#nombre de villes

param n := 4;

#distance entre les villes

param d: 1 2 3 4 :=

      1 0 7 3 5

      2 0 0 0 0

      3 0 8 0 4

      4 0 4 0 0 ;

end;

Dans l'exemple ci-dessus, une cellule dij se définit par sa ligne i et sa colonne j. Ainsi, d12=7 et d43=5.

2.2.2 Le programme linéaire

Le programme linéaire est structuré en plusieurs parties. En premier lieu, on indique l'objectif maximal ou minimal à calculer, par les mots-clés maximize ou minimize. L'objectif est facultatif, suivant le problème.

Puis l'on indique les contraintes linéaires sous la section suject to. Les contraintes commencent par leurs noms, puis optionnellement la suite de conditions de type « pour tout x » entre accolades, suivis de deux-points et de la fonction linéaire. Les conditions peuvent être définies sur un ensemble borné grâce aux opérateurs in et ... De plus, on peut affiner la condition après deux-points. Par exemple, pour indiquer la condition « pour tout i allant de 1 à 10 et tout j allant de 1 à 10, avec i et j différents », on notera : {i in 1..10, j in 1..10 : i!=j}.

Après la section des contraintes, on ajoute la commande de calcul solve et les instructions d'affichage comme printf et display. Le fichier de programme linéaire se termine par une instruction end.

2.2.3 Exemple

Par exemple, avec les données précédentes, pour trouver le plus court qui relie les quatre sommets, on minimise la somme des distances avec comme contraintes de passer au moins une fois par tous les sommets et par des chemins non nuls.

Dans un fichier qui sera utilisé comme modèle pour traiter les données, on indique au début les définitions des données externes à utiliser.

# donnees externes

param n; #nombre de villes

param d{i in 1..n, j in 1..n}; #chemins

La fonction d'optimisation est donc de minimiser la somme des distances. Nous utiliserons donc la matrice d, et une variable binaire xijqui indique si oui ou non on passe par le chemin dij. L'objectif revient ainsi à minimiser la somme (sum) des chemins dijmultipliée par xij, que l'on note :

minimize f: sum {i in 1..n, j in 1..n} d[i,j]*x[i,j];

On ajoute les contraintes sous un bloc Suject to. Une contrainte chemins pour indiquer les n-1 chemins, une contrainte croisements pour n'utiliser qu'un ou deux chemins par ville et une contrainte non_nul pour n'utiliser que des chemins de valeurs non nuls (si d=0 et x(d-1)>=0 alors x=0 oblige).

On obtient le bloc de contraintes suivant :

subject to

#un chemin passant par toutes les villes sans retour initial

chemins : sum {i in 1..n, j in 1..n} x[i,j] = n-1;

#pour chaque ville il ne peut y avoir qu'un ou deux chemins

croisements { i in 1..n} : sum {j in 1..n : i!=j} x[i,j]+sum{j in 1..n : i!=j} x[j,i] <= 2;

#chemin non nul

non_nul {i in 1..n, j in 1..n} : x[i,j]*(d[i,j]-1) >= 0;

Enfin, on lance le calcul avec la fonction solve et l'on affiche les résultats avec display x.

Ce qui nous donne le fichier de programme linéaire program.mod suivant :

#modele du plus court chemin

# donnees externes

param n; #nombre de villes

param d{i in 1..n, j in 1..n}; #chemins

#variable x = si chemin xij est emprunté

var x{i in 1..n, j in 1..n} binary;

# objectif

minimize f: sum {i in 1..n, j in 1..n} d[i,j]*x[i,j];

# contraintes

subject to

#un chemin passant par toutes les villes sans retour initial

chemins : sum {i in 1..n, j in 1..n} x[i,j] = n-1;

#pour chaque ville il ne peut y avoir qu'un ou deux chemins

croisements { i in 1..n} : sum {j in 1..n : i!=j} x[i,j]+sum{j in 1..n : i!=j} x[j,i] <= 2;

#chemin non nul

non_nul {i in 1..n, j in 1..n} : x[i,j]*(d[i,j]-1) >= 0;

printf "---Debut de la resolution ---\n";

solve;

printf "---Fin de la resolution ---\n";

display x;

end;

On lance la commande glpsol :

$ glpsol --model program.mod --data donnees.dat

Et l'on obtient les résultats suivants :

[...]

INTEGER OPTIMAL SOLUTION FOUND

[...]

x[1,1] = 0

x[1,2] = 0

x[1,3] = 1

x[1,4] = 0

x[2,1] = 0

x[2,2] = 0

x[2,3] = 0

x[2,4] = 0

x[3,1] = 0

x[3,2] = 0

x[3,3] = 0

x[3,4] = 1

x[4,1] = 0

x[4,2] = 1

x[4,3] = 0

x[4,4] = 0

Les valeurs de xij égales à 1 indiquent le chemin C={(1;3)(3;4)(4;2)}, qui est bien le plus court chemin qui relie toutes les villes.

Ce problème que l'on vient de résoudre est du type du Voyageur de commerce, un problème NP-Complet dont l'objectif est de savoir s'il existe un plus court chemin qui relie tout les points d'un graphe. Avec 4 villes, ça semble assez simple, mais avec une vingtaine de villes ou plus, ça devient nettement plus difficile si l'on n'a pas ce genre de solveurs sous la main.

3. Problèmes

Voyons à présent deux problèmes classiques. Le premier est celui du sac-à-dos, dont le principe est de pouvoir optimiser l'allocation de la place avec une taille maximale donnée. Plus concrètement, il s'agit de choisir les bonnes boites pour rentabiliser au mieux le port du sac-à-dos. Le second consiste en la résolution d'un Sudoku, qui n'est pas à proprement parler NP-Complet dans sa version 9x9 (il peut être résolu en temps polynomial) mais qui offre cependant un bon exemple de programmation par contraintes.

3.1 Le problème du sac à dos

Le problème qui se pose à Tux est assez simple : il s'agit de remplir son sac à dos avec le plus de blocs de glace pour un maximum de valeur totale, sachant que son sac à dos est limité à 15 litres.

En données, nous avons 7 blocs de glace : 2 blocs de 2 litres à 180 € l'unité, 2 de 3 litres à 200 € l'unité, 2 de 4 litres à 220 € l'unité et un bloc de glace de 6 litres à 420 €. On utilisera donc deux matrices t et v, l'une renseignant sur la taille des blocs, l'autre sur leurs valeurs. Tux dispose d'un sac à dos limité à 15 litres, que nous noterons M, comme volume Maximum.

Ce qui nous donne le fichier de données sac_a_dos.dat suivant :

data;

#volume maximum du sac-à-dos

param M := 15;

#nombre de blocs

param n := 7;

#taille des blocs

param t := 1 2

                   2 2

                   3 3

                   4 3

                   5 4

                   6 4

                   7 6;

#valeur des blocs

param v := 1 180

                    2 180

                    3 200

                    4 200

                    5 220

                    6 220

                    7 420;

end;

Au niveau du programme linéaire, dans un fichier sac_a_dos.mod, on reprend les mêmes principes que dans l'exemple du plus court chemin. En premier lieu, on indique les données externes à utiliser :

#donnees externes

param M; #capacite du sac

param n; #nombre d'objets

param t{i in 1..n}; #taille du bloc i

param v{i in 1..n}; #valeur du bloc i

On utilise également une variable xi binaire qui vaut 1 si Tux prend le bloc i.

#variable xi=1 si bloc i est choisi

var x{1..n} binary;

Le problème est de choisir les blocs de manière à obtenir le maximum de valeur totale, sachant que l'on dispose d'un volume limité. On passera cette limitation dans les contraintes, il nous reste la fonction d'optimisation, qui consiste à maximiser la somme des valeurs (vi) des blocs choisis (xi), c'est-à-dire maximiser la somme des vi facteur xi, pour i allant de 1 à n.

#objectif

maximize f :sum {i in 1..n} v[i]*x[i];

Il nous reste les contraintes. Le sac à dos ne peut contenir plus de 15 livres, donc la somme des tailles (ti) des blocs sélectionnés (xi) doit être inférieure ou égale à ce volume maximum M.

#contraintes

subject to

capacite : sum{i in 1..n} t[i]*x[i] <= M;

Au final, on obtient le PL suivant :

#modele du sac a dos

#donnees

param M; #capacite du sac

param n; #nombre d'objets

param t{i in 1..n}; #taille du bloc i

param v{i in 1..n}; #valeur du bloc i

#variable xi=1 si bloc i est choisi

var x{1..n} binary;

#objectif

maximize f :sum {i in 1..n} v[i]*x[i];

#contraintes

subject to

capacite : sum{i in 1..n} t[i]*x[i] <= M;

printf "---Debut de la resolution ---\n";

solve;

printf "---Fin de la resolution ---\n";

display x;

end;

On lance la commande glpsol :

$ glpsol --model sac_a_dos.mod --data sac_a_dos.dat

On obtient les résultats suivants :

[...]

INTEGER OPTIMAL SOLUTION FOUND

[...]

x[1] = 1

x[2] = 0

x[3] = 0

x[4] = 1

x[5] = 0

x[6] = 1

x[7] = 1

Soit quatre blocs de glace : un bloc de 2 litres, un bloc de 3 litres, un bloc de 4 litres et le bloc de 6 litres, pour un total de 15 litres à 1020 €. Finalement, ce problème fut plus simple que celui du voyageur de commerce. Les problèmes NP-Complet du sac à dos et de sa variante du sac à dos multidimensionnel se posent lorsqu'il s'agit d'optimiser des ressources selon leurs valeurs avec une limite en capacité globale (typiquement dans l'ordonnancement des tâches).

3.2 Le problème du Sudoku

Dans le problème du sudoku, il s'agit de remplir les cases par un chiffre allant de 1 à 9, en respectant plusieurs contraintes : un seul même chiffre par ligne, par colonne et par bloc de 3x3. A priori, il n'y a pas de maximisation ni de minimisation : si l'on respecte les contraintes, toutes les lignes, toutes les colonnes et tous les blocs de 3x3 font la même somme. Il n'y a pas non plus de données externes, si l'on indique les cases sous forme de contraintes. Donc tout le problème réside dans la formulation des contraintes.

Dans notre fichier de programmation linéaire au format GMPL, que l'on nomme sudoku.mod, on commence par définir une variable binaire xijk, qui prend la valeur 1 si la case de ligne i allant de 1 à 9 et de colonne j de 1 à 9 également contient le chiffre k qui va de 1 à 9.

var x{i in 1..9, j in 1..9, k in 1.. 9 } binary;

On définit ensuite les contraintes.

subject to

Un seul chiffre par case, pour éviter d'avoir plusieurs chiffres k pour une même case ij : la somme des valeurs de xijk doit être égale à 1, et ce pour toutes les cases ij :

valeurs {i in 1..9, j in 1.. 9}: sum {k in 1..9} x[i,j,k] = 1;

Pour qu'il n'y ait qu'un même chiffre k sur toute la ligne i, on ajoute une contrainte qui postule que la somme du nombre de fois qu'un chiffre k donné apparaît sur la ligne i est égale à 1. Cela se traduit par :

lignes {i in 1..9, k in 1.. 9}: sum {j in 1..9} x[i,j,k] = 1;

On fait de même pour les colonnes :

colonnes {j in 1..9, k in 1.. 9}: sum {i in 1..9} x[i,j,k] = 1;

Pour éviter que l'on ait plusieurs fois le même chiffre dans un bloc 3x3, on ajoute une contrainte de bloc. Le nombre de fois que le chiffre k apparaît dans un bloc 3x3 doit être égal à 1, ce qui nécessite de parcourir les lignes i et les colonnes j par 3. Pour ce faire, on utilise deux petites variables xx et yy qui prennent successivement les valeurs 0, 3 et 6, ainsi on pourra parcourir les lignes i de xx+1 à xx+3 et les colonnes j de yy+1 à yy+3, pour couvrir les 9 blocs 3x3 du sudoku. Ce qui nous donne :

blocs {xx in {0, 3, 6}, yy in {0, 3, 6}, k in 1..9}:

    sum {i in (xx+1)..(xx+3), j in (yy+1)..(yy+3)} x[i,j,k] = 1;

Il nous reste à indiquer les contraintes des cases données, en forçant xijk pour un case ij. Et après calcul, on affiche le résultat sous forme de grille en s'aidant de boucles for et de fonctions printf (GMLP). On obtient au final le fichier sudoku.mod suivant :

#variable xijk=1 si case ij=k

var x{i in 1..9, j in 1..9, k in 1.. 9 } binary;

#contraintes

subject to

#un seul chiffre par case

valeurs {i in 1..9, j in 1.. 9}: sum {k in 1..9} x[i,j,k] = 1;

#un meme chiffre par ligne

lignes {i in 1..9, k in 1.. 9}: sum {j in 1..9} x[i,j,k] = 1;

#un meme chiffre par colonne

colonnes {j in 1..9, k in 1.. 9}: sum {i in 1..9} x[i,j,k] = 1;

#un meme chiffre par blocs de 3x3

blocs {xx in {0, 3, 6}, yy in {0, 3, 6}, k in 1..9}:

    sum {i in (xx+1)..(xx+3), j in (yy+1)..(yy+3)} x[i,j,k] = 1;

#cases donnees

c1: x[1,2,7]=1; c2: x[1,9,6]=1;

c3: x[2,5,6]=1; c4: x[2,6,9]=1; c5: x[2,8,3]=1;

c6: x[3,1,5]=1; c7: x[3,5,7]=1; c8: x[3,8,4]=1; c9: x[3,9,8]=1;

c10: x[4,1,4]=1; c11: x[4,5,9]=1; c12: x[4,7,7]=1;

c13: x[5,8,6]=1;

c14: x[6,1,8]=1; c15: x[6,2,6]=1; c16: x[6,4,1]=1; c17: x[6,5,3]=1;

c18: x[7,1,6]=1; c19: x[7,3,8]=1; c20: x[7,8,5]=1; c21: x[7,9,3]=1;

c22: x[8,4,9]=1; c23: x[8,7,8]=1;

c24: x[9,1,9]=1; c25: x[9,2,4]=1; c26: x[9,5,1]=1; c27: x[9,6,3]=1; c28: x[9,7,6]=1; c29: x[9,8,7]=1; c30: x[9,9,2]=1;

#calcul

solve;

#affichage

printf "\n\nSolution : \n";

printf "------------------------------------\n";

for {i in 1.. 9} {

    for {j in 1..9} {

         for {k in 1..9 :x[i,j,k] ==1} printf " %d ", k;

         printf "|";   

    }

    printf "\n";

    printf "------------------------------------\n";

}

end;

Je vous laisse découvrir la solution.

Dans ce genre de problème, l'essentiel est de bien vérifier ses contraintes. En effet, la moindre erreur de frappe dans les contraintes de cases données, par exemple en mettant deux fois le même chiffre sur la même ligne, et plus aucune solution n'est possible.

Conclusion

Notre escapade dans le monde de l'optimisation combinatoire, de la recherche opérationnelle et de la programmation linéaire se termine. Nous aurons vu avec GLPK un solveur de PL performant, fiable et open source, capable de résoudre un problème NP-Complet en deux temps trois mouvements à condition d'avoir correctement formulé et modélisé le problème. Dans la pratique, la programmation linéaire est utile pour la résolution des problèmes complexes [5].

Les problèmes NP-Complet se rencontrent fréquemment en informatique. Ils vont de l'ordonnancement des jobs du système d'exploitation jusqu'à l'optimisation de l'espace disque, en passant par de multiples applications (jeux de stratégie, planificateurs, gestion des ressources, ...).

Au niveau de la documentation, sur GLPK, vous trouverez des fichiers au format pdf (graphs.pdf, gmpl.pdf, glpk.pdf) dans le répertoire d'installation des documents (/usr/share/doc/glpk-doc/). Concernant la programmation linéaire, un nombre important de publications sur le sujet existent dans les domaines de l'optimisation combinatoire, de la recherche opérationnelle et des mathématiques appliquées. Un autre logiciel de programmation linéaire de référence, LP-SOLVE, est également disponible librement sous Linux.

Références

[1] E. J. Anderson, P. Nash – Linear Programming in infinite dimensional space – Ed. John Wiley & Sons – 1987

[2] H. S. Wilf – Algorithmes et complexité – Ed. Masson – 1989

[3] Michiels, Aarts, Korts – Theoretical Aspects of Local Search – Ed. Springer – 2007

[4] K. GhediraOptimisation combinatoire par métaheuristiques – Ed. Technip – 2007

[5] R. Faure, B. Lemaire, C. Picouleau – Précis de recherche opérationnelle – Ed. Dunod – 2000