Algorithme du simplexe

L'algorithme du simplexe est un algorithme de résolution des problèmes d'optimisation linéaire. Il a été introduit par George Dantzig à partir de 1947. C'est probablement le premier algorithme permettant de minimiser une fonction sur un ensemble défini par des inégalités^[1]. De ce fait, il a beaucoup contribué au démarrage de l'optimisation numérique. L'algorithme du simplexe a longtemps été la méthode la plus utilisée pour résoudre les problèmes d'optimisation linéaire. Depuis les années 1985-90, il est concurrencé par les méthodes de points intérieurs, mais garde une place de choix dans certaines circonstances (en particulier si l'on a une idée des contraintes d'inégalité actives en la solution).

Le nom de l'algorithme est dérivé de la notion de simplexe et a été suggéré par Motzkin^[2]. En réalité, l'algorithme n'utilise pas de simplexes, mais certaines interprétations de l'ensemble admissible du problème renvoient au concept de simplexe.

Relevant de l'optimisation mathématique, l'algorithme du simplexe fait appel à des notions d'algèbre linéaire et de calcul différentiel.

Les problèmes d'optimisation linéaire que l'algorithme du simplexe résout consistent à minimiser une fonction linéaire de n {\displaystyle n} variables réelles,

x = ( x 1 , … , x n ) ∈ R n ↦ c ⊤ x := ∑ i = 1 n c i x i , {\displaystyle x=(x_{1},\ldots ,x_{n})\in \mathbb {R} ^{n}\mapsto c^{\top }x:=\sum _{i=1}^{n}\,c_{i}x_{i},}

où c = ( c 1 , … , c n ) ∈ R n {\displaystyle c=(c_{1},\ldots ,c_{n})\in \mathbb {R} ^{n}} , sur un ensemble défini au moyen de contraintes affines (ou linéaires) d'égalité et d'inégalité. L'ensemble admissible du problème est donc un polyèdre convexe, que nous noterons A P {\displaystyle {\mathcal {A}}_{P}} (A pour admissible, P pour primal). On écrit le problème sous la forme suivante

