Software 3.6.1 (The tableau method). The functions pivot and simplex, and Example 3.2.1: the basic columns, the vertex \((x_1, x_2)\) and \(Z\) after each pivot (columns \(0, 1\) are \(x_1, x_2\) and \(2, 3, 4\) the slacks):

>>> from fractions import Fraction as F
>>> def pivot(T, r, j):             # make column j a unit vector, 1 in row r
...     T[r] = [v / T[r][j] for v in T[r]]
...     for i in range(len(T)):
...         if i != r and T[i][j] != 0:
...             T[i] = [a - T[i][j] * p for a, p in zip(T[i], T[r])]
...
>>> def simplex(T, basis, n, rule='dantzig', limit=100):
...     """T: rows [A | b], last row [-c | z] (maximise c x); basis[i] is the
...     basic column of row i; n original variables.  Returns the status and
...     the history of (basis, x, z) at each tableau."""
...     def state():
...         x = [F(0)] * n
...         for i, j in enumerate(basis):
...             if j < n: x[j] = T[i][-1]
...         return tuple(basis), tuple(x), T[-1][-1]
...     hist = [state()]
...     for _ in range(limit):
...         z = T[-1][:-1]
...         enter = [j for j, v in enumerate(z) if v < 0]
...         if not enter: return 'optimal', hist
...         if rule == 'bland': j = min(enter)        # smallest index
...         else: j = min(enter, key=lambda k: (z[k], k))   # most negative
...         rows = [i for i in range(len(T) - 1) if T[i][j] > 0]
...         if not rows: return 'unbounded', hist
...         tie = (lambda i: basis[i]) if rule == 'bland' else (lambda i: i)
...         key = lambda i: (T[i][-1] / T[i][j], tie(i))
...         r = min(rows, key=key)
...         pivot(T, r, j); basis[r] = j; hist.append(state())
...         if [h[0] for h in hist].count(tuple(basis)) > 1:
...             return 'cycling', hist
...     return 'iteration limit', hist
...
>>> def tableau(rows): return [[F(v) for v in r] for r in rows]
>>> T = tableau([[2, 1, 1, 0, 0, 26],      # Example 2.3.1, slacks s1, s2, s3
...              [2, 4, 0, 1, 0, 56],
...              [-1, 1, 0, 0, 1, 5],
...              [-10, -15, 0, 0, 0, 0]])  # objective row: max 10 x1 + 15 x2
>>> status, hist = simplex(T, [2, 3, 4], n=2)
>>> status
'optimal'
>>> for B, x, z in hist:                   # basic columns, vertex, objective
...     print(B, [str(v) for v in x], z)
...
(2, 3, 4) ['0', '0'] 0
(2, 3, 1) ['0', '5'] 75
(2, 0, 1) ['6', '11'] 225
(4, 0, 1) ['8', '10'] 230

Reading the output. The method passes through \((0, 0)\), \((0, 5)\), \((6, 11)\) and \((8, 10)\), with \(Z = 0\), \(75\), \(225\), \(230\), the path of Figure 3; each pivot changes one basic column.

MATLAB/Octave: the same row operations, T(r,:) = T(r,:)/T(r,j) and so on; with the Symbolic Math Toolbox, sym entries keep them exact.

Software 3.6.2 (Two phases). Minimise \(20x_1 + 10x_2\) subject to \(x_1 + 2x_2 \leq 40\), \(3x_1 + x_2 \geq 30\), \(4x_1 + 3x_2 \geq 60\) (the problem of Software 2.4.3), with artificial variables in the two \(\geq \) rows. Phase I maximises \(-(a_1 + a_2)\), after subtracting the artificial rows from the objective row so that it is expressed in the non-basic variables; Phase II restores the true objective and drops the artificial columns:

