Software 1.4.1 (Product mix and inspectors). Examples 1.2.1 and 1.2.2:

>>> import numpy as np
>>> from scipy.optimize import linprog
>>> c = [-20, -30]  # maximise 20 x1 + 30 x2: minimise -Z
>>> A = [[10, 40], [3, 2], [2, 3]]         # budget, machine, assembly
>>> b = [4000, 2000, 1400]
>>> r = linprog(c, A_ub=A, b_ub=b, method='highs')
>>> r.x, -r.fun, b - np.array(A) @ r.x     # plan, profit, unused resources
(array([400.,   0.]), 8000.0, array([  0., 800., 600.]))
>>> r = linprog([68.8, 68.0],  # inspectors: daily cost per inspector
...             A_ub=[[-320, -240]], b_ub=[-2000],   # at least 2000 pieces
...             bounds=[(0, 9), (0, 11)], method='highs')
>>> r.x, r.fun
(array([6.25, 0.  ]), 430.0)

Reading the output. The product mix makes \(400\) of \(C_1\) for a profit of \(8000\); the slacks are \(0\) for the budget and \(800\) and \(600\) hours of machine and assembly time. The inspection plan uses \(6.25\) Grade 1 inspectors and no Grade 2, at Rs. 430 a day, a fractional answer that the model’s divisibility assumption allowed.

MATLAB/Octave: linprog(-[20;30], [10 40; 3 2; 2 3], [4000; 2000; 1400], [], [], [0;0]).

Software 1.4.2 (A four-month plan). Example 1.2.3, with the twelve variables ordered as regular production, overtime and stock, and one balance equation per month:

>>> d = np.array([100, 150, 200, 120])     # demand in four months
>>> T = len(d)  # variables: regular, overtime, stock
>>> cost = np.r_[10*np.ones(T), 14*np.ones(T), 2*np.ones(T)]
>>> Aeq = np.zeros((T, 3*T))
>>> for t in range(T):            # stock_{t-1} + R_t + O_t - stock_t = d_t
...     Aeq[t, t] = Aeq[t, T + t] = 1; Aeq[t, 2*T + t] = -1
...     if t > 0: Aeq[t, 2*T + t - 1] = 1
...
>>> bounds = [(0, 140)]*T + [(0, 40)]*T + [(0, None)]*T
>>> r = linprog(cost, A_eq=Aeq, b_eq=d, bounds=bounds, method='highs')
>>> r.x.reshape(3, T), r.fun               # rows: regular, overtime, stock
(array([[140., 140., 140., 120.],
       [  0.,   0.,  30.,   0.],
       [ 40.,  30.,   0.,   0.]]), 5960.0)

Reading the output. Regular production is \(140\) in each of the first three months and \(120\) in the last; overtime of \(30\) is used only in month \(3\); stocks of \(40\) and \(30\) are carried out of months \(1\) and \(2\) into the peak. The total cost is Rs. 5960. Checking month \(3\): \(30 + 140 + 30 - 0 = 200\).

MATLAB/Octave: the same matrices with linprog(cost, [], [], Aeq, d, lb, ub).

Software 1.4.3 (Fitting a line by linear programming). The \(L^\infty \) and \(L^1\) fits of Figure 1 as linear programmes, with the free intercept and slope declared by the bounds (None, None), and least squares for comparison:

>>> xs = np.arange(10.0)
>>> ys = np.array([1.0, 2.1, 2.9, 4.2, 4.8, 6.1, 7.2, 7.9, 9.1, 15.0])
>>> n = len(xs)                   # L-infinity fit: min t, |y - a - bx| <= t
>>> one = np.ones(n)
>>> A = np.r_[np.c_[-one, -xs, -one], np.c_[one, xs, -one]]
>>> r = linprog([0, 0, 1], A_ub=A, b_ub=np.r_[-ys, ys],
...             bounds=[(None, None)]*2 + [(0, None)], method='highs')
>>> cheb = r.x[:2]
>>> I = np.eye(n)  # L1 fit: min sum e_i, |y_i - a - b x_i| <= e_i
>>> A = np.r_[np.c_[-np.ones(n), -xs, -I], np.c_[np.ones(n), xs, -I]]
>>> r = linprog(np.r_[0, 0, np.ones(n)], A_ub=A, b_ub=np.r_[-ys, ys],
...             bounds=[(None, None)]*2 + [(0, None)]*n, method='highs')
>>> l1 = r.x[:2]
>>> ls = np.polyfit(xs, ys, 1)[::-1]       # least squares, for comparison
>>> np.round(np.array([l1, ls, cheb]), 4)  # rows: L1, L2, L-infinity; (a, b)
array([[ 1.    ,  1.02  ],
       [ 0.2945,  1.2745],
       [-1.1722,  1.5556]])

Reading the output. The \(L^1\) line is \(1 + 1.02x\), close to the trend of the first nine points; least squares gives \(0.29 + 1.27x\) and the Chebyshev fit \(-1.17 + 1.56x\). One bad point moves the slope by \(25\%\) and \(53\%\) for the last two. The \(L^1\) fit is the regression analogue of the median, and like it, robust.

MATLAB/Octave: linprog with the same stacked constraint matrices.