simplex.html

On met tout ensemble

On peut avec ce qu'on a résoudre le problème de maximisation de c·x avec les contraintes x≥0 et A·x≤b (x vecteur à n composantes).

On cherche d'abord un sommet du polyèdre des contraintes avec "Sommet(A,b)", puis on utilise le sommet x0 et la base associée S pour "ChercheMax(Aaug, baug, c, x0, S)" où Aaug est A augmentée par l'opposée de la matrice identité dans les n premières lignes et baug est b augmenté par n zéros au début (ceci pour tenir compte de la contrainte de positivité sur x).

> MaxSimplexNonNeg:=proc(A,b,c)
local sommet, m, n, Aaug, baug;
sommet:=Sommet(A,b);
if sommet="pas faisable"
then return "pas faisable"
else m:=RowDimension(A); n:=ColumnDimension(A);
   Aaug:=Matrix(m+n,n); Aaug[1..n,1..n]:=-Matrix(n,n,shape=identity);
   Aaug[n+1..n+m,1..n]:=A;
   baug:=Vector(m+n); baug[n+1..n+m]:=b;
   ChercheMax(Aaug,baug,c,sommet[1],sommet[2])
end if
end proc:

Faisons un essai, en vérifiant avec le paquet "simplex":

> A:=Matrix(4,3,rand(-5..5));b:=Vector(4,rand(-5..5));
c:=Vector[row](3,rand(-5..5));
print(`Notre procédure donne`, MaxSimplexNonNeg(A,b,c));
Inc:=Vector([x,y,z]): Obj:=Multiply(c,Inc):
Contr:={seq(Multiply(Row(A,i),Inc) <= b[i], i=1..4)}:
print(`Le paquet de Maple donne`,simplex[maximize](Obj,Contr,NONNEGATIVE));

Enfin, il reste à traiter le problème général, sans condition de positivité pour les variables. Mais maximiser c·x avec A·x≤b revient (en posant x=y-z) à maximiser c·y - c·z avec y≥0, z≥0 et A·y - A·z ≤ b, problème de maximisation avec variables positives. Il suffit d'utiliser cette transformation pour avoir une procédure traitant le problème général.

> MaxSimplex:=proc(A,b,c)
local m,n,Anonneg,cnonneg,res ;
m:=RowDimension(A); n:=ColumnDimension(A);
Anonneg:= Matrix(m,2*n);
Anonneg[1..m,1..n]:=A ; Anonneg[1..m,n+1..2*n]:=-A;
cnonneg:=Vector[row](2*n); cnonneg[1..n]:=c; cnonneg[n+1..2*n]:=-c;
res:=MaxSimplexNonNeg(Anonneg,b,cnonneg);
[res[1], SubVector(res[2],1..n)-SubVector(res[2],n+1..2*n)]
end proc:

Un essai :

> A:=Matrix(7,3,rand(-5..5));b:=Vector(7,rand(1..5));
c:=Vector[row](3,rand(-5..5));
print(`Notre procédure donne`, MaxSimplex(A,b,c));
Inc:=Vector([x,y,z]): Obj:=Multiply(c,Inc):
Contr:={seq(Multiply(Row(A,i),Inc) <= b[i], i=1..7)}:
print(`Le paquet de Maple donne`, simplex[maximize](Obj,Contr));