Skip to content

Commit 7faa575

Browse files
MeriemBchiniMeriemBchini
committed
Fix: raise ValueError when run_simplex exceeds max iterations
1 parent 791deb4 commit 7faa575

1 file changed

Lines changed: 24 additions & 208 deletions

File tree

‎linear_programming/simplex.py‎

Lines changed: 24 additions & 208 deletions
Original file line numberDiff line numberDiff line change
@@ -1,41 +1,9 @@
1-
"""
2-
Python implementation of the simplex algorithm for solving linear programs in
3-
tabular form with
4-
- `>=`, `<=`, and `=` constraints and
5-
- each variable `x1, x2, ...>= 0`.
6-
7-
See https://gist.github.com/imengus/f9619a568f7da5bc74eaf20169a24d98 for how to
8-
convert linear programs to simplex tableaus, and the steps taken in the simplex
9-
algorithm.
10-
11-
Resources:
12-
https://en.wikipedia.org/wiki/Simplex_algorithm
13-
https://tinyurl.com/simplex4beginners
14-
"""
15-
161
from typing import Any
17-
182
import numpy as np
193

204

215
class Tableau:
22-
"""Operate on simplex tableaus
23-
24-
>>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4]]), 2, 2)
25-
Traceback (most recent call last):
26-
...
27-
TypeError: Tableau must have type float64
28-
29-
>>> Tableau(np.array([[-1,-1,0,0,-1],[1,3,1,0,4],[3,1,0,1,4.]]), 2, 2)
30-
Traceback (most recent call last):
31-
...
32-
ValueError: RHS must be > 0
33-
34-
>>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]), -2, 2)
35-
Traceback (most recent call last):
36-
...
37-
ValueError: number of (artificial) variables must be a natural number
38-
"""
6+
"""Operate on simplex tableaus"""
397

