Recherche Opérationnelle
Workshop

Guide tuteur

Version	Date	Concepteur	Relecteur	Commentaires
1.4	2022-09-14	Benjamin COHEN BOULAKIA		Réorganisation du Workshop sous forme de points de discussion avec les étudiants Quelques explications plus détaillées
1.3	2022-03-15	Benjamin COHEN BOULAKIA	Jean-François DOLLINGER	Simplification de la partie Résolution graphique Quelques remaniements et explications plus détaillées
1.2	2021-05-10	Benjamin COHEN BOULAKIA		La conclusion ouvre vers l'optimisation de programmes linéaires en nombres entiers.
1.1	2020-06-22	Benjamin COHEN BOULAKIA		La partie 1.2 Résolution par la méthode graphique a été simplifiée et raccourcie, le code à compléter est moins important. Certaines parties à compléter sont moins incomplètes.
1.0	2019-25-04	Benjamin COHEN BOULAKIA		Version initiale. Plusieurs exercices repris du Workshop RO de l'ancien bloc Aide à la décision

Compétences du Workshop

• Associer des situations réelles à des problèmes types de la RO
• Modéliser un problème d'optimisation sous forme de programme linéaire
• Identifier les cas d'usage de l'algorithme du Simplexe et de la programmation linéaire en nombres réels
• Expliquer la difficulté des programmes linéaires en nombres entiers/mixtes
• Calculer une solution à un programme linéaire en nombres réels en appliquant une méthode de résolution graphique
• Utiliser une bibliothèque de calcul scientifique pour résoudre un programme linéaire en nombres réels

Objectifs

Le but de ce Workshop est de découvrir la modélisation par Programme Linéaire, le principe de sa résolution en nombres réels, et les limites de cette approche, notamment lorsque les variables de décision sont entières.

Philosophie et recommandations

Il faut veiller, lors du premier exercice, à ce que les étudiants appliquent bien la méthode graphique. Ils pourraient être tentés de se lancer dans une résolution analytique du problème. Dans ce cas particulier, avec peu de variables et de contraintes, cela marcherait, mais ce n'est pas généralisable. S'ils veulent aller dans cette direction, ils doivent comprendre que ce premier exercice n'a pour but que de leur faire comprendre les propriétés des programmes linéaires, sur lesquelles s'appuient les méthodes de résolution qu'ils découvriront plus loin dans le Workshop, et qui, elles, savent gérer un nombre quelconque de dimensions.

