CHAPITRE 1 : INTRODUCTION

Le terme �recherche opérationnelle�, d’origine militaire, n’est pas très suggestif, c’est le moins qu’on puisse dire. La recherche opérationnelle constitue autant une façon d’aborder un problème à l’aide d’outils mathématiques et informatiques qu’une discipline bien définie telle que l’anthropologie ou la linguistique. Ainsi on peut traiter par une approche de RO des problèmes de gestion, de politique électorale, de biologie moléculaire, d’ingénierie, de logistique, d’informatique (vision par ordinateur, infographie, télématique), etc. La RO est une discipline polyvalente, qui emprunte beaucoup aux mathématiques, évidemment, mais également à l’informatique et à l’économie. Elle le leur rend bien.

Dans le cadre de ce cours, nous allons étudier un certains nombres de modèles mathématiques pouvant s’appliquer à des situations concrètes. L’emphase sera mise sur la modélisation et les algorithmes de résolution, parfois au détriment, bien malheureusement, de la rigueur mathématique.

Considérons l’exemple d’un réseau routier urbain fort congestionné. Un candidat aux élections municipales propose, s’il est élu, d’améliorer l’écoulement du trafic par la construction d’une voie rapide. Au cours de l’élection, le débat porte sur la pertinence du projet. Il est admis par tous qu’une voie rapide ne peut que réduire le temps de transport des 6000 automobilistes qui se rendent de leur origine A à leur destination D (voir Figure 1). Par contre, le prix à payer pour la construction de la voie rapide est-il trop élevé pour la diminution de temps de transport appréhendée ?

[figure: 6000; 6000; A; B; C; D; voie; rapide; voie]

Fig. 1 – Le réseau du paradoxe de Braess

Afin de répondre de façon éclairée à cette question, il faut pouvoir évaluer l’impact du projet sur les chemins utilisés par les utilisateurs. Et pour connaître ces chemins, il faut d’abord savoir quels sont les critères de choix de ceux-ci. On suppose que les usagers ne tiennent compte que du temps de parcours et choisirons des chemins de temps de parcours minimum. Si t_{p} représente le temps de parcours du chemin p, t_{min} le temps de parcours minimum et x_{p} le nombre d’utilisateurs du chemin p, ce principe de comportement peut se traduire mathématiquement de la façon suivante :

t_{p} ≥ t_{min} ⟹ x_{p} = 0

ou encore t_{p} ≥ t_{min} p = 1, 2, 3 x_{p}(t_{p} − t_{min}) = 0 p = 1, 2, 3

∑_{p=1}^{3} x_{p} = 6000

x_{p} ≥ 0.

La première équation stipule qu’on ne peut avoir un chemin dont le délai serait inférieur au délai à l’équilibre. La seconde nous assure qu’il ne peut y avoir de flot sur des chemins dont le délai est supérieur au délai à l’équilibre. Les deux suivantes, satisfaction de la demande et non négativité des flots, se justifient d’elles-mêmes. Cette formulation mathématique est loin d’être triviale. Elle implique le bon choix de variables de décision (les flots de véhicules sur les chemins menant de A à D) ainsi que la traduction en équations d’une notion qui semble relever plus de la sociologie que de la mathématique. Ceci constitue la partie la plus ardue de l’exercice de modélisation. Une fois cette étape franchie, il ne reste plus qu’à raffiner un peu notre modèle d’affectation du trafic.

Spécifions d’abord la forme du temps de transport t_{p}. Supposons que ce temps soit égal à la somme des temps de transport sur les arcs formant le chemin p, ce qui paraît assez naturel. Ainsi :

t_{ACD} = t_{AC} + t_{CD}.

Puisque le réseau est congestionné, il est également naturel de supposer que les délais sur les arcs augmentent avec les flots sur ces arcs. On introduit, pour chaque arc du réseau, une fonction de congestion qui spécifie, pour chaque usager de l’arc considéré, le temps de parcours :

t_{AB}(x_{AB}) = \frac{1}{100} x_{AB}
t_{AC}(x_{AC}) = 50 + \frac{1}{1000}x_{AC}
t_{BC}(x_{BC}) = 10 + \frac{1}{1000}x_{BC}
t_{BD}(x_{BD}) = 50 + \frac{1}{1000}x_{BD}
t_{CD}(x_{CD}) = \frac{1}{100} x_{CD}.

Ces temps de parcours peuvent également s’exprimer en fonction des flots sur les chemins. Par exemple, puisque l’arc AB appartient à la fois aux chemins ABD et ABC, on obtient x_{AB} = x_{ABD} + x_{ABCD} et x_{BD} = x_{ABD}, d’où l’on dérive le délai du chemin ABD :

t_{ABD}(x_{ABD}, x_{ACD}, x_{ABCD}) = [\frac{1}{100} (x_{ABD} + x_{ABCD}) + (50 + \frac{1}{1000}x_{BD})].

On dérive de la même façon les délais des chemins ACD et ABCD.

On peut maintenant décrire formellement notre modèle d’affectation du trafic

t_{p} ≥ t_{min} ∀p ∈ {ABD, ACD, ABCD}
x_{ABD}(t_{ABD}(x_{ABD}, x_{ACD}, x_{ABCD}) − t_{min}) = 0
x_{ACD}(t_{ACD}(x_{ABD}, x_{ACD}, x_{ABCD}) − t_{min}) = 0
x_{ABCD}(t_{ABCD}(x_{ABD}, x_{ACD}, x_{ABCD}) − t_{min}) = 0
x_{ABD} + x_{ACD} + x_{ABCD} = 6000
x_{ABD}, x_{ACD}, x_{ABCD} ≥ 0.

Si les trois chemins sont utilisés, il suffit de résoudre le système linéaire

t_{ABD} = t_{ACD}
t_{ABD} = t_{ABCD}
x_{ABD} + x_{ACD} + x_{ABCD} = 6000,

dont l’unique solution est

x_{ABD} = x_{ACD} = x_{ABCD} = 2000.

On vérifie aisément que le délai commun est alors de 92 (minutes). Bien. Mais qu’en était-il du délai avant la construction de la voie rapide BC? Par symétrie, il est facile de voir qu’un équilibre sera atteint lorsque x_{ABD} = x_{ACD} = 3000 et que le délai correspondant est égal à 83 (minutes). La construction de la voie rapide aura donc eu comme conséquence, non seulement de faire élire notre conseiller municipal,^{1} mais également d’augmenter le temps de trajet de tous les automobilistes ! De tels paradoxes ont effectivement étés observés en pratique.

Cet exemple simpliste donne une petit idée de ce que la recherche opérationnelle peut accomplir. Le modèle précédent peut être raffiné pour tenir compte des caractéristiques suivantes d’un réseau urbain :

– la présence de plusieurs modes de transport (automobile, camion, autobus, bicyclette) ;

– les aspects aléatoires du trafic (incidents, accidents, conditions météorologiques) ;

– les variations horaires, journalières et saisonnières de la demande ;

– les modes de contrôle du trafic (feux de circulation, panneaux à messages variables, systèmes d’information) ;

– le caractère dynamique du trafic ;

– la présence de différentes classes d’usagers, chacune avec ses propres critères et comportements.

Plus un modèle est raffiné, plus il risque de se rapprocher d’une situation réelle. Par contre, plus il risque d’être difficile à résoudre. De plus, un problème commun aux modèles détaillés est celui de la collecte de données fiables. Un modèle trop complexe risque d’être lourd et inutilisable.

1Dans une élection subséquente, celui-ci aura bien entendu rejeté le tort sur les ingénieurs responsables du projet.

Les modèles de RO s’expriment souvent (pas toujours, loin de là) comme des problèmes d’optimisation. Dans l’espace euclidien à n dimensions, on obtient typiquement la formulation^{2}

\begin{aligned}max & f(x_{1}, …, x_{n}) \\ & g_{1}(x_{1}, …, x_{n}) & ≤ & 0 \\ & ⋮ & ⋮ & ⋮ \\ & g_{m}(x_{1}, …, x_{n}) & ≤ & 0 \end{aligned}

où x = (x_{1}, …, x_{n}) est le vecteur des variables, f la fonction objectif et g_{1}, …, g_{m} les contraintes du problème. Puisqu’il est inutile de produire une solution qui ne respecte pas les contraintes et que les contraintes s’expriment en fonction des variables de décision, les étapes du processus de construction d’un modèle doivent être, dans l’ordre :

– le choix des variables de décision ;

– la représentation du domaine réalisable (faisabilité d’un vecteur x) à l’aide des fonctions de contrainte g_{i};

– le choix d’un objectif.

Une fois le modèle posé, eh bien, il faut le résoudre. Et ici, on quitte presque le domaine de la RO pour entrer dans celui des mathématiques (programmation mathématique, théorie des graphes, analyse numérique, processus stochastiques) et de l’informatique (structures de données, interfaces usager-machine, programmation). Il n’y a pas de limite aux techniques que peut piller la RO. Le but est de résoudre un problème. La technique de résolution est secondaire.

Dans les chapitres qui suivent, nous étudierons un certain nombre de modèles ainsi que de techniques de résolution de ces modèles. Ces modèles et techniques doivent être perçus comme une boîte à outils. Il ne faut pas hésiter à enrichir celle-ci au besoin.

Maintenant, place aux modèles !

2Noter qu’une contrainte d’égalité g(x) = 0 peut être représentée par les deux contraintes d’inégalité g(x) ≤ 0 et −g(x) ≤ 0.

CHAPITRE 2 : PROGRAMMATION LINÉAIRE

À tout seigneur tout honneur ! La programmation linéaire constitue la pierre angulaire de toute la recherche opérationnelle. Il faut bien sûr éviter de forcer tout modèle à être linéaire. Par contre, 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 s’exprime de façon générique sous la forme

max c_{1}x_{1} + ⋯ + c_{n}x_{n}  objectif
a_{11}x_{1} + ⋯ + a_{1n}x_{n} ≤ b_{1}
⋮  ⋮  contraintes
a_{m1}x_{1} + ⋯ + a_{mn}x_{n} ≤ b_{m}

ou, sous une forme plus compacte :

max ∑_{j=1}^{n} c_{j}x_{j}

∑_{j=1}^{n} a_{ij}x_{j} ≤ b_{i}  i = 1, …, m

et, pour adultes vaccinés avec droit de vote, la forme matricielle :

c : vecteur ligne 1 × n
max cx  x : vecteur colonne n × 1
Ax ≤ b  A : matrice m × n
b : vecteur colonne m × 1

Exemple (voir Figure 1) :

max 3x_{1} + 5x_{2}
x_{1}  ≤ 4
2x_{2} ≤ 12
3x_{1} + 2x_{2} ≤ 18
x_{1} ≥ 0 (−x_{1} ≤ −0)
x_{2} ≥ 0  (−x_{2} ≤ −0)

Un point réalisable (ou admissible) est extrémal (on parle parfois d’un sommet) s’il ne peut se trouver à l’intérieur d’un segment (non nul) entièrement contenu dans la région admissible (polyèdre). Dans l’exemple de la figure 1 il y a cinq points extrémaux : x^{1} = (0, 0), x^{2} = (0, 6), x^{3} = (2, 6), x^{4} = (4, 3) et x^{5} = (4, 0).

Le but du jeu est d’augmenter la valeur de l’objectif. L’observation fondamentale est que l’objectif augmente tout en satisfaisant aux contraintes si et seulement si l’on se dirige dans une direction admissible^{1} d qui fait un angle aigu avec le vecteur c, ce qui advient si et seulement si le produit scalaire des vecteurs c et d est positif. Dans la figure 2, la direction d^{1} est une direction d’amélioration, d^{2} une direction de détérioration et d^{3} une direction de constance. Il est important de noter que

– la droite D est le lieu géométrique des points pour lesquels l’objectif vaut c_{1}x_{1} + c_{2}x_{2}. Ceci correspond à une courbe de niveau (ou droite de constance) sur une carte topographique ;

– les courbes de niveau sont parallèles ;

1Une direction est admissible s’il est possible de s’y engager sans quitter le domaine réalisable. Il se peut que l’extrémité du vecteur définissant cette direction se trouve à l’extérieur du domaine réalisable.

[figure: 1; 2; 3; 4; 5; 6; 1; 2; 3; 4; 5; 6; x_{1}; x_{2}; x^{1}; x^{2}; x^{3}; x^{4}; x^{6}; x^{5}; Région; admissible; (polyèdre); Ax ≤ b; vecteur; (libre); c = (3, 5); 3x_{1} + 2x_{2} = 18]

Fig. 1 – La géométrie d’un programme linéaire

– le vecteur c est le vecteur normal (terme synonyme de perpendiculaire ou orthogonal) associé à la droite D. À chaque point de l’espace correspond une courbe de niveau, c’est-à-dire une droite de constance dans le cas qui nous intéresse. Ainsi la droite D_{1} est une courbe de niveau correspondant à une valeur supérieure de l’objectif ;

– pour un déplacement de longueur unitaire, la plus grande amélioration est obtenue en se dirigeant dans la direction du vecteur c. Cette direction est appelée direction de la plus grande pente et est orthogonale aux droites de constance.

[figure: D; c = (c_{1}, c_{2}); d^{1}; d^{2}; x = (x_{1}, x_{2}); D_{1}; d^{3}]

Fig. 2 – Directions

En suivant la direction admissible de plus grande amélioration (cette direction est unique) et en s’interdisant de quitter la région admissible, on atteint une solution optimale en un nombre fini d’itérations (voir la trajectoire ABCD illustrée à la figure 3). Une solution est optimale si et seulement s’il n’y a pas de direction d’amélioration admissible en ce point.

Pour tout programme linéaire, l’une ou l’autre des quatre situations suivantes doit se présenter :

– il existe une solution optimale et celle-ci est unique ;

– l’objectif est borné et il existe une infinité de solutions optimales qui sont situées sur une face non triviale du polyèdre des contraintes. Dans l’exemple de la figure 1 il existe cinq faces triviales de dimension nulle (points extrémaux), cinq faces de dimension un (côtés du pentagone) et une face de dimension deux (le pentagone lui-même). Dans l’espace, on peut rencontrer des faces de dimensions zéro (points extrémaux), un (arêtes), deux (faces au sens géométrique traditionnel) et trois (le polyèdre lui-même) ;

– il n’existe aucune solution réalisable ;

– le problème est non borné, c’est-à-dire que l’objectif peut augmenter indéfiniment en suivant une demidroite faisant un angle aigu avec le vecteur objectif et entièrement située dans le domaine admissible. Il est clair que cette situation ne peut se produire que si le domaine lui-même n’est pas borné.

[figure: 1; 2; 3; 4; 5; 6; 1; 2; 3; 4; 5; 6; x_{1}; x_{2}; A; B; C; D; c; c; c]

Fig. 3 – Trajectoire menant à une solution optimale

Il n’est pas difficile de démontrer que s’il existe une solution optimale bornée, il existe alors une solution optimale extrémale. Pour résoudre le problème, il suffirait donc d’énumérer les points extrémaux (sommets) du domaine, qui sont en nombre fini. Las ! Ils sont en nombre exponentiel en fonction de la taille (nombre de contraintes) du problème.

Le premier algorithme proposé pour la programmation linéaire visite de façon �efficace� des points extrémaux adjacents, tout en améliorant l’objectif. Dans l’exemple précédent, deux suites partant du point (0, 0) et satisfaisant ces deux critères sont possibles : {x^{1}, x^{2}, x^{3}} et {x^{1}, x^{5}, x^{4}, x^{3}}.

Cet algorithme ressemble singulièrement à l’algorithme d’élimination de Gauss-Jordan pour la résolution de systèmes d’équations linéaires sous-déterminés, c’est-à-dire contenant plus de variables que de contraintes et admettant, en général, une infinité de solutions réalisables.

Appliquons l’algorithme du simplexe à notre exemple. L’algorithme utilise une forme canonique (�standard�) où les variables sont toutes positives ou nulles et les autres contraintes sont toutes des contraintes d’égalité. En ajoutant des variables d’écart (x_{3}, x_{4} et x_{5} dans notre exemple) aux contraintes d’inégalité “≤” et en remplaçant toute variable libre x_{i} par la différence x_{i}^{+} − x_{i}^{−} de deux variables non négatives, il est facile de se ramener à cette forme.^{2}

max 3x_{1} + 5x_{2}
x_{1}  + x_{3}  = 4
2x_{2}  + x_{4}  = 12
3x_{1} + 2x_{2}  + x_{5} = 18
x_{1},  x_{2},  x_{3},  x_{4},  x_{5} ≥ 0

Comme dans la phase finale de la méthode de Gauss-Jordan, on exprime des variables liées (variables de base) en fonction de variables libres (variables hors-base), ce qui nous donne le dictionnaire suivant :

\begin{aligned}max & 3x_{1} & + & 5x_{2} \\ & x_{3} & = & 4 & − & x_{1} \\ & x_{4} & = & 12 & & & − & 2x_{2} \\ & x_{5} & = & 18 & − & 3x_{1} & − & 2x_{2} \\ & x_{1}, & x_{2}, & x_{3}, & x_{4}, & x_{5} & ≥ & 0 \end{aligned}

