Corentin Ducloux | 🚗 Optimisation Linéaire

1 Introduction au problème I

L’entreprise Car-Nivor souhaite se lancer dans la production de voitures et a déjà imaginé plusieurs modèles. Elle commence par analyser le marché avec les modèles Tiger (A) et Piranha (B).

Voici les informations qui ressortent de cette analyse :

Le modèle A lui apporte un bénéfice de 3 000 € par unité vendue, tandis que le modèle B lui apporte un bénéfice de 2 000 € par unité vendue.
La production d’une voiture de modèle A prend 2 heures tandis que celle d’une voiture de modèle B prend seulement 1 heure¹.
Pour des raisons d’absorption de marché, Car-Nivor ne peut pas vendre plus de 100 modèles A, 100 modèles B et 150 voitures (A+B) au total par mois.
Enfin, l’usine de production n’est disponible que 200 heures par mois.

2 Modélisation du problème I

Avant de se lancer tête baissée dans l’utilisation d’un modèle de Boosting, réfléchissons un peu et posons-nous quelques questions.

Quelles quantités peut-on faire varier ? $\Rightarrow$ Les quantités de modèles A et B.
Que cherche-t-on à optimiser ? $\Rightarrow$ Le profit de Car-Nivor.
Quelles sont les contraintes du problème ? $\Rightarrow$ Les contraintes de production et les contraintes de marché

Note

On voit en fait qu’on peut formuler cet énoncé sous la forme d’un problème d’optimisation linéaire sous contraintes !

Plus précisément, on considère les notations suivantes :

$x_{1}$ est la quantité vendue de modèle A,
$x_{2}$ est la quantité vendue de modèle B,
$M_{1}$ est le bénéfice en k€ par unité vendue de modèle A,
$M_{2}$ est le bénéfice en k€ par unité vendue de modèle B,
$z$ est le profit de Car-nivor, qu’on souhaite maximiser – on parle de fonction objectif,
$S$ correspond à l’ensemble des contraintes du problème.

$max_{s l c (S)} z = M_{1} x_{1} + M_{2} x_{2}$

$max_{s l c (S)} z = 3 x_{1} + 2 x_{2}$

$(S) : {\begin{matrix} x_{1} \geq 0; x_{2} \geq 0 \\ x_{1} \leq 100 \\ x_{2} \leq 100 \\ 2 x_{1} + x_{2} \leq 200 \\ x_{1} + x_{2} \leq 150 \end{matrix}$

Interprétation des contraintes

La première contrainte indique qu’on ne peut pas produire de quantité négative de voiture,
La seconde contrainte indique qu’on ne peut pas vendre plus de 100 modèles de type A,
La troisième contrainte indique qu’on ne peut pas vendre plus de 100 modèles de type B,
La quatrième contrainte porte sur le temps de disponibilité de l’usine (200h),
La dernière contrainte indique qu’il est possible de produire au maximum 150 voitures de Type A ou B, et pas plus.

3 Résolution du problème I

3.1 Import des librairies et configuration

import numpy as np
import scipy.optimize as so
import matplotlib.pyplot as plt
from rich import print
from rich.panel import Panel
import sympy as sym
import pandas as pd
from math import pi

plt.rcParams["figure.dpi"] = 600

3.2 Représentation graphique du problème

On peut facilement représenter dans un graph en $2 D$ l’ensemble des contraintes et proposer une solution graphique au problème puisqu’on dispose uniquement de 2 variables.

x = np.linspace(0, 150, 500)
y2 = 200 - 2 * x
y3 = 150 - x
z = (350 - 3 * x) / 2

Code du graphique matplotlib

plt.hlines(
    y=100, xmin=0, xmax=150, linewidth=2, color="r", label="$x_2 ≤ 100$", zorder=1
)
plt.vlines(
    x=100, ymin=0, ymax=150, linewidth=2, color="m", label="$x_1 ≤ 100$", zorder=1
)
plt.plot(x, y2, label="$2x_1 + x_2 \leq 200 $", zorder=1)
plt.plot(x, y3, label="$x_1 + x_2 \leq 150$", zorder=1)
plt.plot(x, z, label="$z = 350k$", zorder=1)

plt.scatter(50, 100, color="black", zorder=2)
plt.fill_between([0, 49.9], [100, 100], [0, 0], hatch="//", color="blue", alpha=0.2)
plt.fill_between(
    x,
    np.min([y2, y3], axis=0),
    where=y2 < 100,
    color="blue",
    alpha=0.2,
    hatch="//",
    interpolate=True,
    label="$admissible$",
)

plt.xlabel("$x_1$", fontsize=15)
plt.ylabel("$x_2$", fontsize=15)
plt.ylim([0, 125])
plt.xlim([0, 125])

plt.title("Ensemble des solutions admissibles")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)

