Skip to content

Instantly share code, notes, and snippets.

@hdgarrood hdgarrood/simplex.py
Last active May 16, 2017

Embed
What would you like to do?
An extremely naive implementation of the simplex algorithm.
from collections import namedtuple
import numpy as np
# This is an implementation of the version of the simplex algorithm we are
# expected to execute by hand in LPMS. There is just one important difference:
# because python and numpy use zero-based indexing, this does too. So for
# example you need to subtract 1 from each element of your basic and nonbasic
# sets.
LPState = namedtuple("LPState", ["basic", "nonbasic"])
LPSolution = namedtuple("LPSolution", ["basic_variables", "basic_variable_values", "objective"])
class LP:
def __init__(self, costs, A, b):
"""
Constructs a linear programming problem of the form
maximize f = (costs)^T . x
subject to A . x <= b
x >= 0
with x being the vector of decision variables.
"""
self._check_dims(costs, A, b)
self.costs = costs
self.A = A
self.b = b
# Computed properties
self.num_constraints = A.shape[0]
self.num_variables = A.shape[1]
self.Abar = np.concatenate([A, np.identity(self.num_constraints)], axis=1)
self.full_costs = np.concatenate([self.costs, np.repeat(0, self.num_constraints)])
def _check_dims(self, costs, A, b):
if len(costs.shape) != 1:
raise ValueError("Expected costs to be a vector")
if len(A.shape) != 2:
raise ValueError("Expected A to be a matrix")
if len(b.shape) != 1:
raise ValueError("Expected b to be a vector")
if A.shape[1] != costs.shape[0]:
raise ValueError("Size mismatch between A and costs")
if A.shape[0] != b.shape[0]:
raise ValueError("Size mismatch between A and b")
def __repr__(self):
def show_constraint(j):
return " " + " + ".join([
repr(coeff) + "*x_" + repr(i) for (i, coeff) in enumerate(self.A[j])
]) + " <= " + repr(self.b[j])
objective_fn = " + ".join([
repr(c) + "*x_" + repr(i) for (i, c) in enumerate(self.costs)
])
constraints = [show_constraint(i) for i in range(0, self.num_constraints)]
return "\n".join([
"maximize f = " + objective_fn,
"subject to:"] +
constraints +
[" x_i >= 0"]
)
def initial_state(self):
# TODO: ensure this is feasible
m = self.num_constraints
n = self.num_variables
return LPState(
basic=range(n, n+m),
nonbasic=range(0, n)
)
def dual(self):
return LP(-self.b,
-np.transpose(self.A),
-self.costs
)
eg2 = LP(
# costs
np.array([1, 2, 3]),
# A
np.array([[1, 0, 2],
[0, 1, 2]]),
# b
np.array([3, 2])
)
eg2_state = LPState(
basic=[0,2],
nonbasic=[1,3,4]
)
assignment2 = LP(
# costs
np.array([1, 10]),
# A
np.array([[1, 20],
[0, 1 ]]),
# b
np.array([100, 1])
)
def iterate(lp, state, verbose=True):
"""Returns an LPSolution if the given state comprises a feasible optimal
basis, a new LPState if it is feasible but non-optimal, and throws an error
if the state is non-feasible or the problem is unbounded."""
def info(msg):
if verbose:
print(msg)
B = lp.Abar[:,state.basic]
N = lp.Abar[:,state.nonbasic]
info("B = {}".format(B))
info("N = {}".format(N))
bhat = np.linalg.solve(B, lp.b)
info("bhat = {}".format(bhat))
if not is_feasible(bhat):
raise ValueError("unfeasible basis; bhat = {}".format(bhat))
basic_costs = lp.full_costs[state.basic]
nonbasic_costs = lp.full_costs[state.nonbasic]
info("basic_costs = {}".format(basic_costs))
info("nonbasic_costs = {}".format(nonbasic_costs))
pi = np.linalg.solve(B.T, basic_costs)
info("pi = {}".format(pi))
reduced_costs = nonbasic_costs - np.matmul(N.T, pi)
info("reduced_costs = {}".format(reduced_costs))
if is_optimal(reduced_costs):
info("Optimal!")
return LPSolution(
basic_variables=state.basic,
basic_variable_values=bhat,
objective=np.dot(basic_costs, bhat)
)
q = np.argmax(reduced_costs)
qp = state.nonbasic[q]
# q-th column of N
a_q = N[:,q]
info("q = {}".format(q))
info("qp = {}".format(qp))
info("a_q = {}".format(a_q))
a_q_hat = np.linalg.solve(B, a_q)
info("a_q_hat = {}".format(a_q_hat))
if is_unbounded(a_q_hat):
raise ValueError("Unbounded!")
possible_ps = np.where(a_q_hat > 0, bhat / a_q_hat, np.infty)
info("possible_ps = {}".format(possible_ps))
p = np.argmin(possible_ps)
pp = state.basic[p]
info("p = {}".format(p))
info("pp = {}".format(pp))
new_nonbasic = list(state.nonbasic)
new_nonbasic[q] = pp
new_basic = list(state.basic)
new_basic[p] = qp
return LPState(
nonbasic=new_nonbasic,
basic=new_basic
)
def is_feasible(bhat):
return (bhat > 0).all()
def is_optimal(reduced_costs):
return (reduced_costs <= 0).all()
def is_unbounded(a_q_hat):
return (a_q_hat <= 0).all()
def main():
"This runs through 1.a.ii of May 2016 paper."
may2016_q1 = LP(
# costs
np.array([3, 1, 2]),
# A
np.array([[1, -2, 4],
[2, -3, -1]]),
# b
np.array([4, 12])
)
initial_state = LPState(
basic=[0,4],
nonbasic=[1,2,3]
)
new_state = iterate(may2016_q1, initial_state)
iterate(may2016_q1, new_state)
if __name__ == '__main__':
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.