Updated date:

Simplex Algorithm in Python

Author:

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

Ravi on April 30, 2020:

can You suggest alternate method for the same algo where we don't need to use tableau

Klaus Scheicher on January 17, 2018:

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 on July 12, 2017:

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 on February 24, 2017:

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

sasan on November 24, 2016:

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

Carlos on May 08, 2016:

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

TomTranter on October 14, 2014:

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 on July 04, 2013:

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

moshahmed on January 02, 2013:

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

taw9 (author) on December 01, 2012:

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

-taw

Karapiperis Al on December 01, 2012:

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 ?