3.3 Mise sous forme normale du problème

On a 4 équations dans le système des contraintes donc le rang est égal à 4.

Les variables positives sont appelées les variables de base :

$β = (\begin{matrix} X_{1} \\ X_{2} \\ X_{3} \\ X_{4} \end{matrix})$

L’idée est ensuite d’exprimer les variables de base en fonction des variables hors base ( $x_{1}$ et $x_{2}$ ) $\Rightarrow$ On réécrit donc le système comme suit :

${\begin{matrix} X_{1} = 100 - x_{1} \\ X_{2} = 100 - x_{2} \\ X_{3} = 200 - 2 x_{1} - x_{2} \\ X_{4} = 150 - x_{1} - x_{2} \end{matrix}$

3.4 Résolution avec `scipy`

Le submodule scipy.optimize fournit un ensemble de fonctions pour minimiser (ou maximiser) des fonctions objectif, possiblement sujettes à certaines contraintes.

Il faut d’abord écrire les matrices de travail avec numpy.

cost_1 = np.array([3, 2])
lhs_1 = np.array([[1, 0], [0, 1], [2, 1], [1, 1]])
rhs_1 = np.array([100, 100, 200, 150])

cost_1 correspond à la matrice de coût, i.e. les coefficients $M_{1}$ et $M_{2}$ ,
lhs_1² correspond au premier membre de la matrice des contraintes,
rhs_1³ correspond au second membre de la matrice des contraintes.

Prévisualisons-les avec sympy.

display(sym.Matrix(cost_1))
display(sym.Matrix(lhs_1))
display(sym.Matrix(rhs_1))

$[\begin{matrix} 3 \\ 2 \end{matrix}]$

$[\begin{matrix} 1 & 0 \\ 0 & 1 \\ 2 & 1 \\ 1 & 1 \end{matrix}]$

$[\begin{matrix} 100 \\ 100 \\ 200 \\ 150 \end{matrix}]$

Pour résoudre le problème, on peut utiliser la fonction linprog du submodule scipy.optimize, dont l’algorithme effectue par défaut une minimisation.

Comme nous souhaitons plutôt maximiser la fonction objectif, il faut inverser le signe de la fonction de coût.

result_1 = so.linprog(-cost_1, A_ub=lhs_1, b_ub=rhs_1, method="highs")

Profitons-en aussi pour créer les fonctions :

total_profit qui retourne le profit total optimal,
details qui affiche le profit ainsi que les quantités de voitures vendues.

Code des fonctions

def total_profit(result: so._optimize.OptimizeResult, cost: np.ndarray) -> np.float64:

    """Retourne le profit total en € lié à la vente de voitures."""

    return np.sum(result["x"] * cost * 1000)


def details(result: so._optimize.OptimizeResult, cost: np.ndarray, types: list) -> str:
    """
    Affiche un message contenant le profit optimal réalisé
    ainsi que le nombre et le type de voiture vendues.
    """
    profit = int(total_profit(result, cost))
    message = f"""
Le profit optimal réalisé par [i]Car-Nivor[/i] est de [u]{profit}[/u] € par mois.
    
=> Pour :
    
"""
    for voiture, type_voiture in zip(result["x"], types):
        message += f"- {voiture.astype(int)} voitures vendues de [red][b]Type {type_voiture}[/b][/red] \n"
    return Panel(message)

Il suffit maintenant d’appeler la fonction details pour obtenir les résultats.

details(result_1, cost_1, ["A", "B"])

╭─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│                                                                                                                 │
│ Le profit optimal réalisé par Car-Nivor est de 350000 € par mois.                                               │
│                                                                                                                 │
│ => Pour :                                                                                                       │
│                                                                                                                 │
│ - 50 voitures vendues de Type A                                                                                 │
│ - 100 voitures vendues de Type B                                                                                │
│                                                                                                                 │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯

3.5 Effondrement de la demande en modèle B

A cause de constructeurs chinois, Car-Nivor ne peut plus écouler qu’un maximum de 50 voitures par mois.

3.5.1 Graphique de la situation

Code du graphique matplotlib

plt.hlines(y=50, xmin=0, xmax=150, linewidth=2, color="r", label="$x_2 ≤ 50$", zorder=1)
plt.vlines(
    x=100, ymin=0, ymax=150, linewidth=2, color="m", label="$x_1 ≤ 100$", zorder=1
)

z = (325 - 3 * x) / 2

plt.plot(x, y2, label="$2x_1 + x_2 \leq 200 $", zorder=1)
plt.plot(x, y3, label="$x_1 + x_2 \leq 150$", zorder=1)
plt.plot(x, z, label="$z = 325k$", zorder=1)

plt.scatter(75, 50, color="black", zorder=2)

plt.fill_between(
    [0, 75],
    [49.9, 49.9],
    [0, 0],
    hatch="//",
    color="blue",
    alpha=0.2,
    interpolate=True,
    label="$admissible$",
)

plt.fill_between(
    x, y2, hatch="//", color="blue", where=y2 < 50, alpha=0.2, interpolate=True
)

plt.xlabel("$x_1$", fontsize=15)
plt.ylabel("$x_2$", fontsize=15)
plt.ylim([0, 125])
plt.xlim([0, 125])
plt.title("Ensemble des solutions admissibles")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)

