Chapitre 4 : Méthodes de résolution de systèmes d’équations et d’inéquations ☝

Avant de s’engouffrer dans les subtilités de l’algorithme du simplexe, il est utile de rappeler ce que peuvent apporter et comment utiliser certaines méthodes de résolution d’équations. Bien entendu, dans la pratique, tout se fait par recours à des logiciels spécialisés. Mais avoir fait tourner les algorithmes de la recherche opérationnelle à la main sur des exemples simples est une condition nécessaire pour comprendre les méthodes qu’ils permettent de développer. Il ne s’agit pourtant pas de proposer une substitution à un cours d’algèbre linéaire. Ce cours repose sur $\texttt{SageMath}$. Les méthodes pour automatiser les calculs seront donc proposées. Il n'en demeure pas moins vrai, qu'il est conseillé de commencer par réellement faire les calculs à la main.

Dans tout ce qui suit, on remarquera qu'en général ne seront utilisés que des nombres appartenant à $\mathbb{N}$ (des entiers) ou des nombres appartenant à $\mathbb{Q}$ (des rationnels). Pour un mathématicien ou un informaticien la raison va de soi, mais il n'est pas sûr qu'un économiste comprenne immédiatement le pourquoi d'un tel usage.

Bien entendu, il est plus simple de manipuler des nombres appartenant à ces deux ensembles dans les exemples réalisés à la main. Mais, il faut se rappeler que la programmation linéaire a un objectif pratique : résoudre des problèmes qui se posent dans la vie réelle. Or, il n'y a aucune raison, a priori, pour que ces problèmes soient formalisables avec des nombres qui appartiennent à ces deux ensembles. Depuis, que l'on sait que $\sqrt{2}$ n'appartient ni à l'un ni à l'autre de ces deux ensembles, on sait qu'en calcul rien n'est simple.

La raison est la suivante : un ordinateur a une capacité mémoire limitée. Or, les réels qui ne sont ni des entiers, ni des rationnels (ceux qui sont donc des irrationnels) nécessitent une capacité mémoire infinie pour être stockés autrement que formellement (on peut stocker un symbole comme $\pi$ mais pas son développement décimal). Il est donc nécessaire de les tronquer. Or, les tronquer revient à utiliser des rationnels. La raison en elle-même est élémentaire : $\mathbb{Q}$ est dense dans $\mathbb{R}$, ce qui est la manière des mathématiciens d'affirmer qu'entre deux réels, il existe toujours une infinité de rationnels. Cela peut sembler étonnant pour des êtres humains qui sont habitués à compter et dénombrer. Pour fixer ce point, imaginez que si vous faîtes une soiréechez vous et que votre appartement ne pouvait recevoir que 10 entiers, il pourrait recevoir aussi une infinité dénombrable de rationnels et une infinité non-dénombrable d'irrationnels. C'est la magie des nombres qui perturbe beaucoup de personnes peut réceptives à la beauté des mathématiques.

Il n'est pas question de donner une démonstration de cette propriété fondamentale de la topologie mathématique. En revanche, il est important de réaliser qu'entre deux rationnels il existe toujours un autre rationnel : on dit que $\mathbb{Q}$ est dense.

Pourquoi ? Rappelons que $x \in \mathbb{Q}$ si et seulement si il existe $p$ et $q$ appartenant tous deux à  $\mathbb{N}$ tels que $ x = p/q$. Considérons $x$ et $y$ dans $\mathbb{Q}$ et prenons en la moyenne arithmétique $z$ (i.e. : $z = (x+y)/2$).

Est-ce un rationnel? Comme $x$ est rationnel, il existe $x_1$ et $x_2$ dans $\mathbb{N}$ tel que $x = x_1/x_2$ et il en va de même pour $y$ (il existe $y_1$ et $y_2$ dans $\mathbb{R}$ tels que $y = y_1/y_2$). Donc :

Comme $x_1$,$x_2$, $y_1$ et $y_2$ sont des entiers, il en va de même pour $x_1 y_2$, $y_1 x_2$ et $x_2 y_2$ sont aussi des entiers. Donc $x_1 y_2 + y_1x_2$ est aussi un entier et $z$ est bien le rapport de deux entiers. Donc $z \in \mathbb{Q}$. Il reste à  démontrer que si $x < y$ alors $x < z < y$.

Mais, si $x < y$, il en va de même pour $x/2 < y/2$. Or : \[ z = \frac{x+y}{2}<\frac{y+y}{2} = y \hspace{.25cm}\text{et}\hspace{.25cm} z = \frac{x+y}{2}>\frac{x+x}{2} = x \]

Par conséquent, on a bien : $x < z < y$. En conséquence, $\mathbb{Q}$ est un ensemble dense. On peut donc toujours trouver un rationnel plus petit que deux rationnels donnés.

Le fait de travailler avec des rationnels est une nécessité quand on pense aux problèmes que peuvent poser les troncations liées à  l'usage pratique des réels et ce particulièrement si ces réels sont petits ou s'ils sont grands. Mais ce n'est pas la panacée. En effet, pour les grands nombres, on est obligé d'employer une notation en puissance de 10 qui est moins précise que la notation habituelle et qui peut concourir à  des erreurs de calcul, surtout quand ces erreurs se composent les unes avec les autres.

Avant de commencer à  aborder le sujet principal, il est important de comprendre un point essentiel pour la compréhension du développement des méthodes qui vont être présentées ci-après. Dans un cours d'initiation à  l'algèbre linéaire, on n'insiste pas trop sur les problèmes pratiques. Par exemple, si on doit s'intéresser au produit suivant : \[ \begin{bmatrix} \alpha & \beta & \gamma \end{bmatrix} \begin{bmatrix} a & b & c\\ d & e & f\\ g & h & i \end{bmatrix} \begin{bmatrix} A & B & C\\ D & E & F\\ G & H & I \end{bmatrix} \]

on n'insiste pas trop sur l'ordre des opérations. Calcule-t-on

puis

ou

puis

Après tout, peu importe. Dans les deux cas on aboutit au même vecteur. Du point de vue de l'algèbre linéaire théorique, c'est une vérité d'église puisque les produits matriciels sont associatifs. Mais, du point de vue du calcul en lui-même est-ce bien la même chose ?

Considérons la première méthode. Quand on effectue le produit matriciel d'un matrice $3\times 3$, on réalise 9 fois : 3 multiplications et deux additions, soit 27 multiplications et 18 additions. Puis, quand on réalise le produit du vecteur ligne par la matrice résultante, on ajoute 3 fois : 3 multiplications et 2 additions, soit 9 multiplications et 6 additions. Au bout du compte, on a réalisé 30 multiplications et 20 additions soit 50 opérations.

En ce qui concerne la seconde méthode, on commence par le produit d'un vecteur avec une matrice et on vient juste de voir que cela correspond à  9 multiplications et 6 additions dont la résultante est un vecteur $1 \times 3$ que l'on doit une seconde fois multiplier par une matrice $3 \times 3$, ce qui fait que l'on réalise le même type d'opérations deux fois. Au bout du compte, on a réalisé 18 multiplications et 12 additions soit 30 opérations. On a fait une économie de 20 opérations avec la seconde méthode soit 33.33\% des opérations. C'est énorme !

Maintenant, adoptons toujours le point de vue de l'économiste sur la résolution d'un système d'équations linéaire du type

On doit donc inverser une matrice $3\times 3$. On commence par calculer son déterminant:

Sachant que

pour calculer le déterminant d'une matrice $3\times 3$, on doit réaliser 3 fois le même groupe d'opérations :

1. \[ \begin{vmatrix} e & f\\ h & i \end{vmatrix} \hspace{.25cm}\text{2 multiplications, 1 soustraction} \]
2. $a (ef-hg)$ rajoute une multiplication.
3. une somme de puissances entières, soit 1 multiplication si la puissance est égale à  2, 2 multiplications si la puissance est égale à  3 et 3 multiplications si la puissance est égale à  4.
4. une élevation à la puissance
5. le produit du terme résultat (i.e. : $1$ ou $-1$ avec le résultat déjà  obtenu).