La deuxième partie du Workshop doit les amener à comprendre pourquoi la résolution par le simplexe (ou toute autre méthode de résolution d'un PLNR) ne fonctionne pas avec des contraintes d'intégrité. Une ouverture vers l'optimisation de programmes linéaires en nombre entiers est présentée en conclusion, mais le but est surtout de connaître l'existence de ces méthodes, et de conclure que leur complexité les exclut de cette séquence pédagogique et du projet sur lequel ils travaillent.

Ce Workshop aborde l'une des méthodes de résolution les plus utilisées dans la Recherche Opérationnelle, la Programmation Linéaire. Nous allons nous concentrer sur la modélisation par programmation linéaire, ses différentes variantes, et différentes méthodes de résolutions avec leurs limites.

Vous aurez un certain nombre d'équations à écrire. Rappel : vous pouvez insérer dans une cellule Markdown des équations au format $\mathrm{\LaTeX}$ en les plaçant entre \$. Par exemple, le code

$$\sum_{i=1}^n n$$

produira l'expression mathématique
$$\sum_{i=1}^n n$$

1. Programme linéaire en nombre réels à 2 dimensions

Une multinationale de bonbons pour les enfants dispose de 2 arômes de base pour la fabrication de ses produits. La saveur de chaque bonbon dépend de la quantité de chaque arôme.
Le premier ingrédient coûte à l'entreprise 10€/kg et le deuxième 20€/kg. Le coût de fabrication d'un kilo de produit (somme du mélange de 2 ingrédients) est de 5€. Les études de gôut auprès des consommateurs a montré que la composition du mélange ne doit pas être supérieure à 50% du premier ingrédient et 80% du deuxième.
Par ailleurs, le mois dernier, l'entreprise a vendu 100 Kg à 50€/kg. L'entreprise ne souhaite pas produire plus que ce qu'elle arrive à vendre. Combien de kilos de bonbons doit-elle fabriquer par mois et les proportions de chaque ingrédient pour maximiser les bénéfices ?

1.1 Modélisation du problème

Commençons par modéliser le problème sous la forme d'un programme linéaire.

Variables de décision du programme

Identifiez les inconnues qui sont à déterminer dans le problème. On les nommera avec des variables $X_1$, $X_2$,…

Solution $X_1$ : quantité d'arôme 1
$X_2$ : quantité d'arôme 2

Contraintes du programme

Quelles contraintes identifiez-vous dans le texte ?

Solution

Les contraintes sont des conditions sur l'usage des ressources, c'est à dire des conditions que doivent respecter les valeurs des variables de décision pour former une solution admissible. Il y en a trois explicites :

• La quantité de demande maximale par mois
• La règle de composition _inférieur à 50% de l'arôme 1
• La règle de composition _inférieur à 80% de l'arôme 2

Il y a aussi une règle implicite, mais qu'il est important d'exprimer :

• Contraintes de non-négativité (nombre de kgs/mois positifs)

Modélisez toutes ces contraintes sous formes d'inéquations sur une ou plusieurs variables de notre problème.

Solution

Dans le cas d'un programme linéaire, ces contraintes sont elle-mêmes linéaires, Il est important de fixer les unités (temps, monétaire, …) et de les respecter dans tout le programme linéaire.

• Quantité de demande maximale par mois :
  $X_1 + X_2 \leq 100$
  
  
• Règle de composition inférieur à 50% de l'arôme 1 :
  $\frac{X_1}{X_1 + X_2} \leq 0.5$, reformulé $0,5X_1 – 0,5X_2 \leq 0$
  
  
• Règle de composition inférieur à 80% de l'arôme 2 :
  $\frac{X_2}{X1 + X2} \leq 0.8$, reformulé $-0,8X_1 + 0,2X_2 \leq 0$
  
  
• Contraintes de non-négativité :
  $X_1, X_2 \geq 0$

Fonction économique

Exprimez la fonction économique du problème.

Solution

La fonction économique décrit ce qu'on veut optimiser lors de l'usage de nos ressources. C'est donc une fonction de nos variables de décision, que l'on doit maximiser ou minimiser. Comme pour les contraintes, dans un programme linéaire, cette fonction économique doit être elle-même linéaire.

Ici, la fonction économique doit maximiser le bénéfice (prix de vente des bonbons - coût des matières premières et de fabrication), soit :
$F (X_1, X_2) = 50\times (X_1+X_2)–10\times X_1–20\times X_2-5\times (X_1+X_2)$
reformulée
$F (X_1, X_2) = 35X_1 + 25X_2$

1.1.4 Programme linéaire complet

Réunissez les réponses du dessus pour exprimer le programme linéaire complet optimisant la composition des bonbons de l'usine. Vous pouvez compléter la formule $\mathrm{\LaTeX}$ ci-dessous ($\text{s.c.}$ veut dire sous contraintes) :

$\begin{array}{rl@{}ll} \text{Minimiser/Maximiser} & … &\\ \text{s.c.}& … & \leq…\\ & … & \geq …\\ & … \end{array}$

Solution

$\begin{array}{rl@{}ll} \text{Maximiser} & 35X_1 + 25X_2 &\\ \text{s.c.}& X_1 + X_2 & \leq 100\\ & 0.5X_1 - 0,5X_2 & \leq 0\\ & -0,8X_1 + 0,2X_2 & \leq 0\\ & X_1 & \geq 0\\ & X_2 & \geq 0 \end{array}$

1.2 Résolution par la méthode graphique

Pour comprendre le comportement de ce genre de système, on va regarder comment il se comporte graphiquement. Chaque variable sera représentée sur un axe du graphique. Ici, le problème a deux variables, il sera représenté sur un plan. Nous allons prendre la convention suivante : $X_1$ sera représenté sur l'axe des abscisses et $X_2$ sera représenté sur l'axe des ordonnées.

Explications Évidemment, cette méthode graphique ne peut être appliquée que sur un programme linéaire comportant au plus trois variables (et encore, cela reste nettement plus lisible avec deux variables). Au-delà, la représentation et la résolution graphiques deviennent impossibles (imaginez la représentation graphique d'un problème à 4 ou 5 dimensions).

Dans le cas qui nous occupe, on pourrait trouver la solution au problème en résolvant ce système d'inéquations, voire en cherchant de manière empirique, mais ce n'est pas le but. On essaye de comprendre ce qu'implique le fait d'avoir des contraintes et une fonction économique linéaires.

Nous allons représenter sur ce plain les inéquations de contraintes et la fonction économique. Pour cela, nous utiliserons matplotlib.

Représentation des contraintes

Commençons par placer les différentes contraintes

Déjà, comment tracer une droite en Python ? C'est très simple avec la fonction plot de la bibliothèque Matplotlib, surtout si on utilise NumPy pour manipuler des intervalles de valeurs linéaires. Il suffit d'utiliser la fonction fill_between.

Explications

Le but ici est de tracer une inéquation, il faudrait donc faire ressortir la zone des valeurs respectant l'inégalité. Pour les inéquations du type $y\geq\dots$, il suffit pour cela d'utiliser fill_between(https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.fill_between.html) (en laissant pour le paramètre y2 la valeur pas défaut 0). La partie respectant l'inégalité reste en blanc.

Avant de l'utiliser pour afficher notre programme linéaire, testons ça sur un petit exemple. Mettons qu'on veuille représenter l'inéquation $2x + y \leq 4$

Explications

Jupyter est parfois un peu capricieux avec Matplotlib, si tout ce que les élèves voient, c'est <Figure size 640x480 with 1 Axes>, il faut réexécuter la cellule.

import numpy as np
import matplotlib.pyplot as plt

xmax = 5
x = np.linspace(0, xmax) # on genere l'intervalle des x
y = 4 - 2*x # on genere l'intervalle des y

ymax = np.amax(y)

plt.plot(x, y)
plt.fill_between(x, y, color='lightgrey')

# on dessine entre 0,0 et les extremas des fonctions
plt.xlim(0, xmax)
plt.ylim(0, ymax)

# on affiche
plt.show()

Si on doit tracer une inégalité du type $y\leq\dots$, il suffit de spécifier une valeur pour le paramètre y2 de la fonction fill_between. La valeur de ymax permettra de couvrir tout le graphique.

ymax = np.amax(y)

plt.plot(x, y)
plt.fill_between(x, y, ymax, color='lightgrey')

# on dessine entre 0,0 et les extremas de la fonction
plt.xlim(0, xmax)
plt.ylim(0, ymax)

# on affiche
plt.show()

On a donc tout ce qu'il nous faut pour afficher les inéquations de nos contraintes. Traçons les droites en rouge, et les zones d'inégalités en gris.

Explications

les contraintes de non négativité pourront être simplement gérées en définissant les limites d'axes, et on ira jusqu'à 100, quitte à ajuster si nécessaire). Rappelons que $X_1$ est sur l'axe des abscisses et $X_2$ sur l'axe des ordonnées.

xmax = 100
x = np.arange(0, xmax) # on genere l'intervalle des x. Pour cela, on utilise la version NumPy de range(),
                       # pour pouvoir effectuer facilement les calculs de Y

y1 = 100-x # opération possible avec la version NumPy de range()
y2 =\
     x   #SOLUTION
y3 =\
     4*x #SOLUTION

ymax = \
      max(max(y1), max(y2), max(y3)) #SOLUTION

plt.plot(x, y1, color='red')
plt.plot(\
          x, y2, color='red') #SOLUTION
plt.plot(\
          x, y3, color='red') #SOLUTION

plt.fill_between(x, y1, ymax, color='lightgrey')
plt.fill_between(\
                  x, y2, color='lightgrey')      #SOLUTION
plt.fill_between(\
                  x, y3, ymax, color='lightgrey')#SOLUTION

# on dessine entre 0,0 et les extrema
plt.ylim(0, ymax)
plt.xlim(0, xmax)

# on affiche
plt.gcf().set_size_inches(18, 10) # cette fois on affiche en grand
plt.show()

Que représente la partie apparaissant en blanc ? Quelle est sa caractéristique mathématique ?

Solution

C'est l'espace des solutions admissibles. Toutes les valeurs de $(X_1, X_2)$ respectant les contraintes sont dans cet espace. C'est une surface convexe.

Représentation de la fonction économique

Nous allons pouvoir utiliser cette propriété pour définir la solution optimale. Il nous faut donc maintenant représenter la fonction économique. Sa forme générale est $ax +by = p$. Cette fonction est la famille des droites correspondant à l’ensemble–solution avec une pente $\dfrac{-a}{b}$. Dans notre cas, que vaut cette pente ?

Solution

La fonction économique est de la forme $F (X_1, X_2) = 35X_1 + 25X_2$. La pente vaut donc $\dfrac{-35}{25}=\dfrac{-7}{5}$.

Le but est de trouver la fonction linéaire de $p$ qui maximise la fonction économique tout en restant dans l'espace des solutions admissibles. Quelle est la forme générale de cette fonction ? Quelle valeur devons-nous donc déterminer ?

Solution

$y=x\times p+n$ avec $p=\dfrac{-7}{5}$, soit $y=x\times \dfrac{-7}{5}+n$
Le but est donc de déterminer le $n$ maximisant la valeur de la fonction économique dans le polygone.

Un Widget intégré automatiquement à notre fonction de calcul nous permettra de faire varier facilement ce paramètre et de déterminer les valeurs admissibles de $X_1$ et $X_2$ maximisant la fonction économique. Le code vous est fourni, il ne vous reste plus qu'à implémenter le calcul et l'affichage de la fonction économique tenant compte de cette valeur réelle !

Explications

Ce code consiste simplement à déplacer les instructions affichant le polygone ainsi que la droite de la fonction économique dans une fonction, et à ajouter un slider permettant à l'utilisateur de choisir une valeur réelle, qui est passée à la fonction en question à chaque modification du slider.

Une fois que vous avez localisé globalement la position de la solution, vous pouvez modifier l'intervalle x pour zoomer sur cette partie et obtenir un résultat plus précis.

from ipywidgets import interact
import ipywidgets as widgets

# cette fonction est appelée à chaque changement d'état du slider.
def trace(n):
    # on reprend le code precedent
    xmax = 60
    x = np.arange(0, xmax) # on genere l'intervalle des x

    y1 = 100-x                                        #SOLUTION
    y2 = x                                            #SOLUTION
    y3 = 4*x                                          #SOLUTION
                                                      #SOLUTION
    ymax = max(np.amax(y1), np.amax(y2), np.amax(y3)) #SOLUTION
                                                      #SOLUTION
    plt.plot(x, y1, color='red')                      #SOLUTION
    plt.plot(x, y2, color='red')                      #SOLUTION
    plt.plot(x, y3, color='red')                      #SOLUTION
    plt.fill_between(x, y1, ymax, color='lightgrey')  #SOLUTION
    plt.fill_between(x, y2, color='lightgrey')        #SOLUTION
    plt.fill_between(x, y3, ymax, color='lightgrey')  #SOLUTION

    # on dessine entre 0 et les extrema
    plt.ylim(0, ymax)
    plt.xlim(0, xmax)
    pente = \
            -7/5      #SOLUTION
    yeco = \
           (x*pente)-n #SOLUTION
    plt.plot(x, yeco)
    
    plt.gcf().set_size_inches(18, 10) # cette fois on affiche en grand
    plt.show()
    
# le slider qui permet de faire varier la valeur passee en parametre de la fonction trace
interact(trace, n=widgets.FloatSlider(min=-200.0, max=0.0, step=1, value=-150.0));

interactive(children=(FloatSlider(value=-150.0, description='n', max=0.0, min=-200.0, step=1.0), Output()), _d…

Résolution

Utilisez le slider pour déplacer la droite sur le polygone. Que peut-on en déduire sur les valeurs de $X_1$ et $X_2$ ? Que vaut la fonction économique ?

Solution

Les élèves devront peut-être ajuster les paramètres d'affichage (les extrema du graphique) et les steps du slider pour obtenir la réponse la plus précise possible.

$X_1=50$
$X_2=50$
valeur optimale : 3000€

Imaginez qu'on ajoute un nombre quelconque de contraintes, donc de droites. Quelle propriété vérifie la solution optimale d'un problème linéaire en nombres réels sur cette partie ? Et que pouvez-vous en déduire sur la solution optimale d'un tel problème ?

Solution

Quel que soit le nombre de droites qu'on ajoute, l'espace des solutions admissible reste un polygone convexe. On voit bien que si on déplace la droite de la fonction économique jusqu'à la limite de cet espace, on trouvera forcément une solution optimale sur l'un des sommets du polygone.

Ce résultat s'étend à un nombre quelconque de dimensions, du fait de la structure linéaire du programme. On parle alors d'un polytope convexe (puisque borné par des hyperplans).

Théorème : Soit $f$ une fonction économique linéaire définie sur un polytope convexe borné. Alors la fonction $f$ atteint sa valeur maximale en au moins un des sommets du polytope convexe.

Ce résultat est très important, c'est lui qui permet de résoudre les problèmes linéaires de manière efficace.

2. Programme linéaire en nombre réels à $n$ dimensions

Arcelor-Mittal, géant mondial de l'acier, a ouvert une nouvelle usine en Grèce. Il fabrique 9 alliages à base de 3 matières premières : le plomb, le zinc et l'étain. Voici un tableau avec les pourcentages de chaque matière dans l'alliage et son coût par unité :

Alliage	1	2	3	4	5	6	7	8	9
% plomb	20	50	30	30	30	60	40	10	10
% zinc	30	40	20	40	30	30	50	30	10
% étain	50	10	50	30	40	10	10	60	80
Coût unit.	7.3	6.9	7.3	7.5	7.6	6.0	5.8	4.3	4.1

Le département de R&D recherche un « nouvel alliage » à prix minimum qui serait un mélange des 9 autres avec les contraintes d'avoir au maximum 30% de plomb, 30% de zinc et 40% d'étain.

2.1 Modélisation du problème

Commençons par modéliser le problème sous la forme d'un programme linéaire. Appliquez la méthode de l'exercice précédent.

Solution

Variables de décision
$X_n$ : Proportion de l'alliage N° $n$ présent dans le nouveau mélange

Contraintes du programme

• Contrainte 1 : maximum de 30% de plomb
  $0,2X_1 + 0,5X_2 + 0,3X_3 + 0,3X_4 + 0,3X_5 + 0,6X_6 + 0,4X_7 + 0,1X_8 + 0,1X_9 \leq 0,3$
  
  

• Contrainte 2 : maximum un 30% de zinc
  $0,3X_1 + 0,4X_2 + 0,2X_3 + 0,4X_4 + 0,3X_5 + 0,3X_6 + 0,5X_7 + 0,3X_8 + 0,1X_9 \leq 0,3$
  
  

• Contrainte 3 : maximum un 40% d’étain
  $0,5X_1 + 0,1X_2 + 0,5X_3 + 0,3X_4 + 0,4X_5 + 0,1X_6 + 0,1X_7 + 0,6X_8 + 0,8X_9 \leq 0,4$
  
  

• Contrainte 4 : la somme de proportions doit être égale à 1
  $\displaystyle\sum_{i=1}^{n}{X_i} = 1$
  
  

• Contraintes de non-négativité
  $X_i \geq 0 ~\forall i\in[1, 9]$

Fonction économique
Minimiser $F(X_1,\dots, X_9 ) = 7,3X_1+ 6_,9X_2 + 7,3X_3 + 7.5X_4 + 7.6X_5 + 6X_6 + 5.8X_7 + 4.3X_8 + 4,1X_9$

Programme linéaire $\begin{array}{rl@{}l} \text{Minimiser} & 7,3X_1+ 6_,9X_2 + 7,3X_3 + 7.5X_4 + 7.6X_5 + 6X_6 + 5.8X_7 + 4.3X_8 + 4,1X_9 &\\ \text{s.c.}& 0,2X_1 + 0,5X_2 + 0,3X_3 + 0,3X_4 + 0,3X_5 + 0,6X_6 + 0,4X_7 + 0,1X_8 + 0,1X_9 & \leq 0,3\\ & 0,3X_1 + 0,4X_2 + 0,2X_3 + 0,4X_4 + 0,3X_5 + 0,3X_6 + 0,5X_7 + 0,3X_8 + 0,1X_9 & \leq 0,3\\ & 0,5X_1 + 0,1X_2 + 0,5X_3 + 0,3X_4 + 0,4X_5 + 0,1X_6 + 0,1X_7 + 0,6X_8 + 0,8X_9 & \leq 0,4\\ & \displaystyle\sum_{i=1}^{n}{X_i} = 1\\ & X_i \geq 0 ~\forall i\in[1, 9]\\ \end{array}$

2.2 Résolution

Nous sommes en présence d'un problème à 9 variables. La méthode de résolution graphique ne peut pas être appliquée ici, nous allons nous tourner vers le solveur PuLP qui calculera lui-même la solution, en exploitant la propriété qu'on a identifiée dans l'exercice précédent.

Explications Il existe de nombreux solveurs sur le marché. Si le solveur d'Excel est rapidement mis à mal, des produits orientés vers le calcul haute performance sont disponibles, tant commerciaux (Cplex, AIMMS…) qu'open source (GLPK, LP_Solve…), dont les performances et les fonctionnalités varient grandement.

PuLP est en fait un framework Python qui permet de s'interfacer avec plusieurs solveurs (il s'interface par défaut avec CoinMP, avec lequel il est livré par défaut). L'intérêt, en dehors de pouvoir utiliser des solveurs commerciaux très performants, c'est que son usage est très intuitif.

On va s'inspirer d'un exemple basique pour implémenter et résoudre le modèle :

Explications

Dans cette solution (et celle de la question suivante), on écrit manuellement les équations. Ça reste encore faisable sur des programmes de cette taille, mais pour des programmes plus volumineux, il vaudra mieux utiliser des matrices de poids pour les contraintes et la fonction objectif, et les combiner avec un dictionnaire de variables grâce à lpSum(). Cette solution est présentée dans le Workshop Performances et optimisation des métaheuristique.

from pulp import *
import numpy as np

# variables du problème
X1 = LpVariable("X_1", 0, None)
X2 = LpVariable(\
    "X_2", 0, None) #SOLUTION
X3 = LpVariable(\
    "X_3", 0, None) #SOLUTION
X4 = LpVariable(\
    "X_4", 0, None) #SOLUTION
X5 = LpVariable(\
    "X_5", 0, None) #SOLUTION
X6 = LpVariable(\
    "X_6", 0, None) #SOLUTION
X7 = LpVariable(\
    "X_7", 0, None) #SOLUTION
X8 = LpVariable(\
    "X_8", 0, None) #SOLUTION
X9 = LpVariable(\
    "X_9", 0, None) #SOLUTION


# probleme
prob = LpProblem("alliage",\
                 LpMinimize) #SOLUTION

# objectif
prob += \
          7.3*X1 + 6.9*X2 + 7.3*X3 + 7.5*X4 + 7.6*X5 + 6.0*X6 + 5.8*X7 + 4.3*X8 + 4.1*X9 #SOLUTION

# contraintes
prob += \
          0.2*X1 + 0.5*X2 + 0.3*X3 + 0.3*X4 + 0.3*X5 + 0.6*X6 + 0.4*X7 + 0.1*X8 + 0.1*X9 <= 0.3 #SOLUTION
prob += \
          0.3*X1 + 0.4*X2 + 0.2*X3 + 0.4*X4 + 0.3*X5 + 0.3*X6 + 0.5*X7 + 0.3*X8 + 0.1*X9 <= 0.3 #SOLUTION
prob += \
          0.5*X1 + 0.1*X2 + 0.5*X3 + 0.3*X4 + 0.4*X5 + 0.1*X6 + 0.1*X7 + 0.6*X8 + 0.8*X9 <= 0.4 #SOLUTION
prob += \
          X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 == 1 #SOLUTION

prob.solve()
print(LpStatus[prob.status])
print("Min=", value(prob.objective))

# variables resultat
for v in prob.variables():
    print("%s=%.2f"%(v.name,v.varValue), end=', ')

Optimal
Min= 4.98
X_1=0.00, X_2=0.00, X_3=0.00, X_4=0.00, X_5=0.00, X_6=0.40, X_7=0.00, X_8=0.60, X_9=0.00, 

Nous avons notre solution optimale ! On peut même pousser le raffinement jusqu'à utiliser $\mathrm{\LaTeX}$ pour améliorer le rendu :

from IPython.display import display, Latex

# variables resultat
result = ""
for v in prob.variables():
    valeur = round(v.varValue*100, 2)
    if valeur > 0:
        result += "$"+v.name+"="+str(valeur)+"\%$, "

display(Latex(result))

$X_6=40.0\%$, $X_8=60.0\%$,

Une remarque : PuLP accepte aussi bien des problèmes de minimisation que de maximisation, et des contraintes d'infériorité que de supériorité. Mais tous les solveurs ne sont pas aussi souples, certains n'acceptent que programmes linéaires sous forme canonique pure, d'autres sous forme standard. Pourquoi ? À votre avis, comment PuLP peut-il faire pour s'interfacer avec ces solveurs solveurs ?

Solution

Un PL sous forme canonique s'écrit :

$ \begin{array}{rl} \text{Maximiser} & \displaystyle\sum_{i=1}^{n}a_nx_n\\ \text{s.c.} & Ax \leq b\\ & x \geq 0 \end{array} $

Un PL sous forme standard s'écrit :

$ \begin{array}{rl} \text{Maximiser} & \displaystyle\sum_{i=1}^{n}a_nx_n\\ \text{s.c.} & Ax = b\\ & x \geq 0 \end{array} $

L'algorithme du Simplexe ne sait résoudre que des programmes linéaires sous forme standard. Certains solveurs appliquent les transformations ci-dessous, d'autres non (Le package de programmation linéaire de SciPy ne le fait pas, par exemple). Dans ce cas, il faut le faire à la main, ou passer par un encapsuleur comme PuLP qui le fait lui-même.

Pour passer d'un PL quelconque à une forme canonique, si on doit passer d'une fonction économique à minimiser $F_{min}(X_1,\dots,X_n)$ à une fonction économique à maximiser $F_{max}(X_1,\dots,X_n)$, il suffit de prendre $F_{max}(X_1,\dots,X_n)=-F_{min}(X_1,\dots,X_n)$. Le même principe s'applique pour les contraintes de supériorité.

Pour passer d'une forme canonique à une forme standard, on peut transformer les inégalités en égalités par introduction de variables supplémentaires positives ou nulles appelées variables d'écart.

$a_1x_1 + \ldots + a_nx_n \leq b$

devient

$a_1x_1 + \ldots + a_nx_n + t = b$ avec $t\geq 0$

Par ailleurs, l'algorithme du Simplexe n'est pas le seul sachant résoudre un programme linéaire. PuLP peut en utiliser un autre. Le connaissez-vous ? Quels sont leurs principes généraux de fonctionnement, et leurs performances ?

Solution

• Simplexe : passe d'un sommet du polytope à l'autre jusqu'à converger vers l'optimal. Dans le pire des cas exponentiel, mais très souvent polynomial
• Méthode de points intérieurs : plus une classe d'algorithmes qu'un algorithme en particulier. Ces algorithmes passent par l'intérieur de l'espace des solutions admissibles, souvent en empruntant un chemin calculé en utilisant le principe du primal/dual. Les premiers algorithmes de ce type étaient très peu performants, mais les développements les plus récents les placent devant Simplexe

3. Programmation linéaire en nombres entiers

Un étudiant possède une bibliothèque musicale contenant 10 titres (certes, c'est peu, et il n'a même pas de compte sur un service de streaming, préférant gérer sa musique localement). Le logiciel qu'il utilise pour gérer sa bibliothèque lui permet de donner une note dans un intervalle de 1 à 5 étoiles à chaque titre, 5 étant la meilleure note. Notre étudiant a un très bon goût musical, et il écoute d'ailleurs sa musique dans un format Lossless, ce qui explique la taille des fichiers. Mais son smartphone est saturé par les milliers de photos, et les jeux inutiles, et sa capacité résiduelle pour stocker la musique est de 500Mo :

Titre	Note	Poids (Mo)
Led Zeppelin - Whole Lotta Love	☆☆	40
Magma - Retrovision (je suis revenu de l'univers)	☆☆☆☆☆	52
The Manhattan Project - Goodbye Pork Pie Hat (live)	☆	27
Jaco Pastorius - The Chicken (live)	☆☆	68
Steve Ray Vaughan - Little Wing	☆☆☆	80
Lynyrd SKynyrd - Free Bird (live)	☆☆☆☆☆	62
Genesis - Back In NYC	☆☆	63
The Who - I Can See For Miles	☆☆☆	53
Jean-Luc Ponty - Jig (live)	☆	56
David Bowie - Space Oddity	☆	59

Il souhaite remplir son smartphone de manière à avoir la meilleure liste de lecture possible à partir de sa bibliothèque, c'est à dire que la somme des notes des titres choisis doit être maximale. À quel problème problème connu cela vous fait-il penser ?

Solution

Il s'agit du problème du sac à dos (knap sack problem en anglais, souvent noté KSP). On considère un sac à dos ayant une capacité maximale, et des objets ayant chacun un encombrement et une valeur. Le but est de déterminer un sous-ensemble des objets disponibles qui puisse tenir dans le sac (c'est-à-dire dont la somme des encombrements ne dépasse pas la capacité du sac) et dont la valeur totale (la somme des valeurs des objets) soit maximale.

3.1 Modélisation du problème

Comment pouvons nous modéliser ce problème sous la forme d'un programme linéaire ?

Solution

On note :

• $w_k\in\mathbb{Z}^{+}$ l'encombrement de l'objet $k$
• $v_k\in\mathbb{Z}^{+}$ la valeur de l'objet $k$
• $x_k\in\{0,1\}$ la variable de décision de l'objet $k$, c'est-à-dire une variable booléenne indiquant si l'objet $k$ est placé dans le sac
• $c$ la capacité du sac

Le problème peut alors être formulé :
$ \begin{array}{rl} \text{Maximiser} & \displaystyle\sum_{i=1}^{n}v_kx_k\\ \text{s.c.} & \displaystyle\sum_{i=1}^{n}w_kx_k \leq c \end{array} $

Un algorithme naïf pour déterminer un résultat optimal serait une approche par force brute, qui teste toutes les solutions (admissibles ou non) et conserve la meilleure solution admissible. Dans le cas de notre étudiant, c'est une solution qui semble envisageable. Mais si l'étudiant supprime plusieurs Giga octets de jeux, et qu'il enrichit sa bibliothèque musicale, cette solution restera-t-elle viable ? Qu'est-ce qui va devenir bloquant ?

Solution

Cela revient à trouver toutes les combinaisons d'un ou plusieurs objets parmi les objets disponibles. Cette énumération considère toutes les solutions, y compris les non admissibles (celles qui ne tiennent pas dans le sac, c'est-à-dire qui ne respectent pas la contrainte $\displaystyle\sum_{i=1}^{n}w_kx_k \leq c$). Cette approche est valide, puisqu'on n'est pas capable de déterminer à l'avance les solutions qui seront admissibles.