Pour une résolution numérique , il suffit simplement de changer le second membre de la matrice des contraintes.

rhs_2 = np.array([100, 50, 200, 150])

result_2 = so.linprog(-cost_1, A_ub=lhs_1, b_ub=rhs_2, method="highs")

details(result_2, cost_1, ["A", "B"])

╭─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│                                                                                                                 │
│ Le profit optimal réalisé par Car-Nivor est de 325000 € par mois.                                               │
│                                                                                                                 │
│ => Pour :                                                                                                       │
│                                                                                                                 │
│ - 75 voitures vendues de Type A                                                                                 │
│ - 50 voitures vendues de Type B                                                                                 │
│                                                                                                                 │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯

On voit qu’un effondrement de la demande en voiture de Type B entraîne une surproduction de voitures de Type A, mais que cette surproduction ne permet pas d’obtenir le profit optimal avant l’effondrement de la demande de Type B.

3.6 Analyse détaillée du problème

On veut étudier graphiquement l’évolution du profit en fonction du nombre maximal de voitures vendues de type B.

Pour faire ceci, je vais utiliser une boucle for et faire varier le second membre de la contrainte pour $x_{2}$ .

list_x2 = [x for x in range(0, 101)]
cost_2 = np.array([3, 2])
lhs_2 = np.array([[1, 0], [0, 1], [2, 1], [1, 1]])

x1_vente = list()
x2_vente = list()
profit = list()

for x2_max in list_x2:
    rhs_2 = np.array([100, x2_max, 200, 150])
    max_bagnole = so.linprog(-cost_2, A_ub=lhs_2, b_ub=rhs_2, method="highs")
    x1_vente.append(max_bagnole["x"][0])
    x2_vente.append(max_bagnole["x"][1])
    profit.append(-max_bagnole["fun"] * 1000)

Code du graphique matplotlib

plt.plot(list_x2, profit)
plt.xlabel("Valeur max de $x_2$")
plt.ylabel("Profit en €")
plt.hlines(y=325000, xmin=0, xmax=50, linewidth=2, linestyles="--", color="m")
plt.vlines(x=50, ymin=300000, ymax=325000, linewidth=2, linestyles="--", color="m")

plt.hlines(y=350000, xmin=0, xmax=100, linewidth=2, linestyles="--", color="m")
plt.vlines(x=100, ymin=300000, ymax=350000, linewidth=2, linestyles="--", color="m")

plt.scatter(50, 325000, color="red", zorder=2, label="π2 = 325k")
plt.scatter(100, 350000, color="green", zorder=2, label="π1 = 350k")

plt.ylim([300000, 351000])
plt.xlim([0, 105])
plt.title(
    f"Profit optimal en fonction de la valeur de $x_2$ compris entre [{list_x2[0]}:{list_x2[-1]}] "
)

plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)

On remarque que le profit est croissant en fonction de la valeur maximale de $x_{2}$ (c’est à dire de voitures vendues de type B), ce dont on pouvait se douter étant donné les résultats précédents.

4 Introduction au problème II

Un manager propose de dessiner deux nouveaux modèles de luxe pour se distinguer des constructeurs chinois : le Lion (C) et le Shark (D). Les bénéfices de ces deux modèles sont de respectivement 5 000 et 10 000 €.