Au bout de ce premier compte, le calcul d'un déterminant d'une matrice $3\times3$ nécessite : 4 multiplications pour le premier bloc, cinq multiplications pour le second, six pour le troisième, 3 soustractions (1 par bloc) et deux additions soit 20 opérations. Maintenant, il faut calculer les matrices des cofacteurs de la matrice initiale soit

Chaque terme de la matrice des cofacteurs rajoute 2 multiplications et une addition et comme il y en a 9, ce calcul revient à  effectuer 18 multiplications et 9 additions. Maintenant, il faut la transposer ce qui revient à  9 opérations et pour finir diviser chaque terme par le déterminant, soit encore 9 divisions. Au bout du compte, on a effectué : 18 multiplications, 3 soustractions, 9 additions, 9 divisions et 1 transpositions qui correspond à  9 permutations.

Mais cette présentation cache la montagne: les multiplications et donc les divisions, prennent beaucoup plus de temps calcul qu'une addition (une soustraction) et une transposition. En réfléchissant bien, une multiplication de deux nombreux à  deux digits (le digit est un symbole, par exemple, un chiffre indo-arabique, utilisé pour l'écriture d'un nombre, en numération de position), par exemple $15 \times 16$, correspond à  7 opérations ($6 \times 15$ correspond à  2 opérations et $1\times 16$ aussi, l'addition de leur résultat correspond à  3 opérations).

Que se passe-t-il si on rajoute un digit à  chaque nombre ? Ce nouveau digit doit être multiplié par tous les digits qui composent l'autre nombre. Pour arriver à  multiplier plusieurs nombres composés de plus de chiffres, nous devons considérer que l'ajout d'un autre chiffre signifie que nous avons maintenant à  multiplier ce chiffre par tous les autres chiffres dans l'autre nombre. Si la multiplication de deux nombres à  deux digits donne $2^2=4$ opérations de multiplication, la multiplication de deux nombres de $n$ digits nécessitera $n^2$ opérations. Une fois les multiplications achevées, on aura deux nombres à  $2 n$ digits que l'on devra additionner, ce qui fait $4 n$ opérations. Finalement, la multiplication de deux nombres à  $n$ digits nécessite $n^2 + 4 n$ opérations (on remarque que pour $n = 2$, la formule donne 8 multiplications alors que dans l'exemple, on en a eu que $7$. C'est parce toutes les multiplications de deux nombres à  deux digits ne donnent pas un nombre à  4 digits. Par exemple, $32 \times 31 = 992$ alors que $32\times 32 = 1024$.). Comparé aux $n$ opérations que prend l'addition, la multiplication est extrêmement coûteuse. Bien entendu, les nombres stockés dans un ordinateur le sont en notation binaire et le coût de la multiplication n'est pas vraiment décrit par l'équation présentée ci-dessus. Mais le message est assez clair : il est nécessaire, quand on le peut, de diminuer ses occurrences. C'est pourquoi depuis deux siècles, mais surtout depuis l'arrivée des ordinateurs, les mathématiciens ont essayé de repenser les opérations de manière à  ce qu'elles soient le moins coûteuse en temps (surtout quand le programme qui doit être calculé en évalue un très grand nombre).

Les efforts des informaticiens ont été concentrés sur cette économie de temps. Par exemple, si on effectue la multiplication suivante $125467\times 256174=32141383258$ d'après la formule précédente, il sera nécessaire de réaliser $6^2 = 36$ multiplications. A. A. Karatsuba et Ofman^* A. A. Karatsuba, & Y. Ofman (1963). Multiplication of multidigit numbers on automata. Soviet Physics Doklady, 7, 595-596. ont trouvé un moyen plus économe. On note tout d'abord que ces deux nombres peuvent s'écrire sous la forme $(a \times 10^3) + b$ (i.e. : $(125\times 1000)+ 467$ pour le premier et $(256 \times 1000) + 174$ pour le second). Leur produit prend donc la forme :

A priori, cette multiplication nécessite 4 multiplications $a_1 a_2$, $a_1b_1$, $b_1 a_2$ et $b_1 b_2$, sachant que le découpage en puissance de 10 est arbitraire et que, de ce fait, les puissance de 10 peuvent être stockées à  l'avance ce qui ne demande aucun temps-machine. En réalité, il est possible de n'avoir recours qu'à  3 multiplications ($a_1 a_2$, $b_1 b_2$ et $(a_1-b_1)(a_2 -b_2)$), si on note que :

De la sorte, pour réaliser la multiplication précédente, on posera : $a_1 = 125$, $a_2 = 256$, $b_1 = 467$ et $b_2 = 174$ et on calculera :

ce qui permet d'écrire :

Tout cela n'est évidemment pas fait pour être réalisé à  la main. Mais de toute évidence, quand on cherche à  gagner du temps machine cela a une importance. Bien entendu, ce qui peut apparaître comme une chicanerie informatique s'est développé à  une époque où les ordinateurs étaient peu puissants et où le principal terminal d'entrée était constitué de cartes perforées. Il était très important de gagner du temps et d'éviter des manipulation Mais, il ne faut pas croire que c'est un passé totalement disparu. Depuis, les demandes de calculs ont été accrues exponentiellement et les calculs sont souvent d'une dimension sans commune mesure avec ce qu'ils étaient dans les années 1960.

D'ailleurs, très peu de temps après la publication de l'algorithme de Karatsuba et Ofman, A. L. Toom ^* A. L. Toom (1963). The complexity of a scheme of functional elements realizing the multiplication of integers. Doklady Akademii Nauk SSSR, 4(3). a encore amélioré la méthode.

Réaliser les multiplications suivantes par la méthode habituelle puis par le méthode de A. A. Karatsuba et Y. Ofman^* A. A. Karatsuba, & Y. Ofman (1963). Multiplication of multidigit numbers on automata. Soviet Physics Doklady, 7, 595-596. .

1. $23456\times 14524$;
2. $34567 \times 12706$;
3. $4987235 \times 2345612$.

Cette méthode dont le nom vérifie la loi d'éponymie de Stigler, puisque sa première référence est donnée au chapitre huit Tableaux Rectangulaires (Fangcheng) du livre Neufs Chapitres sur l'Art des Mathématiques, livres chinois qui aurait été écrit en 150, mais dont la première édition qui nous soit parvenue date de 179. En Europe, on doit la méthode à  Isaac Newton.

Soit à  résoudre un système du type $\boldsymbol{A}\boldsymbol{x}= \boldsymbol{b}$, c'est-à -dire :

$m$ et $n$ ne sont pas nécessairement égaux. La méthode de Gauss-Jordan permet en effectuant 3 types d'opérations élémentaires sur les lignes de la matrice de trouver une solution, de montrer qu'il n'y en a pas ou qu'il en existe une infinité. Il s'agit de transformer la matrice $\boldsymbol{A}$ en matrice triangulaire supérieure en appliquant aussi les mêmes transformations sur $\boldsymbol{b}$ et de résoudre par récurrence le système d'équations qu'elle commande. Les trois types d'opérations sont :

1. permuter des lignes ;
2. ajouter à  une ligne, un multiple scalaire d'une autre ligne ;
3. multiplier une ligne par un scalaire non nul.

L'idée est de réussir à  transformer la matrice de droite en une matrice triangulaire supérieure, c'est-à-dire une matrice telle que tous les éléments situés sous la principale diagonale deviennent nuls, toute opération appliquée sur $\boldsymbol{A}$ devant aussi être appliquée sur $\boldsymbol{b}$. Finalement, à l'arrivée, on doit avoir :

De la sorte, on a directement $x_n = \overline{b}_n$ et comme $\overline{a}_{(n-1)(m-1)}x_{n-1} +\overline{a}_{(n)(m-1)} x_{n} = \overline{b}_{(n-1)}$ et par substitution, on obtient :

Ainsi de suite, en réitérant la substitution arrière, on finit par trouver le vecteur des $x_i$ qui résout le systèmes d'équations. Comme ces trois opérations doivent être menées simultanément sur $\boldsymbol{A}$ et $\boldsymbol{b}$ et que les variables $\boldsymbol{x}$ ne sont pas affectées, on écrit un tableau augmenté $\begin{bmatrix} \boldsymbol{A} | \boldsymbol{b}\end{bmatrix}$ de la forme

