Résolvez un programme linéaire avec le simplexe

La méthode graphique a résolu la brasserie en un coup d’œil. Mais si le brasseur produit huit bières sous quinze contraintes, plus question de dessiner : l’ensemble réalisable vit en dimension 8 et peut compter des milliers de sommets. Les énumérer tous serait beaucoup trop long.

L’algorithme du simplexe, inventé par George Dantzig en 1947, contourne le problème : il part d’un sommet et se déplace de sommet voisin en sommet voisin, en améliorant l’objectif à chaque pas, jusqu’à ce qu’aucun voisin ne fasse mieux. Dans ce chapitre, vous allez le dérouler à la main sur la brasserie, puis le programmer.

La forme standard et les variables d’écart

Le simplexe travaille avec des égalités, plus faciles à manipuler que des inégalités. On transforme chaque contrainte « inférieur ou égal » en égalité en ajoutant une variable d’écart (slack variable) positive, qui mesure ce qui reste de la ressource. Pour la brasserie :

$$\begin{aligned} \max \quad & z = 3x_1 + 2x_2 \\ \text{s.c.} \quad & x_1 + x_2 + e_1 = 4 \\ & x_1 + 3x_2 + e_2 = 9 \\ & x_1 + e_3 = 3 \\ & x_1, x_2, e_1, e_2, e_3 \geq 0 \end{aligned}$$

• $e_1$ est le nombre d’hectolitres de cuve inutilisés ;
• $e_2$ le nombre de sacs de malt restants ;
• $e_3$ le nombre de lots de houblon restants.

Une variable d’écart nulle signifie que la ressource est épuisée : on dit que la contrainte est saturée (ou active). Sous forme matricielle, le problème s’écrit $\max c^T x$ sous $Ax = b$, $x \geq 0$, avec ici $A$ une matrice à 3 lignes (les contraintes) et 5 colonnes (les variables).

Quelques règles pour se ramener à cette forme dans tous les cas :

Situation de départ	Transformation
contrainte $a^T x \leq b$	ajouter une variable d’écart : $a^T x + e = b$, $e \geq 0$
contrainte $a^T x \geq b$	retrancher une variable d’excédent : $a^T x - e = b$, $e \geq 0$ (voir le chapitre 3)
second membre $b$ négatif	multiplier la contrainte par $-1$ (le sens de l’inégalité s’inverse)
variable $x_j$ de signe quelconque	écrire $x_j = x_j^+ - x_j^-$ avec $x_j^+, x_j^- \geq 0$
minimiser $f$	maximiser $-f$ (ou adapter la règle d’entrée)

Solutions de base et sommets

Nous avons 3 équations et 5 inconnues. En fixant 2 variables à zéro, il reste un système de 3 équations à 3 inconnues, qui a (en général) une solution unique. C’est l’idée d’une solution de base :

• on choisit $m$ variables, dites en base (ici $m = 3$, le nombre de contraintes) ;
• on met toutes les autres, dites hors base, à zéro ;
• on résout le système pour obtenir les variables de base.

Si toutes les valeurs obtenues sont positives, la solution de base est réalisable. Et voici le lien avec la géométrie du chapitre 1 : les solutions de base réalisables sont exactement les sommets du polyèdre.

Exemple

Trois choix de base pour la brasserie :

• Base $\lbrace e_1, e_2, e_3 \rbrace$ : on pose $x_1 = x_2 = 0$, d’où $e_1 = 4$, $e_2 = 9$, $e_3 = 3$. C’est le sommet $(0, 0)$ : on ne brasse rien et toutes les ressources restent.
• Base $\lbrace x_1, e_1, e_2 \rbrace$ : on pose $x_2 = e_3 = 0$. La troisième équation donne $x_1 = 3$, puis $e_1 = 1$ et $e_2 = 6$. C’est le sommet $(3, 0)$.
• Base $\lbrace x_1, e_2, e_3 \rbrace$ : on pose $x_2 = e_1 = 0$. La première équation donne $x_1 = 4$, puis $e_2 = 5$ et $e_3 = 3 - 4 = -1$. Cette base n’est pas réalisable : le point $(4, 0)$ consommerait plus de houblon qu’il n’y en a.

