# Simplex Algorithm in Python

Updated on July 24, 2012 ## Introduction

The Simplex algorithm is an optimization procedure for linear programs. As implied by "linear", the objective function for such a problem is a linear combination of the decision variables. Additionally, the region of possible solutions (aka "feasible region") is a convex polyhedron.

Considered one of the most important algorithms of all time, I've always found the existing open source implementations rather frustrating to use. When I was in graduate school, we had to use AMPL, a language specifically for modeling optimization problems (which was horrible and taught me little to nothing).

Later on, through my work, I found myself again dealing with linear optimization problems. By this time, however, I was armed with a relatively strong skillset in Python and could avoid using such a domain-specific language like AMPL. Below is a relatively straightforward implementation of the famous algorithm. As an added bonus to those you who might use this for homework problems, you can print the tableau after each pivot.

## Simplex in Python

```from __future__ import division
from numpy import *

class Tableau:

def __init__(self, obj):
self.obj =  + obj
self.rows = []
self.cons = []

def add_constraint(self, expression, value):
self.rows.append( + expression)
self.cons.append(value)

def _pivot_column(self):
low = 0
idx = 0
for i in range(1, len(self.obj)-1):
if self.obj[i] < low:
low = self.obj[i]
idx = i
if idx == 0: return -1
return idx

def _pivot_row(self, col):
rhs = [self.rows[i][-1] for i in range(len(self.rows))]
lhs = [self.rows[i][col] for i in range(len(self.rows))]
ratio = []
for i in range(len(rhs)):
if lhs[i] == 0:
ratio.append(99999999 * abs(max(rhs)))
continue
ratio.append(rhs[i]/lhs[i])
return argmin(ratio)

def display(self):
print '\n', matrix([self.obj] + self.rows)

def _pivot(self, row, col):
e = self.rows[row][col]
self.rows[row] /= e
for r in range(len(self.rows)):
if r == row: continue
self.rows[r] = self.rows[r] - self.rows[r][col]*self.rows[row]
self.obj = self.obj - self.obj[col]*self.rows[row]

def _check(self):
if min(self.obj[1:-1]) >= 0: return 1
return 0

def solve(self):

# build full tableau
for i in range(len(self.rows)):
self.obj += 
ident = [0 for r in range(len(self.rows))]
ident[i] = 1
self.rows[i] += ident + [self.cons[i]]
self.rows[i] = array(self.rows[i], dtype=float)
self.obj = array(self.obj + , dtype=float)

# solve
self.display()
while not self._check():
c = self._pivot_column()
r = self._pivot_row(c)
self._pivot(r,c)
print '\npivot column: %s\npivot row: %s'%(c+1,r+2)
self.display()

if __name__ == '__main__':

"""
max z = 2x + 3y + 2z
st
2x + y + z <= 4
x + 2y + z <= 7
z          <= 5
x,y,z >= 0
"""

t = Tableau([-2,-3,-2])
t.add_constraint([2, 1, 1], 4)
t.add_constraint([1, 2, 1], 7)
t.add_constraint([0, 0, 1], 5)
t.solve()

```

## Related

0 of 8192 characters used
• Klaus Scheicher

20 months ago

I think your code is missing the first phase in order to find a feasible initial tableau. Your code does not check if there exists a non set of empty solutions to the constraints.

• Jonathan

2 years ago

Thank you for the code.

In the case of rational/integer only linear programs, you could use sympy.Matrix and sympy.Array, which support fractions. This eliminates numerical error.

• Rohith

2 years ago

Thank you very much. This code, which uses classes has helped my to simplify my code.

• sasan

2 years ago

Can you tell me what is array? it is used without definition and python defines that as an error!

• Carlos

3 years ago

Thank you so much for sharing your knowledge, your implementation is very simple and explicit. :)

• TomTranter

4 years ago

Hi,

Thanks for posting this. I've been trying to figure out a problem using linear programming in python. It's possible in matlab using linprog as demonstrated here. http://www.mathworks.co.uk/matlabcentral/fileexcha...

I'm generating a random polyhedra want to know the largest inscribed sphere. Working out the normals of the faces is fine in python but formulating the problem to solve with linear programming is where I'm getting stuck. There are some optimisation functions in scipy but your function is the closest I've seen to something suitable and quick. Do you think it's possible to replicate the matlab code and solve using your function?

Thanks

• Maxim

6 years ago

Sooo .... how do I use it? I should be able to somehow specify a function, and starting parameters?

• moshahmed

6 years ago

Useful and thanks. Bug fix: if 'lhs[i] == 0' should be 'lhs[i] less than equal 1e-5'

• AUTHOR

taw9

6 years ago

You can get Numpy from http://numpy.scipy.org/

-taw

• Karapiperis Al

6 years ago

i copied this whole code in a python 2.7.3 version but there is no library/module with the name numpy. please if numpy is a module of your own can you post it ?

working