Néanmoins, la demande de ces modèles de luxe est assez limitée : 30 unités au total.

La demande mensuelle totale de voiture est toujours de 150.
La disponibilité de l’usine est toujours de 200 heures.
Le modèle C prend au total 5 heures de construction tandis que le modèle D en prend 8.

5 Modélisation du problème II

On considère les notations suivantes :

$x_{1}, x_{2}, x_{3}, x_{4}$ sont les quantités vendues de modèle A, B, C, D,
$M_{1}, M_{2}, M_{3}, M_{4}$ est le bénéfice en k€ par unité vendue de modèle A, B, C, D,
$z$ est le profit de Car-nivor, qu’on souhaite maximiser – on parle de fonction objectif,
$S$ correspond à l’ensemble des contraintes du problème.

$max_{s l c (S)} z = 3 x_{1} + 2 x_{2} + 5 x_{3} + 10 x_{4}$

$(S) : {\begin{matrix} x_{1} \geq 0; x_{2} \geq 0 \\ x_{1} \leq 100 \\ x_{2} \leq 100 \\ x_{3} \leq 30 \\ x_{4} \leq 30 \\ 2 x_{1} + x_{2} + 5 x_{3} + 8 x_{4} \leq 200 \\ x_{1} + x_{2} + x_{3} + x_{4} \leq 150 \end{matrix}$

6 Résolution du problème II

Attention

Nous ne pouvons plus représenter graphiquement la situation, la résolution va donc uniquement se poursuivre en utilisant scipy.

6.1 Résolution avec `scipy`

On procède aux mêmes étapes que dans la section précédente en écrivant les matrices de travail avec numpy.

cost_3 = np.array([3, 2, 5, 10])
lhs_3 = np.array(
    [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [2, 1, 5, 8], [1, 1, 1, 1]]
)
rhs_3 = np.array([100, 100, 30, 30, 200, 150])

result_3 = so.linprog(-cost_3, A_ub=lhs_3, b_ub=rhs_3, method="highs")

details(result_3, cost_3, ["A", "B", "C", "D"])

╭─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│                                                                                                                 │
│ Le profit optimal réalisé par Car-Nivor est de 350000 € par mois.                                               │
│                                                                                                                 │
│ => Pour :                                                                                                       │
│                                                                                                                 │
│ - 50 voitures vendues de Type A                                                                                 │
│ - 100 voitures vendues de Type B                                                                                │
│ - 0 voitures vendues de Type C                                                                                  │
│ - 0 voitures vendues de Type D                                                                                  │
│                                                                                                                 │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯

$\Rightarrow$ Dans cette configuration, Car-Nivor n’a pas intérêt à produire des voitures de luxe de Type C ou D.

6.2 La disponibilité de l’usine double

Dans le cadre du Plan de relance pour l’industrie, des aides conséquentes permettent à Car-Nivor de doubler le temps de disponibilité de l’usine pour atteindre 400 heures.

Pour résoudre le problème il suffit simplement de changer le second membre de la matrice des contraintes en remplaçant 200 par 400.

rhs_4 = np.array([100, 100, 30, 30, 400, 150])

result_4 = so.linprog(-cost_3, A_ub=lhs_3, b_ub=rhs_4, method="highs")

details(result_4, cost_3, ["A", "B", "C", "D"])

╭─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│                                                                                                                 │
│ Le profit optimal réalisé par Car-Nivor est de 580000 € par mois.                                               │
│                                                                                                                 │
│ => Pour :                                                                                                       │
│                                                                                                                 │
│ - 40 voitures vendues de Type A                                                                                 │
│ - 80 voitures vendues de Type B                                                                                 │
│ - 0 voitures vendues de Type C                                                                                  │
│ - 30 voitures vendues de Type D                                                                                 │
│                                                                                                                 │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯

Note

Un doublement de la disponibilité horaire entraîne la production de voitures de luxe de type D, ainsi qu’une augmentation substantielle du profit de Car-Nivor.

6.3 Impact du temps de disponibilité de l’usine sur la stratégie de l’entreprise.

Etudions l’impact du temps de disponibilité en utilisant des tranches de 100h de disponibilité mensuelle.

list_time = [100 * hours for hours in range(1, 11)]
cost_4 = np.array([3, 2, 5, 10])
lhs_4 = np.array(
    [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [2, 1, 5, 8], [1, 1, 1, 1]]
)

x1_vente = list()
x2_vente = list()
x3_vente = list()
x4_vente = list()
profit = list()

for hours in list_time:
    rhs_4 = np.array([100, 100, 30, 30, hours, 150])
    max_bagnole = so.linprog(-cost_4, A_ub=lhs_4, b_ub=rhs_4, method="highs")
    x1_vente.append(max_bagnole["x"][0])
    x2_vente.append(max_bagnole["x"][1])
    x3_vente.append(max_bagnole["x"][2])
    x4_vente.append(max_bagnole["x"][3])
    profit.append(-max_bagnole["fun"] * 1000)

Code du graphique matplotlib

plt.plot(list_time, profit)
plt.scatter(600, profit[5], color="green", zorder=2, label="π seuil = 720k")

plt.xlabel("Disponibilité horaire")
plt.ylabel("Profit en €")
plt.ylim([180000, 775000])
plt.xlim([90, 800])
plt.title(
    f"Profit optimal en fonction de la disponibilité horaire de l'usine \n compris entre [{list_time[0]}:{list_time[-1]}] heures"
)

plt.hlines(
    y=profit[5], xmin=0, xmax=600, linewidth=2, linestyles="--", color="m", zorder=1
)
plt.vlines(
    x=600, ymin=0, ymax=profit[5], linewidth=2, linestyles="--", color="m", zorder=1
)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)