Une solution de base est obtenue en posant x_{j} = 0 pour toutes les variables hors-base, c’est-à-dire les variables qui se trouvent à droite des signes d’égalité. Tout dictionnaire correspond à un point extrémal.^{3} Le dictionnaire précédent correspond au point extrémal (x_{1}, x_{2}, x_{3}, x_{4}, x_{5}) = (0, 0, 4, 12, 18). En deux dimensions, ceci correspond au point x^{1} = (0, 0).

L’algorithme du simplexe substitue à une variable de base (à gauche du dictionnaire) une variable hors-base (à droite du dictionnaire). On choisit évidemment une variable qui fasse augmenter l’objectif. On a le choix entre les variables x_{1} et x_{2}. Choisissons la variable x_{1}. Afin de maintenir la non négativité des variables de base, x_{1} ne peut dépasser la valeur x_{1} = 4. Une fois cette valeur atteinte, la variable de base x_{3} devient nulle et doit quitter la base. En posant x_{1} = 4 − x_{3} (première équation du dictionnaire), on obtient le nouveau dictionnaire :

x_{1} ↑ 4 x_{3} ↓ 0
x_{1} = 4 − x_{3}
max 3(4 − x_{3}) + 5x_{2} = 12 + 5x_{2} − 3x_{3}
x_{1} = 4 − x_{3}  = 4 − x_{3}
x_{4} = 12 − 2x_{2}  = 12 − 2x_{2}
x_{5} = 18 − 3(4 − x_{3}) − 2x_{2} = 6 − 2x_{2} + 3x_{3}

qui correspond au point extrémal x^{5} = (4, 0, 0, 12, 6). À la deuxième itération, x_{2} entre dans la base et x_{5} la quitte.

2N’ayons pas peur d’ajouter des variables. Je rappelle que nous traitons de gros problèmes. Enfin, il faut s’entendre : le terme �nous� désigne un ordinateur qui a été efficacement programmé.

3Cependant, un point extrémal peut avoir plusieurs représentations comme dictionnaire. Noter qu’un dictionnaire comportant des variables de base négatives ne correspond pas à un point extrémal ; voir par exemple le point x^{6} de la figure 1.

x_{2} ↑ 3 x_{5} ↓ 0
x_{2} = 3 + \frac{3}{2}x_{3} − \frac{1}{2}x_{5}
max  27 + \frac{9}{2}x_{3} − \frac{5}{2}x_{5}
x_{2} = 3 + \frac{3}{2}x_{3} − \frac{1}{2}x_{5}
x_{1} = 4 − x_{3}
x_{4} = 6 − 3x_{3} + x_{5}

Ce dictionnaire correspond au point extrémal x^{4} = (4, 3, 0, 6, 0).

On peut également adopter un format de tableau qui évite d’écrire les noms des variables. Plusieurs formats de tableau sont possibles. Vous trouverez ci-dessous le dictionnaire précédent sous forme de tableau ainsi que le tableau optimal obtenu en effectuant une troisième itération où la variable x_{3} réintègre la base. L’élément encerclé du tableau de gauche, situé à l’intersection de la quatrième ligne (variable quittant la base) et de la deuxième colonne (variable entrant dans la base) est appelé pivot et joue un rôle important dans la dérivation de règles arithmétiques permettant de passer d’un tableau à un autre.

[figure: 6; −3; 1; 4; −1; 0; 3; \frac{3}{2}; −\frac{1}{2}; 27; \frac{9}{2}; −\frac{5}{2}; x_{3}; x_{5}; x_{2}; x_{1}; x_{4}; 2; −\frac{1}{3}; \frac{1}{3}; 2; \frac{1}{3}; −\frac{1}{3}; 6; −\frac{1}{2}; 0; 36; −\frac{3}{2}; −1; x_{4}; x_{5}; x_{2}; x_{1}; x_{3}; objectif; objectif; variables hors-base; variables; de base; optimum : x^{3}]

Les coefficients de la première ligne se trouvant sous les variables hors-base sont appelés coûts réduits. Lorsque ceux-ci sont tous négatifs ou nuls, la solution correspondante est bien sûr optimale. Attention : La contraposée de cette affirmation n’est pas nécessairement vraie ! En effet, une solution dont certains coûts réduits sont positifs peut être optimale.

DIVERSES CONSIDÉRATIONS

Comme vous avez pu le constater (enfin je l’espère) à la lecture de ces quelques pages, la programmation linéaire, ce n’est pas vraiment très compliqué. Il y a cependant un certain nombre de points à éclaircir.

1. Initialisation

(a) Une contrainte d’inégalité ax ≥ b peut être mise sous la forme −ax ≤ −b;

(b) Une contrainte d’égalité ax = b peut être remplacé par deux contraintes d’inégalité ax ≤ b et ax ≥ b;

(c) Une variable x_{j} qui n’est pas contrainte en signe peut être remplacée par la différence x_{j} = x′_{j}−x′′_{j} de deux variables non négatives.

Maintenant, que faire si la base correspondant aux variables d’écart n’est pas admissible ? Par exemple soit les contraintes

\begin{aligned}x_{1} & + & x_{2} & ≤ & 2 \\ x_{1} & − & x_{2} & ≥ & 1 \end{aligned}

dont la forme canonique est, après ajout de deux variables d’écart :

\begin{aligned}x_{1} & + & x_{2} & + & x_{3} & & & = & 2 \\ −x_{1} & + & x_{2} & & & + & x_{4} & = & −1. \end{aligned}

Malheureusement, la base constituée des variables d’écart n’est pas réalisable. En changeant le signe de la deuxième contrainte et en ajoutant une variable artificielle (fictive) x_{5} on obtient le système linéaire

\begin{aligned}x_{1} & + & x_{2} & + & x_{3} & & & & & = & 2 \\ x_{1} & − & x_{2} & & & − & x_{4} & + & x_{5} & = & 1 \end{aligned}

qui possède la base admissible x = (0, 0, 2, 0, 1). Une solution x pour le problème modifié sera admissible pour le problème originel si et seulement si la variable artificielle x_{5} est nulle. Il est naturel d’obtenir une telle solution par la résolution du programme linéaire

min  x_{5}
x_{1} + x_{2} + x_{3}  = 2
x_{1} − x_{2}  − x_{4} + x_{5} = 1
x_{1}, x_{2}, x_{3}, x_{4}, x_{5} ≥ 0.

Si la valeur optimale de ce programme est positive, alors le polyèdre originel est vide. Sinon, on élimine les colonnes correspondant aux variables artificielles^{4}, on exprime l’objectif originel en fonction des variables hors-base et l’on entame la �phase II� de la méthode.

Cette technique se généralise à des situations impliquant plusieurs variables artificielles. Il suffit alors de minimiser la somme des variables artificielles.

2. Convergence Toute solution de base réalisable correspond à un point extrémal. Si toutes les solutions de base ont des valeurs (objectif) distinctes, alors l’algorithme du simplexe ne peut visiter une base plus d’une fois et doit donc converger en améliorant l’objectif à chaque itération. Notons que le nombre de solutions de base (réalisables ou non) est énorme : \binom{m}{n} = \frac{m!}{n!(m−n)!}. En pratique cependant le nombre de bases visitées semble linéaire par rapport au nombre de contraintes. Un point extrémal peut admettre plusieurs représentation de base. Il se pourrait alors que que l’on passe d’une base à l’autre sans améliorer l’objectif. Ceci ne peut se produire que si au moins une des variables de base est nulle (une telle solution de base est dite dégénérée). Le danger est alors de visiter de façon cyclique un ensemble de bases dégénérées, ce qui est possible, mais peu probable en pratique. Noter qu’il existe des règles permettant théoriquement d’éviter le phénomène de cyclage.

On trouve à la figure 4 un exemple de point extrémal dégénéré (le point x = (1, 1)) et le programme linéaire correspondant. Soient x_{3}, x_{4} et x_{5} les variables d’écart du problème. Le dictionnaire associé au point (1, 1) est

max  −3 + x_{3} + x_{5}
x_{1} = 1 − x_{3}
x_{2} = 1 + x_{3} − x_{5}
x_{4} = 0 − x_{3} + x_{5} .

4Des ajustements techniques sont nécessaires si une variable artificielle de valeur nulle se trouve dans la base.

Si on introduit la variable x_{3} dans la base, celle-ci prend la valeur zéro. Il y a changement de base (mais pas de point extrémal). À l’itération suivante, il y aura (vérifiez-le) changement de point extrémal et amélioration de l’objectif.

Si on perturbe l’une des contraintes satisfaites à égalité par le point dégénéré, la dégénérescence disparaît. Par exemple, si la contrainte x_{1} + x_{2} ≤ 1 est remplacée par la contrainte x_{1} + x_{2} ≤ 1 − ϵ, soit le point dégénéré (1, 1) est remplacé par les deux points extrémaux (1, 1 − ϵ) et (1 − ϵ, 1) si ϵ est positif (voir figure 4), soit la contrainte ne joue plus aucun rôle dans le problème ainsi que dans l’algorithme. Cette technique de perturbation du vecteur b peut être utilisée en théorie ou en pratique pour éliminer la dégénérescence.

[figure: x_{1}; x_{2}; 1; 1; (1,1); x_{1} = 1; x_{2} = 1; x_{1} + x_{2} = 2; direction de c; max; −2x_{1} − x_{2}; x_{1} ≤ 1; x_{2} ≤ 1; x_{1} + x_{2} ≤ 2; x_{1}, x_{2} ≥ 0; contrainte perturbée]

Fig. 4 – Un exemple de dégénérescence

3. Efficacité

– L’algorithme du simplexe, bien que de complexité théorique exponentielle, est très efficace en pratique ;

– Il existe des algorithmes polynomiaux pour la programmation linéaire, basés sur des principes entièrement différents de ceux de la méthode du simplexe. Entre autres, ces algorithmes �efficaces� ne produisent une solution extrémale qu’à la dernière itération.

4. Implantation numérique

Les problèmes pratiques impliquent fréquemment des milliers, voire des millions de variables et de contraintes. Une implantation naïve peut se révéler désastreuse, suite à des problèmes numériques. Il est en effet important de savoir si une variable, un coût réduit, est positif, nul ou négatif. Les logiciels de programmation linéaire utilisent des techniques sophistiquées pour choisir les variables d’entrée et de sortie ainsi que pour effectuer les changements de base, une opération qui semble pourtant triviale sur de petits exemples.

5. Solution non bornée

Si une variable hors-base de coût réduit positif n’est pas limitée dans son augmentation (ce qui arrive si tous les coefficients de la colonne du dictionnaire correspondant à cette variable sont positifs ou nuls) alors l’objectif peut augmenter indéfiniment le long d’une demi-droite dont l’extrémité est le point extrémal courant. L’algorithme du simplexe reconnaît cette situation et envoie un message approprié à l’utilisateur.

[figure: 2x_{2} = 12; 3x_{1} + 2x_{2} = 18; (0,2); (3,5); (3,2); d = (d_{1}, d_{2}); x^{3} = (2, 6); (3, 5) = 1 × (3, 2) + \frac{3}{2} × (0, 2)]

Fig. 5 – La dualité

6. Dualité

À l’optimum, il n’existe pas de direction d’amélioration admissible, c’est-à-dire :

(3, 5)(d_{1}, d_{2}) ≤ 0 ∀ d = (d_{1}, d_{2}) admissible.

Un théorème stipule (voir l’illustration en Figure 5) que cette condition est satisfaite si et seulement si le vecteur objectif c = (3, 5) peut s’exprimer comme une combinaison linéaire non-négative de vecteurs normaux (orthogonaux) aux contraintes actives. Dans notre cas, on a :

(3, 5) = y_{2}(0, 2) + y_{3}(3, 2),

d’où l’on déduit : y_{2} = 3/2 et y_{3} = 1. En posant y_{1} = y_{4} = y_{5} = 0, on obtient

(3, 5) = y_{1}(1, 0) + y_{2}(0, 2) + y_{3}(3, 2) + y_{4}(−1, 0) + y_{5}(0, −1),

avec y_{2}, y_{3} ≥ 0 et y_{1} = y_{4} = y_{5} = 0 ou, matriciellement :

(3, 5) = ( y_{1} y_{2} y_{3} y_{4} y_{5} ) ( \begin{matrix}1 & 0 \\ 0 & 2 \\ 3 & 2 \\ −1 & 0 \\ 0 & −1 \end{matrix} ) .

Dans le cas général, on obtient (A^{i} désigne la ligne i de la matrice A) :

c = yA  y ≥ 0  et y_{i} = 0 si A^{i}x < b_{i},

ce qu’on peut exprimer synthétiquement comme :

yA = c
y ≥ 0
y(b − Ax) = 0.

Le vecteur de multiplicateurs y est aussi appelé vecteur dual. Un vecteur dual est admissible s’il satisfait aux contraintes

yA = c y ≥ 0.

Un vecteur primal x est alors optimal si et seulement s’il existe un vecteur dual y tel que

Ax ≤ b  primal réalisabilité
yA = c  dual
y ≥ 0  réalisabilité
y(b − Ax) = 0  complémentarité.

Puisque yA = c, la condition de complémentarité peut également s’exprimer comme yb − cx = 0 ou yb = cx. Noter que yb ≥ cx et que l’égalité ne peut être réalisée qu’à l’optimum. Pour tout vecteur dual-réalisable y, la quantité yb représente donc une bonne supérieure sur la valeur optimale de l’objectif (primal) cx. On obtient la meilleure borne en résolvant le programme linéaire dual :

\begin{aligned}min & yb \\ & yA = c \\ & y ≥ 0. \end{aligned}

Puisque le problème dual est un programme linéaire, on peut le mettre sous forme canonique :

max −yb  variables duales du dual
yA ≤ c  u′
−yA ≤ −c  u′′
−y ≤ 0  v

et calculer son programme dual

\begin{aligned}min & cu′ − cu′′ \\ & Au′ − Au′′ − v = −b \\ & u′, u′′, v ≥ 0 \end{aligned}  ou  \begin{aligned}max & c(u′′ − u′) \\ & A(u′′ − u′) + v = b \\ & u′, u′′, v ≥ 0. \end{aligned}

En effectuant la substitution x = u′′ − u′ on élimine les contraintes de signe sur les variables u′ et u′′, ce qui donne :

max  cx
Ax + v = b
v ≥ 0,

qui est équivalent (remarquer que v est une variable d’écart) au programme primal

max  cx
Ax ≤ b.

À partir d’un dictionnaire optimal, il est possible de recouvrer la solution optimale du problème dual. Cette opération est effectuée par tous les logiciels de programmation linéaire.

Puisque le programme dual est un programme linéaire, on peut le résoudre comme le problème primal. Il est parfois avantageux de résoudre le dual plutôt que le primal, par exemple si le dual comporte moins de variables.

7. Le diagramme ci-dessous représente les diverses possiblités pour un couple primal-dual de programmes linéaires.

×  × ×  × fini −∞ ∞ fini ∞ −∞  dual primal  ∞ : primal non borné ou dual non réalisable −∞ : primal non réalisable ou dual non borné

8. Formes non canoniques

Si le programme linéaire n’est pas sous forme canonique, les règles suivantes s’appliquent :

– à une contrainte d’inégalité dans le primal correspond une variable duale non négative ;

– à une contrainte d’égalité dans le primal correspond une variable duale libre ;

– à une variable libre dans le primal correspond une contrainte d’égalité dans le dual ;

– à une variable non négative dans le primal correspond une contrainte d’inégalité dans le dual ;

– les signes des contraintes principales (pas les contraintes x ≥ 0) sont inversés ;

– il y a complémentarité entre les contraintes d’inégalité du primal et les variables duales correspondantes, et vice versa. Ainsi, pour la forme utilisée dans l’algorithme du simplexe, on obtient

primal  dual  complémentarité
max cx  min yb  (c − yA)x = 0
Ax = b  yA ≥ c
x ≥ 0

L’exemple suivant illustre tous les cas de figure.

primal  dual  complémentarité
min 2x_{1} − 3x_{2} + x_{3} − 6x_{4}  max y_{1} + 5y_{2}
x_{1} + x_{2} + x_{3} + x_{4}  = 1  y_{1} + y_{2} = 2  y_{2}(x_{1} − 3x_{2} + 2x_{4} − 5) = 0
x_{1} − 3x_{2} + 2x_{4}  ≥ 5  y_{1} − 3y_{2} ≤ −3  x_{2}(−3 − y_{1} + 3y_{2}) = 0
y_{1}  ≤ 1  x_{3}(1 − y_{1}) = 0
x_{2}, x_{3} ≥ 0  y_{1} + 2y_{2} = −6
x_{1}, x_{4} libres
y_{2} ≥ 0
y_{1} libre

9. Analyse paramétrique

