Software 5.6.1 (A primal–dual pair). Example 5.1.1 and its dual solved separately, the multipliers that HiGHS returns with the primal solution, and the products in the complementary slackness conditions:

>>> 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])
>>> P = linprog(-c, A_ub=A, b_ub=b, method='highs')        # primal: max c x
>>> D = linprog(b, A_ub=-A.T, b_ub=-c, method='highs')     # dual: min b y
>>> P.x, -P.fun, D.x, D.fun                # A^T y >= c, y >= 0
(array([ 8., 10.]), 230.0, array([1.66666667, 3.33333333, 0.        ]), 230.0)
>>> -P.ineqlin.marginals                   # HiGHS's multipliers: the dual
array([1.66666667, 3.33333333, 0.        ])
>>> slack = b - A @ P.x; surplus = A.T @ D.x - c
>>> np.round(slack * D.x, 9) + 0, np.round(surplus * P.x, 9) + 0   # all 0
(array([0., 0., 0.]), array([0., 0.]))

Reading the output. The primal optimum \((8, 10)\) and the dual optimum \((\frac 53, \frac {10}3, 0)\) have the same value \(230\). HiGHS reports the dual solution as the marginals of the inequality constraints (with a minus sign, because it minimised \(-Z\)), so solving the dual separately is never necessary. Every product \(y_i \times \text {slack}_i\) and \(x_j \times \text {surplus}_j\) is zero.

MATLAB/Octave: [x, fval, flag, out, lambda] = linprog(...); the dual is lambda.ineqlin.

Software 5.6.2 (The dual simplex method). A dual simplex method on the tableau, in exact arithmetic, for Example 5.5.1; each line shows the basis (columns \(0, 1\) are \(x_1, x_2\), \(2, 3, 4\) the slacks), the right-hand sides and the objective value:

>>> from fractions import Fraction as F
>>> def dual_simplex(T, basis):
...     """T: rows [A | b], last row [c | -z] with c >= 0 (dual feasible).
...     Pivots until b >= 0 (primal feasible), printing each basis."""
...     while True:
...         rhs = [str(T[i][-1]) for i in range(len(T) - 1)]
...         print(basis, rhs, -T[-1][-1])
...         r = min(range(len(T) - 1), key=lambda i: T[i][-1])
...         if T[r][-1] >= 0: return
...         cols = [j for j in range(len(T[0]) - 1) if T[r][j] < 0]
...         if not cols: raise ValueError('primal infeasible')
...         j = min(cols, key=lambda j: T[-1][j] / -T[r][j])   # ratio test
...         T[r] = [v / T[r][j] for v in T[r]]
...         for i in range(len(T)):
...             if i != r:
...                 T[i] = [a - T[i][j] * p for a, p in zip(T[i], T[r])]
...         basis[r] = j
...
>>> T = [[F(v) for v in row] for row in
...      [[1, 2, 1, 0, 0, 40],              # x1 + 2 x2 <= 40
...       [-3, -1, 0, 1, 0, -30],           # 3 x1 + x2 >= 30, times -1
...       [-4, -3, 0, 0, 1, -60],           # 4 x1 + 3 x2 >= 60, times -1
...       [20, 10, 0, 0, 0, 0]]]            # minimise 20 x1 + 10 x2
>>> dual_simplex(T, [2, 3, 4])             # start from the slack basis
[2, 3, 4] ['40', '-30', '-60'] 0
[2, 3, 1] ['0', '-10', '20'] 200
[2, 0, 1] ['10', '6', '12'] 240

Reading the output. The right-hand sides start at \((40, -30, -60)\) and become non-negative after two pivots. The objective rises from \(0\) to \(200\) to \(240\): in the dual simplex method the value approaches the optimum from below for a minimisation, since each intermediate basis gives a feasible dual solution and hence a lower bound.

MATLAB/Octave: the same row operations; linprog with ’dual-simplex’ uses this method.

Software 5.6.3 (A certificate of infeasibility). For the infeasible system \(x_1 + x_2 \leq 2\), \(-x_1 - x_2 \leq -3\), a Farkas vector found by linear programming, normalised so that its entries sum to \(1\):

>>> A = np.array([[1., 1], [-1, -1]]); b = np.array([2., -3])  # x1 + x2 <= 2,
>>> far = linprog(b, A_eq=np.vstack([A.T, np.ones(2)]),        # x1 + x2 >= 3
...               b_eq=[0, 0, 1], method='highs')  # y >= 0, y A = 0, sum 1
>>> far.x, far.fun                         # y b < 0: a proof of infeasibility
(array([0.5, 0.5]), -0.5)

Reading the output. \(y = (\frac 12, \frac 12)\) satisfies \(A^{\mathsf T}y = 0\) with \(b^{\mathsf T}y = -\frac 12 < 0\): adding half of each constraint gives \(0 \leq -\frac 12\), an impossibility that proves the system has no solution.

MATLAB/Octave: the same linprog call with Aeq = [A’; ones(1,2)].