- HubPages
*»* - Technology
*»* - Computers & Software
*»* - Computer Science & Programming
*»* - Programming Languages

# Simplex Algorithm in Python

## 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 = [1] + obj self.rows = [] self.cons = [] def add_constraint(self, expression, value): self.rows.append([0] + 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 += [0] 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 + [0], 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()

## Comments

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.

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.

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

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

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

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

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

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

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 ?

10