Soit x une solution optimale d’un programme linéaire. Si l’on perturbe légèrement le vecteur c (c devient c′) et que la base optimale n’en est pas modifiée, la valeur optimale du programme perturbé est c′x. Même si la solution x cesse d’être optimale, le nombre c′x nous donne une borne inférieure sur la valeur optimale du problème perturbé.

Symétriquement, si le vecteur de ressources b est perturbé (b devient b′) et que la base demeure réalisable, l’objectif cx = yb devient yb′, où y est une solution optimale du problème dual. Si la base cesse d’être réalisable, alors yb′ fournit une borne supérieure sur la valeur du nouvel objectif. On réalise ici l’importance du vecteur dual-optimal y. La composante y_{i} représente la valeur marginale d’une unité supplémentaire de la iième ressource.

Pour toutes ces raisons, il est important de connaître l’information duale et de l’intégrer au dictionnaire final, pour former un dictionnaire final étendu (voir figure 6). Noter que la numérotation des variables duales doit correspondre à l’ordre des contraintes. Dans notre exemple, l’ordre des contraintes avait été modifié à la deuxième itération, ce qui explique la numérotation des variables duales. Noter que ce n’est pas le fruit du hasard si les composantes positives du vecteur dual optimal sont égales aux coûts réduits des variables hors-base, au signe près, que l’on trouve dans la première ligne du tableau final.

2  −\frac{1}{3}  ^{\frac{1}{3}} 2  ^{\frac{1}{3}}  −\frac{1}{3} 6  −\frac{1}{2}  0 36  −\frac{3}{2}  −1 x_{4}  x_{5} x_{2} x_{1} x_{3} max  optimum : x^{3}  1 0^{\frac{3}{2}}  y_{3} y_{1} y_{2}

Fig. 6 – Un dictionnaire étendu

10. Logiciels

Il existe de nombreux logiciels pour résoudre des problèmes linéaires. Pour vous amuser, vous pouvez tester l’une des nombreuses applets Java, disponibles sur Internet et permettant de résoudre de petits programmes linéaires. Plus sérieusement, vous pouvez utiliser des logiciels spécialisés tels que LINDO, CPLEX ou GAMS. Des versions étudiantes gratuites de GAMS sont disponibles sur Internet.

Lorsque les problèmes deviennent grands, il est nécessaire d’utiliser une syntaxe évoluée. Lorsqu’ils deviennent très grands, il est indispensable d’utiliser un langage spécialisé qui construit, pour une application spécifique, la matrice des contraintes selon la syntaxe du logiciel et les spécifications de l’usager. Ces langages sont appelés �générateurs de matrices�. Enfin, je mentionne que les tableurs les plus populaires, tels que Excel, permettent de résoudre des problèmes linéaires. Avis aux amateurs des produits Microsoft.

11. Forme matricielle

Nous terminons ce chapitre par un exposé de certains résultats basés sur une notation matriciellevectorielle. Pour changer le mal de place, nous considérons un objectif à minimiser plutôt qu’à maximiser.

\begin{aligned}min_{x} & cx \\ & Ax = b & c : & 1 × n \\ & x ≥ 0 & x : & n × 1 \\ & & A : & m × n \\ & & b : & m × 1 \end{aligned} (P)

On suppose que la matrice A (dont la jième colonne sera notée A_{j}) est de plein rang. Soit B une sousmatrice carrée (base) formée de m colonnes indépendantes de A. On associe à B les décompositions

A = [B|N] et x = ( \begin{matrix}x_{B} \\ x_{N} \end{matrix} ). Si l’on fixe le vecteur des variables hors base x_{N} à zéro, on obtient

la solution de base x = ( \begin{matrix}B^{−1}b \\ 0 \end{matrix} ). Une solution de base est admissible si x_{B} ≥ 0. Réécrivons le

programme linéaire (P) en mettant en évidence la décomposition de A et de x :

min_{x_{B}}  c_{B}x_{B}
Bx_{B} + Nx_{N} = b
x_{B}, x_{N} ≥ 0

ou, après quelques manipulations élémentaires :

min_{x_{N}}  c_{B}B^{−1}b + (c_{N} − c_{B}B^{−1}N)x_{N}
B^{−1}b − B^{−1}Nx_{N} ≥ 0.

Puisque x_{B} = B^{−1}b ≥ 0, la solution x_{N} = 0 du programme ci-dessus est réalisable. L’algorithme du simplexe tente d’améliorer la solution en intégrant dans la base une variable hors-base x_{j^{∗}} (j^{∗} ∈ N) dont le coût réduit c_{j^{∗}} − c_{B}B^{−1}A_{j^{∗}} est négatif. La valeur maximale que peut prendre x_{j^{∗}} est donnée par

x_{j^{∗}} = arg inf_{t}{t|(B^{−1}b)_{i} − t(B^{−1}A_{j^{∗}})_{i} ≥ 0 ∀i ∈ B}.

Si t = +∞, le programme linéaire (P) est non borné inférieurement. Sinon, on fait sortir de la base une variable x_{i^{∗}} (i^{∗} ∈ B) telle que

(B^{−1}b)_{i^{∗}} − x_{j^{∗}}(B^{−1}A_{j^{∗}})_{i^{∗}} = 0.

La nouvelle base est : B ← B ∪ A_{j^{∗}} − A_{i^{∗}}. On continue le processus jusqu’à ce que le vecteur des coûts réduits c − c_{B}B^{−1}A soit non négatif. [remarque : la partie de ce vecteur correspondant aux variables de base est nulle puisque c_{B} − c_{B}B^{−1}B = 0.] Dualité : Soit y un vecteur 1 × m. Une borne inférieure sur la valeur du programme (P) s’obtient en substituant aux contraintes Ax = b l’unique contrainte yAx = yb. Si de plus on a que yA ≤ c, on obtient trivialement la borne inférieure cx ≥ yAx = yb. Si l’on veut maintenant obtenir la meilleure borne inférieure à partir de cette technique, on est amené à étudier le programme linéaire dual

\begin{aligned}max_{y} & yb \\ & yA ≤ c \end{aligned} (D)

Remarquons que le dual du programme linéaire (D) est (P). De plus, si le problème primal (P) possède un minimum, on a le résultat :

max_{y} yb = min_{x} cx.

Les conditions primales-duales suivantes caractérisent l’optimalité d’une solution :

\begin{matrix}Ax = b \\ x ≥ 0 \end{matrix} \} primal réalisable

yA ≤ c } dual réalisable

\begin{matrix}yb = cx \\ ou (yA − c)x = 0 \end{matrix} \}

écarts complémentaires (orthogonalité)

Les variables y sont appelées variables duales ou parfois multiplicateurs du simplexe. Si B est une base optimale, c’est-à-dire que les coûts réduits sont tous positifs ou nuls, alors y = c_{B}B^{−1} est une solution optimale du problème dual. Finalement, on dit qu’une base (ou la solution de base correpondante) est dégénérée si au moins une des variables de base est nulle.

CHAPITRE 3 : PROGRAMMATION NON LINÉAIRE

Un programme non linéaire prend la forme générique

\begin{aligned}max & f(x_{1}, …, x_{n}) \\ & g_{1}(x_{1}, …, x_{n}) & ≤ & b_{1} \\ & ⋮ & ⋮ & ⋮ \\ & g_{m}(x_{1}, …, x_{n}) & ≤ & b_{m} \end{aligned}

où au moins une des fonctions f et g_{1}, …, g_{m} est non linéaire.

L’équivalent de la notion de droite de constance en programmation linéaire est la notion de courbe de niveau qui constitue le lieu géométrique des points admettant une valeur commune de l’objectif. Un résultat fondamental est que le gradient de la fonction f au point x = (x_{1}, …, x_{n}), c’est-à-dire le vecteur ligne de ses dérivées partielles

∇f(x) = (\frac{∂f}{∂x_{1}} , …, \frac{∂f}{∂x_{n}} )

est orthogonal à la courbe de niveau α passant par le point x (voir Figure 1), d’équation

{y : f(y) = f(x) = α}.

[figure: x_{1}; x_{2}; x; ∇f(x); 1; 2; 3; 4; ∇f(x); 1; 2; 3; 4; x]

Fig. 1 – Exemples de courbes de niveau

La direction du gradient est la direction de plus grande pente (ascendante). Une direction d est une direction de montée au point x si le produit scalaire ∇f(x)d est positif, et une direction de descente si ∇f(x)d est

[figure: x_{1}; f(x_{1}); max local; max global; min local]

Fig. 2 – Optima locaux et globaux

négatif. L’objet mathématique ∇f(x)d est la dérivée directionnelle de la fonction f dans la direction d au point x.

La géométrie d’un programme non linéaire est beaucoup plus complexe que celle d’un programme linéaire. En particulier :

– la solution n’est plus forcément extrémale ; celle-ci peut se situer à l’intérieur du domaine réalisable, là où aucune contrainte n’est satisfaite à égalité ;

– il peut exister des optima locaux (voir Figure 2). Une condition nécessaire (mais pas nécessairement suffisante) pour qu’un point x^{∗} soit optimal est qu’il n’existe pas de direction d’amélioration admissible en x^{∗} (voir Figure 3) :

∇f(x^{∗})d ≤ 0 ∀d admissible.

Sous certaines hypothèses techniques, une condition nécessaire pour que cette condition soit satisfaite, et donc que x^{∗} soit un candidat à l’optimalité, est que x^{∗} soit une solution optimale du programme linéaire obtenu en linéarisant l’objectif et les contraintes du programme non linéaire originel :

max_{x} f(x^{∗}) + ∇f(x^{∗})(x − x^{∗})
g_{1}(x^{∗}) + ∇g_{1}(x^{∗})(x − x^{∗}) ≤ b_{1}
⋮  ⋮  ⋮
g_{m}(x^{∗}) + ∇g_{m}(x^{∗})(x − x^{∗}) ≤ b_{m}

où, sous forme compacte et canonique :

\begin{aligned}max_{x} & ∇f(x^{∗})x \\ & ∇g_{i}(x^{∗})x & ≤ & b_{i} − g_{i}(x^{∗}) + ∇g_{i}(x^{∗})x^{∗} & i = 1, …, m \end{aligned}

Les conditions d’optimalité de ce programme linéaire sont

∇g_{i}(x^{∗})x ≤ b_{i} − g_{i}(x^{∗}) + ∇g_{i}(x^{∗})x^{∗} i = 1, …, m • primal réalisable
∑_{i=1}^{m} y_{i}∇g_{i}(x^{∗}) = ∇f(x^{∗})  • dual
y_{i} ≥ 0  i = 1, …, m  réalisable
y_{i}(b_{i} − g_{i}(x^{∗})) = 0  • complémentarité

[figure: ∇f(x^{∗}); ∇g_{1}(x^{∗}); ∇g_{2}(x^{∗}); x^{∗}; d; g_{1}(x) = b_{1}; g_{2}(x) = b_{2}; g_{4}(x) = b_{4}; g_{3}(x) = b_{3}]

Fig. 3 – Il n’existe pas de direction de montée admissible

En posant x = x^{∗} dans le système ci-dessus on obtient les conditions nécessaires d’optimalité de Kuhn et Tucker, qui constituent le résultat fondamental de la programmation non linéaire :

\begin{aligned}g_{i}(x^{∗}) & ≤ & b_{i} & i = 1, …, m \\ ∑_{i=1}^{m} y_{i}∇g_{i}(x^{∗}) & = & ∇f(x^{∗}) \\ y_{i} & ≥ & 0 & i = 1, …, m \\ y_{i}(b_{i} − g_{i}(x^{∗})) & = & 0 & i = 1, …, m \end{aligned}

Géométriquement, le gradient ∇f(x^{∗}) de l’objectif s’exprime comme une combinaison linéaire non négative des gradients ∇g_{i}(x^{∗}) des contraintes actives (saturées), comme en programmation linéaire.

Dans le cas où l’objectif est concave et les contraintes sont convexes^{1} , un point satisfaisant aux conditions de Kuhn et Tucker constitue un maximum global du programme non linéaire.

Une fonction d’une seule variable est convexe (voir la figure 4) si sa dérivée est croissante ou si sa dérivée seconde est toujours positive ou nulle. Une fonction de plusieurs variables est convexe si elle est convexe dans toutes les directions. Un critère utile de convexité d’une fonction f est que la matrice hessienne H(x)(matrice des dérivées secondes) soit semidéfinie positive pour toute valeur de x, c’est-à-dire que toutes les valeurs propres de H(x) soient non négatives. Un cas particulier important est celui des formes quadratiques x^{t}Ax où la matrice A est semidéfinie positive, c’est-à-dire que toutes les valeurs propres de A + A^{t} sont non négatives.

Une fonction f est concave si et seulement si son opposée −f est convexe. Noter qu’une fonction linéaire est à la fois convexe et concave.

1On parle alors d’un problème convexe.

[figure: (x^{1}, f(x^{1})); (x^{2}, f(x^{2})); concave; convexe; (x^{1}, f(x^{1})); (x^{2}, f(x^{2})); f (\frac{x^{1}+x^{2}}{2} ) ≥ \frac{1}{2} (f(x^{1}) + f(x^{2})); f (\frac{x^{1}+x^{2}}{2} ) ≤ \frac{1}{2} (f(x^{1}) + f(x^{2}))]

Fig. 4 – Fonction concave, fonction convexe

Appliquons les conditions de Kuhn et Tucker à l’exemple suivant :

max x_{1}^{2} + x_{2}^{2}
x_{1} + 2x_{2} ≤ 5
−x_{1} ≤ 0
−x_{2} ≤ −1

On obtient

x_{1} + 2x_{2} ≤ 5  y_{1} − y_{2} = 2x_{1}  y_{1}(5 − x_{1} − 2x_{2}) = 0
−x_{1} ≤ 0  2y_{1} − y_{3} = 2x_{2}  y_{2}x_{1} = 0
−x_{2} ≤ −1  y_{3}(−1 + x_{2}) = 0
y_{1}, y_{2}, y_{3} ≥ 0.

On trouve facilement la solution de ce système : x^{∗} = (3, 1) et y = (6, 0, 10). Sur la figure 5 on peut vérifier géométriquement l’égalité

∑_{i=1}^{3} y_{i}∇g_{i}(x^{∗}) = ∇f(x^{∗}).

Remarquer que, dans l’exemple précédent, l’optimum global est atteint même si le problème n’est pas convexe. Les choses auraient été différentes si le le domaine avait contenu l’origine ; l’origine, qui est un minimum global de l’objectif, aurait aussi satisfait aux conditions de KKT.

[figure: ∇f(x^{∗}) = (6, 2); ∇g_{1}(x^{∗}) = (1, 2); ∇g_{2}(x^{∗}) = (−1, 0); ∇g_{3}(x^{∗}) = (0, −1); x_{1}; x_{2}; 1; 2; 3; 1; 2; 3]

Fig. 5 – Les conditions de Kuhn et Tucker

ALGORITHMES

Ce n’est certainement pas le but de ce cours de traiter en profondeur des algorithmes de résolution pour les programmes non linéaires, ce qui occuperait le cours dans son entier. Cependant, je vais décrire très brièvement trois stratégies de résolution utilisées en pratique.

1. UNE STRATÉGIE DE PÉNALISATION DES CONTRAINTES

Soit le programme mathématique

max  f(x)
g_{i}(x) ≤ b_{i} i = 1, … m.

Une façon simple de se débarasser des contraintes est d’intégrer dans l’objectif une fonction qui pénalise les vecteurs qui violent les contraintes :

max f(x) − ∑_{i=1}^{m} P_{i}(x),

où P_{i}(x) prend une valeur élevée lorsque g_{i}(x) est strictement positif, ce qui décourage l’adoption de ces vecteurs x. On peut par exemple poser

P_{i}(x) = μ max{0, g_{i}(x) − b_{i}}

où μ est un coefficient positif et � grand �. Si l’on souhaite que l’objectif soit différentiable, on peut élever l’expression précédente au carré. Pour résoudre le problème pénalisé, il suffit alors de trouver les points où le gradient de son objectif s’annule. Mais ceci est une longue histoire qui est du ressort du cours IFT 3512 et non celui d’un premier cours de recherche opérationnelle.

2. UNE PREMIÈRE STRATÉGIE DE LINÉARISATION (FRANK-WOLFE)

Soit le problème

max_{x∈X} f(x).

Si l’ensemble X est un polyèdre, on peut se ramener à un problème de programmation linéaire en faisant une approximation linéaire de l’objectif. Ceci correspond à un développement de Taylor du premier ordre et conduit au sous-problème

max_{y∈X} f(x) + ∇f(x)(y − x).

Soit y(x) une solution de ce programme linéaire. On démontre aisément que la direction d = y(x) − x est une direction d’amélioration. On maximise alors la fonction dans la direction d :

max_{α∈[0,1]} f(x + αd),

ce qui n’est pas difficile car la fonction f(x + αd) ne comporte que la variable α. La solution de cette recherche linéaire conduit au nouveau point x^{+} d’où l’on redémarre le processus (voir figure 6). Noter que si la valeur optimale de α est inférieure à 1, alors le gradient ∇f(x^{+}) doit forcément être orthogonal à la direction d (pourquoi ?).

