Software 6.5.1 (Ranging). From the optimal basis \(\{s_3, x_1, x_2\}\) of the running example: \(x_B\) and \(y\), the range of each \(b_i\) from Theorem 6.1.1, and the range of each \(c_j\) over which the reduced costs keep their sign, found on a fine grid:

>>> import numpy as np
>>> from scipy.optimize import linprog
>>> A = np.array([[2., 1], [2, 4], [-1, 1]]); b = np.array([26., 56, 5])
>>> c = np.array([10., 15])
>>> Af = np.hstack([A, np.eye(3)]); cf = np.r_[c, 0, 0, 0]   # with slacks
>>> basis = [4, 0, 1]                      # s3, x1, x2: the optimal basis
>>> Binv = np.linalg.inv(Af[:, basis]); xB = Binv @ b
>>> y = cf[basis] @ Binv; xB, y          # basic values and shadow prices
(array([ 3.,  8., 10.]), array([1.66666667, 3.33333333, 0.        ]))
>>> for i in range(3):                     # right-hand-side ranges
...     col = Binv[:, i]; t = -xB / np.where(col == 0, np.nan, col)
...     lo = np.nanmax(np.where(col > 0, t, np.nan), initial=-np.inf)
...     hi = np.nanmin(np.where(col < 0, t, np.nan), initial=np.inf)
...     print(f"b{i+1}: {b[i] + lo:6.2f} to {b[i] + hi:6.2f}, y = {y[i]:.4f}")
...
b1:  23.00 to  56.00, y = 1.6667
b2:  26.00 to  62.00, y = 3.3333
b3:   2.00 to    inf, y = 0.0000
>>> def still_optimal(cc):           # all reduced costs <= 0 (maximising)
...     cfull = np.r_[cc, 0, 0, 0]
...     return np.all(cfull - (cfull[basis] @ Binv) @ Af <= 1e-9)
...
>>> grid = np.linspace(0, 60, 60001)       # cost ranges, on a fine grid
>>> for j in range(2):
...     trial = lambda t: np.where(np.arange(2) == j, t, c)
...     ok = [t for t in grid if still_optimal(trial(t))]
...     print(f"c{j+1}: {min(ok):.2f} to {max(ok):.2f}")
...
c1: 7.50 to 30.00
c2: 5.00 to 20.00

Reading the output. The ranges are \(23 \leq b_1 \leq 56\), \(26 \leq b_2 \leq 62\), \(b_3 \geq 2\), and \(7.5 \leq c_1 \leq 30\), \(5 \leq c_2 \leq 20\). The third resource has price \(0\) and can grow without limit, but cannot fall below \(2\), where its slack of \(3\) is used up.

MATLAB/Octave: Binv = inv(Af(:, basis)) and the same ratio formulas.

Software 6.5.2 (The value function). The optimal profit for \(241\) values of \(b_2\) from \(0\) to \(120\), its breakpoints, and the distinct slopes between them:

>>> bs = np.linspace(0, 120, 241); z = []
>>> for b2 in bs:                          # parametric right-hand side b2
...     r = linprog(-c, A_ub=A, b_ub=[26, b2, 5], method='highs')
...     z.append(-r.fun)
...
>>> slopes = np.round(np.diff(z) / np.diff(bs), 4)
>>> breaks = [bs[k] for k in range(1, len(slopes))
...           if slopes[k] != slopes[k - 1]]
>>> breaks, sorted(set(slopes), reverse=True)
([26.0, 62.0], [5.0, 3.3333, 0.0])

Reading the output. The slope changes at \(b_2 = 26\) and \(b_2 = 62\), exactly the ends of the range found by ranging, and the slopes are \(5\), \(\frac {10}3\) and \(0\), decreasing as Theorem 6.1.3 requires. Below \(26\) the constraints \(2x_1 + 4x_2 \leq b_2\) and \(-x_1 + x_2 \leq 5\) bind and an extra unit of \(b_2\) is worth \(5\); above \(62\) the second resource is in surplus.

MATLAB/Octave: a loop of linprog calls.

Software 6.5.3 (A new product and a new constraint). Pricing a new product with column \((1, 1, 1)\) and profit \(12\), then re-solving with it; and re-solving the original problem with the added constraint \(x_1 \leq 6\):

>>> a3, c3 = np.array([1., 1, 1]), 12.0    # a new product: price it out
>>> y @ a3, c3 - y @ a3                    # reduced profit > 0: worth making
(5.0, 7.0)
>>> r = linprog(-np.r_[c, c3], A_ub=np.c_[A, a3], b_ub=b, method='highs')
>>> r.x, -r.fun
(array([ 7., 10.,  2.]), 244.0)
>>> r = linprog(-c, A_ub=np.r_[A, [[1, 0]]], b_ub=np.r_[b, 6],
...             method='highs')            # a new constraint x1 <= 6
>>> r.x, -r.fun
(array([ 6., 11.]), 225.0)

Reading the output. The resources the new product uses are worth \(5\), so a profit of \(12\) gives a reduced profit of \(7 > 0\): it should be made. The new optimum makes \(7\), \(10\) and \(2\) units for a profit of \(244\). The constraint \(x_1 \leq 6\) cuts off the old optimum \((8, 10)\), and the new optimum is the vertex \(C(6, 11)\) with \(Z = 225\): the restriction costs \(5\).

MATLAB/Octave: linprog with the extra column or row.