Le paramètre qui va déterminer le nombre de solutions, c'est donc le nombre d'objets disponibles, ici les chansons de la bibliothèque. Le nombre de solutions en fonction du nombre d'objets disponibles est la combinaison de $i$ objets parmi $n$ pour tout $i\leq n$, avec $n$ le nombre d'objets disponibles, soit :
$\displaystyle \sum_{i=1}^{n}{n\choose i} = 2^{n}$

L'égalité avec $2^{n}$ (d'après l'une des propriétés du triangle de Pascal) s'obtient aussi en raisonnant sur le fait qu'une solution pour le problème à $n$ objets est un mot binaire de $n$ bits, chaque bit est associé à un objet, et indique si l'objet est présent dans le sac. Il y a $2^{n}$ valeurs possibles pour ce mot, donc autant de solutions.

Cette fonction est non polynomiale, mais là encore, cela ne suffit pas pour conclure que le problème est $\mathcal{NP}$-Complet, uniquement qu'on ne peut espérer pouvoir tester toutes les solutions pour un nombre important de chansons disponibles.

3.2 Résolution

Puisque nous avons modélisé le problème sous forme linéaire, pourquoi ne pas tester l'algorithme du Simplexe dessus ? Reprenons les données issues du Smartphone de notre étudiant. Attention aux bornes de nos variables !