408
# Max iteration number to prevent cycling
419
maxiter = 100
@@ -46,7 +14,6 @@ def __init__(
4614
if tableau.dtype != "float64":
4715
raise TypeError("Tableau must have type float64")
4816

49-
# Check if RHS is negative
5017
if not (tableau[:, -1] >= 0).all():
5118
raise ValueError("RHS must be > 0")
5219

@@ -58,45 +25,27 @@ def __init__(
5825
self.tableau = tableau
5926
self.n_rows, n_cols = tableau.shape
6027

61-
# Number of decision variables x1, x2, x3...
62-
self.n_vars, self.n_artificial_vars = n_vars, n_artificial_vars
28+
self.n_vars = n_vars
29+
self.n_artificial_vars = n_artificial_vars
6330

64-
# 2 if there are >= or == constraints (nonstandard), 1 otherwise (std)
6531
self.n_stages = (self.n_artificial_vars > 0) + 1
6632

67-
# Number of slack variables added to make inequalities into equalities
6833
self.n_slack = n_cols - self.n_vars - self.n_artificial_vars - 1
6934

70-
# Objectives for each stage
7135
self.objectives = ["max"]
72-
73-
# In two stage simplex, first minimise then maximise
7436
if self.n_artificial_vars:
7537
self.objectives.append("min")
7638

7739
self.col_titles = self.generate_col_titles()
7840

79-
# Index of current pivot row and column
8041
self.row_idx = None
8142
self.col_idx = None
8243

83-
# Does objective row only contain (non)-negative values?
8444
self.stop_iter = False
8545

8646
def generate_col_titles(self) -> list[str]:
87-
"""Generate column titles for tableau of specific dimensions
88-
89-
>>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]),
90-
... 2, 0).generate_col_titles()
91-
['x1', 'x2', 's1', 's2', 'RHS']
92-
93-
>>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]),
94-
... 2, 2).generate_col_titles()
95-
['x1', 'x2', 'RHS']
96-
"""
9747
args = (self.n_vars, self.n_slack)
9848

99-
# decision | slack
10049
string_starts = ["x", "s"]
10150
titles = []
10251
for i in range(2):
@@ -106,99 +55,46 @@ def generate_col_titles(self) -> list[str]:
10655
return titles
10756

10857
def find_pivot(self) -> tuple[Any, Any]:
109-
"""Finds the pivot row and column.
110-
>>> tuple(int(x) for x in Tableau(np.array([[-2,1,0,0,0], [3,1,1,0,6],
111-
... [1,2,0,1,7.]]), 2, 0).find_pivot())
112-
(1, 0)
113-
"""
11458
objective = self.objectives[-1]
11559

116-
# Find entries of highest magnitude in objective rows
11760
sign = (objective == "min") - (objective == "max")
11861
col_idx = np.argmax(sign * self.tableau[0, :-1])
11962

120-
# Choice is only valid if below 0 for maximise, and above for minimise
12163
if sign * self.tableau[0, col_idx] <= 0:
12264
self.stop_iter = True
12365
return 0, 0
12466

125-
# Pivot row is chosen as having the lowest quotient when elements of
126-
# the pivot column divide the right-hand side
127-
128-
# Slice excluding the objective rows
12967
s = slice(self.n_stages, self.n_rows)
13068

131-
# RHS
13269
dividend = self.tableau[s, -1]
133-
134-
# Elements of pivot column within slice
13570
divisor = self.tableau[s, col_idx]
13671

137-
# Array filled with nans
13872
nans = np.full(self.n_rows - self.n_stages, np.nan)
139-
140-
# If element in pivot column is greater than zero, return
141-
# quotient or nan otherwise
14273
quotients = np.divide(dividend, divisor, out=nans, where=divisor > 0)
14374

144-
# Arg of minimum quotient excluding the nan values. n_stages is added
145-
# to compensate for earlier exclusion of objective columns
14675
row_idx = np.nanargmin(quotients) + self.n_stages
14776
return row_idx, col_idx
14877

14978
def pivot(self, row_idx: int, col_idx: int) -> np.ndarray:
150-
"""Pivots on value on the intersection of pivot row and column.
151-
152-
>>> Tableau(np.array([[-2,-3,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]),
153-
... 2, 2).pivot(1, 0).tolist()
154-
... # doctest: +NORMALIZE_WHITESPACE
155-
[[0.0, 3.0, 2.0, 0.0, 8.0],
156-
[1.0, 3.0, 1.0, 0.0, 4.0],
157-
[0.0, -8.0, -3.0, 1.0, -8.0]]
158-
"""
159-
# Avoid changes to original tableau
16079
piv_row = self.tableau[row_idx].copy()
161-
16280
piv_val = piv_row[col_idx]
16381

164-
# Entry becomes 1
16582
piv_row *= 1 / piv_val
16683

167-
# Variable in pivot column becomes basic, ie the only non-zero entry
16884
for idx, coeff in enumerate(self.tableau[:, col_idx]):
16985
self.tableau[idx] += -coeff * piv_row
17086
self.tableau[row_idx] = piv_row
17187
return self.tableau
17288

17389
def change_stage(self) -> np.ndarray:
174-
"""Exits first phase of the two-stage method by deleting artificial
175-
rows and columns, or completes the algorithm if exiting the standard
176-
case.
177-
178-
>>> Tableau(np.array([
179-
... [3, 3, -1, -1, 0, 0, 4],
180-
... [2, 1, 0, 0, 0, 0, 0.],
181-
... [1, 2, -1, 0, 1, 0, 2],
182-
... [2, 1, 0, -1, 0, 1, 2]
183-
... ]), 2, 2).change_stage().tolist()
184-
... # doctest: +NORMALIZE_WHITESPACE
185-
[[2.0, 1.0, 0.0, 0.0, 0.0],
186-
[1.0, 2.0, -1.0, 0.0, 2.0],
187-
[2.0, 1.0, 0.0, -1.0, 2.0]]
188-
"""
189-
# Objective of original objective row remains
19090
self.objectives.pop()
19191

19292
if not self.objectives:
19393
return self.tableau
19494

195-
# Slice containing ids for artificial columns
19695
s = slice(-self.n_artificial_vars - 1, -1)
19796

198-
# Delete the artificial variable columns
19997
self.tableau = np.delete(self.tableau, s, axis=1)
200-
201-
# Delete the objective row of the first stage
20298
self.tableau = np.delete(self.tableau, 0, axis=0)
20399

204100
self.n_stages = 1
@@ -208,132 +104,52 @@ def change_stage(self) -> np.ndarray:
208104
return self.tableau
209105

210106
def run_simplex(self) -> dict[Any, Any]:
211-
"""Operate on tableau until objective function cannot be
212-
improved further.
213-
214-
# Standard linear program:
215-
Max: x1 + x2
216-
ST: x1 + 3x2 <= 4
217-
3x1 + x2 <= 4
218-
>>> {key: float(value) for key, value in Tableau(np.array([[-1,-1,0,0,0],
219-
... [1,3,1,0,4],[3,1,0,1,4.]]), 2, 0).run_simplex().items()}
220-
{'P': 2.0, 'x1': 1.0, 'x2': 1.0}
221-
222-
# Standard linear program with 3 variables:
223-
Max: 3x1 + x2 + 3x3
224-
ST: 2x1 + x2 + x3 ≤ 2
225-
x1 + 2x2 + 3x3 ≤ 5
226-
2x1 + 2x2 + x3 ≤ 6
227-
>>> {key: float(value) for key, value in Tableau(np.array([
228-
... [-3,-1,-3,0,0,0,0],
229-
... [2,1,1,1,0,0,2],
230-
... [1,2,3,0,1,0,5],
231-
... [2,2,1,0,0,1,6.]
232-
... ]),3,0).run_simplex().items()} # doctest: +ELLIPSIS
233-
{'P': 5.4, 'x1': 0.199..., 'x3': 1.6}
234-
235-
236-
# Optimal tableau input:
237-
>>> {key: float(value) for key, value in Tableau(np.array([
238-
... [0, 0, 0.25, 0.25, 2],
239-
... [0, 1, 0.375, -0.125, 1],
240-
... [1, 0, -0.125, 0.375, 1]
241-
... ]), 2, 0).run_simplex().items()}
242-
{'P': 2.0, 'x1': 1.0, 'x2': 1.0}
243-
244-
# Non-standard: >= constraints
245-
Max: 2x1 + 3x2 + x3
246-
ST: x1 + x2 + x3 <= 40
247-
2x1 + x2 - x3 >= 10
248-
- x2 + x3 >= 10
249-
>>> {key: float(value) for key, value in Tableau(np.array([
250-
... [2, 0, 0, 0, -1, -1, 0, 0, 20],
251-
... [-2, -3, -1, 0, 0, 0, 0, 0, 0],
252-
... [1, 1, 1, 1, 0, 0, 0, 0, 40],
253-
... [2, 1, -1, 0, -1, 0, 1, 0, 10],
254-
... [0, -1, 1, 0, 0, -1, 0, 1, 10.]
255-
... ]), 3, 2).run_simplex().items()}
256-
{'P': 70.0, 'x1': 10.0, 'x2': 10.0, 'x3': 20.0}
257-
258-
# Non standard: minimisation and equalities
259-
Min: x1 + x2
260-
ST: 2x1 + x2 = 12
261-
6x1 + 5x2 = 40
262-
>>> {key: float(value) for key, value in Tableau(np.array([
263-
... [8, 6, 0, 0, 52],
264-
... [1, 1, 0, 0, 0],
265-
... [2, 1, 1, 0, 12],
266-
... [6, 5, 0, 1, 40.],
267-
... ]), 2, 2).run_simplex().items()}
268-
{'P': 7.0, 'x1': 5.0, 'x2': 2.0}
269-
270-
271-
# Pivot on slack variables
272-
Max: 8x1 + 6x2
273-
ST: x1 + 3x2 <= 33
274-
4x1 + 2x2 <= 48
275-
2x1 + 4x2 <= 48
276-
x1 + x2 >= 10
277-
x1 >= 2
278-
>>> {key: float(value) for key, value in Tableau(np.array([
279-
... [2, 1, 0, 0, 0, -1, -1, 0, 0, 12.0],
280-
... [-8, -6, 0, 0, 0, 0, 0, 0, 0, 0.0],
281-
... [1, 3, 1, 0, 0, 0, 0, 0, 0, 33.0],
282-
... [4, 2, 0, 1, 0, 0, 0, 0, 0, 60.0],
283-
... [2, 4, 0, 0, 1, 0, 0, 0, 0, 48.0],
284-
... [1, 1, 0, 0, 0, -1, 0, 1, 0, 10.0],
285-
... [1, 0, 0, 0, 0, 0, -1, 0, 1, 2.0]
286-
... ]), 2, 2).run_simplex().items()} # doctest: +ELLIPSIS
287-
{'P': 132.0, 'x1': 12.000... 'x2': 5.999...}
288-
"""
289-
# Stop simplex algorithm from cycling.
290-
for _ in range(Tableau.maxiter):
291-
# Completion of each stage removes an objective. If both stages
292-
# are complete, then no objectives are left
107+
"""Run simplex algorithm until optimal solution is found."""
108+
109+
for iteration in range(Tableau.maxiter):
110+
293111
if not self.objectives:
294-
# Find the values of each variable at optimal solution
295112
return self.interpret_tableau()
296113

297114
row_idx, col_idx = self.find_pivot()
298115

299-
# If there are no more negative values in objective row
300116
if self.stop_iter:
301-
# Delete artificial variable columns and rows. Update attributes
302117
self.tableau = self.change_stage()
303118
else:
304119
self.tableau = self.pivot(row_idx, col_idx)
305-
return {}
120+
121+
# FIX: raise error instead of returning {}
122+
raise ValueError(
123+
"Simplex algorithm failed to converge.\n"
124+
f"- Iterations performed: {iteration + 1}\n"
125+
f"- Max iterations allowed: {Tableau.maxiter}\n"
126+
f"- Remaining objectives: {self.objectives}\n"
127+
"Possible causes:\n"
128+
"- Cycling due to degeneracy\n"
129+
"- Unbounded feasible region\n"
130+
"- Ill-formed or poorly scaled input\n"
131+
)
306132

307133
def interpret_tableau(self) -> dict[str, float]:
308-
"""Given the final tableau, add the corresponding values of the basic
309-
decision variables to the `output_dict`
310-
>>> {key: float(value) for key, value in Tableau(np.array([
311-
... [0,0,0.875,0.375,5],
312-
... [0,1,0.375,-0.125,1],
313-
... [1,0,-0.125,0.375,1]
314-
... ]),2, 0).interpret_tableau().items()}
315-
{'P': 5.0, 'x1': 1.0, 'x2': 1.0}
316-
"""
317-
# P = RHS of final tableau
318134
output_dict = {"P": abs(self.tableau[0, -1])}
319135

320136
for i in range(self.n_vars):
321-
# Gives indices of nonzero entries in the ith column
322137
nonzero = np.nonzero(self.tableau[:, i])
323138
n_nonzero = len(nonzero[0])
324139

325-
# First entry in the nonzero indices
140+
if n_nonzero == 0:
141+
continue
142+
326143
nonzero_rowidx = nonzero[0][0]
327144
nonzero_val = self.tableau[nonzero_rowidx, i]
328145

329-
# If there is only one nonzero value in column, which is one
330146
if n_nonzero == 1 and nonzero_val == 1:
331147
rhs_val = self.tableau[nonzero_rowidx, -1]
332148
output_dict[self.col_titles[i]] = rhs_val
149+
333150
return output_dict
334151

335152

336153
if __name__ == "__main__":
337154
import doctest
338-
339-
doctest.testmod()
155+
doctest.testmod()

0 commit comments

Comments
 (0)