Chaque intersection de deux droites est une solution de base. Celles qui respectent toutes les contraintes (disques pleins) sont exactement les sommets du polygone ; les autres, comme (4, 0) où e₃ = −1, ne sont pas réalisables.

Deux sommets sont voisins quand leurs bases ne diffèrent que d’une variable. Passer d’un sommet à un voisin, c’est donc faire entrer une variable dans la base et en faire sortir une autre. Toute la question est de choisir lesquelles.

Une itération du simplexe, pas à pas

On part de la base des écarts, toujours réalisable quand $b \geq 0$ : le sommet $(0, 0)$. On écrit chaque variable de base, et l’objectif, en fonction des variables hors base. Cette écriture s’appelle un dictionnaire :

$$\begin{aligned} e_1 &= 4 - x_1 - x_2 \\ e_2 &= 9 - x_1 - 3x_2 \\ e_3 &= 3 - x_1 \\ z &= 0 + 3x_1 + 2x_2 \end{aligned}$$

Choisir la variable entrante

Dans la ligne de $z$, les coefficients des variables hors base s’appellent les coûts réduits. Un coût réduit positif signifie qu’augmenter cette variable fait grimper le bénéfice. La règle usuelle (de Dantzig) fait entrer la variable de plus grand coût réduit positif : ici $x_1$, qui rapporte 3 par hectolitre.

Choisir la variable sortante : le test du rapport

Jusqu’où peut-on augmenter $x_1$, en laissant $x_2$ à zéro ? Tant qu’aucune variable de base ne devient négative :

• $e_1 = 4 - x_1 \geq 0$ impose $x_1 \leq 4$ ;
• $e_2 = 9 - x_1 \geq 0$ impose $x_1 \leq 9$ ;
• $e_3 = 3 - x_1 \geq 0$ impose $x_1 \leq 3$.

La contrainte la plus sévère est la dernière : $x_1$ monte jusqu’à 3, et c’est $e_3$ qui tombe à zéro et sort de la base. En général, on calcule pour chaque ligne le rapport $b_i / a_{ik}$ entre le second membre et le coefficient de la variable entrante, uniquement pour les coefficients $a_{ik}$ strictement positifs, et l’on retient le plus petit.

Le test du rapport : chaque variable de base impose une limite à la variable entrante. On s'arrête à la plus petite, ici x₁ = 3, et c'est e₃ qui tombe à zéro et sort de la base.

Pivoter

On exprime la variable entrante à partir de l’équation de la sortante, $x_1 = 3 - e_3$, et on la remplace partout :

$$\begin{aligned} x_1 &= 3 - e_3 \\ e_1 &= 1 - x_2 + e_3 \\ e_2 &= 6 - 3x_2 + e_3 \\ z &= 9 + 2x_2 - 3e_3 \end{aligned}$$

Nous sommes au sommet $(3, 0)$, avec un bénéfice de 9. Le coût réduit de $x_2$ est encore positif : on continue. $x_2$ entre ; les rapports valent $1/1 = 1$ pour $e_1$ et $6/3 = 2$ pour $e_2$ (la ligne de $x_1$ ne contient pas $x_2$ et ne limite rien). C’est $e_1$ qui sort. Avec $x_2 = 1 - e_1 + e_3$ :

$$\begin{aligned} x_1 &= 3 - e_3 \\ x_2 &= 1 - e_1 + e_3 \\ e_2 &= 3 + 3e_1 - 2e_3 \\ z &= 11 - 2e_1 - e_3 \end{aligned}$$

Tous les coûts réduits sont maintenant négatifs : augmenter $e_1$ ou $e_3$ ne pourrait que faire baisser $z$. On s’arrête. L’optimum est $x_1 = 3$, $x_2 = 1$, $z = 11$, avec $e_2 = 3$ sacs de malt inutilisés. C’est exactement le résultat de la méthode graphique, obtenu en visitant seulement trois sommets : $(0, 0)$, $(3, 0)$ puis $(3, 1)$.

Le simplexe longe les arêtes du polygone, de sommet voisin en sommet voisin, en améliorant z à chaque pivot : deux pivots suffisent pour atteindre l'optimum (3, 1).

À retenir

Une itération du simplexe (pour une maximisation) :