# variables du problème
X1  = LpVariable("X_01", 0, 1)
X2  = LpVariable("X_02", 0, 1) #SOLUTION
X3  = LpVariable("X_03", 0, 1) #SOLUTION
X4  = LpVariable("X_04", 0, 1) #SOLUTION
X5  = LpVariable("X_05", 0, 1) #SOLUTION
X6  = LpVariable("X_06", 0, 1) #SOLUTION
X7  = LpVariable("X_07", 0, 1) #SOLUTION
X8  = LpVariable("X_08", 0, 1) #SOLUTION
X9  = LpVariable("X_09", 0, 1) #SOLUTION
X10 = LpVariable("X_10", 0, 1) #SOLUTION


# probleme
prob = LpProblem("musique", LpMaximize) #SOLUTION

# objectif
prob += 2*X1 + 5*X2 + X3 + 2*X4 + 3*X5 + 5*X6 + 2*X7 + 3*X8 + X9 + X10 #SOLUTION

# contrainte
prob += 40*X1 + 52*X2 + 27*X3 + 68*X4 + 80*X5 + 62*X6 + 63*X7 + 53*X8 + 56*X9 +59*X10 <=500 #SOLUTION

prob.solve()
print(LpStatus[prob.status])
print("Valeur totale = ", value(prob.objective))

