-
Notifications
You must be signed in to change notification settings - Fork 0
/
simplex.py
128 lines (101 loc) · 4.29 KB
/
simplex.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
import numpy as np
def make_tableau(c, b, A):
n = len(c)
m = len(b)
simplex_tableau = np.hstack((A, np.identity(m), np.zeros((m, 1)), np.array([b]).T))
simplex_tableau = np.vstack((simplex_tableau, np.array([np.hstack((-c, np.zeros(m), np.ones(1), np.zeros(1)))])))
return simplex_tableau
def solve(c, b, A, debug_output=False, blands_rule=False):
n = len(c)
m = len(b)
simplex_tableau = make_tableau(c, b, A)
return solve_tableau(n, m, simplex_tableau, debug_output, blands_rule)
def solve_tableau(n, m, simplex_tableau, debug_output=False, blands_rule=False):
if debug_output:
print("Simplex tableau :")
print(simplex_tableau)
while True:
# Column selection
# Equivalent to searching for the variable in the base that would increase the output function the most
optimal = True
column = 0
if blands_rule:
for i, e in enumerate(simplex_tableau[-1]):
if e < 0:
optimal = False
column = i
break
else:
for i, e in enumerate(simplex_tableau[-1]):
if e < 0:
optimal = False
if e < simplex_tableau[-1, column]:
column = i
if optimal:
if debug_output:
print("Optimal simplex tableau found")
break
# Row selection
# Equivalent to searching for the variable outside the base that constraints the augmentation of the output function the most (smallest b_i / leaving variable coefficient)
row = -1
for i in range(m):
# Only positive entries are considered
if simplex_tableau[i, column] <= 0:
continue
if row == -1 or simplex_tableau[row, -1] / simplex_tableau[row, column] > simplex_tableau[i, -1] / simplex_tableau[i, column]:
row = i
if row == -1:
raise Exception(f"Unbounded solutions in tableau (selected column : {column})", simplex_tableau)
pivot = simplex_tableau[row, column]
simplex_tableau[row] /= pivot
for r in range(len(simplex_tableau)):
if r != row:
simplex_tableau[r] -= simplex_tableau[row] * simplex_tableau[r, column]
if debug_output:
print("Current simplex tableau :")
print(simplex_tableau)
return simplex_tableau
# Returns an array of basic variables (row, column, value)
def get_basic_variables(simplex_tableau):
x = []
# We are excluding the last two columns as they do not contain variables
for i in range(len(simplex_tableau.T) - 2):
nonzeros = 0
nonones = 0
oneidx = -1
for cidx, c in enumerate(simplex_tableau.T[i]):
if c != 0:
nonzeros += 1
if c != 1:
nonones += 1
if c == 1:
oneidx = cidx
if nonzeros == 1 and nonones == len(simplex_tableau.T[i]) - 1:
# This is a basic variable
x.append((oneidx, i, simplex_tableau[oneidx, -1]))
return x
def add_constraint(simplex_tableau, constraint_row, constraint_b_value):
rows, columns = simplex_tableau.shape
constraint_row = np.concatenate((constraint_row, [1.0, 0.0, constraint_b_value]))
simplex_tableau = np.insert(simplex_tableau, -2, np.zeros(rows), axis=1)
simplex_tableau = np.insert(simplex_tableau, -1, constraint_row, axis=0)
slack_variable = np.zeros(rows)
slack_variable[-2] = 1
slack_variable = np.array([slack_variable])
return simplex_tableau
if __name__ == '__main__':
print("Simplex algorithm")
print("Maximize <c, x> subject to constraints Ax <= b, x_i >= 0, where A is an m*n matrix")
print("==================")
n = int(input("Size n (size of vector x) : "))
m = int(input("Size m (amount of constraints) : "))
print("c vector :")
c = np.array([float(input("c_" + str(i) + " : ")) for i in range(n)])
print("b vector :")
b = np.array([float(input("b_" + str(i) + " : ")) for i in range(m)])
print("A matrix :")
A = []
for i in range(m):
A.append([float(input("A_" + str(i) + "," + str(j) + " : ")) for j in range(n)])
A = np.array(A)
solve(c, b, A, True)