Creating a Linear Program Solver by Implementing the Simplex Method in Python with NumPy
By Matt Molter — Grad Student at The University of Texas at Austin
I. Introduction
The purpose of this project was to improve my understanding of the simplex method, specifically the tabular method and the 2-phase method. This post is structured as follows. Section II is a brief overview of the simplex method, including the 2-phase variant. Section III discusses my implementation. Section IV shows results. Section V is rough summary. Section VI contains references and links to a Github repo containing the notebooks with code.
For the discussion, we will primarily be talking about minimization problems in standard form as shown below.
II. Simplex Overview
The simplex method is an algorithm for solving linear programming problems. It is based on the fact that an optimal solution to a linear programming problem always lies at an extreme point. This fact leads to a number of ways to find the optimal solution, the most naive of which is to check the cost at each extreme point and whichever is lowest is your optimal solution. However, finding extreme points is computationally difficult, and for large problems, the number of extreme points can be arbitrarily large. This is obviously not ideal.
A more efficient algorithm is the simplex method, which involves finding any initial feasible solution, which can still be computationally difficult, but must be done only once. Once an initial feasible solution is found, you have an initial basis for your solution. In order to determine if there is a more optimal solution, the relative cost increase/decrease for moving along every edge adjacent to your initial point is calculated. If there is no direction in which you can move that involves a decrease in cost, you are at the optimal corner.
Otherwise, the direction in which cost decreases the most is obviously the direction in which you would move. In order to do so, you have to determine how far along the edge you can move, and then calculate your new basis and costs at the new point.
In order to formalize the above actions, we use what is referred to as the Simplex tableau. The tableau is of the format below. B is your basis matrix, how to determine it will be shown below.
The top left quadrant is the current cost of your optimization function. The top right quadrant is the relative costs for each variable. The bottom left quadrant is the value of the variables in the current basis, and the bottom right is your current constraint coefficients.
While this looks like quite a few matrix operations, the tableau can be assembled with less formalism, as will be shown in the implementation below. The matrix operation implementation is faster, but a bit more difficult to understand.
As an example, here is a problem in standard form.
In order to get the tableau, we need to put this in canonical form by adding slack variables, as seen below.
As luck has it, this also provides an initial basic feasible solution, where x = (0, 0, 0, 20, 20, 20). This will allow us to assemble the tableau as below.
It is easy to see how the tableau relates to the problem in canonical form. As can be seen by our initial basic feasible solution, the elements in the initial basis are x4, x5, x6. This means that there is no cost advantage to introducing them into the basis, since they are already in the basis, hence the 0 under these variables in the top row. The decrease in cost of variables outside of the initial basis can be seen to be negative, indicating that the tableau does not show an optimal solution.
In addition, the leftmost column can be seen to be made up of the resources from the righthand side of our constraints. Finally, the lower right section can be seen to be the lefthand side of our constraints in canonical form. Not too much linear algebra needed, right? The question remains, how do we turn this initial feasible solution into the optimal one?
We start by examining the relative cost row, as explained earlier. We notice that there are two variables which have identical decreases in cost. This doesn’t affect correctness, but we will choose the first one that appears to be the variable that we introduce. This is x2 in this case, so we say that x2 enters the basis, and the column it represents is our pivot column.
We can’t have a variable enter the basis without another exiting, so we must now determine which variable exits the basis. We do this by performing what is called the minimum ratio test. What this means is that for each row, we divide the element in the leftmost column by the element in the pivot column. We save these and compare them, and whichever is the lowest is indicative of the pivot row, and that corresponding variable leaves the basis. Again, we choose the first that appears if there are multiple with the same ratio, which does not affect correctness. In this case, this is x4, so we say that x4 leaves the basis.
Before we move on, it is also important to note that the minimum ratio test is performed on a row only if the corresponding element in the pivot column is not 0. Otherwise that corresponding variable will remain in the basis and there is no need to calculate the ratio, and in fact it will break the Simplex algorithm, for what should be an obvious reason
Now that we have our pivot column and row, we can also determine our pivot element, which is the value where the column and row intersect. In this case, it is 2. To begin with, we divide the pivot row by the pivot element in order to normalize it. We then go to each and every other row and perform elementary row operations on it to make the rest of the column values 0.
Once the pivot actions are complete, we are done with this iteration and have a new tableau with which to start the next iteration. It all sounds easy enough, but we’ve ignored one crucial part of the algorithm.
Finding the initial feasible solution is one of the more complicated parts of the Simplex method. In the most trivial case, as shown above, if you have a problem in standard form, and the right-hand side of your constraints are all >= 0, you can introduce nonnegative slack variables and rewrite the constraints in the form Ax + s = b. This is known as canonical form. The solution where x = 0 and s = b is a basic feasible solution and the corresponding basis matrix is the identity matrix.
To generalize better, you need to create an auxiliary problem. We again start by getting the problem into canonical form. Instead of introducing slack variables, we introduce what are referrred to as artificial variables in much the same way, such that Ax + y = b, and the solution where x = 0 and s = b is a basic feasible solution. However, it is not a basic feasible solution to the original problem, but to the new auxiliary problem where our cost function is the sum of the new artificial variables y which we introduced.
In order to move forward, the Simplex method is applied to the new auxiliary problem. The solution to the auxiliary problem will give you information on whether or not the problem is feasible. This is because if x is a feasible solution to the original problem, y = 0, and this should yield a zero cost solution to the original problem. If the cost of the solution to the auxiliary problem is nonzero, then the original problem is infeasible.
Once a problem has been determined as feasible, we move on to the second phase of the 2-phase method. But first, we must ensure no artificial variables remain in the initial basic feasible solution.
First the corresponding columns are removed from the tableau. If any artificial variables remain in the basis, they must also be removed. There are two ways to eliminate artificial variables that remain in the solution. First, if the corresponding row is all zeros, then it is linearly dependent on the other rows, it is unnecessary information and can simply be removed. If any value in the row is nonzero, then that value becomes the pivot value and the pivot routine from above is repeated.
Finally, we can replace the auxiliary cost function with our original cost function. The Simplex method is then applied to the resulting tableau, and when it concludes, we have found our optimal solution .
III. Simplex Implementation
The first step in implementing the Simplex method is to create your tableau. For my initial implementation, I included putting the problem into canonical form from standard form as part of the algorithm. So problems need only be entered in standard form. The example we use is from above, and the initial tableau constructed is seen below.
As above, the top row of the tableau consists of the current cost in the leftmost position, followed by the relative costs of introducing each basis variable.
Below the top row, you have your A matrix from the canonical form of the. problem. The description of generating this matrix is below.
First, we create t1, which is going to be the top row of our tableau. For book keeping reasons, the left most element will be a numpy NaN type and will not be used. The second element will be our current cost, which starts out as 0 before it is calculated. Then comes our initial cost function coefficients for the problem with slack variables. The coefficients for the original variables will be the coefficients in the original problem, while the slack variables will all have 0 coefficients.
Second, we need the indices of our initial basis. Since in this case our initial basis is our slack variables, the initialization is trivial. We take the transposed basis vector, and augment it with our transposed resource vector, original A matrix and basis matrix B. We now have the problem in canonical form. This is stored as t2.
t1 and t2 are then stacked on top of each other to get our starting tableau.
Finally, we convert all values to floating point in order to perform arithmetic later on.
def getTableau(self):
# construct starting tableau
if(self.minmax == "MIN" and self.transform == False):
self.c[0:len(c)] = -1 * self.c[0:len(c)]
self.transform = True
t1 = np.array([None, 0])
numVar = len(self.c)
numSlack = len(self.A)
t1 = np.hstack(([None], [0], self.c, [0] * numSlack))
basis = np.array([0] * numSlack)
for i in range(0, len(basis)):
basis[i] = numVar + i
A = self.A
if(not ((numSlack + numVar) == len(self.A[0]))):
B = np.identity(numSlack)
A = np.hstack((self.A, B))
t2 = np.hstack((np.transpose([basis]), np.transpose([self.b]), A))
tableau = np.vstack((t1, t2))
tableau = np.array(tableau, dtype ='float')
return tableauThe completed simplex method can be seen below, and an explanation follows
def optimize(self):
if(self.minmax == "MIN" and self.transform == False):
for i in range(len(self.c)):
self.c[i] = -1 * self.c[i]
transform = True
tableau = self.getTableau()
if(self.printIter == True):
print("Starting Tableau:")
self.printTableau(tableau)
# assume initial basis is not optimal
optimal = False# keep track of iterations for display
iter = 1while(True):
if(self.printIter == True):
print("----------------------------------")
print("Iteration :", iter)
self.printTableau(tableau)
if(self.minmax == "MAX"):
for profit in tableau[0, 2:]:
if profit > 0:
optimal = False
break
optimal = True
else:
for cost in tableau[0, 2:]:
if cost < 0:
optimal = False
break
optimal = True# if all directions result in decreased profit or increased cost
if optimal == True:
break
# nth variable enters basis, account for tableau indexing
if (self.minmax == "MAX"):
n = tableau[0, 2:].tolist().index(np.amax(tableau[0, 2:])) + 2
else:
n = tableau[0, 2:].tolist().index(np.amin(tableau[0, 2:])) + 2# minimum ratio test, rth variable leaves basis
minimum = 99999
r = -1for i in range(1, len(tableau)):
if(tableau[i, n] > 0):
val = tableau[i, 1]/tableau[i, n]
if val<minimum:
minimum = val
r = i
pivot = tableau[r, n]
print("Pivot Column:", n)
print("Pivot Row:", r)
print("Pivot Element: ", pivot)# perform row operations
# divide the pivot row with the pivot element
tableau[r, 1:] = tableau[r, 1:] / pivot# pivot other rows
for i in range(0, len(tableau)):
if i != r:
mult = tableau[i, n] / tableau[r, n]
tableau[i, 1:] = tableau[i, 1:] - mult * tableau[r, 1:]# new basic variable
tableau[r, 0] = n - 2
iter += 1
if(self.printIter == True):
print("----------------------------------")
print("Final Tableau reached in", iter, "iterations")
self.printTableau(tableau)
else:
print("Solved")
self.x = np.array([0] * len(c), dtype = float)
# save coefficients
for key in range(1, (len(tableau))):
if(tableau[key, 0] < len(c)):
self.x[int(tableau[key, 0])] = tableau[key, 1]
self.optimalValue = -1 * tableau[0,1]
You can see that we start by examining all of the relative costs (or profits in the case of a maximization problem). This allows us to determine which (if any) direction to move in. If there is no direction to move in, the algorithm notes that it is at the optimal solution and ends.
We save the best direction in which to move, and the corresponding variable that will enter the basis to move in this direction as n. This is our pivot column. We then must determine our pivot row. To do this, we perform a minimum ratio test. The minimum ratio is calculated by dividing the rsources in a row by its value in the pivot column. This is done for all rows, and whichever row results in the smallest value is the pivot row, and its associated variable leaves the basis in the next step.
To perform the pivot, the pivot row is first divided by the pivot element. The pivot element is the element whose index is both the pivot column and pivot row. This normalizes that value to 1.
Next, for every other row in the tableau, we perform another row operation and subtract the ratio of its value in the pivot column to the ratio of the pivot value multiplied by the pivot row.
Once this pivot operation is complete, the iteration of the algorithm is complete and it loops around again until optimal.
Finally, there is some misc support code in order to save the coefficients and solution for easy recall.
This method is great, in that it solves problems in standard form with an obvious initial basic feasible solution. There are a few problems that should be apparent though.
- What if the initial problem is infeasible? The solution above cannot handle infeasible problems.
- What if the initial problem is unbounded? This is not a problem with the Simplex method, but was not accounted for in the solution above.
- What if the initial basic feasible solution is not obvious? This method cannot handle this.
This is where the 2-phase Simplex method as stated above comes in. Luckily, the modifications to make a 2-phase solver are relatively straightforward if a Simplex solver is already on hand.
We start by creating a new tableau at the beginning of the problem. Since the problem is required to be in canonical form to begin with, we add artificial variables in order to create our auxiliary problem with which to find a basic feasible solution from.
def getTableauPhase1(self):
# construct starting tableau
numVar = len(self.c)
numArtificial = len(self.A)
t1 = np.hstack(([None], [0], [0] * numVar, [0] * numArtificial))
basis = np.array([0] * numArtificial)
for i in range(0, len(basis)):
basis[i] = numVar + i
A = self.A
if(not ((numArtificial + numVar) == len(self.A[0]))):
B = np.identity(numArtificial)
A = np.hstack((self.A, B))
t2 = np.hstack((np.transpose([basis]), np.transpose([self.b]), A))
tableau = np.vstack((t1, t2))
for i in range(1, len(tableau[0]) - numArtificial):
for j in range(1, len(tableau)):
if(self.minmax == "MIN"):
tableau[0, i] -= tableau[j, i]
else:
tableau[0, i] += tableau[j, i]
tableau = np.array(tableau, dtype ='float')
return tableauThen we break the Simplex routine out of the optimize() method into its own method, simplex(). The new method can be seen below.
def simplex(self, tableau):
if(self.printIter == True):
print("Starting Tableau:")
self.printTableau(tableau)
# assume initial basis is not optimal, problem is feasible, and problem is bounded
optimal = False
feasible = True
bounded = True # keep track of iterations for display
iter = 1 while(True):
if(self.printIter == True):
print("----------------------------------")
print("Iteration :", iter)
self.printTableau(tableau)
if(self.minmax == "MAX"):
for profit in tableau[0, 2:]:
if profit > 0:
optimal = False
break
optimal = True
else:
for cost in tableau[0, 2:]:
if cost < 0:
optimal = False
break
optimal = True # if all directions result in decreased profit or increased cost
if optimal == True:
break
# nth variable enters basis, account for tableau indexing
if (self.minmax == "MAX"):
n = tableau[0, 2:].tolist().index(np.amax(tableau[0, 2:])) + 2
else:
n = tableau[0, 2:].tolist().index(np.amin(tableau[0, 2:])) + 2 # minimum ratio test, rth variable leaves basis
minimum = 99999
r = -1for i in range(1, len(tableau)):
if(tableau[i, n] > 0):
val = tableau[i, 1]/tableau[i, n]
if val<minimum:
minimum = val
r = i
pivotElement = tableau[r, n]
if(self.printIter):
print("Pivot Column:", n)
print("Pivot Row:", r)
print("Pivot Element: ", pivotElement)
print("New Basis: " + str(tableau[r, 0]))for element in tableau[1:, n]:
if (element != 0):
self.bounded = True
if (not self.bounded):
print("Unbounded; No solution.")
return
tableau = self.pivot(tableau, r, n)# new basic variable
tableau[r, 0] = n - 2
iter += 1
if(self.printIter == True):
print("----------------------------------")
print("Final Tableau reached in", iter, "iterations")
self.printTableau(tableau)
else:
print("Solved")
return tableau
It can be seen that this version takes a tableau as an input, rather than creating its own tableau. This is because it will need to be called twice, with two different tableaus in the 2-phase algorithm. It now returns the final tableau for it to be reused as necessary by the caller method. It can also be seen that the pivot code has been removed and added to a method pivot(), which can be seen below. This is also because pivot is used in multiple points of the 2-phase Simplex. It requires a tableau, a pivot row index and a pivot column index as input.
def pivot(self, tableau, r, n):
pivot = tableau[r, n]
# perform row operations
# divide the pivot row with the pivot element
tableau[r, 1:] = tableau[r, 1:] / pivot
# pivot other rows
for i in range(0, len(tableau)):
if i != r:
mult = tableau[i, n] / tableau[r, n]
tableau[i, 1:] = tableau[i, 1:] - mult * tableau[r, 1:]
return tableauThe last element that needs to be added is a method for driving artificial variables out of the basis. This is the reason that the pivot method was best served as its own method. Code can be seen below.
def driveOutAV(self, tableau):
avPresent = False
avInd = []
#remove slack variables
tableau = np.delete(tableau, obj = np.s_[2 + len(self.c):], axis = 1)
#redundant row removal
for i in range(1, len(tableau)):
if tableau[i, 0] > len(self.c) - 1:
redundant = True
for j in range(2, len(self.c) + 2):
if (tableau[i, j] != 0):
redundant = False
avInd.append([i, False])
break
avInd.append([i, redundant])
if(len(avInd) == 0):
print("\nNo artificial variables to drive out\n")
return tableau
for key in avInd:
index, redundant = key
if(redundant):
tableau = np.delete(tableau, obj = index, axis = 0)
avInd.remove(key)
else:
n = -1
for i in range(2, len(tableau[0]) - 1):
if (tableau[r, i] != 0):
n = i
break
tableau = pivot(tableau, index, n)
return tableauThis method begins by first detecting if there are any artificial variables present in the tableau. If there are, it makes a note and adds them to the list of artificial variables. When it does this, it determines if the row for that variable is linearly independent of the other variables. If it is not, then it can simply be removed. Otherwise, we need to pivot around a non-zero element in that row in order to drive it out of the basis and bring one of the original variables back into the basis.
Once we have removed the artificial variables, we can continue with Simplex on the resulting tableau to obtain the solution to our original problem. The new optimize() method can be seen below. It can be seen that there is an additional method, getTableauPhase2(). This method simply replaces the objective function with that from the original problem.
def optimize(self):
#if(self.minmax == "MIN" and self.transform == False):
#for i in range(len(self.c)):
#self.c[i] = -1 * self.c[i]
#transform = True
tableau = self.getTableauPhase1()
print("Phase 1:\n")
tableau = self.simplex(tableau)
if(not self.bounded):
return
if(tableau[0, 1] != 0):
self.feasible = False
print("Problem Infeasible; No Solution")
return
print("Phase 2:\n")
tableau = self.driveOutAV(tableau)
tableau = self.getTableauPhase2(tableau)
tableau = self.simplex(tableau)
if(not self.bounded):
return
self.x = np.array([0] * len(c), dtype = float)
# save coefficients
for key in range(1, (len(tableau))):
if(tableau[key, 0] < len(c)):
self.x[int(tableau[key, 0])] = tableau[key, 1]
self.optimalValue = -1 * tableau[0,1]Lastly, there is the wrapper code that makes this slightly easier to use and a bit more similar to commercial solver packages. You can also select whether or not to display iterations of the Simplex method, or to just display the solution.
class LinearModel:
def __init__(self, A = np.empty([0,0]), b = np.empty([0,0]), c = np.empty([0,0]), minmax = "MAX"):
self.A = A
self.b = b
self.c = c
self.x = [float(0)] * len(c)
self.minmax = minmax
self.printIter = True
self.optimalValue = None
self.transform = False
self.feasible = True
self.bounded = False
def addA(self, A):
self.A = A
def addB(self, b):
self.b = b
def addC(self, c):
self.c = c
self.transform = False
def setObj(self, minmax):
if(minmax == "MIN" or minmax == "MAX"):
self.minmax = minmax
else:
print("Invalid objective.")
self.transform = False
def setPrintIter(self, printIter):
self.printIter = printIter
def printSoln(self):
if(self.feasible):
if(self.bounded):
print("Coefficients: ")
print(self.x)
print("Optimal value: ")
print(self.optimalValue)
else:
print("Problem Unbounded; No Solution")
else:
print("Problem Infeasible; No Solution")
def printTableau(self, tableau):
print("ind \t\t", end = "")
for j in range(0, len(c)):
print("x_" + str(j), end = "\t")
for j in range(0, (len(tableau[0]) - len(c) - 2)):
print("s_" + str(j), end = "\t")
print()
for j in range(0, len(tableau)):
for i in range(0, len(tableau[0])):
if(not np.isnan(tableau[j, i])):
if(i == 0):
print(int(tableau[j, i]), end = "\t")
else:
print(round(tableau[j, i], 2), end = "\t")
else:
print(end = "\t")
print()IV. Testing Implementation
This section will cover some basic tests of the implementation shown previously. It primarily consists of using test problems and comparing the solutions to those obtained from a well known commercial software, Gurobi. Gurobi does not implement Simplex, and instead uses a faster technique called the interior point method. However, that does not change the solutions to these problems.
We will first test the original Simplex implementation.
We begin with the original example as shown above.
The first implementation that I did will put this into canonical form for us assuming the original problem is in standard form, requiring that we enter the problem as below.
model1 = LinearModel()A = np.array([[1, 2, 2],
[2, 1, 2],
[2, 2, 1]])
b = np.array([20, 20, 20])
c = np.array([-10, -12, -12])model1.addA(A)
model1.addB(b)
model1.addC(c)
model1.setObj("MIN")print("A =\n", A, "\n")
print("b =\n", b, "\n")
print("c =\n", c, "\n\n")
model1.optimize()
print("\n")
model1.printSoln()
Executing this will give us the output below.
Compared to the Gurobi code below
m = Model("Example1")x1 = m.addVar(vtype = GRB.CONTINUOUS, name = 'x1')
x2 = m.addVar(vtype = GRB.CONTINUOUS, name = 'x2')
x3 = m.addVar(vtype = GRB.CONTINUOUS, name = 'x3')m.setObjective(-10 * x1 - 12 * x2 - 12 * x3, GRB.MINIMIZE)m.addConstr(1 * x1 + 2 * x2 + 2 * x3 <= 20, 'c0')m.addConstr(2 * x1 + 1 * x2 + 2 * x3 <= 20, 'c1')m.addConstr(2 * x1 + 2 * x2 + 1 * x3 <= 20, 'c2')m.optimize()for v in m.getVars():
print('%s %g' % (v.varName, v.x))
print('Obj: %g' % m.objVal)
Which gives us the below output when executed. It can be seen that the solutions are identical.
Testing the 2-phase Simplex is done in a similar manner. However, for the 2-phase model, input is required in canonical form rather than standard form. Luckily, reference 1 has an example problem that is already in canonical form as below.
This is entered as below.
model1 = LinearModel()A = np.array([[1, 2, 3, 0],
[-1, 2, 6, 0],
[0, 4, 9, 0],
[0, 0, 3, 1]])
b = np.array([3, 2, 5, 1])
c = np.array([1, 1, 1, 0])model1.addA(A)
model1.addB(b)
model1.addC(c)
model1.setObj("MIN")
model1.setPrintIter(True)print("A =\n", A, "\n")
print("b =\n", b, "\n")
print("c =\n", c, "\n\n")
model1.optimize()
print("\n")
model1.printSoln()
Which gives the below output.
Which can then be compared to the Gurobi program below, and it’s output below that.
m1 = Model("Example2")x1 = m1.addVar(vtype = GRB.CONTINUOUS, name = 'x1')
x2 = m1.addVar(vtype = GRB.CONTINUOUS, name = 'x2')
x3 = m1.addVar(vtype = GRB.CONTINUOUS, name = 'x3')
x4 = m1.addVar(vtype = GRB.CONTINUOUS, name = 'x4')m1.setObjective(x1 + x2 + x3, GRB.MAXIMIZE)m1.addConstr(1 * x1 + 2 * x2 + 2 * x3 == 3, 'c0')m1.addConstr(-1 * x1 + 2 * x2 + 6 * x3 == 2, 'c1')m1.addConstr(4 * x2 + 9 * x3 == 5, 'c2')m1.addConstr(3 * x3 + 1 * x4 == 1, 'c3')m1.optimize()for v in m1.getVars():
print('%s %g' % (v.varName, v.x))
print('Obj: %g' % m1.objVal)
Again, results are identical.
V. Takeaways
Despite not being much more complicated than the naive method, the Simplex method offers a significant increase in performance. On my MacBook, the above problems took a fraction of a second. It is also deterministic, and easily traceable should anything go wrong. The 2-phase implementation above also determines if a problem is bounded or infeasible, something the naive approach will not pick up on.
The implementation requires very little math other than some elementary linear algebra, making it simple to understand for those new to linear programming or those with limited math background.
The last takeaway is how difficult it is to create something as user friendly and simple to use as commercial software like Gurobi is. For one, you can enter a problem in a much more intuitive way, exactly as it is stated, rather than wrangling it into some particular form. The amount of background code needed to do this is probably more significant than implementing the Simplex algorithm. And secondly, this implementation is not particularly efficient. I did not test it with any large scale problems, but I suspect that it would quickly become unusable as the number of variables increase.
VI. References and Links
- D. Bertsimas and J. N. Tsitsiklis, Introduction to Linear Optimization. Belmont, MA: Athena Scientific, 1997.
All code featured is hosted in the below Github repository: