Modèles de Recherche Opérationnelle

Fabian Bastin
bastin@[Link]
Département d’Informatique et de Recherche Opérationnelle
Université de Montréal

[Link]
La présente version de ce document s’inspire des notes de Patrice Marcotte, Bernard Gendron, ainsi que des livres
Introduction to Operational Research [2] et The Basics of Practical Optimization [5].

Photographie de couverture: viaduc de Millau, France.

c Fabian Bastin, 2006

Le présent document peut être modifié et redistribué à des fins non commerciales, sous conditions d’être diffusé
sous les même conditions.
c Fabian Bastin, 2015
 Table des matières

1 Introduction 1
1.1 Les origines de la recherche opérationnelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 La nature de la recherche opérationnelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Modélisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 Algorithmes et logiciels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4.1 Un exemple avec GAMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Programmation linéaire 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Modèle général de programmation linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Terminologie de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Hypothèses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 La méthode du simplexe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5.1 Solution de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5.2 Interprétations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5.3 Critère d’optimalité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.5.4 Adaptation à d’autres formes de modèles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.5.5 Obtention d’une base admissible initiale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5.6 Variables à valeurs quelconques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.6 Dualité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.6.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.6.2 Analyse de sensibilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.7 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3 Programmation non linéaire 31

3.1 Fonctions convexes et concaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.1 Ensembles convexes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 Algorithmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.1 L’algorithme du simplexe dans le cas non-linéaire . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.2 Optimisation sans contrainte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.3 Méthode de la bissection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.4 Méthode du gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.5 Optimisation sous contraintes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2.6 Conditions d’optimalité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4 Programmation mixte entière 43

4.1 Présentation générale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Contraintes mutuellement exclusives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.1 Deux contraintes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.2 K contraintes parmi N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.3 Fonction ayant N valeurs possibles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2.4 Objectif avec coûts fixes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

iii
iv TABLE DES MATIÈRES