poids_total = 0.0
poids = [40, 52, 27, 68, 80, 62, 63, 53, 56, 59]
for idx, v in enumerate(prob.variables()):
    print("X%d=%d" % (idx+1, v.varValue), end=', ')
    poids_total += v.varValue*poids[idx]
print("\npoids=%d" % poids_total)

Optimal
Valeur totale =  23.98214286
X1=1, X2=1, X3=1, X4=1, X5=1, X6=1, X7=1, X8=1, X9=0, X10=0, 
poids=500

Ça vous parait correct ? Pourtant, cela ne l'est pas ! Calculez à la main le poids de la solution, en reprenant les valeurs dans le tableau. Qu'obtenez-vous ? Pourquoi ?

Solution

$40+52+27+68+80+62+63+53=445$

On affiche les variables de décision comme des entiers (arrondis à l'inférieur par défaut), alors qu'en fait ce sont des réels. Si on les affiche en tant que tel, on voit bien le problème. La variable $X_9$ est arrondie à l'inférieur lors de l'affichage

Vérifions ça :

poids_total=0.0
for idx, v in enumerate(prob.variables()):
    print("X%d=%f" % (idx+1, v.varValue), end=', ') #SOLUTION
    poids_total += v.varValue*poids[idx]
print("\npoids=%d"%poids_total)

X1=1.000000, X2=1.000000, X3=1.000000, X4=1.000000, X5=1.000000, X6=1.000000, X7=1.000000, X8=1.000000, X9=0.982143, X10=0.000000, 
poids=500

Et là ça ne va plus du tout, La solution n'est pas admissible. On a vu qu'un arrondi inférieur produisait une solution qui, elle, est admissible. Mais avec une valeur de $X_9$ aussi proche de 1, ce n'est peut-être pas la solution optimale. Est-ce qu'un arrondi au plus proche marcherait mieux ?

poids_total=0.0
for idx, v in enumerate(prob.variables()):
    print("X%d=%f" % (idx+1, v.varValue), end=', ')
    poids_total += int(round(v.varValue))*poids[idx]
print("\npoids=%d" % poids_total)

X1=1.000000, X2=1.000000, X3=1.000000, X4=1.000000, X5=1.000000, X6=1.000000, X7=1.000000, X8=1.000000, X9=0.982143, X10=0.000000, 
poids=501

Qu'observe-t-on ?

Solution

La solution générée par arrondi supérieur génère une solution non admissible, elle dépasse la contrainte de contenance du sac.

Parmi toutes les solutions générées, nous en avons une solution qui est admissible. Est-ce une bonne méthode pour trouver une solution optimale ? Prenons cet exemple, dans lequel l'espace des solutions admissibles délimité par les deux contraintes est en blanc. Est-ce que, pour un problème de maximisation, arrondir une solution optimale réelle fonctionnerait ?

Solution

Non. Arrondir les variables de décision comme on a tenté de le faire ici ne permet pas de garantir l'optimalité, même si cela peut arriver. On n'est même pas sûr d'obtenir une solution admissible. Il suffit de considérer la fonction économique (en bleu) suivante :

Dans cet exemple, la variable en ordonnée peut être arrondie à l'inférieur, mais dans ce cas, la variable en abscisse ne peut être arrondie ni au supérieur, ni à l'inférieur, les deux solutions obtenues sont non admissibles.

C'est d'ailleurs un résultat connu. Un programme linéaire en nombre entier ou mixte entiers/réels ne peut être résolu en temps polynomial, à moins que $\mathcal{N}$=$\mathcal{NP}$ . La présence dans le programme linéaire des variables $X_k\in\{0,1\}$ rend donc ce problème $\mathcal{NP}$-Difficile.

Si l'on souhaite traiter ce problème du sac-à-dos par programmation linéaire, il faudra donc s'appuyer sur une approche beaucoup plus complexe, utilisant de la relaxation, des méthodes arborescentes etc.

D'ailleurs, PuLP permet cette approche, il suffit de changer le type des variables du programme (c'est LpContinuous par défaut, mais on peut imposer des variables de type LpInteger ou LpBinary). PuLP utilisera alors si c'est possible les implémentations disponibles dans la bibliothèque avec laquelle il est interfacé (voire des heuristiques si elles sont disponibles, et si le problème est trop complexe pour une méthode exacte). Si les étudiants découvrent cette possibilité en fouillant dans la doc de PuLP, il faut qu'ils comprennent que ces méthodes vont forcément avoir une complexité temporelle très importante (exponentielle dans le pire cas). Elles ne pourront présenter des performances moyennes un tant soit peu compétitives que sur des instances relativement petites, et au prix d'un travail important d'adaptation au problème considéré.

