Software 4.4.1 (The revised simplex method). A revised simplex method with an explicit \(B^{-1}\) updated by eta matrices (Proposition 4.1.1), printing the basis, \(x_B\) and the multipliers \(y\) at each iteration, on Example 4.1.2:

>>> import numpy as np
>>> def revised(A, b, c, basis, tol=1e-9):
...     """Minimise c x, A x = b, x >= 0, from the feasible basis given.
...     Keeps B^{-1} explicitly and updates it by an eta matrix each pivot."""
...     m, n = A.shape
...     Binv = np.linalg.inv(A[:, basis])
...     while True:
...         xB = Binv @ b
...         y = c[basis] @ Binv                        # simplex multipliers
...         rc = c - y @ A                             # reduced costs
...         rc[basis] = 0
...         j = int(np.argmin(rc))
...         print(basis, np.round(xB, 3), 'y =', np.round(y, 3))
...         if rc[j] > -tol: return basis, xB, c[basis] @ xB
...         u = Binv @ A[:, j]                         # the entering column
...         ok = u > tol
...         if not ok.any(): raise ValueError('unbounded')
...         r = np.where(ok)[0][np.argmin(xB[ok] / u[ok])]
...         E = np.eye(m); E[:, r] = -u / u[r]; E[r, r] = 1 / u[r]   # eta
...         Binv = E @ Binv; basis[r] = j
...
>>> A = np.array([[2., 1, 1, 0, 0], [2, 4, 0, 1, 0], [-1, 1, 0, 0, 1]])
>>> b = np.array([26., 56, 5]); c = np.array([-10., -15, 0, 0, 0])
>>> basis, xB, z = revised(A, b, c, [2, 3, 4])
[2, 3, 4] [26. 56.  5.] y = [0. 0. 0.]
[2, 3, 1] [21. 36.  5.] y = [  0.   0. -15.]
[2, 0, 1] [ 3.  6. 11.] y = [ 0.    -4.167  1.667]
[4, 0, 1] [ 3.  8. 10.] y = [-1.667 -3.333  0.   ]
>>> -z
230.0

Reading the output. The bases \(\{s_1, s_2, s_3\}\), \(\{s_1, s_2, x_2\}\), \(\{s_1, x_1, x_2\}\) and \(\{s_3, x_1, x_2\}\) (columns \(2, 3, 4\) are the slacks) are those of the tableau, visiting \(O\), \(D\), \(C\), \(B\). The final multipliers \(y = (-1.667, -3.333, 0)\) are minus the shadow prices, and the optimal value is \(230\).

MATLAB/Octave: y = c(basis)/B solves \(y^{\mathsf T}B = c_B^{\mathsf T}\); u = B∖A(:,j) gives the entering column.

Software 4.4.2 (Bounded variables and flips). Maximise \(3x_1 + 5x_2 + 2x_3\) subject to \(x_1 + x_2 + x_3 \leq 10\), \(2x_1 + x_2 \leq 12\), \(x_2 + 2x_3 \leq 14\), with \(0 \leq x_1 \leq 4\), \(0 \leq x_2 \leq 3\), \(0 \leq x_3 \leq 8\): the bounded-variable method, printing each step’s entering variable and which limit stopped it; then HiGHS with the bounds:

>>> from scipy.optimize import linprog
>>> def bounded(A, b, c, up, tol=1e-9):
...     """Minimise c x, A x + s = b, 0 <= x <= up, s >= 0 (b >= 0): the
...     simplex method with non-basic variables at either bound."""
...     m, n = A.shape
...     A = np.hstack([A, np.eye(m)]); c = np.r_[c, np.zeros(m)]
...     up = np.r_[up, np.full(m, np.inf)]
...     x = np.r_[np.zeros(n), b]; basis = list(range(n, n + m))
...     while True:
...         Binv = np.linalg.inv(A[:, basis])
...         rc = c - (c[basis] @ Binv) @ A
...         free = [j for j in range(n + m) if j not in basis and
...                 ((x[j] < tol and rc[j] < -tol) or
...                  (x[j] > up[j] - tol and rc[j] > tol))]
...         if not free: return x[:n], c @ x
...         j = min(free, key=lambda k: -abs(rc[k]))
...         d = 1.0 if x[j] < tol else -1.0          # raise or lower x_j
...         u = d * (Binv @ A[:, j])                 # basics move by -theta u
...         steps = [(up[j], 'flip', None)]    # x_j reaches its other bound
...         for i, k in enumerate(basis):
...             if u[i] > tol: steps.append((x[k] / u[i], 'to 0', i))
...             if u[i] < -tol:
...                 steps.append(((up[k] - x[k]) / -u[i], 'to up', i))
...         th, kind, i = min(steps, key=lambda t: t[0])
...         x[basis] -= th * u; x[j] += d * th
...         name = f"x{j + 1}" if j < n else f"s{j - n + 1}"
...         if kind != 'flip':
...             out = basis[i]; basis[i] = j
...             who = f"x{out + 1}" if out < n else f"s{out - n + 1}"
...             kind += f", {who} leaves"
...         print(f"{name:3s} {kind:18s} x = {np.round(x[:n], 3)}")
...
>>> c = np.array([-3., -5, -2])            # maximise 3 x1 + 5 x2 + 2 x3
>>> A = np.array([[1., 1, 1], [2, 1, 0], [0, 1, 2]])
>>> b = np.array([10., 12, 14]); up = np.array([4., 3, 8])
>>> x, z = bounded(A, b, c, up)
x2  flip               x = [0. 3. 0.]
x1  flip               x = [4. 3. 0.]
x3  to 0, s1 leaves    x = [4. 3. 3.]
>>> x, -z
(array([4., 3., 3.]), 33.0)
>>> r = linprog(c, A_ub=A, b_ub=b, bounds=[(0, u) for u in up],
...             method='highs')
>>> r.x, -r.fun                            # the same answer from HiGHS
(array([4., 3., 3.]), 33.0)

Reading the output. The first two iterations are bound flips: \(x_2\) rises to its bound \(3\) and \(x_1\) to \(4\), with no change of basis. Only in the third does a constraint stop the step: \(x_3\) rises until the first constraint is tight (\(4 + 3 + 3 = 10\)), and its slack leaves the basis. The optimum is \((4, 3, 3)\), \(Z = 33\). Treating the bounds as rows would have doubled the number of constraints, from \(3\) to \(6\).

MATLAB/Octave: linprog(c, A, b, [], [], [0;0;0], [4;3;8]).

Software 4.4.3 (Simplex against interior points). A random sparse LP with \(400\) constraints and \(1000\) variables bounded in \([0, 1]\), solved by HiGHS’s dual simplex and interior-point methods: status, optimal value, iterations and time:

>>> import time
>>> rng = np.random.default_rng(3)
>>> m, n = 400, 1000                       # a random feasible, bounded LP
>>> A = rng.uniform(0, 1, (m, n)) * (rng.uniform(size=(m, n)) < 0.05)
>>> b = A @ rng.uniform(0, 1, n) + 1
>>> c = -rng.uniform(0, 1, n)
>>> for method in ('highs-ds', 'highs-ipm'):
...     t0 = time.perf_counter()
...     r = linprog(c, A_ub=A, b_ub=b, bounds=(0, 1), method=method)
...     print(method, r.status, round(r.fun, 6), r.nit,
...           f"{time.perf_counter() - t0:.3f}s")
...
highs-ds 0 -366.485484 1043 0.297s
highs-ipm 0 -366.485484 19 0.146s

Reading the output. Both methods reach the same optimum, \(-366.485484\). The simplex method takes about a thousand cheap iterations, roughly \(2.5\) per constraint; the interior-point method takes \(19\) expensive ones. On this small problem both take a fraction of a second, and the times depend on the machine.

MATLAB/Octave: optimoptions(’linprog’,’Algorithm’,’dual-simplex’) or ’interior-point’.