4.2.5 Variables entières en variables 0–1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.2.6 Problème de recouvrement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3 Stratégies de résolutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3.1 Relaxation linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3.2 Approche par énumération . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3.3 Les hyperplans coupants (méthode de coupe) . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.4 Modélisation avec GAMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.6 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5 Réseaux 67
5.1 Graphes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.1.1 Graphe orienté . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.1.2 Graphe non orienté . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1.3 Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1.4 Chemins et circuits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1.5 Connexité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.2 Problème du chemin le plus court - algorithmes de corrections d’étiquettes . . . . . . . . . . . . . . . 69
5.2.1 Algorithme de Dijkstra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.3 Flot dans un réseau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.3.1 Réseau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.3.2 Modèle de flot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.3.3 Algorithme de Bellman-Ford . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.3.4 Modèle du chemin critique (PERT/CPM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.4 Problème de l’arbre partiel minimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.4.1 Algorithme de Prim (1957) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.5 Problème du flot maximum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.5.1 Algorithme de Ford-Fulkerson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.5.2 Flot maximum - Coupe minimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.6 Problème de flot minimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.7 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.8 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

6 Modèles stochastiques 83
6.1 Concepts de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.2 Variable aléatoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.2.1 Variables aléatoires discrètes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2.2 Variables aléatoires continues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2.3 Espérance mathématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2.4 Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.3 Loi de probabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.3.1 Loi de Bernouilli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.3.2 Loi uniforme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.3.3 Loi de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3.4 Loi exponentielle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4 Modèles stochastiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4.1 Processus stochastiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4.2 Chaînes de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.5 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
TABLE DES MATIÈRES v

7 Programmation dynamique 91
7.1 Principe d’optimalité de Bellman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.2 Affectation de ressources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
7.3 Programmation dynamique déterministe et plus court chemin . . . . . . . . . . . . . . . . . . . . . . 96
7.4 Cas probabiliste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
7.5 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

8 Simulation 99
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.2 Files d’attente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.2.1 Concepts de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8.2.2 Modèle M/M/1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
8.3 Simulation à événements discrets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
8.4 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Annexe 105
8.5 Logiciels d’optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.5.1 IOR Tutorial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.5.2 GAMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.5.3 Autres logiciels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
8.6 Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Bibliographie 111
vi TABLE DES MATIÈRES
OR is the art of winning wars without actually fighting.
Aurther Clarke

1
Introduction

1.1 Les origines de la recherche opérationnelle

Les origines de la recherche opérationnelle, en abrégé RO, ne sont pas clairement définies,
bien que Charles Babbage (1791–1871) soit souvent considéré comme le père de la recherche
opérationnelle, en raison de ses travaux sur le coût du transport et le tri du courrier, menant à
l’introduction du Penny Post en Angleterre, en 1840.
Un siècle plus tard, la seconde guerre mondiale, de part son envergure, créa une besoin urgent
d’allouer de manière efficace des ressources limitées aux différentes opérations militaires et aux
activités au sein de chaque opération. En particulier, l’organisation militaire britannique, puis
américaine, mit à contribution un grand nombre de scientifiques pour gérer ces allocations, et
s’occuper d’autres problèmes stratégiques et tactiques. Ce faisant, ils furent appelés à poursuivre
des recherches sur des opérations (militaires), et constituèrent les premières équipes de RO,
notamment celle dirigée par Patrick Blackett, prix Nobel de physique en 1948 [3]. Leurs efforts
furent significatifs dans la marche vers la victoire, par exemple en ce qui touche l’utilisation
du radar, nouvellement développé. Ces succès encouragèrent la poursuite de l’utilisation de la
Figure 1.1 – RO dans d’autres domaines. La croissance importante de l’industrie d’après-guerre entraîna
Patrick Bla- des problèmes, causés par la complexité croissante et la spécialisation dans les organisations,
ckett problèmes en fait proches de ceux présent lors du conflit. Au début des années 1950’s, la RO
avait pénétré une multitude d’organisations commerciales, industrielles, et gouvernementales.
Et ce n’était que le début.
Au moins deux autres facteurs ont joué un rôle clé dans la croissance rapide de la RO. Tout d’abord, des
progrès substantiels ont été obtenus très tôt afin d’améliorer les techniques de RO. Ces techniques, dans leur mise
en pratique, furent soutenues par l’essor des outils informatiques.

1.2 La nature de la recherche opérationnelle

“Rechercher sur des opérations” touche tous les problèmes reliés à la conduite et à la coordination des opérations
(activítés) au sein d’une organisation. Cette organisation peut représenter des domaines très divers: l’industrie
manufacturière, le transport, la construction, les télécommuncations, la finance, les soins de santé,. . . . La RO,
associée à la révolution informatique, pénètre pratiquement tous les secteurs d’activités de la vie courante, même
si sa présence est souvent invisible.
La première étape de la “recherche” est l’observation attentive du problème et sa formulation, ainsi que la collecte
de données associées. Il convient par la suite de construire un modéle scientifique qui tente l’abstraire l’essence du

1
2 CHAPITRE 1. INTRODUCTION

problème réel. Tout modèle est une simplification de la réalité, mais cette représentation doit être suffisamment
précise pour capturer les caractéristiques essentielles de la situation, et de pouvoir tirer des conclusions valides pour
le problème réelle. Il conviendra dès lors de tester ce modèle, et de le modifier au besoin.
Une caractéristique additionnelle est que la RO essaye souvent de trouver la meilleure solution (dite solution
optimale) pour le problème examiné. Cette solution peut ne pas être unique. Cette recherche d’optimalité est un
thème important en RO, même si son interprétation en terme managériels peut être délicate. Il est difficile pour
un individu de pouvoir maîtriser tous les aspects du problèmes à l’étude, de sorte que la RO est généralement
plus un travail d’équipe, avec des experts en mathématiques, statistiques et probabilités, ingénierie, économie,
administration, informatique, physiques, sciences comportementales, et les techniques spécifiques de la RO.

1.3 Modélisation
Un modèle, telle que considéré dans ce cours, est une construction mathématique utilisée pour représenter
certains aspects significatifs de problèmes du monde réel. Il y a beaucoup de types différents de modèles mathéma-
tiques, mais nous nous focaliserons dans un premier temps sur les modèles d’optimisation. Il y a trois composantes
principales dans un modèle d’optimisation:
Variables: elles représentent les composantes du modèle qui peuvent être modifiées pour créer des configurations
différentes.
Contraintes: elles représentent les limitations sur les variables.
Fonction objectif : cette fonction assigne une valeur à chaque configuration différente. Le terme “objectif” vient
du fait que l’objectif est d’optimiser cette fonction.
Exemple 1 (Un exemples de décisions binaires (oui/non)). Un étudiant en quête d’une université projette de
visiter les campus de trois universités du Maine au cours d’un voyage unique, débutant et finissant à l’aéroport de
Portland. Les trois établissements sont dans les villes de Brunswick, Lewiston, et Waterville, et l’étudiant ne veut
visiter chaque ville qu’une seule fois, tout en maintenant le trajet total le plus court possible. Les distances entre ces
villes sont données dans la Table 1.1.

Ville Portland Brunswick Lewiston Waterville

Portland 0 26 34 78
Brunswick 26 0 18 52
Lewiston 34 18 0 51
Waterville 78 52 51 0

Table 1.1 – Distances entre les villes (miles)

L’étape la plus importante dans la construction d’un modèle est le choix des variables qui vont entrer en jeu.
Dans le présent cas, puisque n’importe quel trajet consiste en une série de petits déplacements entre deux villes, il est
raisonnable d’assigner des variables aux décisions de partir ou non d’une ville vers une autre. Pour plus de faciliter,
numérotons les villes comme suit: 1 pour Portland, 2 pour Brunswick, 3 pour Lewiston et 4 pour Waterville. Ainsi,
nous aurons une variable x1,2 égale à 1 si l’étudiant voyage de Portland à Brunswick au cours de son parcours total,
et 0 sinon. Puisqu’il n’y a pas de voyage d’une ville vers cette même ville, nous avons d’ores et déjà les contraintes

xi,i = 0, i = 1, . . . , 4. (1.3.1)

Une fois les variables choisies, nous pouvons essayer de formuler le problème. Ce processus est en fait souvent une
manière utile pour guider le choix des variables.
Chaque ville ne devant être visitée qu’une seule fois, elle ne peut apparaître qu’une seule fois comme ville
d’arrivée. En d’autres termes, pour j fixé, xi,j ne peut être non-nul que pour un i donné, avec i 6= j. Une manière
plus simple d’encoder cette information est d’écrire, pour j = 1, . . . , 4,

x1,j + x2,j + x3,j + x4,j = 1,

1.3. MODÉLISATION 3

ou de manière plus concise.

4
X
xi,j = 1, j = 1, . . . , 4. (1.3.2)
i=1

Les contraintes formulées jusqu’à présent ne garantissent aucune forme de trajet ayant même départ et arrivée.
Par exemple, l’affectation x1,2 = 1, x1,3 = 1, x1,4 = 1, x2,1 = 1, et toutes les autres variables égales à 0, satisfont
les contraintes 1.3.1 et 1.3.2. Cette solution décrit toutefois un schéma de visites impossible puisque Portland est
l’origine de tous les déplacements aux trois autres villes universitaires, mais n’est destination que depuis Brunswick.
Nous avons évidemment aussi besoin des contraintes
4
X
xi,j = 1, i = 1, . . . , 4, (1.3.3)
j=1

afin d’assurer que chaque ville ne serve d’origine que pour exactement un déplacement vers une autre ville. Fi-
nalement, afin d’obtenir un véritable trajet ayant même origine et départ, nous devons rejeter les affectations qui
décrivent des groupes déconnectés de petits déplacements comme x1,2 = x2,1 = 1, x3,4 = x4,3 = 1, avec toutes les
autres variables égales à 0. Nous pouvons forcer ceci avec les contraintes

xi,j + xj,i ≤ 1, i = 1, . . . , 4, et j = 1, . . . , 4.

Cette contrainte exclut tout mini-cycle.

Les contraintes définies, nous devons décrire la distance totale associée à n’importe quel parcours autorisé.
Puisque nos variables ont seulement comme valeurs possibles 0 ou 1, nous pouvons multiplier chacune d’elle par la
distance correspondante entre les deux villes indexées, et les additionner:
4 X
X 4
xi,j ai,j .
i=1 j=1

Notre modèle mathématique consiste à minimiser cette fonction, dite fonction objectif par rapport aux variables
xi,j , tout en satisfaisant les contraintes préalablement décrites:
4 X
X 4
min xi,j ai,j ,
x
i=1 j=1

s.c. xi,i = 0, i = 1, . . . , 4,
4
X
xi,j = 1, j = 1, . . . , 4,
i=1
4
X
xi,j = 1, i = 1, . . . , 4,
j=1

xi,j + xj,i ≤ 1, i = 1, . . . , 4, et j = 1, . . . , 4,
xi,j ∈ {0, 1}, i = 1, . . . , 4, et j = 1, . . . , 4.

Ici, x = (xi,j )i=1,...,4,j=1,...,4 , et s.c., “sous les contraintes”. Le problème d’optimisation ainsi construit constitue un
programme mathématique.
Le problème de visites d’universités est assez petit que pour être résolu explicitement, sans recourir à des mé-
thodes d’optimisation numérique. Puisqu’il y a seulement trois parcours significativement différents, la distance totale
associée à chacun d’eux pourrait être facilement calculée, et nous choisissons le parcours de longueur minimale, qui
est ici
Portland → Brunswick → Waterville → Lewiston → Portland.
avec une distance totale de 163 miles. Il est cependant clair qu’une telle stratégie de résolution ne fonctionne plus
comme le nombre de villes augmente.
4 CHAPITRE 1. INTRODUCTION

Exemple 2 (Un problème de mélange). Un armateur doit construire un navire de guerre à partir de 50 tonnes
d’acier contenant entre 0.5% et 1.25% de carbone (C), entre 0.3% and 0.5% de silicone (Si), pas plus de 0.05% de
sulfure (Su), et pas plus de 0.04% de phosphore (Ph). Un fournisseur produit de l’acier à partir de sept matières
premières dont les qualités, les disponibilités en tonnes, et les coûts en $/tonne sont donnés dans la Table 2. Le
fournisseur veut déterminer la combinaison la moins coûteuse de composants bruts qu’il peut utiliser pour produire
l’acier répondant aux besoins de l’armateur.

Matière première %C % Si % Su % Ph Disponibilité Coût

limonite 3.0 0 0.013 0.015 40 200
taconite 2.5 0 0.008 0.001 30 250
hématite 0 0 0.011 0.05 60 150
magnétite 1.2 0 0.002 0.008 50 220
silicone 1 0 90 0.004 0.002 20 300
silicone 2 0 96 0.012 0.003 30 310
charbon 90 0 0.002 0.01 25 165

Table 1.2 – Données pour le problème de production d’acier

Puisque les fournisseur peut changer les quantités de matéries premières utilisées dans la producton de l’acier,
nous pourrions assigner une variable différente pour représenter la quantiter de chaque matière première:
— x1 = tonnes de limonite,
— x2 = tonnes de taconite,
— x3 = tonnes d’hématite,
— x4 = tonnes de magnétite,
— x5 = tonnes de silicone 1,
— x6 = tonnes de silicone 2,
— x7 = tonnes de charbon.
Notons que les variables sont ici continues, contrairement à l’exemple précédent.
Afin de modéliser les contraintes, observons tout d’abord que les variables dans ce cas sont naturellement bornées
inférieurement par 0 (puisque des quantités négatives ne feraient pas de sens), et bornées supérieurement par leur
quantité disponible, aussi avons-nous:

0 ≤ x1 ≤ 40,
0 ≤ x2 ≤ 30,
0 ≤ x3 ≤ 60,
0 ≤ x4 ≤ 50,
0 ≤ x5 ≤ 20,
0 ≤ x6 ≤ 30,
0 ≤ x7 ≤ 25.

En supposant que n’importe quelle quantité d’une matière première contribue pour la même quantité d’acier, et en
sachant que nous devons produire au moins 50 tonnes, nous avons
7
X
xi ≥ 50.
i=1

Notons que nous ne supposons pas que nous produirons exactement 50 tonnes, puisqu’il peut être nécessaire de
produire d’avantage afin de satisfaire les autres exigences du problème.
L’autre caractéristique contraignante dans ce problème que que l’acier doit contenir un certain pourcentage de
carbone, de silicone, de sulfure et de phosphore. Afin de voir comment ces exigences de composition se traduisent en
contraintes par rapport à nos variables, nous nous concentrerons d’abord sur l’exigence d’avoir entre 0.5% et 1.25%
1.4. ALGORITHMES ET LOGICIELS 5

de carbone, en espérant que les exigences sur le silicone, le sulfure et le phosphore se formulent de manière similaire.
À partir des données, nous connaissons le pourcentage de contribution en carbone de chaque matière première, aussi
nous pouvons facilement calculer la quantité de carbone pour n’importe quel choix de variables comme
0.03x1 + 0.025x2 + 0.012x4 + 0.9x7 .
Cependant, comme nous avons une exigences de proportion de carbone dans l’acier, nous devons diviser cette
quantité de carbone par la quantité d’acier:
 
tonnes de carbone 3.0x1 + 2.5x2 + 1.2x4 + 90x7
%C = 100 = .
tonnes d’acier x1 + x2 + x3 + x4 + x5 + x6 + x7
La contrainte que l’acier contienne entre 0.5% et 1.25% de carbone se traduit dans la paire de contraintes
3.0x1 + 2.5x2 + 1.2x4 + 90x7
0.5 ≤ ≤ 1.25. (1.3.4)
x1 + x2 + x3 + x4 + x5 + x6 + x7
Les contraintes pour les autres composants se formulent de manière similaires.
Puisque ce problème implique de trouver la combinaison la moins coûteuse de matières premières qui rencontre
la demande de 50 tonnes d’acier, la fonction objectif est simplement le coût des matières premières utilisées:
coût = 200x1 + 250x2 + 150x3 + 220x4 + 300x5 + 310x6 + 165x7 ,
où chaque matière première contribue pour son propre coût au total. Le problème d’optimisation est dès lors la
minimisation de cette fonction coût sur tous les choix des variables qui satisfont les contraintes modélisées.
Il n’est plus possible ici d’énumérer les solutions possibles afin de résoudre le modèle, en particuler puisque les
variables considérées sont continues. Nous pouvons néanmoins obtenir une intuition de la solution en considérant le
comportement de la fonction objectif. Ainsi, il est évident que celle-ci décroît quand une des variables diminue, et que
la contribution la plus faible en termes de coût vient des variables avec les plus petits coefficients (i.e. les matières
premières avec les coûts moindres par unité). Si nous ignorons les contraintes de composition, le fournisseur devrait
produire exactement 50 tonnes d’acier à partir des matières premières disponibles les moins chères. Ceci signifierait
utiliser 50 des 60 tonnes disponibles d’hématite (au coût de $150 par tonne), pour un coût total des $7,500.
Avant d’essayer de résoudre le problème d’optimisation complet (à l’aide d’un ordinateur), nous devrions ré-
écrire les contraintes de composition dans une forme plus simple. Par exemple, la contrainte 1.3.4 est nonlinéaire,
mais peut être réexprimée comme deux contraintes linéaires en multipliant chaque terme des inégalités par le dé-
nominateur. Après la simplification de toutes les contraintes de composition de cette manière, et en utilisant un
logiciel d’optimisation, nous obtenons la solution
x1 = 13.7888, x3 = 35.8097, x5 = 0.166667, x7 = 0.234818, x2 = x4 = x6 = 0,
qui se traduit en exactement 50 tonnes d’acier, produit à partir de 13.7888 tonnes de limonite, 35.8097 tonnes
d’hématite, 0.166667 tonnes de silicone 1 et 0.234818 tonnes de charbon. Le coût total de production pour cette
solution est $8,217.96.

1.4 Algorithmes et logiciels

Au-delà de la modélisation, la résolution de problèmes de recherche opérationnelle nécessite de recourir à
des algorithmes adapatés à la nature du problème, et capables de traiter de quelques dizaines à des millions de
variables. Leur étude consistera par conséquent une partie importante du présent document. Même si de nombreux
logiciels mettant en oeuvre ces algorithmes sont commerciaux, le monde de l’open-source n’est pas en reste avec
notamment des projets tels que COIN-OR ([Link] Il existe aussi diverses versions d’évaluations
de solveurs commerciaux, ainsi que l’interface web NEOS ([Link] L’utilisation de tels
outils nécessitent toutefois l’apprentissage de langages de modélisation adaptés ; dans ce cours, nous nous baserons
sur le langage GAMS ([Link] Ces langages servent à décrire dans des termes compréhensibles par
le solveur le problème à résoudre. La dernière version de IOR-Tutorial, développé en Java et proposé en complément
à Hillier et Lieberman [2], peut également se révéler un complément précieux pour se familiariser avec les techniques
de recherche opérationnelle.
6 CHAPITRE 1. INTRODUCTION

1.4.1 Un exemple avec GAMS

GAMS est l’acronyme de General Algebric Modeling System. Il consiste en un langage de description de pro-
blèmes, qui peut être compilé, et en une interface vers différent solveurs. Une version gratuite de démonstra-
tion, permettant de résoudre des problèmes de petite taille, est disponible en téléchargement à l’adress http:
//[Link]. Ce langage peut également être utilisé avec divers solveurs proposés sur NEOS. Nous reviendrons
plus en détail ultérieurement sur la formulation d’un programme mathématique avec GAMS, et nous contenterons
d’un petit exemple introductif à ce stade.
Considérons un brasseur qui disposent des éléments suivants dans son stock:
1. malt (75 kg) ;
2. houblon (60 kg) ;
3. levure (50 kg).
Deux produits peuvent être obtenus: de la bière légère et de la bière noire. Pour un kg de bière légère, il faut 2 kg
de malt, 3 kg de houblon, et 2 kg de levure. Pour un kg de bière noire, il faut 3 kg de malt, 1 kg de houblon, et 5/3
kg de levure. La bière noire est vendue la moitié du prix de la bière légère. En supposant que toute la production
sera vendue, le brasseur souhaiterait décider des quantités de bières noire et légère à produire pour optimiser son
profit. Le programme GAMS résultant est

set b types of beer /light, dark/;

set i inputs /malt, hops, yeast/

parameter r(i) raw supplies /malt 75, hops 60, yeast 50/;

table a(i,b) input requirements

light dark
malt 2 3
hops 3 1
yeast 2 1.67

parameter p(b) selling price / light 2, dark 1/;

variables pi profit (maximand)

x(b) production level;

equations profit defines gross revenue

supply(i) input supply constraint;

profit.. pi =e= sum(b, p(b) * x(b));

supply(i).. sum(b, a(i,b)*x(b)) =L= r(i);

model beer /all/;

[Link](b) = 0;

solve beer using lp maximizing pi;

1.5 Notes
Le contenu des Section 1.1 et 1.2 se base principalement sur Hillier et Lieberman [2], Chapitre 1. La Section 1.3
est construite à partir du Chapitre 1 de Levy [5]. L’exemple GAMS est dû à Thomas F. Rutherford.
“All right,” said Deep Thought. “The Answer to the Great Question. . . ”
The Hitchhiker’s Guide To the Galaxy

2
Programmation linéaire

2.1 Introduction
Bien que la réalité soit souvent loin d’être linéaire, un grand nombre de problèmes peuvent s’écrire sous forme
linéaire, soit directement, soit en première simplification. D’autre part, un très grand nombre de modèles constituent
des extensions de programmes linéaires. Sa compréhension est essentielle à la compréhension de modèles plus
sophistiqués.
Un programme linéaire générique s’écrit sous la forme
maxx c1 x1 + ... + cn x n
t.q. a11 x1 + ... + a1n xn ≤ b1
.. ..
. .
am1 x1 + ... + amn xn ≤ bm ,
ou, sous une forme plus compacte,
n
X
max cj xj
x
j=1
Xn
t.q. aij xj ≤ bi , i = 1, . . . , m.
j=1

La ligne
n
X
cj x j
j=1

représente la fonction objectif, que nous souhaitons maximiser. La maximisation se fait en respectant les m
contraintes
Xn
aij xj ≤ bi , i = 1, . . . , m.
j=1

Sous forme matricielle, le problème se réecrit

max cT x
x
t.q. Ax ≤ b,

7
8 CHAPITRE 2. PROGRAMMATION LINÉAIRE

avec
     
x1 c1 b1
x =  ...  ,
 ..   .. 
c =  . , b=.
 

xn cn bn
 
a11 a12 ... a1n
A =  a21 a22 ... a2n  .
am1 am2 ... amn

La terminologie “linéaire” vient du fait que toutes les fonctions impliquées sont linéaires. Typiquement, nous
ajouterons également des contraintes de non-négativités:

xi ≥ 0, i = 1, . . . , n,

ou, en abrégé,
x ≥ 0.
Nous dirons aussi que le vecteur x appartient à l’orthant positif (i.e. x ∈ Rn+ ).

Exemple 3 (Wyndor Glass). La compagnie Wyndor Glass Co. produit des produits verriers de haute qualité,
incluant des fenêtres et des portes vitrées. Elle dispose à cette fin de trois usines (usine 1, usine 2, usine 3), qui ont
chacune une capacité de production limitée. Les châssis en aluminium et les matériaux sont produits dans l’usine 1,
les châssis en bois sont fabriqués dans l’usine 2, et l’usine 3 produit le verre et assemble les produits. La compagnie
a décidé de mettre en place de ligne de production:
— produit 1: une porte vitrée avec un châssis d’aluminium ;
— produit 2: une fenêtre double-vitrage avec châssis en bois.
Un lot de 20 unités donne lieu à un profit de $3000 et $5000, respectivement pour le produit 1 et le produit 2.
Les données du problème sont synthétisées dans la Table 2.1. Chaque lot d’un produit est le résultat combiné de la

Produit 1 Produit 2 Capacité de

Temps de prod. (h) Temps de prod. (h) production (h)
Usine 1 1 0 4
Usine 2 0 2 12
Usine 3 3 2 18

Table 2.1 – Données du problème Wyndor Glass

production dans les trois usines. Nous souhaitons déterminer le taux de production pour chaque produit (nombre de
lots par semaine) de façon à maximiser le profit total.
Les variables de décision sont
— x1 , le nombre de lots du produit 1 ;
— x2 , le nombre de lots du produit 2.
La fonction objectif est le profit total, qui vaut 3x1 +5x2 , en l’exprimant en miller de dollars. Nous voulons maximiser
ce profit.
Les contraintes concernent tout d’abord les capacítés de production:

x1 ≤4 (usine 1)
2x2 ≤ 12 (usine 2)
3x1 + 2x2 ≤ 18 (usine 3)

Viennent ensuite les contraintes de non-négativité:

x1 ≥ 0, x2 ≥ 0. (nombre positif d’unités produites)

2.1. INTRODUCTION 9

En résumé, nous avons le problème d’optimisation suivant:

max z = 3x1 + 5x2

x

sous les contraintes

x1 ≤4 (usine 1)
2x2 ≤ 12 (usine 3)
3x1 + 2x2 ≤ 18 (usine 3)

x1 ≥ 0, x2 ≥ 0 (non-négativité)

Avant de traiter de méthodes numériques, essayons de visualiser le problème afin de développer une intuition
quant à sa résolution. Le problème peut-être représenté comme sur la Figure 2.1. Pour le réaliser, nous traçons
d’abord les droites correpondant aux contraintes, puis nous déterminons le domaine réalisable en vérifiant le sens
des inégalités pour chacune d’elle. Nous traçons ensuite les droites correspondant à la variation de l’objectif. Dans
l’exemple
3 1
z = 3x1 + 5x2 ⇔ x2 = − x1 + z.
5 5

L’ordonnée à l’origine, dépendant de la valeur de z, est 15 z, et la pente vaut − 35 . Maximiser revient à augmenter z.

x2

9 3x1 + 2x2 ≤ 18
8 x1 ≤ 4

7
6 2x2 ≤ 12

5
4
domaine
3
réalisable
2
1
0 x1
0 1 2 3 4 5 6 7 8 9

Figure 2.1 – exemple Wyndor Glass

La représentation graphique, bien qu’intéressante pour “voir” comment se passent les choses, ne fonctionne plus
dès que nous avons plus de deux variables. Il faut alors passer à l’utilisation de logiciels, utilisant l’algorithme du
simplexe ou de points intérieurs. Nous n’aborderons que le premier dans le cadre de ce cours introductif.
10 CHAPITRE 2. PROGRAMMATION LINÉAIRE

x2

36 = 3x1 + 5x2 8
7
(2, 6)
6
5
20 = 3x1 + 5x2
4
3
10 = 3x1 + 5x2
2
1
0 x1
0 1 2 3 4 5 6 7 8 9

Figure 2.2 – exemple Wyndor Glass: maximisation

2.2 Modèle général de programmation linéaire

Le modèle complet

max cT x
x
Ax ≤ b,
x ≥ 0,

est appelé forme standard. D’autres formes sont possibles et définissent aussi des modèles de programmation linéaire:
— minimiser au lieu de maximiser:
min f (x) = − max f (x);

— les inégalités peuvent être remplacées par des égalités, ou être de sens contraire ;
— certaines variables peuvent ne pas être forcés à être supérieures à 0. Par exemple, considérons la contrainte

x ≥ −4,

qui équivaut à
x + 4 ≥ 0.

Nous pouvons alors définir y = x + 4, de sorte que la contrainte devient

y ≥ 0.

De même, considérons
−10 ≤ x ≤ −2.

En définissant y = x + 10, nous avons les contraintes y ≥ 0, y ≤ 8.

2.3. TERMINOLOGIE DE BASE 11

2.3 Terminologie de base

Nous parlerons de solution réalisable à propos d’une solution (i.e. une instance particulière du vecteur x) pour
laquelle toutes les contraintes sont satisfaites. En d’autres termes cette solution appartient au domaine réalisable.
Par exemple, le point (1, 1) est une solution réalisable pour la Figure 2.1. A contrario, une solution non réalisable
est une solution pour laquelle au moins une contrainte est violée. Elle n’appartient pas au domaine réalisable. Une
solution optimale est une solution donnant la meilleure valeur possible pour l’objectif. Cette valeur est dite valeur
optimale. Un modèle n’a aucune solution optimale si son domaine réalisable est vide, ou si la fonction objectif est
non bornée. Ajoutons par exemple la contrainte

3x1 + 5x2 ≥ 50

dans l’exemple Wyndor Glass. Cette nouvelle contrainte ne peut être satisfaite en même temps que les contraintes
précédentes, comme illustré dans la Figure 2.3. Plus aucun point n’est réalisable. Si nous supprimons plutôt les
contraintes
2x2 ≥ 12, 3x1 + 2x2 ≤ 18,
la fonction objectif devient non bornée dans le domaine réalisable, comme illustré sur la Figure 2.4. Un modèle peut
également présenter une infinité de solutions optimales. Considérons dans l’exemple la fonction objectif

z = 3x1 + 2x2 .

Tout point sur le segment [(2, 6), (4, 3)] est alors solution optimale (voir Figure 2.5).

x2

10
9
8
7
6
5
4
3
2
1
0 x1
0 1 2 3 4 5 6 7 8 9

Figure 2.3 – exemple Wyndor Glass: domaine non réalisable

2.4 Hypothèses
La première hypothèse d’un modèle de programmation linéaire est la proportionnalité:
12 CHAPITRE 2. PROGRAMMATION LINÉAIRE

x2

9 x1 ≤ 4

8
7
6
5
4
domaine
3
réalisable
2
1
0 x1
0 1 2 3 4 5 6 7 8 9

Figure 2.4 – exemple Wyndor Glass: objectif non borné

18 = 3x1 + 3x2
x2

9 3x1 + 2x2 ≤ 18
8 x1 ≤ 4

7
sol

6 2x2 ≤ 12
ut
ion

5
so
pt

4
im
ale

domaine
3
s

réalisable
2
1
0 x1
0 1 2 3 4 5 6 7 8 9

Figure 2.5 – exemple Wyndor Glass: infinité de solutions

2.4. HYPOTHÈSES 13

— la contribution de chaque variable à la valeur de la fonction objectif est proportionnelle à la valeur de cette
variable.
— la contribution de chaque variable au terme de gauche de chaque contrainte fonctionnelle est proportionnelle
à la valeur de cette variable.
Des cas où cette hypothèse n’est pas satisfaite sont
— coût fixe initial (Figure 2.6).
— profit marginal (profit par unité) croissant ou décroissant (Figure 2.7).

x2

4
l
3 n ne
tio
2 or

l
ne
op
pr

n
io
1

t
or
op
pr
0 x1

n
0 1 2 3 4
−1 no
coût fixe initial

Figure 2.6 – coût fixe initial

x2

7
6
5
4
nt
sta

3
on

t
nc

2 tan
s
no

n
1 co

0 x1
0 1 2 3 4

Figure 2.7 – profit marginal

Nous avons d’autre part la propriété d’additivité: la fonction objectif est composée de la somme des contributions
individuelles de chaque variable, et le terme de gauche de chaque contrainte fonctionnelle est composé de la somme
des contributions individuelles de chaque variable. L’additivité interdit les termes de la forme x1 x2 . Si une de ces
hypothèses est invalidée, nous sommes en présence d’un programme non linéaire, étudié au Chapitre 3.
D’autre part, les variables sont continues. Si nous imposons des variables à valeurs entières, nous obtenons
un modèle de programmation en nombres entiers (Chapitre 4). De plus, les valeurs affectées à chaque paramètre
14 CHAPITRE 2. PROGRAMMATION LINÉAIRE

sont des constantes connues avec certitude. Cette hypothèse peut être fort éloignée de la réalité. Nous pouvons
tester cette hypothèse en conduisant une analyse de sensibilité, qui consiste à vérifier la sensibilité du modèle à des
changements de valeurs des paramètres. Nous pouvons également introduire des variables aléatoires, ce qui conduit
typiquement à un problème de programmation stochastique. La programmation stochastique dépasse cependant les
objectifs du présent cours.
Exemple 4 (Horaire de personnel). Nous souhaitons établir un horaire quotidien, sachant que chaque jour est
divisé en périodes et en supposant que nous avons pu estimer un nombre minimum d’employés devant être affecté
durant chaque période. Chaque jour est divisé en quarts de travail de 8 heures. Plusieurs quarts partagent une
même période, mais chaque quart exige un salaire particulier. Nous souhaitons savoir combien d’employés doit-on
affecter à chaque quart de travail de façon à minimiser le total des salaires versés, en respectant le nombre minimum
d’employés pour chaque période. Les données du problème sont données dans la Table 2.2.

Période Quart 1 Quart 2 Quart 3 Quart 4 Quart 5 Minimum employés

6–8 X 48
8–10 X X 79
10–12 X X 65
12–14 X X X 87
14–16 X X 64
16–18 X X 73
18–20 X X 82
20–22 X 43
22–24 X X 52
0–6 X 15
Salaire 170 160 175 180 195

Table 2.2 – données du problème d’horaire de personnel

Comme variables, nous pouvons choisir xj comme le nombre d’employés affectés au quart j. L’objectif est
min z = 170x1 + 160x2 + 175x3 + 180x4 + 195x5 .
Pour chaque période, le nombre d’employés affectés aux différents quarts doit couvrir le minimum d’employés requis
pour cette période. Par exemple, pour la période de 14h à 16h, nous aurons:
x2 + x3 ≥ 64.
Le programme linéaire résultant est
min z = 170x1 + 160x2 + 175x3 + 180x4 + 195x5 ,
s.c. x1 ≥ 48,
x1 + x2 ≥ 79
x1 + x2 ≥ 65
x1 + x2 + x3 ≥ 87
x2 + x3 ≥ 64
x3 + x4 ≥ 73
x3 + x4 ≥ 82
x4 ≥ 43
x4 + x5 ≥ 52
x5 ≥ 15
xj ≥ 0, j = 1, 2, 3, 4, 5.
2.4. HYPOTHÈSES 15

Nous pouvons remarquer la présence de contraintes redondantes:

— si nous avons x1 + x2 ≥ 79, alors x1 + x2 ≥ 65 ; cette dernière contrainte est donc redondante et peut être
éliminée ;
— de même, x3 + x4 ≥ 82 implique x3 + x4 ≥ 73 ;
— les contraintes x1 ≥ 0, x4 ≥ 0, x5 ≥ 0, sont aussi redondantes, mais il n’y a aucun intérêt à les éliminer.
La solution optimale du problème est

x∗ = (x∗1 , x∗2 , x∗3 , x∗4 , x∗5 ) = (48, 31, 39, 43, 15).

Une difficulté est que le nombre d’employés doit toujours être entier, aussi l’hypothèse de continuité n’est pas
satisfaite dans le modèle (bien que la solution optimale dans ce cas particulier soit entière).
Exemple 5 (Réseau de distribution). Considérons deux usines (U 1 et U 2), un centre de distribution (CD), et
deux entrepôts (E1, E2). Chaque usine manufacture un certain nombre d’unités d’un même produit (offre). Chaque
entrepôt requiert un certain nombre d’unités de ce même produit (demande). Sur chaque lien (arc) du réseau, il y
a un coût de transport par unité de produit (coût unitaire). De plus, sur certains arcs, il y a une capactité sur le
nombre d’unités transportées. Le réseau considéré est représenté dans la Figure 2.8. L’objectif est de minimiser le
coût de transport total.

50 unités 900$/unité 30 unités

U1 E1
produites requises
40
0$
/u
nit
é

200$/unité 10 CD 300$/unité 200$/unité

10
ité 0$/
/ un un
$ ité
300 80

40 unités 60 unités
U2 E2
produites requises

Figure 2.8 – Réseau de distribution

Comme d’ordinaire, pour formuler le modèle, identifions en premier lieu les variables d’intérêt. Nous désignerons
par xi,j le nombre d’unités du produit transportées sur l’arc (i, j) (i.e. entre les sommets i et j) La fonction objectif
(en chiffrant le montant total en centaines de dollars) est

min z = 2xU1 ,U2 + 4xU1 ,CD + 9xU1 ,E1 + 3xU2 ,CD + xCD,E2 + 3xE1 ,E2 + 2xE2 ,E1 .

Les contraintes concernent tout d’abord la conservation du flot: en chaque sommet du réseau, le flot sortant moins
le flot entrant est égal au nombre d’unités produites dans le cas d’usine, l’opposé du nombre d’unités requises pour
un entrepôt, et zéro pour le centre de distribution. Nous devons en outre tenir compte de la capacité sur certains
arcs. Par exemple, pour l’arc (U 1, U 2), nous avons

xU1 ,U2 ≤ 10.

Il nous reste à préciser les traditionnelles contraintes de non-négativité:

xi,j ≥ 0.
16 CHAPITRE 2. PROGRAMMATION LINÉAIRE

Le modèle détaillé peut dès lors s’écrire de manière détaillée comme suit:
min z = 2xU1 ,U2 + 4xU1 ,CD + 9xU1 ,E1 + 3xU2 ,CD + xCD,E2 + 3xE1 ,E2 + 2xE2 ,E1 .
s.c. xU1 ,U2 + xU1 ,CD + xU1 ,E1 = 50,
− xU1 ,U2 + xU2 ,CD = 40,
− xU1 ,CD − xU2 ,CD + xCD,E2 = 0,
− xU1 ,E1 + xE1 ,E2 − xE2 ,E1 = −30,
− xCD,E2 − xE1 ,E2 + xE2 ,E1 = −60,
xU1 ,U2 ≤ 10,
xCD,E2 ≤ 80,
xU1 ,U2 ≥ 0, xU1 ,CD ≥ 0, xU1 ,E1 ≥ 0, xU2 ,CD ≥ 0, xCD,E2 ≥ 0, xE1 ,E2 ≥ 0, xE2 ,E1 ≥ 0.
Il s’agit d’un problème de flot à coût minimum, que nous étudierons plus en détails au Chapitre 5. La solution
optimale est
(x∗U1 ,U2 , x∗U1 ,CD , x∗U1 ,E1 , x∗U2 ,CD , x∗CD,E2 , x∗E1 ,E2 , x∗E2 ,E1 ) = (0, 40, 10, 40, 80, 0, 20).
Le nombre d’unités transportées doit toujours être une valeur entière, aussi l’hypothèse de continuité n’est pas
satisfaite dans ce cas. Néanmoins, dans ce cas particulier, la solution trouvée est entière. En fait, pour tout problème
de flot à coût minimum (avec paramètres à valeurs entières), il existe toujours une solution optimale entière.

2.5 La méthode du simplexe

Développée en 1947 par George Dantzig, la méthode du simplexe reste d’actualité pour résoudre des problèmes
de grande taille. Il s’agit d’une méthode algébrique basée sur la résolution de systèmes d’équations linéaires ; nous
nous intéressons ici uniquement aux systèmes d’équations linéaires avec un nombre n de variables supérieur au
nombre m d’équations:
Ax = b,
où x ∈ Rn , b ∈ Rm et A ∈ Rm×n . Dans ce cas, il y a trois possibilités:
1. aucune solution ;
2. une et une seule solution ;
3. une infinité de solutions.
Nous supposerons que toutes les variables sont positives. Le second cas, avec une et seule solution, ne peut survenir
que si n = m, et si la matrice A est inversible, ce qui revient à exiger que nous avons éliminé au préalable
toute équation pour s’écrire comme combinaison linéaire des autres équations. Si nous n’avons qu’une seule solution
admissible, elle est forcément optimale, par conséquent, nous ignorerons ce cas dans le reste du chapitre, et prendrons
n strictement supérieur à m.
Dans le cas où il y a une infinité de solutions, la méthode d’élimination de Gauss-Jordan permet d’identifier
trois types de variables:
— variables fixées ;
— variables dépendantes ;
— variables indépendantes.
Exemple 6 (Préliminaires à l’algorithme du simplexe). Considérons le système
x1 + x2 + x3 + x4 = 4,
x1 + x3 + x4 = 3,
x1 + x2 + 2x4 = 2.
Nous pouvons le réécrire sous la forme
x1 + + 2x4 = 1,
x2 = 1,
x3 − x4 = 2.
2.5. LA MÉTHODE DU SIMPLEXE 17

La variable x2 est fixée, comme elle ne peut prendre que la valeur 1, sans considération pour les autres variables.
A l’opposée, x4 est indépendante, comme elle peut prendre n’importe quelle valeur dans R. x1 et x3 sont quant à
elles dépendantes: le choix de x4 fixe leur valeur, et de plus, il n’est pas possible de les éliminer en combinant des
équations entre elles (au contraire de x4 ).

2.5.1 Solution de base

Il s’agit de la solution obtenue en fixant toutes les variables indépendantes à zéro. Nous qualifierons de variables
hors-base les variables indépendantes fixées à zéro. Les autres variables seront dites variables de base. Une solution
de base est réalisable (ou admissible) lorsque toutes les variables de base ont une valeur positive. Une solution de
base réalisable est dégénérée lorsqu’au moins une variable de base a la valeur 0. Il est possible de montrer qu’une
solution de base réalisable est un sommet du polyèdre.

Exemple 7 (Préliminaires à l’algorithme du simplexe: solution de base). Dans l’exemple, la solution de base est

x1 = 1, x2 = 1, x3 = 2.

Elle est réalisable et non dégénérée.

Pivot
Il est facile de changer le statut des variables par des opérations élémentaires.

Exemple 8 (Préliminaires à l’algorithme du simplexe: pivot).

x1 + + 2x4 = 1,
x2 = 1,
x3 − x4 = 2.

peut se réécrire comme

x1 + 2x3 = 5,
x2 = 1,
−x3 + x4 = −2.
Dans cette nouvelle solution de base, nous avons
— variable hors-base: x3 ;
— variables de base: x1 , x2 , x4 ;
— solution de base non réalisable: x1 = 5, x2 = 1, x4 = −2.

Le pivot est une opération consistant à remplacer une variable de base par une variable hors base pour obtenir
une nouvelle solution de base, dite adjacente.

Exemple 9 (Wyndor Glass). Rappelons que pour cet exemple, les contraintes fonctionnelles sont

x1 ≤4
2x2 ≤ 12
3x1 + 2x2 ≤ 18

Au lieu d’inégalités, nous voudrions des égalités. Pour ce faire, nous ajoutons des variables d’écart (en anglais “slack
variables”) supérieures à 0:
x1 +x3 =4
2x2 +x4 = 12
3x1 +2x2 +x5 = 18
Les variables hors-base sont x1 et x2 . Fixons-les à 0. Nous obtenons comme solution de base

(x1 , x2 , x3 , x4 , x5 ) = (0, 0, 4, 12, 18).

18 CHAPITRE 2. PROGRAMMATION LINÉAIRE

Effectuons un pivot, en remplaçant la variable hors-base par une des variables de base actuelles. Reste à déterminer
comment la choisir. Nous souhaitons en premier lieu que la nouvelle solution de base soit réalisable. Dans cette
solution de base, on aura toujours x1 = 0, et une des variables d’écart deviendra une variable hors-base, donc
prendra la valeur 0. En exploitant le système linéaire précédent et la positivité des variables, nous avons
x1 = 4 −x3 ≥0 4 ≥ x3
x4 = 12 −2x2 ≥0 ⇔ 12 −2x2 ≥0
x5 = 18 −3x1 −2x2 ≥0 18 −2x2 ≥0
En exploitant les inégalités
x4 = 12 − 2x2 ≥ 0
x5 = 18 − 2x2 ≥ 0,
nous obtenons
x2 ≤ 12/2 = 6
x2 ≤ 18/2 = 9.
Par conséquent, en posant x2 = 6, nous obtenons x4 = 0, alors que si nous augmentons d’avantage x2 , la solution
devient non réalisable. Nous effectuons comme pivot le remplacement de la variable de base x4 (qui deviendra hors-
base) par x2 . Nous obtenons alors le système suivant:
x1 +x3 =4
x2 + 21 x4 =6
3x1 −x4 +x5 = 6,
qui donne la solution de base:
(x1 , x2 , x3 , x4 , x5 ) = (0, 6, 4, 0, 6).
Nous effectuons ensuite un pivot pour que la variable x1 entre dans la base (i.e. devienne une variable de base).
Comme x4 = 0,
x3 = 4 −x1 ≥0 x ≤ 4
⇔ 1
x5 = 6 −3x1 ≥ 0 x1 ≤ 2
En posant x1 = 2, nous obtenons x5 = 0. Le pivot revient à remplacer la variable de base x5 par x1 . Le système
obtenu est alors :
x3 + 31 x4 − 31 x5 = 2
x2 + 12 x4 =6
x1 − 13 x4 + 31 x5 = 2
La solution de base correspondante est
(x1 , x2 , x3 , x4 , x5 ) = (2, 6, 2, 0, 0).

2.5.2 Interprétations
Interprétation géométrique
Une solution de base réalisable correspond à un point extrême du domaine réalisable. Un pivot correspond à un
déplacement d’un point extrême à un autre qui lui est adjacente, i.e. toutes les variables hors-base sauf une sont les
mêmes. La méthode du simplexe peut se résumer comme suit.
1. Nous débutons avec une solution de base réalisable initiale (un point extrême).
2. A chaque itération, nous effectuons un pivot, passant ainsi à une solution de base réalisable adjacente (un
point extrême adjacent).
3. L’algorithme s’arrête lorsqu’il identifie une solution de base réalisable optimale (un point extrême correspon-
dant à une solution optimale).
2.5. LA MÉTHODE DU SIMPLEXE 19

x2

3x1 + 2x2 = 18
x1 = 4
x1 = 0

(0, 6) (2, 6) (4, 6) 2x2 = 12

domaine
réalisable

x2 = 0
(0, 0) x1
(4, 0) (6, 0)

Figure 2.9 – exemple Wyndor Glass

Interprétation des variables d’écart

Exemple 10. Dans la solution optimale du problème Wyndor Glass, nous avons x3 = 2, x4 = x5 = 0. Cela indique
que les deux dernières ressources (temps aux usines 2 et 3) sont pleinement utilisées. Une partie de la première
ressource (temps à l’usine 1) n’est pas utilisée: 2 heures

2.5.3 Critère d’optimalité

Exemple 11 (Wyndor Glass: critère d’optimalité). Exprimons l’objectif en fonction des variables hors-base dans
la solution optimale Rappelons que dans cette solution, nous avons

x2 + 21 x4 =6
x1 − 13 x4 + 13 x5 =2

Après substitution dans l’objectif, nous obtenons

max z = 3x1 + 5x2

   
1 1 1
= 3 2 + x4 − x5 + 5 6 − x4
3 3 2
3
= 36 − x4 − x5 .
2
Toute solution réalisable (x1 , x2 , x3 , x4 , x5 ) satisfait
3
z = 36 + 0x1 + 0x2 + 0x3 − x4 − x5 ≤ 36.
2
Le critère d’optimalité s’énonce ainsi comme suit: étant donné que l’objectif s’exprime uniquement en fonction
des variables hors-base de la solution de base réalisable courante, si les coefficients de ces variables dans l’objectif
sont tous négatifs ou nuls, alors la solution de base réalisable courante est optimale. Les coefficients des variables
hors-base dans l’objectif sont appelés coûts réduits (ou coûts relatifs).
20 CHAPITRE 2. PROGRAMMATION LINÉAIRE

Si au moins un coût réduit est positif pour la solution de base courante, nous n’avons pas encore atteint une
solution optimale. Il faut par conséquent effectuer au moins un pivot, mais quelle variable doit-on faire entrer dans
la base. Regardons celle dont le coût réduit est le plus grand parmis toutes les variables hors-base, vu que cette
variable fournit la plus grande augmentation marginale (par unité) de la valeur de l’objectif. Ce qui ne signifie
toutefois pas la plus grande augmentation globale !
Vu que nous faisons entrer une variable dans la base, nous devons également choisir la variable qui va sortir de
la base en tentant de garder toutes les variables non négatives. Supposons que xj est la variable d’entrée. Chaque
variable de base xi s’exprime alors en fonction de la variable d’entrée (puisque les autres variables hors-base sont
nulles):
xi = bi − aij xj .
Dans cette expression, les coefficients bi et aij sont obtenus suite à plusieurs pivots. Mais nous avons nécessairement
bi positif (en partant d’un b positif). Pour que toutes les variables demeurent non négatives suite au pivot, nous
devons avoir
xi = bi − aij xj ≥ 0,
soit, en d’autres termes,
aij xj ≤ bi .
Si aij est négatif, cette inégalité ne limite pas l’augmentation de xj . Si cette condition est satisfaite pour tous les
i, nous pouvons donc augmenter indéfiniment xj ; l’objectif est non borné. Si aij est strictement positif, l’inégalité
limite l’augmentation de xj . Nous prendrons comme variable de sortie celle qui atteint
 
bi
min aij > 0 .
aij

Algorithm 2.1: Simplexe

1. Obtenir une solution de base réalisable.
2. Vérifier le critère d’optimalité: si les coûts réduits de toutes les variables hors-base sont négatifs ou nuls,
stop.
3. Choisir la variable xj , soit celle qui a le coût réduit le plus élevé.
4. Déterminer la variable de sortie:  
bi
min aij > 0 .
aij
5. Effectuer un pivot et déterminer une nouvelle solution de base réalisable. Retour à l’étape 2.

Exemple 12 (Production). Une entreprise fabrique quatre produits. La fabrication de chaque produit nécessite une
certaine quantité de ressources. Les ressources consommées, les stocks des ressources et les bénéfices des produits
sont récapitulés dans la Table 2.3. Nous souhaitons établir un plan de production de façon à maximiser le chiffre

Produit 1 2 3 4 Stock
Ressource A 2 4 5 7 42
Ressource B 1 1 2 2 17
Ressource C 1 2 3 3 24
Profit 7 9 18 17 -

Table 2.3 – Données de production

d’affaires.
2.5. LA MÉTHODE DU SIMPLEXE 21

Soient x1 , x2 , x3 et x4 les quantités respectives de produits 1, 2, 3, 4. Nous avons le problème de programmation

linéaire sous forme standard

max z = 7x1 + 9x2 + 18x3 + 17x4

s.c. 2x1 + 4x2 + 5x3 + 7x4 ≤ 42
x1 + x2 + 2x3 + 2x4 ≤ 17
x1 + 2x2 + 3x3 + 3x4 ≤ 24
x1 , x2 , x3 , x4 ≥ 0.

Introduisons trois variables d’écart x5 , x6 , x7 , qui mesurent pour chaque ressource l’écart entre la quantité initiale-
ment disponible et la quantité consommée par le plan de fabrication donné par les variables de décision x1 , x2 , x3 ,
x4 :

max z = 7x1 + 9x2 + 18x3 + 17x4

s.c. 2x1 + 4x2 + 5x3 + 7x4 + x5 = 42
x1 + x2 + 2x3 + 2x4 + x6 = 17
x1 + 2x2 + 3x3 + 3x4 + x7 = 24
x1 , x2 , x3 , x4 ≥ 0.

Cette formulation permet d’exprimer facilement les variables d’écart comme fonctions affines des variables de déci-
sion:

x5 = 42 − 2x1 − 4x2 − 5x3 − 7x4

x6 = 17 − x1 − x2 − 2x3 − 2x4
x7 = 24 − x1 − 2x2 − 3x3 − 3x4
z = 7x1 + 9x2 + 18x3 + 17x4

Le tableau ci-dessus est appelé un dictionnaire. Les variables x5 , x6 , x7 sont les variables de base et x1 , x2 , x3 ,
x4 sont des variables hors-base. La solution de base associée au dictionnaire est obtenue en donnant la valeur 0 à
toutes les variables hors-base. La solution basique correspondant au dictionnaire ci-dessus est donc

(x1 , x2 , x3 , x4 , x5 , x6 , x7 ) = (0, 0, 0, 0, 42, 17, 24).

Le bénéfice correspondant est z = 0 ; on ne produit rien, on ne consomme aucune ressource et on ne gagne rien.
Néanmoins, c’est une solution réalisable, car elle satisfait toutes les contraintes.
En partant de cette solution basique, nous cherchons à améliorer le bénéfice. Nous sélectionnons la variable
hors-base de coût réduit maximum, i.e. x3 , et gardons les autres variables hors-base à zéro. Il est évident que si
on fait croître x3 à partir de 0, les autres variables hors-base restant nulles, la valeur de la fonction z croît. Nous
cherchons à augmenter au plus la fonction objectif tout en gardant la solution reste réalisable. Les contraintes sur
l’augmentation de x3 sont:

x5 ≥ 0 ⇒ 42 − 5x3 ≥ 0 ⇒ x3 ≤ 8.4
x6 ≥ 0 ⇒ 17 − 3x3 ≥ 0 ⇒ x3 ≤ 8.5
x7 ≥ 0 ⇒ 24 − 3x3 ≥ 0 ⇒ x3 ≤ 8.

La plus restrictive de ces contraintes est x3 ≤ 8. L’interprétation géométrique est la suivante: en partant du sommet
(0, 0, 0, 0) du polyèdre des contraintes, nous nous déplaçons sur une arête de ce polyèdre. Le premier hyperplan
rencontré est x7 = 0. Nous arrivons alors dans un nouveau sommet, à l’intersection des hyperplans x1 = 0, x2 = 0,
x4 = 0, x7 = 0.
Nous allons faire un changement de dictionnaire, i.e. un pivot, en échangeant les rôles de x3 et x7 . Nous utilisons
la troisième équation du premier dictionnaire pour exprimer x3 en fonction de x1 , x2 , x4 et x7 :
1 2 1
x3 = 8 − x1 − x2 − x4 − x7 .
3 3 3
22 CHAPITRE 2. PROGRAMMATION LINÉAIRE

Nous remplaçons ensuite x3 par cette expression dans les autres équations du dictionnaire:
1 2 5
x5 = 2 − x1 − x2 + x7 − 2x4
3 3 3
1 1 2
x6 = 1 − x1 + x2 + x7
3 3 3
1 2 1
x3 = 8 − x1 − x2 − x7 − x4
3 3 3
z = 144 + x1 − 3x2 − 6x7 − x4.
La variable x7 est sortie de la base et la variable x3 est entrée dans la base. La nouvelle base est par conséquent x5 ,
x6 , x3 . La solution basique associée au nouveau dictionnaire est x1 = x2 = x7 = x4 = 0, x5 = 2, x6 = 1, x3 = 8.
Elle correspond au sommet de coordonnées (0, 0, 8, 0) du polyèdre de contraintes. Cette solution définit un plan de
production beaucoup plus intéressant: nous fabriquons 8 unités du produit 3 (x3 = 8) en consommant entièrement
la ressource C (x7 = 0). Il nous reste x5 = 2 unités de la ressource A et x6 = 1 unité de B. Le bénéfice est z = 144.
Dans la nouvelle expression de la fonction z, nous voyons que seule la variable x1 a un coefficient positif. Nous
décidons de faire entrer x1 en base, et ainsi de parcourir une nouvelle arête du polyèdre des contraintes. Nous avons
les limites suivantes sur l’augmentation de x1 :
x5 ≥ 0 ⇒ x1 ≤ 6x6 ≥ 0 ⇒ x1 ≤ 3x3 ≥ 0 ⇒ x1 ≤ 24
d’où x1 ≤ 3 Nous faisons sortir x6 de la base faisons entrer x1 à sa place. Nous obtenons le dictionnaire suivant:
x5 = 1 + x6 − x2 + x7 − 2x4
x1 = 3 − 3x6 + x2 + 2x7
x3 = 7 + x6 − x2 − x7 − x4
z = 147 − 3x6 − 2x2 − 4x7 − x4 .
La solution de base associée est x6 = x2 = x7 = x4 = 0, x5 = 1, x1 = 3, x3 = 7, correspond au sommet (3, 0, 7, 0)
du polyèdre des contraintes. Elle définit le plan de production suivant: nous fabriquons 3 unités du produit 1 et 7
unités du produit 3. Il ne nous reste qu’une unité de la ressource A. Le bénéfice est z = 147. De plus, tous les coûts
réduits sont négatifs, aussi la solution est optimale.

2.5.4 Adaptation à d’autres formes de modèles

Tout modèle de programmation linéaire peut se ramener à la forme suivante:
n
X
max cj x j
j=1
Xn
s.c. aij xj + xn+i = bi , i = 1, 2, . . . , m
j=1

xj ≥ 0, j = 1, 2, . . . , n,
xn+i ≥ 0, i = 1, 2, . . . , m.
Nous supposons de plus que bi est positif, i = 1, 2, . . . , m. Sous cette forme, il est facile d’initialiser la méthode
du simplexe en ajoutant des variables d’écart, et en les prenant comme variables de base. Cela revient de plus
à considérer l’origine comme solution initiale, et il est facile d’effectuer les opérations de pivot. Le situation se
complique avec d’autres formes fonctionnelles pour les contraintes, en particulier dans la recherche d’une solution
de base initiale.

Transformation du min au max

Supposons que nous devions minimiser l’objectif au lieu de le maximiser. Nous utilisons la propriété
n
X n
X
min cj xj = − max − cj xj .
j=1 j=1
2.5. LA MÉTHODE DU SIMPLEXE 23

Nous résolvons le problème de maximisation en changeant les signes des coefficients dans l’objectif. La valeur
optimale du problème de minimisation est l’opposé de celle du problème de maximisation.

Transformation des inégalités en égalités

Pn
Si j=1 aij xj ≤ bi , il y a deux cas:
1. bi ≥ 0: nous ajoutons une variable d’écart positive xni :
n
X
aij xj + xn+i = bi .
j=1

2. bi < 0: nous multiplions l’inégalité par -1, pour se ramener au cas développé ci-dessous.
Pn
Si j=1 aij xj ≥ bi , nous avons à nouveau deux cas possibles:
1. bi ≤ 0: nous multiplions l’inégalité par -1, pour se ramener au cas développé précédemment.
2. bi > 0: nous soustrayons une variable de surplus positive x0i :
n
X
aij xj − x0i = bi .
j=1

Nous nous ramenons alors au cas d’ajout de variables artificielles.

2.5.5 Obtention d’une base admissible initiale

Ajout de variables artificielles
P
Si j=1,...,n aij xj = bi et qu’aucune variable n’est isolée (une variable est isolée si elle à coefficient 1 dans cette
équation et à coefficient 0 dans les autres):
1. nous ajoutons une variable artificielle xn+1 ≥ 0, dont l’unique but est de fournir une variable de base initiale
pour cette contrainte ;
2. nous lui associons un profil très négatif: −M
X
max cj xj − M xn+i
j=1,...,n

s.c. . . .
X
aij xj + xn+i = bi
j=1,...,n

...

Le seul but des variables artificielles est de pouvoir initialiser l’algorithme du simplexe en produisant une solution
de base réalisable initiale. Si le problème est réalisable, nous devrions avoir xn+i = 0. Deux méthodes peuvent être
considérées.
1. La méthode du grand M consiste à optimiser en utilisant une fonction objective formée de la fonction de coût
initiale et de la somme, très fortement pénalisée, des variables artificielles.
2. La méthode à deux phases se déroule comme suit:
Phase 1 Trouver une solution réalisable en minimisant la somme des variables artificielles.
Phase 2 Optimiser en revenant à la fonction de coût initial à partir de la solution initiale trouvée dans la
phase 1.
24 CHAPITRE 2. PROGRAMMATION LINÉAIRE

Exemple 13. Considérons le programme

min z = 0.4x1 + 0.5x2

s.c. 0.3x1 + 0.1x2 ≤ 2.7
0.5x1 + 0.5x2 = 6
0.6x1 + 0.4x2 ≥ 6
x1 ≥ 0, x2 ≥ 0.

Nous devons tout d’abord transformer le système de contraintes afin de pouvoir traiter un système linéaire:

0.3x1 + 0.1x2 + xs1 = 2.7

0.5x1 + 0.5x2 = 6
0.6x1 + 0.4x2 − xs2 = 6
xs1 ≥ 0, xs2 ≥ 0.

En ajoutant les variables artificielles, nous obtenons

0.3x1 + 0.1x2 + xs1 = 2.7

0.5x1 + 0.5x2 + xa2 = 6
0.6x1 + 0.4x2 − xs2 + xa3 = 6
xs1 ≥ 0, xs2 ≥ 0, xa2 ≥ 0, xa3 ≥ 0.

Afin de démarrer la méthode du simplexe, nous appliquons une des deux méthodes précédentes. Méthode à deux
phases:
Phase 1 minimiser xa2 + xa3 jusqu’à obtenir une valeur optimale nulle (si le programme linéaire à une solution
réalisable).
Phase 2 minimiser 0.4x1 + 0.5x2 .
Méthode du grand M :
min 0.4x1 + 0.5x2 + M xa2 + M xa3

2.5.6 Variables à valeurs quelconques

−
Si une variable xj peut prendre des valeurs négatives, nous introduisons deux variables x+
j ≥ 0 et xj ≥ 0. Nous
posons alors
−
xj = x+j − xj .

Une autre possible, si xj ≥ Lj , où Lj est une constante négative, consiste à poser

x+
j = xj − Lj ≥ 0.

2.6 Dualité
2.6.1 Motivation
Exemple 14 (Wyndor Glass: dualité). Supposons qu’une compagnie partenaire de Wyndor Glass, appelée Dual
Glass, aimerait louer du temps aux usines afin de fabriquer des lots de produits. Quel prix horaire pour chaque
usine devrait-elle demander de telle sorte que le résultat soit équitable, soit aucun profit ni perte pour aucun des
deux partenaires ?
Dénotons les nouvelles variables de décision par yi , le prix de location horaire à l’usine i, et reprenons les
données du problèmes telles qu’exprimé précédemment dans la Table 2.1.
2.6. DUALITÉ 25

Produit 1 Produit 2 Capacité de

Temps de prod. (h) Temps de prod. (h) production (h)
Usine 1 1 0 4
Usine 2 0 2 12
Usine 3 3 2 18
Dual Glass cherche a minimiser le prix total qu’elle devra payer pour le temps loué aux trois usines. Le prix total
pour chaque usine peut etre exprimé comme
temps de production maximum (h) × prix
pour louer du temps ($/h). L’objectif est par conséquent
min w = 4y1 + 12y2 + 18y3 .
Les contraintes assurent que le prix total associé à la fabrication d’un lot de chaque produit ne doit pas être inférieur
au profit ($/lot) qu’en retire Wyndor Glass. Le prix total associé à un produit peut être exprimé comme le temps
consacré à la production de chaque lot (h/lot) multiplié par le prix pour louer du temps ($/h). La contrainte associée
au produit 1 peut s’exprimer comme
y1 + 3y3 ≥ 3.
La contrainte associée au produit 2 est
2y2 + 2y3 ≥ 5.
Nous obtenons ainsi le modèle pour Dual Glass, appelé modèle dual:
min w = 4y1 + 12y2 + 18y3
s.c. y1 + 3y3 ≥ 3
2y2 + 2y3 ≥ 5
y1 , y2 , y3 ≥ 0.
Pour rappel, le modèle pour Wyndor Glass, dit modèle primal est
max z = 3x1 + 5x2
s.c. x1 ≤ 4
2x2 ≤ 12
3x1 + 2x2 ≤ 18
Rappelons que pour la solution de base optimale du problème Wyndor Glass, l’objectif s’écrit:
3
z = 36 − x4 − x5 .
2
x4 et x5 sont les variables hors base et les coefficients − 32 et −1 sont leurs coûts réduits. Si on augmente la valeur
de x4 de une unité, le profit diminue de 23 . x4 est également la variable d’écart associée à la contrainte de ressource
pour l’usine 2:
2x2 + x4 = 12.
Augmenter x4 d’une unité signifie par conséquent diminuer le terme de droite correspondant de 1. En effet, si
Wyndor Glass loue a Dual Glass une heure de temps de production à l’usine 2:
— la capacité à l’usine 2 diminue de 1 heure (diminution d’une unité du terme de droite) ;
— la valeur de l’objectif diminue de 23 .
Pour retrouver un profit total égal, il faudra donc demander un prix de 32 (1500$) pour chaque heure de temps louée
à l’usine 2. De manière générale, la solution optimale du dual est donnée par l’opposé des coûts réduits des variables
d’écart (aussi appelés multiplicateurs optimaux, ou shadow prices en anglais). Dans notre exemple:
3
y1 = 0, y2 = , y3 = 1.
2
Le prix de la variable y1 est fixé à 0:
26 CHAPITRE 2. PROGRAMMATION LINÉAIRE

— Wyndor Glass n’exige rien pour une heure louée à l’usine 1, ce qui se justifie par le fait qu’il lui reste du
temps de production non utilisé (la variable d’écart x3 est strictement positive).
— Un prix strictement positif ferait augmenter le profit, et la solution ne serait plus équitable.
Les prix des autres variables est strictement positif:
— puisque le temps de production est utilisé pleinement, louer une heure à Dual Glass revient à perdre une
heure de production, donc à réduire le profit total ;
— pour retrouver le même profit, il faut exiger un prix égal au multiplicateur optimal, égal à l’opposé du coût
réduit.
Nous pouvons constater des écarts complémentaires: pour tout i, nous avons

xn+i .yi = 0.

Tout modèle de programmation linéaire possède un dual. A tout programme linéaire sous forme standard

max z = cT x
s.c. Ax ≤ b,
x ≥ 0,

nous associons le programme dual

min w = bT y
s.c. AT y ≥ c,
y ≥ 0,

Nous obtenons ainsi un couple primal-dual, avec les relations suivantes entre les deux modèles.

Primal Dual
Variable Contrainte
Contrainte Variable
max min
Profit unitaire Terme de droite
Terme de droite Coût unitaire
Ligne Colonne
Colonne Ligne
Contrainte ≤ Contrainte ≥

Si un modele de programmation linéaire possède une solution optimale, il en est de même pour son dual, et les
valeurs optimales des deux modèles sont égales.
La solution optimale du dual correspond aux multiplicateurs optimaux. Nous pouvons les lire directement dans
le tableau optimal du simplexe: ce sont les coefficients dans la ligne correspondant à l’objectif. Le coût réduit
(correspondant à l’opposé du coefficient dans la ligne de l’objectif) mesure la variation de l’objectif entraîné par
une augmentation d’une unité de la valeur de la variable hors-base associée.

2.6.2 Analyse de sensibilité

En général, le coût réduit d’une variable hors-base indique le changement dans l’objectif apporté par une
augmentation d’une unité de la valeur de cette variable. Pour les variables d’écart, ce principe peut se formuler
ainsi: le coût réduit d’une variable d’écart hors-base indique le changement dans l’objectif apporté par une diminution
d’une unité du terme de droite associé. Ceci est un exemple d’analyse de sensibilité: un paramètre (ici, un terme de
droite) est modifié et nous mesurons la sensibilité de la solution optimale à ce changement.
Nous pouvons aussi mesurer la sensibilité de la solution optimale à un changement d’un terme de droite ou d’un
coefficient dans l’objectif en resolvant à nouveau le modèle modifié.

Exemple 15 (Wyndor Glass: sensibilité). Une analyse de sensibilité complète nous apprendrait ce qui suit.
2.6. DUALITÉ 27

— Les multiplicateurs optimaux (Shadow Prices) sont 0, 1500 et 1000.

— La capacité à l’usine 3 peut prendre n’importe quelle valeur entre 12h et 24h sans changer la solution optimale
du dual.
— Le profit unitaire pour le produit 1 peut prendre n’importe quelle valeur entre 0 et 7500$ sans changer la
solution optimale du primal.

Exemple 16. Considérons un fermier cultivant trois types de plantes (maïs, blé, coton), en utilisant des ressources
de terre et de main d’oeuvre. Les revenus nets venant de la culture du maïs, du blé et du coton sont de $109/acre,
$90/acre et $115/acre, respectivement. Un acre de maïs exige 6 heures de travail, un acre de blé, 4 heures, et un
acre de coton, 8 heures. Le fermier a à disposition 100 acres de terre et 500 heures de travail ; il souhaite maximiser
le revenu net de ses plantations.
Nous obtenons le programme mathématique linéaire

max z = 109x1 + 90x2 + 115x3

s.c. x1 + x2 + x3 ≤ 100,
6x1 + 4x2 + 8x3 ≤ 500,
x1 , x2 , x3 ≥ 0,

ce qui peut se traduire en GAMS comme suit:

SETS

Crop Crop production

/ Corn Corn production
Wheat Wheat production
Cotton Cotton production /

Resource Elément de ressources utilisées dans le modèle

/ Land Land used by crop
Labor Labord used by crop /;

VARIABLES
Profit Net income from crops;

POSITIVE VARIABLES
Production(Crop) Production by crop;

SCALAR LandAvailable Total Land / 100 /;

PARAMETER

Revenue(Crop) Revenues from crop production

/ Corn 109
Wheat 90
Cotton 115 /

ResourceAvailable(Resource) Resource availability

/ Land 100
Labor 500 /;

TABLE ResourceUse(Resource, Crop) Resource used in the model

Corn Wheat Cotton
Land 1 1 1
Labor 6 4 8 ;
28 CHAPITRE 2. PROGRAMMATION LINÉAIRE

EQUATIONS
Objective Maximimze farm income
ResourceEq(Resource) Resource Constraint;

Objective..
Profit =E= SUM(Crop, Revenue(Crop)*Production(Crop));

ResourceEq(Resource)..
SUM(Crop, ResourceUse(Resource,Crop)*Production(Crop))
=L= ResourceAvailable(Resource);

MODEL FarmIncome /ALL/;

SOLVE FarmIncome USING LP MAXIMIZING Profit;

L’exécution de GAMS sur ce fichier donne
---- EQU ResourceEq Resource Constraint

LOWER LEVEL UPPER MARGINAL

Land -INF 100.000 100.000 52.000

Labor -INF 500.000 500.000 9.500

LOWER LEVEL UPPER MARGINAL

---- VAR Profit -INF 9950.000 +INF .

Profit Net income from crops

---- VAR Production Production by crop

LOWER LEVEL UPPER MARGINAL

Corn . 50.000 +INF .

Wheat . 50.000 +INF .
Cotton . . +INF -13.000
La première partie de cet extrait indique que les deux contraintes, sur la terre et la main d’oeuvre, sont actives:
toutes les resources sont employées. Les variables duales optimales correspondantes peuvent par conséquent être
non-nulles, et représentent la valeur marginale de ces contraintes (actives). Nous voyons en particulier ici que la
terre a la plus grande valeur marginale.
La seconde partie donnent les quantités de plantes cultivées. Seuls le maïs et le blé sont produits, mais pas le
coton. Le dual contient trois contraintes comme nous avons avons trois variables primales. Il est facile de vérifier
que les deux premières contraintes sont actives, mais pas la troisième, mais dépasse la borne inférieure de 13 unités.
Nous pouvons otenir les valeurs des variables d’écart avec la ligne
OPTION SOLSLACK = 1;
avant la commande SOLVE, ce qui dans cas confirme que les contraintes sont actives.
LOWER SLACK UPPER MARGINAL

Land -INF . 100.000 52.000

Labor -INF . 500.000 9.500
2.7. NOTES 29

Si nous ajoutons la commande

DISPLAY Profit.L, Production.L, ResourceEq.M, ResourceUse;

nous obtenons en sortie:

---- 53 VARIABLE Profit.L = 9950.000 Net income from crops

---- 53 VARIABLE Production.L Production by crop

Corn 50.000, Wheat 50.000

---- 53 EQUATION ResourceEq.M Resource Constraint

Land 52.000, Labor 9.500

---- 53 PARAMETER ResourceUse Resource used in the model

Corn Wheat Cotton

Land 1.000 1.000 1.000

Labor 6.000 4.000 8.000

La syntaxe .L conduit à obtenir les valeurs optimales des variables primales (ou de la solution) quand appliqué à
une variable, et la valeur du terme de gauche d’une équation de contrainte, quand à appliqué à une contrainte. La
syntaxe .M permet d’obtenir le multiplicateur optimal associée à une équation, si appliqué à une équation. Dans le
cas d’une variable, il nous permet de connaître la valeur excédentaire pour la contrainte duale associée.

2.7 Notes
Ce chapitre se base essentiellement sur les notes de cours de Bernard Gendron, 2007. L’exemple 3 est tiré de
Hillier et Lieberman [2], Section 3.1. L’exemple de production est dû à Stefan Balev.
30 CHAPITRE 2. PROGRAMMATION LINÉAIRE
 Programmation non linéaire
3
Nous pouvons généraliser le programme linéaire introduit au chapitre précédent comme suit:

max f (x1 , . . . , xn )
x
s.c.g1 (x1 , . . . , xn ) ≤ b1 ,
..
.
gm (x1 , . . . , xn ) ≤ bm .

Si f et gi , i = 1, . . . , m, sont linéaires, nous retrouvons le programme linéaire précédemment introduit, mais il

n’existe pas de raison de se limiter à de tels cas.
En général, nous supposerons aussi que toutes ces fonctions sont deux fois dérivables.

3.1 Fonctions convexes et concaves

Soient a et b deux points dans Rn . Le segment de droite joignant ces deux points est l’ensemble des points

{x ∈ Rn | ∃ λ ∈ [0, 1] tel que x = a + λ(b − a) = λb + (1 − λ)a}.

Une fonction est convexe si tout point sur le segment de droite joignant f (a) et f (b) se trouve au-dessus du graphe
de f :
f (λb + (1 − λ)a) ≤ λf (b) + (1 − λ)f (a), (3.1.1)
pour toute paire de points (a, b). Si cette inégalité est stricte, la fonction f est dite strictement convexe. En
remplaçant le sens de l’inégalité dans (3.1.1), la fonction f est dite concave:

f (λb + (1 − λ)a) ≥ λf (b) + (1 − λ)f (a),

pour toute paire (a, b). Si cette dernière inégalité est stricte, la fonction f est strictement concave. La Figure 3.1
illustre une fonction strictement convexe et une fonction strictement concave.

Tests de convexité et de concavité

Supposons que f est une fonction deux fois dérivable d’une seule variable. Alors,

31
32 CHAPITRE 3. PROGRAMMATION NON LINÉAIRE

f (x) f (x)
fonction concave

fonction convexe
x x

Figure 3.1 – Fonctions convexe et concave

1. f est convexe si et seulement si

d2 f (x)
∀ x, ≥ 0;
dx2
2. f est strictement convexe si
d2 f (x)
∀ x, > 0;
dx2
3. f est concave si et seulement si
d2 f (x)
∀ x, ≤ 0;
dx2
4. f est strictement concave si
d2 f (x)
∀ x, < 0.
dx2
Supposons que f est une fonction deux dérivable de deux variables. Dans ce cas, nous utilisons la Table 3.1
comme conditions nécessaires.

Convexe Strict. convexe Concave Strict. concave

∂ 2 f (x1 ,x2 )
∂x21
≥0 >0 ≤0 <0
∂ 2 f (x1 ,x2 )
∂x22
≥0 >0 ≤0 <0
 2
∂ 2 f (x1 ,x2 ) ∂ 2 f (x2 ,x2 ) ∂ 2 f (x1 ,x2 )
∂x21
. ∂x2 − ∂x1 x2 ≥0 >0 ≥0 >0
2

Table 3.1 – Test de convexité et de concavité pour une fonction à deux variables deux fois dérivable

Ces tests se généralisent à des fonctions de plus de deux variables en utilisant le Hessien (la matrice des dérivées
secondes). Deux résultats supplémentaires sont utiles pour tester la convexité:
— une fonction qui s’exprime comme la somme de fonctions convexes (concaves) est convexe (concave) ;
— l’opposé d’une fonction concave est convexe et vice-versa.

3.1.1 Ensembles convexes

Un ensemble est convexe si, pour toute paire de points de l’ensemble, le segment de droite joignant ces points
est contenu dans l’ensemble. Il est possible de démontrer qu’un ensemble est convexe grâce aux propriétés suivantes:
— si f est convexe, alors {x | f (x) ≤ b} est convexe ;
3.1. FONCTIONS CONVEXES ET CONCAVES 33

— si f est concave, alors {x | f (x) ≥ b} est convexe ;

— l’intersection d’ensembles convexes est convexe.

Exemple 17 (Wyndor Glass). Modifions l’exemple Wyndor Glass en remplaçant certaines contraintes par une
contrainte non linéaire.

max z = 3x1 + 5x2

s.c. x1 ≤ 4
9x21 + 5x22 ≤ 216
x1 ≥ 0, x2 ≥ 0.

Le problème ainsi modifié peut se représenter comme sur la Figure 3.2. Remarquons que l’objectif est concave (car

x2

8
7
(2, 6)
6
5
4 z = 36 = 3x1 + 5x2
domaine
3
réalisable
2
1
0 x1
0 1 2 3 4 5 6

Figure 3.2 – exemple Wyndor Glass

linéaire) et que le domaine réalisable est convexe: c’est un modèle de programmation convexe. Dans le cas présent,
la solution optimale est sur la frontière du domaine réalisable, mais ne correspond pas à un coin (l’intersection de
deux contraintes).
Modifions encore l’exemple Wyndor Glass, mais cette fois en remplaçant l’objectif lineaire par une fonction non
lineaire, comme illustré sur la Figure 3.3:

max z = 126x1 − 9x21 + 182x2 − 13x22

x1 ≤ 4
2x2 ≤ 12
3x1 + 2x2 ≤ 18

Nous pouvons à nouveau remarquer que l’objectif est concave et que le domaine réalisable est convexe: c’est aussi
un modèle de programmation convexe. Ici, la solution optimale (8/3, 5) est sur la frontière du domaine réalisable,
mais ne correspond pas à un point extrême du domaine réalisable.
34 CHAPITRE 3. PROGRAMMATION NON LINÉAIRE

x2

9
8
7
6
(8/3, 5)
5
4
domaine z = 907
3 z = 857
réalisable z = 807
2
1
0 x1
0 1 2 3 4 5 6 7 8 9

Figure 3.3 – exemple Wyndor Glass - cas 2

Considérons le même exemple, mais avec une fonction objectif différente, tel que représenté sur la Figure 3.4:

max z = 54x1 − 9x21 + 78x2 − 13x22

x1 ≤ 4
2x2 ≤ 12
3x1 + 2x2 ≤ 18

Il s’agit à nouveau d’un modèle de programmation convexe, avec comme solution optimale (3, 3), laquelle se trouve
à l’intérieur du domaine réalisable. La fonction objectif est la somme de deux fonctions d’une seule variable. Si nous
annulons les dérivées de chacune de ces fonctions, nous obtenons la solution unique (3, 3), qui se trouve à l’intérieur
du domaine ; c’est donc nécessairement la solution optimale.
Revenons sur la fonction objectif initiale, mais introduisons une contrainte non-linéaire, qui définit un domaine
réalisable non-convexe, comme illustré sur la Figure 3.5:

max z = 3x1 + 5x2

s.c. x1 ≤ 4
x2 ≤ 7
8x1 − x21 + 14x2 − x22 ≤ 49
x1 ≥ 0, x2 ≥ 0.

Le problème ainsi modifié peut se représenter comme sur la Figure 3.5. Dans ce modèle de programmation non
convexe, nous remarquons la présence de deux maxima locaux:
— x est un maximum local si f (x) ≥ f (a) pour tout a réalisable suffisamment près de a ;
— (4, 3) et (0, 7) sont des maxima locaux, mais x∗ = (0, 7) est la maximum global: f (x) ≥ f (a) pour tout a
réalisable.
Les méthodes classiques de programmation non-linéaire permettent d’identifier un maximum local, mais pas néces-
sairement le maximum global, alors qu’en programmation convexe, tout maximum local est global.
3.1. FONCTIONS CONVEXES ET CONCAVES 35

x2

7
6
z = 117
5
z = 162
4 z = 189
z = 198
3
(3, 3)
2
1
0 x1
0 1 2 3 4 5 6 7

Figure 3.4 – exemple Wyndor Glass - cas 3

x2

8
7 (0.7) = maximum global
6
5
4 z = 35 = 3x1 + 5x2

3 (4.3) = maximum local

domaine
2 z = 27 = 3x1 + 5x2
réalisable
1
0 x1
0 1 2 3 4 5 6

Figure 3.5 – exemple Wyndor Glass - cas 4

36 CHAPITRE 3. PROGRAMMATION NON LINÉAIRE

3.2 Algorithmes
3.2.1 L’algorithme du simplexe dans le cas non-linéaire
Le succès de l’algorithme du simplexe en programmation linéaire a naturellement motivé le développement de
techniques similaires pour l’optimisation non-linéaire. Une telle approche fut proposée en par Nelder et Mead [7].
Néanmoins, la notion de point extrême est beaucoup plus délicate à exploiter en non-linéaire, en particulier en
l’absence de contraintes. Il convient de construire un polyhèdre artificielle, ce qui n’est pas immédiat, et en pratique,
la méthode souffre de l’absence de garantie de convergence, excepté pour des fonctions convexes à une variable [4].
Nous pouvons même observer une convergence vers un point non-stationnaire [6]. Un algorithme d’optimisation
sans garantie de convergence est appelé heuristique.s Cette méthode est pourtant encore fréquemment employée,
en particulier par les ingénieurs. La fonction fminsearch de Matlab 1 implémente ainsi cet algorithme. Malgré ce
succès, aucun algorithme d’optimisation non-linéaire moderne ne se fonde sur le simplexe, et il convient dès lors
d’éviter autant que possible de maintenir un tel usage. De meilleurs résultats seront obtenus en employant des outils
plus appropriés.

3.2.2 Optimisation sans contrainte

Considérons tout d’abord le cas d’un modèle de programmation non-linéaire dans lequel il n’y a aucune
contrainte:
max f (x1 , x2 , . . . , xn ).
Il est alors possible de montrer que si x∗ = (x∗1 , x∗2 , . . . , x∗n ) est un maximum local, alors

∂f (x∗ )
= 0, j = 1, 2, . . . , n.
∂xj

En d’autres mots, lorsque la dérivée s’annule en un point donné, ce point peut être un maximum local, mais l’inverse
n’est pas nécessairement vérifié: ce point peut aussi être un minimum local ou un point de selle. Par contre, si f
est concave, un point où la dérivée s’annule est nécessairement un maximum global. De plus, si f est strictement
concave, un tel maximum global est unique.
Dans le cas d’une fonction non concave, afin de trouver un maximum global, il faut
— identifier tous les maxima locaux ;
— identifier celui de plus grande valeur ;
— vérifier que la fonction est bornée supérieurement (sinon, il n’y a pas de maximum global).
Le problème est qu’identifier tous les maxima locaux peut être extrêmement difficile.

3.2.3 Méthode de la bissection

Considérons en premier lieu le cas le plus simple: la fonction objectif f est concave et comporte une seule variable.
Il est alors possible de montrer que, si x∗ est une solution optimale et il existe xi et xu tels que xi ≤ x∗ ≤ xu et
f (xi ) 6= f (x∗ ), f (xu ) 6= f (x∗ ), alors il existe a et b tels que a ≤ x∗ ≤ b et

df (x)
> 0, si x < a,
dx
df (x)
= 0, si x = x∗ ,
dx
df (x)
< 0, si x > b.
dx

Si f est strictement concave en x∗ , alors a = x∗ = b, comme illustré sur la Figure 3.6.

1. Nous nous référons ici à la version 7.9, Release 2009b, mise en circulation le 4 septembre 2009.
3.2. ALGORITHMES 37

f (x)
df (x)
dx =0

x
x∗

Figure 3.6 – Maximum d’une fonction strictement concave

Algorithm 3.1: Méthode de la bissection

1. Nous fixons d’abord une borne inférieur xi pour laquelle la dérivée en ce point est strictement positive.
2. Nous déterminons également une borne supérieure xu pour laquelle la dérivée en ce point est strictement
négative.
3. Si les deux bornes ne sont pas suffisamment près l’une de l’autre, i.e. la distance xu − xi est trop grande,
nous prenons le point milieu entre les deux bornes comme point candidat xc :

xi + xu
xc = ,
2
et nous passons à l’étape suivante. Sinon, c’est-à-dire si

xu − xi ≤ 2,

avec  > 0 suffisamment petit, nous nous arrêtons: xc est une approximation suffisamment précise de la
solution optimale.
4. Si la dérivée en xc est nulle, arrêt: xc est la solution optimale.
5. Si la dérivée en xc est positive, nous posons xi = xc .
6. Si la dérivée en xc est négative, nous posons xu = xc .
7. Retour à l’étape 3.

3.2.4 Méthode du gradient

Considérons une fonction concave f de plusieurs variables. La méthode de la bissection ne peut plus s’appliquer,
aussi nous devons nous tourner vers une meilleure alternative. Le gradient de f au point x0 est défini comme
 ∂f (x) 0
∂x1
 ∂f (x0 ) 
0
 ∂x2 
∇x f (x ) = 
 ..  .

 . 
∂f (x0 )
∂xn
38 CHAPITRE 3. PROGRAMMATION NON LINÉAIRE

Il est possible de montrer que le gradient correspond a une direction d’augmentation de la valeur de f . Ainsi, dans
la méthode du gradient, nous nous déplaçons dans la direction du gradient en tentant d’augmenter au maximum
la valeur de l’objectif. A partir d’un point initial x0 , nous effectuons un deplacement dans la direction du gradient
vers un nouveau point x:
x = x0 + s,
où s = t∗ ∇x f (x). Dans cette formule, t∗ est la solution du problème de maximisation suivant:

max f (x0 + t∇x f (x0 )).

t≥0

C’est un problème de maximisation d’une fonction concave d’une seule variable: nous pouvons par conséquent le
résoudre, par exemple en employant la méthode de la bissection
Nous posons ensuite x0 = x, et nous itérons ainsi jusqu’à ce que le gradient s’annule (ou presque). Nous avons
par conséquent l’algorithme suivant.

Algorithm 3.2: Méthode du gradient

1. Déterminer une solution initiale x0 . Posons k = 0. Nous fixons également une constante  > 0 suffisamment
petite.
2. Résoudre le problème suivant:
max f (xk + t∇x f (xk )). (3.2.1)
t≥0

3. Soit t∗ la solution optimale de (3.2.1). Posons

xk+1 = xk + t∗ ∇x f (xk ),

et incrémentons k.
4. Si
∂f (xk )
≤ , j = 1, 2, . . . , n,
∂xj
arrêt. Sinon, nous retournons en 2.

3.2.5 Optimisation sous contraintes

Dans le cas d’un modèle de programmation non-linéaire sans contrainte, nous avons vu la condition d’optimalité
suivante: si x∗ = (x∗1 , x∗2 , . . . , x∗n ) est un maximum local, alors

∂f (x∗ )
= 0, j = 1, . . . , n.
∂xj

Cette condition n’est plus valable lorsqu’il y a des contraintes, car une solution optimale peut se trouver sur la
frontière du domaine réalisable

Exemple 18. Nous souhaitons résoudre le programme

max f (x) = 6 − x − x2
x
s.c. x ≥ 0.

Il est facile de se convaincre (voir Figure 3.7) que la solution optimale est x∗ = 0, pourtant

df (0)
= −1 < 0.
dx
3.2. ALGORITHMES 39

f (x)

7
df (x)
6 dx <0
5
4
3
2
1
0 x
0 1 2 3 4

Figure 3.7 – Maximum d’une fonction concave sous contrainte de positivité

3.2.6 Conditions d’optimalité

Avec une contrainte de la forme x ≥ 0, les conditions d’optimalité s’énoncent ainsi: si x∗ = (x∗1 , x∗2 , . . . , x∗n ) est
un maximum local, alors
( ∂f (x∗ )
∂xj = 0, si x∗j > 0,
∂f (x∗ )
∂xj ≤ 0, si x∗j = 0.

Il est plus délicat de dériver les conditions dans le cas général, où nous avons également des contraintes sont de la
forme:
gi (x1 , x2 , . . . , xn ) ≤ bi , i = 1, 2, . . . , m.
Elles s’expriment en fonction des multiplicateurs de Lagrange ui associés à chaque contrainte.

Conditions KKT (Karush-Kuhn-Tucker)

Si x∗ = (x∗1 , x∗2 , . . . , x∗n ) est un maximum local, alors il existe m nombres u1 , u2 , . . . , um tels que
m
∂f (x∗ ) X ∂gi (x∗ )
− ui ≤0
∂xj i=1
∂xj
m
!
∂f (x∗ ) X ∂gi (x∗ )
x∗j − ui =0
∂xj i=1
∂xj
gi (x∗ ) − bi ≤ 0, i = 1, 2, . . . , m,
ui (gi (x∗ ) − bi ) = 0, i = 1, 2, . . . , m,
ui ≥ 0, i = 1, 2, . . . , m,

Les conditions KKT peuvent également s’exprimer en omettant les contraintes de non-négativité, lesquelles
peuvent être réecrites sous la forme
−xi ≤ 0,
ou être tout simplement absentes. La forme qui suit est par conséquent plus générale, et plus communément admise
en programmation non-linéaire (voir par exemple Nocedal et Wright [8]). Si x∗ = (x∗1 , x∗2 , . . . , x∗n ) est un maximum
40 CHAPITRE 3. PROGRAMMATION NON LINÉAIRE

local, alors il existe m nombres u1 , u2 , . . . , um tels que

m
∂f (x∗ ) X ∂gi (x∗ )
− ui =0
∂xj i=1
∂xj
gi (x∗ ) − bi ≤ 0, i = 1, 2, . . . , m,
ui (gi (x∗ ) − bi ) = 0, i = 1, 2, . . . , m,
ui ≥ 0, i = 1, 2, . . . , m,
Exemple 19 (Conditions KKT). Considérons le problème
max f (x1 , x2 ) = ln(x1 + 1) + x2
s.c. 2x1 + x2 ≤ 3,
x1 ≥ 0, x2 ≥ 0.
Nous pouvons donc considérer qu’il n’y a qu’une seule contrainte avec
g1 (x1 , x2 ) = 2x1 + x2 ,
et b1 = 3. Nous associons à cette contrainte un multiplicateur u1 ≥ 0. Outre les contraintes de non-negativité, nous
avons alors les conditions suivantes:
1
∗ − 2u1 ≤ 0,
x1 + 1
1 − u1 ≤ 0,
 
1
x∗1 − 2u1 = 0,
x∗1 +1
x∗2 (1 − u1 ) = 0,
2x∗1 + x∗2 − 3 ≤ 0,
u1 (2x∗1 + x∗2 − 3) = 0.
Nous obtenons u1 ≥ 1. Puisque x∗1 ≥ 0, nous en déduisons
1
− 2u1 < 0.
x∗1 +1
Par conséquent, x∗1 = 0. Puisque u1 6= 0, nous avons
2x∗1 + x∗2 − 3 = 0,
et par conséquent x∗2 = 3. Dès lors, u1 = 1. Les conditions KKT sont donc satisfaites en un seul point: (0, 3). Il
s’agit bien d’un maximum global, car la fonction objectif est concave et le domaine réalisable est convexe (modèle
de programmation convexe).
En utilisant la version plus générale des conditions KKT, nous avons les conditions
1
− 2u1 + u2 ≤ 0,
x∗1 +1
1 − u1 + u3 ≤ 0,
2x∗1 + x∗2 − 3 ≤ 0,
u1 (2x∗1 + x∗2 − 3) = 0,
−x1 ≤ 0,
−x2 ≤ 0,
u2 x1 = 0,
u3 x2 = 0.
Il est facile de vérifier que la solution de ce système est
x1 = 0, x2 = 3, u1 = 1, u2 = 1, u3 = 0.
3.3. NOTES 41

3.3 Notes
Ce chapitre se base en grande partie sur les notes de cours de Bernard Gendron, 2007.
42 CHAPITRE 3. PROGRAMMATION NON LINÉAIRE
The purpose of mathematical programming is insight, not numbers.
Arthur M. Geoffrion

Programmation mixte entière

4
4.1 Présentation générale
Certaines quantités ne peuvent s’écrire sous forme de nombres réels, issus d’un domaine continu. Au contraire,
certaines décisions sont par nature discrètes, et doivent se représenter à l’aide de nombres entiers. Considérons par
exemple une entreprise de transport, qui décide de renouvellement sa flotte de camions. Le nombre de camions à
acheter est un nombre naturel.
La présence de telles variables entières modifie profondément la nature des programmes sous-jacents. Lorsqu’un
problème est linéaire avec des variables entières, nous parlerons de programmation mixte entière. Si toutes les
variables sont entières, nous utiliserons la terminologie de de programmation (pure) en nombres entiers. Un exemple
de programme en nombres entiers a déjà été rencontré dans la Sectionn 1. Si les variables entières sont à valeurs 0
ou 1 (binaires), nous parlerons de programmatión 0–1 (binaire).
Nous pourrions bien entendu considérer le cas non-linéaire, mais les complications sont telles qu’à ce jour,
aucune méthode pleinement satisfaisante n’existe, même si d’intenses recherches sont actuellement conduites à ce
niveau, notamment au sein de l’entreprise IBM.

Exemple 20 (Problème du sac à dos). Considérons un cambrioleur muni d’un sac (unique) pour transporter son
butin. Son problème consiste à maximiser la valeur totale des objets qu’il emporte, sans toutefois dépasser une
limite de poids b correspondant à ses capacités physiques. Supposons qu’il y a n type d’objets que le voleur pourrait
emporter, et que ceux-ci sont en nombre tel que quelle que soit la nature de l’objet considéré, la seule limite au
nombre d’unités que le cambrioleur peut prendre est que le poids total reste inférieur à b. Si l’on associe à l’objet j
une valeur cj et un poids wj , la combinaison optimale d’objets à emporter sera obtenue en résolvant le programme
mathématique
n
X
max cj xj
x
j=1
Xn
s.c. wj xj ≤ b
j=1

xj ∈ N , j = 1, . . . , n.

Ici, xj ∈ N signifie que xj est un naturel, autrement dit un entier non négatif. Intuitivement, nous pourrions penser
que la solution consiste à choisir en premier lieu les objets dont le rapport qualité-poids est le plus avantageux,
quitte à tronquer pour obtenir une solution entière (nous ne pouvons pas diviser un objet). Cette solution peut

43
44 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

malheureusement se révéler sous-optimale, voire mauvaise, comme on le constate sur l’exemple suivant.

max 2x1 + 3x2 + 7x3

x
s.c. 3x1 + 4x2 + 8x3 ≤ 14,
x1 , x2 , x3 ∈ N .

Si nous oublions la contrainte d’intégralité, la solution, de valeur 12.25, est x = (0, 0, 14/8). En tronquant, on obtient
la solution entière x = (0, 0, 1) de valeur 7. Il est facile de trouver de meilleures solutions, par exemple x = (1, 0, 1).
La solution optimale du problème du cambrioleur peut s’obtenir en énumérant toutes les solutions admissibles
et en conservant la meilleure (voir Table 4.1, , où les solutions inefficaces laissant la possibilité d’ajouter un objet
n’ont pas été considérées).

x1 x2 x3 objectif
0 1 1 10
2 0 1 11
0 3 0 9
2 2 0 10
3 1 0 9
4 0 0 8

Table 4.1 – Problème du sac à dos: énumération des solutions

La solution optimale entière, x = (2, 0, 1), diffère passablement de la solution optimale linéaire. Cependant, il
est clair que cette technique d’énumération ne peut s’appliquer à des problèmes de grande taille.

Exemple 21 (California Mfg). Une entreprise doit choisir de nouveaux emplacements pour construire des usines et
des entrepôts. Elle a le choix entre deux emplacements: Los Angeles (LA) et San Fransisco (SF). Nous ne pouvons
construire un entrepôt que dans une ville où nous disposons d’une usine, et nous ne pouvons pas construire plus
d’un entrepôt. Nous associons à chaque construction (d’une usine ou d’un entrepôt dans chacun des lieux envisagés)
sa valeur estimée et son coût de construction. L’objectif est de maximiser la valeur totale estimée, en ne dépassant
pas une limite maximum sur les coûts. Les données du problème sont résumées dans la Table 4.2. Les variables sont
(
1 si la décision j est approuvée ;
xj =
0 si la décision j n’est pas approuvée.

L’objectif est de maximiser la valeur estimée totale:

max z = 9x1 + 5x2 + 6x3 + 4x4 .

Les contraintes fonctionnelles sont

Valeur estimée (millions $) Coût de construction (millions $)

Usine à LA 9 6
Usine à SF 5 3
Entrepôt à LA 6 5
Entrepôt à SF 4 2
Limite maximum - 10

Table 4.2 – Données du problème California Mfg

4.1. PRÉSENTATION GÉNÉRALE 45

1. la limite maximum sur les coûts de construction:

6x1 + 3x2 + 5x3 + 2x4 ≤ 10;

2. nous ne pouvons construire plus d’un entrepôt:

x3 + x4 ≤ 1;

3. l’entrepôt ne sera à LA que si l’usine est à LA:

x3 ≤ x1 ;

4. l’entrepôt ne sera à SF que si l’usine est à SF:

x4 ≤ x2 :

5. contraintes 0–1 (intégralité):

xj ∈ {0, 1}, j = 1, 2, 3, 4;
ou encore
0 ≤ xj ≤ 1 et xj entier, j = 1, 2, 3, 4.
Par conséquent, nous avons le programme mathématique
max z = 9x1 + 5x2 + 6x3 + 4x4
s.c. 6x1 + 3x2 + 5x3 + 2x4 ≤ 10;
x3 + x4 ≤ 1;
− x1 + x3 ≤ 0;
− x2 + x4 ≤ 0;
x1 , x2 , x3 , x4 ≤ 1;
x1 , x2 , x3 , x4 ≥ 0;
x1 , x2 , x3 , x4 entiers.
Le modèle illustre deux cas classiques d’illustrations de variables binaires:
— alternatives mutuellement exclusives: nous ne pouvons construire plus d’un entrepôt, i.e.
x3 + x4 ≤ 1;
— décisions contingentes: nous ne pouvons construire un entrepôt que là où nous avons construit une usine:
x3 ≤ x1 ; x4 ≤ x2 :
Exemple 22 (Localisation). Une entreprise envisage plusieurs sites de construction pour des usines qui serviront à
approvisionner ses clients. À chaque site potentiel i correspond un coût de construction ai , une capacité de production
ui , un coût de production unitaire bi et des coûts de transport cij des usines vers les clients. Soit yi une variable
binaire prenant la valeur 1 si un entrepôt est construit sur le site i, dj la demande de l’usine j et xij la quantité
produite à l’usine i et destinée au marché j (flot de i à j). Un plan de construction optimal est obtenu en résolvant
le programme
X X X XX
min ai yi + bi xij + cij xij
x
i i j i j
X
s.c. xij = dj ,
i
X
xij ≤ ui yi ,
j

xij ≥ 0,
yi ∈ {0, 1}.
46 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

Cette formulation contient deux éléments intéressants: un coût fixe (construction) modélisé par une variable binaire
yi ainsi qu’une contrainte logique forçant les flots provenant d’un site à être nuls si aucune usine n’est construite
en ce site. Notons aussi que certaines variables sont entières alors que d’autres (flots) sont réelles.

Exemple 23 (Contraintes logiques). Des variables binaires peuvent servir à représenter des contraintes logiques. En
voici quelques exemples, où pi représente une proposition logique et xi la variable logique (binaire) correspondante.

contrainte logique forme algébrique

p1 ⊕ p2 = vrai x1 + x2 = 1
p1 ∨ p2 ∨ · · · ∨ pn = vrai x1 + x2 + . . . + xn ≥ 1
p1 ∧ p2 ∧ . . . ∧ pn = vrai x1 + x2 + . . . + xn ≥ n (ou = n)
p1 ⇒ p2 x2 ≥ x1
p1 ⇔ p2 x2 = x1

Exemple 24 (Fonctions linéaires par morceaux). Considérons une fonction objectif à maximiser, pour laquelle
dans chaque intervalle [ai−1 , ai ] la fonction est linéaire, ce que nous pouvons exprimer par:

x = λi−1 ai−1 + λi ai
λi−1 + λi = 1,
λi−1 , λi ≥ 0,
f (x) = λi−1 f (ai−1 ) + λi f (ai ).

Nous pouvons généraliser cette formule sur tout l’intervalle de définition de la fonction f en contraignant les variables
λi à ne prendre que deux valeurs non nulles, et ce pour deux indices consécutifs. Ceci se fait en introduisant des
variables binaires yi associées aux intervalles de linéarité [ai−1 , ai ]:

n
X
x= λi ai ,
i=0
Xn
f (x) = λi f (ai ),
i=0
n
X
λi = 1,
i=0
λi ≥ 0, i = 0, . . . , n
λ0 ≤ y1 ,
λ1 ≤ y1 + y2,
.. .. ..
. . .
λn−1 ≤ yn−1 + yn ,
λ n ≤ yn ,
n
X
yi = 1(un seul intervalle “actif ”)
i=1
yi ∈ {0, 1}, i = 1, . . . , n.

Si la fonction f est concave, les variables binaires yi et les contraintes associées peuvent être éliminées de la
[Link] approche est particulièrement intéressante lorsqu’on optimise des fonctions de plusieurs variables
de la forme i fi (xi ).
4.2. CONTRAINTES MUTUELLEMENT EXCLUSIVES 47

4.2 Contraintes mutuellement exclusives

4.2.1 Deux contraintes
Prenons l’exemple de deux contraintes. L’une ou l’autre des deux contraintes doit être satisfaite, mais pas les
deux simultanément. Par exemple,
— soit 3x1 + 2x2 ≤ 18 ;
— soit x1 + 4x2 ≤ 16.
Soit M un très grand nombre ; le système précédent est équivalent à
— soit

3x1 + 2x2 ≤ 18,

x1 + 4x2 ≤ 16 + M ;

— soit

3x1 + 2x2 ≤ 18 + M,
x1 + 4x2 ≤ 16.

En introduisant une variable binaire y, nous obtenons le système équivalent

3x1 + 2x2 ≤ 18 + M (1 − y),

x1 + 4x2 ≤ 16 + M y.

La signification de cette variable est

— y = 1, si la première contrainte est satisfaite ;
— y = 0, si la seconde contrainte est satisfaite.
Nous avons de la sorte construit deux alternatives mutuellement exclusives.
Nous aurions pu aussi introduire deux variables binaires:
— y1 = 1, si la première contrainte est satisfaite ;
— y2 = 1, si la seconde contrainte est satisfaite.
Nous devons avoir
y1 + y2 = 1.
Afin de se ramener au modèle précédent, il suffit de poser

y1 = y;
y2 = 1 − y.

Il s’agit d’un cas particulier de la situation suivante: K parmi N contraintes doivent être satisfaites. Dans ce cas
plus général, nous introduisons N variables binaires.

4.2.2 K contraintes parmi N

Soit les N contraintes
fj (x1 , x2 , . . . , xn ) ≤ dj , j = 1, 2, . . . , N.
Nous introduisons N variables binaires, avec yj = 1 si la j e contrainte est satisfaite:

fj (x1 , x2 , . . . , xn ) ≤ dj + M (1 − yj ), j = 1, 2, . . . , N.

Il reste à spécifier que seulement K de ces contraintes peuvent êtres satisfaites:

N
X
yj = K.
j=1
48 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

4.2.3 Fonction ayant N valeurs possibles

Soit la contrainte
f (x1 , x2 , . . . , xn ) = d1 , ou d2 ou . . . ou dN .
Nous introduisons N variables binaires, avec yj = 1 si la fonction vaut dj . La contrainte s’écrit alors
N
X
f (x1 , x2 , . . . , xn ) = d j yj ,
j=1

avec
N
X
yj = 1.
j=1

Exemple 25 (Wyndor Glass). Supposons que le temps de production maximum à l’usine 3 n’est pas toujours 18h,
mais pourrait également être 6h ou 12h. Cette contrainte s’écrit alors

3x1 + 2x2 = 6 ou 12 ou 18.

Nous introduisons alors trois variables binaires

3x1 + 2x2 = 6y1 + 12y2 + 18y3 ,

y1 + y2 + y3 = 1.

4.2.4 Objectif avec coûts fixes

Supposons que le coût associé à un produit j est composé de deux parties:
1. un coût fixe initial kj encouru dès qu’une unité de j est produite ;
2. un coût cj proportionnel au nombre d’unités de j produites.
Le coût total associé à la production de xj unités est
(
kj + cj xj si xj > 0,
f (xj ) =
0 si xj = 0.

Supposons de plus que l’objectif consiste à minimiser la somme de n fonctions avec coûts fixes
n
X
min z = fj (xj ).
j=1

Nous introduisons alors n variables binaires:

(
1 si xj > 0,
yj =
0 si xj = 0.

L’objectif s’écrit alors

n
X
min z = cj xj + kj yj .
j=1

Les valeurs de xj et de yj dépendent l’une de l’autre: il s’agit d’un exemple de décisions contingentes. Nous devons
en particulier avoir une contrainte qui précise que xj = 0 si yj = 0. Toutefois, les deux variables ne sont plus
binaires, vu que xj peut être quelconque. Soit Mj une borne supérieure sur la valeur de xj . Nous pouvons écrire la
relation entre les deux variables de cette manière:

xj ≤ M j yj .

Ainsi,
4.2. CONTRAINTES MUTUELLEMENT EXCLUSIVES 49

— si yj = 0, alors xj = 0 ;
— si yj = 1, alors xj ≤ Mj ;
— si xj > 0, alors yj = 1 ;
— si xj = 0, alors toute soluton optimale satisfait yj = 0 lorsque kj > 0 (si kj = 0, la variable yj est inutile).
Nous obtenons par conséquent le programme
n
X
min z = cj x j + k j y j
j=1

s.c. xj ≤ Mj ,
yj ∈ {0, 1}, j = 1, 2, . . . , n.

4.2.5 Variables entières en variables 0–1

Soit x une variable entière générale bornée:
0 ≤ x ≤ u,
N N +1
et soit N l’entier tel que 2 ≤u≤2 . La représentation binaire de x est
N
X
x= 2j yj .
j=0

L’intérêt de cette transformation est que les méthodes de programmation 0–1 sont souvent plus efficaces que les mé-
thodes de programmation en nombres entiers. Elle engendre néanmoins une multiplication du nombre de variables.

Exemple 26. Nous considérons trois types de produits, et deux usines ; nous exprimons le profit par unité de
produit en milliers de dollars. Nous connaissons les ventes potentielles par produit (unités/semaine), et la capacité
de production par usine (h/semaine). Nous avons toutefois comme contrainte que pas plus de deux produits ne
peuvent être fabriqués, et une seule des deux usines doit être exploitée. Les données du problème sont résumées dans
la Table 4.3. Les variables sont xj , le nombre d’unités fabriquées du produit j. Pour représenter la contrainte “pas

Produit 1 Produit 2 Produit 3 Capacité de

temps de production temps de production temps de production production
(h/unité) (h/unité) (h/unité) (h/semaine)
Usine 1 3 4 2 30
Usine 2 4 6 2 40
Profit/unité 5 7 3 –
(1000$)
Ventes potentielles 7 5 9 –
(par semaine)

Table 4.3 – Exemple de production avec variables entières.

plus de deux produits”, nous devons introduire des variables 0–1:

(
1 si xj > 0;
yj =
0 si xj = 0.

Afin de représenter la contrainte “une seule des deux usines”, nous devons ajouter une variables 0–1:
(
1 si l’usine 1 est choisie;
y4 =
0 si sinon.
50 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

L’objectif est
max z = 5x1 + 7x2 + 3x3 .
Les ventes potentielles sont
x1 ≤ 7, x2 ≤ 5, x3 ≤ 9.
L’exigence interdisant d’avoir plus de deux produits se traduit mathématiquement par

y1 + y2 + y3 ≤ 2.

La relation entre les variables continues et les variables 0–1 s’exprime par

x1 ≤ 7y1 , x2 ≤ 5y2 , x3 ≤ 9y3 .

La contrainte portant sur l’utilisation d’une seule usine est

— soit 3x1 + 4x2 + 2x3 ≤ 30 ;
— soit 4x1 + 6x2 + 2x3 ≤ 40.
En utilisant la variable 0–1 (et M très grand), elle se traduit par

3x1 + 4x2 + 2x3 ≤ 30 + M (1 − y4 ),

4x1 + 6x2 + 2x3 ≤ 40 + M y4 .

En résumé, nous avons le modèle

max z = 5x1 + 7x2 + 3x3

s.c. x1 ≤ 7y1 , x2 ≤ 5y2 , x3 ≤ 9y3
y1 + y2 + y3 ≤ 2,
3x1 + 4x2 + 2x3 ≤ 30 + M (1 − y4 ),
4x1 + 6x2 + 2x3 ≤ 40 + M y4 ,
x1 , x2 , x3 ≥ 0,
yj ∈ {0, 1}, j = 1, 2, 3, 4.

Exemple 27. Nous considérons à nouveau trois types de produits, pour lesquels nous pouvons placer cinq annonces
publicitaires, avec un maximum de trois annonces par produit. L’estimation des revenus publicitaires est donnée
dans la Table 4.4, où les profits sont exprimés en millions de dollars. Les variables du problèmes sont le nombre

Nombre d’annonces Produit 1 Produit 2 Produit 3

0 0 0 0
1 1 0 -1
2 3 2 2
3 3 3 4

Table 4.4 – Revenus publicitaires.

d’annonces pour le produit i, dénoté par xi , mais l’hypothèse de proportionnalité est alors violée: nous ne pouvons
représenter l’objectif sous forme linéaire uniquement avec ces variables.
Prenons tout d’abord comme variables
(
1 si xi = j;
yij =
0 sinon.

L’objectif est
max z = y11 + 3y12 + 3y13 + 2y22 + 3y23 − y31 + 2y32 + 4y33 .
4.2. CONTRAINTES MUTUELLEMENT EXCLUSIVES 51

Nous utiliserons les 5 annonces disponibles:

3 X
X 3
jyij = 5.
i=1 j=1

Enfin, la définition des variables 0–1 donne

3
X
yij ≤ 1, i = 1, 2, 3.
j=1

Nous pourrions aussi prendre comme variables

(
1 si xi ≥ j;
yij =
0 sinon.

Autrement dit, nous avons remplacé l’égalité dans la première condition par une inégalité. Cette définition implique

xi = 0 ⇒ yi1 = 0, yi2 = 0, yi3 = 0;

xi = 1 ⇒ yi1 = 1, yi2 = 0, yi3 = 0;
xi = 2 ⇒ yi1 = 1, yi2 = 1, yi3 = 0;
xi = 3 ⇒ yi1 = 1, yi2 = 1, yi3 = 1.

Ce qui peut encore s’énoncer somme

yi(j+1) ≤ yij , i = 1, 2, 3, j = 1, 2.

Supposons que x1 = 3 (il y a trois annonces pour le produit 1). Le profit associé à cette valeur doit être 3. Mais
x1 = 3 veut aussi dire que chaque variable binaire associée au produit 1 vaut 1 ; comment dès lors comptabiliser
correctement la contribution de ces trois variables au profit ? La solution consiste à prendre comme profit associé à
la variable yij la différence cij+1 − cij , où cij est le revenu net si nous plaçons j annonce pour le produit i. Dans
notre exemple, le profit associé à
— y11 est 1-0 = 1 ;
— y12 est 3-1 = 2 ;
— y13 est 3-3 = 0 ;
Nous obtenons ainsi le programme mathématique suivant:

max z = y11 + 2y12 + 2y22 + y23 − y31 + 3y32 + 2y33

s.c. yi(j+1) ≤ yij , i = 1, 2, 3, j = 1, 2,
3 X
X 3
yij = 5,
i=1 j=1

yij ∈ {0, 1}, i = 1, 2, 3, j = 1, 2.

4.2.6 Problème de recouvrement

Exemple 28 (Affectation des équipages). Un problème important des compagnies aériennes consiste à constituer
de façon efficace des équipages pour ses vols. Pour un équipage donné, une rotation consiste en une succession de
services de vol débutant et se terminant en une même ville. Comme il y a un coût associé à chaque séquence de vols,
la compagnie cherche à minimiser les coûts d’affectation des équipages aux séquences tout en assurant le service sur
chacun des vols.
Considérons par exemple un problème avec 11 vols et 12 séquences de vols, dont les données sont décrites dans
la Table 4.5. Les variables sont
(
1 si la séquence de vols j est affectée;
xj =
0 sinon.
52 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

Vol | Séquence 1 2 3 4 5 6 7 8 9 10 11 12
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1 1
5 1 1 1 1
6 1 1 1
7 1 1 1 1 1
8 1 1 1 1
9 1 1 1
10 1 1 1 1
11 1 1 1 1 1
Coût 2 3 4 6 7 5 7 8 9 9 8 9

Table 4.5 – Affectation d’équipages.

L’objectif est

min z = 2x1 + 3x2 + 4x3 + 6x4 + 7x5 + 5x6 + 7x7 + 8x8 + 9x9 + 9x10 + 8x11 + 9x12 .

Nous devons affecter trois équipages

12
X
xj = 3.
j=1

Le service doit être assuré sur chacun des vols:

x1 + x4 + x7 + x10 ≥ 1
x2 + x5 + x8 + x11 ≥ 1
x3 + x6 + x9 + x12 ≥ 1
x4 + x7 + x9 + x10 + x12 ≥ 1
...

En pratique, il faut tenir compte des contraintes imposées par les conventions collectives, ce qui complique singuliè-
rement le problème.

Un problème de recouvrement d’ensemble met en oeuvre

— I: un ensemble d’objets (les vols dans l’exemple précédent) ;
— J : une collection de sous-ensembles de I (e.g. les séquences de vols) ;
— Ji , i ∈ I: les sous-ensembles dans J qui contiennent i.
Nous avons les variables binaires xj , prenant la valeur 1 si le sous-ensemble j est choisi, et 0 sinon. En considérant
un objectif linéaire, avec cj le coût associé au sous-ensemble j. Nous obtenons le programme
X
min cj xj
x
j∈J
X
s.c. xj ≥ 1, i ∈ I;
j∈Ji

xj ∈ {0, 1}, j ∈ J .
4.3. STRATÉGIES DE RÉSOLUTIONS 53

4.3 Stratégies de résolutions

4.3.1 Relaxation linéaire
Il est tentant “d’oublier” les contraintes d’intégralité, et de résoudre le problème en nombres continus ainsi
produit. Nous parlerons alors de relaxation. Ainsi, nous pourrons construire la relaxtion en programme linéraire
d’un programme mixte entier. Une fois le programme relâché résolu, nous pourrions arrondir aux valeurs entières
les plus proches. Dans certains cas, cela peut fonctionner, mais l’exemple du sac à dos montre que ce n’est pas
toujours vérifié. Cette méthode par arrondissement est même parfois désastreuse.

Exemple 29. Considérons le programme

max z = x2
1
s.c. − x1 + x2 ≤ ,
2
7
x1 + x2 ≤ ,
2
x1 , x2 ≥ 0 et entiers.

La relaxation en programme linéaire donne

max z = x2
1
s.c. − x1 + x2 ≤ ,
2
7
x1 + x2 ≤ ,
2
x1 , x2 ≥ 0.

Ce nouveau programme a pour solution 23 , 2 . Que nous arrondissions cette solution à (1, 2) ou (2, 2), nous n’ob-

tenons pas de solution réalisable, comme illustré sur la Figure 4.1

Exemple 30. Considérons le programme

max z = x1 + 5x2
s.c. x1 + 10x2 ≤ 20,
x2 ≤ 2,
x1 , x2 ≥ 0 et entiers.

La version relâchée de programme est

max z = x1 + 5x2
s.c. x1 + 10x2 ≤ 20,
x2 ≤ 2,
x1 , x2 ≥ 0,

qui a pour solution optimale (2, 1.8). En arrondissant à (2, 1) afin de garantir l’admissibilité, nous obtenir la valeur
7 pour la fonction objectif, loin de la valeur optimale du programme mixte entier, avec pour valeur optimale 11, en
(0, 2) (voir Figure 4.2).

4.3.2 Approche par énumération

Un modèle en nombres entiers borné (par exemple, un modèle avec uniquement des variables 0–1) possède un
nombre fini de solutions. Nous pourrions envisager de toutes les énumérer, mais le nombre de solutions explose
rapidement. Pour n = 20 variables 0–1, il y a plus d’un million de solutions possibles. Pour n = 30, c’est plus d’un
54 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

x2

3


2, 2
2

1
domaine réalisable du
modèle de PL
0 x1
0 1 2 3

Figure 4.1 – Exemple de relaxation linéaire: problème d’admissibilité des solutions

milliard. Comme il apparaît qu’il est vite déraisonnable de procéder à une énumération complète des solutions, nous
allons essayer de mettre à profit la relaxation en programme linéaire pour éliminer certaines de ces solutions. Cette
technique d’énumération partielle est connue sous le vocable de branch-and-bound (B&B). Il s’agit d’une approche
diviser-pour-régner:
— décomposition du problème en sous-problèmes plus simples ;
— combinaison de la résolution de ces sous-problèmes pour obtenir la solution du problème original.
Dans l’algorithme de branch-and-bound, chaque sous-problème correspond à un sommet dans l’arbre des solutions.
Nous résolvons la relaxation linéaire de chaque sous-problème. L’information tirée de la relaxation linéaire nous
permettra (peut-être) d’éliminer toutes les solutions pouvant être obtenues à partir de ce sommet.

Algorithme de branch & bound: cas 0–1

Un algorithme simple pour énumérer toutes les solutions d’un modèle 0–1 consiste à:
— choisir un sommet dans l’arbre des solutions ;
— choisir une variable x non encore fixée relativement à ce sommet ;
— générer les deux variables x = 0 et x = 1 (la variable x est dite fixée): chaque alternative correspond à un
sommet de l’arbre des solutions ;
— recommencer à partir d’un sommet pour lequel certaines variables ne sont pas encore fixées.
A la racine de l’arbre, aucune variable n’est encore fixée, tandis qu’aux feuilles de l’arbre, toutes les variables ont
été fixées. Le nombre de feuilles est 2n (pour n variables 0–1).
Le calcul de borne (ou évaluation) consiste à résoudre la relaxation linéaire en chaque sommet. L’élagage (ou
élimination) consiste à utiliser l’information tirée de la résolution de la relaxation linéaire pour éliminer toutes les
solutions émanant du sommet courant. Dès lors, le branch-and-bound est un algorithme de séparation et d’évaluation
successives.
4.3. STRATÉGIES DE RÉSOLUTIONS 55

x2

z ∗ = 10 = x1 + 5x2
1

0 x1
0 1 2 3

Figure 4.2 – Exemple de relaxation linéaire: solution médiocre

Exemple 31 (California Mfg). Rappelons le problème

max z = 9x1 + 5x2 + 6x3 + 4x4

s.c. 6x1 + 3x2 + 5x3 + 2x4 ≤ 10;
x3 + x4 ≤ 1;
− x1 + x3 ≤ 0;
− x2 + x4 ≤ 0;
x1 , x2 , x3 , x4 ≤ 1;
x1 , x2 , x3 , x4 ≥ 0;
x1 , x2 , x3 , x4 entiers.

La relaxation linéaire permet aux variables de prendre les valeurs fractionnaires entre 0 et 1, ce qui conduit à la
solution  
5
, 1, 0, 1 ,
6
avec comme valeur z = 16.5. Branchons sur la variable x1 .
Dénotons sous-problème 1 celui obtenu avec x1 = 0:

max z = 5x2 + 6x3 + 4x4

s.c. 3x2 + 5x3 + 2x4 ≤ 10;
x3 + x4 ≤ 1;
x3 ≤ 0;
− x2 + x4 ≤ 0;
x2 , x3 , x4 binaires.

La solution de la relaxation linéaire est (x1 , x2 , x3 , x4 ) = (0, 1, 0, 1), avec z = 9.

56 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

Le sous-problème 2 est celui obtenu avec x1 = 1:

max z = 5x2 + 6x3 + 4x4 + 9

s.c. 3x2 + 5x3 + 2x4 ≤ 4;
x3 + x4 ≤ 1;
x3 ≤ 1;
− x2 + x4 ≤ 0;
x2 , x3 , x4 binaires.

La solution de la relaxation linéaire est alors

 
4 4
(x1 , x2 , x3 , x4 ) = 1, , 0, ,
5 5

avec z = 16 + 15 .

Nous obtenons dès lors les bornes suivantes:

— sous-problème 1: Z1 ≤ 9 ;
— sous-problème 2: Z2 ≤ 16 + 15 .
Notons que toutes les variables sont binaires et tous les paramètres dans l’objectif sont des valeurs entières. Dès
lors, la borne supérieure pour le sous-problème 2 est 16. Pour le sous-problème 1, la solution obtenue est entiere:
c’est la meilleure solution courante. Nous savons que la valeur optimale cherchée, Z, sera au moins

Z ∗ = 9 : Z ≥ Z ∗.

Quels sous-problèmes pouvons-nous à présent considérer afin de nous approcher de la solution optimale ? Tous
les sous-problèmes actuellement traités doivent-ils donner naissance à d’autres problèmes. Si un sous-problème
ne donne lieu à aucun autre problème, nous parlerons d’élagage, en référence avec l’idée de couper la branche
correspondante dans l’arbre d’exploration.
Considérons tout d’abord le sous-probleme 1: la solution optimale de la relaxation PL est entière. Il ne sert
donc à rien de brancher sur les autres variables, puisque toutes les autres solutions entières (avec x1 = 0) sont
nécessairement de valeur inférieures ou égales à 9. Nous pouvons donc élaguer ce sommet.
Pour le sous-probleme 2, la solution optimale de la relaxation PL n’est pas entière:

Z ∗ = 9 ≤ Z ≤ 16.

La branche (x1 = 1) peut encore contenir une solution optimale. Mais si nous avions eu Z2 ≤ Z∗, nous aurions pu
conclure que la branche ne pouvait améliorer la meilleure solution courante.
Un sous-problème est elagué si une des trois conditions suivantes est satisfaite :
— test 1: sa borne supérieure (valeur optimale de la relaxation PL) est inférieure ou égale à Z ∗ (valeur de la
meilleure solution courante) ;
— test 2: sa relaxation PL n’a pas de solution réalisable ;
— test 3: la solution optimale de sa relaxation PL est entière.
Lorsque le test 3 est verifié, nous testons si la valeur optimale de la relaxation PL du sous-problème, Zi , est supérieure
à Z ∗ . Si Zi > Z ∗ , alors Z ∗ := Zi , et nous conservons la solution, qui devient la meilleure solution courante. En
résumé, nous obtenons l’algorithme ci-dessous.

Algorithm 4.1: Branch-and-Bound: cas binaire

1. Initialisation:
(a) poser Z ∗ = −∞ ;
(b) appliquer le calcul de borne et les critères d’élagage à la racine (aucune variable fixée).
2. Critere d’arrêt : s’il n’y a plus de sous-problemes non élagués, arrêter.
3. Branchement:
4.3. STRATÉGIES DE RÉSOLUTIONS 57

(a) parmi les sous-problèmes non encore élagués, choisir celui qui a été crée le plus récemment (s’il y a
égalité, choisir celui de plus grande borne supérieure) ;
(b) appliquer le Test 1: si le sous-problème est élagué, retourner en 2.
(c) brancher sur la prochaine variable non fixée.
4. Calcul de borne:
(a) résoudre la relaxation PL de chaque sous-problème ;
(b) arrondir la valeur optimale si tous les paramètres de l’objectif sont entiers.
5. Elagage: élaguer un sous-problème si
(a) la borne superieure est inférieure ou égale à Z ∗ ;
(b) la relaxation PL n’a pas de solution réalisable ;
(c) la solution optimale de la relaxation PL est entière: si la borne supérieure est strictement supérieure à
Z ∗ , Z ∗ est mise à jour et la solution de la relaxation PL devient la meilleure solution courante.
6. Retourner en 2.

A partir de quel noeud devrions-nous brancher ? Il y a plusieurs choix possibles ; dans cette version, on propose
comme règle de sélection de choisir le sous-problème le plus récemment créé. L’avantage est que cette approche
facilite la réoptimisation lors du calcul de borne, car il n’y a que peu de changements apportés par rapport au
dernier sous-probleme traité. Le désavantage est que cela peut créer un grand nombre de sous-probèmes. Une autre
option est la règle de la meilleure borne: choisir le sous-problème ayant la plus grande borne supérieure.
Dans cette version, la règle de branchement consiste à choisir la prochaine variable non fixée. Il est souvent plus
intéressant de choisir une variable à valeur fractionnaire. En branchant sur une telle variable, il est certain que les
deux sous-problèmes crées mènent à des solutions differentes de la solution courante. De nombreux critères existent
pour choisir une telle variable de facon à orienter la recherche vers un élagage rapide.

Exemple 32 (suite). Jusqu’à maintenant, voici l’arbre obtenu, en branchant sur la variable x1 .

0 F (3)
9 = Z∗

16
1
16

F (3) indique que le sous-probleme a été élagué (fathomed) en raison du Test 3.

Sélection: nous choisissons le sous-problème 2, le seul qui n’a pas encore été élagué. Nous branchons sur la
prochaine variable, soit x2 . Deux nouveaux sous-problèmes sont crées:
— sous-problème 3: x1 = 1, x2 = 0 ;
— sous-problème 4: x1 = 1, x2 = 1.
Considérons tout d’abord le sous-problème 3 (x1 = 1, x2 = 0). Nous obtenons le problème

max z3 = 6x3 + 4x4 + 9

s.c. 5x3 + 2x4 ≤ 4
x3 + x4 ≤ 1
x3 ≤ 1
x4 ≤ 0
x3 , x4 binaire.

La solution de la relaxation PL est  

4
(x1 , x2 , x3 , x4 ) = 1, 0, , 0 .
5
58 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

et
4
Z = 13 + : Z3 ≤ 13.
5
Le sous-problème 4 (x1 = 1, x2 = 1) devient

max z4 = 6x3 + 4x4 + 14

s.c. 5x3 + 2x4 ≤ 1
x3 + x4 ≤ 1
x3 ≤ 1
x4 ≤ 1
x3 , x4 binaires.

La solution de la relaxation PL est  

1
(x1 , x2 , x3 , x4 ) = 1, 1, 0, .
2
et
Z = 16 : Z4 ≤ 16.
Aucun des tests d’élagage ne s’applique sur ces sous-problèmes. Nous devons dès lors choisir un des deux sous-
problèmes pour effectuer un branchement, puisque ce sont ceux crées le plus récemment. Nous choisissons celui de plus
grande borne supérieure, soit le sous-problème 4. Deux possibilités existent: brancher sur la variable fractionnaire,
x4 , ou choisir la prochaine variable non fixée, x3 . Nous explorons la deuxième approche, i.e. nous branchons sur
x3 , et nous générons deux nouveaux sous-problèmes.
Le sous-problème 5, défini avec x1 = 1, x2 = 1, x3 = 0, s’écrit

max z5 = 4x4 + 14
s.c. 2x4 ≤ 1
x4 ≤ 1
x4 ≤ 1
x4 binaire.

La solution de la relaxation PL est  

1
(x1 , x2 , x3 , x4 ) = 1, 1, 0,
2
et
Z = 16 : Z5 ≤ 16.
Le sous-problème 6, défini avec x1 = 1, x2 = 1, x3 = 1, s’écrit

max z5 = 4x4 + 20
s.c. 2x4 ≤ −4
x4 ≤ 0
x4 ≤ 1
x4 binaire.

La relaxation PL n’a pas de solution realisable: ce sous-probleme est élagué.

Le sous-problème 5 ne peut pas être élagué. Il est crée le plus récemment parmi les sous-problèmes non élagués (3
et 5), aussi choisissons-nous le pour effectuer un branchement. Nous branchons sur x4 et générons les sous-problèmes
suivants:
— sous-problème 7: x1 = 1, x2 = 1, x3 = 0, x4 = 0 ;
— sous-problème 8: x1 = 1, x2 = 1, x3 = 0, x4 = 1.
4.3. STRATÉGIES DE RÉSOLUTIONS 59

Toutes les variables sont fixées, aussi pouvons-nous directement résoudre ces sous-problèmes. Le sous-problème 7 a
pour solution (x1 , x2 , x3 , x4 ) = (1, 1, 0, 0), pour Z7 = 14. La solution est entière, aussi nous élaguons le sous-problème
en vertu du Test 3. Puisque Z7 > Z ∗ , Z ∗ = 14 et la solution du sous-problème devient la meilleure solution courante.
Le sous-problème 8 a pour solution (x1 , x2 , x3 , x4 ) = (1, 1, 0, 1). Cette solution n’est pas réalisable, car la première
contrainte (2x4 ≤ 1) est violée. Le sous-problème est par conséquent élagué par le Test 2.
Le sous-problème 3 est le seul non encore élagué. Appliquons le Test 1: Z3 = 13 ≤ 14 = Z∗. Le sous-problème est
donc elagué. Comme il n’y a plus de sous-problèmes non élagués, nous pouvons nous arrêter. La solution optimale
est:
(x1 , x2 , x3 , x4 ) = (1, 1, 0, 0),
et la valeur optimale est Z ∗ = 14.
L’arbre obtenu suite à l’exécution de l’algorithme se présente comme suit:
x1 x2 x3 x4

0 F (3)
9
0 F (1) 0 F (3)
16 13 14
1 0
16 16
1 1 F (2)
16
1 F (2)

F (j): le sous-problème est élagué par le Test j

Algorithme de branch & bound: cas général

Considérons à présent le cas général d’un modèle de programmation (mixte) en nombres entiers: variables
entières générales et variables continues. Comme précédemment, nous ignorons dans un premier temps les contraintes
d’intégralité (les valeurs des variables entières sont traitées comme variables continues), et résolvons le programme
linéaire résultant. Si la solution optimale de ce programme satisfait aux contraintes d’intégralité, alors cette solution
est aussi solution optimale du programme avec variables entières. Sinon, il doit exister au moins une variable xj
dont la valeur α est fractionnaire. La procédure de branchement se généralise alors comme suit: nous séparons
alors le problème relaxé en deux sous-problèmes ; un sous-problème contiendra la contrainte xj ≤ bαc et le second
la contrainte xj ≥ dαe = bαc + 1. Nous répétons le processus pour chacun des sous-problèmes. Cette procédure
est habituellement représentée sous forme d’un arbre binaire où, à chaque niveau, une partition du sommet père
s’effectue suivant la règle décrite précédemment. Il s’agit alors de parcourir cet arbre d’énumération afin d’y trouver
la solution optimale. L’exploration d’un chemin de l’arbre peut prendre fin pour trois raisons:
— la solution devient entière ;
— le domaine admissible d’un sous-problème devient vide ;
— la valeur de l’objectif correspondant à la solution optimale du problème relaxé est inférieure (moins bonne)
à celle d’une solution admissible connue, possiblement obtenue à un autre sommet de l’arbre.
Dans chacun de ces trois cas on dit que le sommet est sondé, et il est inutile de pousser plus loin dans cette direction.
L’algorithme s’arrête lorsque tous les sommets sont sondés. La meilleure solution obtenue au cours du déroulement
de l’algorithme est alors l’optimum global de notre problème.

Algorithm 4.2: Algorithme de B&B: cas général

1. Initialisation:
(a) Poser Z ∗ = −∞.
(b) Appliquer le calcul de borne et les critères d’élagage à la racine (aucune variable fixée).
(c) Critere d’arrêt : s’il n’y a plus de sous-problèmes non élagués, arrêter.
2. Branchement:
60 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

(a) Parmi les sous-problèmes non encore élagués, choisir celui qui a été crée le plus récemment (s’il y a
égalité, choisir celui de plus grande borne supérieure).
(b) Appliquer le Test 1: si le sous-problème est élagué, retourner en 2.
(c) Brancher sur la prochaine variable entière à valeur non entière dans la relaxation PL.
3. Calcul de borne: résoudre la relaxation PL de chaque sous-problème.
4. Elagage: élaguer un sous-problème si
(a) La borne superieure est inférieure ou égale à Z ∗ .
(b) La relaxation PL n’a pas de solution réalisable.
(c) Dans la solution optimale de la relaxation PL, toutes les variables entières sont à valeurs entières: si la
borne supérieure est strictement supérieure à Z ∗ , Z ∗ est mise à jour et la solution de la relaxation PL
devient la meilleure solution courante.
5. Retourner en 2.

4.3.3 Les hyperplans coupants (méthode de coupe)

Vu que l’énumération des variables entières peut être coûteuse, l’idée est d’ajouter des contraintes redondantes
pour le modèle en nombres entiers, mais non pour la relaxation PL.

Exemple 33. Considérons le programme mathématique

5
max 4x1 + x2
x 2
x1 + x2 ≤ 6
9x1 + 5x2 ≤ 45
x1 , x 2 ∈ N .

Le dictionnaire optimal correspondant à la relaxation linéaire de ce programme contient les deux contraintes

15 5 1
x1 = + x3 − x4 ,
4 4 4
9 9 1
x2 = − x3 + x4 ,
4 4 4

où x3 et x4 sont des variables d’écart. Puisque la variable de base x1 n’est pas entière, cette solution de base n’est
pas admissible. nous pouvons réécrire la première contrainte sous la forme

5 1 15
x1 − x3 + x4 = .
4 4 4
En utilisant l’identité
a = bac + (a − bac),
où (a − bac) représente la partie fractionnaire de a (0 ≤ a − bac < 1), nous obtenons
        
5 3 1 1 15 3
x1 + − + x3 + + x4 = + ,
4 4 4 4 4 4

c’est-à-dire, en mettant tous les coefficients entiers à gauche et les coefficients fractionnaires à droite:

3 3 1
x1 − 2x3 − 3 = − x3 − x4 .
4 4 4
4.3. STRATÉGIES DE RÉSOLUTIONS 61

Puisque les variables x3 et x4 sont non négatives, la partie fractionnaire (constante du membre de droite) est
inférieure à 1, le membre de droite est strictement inférieur à 1. Puisque le membre de gauche est entier, le membre
de droite doit aussi être entier. Or un entier inférieur à 1 doit être inférieur ou égal à zéro. Nous en déduisons une
contrainte additionnelle qui doit être satisfaite par toute solution admissible du problème originel, et que ne satisfait
pas la solution de base courante:
3 3 1
− x3 − x4 ≤ 0.
4 4 4
En utilisant les identités x3 = 6 − x1 − x2 et x4 = 45 − 9x1 − 5x2 , nous obtenons la coupe sous sa forme géométrique:
3x1 + 2x2 ≤ 15.
Cette contrainte linéaire rend inadmissible la solution courante admissible, sans éliminer aucune autre solution
entière. Si la solution du nouveau problème est entière, il s’agit de la solution optimale de notre problème. Sinon,
nous construisons une nouvelle coupe et recommençons.
Exemple 34. Reprenons le problème, illustré sur la Figure 4.3,
max z = 3x1 + 2x2
s.c. 2x1 + 3x2 ≤ 4,
x1 , x2 binaire.
Les solutions réalisables sont (0, 0), (1, 0) et (0, 1).

x2

solution optimale
domaine
réalisable

solution optimale entière

0 x1
0 1

Figure 4.3 – Méthode de coupes

Une contrainte redondante est:
x1 + x2 ≤ 1.
Suite à l’ajout de la contrainte redondante, comme représenté sur la Figure 4.4, le problème est résolu à la racine.
Il y a plusieurs algorithmes permettant de générer de telles inégalités redondantes, appelées coupes. Mais il est
rare que leur ajout permette de résoudre le probleme à la racine. L’ajout de coupes permet toutefois de réduire le
nombre de sous-problèmes traités par l’algorithme de Branch-and-Bound. Nous pouvons même ajouter des coupes
pour chaque sous-problème (pas seulement à la racine): nous obtenons alors un algorithme de branch-and-cut.
62 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

x2

domaine
réalisable solution optimale entière

0 x1
0 1

Figure 4.4 – Méthode de coupes: ajout d’une contrainte

4.4 Modélisation avec GAMS

Exemple 35 (Echangeurs de chaleur). Considérons un problème simplifié d’affection de flux de processus à des
échangeurs de chaleur ; les échangeurs de chaleur ont pour but de transférer la chaleur d’un flux chaud à un flux
froid, afin de limiter la consommation d’énergie nécessaire pour chauffer ou refroidir des flux liquides (par exemple
de l’eau) intervenant dans des processus de production industrielle. Le problème d’optimisation correspondant est:
n X
X n
min z = cij xij
i=1 j=1
n
X
xij = 1, j = 1, . . . , n,
i=1
Xn
xij = 1, i = 1, . . . , n,
j=1

xij ∈ {0, 1}, i = 1, . . . , n, j = 1, . . . , n,

ce qui n’est rien d’autre qu’un problème classique d’affection. i représente ici l’indice des n flux, et j, l’indice des
n échangeurs. La variable binaire xij vaut 1 si le flux i est affecté à l’échangeur j, et 0 sinon. Les contraintes
indiquent que chaque échangeur doit être affecté à un flux, et réciproquement. Les coûts cij d’affectation des flux
aux échangeur sont résumés dans la Table 4.5. Le problème peut se formuler en GAMS comme

[Link]

$TITLE Test Problem

$OFFSYMXREF
$OFFSYMLIST
4.4. MODÉLISATION AVEC GAMS 63

XXX
X Flux
XXX 1 2 3 4
Echangeur XXX
X
A 94 1 54 68
B 74 10 88 82
C 73 88 8 76
D 11 74 81 21

Figure 4.5 – coûts d’affectation des échangeurs de chaleur

SETS
I streams / A, B, C, D /
J exchangers / 1*4 / ;

TABLE C(I,J) Cost of assigning stream i to exchanger j

1 2 3 4
A 94 1 54 68
B 74 10 88 82
C 73 88 8 76
D 11 74 81 21 ;

VARIABLES X(I,J), Z;
BINARY VARIABLES X(I,J);

EQUATIONS ASSI(J), ASSJ(I), OBJ;

ASSI(J).. SUM( I, X(I,J) ) =E= 1;

ASSJ(I).. SUM( J, X(I,J) ) =E= 1;
OBJ.. Z =E= SUM ( (I,J), C(I,J)*X(I,J) ) ;

MODEL HEAT / ALL /;

OPTION LIMROW = 0;
OPTION LIMCOL = 0;
OPTION SOLPRINT = OFF;

SOLVE HEAT USING MIP MINIMIZING Z;

DISPLAY X.L, Z.L ;

La commande OPTION SOLPRINT=OFF supprime l’impression de la solution ; la sortie sera contrôlée par la commande
DISPLAY placé juste à la suite.
Exemple 36 (Problème de production). Nous souhaitons maximiser le profit que nous pourrions obtenir en consi-
dérant 3 produits. Il n’y a pas coût de productions unitaires, mais des coûts fixes sont présents pour le premier et le
deuxième produits. Chaque produit est stocké, et nous ne pouvons pas dépasser le volume disponible. Le code GAMS
du problème est
VARIABLES
z profit;
64 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE

POSITIVE VARIABLES
x1 produit1
x2 produit2
x3 produit3;

BINARY VARIABLES
y1
y2 ;

EQUATIONS
Objective Maximize profit
Ressource
Upper1 limit on product 1
Upper2 limit on product 2
Upper3 limit on product 3;

Objective..
z =E= 2*x1 + 3*x2 + 0.8*x3 - 3*y1 - 2*y2;

Ressource..
0.2*x1 + 0.4*x2 + 0.2*x3 =L= 1;

Upper1..
x1 =L= 3*y1;

Upper2..
x2 =L= 2*y2;

Upper3..
x3 =L= 5;

MODEL example /ALL/;

OPTION optcr = 0.05;

SOLVE example USING mip maximizing z;

Notons la ligne

OPTION optcr = 0.05;

qui nous donne le critère d’optimalité sur la valeur optimale. Dans le cas présent, nous exigeons que l’écart relatif
entre la borne supérieur et la borne inférieure soit inférieur à 5%.
4.5. EXERCICES 65

4.5 Exercices
Exercice 1 (Hillier et Lieberman [2], Exercice 12.7-9). Utiliser l’algorithme de Branch & Bound pour résoudre le
problème

max z = 3x1 + 4x2 + 2x3 + x4 + 2x5 ,

s.c. 2x1 − x2 + x3 + x4 + x5 ≤ 3,
− x1 + 3x2 + x3 − x4 − 2x5 ≤ 2,
2x1 + x2 − x3 + x4 + 3x5 ≤ 1,
xj ≥ 0, j = 1, 2, 3, 4, 5,
xj ∈ {0, 1}, j = 1, 2, 3.

4.6 Notes
Le contenu de ce chapitre est repris des notes de cours de Patrice Marcotte (2001) et de Bernard Gendron
(2007). L’exemple GAMS sur les échangeurs de chaleur est dû à Grossmann.
66 CHAPITRE 4. PROGRAMMATION MIXTE ENTIÈRE
Thank goodness we don’t have only serious problems, but ridiculous
ones as well.
Edsger W. Dijkstra

5
Réseaux

5.1 Graphes
Les graphes sont un outil puissant de modélisation. Ils interviennent naturellement dans la représentation de
réseaux de transport ou de télécommunication, par exemple, mais également dans la représentation de structures
relationnelles abstraites. Mathématiquement, un graphe G prend la forme G = (N, A) où N est un ensemble de
sommets et A ⊆ N × N un ensemble d’arcs. Le terme réseau est un terme générique désignant un graphe dont les
sommets ou les arcs possèdent des attributs: coûts, capacités, longueurs, etc.

5.1.1 Graphe orienté

Exemple 37 (Réseau de distribution). Nous reprenons l’exemple 5, illustré sur la Figure 5.1. Les liens entre les
noeuds du graphe ne peuvent être parcourus que dans un sens précis ; nous parlerons de graphe orienté. Les sommets
du graphe sont A, B, C, D, E, et le graphe posséde les arcs (A, B), (A, C), (A, D), (B, C), (C, E), (D, E), (E, D).

A D

B E

Figure 5.1 – Réseau de distribution

67
68 CHAPITRE 5. RÉSEAUX

5.1.2 Graphe non orienté

Exemple 38 (Parc Seervada). Considérons le graphe représenté sur la Figure 5.2, tiré de Hillier et Lieberman [2],
Section 9.1. Nous avons les sommets O, A, B, C, D, E, T , et les arêtes: {O, A}, {O, B}, {O, C}, {A, B}, {A, D},
{B, C}, {B, D}, {B, E}, {D, E}, {D, T }, {E, T }. Le nombre sur chaque arête représente la distance entre les deux
sommets reliés par cette arête

A
2 7
2
5 T
5
O B D
4
4 7
1 1
3
4
C E

Figure 5.2 – Graphe non orienté

5.1.3 Transformations
Un graphe orienté dérivé d’un graphe non orienté en introduisant deux arcs pour chaque arête, un dans chaque
direction. A l’inverse, nous qualifierons de graphe sous-jacent à un graphe orienté le graphe obtenu en enlevant
l’orientation des arcs. Si G est un graphe orienté, le graphe dérivé du graphe sous-jacent à G n’est pas G.

5.1.4 Chemins et circuits

Un chemin [chaîne] est suite d’arcs [d’arêtes] distinct[e]s reliant deux sommets. Il est non orienté s’il est constitué
d’une suite d’arcs distincts qui relient deux sommets, sans considération de l’orientation des arcs. En d’autres mots,
un chemin non orienté est une chaîne dans le graphe sous-jacent. Un circuit [cycle] est un chemin [chaîne] qui
commence et finit au même sommet. Le circuit est non orienté s’il s’agit d’un cycle dans le graphe sous-jacent.

Exemple 39 (Réseau de distribution). Reprenons l’exemple du réseau de distribution, illustré sur la Figure 5.3.
A → C → E → D décrit un chemin, et en ignorant l’orientation des arcs, nous décrivons également un chemin
non orienté. A → D → E → C → B est chemin non orienté, mais pas un chemin. D → E → D est circuit, et en
omettant l’orientation des arcs, un circuit non orienté. A → B → C → A est un circuit non orienté (mais pas un
circuit).

5.1.5 Connexité
Deux sommets sont connexes s’il existe au moins un chemin non orienté les reliant. Un graphe est connexe si
toute paire de sommets est connexe. Le plus petit graphe connexe à n sommets possède n − 1 arêtes ; il est appelé
un arbre. Une définition alternative consiste à dire qu’un arbre est un graphe connexe sans cycle. Un arbre partiel
est un arbre obtenu à partir d’un graphe connexe en incluant tous les sommets.

Exemple 40. Retournons sur l’exemple 38. Le graphe ci-dessous n’est pas un arbre partiel, comme il n’est pas
connexe.
5.2. PROBLÈME DU CHEMIN LE PLUS COURT - ALGORITHMES DE CORRECTIONS D’ÉTIQUETTES69

A D

B E

Figure 5.3 – Réseau de distribution

T
O B D

C E

Le graphe ci-dessous est connexe, mais présente des cycles. Par conséquent, il ne s’agit pas d’un arbre.

T
O B D

C E

Le graphe ci-dessous est un arbre partiel.

T
O B D

C E

5.2 Problème du chemin le plus court - algorithmes de corrections d’éti-

quettes
Considérons un graphe non orienté et connexe, sur lequel deux sommets ont un rôle particulier:
— source (ou origine): O ;
— puits (ou destination): T .
70 CHAPITRE 5. RÉSEAUX

A chaque arête {i, j}, nous associons une distance cij ≥ 0. Nous cherchons le chemin non orienté (ou chaîne) le plus
court reliant O à T . Par chemin le plus court, nous désignons celui dont la distance totale (somme des distances des
arêtes du chemin) est minimale parmi tous les chemins de O à T . Afin de procéder, nous affectons à chaque noeud
du graphe une étiquette, représentant la meilleure distance connue à un moment donné, de l’origine à ce noeud.

5.2.1 Algorithme de Dijkstra

Il s’agit d’une méthode itérative A chaque itération, nous choisissons le sommet j le plus près de O et nous
fixons d(j), la variable calculant la distance entre 0 et j (le sommet j est dit marqué). Au départ, seul O est marqué
et d(O) = 0. Le sommet le plus près est choisi parmi les sommets non marqués reliés à au moins un sommet marqué.
Le sommet choisi j est celui qui atteint
( )
min min d(i) + cik .
sommets k non marqués sommets i marqués

d(j) est fixée à cette valeur. Nous nous arrêtons lorsque T est marqué si nous cherchons à connaître le chemin allant
de O à T spécifiquement, ou sinon jusqu’à avoir marqué tous les sommets. Par la suite, nous supposerons dans la
description que nous voulons marquer tous les sommets. Plus spécifiquement, nous avons l’algorithme ci-dessous.
L’algorithme de Dijkstra s’applique aussi sur un graphe orienté.

Algorithm 5.1: Algorithme de Dijkstra

Soient G = (N, A) un graphe connexe, de longueur d’arêtes non-négative (la longueur des arêtes est donnée
par une fonction ` : A → R+ ), une source O ∈ N . Désignons par S l’ensemble des sommets marqués.
1. Initialisation. Pour tout u ∈ N , faire
— dist(u) := ∞ ;
— pred(u) := NULL.
où pred(u) indique le prédécesseur de u sur le plus court chemin de O à u. Le noeud origine prend un rôle
particulier, et comme il est évident que le plus court chemin de O à O a une longueur nulle, nous posons
dist(O) = 0. Comme à ce stade, aucun noeud n’a encore été marqué, nous posons S = ∅.
2. Tant que S 6= N , prenons u tel que

dist(u) ≤ dist(v), ∀v ∈ N \S.

Posons S := S ∪ {u}, et, pour tout arc (u, v) existant, si

dist(v) > dist(u) + `(u, v),

posons

dist(v) := dist(u) + `(u, v),

pred(v) = u.

Il est possible de montrer qu’une fois un noeud u marqué, dist(u) donne la longueur du plus court chemin de O
à u.

Exemple 41. L’application de l’algorithme de Dijkstra pour l’exemple 38, repris à la Figure 5.4, donne la série
d’itérations résumée dans la Table 5.1.

Il est également possible de trouver les chemins les plus courts entre la source et tous les autres sommets (il
suffit de regarder la valeur de d(i) en chaque sommet i). Nous pouvons également trouver les chemins les plus courts
entre toutes les paires de sommets, avec n applications de l’algorithme de Dijkstra (mais il est possible de faire
mieux). Si certaines distances sont négatives, l’algorithme de Dijkstra ne s’applique cependant pas.
5.2. PROBLÈME DU CHEMIN LE PLUS COURT - ALGORITHMES DE CORRECTIONS D’ÉTIQUETTES71

A
2 7
2
5 T
5
O B D
4
4 7
1 1
3
4
C E

Figure 5.4 – Plus court chemin

Itération Sommet S N \S Voisins Distance Prédécesseur

marqué non marqués
0 - ∅ N - - -
1 O {O} {A, B, C, D, E, T } A 2 O
B 5 O
C 4 O
2 A {O, A} {B, C, D, E, T } B 4 A
D 9 A
3 B {O, A, B} {C, D, E, T } D 8 B
C 4 O
E 7 B
4 C {O, A, B, C} {D, E, T } E 7 B
5 E {O, A, B, C, E} {D, T } D 8 B
T 14 E
6 D {O, A, B, C, D, E} {T } T 13 D
7 T N ∅ - - -

Table 5.1 – Exécution de l’algorithme de Dijkstra

72 CHAPITRE 5. RÉSEAUX

5.3 Flot dans un réseau

5.3.1 Réseau
Un réseau est un graphe orienté ayant
— des capacités sur les arcs ;
— des sommets d’offre (ou sources) ;
— des sommets de demande (ou puits) ;
— des sommets de transfert.
Le flot dans un réseau est le nombre d’unités circulant sur les arcs du réseau de façon à respecter les capacités et
les contraintes de conservation de flot. En chaque sommet, le flot sortant moins le flot entrant vaut
— l’offre (si le sommet est une source) ;
— la demande (si le sommet est un puits) ;
— 0 (en un sommet de transfert).
Soit A est l’ensemble des arcs du réseau. En dénotant par xij la quantité de flot qui passe sur l’arc (i, j) ∈ A, nous
avons X X
xij − xji = bi , i, j ∈ N
(i,j)∈A+ (i) (j,i)∈A− (i)

où
— bi = 0 (transfert), offre (source), -demande (puits) ;
— N désigne l’ensemble des sommets ;
— A+ (i) est l’ensemble des arcs sortants du sommet i ;
— A− (i) désigne l’ensemble des arcs entrants au sommet i ;
avec 0 ≤ xij ≤ uij .

5.3.2 Modèle de flot

Le chemin le plus court peut être vu comme le flot dans un réseau, lequel est le graphe orienté dérivé. On enlève
les arcs entrant à O et les arcs émanant de T . O est la seule source, avec une offre égale à 1. T est le seul puits,
avec une demande valant 1. Le flot sur chaque arc (i, j) est soit 1, si l’arc appartient au chemin le plus court, soit
0, sinon.

5.3.3 Algorithme de Bellman-Ford

La déficience majeure de l’algorithme de Dijkstra est qu’en présence d’arcs (arêtes) de poids négatif, ajouter un
arc à un chemin peut faire décroître la longueur de ce chemin, de sorte que lors du marquage d’un noeud u, dist(u)
ne correspond plus nécessairement à la longueur du plus court chemin de O à u. En particulier, si considérons le
plus court chemin
O → v1 → v2 → . . . → vk
de O à vk , alors il n’est plus vrai que
d(O, vi ) ≤ d(O, vk ) ∀ i < k.
Il est pourtant fréquent d’avoir des poids négatifs sur certains arcs. L’algorithme de Bellman-Ford permet de traiter
de tels cas, pour autant qu’il n’y ait pas de cycles de longueur négative, c’est-à-dire un cycle dont les sommes des
poids des arêtes est strictement négatif. En effet, s’il est possible d’atteindre de O un cycle C de longueur négative,
et qu’il est possible d’atteindre T depuis C, alors la longueur du plus court chemin de O à T vaut −∞.
Considérons un graphe G = (N, A) connexe, sans cycle négatifs. Il existe alors un chemin

O → v1 → v2 → . . . → vk = T.

Nous allons garder trace des distances à partir de O, en construisant pour chaque noeud des surestimations succes-
sives de la longueur du plus court chemin.

Algorithm 5.2: Algorithme de Bellman-Ford

Soient un graphe G = (N, A) connexe, sans cycle négatifs, un sommet source O et un sommet destination T .
5.3. FLOT DANS UN RÉSEAU 73

1. Initialisation. Pour tout u ∈ N , posons dist(u) := ∞ et pred(u) = NULL. Posons également dist(O) = 0.
Poser i = 0.
2. Si i = |N | − 1, arrêt. Sinon, pour chaque arête {u, v}, posons

dist(v) := min{dist(v), dist(u) + `(u, v)}.

Si dist(v) a été modifié, poser pred(v) = u. Mettre à jour i:

i := i + 1.

Retour en 2.

Considérons le plus court chemin O à T :

O → v1 → v2 → . . . → vk = T.

Pour chaque i ≤ k,
O → v1 → v2 → . . . → vi .
est le plus court chemin de O à vi . En effet
— après l’itération 1, dist(v1 ) = d(O, v1 ),
— après l’itération 2, dist(v2 ) = d(O, v2 ),
et ainsi de suite jusqu’à l’itération i.

5.3.4 Modèle du chemin critique (PERT/CPM)

La théorie des graphes est utilisée également en ordonnancement de tâches, en particulier sur base d’un dia-
gramme de PERT 1 . Considérons un projet divisé en tâches diverses. Un graphe orienté connexe (acyclique) repré-
sente les relations de précédence entre les tâches, avec pour chaque étape du projet un noeud i correspond. L’arc
(i, j) représente une tâche à effectuer, et nous lui associons une pondération tij , représentant la durée de la tâche.
S’il existe également un arc (j, k), cela signifie que la tâche (i, j) doit se terminer avant que (j, k) ne débute. Les
noeuds 1 et N correspondent au début et à la fin du projet. Un chemin dans le graphe représente une séquence de
travaux qui doivent être exécutés dans un certain ordre. Nous souhaitons connaître le temps requis pour terminer
le projet (i.e. avoir effectué toutes les tâches, certaines pouvant s’accomplir en parallèle), ainsi que les activités
critiques, c’est-à-dire celles qui déterminent le temps total du projet.
Un plus long chemin dans le réseau s’appelle un chemin critique et sa longueur correspond à la durée minimale
du projet. Le poids d’un chemin critique est une borne inférieure du temps total nécessaire pour mener à bien le
projet. Il est facile de calculer un tel chemin sur base des algorithmes de plus court chemin afin de trouver un tel
chemin critique. Une première approche consiste à considérer l’opposé de chaque pondération tij , ce qui conduit à
un graphe avec des poids négatifs. Dans ce cas, l’algorithme de Bellman-Ford s’appliquera, pour autant qu’il n’y
ait pas de cycles 2 .

Exemple 42. Nous souhaitons mener à bien le développement d’un nouveau produit, depuis sa production jusqu’au
rapport de projet menant à sa vente. Nous avons la liste des activités reprises dans la Table 5.2, Dénotons l’instant
de départ par O. Nous pouvons représenter les différences tâches, et leurs dépendances, suivant le graphe de la
Figure 5.5. Il est facile de voir que le chemin critique, de longueur 13, vaut

O → A → D → G → J.

Autrement dit, le projet prendra au minimum 13 mois, et les activités critiques sont, dans l’ordre, la conception
du produit, la mise au point du prototype, son test et le rapport de projet. Les autres tâches peuvent s’effectuer en
parallèle (en respectant les contraintes de précédence). Le noeud O ne correspond à aucune activité, mais permet
de définir un point d’entrée dans le graphe. Le point de sortie peut être identifié à J car c’est l’unique activité
1. PERT est un acronyme pour Program Evaluation and Review Technique
2. Pourquoi ?
74 CHAPITRE 5. RÉSEAUX

terminale. Si plusieurs activités ont lieu de manières concomitantes en fin de projet, il n’y a pas d’activité terminale
que nous pourrions identifier. Dans ce cas, nous ajoutons un noeud artificiel, disons T , et traçons un arc depuis
chaque activité de fin vers ce noeud terminal T . Chaque arc portera le pondération à 0, ainsi nous ne modifions pas
la durée du projet de manière arbitraire.

Activité Description Tâches préalables Durée (mois)

A Conception du produit - 5
B Recherche de marché - 1
C Analyse de production A 2
D Prototype A 3
E Brochure de vente A 2
F Anayse de coût C 3
G Test de produit D 4
H Entraînement de vente B, E 2
I Tarification H 1
J Rapport de projet F, G, I 1

Table 5.2 – Tâches d’un projet commercial

C
2 3

3
A D F
5 2 4 1

1
0 E G J
1 2 1

2 2
B H I

Figure 5.5 – Graphe du projet commercial

5.4 Problème de l’arbre partiel minimum

Considérons un graphe non orienté et connexe. A chaque arête {i, j}, nous associons une distance cij ≥ 0. Nous
cherchons à construire un arbre partiel (plus petit graphe connexe contenant tous les sommets) dont la somme
des distances soit minimum parmi tous les arbres partiels du graphe. C’est un exemple simple de problème de
conception de réseaux (network design): choisir une configuration de réseau (sous-ensemble d’arcs) qui optimise un
certain critère.
5.5. PROBLÈME DU FLOT MAXIMUM 75

5.4.1 Algorithme de Prim (1957)

Il s’agit à nouveau d’une approche itérative. L’algorithme est initialisé en choisissant un sommet i (arbitrai-
rement) et en le reliant au sommet j le plus proche. En d’autres termes, nous ajoutons l’arête {i, j}. À chaque
itération, nous choisissons le sommet non relié j le plus près d’un des sommets déjà reliés i et nous ajoutons {i, j}.
Nous nous arrêtons lorsque tous les sommets ont été reliés. En cas d’égalité, nous pouvons choisir arbitrairement.
De telles égalités indiquent qu’il pourrait y avoir plusieurs solutions optimales.

Exemple 43. Considérons à nouveau l’exemple 38.

1. L’initialisation consiste à choisir le sommet O et à le relier au sommet le plus près: on ajoute {O, A}.
2. Le sommet non relié le plus près de O ou de A est B ; comme il est plus près de A, on ajoute {A, B}.
3. Le sommet non relié le plus près de O, de A ou de B est C ; puisqu’il est plus près de B, on ajoute {B, C}.
4. Le sommet non relié le plus près d’un des sommets reliés est E ; nous ajoutons l’arête {B, E},
5. Le sommet non relié le plus près d’un des sommets reliés (E) est D ; nous ajoutons l’arête {E, D}
6. Le sommet non relié le plus près d’un des sommets reliés (D) est T ; nous ajoutons l’arête {D, T }. Tous
les sommets sont à présent reliés, aussi, nous nous arrêtons. La valeur optimale correspond à la somme des
distances des arêtes ajoutées, soit 14.
L’arbre partiel minimum est représenté en rouge dans la Figure 5.6.

A
2 7
2
5 T
5
O B D
4
4 7
1 1
3
4
C E

Figure 5.6 – Arbre partiel minimum

5.5 Problème du flot maximum

Considérons un graphe orienté et connexe. À chaque arc (i, j), nous associons une capacité uij > 0. Il y a deux
sommets spéciaux:
— source (ou origine) O ;
— puits (ou destination) T .
Tous les autres sont des sommets de transfert. L’offre en O et la demande en T sont variables. L’offre en O doit
cependant correspondre à la demande en T , et indique la valeur du flot entre O et T . Nous cherche à maximiser la
valeur du flot entre O et T .
Supposons que nous avons déjà affecté un flot sur les arcs. La capacité résiduelle d’un arc (i, j) est comme

uij − xij

Nous définission le graphe résiduel comme le graphe non orienté sous-jacent, où, sur chaque arête, nous associons
deux valeurs: la capacité résiduelle et le flot déjà affecté.
76 CHAPITRE 5. RÉSEAUX

Exemple 44 (Parc Seervada). Continuons avec l’exemple 38, mais en considérant à présent un graphe orienté,
comme sur la Figure 5.7. Supposons qu’en période de grande affluence, nous disposions d’une flotte d’autobus pour
faire visiter les différents postes d’observation du parc. La réglementation limite le nombre d’autobus pouvant circuler
sur chaque tronçon de route. Comment faire circuler les autobus dans le parc de façon à maximiser le nombre total
d’autobus allant de l’origine (O) à la destination (T ) ?

A
2 7
2
5 T
5
O B D
4
4 7
1 1
3
4
C E

Figure 5.7 – Parc Seervada

Supposons par exemple que l’arc (0, B) a une capacité de 7, ce que nous représentons comme

7
O B

Nous pouvons représenter l’arc comme une arête, en indiquant la capacité disponible au niveau du noeud origine,
et le flot actuel, 0, au niveau de la destination.

O 7 0 B

Si nous affectons cinq unités de flot, nous pouvons représenter la capacité résiduelle au niveau de l’origine, et le
flot au niveau de la destination, comme suit.

O 2 5 B

Interprétation du graphe résiduel: nous avons affecté 5 unités de flot sur l’arc (O, B). Si nous traversons O → B,
2 est la capacité résiduelle et le flot sur (O, B) vaut [Link] nous traverse B → O, la capacité résiduelle est 5 et le flot
sur (B, O) vaut 2.

Un chemin d’augmentation est un chemin allant de la source au puits dans le graphe orienté dérivé du graphe
résiduel. Pour chaque arête {i, j}, l’arc (i, j) possède une capacité résiduelle uij − xij , et l’arc (j, i) possède une
capacité résiduelle xij . Chaque arc du chemin possède une capacité résiduelle strictement positive. La capacité
résiduelle d’un chemin d’augmentation est le minimum des capacités résiduelles de tous les arcs du chemin.

5.5.1 Algorithme de Ford-Fulkerson

Algorithm 5.3: Ford-Fulkerson

1. Initialiser le flot: 0 unité sur chaque arc.
2. Si aucun chemin d’augmentation ne peut être identifié, arrêter: le flot est maximum.
3. Identifier un chemin d’augmentation P ; soit c sa capacité résiduelle.
5.5. PROBLÈME DU FLOT MAXIMUM 77

4. Sur chaque arc de P

(a) augmenter le flot de c ;
(b) diminuer la capacité résiduelle de c.
5. Retourner à l’étape 2.

Algorithm 5.4: Identifier un chemin d’augmentation

1. Marquer la source O (aucun autre sommet n’est marqué) ; tous les sommets sont non visités.
2. S’il n’y a aucun sommet marqué non visité, arrêter: il n’existe aucun chemin d’augmentation.
3. Choisir un sommet marqué non visité i.
4. Visiter i: pour chaque (i, j) de capacité résiduelle strictement positive dans le graphe orienté dérivé du
graphe résiduel, marquer j.
5. Si T est marqué, arrêter: un chemin d’augmentation a été identifié.
6. Retour à l’étape 2.

Exemple 45 (Parc Seervada: suite).

0. Le graphe résiduel initial est représenté sur la Figure 5.8.

A 3
0
1

5 0 0 T
0 9
7 0 4
O B D 0
0
4 2 5 0

0 0 1
6
0 C E
4 0

Figure 5.8 – Parc Seervada: initialisation.

1. Nous pouvons identifier le chemin d’augmentation O → B → E → T . La capacité résiduelle vaut alors

min{7, 5, 6} = 5. En suivant l’algorithme, nous augmentons le flot et diminuons la capacité résiduelle de 5
unités sur tous les arcs de O → B → E → T , ce qui conduit au graphe de la Figure 5.9.
2. Nous pouvons à présent identifier le chemin d’augmentation O → A → D → T , qui donne une capacité
résiduelle de min{5, 3, 9} = 3. Nous augmentons le flot et diminuons la capacité résiduelle de 3 unités sur
tous les arcs de O → A → D → T , comme représenté sur la Figure 5.10.
3. Nous avons le chemin d’augmentation O → A → B → D → T , et la capacité résiduelle min{2, 1, 4, 6} = 1.
Augmentons le flot et diminuons la capacité résiduelle de 1 unité sur tous les arcs de O → A → B → D → T ,
pour obtenir la Figure 5.11.
4. Identifions le chemin d’augmentation O → B → D → T , ce qui donne comme apacité résiduelle min{2, 3, 5} =
2. Nous augmentons le flot et diminuer la capacité résiduelle de 2 unités sur tous les arcs de O → B → D → T ,
ce qui conduite à la Figure 5.12.
5. Prenons le chemin d’augmentation O → C → E → D → T , avec comme capacité résiduelle min{4, 4, 1, 3} =
1. Ceci nous conduit à augmenter le flot et diminuer la capacité résiduelle de 1 unité sur tous les arcs de
O → C → E → D → T . Nous obtenons ainsi le graphe de la Figure 5.13.
78 CHAPITRE 5. RÉSEAUX

A 3
0
1

5 0 0 T 5
0 9
2 5 4
5 O B D 5
0
4 2 0 0

0 5 1
1
0 C E
4 0

Figure 5.9 – Parc Seervada: itération 1.

A 0
3
1

2 0 3 T 8
3 6
2 5 4
8 O B D 5
0
4 2 0 0

0 5 1
1
0 C E
4 0

Figure 5.10 – Parc Seervada: itération 2.

A 0
4
0

1 1 4 T 9
3 5
2 5 3
9 O B D 5
1
4 2 0 0

0 5 1
1
0 C E
4 0

Figure 5.11 – Parc Seervada: itération 3.

5.5. PROBLÈME DU FLOT MAXIMUM 79

A 0
4
0

1 1 6 T 11
3 3
0 7 1
11 O B D 5
3
4 2 0 0

0 5 1
1
0 C E
4 0

Figure 5.12 – Parc Seervada: itération 4.

A 0
4
0

1 1 7 T 12
3 2
0 7 1
12 O B D 5
3
3 2 0 1

0 5 0
1
1 C E
3 1

Figure 5.13 – Parc Seervada: itération 5.

80 CHAPITRE 5. RÉSEAUX

6. Nous pouvons choisir le chemin d’augmentation O → C → E → T . La capacité résiduelle est min{3, 3, 1} = 1.

En augmentant le flot et en diminuant la capacité résiduelle de 1 unité sur tous les arcs de O → C → E → T ,
nous obtenons la Figure 5.14.

A 0
4
0

1 1 7 T 13
3 2
0 7 1
13 O B D 6
3
2 2 0 1

0 5 0
0
2 C E
2 2

Figure 5.14 – Parc Seervada: itération 6.

7. Considérons le chemin d’augmentation: O → C → E → B → D → T , de capacité résiduelle min{2, 2, 5, 1, 2} =

1. Par conséquent, nous augmenter le flot et diminuons la capacité résiduelle de 1 unité sur tous les arcs de
O → C → E → B → D → T , ce qui donne la Figure 5.15.

A 0
4
0

1 1 8 T 14
3 1
0 7 0
14 O B D 6
4
1 2 1 1

0 4 0
0
3 C E
1 3

Figure 5.15 – Parc Seervada: itération 7.

8. Il n’y a plus aucun chemin d’augmentation possible ; nous avons atteint le flot maximum, lequel est réparti
comme représenté sur la Figure 5.16.

5.5.2 Flot maximum - Coupe minimum

Supposons que nous partitionnions l’ensemble des sommets en deux sous-ensembles X, Y . Une coupe est un
sous-ensemble d’arcs allant d’un sommet de X vers un sommet de Y . La capacité d’une coupe est la somme des
capacités des arcs de la coupe. Une coupe minimum est une coupe dont la capacité est minimum parmi toutes les
coupes possibles. Le théorème flot max–coupe min dit que la valeur du flot maximum est égale à la capacité d’une
coupe minimum. Les sommets marqués lors de la dernière itération de l’algorithme de Ford-Fulkerson définissent la
coupe minimum, comme illustré sur la Figure 5.17 pour l’exemple du parc Seervada.
5.5. PROBLÈME DU FLOT MAXIMUM 81

A
4 3
1
8 T 14
7
14 O B D
4
3 6
1
4
3
C E

Figure 5.16 – Parc Seervada: fin.

A 3
0
1

5 0 0 T
0 9
7 0 4
O B D 0
0
4 2 5 0

0 0 1
6
0 C E
4 0

Figure 5.17 – Parc Seervada: coupe minimum.

82 CHAPITRE 5. RÉSEAUX

5.6 Problème de flot minimum

Plutôt que de maximiser le flot sur le réseau, nous pourrions chercher à trouver la façon la plus économique
d’envoyer un certain flot à travers un réseau. Considérons un réseau de flot G(N, A), avec une source s ∈ N et un
puits t ∈ N . On veut envoyer le flux d de s à t. Soient `(i, j) le coût unitaire sur l’arête {i, i}, et f (i, j) le flot {i, j}.
Nous avons le problème linéaire
X
min `(i, j)f (i, j),
i,j∈N

f (i, j) ≤ cij ,
f (i, j) = −f (j, i),
X
f (i, k) = 0, ∀k 6= s, t,
k∈N
X X
f (s, k) = f (k, t) = d.
k∈N k∈N

5.7 Exercices
1. Formulez le problème de plus court chemin comme un problème de programmation linéaire.
Astuce: on peut voir la longueur d’un chemin comme le coût total payé par une unité de flot voyageant de
l’origine à la destination, le coût d’un arc étant sa longueur.
2. Considérons la construction d’une maison, et nous voulons connaître la durée des travaux préléminaires. Cinq
tâches sont considérées.
Tâche Description Durée (jours)
1 Exécution des terrassements 10
2 Mise en place de la grue 2
3 Fondations 5
4 Branchement électrique 3
5 Installation de la fosse septique 6
Nous avons les contraintes suivantes:
(a) la grue ne peut fonctionner que si le branchement électrique est effectué, mais sa mise en place peut se
faire en parallèle du branchement électrique ;
(b) nous avons besoin de la grue pour les fondations ;
(c) l’installation de fosse septique et les fondations ne peuvent être exécutés que si les travaux de terrasse-
ment sont terminés.
Quelle est la durée minimale totale des travaux ? Soyez précis dans le développement de votre réponse.

5.8 Notes
Ce chapitre se base partiellement sur les notes de cours de Bernard Gendron, 2007. La section sur les plus cours
chemins repose également sur du matériel de cours mis en ligne par Otfried Cheong.
 Modèles stochastiques
6
6.1 Concepts de base
Une expérience aléatoire est une expérience dont le résultat n’est pas connu avec certitude. Supposons que tous
les résultats possibles de cette expérience soient connus. L’espace d’échantillonnage Ω est cet ensemble de résultats
possibles d’une expérience aléatoire.
Exemple 46 (Espace d’échantillonnage).

— Tirage d’une pièce de monnaie: Ω = {P, F }.

— Lancement d’un dé: Ω = {1, 2, 3, 4, 5, 6.
— Temps écoulé avant l’arrivée d’un premier client dans un magasin ouvert durant 8 heures, depuis sont
ouverture: Ω = [0, 8].
Un événement E est un sous-ensemble de l’espace d’échantillonnage. Supposons que nous répétions l’expérience
aléatoire un grand nombre de fois (n), et que l’événement E se produise m fois. La probabilité associée à l’événement
E peut être approchée comme
m
P [E] ≈ .
n
Une définition empirique de la probabilité de E est (en version fréquentiste)
m
P [E] = lim ,
n→∞ n
De manière plus formelle
Definition 1. P est une mesure de probabilité si
— 0 ≤ P [E] ≤ 1, pour tout E ⊂ Ω ;
— P [∅] = 0 et P [Ω] = 1 ;
— P [E1 ∪ E2 ] = P [E1 ] + P [E2 ], si E1 et E2 sont disjoints (i.e. E1 ∩ E2 = ∅).
Exemple 47. Tirage d’une pièce de monnaie: P [{P }] = P [{F }] = 0.5.
L’occurence un événement E1 se produit peut influencer la probabilité d’un autre événement E2 . Par exemple,
la probabilité qu’il pleuve demain (E2 ) est plus élevée s’il pleut aujourd’hui (E1 ) que s’il ne pleut pas. Si P [E1 ] > 0,
nous définissons la probabilité conditionnelle associée à l’événement E2 , étant donné E1 , comme:
P [E1 ∩ E2 ]
P [E2 | E1 ] = .
P [E1 ]

83
84 CHAPITRE 6. MODÈLES STOCHASTIQUES

La probabilité conditionnelle jouit des propriétés suivantes:

— 0 ≤ P [E2 | E1 ] ≤ 1 ;
— P [∅ | E1 ] = 0 et P [Ω | E1 ] = 1 ;
— P [E2 U E3 | E1 ) = P [E2 | E1 ] + P [E3 | E1 ], si E2 et E3 sont disjoints.
Deux événements E1 et E2 sont indépendants si

P [E2 | E1 ] = P [E2 ].

De manière alternative, nous pouvons utiliser les définitions suivantes:

P [E1 | E2 ] = P [E1 ]
P [E1 ∩ E2 ] = P [E1 ]P [E2 ].

En général, nous postulons l’indépendance de deux événements pour se servir des définitions ci-dessus, plutôt que de
déduire l’indépendance de deux événements à partir des définitions. K événements E1 , E2 ,. . . , Ek sont indépendants
si
P [E1 ∩ E2 ∩ . . . Ek ] = P [E1 ]P [E2 ] . . . P [Ek ].

6.2 Variable aléatoire

Une variable aléatoire X est une fonction qui associe une valeur numérique X(s) à chaque élément s de l’espace
d’échantillonnage:
X : Ω → R.
Il existe deux types principaux de variables aléatoires:
— variable aléatoire continue: valeurs réelles ;
— variable aléatoire discrète: valeurs entières ou nombre fini de valeurs.
Exemple 48 (Lancer de dés). Considérons l’expérience aléatoire du lancer de deux dés. L’espace d’échantillonnage
est Ω = {(1, 1), (1, 2), . . . , (6, 6)}. Une variable aléatoire X serait la somme des résultats des deux dés. Ainsi, par
exemple,
1
P [X = 2] = P [s ∈ Ω : X(s) = 2] = P [(1, 1)] = ;
36
6 1
P [X ≤ 4] = P [s ∈ Ω : X(s) ≤ 4] = P [(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (3, 1)] = = .
36 6
La fonction de répartition associée à une variable aléatoire X est définie comme

FX (b) = P [X ≤ b] = P [s ∈ Ω | X(s) ≤ b],

Elle a comme propriétés:

— non décroissante ;
— limb→−∞ FX (b) = 0 et limb→∞ FX (b) = 1 ;
— P [a < X ≤ b] = FX (b) − FX (a), vu que

{s ∈ Ω : X(s) ≤ b} = {s ∈ Ω : X(s) ≤ a} ∪ {s ∈ Ω : a < X(s) ≤ b}.

Exemple 49 (Lancer de dés). Soit X la somme des résultats des deux dés. Nous avons

FX (1) = P [X ≤ 1] = 0;
1
FX (2) = P [X ≤ 2] = P [X = 2] = ;
36
6 1
FX (4) = P [X ≤ 4] = = ;
36 6
FX (12) = P [X ≤ 12] = 1.
6.2. VARIABLE ALÉATOIRE 85

6.2.1 Variables aléatoires discrètes

La fonction de masse associée à une variable aléatoire X est définie comme

PX (k) = P [X = k] = P [s ∈ Ω : X(s) = k].

Pour une variable aléatoire discrète, nous avons

X X
FX (b) = P [X ≤ b] = P [X = k] = PX (k).
k≤b X≤b

aussi,
P [a < X < b] = FX (b) − FX (a) − PX (b).

Exemple 50 (Lancer de dés). Soit X la somme des résultats des deux dés. Alors,
1
FX (2) = P [X ≤ 2] = P [X = 2] = PX (2) = ;
36
6 1
FX (4) = P [X ≤ 4] = P [X = 2] + P [X = 3] + P [X = 4] = PX (2) + PX (3) + PX (4) = = .
36 6

6.2.2 Variables aléatoires continues

Une variable aléatoire X est continue si sa fonction de répartition peut être représentée ainsi:
Z b
FX (b) = fx (x)dx.
−∞

La fonction sous l’intégrale est appelée fonction de densité et satisfait les conditions suivantes:

fx (x) ≥ 0, ∀x;
Z ∞
fx (x)dx = 1
−∞

Si la fonction de densité est continue, alors elle est égale à la dérivée de la fonction de répartition:
dFX
fX (x) = (x).
dx
La fonction de masse prend toujours la valeur 0:

PX (x) = 0, ∀x.

Pour tout intervalle de la forme < a, b >, nous avons également

Z b
P [X ∈< a, b >] = f (x)dx = FX (b) − FX (a).
a

6.2.3 Espérance mathématique

L’espérance mathématique (moyenne) associée à une variable aléatoire X est dénotée E[X] et définie comme
suit: (P P
k kPX (k) = k kP [X = k] si X est discrète,
E[X] = R ∞
−∞
xfX (x)dx si X est continue.

Exemple 51 (Lancer de deux dés). Soit X la somme des résultats des deux dés.
X X
E[X] = kP [X = k] = kP [X = k] = 7.
k k=2,...,12
86 CHAPITRE 6. MODÈLES STOCHASTIQUES

Nous pouvons également considérer l’espérance d’une fonction g(X). Si X est une variable aléatoire discrète,
cela donne X
E[g(X)] = g(k)PX (k).
k
Pour une variable continue, nous aurons
Z +∞
E[g(X)] = g(X)fX (x)dx.
−∞

En particulier, nous avons, si a ∈ R et X, Y sont deux variables aléatoires,

E[aX] = aE[X],
E[X + Y ] = E[X] + E[Y ].

6.2.4 Variance
La variance associée à une variable aléatoire X, dénotée σ 2 (X), est définie par la formule
σ 2 (X) = E[X 2 ] − E[X]2 = E[(X − E[X])2 ].
Exemple 52 (Lancer de deux dés). Soit X la somme des résultats des deux dés.
X
σ 2 (X) = E[X 2 ] − E[X]2 = k 2 P [X = k] − E[X]2 ≈ 5.833.
k

6.3 Loi de probabilité

Une loi de probabilité est un modèle d’une expérience aléatoire. Elle est habituellement représentée par la
fonction de répartition d’une variable aléatoire. Si cette dernière est discrète, la loi de probabilité est dite discrète.
Une loi de probabilité discrète peut être représentée par sa fonction de masse. Si la variable aléatoire est continue,
la loi de probabilité est dite continue, et peut être représentée par sa fonction de densité.

6.3.1 Loi de Bernouilli

L’espace d’échantillonnage se limite à deux éléments, dénotés par exemple par Ω = {S, E}. Nous définissons la
variable aléatoire X comme suit:
X(S) = 1 et X(E) = 0.
La fonction de masse est
PX (1) = p et PX (0) = 1 − p,
où p est un paramètre. De manière équivalente, nous pouvons l’écrire comme
PX (x) = px (1 − p)1−x .
Nous avons de plus
E[X] = p et σ 2 (X) = p(1 − p).
Par exemple, le tirage d’une pièce de monnaie suit une loi de Bernouilli, avec p=1/2.

6.3.2 Loi uniforme

Une variable aléatoire continue X (qui prend ses valeurs dans l’intervalle [a, b]) suit une loi uniforme (notée
U [a, b]) si sa fonction de densité est:
1
fX (x) = , ∀x ∈ [a, b].
b−a
Celle loi modélise l’expérience aléatoire consistant à choisir au hasard un point de [a, b] (la probabilité de choisir un
point dans un sous-intervalle est proportionnelle à la longueur de ce sous-intervalle).
Si X est une variable aléatoire continue (quelconque), nous avons la propriété suivante:
FX (x) ∼ U [0, 1].
6.4. MODÈLES STOCHASTIQUES 87

6.3.3 Loi de Poisson

Une variable aléatoire X suivant une loi de Poisson est une variable aléatoire qui sert à compter le nombre
d’apparitions d’un phénomène aléatoire durant un intervalle de temps de longueur t. Il pourrait s’ agit par exemple
du nombre d’appels reçus par un téléphoniste. X a alors pour fonction de masse

(θt)k e−θt
PX (k) = P [X = k] = , k = 0, 1, 2, . . . ,
k!
où θ représente le taux moyen.

Exemple 53 (Téléphoniste). Un téléphoniste reçoit en moyenne 2 appels par minute ; quelle est la probabilité de
recevoir moins de 5 appels en 4 minutes ? Nous avons

4 4
X X 8k e−8
P [X < 5] = PX (k) = ≈ 0.1.
k!
k=0 k=0

6.3.4 Loi exponentielle

Soit une variable aléatoire X représentant le temps d’attente entre deux apparitions du phénomène aléatoire en
supposant que le nombre d’apparitions durant un intervalle t suit une loi de Poisson de paramètre θ. La fonction
de répartition vérifie alors
1 − FX (x) = P [X > x] = e−θx , x ≥ 0.

Il s’agit de la loi exponentielle de fonction de densité:

(
θe−θx si x > 0,
fX (x) =
0 sinon.

L’espérance mathématique est:

1
E[X] = .
θ
C’est le taux moyen entre deux apparitions du phénomène aléatoire.

6.4 Modèles stochastiques

Un système stochastique est un système évoluant de manière probabiliste dans le temps. Les exemples sont
nombreux, avec par exemple la température quotidienne ou un centre d’appels téléphoniques. Un modèle stochas-
tique est une représentation mathématique d’un système stochastique. Nous verrons brièvement deux cas classiques
de modèles stochastiques: les processus stochastiques et les files d’attente.

6.4.1 Processus stochastiques

Un processus stochastique est une suite de variables aléatoires évoluant dans le temps, que nous dénoterons
{Xt }, t ∈ T . En général, T est un ensemble discret: T = {0, 1, 2, . . .}. De plus, chaque variable aléatoire peut prendre
une valeur parmi M + 1 états: Xt ∈ {0, 1, . . . , M }.

Exemple 54 (Précipitations quotidiennes).

(
1 s’il y a des précipitations,
Xt =
0 s’il n’y a pas de précipitations.
88 CHAPITRE 6. MODÈLES STOCHASTIQUES

6.4.2 Chaînes de Markov

Un processus stochastique est une chaîne de Markov s’il possède la propriété markovienne:

P [Xt+1 = j | X0 = k0 , X1 = k1 , . . . , Xt−1 = kt−1 , Xt = i] = P [Xt+1 = j | Xt = i].

Cette propriété signifie que la probabilité d’un événement futur, étant donné des événements passés et un état au
temps présent, ne dépend pas du passé, mais uniquement de l’état actuel. La probabilité de transition entre les
états i et j est définie comme
pij = P [Xt+1 = j | Xt = i].
Cette probabilité de transition est dite stationnaire si:

P [Xt+1 = j | Xt = i] = P [X1 = j | X0 = i], t = 1, 2, . . .

Puisqu’il s’agit de probabilité, nous devons en outre avoir

pij ≥ 0, i, j ∈ {0, 1, . . . , M };
M
X
pij = 1 ≥ 0, i ∈ {0, 1, . . . , M }.
j=0

À partir des probabilités de transition, nous pouvons construire

— La matrice des transitions, ayant M + 1 rangées (les états présents) et M + 1 colonnes (les états futurs),
chaque entrée (i, j) de la matrice correspondant à pij .
— Le graphe (ou diagramme) des transitions, ayant M + 1 sommets et tel qu’il y a un arc entre les états i et
j si pij > 0.

Exemple 55 (Précipitations). Supposons que la probabilité qu’il n’y ait pas de précipitations à Montréal demain,
étant donné:
— qu’il n’y en a pas aujourd’hui est 0.8 ;
— qu’il y en a aujourd’hui : 0.6.
Ces probabilités ne changent pas, même si nous tenons compte de ce qui se passe avant aujourd’hui. La propriété
markovienne est satisfaite:

P [Xt+1 = 0 | X0 = k0 , X1 = k1 , . . . , Xt−1 = kt−1 , Xt = 0] = P [Xt+1 = 0 | Xt = 0]

P [Xt+1 = 0 | X0 = k0 , X1 = k1 , . . . , Xt−1 = kt−1 , Xt = 1] = P [Xt+1 = 0 | Xt = 1]

Nous avons donc une chaîne de Markov dont les probabilités de transition sont:

p00 = P [Xt+1 = 0 | Xt = 0] = 0.8,

p10 = P [Xt+1 = 0 | Xt = 1] = 0.6.

Grâce aux propriétés des probabilités de transition, nous pouvons déduire celles qui manquent :

p01 = P [Xt+1 = 1 | Xt = 0] = 1 − 0.8 = 0.2.

p11 = P [Xt+1 = 1 | Xt = 1] = 1 − 0.6 = 0.4.

Ceci donne la matrice de transition:  

0.8 0.2
[P ] =
0.6 0.4
Le graphe de transition est quant à lui représenté sur la Figure 6.1.
6.5. NOTES 89

0.8 0.6 0.4

0 1
0.2

Figure 6.1 – Graphe de transition

Exemple 56 (Marché boursier). À la fin de chaque jour, nous enregistrons le prix de l’action de Google au marché
de Wall Street: (
1 si le prix de l’action n’a pas augmenté à la fin du jour t;
Xt =
0 si le prix de l’action a augmenté à la fin du jour t.
Nous supposons de plus que la probabilité que le prix augmente demain étant donné qu’il a augmenté aujourd’hui
est de 0.7, et qu’il n’a pas augmenté aujourd’hui, 0.5. Nous avons une chaîne de Markov avec comme matrice de
transition:    
p00 p01 0.7 0.3
[P ] = =
p10 p11 0.5 0.5
Supposons maintenant que la probabilité que le prix de l’action augmente demain dépend non seulement de ce
qui est arrivé aujourd’hui, mais également de ce qui est arrivé hier. Le processus stochastique défini précédemment
n’est alors plus une chaîne de Markov. Nous pouvons néanmoins nous en sortir en introduisant un état pour chaque
combinaison d’états possibles sur deux jours consécutifs. Nous définissons alors le processus stochastique suivant,
où l’indice t représente deux jours consécutifs:


 0 si le prix de l’action a augmenté hier et aujourd’hui;

3 si le prix de l’action n’a pas augmenté, ni hier, ni aujourd’hui;
Xt =


 2 si le prix de l’action a augmenté hier, mais pas aujourd’hui;
1 si le prix de l’action a augmenté aujourd’hui, mais pas hier.


Remarquons qu’il est impossible de passer de l’état 0 au temps t à l’état 1 au temps t + 1, car Xt = 0, si le prix
augmente hier et aujourd’hui, et Xt+1 = 1, si le prix augmente demain, mais pas aujourd’hui. La probabilité que le
prix de l’action augmente demain vaut
— s’il a augmenté hier et aujourd’hui: 0.9 ;
— s’il a augmenté aujourd’hui, mais pas hier: 0.6 ;
— s’il a augmenté hier, mais pas aujourd’hui: 0.5 ;
— s’il n’a pas augmenté, ni hier, ni aujourd’hui: 0.3.
La matrice de transition est  
0.9 0 0.1 0
0.6 0 0.4 0 
P =  0 0.5 0 0.5


0 0.3 0 0.7

6.5 Notes
Ce chapitre se base essentiellement sur les notes de cours de Bernard Gendron, 2007.
90 CHAPITRE 6. MODÈLES STOCHASTIQUES
 Programmation dynamique
7
Supposons un problème décomposé en étapes. Chaque étape possède un certain nombre d’états correspondant
aux conditions possibles au début d une étape. À chaque étape, nous devons prendre une décision qui associe à
l’état courant l’état au début de la prochaine étape. Dans le cas d’un problème d optimisation, la programmation
dynamique identifie une politique optimale: une décision optimale à chaque étape pour chaque état possible.

7.1 Principe d’optimalité de Bellman

Étant donné un état courant, la politique optimale pour les étapes à venir est indépendante des décisions prises
aux étapes précédentes. Lorsque ce principe est satisfait, nous pouvons recourir à l’approche suivante pour résoudre
un problàme ayant N étapes:
— identifier une décision optimale pour la dernière étape (N ) ;
— pour chaque étape n, de N − 1 à 1: étant donné une politique optimale pour l’étape n + 1, identifier une
politique optimale pour l’étape n au moyen d une relation de récurrence.
Lorsque n = 1, nous obtenons une politique optimale pour le problème.
Nous utiliserons les notations suivantes:
— N : nombre d’étapes ;
— n: indice de l’étape courante ;
— sn : état au début de l’étape n ;
— xn : variable de décision associée à l’étape n ;
— fn (sn , xn ): contribution des étapes n, n + 1,. . . , N , à la valeur de l objectif lorsqu’on se trouve à l’état sn
au début de l’étape n, que la décision est xn , et que des décisions optimales sont prises par la suite ;
— x∗n : solution optimale à l’étape n ;
— fn∗ (sn ) = fn (sn , x∗n ) = min{fn (sn , xn )} ou max{fn (sn , xn )}.
À chaque itération de la méthode, on construira un tableau représentant
— la liste des états possibles à l’étape n ;
— les valeurs possibles de la variable xn ;
— la contribution à la valeur de l objectif fn (sn , xn ) ;
— la valeur optimale fn∗ (sn ) et la solution optimale x∗n .

sn xn = a xn = b xn = c fn∗ (sn ) x∗n

1 fn (1, a) fn (1, b) fn (1, c) fn∗ (1)
2 fn (2, a) fn (2, b) fn (2, c) fn∗ (2)

91
92 CHAPITRE 7. PROGRAMMATION DYNAMIQUE

Semaine Production max, Production max, Coût unitaire,

temps régulier temps supp. temps régulier
1 2 2 300$
2 3 2 500$
3 1 2 400$

Table 7.1 – Coûts de production

Exemple 57 (Production-inventaire). Une compagnie doit fournir à son meilleur client trois unités du produit P
durant chacune des trois prochaines semaines. Les coûts de production sont donnés dans la Table 7.1. Le coût pour
chaque unité produite en temps supplémentaire est 100$ de plus que le coût par unité produite en temps régulier. Le
coût unitaire d’inventaire est de 50$ par semaine. Au début de la première semaine, il y a deux unités du produit
dans l’inventaire. La compagnie ne veut plus rien avoir dans son inventaire au bout des trois semaines. Combien
doit-on produire d’unités chaque semaine afin de rencontrer la demande du client, tout en minimisant les coûts ?
Une étape correspond ici à une semaine, N = 3, n = 1, 2, 3, et
— sn , l’état au début de la semaine n, est le nombre d’unités de produit dans l inventaire ;
— xn : nombre d’unités produites lors de la semaine n ;
— sn+1 = sn + xn − 3 (puisque nous devons livrer 3 unités au client à chaque semaine) ;
— s1 = 2 (puisqu’il y a 2 unités au début).
Soient cn , le coût unitaire de production au cours de la semaine n, rn la production maximale en temps régulier
pendant la semaine n, et mn la production maximale durant la semaine n. Le coût au cours de la semaine n est
pn (sn , xn ) = cn xn + 100 max(0, xn − rn ) + 50 max(0, sn + xn − 3).
Le coût total vaut dès lors
∗
fn (sn , xn ) = pn (sn , xn ) + fn+1 (sn+1 ),
et le coût optimal répond à la récurrence
fn∗ (sn ) = min{pn (sn , xn ) + fn+1
∗
(sn+1 ) | 3 − sn ≤ xn ≤ mn }.
Si n = 3, nous devons poser
f4∗ (s4 ) = 0.
De plus,
s4 = 0 = s3 + x3 − 3,
aussi
x3 = 3 − s3 .
Calculons d abord les valeurs f3∗ (s3 ) et x∗3 . Nous obtenons le tableau
s3 f3 (s3 , 3 − s3 ) f3∗ (s3 ) x∗3
0 1400 1400 3
1 900 900 2
2 400 400 1
≥3 0 0 0
Voyons maintenant comment nous pouvons calculer les valeurs f2∗ (s2 ) et x∗2 , lorsque s2 = 0. Nous devons avoir
x2 ≥ 3 (car nous devons livrer au moins 3 unités du produit). D’autre part,
f2 (0, 3) = p2 (0, 3) + f3∗ (0) = 1500 + 1400 = 2900,
f2 (0, 4) = p2 (0, 4) + f3∗ (1) = 2150 + 900 = 3050,
f2 (0, 5) = p2 (0, 5) + f3∗ (2) = 2800 + 400 = 3200,
f2∗ (0) = min{f2 (0, 3), f2 (0, 4), f2 (0, 5)} = f2 (0, 3).
Par conséquent, x∗2 = 3. Nous procédons de la même manière pour s2 = 1, 2, 3.
7.1. PRINCIPE D’OPTIMALITÉ DE BELLMAN 93

s2 x2 = 0 x2 = 1 x2 = 2 x2 = 3 x2 = 4 x2 = 5 f2∗ (s2 ) x∗2

0 - - - 2900 3050 3200 2900 3
1 - - 2400 2450 2600 2850 2400 2
2 - 1900 1950 2000 2250 2900 1900 1
≥3 1400 1450 1500 1650 2300 2950 1400 0

Pour la première étape (n = 1), nous avons s1 = 2 (il y a 2 unités au départ dans l inventaire). Nous devons donc
avoir x1 ≥ 1. De plus,

f1 (2, 1) = p1 (2, 1) + f2∗ (0) = 300 + 2900 = 3200,

f1 (2, 2) = p1 (2, 2) + f2∗ (1) = 650 + 2400 = 3050,
f1 (2, 3) = p1 (2, 3) + f2∗ (2) = 1100 + 1900 = 3000,
f1 (2, 4) = p1 (2, 4) + f2∗ (3) = 1550 + 1400 = 2950.

Par conséquent, f1∗ (2) = f1 (2, 4) et x1 ∗ = 4. Sous forme de tableau, cela donne:

s1 x1 = 1 x1 = 2 x1 = 3 x1 = 4 f1∗ (s1 ) x∗1

2 3200 3050 3000 2950 2950 4

La politique optimale est donc:

x∗1 = 4, x∗2 = 0, x∗3 = 3
donnant lieu à s2 = 3, s3 = 0. Le coût total est f1∗ (2) = 2950$.

Exemple 58 (Affectation de ressources). Trois équipes de recherche travaillent sur un même problème avec chacune
une approche différente. La probabilité que l’équipe Ei échoue est P [E1 ] = 0.4, P [E2 ] = 0.6, P [E3 ] = 0.8. Ainsi,
la probabilité que les trois équipes échouent simultanément est: (0.4)(0.6)(0.8) = 0.192. Nous voulons ajouter deux
autres scientifiques afin de minimiser la probabilité d’échec. À quelles équipes faut-il affecter ces deux scientifiques
afin de minimiser la probabilité d’échec ? La Table 58 représente les probabilités d’échec pour chaque équipe en
fonction du nombre de scientifiques ajoutés.

Scientifiques ajoutés E1 E2 E3
0 0.4 0.6 0.8
1 0.2 0.4 0.5
2 0.15 0.2 0.3

Table 7.2 – Affectation de ressources.

Nous faisons correspondre l’équipe n à l’étape n, pour n = 1, 2, 3 (N = 3). L’état sn est le nombre de scientifiques
non encore affectés. La décision xn correspond au nombre de scientifiques affectés à l’équipe En . Soit pn (xn ) la
probabilité d’échec de l’équipe En si xn scientifiques lui sont affectés. Alors,
∗
fn (sn , xn ) = pn (xn ).fn+1 (sn − xn )
fn∗ (sn ) = min{fn (sn , xn ) | xn = 0, 1, . . . , sn },

avec f4∗ (0) = 1. Calculons d abord f3∗ (s3 ) et x∗3 . Nous devons avoir s3 = x3 (nous affectons tous les scientifiques
non encore affectés), de sorte que nous avons le tableau

s3 f3 (s3 , s3 ) f3∗ (s3 ) x∗3

0 0.8 0.8 0
1 0.5 0.5 1
2 0.3 0.3 2
94 CHAPITRE 7. PROGRAMMATION DYNAMIQUE

Pour calculer les valeurs f2∗ (s2 ) et x∗2 , lorsque s2 = 1, observons d’abord que nous devons avoir x2 ≤ 1. Nous avons

f2 (1, 0) = p2 (1, 0) × f3∗ (1) = 0.6 × 0.5 = 0.3,

f2 (1, 1) = p2 (1, 1) × f3∗ (0) = 0.4 × 0.8 = 0.32,
f2 ∗ (1) = min{f2 (1, 0), f2 (1, 1)} = f2 (1, 0).

Dès lors, x∗2 = 0. Nous pouvons procéder de même pour s2 = 0 et 2. Ceci donne le tableau ci-dessous.

s2 x2 = 0 x2 = 1 x2 = 2 f2∗ (s2 ) x∗2

0 0.48 - - 0.48 0
1 0.3 0.32 - 0.3 0
2 0.18 0.2 0.16 0.16 2

Pour la dernière étape (n = 1), nous avons s1 = 2 (il y a 2 scientifiques à affecter). Nous avons

f1 (2, 0) = p1 (2, 0).f2∗ (2) = 0.4 × 0.16 = 0.064,

f1 (2, 1) = p1 (2, 1).f2∗ (1) = 0.2 × 0.3 = 0.06,
f1 (2, 2) = p1 (2, 2).f2∗ (0) = 0.15 × 0.48 = 0.072.

Aussi, f1∗ (2) = f1 (2, 1) et x∗1 = 1.

s1 x1 = 0 x1 = 1 x1 = 2 f1∗ (s1 ) x∗1

0 0.064 0.06 0.072 0.06 1

La politique optimale est donc: x∗1 = 1 (d’où s2 = 1), x∗2 = 0 (d’où s3 = 1), x∗3 = 1. La probabilité d’échec est:
f1∗ (1) = 0.06.

7.2 Affectation de ressources

La programmation dynamique est souvent utilisée pour résoudre des problèmes d affectation de ressources.
C’était le cas dans notre dernier exemple, où il y avait une ressource à affecter (les scientifiques). Que faire lorsqu’il
y a plus d une ressource ? Dans ce cas, chaque état est un vecteur dont chaque coordonnée i représente la quantité
de la ressource i non encore affectée.

Exemple 59 (problème Wyndor Glass). Nous considérons deux types de produits (produit 1, produit 2) et trois
usines (usine 1, usine 2, usine 3). La ressource i est la capacité de production pour l’usine i, exprimée en heures
par semaine. Nous connaissons le profit par lot (20 unités) de chaque produit, et chaque lot du produit 1 (2) est
le résultat combiné de la production aux usines 1 et 3 (2 et 3). Nous souhaitons déterminer le taux de production
pour chaque produit (nombre de lots par semaine) de façon à maximiser le profit total. Les données sont reproduites
dans la Table 59. L’étape n sera le produit n, n = 1, 2 (N = 2). Nous avons comme état sn (R1 , R2 , R3 )n , où Ri

Produit 1 (temps de Produit 2 (temps de Capacité de

production, h/lot) production, h/lot) production (h)
Usine 1 1 0 4
Usine 2 0 2 12
Usine 3 3 2 18
Profit($)/lot 3000 5000

Table 7.3 – Allocation de ressources pour le problème Wyndor Glass

est la quantité de la ressource i (temps de production à l usine i) non encore affectée. La décision xn est le nombre
7.2. AFFECTATION DE RESSOURCES 95

de lots du produit n, et
∗
fn (sn , xn ) = cn xn + fn+1 (sn+1 ) (f3∗ (s3 ) = 0),
f2∗ (s2 ) = max{5x2 | 2x2 ≤ 12, 2x2 ≤ 18 − 3x1 , x2 ≥ 0}
f1∗ (s1 ) = max{3x1 + f2∗ (s2 ) | x1 ≤ 4, 3x1 ≤ 18, x1 ≥ 0}
De plus,
s1 = (4, 12, 18), s2 = (4 − x1 , 12, 18 − 3x1 ).
Afin de résoudre le problème, nous cherchons d abord f2∗ (s2 ) et x∗2 . Nous avons
f2∗ (s2 ) = max{5x2 | 2x2 ≤ 12, 2x2 ≤ 18 − 3x1 , x2 ≥ 0}.
Les contraintes
2x2 ≤ 12 et 2x2 ≤ 18 − 3x1
impliquent  
18 − 3x1
x2 ≤ min 6, .
2
Par conséquent   
18 − 3x1
f2∗ (s2 ) = max 5x2 | 0 ≤ x2 ≤ min 6, ,
2
et  
18 − 3x1
x∗2 = min 6, . (7.2.1)
2
Nous avons donc  
18 − 3x1
f2∗ (s2 )
= 5 min 6, | 3x1 ≤ 18 .
2
Calculons maintenant f1∗ (s1 ) et x∗1 . Nous avons
f1∗ (s1 ) = max{3x1 + f2∗ (s2 ) | x1 ≤ 4, 3x1 ≤ 18, x1 ≥ 0}.
Or,
x1 ≤ 4
implique
3x1 ≤ 18.
Ainsi    
18 − 3x1
f1∗ (s1 ) = max 3x1 + 5 min 6, | 0 ≤ x1 ≤ 4 .
2
Or, (
 
18 − 3x1 6 si 0 ≤ x1 ≤ 2,
min 6, = 18−3x1
2 2 si 2 ≤ x1 ≤ 4.
Par conséquent,
f1∗ (s1 ) = max{f1 (s1 , x1 ) | 0 ≤ x1 ≤ 4},
avec (
3x1 + 30 si 0 ≤ x1 ≤ 2,
f1 (s1 , x1 ) =
45 − 29 x1 si 2 ≤ x1 ≤ 4.
Le maximum est atteint en x∗1 = 2 et f1∗ (s1 ) = 36.
Comme, de (7.2.1),  
18 − 3x1
x∗2 = min 6, ,
2
et x∗1 = 2, nous avons x∗2 = 6. La politique optimale est donc: x∗1 = 2, x∗2 = 6 de valeur f1∗ (s1 ) = 36.
En utilisant la même approche, il est possible de résoudre par la programmation dynamique des modèles où
— l’objectif ou les contraintes sont non linéaires ;
— certaines variables sont entières.
96 CHAPITRE 7. PROGRAMMATION DYNAMIQUE

7.3 Programmation dynamique déterministe et plus court chemin

Supposons que pour s1 fixé, nous voulions résoudre
N
X −1
min gN (xN ) + gk (sk , xk ).
µ0 ,...,µN −1
k=1

Nous pouvons calculer les décisions optimales x1 , . . . , xN dès le départ, car aucune nouvelle information n’est
obtenue en cours de route. Si les ensembles d’états et de décisions sont finis à chaque étape, résoudre ce problème
équivaut à trouver un plus court chemin dans un réseau, où les noeuds sont tous les états (k, sk ) possibles, pour
1 ≤ k ≤ N , auxquels nous ajoutons un noeud artificiel t qui correspond à l’état où tout est terminé (étape N + 1).
Pour chaque noeud (k, sk ), k < N , et chaque décision xk , nous introduisons un arc de longueur g(sk , xk ) partant
du noeud (k, sk ) et allant au noeud (k + 1, sk+1 ), avec sk+1 = gk (sk , xk ). Chaque noeud (N, sN ) est relié au noeud
t par un arc de longueur gN (xN ). Nous cherchons alors un plus court chemin de (0, s0 ) à t.
...

gN
(x
N
... gN

)
(x
N )
État initial ) Noeud terminal artificiel
(x N
gN
...

N
)
(x
gN
...

Étape 0 Étape 1 Étape 2 Étape N − 1 Étape N

Si on numérote les noeuds couche par couche, par ordre croissant de valeur de k, on obtient un réseau sans cycle
et ordonné topologiquement (i.e., un arc (i, j) ne peut exister que si i < j). Dans le cas où il n’est pas nécessaire
de mémoriser le numéro d’étape, on peut simplifier le réseau en agrégeant des noeuds. Le réseau résultant peut
ne plus être ordonné topologiquement. Inversement, tout problème de recherche d’un plus court chemin dans un
réseau peut se formuler comme un problème de programmation dynamique déterministe (avec des coûts additifs
entre étapes).

7.4 Cas probabiliste

Dans la programmation dynamique probabiliste, la décision optimale à l’étape n dépend de la loi de probabilité
de l’état à l’étape n + 1. Ainsi, nous passerons de l’étape n à l’étape n + 1, se retrouvant ainsi dans l’état i, avec
une probabilité pi . Étant donné S états à l’étape n + 1, nous aurons pi ≥ 0 et
S
X
pi = 1.
i=1

La relation de récurrence dépend de cette loi de probabilité et de l’objectif à optimiser.

Exemple 60 (Jeu de hasard). Le jeu consiste à miser un nombre quelconque de jetons. Si nous gagnons, nous
remportons le nombre de jetons misés. Si nous perdons, nous perdons le nombre de jetons misés. Un statisticien
croit pouvoir gagner chaque jeu avec une probabilité égale à 2/3. Ses collègues parient avec lui qu’en misant au
départ 3 jetons, il aura moins de 5 jetons après 3 parties. Combien de jetons miser à chacune des 3 parties ?
Le modèle peut se formuler comme suit: l’étape n correspond à la partie n, n = 1, 2, 3 (N = 3). L’état sn est
le nombre de jetons au début de la partie n, et la décision xn correspond au nombre de jetons à parier à la partie
n. Nous souhaitons maximiser la probabilité d’avoir au moins 5 jetons après 3 parties. Soit fn (sn , xn ) la probabilité
7.4. CAS PROBABILISTE 97

de terminer avec au moins 5 jetons, étant donné que nous sommes dans l’état sn à l’étape n, que nous misons xn
jetons et que nous effectuons des décisions optimales aux étapes n + 1, . . . , N . Nous avons

fn∗ (sn ) = max{fn (sn , xn ) | xn = 0, 1, . . . , sn }.

étant donné que nous dans l’état sn à l’étape n et que nous misons xn jetons, nous pouvons:
— perdre et se retrouver dans l’état sn − xn avec une probabilité 1/3 ;
— gagner et se retrouver à l’état sn + xn avec une probabilité 2/3.
Nous avons la fonction de récurrence
1 ∗ 2 ∗
fn (sn , xn ) = fn+1 (sn − xn ) + fn+1 (sn + xn ).
3 3
De plus, (
0 si s4 < 5,
f4∗ (s4 ) =
1 si s4 ≥ 5.
Calculons d abord f3∗ (s3 ) et x∗3 . Nous avons le tableau
s3 x3 = 0 x3 = 1 x3 = 2 x3 = 3 x3 = 4 f3∗ (s3 ) x∗3
0 0 - - - - 0 0
1 0 0 - - - 0 ≥0
2 0 0 0 - - 0 ≥0
2 2 2
3 0 0 3 3 - 3 ≥2
2 2 2 2 2
4 0 3 3 3 3 3 ≥1
≥5 1 1 0

Voyons maintenant comment nous pouvons calculer les valeurs f2∗ (s2 ) et x∗2 , lorsque s2 = 3. Il est clair que
nous devons avoir x2 ≤ 3. De plus,
1 ∗ 2 ∗ 2
f2 (3, 0) = f (3) + f (3) = ;
3 3 3 3 3
1 2 ∗ 4
f2 (3, 1) = f3∗ (2) + f (4) = ;
3 3 3 9
1 ∗ 2 ∗ 2
f2 (3, 2) = f3 (1) + f (5) = ;
3 3 3 3
1 ∗ 2 ∗ 2
f2 (3, 3) = f3 (0) + f (6) = .
3 3 3 3

Ainsi,
f2∗ (3) = max{f2 (3, 0), f2 (3, 1), f2 (3, 2), f2 (3, 3)} = f2 (3, 0).
Dès lors, x∗2 = 0, 2 ou 3. De la même façon, pour s2 , nous avons le tableau,
s2 x2 = 0 x2 = 1 x2 = 2 x2 = 3 x2 = 4 f2∗ (s2 ) x∗2
0 0 - - - - 0 0
1 0 0 - - - 0 ≥0
4 4 4
2 0 9 9 - - 9 1, 2
2 4 2 2 2
3 3 9 3 3 - 3 0, 2, 3
2 8 2 2 2 8
4 3 9 3 3 3 9 ≥1
≥5 1 1 0

Pour la première étape (n = 1), nous avons s1 = 3, comme nous pouvons affecter 3 jetons au départ du jeu. Ceci
donne le tableau
s1 x1 = 0 x1 = 1 x1 = 2 x1 = 3 f1∗ (s2 ) x∗1
2 20 2 2 20
3 3 27 3 3 27 1
98 CHAPITRE 7. PROGRAMMATION DYNAMIQUE

Dès lors
1 ∗ 2 ∗ 2
f1 (3, 0) = f (3) + f (3) = ,
3 2 3 2 3
1 2 ∗ 20
f1 (3, 1) = f2∗ (2) + f (4) = ,
3 3 2 27
1 ∗ 2 ∗ 2
f1 (3, 2) = f2 (1) + f (5) = ,
3 3 2 3
1 ∗ 2 ∗ 2
f1 (3, 3) = f2 (0) + f (6) = ,
3 3 2 3
et par conséquent,
f1∗ (3) = max{f1 (3, 0), f1 (3, 1), f1 (3, 2), f1 (3, 3)} = f1 (3, 1).
Nous en déduisons
x∗1 = 1.
La politique optimale est donc :
— x∗1 = 1 ;
— Si nous gagnons (s2 = 4), x∗2 = 1.
— Si nous gagnons (s3 = 5), x∗3 = 0.
— Si nous perdons (s3 = 3), x∗3 = 2 ou 3.
— Si nous perdons (s2 = 2), x∗2 = 1 ou 2.
— Si nous gagnons (s3 = 3 ou 4), x∗3 = 2 ou 3 (si x∗2 = 1) ou 1 ≤ x∗3 ≤ 4 (si x∗2 = 2).
— Si nous perdons (s3 = 1 ou 0), x∗3 ≥ 0 (mais le pari est perdu !).

7.5 Notes
Ce chapitre se base essentiellement sur les notes de cours de Bernard Gendron, 2007. Le lien avec les plus courts
chemins est développé notamment dans [1].
 8
Simulation

8.1 Introduction
Simuler un système stochastique consiste à imiter son comportement pour estimer sa performance Un modèle
de simulation est une représentation du système stochastique permettant de générer un grand nombre d’événements
aléatoires et d’en tirer des observations statistiques.
Exemple 61 (Jeu de hasard). Chaque partie consiste à tirer une pièce de monnaie jusqu’à ce que la différence
entre le nombre de faces et le nombre de piles soit égale à 3. Chaque tirage coûte 1$, et chaque partie jouée rapporte
8$ au joueur. Nous aurons par exemple les jeux
— FFF: gain de 8$-3$=5$ ;
— PFPPP: gain de 8$-5$=3$ ;
— PFFPFPFPPPP: perte de 8$-11$=3$.
Nous nous demandons par conséquent s’il est intéressant de jouer. Pour répondre à cette question, nous allons
simuler le jeu. Il y a deux façons de le faire:
— nous pouvons jouer pendant un certain temps sans miser d’argent ;
— nous pouvons simuler le jeu par ordinateur.
Néanmoins, simuler une seule partie ne nous aide pas à prendre une décision. Pour cela, il faudrait voir ce qui se
passe sur un grand nombre de parties et mesurer le gain moyen (ou la perte moyenne).
Dans cet exemple, nous pouvons définir les éléments d’un modèle de simulation comme suit
système stochastique : tirages successifs ;
horloge : nombre de tirages ;
définition de l’état du système : N (t), le nombre de faces moins le nombre de piles après t tirages ;
événements modifiant l’état du système : tirage de pile ou de face ;
méthode de génération d’événements : génération d’un nombre aléatoire uniforme ;
formule de changement d’état : N (t + 1) = N (t) + 1, si F est tirée ; N (t) − 1, si P est tirée ;
performance : 8 − t, lorsque N (t) atteint +3 ou -3.

8.2 Files d’attente

8.2.1 Concepts de base
Un système de file d’attente consiste d’un ou de plusieurs serveurs qui fournissent un service d’une certaine
nature à des clients qui se présentent, et qui sont issus d’une population donnée. Les clients qui arrivent et trouvent

99
100 CHAPITRE 8. SIMULATION

tous les serveurs occupés rejoindront (généralement) une ou plusieurs files (ou lignes) devant les serveurs, d’où la
qualification de file d’attente. . Quelques exemples de files d’attente sont repris dans la Table 8.1. Le système de file

système Serveurs Clients

Banque Guichets Clients
Hôpital Docteurs, infirmiéres, lits Patients
Système informatique Unités centrales, dispositifs I/O Travaux
Manufacture Machines, travailleurs Composants
Centre d’appel Serveurs Appels entrants

Table 8.1 – Exemples de files d’attente

d’attente est caractérisé par trois composants:

1. le processus d’arrivée ;
2. le mécanisme de service ;
3. la discipline de file.
Spécifier le processus d’arrivée pour un système de file d’attente revient à décrire comment les clients arrivent
dans le système. Soit Ai le temps d’inter-arrivées entre les arrivées des (i − 1)e et ie clients. Si A1 , A2 , . . ., sont
supposés être des variables i.i.d., nous dénoterons le temps d’interarrivée moyen (ou espéré) par E[A] et appellerons
λ = 1/E[A] le taux d’arrivée des clients.
Le mécanisme de service pour un système de file d’attente est articulé en spécifiant le nombre de serveurs
(dénoté par s), où chaque serveur a sa propre file ou au contraire une seule file alimente tous les serveurs, et la
distribution de probabilité des temps de service des clients. Soit Si le temps de service du ie client. Si S1 , S2 , . . .,
sont des variables aléatoires i.i.d., nous dénoterons le temps de service moyen d’un client par E[S] et appelleront
ω = 1/E[S] le taux de service d’un serveur.
Nous parlerons de situation transitoire lorsque l’état du système dépend grandement de la situation initiale et
du temps écoulé, et de situation d’équilibre lorsque l’état du système peut être considéré indépendant de la situation
initiale et du temps écoulé. En situation d’équilibre, nous avons la formule de Little:

L = λW,

où L est le nombre moyen de clients dans le système, λ, le taux moyen d’arrivée des nouveaux clients, et W le temps
moyen dans le système.

8.2.2 Modèle M/M/1

Il s’agit du modèle de file d’attente le plus courant, avec les caractéristiques suivantes:
— file d’attente: nombre infini de clients ;
— stratégie de service: premier arrivé, premier servi (FIFO) ;
— un seul serveur ;
— les taux d’arrivée et de service obéissent à des lois de Poisson. De manière équivalente, le temps entre l’arrivée
de deux clients successifs et le temps de service obéissent à des lois exponentielles: on parle de processus
Markoviens.
La notation générale pour les modèles de file d’attente est X/Y /s, où X représente la loi du temps d’interarrivée,
Y est la loi du temps de service, et s donne le nombre de serveurs.

Exemple 62 (File d’attente M/M/1). En situation d’équilibre, plusieurs résultats analytiques (obtenus par analyse
du modèle mathématique) existent. Soit λ le taux moyen d’arrivée et µ le taux moyen de service. Supposons que
λ < µ Le nombre moyen de clients dans le système vaut
λ
L= .
µ−λ
8.2. FILES D’ATTENTE 101

Le temps moyen d’attente dans le système vaut quant à lui

1
W = .
µ−λ
Il est aisé de vérifier ces résultats par simulation. Dans le cas présent, les éléments de la simulation sont:
système stochastique : file d’attente M/M/1 ;
horloge : temps écoulé ;
définition de l’état du système : N (t), le nombre de clients dans le système au temps t ;
événements modifiant l’état du système : arrivée ou fin de service d’un client ;
formule de changement d’état : N (t + 1) = N (t) + 1, si arrivée ; N (t) − 1, si fin de service.
Nous allons voir deux méthodes pour étudier l’evolution du système dans le temps:
— par intervalles de temps fixe ;
— par génération d’événement.
Nous supposerons que les valeurs des paramètres de notre système sont λ = 3 clients par heure, et µ = 5 clients
par heure.
Dans le cas d’intervalles de temps fixe, l’idée est
1. faire écouler le temps d’un petit intervalle ∆t ;
2. mettre à jour le système en déterminant les événements qui ont pu se produire durant l’intervalle ∆t, et
recueillir l’information sur la performance du système ;
3. retourner en 1.
Ici, les événements sont soit des arrivées, soit des départs (fins de service). Si ∆t est suffisamment petit, nous
pouvons considérer qu’il ne se produira qu’un seul événement (arrivée ou départ) durant cet intervalle de temps.
Prenons ∆t égal à 0.1 heure (6 minutes). La probabilité qu’il y ait une arrivée durant cet intervalle de temps
vaut:
PA = 1 − e−λ∆t = 1 − e−3/10 ≈ 0.259.
La probabilité qu’il y ait un départ durant cet intervalle de temps est:

PD = 1 − e−µ∆t = 1 − e−5/10 ≈ 0.393.

Pour générer un événement, nous pouvons simplement tirer deux nombres aléatoires selon une loi U [0, 1]. Si le
premier nombre est strictement plus petit que 0.259, nous aurons une arrivée. Si le deuxième nombre est strictement
inférieur à 0.393, nous considérerons un départ (si un client était en cours de service). Ceci pourrait par exemple
donner les chiffres de la Table 8.2.
D’après cet exemple, il est possible d’estimer les performances du système. Que se passe-t-il si nous voulons
mesurer W , le temps moyen passé dans le système ? Nous avons deux clients qui sont entrés dans le système
et chacun y est resté 18 minutes ou 0.3 heures, aussi peut-on estimer W par 0.3. Pourtant, la vraie valeur est
W = 1/(µ − λ) = 0.5. Il faudrait un échantillon beaucoup plus grand, et ce d’autant plus que nous devrions simuler
le système en état d’équilibre.
La génération d’événement peut se résumer comme suit:
1. faire écouler le temps jusqu’au prochain événement ;
2. mettre à jour le système en fonction de l’événement qui vient de se produire et générer aléatoirement le temps
jusqu’au prochain événement ; recueillir l’information sur la performance du système ;
3. retourner à 1.
Nous pourrions ainsi obtenir la suite d’événements reprise dans la Table 8.3.
Si nous considérons une simulation comportant l’arrivée de 10000 clients, les résultats montrent que:
— le nombre moyens de clients dans le système est L ≈ 1.5 ;
— le temps moyen dans le système est W ≈ 0.5.
Remarque: des résultats analytiques existent pour des modèles simples (comme M/M/1), mais pas pour des files
d’attente plus complexes.
102 CHAPITRE 8. SIMULATION

t (min) N (t) Nombre 1 Arrivée Nombre 2 Départ

0 0
6 1 0.096 Oui -
12 1 0.569 Non 0.665 Non
18 1 0.764 Non 0.842 Non
24 1 0.492 Non 0.224 Oui
30 0 0.950 Non -
36 0 0.610 Non -
42 1 0.145 Oui -
48 1 0.484 Non 0.552 Non
54 1 0.350 Non 0.590 Non
60 0 0.430 Non 0.041 Oui

Table 8.2 – Simulation par intervalles de temps fixe.

t (min) N (T ) Temps Temps Prochaine Prochain Prochain

d’interarrivée de service arrivée départ événement
0 0 2.019 - 2.019 - Arrivée
2.019 1 16.833 13.123 18.852 15.142 Départ
15.142 0 - - 18.852 - Arrivée
18.852 1 28.878 22.142 47.730 40.994 Départ
40.994 0 - - 47.730 - Arrivée
47.730 1

Table 8.3 – Simulation par événements.

8.3. SIMULATION À ÉVÉNEMENTS DISCRETS 103

8.3 Simulation à événements discrets

L’approche s’est imposée dans nombre de problèmes de simulations. Le principe général est résumé dans l’ordi-
nogramme de la Figure 8.1.

Planification
Horloge de du ou des pre-
simulation Initialisation miers(s) évé-
mise à 0 nement(s)

Gestion de
l’événement
actuel

Dépilement
du prochain Planification
événement et du ou des
avancement événements
de l’horloge en découlant.
de simulation

Liste
non d’événements
vide ?

oui

Arrêt

Figure 8.1 – Schéma d’une simulation à événements discrets.

8.4 Notes
Ce chapitre se base essentiellement sur les notes de cours de Bernard Gendron, 2007. Pour plus de détails sur
la simulation à événements discrets, nous renvoyons le lecteur aux notes du cours IFT3245, du même auteur: .
104 CHAPITRE 8. SIMULATION
 Annexe

8.5 Logiciels d’optimisation

8.5.1 IOR Tutorial
Le logiciel IOR Tutorial, qui accompagne le livre Introduction to Operational Research [2] est disponible gratuite-
ment à l’adresse [Link]
Le logiciel est fourni sous la forme d’un exécutable sous l’environnement Microsoft Windows, et comme archive
Java pour les autres plate-formes. Il est par conséquent utilisable sur la majorité des systèmes actuels. Le logiciel
est intéressant pour rapidement tester des programmes d’optimisation simples, mais ne peut être considéré comme
un solveur à part entière. Son interface graphique le rend néanmoins adapté pour une expérimentation rapide avec
les exemples introduits dans lecours.

8.5.2 GAMS
GAMS est l’acronyme de Generalized Algebric Modeling System. Il s’agit avant tout d’un langage servant à
la formulation et à la résolution de modèle de programmation mathématique. En pratique, il s’agit d’un paquetage
intégré qui permet de
— spécifier la structure du modèle d’optimisation ;
— spécifier et calculer les données qui entrent dans le modèle ;
— résoudre ce modèle ;
— produire un rapport sur un modèle ;
— conduire une analyse comparative.
GAMS peut être téléchargé à l’adresse [Link] La licence de GAMS et des solveurs académique est
passablement élevés, toutefois il est utilisable de manière gratuite en version de démonstration, ce qui permettra
d’illustrer certains concepts de base du cours.
Sous des environnements de type UNIX (incluant Mac OS X et la majorité des distributions Linux), GAMS
s’utilise uniquement en ligne de commande. Sous Microsoft Windows, il dispose d’une interface graphique.

Formulation d’un problème simple

Exemple 63. Considérons le problème d’optimisation linéaire

max 109x1 + 90x2 + 115x3

x
s.c. x1 + x2 + x3 ≤ 100;
6x1 + 4x2 + 8x3 ≤ 500;
x1 , x2 , x3 ≥ 0.

Le fichier GAMS décrivant se problème est constitué des parties suivantes:

1. spécification des variables ;
2. spécification des équations:
— déclaration ;
— spécification de la structure algébrique ;
3. définition du modèle ;

105
106 ANNEXE

4. définition de la méthode de résolution.

Ce problème peut être formulé sous GAMS comme suit:

VARIABLES
Z Variable Z;

POSITIVE VARIABLES
X1 Variable X1
X2 Variable X2
X3 Variable X3;

EQUATIONS
Equation1 Equation 1
Equation2 Equation 2
Equation3 Equation 3;

Equation1..
Z =E= 109*X1 + 90*X2 + 115*X3;

Equation2..
X1 + X2 + X3 =L= 100;

Equation3..
6*X1 + 4*X2 + 8*X3 =L= 500;

MODEL Example1 /ALL/;

SOLVE Example1 USING LP MAXIMIZING Z;

Remarquons d’ores et déjà que chaque instruction se termine par un point-virgule. Leur omission produit une erreur
de compilation. Détaillons les instructions.

Spécification des variables

GAMS exige que les variables soient identifiées. Dans l’exemple précédent, nous avons les variables Z, X1, X2,
X3:

VARIABLES
Z Variable Z;

POSITIVE VARIABLES
X1 Variable X1 (Optional Text)
X2 Variable X2
X3 Variable X3;

Pour tout modèle construit avec GAMS, il convient d’identifier l’objectif à l’aide d’une variable, ici Z. Autrement
dit,
max cT x
x

devient

max z
s.c. z = cT x,
Z = CX.
8.5. LOGICIELS D’OPTIMISATION 107

Les noms de variables peuvent avoir jusqu’à 31 caractères. Il est également possible de spécifier plus précisément le
type des variables:
VARIABLE variables non restreintes
POSITIVE VARIABLE variables non négatives
NEGATIVE VARIABLE variables non positives
BINARY VARIABLE variables binaires (dans {0, 1})
INTEGER VARIABLE variables entières (naturelles)

Dénomination des équations

La spécification d’une équation consiste en deux parties. Il convient tout d’abord de déclarer ces équations:
GAMS exige du modélisateur de nommer chaque équation à l’oeuvre dans le modèle. Dans l’exemple, les équations
sont déclarées après le mot-clé EQUATIONS

EQUATIONS
Equation1 Equation 1
Equation2 Equation 2
Equation3 Equation 3;

Comme pour les variables, le nom d’une équation peut prendre jusqu’à 31 caractères. A droite du nom de l’équation
figure un (bref) texte d’explications.
Il nous faut de plus spécifier la structure algébrique: après avoir nommé les équations, la structure algébrique
exacte doit être spécifiée en utilisant la notation ..:

Equation1..
Z =E= 109*X1 + 90*X2 + 115*X3;

Equation2..
X1 + X2 + X3 =L= 100;

Equation3..
6*X1 + 4*X2 + 8*X3 =L= 500;

La forme algébrique exige l’utilisation d’une syntaxe speciale afin de définir la forme exacte de l’équation:
=E= contrainte d’égalité
=L= inférieur ou égal
=G= supérieur ou égal

Spécification du modèle
Le mot-clé MODEL est utilisé pour identifier les modèles à résoudre. Il convient de
1. donner un nom au modèle (par exemple Example1) ;
2. spécifier les équations à inclure, entre des barres italiques / /.
Ceci donne pour notre exemple

MODEL Example1 /ALL/;

Nous pourrions aussi avoir

MODEL Example1 /Equation1, Equation2/;

108 ANNEXE

Spécification du solveur
Le mot-clé SOLVE indique à GAMS d’appliquer un solveur au modèle nommé, en utilisant les données définies
juste avant l’instruction SOLVE. Ainsi, dans notre exemple, nous avions:

SOLVE Example1 USING LP MAXIMIZING Z;

Si nous avions affaire à un problème linéaire de minimisation, nous pourrions écrire

SOLVE Example1 USING LP MINIMIZING Z;

En place de LP, nous devrions écrire MIP pour traiter un problème de programmation entière mixte:

SOLVE Example1 USING MIP MAXIMIZING Z;

De même, nous spécifierons un problème non-linéaire avec le mot-clé NLP:

SOLVE Example1 USING NLP MAXIMIZING Z;

Rapport de solution
A la fin de l’exécution, GAMS produit un rapport indiquant la solution trouvée, la valeur de la fonction
objectif en cette solution, ainsi que différentes informations permettant d’analyser le comportement de l’algorithme
d’optimisation, et diverses propriétés du problème en cours d’étude. En particulier, le résumé de rapport donne le
nombre total de cas non optimaux, non réalisables et non bornés rencontrés.

**** REPORT SUMMARY : 0 NONOPT

0 INFEASIBLE
0 UNBOUNDED

L’information sur les solutions peut être affichée de différentes manières:

1. sortie standard de GAMS ;
2. utilisation des commandes DISPLAY ;
3. rappports additionels sur base des valeurs des solutions.
La sortie standard de GAMS présente la solution sous forme d’un tableau, qui dans le cas de l’exemple est:

LOWER LEVEL UPPER MARGINAL

---- EQU Equation1 . . . 1.000

---- EQU Equation2 -INF 100.000 100.000 52.000
---- EQU Equation3 -INF 500.000 500.000 9.500

Equation1 Equation 1
Equation2 Equation 2
Equation3 Equation 3

LOWER LEVEL UPPER MARGINAL

---- VAR Z -INF 9950.000 +INF .

---- VAR X1 . 50.000 +INF .
---- VAR X2 . 50.000 +INF .
---- VAR X3 . . +INF -13.000

Le point simple “.” représente un zéro, tandis que INF représente l’infini.
8.5. LOGICIELS D’OPTIMISATION 109

Sommations
P
La notation mathématique j xj se traduira par

SUM(j,x(j))

En d’autres termes, nous avons la syntaxe

SUM( index of summation, names(index))

P P
Il est possible d’imbriquer les sommations. Ainsi, j i xji donnera

SUM(j,SUM(i,x(j,i))

ou encore

SUM((j,i),x(j,i))

Définition d’ensemble
Spécifier les variables une à une est vite fastidieux, pour ne pas dire irréaliste (par exemple si nous avons un
million de variables), c’est pourquoi en modélisation algébrique, nous utilisons des indices. GAMS met à disposition
le mot clé SET pour définir des ensembles, parcouru par un indice spécifié juste avant:

SET ItemName optional explanatory text for item

/ element1 optional explanatory text for element,
element2 optional explanatory text for element /

Par exemple,

SET i index /1*10/

défini l’ensemble {1, 2, . . . , 10}, qui peut être parcouru avec l’indice i. Il est possible d’associer un autre indice au
même ensemble au moyen de la commande ALIAS, par exemple

Exemple 64. ALIAS (i,j) ;

permet d’utiliser j au lieu de i. Cette commande est utile par exemple pour traduire une contrainte de la forme

xij + xji = 1, i = 1, . . . , n, j = 1, . . . , n,

j et i doivent indicer le même ensemble, mais écrire

SET i / 1*10 /
j / 1*10 /

mènerait GAMS à considérer i et j comme indexant deux ensembles.

Entrée de données
Les données sont entrées au moyen de quatre différents types de commandes GAMS:
SCALAR , pour les éléments ne dépendant pas d’ensembles ;
PARAMETERS , pour les éléments qui sont des vecteurs ;
TABLES , pour des éléments de deux dimensions ou plus ;
PARAMETERS , en affectation directe.
La manière la plus simple d’entrée des données est l’emploi de la commande SCALAR, qui prend la syntaxe suivante
dans le format basique:
110 ANNEXE

SCALAR ItemName optional text / value /;

Dans le cas de vecteurs, nous utiliserons la syntaxe
PARAMETER ItemName(Set) optional text;
/ element1 value,
element2 value /;
Pour les entrées multidimensionnelles, nous aurons
TABLE ItemName(set1dep,set2dep) optional text
set2elem1 set2elem2
set1element1 value11 value12
set1element2 value12 value22 ;
Plutôt que d’utiliser des constantes, nous pouvons affecter les valeurs au moyen d’expressions mathématiques:
PARAMETER ItemName(set1dep,set2dep) optional text;
ItemName(set1dep,set2dep) = some expression;
L’indice set2dep est optionnel.
L’exemple 63 peut ainsi se reformuler comme
SETS
i Variables / X1, X2, X3 /;

PARAMETERS
c(i) Costs / X1 109
X2 90
X3 115 /

a(i) Coeff / X1 6
X2 4
X3 8 /;

VARIABLES
Z Variable Z;

POSITIVE VARIABLES
x(i) Variables;

EQUATIONS
Equation1 Equation 1
Equation2 Equation 2
Equation3 Equation 3;

Equation1..
Z =E= sum(i, c(i)*x(i));

Equation2..
sum(i, x(i)) =L= 100;

Equation3..
sum(i, a(i)*x(i)) =L= 500;

MODEL Example1 /ALL/;

SOLVE Example1 USING LP MAXIMIZING Z;

8.6. NOTES 111

8.5.3 Autres logiciels

Les solutions logiciels retenues sont loin d’être les seules envisageables, et les choix sont volontairement orientés
afin de répondre aux objectifs d’un cours d’introduction, mais ne représentent pas forcément un choix optimal.
Il existe d’autres langages de modélisation intéressant (on pensera par exemple à AMPL), d’autres solveurs com-
merciaux ou gratuits. Nous nous désinteresserons cependant des solveurs intégrés dans les tableurs, les algorithmes
proposés par défaut dans ceux-ci étant en général de faible qualité, en particulier en optimisation non-linéaire, et
limités dans le nombre de variables qui peuvent être traités. Il est toutefois bon de noter que des outils existent
pour importer/exporter des données entre un tableur et des outils tels que GAMS.

8.6 Notes
La section sur GAMS s’inspire partiellement du cours Applied Mathematical Programming donné par Bruce A.
McCarl à l’Université Texas A&M.
112 ANNEXE
 Bibliographie

[1] Dimitri P. Bertsekas. Dynamic Programming and Optimal Control, volume I. Athena Scientific, Belmont, MA,
USA, 3rd édition, 2005.
[2] Frederick S. Hillier et Gerald J. Lieberman. Introduction to Operations Research. McGraw-Hill, New York,
USA, seventh édition, 2001.
[3] Maurice Kirby. Patrick maynard stuart blackett. International Transactions in Operational Research, 10(4):405–
407, 2003.
[4] Jeffrey C. Lagarias, James A. Reeds, Margaret H. Wright et Paul E. Wright. Convergence properties of the
Nelder-Mead simplex method in low dimensions. SIAM Journal on Optimization, 9(1):112–147, 1998.
[5] Adam B. Levy. The Basics of Practical Optimization. SIAM, Philadelphia, USA, 2009.
[6] Ken I. M. McKinnon. Convergence of the Nelder-Mead simplex method to a nonstationary point. SIAM Journal
on Optimization, 9(1):148–158, 1998.
[7] John A. Nelder et Roger Mead. A simplex method for function minimization. Computer Journal, 7:308–313,
1965.
[8] Jorge Nocedal et Stephen J. Wright. Numerical Optimization. Springer, New York, NY, USA, 1999.

113