(la barre de séparation entre $\boldsymbol{A}$ et $\boldsymbol{b}$ étant facultative) et c'est sur ce tableau qu'on applique les opérations. Considérons le système

l'on transforme en

Éliminons $a_{21} = 2$ par l'opération $L_2\ - 2\ L_1$, ce qui donne :

Puis, éliminons $a_{31}=1$ par l'opération $L_3 - L_1$ :

Pour finir, on peut éliminer $a_{32}=1$ par l'opération $L_3 - L_2$, ce qui donne

En opérant ainsi, on est passé d'un système $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ à  un système $\boldsymbol{U}\boldsymbol{x}=\boldsymbol{c}$. On a alors :

Ainsi, on a $x_3 = -3/4$. En substituant ce résultat dans $x_2 - 3 x_3 = - 5$, ce qui donne $x_2 = -5+9/4 = -11/4$ et $x_1+x_2 + 2 x_3 = 3$ et $x_1 = 12/4 +11/4+ 6/4 = 29/4$. Il est possible de vérifier directement que ce résultat est le bon en substituant $x_1$, $x_2$ et $x_3$ dans la seconde équation du système initial, c'est-à-dire :

Bien entendu, si le système est irréalisable, il y aura une incohérence dans les équations qui aboutira, par exemple, à  nous forcer à  résoudre une équation du type $x_3 = \alpha \not= 0$. Ce point sera abordé à  la section suivante.

Il est facile de demander à $\texttt{SageMath}$ de faire tourner la méthode de Gauss sur une matrice. Dans une discution du forum $\texttt{Ask SageMath}$ (pour avoir accès à ce forum et/ou poser une question, il suffit de passer par une barre de recherche), le code qui suit a été proposé. Comme, ce code est utilisé uniquement pour des raisons pédagogiques, on ne s'intéresse qu'à des matrices de rationnels.

$\texttt{Code}$

Les cellules qui se suivent avec un cli (pour i = 1,2,...) identique sont liées et doivent être calculées dans l'ordre. Pour éviter que le serveur ne se déconnecte, il est conseillé des les lancer assez rapidement les unes après les autres.

$\texttt{Réduction Gaussienne}$

La matrice qui est utilisée dans la cellule qui suit peut être modifiée.

Le calculs précédents ont une portée essentiellement pédagogique. $\texttt{Ask SageMath}$ offre une commande qui fait directement le travail.

Déterminer la valeur de chaque élément pour trouver le nombre qui remplace le point d'interrogation (inspiré de R. Galland^* R. W. Galland (1997). Le Grand livre des énigmes Nikola Tesla – 100 énigmes de physique, challenges mathématiques et casse-têtes logiques ! Marabout. ).

				15
				13
				?
				13
13	16	22	13

Il s'agit comme son nom l'indique d'un développement de la méthode de Gauss. Soit toujours à  résoudre un système du type $\boldsymbol{A}\boldsymbol{x}= \boldsymbol{b}$, c'est-à -dire :

$m$ et $n$ ne sont pas nécessairement égaux. La méthode de Gauss-Jordan permet en effectuant 3 types d'opérations élémentaires sur les lignes de la matrice de trouver une solution, de montrer qu'il n'y en a pas ou qu'il en existe une infinité. Contrairement à  la méthode de Gauss, on recherche à  diagonaliser $\boldsymbol{A}$. Les trois types d'opérations sont encore:

1. Permuter des lignes ;
2. Ajouter à  une ligne, un multiple scalaire d'une autre ligne ;
3. Multiplier une ligne par un scalaire non nul.

Considérons l'exemple précédent :

que l'on transforme en

