The 4565pt2RevisedSimplexMethod.nb downloadable on internet is a step by step implementation of the revised simplex method very usefull for teaching. I am trying to typeset it in french. Until Iteration 2 all is fine, but from In[32] the behavior become weird finishing in an Out[39]= {x2, X2} in place of the true Out[39]= {x2, X1}. Here is my code (in french). I cannot understand why it doesn't work.
Thanks for the help
`Méthode du simplex révisée (dim2)
Initalisation
In[1]:= Off[General::spell1] Off[General::spell]
Données du Programme
In[3]:= A = {{40, 80}, {6, 8}}; x = {x1, x2}; b = {560, 72}; xs = {x3, x4}; c = {24, 36}; X = Join[x, xs]
Out[8]= {x1, x2, x3, x4}
Itération 0
In[9]:= xB = {x3, x4}; B = IdentityMatrix[Length[xB]]; InverseB = Inverse[B]; InverseB.b
Out[12]= {560, 72}
In[13]:= cB = {0, 0}; Z = cB.InverseB.b
Out[14]= 0
Maintenant, on teste les coefficients des variables non-basiques dans le rang de la fonction objectif de manière à vérifier que l'on est pas arrivé à la (une) solution optimale --- auquel cas tous ces coefficients seront non-négatifs. Les coefficients de x1 et x2 sont donnés par :
In[15]:= cB.InverseB.A - c
Out[15]= {-24, -36}
La négativité de ces deux coefficients montre que la solution n'est pas optimale.
Itération 1
La base actuelle est
In[16]:= xB
Out[16]= {x3, x4}
Du text d'optimalité réalisé dans la première itération, il vient que les coefficients des variables non-basique dans la fonction objectif ne sont pas tous non-négatifs. On doit choisir la variable ayant le coefficient le plus négatif comme variable entrante. Dans notre cas -36 < -24 dont c'es la variable x2. Soit k l'indice de la variable entrante.
In[17]:= k = 2;
On applique le test du rapport minimum pour les valeurs courantes des coefficients de la variable entrante.
In[18]:= (InverseB.b/InverseB.A)[[All, k]]
Out[18]= {7, 9}
Le rapport le plus petit est 7, donc c'est x3 la variable sortante. Soit r sa ligne.
In[19]:= r := 1;
On permute la variable entrante et la variable sortante dans le vecteur xB.
In[20]:= xB[[r]] = x[[k]]
Out[20]= x2
La nouvelle base est :
In[21]:= xB
Out[21]= {x2, x4}
Pour obtenir la nouvelle B^-1, on doit construire la matrice E --- on note que E est un symbole réservé pour l'exponentielle dans Mathématica et qu' on doit le remplacer par e ---
In[22]:= Clear[[Eta], eta, e]; eta[i_] = If[i != r, -(InverseB.A)[[i, k]]/(InverseB.A)[[r, k]], 1/(InverseB.A)[[r, k]]]; [Eta] = Table[eta[i], {i, 1, Length[xB]}]
Out[24]= {1/80, -(1/10)}
In[25]:= e = IdentityMatrix[Length[xB]]; For[i = 1, i < Length[xB] + 1, e[[i, r]] = [Eta][[i]]; i++]; e
Out[27]= {{1/80, 0}, {-(1/10), 1}}
On peut évaleur la nouvelle valeur de B^-1 à partir de l'ancienne
In[28]:= InverseB = e.InverseB
Out[28]= {{1/80, 0}, {-(1/10), 1}}
Puis trouver la nouvelle valeur de xB
In[29]:= InverseB.b
Out[29]= {7, 16}
On trouve les coefficients des nouvelles variables de base dans la fonction objectif
In[30]:= cB = {36, 0}
Out[30]= {36, 0}
On trouve la nouvelle valeur de z
In[31]:= cB.InverseB.b
Out[31]= 252
On vérifie si la solution est optimale en regardant les coefficients des variables non-basiques, x1 et x3 dans la fonction objectif.
Le coefficient de x1 est
In[32]:= (cB.InverseB.A - c)[[1]]
Out[32]= -6
Le coefficient de x3 est
In[33]:= (cB.InverseB)[[1]]
Out[33]= 9/20
Puisque le coefficient de x1 n'est pas négatif, z n'est pas pas maximisé.
Itération 2
Les variables de base sont les suivantes :
In[34]:= xB
Out[34]= {x2, x4}
Du fait du test d'optimalité réalisé à l'itération précédente, les coefficients des variables hors base ne sont pas tous non-négatifs. On choisit la variable hors-base ayant le coefficient le plus négatif dans l'objectif comme variable entrante. Cette fois-ci il s'agit de x1. Soit k l'indice de la variable entrante.
In[35]:= k = 1;
On applique le test du rapport minimal en divisant l'entrée de droite par les coefficients correspondant de la variable entrante :
In[36]:= (InverseB.b/InverseB.A)[[All, k]]
During evaluation of In[36]:= Power::infy: Infinite expression 1/0 encountered. >>
Out[36]= {14, 8}
{14, 8}
{14, 8}
Le plus petit des rapports est 8 ce qui correspond à la variable x4 qui est de ce fait la variable sortante. Soit r la ligne correspondante.
In[37]:= r = 2;
On permute la variable entrante et la variable sortante
In[38]:= xB[[r]] = xB[[k]];
In[39]:= xB
Out[39]= {x2, x2}
Il y a une erreur que je ne comprends pas !!!
On cherche la nouvelle matrice B^-1
In[40]:= Clear[[Eta], eta, e];
In[41]:= eta[i_] = If[i != r, -(InverseB.A)[[i, k]]/(InverseB.A)[[r, k]], 1/(InverseB.A)[[r, k]]]; [Eta] = Table[eta[i], {i, 1, Length[xB]}]
During evaluation of In[41]:= Part::partw: Part 3 of {{1/2,1},{2,0}} does not exist. >>
Out[42]= {-(1/2) {{1/2, 1}, {2, 0}}[[3, 1]], -(1/2) {{1/2, 1}, {2, 0}}[[3, 1]]}`