1. Entrée : la variable hors base de plus grand coût réduit positif. S’il n’y en a pas, la solution est optimale.
2. Sortie : la variable de base de plus petit rapport $b_i / a_{ik}$, calculé seulement pour $a_{ik} > 0$.
3. Pivot : exprimer la variable entrante avec l’équation de la sortante et la substituer partout.

Le tableau du simplexe

À la main, on range les mêmes calculs dans un tableau : une ligne par variable de base, une colonne par variable, le second membre à droite, et une dernière ligne pour les coûts réduits. Voici les trois tableaux de la brasserie (la case en bas à droite donne la valeur de $z$).

Tableau initial, base $\lbrace e_1, e_2, e_3 \rbrace$ :

Base	$x_1$	$x_2$	$e_1$	$e_2$	$e_3$	$b$
$e_1$	1	1	1	0	0	4
$e_2$	1	3	0	1	0	9
$e_3$	1	0	0	0	1	3
coûts réduits	3	2	0	0	0	$z = 0$

Après le premier pivot ($x_1$ entre, $e_3$ sort) :

Base	$x_1$	$x_2$	$e_1$	$e_2$	$e_3$	$b$
$e_1$	0	1	1	0	$-1$	1
$e_2$	0	3	0	1	$-1$	6
$x_1$	1	0	0	0	1	3
coûts réduits	0	2	0	0	$-3$	$z = 9$

Tableau final ($x_2$ entre, $e_1$ sort) :

Base	$x_1$	$x_2$	$e_1$	$e_2$	$e_3$	$b$
$x_2$	0	1	1	0	$-1$	1
$e_2$	0	0	$-3$	1	2	3
$x_1$	1	0	0	0	1	3
coûts réduits	0	0	$-2$	0	$-1$	$z = 11$

Le pivot (en gras) se fait par opérations sur les lignes : on divise la ligne du pivot par le pivot, puis on retranche à chaque autre ligne (y compris celle des coûts réduits) le bon multiple de la ligne du pivot pour annuler la colonne entrante. Gardez un œil sur les coûts réduits finaux des écarts, $-2$ et $-1$ : ils ont un sens économique précis, que vous découvrirez au chapitre 4.

Attention

Les conventions de signe varient d’un livre à l’autre. Certains écrivent la ligne de $z$ sous la forme $z - 3x_1 - 2x_2 = 0$ et cherchent alors le coefficient le plus négatif. Le raisonnement est le même ; fixez une convention et tenez-vous-y pendant tout le calcul.

Arrêt, non-bornitude et dégénérescence

Trois situations particulières se lisent directement dans le tableau.

• Problème non borné : si la colonne de la variable entrante n’a aucun coefficient strictement positif, rien ne limite sa croissance, et $z$ peut augmenter indéfiniment. On s’arrête en le signalant.
• Optimums multiples : si, à l’optimum, une variable hors base a un coût réduit nul, on peut la faire entrer sans changer $z$ : on obtient un autre sommet optimal.
• Dégénérescence : si une variable de base vaut 0, un pivot peut changer de base sans changer de sommet. Dans des cas rares, l’algorithme peut alors tourner en rond (on parle de cyclage). La règle de Bland (choisir, parmi les candidates, la variable de plus petit indice) l’évite à coup sûr.

Côté performance, le simplexe peut, dans le pire des cas, visiter un nombre exponentiel de sommets. En pratique, il en visite très peu et reste l’une des méthodes les plus utilisées au monde, à côté des méthodes de points intérieurs qui traversent le polyèdre au lieu d’en longer les arêtes.

Voici une implémentation compacte, pour le cas $\max c^T x$ sous $Ax \leq b$, $x \geq 0$, avec $b \geq 0$ :

import numpy as np

def simplexe(c, A, b):
    """Maximise c @ x sous A @ x <= b, x >= 0, avec b >= 0."""
    m, n = A.shape
    T = np.hstack([A, np.eye(m), b.reshape(-1, 1)]).astype(float)  # contraintes et écarts
    z = np.concatenate([c, np.zeros(m + 1)]).astype(float)          # coûts réduits, puis -valeur
    base = list(range(n, n + m))                                     # départ : les écarts en base
    while True:
        k = int(np.argmax(z[:-1]))              # variable entrante
        if z[k] <= 1e-12:
            break                               # plus de coût réduit positif : optimum
        col = T[:, k]
        if np.all(col <= 1e-12):
            raise ValueError("problème non borné")
        rapports = np.full(m, np.inf)
        rapports[col > 1e-12] = T[col > 1e-12, -1] / col[col > 1e-12]
        r = int(np.argmin(rapports))            # ligne de la variable sortante
        T[r] /= T[r, k]                         # pivot
        for i in range(m):
            if i != r:
                T[i] -= T[i, k] * T[r]
        z -= z[k] * T[r]
        base[r] = k
    x = np.zeros(n + m)
    x[base] = T[:, -1]
    return x[:n], -z[-1]