La première étape consiste à  trouver un pivot, c'est-à-dire une entrée de la matrice de gauche telle que l'on désire que tous les éléments qui sont dans la même colonne puissent être rendus nuls par une opération de type (Cette étape sera réalisée autant de fois qu'il y a de colonnes. Il est évident que le nombre d'étapes nécessaire à  la complétion d'une élimination de Gauss-Jordan est fini.}. Cette étape doit être réalisée en gardant à  l'esprit les règles suivantes :

1. toujours commencer par la colonne qui possède le plus de zéros.
2. N'utiliser une colonne (une ligne) qu'une seule fois et si possible utiliser un pivot situé sur la principale diagonale.
3. Ne jamais pivoter sur un zéro ni sur un élément de droite du tableau.

En général, on choisit une entrée égale à  l'unité (ici, $a_{11}= 1$). On réalise l'opération suivante : $\ L_2 - 2 L_1$. On obtient donc :

Puis, on effectue de $L_3 - L_1$

On recherche maintenant le pivot dans la seconde colonne ($a_{22} = 1$. On opère $L_1 - L_2$) :

$L_3 - L_2$, i.e. :

On n'a plus de $1$ dans la troisième colonne, mais on peut en réintroduire un en effectuant l'opération $-L_3/8$

On peut alors réaliser l'opération : $L_1 - L_3$

Alors, on réalise l'opération $L_2 - 3 L_3$, soit :

Finalement, on a transformé le système d'équations

ce qui montre que la solution est $\boldsymbol{x}^\top = \begin{bmatrix} 29/4 & -11/4 & -3/4 \end{bmatrix} $.

Considérons maintenant le système suivant :

que l'on transforme en

Répétons les étapes de l'algorithme de Gauss-Jordan, mais cette fois en regroupant le plus grand nombre d'opérations possible. La première étape permet de prendre un pivot sur la première ligne ($-L_1$) soit :

Puis, en une unique étape, on réalise les opérations suivantes : $L_2 - L_1$, $L_3 - 2L_1$, $L_4 - 12\ L_1$ et $L_5 - 3\ L_1$, ce qui donne :

On trouve le pivot de la seconde colonne par $-L_2$, i.e. :

On élimine la seconde colonne par la suite d'opérations suivante : $L_1 + 7 L_2$, $L_3 - 3L_2$, $L_4-72 L_2$ et $L_5 - 2 L_2$, ce qui donne :

Pour obtenir le pivot de la troisième colonne, on commence par l'opération $L_4/88$, pour obtenir :

On élimine la troisième colonne par la suite d'opérations suivante : $L_1 - 156 L_3$, $L_2 + 21L_3$, $L_4-1630 L_3$ et $L_5 - 88 L_2$, ce qui donne

On s'aperçoit alors que le système est inconsistant puisque la ligne $L_4$ signifie que :

est rigoureusement impossible. Ce système n'a donc pas de solution. Considérons pour terminer cette suite d'exemples, le système suivant :

que l'on transforme en :

On sélectionne le pivot par la première opération ($L_1/2$), ce qui donne

puis on réalise $L_2 + 6 L_1$ pour obtenir

$L_2$ est parfaitement compatible ($0 = 0$) et laisse une seule équation :

équation qui peut être vérifiée par une infinité de couples $(x_1, x_2)$.

La même méthode peut servir aussi à  inverser une matrice (si tant est qu'elle soit inversible). Cette fois, on applique la méthode à  partir du tableau \[ \left[\begin{array}{cccc|cccc} a_{11} & a_{12} & \cdots & a_{1m}& 1 & 0 & \cdots & 0\\ a_{21} & a_{22} & \cdots & a_{2m}& 0 & 1 & \cdots& 0\\ \vdots & \vdots & \ddots &\vdots& \vdots& \vdots& \ddots\\ a_{n1} & a_{n2} & \cdots & a_{nn}& 0& 0& \cdots & 1 \end{array} \right] \]

Considérons la matrice \[ \begin{bmatrix} 2 & 0 & 2\\ 0 & 2 & 0\\ 3 & 0 & 1 \end{bmatrix} \]

On part du tableau

On construit le pivot de la première colonne avec $L_1/2$ et celui de la seconde avec $L_2/2$ ce qui donne

Si on applique $L_3 - 3 L_1$, on obtient

On applique $-L_3/2$

Pour finir, on soustrait $L_3$ à  $L_1$ (on calcule $L_1-L_3$) : \[ \begin{array}{cc} L_1 & \\ L_2 &\\ L_3 \end{array} \left[ \begin{array}{ccc|ccc} 1 & 0 & 0 & -1/4 & 0 & 1/2\\ 0 & 1 & 0 & 0 & 1/2 & 0 \\ 0 & 0 & 1 & 3/4 & 0 & -1/2 \end{array}\right] \]

Il est alors élémentaire de vérifier que la matrice

est bien la matrice inverse de la matrice initiale puisque :

Il existe une commande dans $\texttt{SageMath}$ qui permet d'obtenir la solution d'un système matriciel par l'élimination de Gauss-Jordan (on notera la manière alternative simple dont on définit une matrice : on définit une liste puis on donne le nombre de liste et de colonnes et, éventuellement, l'anneau dans lequel les calculs seront effectués).

On notera que l'on peut obtenir le même résultat avec la commande $\mathtt{M.echelon\_form('row\_reduction')}$.

La commande a donc permis de résoudre le système d'équations linéaire suivant : \[ \begin{bmatrix} 1&2&3&4&5 \\ 1/2&2/3&3&4&5 \\ 1&2/7&3/5&4&5 \\ 6&2&3/5&4&5/11 \\ 6/5&1/9&1&2/13&5 \\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4\\ x_5 \end{bmatrix}= \begin{bmatrix} 2/3\\ 5/3\\ 3/7\\ 4\\ 7 \end{bmatrix} \]

Pour obtenir la solution :

$\texttt{SageMath}$ offre une autre méthode, plus directe, pour obtenir la solution d'un système d'équations linéaires.

Les vrais flemmards peuvent aussi remplacer $\mathtt{A.solve\_right(b)}$ par $\mathtt{A\backslash b}$.

On notera qu'il existe d'autres façons de résoudre un système d'équations linéaires. Celà provient du fait que $\texttt{SageMath}$ est un compendium de logiciels qui ont leurs propres solveurs. Les personnes intéressées peuvent consulter la page Linear algebra.

On pourra se demander pourquoi présenter l'algorithme d'élimination de Gauss alors que l'algorithme de Gauss-Jordan semble une amélioration du premier. Simplement, parce que l'on peut démontrer que l'algorithme de Gauss est plus rapide que celui de Gauss-Jordan.

L'idée de la décomposition d'une matrice $\boldsymbol{A}$ comme le produit de deux matrice $ \boldsymbol{L}$ et $\boldsymbol{U}$ remonte à  \citet{Turing48}\sindex[aut]{Turing A.}. Elle a été anticipé sur un cas particulier par la méthode de Cholesky (voir \citet{CommandantBenoit24}\sindex[aut]{Commandant Benoit}). Elle consiste à  décomposer une matrice en deux sous matrices triangulaires, la première ($\boldsymbol{L}$ comme lower) étant triangulaire inférieure et la seconde ($\boldsymbol{U}$ comme upper) étant triangulaire supérieure. Par exemple, pour une matrice $4\times3$, on aura

Il n'est pas toujours possible de décomposer directement la matrice $\boldsymbol{A}$ sous la forme $\boldsymbol{L}\boldsymbol{U}$. C'est, par exemple, le cas de la matrice :

simplement du fait de la présence du zéro. Cependant, il est possible, par permutation, de transformer $\boldsymbol{A}$ de telle sorte que l'opération puisse être validée. Pour ce faire, on introduit une matrice de permutation $\boldsymbol{P}$ telle que :

où $\boldsymbol{P}$ est une matrice qui échange les lignes d'une autre matrice. Par exemple :

est une matrice de permutation^! Précisément, une matrice de permutation est composée de $0$ et de $1$ ayant un unique $1$ par ligne et par colonne..

Si on applique cette matrice à  une matrice $3\times 3$ quelconque, par exemple la matrice

on obtient

L'intérêt de la décomposition^!. On parle aussi de factorisation. $\boldsymbol{L} \boldsymbol{U}$ est encore une fois de diminuer la complexité liée à  l'inversion d'une matrice nécessaire pour résoudre un système d'équations. Supposons que $\boldsymbol{P} = \boldsymbol{I}$, alors, si le système :

est décomposé, on a :

on peut, dans un premier temps, le résoudre en introduisant le vecteur $\boldsymbol{y} =\boldsymbol{U}\boldsymbol{x}$, ce qui donne :

nouveau système qui est solvable de manière récursive vers le bas (on parle de descente) étant donné la structure de $\boldsymbol{L}$ ce qui est moins coûteux que d'inverser $\boldsymbol{L}$. Enfin, de $\boldsymbol{y}$, il est élémentaire de trouver $\boldsymbol{x}$ par des opérations du même type. On notera que cette seconde étape se fait par remontée, étant donnée la structure de $\boldsymbol{U}$. Objectivement, la décomposition $\boldsymbol{L} \boldsymbol{U}$ n'est autre qu'une forme de l'élimination gaussienne. Il demeure néanmoins à  exposer un algorithme qui permet de l'obtenir directement. Voici une méthode qui évite de réaliser un pivotage.

On écrit les matrices comme des matrices par block, i.e. :

où $a_{11}$ sont des scalaires et les $\boldsymbol{0}_k$ ($k= m \text{ ou }n$ des vecteurs de $0$). On détermine $\boldsymbol{L}$ et $\boldsymbol{U}$ de l'égalité $\boldsymbol{A} = \boldsymbol{L}\boldsymbol{U}$, c'est-à -dire de :

On commence par déterminer tout d'abord les lignes de $\boldsymbol{U}$ et les colonnes de $\boldsymbol{L}$ par identification : \[ u_{11} = a_{11}, \hspace{.25cm}\boldsymbol{U}_{12} = \boldsymbol{A}_{12}\hspace{.25cm}\text{et}\hspace{.25cm}\boldsymbol{L}_{21} = \frac{\boldsymbol{A}_{21}}{a_{11}} \]

Il ne reste plus qu'à  factoriser la matrice $(n-1) \times (n-1)$ :

Un exemple permet d'apprécier la démarche. Considérons la matrice :

On applique les différentes formules pour obtenir :

ce qui signifie que :

Pour finir

La seconde ligne donne naissance à  trois égalités. La première $2/3 = l_{32}(8/3)$ donne $l_{32} =1/4$. Après substitution dans les deux autres équations, on a :

ce qui montre que la méthode est récursive et enfin, comme $0 = l_{42}(8/3)$, ce qui implique que $l_{42} = 0$ et :

Ainsi :

Un dernier point : la décomposition $\boldsymbol{L}\boldsymbol{U}$ des matrices n'est pas nécessairement unique.

$\texttt{SageMath}$ a une commande qui permet d'obtenir la décomposition $\boldsymbol{L}\boldsymbol{U}$ d'une matrice. Comme cette commande possède des options, il est nécessaire d'être précis.

On notera que les matrices peuvent être différentes selon la méthode retenue comme le montrent les exemples d'application suivants (empruntés à la documentation). On commence ave l'option $\mathtt{partial}$ (ne pas oublier les guillemets pour passer une chaîne de caractères).

Puis, on continue avec l'option $\mathtt{nonzero}$

Puis, on continue avec l'option $\mathtt{compact}$

On note que le format $\mathtt{compact}$ n'est pas très utile pour la suite du cours. Néanmoins, si on veut comprendre le lien entre les formats, on peut consulter la doc à la page Base class for matrices, part 2

La question ici est différente de celle de la sous-section précédente qui était consacrée aux solutions des systèmes d'équations (qui se présentent avec des égalités). Ici, il ne s'agit pas de trouver une solution pour elle-même, mais de montrer que le système en admet au moins une ou n'en admet pas, i.e. : $\boldsymbol{A} \boldsymbol{x}\leq\boldsymbol{b}$ pour $\boldsymbol{x}\geq \boldsymbol{0}$. Évidement, si $\boldsymbol{b} \geq 0$ et que toutes les entrées de $\boldsymbol{A}$ sont positives, le problème est immédiatement résolu puisque $\boldsymbol{x} = \boldsymbol{0}$ est admissible. Comment faire pour trouver une solution?

Le premier auteur a s'être inquiété de ce problème semble avoir été Joseph Fourier^* J. Fourier (1826). Solution d'une question particuli\`ere du calcul des inégalités. Nouveau Bulletin des Sciences par la Société Philomatique de Paris, 99-100. . Il fut relayé par George Boole^* G. Boole (1854). On the conditions by wich the solutions of questions in the theory of probabilities are limited. The Philosophical Magazine, VIII. et au XX^e siècle par par Loyd L. Dines^* L. L. Dines (1918). System of linear inequalities. Annals of Mathematics, 20, 191-199. , Theodore Motzkin^* T. Motzkin (1936). Beiträge zur Theorie der linearen Ungleichungen. Université de Bale. puis, conjointement, par David Avis et Bohdan Kaluzni^* D. Avis, & B. Kaluzny (2004). Solving inequalities and proving Farkas' lemma made easy. The American Mathematical Monthly, 2, 152-157. .

On doit tout d'abord réaliser la principale différence entre une égalité et une inégalité. Dans le cas de la première, multiplier les deux membres par un nombre négatif ne change pas la nature de l'inégalité. En revanche, en présence d'une inégalité du type $\geq$, la multiplication par un nombre négatif inverse le sens de l'inégalité vers $\leq$ (et réciproquement si l'on part de l'inégalité opposée).

inégalités qui peuvent s'écrire de manière plus compacte comme :

ou encore $0\leq x \leq 5/3$. On a ainsi réussi à  remplacer 4 inégalités par deux inégalités qui caractérisent le même problème. Globalement, l'idée de l'élimination de Fourier-Motzkin consiste à  dériver du système initial un système d'inégalités dans lequel une variable a été éliminée (le choix de cette variable est arbitraire. En itérant cette démarche un nombre suffisant de fois, on peut aboutir sur un système à  une, voir deux variables et en déduire la faisabilité du système initial).

Avant de développer la méthode pour plusieurs variables, il est intéressant de développer une fonction qui permette de retrouver ce résultat dans $\texttt{SageMath}$.

Il est facile de demander à $\texttt{SageMath}$ de faire tourner la méthode de Gauss sur une matrice. Dans une discution du forum $\texttt{Ask SageMath}$ (pour avoir accès à ce forum et/ou poser une question, il suffit de passer par une barre de recherche), le code qui suit a été proposé. Comme, ce code est utilisé uniquement pour des raisons pédagogiques, on ne s'intéresse qu'à des matrices de rationnels.

$\texttt{Code}$

On peut donc appliquer cette fonction sur l'exemple suivant :

De manière générale, dans $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$, on peut avoir les trois types de lignes suivantes :

Oublions pour l'instant le dernier type de lignes (en fait, il est aisé de les remplacer par deux lignes identiques à  l'inégalité près). On peut alors isoler $x_l$ dans les deux premiers types d'inégalités pour obtenir :

ou, en condensant les notations :

Supposons maintenant que, parmi les $m$ équations qui composent le système, il y en ait $m_1$ qui ont $a_{kl} >0$ et $m_2$ avec $a_{kl} <0$ (comme on l'a signalé précédemment on ne tient pas compte de $a_{kl} =0$). Il est alors évident que $x_k$ appartient à  l'intervalle :

On peut éliminer $x_k$ et conserver simplement l'inégalité

d'où l'on tire une nouvelle inégalité concernant toutes les autres variables $x_p$ ($x_l$ étant exclue) que l'on rajoute aux $m$ équations pour lesquelles $a_{kl} = 0$.

Notons que si $\max\{a_1, a_2, \ldots,a_{m_1}\} \leq \min\{b_1, b_2, \ldots,b_{m_2}\}$ alors :

Donc, l'équation~\ref{fme1} permet de définir un ensemble de $m_1\times m_2$ inégalités pour $n-1$ variables. Notons qu'un polyèdre $\mathcal{P}\in \mathbb{R}^n$ est l'intersection d'un nombre fini de demi-hyperplans $\boldsymbol{a}_i^\top\boldsymbol{x}\leq \boldsymbol{b}_i$ ($\boldsymbol{a}_i$, $\boldsymbol{x}$ et $\boldsymbol{b}_i$ étant tous des vecteurs de $\mathbb{R}^n$). En d'autres termes :

L'opération d'élimination de Fourier-Motzkin part d'un polyèdre dans $\mathbb{R}^n$ puisqu'un ensemble de demi-hyperplans est donné et par élimination d'une variable aboutit à  un autre ensemble d'hyper-plans donc à  un autre polyèdre dans $\mathbb{R}^{n-1}$. En opérant ainsi, on réalise une projection de $\mathbb{R}^n$ sur $\mathbb{R}^{n-1}$. Or, il est possible de démontrer que l'objet projeté est lui-même un polyèdre. On peut donc recommencer sur une autre variable (en général, pour simplifier, on opère dans l'ordre croissant des indices) et réaliser une projection de $\mathbb{R}^{n-1}$ sur $\mathbb{R}^{n-2}$ (voir la figure suivante pour les diverses projections de $\mathbb{R}^3$ dans $\mathbb{R}^2$). On opère ainsi de suite jusqu'à  ce qu'il ne demeure qu'une seule variable.

Tuto : Polyhèdres

Le graphique ci-dessus représente un cube. C'est un volume assez simple, car en fonction de sa position dans l'espace, il est assez facile de repérer ses coordonnées de ses sommets. Par contre, il est peut-être plus difficile de déterminer les demi-plans dont il est l'intersection. On reviendra plus loin sur ce point. Il est ici assez évident que son intérieur n'est pas vide. Appliquer la méthode de Fourier-Motzkin sur ce cube peut se faire de différentes manières (certaines flèches ne sont visibles que lorsqu'on effectue une rotation du graphique à l'aide de la souris) :

1. On peut éliminer $x_3$ (flèche bleue). Après quoi, soit on élimine $x_2$ (flèche rouge), soit $x_1$ (flèche jaune) ;
2. On peut éliminer $x_2$ (flêche magenta). Après quoi, soit on élimine $x_3$ (flèche verte), soit $x_1$ (flèche vert olive) ;
3. On peut éliminer $x_1$ (flêche grise). Après quoi, soit on élimine $x_2$ (flèche cyan foncée), soit $x_3$ (flèche noire).

Dans tous les cas, on aboutit sur un un intervale final qui est clairement défini puisque les inégalités qui caractérisent le cube sont compatibles les unes avec les autres (le cube a un intérieur et une frontière clairement établis). Cependant, on notera que des inégalités ne déterminent pas toutes des polyèdres (solides délimités de toute part par des polygones plans). On peut commencer par le cas simple suivant :

On isole $x_1$ :

ce qui implique que

inégalité qui donne naissance aux deux inéquations :

soit

Ceci signifie que pour tout $x_2 \in [0, 12/5]$, il existe $x_1$ qui satisfait le système initial des trois inégalités comme cela est représenté à  la figure suivante (ici encore, il est important de lire le code pour savoir comment obtenir une représentation avec $\texttt{SageMath}$).

$\texttt{Code}$

On peut donc appliquer cette fonction sur l'exemple suivant :

On peut maintenant passer à des systèmes de plus grande dimension. On peut commencer par le cas simple suivant :

On isole $x_1$ :

ce qui implique que

inégalité qui donne naissance aux deux inéquations :

soit

Ceci signifie que pour tout $x_2 \in [0, 12/5]$, il existe $x_1$ qui satisfait le système initial des trois inégalités comme cela est représenté à  la figure suivante

Pour simplifier, on part du dodécahèdre (voir le tuto sur les polyèdres, plus haut), mais on pourrait partir de n'importe lequel des polyèdres. La fonction suivante permet d'éliminer une dimension au choix ($x_0$, $x_1$ ou $x_2$).

Après quelques instants de réflexion, il est évident que la méthode peut être appliquée pour déterminer la solution d'un problème linéaire. En effet, supposons qu'il s'agisse d'un problème de maximisation du type :

Ceci signifie que $z\leq \sum_{i=1}^n a_i x_i$ ou, plus simplement, que $z- \sum_{i=1}^n a_i x_i\leq 0$. Inversement, si on est face à un problème de minimisation

on aura $z\geq \sum_{i=1}^n a_i x_i$ ou, plus simplement, que $z- \sum_{i=1}^n a_i x_i\geq 0$. On peut dans les deux cas traiter l'objectif comme une simple inégalité supplémentaire et si on prend soin à conserver $z$ jusqu'à la dernière étape, on obtiendra une borne de faisabilité qui, à l'égalité, donnera la valeur souhaitée de l'objectif.

Considérer les inégalités $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$ définies par :

\[ \boldsymbol{A} = \begin{bmatrix}0 & 4\\0 & -35\\0 & 12\\0 & -1\end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm} \boldsymbol{b}=\begin{bmatrix}100\\-3\\72\\40\end{bmatrix}\hspace{.25cm}\text{pour}\hspace{.25cm}\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\end{bmatrix} \]

Existe-t-il des valeurs de $x_2$ qui résolvent ces inégalités ?

☕_Réponse 🕵️‍

Considérer les inégalités $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$ définies par :

\[ \boldsymbol{A} = \begin{bmatrix}3 & 4\\5 & 6\\4 & -7\\2 & -3\end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm} \boldsymbol{b}=\begin{bmatrix}14\\50\\18\\-8\end{bmatrix}\hspace{.25cm}\text{pour}\hspace{.25cm}\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\end{bmatrix} \]

1. Dans le plan $(x_1, x_2)$ donner une représentation graphique du domaine définit par $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$.
2. En éliminant $x_1$ par la méthode de Fourier-Motzkin, donner le domaine de validité de $x_2$.
3. En éliminant $x_2$ par la méthode de Fourier-Motzkin, donner le domaine de validité de $x_1$.

☕_Réponse 🕵️‍

Considérer les inégalités $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$ définies par :

\[ \boldsymbol{A} = \begin{bmatrix}1&-2 & -3\\0 & 1 & -2\\0&2&1\\0& 0 & 1\\0 &-1&0\\0 & 0 & -1\end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm} \boldsymbol{b}=\begin{bmatrix}0\\4\\18\\10\\0 \\0\end{bmatrix}\hspace{.25cm}\text{pour}\hspace{.25cm}\boldsymbol{x} = \begin{bmatrix}z\\x_1 \\x_2\end{bmatrix} \]

1. En éliminant $x_2$ et $x_1$ dans un ordre quelconque, montrer qu'il existe une borne qui limite $z$.
2. Si $z$ prend cette valeur limite, la dernière inégalité doit être une égalité. Mais, de ce fait, les équations qui l'engendrent doivent-elles même être des équations. En déduire la valeur optimale de $x_1$ et $x_2$.
3. Considérer le programme suivant :
\[ \max_{\{x_1,x_2\}} \left\{\left. 2x_1+ 3 x_2\right| x_1- 2x_2\leq 4, 2x_1+x_2\leq 18,x_2 \leq 10, x_1 \geq 0,x_2\geq0\right\} \]
1. donner une représentation graphique de l'ensemble des contraintes.
2. Donner une représentation graphique de l'objectif.
3. Quelle est la valeur optimale de $z$ et de ses arguments $x_1$ et $x_2$ ?
4. Qu'en conclure ?

☕_Réponse

Considérer les inégalités $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$ définies par :

\[ \boldsymbol{A} = \begin{bmatrix}1&-5 & -1\\0 & 2 & 1\\0&0&1\\0& 2 & 3\\0 &0&-1\\0 & 0 & -1\end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm} \boldsymbol{b}=\begin{bmatrix}0\\4\\18\\10\\0 \\0\end{bmatrix}\hspace{.25cm}\text{pour}\hspace{.25cm}\boldsymbol{x} = \begin{bmatrix}z\\x_1 \\x_2\end{bmatrix} \]

1. En éliminant $x_2$ et $x_1$ dans un ordre quelconque, montrer qu'il existe une borne qui limite $z$.
2. Si $z$ prend cette valeur limite, la dernière inégalité doit être une égalité. Mais de ce fait, les équations qui l'engendrent doivent-elles même être des équations. En déduire la valeur optimale de $x_1$ et $x_2$.
3. Considérer le programme suivant :
\[ \max_{\{x_1,x_2\}} \left\{\left. 2x_1+ 3 x_2\right| x_1- 2x_2\leq 4, 2x_1+x_2\leq 18,x_2 \leq 10, x_1 \geq 0,x_2\geq0\right\} \]
1. Donner une représentation graphique de l'ensemble des contraintes.
2. Donner une représentation graphique de l'objectif.
3. Quelle est la valeur optimale de $z$ et de ses arguments $x_1$ et $x_2$ ?
4. Qu'en conclure ?

☕_Réponse

Trouver l'intervalle de validité de la fonction $u = x_1+x_2$ quand $\boldsymbol{x}= \begin{bmatrix}x_1 & x_2\end{bmatrix}^\top$ est contraint par le système $\boldsymbol{A}\boldsymbol{x}\leq \boldsymbol{b}$ où les constantes sont définie par

\[ \boldsymbol{A} = \begin{bmatrix}-2 & -1\\-1 & 2\\3 & -1\\3 & 1\end{bmatrix}\hspace{.25cm}\text{et}\hspace{.25cm} \boldsymbol{b}=\begin{bmatrix}20\\-10\\60\\80\end{bmatrix}\hspace{.25cm}\text{pour}\hspace{.25cm}\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\end{bmatrix} \] ☕_Réponse 🕵️‍

4.7.2 La méthode d'Avis et kaluzni ☝

D. Avis et B. Kaluzny^* D. Avis, & B. Kaluzny (2004). Solving inequalities and proving Farkas' lemma made easy. The American Mathematical Monthly, 2, 152-157. ont proposé une autre méthode qui se rapproche de la méthode canonique de la programmation linéaire. Soit le système d'inégalités :

\[ \boldsymbol{A} \boldsymbol{x} \leq \boldsymbol{b}, \hspace{.25cm}\boldsymbol{x}\geq \boldsymbol{0} \]

où $\boldsymbol{A}$ est une matrice $m\times n$, $\boldsymbol{x}$ et $\boldsymbol{b}$ des vecteurs $n\times 1$. On introduit un vecteur $\boldsymbol{e}$ de variables d'écart non négatives ---~\textit{i.e.} : $\boldsymbol{e}\geq \boldsymbol{0}$~--- qui permettent de transformer l'inégalité en égalité, i.e. : \[ \boldsymbol{A} \boldsymbol{x} +\boldsymbol{e} = \boldsymbol{b} \hspace{.25cm}\Longleftrightarrow\hspace{.25cm}\boldsymbol{e} = \boldsymbol{b}- \boldsymbol{A} \boldsymbol{x} \]

Un tel système s'appelle un dictionnaire. Les variables qui sont dans le membre de gauche de l'égalité (pour l'instant les $e_i$) sont appelées variables basiques. Les variables qui sont dans le membre de droite sont appelées {variables hors-base. On remarque que toute solution positive du système d'égalités s'étend aisément au système initial d'inégalités.

Maintenant, de deux choses l'une, ou bien $\boldsymbol{b}> \boldsymbol{0}$, ce qui signifie que tous les éléments qui composent $\boldsymbol{b}$ sont positif, ou $\boldsymbol{b}$ est constitué d'éléments positifs ou nuls. Dans le premier cas, si on pose $\boldsymbol{x} = \boldsymbol{e} = \boldsymbol{0}$, l'inégalité est vérifiée trivialement puisque $\boldsymbol{0}= \boldsymbol{A}\boldsymbol{x} \leq \boldsymbol{b}$ et $\boldsymbol{x} = \boldsymbol{0}\geq \boldsymbol{0}$. Dans le second cas, on a une solution immédiate : $\boldsymbol{e} = \boldsymbol{b}$ et $\boldsymbol{x} = \boldsymbol{0}$.\ Elle n'est pas admissible puisque dans ce cas, tous les $e_i$ ne sont pas positifs. L'algorithme suggérée par \citet{AvisKaluzny04}\sindex[aut]{Avis D.}\sindex[aut]{Kaluzny B.} consiste à  sauter de solutions irréalisables voisines jusqu'à  ce qu'éventuellement, si le système est solvable, on aboutisse sur une solution réalisable ---~on dit encore \textit{faisable}. On commence par unifier les indices de telle sorte à  ce que \[ i = \underbrace{1, 2, \ldots, n}_{x_i}, \underbrace{n+1,n+2,\ldots, 2n}_{e_i} \]

Ainsi, quand on parle des indices, toutes les variables sont confondues ---~\textit{i.e.} : quand on parle du plus petit indice, il s'agit conjointement des indices des $x_i$ et des $e_i$. L'algorithme proposé est le suivant (nommé la Règle-b) :

  Algorithme d'Avis & Kaluzni

1. On part de $\boldsymbol{e} = \boldsymbol{b}$ (ce qui implique que $\boldsymbol{x} = \boldsymbol{0}$).
2. On choisit le $e_i$ ayant le plus petit indice $i$ tel que $e_i<0$.
3. Dans l'équation qui détermine $e_i$, on choisit le $x_i$ ayant le plus petit indice associé à  un coefficient positif. Notons le $x_{\underline{i}}.$ S'il n'y en a pas, le système est infaisable. S'il y en a au moins un :
4. On résout cette équation en $x_{\underline{i}}$.
5. On substitue à  $x_{\underline{i}}$ dans toutes les autres équations sa nouvelle définition de telle sorte que toutes ces équations ne dépendent plus de $x_{\underline{i}}$. On obtient alors un nouveau dictionnaire hybride dans lequel les variables de base et hors base sont composées de variables d'écart et de variables initiales.
6. On recommence jusqu'à  ce qu'éventuellement, en fixant toutes les variables hors base à  $0$, toutes les variables de base puisent être positives, auquel cas on a montré que le système est faisable.

Considérons l'exemple suivant emprunté à 

D. Avis et B. Kaluzny^* D. Avis & B. Kaluzny (2004). Solving inequalities and proving Farkas' lemma made easy. The American Mathematical Monthly, 2, 152-157. :

\[ \begin{array}{l} -x_1 - 2x_2 + x_3 \leq -1\\ x_1 - 3x_2 - x_3 \leq 2\\ -x_1 - 2x_2 + 2 x_3 \leq -2 \end{array} \]

On commence par transformer ce système en égalités^!On doit faire attention à  l'effet de présentation qui est lié à  l'isolation des variables et qui force à  inverser le signe des coefficients associés aux variables hors base. :

\[ \begin{array}{l} {\color{red}{e_4}} = -1 + {\color{cyan}{x_1}} + 2 x_2 -x_3\\ e_5 = 2 -x_1+\ 3 x_2 + x_3\\ e_6 = -2 +x_1 + 2\ x_2 - 2x_3 \end{array} \]

La solution $\begin{bmatrix}e_4 & e_5 & e_6\end{bmatrix}^\top = \begin{bmatrix}-1 & 2 & -2\end{bmatrix}^\top, \begin{bmatrix}x_1 & x_2 & x_3\end{bmatrix}^\top = \begin{bmatrix}0 & 0 & 0\end{bmatrix}^\top$ n'est pas faisable puisque $e_1$ et $e_3$ y prennent des valeurs négatives. Considérons la quand même. Comme dans ce cas, $e_4$ est la variable de base ayant le plus petit indice prenant une valeur négative et comme, dans cette équation, $x_1$ est la variable hors base ayant le plus petit indice associé à  un coefficient positif on résout cette équation en $x_1$, ce qui donne :

\[ x_1 = 1 -2x_2 + x_3 + e_4 \]

Par substitution dans la seconde équation, on obtient :

\[ e_5 = 2 -\left(1 -2x_2 + x_3 + e_4\right)+\ 3 x_2 + x_3 = 1+5x_2 -e_4 \]

et

\[ e_6 = -2 + \left(1 -2x_2 + x_3 + e_4\right) + 2\ x_2 - 2x_3 =-1 -x_3+e_4 \]

ce qui donne le système :

\[ \begin{array}{l} x_1 = 1 -2x_2 + x_3 + e_4\\ e_5 = 1+5x_2 -e_4\\ {\color{red}{e_6}} = -1 -x_3+{\color{cyan}{e_4}} \end{array} \]

Il n'est pas possible de s'arrêter là, car la solution $\begin{bmatrix}e_4 & e_5 & e_6\end{bmatrix}^\top = \begin{bmatrix}0 & 1 & -1\end{bmatrix}^\top$, $\begin{bmatrix}x_1 & x_2 & x_3\end{bmatrix}^\top = \begin{bmatrix}1 & 0 & 0\end{bmatrix}^\top$ n'est pas faisable puisque $e_6 = -1.$ Donc, on réitère. La variable de base prenant une valeur négative ayant le plus petit indice est $e_6$. Dans l'équation qui définit $e_6$, la variable hors base dont le coefficient est positif ayant le plus petit indice est $e_4$. On fait entrer $e_4$ dans la base et sortir $e_6$, soit :

\[ e_4 = 1+ x_3+ e_6\]

Par substitution, on obtient :

\[ x_1 = 1 -2x_2 + x_3 +\left( e_6 +\ 1 + x_3\right) = 2 \ -2 x_2+2x_3+e_6 \]

et

\[ e_5 = 1+5x_2 -\left(e_6 +\ 1 + x_3\right) =0 +5 x_2-x_3-e_6 \]

ce qui donne le système suivant :

\[ \begin{array}{l} x_1 = 2 \ -2 x_2+2x_3+e_6\\ e_5=0 +5 x_2-x_3-e_6\\ e_4 = 1+ x_3+ e_6\end{array} \]

Nous pouvons nous arrêter là. En effet, en fixant les variables de bases à  $x_1$, $e_5$ et $e_4$, l'attribution des valeurs $\begin{bmatrix}e_4 & e_5 & e_6\end{bmatrix}^\top = \begin{bmatrix}1 & 0 & 0\end{bmatrix}^\top$, $\begin{bmatrix}x_1 & x_2 & x_3\end{bmatrix}^\top = \begin{bmatrix}2 & 0 & 0\end{bmatrix}^\top$ donne une solution faisable. Il existe donc bien une solution au problème. Appliquons la méthode au système :

\[ \begin{array}{l} -x_1 + 2x_2 + x_3 \leq 3\\ 3 x_1 - 2x_2 + x_3 \leq -17\\ -x_1 - 6x_2 + 23 x_3 \leq 19 \end{array} \]

encore une fois emprunté à  D. Avis & B. Kaluzny. Commençons par introduire les variables d'écart $e_4$, $e_5$ et $e_7$. \[ \begin{array}{l} -x_1 + 2x_2 + x_3 + e_4 = 3\\ 3 x_1 - 2x_2 + x_3 + e_5 = -17\\ -x_1 - 6x_2 + 23 x_3 + e_6 = 19 \end{array} \,\,\,\,\,\Longleftrightarrow\,\,\,\,\, \begin{array}{l} e_4 = 3+x_1 - 2x_2 - x_3 \\ {\color{red}{e_5}} = -17-3 x_1 + 2{\color{red}{x_2}} - x_3 \\ e_6 = 19 +x_1 + 6x_2 - 23 x_3 \end{array} \]

Comme précédemment des valeurs $\begin{bmatrix}e_4 & e_5 & e_6\end{bmatrix}^\top = \begin{bmatrix}3& -17 & 19\end{bmatrix}^\top$, $\begin{bmatrix}x_1 & x_2 & x_3\end{bmatrix}^\top = \begin{bmatrix}0 & 0 & 0\end{bmatrix}^\top$ ne constituent pas une solution faisable du fait que $e_5$ prend une valeur négative. Comme $x_2$ a un coefficient positif, il rentre dans la base, ce qui donne :

\[ x_2 = \frac{1}{2}\left(e_5 +\ 17 +\ 3x_1 + x_3\right) \]

Par substitution, on a :

\[ e_4 = 3+x_1 - 2\left( \frac{1}{2}\left(e_5 +\ 17 +\ 3x_1 - x_3\right)\right) - x_3 =-14- 2x_1-2\ x_3-e_5 \]

et

\[ e_6 = 19 +x_1 + 6\left(\frac{1}{2}\left(e_5 +\ 17 +\ 3x_1 - x_3\right)\right) - 23 x_3=70 +4x_1-3 x_3 + 3 e_5 \]

Ainsi, on a le système :

\[ \begin{array}{c} e_4=-14- 2x_1-2\ x_3-e_5\\ x_2 = \displaystyle{\frac{17}{2}+\frac{3}{2}x_1+\frac{3}{2}x_3+\frac{1}{2}e_5}\\ e_6 = 70 +4x_1-3 x_3 + 3 e_5 \end{array} \]

mais la solution $\begin{bmatrix}e_4 & e_5 & e_6\end{bmatrix}^\top = \begin{bmatrix}-14& 0 & 70\end{bmatrix}^\top$, $\begin{bmatrix}x_1 & x_2 & x_3\end{bmatrix}^\top = \begin{bmatrix}0 & 17/2 & 0\end{bmatrix}^\top$ n'est pas plus acceptable puisque $e_4 \leq 0$. Il serait donc bien venu de sortir $e_4$ de la base mais, cette fois, on ne sait comment faire puisque tous les coefficients du membre de droite de l'équation de $e_4$ sont négatifs. Cette équation est inconsistante dans la mesure où $e_4 \geq 0$, $x_1 \geq0$, $x_3 \geq 0$ et $e_5 \geq0$. Elle ne peut donc être satisfaite: \[ 0\geq e_4=-14- 2x_1-2\ x_3-e_5\leq -14 \]

Étant donnée la simplicité de cet algorithme, il est particulièrement intéressant de montrer ses propriétés :

1. il s'arrête après un nombre fini de pas.
2. S'il termine à  l'étape 6 de sa description page (voir l'encadré ci-dessus), la solution de base proposée est faisable.
3. S'il termine à  l'étape 3, elle est infaisable.

Pour discuter de ces propriétés, nous allons simplifier les notations en notant les variables d'écart comme les variables initiales (i.e. : $e_i= x_i$, $i = n+1, \ldots, n+m$). Par définition, le système qui doit être analysé comporte $m$ inégalités et $n$ variables. Après transformation en système d'égalités, il comporte toujours $m$ inégalités mais, cette fois, il a $n+m$ variables. Au plus a-t-il donc $\complement_{n+m}^n$ vecteurs candidats pour constituer une base. Si l'algorithme ne s'arrête pas en un nombre fini de pas, c'est donc qu'après avoir exploré toutes les bases potentielles, il revient sur une base déjà  explorée. En d'autres termes, il cycle.

Supposons que le système analysé cycle. Puisque toutes les variables entrent un moment donné dans la base, considérons la dernière $x_{n+m}$. Au moment où elle entre dans la base, une des équations doit être du type :

\begin{equation} \color{cyan}{x_k = -b^\prime_k -\sum_{j\in \mathcal{N}\setminus\{n+m\}} a^\prime_{jk} x_k + a^\prime_{j{n+m}} x_{n+m}} \end{equation}

Précisons : $k\in \mathcal{B}$, où $\mathcal{B}$ est l'ensemble des variables de base au moment du choix de $x_{n+m}$ comme variable entrante et $\mathcal{N}$ est l'ensemble des variables qui ne sont pas dans la base. Puisque $x_k$ va sortir, $-b^\prime_k< 0$. Puisque $x_{n+m}$ entre, elle doit posséder l'indice le plus petit associé à  un coefficient strictement positif. Donc, pour tout $j\in \mathcal{N}/\{n+m\}$, on doit avoir $-a^\prime_{jk}\leq0$ et $a^\prime_{j{n+m}} >0$. En particulier, ceci montre que toute solution du système d'équations pour laquelle $x_1\geq 0,\ldots x_{n+m-1} \geq 0$ doit vérifier $x_{n+m} > 0$.

Maintenant, considérons le cas où $x_{n+m}$ vas être choisit pour quitter la base. On doit avoir :

\begin{eqnarray} x_i &=& b^\prime_i + \sum_{j \in \mathcal{N}} a_{ij} x_j, \hspace{.25cm}i \in \mathcal{B}\setminus \{n+m\}\label{ak2}\\ x_{n+m} &=& -b^\prime_{n+m} +\sum_{j \in \mathcal{N}} a_{(n+m)j} x_j\label{ak3} \end{eqnarray}

Puisque $x_{n+m}$ doit sortir, on a nécessairement $b^\prime _i \geq 0$ pour tout $i \in \mathcal{B}\setminus \{n+m\}$ et $-b^\prime_{n+m}< 0$. Si l'on fixe les variables hors base à  $0$ (pour tout $i \in \mathcal{B}\setminus \{n+m\}$, $x_i = 0$, le système d'équations précédent montre qu'il existe une solution telle que $x_1\geq 0, \ldots x_{n+m -1}\geq 0$ et $x_{n+m} < 0$). Mais, ceci contredit ce qui a été montré à  partir de l'équation écrite en couleur. Il n'est donc pas possible que l'algorithme cycle sur la variable $x_{n+m}$.

Maintenant, supposons qu'il existe un cycle dans lequel $x_{n+m}$ demeure toujours dans la base. On peut donc retirer l'équation correspondant à  $x_{n+m}$ sans changer les valeurs des pivots durant le cycle et il en va de même si pendant le cycle $x_{n+m}$ demeure toujours une variable hors base. Dans les deux cas, on peut réduire le système en éliminant une contrainte et redémarrer le raisonnement initial sur le nouveau système et ainsi de suite. Ceci montre que l'algorithme est bien fini (i.e. : stoppe après un nombre de pas fini). Par conséquent, l'algorithme soit s'arrête par ce qu'il a trouvé une solution, soit parce qu'il donne un certificat d'infaisabilité.

En terme générique, la finitude de l'algorithme donne une démonstration simple du lemme de Farkas.

Avec la méthode d'Avis et Kaluzni, vérifier s'il existe au moins une solution au système d'inéquations $\boldsymbol{A}\boldsymbol{x} \leq \boldsymbol{b }$ défini par les matrices et vecteurs suivants :

1. $\boldsymbol{A} = \begin{bmatrix}1 & 1&-2\\2 & -1&-1\\-1 & 2& 1\\-1 & -1 & 1\end{bmatrix}\,\,\,\,\,\text{et}\,\,\,\,\, \boldsymbol{b}=\begin{bmatrix}0\\2\\1\\-1\end{bmatrix}\,\,\,\,\,\text{pour}\,\,\,\,\,\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\\x_3\end{bmatrix}$
2. $\boldsymbol{A} = \begin{bmatrix}1 & -2\\1 & 1\\1 & 0\\-2 & 1\\-1 & 0\\ 0& 8\end{bmatrix}\,\,\,\,\,\text{et}\,\,\,\,\, \boldsymbol{b}=\begin{bmatrix}-2\\3\\2\\0\\-1\\11\end{bmatrix}\,\,\,\,\,\text{pour}\,\,\,\,\,\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\end{bmatrix}$
3. $\boldsymbol{A} = \begin{bmatrix}3 & 4\\5 & 6\\4 & -7\\2 & -3\end{bmatrix}\,\,\,\,\,\text{et}\,\,\,\,\, \boldsymbol{b}=\begin{bmatrix}14\\50\\18\\-8\end{bmatrix}\,\,\,\,\,\text{pour}\,\,\,\,\,\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\end{bmatrix}$
4. $\boldsymbol{A} = \begin{bmatrix}-1 & 3\\1 & 1\\2 & -1\end{bmatrix}\,\,\,\,\,\text{et}\,\,\,\,\, \boldsymbol{b}=\begin{bmatrix}3\\3\\6\end{bmatrix}\,\,\,\,\,\text{pour}\,\,\,\,\,\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\end{bmatrix}$
5. $\boldsymbol{A} = \begin{bmatrix}2 & 3&-4\\-1 & -2&1\\1 & 4&2\\0&3&5\end{bmatrix}\,\,\,\,\,\text{et}\,\,\,\,\, \boldsymbol{b}=\begin{bmatrix}5\\-3\\6\\2\end{bmatrix}\,\,\,\,\,\text{pour}\,\,\,\,\,\boldsymbol{x} = \begin{bmatrix}x_1 \\x_2\\x_3\end{bmatrix}$

☕_Réponse

Montrer, en appliquant la méthode d'Avis et Kaluzni que le système suivant est infaisable :

\[ \begin{bmatrix}-2 & 3 & -5\\-1 & -3 & 2\\1&1 &4\\2 & -2 & -2\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix} =\begin{bmatrix}5\\5\\4\\-10\end{bmatrix} \] ☕_Réponse

👈 ☝ 👉