{ inf x c ⊤ x x ∈ A P {\displaystyle \left\{{\begin{array}{l}{\inf }_{x}\;c^{\top }x\\x\in {\mathcal {A}}_{P}\end{array}}\right.} .

L'ensemble admissible A P {\displaystyle {\mathcal {A}}_{P}} peut être défini de manières variées.

• La forme la plus couramment utilisée pour présenter et étudier les algorithmes, qui est celle supposée par l'algorithme du simplexe révisé, est la forme standard :
  A P = { x ∈ R n : A x = b ,   x ⩾ 0 } , {\displaystyle {\mathcal {A}}_{P}=\{x\in \mathbb {R} ^{n}:Ax=b,~x\geqslant 0\},}
  où A ∈ R m × n {\displaystyle A\in \mathbb {R} ^{m\times n}} est une matrice réelle de type m × n {\displaystyle m\times n} et b ∈ R m {\displaystyle b\in \mathbb {R} ^{m}} . L'écriture x ⩾ 0 {\displaystyle x\geqslant 0} signifie que toutes les variables x i {\displaystyle x_{i}} doivent être positives. Sous cette forme, A P {\displaystyle {\mathcal {A}}_{P}} apparaît comme l'intersection d'un sous-espace affine et de l'orthant positif.
• On rencontre aussi la forme canonique :
  A P = { x ∈ R n : A x ⩽ b } , {\displaystyle {\mathcal {A}}_{P}=\{x\in \mathbb {R} ^{n}:Ax\leqslant b\},}
  où A ∈ R m × n {\displaystyle A\in \mathbb {R} ^{m\times n}} est une matrice réelle de type m × n {\displaystyle m\times n} , b ∈ R m {\displaystyle b\in \mathbb {R} ^{m}} et l'inégalité A x ⩽ b {\displaystyle Ax\leqslant b} doit se comprendre composante par composante : ( A x − b ) i ⩽ 0 {\displaystyle (Ax-b)_{i}\leqslant 0} , pour i ∈ [ [ 1 , m ] ] ≡ { 1 , … , m } {\displaystyle i\in [\![1,m]\!]\equiv \{1,\ldots ,m\}} . Sous cette forme, A P {\displaystyle {\mathcal {A}}_{P}} apparaît comme l'intersection de m {\displaystyle m} demi-espaces de R n {\displaystyle \mathbb {R} ^{n}} . C'est souvent la forme utilisée pour illustrer le comportement de l'algorithme du simplexe.

Ces deux manières de représenter un polyèdre convexe sont équivalentes dans le sens où l'on peut passer d'une expression à l'autre. Ainsi, un polyèdre convexe représenté par la forme standard { x ∈ R n : A x = b ,   x ⩾ 0 } {\displaystyle \{x\in \mathbb {R} ^{n}:Ax=b,~x\geqslant 0\}} , peut s'exprimer sous la forme canonique par

{ x ∈ R n : A x ⩽ b ,   − A x ⩽ − b ,   − x ⩽ 0 } . {\displaystyle \{x\in \mathbb {R} ^{n}:Ax\leqslant b,~-Ax\leqslant -b,~-x\leqslant 0\}.}

De même, un polyèdre convexe représenté par la forme canonique { x ∈ R n : A x ⩽ b } {\displaystyle \{x\in \mathbb {R} ^{n}:Ax\leqslant b\}} , peut s'exprimer sous la forme standard par

{ ( u , v , s ) ∈ R n × R n × R m : A u − A v + s = b ,   u ⩾ 0 ,   v ⩾ 0 ,   s ⩾ 0 } , {\displaystyle \{(u,v,s)\in \mathbb {R} ^{n}\times \mathbb {R} ^{n}\times \mathbb {R} ^{m}:Au-Av+s=b,~u\geqslant 0,~v\geqslant 0,~s\geqslant 0\},}

où x {\displaystyle x} a été décomposé en u − v {\displaystyle u-v} avec u ⩾ 0 {\displaystyle u\geqslant 0} et v ⩾ 0 {\displaystyle v\geqslant 0} et les composantes de s {\displaystyle s} sont appelées des variables d'écart.

La version de l'algorithme du simplexe que nous présentons dans cette section est celle connue sous le nom d'algorithme du simplexe révisé. On suppose que le problème à résoudre est écrit sous la forme standard

( P L ) { inf x c ⊤ x A x = b x ⩾ 0 {\displaystyle (P_{L})\qquad \left\{{\begin{array}{l}{\inf }_{x}\;c^{\top }x\\Ax=b\\x\geqslant 0\end{array}}\right.} .

Il consiste donc à minimiser la fonction linéaire x ∈ R n ↦ c ⊤ x {\displaystyle x\in \mathbb {R} ^{n}\mapsto c^{\top }x} , aussi appelée fonction-coût, sur l'ensemble admissible

A P := { x ∈ R n : A x = b ,   x ⩾ 0 } . {\displaystyle {\mathcal {A}}_{P}:=\{x\in \mathbb {R} ^{n}:Ax=b,~x\geqslant 0\}.}

Cet ensemble est un polyèdre convexe, que l'on supposera non vide. On supposera également que

A {\displaystyle A} est surjective,

c'est-à-dire que ses lignes sont linéairement indépendantes.

On peut montrer que lorsque le problème ( P L ) {\displaystyle (P_{L})} a une solution, il a une solution sur un sommet du polyèdre convexe A P {\displaystyle {\mathcal {A}}_{P}} . On sait comment calculer tous ces sommets, qui sont en nombre fini, si bien que le problème de résoudre ( P L ) {\displaystyle (P_{L})} pourrait être de sélectionner le sommet qui donne à la fonction à minimiser sa plus petite valeur. Cependant, le nombre fini de sommets est en général très grand et dépend souvent exponentiellement des dimensions n {\displaystyle n} et m {\displaystyle m} du problème, si bien que cette approche ne pourrait être utilisée que pour résoudre des problèmes de petite dimension. L'algorithme du simplexe va rechercher une solution parmi les sommets de A P {\displaystyle {\mathcal {A}}_{P}} , mais en n'en calculant qu'une partie d'entre eux, en éliminant séquentiellement les sommets donnant à la fonction-coût une valeur supérieure à celle obtenue à l'itéré courant.

L'algorithme du simplexe est géométriquement très simple : chaque itération consiste à passer d'un sommet (face de dimension 0) du polyèdre convexe A P {\displaystyle {\mathcal {A}}_{P}} à un sommet adjacent en suivant une arête (face de dimension 1) particulière de ce polyèdre, de manière à faire décroître la fonction-coût. S'il n'y a pas de sommet le long de l'arête sélectionnée (parce que cette arête a une longueur infinie), le problème sera non borné (la valeur minimale de la fonction-coût sera − ∞ {\displaystyle -\infty } ). Il est bien de garder cette idée générale à l'esprit car la description algébrique d'une itération est relativement longue et doit prendre en compte quelques cas particuliers (problème non borné et pas nul) qui distraient.

Si x ∈ R + n {\displaystyle x\in \mathbb {R} _{+}^{n}} , on note

I + ( x ) := { i : x i > 0 } et I 0 ( x ) := { i : x i = 0 } . {\displaystyle I^{+}(x):=\{i:x_{i}>0\}\quad {\mbox{et}}\quad I^{0}(x):=\{i:x_{i}=0\}.}

Par propriété, x ∈ A P {\displaystyle x\in {\mathcal {A}}_{P}} est un sommet de A P {\displaystyle {\mathcal {A}}_{P}} si la sous-matrice A ∙ I + ( x ) {\displaystyle A_{\bullet I^{+}(x)}} , formée des colonnes de A {\displaystyle A} avec indices dans I + ( x ) {\displaystyle I^{+}(x)} , est injective. On a donc nécessairement | I + ( x ) | ⩽ m {\displaystyle |I^{+}(x)|\leqslant m} et on peut bien sûr avoir | I + ( x ) | < m {\displaystyle |I^{+}(x)|<m} . Un sommet x ∈ A P {\displaystyle x\in {\mathcal {A}}_{P}} est dit « dégénéré » si | I + ( x ) | < m {\displaystyle |I^{+}(x)|<m} et il est « non dégénéré » si | I + ( x ) | = m {\displaystyle |I^{+}(x)|=m} .

Une « base d'indices » est un ensemble B {\displaystyle B} de m {\displaystyle m} indices pris dans [ [ 1 , n ] ] {\displaystyle [\![1,n]\!]} tels que la sous-matrice A ∙ B {\displaystyle A_{\bullet B}} formée des m {\displaystyle m} colonnes correspondantes de A {\displaystyle A} soit inversible. Si x ∈ R n {\displaystyle x\in \mathbb {R} ^{n}} , les composantes x i {\displaystyle x_{i}} avec i ∈ B {\displaystyle i\in B} sont alors dites « basiques » et celles avec i ∉ B {\displaystyle i\notin B} sont dites « non basiques ». Une base d'indices B {\displaystyle B} est dite « associée » à un sommet x ∈ A P {\displaystyle x\in {\mathcal {A}}_{P}} si I + ( x ) ⊂ B {\displaystyle I^{+}(x)\subset B} .

La littérature scientifique en optimisation linéaire utilise un jargon ayant évolué depuis les contributions fondatrices. Certains termes peuvent porter à confusion. Il en va ainsi de « solution » pour désigner un point (par la suite,ce terme désignera une solution du problème d'optimisation linéaire), de « solution admissible » pour désigner un point admissible ou encore de « solution basique admissible » pour désigner un sommet de A P {\displaystyle {\mathcal {A}}_{P}} . Par suite, la terminologie sera conforme à celle utilisée en analyse convexe et en optimisation non linéaire.

Jargon de l'optimisation linéaire	Terminologie adoptée
base	base d'indices
solution	point
solution admissible	point admissible
solution basique	-
solution basique admissible	sommet
solution admissible optimale	solution
solution basique admissible optimale	solution-sommet

Une itération de l'algorithme démarre donc en un sommet x ^ {\displaystyle {\hat {x}}} de A P {\displaystyle {\mathcal {A}}_{P}} . Le calcul d'un tel sommet n'est pas une opération triviale, mais nous verrons plus loin comment on peut la réaliser. On note B {\displaystyle B} une base d'indices associée à ce sommet et N {\displaystyle N} le complémentaire de B {\displaystyle B} dans [ [ 1 , n ] ] {\displaystyle [\![1,n]\!]} . Au départ de l'itération, on a donc

x ^ B ⩾ 0 , x ^ N = 0 et A ∙ B x ^ B = b (donc   x ^ B = A ∙ B − 1 b ) , {\displaystyle {\hat {x}}_{B}\geqslant 0,\quad {\hat {x}}_{N}=0\quad {\mbox{et}}\quad A_{\bullet B}{\hat {x}}_{B}=b\quad {\mbox{(donc}}~{\hat {x}}_{B}=A_{\bullet B}^{-1}b),}

où, comme précédemment, A ∙ B {\displaystyle A_{\bullet B}} (resp. A ∙ N {\displaystyle A_{\bullet N}} ) désigne la sous-matrice de A {\displaystyle A} formée de ses colonnes avec indices dans B {\displaystyle B} (resp. N {\displaystyle N} ).

L'algorithme du simplexe génère en réalité une suite de bases d'indices plutôt qu'une suite de sommets. Il y a une distinction entre les deux notions lorsque le sommet est dégénéré, auquel cas il peut y avoir plusieurs bases d'indices correspondant à un même sommet. Si l'algorithme du simplexe visite un sommet dégénéré, il est possible qu'il ne change pas de sommet à l'itération suivante, mais il changera alors de base d'indices. Cependant, décrire l'algorithme en matière de bases d'indices fait perdre l'aspect géométrique de l'algorithme, qu'il nous semble précieux de conserver. Dès lors, nous considérerons que l'itéré de l'algorithme est un sommet, mais que certaines itérations font un déplacement nul.

Reconnaître l'optimalité

Soit x {\displaystyle x} un point de l'ensemble admissible A P {\displaystyle {\mathcal {A}}_{P}} , qui vérifie donc A x = b {\displaystyle Ax=b} . Comme B {\displaystyle B} est une base d'indices, on peut exprimer les composantes basiques x B {\displaystyle x_{B}} de x {\displaystyle x} en fonction de b {\displaystyle b} et de ses composantes non basiques x N {\displaystyle x_{N}}  :

x B = A ∙ B − 1 ( b − A ∙ N x N ) . {\displaystyle x_{B}=A_{\bullet B}^{-1}(b-A_{\bullet N}x_{N}).}

On peut également exprimer le coût c ⊤ x {\displaystyle c^{\top }x} en fonction de x N {\displaystyle x_{N}}  :

c ⊤ x = c B ⊤ A ∙ B − 1 ( b − A ∙ N x N ) + c N ⊤ x N . {\displaystyle c^{\top }x=c_{B}^{\top }A_{\bullet B}^{-1}(b-A_{\bullet N}x_{N})+c_{N}^{\top }x_{N}.}

Son gradient par rapport à x N {\displaystyle x_{N}} est appelé le coût réduit. Il s'écrit

r = c N − A N ⊤ A ∙ B − ⊤ c B . {\displaystyle r=c_{N}-A_{N}^{\top }A_{\bullet B}^{-\top }c_{B}.}

Dans l'algorithme du simplexe, ce coût réduit sert, d'une part, à détecter l'optimalité éventuelle de l'itéré courant x ^ {\displaystyle {\hat {x}}} et, d'autre part, à sélectionner une arête de A P {\displaystyle {\mathcal {A}}_{P}} le long de laquelle la fonction-coût décroît lorsque x ^ {\displaystyle {\hat {x}}} n'est pas optimal.

Si un sommet optimal est dégénéré, il peut y avoir un coût réduit r ⪈ 0 {\displaystyle r\ngeqslant 0} pour une base d'indices arbitraire associée à ce sommet, comme le montre l'exemple suivant :

{ min 2 x 1 + x 2 x 1 + x 2 = 0 x ⩾ 0. {\displaystyle \left\{{\begin{array}{l}\min \,2x_{1}+x_{2}\\x_{1}+x_{2}=0\\x\geqslant 0.\end{array}}\right.}

L'ensemble admissible est réduit au singleton { ( 0 , 0 ) } {\displaystyle \{(0,0)\}} . La solution du problème est donc forcément x ^ = ( 0 , 0 ) {\displaystyle {\hat {x}}=(0,0)} , qui est un sommet dégénéré. Si l'on prend pour base d'indices B = { 1 } {\displaystyle B=\{1\}} , le coût réduit r = − 1 {\displaystyle r=-1} est strictement négatif. Ceci signifie que l'on peut faire décroître le critère en augmentant x N = x 2 {\displaystyle x_{N}=x_{2}} tout en satisfaisant la contrainte d'égalité (c'est le sens du coût réduit). Mais ici on ne peut pas augmenter x 2 {\displaystyle x_{2}} sans sortir de l'ensemble admissible (ce ne serait pas le cas si le sommet était non dégénéré). À l'inverse, si l'on prend pour base d'indices B = { 2 } {\displaystyle B=\{2\}} , le coût réduit r = 1 {\displaystyle r=1} est positif, ce qui est annoncé par la proposition (il n'y a pas d'autre base d'indices).

La technique utilisée par l'algorithme du simplexe pour détecter l'optimalité éventuelle du sommet courant x ^ {\displaystyle {\hat {x}}} est la positivité du coût réduit, calculé en utilisant la base d'indices B {\displaystyle B} courante. Il s'agit d'un critère essentiellement primal (il ne fait pas intervenir de multiplicateur ou variable duale). Lorsque l'itéré courant est un sommet-solution non dégénéré, il n'y a qu'une seule base d'indices associée à ce sommet, si bien que le coût réduit est positif et l'algorithme s'interrompt. À l'inverse, lorsque l'itéré courant est un sommet-solution dégénéré, il se peut que la base d'indices courante ne permette pas d'avoir un coût réduit positif. Il est donc important que l'algorithme dispose d'un mécanisme lui permettant de changer de base d'indices si cela est nécessaire jusqu'à en trouver une permettant d'avoir un coût réduit positif (comme dans l'exemple ci-dessus). Un mécanisme permettant d'obtenir la convergence de l'algorithme du simplexe (c'est-à-dire d'éviter son cyclage) est appelé règle d'anti-cyclage ; les principales règles d'anti-cyclage seront vues ci-dessous.

Déplacement le long d'une arête

Si r {\displaystyle r} a une composante strictement négative, disons r j < 0 {\displaystyle r_{j}<0} , cela veut dire que l'on peut faire décroître le coût en augmentant la composante j {\displaystyle j} de x ^ N {\displaystyle {\hat {x}}_{N}} . On est donc tenté de chercher un nouveau point admissible en faisant un déplacement suivant une direction d {\displaystyle d} , c'est-à-dire

x ( α ) = x ^ + α d , {\displaystyle x(\alpha )={\hat {x}}+\alpha d,}

telle que la composante non basique de d {\displaystyle d} soit

d N = e N j . {\displaystyle d_{N}=e_{N}^{j}.}

On a noté e j {\displaystyle e^{j}} le j {\displaystyle j} -ième vecteur de base de R n {\displaystyle \mathbb {R} ^{n}} et e N j = ( e j ) N {\displaystyle e_{N}^{j}=(e^{j})_{N}} . Pour que ce déplacement d = ( d B , d N ) {\displaystyle d=(d_{B},d_{N})} soit acceptable, il faut d'abord que l'on ait A x ( α ) = b {\displaystyle Ax(\alpha )=b} , donc A d = 0 {\displaystyle Ad=0} , ce qui détermine sa composante basique :

d B = − A ∙ B − 1 A ∙ N e N j . {\displaystyle d_{B}=-A_{\bullet B}^{-1}A_{\bullet N}e_{N}^{j}.}

Sur le choix de l'indice j {\displaystyle j} . Remarquons que le coût décroît bien le long de d {\displaystyle d} puisque l'on a

c ⊤ d = r j < 0. {\displaystyle c^{\top }d=r_{j}<0.}

Si r {\displaystyle r} a plusieurs composantes strictement négatives, il semble donc raisonnable de choisir l'indice j {\displaystyle j} parmi ceux donnant la composante de r {\displaystyle r} la plus négative. C'est ce que l'on appelle la règle du coût réduit minimal. Cette règle ne garantit cependant pas l'efficacité globale de l'algorithme qui est principalement liée au nombre total d'itérations, c'est-à-dire au nombre de sommets visités (aspect global), ce qui ne peut se déduire d'un calcul de dérivée (aspect local). D'autres règles existent (comme celles introduites par les règles d'anti-cyclage décrites à la section Règles d'anti-cyclage) et les algorithmes du simplexe diffèrent en particulier par l'heuristique adoptée à cette étape.

Il est intéressant d'observer que le déplacement porté par la direction d {\displaystyle d} se fait le long d'une arête de A P {\displaystyle {\mathcal {A}}_{P}} .

Détection d'un problème non borné

Si d B ⩾ 0 {\displaystyle d_{B}\geqslant 0} , alors

∀ α > 0 : x ( α ) ⩾ 0 , {\displaystyle \forall \,\alpha >0:\qquad x(\alpha )\geqslant 0,}

et comme le coût c ⊤ x ( α ) = c ⊤ x ^ + α c ⊤ d {\displaystyle c^{\top }x(\alpha )=c^{\top }{\hat {x}}+\alpha c^{\top }d} décroît strictement le long de la direction de descente d {\displaystyle d} , le problème n'est pas borné.

Nouvelle base d'indices Si d B ⪈ 0 {\displaystyle d_{B}\ngeqslant 0} , on ne peut plus prendre un pas arbitrairement grand. Pour que l'on ait x ( α ) B ⩾ 0 {\displaystyle x(\alpha )_{B}\geqslant 0} , il faut que α ⩽ α ^ {\displaystyle \alpha \leqslant {\hat {\alpha }}} , où

α ^ = min { − x ^ i d i : i ∈ B ,   d i < 0 } . {\displaystyle {\hat {\alpha }}=\min \left\{-{\frac {{\hat {x}}_{i}}{d_{i}}}:i\in B,~d_{i}<0\right\}.}

Lorsque le sommet x ^ {\displaystyle {\hat {x}}} est dégénéré (il y a des x ^ i = 0 {\displaystyle {\hat {x}}_{i}=0} pour i ∈ B {\displaystyle i\in B} ), ce pas maximal α ^ {\displaystyle {\hat {\alpha }}} peut être nul (voir ci-après). Soit k {\displaystyle k} un indice donnant le min {\displaystyle \min } ci-dessus. Alors, x ^ k + α ^ d k = 0 {\displaystyle {\hat {x}}_{k}+{\hat {\alpha }}d_{k}=0} et on peut faire sortir l'indice k {\displaystyle k} de la base d'indices B {\displaystyle B} , et y faire entrer l'indice j {\displaystyle j} . La nouvelle base d'indices s'écrit

B + = ( B ∪ { j } ) ∖ { k } . {\displaystyle B_{+}=\left(B\cup \{j\}\right)\setminus \{k\}.}

L'opération de mise à jour de la base d'indices B {\displaystyle B} en B + {\displaystyle B_{+}} , qui consiste à lui adjoindre l'indice j {\displaystyle j} et à lui ôter l'indice k {\displaystyle k} , est parfois appelée pivotage et la règle déterminant le choix des indices j {\displaystyle j} et k {\displaystyle k} est alors appelée règle de pivotage.

Progrès ou stagnation

Deux situations peuvent maintenant se présenter.

• Si α ^ > 0 {\displaystyle {\hat {\alpha }}>0} , le coût décroît strictement et on peut passer à l'itération suivante avec x ^ + := x ^ + α ^ d {\displaystyle {\hat {x}}_{+}:={\hat {x}}+{\hat {\alpha }}d} comme nouveau sommet et B + {\displaystyle B_{+}} comme nouvelle base d'indices.
• Si α ^ = 0 {\displaystyle {\hat {\alpha }}=0} (ceci ne peut se produire que si le sommet x ^ {\displaystyle {\hat {x}}} est dégénéré), il y a un changement de base d'indices sans changer de sommet (le pas α ^ {\displaystyle {\hat {\alpha }}} est nul). Si l'on ne prend pas quelques précautions, l'algorithme peut cycler (par exemple en faisant entrer k {\displaystyle k} dans la base et en faisant sortir j {\displaystyle j} à l'itération suivante). On a mis au point des règles d'anti-cyclage pour faire face à cette situation. Certaines d'entre elles sont présentées dans la section suivante.

Nous énonçons ci-dessous quelques règles d'anti-cyclage et renvoyons le lecteur aux articles qui les introduisent pour une démonstration de leur propriété d'anti-cyclage. Ces articles sont souvent difficiles à comprendre, si l'on n'est pas familier avec le jargon développé par les spécialistes de l'algorithme du simplexe, en particulier avec la description de l'algorithme sous forme de tableau.

La règle de Bland (en)^[3] consiste à faire entrer dans la base B + {\displaystyle B_{+}} le plus petit indice j ∈ N {\displaystyle j\in N} tel que le coût réduit r j < 0 {\displaystyle r_{j}<0} (voir ci-dessus) et à en faire sortir le plus petit indice k ∈ a r g m i n ⁡ { − x ^ i / d i : i ∈ B {\displaystyle k\in \operatorname {arg\,min} \{-{\hat {x}}_{i}/d_{i}:i\in B} , d i < 0 } {\displaystyle d_{i}<0\}} (voir ci-dessus).

On peut résumer l'algorithme décrit ci-dessus comme suit.

On a le résultat de convergence finie suivant.

Pour utiliser l'algorithme du simplexe, il faut disposer d'un itéré initial qui est un sommet de l'ensemble admissible A P := { x ∈ R n : A x = b ,   x ⩾ 0 } {\displaystyle {\mathcal {A}}_{P}:=\{x\in \mathbb {R} ^{n}:Ax=b,~x\geqslant 0\}} . Nous présentons dans cette section plusieurs manières de faire face à cette exigence.

Comme son nom l'indique, la technique des deux phases décompose la résolution d'un problème d'optimisation linéaire en deux étapes. La phase/étape I consiste à résoudre un problème d'optimisation linéaire auxiliaire, dont on connait un sommet, par l'algorithme du simplexe. La solution de ce problème auxiliaire fournit un sommet du problème ( P L ) {\displaystyle (P_{L})} (si A P ≠ ∅ {\displaystyle {\mathcal {A}}_{P}\neq \varnothing } ) ou indique que ce problème n'a pas de point admissible. Dans la phase II, on résout le problème ( P L ) {\displaystyle (P_{L})} par l'algorithme du simplexe, à partir du sommet obtenu dans la première phase.

La phase I consiste à résoudre le problème d'optimisation linéaire auxiliaire suivant :

{ min ∑ i = 1 m z i A x + D z = b x ⩾ 0 , z ⩾ 0 , {\displaystyle \left\{{\begin{array}{l}\textstyle \min \,\sum _{i=1}^{m}z_{i}\\Ax+Dz=b\\x\geqslant 0,\quad z\geqslant 0,\end{array}}\right.}

où D {\displaystyle D} est la matrice diagonale définie par D i i = 1 {\displaystyle D_{ii}=1} si b i ⩾ 0 {\displaystyle b_{i}\geqslant 0} et D i i = − 1 {\displaystyle D_{ii}=-1} sinon. On peut utiliser pour cela l'algorithme du simplexe, démarrant en ( 0 , D b ) {\displaystyle (0,Db)} , qui est un sommet de son ensemble admissible

Cette technique combine les phases I et II de manière à ne devoir résoudre qu'un seul problème d'optimisation linéaire, à savoir le problème

( P L , M ) { min c ⊤ x + M ∑ i = 1 m z i A x + D z = b x ⩾ 0 , z ⩾ 0 , {\displaystyle (P_{L,M})\qquad \left\{{\begin{array}{l}\textstyle \min \,c^{\top }x+M\,\sum _{i=1}^{m}z_{i}\\Ax+Dz=b\\x\geqslant 0,\quad z\geqslant 0,\end{array}}\right.}

où D {\displaystyle D} est pris comme dans la technique des deux phases et M > 0 {\displaystyle M>0} est une constante choisie « suffisamment grande ». On ne connait pas a priori la valeur qu'il faut donner à M {\displaystyle M} pour que le problème ( P L , M ) {\displaystyle (P_{L,M})} soit équivalent au problème ( P L ) {\displaystyle (P_{L})} , c'est-à-dire pour qu'en une solution ( x ^ , z ^ ) {\displaystyle ({\hat {x}},{\hat {z}})} du problème ( P L , M ) {\displaystyle (P_{L,M})} on ait z ^ = 0 {\displaystyle {\hat {z}}=0} . Le raisonnement suivant montre pourquoi ils le seront si M {\displaystyle M} est « suffisamment grand ».

On peut en effet voir le problème ( P L , M ) {\displaystyle (P_{L,M})} comme la pénalisation exacte, au moyen de la norme ℓ 1 {\displaystyle \ell _{1}} , de la contrainte z = 0 {\displaystyle z=0} dans

{ min c ⊤ x A x + D z = b z = 0 x ⩾ 0 , z ⩾ 0. {\displaystyle \left\{{\begin{array}{l}\min \,c^{\top }x\\Ax+Dz=b\\z=0\\x\geqslant 0,\quad z\geqslant 0.\end{array}}\right.}

Ce dernier problème est équivalent à ( P L ) {\displaystyle (P_{L})} . Dès lors, selon la théorie de la pénalisation exacte, si M {\displaystyle M} est supérieur à la norme ℓ ∞ {\displaystyle \ell _{\infty }} (norme duale de la norme ℓ 1 {\displaystyle \ell _{1}} pour le produit scalaire euclidien) de tout multiplicateur optimal associé à la contrainte z = 0 {\displaystyle z=0} dans ce dernier problème, et si A P ≠ ∅ {\displaystyle {\mathcal {A}}_{P}\neq \varnothing } , alors toute solution ( x ^ , z ^ ) {\displaystyle ({\hat {x}},{\hat {z}})} du problème ( P L , M ) {\displaystyle (P_{L,M})} sera telle que z ^ = 0 {\displaystyle {\hat {z}}=0} et donc x ^ {\displaystyle {\hat {x}}} sera solution de ( P L ) {\displaystyle (P_{L})} .

En termes géométriques, l'ensemble des inégalités linéaires définit un polytope dans l'espace à n {\displaystyle n} dimensions (polygone en 2 dimensions et polyèdre en 3 dimensions) et il s'agit de trouver le sommet optimal pour la fonction de coût donnée. En effet, la fonction que l'on cherche à minimiser étant linéaire sur le polytope, elle y est en particulier concave. Or une fonction concave et minorée sur un polytope admet un minimum en un des sommets du polytope. La recherche d'un point de minimum peut donc se restreindre aux sommets du polytope (qui peuvent être très nombreux néanmoins).

L'idée de l'algorithme consiste à partir d'un sommet quelconque du polytope et, à chaque itération, d'aller à un sommet adjacent s'il est possible d'en trouver un meilleur pour la fonction objectif. S'il n'y en a pas, l'algorithme s'arrête en concluant que le sommet courant est optimal. En général, il y a plusieurs sommets adjacents au sommet courant qui sont meilleurs pour l'objectif. Il faut en sélectionner un seul, la règle de sélection est appelée règle de pivotage.

Le polytope est donné par f contraintes linéaires : une matrice A de taille (f, n) et un vecteur b de taille f. On cherche à trouver

min x ∈ R n , A x ≤ b φ ( x ) {\displaystyle \min _{x\in \mathbb {R} ^{n},\,Ax\leq b}\;\varphi (x)}

avec φ {\displaystyle \varphi } une forme linéaire. Un sommet du polytope se définit comme un sous-ensemble I ⊂ {1, …, f} tel que

• Card I = n
• la matrice carrée extraite A_I est inversible
• le point x_I = A_I⁻¹b_I vérifie A x_I ≤ b

La colonne j de A_I⁻¹ est le déplacement qui ajoute 1 à la I_j-ème forme linéaire de A, en gardant les formes I_k, k ≠ j, constantes. C'est un déplacement sur une arête du polytope. Une version de l'algorithme peut alors s'écrire comme suit :

• Trouver un sommet I ⊂ {1, …, f}
• n arêtes arrivent à I. Calculer la matrice A_I⁻¹. Pour k=1 à n, regarder le signe de φ.Col_k[A_I⁻¹].
  • Si ce signe est négatif ou nul pour tout k, Terminer. φ a son minimum en ce sommet I, sous les contraintes A x ≤ b.
  • Si ce signe est strictement positif pour un k, soit u = Col_k[A_I⁻¹] et aller au point suivant.
• Chercher le sommet le plus proche sur la demi-droite x_I + R^*₊u. Pour tout h ∈ {1, …, f} \ I, quand A_hu ≠ 0, calculer λ h = b h − A h x I A h u ∈ R {\displaystyle \lambda _{h}={\frac {b_{h}-A_{h}x_{I}}{A_{h}u}}\in \mathbb {R} } .
  • S'il y a un h avec λ_h > 0, prendre le plus petit de ces λ_h > 0. Assigner {h} ∪ I \ {k} à I et aller au point précédent.
  • Sinon terminer. Cette arête du polytope est infinie ; en la suivant on trouve min x ∈ R n , A x ≤ b φ ( x ) = − ∞ {\displaystyle \min _{x\in \mathbb {R} ^{n},\,Ax\leq b}\;\varphi (x)=-\infty } .

Il a été montré pour les principales règles de pivotage employées que l'algorithme du simplexe pouvait prendre un temps de calcul exponentiel. En particulier, on ne sait pas s'il existe une règle de pivotage qui assurerait que l'algorithme se termine après un nombre polynomial d'étapes.

On peut montrer que le nombre d'itérations de l'algorithme est majoré par : N − 2 ν − 1 {\displaystyle {\frac {N-2}{\nu -1}}} où ν {\displaystyle \nu } est le plus petit nombre d'arêtes reliées à un même sommet du polytope parcouru par le simplexe et N {\displaystyle N} est le nombre de sommets. On remarquera que ν {\displaystyle \nu } est minoré par la dimension de l'espace dans lequel vit le polytope.

Néanmoins, l'algorithme du simplexe est très efficace en pratique et il est implémenté dans tous les solveurs d'optimisation linéaire.

Une analyse un peu différente, l'analyse lisse permet d'expliquer l'efficacité du simplexe en pratique, en calculant des performances sur des entrées légèrement perturbées (Spielman et Teng 2004).

On trouve d'autres algorithmes de résolution de problèmes d'optimisation linéaire : l'algorithme en croix (en), la méthode du gradient projeté, la méthode du lagrangien augmenté, la méthode de l'ellipsoïde, les méthodes affines, les méthodes de points intérieurs, etc.

L'algorithme du simplexe est le plus souvent attribué à George Dantzig qui l'a découvert en 1947. Des techniques similaires avaient été découvertes par d'autres mathématiciens précédemment, mais sans parvenir à attirer l'attention^[4].

Cependant, Dantzig propose Fourier et son travail précurseur, en 1823-1830, sur l’algorithmique du simplexe. « De manière assez intéressante, en dépit de sa vaste applicabilité aux problèmes de tous les jours, la programmation linéaire resta inconnue avant 1947. Fourier était conscient de son potentiel dès 1823. » G.B. Dantzig, Linear Programming, Theory and Extensions ^[5]

• Méthode de Nelder-Mead
• Optimisation linéaire

• (en) Simplex Algorithm. Illustration détaillée de l'exécution de l'algorithme (version « tableau »).
• [PDF] Exemple de l'algorithme du Simplexe. L’algorithme du simplexe appliqué de manière très didactique à un exemple (version « tableau »).
• PHPSimplex: Outil en ligne pour résoudre les problèmes d'optimisation linéaire développé par Daniel Izquierdo et Juan José Ruiz, Université de Málaga (UMA, Espagne).
• J. Ch. Gilbert, Éléments d'Optimisation Différentiable — Théorie et Algorithmes, syllabus de cours à l'ENSTA ParisTech, Paris.

• (en) G.B. Dantzig (1990). Origins of the simplex method. In G. Nash, éditeur, A History of Scientific Computing, ACM Press Hist. Ser., pages 141–151. ACM Press, Reading, MA, États-Unis.
• (en) K.G. Murty (1983). Linear programming. John Wiley & Sons Inc., New York. (ISBN 0-471-09725-X). MR 720547.
• (en) M. Padberg (1999). Linear Optimization and Extensions, deuxième édition, Springer-Verlag.
• Daniel A. Spielman et Shang-Hua Teng, « Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time », Journal of the ACM, vol. 51, n^o 3,‎ 2004, p. 385–463 (lire en ligne)

• Portail de l'informatique théorique