C'est finalement assez logique, puisque le problème de décision est $\mathcal{NP}$-Complet (on le sait depuis 1972). Le nombre exponentiel de solution, qu'on a démontré tout à l'heure, était-il un indice ?

Solution

Non. Ce n'est pas parce qu'un problème d'optimisation a un nombre exponentiel de solutions qu'il est $\mathcal{NP}$-Complet (plus précisément, que son problème de décision associé est $\mathcal{NP}$-Complet). De même, un problème que l'on a modélisé par un PLNE n'est pas nécessairement difficile. Il peut très bien exister un algorithme polynomial qui résout ce problème.

Par contre, si un problème est prouvé $\mathcal{NP}$-Complet, cela implique nécessairement qu'il admet un nombre exponentiel de solutions, et que sa modélisation par programmation linéaire nécessite forcément des variables en nombre entier.

4 Conclusion

Ce Workshop est maintenant fini, vous avez utilisé la programmation linéaire, vous en avez découvert les principes mathématiques de base, le fonctionnement général, et vous en avez testé les limites.

Explications

On n'a pourtant qu'effleuré le sujet, la programmation linéaire va bien plus loin, justement pour contourner ces limites. On pourrait d'ailleurs essayer de continuer à suivre cette piste pour le problème du sac à dos. Les méthodes d'optimisation de programme linéaire en nombre entier utilisent notamment les méthodes arborescentes pour traiter ce genre de problème. Mais il s'agit d'approches assez compliquées, et pour lesquelles les élèves n'ont pas le temps.

Et il y a la programmation mathématique en général, sujet encore plus vaste. Tous les problèmes ne sont pas linéaires… Le problème d'affectation quadratique, par exemple, est extrêmement connu et étudié. Il est particulièrement difficile, du fait de la présence d'inégalités quadratiques dans la fonction économique. Dans ce cas, une des approches qui donne de bons résultats est la programmation par contraintes, qui peut considérer des contraintes non linéaires. Là encore, c'est hors du périmètre de ce bloc.

Les étudiants peuvent s'informer sur ces méthodes, mais il est important qu'ils n'y passent pas trop de temps. Le but est de leur faire comprendre que, dans le cadre de cette séquence, ces méthodes, aussi efficaces soient-elles, ne sont pas envisageables, car trop compliquées à mettre en œuvre. Il faut surtout qu'ils prennent conscience de l'existence de ces approches, et de l'investissement que cela représente.

Et si vous êtes déçus de ne pas être venus à bout du problème du sac à dos, rassurez-vous, on ne s'avoue pas vaincus, dès la semaine prochaine, on s'y attaque à nouveau. Et cette fois-ci, avec beaucoup plus de succès !