Simplex method: doSimplex[A, B]

Function to to max a functions subject to constraints. Assumptions: 

   * all functions are linear
   * all inequalities are of the form “lin_fun <= const” with a positive const 
   * all variables are constrained to be non-negative (but coefficients could be negative)
   
The coefficient matrix A and the vector of constants B are input to the doSimplex method. Note: the variable “tableau” will be created and overwritten if it already exists - this variable will contain the final answer.

First method: doPivot. It defines the basic pivot function to pivot around the element at [row, col]. Note: the matrix must be defined in a variable named tableau. That variable will be defined and the original values, if any, will be overwritten in the setup function, called by the doSimplex function.

doPivot[row_, col_] := 
 Module[{i}, 
  tableau[[row]] =tableau[[row]] / tableau[[row,col]];
  Do[
   If[i≠ row,
    If[tableau[[i,col]]≠0,
	tableau[[i]] = tableau[[i]] - tableau[[i,col]]*tableau[[row]]]],
 {i,1,Length[tableau]}]]

Now defining functions getPivotCol and getPivotRow to find the row and column of the pivot element. Finding the pivot row is more complicated, since you are dividing, which could result in a division by zero.

getPivotCol[A_]:=
Module[{lastRow},
lastRow =Drop[A[[-1]],-1]; 
Position[lastRow,Min[lastRow]][[1,1]]]

getPivotRow[A_, pivotCol_]:= 
Module[{quotient, denom, numer},
denom = Drop[A[[All,pivotCol]],-1];
numer = Drop[A[[All,-1]],-1];
quotient = numer;
Do[
If[(denom[[i]]≤0) || (numer[[i]]≤ 0),
quotient[[i]]=Infinity, 
quotient[[i]] = numer[[i]] / denom[[i]]],
{i,1,Length[denom]}];
Position[quotient,Min[quotient]][[1,1]]
]

The setup method setupTableau appends the proper identity matrix for the slack variables and defines the last column as the (positive) constants. For example, if the input matrix was

      1  2  3  4
      4  5  6  7
     -1 -2 -3 -4
     
 and the constant vector <10, 20, 0>, this would describe the Linear Optimization problem with 4 variables and 2 inequalities:
 
    Maximize x + 2y + 3z + 4w subject to
    
        x + 2y + 3z + 4w <= 10
       4x + 5y + 6z + 7w <= 20

setupTableau[a_, b_]:=
( nVars = Length[a[[1]]];
nEqs = Length[a] - 1;
tableau=Transpose[Join[Join[Transpose[a],Join[IdentityMatrix[nEqs+1],{b}]]]];
)

Defining functions checkColForAnswers and checkForAnswers to check the tableau for the answers.

checkColForAnswer[col_]:=
Module[{values},
values = tableau[[All, col]];
If[(Total[values]==1) && (Max[values]==1),
Position[values,1][[1,1]]]];

checkForAnswers[]:=
Module[{pos},
Do[
pos =checkColForAnswer[i];If [pos >0,Print["x",i," = ",tableau[[pos,-1]]],Print["x",i," = 0"]],
{i,1,nVars}];
Print["Max = ", tableau[[-1,-1]]]]

Now defining the main simplex method doSimplex: It takes a simplex matrix as input and creates the global variable tableau. That matrix will contain the final matrix and you can read off the answers. The function then calls “doPivot” as long as there’s a negative value in the bottom row.

doSimplex[a_, b_] :=
Module[{lastRow,pivotRow,pivotCol},
setupTableau[a,b];
While[ Min[tableau[[Length[tableau]]]] < 0,
pivotCol = getPivotCol[tableau];
pivotRow = getPivotRow[tableau, pivotCol];
doPivot[pivotRow,pivotCol]];
  checkForAnswers[]; ]

Let’s check our doSimplex function to find the maximum of 4x + 2y subject to 

   10x+ 5y <= 50
   12 x + y <= 48, 
    3x + 2y <= 18
    
(and of course x >= 0 and y >= 0)

Maximize[2x-3y+4z, 4x-3y+z≤3 && x+y+z ≤ 10 && 2x + y-z≤ 10 && x ≥0&&y≥0&&z≥0,{x,y,z}]

doSimplex[{{4,-3,1},{1,1,1},{2,1,-1},{-2,3,-4}}, {3, 10,10, 0}]

Maximize[{4x + 2y,
10x+5y≤ 50 && 12x +y≤48 && 3x+2y ≤18&& 
x ≥ 0 && y ≥ 0 },{x,y}]

doSimplex[{{10,5},{12,1},{3,2},{-4,-2}}, {50, 48, 18, 0}]

Finally we can test our program by generating random problems with nVars variables and nVars/2 equations. Each time you run this code, a new problem is generated and solved. In addition, our Simplex method is timed and the total time it took to computw the answer is shown.

nVars = 400; nEqs = nVars / 2;
randA=RandomInteger[10, {nVars,nEqs}];
randA[[-1]] = -randA[[-1]];
randB = Table[RandomInteger[50],{i,nVars-1}];
randB = Append[randB,0];
Timing[doSimplex[randA, randB]]