x, valeur = simplexe(np.array([3, 2]), np.array([[1, 1], [1, 3], [1, 0]]), np.array([4, 9, 3]))
print(x, valeur)   # [3. 1.] 11.0

En entretien

« Le simplexe se déplace de sommet en sommet du polyèdre réalisable, en faisant entrer dans la base une variable qui améliore l’objectif et sortir celle qui atteint zéro la première. Il s’arrête quand plus aucun coût réduit n’est favorable. »

En résumé

• La forme standard transforme les inégalités en égalités grâce aux variables d’écart, qui mesurent les ressources inutilisées.
• Une solution de base fixe à zéro les variables hors base et résout pour les $m$ variables de base ; les solutions de base réalisables sont les sommets du polyèdre.
• Le simplexe fait entrer la variable de meilleur coût réduit, fait sortir celle de plus petit rapport (sur les coefficients strictement positifs), puis pivote.
• Il s’arrête quand aucun coût réduit n’est favorable ; une colonne entrante sans coefficient positif signale un problème non borné.
• Un coût réduit nul à l’optimum révèle des optimums multiples ; la dégénérescence peut faire cycler l’algorithme, ce qu’évite la règle de Bland.

Exercices

Exercice · Identifier des solutions de base

facile

Pour la brasserie sous forme standard, calculez la solution de base associée à $\lbrace x_1, x_2, e_3 \rbrace$ et à $\lbrace x_2, e_1, e_3 \rbrace$. Sont-elles réalisables ? À quels sommets correspondent-elles ?

Voir la solution

Base $\lbrace x_1, x_2, e_3 \rbrace$ : on pose $e_1 = e_2 = 0$. Les deux premières équations deviennent $x_1 + x_2 = 4$ et $x_1 + 3x_2 = 9$, d’où $x_2 = 2{,}5$ et $x_1 = 1{,}5$. La troisième donne $e_3 = 3 - 1{,}5 = 1{,}5$. Tout est positif : c’est le sommet réalisable $(1{,}5 ; 2{,}5)$, où cuves et malt sont saturés.

Base $\lbrace x_2, e_1, e_3 \rbrace$ : on pose $x_1 = e_2 = 0$. La deuxième équation donne $x_2 = 3$, puis $e_1 = 4 - 3 = 1$ et $e_3 = 3$. Tout est positif : c’est le sommet $(0, 3)$, où seul le malt est saturé.

En code

Python

import numpy as np

# Forme standard : colonnes x1, x2, e1, e2, e3
noms = ["x1", "x2", "e1", "e2", "e3"]
A = np.array([[1, 1, 1, 0, 0],
              [1, 3, 0, 1, 0],
              [1, 0, 0, 0, 1]], dtype=float)
b = np.array([4, 9, 3], dtype=float)

for base in (["x1", "x2", "e3"], ["x2", "e1", "e3"]):
    idx = [noms.index(v) for v in base]
    valeurs = np.linalg.solve(A[:, idx], b)        # hors base à zéro
    print(base, valeurs, "réalisable" if np.all(valeurs >= 0) else "non réalisable")
# ['x1', 'x2', 'e3'] [1.5 2.5 1.5] réalisable
# ['x2', 'e1', 'e3'] [3. 1. 3.] réalisable

R

# Forme standard : colonnes x1, x2, e1, e2, e3
A <- rbind(c(1, 1, 1, 0, 0),
           c(1, 3, 0, 1, 0),
           c(1, 0, 0, 0, 1))
colnames(A) <- c("x1", "x2", "e1", "e2", "e3")
b <- c(4, 9, 3)