Cette méthode est très appréciée dans le domaine du transport où le sous-problème linéaire est fréquemment un problème à la structure très simple qui peut être résolu par des algorithmes spécialisés.

[figure: X; x; y(x); ∇f(x); ∇f(x^{+}); x^{+}; d]

Fig. 6 – Un premier algorithme de linéarisation

3. UNE SECONDE STRATÉGIE DE LINÉARISATION

Une fonction f est additive si elle peut se mettre sous la forme

f(x_{1}, …, x_{n}) = f_{1}(x_{1}) + … + f(x_{n}).

Considérons le programme mathématique

max  f(x)
g_{i}(x) ≤ b_{i} i = 1, …, m

où toutes les fonctions impliquées sont additives. Il est alors possible d’approximer chaque fonction, sur un intervalle [a, b], par une fonction linéaire par morceaux et de substituer au problème de départ le programme linéaire (voir figure 7) :

max  ∑_{k=1}^{K} λ_{k}f(x^{k})

∑_{k=1}^{K} λ_{k}g_{i}(x^{k}) ≤ b_{i} i = 1, …, m

∑_{k=1}^{K} λ_{k} = 1

λ_{k} ≥ 0 k = 1, … K.

[figure: a = x^{1}; b = x^{5}; x; f(x); x^{2}; x^{3}; x^{4}]

Fig. 7 – Approximation linéaire par morceaux d’une fonction concave

Les nouvelles variables sont maintenant les λ_{k}. On démontre, si le problème est convexe, qu’au plus deux variables λ_{k} peuvent prendre des valeurs non nulles, et que ces variables doivent avoir des indices consécutifs. Si ce n’était pas le cas, l’approximation linéaire pourrait être invalide. Noter enfin que, plus le nombre de points d’évaluation x^{k} est élevé, plus l’approximation sera fine et plus le programme linéaire à résoudre sera grand.

4. AUTRES CONSIDÉRATIONS Il n’y a pas de stratégie universelle qui permette de résoudre efficacement tous les problèmes non linéaires. Il est important de tenir compte de la structure particulière que l’on traite. Par exemple, il existe des algorithmes spécialisés pour résoudre des programmes quadratiques, c’est-à-dire des programmes où l’objectif est quadratique et où les contraintes sont linéaires. Les conditions de Kuhn et Tucker de ces programmes sont très semblables à celles de programmes linéaires et il est possible de les résoudre par des variantes de la méthode du simplexe. Pour en savoir plus, vous pouvez consulter l’excellent livre “Nonlinear Programming” de Bazaraa, Sherali et Shetti.

CHAPITRE 4 : THÉORIE DES JEUX 7 octobre 2003

La théorie des jeux constitue un domaine fascinant où l’on cherche à décrire par des équations le comportement stratégique d’êtres humains. On distingue les jeux non coopératifs où les coalitions ne sont pas admises et les jeux coopératifs où les joueurs peuvent former des coalitions afin de maximiser un revenu commun à la coalition et qu’ils se partagent par la suite, à l’intérieur de chaque coalition. Dans cette section, nous ferons l’hypothèse que les joueurs annoncent simultanément leur stratégie.

Jeux non coopératifs

Dans cette section, on ne donne pas de description détaillée du déroulement du jeu, ce qui rendrait impossible son analyse mathématique. On se contente de donner le résultat du jeu (gains) correspondant aux choix stratégiques des joueurs. Ainsi, dans le jeu d’échecs, une stratégie correspond à une règle qui détermine, pour chaque position, le mouvement à effectuer. D’un point de vue mathématique, ce jeu n’a que peu d’intérêt et n’a que 3 issues possibles : soit il existe une stratégie gagnante pour les blancs, soit il existe une stratégie gagnante pour les noirs, soit toutes les parties entre deux joueurs �parfaits� se terminent par une nulle.

Considérons un jeu impliquant n joueurs cherchant chacun à maximiser un objectif qui lui est propre. L’objectif f_{i}(x_{1}, …, x_{n}) du joueur i dépend à la fois de sa stratégie x_{i} ainsi que des stratégies de ses adversaires. Si l’on suppose que le domaine X_{i} des stratégies admissibles du joueur i ne dépend pas des stratégies de ses adversaires, celui-ci est confronté au programme mathématique

max_{x_{i}∈X_{i}} f_{i}(x_{1}, …, x_{i}, …, x_{n}).

On fera l’hypothèse que ce programme admet une solution unique. Un ensemble de stratégies (x_{1}^{∗}, …, x_{n}^{∗}) est un équilibre de Nash si aucun joueur ne peut améliorer son sort en modifiant unilatéralement sa stratégie, c’est-à-dire, pour tout i ∈ {1, …, n}, on a x_{i}^{∗} ∈ X_{i} et

f_{i}(x_{1}^{∗}, …, x_{i−1}^{∗}, x_{i}^{∗}, x_{i+1}^{∗}, …, x_{n}^{∗}) ≥ f_{i}(x_{1}^{∗}, …, x_{i−1}^{∗}, x_{i}, x_{i+1}^{∗}, …, x_{n}^{∗}) ∀x_{i} ∈ X_{i}.

Considérons par exemple un duopole où deux firmes tentent d’écouler un produit sur un marché. Le prix de vente du produit sur le marché est spécifié par une fonction décroissante de la demande totale. Cette fonction est appelée loi de demande inverse. Si c_{1}(x_{1}) et c_{2}(x_{2}) représentent, respectivement, les coûts de productions des deux firmes, les objectifs des firmes s’expriment comme

f_{i}(x_{1}, x_{2}) = x_{i}p(x_{1} + x_{2}) − c_{i}(x_{i}).

Ainsi, si p(x) = 200−x, c_{1}(x_{1}) = x_{1}^{2}+20 et c_{2}(x_{2}) = 2x_{2}^{2}+10, on obtient les deux problèmes d’optimisation :

max_{x_{1}≥0}  x_{1}(200 − (x_{1} + x_{2})) − (x_{1}^{2} + 20)
max_{x_{2}≥0}  x_{2}(200 − (x_{1} + x_{2})) − (2x_{2}^{2} + 10).

En cherchant à annuler les dérivées de chacun des objectifs, on obtient le système linéaire

\begin{aligned}200 & − & 2x_{1} & − & x_{2} & − & 2x_{1} & = & 0 \\ 200 & − & x_{1} & − & 2x_{2} & − & 4x_{2} & = & 0 \end{aligned}

dont la solution est

(x_{1}, x_{2}) = (1000/23, 600/23).

Les profits respectifs des deux firmes sont approximativement 3760 et 2031. On vérifie que les profits des deux joueurs sont non négatifs. Si ce n’était pas le cas, il faudrait reprendre l’analyse en éliminant l’une des deux firmes car, si un profit est négatif, il est préférable de ne rien produire^{1}. Il se peut même qu’à l’équilibre, la production des deux firmes soit nulle.

1Paradoxe : ce n’est pas nécessairement la firme de profit négatif qui sera éventuellement sortie du marché !

Il est à noter qu’en fusionnant, les deux firmes augmenteraient leur profit global en résolvant le programme

max_{x_{1}≥0,x_{2}≥0} x_{1}(200 − (x_{1} + x_{2})) − (x_{1}^{2} + 20) + x_{2}(200 − (x_{1} + x_{2})) − (2x_{2}^{2} + 10)

dont la solution (x_{1}, x_{2}) = (40, 20) fournit un profit total de 5970 (3980 pour la première firme et de 1990 pour la seconde). Il existe des situations où le profit de chacune des firmes augmente après fusion. Ceci nous rappelle le paradoxe de Braess en équilibre du trafic.

Il n’est pas dans notre propos de traiter des algorithmes permettant de trouver un équilibre de Nash, sinon pour indiquer qu’il est possible de procéder par tâtonnement, en résolvant successivement les problèmes des deux joueurs, en espérant que le processus converge !

Jeux matriciels

[Dans cette section on adopte la convention que tous les vecteurs sont des vecteurs colonne et l’on utilisera l’opérateur de transposition]

Peut-être la forme la plus simple d’un jeu est celle où chaque joueur possède un nombre fini de stratégies et où les gains associés à chaque combinaison de stratégies est connue des deux joueurs. Si le gain du premier joueur correspond au déboursé du second, on dit que le jeu est à somme nulle. (Noter que le duopole de la section précédente n’est pas à somme nulle.) Un tel jeu peu se représenter sous forme matricielle, le premier joueur étant associé aux lignes de la matrice et le second à ses colonnes. Les coefficients de la matrice représentent les gains du premier joueur, qui cherche à maximiser ses gains, alors que le second joueur tente de minimiser ses pertes, c’est-à-dire les gains du joueur ligne.

Considérons le jeu matriciel associé à la matrice

[ \begin{matrix}4 & 3 & 10 \\ 6 & 5 & 7 \end{matrix} ] .

Il est facile de vérifier que le couple stratégique (2, 2) constitue un équilibre de Nash, et que le profit du joueur ligne, qu’on appelle la valeur du jeu, est égal à 5. Le coefficient a_{22} = 5 de la matrice est un point de selle, c’est-à-dire qu’il constitue le minimum de sa ligne et le maximum de sa colonne.

Considérons maintenant la matrice

[ \begin{matrix}−1 & 1 \\ 3 & −2 \end{matrix} ] ,

qui ne possède pas de point de selle. Le jeu correspondant ne possède donc pas d’équilibre de Nash. Une façon élégante de corriger la situation est d’élargir l’ensemble des stratégies en introduisant un élément aléatoire dans le jeu. Une stratégie mixte est un vecteur de probabilités dont la dimension est égale au nombre de stratégies. L’ensemble des stratégies mixtes du premier joueur est

X = {x = (x_{1}, …, x_{m}) : ∑_{i=1}^{m} x_{i} = 1, x ≥ 0}.

Ainsi, la probabilité que le joueur ligne adopte la stratégie i est égale à x_{i}. Symétriquement, l’ensemble des stratégies mixtes du joueur colonne est

Y = {y = (y_{1}, …, y_{m}) : ∑_{j=1}^{n} y_{j} = 1, y ≥ 0}.

La probabilité que le couple stratégique (i, j) soit observé est égale à x_{i}y_{j}, et le gain moyen (espérance de gain) du joueur ligne est donnée par l’expression

∑_{i=1}^{m} ∑_{j=1}^{n} a_{ij}x_{i}y_{j} = x^{t}Ay.

Un équilibre de Nash (x, y) devra satisfaire aux conditions d’optimalité des programmes

joueur 1  joueur 2
max_{x} x^{t}Ay  min_{y} x^{t}Ay
∑_{i=1}^{m} x_{i} = 1  ∑_{j=1}^{n} y_{j} = 1.
x ≥ 0  y ≥ 0

Les conditions d’optimalité de ces programmes linéaires sont

joueur 1  joueur 2
∑_{i=1}^{m} x_{i} = 1  ∑_{j=1}^{n} y_{j} = 1
x ≥ 0  y ≥ 0

w ( 1 ⋮ 1 ) ≥ Ay  v ( 1 ⋮ 1 ) ≤ A^{t}x

x^{t}Ay = w  x^{t}Ay = v.

On en déduit que la valeur du jeu est égale à v = w = x^{t}Ay. Par une coïncidence extraordinaire, en interchangeant le deuxième groupe de contraintes de chacun de ces systèmes et en oubliant la dernière ligne, on obtient les programmes linéaires primal et dual suivants :

max v  min w
∑_{i=1}^{m} x_{i} = 1  ∑_{j=1}^{n} y_{j} = 1
x ≥ 0  y ≥ 0

v ( 1 ⋮ 1 ) ≤ A^{t}x  w ( 1 ⋮ 1 ) ≥ Ay.

Ces deux programmes ont une interprétation intuitive. En effet, chaque joueur y détermine une stratégie visant à maximiser son gain minimal ou à minimiser sa perte maximale. De telles stratégies sont appelées, respectivement, maximin et minimax. Chaque joueur se prémunit contre la pire situation possible. Ces programmes linéaires sont équivalents aux programmes

max_{x∈X} min {A^{t}x}_{i}  min_{y∈Y} max {Ay}_{i}.

Dans le petit exemple précédent, on obtient

max_{x∈X} min {−x_{1} + 3x_{2}, x_{1} − 2x_{2}}  min_{y∈Y} max {−y_{1} + y_{2}, 3y_{1} − 2y_{2}}.

En posant x_{2} = 1 − x_{1} on obtient le problème unidimensionnel

max_{x∈[0,1]} min {3 − 4x_{1}, 3x_{1} − 2},

dont la solution est x_{1} = 5/7. Le vecteur stratégique optimal est donc x = (5/7, 2/7) et la valeur du jeu est égale à 1/7 (voir figure 1). De la même façon, on obtient le vecteur stratégique y = (3/7, 4/7) pour le second joueur.

[figure: x_{1}; 3 − 4x_{1}; 3x_{1} − 2; 0; 1; 2; 3; −1; −2; 1; \frac{5}{7}; \frac{2}{7}; (\frac{5}{7}, \frac{1}{7})]

Fig. 1 – Calcul graphique de l’équilibre

Jeux coopératifs

Le concept d’équilibre de Nash ne permet pas de modéliser les situations où des coalitions de joueur peuvent se former. La forme caractéristique d’un jeu a été introduite pour traiter correctement de telles situations. Un tel jeu est spécifié entièrement par sa fonction caractéristique v qui associe à chaque coalition de joueurs S ⊆ N = {1, …, n} la valeur minimale qu’elle peut obtenir indépendamment du comportement des joueurs qui ne font pas partie de la coalition.

On fait l’hypothèse raisonnable que la fonction caractéristique est superadditive, c’est-à-dire que, pour tout couple de coalitions A et B ne partageant aucun membre commun, on a

v(A ∪ B) ≥ v(A) + v(B).

En d’autres mots, l’union fait la force ! Une imputation x = (x_{1}, …, x_{n}) est une répartition du gain maximal v(N) parmi les joueurs, qui satisfait la condition supplémentaire

x_{i} ≥ v({i}).

Cette dernière condition est tout à fait raisonnable. En effet, aucun joueur rationnel n’accepterait de recevoir moins que ce qu’il peut obtenir en ne coopérant pas. On peut généraliser cette condition à toutes les coalitions, c’est-à-dire

∑_{i∈S} x_{i} ≥ v(S)

pour toutes les coalitions S (on a l’égalité si S = N). L’ensemble des imputations satisfaisant cette dernière condition est appelé le cœur du jeu.

Exemple

Un inventeur (joueur 1) cherche une compagnie (joueur 2 ou joueur 3) pour commercialiser son produit. Si le produit est commercialisé, il rapportera 10^{6} euros. La fonction caractéristique du jeu est, naturellement :

v(∅) = v({1}) = v({2}) = v({3}) = v({2, 3}) = 0
v({1, 2}) = v({1, 3}) = v({1, 2, 3}) = 10^{6}.

Les équations définissant le cœur du jeu sont

x_{1} + x_{2} ≥ 10^{6}
x_{1} + x_{3} ≥ 10^{6}
x_{2} + x_{3} ≥ 0
x_{1} + x_{2} + x_{3} = 10^{6}
x_{1}, x_{2}, x_{3} ≥ 0

dont l’unique solution est x = (10^{6}, 0, 0).

Dans ce cas, le cœur contient un seul élément, ce qui est rassurant. Malheureusement, il est facile de construire des exemples où le cœur du jeu (un polyèdre convexe) contient une nombre infini d’éléments ou, pire, est vide. Plusieur propositions ont été faites afin de remédier à cette situation. Elles ont toutes leurs qualités et leurs défauts.

Extensions

Ce qui précède n’est évidemment qu’un bref survol d’un domaine très riche. En vrac, voici quelques notions qui ont fait l’objet d’études sérieuses :

jeux répétés Jeux où les joueurs peuvent adapter leurs stratégies en observant le comportement des

adversaires. jeux dynamiques Jeux qui se déroulent dans le temps. jeux différentiels Jeux qui se déroulent en temps réels et sont décrits par des systèmes d’équations

différentiels. menaces Modélisation d’un élément important ! négociation Important dans les jeux coopératifs. indices de pouvoir Concepts alternatifs à celui de cœur. propension au risque Certains joueurs ont une aversion au risque et préfère �jouer sûr�, c’est-à-dire

accepter un gain moins élevé associé à un risque (variance du gain) plus faible.

CHAPITRE 5 : PROGRAMMATION LINÉAIRE AVEC VARIABLES ENTIÈRES dernières modifications : 28 février 2001

Lorsque les variables d’un programme mathématique représentent des décisions stratégiques, ou simplement lorsque certaines variables ne peuvent prendre des valeurs fractionnaires, il faut intégrer dans le modèle des contraintes d’intégralité. Ceci modifie profondément la nature des programmes linéaires.

Prenons l’exemple, peu édifiant je l’avoue, du cambrioleur muni d’un sac qui lui sert à transporter le butin tiré d’un coffre-fort^{1} : bijoux, argent, caviar, portos anciens, objets divers. 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. Si l’on associe à l’objet j une valeur c_{j} et un poids w_{j}, la combinaison optimale d’objets à emporter sera obtenue en résolvant le programme mathématique

