Last active
April 6, 2018 20:05
-
-
Save thien/3df4fb97c7383244d5ccc7fd56b3b8bf to your computer and use it in GitHub Desktop.
Simplex Algorithm implementation in Python (3.x)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
def simplex(c, A, b, basis): | |
z = 0 | |
tableau = None | |
foundOptimal = False | |
numRound = 1 | |
# we'll be using tableau because it's easier and makes more sense. | |
# using the basis, check if A is in canonical form. | |
# i.e: the basis variables used create the identity matrix. | |
if isCanonicalForm(A,basis): | |
# convert c, A, and b to tableau form. | |
tableau = componmentsToTableau(-c,A,b) | |
while not foundOptimal: | |
print("Round:", numRound) | |
numRound += 1 | |
# to get the entering variable, look @ the negative entries in | |
# the top row (which is c). Use Blands rule to get the entering | |
# variable. | |
enteringVariable = getEnteringVariable(tableau) | |
# now we know the column to go to, we need to look into what | |
# row to pivot from. | |
pivotLocation = findMin(A, enteringVariable, b) | |
if pivotLocation == -1: | |
print("The solution is unbounded..") | |
return -1 | |
print("Pivot Location [regards to A]: (" + str(enteringVariable+1)+","+str(pivotLocation+1)+")") | |
# now we can pivot. | |
tableau = pivot(tableau,pivotLocation+2,enteringVariable+2) | |
print("new Tableau:\n",tableau, "\n") | |
# extract the components of the equation | |
(A,b,c,z) = tableauToComponments(tableau) | |
# find the leaving column. We can tell which one is the | |
# leaving column by looking at the columns that make up A | |
# and removing ones that don't make the basis. | |
leavingVariable, basis = findLeavingColumn(A,basis) | |
print("Entering, Leaving Columns (for Bases):", enteringVariable+1, leavingVariable) | |
# at this stage we've just finished a round of simplex. | |
# check if C is all positive. if they are, then we've found the optimal | |
# solution. | |
if min(c) >= 0: foundOptimal = True | |
print("---------------------------------Best: ", z) | |
print() | |
print("Found the (optimal?) solution:") | |
print("z = ", z) | |
print("Where x = ", produceX(A,basis,b)) | |
return False | |
def produceX(A,basis,b): | |
# use the basis to produce the columns in A | |
# iterate through the basis (it starts from 1 so -1 to them) | |
basis = [x-1 for x in basis] | |
# find how many items go in x. | |
numberOfXValues = A.shape[1] | |
result = [] | |
# iterate through the indexes | |
for x in range(numberOfXValues): | |
if x in basis: | |
pos = A[:,x].T.tolist()[0].index(1) | |
# add b's value to x | |
result.append(b.tolist()[0][pos]) | |
else: | |
# fill it with a zero. | |
result.append(0.0) | |
return result | |
def getEnteringVariable(t): | |
# repull c. | |
c = t[0].tolist()[0][1:-1] | |
# using bland's rule (Choose the entering variable k with | |
# the SMALLEST INDEX) and the negative entries, get the | |
# first one that pops up. | |
for i in range(len(c)): | |
if c[i] < 0: | |
return i | |
# otherwise we can't find a solution. | |
return -1 | |
def isCanonicalForm(A,basis): | |
""" | |
Check if the basis variables actually create | |
the identity matrix. | |
""" | |
# get the basic variables X_b | |
X_b = [] | |
for b in basis: | |
X_b.append(np.asarray(np.transpose(A[:,b-1])).tolist()[0]) | |
# convert X_b into a numpy matrix | |
X_b = np.matrix(X_b) | |
# get the size of the basis to create an identity matrix. | |
dim = X_b.shape[0] | |
I = np.eye(dim) | |
# check if they're the same. | |
if np.array_equal(X_b, I): | |
return True | |
else: | |
return False | |
def findMin(A, ent, b): | |
# convert those subsections into lists | |
# so we can traverse them easier. | |
a_col = A[:,ent].T.tolist()[0] | |
b_cop = b.copy().tolist()[0] | |
negatives = [] | |
contenders = [] | |
for i in range(len(a)): | |
division = b_cop[i]/a_col[i] | |
# if it's a negative value then we mess up the number | |
# (we want the smallest positive) | |
if division < 0: | |
division = 2036854775807 | |
negatives.append(division) | |
contenders.append(division) | |
# if they're all negative then our problem is unbounded. | |
if len(negatives) == len(contenders): | |
return -1 | |
# now we want the index of the smallest positive number. | |
index = contenders.index(min(contenders)) | |
return index | |
def componmentsToTableau(c,A,b,z=None): | |
""" | |
Converts the matrix into an augmented matrix | |
based on the linear system z-c^(T)x = z' | |
""" | |
tableau = None | |
# create top | |
top = [1.0]+np.transpose(c).tolist()[0]+[0.0] | |
top = np.matrix(top) | |
# create zeros | |
lhs_zero = np.array([np.zeros(A.shape[0])]).T | |
# transpose b | |
b_t = b.T | |
# initiate bottom variable | |
bottom = np.asarray(A.copy()) | |
# concatenate those zeros and bottom | |
bottom = np.hstack((lhs_zero, bottom)) | |
# concatenate bottom and b | |
bottom = np.hstack((bottom,b_t)) | |
# now we concatenate top and bottom | |
tableau = np.vstack((top, bottom)) | |
return tableau | |
def tableauToComponments(tableau): | |
""" | |
Converts the tableau to regular components. | |
""" | |
# let's create A | |
A = tableau[1:] | |
A = np.delete(A, 0, 1) | |
A = np.delete(A, A.shape[1]-1, 1) | |
# creating b | |
b = tableau[:,-1][1:].T | |
# creating c | |
c = tableau[0].tolist()[0][1:-1] | |
# creating z | |
z = tableau[0,-1] | |
return (A,b,c,z) | |
def findLeavingColumn(A, basis): | |
""" | |
We'll use bland's other rule here | |
i.e we get the smallest index first; or rather | |
the first thing we find that isn't to scratch will suffice. | |
""" | |
# iterate through the columns and see if we can | |
# find the identity matrix. | |
new_basis = [] | |
for i in range(A.shape[1]): | |
# get the column | |
column = A[:,i].T.tolist()[0] | |
count_zero = 0 | |
int_is_one = False | |
for j in column: | |
if j == 0: count_zero += 1 | |
else: | |
if j == 1: int_is_one = True | |
if int_is_one and count_zero == len(column)-1: | |
new_basis.append(i+1) | |
# create placeholder variable for our leaving | |
# column | |
odd_one_out = -1 | |
for i in basis: | |
# this oneliner finds the first leaving value. | |
if i not in new_basis: odd_one_out = i; break | |
return (odd_one_out, new_basis) | |
def pivot(a, i, j): | |
""" | |
Remember with elementary row operations you can do the following: | |
1. Multiply a row with a constant. | |
2. Add a row * constant to another row. | |
3. Swap the positions of two rows. | |
Pivoting on an entry T_(i,j) != 0 means transforming the matrix T into | |
a matrix T' by performing the following operations: | |
- multiply the ith row of T by 1/T_(i,j) (to make T_(i,j) == 1) | |
- For each k != i, add the ith row * constant to the kth row. (to | |
make T'_(k,j) = 0) | |
""" | |
# make sure we simplify the code s.t i and j refer to np indexes. | |
i, j = i - 1, j - 1 | |
# make a copy of the matrix we're going to manipulate | |
T = a.copy() | |
# get the pivot | |
pivot = T[i,j] | |
# 1. multiply the ith row of T by 1/T_(i,j) (to make T_(i,j) == 1) | |
divider = 1/pivot | |
T[i] = np.multiply(T[i], divider) | |
# 2. For each k != i, add the ith row * -constant to the kth row. | |
# (to make T'_(k,j) = 0) | |
for k in range(T.shape[0]): | |
if k != i: | |
# find the negative constant | |
constant = -T[k,j] | |
multiplier = np.multiply(T[i],constant) | |
T[k] = np.add(T[k], multiplier) | |
return T | |
if __name__ == "__main__": | |
# Bounded Example: | |
c = np.matrix([[2],[3],[0],[0],[0]], dtype=float) | |
a = np.matrix([ | |
[1,1,1,0,0], | |
[2,1,0,1,0], | |
[-1,1,0,0,1] | |
], dtype=float) | |
b = np.matrix([6,10,4]) | |
basis = [3,4,5] | |
# ----------- | |
# Unbounded Example: | |
# c = np.matrix([[-1],[3],[0],[0],[1]],dtype=float) | |
# a = np.matrix([ | |
# [-2,4,1,0,1], | |
# [-3,7,0,1,1] | |
# ], dtype=float) | |
# b = np.matrix([1,3]) | |
# basis = [3,4] | |
simplex(c, a, b, basis) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I've verified the program to detect unbounded examples; that being said I've only tested it on two LPs. I'll update it in the event that other LPs break it.