for (base in list(c("x1", "x2", "e3"), c("x2", "e1", "e3"))) {
  valeurs <- solve(A[, base], b)                  # hors base à zéro
  print(valeurs)
}
# x1 = 1.5, x2 = 2.5, e3 = 1.5 : sommet (1,5 ; 2,5)
# x2 = 3, e1 = 1, e3 = 3 : sommet (0, 3)

Exercice · Dérouler le simplexe

moyen

Résolvez par le simplexe le programme de l’exercice 1 du chapitre précédent : $\max z = x_1 + x_2$ sous $2x_1 + x_2 \leq 8$, $x_1 + 2x_2 \leq 7$, $x_1, x_2 \geq 0$. En cas d’égalité des coûts réduits, faites entrer la variable de plus petit indice.

Voir la solution

Dictionnaire initial : $e_1 = 8 - 2x_1 - x_2$, $e_2 = 7 - x_1 - 2x_2$, $z = x_1 + x_2$.

Itération 1 : $x_1$ entre. Rapports $8/2 = 4$ et $7/1 = 7$ : $e_1$ sort. On obtient $x_1 = 4 - \tfrac12 x_2 - \tfrac12 e_1$, puis $e_2 = 3 - \tfrac32 x_2 + \tfrac12 e_1$ et $z = 4 + \tfrac12 x_2 - \tfrac12 e_1$.

Itération 2 : $x_2$ entre (coût réduit $\tfrac12$). Rapports $4 / \tfrac12 = 8$ pour $x_1$ et $3 / \tfrac32 = 2$ pour $e_2$ : $e_2$ sort. Alors $x_2 = 2 + \tfrac13 e_1 - \tfrac23 e_2$, $x_1 = 3 - \tfrac23 e_1 + \tfrac13 e_2$ et $z = 5 - \tfrac13 e_1 - \tfrac13 e_2$.

Les coûts réduits sont négatifs : l’optimum est $(3, 2)$ avec $z = 5$, comme par la méthode graphique.

En code

Python

import numpy as np

A = np.array([[2, 1], [1, 2]], dtype=float)
b = np.array([8, 7], dtype=float)
c = np.array([1, 1], dtype=float)
m, n = A.shape
T = np.hstack([A, np.eye(m), b[:, None]])       # tableau : x1, x2, e1, e2 | b
z = np.concatenate([c, np.zeros(m + 1)])        # coûts réduits, puis -valeur
base = [n, n + 1]                               # départ : les écarts
while z[:-1].max() > 1e-12:
    k = int(np.argmax(z[:-1]))                  # plus grand coût réduit (plus petit indice)
    col = T[:, k]
    rapports = np.full(m, np.inf)
    rapports[col > 1e-12] = T[col > 1e-12, -1] / col[col > 1e-12]
    r = int(np.argmin(rapports))                # test du rapport
    T[r] /= T[r, k]
    for i in range(m):
        if i != r:
            T[i] -= T[i, k] * T[r]
    z -= z[k] * T[r]
    base[r] = k
    print(f"x{k + 1} entre, rapports {rapports}, z = {-z[-1]}, coûts réduits {z[:-1].round(3)}")
# x1 entre, rapports [4. 7.], z = 4.0, coûts réduits [ 0.   0.5 -0.5  0. ]
# x2 entre, rapports [8. 2.], z = 5.0, coûts réduits [ 0. 0. -0.333 -0.333]
x = np.zeros(n + m)
x[base] = T[:, -1]
print(x[:n])                                    # [3. 2.]

R

A <- rbind(c(2, 1), c(1, 2)); b <- c(8, 7); cc <- c(1, 1)
m <- nrow(A); n <- ncol(A)
T <- cbind(A, diag(m), b)            # tableau : x1, x2, e1, e2 | b
z <- c(cc, rep(0, m + 1))            # coûts réduits, puis -valeur
base <- n + 1:m                      # départ : les écarts
while (max(z[1:(n + m)]) > 1e-12) {
  k <- which.max(z[1:(n + m)])       # plus grand coût réduit (plus petit indice)
  rapports <- ifelse(T[, k] > 1e-12, T[, n + m + 1] / T[, k], Inf)
  r <- which.min(rapports)           # test du rapport
  T[r, ] <- T[r, ] / T[r, k]
  for (i in setdiff(1:m, r)) T[i, ] <- T[i, ] - T[i, k] * T[r, ]
  z <- z - z[k] * T[r, ]
  base[r] <- k
  cat(paste0("x", k, " entre, rapports ", rapports[1], " et ", rapports[2],
             ", z = ", -z[n + m + 1]), "\n")
}
# x1 entre, rapports 4 et 7, z = 4
# x2 entre, rapports 8 et 2, z = 5
round(unname(z[1:(n + m)]), 3)       # 0 0 -0.333 -0.333
T[order(base), n + m + 1]            # x1 = 3, x2 = 2