max  ∑_{j=1}^{n} c_{j}x_{j}

∑_{j=1}^{n} w_{j}x_{j} ≤ b

x_{j} entier non négatif.

La solution de ce problème semble triviale : pourquoi ne pas choisir de façon prioritaire les objets dont le rapport qualité-poids est le plus avantageux, quitte à tronquer pour obtenir une solution entière (il n’y a pas de demi-bouteilles de porto dans le coffre-fort) ? Et bien, tout simplement parce que cette solution peut se révéler sous-optimale, voire mauvaise, comme on le constate sur l’exemple suivant.

max 2x_{1} + 3x_{2} + 7x_{3}
3x_{1} + 4x_{2} + 8x_{3} ≤ 14
x_{1},  x_{2},  x_{3} entiers non négatifs.

[figure: x_{1}; x_{2}; c; tronqué; optimal; 1; 2; 1; 2]

Fig. 1 – Un programme en variables entières Si l’on oublie la contrainte d’intégralité, la solution, de valeur 12\frac{1}{4}, 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. Une représentation géométrique d’un autre exemple, impliquant 2 variables, est donné à la figure 1.

1Le coffre-fort, immense, contient un nombre incalculable d’objets de chaque type.

La solution optimale du problème du cambrioleur peut s’obtenir en énumérant toutes les solutions admissibles et en conservant la meilleure (voir tableau ci-dessous, où les solutions inefficaces laissant la possibilité d’ajouter un objet n’ont pas été considérées).

x_{1} x_{2} x_{3} 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

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. J’y reviendrai après avoir formulé d’autres problèmes où certaines variables ne peuvent prendre des valeurs fractionnaires.

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. Chaque service de vol est constitué d’une séquence de segments de vol (voir figure 2).

[figure: domicile; s_{1}; s_{3}; s_{2}; s_{5}; s_{9}; s_{2}; s_{6}; segment de vol; quart de travail]

Fig. 2 – Un exemple de rotation

Soit T l’ensemble des rotations, T (i) l’ensemble des rotations contenant le segment de vol i et c_{j} le coût associé à la rotation j. Soit x_{j} une variable binaire prenant la valeur 1 si la rotation j est utilisée. Le problème d’affectation des équipages aux rotations peut s’exprimer comme le problème de recouvrement

min  ∑_{j∈T} c_{j}x_{j}

∑_{j∈T (i)} x_{j} ≥ 1

x_{j} ∈ {0, 1}.

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

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 a_{i}, une capacité de production u_{i}, un coût de production unitaire b_{i} et des coûts de transport c_{ij} des usines vers les clients.

Soit y_{i} une variable binaire prenant la valeur 1 si un entrepôt est construit sur le site i, d_{j} la demande de l’usine j et x_{ij} 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

min  ∑_{i} a_{i}y_{i} + ∑_{i} b_{i} ∑_{j} x_{ij} + ∑_{i} ∑_{j} c_{ij}x_{ij}

∑_{i} x_{ij} = d_{j}

∑_{j} x_{ij} ≤ u_{i}y_{i}

x_{ij} ≥ 0 y_{i} ∈ {0, 1}.

Cette formulation contient deux éléments intéressants : un coût fixe (construction) modélisé par une variable binaire y_{i} 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. Un tel programme linéaire est appelé mixte.

Contraintes logiques

Des variables binaires peuvent servir à représenter des contraintes logiques. En voici quelques exemples, où p_{i} représente une proposition logique et x_{i} la variable logique (binaire) correspondante. Vous pouvez agrandir cette liste par vous-même.

contrainte logique  forme algébrique
p_{1} ⊕ p_{2} = vrai  x_{1} + x_{2} = 1
p_{1} ∨ p_{2} ∨ … ∨ p_{n} = vrai  x_{1} + x_{2} + … + x_{n} ≥ 1
p_{1} ∧ p_{2} ∧ … ∧ p_{n} = vrai  x_{1} + x_{2} + … + x_{n} ≥ n (ou = n)
p_{1} ⇒ p_{2}  x_{2} ≥ x_{1}
p_{1} ⇔ p_{2}  x_{2} = x_{1}

Fonctions linéaires par morceaux

Considérons une fonction objectif à maximiser de la forme illustrée à la figure ci-contre. Dans chaque intervalle [a_{i−1}, a_{i}] la fonction est linéaire, ce qu’on peut exprimer par :

x = λ_{i−1}a_{i−1} + λ_{i}a_{i}
λ_{i−1} + λ_{i} = 1
λ_{i−1}, λ_{i} ≥ 0
f(x) = λ_{i−1}f(a_{i−1}) + λ_{i}f(a_{i}).

On peut 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 y_{i} associées aux intervalles de linéarité [a_{i−1}, a_{i}] (voir figure 3) :

[figure: a_{0}; a_{1}; a_{2}; a_{3}; a_{4}; f(a_{0}); f(a_{1}); f(a_{2}); f(a_{3}); f(a_{4}); x; f(x); y_{1}; y_{2}; y_{3}; y_{4}]

Fig. 3 – Fonction linéaire par morceaux

x = ∑_{i=0}^{n} λ_{i}a_{i}

f(x) = ∑_{i=0}^{n} λ_{i}f(a_{i})

∑_{i=0}^{n} λ_{i} = 1
λ_{i} ≥ 0  i = 0, …, n
λ_{0} ≤ y_{1}
λ_{1} ≤ y_{1} + y_{2}
⋮ ⋮ ⋮
λ_{n−1} ≤ y_{n−1} + y_{n}
λ_{n} ≤ y_{n}
∑_{i=1}^{n} y_{i} = 1  (un seul intervalle � actif �)
y_{i} ∈ {0, 1}  i = 1, …, n.

Si la fonction f est concave, les variables binaires y_{i} et les contraintes associées peuvent être éliminées de la formulation. Cette approche est particulièrement intéressante lorsqu’on optimise des fonctions de plusieurs variables de la forme ∑_{i} f_{i}(x_{i}).

Une première stratégie de résolution : l’énumération implicite

Considérons un programme mathématique mixte. Par paresse, il est naturel de laisser de côté, dans un premier temps, les contraintes d’intégralité, et de se contenter de résoudre un programme linéaire (relaxation linéaire notée RL) du programme en variables entières). 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 x_{j} dont la valeur α est fractionnaire. On sépare alors le problème en deux : un sous-problème contiendra la contrainte x_{j} ≤ ⌊α⌋ et le second la contrainte x_{j} ≥ ⌈α⌉ = ⌊α⌋ + 1. Il est clair que ceci crée une partition du problème relaxé. ^{2} On répète 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.

2Bon, je sais que cette description n’est pas très rigoureuse, mais je pense qu’elle est compréhensible.

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.

Cette technique d’énumération partielle est connue sous le vocable de Branch-and-Bound. Appliquons-la au problème du cambrioleur (voir figure 4). Dans la figure, la solution optimale de chaque relaxation linéaire est indiquée sous le nœud correspondant. Cette solution est obtenue en favorisant les variables de rapport qualité/poids élevé tout en respectant les contraintes rajoutées sur les branches de l’arbre d’énumération.

[figure: RL_{0}; RL_{1}; RL_{6}; (0,0,14/8); obj=12,25; (0,6/4,1); obj=11,5; sondé; (non réalisable); x_{3} ≤ 1; x_{3} ≥ 2; RL_{2}; x_{2} ≤ 1; RL_{3}; RL_{4}; (2/3,1,1); obj=11\frac{1}{3}; (0,1,1); obj=10; sondé; (réalisable); (1,3/4,1); obj=11\frac{1}{4}; RL_{5}; (2,0,1); obj=11; sondé; (réalisable); sondé (pas meilleur que RL_{5}); sondé (pas meilleur que RL_{5}); x_{1} = 0; x_{1} ≥ 1; x_{2} = 0]

Fig. 4 – La procédure de Branch-and-Bound

Plusieurs remarques sont de rigueur :

1. Dans l’exemple, l’algorithme ne semble pas plus efficace que l’énumération complète. En effet il y a 17 solutions admissibles (vérifiez-le !), dont seulement 6 sont efficaces, alors que l’arbre de Branchand-Bound possède 7 sommets et exige une certaine quantité de travail pour traiter chaque sommet. Sur des problèmes de grande taille, par contre, le nombre de sommets visités est en général beaucoup plus faible que le nombre de solutions réalisables.

2. On peut tirer profit de la structure du problème pour obtenir des informations additionnelles sur la valeur des variables. Ainsi, il est clair que x_{3} ne peut excéder 2. On peut alors remplacer une contrainte de la forme x_{3} ≥ 1 par x_{3} = 1, ce qui fixe la valeur de x_{3} et l’élimine du sous-problème.

3. Le programme linéaire associé à un sommet ne diffère du problème père que par une contrainte. Il est alors judicieux de résoudre ce problème par l’algorithme dual du simplexe, c’est-à-dire l’algorithme du simplexe appliqué au problème dual.

4. L’ordre dans lequel l’arbre est exploré est important. Si RL5 n’avait pas été visité avant RL_{6}, celui-ci n’aurait pu être sondé aussi tôt.

5. La solution optimale est obtenue au nœud RL_{5}. Il faut cependant pousser l’exploration un plus loin pour démontrer que cette solution est bien optimale.^{3} Dans plusieurs applications, on perd beaucoup plus de temps à démontrer qu’une solution est optimale qu’à la trouver !

6. Si l’objectif du problème linéaire relaxé vaut z, la solution entière ne peut valoir plus que ⌊z⌋ puisque les coefficients de l’objectif sont entiers. Ceci ne serait pas le cas si les coefficients étaient fractionnaires.

7. Deux stratégies d’exploration sont utilisées en pratique : fouille en profondeur (c’est celle utilisée dans l’exemple) ou exploration prioritaire des sommets correspondant aux meilleures valeurs de l’objectif. Cette dernière stratégie a l’avantage de limiter le nombre de sommets visités. Par contre, elle ne permet pas de tirer profit de la structure du problème père, à moins de garder en mémoire les dictionnaires optimaux de tous les sommets, ce qui est coûteux.

8. Au sommet initial, aucune solution admissible n’est encore disponible. On pourrait en obtenir à l’aide d’algorithmes heuristiques. Dans le problème du cambrioleur, on pourrait par exemple choisir les variables dans l’ordre décroissant de leur rapport qualité-prix, en exigeant que les variables soient entières. Ceci donnerait la solution x = (0, 1, 1) de valeur 10. Cette valeur pourrait servir à sonder de futurs sommets. Puisque les sous-problèmes de chaque sommet possèdent la structure du problème de départ, cette technique pourrait être appliquée à chaque sommet de l’arbre.

9. De façon surprenante, la méthode de Branch-and-Bound fonctionne souvent bien. Malgré la taille théoriquement énorme de l’arbre d’énumération, on parvient fréquemment à une solution optimale ou quasi-optimale après avoir exploré quelques dizaines de sommets.

10. La méthode peut facilement être adaptée pour trouver toutes les solutions optimales du problème.

Une seconde stratégie : les hyperplans coupants

Considérons le programme mathématique

max  4x_{1} + \frac{5}{2}x_{2}
x_{1} + x_{2} ≤ 6
9x_{1} + 5x_{2} ≤ 45
x_{1}, x_{2} entiers non négatifs.

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

x_{1} = \frac{15}{4} + \frac{5}{4}x_{3} − \frac{1}{4}x_{4}
x_{2} = \frac{9}{4} − \frac{9}{4}x_{3} + \frac{1}{4}x_{4},

où x_{3} et x_{4} sont des variables d’écart. Puisque la variable de base x_{1} n’est pas entière, cette solution de base n’est pas admissible. On réécrit la première contrainte sous la forme

x_{1} − \frac{5}{4}x_{3} + \frac{1}{4}x_{4} = \frac{15}{4}.

3Encore que, dans le problème du cambrioleur, il est clair que la partie de l’arbre située à droite du sommet initial RL_{0} ne sera jamais admissible.

En utilisant l’identité

a = ⌊a⌋ + (a − ⌊a⌋),

où a − ⌊a⌋ représente la partie fractionnaire de a (0 ≤ a − ⌊a⌋ < 1), on obtient

x_{1} + (⌊−\frac{5}{4}⌋ + \frac{3}{4}) x_{3} + (⌊\frac{1}{4}⌋ + \frac{1}{4}) x_{4} = (⌊\frac{15}{4}⌋ + \frac{3}{4}) ,

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

x_{1} − 2x_{3} − 3 = \frac{3}{4} − \frac{3}{4}x_{3} − \frac{1}{4}x_{4}.

Puisque les variables x_{3} et x_{4} 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. On en déduit 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 :

\frac{3}{4} − \frac{3}{4}x_{3} − \frac{1}{4}x_{4} ≤ 0.

En utilisant les identités x_{3} = 6 − x_{1} − x_{2} et x_{4} = 45 − 9x_{1} − 5x_{2} on obtient la coupe sous sa forme géométrique :

3x_{1} + 2x_{2} ≤ 15.

Cette contrainte linéaire rend inadmissible la solution courante inadmissible ^{4} , sans éliminer aucune autre solution entière (voir figure 5). Si la solution du nouveau problème est entière, c’est la solution optimale de notre problème. Sinon, on construit une nouvelle coupe et on recommence.

La résolution du nouveau programme linéaire s’effectue, comme dans l’algorithme de Branch-and-Bound, à l’aide de l’algorithme du simplexe appliqué au problème dual, afin de tirer profit de la solution duale optimale obtenue précédemment.

Cette technique seule ne permet pas de résoudre des problèmes de grande taille. Cependant elle peut être combinée de façon efficace avec la technique de Branch-and-Bound.

4Dans la solution courante, les variables hors-base x_{3} et x_{4} sont nulles, ce qui donne 3/4 ≤ 0 !

[figure: x_{1}; x_{2}; x_{1} + x_{2} = 6; 3x_{1} + 2x_{2} = 15 (coupe); 9x_{1} + 5x_{2} = 45; (3,3); (\frac{15}{4}, \frac{9}{4}); (4,5/2); (4,5/2)]

Fig. 5 – Une droite coupante CHAPITRE 6 : RÉSEAUX 19 novembre 2003

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 arcs possèdent des attributs : coûts, capacités, longueurs, etc. Ce qui est amusant avec les réseaux, c’est qu’ils permettent de visualiser des situations abstraites.^{1}

Dans ce chapitre nous allons étudier quatre problèmes classiques de la théorie des réseaux. Dans le texte, les termes �distance�, �coût� et �longueur� sont synonymes.

1. Le problème de flot à coût minimum

Un très grand nombre de programmes linéaires possèdent une composante �réseau� très importante. Nous allons étudier dans cette section le problème consistant à approvisionner, au moindre coût, des entrepôts à partir des unités de production (usines) d’une compagnie. Les données du problème sont :

– le réseau de transport ;

– la localisation et la capacité de production des usines (sources) ;

– la localisation et la demande des entrepôts (destinations) ;

– les coûts de transport c_{ij} sur les arcs du réseau. Il reste à déterminer les productions des usines ainsi que les chemins utilisés pour approvisionner les entrepôts. Désignons par b_{i} la demande associée au sommet i. Par convention, un sommet d’offre (usine) possède une demande négative et un sommet de passage (ni usine ni entrepôt) possède une demande nulle. Pour simplifier, nous supposerons que l’offre est égale à la demande, c’est-à-dire : ∑_{i} b_{i} = 0.

Il semble naturel de formuler ce problème en terme de flots sur les chemins reliant les sources et les destinations du résau. Or, le nombre de chemins étant astronomique, il est plus simple d’utiliser les flots sur les arcs comme variables, quitte à récupérer les flots de chemins par la suite.

Soit x_{ij} le flot sur l’arc (i, j). De façon naturelle, les flots doivent être non négatifs. Ils doivent également satisfaire aux conditions de conservation qui stipulent que le flot ne peut être créé qu’aux sources et absorbé aux destinations. En tout autre sommet du réseau, la flot entrant doit correspondre au flot sortant. (voir figure 1). Le PFCM prend ainsi la forme du programme linéaire

min_{x∈X} cx

où

X = {x ≥ 0 : ∑_{j} x_{ji} − ∑_{j} x_{ij} = b_{i}, ∀i ∈ N} = {x ≥ 0 : Ax = b}.

[figure: 1; 2; 4; 5; 3; 0; 6; −b_{2} = 3; x_{12} − (x_{24} + x_{25}) = b_{2}; 1; 3; 4; 5; b_{3} = 2; (x_{13} + x_{53}) − x_{34} = b_{3}; 5; 1; 2; 3; 4; x_{15} + x_{25} = x_{53} + x_{54}; 4; 1; 3; 2; 1; 3; 4]

Fig. 1 – Conservation du flot (offre, demande, passage)

1Pour obtenir une liste des termes les plus fréquents de la théorie des graphes, vous pouvez consulter mes notes du cours IFT 1063.