Graphiquement, on remarque qu’il y n’y a plus intérêt à augmenter la disponibilité horaire à partir de 600 heures puisqu’on satisfait la contrainte de la demande de voiture $D_{x_{i}} = 150$ et que le profit n’augmente plus à partir de ce point.

On peut, en plus, représenter le nombre de voitures vendues selon leur type en faisant varier la disponibilité horaire.

Code du graphique matplotlib

plt.plot(list_time, x1_vente, label="Type A")
plt.plot(list_time, x2_vente, label="Type B")
plt.plot(list_time, x3_vente, label="Type C")
plt.plot(list_time, x4_vente, label="Type D")
plt.xlabel("Disponibilité horaire")
plt.ylabel("Nombre de voitures vendues")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)

Note

Si initialement les voitures de type B sont favorisées, on voit que Car-Nivor en produit de moins en moins pour focaliser sa production sur les type A, et dès qu’elle le peut, elle commence à produire des voitures de Type D. Pour finir, elle choisit de compléter sa production avec des voitures de Type C.

7 Introduction au problème III

L’usine de production fait remarquer à la direction que les flux de matières premières sont tendus et demande de les prendre en compte pour faire des modèles plus réalistes. Pour ce faire, elle crée le tableau suivant :

	Type A	Type B	Type C	Type D	Limite des stocks mensuels
Fer	2.0	1.0	4.0	2.0	125
Cuivre	0.2	0.2	0.5	1.0	55
Aluminium	1.0	1.0	2.0	2.0	300
Acier	2.0	2.0	1.0	2.0	100
Caoutchouc	1.0	1.0	1.5	0.8	60
Terres rares	0.0	0.0	2.0	2.0	50
Verre	1.0	1.0	1.0	1.2	100
Tissu	2.0	2.0	1.0	1.0	220
Cuir	0.0	0.0	1.0	1.0	30

Section en travaux

🐱‍🚀 Patience, ça arrive…

Retour au sommet

Notes de bas de page

Ces valeurs sont clairement fantaisistes et ne reflètent pas du tout le temps moyen d’assemblage réel d’une voiture, qui dépend de nombreux facteurs.↩︎
LHS = Left Hand Side / Premier membre d’une équation.↩︎
RHS = Right Hand Side / Second membre d’une équation.↩︎

1 Introduction au problème I

2 Modélisation du problème I

3 Résolution du problème I

3.1 Import des librairies et configuration

3.2 Représentation graphique du problème

3.3 Mise sous forme normale du problème

3.4 Résolution avec scipy

3.5 Effondrement de la demande en modèle B

3.5.1 Graphique de la situation

3.6 Analyse détaillée du problème

4 Introduction au problème II

5 Modélisation du problème II

6 Résolution du problème II

6.1 Résolution avec scipy

6.2 La disponibilité de l’usine double

6.3 Impact du temps de disponibilité de l’usine sur la stratégie de l’entreprise.

7 Introduction au problème III

Notes de bas de page

3.4 Résolution avec `scipy`

6.1 Résolution avec `scipy`