>>> T = tableau([[1, 2, 1, 0, 0, 0, 0, 40],    # Phase I: x1 + 2 x2 + s1 = 40
...              [3, 1, 0, -1, 0, 1, 0, 30],   # 3 x1 + x2 - s2 + a1 = 30
...              [4, 3, 0, 0, -1, 0, 1, 60],   # 4 x1 + 3 x2 - s3 + a2 = 60
...              [0, 0, 0, 0, 0, 1, 1, 0]])    # max -(a1 + a2)
>>> for k in (1, 2): T[-1] = [u - v for u, v in zip(T[-1], T[k])]  # price out
>>> basis = [2, 5, 6]
>>> status, hist = simplex(T, basis, n=2)
>>> status, -T[-1][-1], basis              # Phase I value 0: feasible
('optimal', Fraction(0, 1), [2, 0, 1])
>>> T[-1] = [F(20), F(10)] + [F(0)] * 6   # Phase II: maximise -cost
>>> for i, j in enumerate(basis):
...     if T[-1][j] != 0:                  # price out the basic columns
...         T[-1] = [u - T[-1][j] * v for u, v in zip(T[-1], T[i])]
...
>>> for row in T: row[5] = row[6] = F(0)   # drop the artificial columns
>>> status, hist = simplex(T, basis, n=2)
>>> status, [str(v) for v in hist[-1][1]], -hist[-1][2]
('optimal', ['6', '12'], Fraction(240, 1))

Reading the output. Phase I reaches \(a_1 + a_2 = 0\) with basis \(\{s_1, x_1, x_2\}\) (columns \(2, 0, 1\)), no artificial variable remaining, so the problem is feasible. Phase II finds it already optimal there: \((6, 12)\) at cost \(240\), as enumeration gave.

MATLAB/Octave: linprog handles the start internally; the two phases are for understanding.

Software 3.6.3 (Cycling and the Klee–Minty cube). Beale’s example (maximise \(10x_1 - 57x_2 - 9x_3 - 24x_4\) subject to three constraints with right-hand sides \(0\), \(0\) and \(1\)) under Dantzig’s rule and Bland’s rule; then the number of pivots Dantzig’s rule takes on the Klee–Minty cube in dimensions \(2\) to \(8\):

>>> beale = [[F(1, 2), F(-11, 2), F(-5, 2), 9, 1, 0, 0, 0],  # Beale, as in
...          [F(1, 2), F(-3, 2), F(-1, 2), 1, 0, 1, 0, 0],    # Chvatal
...          [1, 0, 0, 0, 0, 0, 1, 1],
...          [-10, 57, 9, 24, 0, 0, 0, 0]]
>>> for rule in ('dantzig', 'bland'):
...     status, hist = simplex(tableau(beale), [4, 5, 6], n=4, rule=rule)
...     print(rule, status, len(hist) - 1, 'pivots, z =', hist[-1][2])
...
dantzig cycling 6 pivots, z = 0
bland optimal 7 pivots, z = 1
>>> def klee_minty(n):                     # Klee and Minty's cube
...     rows = []
...     for i in range(n):
...         r = [2**(i - j + 1) if j < i else int(j == i) for j in range(n)]
...         slack = [int(k == i) for k in range(n)]
...         rows.append(r + slack + [5**(i + 1)])
...     rows.append([-2**(n - 1 - j) for j in range(n)] + [0] * n + [0])
...     return tableau(rows)
...
>>> [len(simplex(klee_minty(n), list(range(n, 2*n)), n, limit=10**4)[1]) - 1
...  for n in range(2, 9)]                 # pivots taken by Dantzig's rule
[3, 7, 15, 31, 63, 127, 255]

Reading the output. With Dantzig’s rule the method returns to its starting basis after six degenerate pivots, with \(Z\) still \(0\): it would cycle forever. Bland’s rule escapes after seven pivots and finds the optimum \(Z = 1\). On the Klee–Minty cubes the pivot counts are \(3, 7, 15, \ldots , 255\), exactly \(2^n - 1\): the method visits every vertex.

MATLAB/Octave: the same loops; linprog uses anti-cycling safeguards and never cycles.