[figure: 1; 2; 5; 3; 4; 5; 3; 2; 6; 10$; 1$; 3$; 4$; −3$; 9$; 2$; 7$]

Fig. 2 – Exemple de PFCM

Ci-dessous vous trouverez toutes les données de l’exemple que nous allons tenter de résoudre (voir figure 2). La matrice A représente la matrice d’incidence sommets-arcs du graphe. Chaque colonne décrit un arc alors que la ligne i contient l’information pertinente au sommet i.

A = [ \begin{matrix}−1 & −1 & −1 \\ +1 & & & −1 & −1 \\ & +1 & & & & −1 & +1 \\ & & & +1 & & +1 & & +1 \\ & & +1 & & +1 & & −1 & −1 \end{matrix} ]

b = [ \begin{matrix}−5 \\ −3 \\ +2 \\ +6 \\ 0 \end{matrix} ] c = [1, 10, 3, 7, 9, −3, 4, 2].

Ce programme étant linéaire, il est naturel de lui appliquer l’algorithme du simplexe. Celui-ci prend une forme très particulière sur les réseaux. Sans entrer dans les détails techniques, mentionnons qu’un point extrémal (solution de base) est obtenu en affectant le flot uniquement sur un arbre du réseau touchant à tous les sommets. À chaque arbre correspond une et une seule affectation du flot qui satisfait aux offres et aux demandes du problème. Par exemple, sur l’arbre de la figure 3, le flot sur les arcs est obtenu en parcourant les arcs dans l’ordre inverse de leur niveau dans l’arbre (ordre topologique inverse, c’est-à-dire en allant des feuilles vers la racine). ^{2} Le coût de cette solution initiale est de 65$.

On introduit, à chaque sommet i, un nombre y_{i} représentant la longueur de l’unique chemin allant de la racine au sommet i dans l’arbre. Si un arc est utilisé à rebours, son coût change de signe. La distance entre deux sommets i et j du réseau, en passant par l’arbre, est alors y_{j} − y_{i}. Par exemple, la distance de 5 à 3 peut être obtenue de 2 façons (voir figure 4) :

distance = y_{3} − y_{5}  calcul direct
= −c_{51} + c_{13} en passant par le chemin 5 − 1 − 3.

Les variables associées aux arcs de l’arbre sont les variables de base. Suivant la philosophie de l’algorithme du simplexe, l’on se pose la question suivante : que se passe-t-il si l’on augmente la valeur d’une variable hors-base, par exemple la variable x_{53} ? L’objectif diminuera si le coût de l’arc (5, 3) (coût direct) est inférieur au coût obtenu en passant par l’arbre pour se rendre de 5 à 3. Le coût �dans l’arbre� étant égal à y_{3} − y_{5} = 10 − 3 = 7, il est avantageux d’utiliser l’arc (5,3) dont le coût est 4 < 7. De façon générale, l’augmentation d’une variable hors-base fera diminuer l’objectif si et seulement si

c_{ij} < y_{j} − y_{i},

2La racine de l’arbre est choisie de façon arbitraire (sommet 1 dans l’exemple).

[figure: 1; 2; 5; 3; 4; 5; 3; 2; 6; 2; 3; 0; 6]

Fig. 3 – Solution de base (point extrémal)

[figure: 1; 2; 5; 3; 4; 10$; 1$; 3$; 7$; 0; 10; 1; 3; 8; 4$; y_{i}]

Fig. 4 – Distance entre deux sommets dans l’arbre

c’est-à-dire

c̄_{ij} = c_{ij} + y_{i} − y_{j} < 0.

Le nombre c̄_{ij} est le coût réduit associé à l’arc (variable) hors-base (i, j). Dans notre cas :

c̄_{53} = 4$ + 3$ − 10$ = −3$.

En augmentant le flot d’une quantité ϵ sur l’arc (5,3), il est nécessaire d’effectuer certains ajustements pour maintenir la conservation du flot. Ceci ne peut se faire que d’une seule façon : en diminuant le flot d’une quantité ϵ sur le chemin allant de 5 à 3 dans l’arbre, tout en tenant compte du sens des arcs (une diminution du flot d’un arc parcouru dans le mauvais sens correspond à une augmentation du flot).

Mais voyons les choses de façon légèrement différentes. En introduisant l’arc (5,3), on crée un cycle unique (5 − 3 − 1 − 5) dans le réseau. En modifiant le flot par une même quantité (positive pour les arcs dans le même sens que (5,3) et négative pour les autres) sur tous les arcs du cycle, on obtient un flot admissible, à la condition que tous les flots demeurent non négatifs. On augmente donc le flot x_{53} de la plus grande quantité ϵ qui satisfasse à cette dernière contrainte. Lorsque x_{53} prend la valeur ϵ, les flots du circuit prennent les valeurs respectives x_{13} = 2 − ϵ (sens contraire) et x_{15} = 0 + ϵ (bon sens). Le nombre ϵ ne peut donc dépasser la valeur 2. Pour ϵ = 2, la variable x_{53} intègre la base où elle prend la valeur 2, la variable x_{13} quitte la base, et la variable de base x_{15} voit sa valeur passer de 0 à 2 (voir figure 5). Le résultat, comme vous pouvez le constater, est une nouvelle solution de base (arbre).

La règle générale déterminant simultanément la valeur maximale de ϵ et la variable quittant la base est :

[figure: 1; 5; 3; 2 − ϵ; 0 + ϵ; 0 + ϵ; arc entrant]

Fig. 5 – Modification du flot dans le cycle

ϵ = min_{(k,l)∈cycle}{x_{kl} : (k, l) est dans le sens inverse de l’arc entrant}.

L’opération décrite ci-dessus correspond à un changement de base dans l’algorithme du simplexe. Le nouvel objectif est obtenu en additionnant à l’ancien objectif le produit du coût réduit par la modification du flot sur l’arc entrant (ϵ) :

nouvel objectif = 65$ + (−3$) × 2 = 59$.

Une fois le pivot effectué, on met à jour les nombres y_{i}^{3} et on se remet à la recherche d’une variable de coût réduit c̄_{ij} < 0. On découvre (voir figure 6) :

c̄_{34} = c_{34} + y_{3} − y_{4} = −3$ + 7$ − 8$ = −4$ < 0.

En parcourant le cycle dont le premier arc est (3,4), on réalise (voir figure 6) que la valeur maximale de ϵ est

min{6 − ϵ, 3 − ϵ} = 3.

[figure: 1; 2; 5; 3; 4; 3 − ϵ; 2 + ϵ; 6 − ϵ; 2 + ϵ; 0 + ϵ; 0; 7; 1; 8; 3]

Fig. 6 – Deuxième itération

Le nouvel objectif est égal à

59$ − c̄_{34}ϵ = 59$ − 4$ × 3 = 47$.

On met à jour l’arbre, on réévalue les nombres y_{i} et on obtient l’arbre de la figure 7. Les coûts réduits des variables hors base sont maintenant :

\begin{aligned}c̄_{12} & = & 1$ + 0$ − (−3$) & = & 4$ & ≥ & 0 \\ c̄_{13} & = & 10$ + 0$ − 7$ & = & 3$ & ≥ & 0 \\ c̄_{25} & = & 9$ + (−3$) − 3$ & = & 3$ & ≥ & 0. \end{aligned}

Puisque ceux-ci sont tous positifs ou nuls, la solution est optimale.

3Si on utilise le sommet 5 plutôt que le sommet 1 comme racine on obtient : y_{5} = 0, y_{1} = y_{5} − c_{15} (l’arc (1,5) est emprunté à contre-sens), y_{2} = y_{1} + c_{12} = −2, y_{3} = 4 et y_{4} = 5. On vérifie que les coûts réduits des variables hors-base ne sont pas modifiés car tous les y_{i} ont diminué d’une même quantité 3.

[figure: 1; 2; 5; 3; 4; 5; 3; 5; 3; 0; −3; 3; 7; 4]

Fig. 7 – Solution optimale

Remarques

1. Uniquement la différence y_{i} − y_{j} importe. On pourrait ajouter une constante à tous les nombres y_{i} sans changer quoi que ce soit à l’algorithme. C’est ce qui se passe lorsqu’on change de racine.

2. Le choix de la racine d’un arbre est arbitraire. Les seules considérations sont d’ordre informatique. En effet, il est parfois utile, pour minimiser les modifications lors du pivot, de changer de racine. Ceci implique des manipulations non triviales des structures de données associées au graphe

3. En cas d’égalité, le choix de la variable sortante est arbitraire. La solution de base subséquente est alors dégénérée.

4. Si, à une itération donnée, les variables du cycle augmentent toutes, alors le problème est non borné.

5. Pour trouver une solution de base initiale, on peut utiliser la méthode dite de la phase I, qui consiste à mettre tout le flot sur des arcs artificiels liant les sources aux destinations. Les arcs artificiels ont un coût de 1 alors que les autres arcs ont un coût nul. La technique dite du �grand M� associerait aux arcs artificiels un coût de M alors que les autres arcs conserveraient leur coût originel.

6. En pratique, les arcs possèdent fréquemment des bornes supérieures sur le flot qu’ils peuvent supporter. Dans ce cas, on utilise une variante de la méthode du simplexe, dite �simplexe avec bornes sur les variables� où une variable est hors-base lorsque sa valeur est égale soit à zéro, soit à sa borne supérieure. Le choix des variables d’entrée est légèrement modifié pour tenir compte de cette nouvelle définition.

7. Il est facile de trouver des itinéraires des flots (chemins) qui soient compatibles avec les flots sur les arcs obtenus par l’algorithme (exercice).

2. Flot maximum dans un réseau

Ce problème ultra classique consiste à maximiser le flot d’une source (origine) s à une destination t dans un réseau dont les arcs sont munis de capacités u_{ij} (voir figure 8). Il s’exprime comme un programme linéaire similaire à celui du PFCM :

max_{v,x}  v