Exercice · Détecter un problème non borné

difficile

Appliquez le simplexe à $\max z = x_1 + x_2$ sous $-x_1 + x_2 \leq 1$, $x_1 - x_2 \leq 2$, $x_1, x_2 \geq 0$. Que se passe-t-il ? Donnez une famille de solutions réalisables dont la valeur tend vers l’infini.

Voir la solution

Dictionnaire initial : $e_1 = 1 + x_1 - x_2$, $e_2 = 2 - x_1 + x_2$, $z = x_1 + x_2$.

Itération 1 : $x_1$ entre. Dans la ligne de $e_1$, le coefficient de $x_1$ dans la contrainte est $-1$ (augmenter $x_1$ augmente $e_1$) : pas de limite. La ligne de $e_2$ impose $x_1 \leq 2$ : $e_2$ sort. On obtient $x_1 = 2 + x_2 - e_2$, $e_1 = 3 - e_2$ et $z = 2 + 2x_2 - e_2$.

Itération 2 : $x_2$ entre (coût réduit 2). Mais $x_2$ n’apparaît pas dans la ligne de $e_1$, et il apparaît avec un signe $+$ dans celle de $x_1$ : aucune variable de base ne diminue quand $x_2$ augmente. La colonne entrante n’a aucun coefficient strictement positif, le problème est non borné.

Concrètement, pour $t \geq 0$, le point $x_1 = 2 + t$, $x_2 = t$ est réalisable ($-x_1 + x_2 = -2 \leq 1$ et $x_1 - x_2 = 2 \leq 2$) et vaut $z = 2 + 2t$, qui tend vers l’infini.

En code

Python

import numpy as np
from scipy.optimize import linprog

A_ub, b_ub = [[-1, 1], [1, -1]], [1, 2]
res = linprog([-1, -1], A_ub=A_ub, b_ub=b_ub, bounds=[(0, None)] * 2)
print(res.status, res.message)       # 3 : le problème est non borné

# La famille x = (2 + t, t) reste réalisable et z = 2 + 2t grandit sans fin
for t in [0, 10, 1000]:
    x = np.array([2 + t, t])
    print(t, np.all(np.array(A_ub) @ x <= b_ub), x.sum())   # True, z = 2, 22, 2002

R

simplexe <- function(cc, A, b) {     # max cc'x sous A x <= b, x >= 0, b >= 0
  m <- nrow(A); n <- ncol(A)
  T <- cbind(A, diag(m), b); z <- c(cc, rep(0, m + 1)); base <- n + 1:m
  while (max(z[1:(n + m)]) > 1e-12) {
    k <- which.max(z[1:(n + m)])
    if (all(T[, k] <= 1e-12)) stop("colonne de x", k, " sans coefficient positif : non borné")
    r <- which.min(ifelse(T[, k] > 1e-12, T[, n + m + 1] / T[, k], Inf))
    T[r, ] <- T[r, ] / T[r, k]
    for (i in setdiff(1:m, r)) T[i, ] <- T[i, ] - T[i, k] * T[r, ]
    z <- z - z[k] * T[r, ]; base[r] <- k
    cat(paste0("x", k, " entre, z = ", -z[n + m + 1]), "\n")
  }
}
A <- rbind(c(-1, 1), c(1, -1)); b <- c(1, 2)
tryCatch(simplexe(c(1, 1), A, b), error = function(e) cat(conditionMessage(e), "\n"))
# x1 entre, z = 2
# colonne de x2 sans coefficient positif : non borné
t <- c(0, 10, 1000); x1 <- 2 + t; x2 <- t   # la famille (2 + t, t)
all(-x1 + x2 <= 1 & x1 - x2 <= 2)          # TRUE : toujours réalisable
x1 + x2                                    # 2 22 2002 : z tend vers l'infini