introduction_metaheuristiques_Taillard
introduction_metaheuristiques_Taillard
introduction_metaheuristiques_Taillard
métaheuristiques
17
20
c É.D. Taillard 25 septembre 2017
e
br
em
pt
se
25
g
Al
T_
T_ ii
Al
g
25
se
pt
em
br
e
20
17
Introduction
17
Si le lecteur sait ce qu’est une heuristique — un algorithme permettant de
20
trouver avec un effort raisonnable une solution pas forcément la meilleure à un
problème donné, algorithme se basant sur des connaissances acquises par
l’expérience — la définition d’une métaheuristique est plus discutée.
Beaucoup d’auteurs ont adopté celle d’un processus maître itératif guidant
e
et modifiant les opérations d’une heuristique subordonnée. En conséquence,
les ouvrages dans le domaine sont décomposés en chapitres présentant cha-
br
cun une métaheuristique, comme par exemple le recuit simulé, la recherche
avec tabous, les colonies de fourmi artificielles ou les algorithmes génétiques,
em
pour ne citer que les plus connues.
Le présent livre aborde les métaheuristiques sous un angle différent. Il les
présente la forme d’un ensemble de méthodes de base qui sont combinées
entre elles pour constituer une heuristique. Ces principes de base sont suf-
pt
ment compromise.
La décomposition en sous-problèmes Lorsque l’on doit résoudre un pro-
blème de grande taille, il est nécessaire de le décomposer en sous-
g
problèmes plus petits. Ces derniers peuvent être eux-mêmes des pro-
blèmes difficiles, qu’il faut approcher par une technique appropriée, par
Al
Heuristique Métaheuristique
Domaine Un problème générique Optimisation combinatoire
d’application d’optimisation
Connaissance Spécifiques au problème Méthodes heuristiques
à inclure d’optimisation
Données à Les valeurs numériques Un problème générique
fournir d’un exemple de problème d’optimisation
Résultat Solution heuristique de Un algorithme heuristique
l’exemple de problème
17
duction au monde discret des méthodes basées sur les gradients pour
l’optimisation différentiable.
Randomisation et apprentissage Finalement, la répétition de constructions
20
ou de modifications permet d’améliorer la qualité des solutions produites,
pour peu que l’on fasse intervenir une composante aléatoire et/ou un
processus d’apprentissage.
e
Les métaheuristiques sont devenues un outil incontournable lorsqu’il s’agit
de s’attaquer à la résolution de problèmes d’optimisation difficiles, même si
br
elles ont parfois été décriées, surtout dans les années 1980 par des personnes
qui opposaient méthodes « exactes » et heuristiques. Depuis, on s’est rendu
compte que bien des méthodes exactes, donc sensées fournir une solution la
em
meilleure possible, ne garantissaient l’optimalité qu’avec une précision limitée
et qu’elles utilisaient abondamment des heuristiques !
pt
Public-cible
se
ambiguïté sur les inévitables interprétations qui peuvent être faites lors de la
présentation d’une trame de métaheuristique.
Comme un code source n’est utile que lorsque l’on veut connaître tous les
détails d’une méthode mais qu’il est d’une lecture peu attractive, nous avons
dépouillés les codes proposés au maximum, en essayant de les limiter à une
seule page. Certes, ces codes ne sont pas des « chevaux de course », on en
trouvera de plus complets et de plus efficaces par exemple dans la « Source
Code Library », mais ils ne contiennent que la quintessence des méthodes et
leur extrême concision devrait permettre au lecteur de les comprendre. Ainsi,
plus d’une douzaine de méthodes différentes ont été implantées tout en pre-
nant ensemble moins du dixième du nombre de lignes de code d’une des im-
3
plantations les plus rapide. Nous avons cependant dû, pour respecter cette
concision, compacter quelque peu les codes notamment en ne les commentant
que succinctement et parfois en plaçant plusieurs instructions sur une même
ligne ; nous prions donc le lecteur habitué à des codes plus aérés d’être indul-
gent.
Le choix des exercices a toujours été guidé par les solutions auxquelles on
s’attend, qui doivent être aussi univoques que possible. En effet, lorsque l’on
conçoit une nouvelle heuristique, il n’y a plus de solution juste ou fausse. Il y
a seulement des heuristiques qui fonctionnent bien pour certaines données du
17
problème considéré et d’autres qui ne produisent pas de résultats satisfaisants
— piètre qualité des solutions, temps de calcul prohibitif, etc.
Rien n’est plus déstabilisant pour un étudiant que des réponses ambiguës,
20
voire philosophiques. C’est pour cela les solutions pour l’ensemble des exer-
cices sont données en annexe.
Organisation du livre
e
br
Le présent ouvrage est découpé en 3 parties. La première (chapitres 1 et 2)
rappelle quelques bases de programmation linaire, de la théorie des graphes,
em
de théorie de la complexité et présente quelques problèmes, simples et diffi-
ciles, d’optimisation combinatoire. Le but de cette première partie est de rendre
le domaine intelligible par un lecteur ne disposant pas de connaissances parti-
culières en optimisation combinatoire.
pt
tives assemblant ces ingrédients et les utilisant dans des processus itératifs
fonctionnant sans mémoire sont y sont également incorporées.
La troisième partie présente des métaheuristiques plus évoluées qui uti-
25
un problème difficile qui ne diffère que par un petit détail pouvant paraître
anodin à première vue. Il met également en évidence le fait qu’un pro-
blème peut avoir diverses modélisations, même si sa formulation reste
très simple.
3. Méthodes constructives Ce chapitre présente des méthodes de construc-
tion de solutions, en commençant par un rappel sur les méthodes de
séparation et évaluation, largement utilisées pour la conception d’algo-
rithmes exacts. Ensuite deux méthodes élémentaires sont présentées, la
construction aléatoire et la construction gloutonne. Cette dernière sélec-
17
tionne séquentiellement les éléments à ajouter à une solution partielle, en
ne remettant jamais en cause les choix qui ont été faits. Cette méthode
peut être améliorée en évaluant plus profondément les conséquences du
20
choix d’un éléments. La recherche en faisceaux et la méthode pilote en
font partie.
4. Recherches locales Les méthodes d’amélioration d’une solution forment
l’épine dorsale de la plupart des métaheuristiques. Ces méthodes sont
e
basées sur la définition d’un ensemble de solutions voisines, pour toute
br
solution du problème. La définition de cet ensemble dépend naturelle-
ment de la modélisation du problème. Selon cette dernière, un voisinage
s’exprimant naturellement peut s’avérer trop petit pour mener à des solu-
em
tions de qualité ou au contraire trop grand, induisant des temps de calcul
prohibitifs. Diverses méthodes ont été proposées pour agrandir le voisi-
nage, comme la recherche en éventail ou les chaînes d’éjections, ou pour
le réduire, comme la recherche granulaire ou la liste de candidats.
pt
nages.
Al
17
tions, on peut chercher à apprendre comment les combiner et comment
gérer cette population. La méthode la plus populaire dans ce domaine est
sans conteste les algorithmes génétiques. Ces derniers sont cependant
20
une métaheuristique moins évoluée que la recherche par dispersion. Fi-
nalement, parmi les dernières-nées des métaheuristiques on trouve les
méthodes à essaim particulaires, qui semblent adaptés à l’optimisation
continue.
e
br
em
pt
se
25
g
Al
T_
T_ 6
Al
g
25
se
pt
em
br
e
20
17
17
20
Première partie
e
br
Concepts de base
em
pt
se
25
g
Al
T_
T_
Al
g
25
se
pt
em
br
e
20
17
Chapitre 1
17
Éléments de théories des
20
graphes et de la complexité
e
1.1 Optimisation combinatoire
br
em
Le domaine d’application privilégié des métaheuristiques est l’optimisation
combinatoire. Pour fixer les idées, il n’est donc pas inutile de brièvement rap-
peler en quoi consiste ce domaine. Prenons un petit exemple de problème
pt
qui ont une frontière commune ne reçoivent pas la même couleur. Dans la fi-
gure 1.1 on a utilisé 5 couleurs différentes, sans se préoccuper de l’attribution
politique des îles ou enclaves.
Ce problème est combinatoire. S’il y a n divisions à colorer et que l’on désire
25
le petit exemple suivant : peut-on colorer la Suisse (s) et les pays voisins
(Allemagne d, France f , Italie i, Autriche a et Liechtenstein l) avec 3 cou-
leurs (voir figure 1.2) ? Pour cela, on peut introduire 18 variables booléennes
bs1 , bs2 , bs3 , bd1 , bd2 , . . . , bl3 qui indiqueront pour chaque pays la couleur 1, 2 ou
3 qu’il reçoit. Il s’agit d’écrire une longue formule booléenne qui prendra la va-
leur vraie si, et seulement s’il existe une coloration admissible en 3 couleurs.
Tout d’abord, on peut imposer que la Suisse soit colorée dans au moins une
couleur : bs1 ∨ bs2 ∨ bs3 . Mais il ne faut pas qu’elle reçoive à la fois la couleur 1
et la couleur 2, ce qui peut s’écrire : bs1 ∧ bs2 , ce qui est équivalent à : bs1 ∨ bs2 .
Elle ne devra également pas prendre à la fois la couleur 1 et la couleur 3 ou la
couleur 2 et la couleur 3. Ainsi, pour imposer que la Suisse soit colorée avec
10 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
17
20
e
F IGURE 1.1 – Une ancienne carte de l’Europe colorée en 5 couleurs (en comptant celle
du fond).
br
em
exactement 1 couleur, on aura la conjonction de 4 clauses :
Pour chacun des pays concernés, il faudra également écrire une conjonction de
4 clauses, mais avec les variables correspondant aux autres pays. Finalement,
se
pour chaque frontière, il faudra imposer que les couleurs de part et d’autre
soient différentes. Par exemple, pour la frontière entre la Suisse et la France,
on devra avoir :
(bs1 ∨ bf 1 ) ∧ (bs2 ∨ bf 2 ) ∧ (bs3 ∨ bf 3 )
25
De manière générale, si on désire colorer une carte avec n pays qui ont,
en tout, m frontières communes en k couleurs, on peut se poser la question
du nombre de variables que le problème aura ainsi que du nombre de clauses
qui devront être toutes vraies simultanément. Tout d’abord, il faudra introduire
g
k·(k−1)
2 clauses pour s’assurer que chaque pays reçoive exactement 1 couleur.
Finalement, pour chaque frontière, il faudra écrire m·k clauses. Le problème de
la coloration d’une telle carte en k couleurs aura une solution si, et seulement
T_
s’il existe des valeurs vrai ou f aux pour chacune des variables telles que
l’ensemble de toutes ces clauses soient vraies simultanément. Ce problème
est connu comme celui de l’acceptabilité et joue un rôle central en théorie de
la complexité. Tout ce développement a permis de formaliser le problème par
une seule, longue, formule booléenne, mais ne nous dit pas encore comment
trouver une solution à ce problème !
Un algorithme, très primitif, permettant de trouver une solution est d’es-
sayer toutes les valeurs possibles pour les variables (il y a 2nk jeux de valeurs
différentes), et pour chaque jeu, de vérifier si la formule est juste.
Tel que modélisé par une formule booléenne, le problème de la coloration
d’une carte est un problème de décision, dont la réponse est vrai (on arrive à
1.1 O PTIMISATION COMBINATOIRE 11
colorer la carte en k couleurs) ou f aux (on n’y arrive pas). En supposant que
l’on dispose d’un algorithme A qui permette de trouver les valeurs à donner aux
variables booléennes pour que la formule soit vraie — ou qui dise que de telles
valeurs n’existent pas — est-il possible de résoudre le problème d’optimisation,
c’est-à-dire de trouver la valeur minimale de k telle qu’il existe une coloration
admissible ?
Une façon de répondre à cette question est de constater que l’on a be-
soin, au maximum, de n couleurs pour colorer n surfaces, en attribuant une
couleur différente à chacune. Par conséquent, on sait qu’une coloration en n
17
couleurs existe. On peut appliquer l’algorithme A en lui demandant de trouver
une coloration en n − 1 couleurs, puis n − 2, etc. jusqu’à ce qu’il réponde qu’il
n’existe pas de coloration. On saura alors que la dernière valeur pour laquelle
20
l’algorithme a trouvé une solution correspond à une coloration optimale.
Une technique plus rapide est de procéder par dichotomie : plutôt que de
diminuer d’une unité à chaque appel de A, on maintient deux valeurs kmin
et kmax pour lesquelles on sait, respectivement, qu’il n’y a pas de coloration
e
admissible et qu’il en existe une. En éliminant le cas de la carte triviale qui n’a
pas de frontière, on sait que l’on peut démarrer avec kmin = 1 et kmax = n. On
réponse
c couleurs.
est non,
em
kmin +kmax
on modifie kmin ← b 2 c. On répète le procédé jusqu’à ce que kmax =
kmin + 1, valeur qui correspond au nombre optimum de couleurs.
fonction linéaire et que l’ensemble peut être décrit au moyen de fonctions li-
néaires, on parle de programmation linéaire.
Un programme linéaire sous forme canonique s’écrit mathématiquement de
la manière suivante :
g
Maximiser z = c1 x1 + c2 x2 + . . .+ cn xn (1.1)
Al
contraintes ...
am1 x1 + am2 x2 + . . .+ amn xn 6 bm
xi > 0(i = 1, . . . , n) (1.3)
Dans le cas particulier où tous les coefficients cj , aij et bi sont positifs, une
interprétation du problème consiste à chercher des quantités xj à produire,
non négatives (1.3), chaque unité du produit j rapportant un revenu cj , de
sorte à maximiser le revenu global z (1.1). La production d’une unité du produit
j consomme des quantités aij de chaque ressource i, disponible en quantité
bi (contraintes 1.2).
12 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
17
indiquent si la couleur k (pour k = 1, . . . , n) est utilisée (yk = 1) ou non (yk = 0)
ainsi que des variables xik qui indiquent si la surface i reçoit la couleur k. Le
programme linéaire en nombres entiers suivant permet de formaliser le pro-
20
blème de la coloration d’une carte :
n
X
Minimiser yk (1.4)
e
k=1
Sous contraintes
Xn
xik =
k=1 br
1 i = 1, . . . , n (1.5)
em
xik − yk 6 1 i, k = 1, . . . , n (1.6)
xik + xjk 6 1 (i, j) ont une frontière commune, (1.7)
k = 1, . . . , n
pt
mier ensemble de contraintes (1.5) impose que chaque sommet reçoive une
couleur exactement ; le second ensemble (1.6) que l’on n’attribue pas à un
sommet une couleur non utilisée ; et l’ensemble (1.7) empêche de donner
la même couleur à deux surface contigües. Les contraintes d’intégrité (1.8)
25
peuvent également s’écrire en faisant ressortir des inégalités linéaires (yk >
0, yk 6 1, yk ∈ Z).
La programmation linéaire est un outil de modélisation et de formalisa-
tion de problèmes très puissant. Lorsqu’il n’y a pas de contraintes d’intégrité
g
riables et des milliers de contraintes. Dans ce cas, la résolution est à peine plus
complexe que la résolution d’un système d’équations linéaires. La principale li-
mitation vient essentiellement de la place mémoire nécessaire au stockage des
T_
17
20
F IGURE 1.2 – La Suisse et ses pays limitrophes que l’on désire colorer. Chaque pays est
symbolisé par un disque et un frontière communes est symbolisée par un trait reliant
e
entre les pays-disques correspondants. On obtient ainsi une représentation schéma-
tique du problème de coloration sous forme d’un graphe.
br
présenter d’autres modélisations du problème de coloration (voir section 2.9),
em
nous rappellerons très brièvement quelques définitions de la théorie des gra-
phes sur lesquelles le reste de l’ouvrage s’appuiera.
pt
a des éléments qui présentent des relations entre eux. On représente les élé-
ments par un point et si deux éléments sont en relation, on les relie par un seg-
ment. Ainsi, le problème de coloration d’une carte vu précédemment pourra se
25
Graphe non orienté, sommet, arête Un graphe non orienté G est une col-
lection d’éléments V appelés également sommets (« vertex » en globish) ou
g
G = (V, E). Si v et w sont deux sommets reliés par une arête, on dit qu’ils
sont adjacents ou encore que cette arête est incidente avec v et w. Lorsque
T_
Graphe orienté, arc Dans certains cas, les relations entre les couples d’élé-
ments sont ordonnées. On parle alors de graphe orienté et les « arêtes orien-
tées » sont appelées des arcs. Il convient dès lors de distinguer les deux extré-
mités d’un arc en une extrémité initiale et une extrémité terminale ou finale. Un
arc dont l’extrémité initiale et terminale est confondue est également appelé
14 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
boucle, comme pour le cas non-orienté, de même que deux arcs qui ont les
mêmes extrémités initiales et terminales sont des arcs parallèles ou multiples.
Deux arcs a et b dont l’extrémité initiale de l’un est l’extrémité terminale de
l’autre et vice versa sont dit tête bêche ; dans ce cas, les arcs a et b ne forment
pas des arcs multiples.
Degré Le degré d’un sommet v d’un graphe non orienté, noté deg(v) est le
nombre d’arêtes incidentes à v. Une boucle augmente le degré d’un sommet
de 2. Un sommet de degré 1 est dit pendant. Un graphe est régulier si tous
17
ses sommets ont le même degré. Pour un graphe orienté, on parle de degré
sortant ou degré extérieur, noté deg + (v) le nombre d’arcs sortant du sommet
v. Le degré entrant ou degré intérieur, noté deg − (v) est le nombre d’arcs dont
20
l’extrémité finale est v.
Chaîne, chemin, cycle, circuit Une chaîne est une suite alternée de som-
mets et d’arêtes, débutant et finissant par un sommet, telle que chaque arête
e
est encadrée par ses deux extrémités. Un cycle est une chaîne dont les deux
br
extrémités sont identiques. Pour un graphe orienté, l’équivalent d’une « chaîne
orientée » est appelé chemin. Dans ce cas, chaque arc est encadré d’abord
par son extrémité initiale puis par son extrémité terminale. Un « cycle orienté »
em
est appelé un circuit.
La longueur d’une chaîne (ou d’un chemin) est le nombre d’arêtes (d’arcs)
qui la compose. Une chaîne (ou un chemin) est simple si une arête (un arc) n’y
apparaît pas plus d’une fois. Une chaîne (ou un chemin) est élémentaire si un
pt
sommet n’y apparaît pas plus d’une fois (à l’exception des extrémités).
se
sommets. Un graphe orienté est fortement connexe s’il existe un chemin entre
tout couple de sommet (de l’un vers l’autre et vice versa).
cycle. Les sommets pendants d’un arbre sont également appelés feuilles. Un
graphe sans cycle est une forêt (qui peut donc être composée de plusieurs
Al
arbres). Un « arbre orienté » tel qu’il existe un chemin d’un sommet r à tous les
autres est une arborescence de racine r.
T_
cycle qui passe par toute les arêtes ou tous les sommets exactement une fois
(à l’exception du sommet de départ et d’arrivée, qui est confondu).
Graphes particuliers Un graphe simple non orienté est complet s’il existe
une arête entre toute paire de sommet. Un graphe G = (V, E) tel que V =
V1 ∪ V2 , V1 ∩ V2 = ∅ et toute arête de E a l’une de ses extrémités dans V1 et
l’autre dans V2 est biparti. Un graphe G = (V, E) tel que V = V1 ∪ V2 ∪ · · · ∪ Vm ,
Vi ∩ Vj = ∅ ∀i, j ∈ V et que pour tout arête de E, il existe un i tel que l’arête a
l’une de ses extrémités dans Vi et l’autre dans Vi+1 est un graphe en couches.
17
Une clique est un sous-ensemble des sommets d’un graphe qui induit un sous-
graphe complet. Un ensemble stable est un sous-ensemble de sommets non
adjacents deux à deux.
20
Coloration, couplage Le problème de la coloration des sommets d’un gra-
phe qui a été utilisé comme exemple introductif à la section 1.1 sur l’optimi-
e
sation combinatoire permet de définir le nombre chromatique d’un graphe G,
parfois noté χ(G), qui est le nombre minimum de couleurs nécessaires pour
br
colorer les sommets de G sans que deux sommets adjacents ne reçoivent la
même couleur. Le nombre de stabilité d’un graphe G, parfois noté α(G) est le
em
cardinal maximum d’un ensemble stable de G. L’indice chromatique d’un gra-
phe G, parfois noté q(G) est le nombre minimum de couleurs nécessaires pour
colorer les arêtes de G de sorte que deux arêtes incidentes ne reçoivent pas
la même couleur.
pt
d’un graphe est un transversal s’il recouvre toutes les arêtes du graphe, c’est-
à-dire si chaque arête est incidente à au moins un sommet du transversal.
25
Réseau
Dans bien des cas, on veut associer une valeur à chaque arête ou à chaque
arc. La valeur w(e) ou poids de l’arête ou de l’arc e représentera typiquement
g
Flots Pour résoudre des problèmes de flot, on associe une limitation de débit
aux arcs d’un réseau, ou capacité wij > 0 entre les nœuds i et j. Un flot
d’un sommet-source s à un sommets-puits t dans un réseau R = (V, E, w) est
l’attribution d’un flux, ou quantité de flot xij telle que la somme des flux sortant
d’un
P sommet soitPégale à la somme des flux entrants (loi de conservation) :
i∈Succ(j) xij = k∈P red(i) xki ∀i ∈ V, i 6= s, t. La valeur du flot de s à t est
la somme des flux qui sortent de s (qui est identique à celle qui entre en t, en
raison de la loi de conservation). Lorsque 0 6 xij 6 wij ∀(i, j) ∈ E, le flot est
dit compatible ou admissible.
La séparation des sommets d’un réseau R = (V, E, w) en A ⊂ V et son
complémentaire A ⊂ V permet de définir une coupe ⊆ E de A vers A, qui est
16 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
l’ensemble des arcs qui ont leur extrémité initiale dans A et terminale dans A.
La capacité d’une coupe est la somme des poids des arcs de la coupe.
Les flots permettent de modéliser des problèmes qui n’ont, à première vue,
rien en commun avec des débits, comme par exemple les problèmes d’affec-
tation de ressources (voir par exemple la section 2.4.2). Plus loin dans ce cha-
pitre, nous passerons en revue quelques algorithmes bien connus et efficaces
pour la résolution de problèmes, comme la recherche d’un arbre, d’un chemin
ou d’un flot optimum dans un réseau. D’autres problèmes, comme celui de la
coloration de graphes, sont difficiles, et les seuls algorithmes connus pour les
17
résoudre prennent un temps qui peut croître exponentiellement avec la taille
du graphe.
La classification en problèmes faciles et difficiles ressort de la théorie de la
20
complexité. Les métaheuristiques ont été conçues pour trouver des solutions
acceptables avec un effort de calcul limité pour des problèmes en principe
difficiles. Avant d’élaborer un nouvel algorithme sur la base des principes des
métaheuristiques, il est essentiel de s’assurer que le problème traité est bien
e
difficile et qu’il n’existe pas déjà un algorithme efficace pour le résoudre. La
suite de ce chapitre expose quelques bases théoriques dans le domaine de la
br
classification des problèmes selon leur difficulté.
em
1.2 Éléments de théorie de la complexité
La théorie de la complexité a pour but de classer les problèmes afin de
pt
La raison est très simple : on conçoit parfaitement que l’on ait besoin de
plus d’opérations pour traiter un plus gros volume de données, ce qui élimine
les fonctions non croissantes, comme les fonctions trigonométriques. Si l’on
se limite aux méthodes séquentielles, on conçoit aisément que l’on doive lire
g
chaque donnée au moins une fois, ce qui implique une croissance au moins
Al
17
nombre exponentiel) et en isolant le plus court, mais à l’aide d’algorithmes
faisant appel à des propriétés mathématiques du plus court chemin. Ces al-
gorithmes effectuent un nombre d’opérations polynomial en la taille du réseau.
20
Le problème de la recherche d’un plus court chemin est donc un problème
simple. Par contre, la recherche d’un plus long (ou d’un plus court) chemin
simple ou élémentaire (sans boucles ou sans passer deux fois par le même
sommet) entre deux sommets est un problème difficile, car on ne connaît pas
e
d’algorithme polynomial pour le résoudre.
Pour terminer, il faut mentionner que la classe des polynômes présente
br
une propriété mathématique intéressante : elle est fermée. La composition de
deux polynômes est encore un polynôme. Pratiquement, dans le contexte de
em
la programmation, cela signifie que si à l’intérieur d’un algorithme, on utilise un
nombre polynomial de fois un sous-algorithme, lui-même polynomial, le nom-
bre total d’instructions effectuées par tous les appels du sous-algorithme est
encore polynomial.
pt
cet algorithme n’est pas vraiment élémentaire ; son analyse et son élaboration
constituent même des tâches aux limites des capacités humaines. Donc, tester
si un nombre est premier est un problème simple (car il existe un algorithme
polynomial pour le résoudre), mais l’algorithme en question est toutefois dif-
ficile à implanter et requerrait un temps de calcul prohibitif pour montrer que
243112609 − 1 est premier. Inversement, il peut exister des algorithmes qui en
théorie pourraient dégénérer mais qui se comportent tout-à-fait bien en pra-
tique, comme l’algorithme du simplexe en programmation linéaire.
Les ressources utilisées pour durant l’exécution d’un algorithme sont limi-
17
tées. Elles sont de plusieurs types : nombre de processeurs, place mémoire,
temps. En se restreignant à cette dernière ressource, si l’on désire mesurer
l’efficacité d’un algorithme, on pourrait évaluer sa durée d’exécution sur une
20
machine donnée. Malheureusement, cette mesure présente de nombreuses
faiblesses. Tout d’abord, elle est relative à une machine particulière, dont la
durée de vie est limitée à quelques années. Ensuite, la manière dont l’algo-
rithme a été implanté (langage de programmation, compilateur, options, sys-
e
tème d’exploitation) peut influencer énormément son temps d’exécution. Par
conséquent, on préfère mesurer le nombre d’opérations caractéristiques qu’un
br
algorithme va effectuer. En effet, ce nombre ne dépend pas de la machine ou
du langage et peut parfaitement s’évaluer théoriquement.
em
On appelle complexité d’un algorithme une fonction f (n) qui donne le nom-
bre d’opérations caractéristiques exécutées dans le pire des cas lorsqu’il s’exé-
cute sur un problème dont le volume des données est n. Il faut mentionner que
cette complexité n’a rien à voir avec la longueur du code ou la difficulté de
pt
luation du pire des cas est très importante dans les applications où le temps
de réponse est critique.
25
g(n) = 0, 2 · n3 opérations. Il est clair que pour un jeu de données avec n = 10,
A1 effectue 5 fois plus d’opérations que A2 . Par contre, dès que n > 50, c’est
T_
A2 qui en fera le plus. En réalité, quels que soient les cœfficients positifs de
n2 et n3 dans f (n) et g(n), la fonction g sera plus grande que f à partir d’un
certain n. On appelle ordre la croissance asymptotique d’une fonction.
17
des cas au moins autant d’opérations que A n’en effectue dans le pire des cas)
ou qu’un algorithme C est optimal (e.g. C résout un problème en effectuant dans
le pire des cas un nombre d’opérations pas plus grand que le minimum qui est
20
requis par n’importe quel algorithme pour résoudre ce problème).
Si le meilleur et le pire des cas sont les mêmes, c’est-à-dire si ∃0 < c1 < c2
tel que c1 · g(n) 6 f (n) 6 c2 · g(n) alors on notef ∈ Θ(g). Il faut bien différencier
la notation Θ(·) d’une notion (souvent mal définie) de complexité moyenne : en
e
effet, si l’on prend l’algorithme de tri rapide de n éléments, on peut dire qu’il
sera en Ω(n) et en O(n2 ), mais pas en Θ(n · log(n)), même si en moyenne il
prend un temps proportionnel à n · log(n).
br
Cependant, on peut montrer que l’espérance du temps de calcul de l’algo-
em
rithme de tri rapide pour un ensemble de n éléments aléatoirement mélangés
est proportionnelle à (n · log(n)). Les notations O(·) (espérance ou moyenne
théorique) et Ô(·) (moyenne empirique) seront utilisées plus loin, tout en souli-
gnant qu’elles ne sont pas couramment utilisées car on doit spécifier sur quel
pt
avantages :
Al
complexe.
— 25n3 = O(3n3 ) et 3n3 = O(25n3 ), ce qui signifie que deux fonctions
différant uniquement par un facteur sont du même ordre, ce qui permet
de s’affranchir des vitesses relatives des calculateurs. Plutôt que d’écrire
O(25n3 ), on écrira donc O(n3 ), ce qui est équivalent et plus simple.
— 3n3 + 55n = O(n3 ), ce qui implique que l’on peut négliger les termes
d’ordre inférieur et ne retenir que la puissance la plus élevée.
Il est important d’insister sur le fait que la complexité d’un algorithme est
une notion théorique, qui se dérive par réflexion et calculs qui peuvent être
faits avec une feuille et un crayon. Cette complexité s’exprime généralement
20 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
αΩ
O(n!)
O(2n )
O(n3 )
S˚ O(n2 )
O(n log n)
h O(n)
17
s
O(log n)
20
µs O(1)
e
ns 1 2 3 4 5 6 7 8 9 10 11 12 13 14
10 10 10 10 10 10 10 10 10 10 10 10 10 10 1015
br
F IGURE 1.3 – Illustration de la croissance de quelques fonctions souvent utilisées pour
em
exprimer la complexité d’un algorithme. L’axe horizontal donne la taille du problème
(avec une croissance exponentielle) et l’axe vertical donne l’ordre de grandeur du temps
de calcul (avec une croissance bi-exponentielle).
pt
garithme, car O(loga (n)) = O(logb (n)) ), O(nc ) (puissance fractionnaire ; avec
0 < c < 1), O(n) (linéaire), O(n · log(n)) (log-linéaire ou quasi-linéaire), O(n2 )
(quadratique), O(n3 ) (cubique), O(nc ) (polynomial ; avec c > 1 constante),
O(nlog(n) ) (quasi-polynomial ; super-polynomial, sous-exponentiel), O(cn ) (ex-
g
ponentiel, avec c > 1 constante), O(n!) (factoriel). La figure 1.3 illustre la crois-
sance de quelques-unes de ces fonctions.
Al
17
Cependant, comme on veut se restreindre à un problème de décision, la ques-
tion posée ne sera pas « quelle est la tournée de longueur minimum », mais
20
Question : Existe-t-il un tour passant par toutes les villes de C et n’ayant
pas une longueur supérieure à P B ? (i.e. trouver σ, une permutation des élé-
n−1
ments 1, 2, ?, n telle que dσn ,1 + i=1 dσi , dσi+1 6 B. )
e
br
Schéma de codage, langage et machine de Turing déterministe
l’automate passe d’un état à un autre selon la lettre du mot en cours de lecture
et des transitions associées. Comme un tel automate a un nombre fini d’états,
cette machine a donc une mémoire bornée.
Un modèle un peu plus complexe est un automate à pile, fonctionnant de fa-
çon semblable à une machine à états finis mais disposant d’une pile. À chaque
étape on lit une lettre du mot ainsi que le symbole en sommet de pile (si elle
n’est pas vide) et on se met dans un nouvel état, en empilant au besoin un
nouveau symbole en sommet de pile. Ce type d’automate permet de faire des
calculs un peu plus complexes, par exemple reconnaître les mots d’un langage
17
non contextuel (e.g. effectuer l’analyse syntaxique d’un programme décrit par
une grammaire de type 2). Un modèle de calculateur plus puissant qu’un auto-
mate à pile est la machine de Turing.
20
Machine de Turing déterministe Pour pouvoir représenter mathématique-
ment le fonctionnement d’un ordinateur, Allan Turing a imaginé en 1936 une
e
machine fictive (il n’existait pas d’ordinateurs à cette époque) dont le fonction-
nement peut être modélisé par une fonction de transition. Avec cette machine,
br
il est possible de mettre en œuvre tous les algorithmes usuels ou de recon-
naître en un temps fini un mot généré par une grammaire générale de type 0.
La figure 1.4 illustre une telle machine, composée d’un programme qui com-
em
mande le défilement d’une bande magnétique ainsi qu’une tête de lecture et
d’écriture.
pt
Programme δ
se
État
25
Tête de lecture-écriture
Bande infinie
··· b b b b 1 1 0 0 0 0 1 0 1 1 0 0 1 0 1 1 0 b b ···
g
-2 -1 0 1 2 3 4
Al
Cases
F IGURE 1.4 – Schéma de principe d’une machine de Turing, qui permet de modéliser et
T_
de formaliser un calculateur
17
Soit q, l’état courant de la machine, σ le symbole lu sur la bande ainsi que
(q 0 , σ 0 , ∆) = δ(σ, q). Un pas de la machine consiste à :
— Effacer σ sur la bande et écrire σ 0 à la place.
20
— Déplacer la tête de lecture d’une case à gauche si ∆ = −1 et d’une case
à droite si ∆ = 1.
— Se mettre dans l’état q 0 .
e
La machine s’arrête lorsqu’on se trouve dans l’un des deux états qY ou qN .
C’est pour cette raison que la fonction de transition δ n’est définie que pour les
états non finaux de la machine.
br
Bien que très simple au niveau conceptuel, une machine de Turing permet
em
de représenter schématiquement tout ce qui se passe dans un ordinateur cou-
rant, contrairement à d’autres modélisations de machines plus simples, comme
les automates à états finis (qui n’ont qu’une tête de lecture se déplaçant tou-
jours dans le même sens) ou les automates à pile. Ces machines permettent
pt
Q = q0 , q1 , q2 , q3 , qY , qN
État 0 1 b
(q1 , b, −1)
Al
q0 (q0 , 0, 1) (q0 , 1, 1)
q1 (q2 , b, −1) (q3 , b, −1) (qN , b, −1)
q2 (qY , b, −1) (qN , b, −1) (qN , b, −1)
T_
17
TM (n) que la machine a besoin pour s’arrêter, quelle que soit la chaîne x de
longueur n écrite sur la bande initialement. On dit qu’un programme pour ma-
chine de Turing déterministe est en temps polynomial s’il existe un polynôme p
20
tel que TM (n) 6 p(n)
La classe P des langages comprend tous les langages L tels qu’il existe
un programme pour machine de Turing déterministe reconnaissant L en temps
polynomial. Par abus de langage, on dira que le problème π appartient à la
e
classe P si le langage associé à π et à un schéma de codage e (non spécifié
mais supposé raisonnable) appartient à P . Lorsque l’on utilise l’expression « il
br
existe un programme », on fait en quelque sorte l’hypothèse que l’on dispose du
meilleur programme possible pour ce problème, mais sans forcément connaître
em
ce programme. Inversement, si l’on connaît un programme — pas forcément le
meilleur — s’exécutant en temps polynomial pour ce problème, alors on peut
affirmer que le problème fait partie de la classe de complexité P .
pt
vivons. Conceptuellement, cette machine est composée d’un module qui de-
vine la solution du problème et l’écrit sur les cases de la bande d’indices né-
gatifs (voir la figure 1.5). Cet artifice permet de s’affranchir de notre ignorance
d’un algorithme efficace pour résoudre un problème : la machine le devine.
La spécification d’un programme pour machine de Turing non déterministe
est identique à celle d’une machine déterministe. Initialement la machine se
trouve dans l’état q0 , la bande contient la chaîne x codant le problème dans les
cases 1 à |x| de la bande et le programme est inactif. Suit alors une phase divi-
natoire durant laquelle le module écrit des symboles aléatoires dans les cases
négatives et s’arrête arbitrairement. Ensuite, le programme de la machine s’ac-
tive et elle fonctionne comme une machine de Turing déterministe.
1.2 É LÉMENTS DE THÉORIE DE LA COMPLEXITÉ 25
Programme δ
Module divinatoire
État
17
Bande infinie
··· b b 1 0 1 1 0 0 0 0 1 0 1 1 0 0 1 0 1 1 0 b b ···
-2 -1 0 1 2 3 4
20
Cases
F IGURE 1.5 – Schéma de principe d’une machine de Turing non déterministe, qui per-
met de de formaliser la classe N P mais n’existe pas dans la réalité
e
br
Avec une telle machine, il est clair qu’une chaîne donnée x peut engendrer
divers calculs, à cause du caractère non déterministe de la phase divinatoire.
em
Un problème ayant une solution peut engendrer des exécutions se terminant
par l’état qN de la machine, ou des exécutions de diverses durées se terminant
par l’état qY (mais la machine ne peut pas terminer dans l’état qY pour un
problème n’ayant pas de solution). Le langage LM reconnu par la machine M
pt
est par définition l’ensemble des chaînes x ∈ Σ∗ telles qu’il existe au moins un
calcul pour lequel la chaîne x est acceptée. Le temps de calcul TM (n) est le
se
Théorème : Si π est dans N P , il existe un polynôme p(n) tel que π peut être
T_
∝T (n)
Problème 1 Problème 2
Codage
f (Donnees1)
Données 1 Données 2
Résolution
Résultat 1 Résultat 2
17
Reconstruction
F IGURE 1.6 – Transformation en temps T (n) d’un problème 1 vers un problème 2. Dans
la pratique seule sont réalisées les opérations marquées par les flèches pleines.
20
0 0
on a |Γ|q (n) = 2q (n)·log2 (|Γ|) , d’où le résultat.
e
br
Transformation polynomiale : La notion de transformation polynomiale d’un
premier problème en un second problème est fondamentale en théorie de la
complexité, car elle est d’une grande utilité pour la classification des problèmes.
em
En effet, si l’on sait — ou, pour les problèmes difficiles : si l’on savait — ré-
soudre efficacement le second problème et que le premier peut se mettre sous
la forme du second par une transformation peu coûteuse en ressources de
calcul, alors on peut aussi résoudre efficacement le premier problème.
pt
que si, et seulement s’il existe un parcours hamiltonien dans le graphe initial.
On en déduit que le problème du parcours hamiltonien peut se transformer en
un problème de voyageur de commerce. Il faut noter que l’inverse n’est pas
forcément vrai.
Classe NP-Complet
17
peut se transformer polynomialement en π.
20
en constatant que la composition de deux polynômes est encore un polynôme,
on peut énoncer quelques théorèmes dont la démonstration est évidente.
— Si π est NP-complet et si π peut être résolus en temps polynomial, alors
e
P = NP.
— Si π est NP-complet et si π n’appartient pas à la classe P , alors P 6= N P .
br
— Si π1 se transforme polynomialement en π2 et π2 se transforme polyno-
mialement en π3 , alors π1 se transforme polynomialement en π3 .
em
— Si π1 est NP-complet, π2 appartient à la classe N P et π1 se transforme
polynomialement en π2 , alors π2 est NP-complet.
On ne connaît pas de problème NP-complet qui se résolve en temps po-
pt
Il n’est pas possible d’imaginer quels sont tous les problèmes de N P , et en-
core moins de trouver une transformation pour chacun d’eux vers le problème
universel. Or, un tel problème existe, et le premier pour lequel on a démontré
qu’il était N P -complet est le problème d’acceptabilité. À partir de ce résultat,
on a réussi à montrer que tout un ensemble de problèmes faisaient aussi par-
tie de la classe NP-complet, en utilisant le principe énoncé dans la remarque
ci-dessus.
fausse si, et seulement si tous ses littéraux sont faux. Un problème d’accepta-
bilité est un ensemble de clauses reliées entre elles par des « et » logiques. Un
exemple de problème d’acceptabilité aura pour réponse « oui » s’il existe une
affectation de valeurs aux variables booléennes telles que toutes les clauses
du problème soient simultanément vraies. Par exemple (u1 ∨ u2 ) ∧ (u1 ∨ u2 ) est
un exemple de problème acceptable. Par contre, (u1 ∨u3 )∧(u1 ∨u3 )∧(u1 )∧(u2 )
n’est pas un exemple acceptable.
17
ceptabilité est NP-complet, il faut tout d’abord vérifier que ce problème soit
dans N P , ce qui est évident. Ensuite, il faut montrer que tout problème de
N P peut se transformer polynomialement en un problème d’acceptabilité. Une
20
telle transformation ne peut se faire explicitement, car on ne peut pas prendre
individuellement tous les problèmes de N P et leur trouver une transformation
adéquate, premièrement parce que le nombre de problèmes de N P est infini et
deuxièmement parce que trouver une transformation directe pourrait être très
e
difficile. Au lieu de cela, on montre implicitement qu’une telle transformation
br
existe, cette transformation pouvant être inconnue.
Pour prouver que cette transformation existe, on part de la définition de la
classe N P . Si un problème π appartient à cette classe de complexité, c’est
em
qu’il existe un programme en temps polynomial pour machine de Turing non
déterministe reconnaissant π. Cette machine est spécifiée par son alphabet de
bande, son ensemble d’états et la fonction δ de transition. Si l’on veut montrer
que π peut se transformer polynomialement en un problème d’acceptabilité, il
pt
complet.
La construction de cette fonction f se fait par six groupes de clauses qui
ne seront vraies que pour une exécution du programme conforme aux spécifi-
cations de fonctionnement de la machine de Turing. Les variables booléennes
g
chine est dans l’état q), Htj (au temps t, la tête examine la case j) et Stjσ
(au temps t, la case j de la bande contient le symbole σ). Comme le temps
de calcul est limité par un polynôme p(n) en la taille n du problème, il y a un
T_
17
4. Au temps 0 (début de la phase de vérification), on est dans l’état q0 ,
20
la tête est sur la case 1, la bande contient dans les cases 1 à n les
symboles codant l’exemple de problème de π et il y a le symbole blanc
dans les autres cases non négatives (par contre, on ne sait pas ce qu’il
y a dans les cases négatives qui devraient contenir une solution avec
e
réponse « oui » au problème).
br
Q0q0 ∧ H01 ∧ S00b ∧ S01x1 ∧ S02x2 ∧ · · · ∧ S0nxn ∧ S0n+1b ∧ · · · ∧ S0,p(n)+1,b
em
5. À un temps borné par un polynôme en n, la machine doit se trouver dans
l’état qY .
Qp(n)qY
pt
se
la machine est dans l’état q et qu’elle lit le symbole σ, alors, au temps t+1,
la machine doit se trouver dans l’état q 0 , le symbole s0 doit être écrit dans
la case j et la tête se trouver en position j + ∆, avec (q 0 , s0 , ∆) = δ(s, q).
g
(Htj ∨Qtq ∨Stjσ ∨Ht+1,j+∆ )∧(Htj ∨Qtq ∨Stjσ ∨Qt+1q0 )∧(Htj ∨Qtq ∨Stjσ ∨St+1jσ0 )
T_
De plus, seule la case j peut être modifiée. Toutes les autres doivent
conserver le symbole qu’elles avaient au temps t.
x̄ ∨ y ∨ z̄ x ∨ ȳ y∨z
x̄ x z
z̄
y ȳ y
17
F IGURE 1.7 – Transformation de l’exemple de problème d’acceptabilité : (x̄ ∨ y ∨ z̄) ∧
(x ∨ ȳ) ∧ (y ∨ z) en un problème de stable.
20
Problème de l’ensemble stable : Données : G = (X, E), un graphe et k,
un entier Question : Existe-t-il un sous-ensemble X 0 ⊆ X, |X 0 | = k tel que :
e
∀i, j ∈ X 0 , (i, j) ∈
/ E (sous-ensemble de k sommets sans connexion) ?
br
La transformation d’un exemple de problème d’acceptabilité vers un en-
semble stable se fait de la façon suivante :
em
— On associe un sommet à tous les littéraux de toutes les clauses
— On relie en une clique (entièrement) tous les sommets associés à une
clause
— On relie les sommets-littéraux incompatibles (une variable et sa négation)
pt
Une telle transformation est illustrée en figure 1.7 pour un petit problème à trois
littéraux et trois clauses.
25
En cas d’égalité, chaque équipe reçoit un point. Étant donné une suite de
Al
scores pour chaque équipe, cette suite peut-elle être le résultat obtenu à la
fin d’un championnat ?
Remarque : si le gagnant reçoit seulement deux points, alors il existe un algo-
T_
17
n’est pas polynomial.
Par contre, les problèmes d’acceptabilité, du stable ou du cycle hamiltonien
ne sont pas des problèmes sur des nombres (car ils ne font simplement pas
20
intervenir de nombres).
On dit qu’un algorithme est pseudo-polynomial s’il s’exécute en un temps
borné par un polynôme dépendant de la taille n des données et du plus grand
nombre M intervenant dans le problème. Le problème de bipartition d’un en-
e
semble en deux sous-ensembles de poids égaux est un problème NP-complet
pour lequel il existe un algorithme pseudo-polynomial simple.
br
Exemple de problème de bipartition d’un ensemble : Est-il possible de
em
diviser l’ensemble {5, 2, 1, 6, 4} en deux sous-ensembles de poids égaux ? La
somme des poids pour cet exemple est de 18. On cherchera donc deux sous-
ensembles de poids 9. Pour résoudre ce problème, on crée un tableau de n
lignes, où n est le cardinal de l’ensemble, et M = 9 colonnes, où M est la
pt
ligne, mais décalé du poids du deuxième élément (ici : 2). On répète ensuite
Al
le procédé jusqu’à ce que tous les éléments aient été considérés. S’il y a un
× dans une case le la dernière colonne, c’est qu’il est possible de créer un
sous-ensemble de poids M , ce qui est le cas pour cet exemple. Une solution
T_
est : {2, 1, 6}{5, 4}. La complexité de l’algorithme est O(M · n), ce qui est bien
polynomial en n et en M .
Somme des poids
Élément 0 1 2 3 4 5 6 7 8 9
5 × ×
2 × × × ×
1 × × × × × × × ×
6 × × × × × × × × ×
4 × × × × × × × × × ×
Soit π, un problème sur des nombres et πp(n) les sous-problèmes de π res-
treints aux exemples pour lesquels M 6 p(n). L’ensemble πp(n) ne fait donc
32 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
17
est NP-complet, les sous-problèmes de voyageur de commerce ne faisant in-
tervenir que des petits nombres sont également NP-complets.
Par contre, les problèmes que l’on peut résoudre par programmation dyna-
20
mique, comme le problème du sac de montagne ou celui de la partition d’un
ensemble ne sont pas fortement NP-complets. En effet, si pour ce dernier pro-
blème la somme des poids des n éléments est bornée par un polynôme p(n),
l’algorithme ci-dessus aura une complexité polynomiale O(n · p(n)).
e
1.2.4
br
Autres classes de complexité
D’innombrables autres classes de complexité ont été proposées. Parmi
em
celles qui sont le plus fréquemment rencontrées dans la littérature et qui peu-
vent être décrites de façon intuitive, on peut citer :
Classe NP-difficile Les problèmes considérés ci-dessus sont des problèmes
pt
données dont la taille ne tiendrait pas dans la mémoire vive des ordina-
teurs.
Classe N C La classe N C est celle des problèmes pouvant se résoudre en
temps poly-logarithmique sur une machine disposant d’un nombre po-
lynomial de processeurs. Les problèmes de cette classe peuvent donc
résoudre en parallèle en un temps plus faible que celui dont on a be-
soin pour lire les données en séquentiel. Le tri des éléments d’un tableau
entre dans la catégorie N C.
17
Peu de résultats ont été établi en ce qui concerne les relations existantes
entre ces diverses classes de complexité. À l’exception des évidentes inclu-
sions au sens large L ⊆ P ⊆ N P ⊆ NP-complet ⊆ P-SPACE et N C ⊆ P , il
n’a été possible d’établir que L 6= P-SPACE. Il est communément conjecturé
20
que P 6= N P , mais cette conjecture fait partie des quelques problèmes du
millénaire.
e
br
em
pt
se
25
g
Al
T_
34 É LÉMENTS DE THÉORIES DES GRAPHES ET DE LA COMPLEXITÉ
Exercices
Exercice 1.1. On cherche à dessiner 5 segments sur le plan de manière à ce
que chaque segment en coupe exactement 3 autres. Formaliser ce problème
en termes de graphes.
n
Exercice 1.2. Simplifier les expressions suivantes : O(n5 + 2n ), O(5n + 22 ),
Ω(n2 · n! + (n + 2)!), Ω(nlog(log(n)) + 23n), O(nlog(n) + n5+cos(n) )
Exercice 1.3. Écrire un programme pour machine de Turing déterministe qui
17
reconnaît si la sous-chaîne ane est contenue sur la bande. L’alphabet d’entrée
est : Σ = {a, c, e, n}. Spécifier l’alphabet de bande Γ, l’ensemble Q des états
possibles de la machine et la fonction de transition δ.
20
Exercice 1.4. Démontrer que la recherche d’une clique de taille donnée dans
un graphe est un problème NP-complet.
e
commerce asymétrique en un problème symétrique en doublant le nombre de
villes.
br
em
pt
se
25
g
Al
T_
Chapitre 2
17
Exemples et modélisation de
20
problèmes combinatoires
e
Maintenant que nous avons passé en revue les principales définitions de
br
théorie des graphes et de théorie de la complexité, nous pouvons donner quel-
ques exemples de problèmes d’optimisation combinatoire. Certains de ces pro-
em
blèmes sont faciles, mais l’ajout d’une contrainte qui semble anodine peut les
rendre difficiles. Nous rappellerons également brièvement le principe de fonc-
tionnement de quelques algorithmes simples qui permettent de résoudre cer-
tains de ces problèmes.
pt
Les deux algorithmes les plus connus pour résoudre ce problème sont celui de
Kruskal et celui de Prim. Ils sont basés tous deux sur une méthode gloutonne :
à chaque étape on ajoute un élément à la structure que l’on construit, sans
jamais remettre en question ultérieurement le choix de cet élément.
L’algorithme de Kruskal (2.1) démarre avec un graphe T = (V, ET = ∅) et
ajoute à ET successivement une arête de E de poids aussi faible que possible
tout en veillant à ne pas créer de cycle. L’algorithme de Prim démarre avec un
graphe T = (V 0 = {s}, ET = ∅) et ajoute successivement un sommet v à V 0 et
une arête e à ET , telle que le poids de e soit aussi faible que possible et qu’une
de ses extrémité fasse partie de V 0 et l’autre non. Autrement dit, Kruskal part
d’une forêt constituée d’autant d’arbres qu’il y a de sommets et cherche à réunir
36 E XEMPLES ET MODÉLISATION DE PROBLÈMES COMBINATOIRES
Algorithme 2.1 : (Kruskal) Recherche d’un arbre de poids minimum. Pour une
implantation efficace, il faut utiliser une structure de données adéquate d’en-
sembles disjoints pour pouvoir tester si l’arête que l’on tente d’ajouter fait partie de
la même composante connexe ou non. Dans ce cas, la complexité de l’algorithme
est en O(|E| log |E|).
Entrées : R = (V, E, w) connexe
17
Résultat : Arbre recouvrant T = (V, ET ) de poids minimum
1 Trier et renuméroter les arêtes par poids non décroissant
w(e1 ) 6 w(e2 ) 6 · · · 6 w(e|E| );
20
2 ET = ∅;
3 pour k = 1 . . . |E| faire
4 si ET ∪ {ek } ne forme pas un cycle alors
5 ET ← ET ∪ {ek }
e
br
em
pt
Algorithme 2.2 : (Prim) Recherche d’un arbre de poids minimum. Pour une im-
plantation efficace, il faut utiliser une structure de données adéquate pour extraire
se
le sommet de L avec le plus petit poids (ligne 8) et pour modifier les valeurs des
poids (ligne 14). Un tas de Fibonacci ou une queue de Brodal permettent une
implantation de l’algorithme en O(|E| + |V | log |V |).
Entrées : R = (V, E, w) connexe, un sommet particulier s ∈ V
25
4 λs = 0 ; ET ← ∅;
Al
9 si i 6= s alors
10 ET ← ET ∪ {pi , i}
11 pour tous les Sommets v adjacent à i faire
12 si j ∈ L et λj > w({i, j}) alors
14 λj ← w({i, j});
15 pj ← i;
2.1 A RBRES OPTIMAUX 37
tous ces arbres en un seul alors que Prim part d’un arbre constitué d’un seul
sommet et cherche à le faire croître jusqu’à tous les intégrer.
17
ner consiste à relier un ensemble donné de points du plan par des lignes
dont la longueur est la plus faible possible. La figure 2.1 montre que l’arbre
de poids minimum qui n’utilise que des arêtes reliant directement les sommets
20
qui doivent figurer dans l’arbre peut avoir un poids supérieur à un arbre où l’on
a ajouté des nœuds de Steiner judicieusement choisis. Le choix des sommets
à ajouter engendre une combinatoire qui rend le problème NP-difficile.
e
br Nœuds de Steiner
em
Nœuds ordinaires
pt
F IGURE 2.1 – Arbre de poids minimum 2.1(a) sur un ensemble n’utilisant que les nœuds
ordinaires, qui sont donc directement reliés entre eux et arbre de Steiner de poids mini-
mum 2.1(b) où l’on peut utiliser des nœuds supplémentaires.
17
La recherche de chemins optimaux est un problème vieux comme le monde
et bien connu du grand public, en particulier depuis que les automobiles sont
20
équipées en série de systèmes de navigation automatisés. Connaissant sa po-
sition actuelle sur un réseau de transport, il s’agit de trouver le meilleur chemin
vers une destination donnée. Le critère usuellement retenu pour l’optimalité du
chemin est le temps, mais il peut également s’agir de la distance, notamment
e
s’il s’agit d’un trajet à faire à pied.
sive le qualificatif de longueur du chemin pour désigner cette somme des poids,
étant sous-entendu que les poids peuvent représenter autre chose qu’une dis-
tance, comme par exemple un temps, une consommation d’énergie, etc. Au
se
tra. Il est présenté dans l’algorithme 2.3. Pour une implantation efficace, il faut
Al
Algorithme 2.3 : (Dijkstra) Recherche d’un plus court chemin de s à tous les
autres sommets dans un réseau à pondération non négative. On a indiqué en
couleur les deux différences qu’il y a entre l’algorithme de Dijkstra et celui de
Prim.
Entrées : Réseau R = (V, E, w) orienté avec w(e) > 0 ∀e ∈ E, donné
sous la forme de listes des successeurs succ(i) de chaque
sommet i ∈ V , un sommet de départ particulier s
Résultat : Prédécesseur predj immédiat de j sur un plus court chemin
de s à j, ∀j ∈ V et longueur λj du plus court chemin de s à j.
17
1 pour tous les Sommets i ∈ V faire
2 λi ← ∞ ; // Longueur connue du plus court chemin de s à i
3 predi ← ∅ ;
20
4 λs = 0 ; L ← V ; // Sommets pour lesquels le plus court chemin n’est pas
définitif
5 tant que L 6= ∅ faire
7 Retirer de L le sommet i ayant le plus petit λi ;
e
8 pour tous les Sommets j ∈ succ(i) faire
11
12
9 si j ∈ L et λj > λi + w(i, j) alors
λj ← λi + w(i, j));
predj ← i br
em
f o r ( i = 0 ; i < n ; i ++)
{ lambda [ i ] = i n f i n i t e ; known [ i ] = i ; }
Al
lambda [ s ] = 0 ;
swap ( known +s , known + 0 ) ; / ∗ Only s h o r t e s t path t o known [ 0 ] = s a l r e a d y known ∗ /
{ f o r ( j = i ; j < n ; j = j +1)
i f ( lambda [ known [ i −1]] + d [ known [ i − 1 ] ] [ known [ j ] ] < lambda [ known [ j ] ] )
{ lambda [ known [ j ] ] = lambda [ known [ i −1]] + d [ known [ i − 1 ] ] [ known [ j ] ] ;
pred [ known [ j ] ] = known [ i −1];
j_kept = j ;
}
swap ( known + i , known + j _ k e p t ) ; / ∗ s h o r t e s t path i d e n t i f i e d t o known [ i ] ∗/
}
f r e e ( known ) ;
} /∗ d i j k s t r a ∗/
Lorsque l’on a des pondérations négatives, les plus courts chemins n’existent
que s’il n’y a pas de circuit de longueur négative dans le réseau. Un algorithme
plus général pour trouver ces plus courts chemins a été proposé par Bellman
40 E XEMPLES ET MODÉLISATION DE PROBLÈMES COMBINATOIRES
et Ford (voir l’algorithme 2.4). Il est basé sur la vérification, pour chaque arc,
des condition de Bellman : λj 6 λi + w(i, j). Autrement dit, la longueur du plus
court chemin de s à j ne doit pas dépasser celle de s à i à laquelle on ajoute la
longueur de l’arc (i, j). En effet, si c’était le cas, il existerait un chemin encore
plus court jusqu’à j, passant par i.
17
ce qui signifie que la longueur (négative) du plus court chemin n’est pas bornée.
Cet algorithme est excessivement simple à programmer (le code n’est guère plus
long que le pseudo-code donné ici) ; sa complexité est en O(|E||V |).
20
Entrées : Réseau R = (V, E, w) orienté, donné sous la forme de listes
d’arcs, un sommet de départ particulier s
Résultat : Prédécesseur predj immédiat de t sur un plus court chemin
de s à j avec sa longueur λj , ∀j ∈ V , ou indication de
e
l’existence d’un circuit de longueur négative
1
2
pour tous les i ∈ V faire
λi ← ∞ ; predi ← ∅
λs ← 0;
br
em
3
4 k←0; // Compteur du nombre d’étapes
5 Continuer ← vrai ; // Au moins un λ a été modifié à la dernière étape
6 tant que k < |V | et Continuer faire
Continuer ← f aux;
pt
7
8 k ← k + 1;
9 pour tous les arc (i, j) ∈ E faire
se
14 si k = |V | alors
15 Il y a un circuit de longueur négative accessible depuis s
g
Al
de partir d’une solution de départ complète, même très mauvaise mais facile à
construire, et de tenter de l’améliorer. La trame générale de cet algorithme est
celle d’une méthode d’amélioration locale. À chaque étape de l’algorithme, on
vérifie si les conditions de Bellman sont satisfaites pour l’ensemble des arcs. Si
c’est le cas, on a l’ensemble des plus courts chemins. Si on trouve un sommet
j pour lequel λj < λi + w(i, j), on met à jour le meilleur chemin connu jusqu’au
sommet j en lui fixant le sommet i comme sommet prédécesseur. Lorsque l’on
fait une telle modification, cela peut invalider la condition de Bellman pour des
arcs qui la satisfaisaient. Il faut donc revérifier pour l’ensemble des arcs si une
modification n’en implique pas d’autres.
Se pose alors la question : Sans autre précaution, un algorithme basé sur
2.2 C HEMINS OPTIMAUX 41
17
tion avec un critère d’arrêt bien défini : Si l’on a fait |V | étapes et qu’il y a
encore eu des modifications, alors on peut s’arrêter en indiquant que le ré-
seau possède un circuit de longueur négative. Si à une étape on vérifie que les
20
conditions de Bellman sont satisfaites pour tous les arcs, alors on peut s’arrêter
car on a l’ensemble des plus courts chemins.
La recherche de chemins optimaux apparaît dans de nombreuses appli-
cations, en particulier en planification de projet et en ordonnancement. Les
e
problèmes qui peuvent être résolus par une technique connue sous le nom de
programmation dynamique peuvent se formuler sous la forme de la recherche
br
d’un chemin optimal dans un réseau en couches. Cette technique utilise la to-
pologie particulière du réseau pour trouver la solution du problème sans avoir
à construire explicitement le réseau.
em
Formulation du plus court chemin en termes de PL
pt
est utilisé par le plus court chemin. La formulation ci-dessous peut paraître
incomplète : en effet, on voudrait que les variables xij prennent soit la valeur
0 (auquel cas l’arc (i, j) ne fait pas partie du plus court chemin), soit la valeur
1 (l’arc en fait partie). Les contraintes 2.4 sont suffisantes, car si une variable
25
prend une valeur fractionnaire, cela signifie qu’il existe plusieurs plus courts
chemins de s à t. La contrainte 2.3 impose en effet qu’il y ait une « quantité »
unité qui arrive en t. Cette quantité peut être fractionnée, mais chaque fraction
doit utiliser un plus court chemin, à l’optimum. Les contraintes 2.2 imposent
g
que pour tout sommet intermédiaire j, la quantité arrivant dans ce sommet doit
Al
X
Minimiser z= w(i, j)xij (2.1)
i,j
n
X n
X
xij − xjk = 0 ∀j 6= s, j 6= st (2.2)
i=1 k=1
Xn Xn
xit − xtk = 1 (2.3)
i=1 k=1
xij > 0 ∀i, j (2.4)
42 E XEMPLES ET MODÉLISATION DE PROBLÈMES COMBINATOIRES
Maximiser λt (2.5)
Sous λj − λi 6 w(i, j) ∀i, j (2.6)
Contraintes λs = 0 (2.7)
17
La dualité joue un rôle important en programmation linéaire. On peut en
effet montrer que toute solution admissible du premier problème a une valeur
20
qui ne peut être inférieure à une solution admissible du second. Si l’on trouve
une solution au premier problème qui vaut exactement une solution du second,
alors ces deux solutions sont optimales. Pour le problème du plus court che-
min, la valeur optimale de λt correspond à la somme des longueurs des arcs
e
qu’il faut utiliser dans un chemin optimum de s à t.
2.2.2 br
Plus court chemin élémentaire
em
Le problème du plus court chemin est mal défini, en raison des circuits de
longueur négative. Par contre, on pourrait ajouter une contrainte très naturelle,
qui le rend parfaitement défini : chercher le plus court chemin élémentaire d’un
pt
sommet particulier s à tous les autres. Dans ce cas, même s’il y a des circuits
de longueur négative, le problème a bien une solution finie, car le fait d’exiger
que le chemin soit élémentaire empêche de passer deux fois pas le même
se
un réseau R = (V, E, w), dont l’ensemble V peut être interprété comme des
villes, on cherche un circuit de longueur minimale passant une et une seule fois
par chacune des villes. Ce problème est l’archétype de l’optimisation combina-
T_
17
optimal l’est également.
1 s t
20
10 15
−40 −35
2 5 2 5
32 25 −18
12 −23 −37
e
18 29 −32 −21
3
27 13
4 br −38
3
−25
4
em
9 −41
F IGURE 2.2 – Transformation polynomiale d’un problème de voyageur de commerce
vers un problème de plus court chemin élémentaire.
pt
X
Minimiser z= w(i, j)xij (2.8)
g
(i,j)
Al
n
X
xij = 1 ∀j (2.9)
i=1
T_
Xn
xij = 1 ∀i (2.10)
j=1
Xn
y1j = n − 1 (2.11)
j=2
n
X Xn
yik − xkj = 1 ∀k 6= 1 (2.12)
i=1 j=1
La contrainte 2.9 impose que l’on entre exactement une fois dans chaque
ville et la contrainte 2.10 qu’on en ressorte exactement une fois. Le flux à tra-
vers l’arc (i, j) est donné par les variables yij . On impose que la somme des
flux sortant du sommet 1 par 2.11 et la consommation unitaire de chaque som-
met par 2.12. Finalement, on empêche de faire passer un flux à travers un
arc non utilisé par 2.13. Cette formulation n’est sans doute pas la plus effi-
cace pour résoudre pratiquement des problèmes de voyageur de commerce,
mais elle présente relativement peu de variables (2 par arc) et un nombre de
contraintes proportionnel au nombre d’arcs du problème.
17
2.2.4 Tournées de véhicules
20
Des problèmes utilisant celui du voyageur de commerce comme sous-
problème apparaissent naturellement dans l’élaboration de tournées de vé-
hicule. Dans sa forme la plus simple, celui-ci peut se formuler ainsi : Un en-
semble V de clients demandant des quantités qi de biens (i = 1, . . . , |V |)
e
doivent être desservis avec un véhicule de capacité Q à partir d’un entre-
pôt d. Il faut déterminer une décomposition de l’ensemble des clients en m
br
P
sous-ensembles V1 , . . . Vm tels que i∈Vi qi 6 Q et trouver pour chaque sous-
ensemble Vj ∪ {d}, (j = 1, . . . , m) une tournée de voyageur de commerce de
em
sorte que la distance totale parcourue par le véhicule soit aussi faible que pos-
sible.
pt
se
25
Dépôt
g
Al
T_
17
— Une tournée mêle collecte et distribution
— On a plus d’un entrepôt
20
— Les entrepôts abritent une flotte de véhicules hétérogènes ou non
— La position des entrepôts peuvent être choisies
— Etc.
Comme il s’agit de trouver dans quel ordre desservir les clients, les pro-
e
blèmes d’élaboration de tournées font également partie des problèmes d’or-
br
donnancement. Du reste, dans la littérature anglo-saxonne, on rencontre aussi
bien le terme de « vehicle routing » que celui de « vehicle scheduling ».
em
2.3 Ordonnancement
Les problèmes d’ordonnancement regroupent tous ceux dont le but est de
pt
sur une machine. Des opérations ayant une forte interdépendance sont regrou-
pées en tâches. Par exemple, une tâche peut représenter l’ensemble des opé-
rations qu’un objet doit subir. Souvent, ces dernières ne peuvent être réalisées
dans n’importe quel ordre ; on parle alors de gamme opératoire.
25
etc.
Al
17
blèmes de collision. Il peut y avoir plusieurs machines du même type, des ma-
chines pouvant réaliser diverses opérations, etc.
20
Exemple : atelier de peinture Un des problèmes d’ordonnancement les plus
simples à formule est celui de la peinture d’objets. Il n’y a ici qu’une seule res-
source : la machine qui permet de déposer la peinture sur les objets. Chaque
tâche ne comporte qu’une opération : peindre un objet i (i = 1, . . . , n) dans
e
une couleur donnée, ce qui prend un temps ti . Après avoir peint un objet i, il
br
faut nettoyer la machine pour pouvoir peindre le suivant, j, dans une autre cou-
leur, ce qui prend un temps sij . Notons que généralement sij 6= sji : en effet,
des traces d’une teinte foncée dans une couleur claire ont plus d’impact que
em
le contraire. Après avoir peint tous les objets, il faut nettoyer la machine pour
qu’elle soit prête pour le jour suivant, ce qui prend un temps r. L’objectif est de
trouver l’ordre dans lequel faire passer les objets sur la machine pour minimi-
ser l’heure à laquelle on P termine le travail. On cherche donc une permutation p
pt
n−1
des n objets minimisant i=1 (tpi + spi pi+1 ) + tpn + r.
Ce problème d’ordonnancement se ramène à celui d’un voyageur de com-
se
merce à n + 1 villes. Pour cela, on pose wij = ti + sij (i, j = 1, . . . , n), wi0 = r
et w0i = 0, (i = 1, . . . , n). On vérifie que la tournée la plus courte sur les villes
0, . . . , n donne l’ordre optimal dans lequel il faut peindre les objets, en interpré-
tant l’« objet » 0 comme le début d’une journée de travail.
25
dans cet ordre. Un objet j doit donc subir m opérations qui prennent un temps
tij , (i = 1, . . . , m, j = 1, . . . , n). Il s’agit de déterminer dans quel ordre faire
passer les objets sur la chaîne de traitement, autrement dit, trouver une per-
mutation des objets telle que la durée des travaux soit minimisée. Il faut noter
qu’il y a un espace tampon entre chaque poste, où l’on peut au besoin stocker
des objets en attendant que le poste termine le traitement d’un objet arrivé
plut tôt. Une façon commode de représenter une solution pour la planification
de tâches est le diagramme de Gannt dans lequel l’axe des abscisses repré-
sente le temps et celui des ordonnées les ressources. La figure 2.4 donne le
diagramme de Gannt d’une solution quelconque, où chaque opération a été
planifiée le plus tôt possible ainsi qu’un ordonnancement optimal au plus tard.
2.3 O RDONNANCEMENT 47
Machine
Temps
Étendue des travaux
(a) Ordonnancement quelconque au plus tôt
17
Machine
20
Temps
Étendue des travaux
(b) Ordonnancement optimal au plus tard
e
br
F IGURE 2.4 – Chaîne de traitement à gamme opératoire unique. Diagramme de Gannt
d’un ordonnancement quelconque au plus tôt pour un petit exemple à 4 ressources et
em
5 tâches 2.4(a) et ordonnancement optimal au plus tard 2.4(b) pour le même exemple.
dans une liste toutes les opérations par durée croissante. On sélectionne de
la liste l’opération dont la durée est la plus petite, si cette opération a lieu sur
se
Minimiser dω (2.15)
dmj + tmj 6 dω (j = 1 . . . n) (2.16)
T_
terminé son traitement sur une machine i avant de passer sur la machine i + 1
(2.17). Les variables yjk , indiquent si l’objet j doit être traité avant l’objet k ; il
est nécessaire de n’introduire que n·(n−1)/2 de ces variables, étant donné que
ykj devrait prendre la valeur complémentaire 1 − yjk . Les deux contraintes 2.18
et 2.19 font intervenir une grande constante M et permettent l’introduction de
contraintes disjonctives : soit l’objet j passe avant l’objet k, soit c’est l’inverse.
Si yjk = 1, j passe avant k et les contraintes 2.19 sont trivialement satisfaites
pour toute machine i, si M est suffisamment grand. Inversement, si yjk =
0, ce sont les contraintes 2.18 qui sont trivialement satisfaites alors que les
17
contraintes 2.19 imposent de terminer le traitement de l’objet k sur la machine
i avant que cette dernière ne puisse commencer le traitement de l’objet j.
20
2.3.2 Gamme opératoire variable
Un problème un peu plus général est celui de l’ordonnancement à gamme
opératoire variable, connu dans la littérature anglo-saxonne sous le terme de
e
« jobshop ». Ici, chaque tâche doit subir un certain nombre d’opérations sur
br
des machines données. La séquence des opérations pour une tâche est fixée,
mais deux tâches n’ont pas forcément la même séquence et toutes les tâches
ne passent pas forcément sur toutes les machines.
em
Machine 3
6 4 3
3
pt
se
0 4
α 4 7 ω
0 2
2 2 3
Machine 2
25
4
3
1 5
Machine 1
g
à gamme opératoire variable. Une tâche est composée de 3 opérations et deux autres
de 2. Trois machines sont impliquées. La pondération des arcs correspond à la durée
de l’opération correspondante. Les arcs représentant les relations de précédence des
T_
opérations d’une même tâche sont en pointillés. Le chemin critique, c’est-à-dire le plus
long chemin de α à ω est indiqué en gras.
Minimiser dω (2.22)
17
di + ti 6 dj ∀(i, j) (2.23)
di + ti 6 dω ∀i (2.24)
di + ti 6 dk + M · (1 − yik ) ∀i, k sur une même machine (2.25)
20
dk + tk 6 di + M · yik ∀i, k sur une même machine (2.26)
di > 0 ∀i (2.27)
yik ∈ {0, 1} ∀i, k sur une même machine (2.28)
e
Les contraintes 2.23 imposent de terminer l’opération i qui précède l’opé-
br
ration j d’une tâche donnée avant de pouvoir débuter le traitement de cette
dernière. Les contraintes 2.24 imposent que les heures de fin de traitement de
em
toutes les opérations précèdent la fin du traitement. Les variables yik associées
aux contraintes disjonctives 2.25, 2.26 déterminent si l’opération i précède ou
non l’opération k qui a lieu sur la même machine.
pt
clavier, des postes de travail à des employés. Certains de ces problèmes peu-
vent être résolus efficacement en les modélisant sous la forme de problèmes
de flots dans des réseaux.
g
(i, j) ∈ E de sorte à maximiser la somme des flux qui sortent d’un sommet
particulier s pour aller vers un autre sommet t tout en respectant les contraintes
de conservation pour les autres sommets — la somme des flux entrants est
égale à celle des flux sortants — et de capacité — les flux xij doivent être
positifs, sans dépasser la valeur w(i, j) associés aux arcs.
Pour résoudre ce problème, Ford et Fulkerson ont proposé un algorithme
assez simple ( 2.5), basé sur une méthode d’amélioration : On part d’un flot nul
(qui est admissible) et on l’augmente à chaque étape, le long d’un chemin de
s à t jusqu’à obtenir le flot optimum. La première étape de cet algorithme est
illustrée dans la figure 2.6. Pour ne pas se retrouver bloqué dans une situation
où il n’existe plus de chemin de s à t mais sans que le flot soit maximal, comme
50 E XEMPLES ET MODÉLISATION DE PROBLÈMES COMBINATOIRES
Capacité = 2
s 1 2 t
Capacité = 1
(a) Flot de départ nul
17
Flux = 0
20
s 1 2 t
e
Flux = 1
(b) Flot bloquant
br
em
F IGURE 2.6 – Représentation graphique d’un problème de flot. Un triangle vide indique
une unité de capacité non utilisée alors qu’elle l’est pour un triangle rempli. L’algorithme
de Ford et Fulkerson part d’un flot nul 2.6(a) et l’augmente au maximum sur un chemin
de s à t. Si le premier chemin trouvé est s − 1 − 2 − t on ne peut augmenter le flot que
d’une unité le long de ce trajet. On obtient un flot bloquant 2.6(b), car il n’existe plus de
pt
s 1 2 t
g
Al
s 1 2 t
V \A
17
bloquant de la figure 2.6(b), ce qui permet d’obtenir le flot optimum. La preuve de l’op-
timalité découle du fait que la capacité de la coupe A séparant s de t indiquée en ligne
de traits est égale à la valeur du flot.
20
Algorithme 2.5 : (Ford-Fulkerson) Recherche d’un flot maximum de s à t.
Entrées : Réseau R = (V, E, w) orienté, un sommet-source s et un
e
sommet-puits t
Résultat : Flot optimum de s à t
1 Partir d’un flot nul dans tous les arcs;
2 répéter br
em
3 Construire le réseau des augmentations R∗ correspondant au flot
courant;
4 si Un chemin de s à t existe dans R∗ alors
5 Chercher le flux maximum que l’on peut faire passer le long du
pt
chemin de s à t dans R∗ ;
6 Superposer ce flux au flot courant (diminuer le flot dans les arcs
(i, j) de R qui apparaissent inversés (j, i) sur le chemin de R∗ )
se
d’une unité, pour des réseaux à capacités entières. Si le flot maximum dans un
Al
Dans le cas d’un réseau dense (|E| ∈ O(|V 2 |)), cette complexité se simplifie
en O(|V |5 ).
De nombreux algorithmes ont été proposés pour la résolution du problème
du flot maximum. Au niveau de la complexité, celui de Malhotra, Kumar, Mahe-
shwari, en O(|V |3 ) est un des meilleur.
Dans de nombreuses applications, faire passer une unité de flux à travers
un arc (i, j) coûte c(i, j). On a donc un réseau R = (V, E, w, c), où w(i, j) est la
capacité de l’arc (i, j) et c(i, j) le coût unitaire d’une unité de flux traversant cet
17
arc. Se pose alors le problème du flot maximum à coût minimum. Ce problème
peut être résolu avec l’algorithme de Busacker et Gowen 2.6 lorsque le réseau
de départ ne possède pas de circuit à coût négatif.
20
Algorithme 2.6 : (Busacker et Gowen) Recherche d’un flot maximum à coût
minimum de s à t.
e
Entrées : Réseau R = (V, E, w, c) orienté, sans circuit à coût négatif, un
1 br
sommet-source s et un sommet-puits t
Résultat : Flot maximum à coût minimum de s à t
Partir d’un flot nul dans tous les arcs;
em
2 répéter
3 Construire le réseau des augmentations R∗ correspondant au flot
courant;
4 si Un chemin de s à t existe dans R∗ alors
pt
Comme ce qui avait été constaté pour les algorithmes de Prim et Dijks-
tra, on remarque une différence minime entre les algorithmes de Ford & Ful-
kerson 2.5 (ou plutôt son amélioration par Edmons & Karp) et Busacker &
g
cité soit positive) avec coût c(i, j) inchangé et un arc inversé (j, i) de capacité
xij et de coût −c(i, j).
Dans le cas général, le problème du flot maximum à coût minimum est N P -
difficile. On peut en effet transformer polynomialement le problème du voyageur
de commerce dans ce problème, de façon similaire à la transformation du plus
court chemin élémentaire (voir figure 2.2).
Si les algorithmes pour trouver des flots optimaux présentés ci-dessus peu-
vent résoudre de nombreux problème directement liés à la gestion de flux,
comme l’acheminement d’énergie électrique ou certains problèmes de trans-
port, leur principale application est le traitement de problèmes d’affectation.
2.4 P ROBLÈMES D ’ AFFECTATION 53
n X
n
17
X
Minimiser ciu xiu (2.29)
i=1 u=1
Sous contraintes
20
Xn
xiu = 1 u = 1, . . . , n (2.30)
i=1
Xn
xiu = 1 i = 1, . . . , n (2.31)
e
u=1
xiu ∈ {0, 1}
br (i, u = 1, . . . , n)
17
pertes de temps peut se modéliser par le programme quadratique en variables
0 − 1 suivant, dans lequel la variable xiu prend la valeur 1 si l’employé i occupe
le bureau u et la valeur 0 sinon :
20
n X
X n X
n X
n
Minimiser aiu bjv xiu xjv (2.33)
i=1 j=1 u=1 v=1
e
Sous contraintes
n
br
X
xiu = 1 u = 1, . . . , n (2.34)
i=1
n
em
X
xiu = 1 i = 1, . . . , n (2.35)
u=1
xiu ∈ {0, 1} (i, u = 1, . . . , n) (2.36)
pt
la recherche d’une
Pn permutation
Pn p se formule :
Minimiser i=1 j=1 aij · bpi pj .
De nombreuses applications peuvent le formuler comme un problème d’af-
25
fectation quadratique :
Attribution de bureaux à des employés C’est l’exemple qui vient d’être donné.
Attribution de blocs dans un FPGA La programmation d’un réseau de
portes programmables in situ (Field Programmable Gate Array en an-
g
glais) nécessite de connecter entre eux des blocs logiques présent sur
Al
une puce de silicium. Ces blocs permettent d’implanter des équations lo-
giques, des multiplexeurs ou des registres à décalage. Lorsque l’on veut
programmer un FPGA pour qu’il réalise une fonctionnalité donnée, on
T_
commence par établir un schéma dans lequel des on connecte des mo-
dules. On peut décrire ce schéma au moyen d’une matrice de routage
A = (aij ) qui donne le nombre de connexions entre les modules i et j.
Ensuite, il faut attribuer à chaque module i un bloc logique pi sur la puce.
Cependant, on doit prendre soin de ne pas réaliser cette affectation n’im-
porte comment car le temps de propagation des signaux dépend de la
longueur des liaisons. Connaissant la longueur buv de la liaison entre les
blocs logiques u et v, le problème de la minimisation de la somme des
temps de propagation est donc un problème d’affectation quadratique.
Configuration d’un clavier Pour saisir un texte sur le clavier d’un téléphone,
les 26 lettres de l’alphabet ainsi que l’espace on été attribuées aux touches
2.5 C OLORATION DE GRAPHES 55
17
(a) Clavier standard (b) Clavier optimisé
20
F IGURE 2.9 – Clavier de téléphone standard 2.9(a) et clavier optimisé pour la saisie de
texte français 2.9(b)
e
En supposant que la frappe d’une touche prenne une unité de temps, que
le déplacement d’une touche à une autre prenne deux unités de temps
br
et qu’il faille attendre 6 unités de temps avant de pouvoir commencer à
saisir un nouveau symbole placé sur la même touche que le précédent,
em
on peut calculer qu’il faut 70 unités de temps pour taper le texte « a ce
soir bisous » : en effet, il faut 1 unité pour saisir le « a » sur la touche 2, se
déplacer sur la touche 0, ce qui prend deux unités, appuyer une fois sur
cette touche pour l’espace, 1 unité, se déplacer à nouveau sur la touche
pt
qu’il faut seulement 51 unités de temps, soit près d’un tiers de moins. Ce
clavier optimisé a été obtenu en résolvant un problème d’affectation qua-
dratique pour lequel les cœfficients aij représente la fréquence d’appari-
tion du symbole j après le symbole i dans un texte-type et buv représente
25
3 e 2 e
a b
4
d f a b
c
1 G → L(G) d c f
g
17
5 g
F IGURE 2.10 – Graphe et son graphe représentatif des arêtes, montrant que les pro-
20
blème de la coloration des sommets et des arêtes sont des problèmes équivalents
arêtes d’un graphe, deux arêtes ayant un sommet incident commun doivent
e
recevoir des couleurs différentes. Pour la coloration des sommets, on veut que
br
deux sommets adjacents reçoivent une couleur différente.
Les deux problèmes sont équivalents : À partir d’un graphe G = (V, E), on
peut construire le graphe représentatif des arêtes L(G) = (E, F ), où il existe
em
une arête f ∈ F entre les sommets e et e0 ∈ E de L(G) si ces deux arêtes de G
ont un sommet en commun. Cette construction est illustrée dans la figure 2.10.
Ainsi à une coloration des arêtes de L(G) correspond une coloration des som-
mets de G et vice-versa.
pt
17
lorsque l’on associe un poids à chaque sommet. On cherche dans ce cas un
sous-ensemble de sommets indépendants dont la somme des poids est la plus
élevée possible. Naturellement, si le poids associé à chaque sommet est uni-
20
taire, cette variante est équivalente au stable maximal.
e
br
em
pt
F IGURE 2.11 – Problème de l’étiquetage d’un plan : On désire placer le nom de trois
se
objets sur un plan de sorte que les textes ne se recouvrent pas, pour en conserver la li-
sibilité. Ici, on retient 4 possibilités pour le placement de l’étiquette d’un objet. En termes
de graphe, ce problème peut se modéliser sous la forme d’un stable maximum, les som-
mets du stable correspondant aux places finalement choisies pour les étiquettes. Dans
cet exemple, le nom de chaque ville pourrait être placé en haut et à droite.
25
gauche. Les préférences peuvent être modélisées sous la forme de poids as-
sociées aux position. Le problème de l’étiquetage consiste alors à trouver un
stable de poids maximum.
D’autres problèmes peuvent se modéliser exactement de la même ma-
nière : on mentionnera tout d’abord celui de l’attribution de la position sur un
quai pour un bateau faisant escale dans un port. Traduit en termes d’étique-
tage, une étiquette aura comme largeur la durée prévue de l’escale et comme
hauteur la longueur de quai nécessaire pour ce bateau. Les positions possibles
pour cette étiquette sont déterminées par l’heure à laquelle le bateau arrive et
17
les emplacements du ou des quais pouvant recevoir le bateau.
Une autre application de ce problème est l’attribution de niveaux de vol
pour des avions commerciaux. Connaissant les heures de départ prévues de
20
chaque avion et les trajets qu’ils empruntent, on peut déterminer des zones de
collision potentielle entre avion. La dimension de ces zone dépend de l’incer-
titude sur les heures de départ effectives et les trajets réellements suivis. Une
étiquette correspondra donc à une zone. Les positions possibles des étiquettes
e
seront les divers niveaux de vol que l’avion pourrait utiliser.
2.7 Classification br
em
pt
se
25
g
Al
T_
17
construction d’arbres phylogénétiques, popularisé par Darwin au XIXe siècle,
qui consiste à trouver des relations de « parenté » entre les éléments, chaque
nœud de l’arbre représentant un « ancêtre » commun.
20
En toute généralité, la classification consiste à regrouper les éléments qui
sont similaires et à séparer ceux qui ne le sont pas. Cela suppose que l’on
puisse quantifier la « dissimilarité » d(i, j) 6 0 qu’il y a entre deux éléments i
et j appartenant à l’ensemble E que l’on cherche à classer. Souvent, la fonc-
e
tion d(i, j) est une distance, (si l’on a symétrie : d(i, j) = d(j, i), séparation :
d(i, j) = 0 ⇐⇒ i = j et inégalité triangulaire : d(i, k) 6 d(i, j) + d(j, k)), mais
br
pas nécessairement. Pour garantir la stabilité des algorithmes, on supposera
toutefois que d(i, j) > 0 et d(i, i) = 0, ∀i, j ∈ E.
Sitôt que l’on dispose une telle fonction, on peut mesurer l’homogénéité
em
d’un groupe G ⊂ E. Plusieurs définitions ont été proposées pour cette mesure :
Diamètre d’un groupe Valeur maximale de la fonction d(i, j) pour deux enti-
tés i et j appartenant à G : maxi,j∈G d(i, j)
pt
Clique
PSommeP des dissimilarités entre toutes les paires d’éléments de G :
i∈G j∈G d(i, j)
2.7.1 P-médiane
Le problème de la p-médiane est un des plus connus en classification non
supervisée. En reprenant les définitions présentées ci-dessus, il s’agit de mini-
miser la somme des étoiles. EnP d’autres termes, il faut trouver les p éléments
c1 , . . . , cp de E qui minimisent : i∈E minr=1,...,p d(i, cr ) Ce problème est N P -
17
difficile.
Un algorithmes heuristiques pour ce problème est la méthode d’améliora-
tion PAM (Partition Around Medoids) 2.7.
20
Algorithme 2.7 : (PAM) Méthode d’amélioration pour le partitionnement autours
de médians.
Entrées : Ensemble E d’éléments avec une fonction mesurant la
e
dissimilarité d(i, j) entre éléments i et j ; p éléments médians
c1 , . . . , cp de E
1 répéter br
Résultat : Groupes G1 , . . . , Gp ⊂ E
em
2 Associer chaque élément i ∈ E au médian le plus proche, et créer
les groupes G1 , . . . , Gp ⊂ E;
3 pour chaque Médian cj faire
4 pour chaque Élément i ∈ E faire
pt
médians;
7 jusqu’à ce que les p médians ne bougent plus;
25
17
3 pour chaque i ∈ 1, . . . , n faire
4 ci = centre de gravité de Gi
jusqu’à ce que les p centres ne bougent plus;
20
5
d’une méthode d’amélioration qui démarre avec des centres placés généra-
e
lement aléatoirement. Il alterne une phase d’association des éléments à leur
br
centre le plus proche (ligne 2) et une phase de repositionnement optimal des
centres pour chaque groupe précédemment construit (ligne 4). L’algorithme
s’arrête lorsqu’il n’y a plus de modification dans la position des centres. Cet
em
algorithme est relativement rapide, en Ω(p · n), mais il est très sensible à la
solution initiale fournie ainsi qu’aux mesures extrêmes. Si la dispersion des
éléments est élevée, on a intérêt à utiliser une autre mesure de la dissimilarité,
par exemple la distance ordinaire (norme l1 ), ce qui implique de modifier la
pt
n
X
Maximiser ci · x i (2.37)
i=1
T_
Sous contraintes
n
X
vi · xi 6 V (2.38)
i=1
xi ∈ {0, 1} (i = 1, . . . , n) (2.39)
17
sation d’un même problème, celui de la coloration des sommets d’un graphe.
20
e
br
F IGURE 2.13 – Coloration d’un graphe avec un nombre minimum de couleur. On peut
em
montrer que cette coloration est optimale, car le graphe contient une clique de quatre
éléments.
pt
néaire présenté en page 12 sont les suivantes : l’objectif (2.40) consiste à mini-
miser (1.4), qui est équivalente à minimiser l’indice de la couleur la plus élevée
attribué à un sommet du graphe. La contrainte (2.41) résume l’ensemble des
contraintes (1.5) (1.6) (1.7) (1.8).
De manière un peu moins formelle, on peut exprimer le problème de la
coloration de graphe comme :
Minimiser : c = f1 (s) (2.42)
Sous contraintes : 2 sommets adjacents ont des couleurs différentes
(2.43)
et : Nombre de couleurs utilisées par s 6 c (2.44)
2.9 M ODÉLISATION : F ONCTION - OBJECTIF ET FONCTION - UTILITÉ 63
F IGURE 2.14 – Coloration d’un graphe avec une couleur de trop. Un certain nombre de
17
changements sont nécessaires pour la supprimer. Pour cette solution, on a f1 = 5 alors
que pour celle de la figure 2.13 on a f1 = 4.
20
Un objectif consistant à minimiser un maximum — ou à maximiser un mini-
mum — comme (2.42) n’est pas propice à la découverte ou à l’amélioration de
solutions satisfaisant toutes les contraintes. Que faire si une solution que l’on
e
a découverte utilise une couleur de plus que l’optimum ? On ne sait pas s’il y
a juste un sommet qui utilise la couleur supplémentaire ou s’il y en a un grand
nombre.
br
em
pt
se
F IGURE 2.15 – Coloration non admissible d’un graphe avec un nombre donné de cou-
leurs. Pour cette solution, on a f2 = 4 + 2λ alors que pour celle de la figure 2.13 on a
25
grant des contraintes qui ont été relaxées. En l’occurrence, on peut relaxer (2.43)
et formuler le problème, où λ est une constante suffisamment grande :
Al
Il est facile de voir que pour une valeur de λ suffisemment grande (par
exemple en posant λ = nombre chromatique), une solution optimisant (2.45)
est également optimale pour (2.42). Ainsi, un triangle coloré avec une seule
couleur aura une valeur de 1 + 3λ, avec deux couleurs on aura une valeur de
2 + λ alors que la coloration optimale a une valeur de 3. Pour 0 < λ < 1/2,
l’optimum de f2 a lieu pour une couleur, pour 1/2 < λ < 1 pour deux couleurs
et pour λ > 1 pour trois couleurs. Par exemple, la solution de la figure 2.13
a une valeur de f2 = 4, celle de la figure 2.14 est de f2 = 5 et celle de la
figure 2.15 est de f2 = 4 + 2λ.
64 E XEMPLES ET MODÉLISATION DE PROBLÈMES COMBINATOIRES
17
20
F IGURE 2.16 – Coloration partielle d’un graphe avec un nombre donné de couleurs.
e
Pour cette solution, on a f3 = 1.
br
Au lieu de relaxer la contrainte stipulant que deux sommets adjacents ne
doivent pas recevoir la même couleur, on aurait pu relaxer celle qui impose
em
d’attribuer une couleur à chaque sommet :
F IGURE 2.17 – Coloration d’un graphe obtenue par orientation sans circuit de ses
arêtes. Le plus long chemin est indiqué en gras. Pour cette solution, on a f4 = 4 (pour
une coloration en cinq couleurs).
En effet, une fois que l’on a une telle orientation, on peut facilement trou-
ver une coloration admissible : Les sommets sans prédécesseur reçoivent la
couleur 1 ; ils ne peuvent être connectés par une arête, donc il n’y a pas de vio-
17
lation de contrainte. On peut les retirer du graphe avant d’attribuer la couleur 2
à ceux qui restent et qui n’ont pas de prédécesseur, etc. Ainsi, le nombre de
couleurs obtenues vaut un de plus que la longueur du plus long chemin.
20
Le choix d’une modélisation a une influence majeure sur la capacité de ré-
solution du problème par une méthode donnée. Il convient donc de prêter une
très grande attention à la phase d’analyse du problème. Par exemple, si on
choisit de résoudre un problème d’établissement d’horaire à l’aide de la colo-
e
ration de graphe, l’utilisateur ne cherchera pas nécessairement à trouver un
horaire utilisant la plus petite grille possible — ce qui se traduit par la recher-
br
che d’une coloration avec un nombre minimum de couleur. Il spécifiera plutôt
que la grille comportera un nombre donné de périodes et que le problème con-
siste à trouver un horaire pouvant s’intégrer dans cette grille. Dans ce cas, des
em
fonctions-utilité du type f2 (avec λ proche de 0 et c fixé au nombre de périodes
désirées) ou f3 conviendront certainement mieux que f1 ou f4 .
La technique de relaxation de contraintes par l’intégration d’une mesure
de leurs violations dans une fonction-utilité ne peut être envisagée que pour
pt
férences plutôt que des conditions rédhibitoires. Une approche plus souple et
fréquemment employée est de considérer plusieurs objectifs simultanément et
de laisser à l’utilisateur le choix du compromis — choisir une solution privi-
légiant tel objectif plutôt que tel autre. Mais pour cela, les méthodes doivent
25
pouvoir proposer plusieurs solutions variées en lieu et place de celle qui opti-
miserait un unique objectif.
g
17
2.10.1 Ensemble Pareto
20
Aéroport 1
30/120 200/45
e
15/12 35/15
Aéroport 2 250/55
110/50
pt
F IGURE 2.18 – Exemple de problème multi-objectifs : Pour se rendre d’un lieu de départ
à une destination, on a le choix entre plusieurs itinéraires. Chaque itinéraire a un coût
et une durée donnés, indiqués en regard des arcs. Depuis le lieu de départ, on peut se
se
rendre à la gare, soit en bus, soit en taxi, ou encore aller directement en taxi à l’aéroport
le plus proche. Depuis la gare, on peut rejoindre soit l’aéroport le plus proche, soit un
autre aéroport, un peu plus éloigné mais mieux desservi et moins cher. Ensuite, on
prend l’avion pour arriver à l’aéroport le plus proche de la destination, où l’on prend à
25
17
20
e
br
F IGURE 2.19 – Représentation des solutions du problème de la figure 2.18 dans un
diagramme coût/temps. Les 9 solutions non-dominées sont mises en évidence. Seules
les 5 solutions supportées pourront être découvertes par une méthode scalaire exacte.
em
2.10.2 Scalarisation
pt
K
X
Minimiser : wi · fi (s) (2.49)
g
i=1
Sous contraintes : s∈S
Al
P
Généralement on utilise 0 6 wi 6 1 et i wi = 1. Cette technique est
T_
17
Minimiser : f~(s) = (f2 (s), . . . , fK (s)) (2.50)
Sous contrainte : f1 (s) 6 v1
20
avec : s ∈ S
e
supportée.
br
em
pt
se
25
g
Al
T_
2.10 O PTIMISATION MULTIOBJECTIFS 69
Exercices
Exercice 2.1 (Connexion de points du plan). Un ensemble V de points du plan
doivent être connectés. Comment réaliser les connexions pour que la longueur
totale de ces dernières soit aussi faible que possible ? Application : considérer
les 3 points (0, 0), (30, 57), et (66, 0).
Exercice 2.2 (Accessibilité par des poids lourds). Dans le réseau routier de la
figure 2.20, la charge maximale (en tonnes) est donnée en regard de chaque
tronçon. Quel est le poids du véhicule le plus lourd pouvant aller de A à B ?
17
33 87
1 5 47
6 B
20
49 55
28 42
A 3
49
50
e
32
2
48 85
br
4 8
em
59
65 7 34
pt
3
0,729 0,81
1
g
1 6
0,81 0,9 0,729 0,81
Al
s 0,9 4 0,9 8
0,81
0,729
T_
1 0,9
2 0,9 7
0,729 0,81
5
17
C 4.5 6 5.4 4
D 5.5 4.5 5 3.8
20
trois unités de production notées 1, 2 et 3 et veut créer trois nouvelles unités,
notées 4, 5 et 6. Pour cela, elle a retenu trois emplacements notés a, b et c. Une
représentation de la disposition des unités de production est donnée dans la
figure 2.22. Les pièces produites doivent être transférées d’une unité existante
e
vers une nouvelle unité en utilisant uniquement les connexions représentées
6 a
br
sur la figure. Par exemple, la distance entre les unités 1 et b est de 3 + 2 = 5.
em
5 3
4 2
pt
3 b
se
1 1 c
25
1 2 3 4 5 6
4 5 6
1 9 4 0
2 2 7 5
3 3 0 10
Quelle nouvelle unité (4, 5 ou 6) placer en a, (respectivement en b et en
c) de manière à minimiser la distance totale à parcourir chaque jour lors des
transferts de pièces entre les unités de production existantes et les nouvelles ?
Si les nouvelles unités doivent également transférer des pièces entre elles,
peut-on résoudre le problème avec la même technique ?
Exercice 2.7 (Examen oraux). 6 étudiants (A, . . . , F ) doivent passer des exa-
mens oraux pour différents modules (1, . . . , 5). La durée de chaque examen est
2.10 O PTIMISATION MULTIOBJECTIFS 71
d’une heure. Comment construire un horaire aussi court que possible pour le
passage des étudiants, sachant qu’un étudiant ne peut faire qu’un examen à
la fois et un professeur ne peut évaluer qu’un étudiant à la fois ?
Étudiant
A B C D E F
1 × × ×
2 × × × × ×
Module 3 × ×
17
4 × × ×
5 × ×
Exercice 2.8 (Examen écrit). Le tableau suivant résume les inscriptions aux
20
examens de divers modules par les étudiants A, . . . , M . Si chaque candidat
pas au plus un examen par jour, combien de jours sont nécessaires, au mini-
mum, pour organiser la session d’examens ? Tous les étudiants devant passer
le même module sont examinés en même temps.
e
Étudiant
1
A B C D
×
E
× ×br
F G H I
×
J
×
K L M
em
2 × × × ×
3 × × × ×
4 × ×
Module 5 × × × ×
pt
6 × ×
7 × × × ×
8 × × ×
se
tion r.
Exercice 2.11 (Configuration des touches d’un téléphone portable). On dé-
Al
sire reconfigurer les touches d’un téléphone portable. On désire placer unique-
ment les 26 lettres minuscules ainsi que l’espace sur les touches 0, 2, 3, . . . , 9.
T_
de sorte que le nombre d’arêtes ayant une extrémité dans X et l’autre dans
Y soit aussi faible que possible. Comment construire un exemple de problème
d’affectation quadratique pour la bipartition d’un graphe ? Même question pour
la transformation du TSP en QAP.
Exercice 2.13 (Bipartition spéciale). On considère un jeu de cartes numéro-
tées de 1 à 50. On désire former deux tas avec ces cartes, le premier avec une
somme des valeurs numérotées de 1170 et l’autre avec un produit de 36000
Comment coder une tentative de solution à ce problème ? Comment évaluer la
qualité d’une tentative de solution ?
17
Exercice 2.14 (Carré magique). On désire créer un carré magique d’ordre n.
Ce carré comporte n × n cases devant être remplies avec les nombres de 1 à
20
n2 . La somme des nombres dans chaque ligne, colonne et diagonale doit valoir
(n3 − n)/2. On donne ci-dessous un exemple de carré magique d’ordre 4.
34
e
%
4 14 15 1 → 34
9
5 br
7 6 12 →
11 10 8 →
34
34
em
16 2 3 13 → 34
↓ ↓ ↓ ↓ &
34 34 34 34 34
Comment coder une tentative de solution à ce problème ? Comment éva-
pt
traitement à séquence fixée, comment calculer les heures de fin au plus tôt fij
et de début au plus tard dij de chaque opération si les objets passent dans un
ordre donné par une permutation p ?
g
différentes à produire (dans un ordre qui peut être choisi par le gestionnaire de
la chaîne) et que l’on connaît le temps tij du traitement de la plaque i sur la
machine j. De plus, lorsqu’une machine a terminé le traitement d’une plaque,
cette dernière doit immédiatement passer sur la machine suivante, qui doit la
traiter sans attendre pour qu’elle ne refroidisse pas. Une machine ne traite
qu’une plaque à la fois. Le gestionnaire de la chaîne doit déterminer dans quel
ordre produire les n plaques pour terminer la production aussi rapidement que
possible. Comment modéliser ce problème sous la forme d’un voyageur de
commerce ?
17
20
Deuxième partie
e
br
Techniques heuristiques de
em
base
pt
se
25
g
Al
T_
T_
Al
g
25
se
pt
em
br
e
20
17
Chapitre 3
17
Méthodes constructives
20
3.1 Construction systématique
e
br
Lorsque l’on doit trouver la meilleure solution possible à un problème d’op-
timisation combinatoire, la première idée qui vient est d’essayer de construire
toutes les solutions du problème, d’évaluer leur faisabilité et leur qualité et de
em
retourner la meilleure qui satisfait l’ensemble des contraintes. Évidemment,
cette approche ne peut se prêter qu’à des problèmes de taille limitée. Prenons
l’exemple d’un tout petit problème de sac de montagne en variables 0 − 1 à
deux contraintes :
pt
(3.1)
contraintes 4x1 + 2x2 + 3x3 + 2x4 + x5 6 7
xi ∈ {0, 1}(i = 1, . . . , 5)
25
17
sommairement si un sous-problème ne peut mener à une solution meilleure
que la meilleure connue jusque là. Il s’agit de la méthode d’évaluation et sépa-
ration.
20
3.1.1 Évaluation et séparation
e
Pour estimer rapidement si un sous-problème peut avoir une solution, et si
cette dernière est prometteuse, une technique consiste à relaxer une ou plusi-
br
eurs contraintes. La solution optimale du problème relaxé n’est pas forcément
admissible pour le problème initial. Par contre, il est possible de faire quelques
em
déductions à partir de la résolution du problème relaxé : si ce dernier n’a pas
de solutions ou qu’il a une solution optimale plus mauvaise que la meilleure
solution admissible connue du problème initial, il est inutile de continuer à dé-
velopper la branche partant de ce sous-problème. Si le sous-problème relaxé
pt
a une solution optimale qui est admissible pour le problème initial, il est éga-
lement inutile de développer la branche, puisqu’on connaît déjà comment fixer
optimalement les variables encore libres de cette branche.
se
Intégralité des variables Ce qui rend le problème (3.1) difficile, c’est que l’on
25
impose aux variables de prendre des valeurs entières. On peut donc re-
laxer cette contrainte et résoudre le problème :
17
sa résolution pourrait parfois soulever certaines difficultés.
Combinaison de relaxations Il est évidemment possible de combiner plusi-
eurs types de relaxation, par exemple l’agrégation de contraintes et l’in-
20
tégralité des variables. Dans l’exemple, cela donne :
max S = 9x1 + 5x2 + 7x3 + 3x4 + x5
sous 8x1 + 5x2 + 8x3 + 4x4 + 2x5 6 17 (3.4)
e
contraintes 0 6 xi 6 1(i = 1, . . . , 5)
br
Ce problème peut être résolu en O(n log n) de la façon suivante : Les
variables sont triées dans l’ordre des ri /vi décroissants (où ri est le re-
venu de l’objet i et vi son volume (agrégé) ; dans l’exemple, les indices
em
sont déjà triés) puis on prend les objets les uns après les autres dans
cet ordre jusqu’à ce qu’il soit impossible d’ajouter entièrement un nouvel
objet (ce qui donne x1 = x2 = 1). L’objet suivant est fractionné de façon
à remplir totalement le sac ( =⇒ x3 = 4/8 pour une valeur totale du sac
pt
dant souvent géré en queue de priorité, ce qui suppose le calcul d’une valeur
Al
pour chaque sous-problème, valeur que devrait être fortement corrélée avec
celle de la meilleure solution pouvant être obtenue en développant la branche
(par exemple, la valeur S du problème relaxé). Souvent, le choix d’une méthode
T_
17
Algorithme 3.1 : Trame d’une méthode d’évaluation et séparation pour un ob-
jectif f à maximiser. Il est nécessaire de fournir trois méthodes : une pour la ges-
tion des sous-problèmes à résoudre (généralement, une queue de priorité (basée
20
sur un critère heuristique) ou une pile), une méthode de calcul de relaxation des
sous-problèmes et une heuristique choisissant la prochaine variable sur laquelle
séparer un problème en plusieurs sous-problèmes.
Entrées : Un problème à n variables x1 , . . . , xn , politique α de gestion
e
des sous-problèmes, méthode β de relaxation, politique γ de
1 f ∗ ← −∞;
branchement
br
Résultat : Une solution optimale x∗ de valeur f ∗
// Valeur de la meilleure solution trouvée
em
2 F ←∅; // Ensemble des variables fixées
3 L ← {x1 , . . . , xn } ; // Ensemble des variables libres
4 Q ← {(F, L)} ; // Ensemble des sous-problèmes à résoudre
5 tant que Q 6= ∅ faire
pt
10 x∗ ← x;
11 f ∗ ← f (x)
12 sinon si f (x) > f ∗ alors
13 Choisir un xk ∈ L selon la politique γ; // Développer la branche
g
16
17 sinon
Élaguer la branche ; // Elle n’a pas de solution meilleure que x∗
T_
18
3.2 C ONSTRUCTION ALÉATOIRE 79
17
riable xk à fixer est très grand ; notamment lorsque xk peut prendre n’im-
porte quelle valeur entière. Une technique de branchement est de considérer
la valeur y non entière d’une variable xk pour le problème relaxé et de créer
20
deux branches, l’une avec la contrainte supplémentaire xk 6 byc et l’autre
xk > byc + 1. Dans ce cas, les ensembles de variables fixées et libres ne sont
pas modifiés à la ligne 16. Cette technique a été proposée par [?].
Ces dernières années, les méthodes dites exactes pour la résolution de
e
programmes linéaires en nombres entiers ont fait de gros progrès, notamment
grâce au développement de relaxations, de politiques de séparation et de bran-
br
chement de plus en plus élaborées. Certains logiciels, comme CPLEX ou Gu-
robi, intègrent des méthodes basées sur des métaheuritiques pour le calcul
em
de bornes ou pour l’obtention de solutions de bonne qualité, ce qui permet
d’élaguer plus rapidement l’arbre d’énumération. Malgré cela, on n’échappe
pas à une croissance exponentielle du temps de résolution lorsque la taille des
problèmes augmente.
pt
Une méthode simple et rapide pour obtenir une solution est de la générer
aléatoirement parmi l’ensemble des solutions admissibles. On ne peut évidem-
25
ment espérer trouver ainsi de manière fiable une excellente solution, mais cette
méthode est largement employée dans des trames itératives où l’on répète une
phase constructive suivie d’une phase d’amélioration avec une recherche lo-
cale.
g
Notons qu’il n’est pas forcément évident d’écrire une procédure qui per-
mette de générer chaque solution avec la même probabilité. L’exercice 3.1
Al
Algorithme 3.2 : Trame d’une méthode gloutonne. Il ne s’agit pas d’un algo-
rithme à proprement parler, car différentes implantations sont possibles, selon la
définition de l’ensemble E des éléments constituant les solutions et la fonction de
coût incrémental.
Entrées : s solution partielle triviale (généralement ∅) ; ensemble E
17
d’éléments dont sont constituées les solutions du problème ;
fonction de coût incrémental c(s, e)
Résultat : Solution complète s
20
1 R ← E; // Éléments pouvant encore être ajoutés à s
2 tant que R 6= ∅ faire
3 ∀e ∈ R, évaluer c(s, e);
4 Choisir un e0 optimisant c(s, e0 );
e
5 s ← s ∪ e0 ; // Ajouter e0 à la solution partielle s
6
compléter s;
br
Retirer de R les éléments qui ne peuvent plus être utilisés pour
em
Des algorithmes aux performances très variées peuvent être obtenus selon
la définition de E et c(s, e). Si l’on prend l’exemple de l’arbre de Steiner, on
pourrait considérer que E est l’ensemble des arêtes du problème et la fonction
pt
de coût incrémental le poids de chaque arête. Une solution partielle est alors
constituée d’une forêt.
se
F IGURE 3.1 – Exécution d’un algorithme glouton basé sur le poids des arêtes pour le
voyageur de commerce.
17
Plus proche voisin
Une méthode gloutonne des plus simples à programmer pour le problème
20
du voyageur de commerce est celle connue sous le nom de plus proche voisin.
Les éléments à insérer sont les villes plutôt que les arêtes. Une solution par-
tielle s est donc un chemin dans lequel les villes sont visitées dans l’ordre de
leur insertion. Le coût incrémental correspond au poids de l’arête qui permet
e
de relier la prochaine ville. La figure 3.2 illustre l’exécution de cette heuristique
sur le même exemple de problème que précédemment. Le fait d’obtenir une
br
solution identique à la méthode précédente est une coïncidence. Elle peut être
programmée de façon très concise en Θ(n2 ), où n est le nombre de villes,
em
comme le montre le code 3.1.
pt
se
F IGURE 3.2 – Exécution du plus proche voisin pour le problème du voyageur de com-
merce.
25
double t s p _ n e a r e s t _ n e i g h b o u r ( i n t n , / ∗ Number o f c i t i e s ∗/
double ∗∗d , /∗ Distance matrix ∗/
int s, /∗ S t a r t i n g c i t y ∗/
Al
f o r ( i = 0 ; i < n ; i ++)
solution [ i ] = i ;
swap ( s o l u t i o n + s , s o l u t i o n + 0 ) ; / ∗ Tour s t a r t s from c i t y s ∗ /
17
L’heuristique du regret le plus grand choisit donc la ville e qui maximise
c(s, e) = minj,k∈R dje + dek − minr∈R die + der .
20
Insertion la moins coûteuse
e
dans une tournée partielle. L’ensemble E est donc celui des villes et la solution
triviale de départ un cycle sur les deux villes les plus proches. Le coût incré-
br
mental c(s, e) d’une ville est le détour minimal qu’il faut consentir pour insérer
la ville e dans la tournée partielle s entre deux villes successives de s.
em
pt
Nous n’avons donné ici qu’une toute petite palette des méthodes construc-
tives gloutonnes qui ont été proposées pour le voyageur de commerce. La
3.4 A MÉLIORATION DE LA CONSTRUCTION GLOUTONNE 83
qualité des solutions qu’elles produisent est très variable. Il n’est généralement
pas très difficile de trouver des exemples de problèmes pour lesquels l’heuris-
tique gloutonne se fourvoie et fait des choix de plus en plus mauvais. Sur des
semis de point uniformément distribués dans un carré du plan euclidien, elles
fournissent des solutions dont la longueur dépasse de quelques dizaines de
pour-cent celle de l’optimum.
17
Méthode DSATUR
20
mets d’un graphe essaie de déterminer le sommet pour lequel l’attribution
d’une couleur risque d’être la plus problématique. La méthode DSATUR [?]
fait l’hypothèse que ce sera le sommet qui a déjà de nombreux voisins colorés
avec une large palette de couleur. On définit donc le degré de saturation d’un
e
sommet v, noté DS(v), le nombre de couleurs différentes employées par les
sommets adjacents de v. À degré de saturation égal — notamment au départ,
br
lorsqu’aucun sommet n’est coloré — on choisira le sommet de degré le plus
élevé. À degré de saturation et degré égaux, les sommets sont départagés
em
arbitrairement.
11 R ← R \ v0 ;
17
F IGURE 3.5 – Recherche en faisceaux avec p = 3 et k = 3. Avant de choisir l’élément à
insérer dans la solution partielle, on fait une énumération en largeur jusqu’à une profon-
deur de k, en ne retenant à chaque profondeur que les p candidats jugés les meilleurs.
20
donc limiter le développement de certaines branches, ce qui se fait générale-
ment sur la base de critères gloutons.
e
Dans cette section, nous présentons deux techniques d’énumération par-
br
tielles qui ont été proposées pour améliorer un algorithme glouton. La pre-
mière, la recherche en faisceaux (Beam Search dans la littérature anglo-
saxonne) a été proposée dans le cadre d’une application en reconnaissance
em
de la parole [?]. La seconde, proposée par [?] est la méthode pilote, plus ré-
cente, a été présentée directement comme la trame d’une méta-heuristique.
D’autres trames en ont été dérivées [?].
pt
mais en faisant une énumération complète jusqu’à une profondeur telle que
l’on ait plus de p nœuds au dernier niveau. De l’ensemble de ces nœuds, on
T_
en retient les p meilleurs qui sont utilisés pour générer les candidats au niveau
suivant.
17
20
Meilleure solution complète
F IGURE 3.6 – Méthode pilote. On ajoute à la solution de départ un élément, puis on
applique une heuristique pilote qui la complète totalement avant de répéter le processus
avec un autre élément ajouté à la solution de départ. L’élément finalement inséré est
e
celui ayant mené à la meilleure des solutions complètes.
br
em
Cette trame énumère toutes les solutions partielles qu’il est possible d’ob-
tenir en ajoutant un élément à la solution de départ. On applique ensuite à
toute ces solutions partielles l’heuristique pilote pour aboutir à autant de solu-
tion complètes. La solution partielle qui a été à l’origine de la meilleure solution
pt
complète est utilisée comme nouvelle solution de départ, jusqu’à ce qu’il n’y ait
plus d’élément à ajouter. La figure 3.6 illustre une étape de la méthode.
se
25
g
Al
T_
86 M ÉTHODES CONSTRUCTIVES
17
4 pour tous les e ∈ R faire
5 Compléter sp avec e pour obtenir se ;
6 Appliquer h(se ) pour obtenir une solution complète s;
20
7 si f (s) 6 v alors
8 v ← f (s);
9 sr ← se ;
10 si s est meilleure que s∗ alors
e
11 s∗ ← s; // Mémoriser la solution complète s
12
13
sp ← sr ;
br
// Ajouter un élément à la solution partielle sp
Retirer de R les éléments qui ne peuvent plus être utilisés pour
compléter sp ;
em
La trame 3.4 précise le fonctionnement de cette méta-heuristique. Notons
que dans cette trame, la dernière des solutions «partielles» est une solution
pt
complète admissible qui n’est pas forcément la solution retournée par cet algo-
rithme. En effet, l’heuristique pilote peut générer une solution complète qui ne
se
17
f o r ( q = 0 ; q < n−1; q++) / ∗ c i t i e s up t o q a t t h e i r f i n a l pos . i n s o l u t i o n ∗ /
{ / ∗ Choosing t h e n e x t c i t y r t o i n s e r t a t p o s i t i o n q ∗ /
cost_ins_r = i n f i n i t e ;
f o r ( r = q ; r < n ; r ++)
20
{ f o r ( i = 0 ; i < n ; i ++)
s [ i ] = solution [ i ] ;
swap (& s [ q ] , &s [ r ] ) ;
/ ∗ n e a r e s t neighbour on c i t i e s from q+1 ∗ /
f o r ( i = q +1; i < n−1; i = i + 1 )
{ dist_ins = i n f i n i t e ;
e
f o r ( j = i ; j < n ; j = j +1)
i f ( d [ s [ i −1]][ s [ j ] ] < d i s t _ i n s )
br
{ d i s t _ i n s = d [ s [ i −1]][ s [ j ] ] ;
j_kept = j ;
}
swap(& s [ i ] , &s [ j _ k e p t ] ) ;
em
} /∗ f o r i ∗/
cost = tsp_length ( n , d , s ) ;
t s p _ 2 o p t _ f i r s t ( n , d , s , &c o s t ) ;
i f ( cost_ins_r > cost )
{ cost_ins_r = cost ;
r_kept = r ;
pt
}
} /∗ f o r r ∗/
swap(& s o l u t i o n [ r _ k e p t ] , & s o l u t i o n [ q ] ) ;
} /∗ f o r q ∗/
se
free ( s ) ;
cost = tsp_length ( n , d , s o l u t i o n ) ;
t s p _ 2 o p t _ f i r s t ( n , d , s o l u t i o n , &c o s t ) ;
return cost ;
}
25
g
Al
T_
88 M ÉTHODES CONSTRUCTIVES
3.5 Exercices
Exercice 3.1. Permutation aléatoire Écrire une procédure pour générer une
permutation aléatoire de n éléments contenus dans un tableau p. On désire
que la probabilité de trouver l’élément pi en j ième position dans p après permu-
tation soit 1/n. Décrire les défaut des algorithmes 3.5 et 3.6.
17
Entrées : Un ensemble de n éléments e1 , . . . , en
Résultat : Une permutation p des éléments
1 i ← 0; // Nombre d’éléments choisis
20
2 tant que i 6= n faire
3 Tirer un nombre aléatoire u uniformément entre 1 et n;
4 si eu n’a pas déjà été choisi alors
5 i ← i + 1;
e
6 p i ← eu
br
em
Algorithme 3.6 : Autre mauvais algorithme pour générer une permutation aléa-
toire de n éléments
Entrées : Un ensemble de n éléments e1 , . . . , en
Résultat : Une permutation p des éléments
pt
7 pi ← eu0
8 sinon
9 p i ← eu
g
Al
Exercice 3.4. Voyageur de commerce avec limitation des voisins Pour ac-
célérer une méthode gloutonne pour le TSP, on ne considère que les arêtes
des 40 plus proches voisins de chaque sommets. Est-ce susceptible de dimi-
nuer la complexité algorithmique de la méthode ? Est-ce que cela peut soulever
certains problèmes ?
3.5 E XERCICES 89
17
blèmes de la littérature.
Exercice 3.7. Méthodes gloutonnes pour l’élaboration de tournées Propo-
20
ser deux heuristiques gloutonnes pour le problème de l’élaboration de tournées
de véhicules.
Exercice 3.8. Heuristique pilote avec recherche locale Le code 3.2 fait in-
e
tervenir comme heuristique pilote une complétion de la solution partielle avec
l’heuristique du plus proche voisin (voir l’algorithme 3.1, suivie d’une amélio-
br
ration de la solution complète avec une recherche locale tsp_2opt_first.
Est-ce que la suppression de l’optimisation de la solution complète qui a été
construite, juste avant le return final a une grande incidence sur la qualité de
em
la solution retournée par la méthode ? Est-ce qu’une méthode pilote où l’heu-
ristique pilote est simplement l’heuristique du plus proche voisin (qui peut être
obtenue par simple suppression des deux appels à tsp_2opt_first, bien que
pt
17
20
e
br
em
pt
se
25
g
Al
T_
Chapitre 4
17
Améliorations locales
20
En examinant les solutions produites par des heuristiques constructives
e
gloutonnes comme celles présentées au chapitre précédent, on remarque im-
médiatement qu’elles ne sont pas optimales. Par exemple, la solution du pro-
br
blème du voyageur de commerce euclidien obtenue par l’heuristique du plus
proche voisin présentée en figure 3.2 comporte des arêtes qui se croisent,
ce qui est évidemment sous-optimal. En effet, il est possible de remplacer les
em
deux arêtes qui se croisent par deux autres dont la somme des longueurs est
moindre tout en conservant une tournée. Ce processus peut ensuite être ré-
itéré jusqu’à l’obtention d’une solution qui ne peut plus être améliorée par le
même procédé, comme illustré en figure 4.1
pt
se
4.1 Trame
L’idée générale est donc de partir d’une solution obtenue à l’aide d’une mé-
25
mettant de trouver une direction d’amélioration n’existe pas ; elle est remplacée
Al
T_
17
2 Chercher une modification de s en s0 qui l’améliore;
3 si s0 est meilleure que s alors
4 s ← s0 ;
20
5 jusqu’à ce qu’on ne trouve plus de modifications améliorant s;
e
4.2 Voisinage
br
Pour pouvoir appliquer la trame 4.1 d’amélioration, il faut définir comment
em
obtenir les solutions voisines. Formellement, on doit définir, pour toute solution
s ∈ S, un ensemble V (s) ⊂ S. Dans la trame, la recherche d’une modification
de s passe donc par l’énumération des solutions de V (s) pour en extraire une,
s0 , qui soit meilleure que s.
pt
Une manière commode pour définir un voisinage V (s) est de spécifier les
modifications, appelés communément les mouvements que l’on peut apporter
se
Ce voisinage, consistant à remplacer deux arêtes par deux autres est connu
dans la littérature sous le nom de 2-opt. L’ensemble M (s) des mouvements de
type 2-opt que l’on peut appliquer à la solution s peut être défini formellement
par M (s) = {[i, j], i, j ∈ s, i 6= j, j 6= si , i 6= sj }.
g
17
4.2.2 Politique du meilleur mouvement
20
Dans cette politique, on examine entièrement le voisinage à chaque itéra-
tion. La solution voisine choisie comme solution courante pour l’itération sui-
vante est la meilleure identifiée.
e
Algorithme 4.3 : Trame d’une méthode d’amélioration avec politique du meilleur
mouvement.
br
Entrées : Solution s, spécification du voisinage V (·), et de la
em
fonction-objectif f (·) à minimiser.
Résultat : Solution s améliorée
1 répéter
2 f in ← vrai;
pt
3 valeur_meilleur_voisin ← ∞;
4 pour tous les s0 ∈ V (s) faire
si f (s0 ) < valeur_meilleur_voisin alors
se
5
6 valeur_meilleur_voisin ← f (s0 );
7 meilleur_voisin ← s0
si valeur_meilleur_voisin < f (s) alors
25
8
9 s ← s0 ;
10 f in ← faux
11 jusqu’à f in;
g
Al
Parfois, l’ensemble V (s) est si grand que son énumération se fait soit im-
plicitement pour en extraire la meilleure solution voisine, soit de façon heuris-
T_
104
Temps de calcul
103
102
10
1
10−1 2-opt best
10−2 2-opt first
17
10−3 Θ(n3.25 )
- −4
10 Θ(n2.33 )
−5
10
20
10 102 103 104 105
Taille du problème
e
F IGURE 4.2 – Évolution du temps de calcul en fonction du nombre de villes d’un pro-
blème de voyageur de commerce pour deux variantes de méthodes d’amélioration.
br
em
chemin optimum avec l’algorithme 2.4 de Bellman-Ford et la programmation
linéaire avec l’algorithme du Simplexe.
Il faut insister sur le fait qu’un optimum local relativement à un voisinage
n’en est pas forcément un pour un autre voisinage.
pt
Comme l’ensemble des solutions d’un problème combinatoire est fini et que
les trames 4.2 et 4.3 n’effectuent une modification que si la solution est stric-
se
ment et les exemples de problèmes euclidiens ont été tirés de la TSPLIB [?].
Comme les échelles sont logarithmiques, des temps de calculs approximative-
Al
TSP 3-opt
17
effet deux solutions voisines de même coût pourraient être considérées toutes
deux meilleures l’une que l’autre, suite à des imprécisions numériques. En ac-
ceptant toutes les solutions dont la différence de coût est fautivement négative,
20
on s’exposerait à un code sans fin. En choisissant correctement la valeur de ,
on s’affranchit de ce problème. En pratique, on peut choisir = 10−14 · max dij ,
la limite de la précision relative du type double étant approximativement 10−15
(voir code A.4).
e
br
Listing 4.1 – tsp_ 3opt_ first.c Implantation en C d’une méthode d’amélioration basée
sur le voisinage 3-opt. En pratique sa complexité en Ω(n3 ) la rend difficilement utilisable
pour des problèmes présentant plusieurs centaines de villes.
em
void t s p _ 3 o p t _ f i r s t ( i n t n , / ∗ Number o f c i t i e s ∗/
double ∗∗d , / ∗ D i s t a n c e m a t r i x , symmetry n o t r e q u i r e d ∗/
i n t best_sol [ ] , / ∗ InOut S o l u t i o n p r o v i d e d and r e t u r n e d ∗/
double ∗ b e s t _ c o s t ) / ∗ InOut Length o f t h e b e s t t o u r ∗/
pt
{
i n t ∗s = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ; / ∗ Successors o f c i t i e s ∗/
i n t i , j , k , temp , / ∗ indexes ∗/
/ ∗ L a s t v a l u e s o f indexes ∗/
se
i = 0; j = s [ i ] ; k = s [ j ] ;
i_end = i ; j_end = j ; k_end = k ;
do
{ delta = d [ i ] [ s [ j ] ] + d [ j ] [ s [ k ] ] + d [ k ] [ s [ i ] ] / ∗ Computation o f move c o s t ∗ /
−d [ i ] [ s [ i ] ] − d [ j ] [ s [ j ] ] − d [ k ] [ s [ k ] ] ;
g
i f ( d e l t a < −e p s i l o n ) / ∗ I m p r o v i n g move ? ∗ /
{ temp = s [ i ] ; s [ i ] = s [ j ] ; s [ j ] = s [ k ] ; s [ k ] = temp ; / ∗ Perform move ∗ /
temp = j ; j = k ; k = temp ; / ∗ Replace j between i and k ∗ /
Al
k = s[k]; / ∗ Next k ∗ /
i f ( k == i )
{ j = s[ j ]; k = s[ j ];} / ∗ k a t i t s l a s t value , n e x t j ∗ /
i f ( k == i )
{ i = s[ i ]; j = s[ i ]; k = s[ j ];} / ∗ j a t i t s l a s t value , n e x t i ∗ /
} / ∗ A complete t o u r performed w i t h o u t improvement ∗ /
w h i l e ( i ! = i_end | | j ! = j_end | | k ! = k_end ) ;
/ ∗ R e b u i l d s o l u t i o n under t h e form o f a p e r m u t a t i o n o f c i t i e s ∗ /
f o r ( j = 0 , i = 0 ; i < n ; j = s [ j ] , i ++)
best_sol [ i ] = j ;
free ( s ) ;
}
96 A MÉLIORATIONS LOCALES
17
F IGURE 4.3 – Mouvement de type 3-opt où trois arcs sont remplacés par trois autres
arcs. Les autres chemins de la solution sont parcourus dans le même sens, avant et
20
après modification. Une autre façon d’appréhender le voisinage 3-opt est le déplace-
ment d’un sous-chemin ailleurs dans la tournée.
e
br
em
pt
se
F IGURE 4.4 – Mouvement de type Or-opt où trois sommets successifs sont déplacés à
l’intérieur de la tournée. Dans certains cas, il est souhaitable de modifier l’ordre respectif
de visite des trois sommets.
25
g
TSP Or-opt
Al
Un autre type de voisinage proposé par [?] consiste à raisonner sur le dé-
placement de villes à l’intérieur de la tournée. L’idée est d’examiner s’il est
T_
possible de placer ailleurs trois villes qui sont visitées successivement dans la
solution courante. L’originalité de la méthode proposée par Or est de ne pas
se contenter d’un seul voisinage : Une fois qu’il n’est plus possible d’améliorer
la solution en déplaçant trois villes, on cherche à n’en déplacer que deux. Dès
qu’un déplacement de deux villes permet d’améliorer la solution, on revient
au premier voisinage. Ensuite, lorsqu’on a une solution localement optimale
relativement à ces deux voisinages, on essaie de ne déplacer qu’une ville.
Le voisinage est différent d’une limitation du voisinage 3-opt dans lequel un
sous-chemin serait limité à 3, 2 ou 1 ville. En effet, il est possible d’inverser
l’ordre de parcours du chemin avec ce voisinage.
Il est possible de tester en Θ(n2 ) si une solution est Or-optimale.
4.2 VOISINAGE 97
17
20
F IGURE 4.5 – Structure de données permettant d’effectuer un mouvement de type 2-
opt en temps constant. Les entrées 2i et 2i + 1 d’un tableau donnent les indices du
même tableau identifiant les deux villes adjacentes à la ville i. En partant de l’indice
e
0, on peut reconstituer la tournée en allant de ville adjacente en ville adjacente et en
partant de l’indice 1, on peut reconstituer la tournée dans l’autre sens. La réalisation
d’un mouvement 2-opt s’effectue par modification de 4 entrées du tableau.
br
em
Structure de données pour TSP 2-opt
villes à visiter.
void build_2opt_data_structure ( i n t n ,
const i n t s o l u t i o n [ ] ,
int t [])
{ int i ;
/ ∗ Forward t o u r ∗ /
f o r ( i = 0 ; i < n−1; i = i + 1 )
t [2∗ s o l u t i o n [ i ] ] = 2∗ s o l u t i o n [ i + 1 ] ;
t [2∗ s o l u t i o n [ n −1]] = 2∗ s o l u t i o n [ 0 ] ;
/ ∗ Backward t o u r ∗ /
f o r ( i = 1; i < n ; i = i + 1)
t [2∗ s o l u t i o n [ i ] + 1 ] = 2∗ s o l u t i o n [ i −1]+1;
t [2∗ s o l u t i o n [ 0 ] + 1 ] = 2∗ s o l u t i o n [ n−1]+1;
}
98 A MÉLIORATIONS LOCALES
17
Listing 4.3 – tsp_ 2opt_ first.c Implantation en C d’une méthode d’amélioration pour
le voyageur de commerce basée sur le voisinage 2-opt. Cette implantation utilise la
structure de donnée présentée en figure 4.5 permettant d’effectuer un mouvement en
20
temps constant.
void t s p _ 2 o p t _ f i r s t ( i n t n , / ∗ Number o f c i t i e s ∗/
double ∗∗d , / ∗ D i s t a n c e m a t r i x , must be s y m m e t r i c a l ∗/
i n t best_sol [ ] , / ∗ InOut S o l u t i o n p r o v i d e d and r e t u r n e d ∗/
double ∗ b e s t _ c o s t ) / ∗ InOut Length o f t h e b e s t t o u r ∗/
e
{ int ∗t ; / ∗ Doubly connected edge l i s t f o r e f f i c i e n t l y p e r f o r m i n g moves ∗ /
i n t i , j , next_i , next_j , l a s t _ i ;
double d e l t a ;
br /∗ Indices ∗/
/ ∗ Cost d i f f e r e n c e ∗ /
em
t = ( i n t ∗) m a l l o c ( ( s i z e _ t ) ( 2 ∗ n ) ∗ s i z e o f ( i n t ) ) ;
build_2opt_data_structure ( n , best_sol , t ) ;
i = l a s t _ i = 0; / ∗ S t a r t moves e v a l u a t i o n from c i t y 0 ∗ /
w h i l e ( t [ t [ i ]] > >1 ! = l a s t _ i ) / ∗ Index i has made 1 t u r n w i t h o u t impovement ∗ /
{ j = t [ t [ i ]];
pt
} /∗ while i ∗/
/ ∗ R e b u i l d s o l u t i o n : number o f c i t y i n j t h p o s i t i o n ∗ /
Al
f o r ( i = 0 , j = 0 ; j < n ; i = t [ i ] , j ++)
b e s t _ s o l [ j ] = i > >1;
T_
free ( t ) ;
} / ∗ t s p 2−o p t ∗ /
Connexité
Faible diamètre
17
Un voisinage devrait permettre d’atteindre une solution optimale en peu
d’étapes. Par définition, le diamètre d’un graphe est la longueur maximale d’un
20
plus court chemin reliant deux sommets.
Faible rugosité
e
Un voisinage devrait posséder un nombre d’optima locaux aussi faible que
br
possible et une forte corrélation entre les valeur des solutions voisines. L’idéal
serait qu’il n’en possède qu’un, qui serait alors l’optimum global, atteignable
em
à tous les coups en partant de n’importe quelle solution. Évidemment, cette
propriété n’est pas remplie pour les problèmes difficiles. Néanmoins, trouver
un voisinage adéquat pour le problème considéré est essentiel pour le succès
d’une méthode d’amélioration.
pt
2 2
(x+ x + τ )/2, où τ est un paramètre que l’on fait tendre vers 0 si l’on veut se
Al
Faible taille
Facilité d’évaluation
17
4.3 Limitation du voisinage
20
Typiquement, la taille d’un voisinage croît de façon quadratique ou cubique
avec la taille du problème. Par conséquent, une recherche locale nécessitant
e
d’examiner de nombreuses fois l’ensemble du voisinage requiert un effort de
calcul prohibitif lorsque la taille du problème augmente. Diverses techniques
br
ont été proposées pour limiter la croissance des calculs.
em
4.3.1 Liste de candidats
Une première idée consiste à faire l’hypothèse qu’un mouvement favorable
pour une solution restera bon pour des solutions pas trop différentes. Une
pt
été passablement modifiée ; des mouvements qui n’étaient pas favorables peu-
vent donc le devenir et vice versa.
g
17
20
e
br
em
pt
se
25
g
Al
T_
triangulation. Ceci suggère qu’une heuristique travaillant sur des solutions ba-
sées essentiellement sur les arêtes de la triangulation de Delaunay peut don-
ner d’excellent résultats tout en nécessitant un effort de calcul très limité. En
particulier, une recherche locale basée sur un voisinage de type k-opt pourra
fonctionner même pour des valeurs de k valant 5 ou 6, s’il est restreint aux
arêtes du Delaunay, puisque le nombre de solutions voisines d’un problème de
taille n sera en O(n · 6k ).
17
Limitation du voisinage pour le voyageur de commerce avec 1-arbres
20
technique plus générale de limitation du voisinage pour le problème du voya-
geur de commerce, proposée par [?], est basée sur le concept de 1-arbre, qui
n’est rien d’autre qu’une forme de relaxation lagrangienne.
Un 1-arbre dans un réseau comportant des sommets 1, 2, . . . , n est un arbre
e
de poids minimum sur les sommets 2, . . . , n auquel on ajoute les deux arêtes
br
adjacentes au sommet 1 qui ont le poids le plus faible. Le nœud 1 joue ici
un rôle particulier, d’où le terme de 1-arbre. On peut reformuler le problème
du voyageur de commerce en spécifiant que l’on cherche un 1-arbre de poids
em
minimum pour lequel le degré de tous les sommets est 2.
P
min z = d x
Pn ij ij
(i,j)∈H
pt
P Pn Pn
min z(λ) = (i,j)∈H dij xij + i=1 λi ( j=1 xij − 2)
(4.2)
sous contrainte H est un 1-arbre
minimum dans un réseau dont le poids de l’arête (i, j) est modifiée en dij +
Al
17
rieur à 1 contrôlant ce qui a été appelé la granularité de la recherche locale.
Pour le problème de distribution, il est toutefois nécessaire de conserver les
longs trajets directs que l’on pourrait avoir entre clients et dépôt.
20
Une technique similaire a été abondamment utilisée pour le problème du
voyageur de commerce : au lieu de considérer un graphe complet où chaque
ville est reliée à toutes les autres, on se contente d’un graphe où chaque ville
n’est connectée qu’à ses p plus proches voisines, avec p limité à quelques
e
dizaines. Ainsi, la taille d’un voisinage 2-opt est en n · p2 au lieu de n2 . La
br
perte sur la qualité des solutions obtenues avec un tel voisinage réduit est
souvent négligeable. Par contre, l’implantation de cette idée n’est pas triviale :
tout d’abord, il se peut que le graphe où chaque nœud n’est relié qu’à ses p
em
plus proches voisins ne soit pas connexe. Il faut donc ajouter des connexions
plus longues pour qu’il contienne un cycle passant par tous les nœuds. Ensuite,
l’implantation de la recherche locale est plus complexe. Par exemple, pour le
voisinage 2-opt, on ne peut plus exploiter directement la structure de données
pt
i et une ville j proche de i (il n’y en a plus que p), on ne peut plus identifier
la ville sj qui succède à j en parcourant la tournées dans le sens i → si . De
la même manière, pour le voisinage 3-opt, on peut rapidement détecter trois
villes i, j et k proches qui seraient candidates pour un mouvement 3-opt, mais
25
17
Meilleure solution visitée : solution de départ pour la prochaine itération
20
celles énumérées.
e
k modifications successives : Ve (s) = {s0 |s0 = s ⊕ m1 ⊕ · · · ⊕ mk , m1 , . . . mk ∈
M (s)} La taille de Ve croît exponentiellement avec k. Pour éviter une telle crois-
br
sance, on peut utiliser la stratégie de la recherche en faisceaux présentée à la
section 3.4.1, mais adaptée à une recherche locale plutôt qu’à une méthode
em
constructive. Cette technique a été appelée recherche en éventail avec filtre. À
chaque niveau, on ne retient que les p meilleures solutions voisines (qui peu-
vent être moins bonnes que la solution de départ) et on évalue les voisins de
ces dernière avant de réitérer le processus jusqu’au niveau k. Ainsi, à chaque
pt
Tournée initiale
a d c b
Chaîne
Structure de référence
17
20
Nouvelle chaîne
a0 d0 c0 b0
e
Nouvelle structure
br
em
Nouvelle tournée
pt
arête [b, d] pour obtenir une structure de référence qui peut être transformée soit en une
autre structure de référence par éjection de l’arête [c, d] et ajout d’une autre arête, soit
en une tournée par éjection de la même arête et ajout de l’arête [a, c].
25
de référence.
T_
17
Solution de départ Chaîne Structure 1
20
e
br
em
Structure 2 Structure 3 Solution acceptée
pt
F IGURE 4.9 – Application d’une chaîne d’éjection sur un petit problème de voyageur de
commerce. Après 4 éjections, on réussit à améliorer la solution de départ.
se
mécanisme peut sembler artificiel dans cette figure, mais avec des villes bien
25
réparties sur le plan, on se rend compte qu’il est possible d’obtenir des modi-
fications relativement complexes et améliorant des solutions pas franchement
mauvaises, comme illustré en figure 4.9.
g
Al
T_
4.4 E XTENSION DU VOISINAGE 107
17
Algorithme 4.4 : Trame d’une méthode d’amélioration par chaîne d’éjections
20
pour le problème du voyageur de commerce.
Entrées : Solution p d’un voyageur de commerce
Résultat : Solution p améliorée
1 répéter
e
2 Initier la chaîne par l’éjection d’une arête [a, b];
3 répéter
4
br
Trouver l’arête [b, d] qui, une fois ajoutée, minimise le poids de la
structure de référence et telle que l’arête [c, d] n’ait pas encore été
em
supprimée dans la chaîne d’éjections courante;
5 si L’arête [b, d] a été trouvée et Le poids de la structure de
référence est inférieure au poids de la tournée p alors
6 si L’ajout de [a, c] et la suppression de [c, d] améliore la
pt
solution alors
7 Succès : Remplacer la solution p par celle qui a été
trouvée;
se
8 sinon
9 Ajouter l’arête [b, d];
10 Supprimer l’arête [c, d];
25
11 b ← c;
12 sinon
13 Échec de la chaîne : revenir à la solution p et tenter une autre
g
éjection
14 jusqu’à ce qu’aucune arête [b, d] convenable existe;
Al
15 jusqu’à ce que toutes les arêtes de la solution p ont initié une chaîne
sans succès;
T_
108 A MÉLIORATIONS LOCALES
17
f o r ( i = 0 ; i < n ; i ++) tabu [ i ] = ( i n t ∗) c a l l o c ( ( s i z e _ t ) n , s i z e o f ( i n t ) ) ;
f o r ( i = 0 ; i < n ; i ++) s [ b e s t _ s o l [ i ] ] = b e s t _ s o l [ ( i +1)%n ] ;
no_iter = 0;
a = end_a = 0 ; / ∗ I n i t i a t e e j e c t i o n c h a i n from c i t y 0 ∗/
20
do
{ b = s[a];
length_path = ∗best_cost − D[ a ] [ b ] ;
n o _ i t e r ++;
improved = 0 ;
do / ∗ I d e n t i f y b e s t r e f e r e n c e s t r u c t u r e w i t h edge ( a−b ) removed ∗/
e
{ path_modified = 0;
cost_ref_struct = i n f i n i t e ;
best_c = c = s [ b ] ;
while ( s [ c ] != a )
{ d = s[c];
{
br / ∗ E j e c t i o n can be propagated
∗/
em
best_c = c ;
c o s t _ r e f _ s t r u c t = length_path − D[ c ] [ d ] + D[ c ] [ a ] + D[ b ] [ d ] ;
break ;
}
i f ( tabu [ c ] [ d ] ! = n o _ i t e r && / ∗ edge c − d has n o t been added ∗/
length_path + D[ b ] [ d ] < c o s t _ r e f _ s t r u c t )
pt
{ c o s t _ r e f _ s t r u c t = length_path + D[ b ] [ d ] ;
best_c = c ;
}
se
c = d; / ∗ Next v a l u e f o r d ( and c ) ∗/
}
i f ( c o s t _ r e f _ s t r u c t + e p s i l o n < ∗ b e s t _ c o s t ) / ∗ Adm. r e f . s t r u c t . found ∗/
{ path_modified = 1;
c = best_c ;
d = s[c]; / ∗ Update r e f e r e n c e s t r u c t u r e ∗/
25
{ ssi = s [ si ] ; s [ si ] = i ; i = si ; si = ssi ; }
b = c;
Al
improved = 1 ;
f o r ( i = 0 , j = 0 ; j < n ; j ++) / ∗ S t o r e improved s o l u t i o n ∗/
{ best_sol [ j ] = i ; i = s [ i ] ; }
}
}
}
while ( path_modified ) ;
f o r ( i = 0 ; i < n ; i ++) s [ b e s t _ s o l [ i ] ] = b e s t _ s o l [ ( i +1)%n ] ;
a = s[a];
}
w h i l e ( a ! = end_a | | improved ) ; / ∗ a has made 1 t u r n w i t h o u t improvement ∗/
f o r ( i = 0 ; i < n ; i ++) f r e e ( tabu [ i ] ) ;
f r e e ( tabu ) ;
free ( s ) ;
} / ∗ t s p LK ∗ /
4.5 U TILISATION DE PLUSIEURS VOISINAGES OU PLUSIEURS MODÈLES 109
17
péter alternativement des processus d’amélioration tant que la solution trouvée
n’est pas un optimum local relativement à toutes les structures de voisinage
considérées.
20
Finalement, mentionnons que l’on peut encore passer d’une modélisation
du problème à une autre. Dans la mesure où la structure de voisinage n’est
pas la même pour les diverses modélisations, il est également possible d’itérer
des méthodes d’améliorations utilisant divers modèles. Cette technique n’est
e
pas forcément applicable telle quelle, puisqu’une solution admissible pour un
br
modèle ne l’est pas forcément pour un autre. Dans ce cas, il faut prévoir des
méthodes de réparation lorsqu’on change de modèle, ce qui peut impliquer que
le processus n’est plus une méthode d’amélioration stricte. Un corollaire est
em
que la recherche pourrait « tourner en rond », la réparation suivie de « l’amé-
lioration » avec un modèle pouvant défaire les modifications réalisées par un
autre modèle.
pt
cales
Diverses techniques relativement simple ont été proposées pour réaliser
25
17
Résultat : Approximation P de l’ensemble Pareto
1 P = s;
2 pour Imax itération faire
20
3 Générer une vecteur de pondérations → −
w;
4 répéter
5 f in ← vrai;
6 valeur_meilleur_voisin ← ∞;
e
7 pour tous les s0 ∈ V (s) faire
si s0 n’est pas dominée par les solutions de P alors
8
9
br
Insérer s0 dans P et retirer de P les solutions dominées
par s0
em
→
−
10 si →
−
w · f (s0 ) < valeur_meilleur_voisin alors
→
−
11 valeur_meilleur_voisin ← → −w · f (s0 );
12 meilleur_voisin ← s0
pt
→
−
13 si valeur_meilleur_voisin < →
−
w · f (s) alors
14 s ← s0 ;
se
15 f in ← faux
16 jusqu’à f in;
25
Une autre approche basée sur une recherche locale pour estimer l’en-
semble Pareto ne passe pas par une scalarisation. L’idée de [?] est de par-
Al
Pareto n’est pas stabilisé, on considère toutes les solutions voisines de celles
de l’ensemble Pareto et on met à jour ce dernier avec celles-ci. La méthode
se programme récursivement de façon très concise, comme le montre l’algo-
rithme 4.6.
Le code 4.5 implante une recherche locale Pareto pour le problème du
voyageur de commerce.
4.6 O PTIMISATION MULTIOBJECTIF PAR AMÉLIORATIONS LOCALES 111
Solution de départ
3e optimum local
17
Trajectoire de la recherche locale
20
e
Solutions ayant fait avancer le front Pareto
br
em
Front Pareto après 200 scalarisations
pt
25
2e optimum local
g
Objectif 1
Al
la technique d’agrégation des objectifs. La méthode d’amélioration est basée sur les
chaînes d’éjections et les pondérations des objectifs sont tirées aléatoirement chaque
fois que la recherche atteint un optimum local. La solution de départ est celle que l’on
obtient à l’aide d’une méthode constructive gloutonne ne travaillant que sur le premier
objectif. Il est intéressant de constater que le troisième optimum local est meilleur que
la solution de départ, pour les deux objectifs.
112 A MÉLIORATIONS LOCALES
17
4 Mise_a_jour_Pareto;
Entrées : Solution s, valeurs des objectifs →
−
v
Résultat : Ensemble Pareto P mis à jour
si (s, →
−
20
5 v ) domine une solution de P ou si P = ∅ alors
6 Retirer de P toutes les solutions dominées par (s, →
−
v );
7
→
−
P ← P ∪ (s, v );
8 Énumération_voisinage(s);
e
br
Listing 4.5 – tsp_3opt_pareto.c Code en C d’une recherche locale Pareto basée sur
les mouvements de type 3-opt pour le problème du voyageur de commerce. Cette fonc-
em
tion en utilise une autre qui met à jour l’ensemble Pareto pour chaque solution voisine
de la solution qui lui est fournie. La structure de donnée utilisée pour stocker l’ensemble
Pareto est discutée en section 4.6.3.
/ ∗ F u n c t i o n update_pareto_3opt w i l l b u i l d t h e neighbour s o l u t i o n i f necessary ∗ /
pt
/ ∗ I t r e c u r s i v e l y c a l l s t s p _ 3 o p t _ p a r e t o w i t h t h e neighbour s o l u t i o n ∗/
v o i d update_pareto_3opt (NODE ∗∗ pareto , POINT p o i n t ,
i n t n , i n t ∗s , i n t i , i n t j , i n t k , double ∗∗∗d ) ;
se
{
i n t i , j , k , dim ; / ∗ indexes ∗ /
POINT c o s t s _ n e i g h b o u r ; / ∗ Costs o f a neighbour s o l u t i o n ∗ /
i n t ∗s = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n∗ s i z e o f ( i n t ) ) ; / ∗ update_pareto a l t e r s succ ∗ /
memcpy ( s , succ , ( s i z e _ t ) n∗ s i z e o f ( i n t ) ) ; / ∗ make a copy ∗/
i = 0; j = s [ i ] ; k = s [ j ] ;
g
do
{
Al
17
Les deux techniques présentées ci-dessus pour l’optimisation multi-objectifs
peuvent être relativement peu efficaces si l’on n’utilise pas une structure de
données adéquate pour stocker l’ensemble Pareto. En effet, les évaluations
20
des objectifs d’une solution voisine peuvent souvent se faire en temps cons-
tant, soit un temps se mesurant en nano-secondes sur les ordinateurs actuels.
Si l’ensemble Pareto contient p solutions et en utilisant une simple liste pour les
stocker, vérifier si une solution voisine n’est pas dominée va prendre beaucoup
e
plus de temps et la recherche locale va en être ralentie d’un facteur proportion-
nel à p.
br
Il n’est pas rare que l’ensemble Pareto contienne des milliers de solutions.
Il est donc nécessaire d’utiliser une structure de données permettant de tester
em
la dominance d’une solution en un temps ne croissant pas linéairement avec
la taille de l’ensemble Pareto.
Casiers
pt
Arbre KD
T_
Dans le cas général, une structure de données dont le temps pour faire une
requête croît faiblement avec p est l’arbre KD. Il s’agit d’un arbre de recherche
binaire, où un nœud à la profondeur d discrimine les autres éléments de l’arbre
sur la dimension d modulo K. Le code 4.6 donne les procédures élémentaires
d’ajout d’un nouvel élément dans un arbre KD et le parcours de ce dernier.
114 A MÉLIORATIONS LOCALES
17
i f (∗ t r e e == NULL )
∗ t r e e = element ;
e l s e i f ( ( ∗ t r e e )−>key [ deepness%K ] < element−>key [ deepness%K ] )
KD_add(&(∗ t r e e )−>r , element , deepness + 1 ) ;
else
20
KD_add(&(∗ t r e e )−>L , element , deepness + 1 ) ;
}
e
int i ;
i f ( tree )
br
{ KD_scan_and_delete ( t r e e −> L , n ) ;
fo r ( i = 0; i < K; i ++)
p r i n t f ( "%g ", t r e e −>key [ i ] ) ;
em
fo r ( i = 0; i < n; i ++)
p r i n t f ( "%d ", t r e e −>i n f o [ i ] ) ;
putchar ( ’ \ n ’ ) ;
KD_scan_and_delete ( t r e e −> r , n ) ;
f r e e ( t r e e −>i n f o ) ;
pt
free ( tree ) ;
}
}
se
tombe sur une feuille, qui est simplement éliminée. Le code 4.7 implante une
Al
17
o p t = KD_find_opt (&(∗∗ t r e e ) . L , dim , deep +1 , opt , min , value , deep_opt ) ;
i f ((∗∗ tree ) . r )
o p t = KD_find_opt (&(∗∗ t r e e ) . r , dim , deep +1 , opt , min , value , deep_opt ) ;
r e t u r n opt ;
20
}
e
double v a l _ r e p l ;
i n t deep_repl ;
{
val_repl = −i n f i n i t e ;
r e p l a c i n g = KD_find_opt (&(∗∗ node ) . L , deep%K , deep +1 , NULL ,
0 , &v a l _ r e p l , &deep_repl ) ;
se
}
else / ∗ Find t h e r e p l a c i n g node w i t h minimal c o o r d i n a t e a t r i g h t ∗ /
{
val_repl = i n f i n i t e ;
r e p l a c i n g = KD_find_opt (&(∗∗ node ) . r , deep%K , deep +1 , NULL ,
1 , &v a l _ r e p l , &deep_repl ) ;
25
f o r ( i n t i = 0 ; i < K ; i ++)
(∗∗ node ) . key [ i ] = (∗∗ r e p l a c i n g ) . key [ i ] ;
KD_delete ( r e p l a c i n g , deep_repl ) ;
}
g
Al
un point situé entre le point idéal, qui est celui dont les valeurs sur toutes les di-
mensions sont celles de l’optimum des objectifs pris séparément. En pratique,
le point idéal n’est pas connu, mais on peut se contenter d’une approximation,
qui peut être aussi grossière que l’on veut, par exemple (−∞, . . . , −∞) si tous
les objectifs doivent être minimisés. Si l’arbre KD comporte un tel point, alors
celui-ci domine la solution essayée, qui est donc ignorée. Sinon, la nouvelle
solution en domine d’autres, qui doivent être éliminées, et elle est ajoutée à
l’arbre KD.
Le code 4.8 permet d’effectuer la mise à jour d’un ensemble Pareto lors-
qu’on cherche à y inclure un nouveau point.
116 A MÉLIORATIONS LOCALES
17
NODE ∗∗ r e s ;
i f (!(∗ tree ) ) / ∗ t r e e empty : no node found ∗ /
r e t u r n NULL ;
i f ( KD_in (∗ t r e e , min , max ) ) / ∗ The r o o t o f t r e e i s i n t h e box ∗ /
20
{
∗deep_found = deep ;
return tree ;
}
i f ( max [ deep%K ] > (∗ t r e e )−>key [ deep%K ] && / ∗ Look a t r i g h t sub−t r e e ∗ /
( r e s = KD_find (&(∗ t r e e )−>r , min , max , deep +1 , deep_found ) ) )
e
r e t u r n res ;
i f ( min [ deep%K ] < (∗ t r e e )−>key [ deep%K ] ) / ∗ Look a t t h e l e f t ∗ /
r e t u r n KD_find (&(∗ t r e e )−>L , min , max , deep +1 , deep_found ) ;
}
r e t u r n NULL ;
br / ∗ No node found ∗ /
element−>i n f o = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
memcpy ( element−>i n f o , s , ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
Al
t s p _ 3 o p t _ p a r e t o ( n , d , element−>i n f o , p o i n t , p a r e t o ) ;
}
}
4.7 E XERCICES 117
4.7 Exercices
Exercice 4.1. Propriétés des voisinges 2-opt et 3-opt Montrer les propriétés
suivantes des voisinages 2-opt et 3-opt :
— L’inversion de deux villes ou un mouvement de type 3-opt peut être ob-
tenu par une succession de mouvements 2-opt
— 2-opt et 3-opt sont des voisinages possédant la propriété de connectivité.
Donner une borne supérieur au diamètre de ces voisinages
17
Exercice 4.2. Limitation du voisinage 3-opt Pour limiter la taille du voisinage
d’un problème de voyageur de commerce, seules les 40 villes les plus proches
de chaque ville sont considérées. Avec une telle limitation, quelle est la com-
20
plexité de la vérification qu’une solution est 3-optimale ? A-t-on besoin d’une
structure de données particulière pour atteindre une complexité algorithmique
minimale ? Le voisinage engendré par une telle limitation est-il connexe.
e
Exercice 4.3. Voisinages pour l’élaboration de tournées de véhicules Pro-
poser 4 voisinages différents pour le problème de distribution de biens. Donner
br
la taille de ces voisinages en fonction du nombre n de clients du problème
ainsi que du nombre m de tournées d’une solution. Préciser si ces voisinages
em
disposent de la propriété de connectivité, selon la modélisation du problème.
Exercice 4.4. Voisinages pour l’arbre de Steiner Deux représentations de
solutions ont été proposées à la section 2.1.2 pour l’arbre de Steiner. Donner
des voisinages adaptés à chacune de ces modélisations du problème.
pt
cules Proposer une technique basée sur les chaînes d’éjections pour le pro-
blème de distribution de bien. Préciser comment initialiser la chaîne, comment
la propager et comment la terminer. Estimer la complexité algorithmique de
l’évaluation d’une chaîne pour une solution avec n clients et m tournées.
25
g
Al
T_
118 A MÉLIORATIONS LOCALES
17
20
e
br
em
pt
se
25
g
Al
T_
Chapitre 5
17
Méthodes de décomposition
20
5.1 Considérations sur la taille des problèmes
e
La complexité algorithmique très brièvement exposée en section 1.2, a pour
br
but d’évaluer les ressources de calcul nécessaire à l’exécution d’un algorithme
en fonction de la taille des données qu’on lui fournit. On ne peut pas classer
em
les problèmes — petits ou grands — uniquement en fonction de leur taille ab-
solue : trier un tableau de 1000 éléments est assurément beaucoup plus facile
que de trouver la tournée optimale d’un voyageur de commerce sur 100 villes.
Le temps à disposition pour trouver une solution entre évidemment en ligne de
pt
17
tion et la calibration des méthodes heuristiques. En effet, la connaissance
des solutions optimales permet de déterminer la qualité de l’heuristique
et d’optimiser la valeur de ses paramètres en conservant des temps de
20
calcul raisonnables
e
heuristiques. Ces problèmes, fréquemment rencontrés dans les applica-
tions réelles, sont de taille trop importante pour pouvoir être résolus ef-
br
ficacement par des méthodes exactes ou pour qu’un humain réussisse
à deviner des solutions de qualité suffisante. La taille maximale des pro-
em
blèmes traitables par une métaheuristique n’utilisant pas de méthodes de
décomposition est lié à sa complexité algorithmique que ce soit au niveau
du temps de calcul ou de la mémoire. Avec plus de 104 éléments, il de-
vient en effet difficilement envisageable d’utiliser une méthode construc-
pt
tive ou une taille de voisinage en O(n2 ), et ce, d’autant plus si l’on doit
mémoriser, pour des raisons d’efficacité, une matrice n × n. La com-
plexité algorithmique d’un programme basé sur une métaheuristique est
se
une taille de 105 n’est pas exceptionnelle. Dans ce cas, on doit utiliser des
méthodes de décomposition. Ce chapitre présente quelques techniques
générales permettant d’aborder les problèmes de grande taille. Il faut ce-
T_
Les très grands problèmes Taille : n > 107 éléments. Lorsque la taille du
problème dépasse 107 à 109 éléments, il n’est plus possible de stocker
les données entièrement en mémoire vive. Il est nécessaire dans ce cas
de travailler sur des portions du problème, généralement à l’aide d’algo-
rithmes parallèles pour conserver des temps de traitement convenables.
Le traitement de ce type de problèmes soulève des questions d’ordre
essentiellement techniques et sort du cadre de cet ouvrage.
5.2 A LGORITHMES RÉCURSIFS 121
17
— Les morceaux doivent être indépendants les uns des autres
— L’assemblage des morceaux doit être d’une complexité inférieure à celle
20
de la résolution directe du problème
Toute la difficulté vient de la manière dont on définit les morceaux : ils doivent
représenter une portion logique du problème pour que leur assemblage, une
fois résolus, soit simple. Les algorithmes récursifs pour le tri d’un tableau consti-
e
tuent un exemple typique de technique de décomposition : on découpe le ta-
br
bleau en deux parties à peu près égales, celles-ci sont triées par deux appels
récursifs, si elles contiennent plus d’un élément ; finalement les deux parties
localement triées sont parcourues pour reconstituer l’ensemble du tableau.
em
5.2.1 Théorème général de récurrence
L’efficacité d’un algorithme récursif peut s’évaluer dans un certain nombre
pt
tion à être soit plus petite, soit plus grande que nlogb (a) .
Souvent, a = b : on a un appel récursif pour toutes les parties. Dans ce
cas, le théorème établit que si la reconstitution peut se faire en un temps sub-
linéaire, alors on peut traiter le problème en temps linéaire ; si la reconstitution
prend un temps linéaire — ce qui est typiquement le cas pour les algorithmes
de tri — le problème peut se résoudre en O(n log n). Le dernier cas indique
simplement que toute la difficulté de l’algorithme est concentré dans les opéra-
tions de reconstitution. Finalement, il faut mentionner que le théorème ne traite
pas tous les cas pour la fonction f (n).
On rencontre également des cas où a 6= b. Par exemple, la recherche par di-
chotomie dans un tableau trié procède par division du tableau en b = 2 parties
122 M ÉTHODES DE DÉCOMPOSITION
dont une seule (a = 1) devra être traitée récursivement. Pour ce problème, il n’y
a pas de « reconstitution », mais l’effort de calcul à réaliser entre deux appels
récursifs est constant (f (n) = Θ(1)). En appliquant le second cas ci-dessus,
on trouve que T (n) = Θ(nlog2 1 · log n) = Θ(log n).
Un autre exemple est celui de la recherche d’un point du plan parmi un
ensemble de n points stockés sous la forme d’un arbre 2D équilibré (voir la
structure de données discutée en section 4.6.3). Avec une telle structure de
données, on peut diviser par deux le nombre de points restant à examiner
en traitant au maximum a = 2 parties parmi b = 4. En effet, contrairement
17
à un arbre binaire en une dimension, on ne peut pas assurer de diviser par
deux ce nombre à chaque étage de l’arbre, mais seulement tous les deux
étages. Le problème étant celui d’une recherche, il n’y a pas de reconstruction
20
et f (n) = O(1). En remarquant que log4 (2) = 1/2 et en choisissant = 1/2,
on remarque que l’on est dans le premier cas. On en déduit que la complexité
d’une recherche dans un arbre 2D est en Θ(n1/2 ). En pratique cependant, si
les point sont bien répartis, on observe un comportement nettement plus favo-
e
rable, se rapprochant de log n.
Dans le cas d’algorithmes heuristiques procédant par récursion, on arrête
br
communément cette dernière prématurément, avant que les morceaux soient
d’une taille si faible que leur résolution devienne triviale. En effet, la procédure
de découpe d’un problème et celle de reconstitution d’une solution à partir des
em
morceaux se font de manière heuristique, car on ne peut en général garantir
l’optimalité de la reconstitution à partir de morceaux, même si ces derniers ont
été résolus exactement. Cela signifie que les « zones frontières » entre deux
pt
effet, s’ils sont trop gros, leur résolution exacte prend trop de temps ou l’appli-
cation d’une heuristique risque de produire des morceaux de basse qualité.
25
nible. La résolution d’un problème de grande taille implique une limitation sur
Al
17
(a) Ensemble des éléments (b) Échantillon
20
e
br
em
(c) Centres (d) Groupes
pt
L’ensemble des éléments d’un problème. 5.1(b) Échantillonnage des éléments. 5.1(c)
Centres des groupes trouvés avec les éléments de l’échantillon. 5.1(d) Affectation de
l’ensemble des éléments aux centres les plus proches.
25
1. On ne peut utiliser la notation O(·) ici car on fait l’hypothèse que les p groupes contiennent
approximativement le même nombre d’éléments. Dans le pire des cas, quelques groupes pour-
raient en contenir un nombre en O(n), et O(n) groupes pourraient en contenir un nombre en
O(1), ce qui impliquerait une complexité théorique en O(n2 ).
124 M ÉTHODES DE DÉCOMPOSITION
17
20
e
br
em
pt
se
25
F IGURE 5.2 – Décomposition à deux niveaux d’un problème. Les éléments sont re-
√
g
17
√
√En procédant par échantillonnage, et en choisissant a = Θ( n) et b =
Θ( n), on obtient en O(n3/2 ) une décomposition en un nombre de petits mor-
20
ceaux proportionnel à n dont la taille est approximativement identique ce qui,
pour certains problème peut avoir un sens : Par exemple, pour l’élaboration de
tournées de véhicules, le nombre maximum de clients que l’on peut placer sur
une tournée dépend de l’application (réparateur, distribution de colis, ramas-
e
sage d’ordures) et non du nombre total de clients que l’on a dans le problème.
Un exemple de gain possible de complexité par
br
√ rapport à une méthode directe
qui serait en O(n2 ) est de choisir a = b = Θ( n).
em
Cette technique de décomposition est illustrée en figure 5.2, où l’on a super-
posé les relations de proximités entre les gros morceaux et les groupes d’envi-
ron 15 éléments obtenus par décomposition des gros morceaux. Les éléments
pt
d’une procédure en O(n2 ) pour construire une tournée, une construction récur- √
sive peut se faire en O(n3/2 ) : Tout d’abord, on cherche une tournée sur les n
Al
centres correspondant
√ aux centres des gros morceaux. Ceci nécessite un tra-
vail en O(n · n). Ensuite, pour chaque gros morceau, on cherche une tournée
T_
sur les entités du morceau, auxquelles on ajoutera deux unités fictives, cor-
respondant aux connexions qu’il faudra réaliser avec les deux gros morceaux
adjacents.
17
phase de destruction consiste à sélectionner un sous-ensemble de variables,
en incorporant une certaine stochasticité dans ce processus. Dans sa forme la
plus simple, celle-ci consiste à sélectionner un nombre constant de variables,
20
de façon totalement aléatoire. Un forme plus élaborée est de sélectionner aléa-
toirement une variable-germe et un certain nombre d’autres, qui sont le plus
en relation avec la variable-germe. La phase de réparation consiste à tenter
d’améliorer la solution en résolvant un sous-problème sur les variables qui ont
e
été sélectionnées, tout en fixant la valeur des autres variables à celle prise
dans la solution de départ.
br
Le nom de cette technique vient du fait qu’un très grand nombre de possibi-
lités existe pour reconstruire une solution, nombre croissant exponentiellement
em
avec la taille du problème ou du moins plus élevé que celui que l’on pourrait rai-
sonnablement énumérer extensivement. Ainsi la phase de reconstruction con-
siste à choisir une solution parmi un grand nombre de possibilités. Comme la
majeur partie des variables conservent leur valeur d’une solution à la suivante,
pt
il s’agit bien, conceptuellement, d’une recherche locale, mais avec une taille de
voisinage élevée. La trame de LNS est donnée par l’algorithme 5.1.
se
3 s0 ← r(d(s));
4 si a(s, s0 ) alors
5 s ← s0
T_
On remarque que cette trame laisse une très grande liberté au program-
meur pour choisir diverses options :
Méthode de destruction d(·) Cette méthode est censée détruire une partie
de la solution courante. Les auteurs préconisent qu’elle ne soit pas déter-
ministe, de sorte que deux appels successifs détruisent des portions dif-
férentes. Une autre vision de cette méthode est de fixer un certain nom-
5.4 M ÉTHODES D ’ AMÉLIORATION POUR PROBLÈMES DE GRANDE TAILLE 127
17
bitif.
Méthode de reconstruction r(·) Cette méthode est censée reconstruire la
partie d’une solution qui serait détruite. Une autre vision de cette mé-
20
thode est de réoptimiser la portion du problème correspondant aux va-
riables qui ont été libérées par la méthode de destruction. Une option
possible pour la méthode de reconstruction est d’utiliser une méthode
exacte, notamment de la programmation par contrainte. Une autre option
e
est d’utiliser une méthode heuristique, qu’elle soit simple, comme un al-
gorithme glouton, ou plus évoluée, comme une recherche avec tabous,
une recherche à voisinage variable ou autre.
br
Critère d’acceptation a(·, ·) Le plus simple des critères d’acceptation est d’uti-
em
liser la valeur de la fonction-utilité des deux solutions passées en para-
mètre :
0
a(s, s ) =
F aux Sinon
se
D’autres critères ont été proposés, notamment ceux utilisés par des mé-
thodes de type recuit simulé (c.f. section 6.1).
Critère d’arrêt La trame ne fournit aucune suggestion concernant le critère
25
à choisir un client-germe au hasard. Les autres clients sont triés à l’aide d’une
fonction mesurant la relation existante avec le client-germe. Cette fonction est
inversement proportionnelle à la distance entre clients, et au fait que les clients
font partie de la même tournée. L’idée est de sélectionner un sous-ensemble
de clients proches, faisant partie de tournées différentes. Ces clients sont choi-
sis aléatoirement, mais avec un biais pour favoriser le choix de ceux les plus
fortement en relation avec le client-germe.
La méthode de reconstruction est basée sur un modèle de programme li-
néaire en nombres entiers qui est résolu avec une technique de séparation
et évaluation avec propagation de contraintes. On contraint cette méthode à
ne modifier que les variables associées aux clients choisis par la méthode de
128 M ÉTHODES DE DÉCOMPOSITION
17
20
(a) Solution initiale
e
br
em
pt
se
destruction. De plus, pour éviter que le temps de calcul n’explose, ce qui sur-
vient fréquemment avec les méthodes exactes, l’arbre d’énumération n’est pas
examiné complètement, mais élagué heuristiquement. Un cycle de destruction-
reconstruction est illustré en figure 5.3.
Il existe des algorithmes basés sur la trame de LNS qui ont été proposés
bien avant cette dernière. On peut citer notamment l’algorithme de déplace-
ment du goulot d’étranglement [?] pour un problème d’ordonnancement. Dans
cet article, la méthode de destruction sélectionne la machine qui constitue le
goulot d’étranglement. La méthode de reconstruction réordonne les opérations
17
de cette machine, en considérant que les séquences sur les autres machines
ne sont pas modifiées. Cela signifie que pour chaque opération sur la machine-
goulot, on a un délai avant lequel l’opération ne peut avoir lieu (parce que l’élé-
20
ment sur lequel elle a lieu doit subir d’autres opérations sur d’autres machines
préalablement) et un délai après lequel l’opération ne doit pas être faite si l’on
s’interdit de dégrader la solution (parce que l’élément sur lequel elle a lieu doit
subir d’autres opérations sur d’autres machines ultérieurement). Comme tous
e
les choix sont déterministes, que les réoptimisations se font de manière exacte,
et que la solution courante n’est modifiée que si la réoptimisation l’améliore
br
strictement, la méthode a un critère d’arrêt naturel.
La méthode POPMUSIC présentée à la section suivante a été développée
em
indépendamment de LNS. Elle peut être vue comme une méthode de type LNS
un peu plus rigide, dans le sens où elle dirige mieux le programmeur dans le
choix des options, notamment en ce qui concerne le critère d’arrêt.
pt
5.4.2 POPMUSIC
se
lement cette méthode a reçu l’acronyme moins attractif de LOPT (pour Local
Optimizations) [?, ?].
Pour les problèmes de grande taille, on peut considérer qu’une solution
est composée d’un certain nombre de parties, parties qui sont elles-mêmes
g
Parties indépendantes
17
20
Sous-problème
e
Partie-Germe
br
em
F IGURE 5.4 – Pour appliquer la trame POPMUSIC à un problème de classification
pt
avec centres, on peut définir une partie comme l’ensemble des éléments rattachés à
un même centre. Les éléments du groupe de la partie-germe, auxquels on adjoint les
parties les plus proches forment un sous-problème que l’on tente d’optimiser indépen-
se
damment. L’optimisation de groupes très distants ne peut amener de bénéfice, car ces
parties sont de facto indépendantes.
25
1 U = {s1 , . . . , sq };
2 tant que U 6= ∅ faire
3 Sélectionner sg ∈ U // sg : partie-germe ;
T_
17
20
e
br
em
pt
se
25
g
Al
17
l’origine de sa définition est retirée de l’ensemble U . Une fois que U est vide,
le processus s’arrête. Si un sous-problème R a été amélioré avec succès, cela
signifie qu’un certain nombre de parties ont été modifiées. De nouvelles amé-
20
liorations deviennent possible dans leur voisinage. Dans ce cas, toutes les par-
ties de U qui n’existent plus dans la solution améliorée en sont retirées avant
d’y incorporer l’ensemble des parties de R. La trame 5.2 formalise la méthode
POPMUSIC. Pour transcrire cette trame en un code pour un problème donné,
e
plusieurs choix doivent être pris :
br
Obtention de la solution initiale La méthode a besoin d’une solution avant
de démarrer. La technique présentée dans la section 5.3 permet d’obtenir
une solution initiale appropriée pour la trame POPMUSIC avec un effort
em
de calcul limité. Cependant, une méthode basée sur POPMUSIC peut
très bien s’avérer efficace pour des problème de taille limitée et on peut
démarrer avec un solution obtenue avec un algorithme plus gourmand en
ressources de calcul.
pt
Définition d’une partie La définition de ce qu’est une partie n’est de loin pas
unique pour un problème donné. Par exemple, pour les problèmes d’éla-
se
17
20
F IGURE 5.6 – Pour l’étiquetage de plans, une partie peut être un objet devant être
étiqueté. Ici, on considère 4 positions possibles pour l’étiquette de chaque objet. Deux
e
objets seront à une distance de 1 si leurs étiquettes peuvent se chevaucher. Le chiffre
br
à l’intérieur de chaque disque représente sa distance à l’objet-germe, noté 0. Un sous-
problème sera constitué des r = 25 objets les plus proches de l’objet-germe (dont
la distance est ici au maximum 4) auxquels ont adjoint les objets dont les étiquettes
em
pourraient entrer en collisions avec ces r objets. Seules les positions des étiquettes de
ces derniers peuvent être modifiées lors de l’optimisation d’un sous-problème.
pt
de préciser ici que l’implantation de celle-ci est facilitée par le fait que
la taille des sous-problèmes peut être déterminée, notamment en ce qui
T_
Avec un critère d’arrêt basé sur le fait que l’ensemble U est vide, l’effort de
calcul pourrait potentiellement devenir prohibitif pour des problèmes de grande
taille, sachant que plusieurs parties sont susceptibles d’y être introduites pour
chaque optimisation de sous-problème. En pratique, on observe qu’il n’en est
rien : le nombre de sous-problèmes à résoudre croît quasi-linéairement avec la
taille du problème. Ceci est illustré en figure 5.7 pour un problème de localisation-
routage [?].
134 M ÉTHODES DE DÉCOMPOSITION
17
Temps de calcul [s]
105
20
104
e
103
br
em
102
10 Θ(n)
Solution initiale
se
Θ(n3/2 )
1
104 105 106
25
Taille du problème
F IGURE 5.7 – Temps de calcul observés pour la création d’une solution initiale d’un
problème de localisation-routage avec la technique présentée en section 5.3 et temps
global d’optimisation des sous-problèmes avec la trame POPMUSIC. On remarque que
g
la croissance du temps de calcul semble inférieure à l’analyse de Θ(n3/2 ) qui a été faite
Al
en section 5.3 et que le temps pour l’optimisation des sous-problèmes est pratiquement
linéaire.
T_
5.4 M ÉTHODES D ’ AMÉLIORATION POUR PROBLÈMES DE GRANDE TAILLE 135
17
chaque fois décalés d’une ville. Pour déterminer s’il faut continuer de tenter des
optimisations, on mémorise la ville initiale de la dernière portion de chemin qui
a été optimisée avec succès. Si l’on fait un tour complet sans amélioration, le
20
processus peut s’arrêter.
e
void tsp_3opt_limited ( i n t n , / ∗ Number o f c i t i e s ∗/
double ∗∗d , /∗ Distance matrix ∗/
int r ,
i n t best_sol [ ] ,
double ∗ b e s t _ c o s t )
br
/ ∗ Maximum n r o f nodes o f subpath o p t i m i z e d
/ ∗ InOut S o l u t i o n p r o v i d e d and r e t u r n e d
/ ∗ InOut Length o f t h e b e s t t o u r
∗/
∗/
∗/
em
{
i n t ∗s = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ; / ∗ Successors o f c i t i e s ∗/
i n t i , j , k , t , u , temp , / ∗ indexes ∗/
i_end ; /∗ Last values of index i ∗/
double d e l t a ; / ∗ Cost o f a move ∗/
pt
s [ b e s t _ s o l [ i −1]] = b e s t _ s o l [ i ] ;
s [ b e s t _ s o l [ n −1]] = b e s t _ s o l [ 0 ] ;
i_end = i ; / ∗ Update l a s t v a l u e s o f i n d e x i ∗ /
}
}
i = s[ i ];
T_
}
w h i l e ( i ! = i_end ) ; / ∗ A complete t o u r scanned w i t h o u t improvement ∗ /
j = 0; / ∗ R e b u i l d s o l u t i o n under t h e form o f a p e r m u t a t i o n o f c i t i e s ∗ /
f o r ( i = 0 ; i < n ; i ++)
{ best_sol [ i ] = s [ j ] ; j = s [ j ] ; }
free ( s ) ;
}
— Plutôt que de mettre au point une recherche locale ad-hoc comme celle
du code 5.1 pour optimiser des sous-chemins, il serait préférable d’utili-
ser une méthode générale de résolution d’un problème de voyageur de
commerce, comme par exemple le code 4.4.
17
Pour démarrer avec une solution avec une bonne structure sans utiliser
un algorithme de complexité trop élevée, on peut s’inspirer de la technique
20
présentée en section 5.3. On commence par échantillonner aléatoirement k
des n villes, avec k < n, un paramètre à choisir en fonction de n de telle sorte
que la complexité algorithmique de la méthode soit aussi faible que possible ;
on y reviendra. Un problème de voyageur de voyageur de commerce est résolu
e
(éventuellement de façon heuristique) sur les villes échantillonnées à l’aide
br
d’une méthode d’optimisation de complexité c(k). On trouve pour chacune des
n − k villes restantes la ville de l’échantillon qui lui est la plus proche. Ceci peut
se faire trivialement en Θ(k(n − k)). Chaque ville restante est ensuite insérée
em
(dans un ordre quelconque) dans la tournée partielle, juste après la ville de
l’échantillon qui lui est le plus proche. Ceci peut se faire également trivialement
en Θ(n − k).
On obtient ainsi une tournée mauvaise, mais qui a de bonne propriétés
pt
Une très bonne tournée sur un échantillon de k villes peut être obtenue avec
le code 4.4 implantant des chaînes d’éjections. On observe empiriquement
T_
17
20
e
br
em
pt
se
(a) Échantillon et tournée complète en cours (b) Tournée optimisée avec POPMUSIC et,
d’optimisation en superposition, la tournée optimale
25
sur cette tournée partielle juste après la ville de l’échantillon qui lui est le plus proche.
Ensuite, on améliore avec le voisinage 2-opt (Code 4.3) des sous-chemins compor-
Al
tant deux groupes de villes insérées après deux villes successives de l’échantillon.
La figure 5.8(a) montre la tournée sur l’échantillon ainsi que la tournée complète en
cours d’optimisation. Finalement, des sous-chemins de 50 villes sont optimisés à l’aide
T_
de POPMUSIC en utilisant une recherche locale basée sur les chaînes d’éjection
(Code 4.4). La tournée optimale est superposée en trait fin à celle obtenue par POP-
MUSIC.
138 M ÉTHODES DE DÉCOMPOSITION
tances D suivante :
0 d12 + L d13 + L ... 0
d21 + L 0 d23 ... d2r + L
D=
d31 + L d32 0 ... d3r + L
.. .. .. .. ..
. . . . .
0 dr2 + L dr3 + L ... 0
17
chemin. Pour éviter cela, on peut associer à chaque ville un drapeau indiquant
si elle a déjà été la ville de départ d’un sous-chemin optimisé. Lors d’une opti-
misation réussie, on peut abaisser tous les drapeaux des villes qui se trouvent
20
jusqu’à 50 position avant ou après une séquence de ville qui a été modifiée
durant l’optimisation (tout en laissant levé le drapeau de la ville de départ du
sous-chemin).
e
5.4.3 Commentaires
br
La principale différence entre LNS et POPMUSIC est que cette dernière
fixe le critère d’arrêt et celui d’acceptation d’une solution du large voisinage.
em
En effet, dans la trame de POPMUSIC on accepte de modifier la solution que
si on a une amélioration stricte. Pour plusieurs problèmes, cette trame semble
suffisante pour obtenir des solution de bonne qualité, cette dernière étant for-
tement conditionnée à la capacité de la méthode d’optimisation utilisée. La phi-
pt
losophie de POPMUSIC est de conserver une trame aussi simple que possible
et d’améliorer au besoin la méthode d’optimisation pour qu’elle puisse mieux
traiter des sous-problèmes de plus grande taille plutôt que de complexifier la
se
contraintes que l’on ajoute au problème sur la base d’une solution existante.
Ces contraintes permettent d’utiliser une méthode d’optimisation qui serait in-
applicable au problème entier. Du reste, [?] dans leur Méthode du Corridor
prennent la problématique par l’autre bout : étant donné une méthode d’opti-
g
toutes interdépendantes. Le choix d’une option ayant une incidence sur les
autres explique en partie pourquoi des méthodes en réalité très similaires sont
présentées par des trames et des noms différents.
5.5 E XERCICES 139
5.5 Exercices
Exercice 5.1. Complexité de la recherche dichotomique En appliquant le
théorème général de récurrence de la section 5.2.1, déterminer la complexité
algorithmique de la recherche d’un élément dans un tableau trié au moyen
d’une recherche dichotomique
Exercice 5.2. POPMUSIC pour le problème de la chaîne de traitement
Dans le cadre de la méthode de décomposition POPMUSIC, comment définir
une partie et un sous-problème pour le problème de la chaîne de traitement ?
17
Comment tenir compte de l’interaction entre le sous-problème et les parties qui
ne doivent pas être optimisées ?
20
Exercice 5.3. Complexité algorithmique de POPMUSIC Si la taille des sous-
problèmes est indépendante de la taille du problème, on peut considérer qu’un
sous-problème peut être résolu en temps constant. Des observations empi-
riques, comme celles présentées dans la figure 5.7 montrent que le nombre de
e
fois qu’une partie est insérée dans U est également indépendant de la taille du
problème, en moyenne. Si ces hypothèses sont satisfaites, quelles sont les par-
br
ties les plus complexes de la trame de POPMUSIC, en termes de complexité
algorithmique ?
em
Exercice 5.4. Complexité minimale de POPMUSIC pour le voyageur de
commerce Justifier le choix de la taille de l’échantillon en Θ(n0.56 ) suggéré en
page 136 pour une implantation de complexité empirique minimale de POP-
pt
17
20
e
br
em
pt
se
25
g
Al
T_
17
20
Troisième partie
e
br
Techniques populaires
em
pt
se
25
g
Al
T_
T_
Al
g
25
se
pt
em
br
e
20
17
Chapitre 6
17
Méthodes randomisées
20
En appliquant les principes présentés dans les chapitres précédents, nous
arrivons à construire une solution et à l’améliorer jusqu’à trouver un optimum
e
local. De plus, si le problème est compliqué ou de grande taille, il faut le dé-
br
composer en sous-problèmes plus simples à résoudre. Que faire si en utilisant
ces techniques on obtient des solutions dont la qualité n’est pas suffisante et
que l’on a la possibilité d’investir plus d’effort de calcul pour tenter d’autres
em
améliorations ?
Une première option venant à l’esprit est d’essayer de mettre en place un
processus d’apprentissage, ce que nous examinerons dans des chapitres ul-
térieurs.
pt
Une deuxième option — plus simple à mettre en œuvre a priori — est d’in-
corporer une composante aléatoire dans une méthode « d’amélioration » où
se
la méthode GRASP.
Finalement, on peut mettre en œuvre une stratégie dans laquelle on fait
alterner des phases d’intensification et de diversification de la recherche. La
T_
recherche à voisinage variable fait usage de cette stratégie, en itérant des mé-
thodes d’améliorations utilisant un voisinage de base après avoir dégradé la
solution par des mouvements choisis aléatoirement dans d’autres voisinages.
interne.
État visqueux
17
Refroidissement lent Refroidissement rapide
20
Recuit Trempe
e
br
em
pt
plus rapide, on parle de trempe et les cristaux se forment de façon désordonnée. Dans
certaines conditions, un refroidissement très rapide ne laisse pas le temps aux cristaux
de se former et la matière reste dans un état amorphe.
g
grains qui sont de petits cristaux dont l’orientation est différente pour chaque
grain. Il s’agit du processus de trempe, qui est utilisé notamment pour durcir
certains aciers.
Au contraire, si le refroidissement est lent, les molécules arriveront à former
des cristaux beaucoup plus grands, correspondant à leur état d’énergie mini-
male. En répétant le procédé, on peut encore augmenter la taille des cristaux,
voire à obtenir un mono-cristal, c’est le procédé de recuit.
[?] et [?] ont eu indépendamment l’idée de simuler ce processus en optimi-
sation combinatoire, en faisant l’analogie entre une fonction-objectif à minimi-
ser et l’énergie des molécules.
Si à haute température une molécule a une énergie suffisante pour aller
6.1 R ECUIT SIMULÉ 145
17
fie que le mouvement m améliore la solution s, alors on l’accepte et la nouvelle
solution devient s⊕m. Sinon, le mouvement m est tout de même accepté, mais
avec une probabilité proportionnelle à e−∆/T , où T est un paramètre simulant
20
la température.
À chaque étape, la température T est diminuée. Plusieurs formules ont été
proposées pour ajuster la température. Parmi les plus fréquemment rencon-
T
trées, citons T ← α · T et T ← 1+αT , où 0 < α < 1 est un paramètre réglant la
e
vitesse de décroissance de la température.
La méthode comporte encore deux paramètres, au minimum : Tinit et Tf in ,
br
les températures initiales et finales. L’algorithme 6.1 donne la trame d’un recuit
simulé de base.
em
Algorithme 6.1 : Trame d’un recuit simulé élémentaire. D’innombrables va-
riantes d’algorithmes basés sur cette trame ont été proposés.
Entrées : Solution initiale s ; fonction utilité f à minimiser ; structure de
pt
1 T ← Tinit ;
2 tant que T > tf in faire
3 Générer un mouvement m aléatoirement dans M ;
4 ∆ = f (s ⊕ m) − f (s);
25
nales ont une influence très différente selon l’unité de mesure de la fonction
utilité f . En effet, si f mesure la longueur d’une tournée de voyageur de com-
merce, il convient d’adapter ces paramètres selon que cette mesure est expri-
mée en mètres ou en kilomètres.
Pour rendre l’algorithme plus robuste, on ne demande pas à l’utilisateur de
fournir directement des températures, mais par exemple des taux d’acceptation
de mouvements dégradants τinit et τf in , ce qui est beaucoup plus intuitif. Les
températures sont alors calculées automatiquement en fonction de ces taux,
en effectuant quelques centaines ou milliers de pas d’une marche aléatoire,
par enregistrement de statistiques sur les valeurs moyennes des ∆.
La figure 6.2 illustre l’évolution de la longueur de la tournée d’un voyageur
146 M ÉTHODES RANDOMISÉES
17
20
Longueur de la tournée
e
br
em
Amélioration de la solution
pt
Température (x100)
se
Itérations
de commerce sur 225 villes en fonction du nombre d’itérations pour une exé-
cution de l’algorithme 6.1 avec comme température initiale 5 · dmax /n, comme
Al
température finale 20 · dmax /n2 et α = 0.99. Il est fourni à l’algorithme une so-
lution de relativement bonne qualité. On remarque une importante dégradation
T_
Listing 6.1 – tsp_SA.c Code C implantant une variante de recuit simulé pour le pro-
blème du voyageur de commerce symétrique. La fonction build_2opt_data_structure
qui initialise la structure de donnée décrite en section 4.2.3 pour effectuer un mouve-
ment de type 2opt en temps constant est données par le code 4.2.
v o i d tsp_SA ( i n t n , / ∗ Number o f c i t i e s ∗/
double ∗∗d , / ∗ D i s t a n c e m a t r i x , must be s y m m e t r i c a l ∗/
i n t best_sol [ ] , / ∗ InOut S o l u t i o n p r o v i d e d and r e t u r n e d ∗/
double ∗ b e s t _ c o s t , / ∗ InOut Length o f t h e b e s t t o u r ∗/
double i n i t i a l _ t e m p e r a t u r e , / ∗ SA Parameters ∗/
double f i n a l _ t e m p e r a t u r e ,
double alpha )
17
{ int∗ t ; / ∗ Successor and predecessor o f each c i t y [ see tsp −2o p t ] ∗/
i n t i , j , next_i , next_j , r , s ; /∗ indices ∗/
long i t e r a t i o n ; /∗ I t e r a t i o n counter ∗/
double l e n g t h = ∗ b e s t _ c o s t ; / ∗ Length o f c u r r e n t t o u r ∗/
double T = i n i t i a l _ t e m p e r a t u r e ; /∗ Current temperature ∗/
20
double d e l t a ; / ∗ Cost o f a t r i a l move ∗/
t = ( i n t ∗) m a l l o c ( ( s i z e _ t ) ( 2 ∗ n ) ∗ s i z e o f ( i n t ) ) ;
build_2opt_data_structure ( n , best_sol , t ) ;
i t e r a t i o n = 0; / ∗ Number o f i t e r a t i o n s performed ∗/
e
while (T > final_temperature ) / ∗ Repeat , w h i l e t e m p e r a t u r e h i g h enough ∗/
{ i = 0; / ∗ S t a r t moves e v a l u a t i o n from c i t y 0 ∗/
{ j = t [ t [ i ]];
w h i l e ( j >>1 && ( t [ j ] > >1 | | i > >1) )
{ delta = d [ i > > 1 ] [ j > >1] br
w h i l e ( t [ t [ i ]] > >1 && t [ i ] > >1) / ∗ Index i a r r i v e d a t t h e l a s t t o c o n s i d e r ?
i = t [ i ]; /∗ Explanation l e f t in exercise ∗/
j = t [ i ];
se
p r i n t f ( "SA %l d %f \ n " , i t e r a t i o n , l e n g t h ) ;
}
} / ∗ Move accepted ∗ /
i t e r a t i o n ++;
i f ( i t e r a t i o n % ( n∗n ) == 0 ) / ∗ Temperature decrease every n∗n i t e r . ∗ /
g
T ∗= alpha ;
j = t [ j ]; / ∗ Next j ∗ /
Al
}
i = t [ i ] ; / ∗ Next i ∗ /
}
} /∗ while (T > final_temperature ) ∗/
T_
free ( t ) ;
}
Le code 6.1 implante une sorte de recuit simulé pour le problème du voya-
geur de commerce. Il est basé sur une structure de voisinage 2-opt. Afin de
pouvoir évaluer un mouvement et le réaliser en temps constant, la structure de
données introduite en section 4.2.3 a été utilisée. Cependant, cette structure
de données ne permet pas d’évaluer en un temps constant un mouvement
(i, j) tiré aléatoirement, car, étant donné une ville i et celle qui lui succède
si lorsqu’on fait le parcours dans un sens donné, on ne peut pas identifier la
ville sj qui succède à j autrement qu’en suivant les successeur de i jusqu’à
atteindre j puis sj .
148 M ÉTHODES RANDOMISÉES
L’artifice utilisé dans ce code est de ne pas tirer les mouvement aléatoire-
ment, mais d’examiner systématiquement tout le voisinage. L’autre adaptation
de la trame 6.1 est de ne pas diminuer la température à chaque itération, mais
seulement toutes les n2 itérations. Ainsi, une valeur de α entre 0.8 et 0.99
donne des résultats acceptables, indépendamment de la taille du problème, et
donc du nombre d’itérations nécessaires à l’obtention de bonnes solutions. Fi-
nalement, mentionnons que l’utilisateur doit fournir les températures initiales et
finales de façon absolue. En pratique, Tinit = 5 · dmax /n et Tf in = 20 · dmax /n2
sont des valeurs qui peuvent convenir pour démarrer un ajustement de ces
17
paramètres.
20
6.2 Acceptation à seuil
La méthode de l’acceptation à seuil est une recherche locale pure, dans la
mesure où l’on ne fait que se déplacer d’une solution à l’une de ses voisines.
e
Tout comme les méthodes du démon, du grand déluge ou du recuit simulé,
les mouvements sont choisis aléatoirement, mais ne sont pas systématique-
br
ment appliqués à la solution courante. Dans le cas où le mouvement permet
une amélioration, la solution voisine est acceptée. Si par contre le mouvement
em
détériore la solution, le mouvement n’est accepté que si la détérioration est
inférieure à un certain seuil. Ce dernier est progressivement diminué jusqu’à
prendre une valeur nulle, de sorte que la méthode s’arrête dans un optimum
local.
pt
Algorithme 6.2 : Trame d’une acceptation à seuil. Les valeurs des seuils T1 à
se
TImax ne sont pas forcément fournis explicitement mais calculés, par exemple en
fournissant uniquement le seuil initial et en le multipliant par un autre paramètre α
à chaque ronde de R itérations.
Entrées : Solution initiale s ; fonction utilité f à minimiser ; structure de
25
voisinage M , paramètres T , R, τ1 , . . . τT
Résultat : Solution s∗
1 s∗ ← s;
2 pour t allant de 1 à T faire
g
[?] ont proposé une méthode automatique pour fixer les seuils : on com-
mence par effectuer des mouvements aléatoires en enregistrant les variations
de la valeur de la fonction-objectif (en valeur absolue). Cela permet de déter-
miner la fonction de distribution empirique F de l’amplitude des mouvements.
On fixe les seuils en utilisant l’inverse de cette fonction (voir illustration en fi-
gure 6.3) : Ti = F −1 (0.8·(Imax −i)/Imax ). Autrement dit, le premier seuil permet
6.3 G RAND DÉLUGE 149
0.8
17
20
x
e
0
τT = 0 τ2 τ1 ∆max
br
F IGURE 6.3 – Technique de détermination des seuils τ1 , . . . τT pour l’acceptation à
em
seuils. La fonction de répartition empirique F est obtenue en effectuant un certain nom-
bre de mouvements et en enregistrant leur amplitude en valeur absolue.
pt
terrestres finirent tous noyé, sauf ceux ayant été embarqués sur l’arche de
Noé. On imagine que les animaux restés à terre, pris de panique, couraient
Al
qui sont acceptés tant qu’ils ne mènent pas à une solution dont la qualité est
inférieur à un seuil L — le niveau de l’eau — que l’on augmente de P — le
débit de la pluie — à chaque itération. Le processus s’arrête lorsque la valeur
de la solution courante est inférieure à L. L’algorithme 6.3 donne la trame de
cette méthode.
150 M ÉTHODES RANDOMISÉES
17
voisinage M , paramètres L et P
Résultat : s∗
1 s∗ ← s;
20
2 tant que f (s) < L faire
3 Tirer un mouvement m aléatoirement dans M ;
4 si f (s ⊕ m) > L alors
5 s ← s ⊕ m;
e
6 si f (s) < f (s∗ ) alors
s∗ ← s
7
8 L ← L + P;
br
em
pt
se
25
g
Al
T_
17
Traduit en termes de recherche locale, une mise consiste à autoriser un
mouvement dégradant la solution d’une quantité au plus égale au montant à
disposition. Si la modification améliore la solution, le budget est augmenté d’au-
20
tant, jusqu’à concurrence du seuil maximum. L’algorithme 6.4 donne la trame
de cette méthode.
e
à implanter, mais, comme pour l’acceptation à seuil, ses paramètres doivent être
ajustés en fonction de la valeur numérique des données.
br
Entrées : Solution s, fonction-utilité à minimser f , structure de voisinage
M , paramètres Imax , D, Dmax
em
Résultat : s∗
1 s∗ ← s;
2 pour Imax itérations faire
3 Tirer un mouvement m aléatoirement dans M ;
pt
4 ∆ = f (s ⊕ m) − f (s);
5 si ∆ 6 D alors
se
10 s ← s ⊕ m;
g
Al
17
4 si f (s ⊕ m) + bruit(i) < f (s) alors
5 s ← s ⊕ m;
6 si f (s) < f (s∗ ) alors
20
7 s∗ ← s
e
de bruit aléatoire. Généralement, l’espérance de la loi est nulle et sa variance
br
diminue avec le numéro d’itération, ce qui fait qu’à la fin de l’algorithme, la
méthode se rapproche de plus en plus d’une méthode d’amélioration.
em
Répartition des niveaux de bruit
pt
se
0
25
Itérations
g
Al
17
i n t i , j , next_i , next_j , r , s ; /∗ indices ∗/
long i t e r a t i o n ; /∗ I t e r a t i o n counter ∗/
double l e n g t h = ∗ b e s t _ c o s t ; / ∗ Length o f c u r r e n t t o u r ∗/
20
double d e l t a ; / ∗ Cost o f a t r i a l move ∗ /
t = ( i n t ∗) m a l l o c ( ( s i z e _ t ) ( 2 ∗ n ) ∗ s i z e o f ( i n t ) ) ;
build_2opt_data_structure ( n , best_sol , t ) ;
i t e r a t i o n = 0; / ∗ Number o f i t e r a t i o n s performed ∗/
while ( current_noise > final_noise ) / ∗ Repeat , w h i l e n o i s e h i g h enough ∗/
e
{ i = 0; / ∗ S t a r t moves e v a l u a t i o n from c i t y 0 ∗/
w h i l e ( t [ t [ i ]] > >1 && t [ i ] > >1) / ∗ Index i a r r i v e d a t t h e l a s t t o c o n s i d e r ? ∗/
br
{ j = t [ t [ i ]];
w h i l e ( j >>1 && ( t [ j ] > >1 | | i > >1) )
{ delta = d [ i > > 1 ] [ j > >1] + d [ t [ i ] > > 1 ] [ t [ j ] > >1]
− d [ i > > 1 ] [ t [ i ] > >1] − d [ j > > 1 ] [ t [ j ] > > 1 ] ;
em
n o i s e = ( rando ( ) −0.5) ∗ c u r r e n t _ n o i s e ;
i f ( d e l t a + noise < 0) / ∗ Move accepted ∗ /
{ next_i = t [ i ] ; next_j = t [ j ] ;
t [ i ] = j ^1; t [ next_i ^1] = next_j ;
t [ j ] = i ^1; t [ next_j ^1] = next_i ;
length = length + delta ;
pt
i = t [ i ]; / ∗ To a v o i d r e v e r s i n g i m m e d i a t e l y a degrading move ∗ /
j = t [ i ];
se
{ b e s t _ s o l [ s ] = r > >1; r = t [ r ] ; }
p r i n t f ( " N o i s i n g method improvement %l d : %f \ n " , i t e r a t i o n , l e n g t h ) ;
}
} / ∗ Move accepted ∗ /
i t e r a t i o n ++;
i f ( i t e r a t i o n % ( n∗n ) == 0 ) / ∗ Noise decrease every n∗n i t e r . ∗ /
g
c u r r e n t _ n o i s e ∗= alpha ;
j = t [ j ]; / ∗ Next j ∗ /
}
Al
i = t [ i ] ; / ∗ Next i ∗ /
}
} /∗ while ( current_noise > f i n a l _ n o i s e ) ∗/
T_
free ( t ) ;
}
6.6 GRASP
La méthode GRASP, dont le nom ne vient pas de l’anglais « empoigner »,
mais qui est l’acronyme de Greedy Randomized Adaptive Search Procedure,
soit une recherche gloutonne randomisée suivie d’une méthode d’amélioration.
Cette méthode comporte deux paramètres, Imax , le nombre de répétitions de
la boucle la plus externe de l’algorithme et α, permettant de régler le degré de
154 M ÉTHODES RANDOMISÉES
17
f à minimiser, paramètres Imax et 0 6 α 6 1, méthode
d’amélioration a
Résultat : Solution complète s∗
20
1 f ∗ ← ∞;
2 pour Imax itérations faire
3 Initialiser s à une solution partielle triviale;
4 R ← E; // Éléments pouvant encore être ajoutés à s
e
5 tant que R 6= ∅ faire
Trouver cmin = mine∈R c(s, e) et cmax = maxe∈R c(s, e);
br
6
7 Choisir aléatoirement, uniformément un e0 ∈ R tel que
cmin 6 c(s, e0 ) 6 cmin + α(cmax − cmin );
em
8 s ← s ∪ e0 ; // Ajouter e0 à la solution partielle s
9 Retirer de R les éléments qui ne peuvent plus être utilisés pour
compléter s;
10 s0 ← a(s); // Trouver l’optimum local associé à s
pt
13 s∗ ← s0
à moins que plusieurs éléments aient les mêmes coûts incrémentaux, la ré-
pétition de constructions n’a pas de sens. Pour un paramètre α = 1, on a
une construction purement aléatoire. La méthode se transforme en répétitions
d’une méthode d’amélioration partant de solutions aléatoires ; cette technique
g
est souvent pratiquée car elle permet d’obtenir de bien meilleures solutions
Al
17
double l e n g t h ; / ∗ Cost o f c u r r e n t s o l u t i o n ∗/
int nr_iteration ;
double cmin , cmax ;
s o l u t i o n = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
20
f o r ( n r _ i t e r a t i o n = 0 ; n r _ i t e r a t i o n < i t e r a t i o n s ; n r _ i t e r a t i o n ++)
{ generate_random_permutation ( n , s o l u t i o n ) ;
f o r ( i = 0 ; i < n−1; i ++) /∗ s o l u t i o n [ i ] i s the l a s t c i t y i n s e r t e d ∗/
{ / ∗ Determine min and max i n c r e m e n t a l c o s t s ∗ /
cmin = i n f i n i t e ; cmax = − i n f i n i t e ;
f o r ( j = i +1; j < n ; j ++)
e
{ i f ( cmin > d [ s o l u t i o n [ i ] ] [ s o l u t i o n [ j ] ] )
cmin = d [ s o l u t i o n [ i ] ] [ s o l u t i o n [ j ] ] ;
}
br
i f ( cmax < d [ s o l u t i o n [ i ] ] [ s o l u t i o n [ j ] ] )
cmax = d [ s o l u t i o n [ i ] ] [ s o l u t i o n [ j ] ] ;
em
/ ∗ Find t h e n e x t c i t y t o i n s e r t , biased on l o w e r c o s t ∗ /
f o r ( j = i +1; d [ s o l u t i o n [ i ] ] [ s o l u t i o n [ j ] ] > cmin + alpha ∗(cmax−cmin ) ; j ++)
{ } / ∗ C i t i e s t o i n s e r t a l r e a d y i n random o r d e r ; f i r s t c h o i c e i s random ∗ /
swap(& s o l u t i o n [ i + 1 ] , & s o l u t i o n [ j ] ) ;
} / ∗ s o l u t i o n i s complete now ∗ /
pt
}
free ( solution ) ;
return best_cost ;
} / ∗ tsp_GRASP ∗ /
g
Al
che appelée parfois oscillations stratégiques qui consiste à alterner des phase
d’intensification et de diversification de la recherche. L’idée à la base de cette
méthode est de s’appuyer sur plusieurs voisinages M1 . . . Mp . Le voisinage
M1 est utilisé de façon standard pour trouver des optima locaux ; il s’agit d’une
des manière les plus simple d’intensifier la recherche. Les autres voisinages
permettent de s’échapper de ces optima en effectuant des mouvements aléa-
toires qui, généralement, détruisent de plus en plus leur structure. Le tirage
de mouvement aléatoires est également une des manière les plus simples de
diversifier la recherche, pour explorer d’autres portions de l’espace des solu-
tions. La trame d’une méthode à voisinage variable de base est donnée par
l’algorithme 6.7.
156 M ÉTHODES RANDOMISÉES
Algorithme 6.7 : Méthode à voisinage variable. Lorsque l’on dispose d’un faible
nombre p de voisinage, on répète plusieurs fois l’algorithme.
Entrées : Solution s, fonction-utilité f à minimiser, structures de
voisinage M1 . . . Mp
Résultat : s∗
1 s∗ ← s;
2 k ← 1;
3 tant que k 6 p faire
4 Tirer un mouvement aléatoire m ∈ Mk ;
17
5 s ← s ⊕ m;
6 Trouver l’optimum local s0 associé à s dans le voisinage M1 ;
7 si f (s0 ) < f (s∗ ) alors
20
8 s∗ ← s0 ;
9 k←1
10 sinon
s ← s∗ ;
e
11
12 k ←k+1
br
em
Listing 6.4 – tsp_VNS.c Implantation d’une méthode à voisinage variable pour le pro-
blème du voyageur de commerce. Le voisinage Mk consiste à échanger deux villes k
fois. La méthode permettant de réparer une solution perturbée est une chaîne d’éjec-
tions. La particularité de cette implantation, outre son extrême simplicité, est de ne
pt
s o l u t i o n = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
k = 1;
while ( k < n)
{ f o r ( i = 0 ; i < n ; i ++)
solution [ i ] = best_sol [ i ] ;
g
f o r ( i = 0 ; i < k ; i ++)
swap(& s o l u t i o n [ u n i f ( 0 , n −1)] , & s o l u t i o n [ u n i f ( 0 , n − 1 ) ] ) ;
Al
17
20
e
br
em
pt
se
25
g
Al
T_
158 M ÉTHODES RANDOMISÉES
17
20
e
br
em
pt
se
25
g
Al
T_
Chapitre 7
17
Apprendre à construire
20
Après les quatre principes de base vus précédemment — modélisation,
e
décomposition, construction et amélioration — nous abordons dans ce cha-
pitre le cinquième principe : celui de l’apprentissage. Les algorithmes vus dans
br
le chapitre précédent se basent uniquement sur le hasard pour tenter d’ob-
tenir des solutions meilleures que celles que pourraient fournir des méthodes
em
constructives gloutonnes ou des recherches locales. Ceci n’est sans doute pas
très satisfaisant du point de vue intellectuel, bien que, comme dit l’adage, « le
hasard fait bien les choses ». Sans y renoncer totalement, nous examinerons
dans ce chapitre comment mettre en œuvre des techniques d’apprentissage
pt
— Répéter des expériences et analyser les succès (ou les échecs, car on
apprend peut-être plus en faisant des erreurs !)
25
17
20
e
br
em
pt
se
F IGURE 7.1 – Comportement d’une colonie de fourmi séparée d’une source de nour-
riture par un parcours qui bifurque. Initialement les fourmis se répartissent également
dans les deux branches 7.1(a). Les fourmis ayant emprunté le chemin le plus court
g
vont arriver avant à la source de nourriture et vont donc déposer plus rapidement une
quantité additionnelle de phéromones sur le chemin du retour. La quantité de phéro-
Al
mones déposée sur le plus court chemin croît plus rapidement. Après un certain temps,
pratiquement toutes les fourmis vont utiliser la branche la plus courte 7.1(b)
T_
7.1 F OURMIS ARTIFICIELLES 161
17
éclaireurs vont se répartir également dans les deux branches 7.1(a). Tout en
explorant, chaque fourmi dépose une trace chimique qu’elle est capable de dé-
tecter avec ses antennes, ce qui l’aidera lors de son retour dans la fourmilière.
20
Une telle substance chimique porteuse d’information est appelée phéromone.
Une fourmi ayant découvert une source de nourriture va déposer une quan-
tité d’autant plus grande de phéromones sur le chemin du retour que la source
est importante. Naturellement, une fourmi ayant découvert un chemin court
e
va pouvoir retourner plus rapidement que celle qui a emprunté la mauvaise
br
branche. Par conséquent, la quantité de phéromones déposées sur le plus
court chemin croît plus rapidement. Le corollaire est qu’une nouvelle fourmi
disposera d’une information sur le chemin à prendre et biaisera son choix
em
en faveur de la branche la plus courte. Après un certain temps, on observe
que pratiquement toutes les fourmis utilisent la branche la plus courte 7.1(b).
Ainsi, la colonie arrive à déterminer collectivement un chemin optimal, alors
que chaque individu ne voit pas plus loin que le bout de ses antennes.
pt
d’une solution, les valeurs des éléments constituant cette dernière seront
augmentées d’une quantité d’autant plus grande que la solution est de
bonne qualité.
— Le phénomène d’oubli sera simulé par l’évaporation des traces au cours
du temps.
Reste ensuite à préciser comment toutes ces composantes peuvent être
mises en place. Le processus de construction va utiliser une technique de
construction randomisée, presque similaire à la méthode GRASP, si ce n’est
que la composante aléatoire doit être biaisée non seulement par la fonction
ce coût incrémental c(s, e), qui représente l’intérêt a priori d’inclure l’élément
162 A PPRENDRE À CONSTRUIRE
e dans la solution partielle, mais également par la valeur τe qui sera l’intérêt
a posteriori de cet élément, intérêt qui ne sera connu qu’après avoir construit
une multitude de solutions. Le mariage de ces deux formes d’intérêts se réa-
lise en choisissant le prochain élément e à inclure dans la solution partielle s
avec une probabilité proportionnelle à τeα · c(s, e)β , ou α > 0 et β < 0 sont deux
paramètres de la méthode permettant de balancer l’importance que l’on donne
respectivement à la mémoire et au coût incrémental.
La mise à jour des phéromones artificielles se fait en deux étapes, néces-
sitant chacune un paramètre. Tout d’abord, on simule l’évaporation des phéro-
17
mones en multipliant toutes les valeurs τe par 1 − ρ, où 0 6 ρ 6 1 représente
le taux d’évaporation. Ensuite, chaque élément e constituant une solution s ve-
nant d’être construite verra sa valeur τe renforcée d’une quantité 1/f (s), où
20
f (s) est la valeur de la solution, que l’on suppose supérieure à 0.
e
Les premières applications de colonies de fourmis artificielles ne compor-
br
taient que les composantes décrites plus haut. Il a été observé que la rétro-
action positive engendrée par la mise à jour des traces de phéromones artifi-
cielles était difficile à contrôler à l’aide des trois paramètres α, β et ρ. En effet,
em
le point de basculement entre un processus presque aléatoire dépourvu d’ap-
prentissage et un processus qui apprend trop rapidement à répéter la construc-
tion d’une même solution est difficile à trouver.
Pour y remédier, [?] ont suggéré de limiter les valeurs des traces entre deux
pt
éléments voient leur valeur grandir tellement qu’il ne deviendrait plus possible
de construire une solution sans eux.
On aboutit ainsi à la trame présentée plus loin, qui s’est révélée beau-
coup plus efficace que de nombreuses autres trames proposées antérieure-
25
ment. Cette trame comporte une méthode d’amélioration. En effet, les implan-
tations de colonies de fourmis artificielles « pures » basées uniquement sur
la construction de solutions se sont révélées peu performantes et difficiles à
mettre au point. Il peut exister des exceptions, notamment pour le traitement
g
17
4 pour Imax itérations faire
5 pour k = 1 . . . m faire
Initialiser s à une solution partielle triviale;
20
6
7 R ← E; // Éléments pouvant encore être ajoutés à s
8 tant que R 6= ∅ faire
9 Choisir aléatoirement e ∈ R avec probabilité proportionnelle à
τeα · c(s, e)β ;
e
// Formule des colonies de fourmis
10 s ← s ∪ e;
11
pour compléter s;
br
Retirer de R les éléments qui ne peuvent plus être utilisés
em
12 sk ← a(s); // Trouver l’optimum local sk associé à s
13 si f ∗ > f (sk ) alors
14 f ∗ ← f (sk );
15 s∗ ← sk ; // Mettre à jour la meilleure solution connue
pt
16 pour ∀e ∈ E faire
17 τe ← (1 − ρ) · τe ; // Évaporation des traces de phéromones
se
peut être délicat pour certains problèmes de concevoir une fonction de coût
g
C’est pourquoi une trame simplifiée appelée FANT (pour Fast ANT system)
a été proposée. En plus du nombre d’itérations Imax , cette trame ne comporte
T_
225
175
Avant optimisation
150
Après optimisation
125 Amélioration
100
17
75
50
20
25
0
0 100 200 300 400 500 600 700 800
e
Itérations
br
em
F IGURE 7.2 – Comportement de FANT sur un problème de voyageur de commerce à
225 villes. Ce diagramme donne, en fonction du nombre de solutions construites par
la méthode, le nombre d’arêtes qui diffèrent de la meilleure solution trouvée jusque là,
juste après la construction au moyen des traces de phéromones artificielles et après
amélioration à l’aide de la technique de chaînes d’éjections. Les lignes verticales in-
pt
rée. Finalement, FANT incorpore une méthode d’amélioration locales sur les
solutions construites au moyen des traces. Il a en effet été remarqué, pour de
25
Algorithme 7.2 : Trame de FANT. Bien que pouvant paraître complexe, la plus
grande partie de cette trame consiste à adapter automatiquement le poids res-
pectif τc que l’on donne à la solution nouvellement construite en regard du poids
τb associé à la meilleure solution construite et à réinitialiser les traces si l’on a trop
17
appris ou si on vient d’améliorer la meilleure solution connue.
Entrées : ensemble E d’éléments dont sont constituées les solutions du
problème ; fonction utilité f à minimiser, paramètres Imax , τb et
20
méthode d’amélioration a(·)
Résultat : Solution complète s∗
1 f ∗ ← ∞;
τc ← 1;
e
2
3 pour ∀e ∈ E faire
4
5
τe ← τc
pour Imax itérations faire br
em
6 Initialiser s à une solution partielle triviale;
7 R ← E; // Éléments pouvant encore être ajoutés à s
8 tant que R 6= ∅ faire
9 Choisir aléatoirement e ∈ R avec probabilité proportionnelle à τe ;
pt
10 s ← s ∪ e;
11 Retirer de R les éléments incompatibles pour compléter s;
se
21 pour ∀e ∈ E faire
22 τe ← τc ; // Effacer toutes les traces
T_
23 pour ∀e ∈ s0 faire
24 τe ← τe + τc ; // Renforcement des traces de la solution courante
∗
25 pour ∀e ∈ s faire
26 τe ← τe + τb ; // Renforcement des traces de la meilleure solution
166 A PPRENDRE À CONSTRUIRE
17
i n t increment , / ∗ Reinforcement o f l o c a l s o l u t i o n ∗ /
int∗ trail []) / ∗ Pheromone t r a i l s ∗ /
{ int i , j ;
f o r ( i = 0 ; i < n ; i ++)
f o r ( j = 0 ; j < n ; j ++)
20
t r a i l [ i ] [ j ] = increment ;
f o r ( i = 0 ; i < n ; i ++)
t r a i l [ i ] [ i ] = 0;
}
e
void u p d a t e _ t r a i l ( i n t n , / ∗ Problem s i z e ∗/
const i n t p [ ] , /∗ Local s o l u t i o n ∗/
c o n s t i n t best_p [ ] , /∗ Global best s o l u t i o n ∗/
{ i n t i = 0;
i n t ∗increment ,
i n t R,
int∗ trail [])
br / ∗ Reinforcement o f l o c a l s o l u t i o n
/ ∗ Reinforcement o f g l o b a l b e s t s o l u t i o n
∗/
∗/
em
w h i l e ( i < n && p [ i ] == best_p [ i ] ) i ++;
i f ( i == n ) / ∗ New s o l u t i o n p same as b e s t s o l u t i o n found ∗/
{
(∗ i n c r e m e n t ) + + ;
i n i t _ t r a i l ( n , ∗increment , t r a i l ) ;
}
pt
else
f o r ( i = 0 ; i < n ; i ++)
{ t r a i l [ p [ i ] ] [ p [ ( i +1)%n ] ] += ∗ i n c r e m e n t ;
se
Une fois que les trois procédures données dans les codes 7.1 et 7.2 ainsi
25
17
{ sum = 0 ;
f o r ( j = i + 1 ; j < n ; j ++)
sum += t r a i l [ p [ i − 1 ] ] [ p [ j ] ] ;
t a r g e t = u n i f ( 0 , sum−1);
f o r ( j = i , sum = t r a i l [ p [ i − 1 ] ] [ p [ j ] ] ; sum < t a r g e t ; j ++)
20
sum += t r a i l [ p [ i − 1 ] ] [ p [ j + 1 ] ] ;
swap(&p [ j ] , &p [ i ] ) ;
}
return tsp_length (n , d , p ) ;
}
e
Listing 7.3 – tsp_FANT Implantation de la méthode FANT pour le problème du voya-
double tsp_FANT ( i n t n ,
double∗∗ d , br
geur de commerce. La recherche locale utilisée est celle donnée dans le code 4.4
/ ∗ Number o f c i t i e s
/∗ Distance matrix
∗/
∗/
em
i n t best_sol [ ] , / ∗ Out S o l u t i o n r e t u r n e d ∗/
i n t R, / ∗ FANT Parameters ∗/
int iterations ) / ∗ Number o f s o l u t i o n s generated ∗/
{ int i , j , nr_iteration , /∗ Indices ∗/
increment = 1; / ∗ Reinforcement o f elements o f t h e s o l u t i o n generated ∗/
double c o s t ; / ∗ Cost o f s o l u t i o n generated by t h e a n t ∗/
pt
double b e s t _ c o s t = i n f i n i t e ; / ∗ Length o f t h e b e s t t o u r ∗/
int∗ p; / ∗ S o l u t i o n generated by an a n t process ∗/
i n t ∗∗ t r a i l ; / ∗ Pheromone t r a i l m a t r i x ∗/
se
p = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
t r a i l = ( i n t ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ∗ ) ) ;
f o r ( i = 0 ; i < n ; i ++)
t r a i l [ i ] = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
25
i n i t _ t r a i l ( n , increment , t r a i l ) ;
/ ∗ FANT i t e r a t i o n s ∗ /
f o r ( n r _ i t e r a t i o n = 0 ; n r _ i t e r a t i o n < i t e r a t i o n s ; n r _ i t e r a t i o n ++)
{ / ∗ B u i l d a new s o l u t i o n ∗ /
cost = g e n e r a t e _ s o l u t i o n _ t r a i l ( n , d , p , t r a i l ) ;
g
/ ∗ Improve s o l u t i o n w i t h a l o c a l search ∗ /
tsp_LK ( n , d , p , &c o s t ) ;
Al
/ ∗ Best s o l u t i o n improved ? ∗ /
i f ( cost < best_cost − epsilon )
{ best_cost = cost ;
T_
p r i n t f ( "FANT %d %f \ n " , n r _ i t e r a t i o n , c o s t ) ;
f o r ( j = 0 ; j < n ; j ++) b e s t _ s o l [ j ] = p [ j ] ;
increment = 1;
i n i t _ t r a i l ( n , increment , t r a i l ) ;
}
else
/ ∗ Memory update ∗ /
u p d a t e _ t r a i l ( n , p , b e s t _ s o l , &increment , R, t r a i l ) ;
};
free (p ) ;
f o r ( i = 0 ; i < n ; i ++)
free ( t r a i l [ i ] ) ;
free ( t r a i l ) ;
return best_cost ;
} / ∗ tsp_FANT ∗ /
168 A PPRENDRE À CONSTRUIRE
17
20
(a) Une excellente solution (b) Tournées rapidement obtenues
e
F IGURE 7.3 – Une très bonne solution à un problème d’élaboration de tournées de
véhicules 7.3(a) et quelques tournées obtenues rapidement à l’aide de recherches avec
br
tabou sommaires 7.3(b). On remarque une grande similitudes entre ces dernières et
celles de la bonne solution.
em
7.2 Construction de vocabulaire
pt
17
20
e
br
em
pt
se
25
figure 7.6, on a représenté toutes les arêtes présentes dans plus de deux tiers
des 100 tournées obtenues en appliquant une méthode d’amélioration basée
sur un voisinage de chaînes d’éjections démarrant avec une solution aléatoire.
La solution optimale de ce problème étant connue, il a été possible de
mettre en évidence la quinzaine d’arêtes fréquemment obtenues qui ne font
pas partie de la solution optimale. Il est intéressant de remarquer que près
des 80% des arêtes de la solution optimale ont été identifiées en initialisant le
dictionnaire avec une simple méthode d’amélioration.
170 A PPRENDRE À CONSTRUIRE
17
(a) Tentative de solution (b) Solution complétée et améliorée
20
F IGURE 7.5 – Une tentative de phrase est construite en choisissant de façon randomi-
sée des mots du dictionnaire 7.5(a). Cette tentative est complétée et améliorée 7.5(b).
e
Exercices
br
Exercice 7.1. Colonie de fourmis pour l’arbre de Steiner Pour le problème
em
de l’arbre de Steiner, comment définir les traces d’une colonie de fourmis arti-
ficielles ? Décrire comment ces traces sont exploitées.
Exercice 7.2. Ajustement du paramètre de FANT Trouver de bonnes valeurs
pour le paramètre R de la méthode tsp_FANT donnée par le code 7.3 lorsque
pt
17
20
e
br
em
pt
se
25
g
Al
T_
17
20
e
br
em
pt
se
25
g
Al
T_
Chapitre 8
17
Apprendre à modifier
20
Les recherches locales jouent un rôle central dans les méta-heuristiques.
e
Pratiquement toutes les méthodes heuristiques efficaces incorporent ces tech-
niques. Du reste, les méta-heuristiques sont parfois définies comme un pro-
cessus maître guidant une recherche locale.
br
Dans le chapitre 4, nous avons déjà vu quelques techniques d’adaptation
em
d’un voisinage de base, en particulier sa limitation au moyen de liste de mou-
vements candidats ou la recherche granulaire et son extension par une re-
cherche en éventail ou une chaîne d’éjections. Le chapitre 6 consacré aux
méthodes randomisées présente de nombreuses extensions de recherches lo-
pt
les solutions visitées ou de l’effort de calcul pour les comparer aux solutions
voisines.
17
T de m éléments la valeur i dans son entrée h(si ) mod m. Ainsi, la valeur tk
de l’entrée k du tableau T indique à quelle itération une solution dont la valeur
de hachage (modulo m) a été visitée.
20
La fonction h n’est généralement pas bijective sur l’ensemble des solutions
du problème, ce qui fait que plusieurs solutions différentes peuvent avoir la
même valeur de hachage et la taille m du tableau T doit être limitée en rai-
son de la mémoire disponible. Cette technique est donc une approximation du
e
concept d’interdiction de solutions déjà visitées. En effet, non seulement ces
leurs de hachage.
br
dernières sont interdites, mais également toutes celles qui ont les mêmes va-
che avec tabous : la durée des interdiction, parfois appelé longueur de la liste
de tabous.
se
Fonctions de hachage
Le choix d’une fonction de hachage pour implanter une recherche avec ta-
25
bous n’est pas très difficile. Dans certains cas, la valeur de la fonction-objectif
convient parfaitement, notamment lorsque le voisinage comporte de nombreux
mouvements à coût nul. En effet, les modifications neutres rendent l’appren-
tissage difficile, dans la mesure où l’on choisit à chaque itération le meilleur
g
d’itérations de revenir à une valeur donnée permet dans bien des cas de casser
la structure d’un optimum local et d’en découvrir un autre.
Une fonction de hachage générale est la suivante : en considérant qu’une
solution est composée d’éléments e ∈ E, pour reprendre la notation introduite
dans le chapitre 3 sur les méthodes constructives, on peut associer à tous les
éléments e des valeurs entières ze . Ces valeurs sont générées aléatoirement
au début Pde l’algorithme. La valeur de hachage d’une solution s est donnée par
h(s) = e∈s ze .
Une technique de hachage plus sophistiquée, utilisant de multiples tables,
est discutée dans [?]. Elle permet d’obtenir l’équivalent d’une table très grande,
tout en limitant l’espace-mémoire.
8.1 R ECHERCHE AVEC TABOUS 175
17
Pour être concret, prenons l’exemple du voyageur de commerce symé-
trique. Un mouvement 2-opt peut être caractérisé par le couple [i, j] qui con-
siste à remplacer les arêtes [i, si ] et [j, sj ] de la solution courante s par les
arêtes [i, j] et [si , sj ] (on suppose ici que la solution s donne le « successeur »
20
de chaque ville et que la ville j vient « après » la ville i lorsque l’on parcours les
villes dans l’ordre donné par s). Si le mouvement [i, j] est réalisé à une itération,
on peut interdire le mouvement inverse [i, si ] durant les itérations suivantes. Il
s’agit là d’une interdiction directe basée sur l’inverse d’un mouvement.
e
Une autre possibilité est d’interdire de façon indirecte, après avoir effectué
br
le mouvement [i, j], les mouvements qui mèneraient à une solution comportant
à la fois les arêtes [i, si ] et [j, sj ].
em
Par abus de langage, on notera m−1 l’inverse d’un mouvement, ou les ca-
ractéristiques d’une solution que l’on ne voudrait pas retrouver après avoir ef-
fectué le mouvement m d’un voisinage caractérisé par un ensemble M de mou-
vements. Bien que l’on ait (s ⊕ m) ⊕ m−1 = s, il peut y avoir diverses manières
pt
de définir m−1 . La trame de la plus élémentaire des recherche avec tabous est
donnée par l’algorithme 8.1.
se
17
9 s4 (1, 0, 1, 1, 0, 0, 1, 1, 0) 45 40 (5, 11, 9, 12, 0, 0, 6, 10, 8)
10 s5 (1, 0, 1, 1, 1, 0, 1, 1, 0) 49 45 (5, 11, 9, 12, 13, 0, 6, 10, 8)
Tableau 8.1 – Évolution sur 10 itérations d’une recherche avec tabous élémentaire
20
pour le problème du sac de montagne 8.1. La première colonne donne le numéro de
l’itération, la seconde la variable qui est modifiée, la suivante la solution, r donne le
revenu du sac, v le volume occupé et la dernière colonne donne l’état de la liste des
interdictions. Cette recherche interdit de modifier à nouveau une variable pendant 3
e
itérations.
br
Mise en œuvre de l’interdiction sur les mouvements
em
Dans la mesure où la taille du voisinage n’est pas trop importante, il est pos-
sible de mémoriser, pour chaque caractéristique de mouvements, le numéro de
l’itération de la recherche avec tabou à partir de laquelle cette caractéristique
peut à nouveau être utilisée. Donnons tout de suite un petit exemple d’une telle
pt
max r = 12s1 + 10s2 + 9s3 + 7s4 + 4s5 + 8s6 + 11s7 + 6s8 + 13s9
sous 10s1 + 12s2 + 8s3 + 7s4 + 5s5 + 13s6 + 9s7 + 6s8 + 14s9 6 45
contraintes si ∈ {0, 1} (i = 1, . . . , 9)
(8.1)
25
signifie qu’à la première itération, toutes les variables peuvent être modifiées.
Pour ce petit problème, on va supposer que la durée des interdictions est de
3. La solution initiale peut également être initialisée à s = 0, ce qui constitue la
plus mauvaise des solutions admissibles du problème.
Le tableau 8.1 donne l’évolution d’une recherche avec tabous sur ce petit
exemple de problème. Sans surprise, à la première itération, on insère l’objet
9 avec le revenu le plus élevé. À la fin de l’itération 1, il sera interdit jusqu’à
l’itération t9 = 4 = 1 + 3 de poser à nouveau s9 = 0.
Tant que l’on arrive à insérer des objets dans le sac, la recherche avec
tabous se comporte comme un algorithme constructif glouton et arrive à la
solution localement optimale s = (1, 1, 0, 0, 0, 0, 1, 0, 1) de valeur 46 à l’itéra-
8.1 R ECHERCHE AVEC TABOUS 177
tion 4. À l’itération 5, il faut retirer un objet, car le sac est entièrement plein.
Seul l’objet 9 peut être retiré en raison des conditions tabous. La fonction-
objectif diminue donc de 46 à 33, mais de la place est libérée. À l’itération 6,
le meilleur mouvement consisterait à ajouter l’objet 9, mais ce mouvement est
tabou. Il correspondrait du reste à revenir à la solution visitée à l’itération 4.
Le meilleur mouvement autorisé est donc d’ajouter l’objet 3, puis à l’itération
suivante l’objet 8, menant à un nouvel optimum local s = (1, 1, 1, 0, 0, 0, 1, 1, 0)
de valeur r = 48, avec un sac à nouveau entièrement plein avec v = 45. À
l’itération 8, il faut à nouveau retirer un objet et on pose s2 = 0. La place ainsi
17
libérée permet d’ajouter les objets 4 et 5 ce qui permet d’obtenir une solution
s = (1, 0, 1, 1, 1, 0, 1, 1, 0) encore meilleure que les deux optima locaux précé-
demment obtenus.
20
Dans l’exemple du voyageur de commerce, le type d’interdictions décrites
ci-dessus peuvent être implantées au moyen d’une matrice T dont l’entrée tij
donne l’itération à partir de laquelle on peut à nouveau effectuer un mouvement
qui ferait apparaître l’arête [i, j] dans la solution. Ce principe s’étend à tout
e
problème combinatoire pour lequel on cherche une permutation optimale.
tion suivante. Il est clair que la durée maximale des interdictions est limitée par
la taille du voisinage : Pour autant qu’il existe à chaque itération un mouvement
se
Il faut donc apprendre quelle est cette durée pour l’exemple de problème
Al
2650
17
2575
20
2550
0 5 10 15 20 25 30 35 40 45 50
Durée des interdictions
e
F IGURE 8.1 – Influence de la durée des interdictions d’une recherche avec tabous ap-
br
pliquée au problème de l’affectation quadratique où il faut trouver une permutation de
coût aussi faible que possible. Une faible durée permet de visiter en moyenne des solu-
tions de meilleure qualité, mais la recherche n’arrive pas à s’échapper d’optima locaux,
em
ce qui fait que la qualité de la meilleure solution trouvée n’est pas très bonne. Inver-
sement, si la durée des interdictions est très élevée, la qualité moyenne des solutions
visitée diminue, de même que celle des meilleures solutions trouvées. Dans ce cas, un
bon compromis semble être une durée d’interdictions comprise entre 15 et 20, soit un
pt
peu plus grand que la taille des problèmes traités (ou de la racine carrée de la taille du
voisinage).
se
recherche a été capable de découvrir. Il faut donc trouver une durée d’interdic-
tions suffisamment grande pour ne pas cycler, mais aussi faible que possible
25
25
20
15
17
10
20
5
e
0
0 5 10 15 20
br
25 30 35 40
dmin
em
F IGURE 8.2 – Technique du choix aléatoire de la durée des interdictions entre dmin et
dmin +∆ à chaque itération. Un cercle vide indique que la recherche avec tabou n’a pas
été en mesure de trouver systématiquement l’optimum de 500 exemples de problèmes
d’affectation quadratique. La dimension du cercle est proportionnelle au nombre d’op-
pt
tima trouvés. Un disque rempli indique que l’optimum a été systématiquement trouvé.
La dimension du disque est proportionnelle au nombre moyen d’itérations qu’il a fallu
se
teindre l’optimum.
On remarque qu’avec une durée non aléatoire des interdictions (∆ = 0), il
n’a jamais été possible de trouver toutes les solutions optimales, même avec
des durées relativement importantes. Inversement, on remarque qu’avec des
g
de voyageur de commerce pour tout couple (dmin , ∆). On remarque des si-
militudes avec la figure 8.2 ; dans ce cas il semble également qu’une durée
d’interdiction aléatoire mais proportionnelle à la taille du problème soit un bon
compromis.
Condition d’aspiration
∆
2n
3n
2 Qualité
+4%
17
+3%
n +2%
20
+1%
Optimum
n
2
e
0 br dmin
em
0 n 3n
2 n 2 2n
F IGURE 8.3 – Qualité des solutions obtenues avec une recherche tabou où la durée
des interdictions est choisie aléatoirement entre dmin et dmin + ∆ pour un problème de
pt
il sera retenu. Dans la littérature, ceci est appelé une condition d’aspiration.
Il est possible de définir d’autres conditions d’aspiration moins triviales, no-
tamment pour mettre en œuvre une mémoire à long terme.
g
une mémoire à court terme qui peut se révéler très efficace pour des problèmes
de taille modérée. Par contre, dès que l’on s’attaque à des problèmes plus
difficiles, ce seul mécanisme n’est pas suffisant. Une stratégie de recherche
qui a été proposée dans le cadre de la recherche avec tabous est d’alterner
des phases d’intensification et de diversification. Le but de l’intensification est
d’examiner en profondeur une portion limitée de l’espace de recherche, en
conservant des solutions qui ont une structure globalement similaire. Une fois
que l’on pense avoir exploré toutes les solutions intéressantes de cette portion,
il faut aller voir ailleurs. En d’autres termes, on diversifie la recherche en brisant
la structure de la solution sur laquelle on travaille.
L’intensification de la recherche peut être faite à l’aide du simple mécanisme
8.2 O SCILLATIONS STRATÉGIQUES 181
17
la structure d’une solution est d’effectuer des modifications qui n’ont jamais
été retenues durant de nombreuses itérations. Si l’on dispose d’une recherche
avec tabous de base, où l’on a mémorisé pour chaque mouvement l’itération
à partir de laquelle il peut à nouveau être effectué, la mise en œuvre de cette
20
forme de mémoire à long terme est pratiquement gratuite : en effet si le nu-
méro d’itération mémorisé pour un mouvement est bien plus petit que celui
de l’itération courante, c’est que le mouvement en question n’a pas été effec-
e
tué pendant longtemps. On peut ainsi forcer l’utilisation de cette modification,
quelle que soit la qualité de la solution vers laquelle elle conduit.
br
Ce mécanisme requiert un nouveau paramètre, K représentant le nombre
d’itérations à partir duquel un mouvement jamais choisi sera forcé. Naturel-
lement, ce paramètre doit être plus grand que la taille du voisinage, sinon la
em
recherche dégénérera en n’effectuant plus que des mouvements forcés. Cet
type de mémoire à long terme est une autre forme de mise en œuvre des
conditions d’aspiration que l’on a introduite dans la section précédente.
pt
puis un autre ailleurs et ainsi de suite jusqu’à ce qu’il soit à nouveau autorisé
de défaire le premier nœud.
Pour éviter ce comportement, il faut mémoriser le nombre de fois fm qu’un
mouvement m a été choisi et en limiter son usage. Un mécanisme est de pé-
g
∆
2n
3n
2 Qualité
+4%
17
+3%
n +2%
20
+1%
Optimum
n
2
e
0 br dmin
em
0 n 3n
2 n 2 2n
F IGURE 8.4 – Même diagramme que pour la figure 8.3, mais la recherche avec tabous
a été dotée du mécanisme de mémoire à long terme de pénalisation des mouvements
pt
problème.
25
chaque mouvement peut être implanté. Ici, la liste des interdictions est une
matrice dont l’entrée (i, j) donne le numéro de l’itération à partir de laquelle
on peut à nouveau utiliser l’arête [i, j] dans un mouvement. Le comptage des
mouvements est implanté de façon similaire. Le code 8.1 donne la procédure
g
17
{ / ∗ F o r b i d t o use again edges ( i , i +1) and ( j , j +1) ∗ /
t a b u _ l i s t [ s o l [ i ] ] [ s o l [ i + 1 ] ] = t a b u _ l i s t [ s o l [ j ] ] [ s o l [ ( j +1)%n ] ] =
t a b u _ l i s t [ s o l [ i + 1 ] ] [ s o l [ i ] ] = t a b u _ l i s t [ s o l [ ( j +1)%n ] ] [ s o l [ j ] ]
= tabu_duration ;
20
frequency [ s o l [ i ] ] [ sol [ j ]]++;
frequency [ s o l [ j ] ] [ sol [ i ]]++;
frequency [ s o l [ i + 1 ] ] [ s o l [ ( j +1)%n ] ] + + ;
frequency [ s o l [ ( j +1)%n ] ] [ s o l [ i + 1 ] ] + + ;
/ ∗ Update s o l u t i o n c o s t ∗ /
e
∗length = ∗length + d [ sol [ i ] ] [ sol [ j ] ] + d [ s o l [ i + 1 ] ] [ s o l [ ( j +1)%n ] ]
− d [ s o l [ i ] ] [ s o l [ i + 1 ] ] − d [ s o l [ j ] ] [ s o l [ ( j +1)%n ] ] ;
}
f o r ( i n t r = 0 ; r < ( j −i ) / 2 ; r ++)
swap(& s o l [ r + i + 1 ] , &s o l [ j −r ] ) ;
br
/ ∗ Reverse path from i + 1 t o j ∗ /
em
Le code 8.2 donne une implantation d’une recherche avec tabous pour le
problème du voyageur de commerce basé sur un voisinage de type 2-opt. Deux
types de mémoires sont utilisées : une mémoire classique à court terme qui
pt
Listing 8.2 – tsp_TS.c Implantation d’une recherche avec tabous pour le problème du
voyageur de commerce.
v o i d tsp_TS ( i n t n , / ∗ Number o f c i t i e s ∗/
double ∗∗d , / ∗ D i s t a n c e m a t r i x , must be s y m m e t r i c a l ∗/
i n t best_sol [ ] , / ∗ InOut S o l u t i o n p r o v i d e d and r e t u r n e d ∗/
double ∗ b e s t _ c o s t , / ∗ InOut Length o f t h e b e s t t o u r ∗/
int iterations , / ∗ Number o f TS i t e r a t i o n s ∗/
i n t min_tabu_duration , / ∗ Minimal f o r b i d d e n d u r a t i o n ∗/
i n t max_tabu_duration ,
double F )
{ int i , j , /∗ Indices ∗/
nr_iteration , /∗ I t e r a t i o n counter ∗/
i_kept , j_kept ; / ∗ Move t o perform ∗/
17
double d e l t a , / ∗ Cost o f c u r r e n t move ∗/
delta_penalty , / ∗ Cost o f move k e p t w i t h p e n a l t y ∗/
length = ∗best_cost ; / ∗ Length o f c u r r e n t s o l u t i o n ∗/
i n t ∗ s o l = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ; /∗ Current s o l u t i o n ∗/
20
i n t ∗∗ t a b u _ l i s t ; / ∗ I t e r a t i o n from which an a r c can be moved again ∗/
i n t ∗∗ f r e q u e n c y ;
double p e n a l t y ;
f o r ( i = 0 ; i < n ; i ++) s o l [ i ] = b e s t _ s o l [ i ] ;
e
/ ∗ I n i t i a l i z a t i o n o f tabu l i s t : any move i n i t i a l l y allowed ∗/
t a b u _ l i s t = ( i n t ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ∗));
f o r ( i = 0 ; i < n ; i ++)
f o r ( i = 0 ; i < n ; i ++)
br
t a b u _ l i s t [ i ] = ( i n t ∗) c a l l o c ( ( s i z e _ t ) n , s i z e o f (
f r e q u e n c y = ( i n t ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t
f r e q u e n c y [ i ] = ( i n t ∗) c a l l o c ( ( s i z e _ t ) n , s i z e o f (
int ));
∗));
int ));
em
/ ∗ Main tabu search l o o p ∗ /
f o r ( n r _ i t e r a t i o n = 1 ; n r _ i t e r a t i o n <= i t e r a t i o n s ; n r _ i t e r a t i o n ++)
{ delta_penalty = i n f i n i t e ;
i _ k e p t = j _ k e p t = −1; / ∗ Dummy v a l u e t o a v o i d c o m p i l e r warnings ∗ /
pt
t a b u _ l i s t [ s o l [ i + 1 ] ] [ s o l [ ( j +1)%n ] ] <= n r _ i t e r a t i o n | |
length + delta < ∗best_cost − epsilon ) ) /∗ Aspirated ? ∗/
{ delta_penalty = delta + penalty ;
i_kept = i ; j_kept = j ;
}
} ; /∗ f o r j ∗/
g
} ; /∗ f o r i ∗/
Al
s o l , &l e n g t h , t a b u _ l i s t , f r e q u e n c y ) ;
else
p r i n t f ( " A l l moves are f o r b i d d e n : tabu l i s t t o o l o n g \ n " ) ;
f o r ( i = 0 ; i < n ; i ++) f r e e ( f r e q u e n c y [ i ] ) ;
f r e e ( frequency ) ;
f o r ( i = 0 ; i < n ; i ++) f r e e ( t a b u _ l i s t [ i ] ) ;
free ( tabu_list ) ; free ( sol ) ;
} / ∗ tsp_TS ∗ /
8.2 O SCILLATIONS STRATÉGIQUES 185
Redéparts
Une technique fréquemment employée pour intensifier une recherche avec
tabous est de repartir avec la meilleure solution trouvée jusque là si la re-
cherche semble stagner, par exemple s’il n’y a pas eu d’amélioration de cette
meilleures solution durant un relativement grand nombre d’itération. Afin de ne
pas perdre tout le travail fait depuis la dernière amélioration de la meilleure
solution, le redépart se fait en conservant les informations collectées durant
ces itérations, notamment la liste des interdictions et les statistiques sur les
fréquences d’utilisation des mouvement, si ce mécanisme a été utilisé.
17
Ainsi, les structures de données guidant la recherche étant dans un état
différent après redémarrage, la trajectoire suivie par la recherche le sera éga-
lement. Ce mécanisme peut être vu comme l’opposé de celui présenté plus
20
haut où l’on force l’utilisation d’attribut négligés pendant de nombreuses itéra-
tions. Son but est de réaliser une intensification de la recherche et non une
diversification. Naturellement, l’implantation de ce mécanisme implique l’intro-
duction de nouveaux paramètres qu’il faut ajuster, tel que le nombre d’itérations
e
à effectuer avant un redépart et une éventuelle adaptation de la valeur d’autres
br
paramètres (durée des interdictions, pénalité sur la fréquence) pour inciter la
recherche à suivre des trajectoires diversifiées.
em
pt
se
25
g
Al
T_
186 A PPRENDRE À MODIFIER
Exercices
Exercice 8.1. Recherche tabou pour une fonction explicite Une fonction
entière [−7, −6] × [6, 7] → [−10, 650] est donnée explicitement dans le tableau
suivant :
x
y -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
-6 248 216 210 222 230 234 256 304 336 372 428 495 585 650
-5 193 175 157 166 174 181 215 249 295 329 382 454 539 597
-4 138 144 126 116 124 150 184 194 250 305 361 425 480 566
-3 123 89 85 97 105 109 129 179 209 246 302 368 458 525
17
-2 92 58 70 70 78 94 98 148 168 223 282 339 413 510
-1 68 34 46 46 54 70 74 124 144 199 258 315 388 486
0 51 17 14 25 33 38 57 107 136 174 230 296 386 454
1 18 25 5 -4 3 29 65 74 131 185 240 305 361 445
2 27 6 -10 0 8 13 46 83 126 160 213 284 371 429
20
3 33 0 -3 7 15 20 39 89 118 156 212 278 368 436
4 33 12 -4 6 14 19 52 89 132 166 219 290 377 435
5 30 37 17 7 15 41 77 86 143 197 252 317 373 457
6 69 35 32 43 51 56 75 125 154 192 248 314 404 472
7 92 58 70 70 78 94 98 148 168 223 282 339 412 510
Pour trouver le minimum de cette fonction, on considère une recherche avec
e
tabous basée sur le voisinage où l’on modifie d’une unité la valeur d’une va-
br
riable. Les conditions tabous consistent à interdire pendant d itérations d’incré-
menter (respectivement : de décrémenter) une variable qui a été décrémentée
(respectivement : incrémentée). Considérer tout d’abord une durée d’interdic-
em
tions de d = 3 et (−7, 7) comme solution de départ. Ensuite, partir de (−7, −6)
et utiliser t = 1. La recherche s’arrête s’il n’y a plus de mouvement autorisé ou
si 25 itérations ont été effectuées.
pt
deux clients appartenant à des tournées différentes. Proposer des critères ta-
bous pour ce voisinage.
Exercice 8.3. Application de recherche avec tabous au QAP Soit l’exemple
25
F = 1 2 0 1 2
D= 2 3 0 0 0
2 1 1 0 1 4 0 0 0 5
Al
3 2 2 1 0 1 2 0 5 0
T_
17
Apprendre par les solutions
20
Avec une modélisation adéquate d’un problème d’optimisation, il est très
e
facile de construire de nombreuses solutions différentes notamment au moyen
de méthodes randomisées. Dès lors, on peut chercher à apprendre à élabo-
br
rer de nouvelles solutions à partir celles précédemment construites. Dans ce
chapitre, nous allons donc étudier comment exploiter une population de so-
lutions, ce qui nous amènera à rassembler dans une méthode de résolution
em
les diverses briques des méta-heuristiques étudiées plus haut (construction et
améliorations de solutions).
pt
permanent. Cela signifie que les êtres vivants peuvent optimiser leur chances
de survie, donc de résoudre des problèmes d’une extrême complexité. Dès
lors, pourquoi ne pas tenter de reproduire artificiellement cette évolution pour
résoudre des problèmes d’optimisation combinatoire.
g
Population
17
20
Évaluation
e
Mutation Croisement
br
F IGURE 9.1 – Boucle générationnelle dans un algorithme évolutionnaire. À partir d’une
em
population de solutions, symbolisées ici par des bâtonnets de couleur, on sélectionne
des individus qui vont se reproduire par croisement et mutation. Les solutions-filles
ainsi générées sont évaluées et incorporées à la population. Finalement, des individus
de cette dernière sont éliminés par un opérateur de sélection pour la survie afin de
pt
2 répéter
3 Sélectionner des individus de P pour la reproduction;
4 Combiner ces individus et leur faire subir des mutations pour obtenir
λ nouvelle solutions;
5 Sélectionner µ individus pour la survie, parmi les µ + λ ; ces µ
individus formeront la population P de la prochaine génération
6 jusqu’à ce qu’un critère d’arrêt soit satisfait;
9.2 A LGORITHMES GÉNÉTIQUES 189
la meilleure) au détriment de ceux qui sont plus faibles, malades, mal adaptés,
à l’image de ce qui se passe dans la nature.
Les individus sélectionnés sont ensuite mélangés entre eux (par exemple
par paire) à l’aide d’un opérateur de croisement pour former λ nouvelles so-
lutions appelées solutions-enfants qui subissent au passage des modifications
aléatoires au moyen d’un opérateur de mutation. Ces deux opérateurs simulent
la reproduction sexuée des espèces vivantes, en faisant l’hypothèse qu’avec
un peu de chance les bonnes caractéristiques (les bons gènes contenus dans
l’ADN) des solutions-parents seront transmises à leurs enfants et que les mu-
17
tations aléatoires feront apparaître de nouveaux gènes favorables.
Finalement, les nouvelles solutions sont évaluées et un opérateur de sé-
lection pour la survie élimine λ solutions parmi les µ + λ disponible pour se
20
ramener à une nouvelle population de µ individus. La figure 9.1 illustre le pro-
cessus d’une boucle générationnelle.
La trame des algorithmes évolutionnaires laisse une grande liberté dans
les choix à opérer pour l’implantation des divers opérateurs et paramètres. Par
e
exemple [?] dans sa « stratégie d’évolution » n’utilise pas d’opérateur de croise-
ment entre deux solutions. Dans cette technique, les solutions de la population
br
sont modifiées avec un opérateur de mutation et entrent en compétition les
unes avec les autres, à l’instar de la reproduction parthénogénétique.
em
9.2 Algorithmes génétiques
pt
paradoxal, puisque dans cette référence, le cœur de l’étude portait sur la com-
préhension des mécanismes de convergence de ces algorithmes, et non sur
leur capacité à optimiser des problèmes difficiles. Pendant longtemps, la com-
munauté gravitant dans ce domaine a continué à travailler sur la théorie de
25
17
individu aura une mesure de qualité fi = (1 − rµi )p , où p > 0 est un paramètre
permettant de moduler la pression sélective. Une pression p = 0 implique un
tirage uniforme parmi la population (pas de pression sélective) alors que p = 2
20
représente une assez forte pression. Le code 9.1 donne une implantation de
cet opérateur pour une pression sélective p = 1.
Listing 9.1 – rank_based_selection.c Implantation d’un opérateur de sélection pour
la reproduction basé sur le rang, avec une pression sélective p = 1. Le meilleur des µ
e
individus a une probabilité de 2µ/(µ · (µ + 1))) d’être choisi, alors que le plus mauvais
a une probabilité de 2/(µ · (µ + 1))).
br
i n t r a n k _ b a s e d _ s e l e c t i o n ( i n t mu)
{ r e t u r n mu − ( i n t ) c e i l ( s q r t ( 0 . 2 5 + 2∗ u n i f ( 1 , mu∗(mu+ 1 ) / 2 ) ) − 0 . 5 ) ; }
em
Sélection proportionnelle
L’opérateur de sélection le plus simple est de tirer aléatoirement un individu
pt
proportionnellement
P à sa mesure de qualité. L’individu i aura donc une proba-
bilité fi / fi d’être choisi. En principe, on ne sélectionne pas qu’un individu à
se
mutation.
T_
Sélection naturelle
Il est également envisageable de procéder à une sélection totalement aléa-
toire et uniforme pour la reproduction, à l’image de ce qui se passe pour les
êtres vivants. La convergence de l’algorithme doit alors être guidée par l’opé-
rateur de sélection pour la survie, qui assure un biais vers les meilleures indi-
vidus.
Sélection complète
Si l’on ne choisit pas une taille de population trop grande, il est également
possible de faire participer tous les individus, de manière systématique, pour
9.2 A LGORITHMES GÉNÉTIQUES 191
17
20
e
F IGURE 9.2 – Reproduction de coccinelles. On peut imaginer que les couples du haut
vont produire des enfants très semblables à eux-mêmes, alors que ceux du bas vont
br
conserver une certaine diversité génétique dans la population, et, avec un peu de
chance, produire quelques enfants nettement mieux adaptés qu’eux à leur environne-
ment.
em
la reproduction. Comme pour la sélection naturelle, l’évolution de la popula-
tion vers de bonnes solutions dépend alors de l’opérateur de sélection pour la
pt
façon aléatoire.
Le but de cet opérateur est de créer un nouvel individu, différent de ses pa-
rents, mais ayant hérité de certaines de ses caractéristiques. Avec un peu de
T_
Parent 1
A C A G T T T G G C A G
Parent 2
C A T G C T G A G T C G
Enfant 1
C A A G T T G A G C C G
Enfant 2
17
A C T G C T T G G T A G
20
enfant est tiré aléatoirement du premier ou du second parent, le deuxième enfant rece-
vant l’élément complémentaire.
e
ou plusieurs enfants presque identique à lui-même, où seules des mutations
br
spontanées font évoluer le patrimoine génétique de la population.
em
Croisement uniforme
sante, et qu’à tout vecteur de cette taille on peut faire correspondre une solution
admissible du problème.
Ceci n’est pas le cas pour un problème où l’on cherche une permutation
de n éléments. Une technique d’adaptation du croisement uniforme pour cette
g
parents, pour autant que l’élément en question n’ait pas encore été choisi. Si
les deux parents ont des éléments déjà choisis à la position à remplir, cette
T_
Croisement 1-point
Parent 1
9 5 12 8 2
6 1 11 3 10 4 7
Parent 2
3 11 4 8 12 5 1 2 10 6 9 7
Phase 1
3 11 12 8 2 6 1 10 9 7
Enfant
17
3 11 12 8 2 6 1 5 10 4 9 7
F IGURE 9.4 – Croisement uniforme sur une permutation. Un enfant est produit en deux
passes. On choisit d’abord successivement les éléments de l’un ou l’autre des parents,
20
pour autant que ces éléments n’aient pas été choisis (auquel cas, on prend l’unique élé-
ment disponible, ou on laisse la case vide si les deux éléments des parents figurent déjà
dans l’enfant. La seconde passe complète aléatoirement l’enfant à l’aide des éléments
encore disponibles.
e
Parent 1 br
em
A C A G T T T G G C A G
Parent 2
C A T G C T G A G T C G
pt
Enfant 1
se
A C A G T T T G G T C G
Enfant 2
C A T G C T G A G C A G
25
Croisement 2-points
Parent 1
A C A G T T T G G C A G
Parent 2
C A T G C T G A G T C G
Enfant 1
A C A G C T G A G C A G
17
Enfant 2
C A T G T T T G G T C G
20
jumeaux » par tirage aléatoire de deux points points de croisement (ici, 4 et 8). Copie
des éléments d’un parent jusqu’au premier point de croisement et à partir du second
point de croisement et ceux de l’autre pour la partie intermédiaire.
e
Croisement OX
br
Pour chaque problème, on peut inventer un opérateur de croisement spé-
em
cifique. Par exemple, pour celui du voyageur de commerce, on peut faire le
raisonnement que des portions de chemins devraient être copiées dans la
solution-enfant à partir des parents. Si l’on considère une solution comme étant
une permutation des villes, on se rend compte que le croisement uniforme vu
pt
plus haut et adapté au cas des permutations n’a pas vraiment de sens : la ville
de départ n’est pas déterminante, donc ce n’est pas la position d’une ville dans
se
la tournée qui est importante, mais les villes qui lui précèdent et succèdent.
Dans ce cas, on aura intérêt à adapter l’opérateur de croisement 2-points pour
les problèmes où ce sont les séquences qui sont importantes.
L’opérateur de croisement OX imaginé pour le problème du voyageur de
25
Parent 1
9 5 12 8 2
6 1 11 3 10 4 7
Parent 2
3 11 4 8 12 5 1 2 10 6 9 7
Enfant
2 6 1 11 4 8 12 5 10 9 7 3
17
F IGURE 9.7 – Croisement OX, spécifiquement développé pour le problème du voyageur
de commerce. On commence par tirer aléatoirement deux points de croisement. La
portion intermédiaire du premier parent est copiée chez l’enfant. Dans cet exemple
cette portion se termine par la ville 11. On repère cette ville dans le second parent
20
et on complète l’enfant à partir de là, en l’occurrence par la ville 4. Les villes figurant
déjà dans l’enfant (1, 2, puis 6) sont sautées. Lorsqu’on arrive à la « dernière » ville du
second parent (7) on revient à la première (3).
e
br
Listing 9.2 – OX_crossover.c Implantation de l’opérateur de croisement OX préservant
em
l’ordre d’une sous-séquence
v o i d OX_crossover ( i n t n , / ∗ Number o f c i t i e s ∗/
c o n s t i n t p1 [ ] , c o n s t i n t p2 [ ] , / ∗ Parent s o l u t i o n s ∗/
int child [ ] )
{ i n t i , j , temp ,
pt
a l r e a d y _ i n s e r t e d = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
f o r ( i = 0 ; i < n ; i ++)
already_inserted [ i ] = 0;
p o i n t 1 = u n i f ( 1 , n−3);
i f ( p o i n t 1 >= p o i n t 2 )
{ temp = p o i n t 2 ;
p o i n t 2 = ++ p o i n t 1 ;
p o i n t 1 = temp ;
}
g
/ ∗ Copy t h e p o r t i o n o f p1 a t t h e b e g i n n i n g o f c h i l d ∗ /
Al
/ ∗ Find i n p2 t h e l a s t element i n s e r t e d i n c h i l d ∗ /
f o r ( j = 0 ; p2 [ j ] ! = c h i l d [ p o i n t 2−p o i n t 1 ] ; j + + ) { } /∗ j u s t increment j ∗/
/ ∗ I n s e r t t h e r e m a i n i n g elements i n c h i l d , i n t h e o r d e r o f appearance i n p2 ∗ /
n r _ i n s e r t e d = p o i n t 2−p o i n t 1 + 1 ;
f o r ( i = j ; n r _ i n s e r t e d < n ; i ++)
i f ( ! a l r e a d y _ i n s e r t e d [ p2 [ i%n ] ] )
{ c h i l d [ n r _ i n s e r t e d ] = p2 [ i%n ] ;
a l r e a d y _ i n s e r t e d [ p2 [ i%n ] ] = 1 ;
n r _ i n s e r t e d ++;
}
free ( already_inserted ) ;
} / ∗ OX_crossover ∗ /
196 A PPRENDRE PAR LES SOLUTIONS
Coût
135
17
115
110
20
105
e
100
100 1’000 10’000 100’000
br Générations
em
F IGURE 9.8 – Influence du taux de mutation sur la qualité des solutions produites
en fonction du nombre de boucles générationnelles réalisées. Seule la valeur de la
meilleure solution de la population est représentée, en fonction du nombre d’itérations
de la boucle générationnelle. L’algorithme favorise les meilleures solutions de la po-
pulations au moyen d’opérateurs de sélection pour la reproduction et de survie. Sans
pt
mutation, la population converge relativement rapidement vers des individus tous sem-
blables et de mauvaise qualité. Plus le taux de mutation est élevé, plus la convergence
se
consacré.
L’opérateur de mutation a deux rôles : tout d’abord, la modification locale
Al
peut améliorer la solution, ensuite, même si cette dernière n’est pas amélio-
rée, elle ralentit la convergence globale de l’algorithme par renforcement de
T_
Listing 9.3 – mutate.c Implantation d’un opérateur de mutation pour des problèmes
sur des permutations.
/ ∗ M u t a t i o n o p e r a t o r f o r p e r m u t a t i o n perform m u t a t i o n _ r a t e ∗n / 2 random swaps ∗/
v o i d mutate ( i n t n , /∗ Permutation size ∗/
double m u t a t i o n _ r a t e ,
int p []) / ∗ P e r m u t a t i o n t o mutate ∗ /
{ i n t i , nr_mutations ;
nr_mutations = ( i n t ) ( mutation_rate ∗ n / 2 . 0 ) ;
f o r ( i = 0 ; i < n r _ m u t a t i o n s ; i ++)
swap(&p [ u n i f ( 0 , n −1)] , & p [ u n i f ( 0 , n − 1 ) ] ) ;
}
17
Le code 9.3 donne une implantation d’un opérateur de mutation pour des
problèmes sur des permutations.
20
9.2.4 Sélection pour la survie
Le dernier opérateur-clé des algorithmes est celui de la sélection pour la
e
survie, qui a pour but de ramener la population à sa taille initiale de µ indivi-
dus, après que l’on ait généré λ nouvelles solutions. Plusieurs politiques de
br
sélection ont été imaginées, selon les valeurs choisies pour les paramètres µ
et λ.
em
Remplacement générationnel
Stratégie d’évolution
g
Remplacement stationnaire
140
Coût
µ=5
135
µ = 10
130
125
17
120
115
µ = 20
20
110
µ = 50
105 µ = 100
e
µ = 200
µ = 500
100
10 100
br 1’000 10’000 100’000
em
Générations
F IGURE 9.9 – Influence de la taille de la population sur la qualité des solutions. Lorsque
la population est trop petite, elle converge très rapidement avec peu de chances d’at-
pt
teindre de bonnes solutions. Inversement, une grande population convergera très len-
tement, mais vers de meilleures solutions.
se
Remplacement élitiste
25
17
i n t child_rank = 0;
int i ;
20
/ ∗ Find t h e rank o f t h e c h i l d ∗ /
child_rank = 0;
f o r ( i = 0 ; i < pop_size ; i ++)
i f ( length [ i ] < length_child )
c h i l d _ r a n k ++;
e
i f ( c h i l d _ r a n k < pop_size −1) / ∗ The c h i l d i s n o t dead−born ∗/
br
{ / ∗ May t h e c h i l d be i d e n t i c a l t o an i n d i v i d u a l o f t h e p o p u l a t i o n ?
i f ( f a b s ( l e n g t h [ o r d e r [ c h i l d _ r a n k ] ] − l e n g t h _ c h i l d ) > e p s i l o n &&
( c h i l d _ r a n k == 0 | |
f a b s ( l e n g t h [ o r d e r [ c h i l d _ r a n k −1]] − l e n g t h _ c h i l d ) > e p s i l o n ) )
∗/
em
{ / ∗ The c h i l d n o t p r e s e n t i n t h e p o p u l a t i o n , i t r e p l a c e w o r s t i n d i v i d u a l ∗/
f o r ( i = 0 ; i < n ; i ++)
p o p u l a t i o n [ o r d e r [ pop_size − 1 ] ] [ i ] = c h i l d [ i ] ;
l e n g t h [ o r d e r [ pop_size −1]] = l e n g t h _ c h i l d ;
f o r ( i = 0 ; i < pop_size ; i ++)
i f ( rank [ i ] >= c h i l d _ r a n k )
pt
rank [ i ] + + ;
rank [ o r d e r [ pop_size −1]] = c h i l d _ r a n k ;
f o r ( i = 0 ; i < pop_size ; i ++)
se
o r d e r [ rank [ i ] ] = i ;
}
else
c h i l d _ r a n k = pop_size ; / ∗ C h i l d a l r e a d y p r e s e n t i n p o p u l a t i o n , i g n o r e i t ∗/
}
25
return child_rank ;
}
g
être améliorée par une simple modification locale, telle que vue au chapitre 4.
Secondement, la diversité de la population décroît au fur et à mesure des ité-
rations de la boucle générationnelle pour finalement ne comporter que des
clones du même individu.
Pour contrer ces deux inconvénients, [?] a imaginé ce qu’il a appelé les
algorithmes mémétiques. Le premier de ces défaut est résolu par l’application
d’une recherche locale après avoir produit une solution-enfant. La manière la
plus simple d’éviter la duplication d’individus dans la population est d’éliminer
immédiatement les doublons.
Le code 9.5 illustre une implantation simple d’un algorithme mémétique
pour le problème du voyageur de commerce où les enfants sont améliorés à
200 A PPRENDRE PAR LES SOLUTIONS
l’aide d’une recherche locale basée sur les chaînes d’éjections et où ils ne rem-
placent la plus mauvaise solution de la population que s’ils sont de meilleure
qualité que cette dernière et que leur évaluation est différente de toutes celles
de la population, ce qui assure de ne pas créer de doublons. Cet algorithme
n’implante qu’une version élémentaire d’un algorithme mémétique. Une ges-
tion de la population proposée par [?] est un peu plus évoluée. Ces auteurs
suggèrent de d’évaluer, pour chaque solution produite, une mesure de simila-
rité avec les solutions contenues dans la population. Les solution trop similaires
sont écartées pour maintenir une diversité suffisante, pour que l’algorithme ne
17
converge pas prématurément.
20
plus sophistiquées ont été imaginées, notamment celle basée sur des îlots.
L’idée est de faire évoluer plusieurs populations de faible taille indépendam-
ment pendant un nombre restreint d’itérations. Du fait de leur faible taille, la
qualité des solutions de chaque population s’améliore rapidement. Pour ne pas
e
perdre de diversité, des individus sont transférés régulièrement d’une popula-
tion à une autre, ce qui atténue les problèmes de consanguinité.
br
em
pt
se
25
g
Al
T_
9.3 A LGORITHMES MÉMÉTIQUES 201
17
i n t generations , / ∗ Number o f c h i l d r e n generated ∗/
double m u t a t i o n _ r a t e )
{ i n t i , j , nr_generation , /∗ Indices ∗/
p1 , p2 ; / ∗ Rank o f p a r e n t s chosen ∗/
double l e n g t h _ c h i l d ; /∗ Fitness of c h i l d ∗/
20
i n t ∗ c h i l d = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ; / ∗ c h i l d −s o l u t i o n ∗/
i n t ∗ order ; / ∗ I n v e r s e o f t h e rank ∗/
double∗ l e n g t h = ( double ∗) m a l l o c ( ( s i z e _ t )mu ∗ s i z e o f ( double ) ) ; / ∗ Sol . f i t . ∗/
i n t ∗∗ p o p u l a t i o n ; /∗ Population of sol ut ion s ∗/
o r d e r = ( i n t ∗) m a l l o c ( ( s i z e _ t )mu ∗ s i z e o f ( i n t ) ) ;
e
p o p u l a t i o n = ( i n t ∗∗) m a l l o c ( ( s i z e _ t )mu ∗ s i z e o f ( i n t ∗ ) ) ;
f o r ( i = 0 ; i < mu; i ++)
/ ∗ Generate a p o p u l a t i o n o f random s o l u t i o n s ∗ /
f o r ( i =0; i < mu; i ++) br
p o p u l a t i o n [ i ] = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
em
{ generate_random_permutation ( n , p o p u l a t i o n [ i ] ) ;
length [ i ] = tsp_length (n , d , population [ i ] ) ;
}
/ ∗ Order t h e i n d i v i d u a l by d e c r e a s i n g f i t n e s s and a s s i g n them a rank ∗ /
f o r ( i = 0 ; i < mu; i ++)
order [ i ] = i ;
pt
{ / ∗ Generate a new c h i l d ∗ /
p1 = r a n k _ b a s e d _ s e l e c t i o n (mu ) ;
p2 = r a n k _ b a s e d _ s e l e c t i o n (mu ) ;
OX_crossover ( n , p o p u l a t i o n [ o r d e r [ p1 ] ] , p o p u l a t i o n [ o r d e r [ p2 ] ] , c h i l d ) ;
mutate ( n , m u t a t i o n _ r a t e , c h i l d ) ;
length_child = tsp_length (n , d , child ) ;
g
tsp_LK ( n , d , c h i l d , & l e n g t h _ c h i l d ) ;
i f ( i n s e r t _ c h i l d ( n , c h i l d , l e n g t h _ c h i l d , mu, p o p u l a t i o n , l e n g t h , o r d e r ) == 0 )
Al
p r i n t f ( "GA %d %f \ n " , n r _ g e n e r a t i o n , l e n g t h _ c h i l d ) ;
}
/ ∗ Return t h e b e s t s o l u t i o n o f t h e p o p u l a t i o n ∗ /
f o r ( i = 0 ; i < n ; i ++)
T_
free ( child ) ;
free ( length ) ;
free ( order ) ;
f o r ( i =0; i < mu; i ++)
free ( population [ i ] ) ;
free ( population ) ;
return tsp_length (n , d , best_sol ) ;
}
202 A PPRENDRE PAR LES SOLUTIONS
17
F IGURE 9.10 – La recherche par dispersion brise le tabou d’une reproduction limitée au
20
croisement de deux solutions.
e
La recherche par dispersion est presque aussi ancienne que les algo-
br
rithmes génétiques. [?] a proposé cette technique dans le cadre de la pro-
grammation linéaire en nombres entiers. À l’époque, il brisait certains tabous,
em
comme le fait de pouvoir représenter une solution sous une forme naturelle et
non codée par un vecteur binaire ou de croiser plus de deux solutions entre
elles, comme illustré en figure 9.10.
Les principales idées de la recherche par dispersion sont les suivantes,
pt
sons possibles des individus de la population, qui doit donc être limitée à
quelques dizaines de solutions.
Opérateur de réparation/amélioration Du fait de la représentation naturelle
des solutions, le « croisement » simultané de combinaisons de plusieurs
individus ne produit pas forcément une solution réalisable. Un opérateur
de réparation projetant une solution potentielle irréalisable dans l’espace
des solutions admissibles est donc prévu. Cet opérateur peut également
au passage améliorer une solution réalisable, notamment à l’aide d’une
recherche locale.
Gestion de la population Une population de référence, de faible taille, est dé-
composée en un sous-ensemble de solutions-élites et d’autres solutions
9.4 R ECHERCHE PAR DISPERSION 203
17
ment générer des solutions aussi dispersées que possible, ni comment com-
biner plusieurs solutions entre elles ou encore comment implanter la répara-
tion et l’amélioration des solutions potentielles. Toutes ces choses sont dépen-
20
dantes du problème traité et doivent être adaptées au cas par cas. Le para-
graphe suivant donne un exemple d’adaptation pour le problème du sac de
montagne.
e
Algorithme 9.2 : Trame d’une recherche par dispersion.
des solutions-élites
Résultat : Meilleure solution trouvée
br
Données : Taille µ de la population complète, E taille du sous-ensemble
em
1 Générer de façon systématique une (grande) population P de solutions
potentielles aussi dispersées que possible;
2 répéter
pt
potentielles;
8 Joindre les solutions potentielles à l’ensemble de référence pour
Al
max r = 11s1 + 10s2 + 9s3 + 12s4 + 10s5 + 6s6 + 7s7 + 5s8 + 3s9 + 8s10
sous 33s1 + 27s2 + 16s3 + 14s4 + 29s5 + 30s6 + 31s7 + 33s8 + 14s9 + 18s10 6 100
contraintes si ∈ {0, 1}(i = 1, . . . , 10)
(9.1)
204 A PPRENDRE PAR LES SOLUTIONS
17
8 (0,1,1,0,1,1,0,1,1,0) 43 (0,1,1,1,1,0,0,0,1,0) 44
9 (0,1,1,1,0,1,1,1,0,1) 57 (0,1,1,1,0,0,0,0,1,1) 42 = solution 1
10 (0,1,1,1,1,0,1,1,1,1) 64 (0,1,1,1,0,0,0,0,1,1) 42 = solution 1
20
Tableau 9.1 – Population initiale dispersée P pour le problème du sac de montagne 9.1
et résultat de l’application de l’opérateur de réparation/amélioration sur les solutions
potentielles. Celles qui ne sont pas réalisables sont mises en évidence, de même que
les E = 3 solutions-élites.
e
br
Population initiale Les solutions de ce problèmes sont donc des vecteurs bi-
naires de 10 composantes. Pour générer un ensemble de solutions potentielles
em
aussi dispersées que possible, on peut choisir de mettre tous les objets dans
le sac, ou un sur deux, ou un sur trois, etc. et, pour chaque solution potentielle
ainsi générée, la solution complémentaire. Naturellement, toutes ces solutions
ne sont pas admissibles. En particulier, la solutions comportant tous les objets
pt
tions potentielles. On peut procéder ainsi : tant que la solution n’est pas admis-
sible, retirer l’objet ayant le plus mauvais rapport revenu/volume. Une solution
admissible peut être améliorée gloutonnement, en lui ajoutant l’objet avec le
25
Tableau 9.2 – Détermination des solutions de la population qui sont aussi différentes
que possible des solutions-élites. Si l’on veut un ensemble de référence de µ = 5
17
solutions, on retiendra les solutions 3 et 7 en plus des 3 solutions-élites car ce sont
celles qui maximisent la plus petite distance avec une de ces dernières.
20
38
qu’elle est bonne. Une idée est donc d’attribuer un poids de 38+36+44 à la so-
36 44
lution 3, de 38+36+44 à la solution 7 et de 38+36+44 à la solution 8. Le vecteur
ainsi obtenu est arrondi pour le projeter vers des valeurs binaires :
e
0.322·(1, 0, 0, 1, 0, 0, 1, 0, 0, 1)+
0.305·(0, 1, 0, 1, 0, 1, 0, 0, 0, 1)+
0.373·(0, 1, 1, 1, 1, 0, 0, 0, 1, 0) br
em
=(0.322, 0.678, 0.373, 1.000, 0.373, 0.305, 0.322, 0.000, 0.373, 0.627)
Arrondi(0, 1, 0, 1, 0, 0, 0, 0, 0, 1)
pt
trouvées par une recherche avec tabous. Parmi ces solutions, on en choisit
deux, qui ont été reliées entre elles avec la recherche tabou par un chemin
allant de solution en solution voisine et de les relier à nouveau par un nouveau
chemin, plus court.
g
ture de voisinage pour l’implanter. On choisit donc dans la population une solu-
tion de départ et une solution d’arrivée. On évalue tous les voisins de la solution
T_
17
20
Solution de départ
e
br
em
pt
se
Solution d’arrivée
25
17
double c o s t = ∗ b e s t _ c o s t ; / ∗ Cost o f c u r r e n t s o l u t i o n ∗/
i n t ∗ s u c c _ t a r g e t , ∗succ , ∗pred ; /∗ Alternate sol ut io ns representation ∗/
s u c c _ t a r g e t = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
succ = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
20
pred = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
e
}
do
{ i = 0;
best_delta = i n f i n i t e ;
do br
em
{ i f ( succ [ i ] ! = s u c c _ t a r g e t [ i ] ) / ∗ t r y i −succ [ i ] −> i −s u c c _ t a r g e t [ i ] ∗ /
{ j = pred [ s u c c _ t a r g e t [ i ] ] ; / ∗ j −succ [ j ] n o t i n t a r g e t s o l u t i o n ∗ /
f o r ( k = s u c c _ t a r g e t [ i ] ; k ! = i ; k = succ [ k ] )
i f ( succ [ k ] ! = s u c c _ t a r g e t [ k ] )
{ d e l t a = d [ i ] [ succ [ j ] ] + d [ j ] [ succ [ k ] ] + d [ k ] [ succ [ i ] ]
−d [ i ] [ succ [ i ] ] − d [ j ] [ succ [ j ] ] − d [ k ] [ succ [ k ] ] ;
pt
}
}
}
i = succ [ i ] ;
}
w h i l e ( b e s t _ d e l t a >= 0 . 0 && i ) ; / ∗ Impr . move found o r a l l moves e v a l u a t e d ∗ /
25
sol [ 0 ] = 0;
f o r ( i = 0 ; i < n−1; i ++)
s o l [ i +1] = succ [ s o l [ i ] ] ;
}
T_
}
}
while ( best_delta < i n f i n i t e ) ;
f r e e ( s u c c _ t a r g e t ) ; f r e e ( succ ) ; f r e e ( pred ) ;
} /∗ tsp_path_relinking ∗/
9.5.1 GRASP-PR
Une méthode utilisant les principales composantes des méta-heuristiques
(construction, recherche locale et gestion d’une population de solution) tout en
restant relativement simple et avec peu de paramètres est la méthode GRASP-
PR (procédure de recherche gloutonne adaptative avec chemin de liaison) de
208 A PPRENDRE PAR LES SOLUTIONS
17
tion de P qui est la plus différente d’elle-même tout en étant plus mauvaise.
L’algorithme 9.3 donne la trame de GRASP-PR.
20
Algorithme 9.3 : Trame de GRASP-PR.
Données : Structure de voisinage N , paramètres Imax , 0 6 α 6 1, µ
Résultat : Meilleure solution trouvée
1 P ← ∅;
e
2 tant que |P | < µ faire
3
4
5 si s0 ∈
/ P alors br
Générer une solution s avec une méthode gloutonne biaisée avec α;
s0 ← minimum local associé à s et N ;
em
6 P ← P ∪ s0
7 pour Imax itérations faire
8 Générer une solution s avec une méthode gloutonne biaisée avec α;
pt
11
meilleure solution s∗ du chemin;
12 si s∗ ∈/ P et s∗ strictement meilleur qu’une solution de P alors
13 s remplace la solutions de P la plus différente de s∗ et plus
∗
mauvaise que s∗
25
g
Les essaims particulaires sont un peu particuliers car ils ont été conçus
T_
17
renciant par le graphe d’influence et les formules utilisées pour calculer les
déviations de la vitesse des particules. Dans sa version la plus simple, une
particule p d’un ensemble en comportant P est influencée par seulement deux
solutions : la meilleure solution →
−
20
g trouvée par l’ensemble des particules et la
−
→
meilleure solution mp qu’elle a elle-même trouvée. La nouvelle vitesse de la
particule est un vecteur dont chaque composante est modifiée avec des poids
aléatoirement tirés entre 0 et φ1 dans la direction de −m→ et tirés entre 0 et φ
p 2
→
−
e
dans la direction de g , où φ1 et φ2 sont des paramètres de la méthode. De
plus, une particule a une inertie ω également donnée en paramètre. L’algo-
br
rithme 9.4 donne une trame simple d’essaim particulaire.
em
Algorithme 9.4 : Trame d’une méthode à essaim particulaires.
Données : Fonction f : [→ −
x min , →
−
x max ] ∈ Rn → R à minimiser,
paramètres P, ω, φ1 , φ2 , Imax
Résultat : →−
g
pt
∗
1 f = ∞;
2 pour p = 1 . . . P faire
→
− −−−→ −
sp ← unif (→ x min , →
−
se
3 x max );
4
−→ →
−
mp ← sp ;
→
− −−−→ −
5 vp ← unif (→ x min − → −
x max , →
−
x max − →−x min );
6
∗ →
−
si f > f ( sp ) alors
25
7 f ∗ ← f (→ −
sp );
8
→
−
g ←s ; →
−
p
10
→←−
−
u
−−→ → − → −
unif ( 0 , 1 ) ;
1
− −−− → →
− →
−
→ ← unif ( 0 , 1 ) ;
Al
11 u2
→
−
vp ← ω → −
vp + φ 1 − →·I·− →+φ → − −
→
12 m p u1 2 g · I · u2 ;
→
− −− →
sp ← − −→ min(→ −
sp + →−
vp , →
−
x max ), →
−
T_
13 max( x min );
14
−
→ →
−
si f (mp ) > f ( sp ) alors
15
−
m→←→ −
sp ;
p
16 si f > f (→
∗ −
sp ) alors
17 f ← f (→
∗ −
sp );
18
→
−g ←s ; →
−
p
17
20
e
br
em
pt
se
25
g
Al
T_
9.6 E SSAIMS PARTICULAIRES 211
Exercies
Exercice 9.1. Algorithme génétique pour fonction à une dimension On
doit optimiser une fonction f d’une variable entière x, 0 6 x < 2n . Dans le
cadre d’un algorithme génétique doté d’un opérateur de croisement de deux
solutions, comment coder x sous la forme d’un vecteur binaire ?
Exercice 9.2. Séquences d’inversions Une permutation p des éléments de
1 à n peut être représentée par une séquence d’inversions s, où si compte le
nombre d’éléments de p1 , . . . , pk = i qui sont plus grand que i. Par exemple, la
17
permutation p = (2, 4, 6, 1, 5, 3) a pour séquence d’inversion s = (3, 0, 3, 0, 1, 0).
À quelles permutations correspondent les séquences d’inversions (4, 2, 3, 0, 1, 0)
et (0, 0, 3, 1, 2, 0) ? Donner des conditions nécessaires et suffisantes pour qu’un
20
vecteur s soit une séquence d’inversions d’une permutation. Est-ce que les
opérateurs de croisement standards 1-point, 2-points et uniformes peuvent être
appliqués aux séquences d’inversions ? Comment les séquences d’inversions
peuvent être utilisées dans le contexte de la recherche par dispersion ?
e
Exercice 9.3. Sélection basée sur le rang Quelle est la probabilité de la
br
fonction rank_based_selection(m) donnée dans l’algorithme 9.1 de retourner
une valeur v donnée ?
em
Exercice 9.4. Ajustement de paramètre d’un algorithme génétique Ajuster
la taille de la population et le taux de mutation de la procédure tsp_GA donnée
par le code 9.5 si cette dernière génère en tout 5n enfants.
pt
17
20
e
br
em
pt
se
25
g
Al
T_
Annexe A
17
Codes
20
Cette annexe donne les codes C des quelques procédures utilitaires appa-
raissant dans divers algorithmes présentés dans ce livre.A.1 A.2
e
Listing A.1 – random_generators.c Génération de nombres pseudo-aléatoires.
br
rando() n’a pas de paramètres et retourne un double uniformément réparti entre 0
et 1 et unif retourne un entier entre deux bornes données en paramètres. Génération
em
d’un tableau de n entiers rempli avec une permutation aléatoires des nombres 0 à n − 1.
/ ∗ Random number ] 0 , 1 ] g e n e r a t o r , proposed by Lecuyer , r e p l a c e bad rand ( ) ∗/
double rando ( v o i d )
{ s t a t i c i n t x10 = 12345 , x11 = 67890 , x12 = 13579 , /∗ i n i t i a l value ∗/
x20 = 24680 , x21 = 98765 , x22 = 43210; / ∗ o f seeds ∗ /
pt
c o n s t i n t m = 2147483647; c o n s t i n t m2 = 2145483479;
c o n s t i n t a12= 63308; c o n s t i n t q12 =33921; c o n s t i n t r12 =12979;
c o n s t i n t a13=−183326; c o n s t i n t q13 =11714; c o n s t i n t r13 =2883;
c o n s t i n t a21= 86098; c o n s t i n t q21 =24919; c o n s t i n t r21= 7417;
se
x10 = x11 ; x11 = x12 ; x12 = p12−p13 ; i f ( x12 < 0 ) x12 = x12 + m;
h = x20 / q23 ; p23 = −a23 ∗( x20−h∗q23)−h∗r23 ;
h = x22 / q21 ; p21 = a21 ∗( x22−h∗q21)−h∗r21 ;
i f ( p23 < 0 ) p23 = p23 + m2; i f ( p21 < 0 ) p21 = p21 + m2;
x20 = x21 ; x21 = x22 ; x22 = p21−p23 ; i f ( x22 < 0 ) x22 = x22 + m2;
i f ( x12 < x22 ) h = x12 − x22 + m; e l s e h = x12 − x22 ;
g
i f ( h == 0 ) r e t u r n ( 0 . 5 ) ;
e l s e r e t u r n ( h∗invm ) ;
Al
v o i d swap ( i n t ∗ a , i n t ∗ b )
{ i n t temp = ∗a ; ∗a = ∗b ; ∗b = temp ; }
Listing A.2 – tsp_length.c Calcul de la longueur d’une tournée d’un voyageur de com-
merce. Une solution est donnée sous la forme d’un tableau comportant l’ordre dans le-
quel les villes sont parcourues. Parfois une autre représentation d’une solution, comme
un tableau qui donne les villes qui succèdent, serait meilleure.
double t s p _ l e n g t h ( i n t n , / ∗ Number o f c i t i e s ∗ /
double ∗∗d , /∗ Distance matrix ∗/
const i n t s o l u t i o n [ ] ) / ∗ Order o f t h e c i t i e s ∗ /
{ int i ;
double sum = d [ s o l u t i o n [ n − 1 ] ] [ s o l u t i o n [ 0 ] ] ;
f o r ( i = 0 ; i < n−1; i ++)
sum += d [ s o l u t i o n [ i ] ] [ s o l u t i o n [ i + 1 ] ] ;
r e t u r n sum ;
17
}
20
de commerce. La méthode de résolution ici utilisée est GRASP, qui utilise elle-
même une recherche locale basée sur les listes d’éjections, ainsi que d’autres
fonctions utilitaires. Une valeur de eplsilon permettant d’éviter des problèmes
numériques est de prendre 10−14 fois la plus grande distance séparant 2 villes.
e
Le code A.4 permet quant à lui de tester l’algorithme FANT.
br
Le code A.5 permet quant à lui de tester un algorithme mémétique.
Le code A.6 permet de tester une recherche avec tabous.
Le code A.7 permet de tester une recherche locale multi-objectifs. Ces
em
codes ont été simplifiés de sorte qu’un utilisateur peu habitué à la program-
mation en C, notamment en ce qui concerne la compilation séparée, n’ait qu’à
compiler le programme de test qu’il désire pour obtenir un exécutable fonction-
nel. Leur style de programmation n’est donc pas exemplaire !
pt
se
25
g
Al
T_
215
17
# d e f i n e i n f i n i t e 9.99 e99
double e p s i l o n ;
# i n c l u d e <math . h>
20
# i n c l u d e " random_generators . c " /∗ Listing A.1 ∗/
#include " tsp_length . c " /∗ Listing A.2 ∗/
# i n c l u d e " tsp_LK . c " /∗ Listing 4.4 ∗/
# i n c l u d e " tsp_GRASP . c " /∗ Listing 6.3 ∗/
i n t main ( v o i d )
e
{
int n, i , j ;
int ∗solution ;
double ∗∗ d i s t a n c e ;
double ∗x , ∗y ; br
em
p r i n t f ( " Number o f c i t i e s : \ n " ) ;
s c a n f ( "%d " ,&n ) ;
s o l u t i o n = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
d i s t a n c e = ( double ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ∗ ) ) ;
pt
f o r ( i = 0 ; i < n ; i ++)
d i s t a n c e [ i ] = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
x = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
se
y = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
f o r ( i = 0 ; i < n ; i ++)
{ x [ i ] = rando ( ) ; y [ i ] = rando ( ) ; }
epsilon = 0.0;
25
}
e p s i l o n ∗= 1 . 0 e−14;
Al
tsp_GRASP ( n , d i s t a n c e , s o l u t i o n , 100 , 0 . 8 ) ;
f o r ( i = 0 ; i < n ; i ++)
T_
free ( distance [ i ] ) ;
free ( solution ) ;
r e t u r n EXIT_SUCCESS ;
}
216 C ODES
Listing A.4 – test_tsp_FANT.c Programme de test d’une méthode inspirée des colo-
nies de fourmi artificielles.
/ ∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
Programme f o r t e s t i n g procedures presented i n :
E . T a i l l a r d , I n t r o d u c t i o n aux m é t a h e u r i s t i q u e s , 2015
Compile : gcc −O2 test_tsp_FANT . c
Example o f e x e c u t i o n r e s u l t :
Number o f c i t i e s :
30
Number o f FANT i t e r a t i o n s , FANT parameter :
50 30
FANT 0 1.761522
17
FANT 1 1.660644
FANT 2 1.646723
FANT 3 1.616719
FANT 4 1.614217
20
FANT 28 1.607204
Cost o f s o l u t i o n found w i t h FANT 1.607204e+00
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗ /
# i n c l u d e < s t d i o . h>
# i n c l u d e < s t d l i b . h>
# i n c l u d e <math . h>
e
# d e f i n e i n f i n i t e 9.99 e99
double e p s i l o n ;
# i n c l u d e " random_generators . c "
#include " tsp_length . c "
# i n c l u d e " tsp_LK . c "
#include " i n i t _ u p d a t e _ t r a i l . c "
br /∗
/∗
/∗
/∗
Listing
Listing
Listing
Listing
A.1
A.2
4.4
7.1
∗/
∗/
∗/
∗/
em
#include " generate_solution_trail . c " /∗ Listing 7.2 ∗/
# i n c l u d e " tsp_FANT . c " /∗ Listing 7.3 ∗/
i n t main ( v o i d )
{
pt
i n t n , i , j , f a n t _ i t e r a t i o n s , fant_parameter ;
int ∗solution ;
double ∗∗ d i s t a n c e ;
double l e n g t h ;
se
s o l u t i o n = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
d i s t a n c e = ( double ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ∗ ) ) ;
f o r ( i = 0 ; i < n ; i ++)
d i s t a n c e [ i ] = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
epsilon = 0.0;
g
{ d i s t a n c e [ i ] [ j ] = d i s t a n c e [ j ] [ i ] = rando ( ) ;
i f ( epsilon < distance [ i ] [ j ] )
epsilon = distance [ i ] [ j ] ;
}
T_
e p s i l o n ∗= 1 . 0 e−14;
l e n g t h = tsp_FANT ( n , d i s t a n c e , s o l u t i o n , f a n t _ p a r a m e t e r , f a n t _ i t e r a t i o n s ) ;
p r i n t f ( " Cost o f s o l u t i o n found w i t h FANT %e \ n " , l e n g t h ) ;
f o r ( i = 0 ; i < n ; i ++)
free ( distance [ i ] ) ;
free ( distance ) ;
free ( solution ) ;
r e t u r n EXIT_SUCCESS ;
}
217
Listing A.5 – test_tsp_GA.c Programme de test d’un algorithme mémétique pour pro-
blème de voyageur de commerce dont la matrice des distance est symétrique et géné-
rée aléatoirement.
/ ∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
Programme f o r t e s t i n g procedures presented i n :
E . T a i l l a r d , I n t r o d u c t i o n aux m é t a h e u r i s t i q u e s , 2015
Compile : gcc −O2 test_tsp_GA . c −lm
Example o f e x e c u t i o n r e s u l t :
Number o f c i t i e s :
30
Size o f t h e p o p u l a t i o n , m u t a t i o n r a t e , number o f g e n e r a t i o n s :
17
10 0.01 60
GA i n i t i a l p o p u l a t i o n 11.592710
GA 0 1.707681
GA 2 1.696093
GA 4 1.614217
20
GA 59 1.607204
Cost o f s o l u t i o n found w i t h GA: 1.607204e+00
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗ /
# i n c l u d e < s t d i o . h>
# i n c l u d e < s t d l i b . h>
# i n c l u d e <math . h>
e
# d e f i n e i n f i n i t e 9.99 e99
br
double e p s i l o n ;
# i n c l u d e " random_generators . c " /∗ Listing A.1 ∗/
#include " tsp_length . c " /∗ Listing A.2 ∗/
# i n c l u d e " tsp_LK . c " /∗ Listing 4.4 ∗/
em
# include " rank_based_selection . c " /∗ Listing 9.1 ∗/
# i n c l u d e " OX_crossover . c " /∗ Listing 9.2 ∗/
# i n c l u d e " mutate . c " /∗ Listing 9.3 ∗/
#include " insert_child . c " /∗ Listing 9.4 ∗/
# i n c l u d e " tsp_GA . c " /∗ Listing 9.5 ∗/
pt
i n t main ( v o i d )
{
i n t n , i , j , population_size , nr_generations ;
int ∗solution ;
se
double ∗∗ d i s t a n c e ;
double l e n g t h , m u t a t i o n _ r a t e ;
s o l u t i o n = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
d i s t a n c e = ( double ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ∗ ) ) ;
f o r ( i = 0 ; i < n ; i ++)
g
d i s t a n c e [ i ] = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
epsilon = 0.0;
Al
l e n g t h = tsp_GA ( n , d i s t a n c e , s o l u t i o n , p o p u l a t i o n _ s i z e ,
nr_generations , mutation_rate ) ;
p r i n t f ( " Cost o f s o l u t i o n found w i t h GA: %e \ n " , l e n g t h ) ;
f o r ( i = 0 ; i < n ; i ++)
free ( distance [ i ] ) ;
free ( distance ) ;
free ( solution ) ;
r e t u r n EXIT_SUCCESS ;
}
218 C ODES
Listing A.6 – test_tsp_TS.c Programme de test d’une recherche avec tabous pour
problème de voyageur de commerce dont la matrice des distance est symétrique et
générée aléatoirement.
/ ∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
Programme f o r t e s t i n g procedures presented i n :
E . T a i l l a r d , I n t r o d u c t i o n aux m é t a h e u r i s t i q u e s , 2015
Compile : gcc −O2 t e s t _ t s p _ T S . c
Example o f e x e c u t i o n :
Number o f c i t i e s :
30
Number o f tabu i t e r a t i o n s , min . and max . t a b u _ d u r a t i o n , p e n a l t y :
200 4 20 0.0005
17
TS 1 11.855969
TS 2 10.474485
TS 3 8.993716
...
TS 17 1.828523
20
TS 21 1.698364
TS 26 1.650980
TS 156 1.629700
TS 157 1.610270
TS 158 1.607204
Cost o f s o l u t i o n found w i t h TS 1.607204e+00
e
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗ /
# i n c l u d e < s t d i o . h>
br
# i n c l u d e < s t d l i b . h>
# d e f i n e i n f i n i t e 9.99 e99
double e p s i l o n ;
em
# i n c l u d e " random_generators . c " /∗ Listing A.1 ∗/
#include " tsp_length . c " /∗ Listing A.2 ∗/
# i n c l u d e " perform_move_update . c " /∗ Listing 8.1 ∗/
# i n c l u d e " tsp_TS . c " /∗ Listing 8.2 ∗/
i n t main ( v o i d )
pt
{
i n t n , i , j , i t e r a t i o n s , min_tabu , max_tabu ;
int ∗solution ;
double ∗∗ d i s t a n c e ;
se
double f r e q _ p e n a l t y , l e n g t h ;
s o l u t i o n = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
d i s t a n c e = ( double ∗∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ∗ ) ) ;
f o r ( i = 0 ; i < n ; i ++)
d i s t a n c e [ i ] = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
g
epsilon = 0.0;
f o r ( i = 0 ; i < n−1; i ++)
Al
epsilon = distance [ i ] [ j ] ;
}
e p s i l o n ∗= 1 . 0 e−14;
generate_random_permutation ( n , s o l u t i o n ) ;
length = tsp_length ( n , distance , s o l u t i o n ) ;
tsp_TS ( n , d i s t a n c e , s o l u t i o n , &l e n g t h ,
i t e r a t i o n s , min_tabu , max_tabu , f r e q _ p e n a l t y ) ;
p r i n t f ( " Cost o f s o l u t i o n found w i t h TS %e \ n " , l e n g t h ) ;
f o r ( i = 0 ; i < n ; i ++)
free ( distance [ i ] ) ;
free ( distance ) ;
free ( solution ) ;
r e t u r n EXIT_SUCCESS ;
}
219
17
# i n c l u d e < s t r i n g . h>
# d e f i n e i n f i n i t e 9.99 e99
#define K 3 / ∗ Number o f o b j e c t i v e s ∗ /
double e p s i l o n ;
20
# i n c l u d e " random_generators . c " / ∗ L i s t i n g A.1 ∗ /
# i n c l u d e " KD_tree_add_scan . c " / ∗ L i s t i n g 4.6 ∗ /
# i n c l u d e " KD_tree_delete . c " / ∗ L i s t i n g 4.7 ∗ /
# include " tsp_3opt_pareto . c " / ∗ L i s t i n g 4.5 ∗ /
# i n c l u d e " KD_tree_update_pareto . c " / ∗ L i s t i n g 4.8 ∗ /
e
i n t main ( v o i d )
{
int n, i , j , k;
i n t ∗succ ;
double ∗∗∗ d i s t a n c e ;
POINT l e n g t h s ; br
em
NODE ∗ p a r e t o _ f r o n t = NULL ;
succ = ( i n t ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( i n t ) ) ;
pt
f o r ( i = 0 ; i < n ; i ++)
d i s t a n c e [ k ] [ i ] = ( double ∗) m a l l o c ( ( s i z e _ t ) n ∗ s i z e o f ( double ) ) ;
}
f o r ( j = 0 ; j < n ; j ++)
d i s t a n c e [ k ] [ i ] [ j ] = rando ( ) ;
e p s i l o n = 1 . 0 e−13;
lengths [ k ] = 0;
f o r ( i = 0 ; i < n ; i ++)
l e n g t h s [ k ] += d i s t a n c e [ k ] [ i ] [ succ [ i ] ] ;
}
T_
t s p _ 3 o p t _ p a r e t o ( n , d i s t a n c e , succ , l e n g t h s , & p a r e t o _ f r o n t ) ;
KD_scan_and_delete ( p a r e t o _ f r o n t , n ) ;
f o r ( k = 0 ; k < K ; k ++)
{
f o r ( i = 0 ; i < n ; i ++)
free ( distance [ k ] [ i ] ) ;
free ( distance [ k ] ) ;
}
free ( distance ) ;
f r e e ( succ ) ;
r e t u r n EXIT_SUCCESS ;
}
T_ 220
Al
g
25
se
pt
em
br
e
20
17
C ODES