∑_{j} x_{ji} − ∑_{j} x_{ij} = \{ \begin{matrix}−v & si i = s \\ v & si i = t \\ 0 & sinon \end{matrix}

0 ≤ x_{ij} ≤ u_{ij}.

En fait, il est possible de formuler ce problème comme un PFCM. Cependant, nous allons présenter un algorithme spécialement adapté au problème du flot maximum.

[figure: 1; 2; 3; 4; 5; 6; source; destination; 7; 3; 1; 6; 3; 8; 2; 2; 8; s = 1; t = 6; capacité]

Fig. 8 – Un problème de flot maximal dans un réseau

Pour résoudre le problème, il est naturel de trouver un chemin non saturé et d’y pousser le plus grand flot possible. Ceci s’effectue par une procédure d’étiquetage. Un sommet i est étiqueté s’il existe un chemin non saturé menant de la source s au sommet i. Soit S l’ensemble des sommets étiquetés (s ∈ S par définition) et T son complément. Si (i, j) ∈ S × T et que x_{ij} < u_{ij}, on retire le sommet j de T pour l’intégrer à S, en conservant l’information sur le prédécesseur (sommet i) afin de pouvoir retracer le chemin. La procédure d’étiquetage se termine si soit

– t est étiqueté ;

– il n’y a plus aucun sommet à étiqueter et t ∈ T. Dans le deuxième cas, l’algorithme se termine. Dans le premier cas, on a découvert un chemin

p = (i_{0} = s, i_{1}, …, i_{l−1}, i_{l} = t)

de capacité résiduelle positive. L’augmentation maximale de flot sur ce chemin est égale à

ϵ = min_{1≤k≤l}{u_{i_{k−1}i_{k}} − x_{i_{k−1}i_{k}}}.

La figure 9 illustre l’algorithme. Au départ, le flot est nul. Puis, par la procédure d’étiquetage, on trouve le chemin d’augmentation 1 − 2 − 5 − 6 dont la capacité est 4. On augmente le flot de 4 unités le long de ce chemin et on recommence. À la deuxième itération, on trouve le chemin 1−4−5−6 de capacité résiduelle 1. On augmente le flot de une unité et on s’arrête car plus aucun ne chemin ne posséde de capacité résiduelle positive.

Malheureusement, la solution obtenue par cette procédure n’est pas optimale. Il est facile de vérifier qu’en affectant un flot de 2 sur chacun des chemins 1 − 2 − 3 − 6, 1 − 4 − 5 − 6 et 1 − 2 − 5 − 6 on obtient un flot de 6 unités. Une modification de l’algorithme s’impose donc. Nous allons permettre d’étiqueter un sommet en utilisant un arc dans le sens inverse, sous la condition que le flot sur cet arc est strictement positif. On peut en effet, comme dans la méthode du simplexe de la section précédente, maintenir la conservation de flot en diminuant le flot des arcs utilisés à rebours. Grâce à cette astuce, on identifie le chemin d’augmentation 1 − 4 − 5 − 2 − 3 − 6 dans la figure 9. Le sommet 2 est étiqueté à partir du sommet 5 et de l’arc (2,5). La règle permettant de déterminer l’augmentation maximale ϵ doit être modifiée pour interdire à un flot qui diminue de devenir négatif. Soit p^{+} l’ensemble des arcs du chemin parcourus dans le sens naturel et p^{−} l’ensemble des arcs du chemin parcourus dans le sens inverse. On a

ϵ_{1} = min_{(i,j)∈p^{+}}{u_{ij} − x_{ij}}
ϵ_{2} = min_{(i,j)∈p^{−}}{x_{ij}}
ϵ = min{ϵ_{1}, ϵ_{2}}.

Le flot optimal est maintenant tel qu’illustré à la figure 10.

[figure: 1; 2; 3; 4; 5; 6; s; t; 4 ≤ 4; 0 ≤ 3; 0 ≤ 2; 1 ≤ 3; 4 ≤ 6; 1 ≤ 2; 5 ≤ 5; v = 4 + 1 = 5; 1; 2; 3; 4; 5; 6; s; t; 4 ≤ 4; 0 ≤ 3; 0 ≤ 2; 0 ≤ 3; 4 ≤ 6; 0 ≤ 2; 4 ≤ 5; v = 4; 1; 2; 3; 4; 5; 6; s; t; 4; 3; 2; 3; 6; 2; 5; v = 0; capacité; flot, capacité; +ϵ; +ϵ; −ϵ; +ϵ; ϵ]

Fig. 9 – Un problème de flot maximal dans un réseau

Remarques

1. On peut associer à chaque sommet j une étiquette f_{j} représentant une augmentation possible du flot jusqu’à j. On initialise y_{s} à = +∞ et y_{j} à −∞ j ≠ s. Lors du balayage du sommet i, on révise les étiquettes des voisins de i :

– y_{j} = max{y_{j}, min{y_{i}, u_{ij} − x_{ij}}} si l’arc (i, j) existe ;

– y_{j} = max{y_{j}, min{y_{i}, x_{ji}}} si l’arc (j, i) existe. Si l’étiquette est modifiée, on note le sommet prédécesseur i afin de pouvoir retracer le chemin.

Ce procédé permet d’identifier le chemin d’augmentation ayant la plus grande capacité résiduelle. Cette stratégie gloutonne n’est pas toujours la meilleure et peut même ne pas converger si les capacités ne sont pas entières ! Une stratégie qui converge dans tous les cas consiste à utiliser un chemin d’augmentation possédant le plus petit nombre d’arcs.

2. Soit S l’ensemble des sommets étiquetés lorsque l’algorithme se termine. Puisque la solution est optimale, t ∉ S, c’est-à-dire : t ∈ T. Définissons :

C = A ∩ (S × T).

Puisque s ∈ S et t ∈ T, l’ensemble C n’est pas vide. Cet ensemble, que l’on appelle une coupe, déconnecte la source s de la destination t. Il s’ensuit que la somme des capacités des arcs de C est une borne supérieure sur le flot maximal. À l’optimum, cette borne est égale à la valeur maximale du flot, et on nomme la coupe correspondante coupe minimale (voir figure 11). Un théorème célèbre, dû à Ford et Fulkerson, affirme que le valeur maximale du flot est égale à la capacité minimale d’une coupe (MAX FLOT = MIN COUPE). Ce résultat établit l’égalité des objectifs primal et dual lorsqu’il existe une solution optimale.

[figure: 1; 2; 3; 4; 5; 6; s; t; 4 ≤ 4; 1 ≤ 3; 1 ≤ 2; 2 ≤ 3; 3 ≤ 6; 2 ≤ 2; 5 ≤ 5; v = 5 + 1 = 6]

Fig. 10 – Le flot optimal

[figure: 1; 4; 2; 3; 5; 6; 4; 2; S; T; C = {(1, 2), (4, 5)}]

Fig. 11 – Une coupe minimale

3. Les plus courts chemins

Le problème de la recherche d’un plus court chemin entre deux sommets s et t d’un réseau est l’un des plus étudiés en recherche opérationnelle. Il peut s’exprimer comme un programme linéaire en variables entières :

min_{x}  ∑_{(i,j)∈A} c_{ij}x_{ij}

∑_{j} x_{ji} − ∑_{j} x_{ij} = \{ \begin{matrix}−1 & si i = s \\ +1 & si i = t \\ 0 & sinon \end{matrix}

x_{ij} ∈ {0, 1}.

Dans le cas où il n’existe pas de circuit dont la somme des longueurs d’arc est négative, on peut remplacer les contraintes x_{ij} ∈ {0, 1} par x_{ij} ≥ 0, tout simplement.

En pratique, il est (presque) aussi simple de calculer les plus courts chemins de s vers tous les autres sommets de N que de s vers t. De plus, comme la solution de ce problème plus général est un arbre de plus courts chemins et donc une solution de base, les contraintes d’intégralité sur les variables ne sont pas requises. On peut donc appliquer l’algorithme du simplexe à ce cas particulier de PFCM. Nous allons plutôt présenter un schéma algorithmique combinatoire très général qui partage certains points communs avec la méthode du simplexe. Ce schéma attribue à chaque sommet i une étiquette y_{i} qui, à la fin de l’algorithme, correspond à la longueur d’un plus court chemin de s à i. Afin de retracer les plus courts chemins, on met à jour le vecteur des prédécesseurs pred.

ALGORITHME GÉNÉRIQUE

1. B ← {s} y_{s} ← 0 y_{i} ← +∞ ∀i ≠ s

2. si B = ∅ STOP sinon choisir i ∈ B

3. BALAYER LE SOMMET i. ∀(i, j) ∈ A faire : si y_{i} + c_{ij} < y_{j} alors

– y_{j} ← y_{i} + c_{ij}
– B ← B ∪ {j}
– pred(j) ← i

4. Retirer i de B et retourner à 2.

Appliquons l’algorithme à l’exemple de la figure 12.

[figure: 5; 2; 2; −2; 1; 3; 3; 1; 2; 3; 4; 5; 6; s; longueur]

Fig. 12 – Réseau (plus court chemins)

• choisir le sommet 1 et le balayer B ← {1} y ← {0, ∞, ∞, ∞, ∞, ∞}

• balayer 1 y_{1} + c_{12} = 0 + 5 < y_{2} = ∞ ⇒ B ← {1, 2}, y_{2} ← 5, pred(2) ← 1 y_{1} + c_{13} = 0 + 1 < y_{3} = ∞ ⇒ B ← {1, 2, 3}, y_{3} ← 1, pred(3) ← 1 B ← B − {1} = {2, 3}

• choisir le sommet 3 et le balayer y ← (0, 5, 1, 4, ∞, ∞), pred ← (−, 1, 1, 3, −, −), B ← {2, 4}

• choisir le sommet 4 et le balayer y ← (0, 5, 1, 4, ∞, 7), pred ← (−, 1, 1, 3, −, 4), B ← {2, 6}

• choisir le sommet 2 et le balayer y ← (0, 5, 1, 3, 7, 7), pred ← (−, 1, 1, 2, 2, 4), B ← {4, 5, 6}

• choisir le sommet 4 et le balayer y ← (0, 5, 1, 3, 7, 6), pred ← (−, 1, 1, 2, 2, 4), B ← {5, 6} (réétiquetage du sommet 6)

• choisir le sommet 5 et le balayer y ← (0, 5, 1, 3, 7, 6), pred ← (−, 1, 1, 2, 2, 4), B ← {6} (aucune modification d’étiquettes ou de prédécesseurs)

• choisir le sommet 6 et le balayer y ← (0, 5, 1, 3, 7, 6), pred ← (−, 1, 1, 2, 2, 4), B ← ∅ (fin de l’algorithme)

Remarques

1. Le choix du sommet à balayer influence l’efficacité de l’algorithme. Le choix du sommet d’étiquette minimale est celui adopté dans l’algorithme de Dijkstra. Si toutes les longueurs d’arc sont non négatives, alors l’étiquette d’un sommet balayé devient permanente. Ce résultat permet de démontrer facilement que la complexité de la méthode est O(n^{2}), où n est le nombre de sommets du graphe.

Si les sommets sont choisis suivant l’ordre de leur apparition dans l’ensemble B, alors on obtient l’algorithme de Ford-Bellman-Moore dont la complexité est O(n^{3}) mais qui est souvent plus efficace que la méthode de Dijkstra en pratique.

2. S’il existe un circuit de longueur négative dans le réseau, alors il n’existe pas d’arbre de plus courts chemins (solution non bornée inférieurement). L’algorithme peut détecter cette situation : si l’étiquette du sommet i est améliorée à une nième itération, ce sommet fait partie d’un circuit négatif qu’on peut retracer à l’aide de la liste des prédécesseurs.

3. Un algorithme permet de résoudre élégamment le problème qui consiste à trouver les plus courts chemins entre tous les couples de sommets d’un réseau. Cet algorithme est inspiré de l’algorithme de Warshall permettant de déterminer la fermeture transitive d’un graphe (revoir les notes du cours IFT 1063). L’algorithme de Floyd-Warshall construit les longueurs w_{ij}^{[k]} des plus courts chemins du sommet i au sommet j n’utilisant comme sommets intermédiaires que des sommets de l’ensemble {1, 2, …, k}. On pose : w_{ij}^{[0]} = c_{ij} et l’on détermine itérativement

w_{ij}^{[k]} = min{w_{ij}^{[k−1]}, w_{ik}^{[k−1]} + w_{kj}^{[k−1]}}.

Par construction, l’algorithme se termine à l’itération n avec la matrice W des plus courtes distances entre tous les couples de sommets. Si l’on souhaite retracer les chemins, il suffit de mettre à jour simultanément une matrice de prédécesseurs. La complexité de cet algorithme est O(n^{3}).

4. Tous les livres de recherche opérationnelle traitent de gestion de projet par la technique du �chemin critique�. Ce triste chapitre de la recherche opérationnelle m’ennuie profondément. C’est pourquoi je n’en parlerai plus...

4. Arbre sous-tendant de poids minimal

Soit le problème consistant à relier au moindre coût un certain nombre de villes par un réseau de fibre optique, sachant que l’installation des câbles se fera le long des routes du réseau routier. Soit c_{ij} le coût de construction unitaire sur l’arc (i, j). La solution du problème sera certainement un arbre, c’est-à-dire un sous-graphe acyclique comportant n−1 arcs. Le �meilleur� arbre sera solution du programme mathématique

min_{x}  ∑_{i,j} c_{ij}x_{ij}

∑_{i} ∑_{j} x_{ij} = n − 1

∑_{(i,j)∈S×S} x_{ij} ≤ |S| − 1, ∀S ⊆ N

x_{ij} ∈ {0, 1},

où |S| représente la cardinalité (nombre d’éléments) de l’ensemble S.

Supposons que nous soyions en possession d’une solution partielle du problème (voir figure 13) qui forme une forêt^{4} .

On appelle admissible un arc qui ne fait pas partie de la forêt et dont l’ajout ne crée pas de cycle. Soit T un arbre de la forêt. Puisque T doit être relié au reste du graphe, il est naturel de considérer l’arc admissible adjacent à T qui soit de coût minimum. À chaque arbre T est ainsi associé un arc intéressant indiqué par un arc pointillé dans la figure 13.

Et maintenant, pourquoi ne pas insérer les arcs intéressants dans la solution ? De façon surprenante, on peut démontrer que cette stratégie �gloutonne� conduit à la solution optimale. Nous allons en considérer trois variantes.

4Une forêt est un graphe acyclique. Une forêt est l’union d’arbres dont certains peuvent ne contenir qu’un seul point et aucun arc.

[figure: 1; 2; 3; 4; 5; 6; 7; 1; 2; 3; 4; 5; 6; 7; 1; 4; 1; 2; 4; 6; 2; 9; 5; 4; 2; T_{1}; T_{2}; T_{3}; coût; 7; forêt {T_{1}, T_{2}, T_{3}}; T_{1} = {1, 2, 3, 5} T_{2} = {4, 6} T_{3} = {7}]

Fig. 13 – Graphe pour le problème de l’arbre sous-tendant minimal

1. À la première itération, insérer dans la solution l’arc de coût minimum. Aux itérations subséquentes, insérer l’arc intéressant adjacent au seul arbre non trivial. Cette stratégie, attribuée à Prim, fait croître un seul arbre dans la forêt.

2. À chaque itération, insérer l’arc admissible de coût minimum. Cette stratégie est attribuée à Kruskal.

3. Insérer simultanément tous les arcs intéressants.^{5} Cet algorithme, proposé en 1926 par Borůvka est en quelque sorte le plus moderne : son auteur aurait-il eu la prémonition du parallélisme ?

Appliquons ces trois algorithmes au problème de la figure 13. Les itérés successifs sont présentés dans les tableaux ci-dessous. On note que les solutions, bien que de même coût total (14), sont différentes. Ceci peut se produire même si les coûts des arcs sont tous distincts.

[figure: 1; 2; 3; T_{1} = {1, 2}; T_{2} = {3}; 1; 1]

Fig. 14 – Mauvais choix dans l’algorithme de Borůvka

5Attention à ne pas créer de cycle en cas d’égalité de coût (voir figure 14 où l’insertion simultanée des arcs (1,3) et (2,3)

crée un cycle.)

itération 1  2  3  4  5  6
arc inséré (1,2) (1,7) (2,3) (2,5) (5,6) (5,4)
alternative (5,6) (2,5) (7,6) (4,6)

(7,6)

coût 1 2 4 4 1 2

Algorithme de Prim

itération 1 2 3 4 5 6
arc inséré (1,2) (5,6) (4,6) (1,7) (6,7) (2,3)
alternative (5,6) (5,4) (2,3)
(1,7)  (2,5)
coût 1 1 2 2 4 4

Algorithme de Kruskal

itération 1  2
arcs insérés (1,2) (2,1) (3,2) (4,5) (5,6) (6,5) (7,1) (2,5)

alternative

alternative  (1,7)  (2,5)
coût 1  –  4  2  1  –  2  4

Algorithme de Borůvka CHAPITRE 7 : PROGRAMMATION DYNAMIQUE 26 novembre 2003

La programmation dynamique est une technique d’exploration des solutions d’un programme mathématique qui, comme un chat, est plus facile à reconnaître qu’à décrire. Je vais en donner quatre exemples classiques.

1. Le problème du sac-à-dos

Dans cette section, le problème du sac-à-dos

max  ∑_{j=1}^{n} c_{j}x_{j}

∑_{j=1}^{n} w_{j}x_{j} ≤ b

x_{j} entier non négatif

est résolu par une technique différente de celle (Branch-and-Bound) étudiée précédemment. À l’itération k du processus, on détermine une solution optimale correspondant à un sac-à-dos de capacité k. Si z(k) est la valeur de la solution optimale, on obtient trivialement l’équation de récurrence

z(k) = max_{j:w_{j}≤k}{c_{j} + z(k − w_{j})}

avec, comme conditions initiales :

z(k) = 0 si k < min_{j}{w_{j}}.

À chaque itération, on note j(k) un objet pour lequel le maximum a été atteint. Appliquons la technique au problème

max 2x_{1} + 3x_{2} + 7x_{3}
3x_{1} + 4x_{2} + 8x_{3} ≤ 14
x_{1},  x_{2},  x_{3} entiers non négatifs.
k z(k)  j(k)
0 0  –
1 0  –
2 0  –
3 2  1
4 max{2 + z(4 − 3), 3 + z(4 − 4)} = 3  2
5 max{2 + z(5 − 3), 3 + z(5 − 4)} = 3  2
6 max{2 + z(6 − 3), 3 + z(6 − 4)} = 4  1
7 max{2 + z(4), 3 + z(3)} = 5  1,2
8 max{2 + z(5), 3 + z(4), 7 + z(0)} = 7  3
9 max{2 + z(6), 3 + z(5), 7 + z(1)} = 7  3
10 max{2 + z(7), 3 + z(6), 7 + z(2)} = 7  1,2,3
11 max{2 + z(8), 3 + z(7), 7 + z(3)} = 9  1,3
12 max{2 + z(9), 3 + z(8), 7 + z(4)} = 10  2,3
13 max{2 + z(10), 3 + z(9), 7 + z(5)} =10  2,3
14 max{2 + z(11), 3 + z(10), 7 + z(6)} = 11 1,3

On retrace la solution optimale en procédant à rebours : j(14) = 1, j(14 − w_{1}) = j(11) = 1, j(11 − w_{1}) = j(8) = 3. La solution optimale est donc x^{∗} = (2, 0, 1).

Cet algorithme permet non seulement d’obtenir la solution optimale du problème de départ, mais également les solutions optimales correspondant à tous les sacs-à-dos de capacité inférieure ou égale à 14. Notons qu’il est possible d’améliorer l’efficacité de l’algorithme.

2. Les plus courts chemins

Soit t un sommet destination et y_{i}^{[k]} la distance d’un plus court chemin du sommet i au sommet t qui utilise au plus k arcs. La formule récursive suivante permet de calculer y_{i}^{[k]}^{1}

y_{i}^{[k]} = min_{j}{c_{ij} + y_{j}^{[k−1]}}

et on note j_{i}^{[k]} un sommet pour lequel le minimum est atteint. Les condition initiales sont

y_{i}^{[0]} = \{ \begin{matrix}0 si i = t \\ ∞ sinon. \end{matrix}

À moins qu’il n’existe des circuits négatifs, l’algorithme basé sur cette équation récursive trouvera après au plus n − 1 itérations (n désigne le nombre de sommets dans le réseau) un arbre de plus courts chemins de tout sommet i vers la destination t. En effet, s’il existe un plus court chemin, il doit exister un plus court chemin simple (sans circuit) et un tel chemin ne peut comporter plus de n − 1 arcs.

Appliquons cette technique au réseau ci-dessous, avec t = 6.

[figure: 4; 1; 6; 3; 2; 5; 1; 2; 3; 4; 5; 6; longueur; −1; t; 2; 0; −3]

Fig. 1 – Le réseau de l’exemple de plus courts chemins

Les résultats sont indiqués au tableau 1. Le premier nombre indique la longueur d’un chemin partiel. Le nombre entre parenthèses correspond au successeur succ^{[k]}(i) du sommet i dans le chemin obtenu à l’itération k, et permet de retracer le chemin optimal (1,3,2,4,5,6). Je donne au-dessous du tableau deux exemples de calcul de y_{i}^{[k]} et succ^{[k]}(i) :

5  4  3  2  1  0 ← itérations
1 7(3) 8(3) 9(3) ∞ ∞ ∞
2 5(4) 5(4) 5(4) 7(4) ∞ ∞
sommets → 3 4(2) 4(2) 5(4) 6(5) ∞ ∞
4 3(5) 3(5) 3(5) 3(5) 5(6) ∞
5 6(6) 6(6) 6(6) 6(6) 6(6) ∞
6 0(6) 0(6) 0(6) 0(6) 0(6) 0(6)

Tab. 1 – Programmation dynamique pour les plus courts chemins : y_{i}^{[k]} (succ^{[k]}(i))

1Le minimum est pris sur tous les indices de sommets adjacents à i, c’est-à-dire pour lesquels il existe un arc de i à j. Par convention, on suppose que i est relié à lui-même par un arc de longueur nulle.

y_{2}^{[2]} = min{c_{22} + y_{2}^{[1]}, c_{24} + y_{4}^{[1]}, c_{25} + y_{5}^{[1]}}
= min{0 + ∞, 2 + 5, 1 + 6} = 7
succ^{[2]}(2) = 4 ou 5
y_{1}^{[4]} = min{c_{11} + y_{1}^{[3]}, c_{12} + y_{2}^{[3]}, c_{13} + y_{3}^{[3]}}
= min{0 + 9, 4 + 5, 3 + 5} = 8
succ^{[4]}(1) = 3

Remarques

– Si, lors d’une itération complète, il n’y a aucune amélioration, alors on peut arrêter l’algorithme.

– Si on poursuit l’algorithme après la (n − 1)ième itération et qu’il y a amélioration, alors le réseau doit forcément comporter un circuit négatif, c’est-à-dire un circuit dont la somme des coûts des arcs est négative. L’algorithme détectera ce circuit.

– La complexité de l’algorithme est O(n^{3}).

– Noter la similarité de cet algorithme avec l’algorithme générique du chapitre 6 où l’on choisit les sommets à balayer dans l’ordre de leur insertion dans l’ensemble B. En fait, si on autorise l’ensemble B à contenir plusieurs copies d’un même sommet^{2} alors les deux algorithmes sont identiques.

– Le problème du sac-à-dos peut se formuler sous forme d’un problème de plus court chemin en introduisant b + 1 sommets indexés par les entiers de 0 à b et des arcs (k, k + w_{j}) de coût c_{j}. L’algorithme de programmation dynamique tient compte de la structure particulière de ce graphe (graphe topologique, ça devrait vous rappeler des bons souvenirs du cours IFT 1063 : ordres partiels, diagrammes de Hasse, etc.). Ainsi, à l’itération k, la quantité z(k) (qui joue le rôle de y_{i}^{[k]}) est-elle fixée de façon définitive.

3. Gestion de stocks

Considérons un problème de gestion de stocks sur horizon infini où tout (demande par unité de temps d, coût unitaire d’inventaire h, coût fixe de commande c) est constant. On souhaite déterminer le niveau et l’intervalle optimal des commandes. Il est facile de démontrer qu’une stratégie optimale consiste à commander, à intervalle fixe de durée I, une quantité x du produit. Le niveau de l’inventaire est alors donné par la fonction linéaire par morceaux illustrée à la figure 2.

[figure: temps; stock; x; I = x/d]

Fig. 2 – Solution d’un problème simple de gestion de stocks

La demande étant uniforme, le niveau de stock s’annule x/d unités de temps après la commande. Puisque le niveau moyen du stock est égal à x/2, le coût moyen unitaire d’inventaire est égal à hx/2. Le coût unitaire moyen du système prend la forme

2... ce qui est inefficace, car un sommet i pourrait être balayé avec une valeur y_{i} qui est périmée.

f(x) = \frac{c}{x/d} + h\frac{x}{2} = \frac{cd}{x} + \frac{h}{2}x.

Le minimum de cette fonction est atteint lorsque la dérivée f′(x) = −\frac{cd}{x^{2}} − \frac{h}{2} s’annule, c’est-à-dire

x = \sqrt{\frac{2cd}{h}}.

Pas grand-chose à voir avec la programmation dynamique. Cependant, les choses deviennent plus intéressantes et plus réalistes lorsque les paramètres varient, par exemple s’il y a des économies d’échelle ou des variations saisonnières de la demande. Pour simplifier un peu, on suppose que les commandes sont passées à intervalles réguliers et que l’horizon de planification est fini. On utilise les notations suivantes :

[k − 1, k]  période k (k = 1, …, n)
s_{k}  niveau de l’inventaire à la fin de la période k (état du système )
c_{k}(x)  coût de commande pour une quantité x
h_{k}(s)  coût d’inventaire pour un niveau moyen s
d_{k}  demande à la période k
f^{[k]}(s)  coût d’une politique optimale de commande débutant en début de période k (instant k − 1) à l’état s
x^{[k]}(s)  commande optimale au début de la période k

L’évolution de l’inventaire est illustrée à la figure 3.

[figure: 1; 0; 2; 3; 4; k; stock; s_{0}; s_{1}; s_{2}; s_{3}; 1; 2; 3; 4; 5; s_{4}; s_{0} + x_{1}; s_{1} + x_{2}; s_{2} + x_{3}; s_{3} + x_{4}]

Fig. 3 – Évolution typique du niveau d’inventaire

L’équation de la programmation dynamique pour ce problème prend la forme^{3}, pour toute valeur possible du niveau de stock s :

f^{[n+1]}(s) = 0 ∀s

f^{[k]}(s) = min_{x admissible} \{c_{k}(x) + h_{k} × \frac{s + x + s_{k}}{2} + f^{[k+1]}(s_{k})\} k = n, n − 1, …, 2, 1

avec l’équation d’évolution s_{k} = s + x − d_{k}. Dans l’équation ci-dessus, (s + x + s_{k})/2 représente le niveau de stock moyen en période k.

3Une commande x est admissible si elle respecte la contrainte de capacité et si elle permet de satisfaire la demande pour la période k.

Ces équations de récurrence donnent, pour chaque période k et chaque niveau d’inventaire s, le coût d’une politique de commande optimale du début de la période k jusqu’à la fin de l’horizon de planification. Dans le cas où tous les paramètres et variables sont entiers, ces équations sont faciles à résoudre.

Considérons l’exemple sur 3 périodes correspondant aux données suivantes : h_{k}(s) = 2s, vecteur de demande= (2,1,4), capacité de stockage = 5 et c_{k}(x) = c(x) avec comme fonction de coût de commande

\begin{matrix}x & 0 & 1 & 2 & 3 & 4 & 5 \\ c(x) & 0 & 10 & 16 & 20 & 22 & 23 & . \end{matrix}

On suppose que le stock initial est nul. Une solution optimale est donnée au tableau 2. Le premier chiffre indique le coût de la solution partielle. Il est suivi, entre parenthèses, du niveau de commande optimal x^{[k]}(s). Il est facile de démontrer que le niveau de stock terminal est nul. Le niveau de commande de la dernière période (période 3) s’obtient donc trivialement. Aux autres périodes (2 et 1) on obtient successivement :

s f^{[1]}(s) f^{[2]}(s) f^{[3]}(s) f^{[4]}(s)
5  13(0)  6(0)  –
4  21(0)  4(0)  0
3  25(0) 14(1)  0
2  27(0) 20(2)  0
1  27(0) 24(3)  0
0 51(3) 36(5) 26(4)  0

Tab. 2 – Solution optimale du problème de gestion de stock : f^{[k]}(s_{k}) (x^{[k]}(s))

Voici quelques exemples de calcul des valeurs optimales :

f^{[2]}(0) = min\{c(1) + h (\frac{1+0}{2}) + f^{[3]}(0 + 1 − 1),  x optimal
c(2) + h (\frac{2+1}{2}) + f^{[3]}(0 + 2 − 1),
c(3) + h (\frac{3+2}{2}) + f^{[3]}(0 + 3 − 1),
c(4) + h (\frac{4+3}{2}) + f^{[3]}(0 + 4 − 1),
c(5) + h (\frac{5+4}{2}) + f^{[3]}(0 + 5 − 1)\}
= min{10 + 1 + 26, 16 + 3 + 24, 20 + 5 + 20, 22 + 7 + 14, 23 + 9 + 4} = 36  5
f^{[2]}(1) = min{0 + 1 + 26, 10 + 3 + 24, 16 + 5 + 20, 20 + 7 + 14, 22 + 9 + 4} = 27  0
f^{[2]}(2) = min{0 + 3 + 24, 10 + 5 + 20, 16 + 7 + 14, 20 + 9 + 4} = 27  0
f^{[2]}(3) = min{0 + 5 + 20, 10 + 7 + 14, 16 + 9 + 4} = 25  0
f^{[2]}(4) = min{0 + 7 + 14, 10 + 9 + 4} = 21  0
f^{[2]}(5) = min{0 + 9 + 4} = 13  0
f^{[1]}(0) = min{16 + 2 + 36, 20 + 4 + 27, 22 + 6 + 27, 23 + 8 + 25} = 51  3

Dans ces équations, le premier nombre représente le coût de commande, le second le coût d’inventaire et le troisième le coût futur. La somme des deux premiers coûts représente le coût immédiat. Noter que le minimum n’est pris que sur les valeurs admissibles du niveau de commande : il faut satisfaire la demande à venir sans dépasser la capacité de l’entrepôt (5). Ainsi, si le niveau d’inventaire est nul au début de la deuxième période, il faut commander au moins une unité pour satisfaire à la demande de cette période.

On trouve la solution optimale en partant de la cellule de gauche, en se dirigeant vers la droite et en respectant l’équation d’évolution de l’inventaire. Ainsi, x_{1} = 3, x_{2} = 0 et x_{3} = 4. La mise en œuvre informatique de cet algorithme est d’une simplicité déconcertante.

On remarque une similitude entre cet algorithme et l’algorithme de programmation dynamique présenté à la section précédente (plus courts chemins). En fait, ces algorithmes sont identiques. On peut en effet formuler le problème d’inventaire comme un problème de plus court chemin en en créant, à chaque début de période, des nœuds correspondant aux niveaux d’inventaire. À chaque niveau de commande x correspond un arc liant un sommet (k, s) au sommet (k + 1, s_{k} = s + x − d_{k}) de coût c_{k}(x_{k}) + h_{k}((s + x + s_{k})/2).

4. Un modèle d’achats stochastique

Les situations où tous les paramètres d’un problème sont connus avec certitude sont peu fréquentes. Par exemple, dans le modèle précédent, il est sans doute irréaliste de faire l’hypothèse que la demande est une constante immuable. Dans cette section, nous considérons la situation où un détaillant achète un produit d’un grossiste pour le revendre. Malheureusement, le commerçant doit passer ses commandes au début du mois, avant de connaître le niveau de la demande mensuelle. Par expérience, il en connaît cependant la distribution de probabilité. S’il achète une trop grande quantité du produit et que la demande est faible, les coûts d’inventaire seront élevés et il risque de se retrouver avec des invendus. S’il achète peu du produit, il risque de perdre des ventes. Notre algorithme de programmation dynamique devra donc lui suggérer une politique d’achat qui maximise son gain moyen, étant donné un niveau de stock au début de chaque mois. On introduit les notations suivantes :

r_{k}(v)  revenu au mois k pour un niveau de vente v
Prob (d_{k} = d) probabilité que le niveau de demande soit égal à d
v_{k}  niveau de vente au mois k

On suppose que toutes les quantités sont entières et non négatives. On a

Prob (d_{k} = d) ≥ 0

∑_{d} Prob (d_{k} = d) = 1.

Il faut modifier l’équation d’évolution pour interdire au niveau de stock de devenir négatif :

s_{k} = max{s_{k−1} + x_{k} − d_{k}, 0}

et introduire une équation spécifiant le niveau de vente v_{k} au mois k correspondant à une demande d :

v_{k} = min{d_{k}, s_{k−1} + x_{k}}.

Au début du mois k, le gain maximal s’obtient en résolvant l’équation de récurrence probabiliste

f^{[k]}(s) = max_{x admissible} ∑_{d} Prob (d_{k} = d) [r_{k}(v_{k}) − c_{k}(x) − h_{k} × (\frac{s + x + s_{k}}{2} ) + f^{[k+1]}(s_{k})] .

Une commande x est admissible si elle respecte la capacité de l’entrepôt (s_{k−1} + x ≤ capacité).

À chaque itération, on maximise le gain moyen jusqu’à la fin de la période de planification. Il est important de noter que la politique optimale n’est pas fixée au début de la période de planification mais s’ajuste, mois après mois, aux aléas du marché de la demande. Noter que v_{k} et s_{k} sont des fonctions de d.

Il est maintenant plus que temps de considérer un exemple concret dont les paramètres sont

– Stock initial nul.

– Prix de vente unitaire : 4.

– Horizon de planification : 3 mois.

– Capacité de l’entrepôt : 3.

– Coût unitaire d’inventaire : 2.

– Répartition de la demande et coûts de commande (c_{k}(x) = c(x) pour tout k)

\begin{matrix}d & 0 & 1 & 2 & 3 \\ Prob (d_{k} = d) & \frac{1}{8} & \frac{2}{8} & \frac{3}{8} & \frac{2}{8} \\ x & 0 & 1 & 2 & 3 \\ c(x) & 0 & 2 & 3 & 3 \end{matrix}

Comme les calculs sont fastidieux, je me contenterai de donner quelques éléments du tableau contenant les f^{[k−1]}(s_{k−1}) (voir tableau 3) et d’illustrer le calcul de certains de ses éléments. Je vous encourage fortement à écrire un petit programme qui vous permette de compléter ce tableau.

s f^{[1]}(s) f^{[2]}(s) f^{[3]}(s) f^{[4]}(s)
3  ? \frac{22}{8}(0)  0
2  ? \frac{28}{8}(0)  0
1  \frac{177}{64}(1) \frac{19}{8}(0)  0
0  ? \frac{64}{64}(1)  \frac{3}{8}(1)  0

Tab. 3 – Solution optimale partielle du problème de gestion des achats stochastique

Le calcul de l’élément f^{[3]}(1) s’effectue de la façon suivante, en notant qu’on ne peut satisfaire une demande supérieure à 1 si le niveau de stock est nul et que l’on ne commande rien. Il suffit donc de ne considérer, dans le calcul du premier élément du gain moyen, que deux cas : soit la demande est nulle (avec probabilité 1/8) soit la demande est supérieure ou égale à 1 (avec probabilité 7/8). Si l’on commande une unité, il y a trois possibilités de vente (0, 1 ou 2). Si l’on commande 2 unités, on pourra répondre à toute demande pendant le troisième mois. Comme dans le cas déterministe, le gain comporte deux composantes : un gain immédiat (revenu de la vente − coût de commande − coût d’inventaire) et un gain à venir. Ceci nous donne :

f^{[3]}(1) = max\{\frac{1}{8}(0 − 0 − 2) + \frac{7}{8}(4 − 0 − 1),  x = 0
\frac{1}{8}(0 − 2 − 4) + \frac{2}{8}(4 − 2 − 3) + \frac{5}{8}(8 − 2 − 2),  x = 1
\frac{1}{8}(0 − 3 − 6) + \frac{2}{8}(4 − 3 − 5) + \frac{3}{8}(8 − 3 − 4), +\frac{2}{8}(12 − 3 − 3)\} x = 2
= max\{\frac{19}{8}, \frac{12}{8}, \frac{4}{8}\} = \frac{19}{8}  x^{[3]}(1) = 0.

Si le niveau de stock au début du troisième mois est égal à 2, il est inutile de considérer l’achat de plus d’une unité du produit :

f^{[3]}(2) = max\{\frac{1}{8}(0 − 0 − 4) + \frac{2}{8}(4 − 0 − 3) + \frac{5}{8}(8 − 0 − 2),  x = 0
\frac{1}{8}(0 − 2 − 6) + \frac{2}{8}(4 − 2 − 5) + \frac{3}{8}(8 − 2 − 4), +\frac{2}{8}(12 − 2 − 3)\} x = 1
= max\{\frac{28}{8}, \frac{6}{8}\} = \frac{28}{8}  x^{[3]}(2) = 0.

Et ainsi de suite :

f^{[3]}(3) = \frac{1}{8}(0 − 0 − 6) + \frac{2}{8}(4 − 0 − 5) + \frac{3}{8}(8 − 0 − 4) + \frac{2}{8}(12 − 0 − 3)\} = \frac{22}{8} x = 0  x^{[3]}(3) = 0.

f^{[2]}(1) = max\{\frac{1}{8}(0 − 0 − 2 + \frac{19}{8}) + \frac{7}{8}(4 − 0 − 1 + \frac{3}{8}),  x = 0
\frac{1}{8}(0 − 2 − 4 + \frac{28}{8}) + \frac{2}{8}(4 − 2 − 3 + \frac{19}{8}) + \frac{5}{8}(8 − 2 − 2 + \frac{3}{8}),  x = 1
\frac{1}{8}(0 − 3 − 6 + \frac{22}{8}) + \frac{2}{8}(4 − 3 − 5 + \frac{28}{8}) + \frac{3}{8}(8 − 3 − 4 + \frac{19}{8}) + \frac{2}{8}(12 − 3 − 3 + \frac{3}{8})\} x = 2
= max\{\frac{172}{64}, \frac{177}{64}, \frac{125}{64}\} = \frac{177}{64}  x^{[2]}(1) = 1.

Il existe plusieurs variantes et extensions de ce modèle stochastique. Par exemple, on pourrait attribuer une valeur aux invendus à la fin de l’horizon de planification. On pourrait également pénaliser l’entreprise en cas de rupture de stock. Ou considérer une demande qui varie dans le temps (demande saisonnière).

La programmation dynamique

Tel que promis en début de section, je n’ai pas donné de définition formelle de la programmation dynamique, comme je l’avais fait, par exemple, pour la programmation linéaire. On peut cependant identifier deux points fondamentaux communs aux situations considérées, qui permettent d’utiliser un schéma algorithmique relativement général :

– une division de l’horizon de planification en périodes ;

– une fonction objectif additive qui, au début d’une période donnée, est la somme de deux termes : un gain passé et un gain à venir. Cette forme fonctionnelle permet d’optimiser en ne tenant compte que de l’état présent, indépendamment des décisions qui ont mené à cet état.

La définition des périodes et états est souvent naturelle. Parfois, cependant, un effort d’imagination est nécessaire pour mettre le problème sous une forme qui se prête à la programmation dynamique. C’est le cas du problème de sac-à-dos où il n’y a pas d’aspect dynamique � naturel �. Dans cet exemple � dégénéré �, les périodes correspondent à la capacité résiduelle du sac et il n’y a pas vraiment d’états. On peut en fait considérer qu’au début de la période, il n’existe qu’un seul état k correspondant à un sac-à-dos vide, et qui ne joue donc aucun rôle.

En conclusion, je vais élaborer brièvement sur les avantages et inconvénients de la programmation dynamique.

Avantages

– Très flexible. Il est facile d’ajouter des contraintes, d’incorporer des fonctions quelconques.

– Fournit une solution globalement optimale.

– Fournit non seulement une solution optimale, mais un ensemble de solutions optimales à chaque étape intermédiaire et pour tous les états possibles. Ainsi, il est possible d’utiliser les résultats dans les situations où des modifications du système mènent à des changements d’état impromptus.

– Permet de traiter des problèmes stochastiques.

– S’adapte à des situations où le temps, les états et les décisions varient continûment. Mais là, c’est pas mal plus complexe.

Inconvénients

– Ne s’applique que sous les hypothèses d’additivité et d’amnésie (le passé n’influence pas le futur).

– Parfois lent et donc inefficace.

– Ne se prête pas bien à la réoptimisation, contrairement à la programmation linéaire.