Software 7.6.1 (Three starts and the optimum). The north-west corner, least-cost and Vogel rules for Example 7.1.2, and the LP with one equation per row and column, built with Kronecker products:

>>> import numpy as np
>>> from scipy.optimize import linprog
>>> C = np.array([[19, 30, 50, 10], [70, 30, 40, 60], [40, 8, 70, 20]])
>>> s = np.array([7, 9, 18]); d = np.array([5, 8, 7, 14])  # balanced: 34 = 34
>>> def northwest(s, d):
...     s, d, X = s.copy(), d.copy(), np.zeros((len(s), len(d)), int)
...     i = j = 0
...     while i < len(s) and j < len(d):
...         q = min(s[i], d[j]); X[i, j] = q; s[i] -= q; d[j] -= q
...         if s[i] == 0: i += 1
...         else: j += 1
...     return X
...
>>> def least_cost(C, s, d):
...     s, d, X = s.copy(), d.copy(), np.zeros(C.shape, int)
...     for k in np.argsort(C, axis=None, kind='stable'):
...         i, j = divmod(k, C.shape[1]); q = min(s[i], d[j])
...         X[i, j] += q; s[i] -= q; d[j] -= q
...     return X
...
>>> def vogel(C, s, d):
...     s, d, X = s.copy(), d.copy(), np.zeros(C.shape, int)
...     rows, cols = set(range(len(s))), set(range(len(d)))
...     while rows and cols:
...         def pen(costs):                # difference of the two smallest
...             c = sorted(costs); return c[1] - c[0] if len(c) > 1 else c[0]
...         cand = [(pen([C[i, j] for j in cols]), 'r', i) for i in rows]
...         cand += [(pen([C[i, j] for i in rows]), 'c', j) for j in cols]
...         p, kind, k = max(cand, key=lambda t: t[0])
...         if kind == 'r': i = k; j = min(cols, key=lambda j: C[i, j])
...         else: j = k; i = min(rows, key=lambda i: C[i, j])
...         q = min(s[i], d[j]); X[i, j] = q; s[i] -= q; d[j] -= q
...         if s[i] == 0: rows.discard(i)
...         elif d[j] == 0: cols.discard(j)
...     return X
...
>>> starts = {'northwest': northwest(s, d), 'least cost': least_cost(C, s, d),
...           'Vogel': vogel(C, s, d)}
>>> for name, X in starts.items():
...     cost = int((C * X).sum())
...     print(f"{name:10s} cost {cost:4d}, cells {(X > 0).sum()}")
...
northwest  cost 1015, cells 6
least cost cost  814, cells 6
Vogel      cost  779, cells 6
>>> m, n = C.shape                   # the LP: one equation per row and column
>>> rowsum = np.kron(np.eye(m), np.ones(n))
>>> colsum = np.kron(np.ones(m), np.eye(n))
>>> Aeq = np.vstack([rowsum, colsum])
>>> r = linprog(C.ravel(), A_eq=Aeq, b_eq=np.r_[s, d], method='highs')
>>> X = r.x.reshape(m, n) + 0.0; X, r.fun
(array([[ 5.,  0.,  0.,  2.],
                                                                                      
                                                                                      
       [ 0.,  2.,  7.,  0.],
       [ 0.,  6.,  0., 12.]]), 743.0)

Reading the output. The three starts cost \(1015\), \(814\) and \(779\), each using \(6 = m + n - 1\) cells. The optimum is \(743\): ship \(5\) and \(2\) from \(S_1\) to \(D_1\) and \(D_4\), \(2\) and \(7\) from \(S_2\) to \(D_2\) and \(D_3\), \(6\) and \(12\) from \(S_3\) to \(D_2\) and \(D_4\). The LP solution is integral, as total unimodularity guarantees.

MATLAB/Octave: Aeq = [kron(eye(m), ones(1,n)); kron(ones(1,m), eye(n))].

Software 7.6.2 (MODI multipliers). The dual values of the supply and demand equations from HiGHS, shifted so that \(u_1 = 0\), the reduced costs \(c_{ij} - u_i - v_j\), and the rank of the constraint matrix:

>>> u = r.eqlin.marginals[:m]; v = r.eqlin.marginals[m:]
>>> u, v = u - u[0], v + u[0]        # normalise u1 = 0; u + v is unchanged
>>> red = C - u[:, None] - v[None, :]      # c_ij - u_i - v_j
>>> u, v
(array([ 0., 32., 10.]), array([19., -2.,  8., 10.]))
>>> np.round(red, 6) + 0.0
array([[ 0., 32., 42.,  0.],
       [19.,  0.,  0., 18.],
       [11.,  0., 52.,  0.]])
>>> np.linalg.matrix_rank(Aeq), m + n - 1  # one equation is redundant
(6, 6)

Reading the output. \(u = (0, 32, 10)\) and \(v = (19, -2, 8, 10)\). The reduced costs are \(0\) on the six basic cells and positive elsewhere, so the plan is optimal and, since no non-basic cell has \(d_{ij} = 0\), it is the unique optimum. The rank is \(6 = m + n - 1\): one of the seven equations is redundant, which is why one multiplier can be set freely.

MATLAB/Octave: lambda.eqlin from linprog.

Software 7.6.3 (Total unimodularity and assignment). A brute-force check of total unimodularity on the \(2 \times 3\) transportation matrix and on the incidence matrix of a triangle; then Example 7.5.1 by SciPy’s assignment solver and by its LP relaxation:

>>> from itertools import combinations
>>> from scipy.optimize import linear_sum_assignment
>>> def totally_unimodular(M):             # every square minor in {-1, 0, 1}
...     r, c = M.shape
...     for k in range(1, min(r, c) + 1):
...         for R in combinations(range(r), k):
...             for K in combinations(range(c), k):
...                 if round(abs(np.linalg.det(M[np.ix_(R, K)]))) > 1:
...                     return False
...     return True
...
>>> small = np.vstack([np.kron(np.eye(2), np.ones(3)),   # a 2 x 3 transport
...                    np.kron(np.ones(2), np.eye(3))])
>>> triangle = np.array([[1, 1, 0], [0, 1, 1], [1, 0, 1]])  # an odd cycle
>>> totally_unimodular(small), totally_unimodular(triangle)
(True, False)
>>> W = np.array([[9, 11, 14, 11], [6, 15, 13, 13],
...               [12, 13, 6, 8], [11, 9, 10, 9]])
>>> rows, cols = linear_sum_assignment(W)  # Hungarian-type algorithm
>>> cols, W[rows, cols].sum()
(array([1, 0, 2, 3]), 32)
>>> I4, e4 = np.eye(4), np.ones(4)
>>> A2 = np.vstack([np.kron(I4, e4), np.kron(e4, I4)])
>>> r = linprog(W.ravel(), A_eq=A2, b_eq=np.ones(8), method='highs')
>>> r.x.reshape(4, 4) + 0.0, r.fun   # the LP relaxation is already 0-1
(array([[0., 1., 0., 0.],
       [1., 0., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]]), 32.0)

Reading the output. The transportation matrix is TU; the triangle’s is not (its determinant is \(2\)). The assignment solver sends workers \(1, 2, 3, 4\) to jobs \(2, 1, 3, 4\) (Python counts from \(0\)) at cost \(32\), as the Hungarian method found, and the LP relaxation returns a permutation matrix with the same cost: the LP needs no integrality constraints.

MATLAB/Octave: matchpairs(W, 1e6) solves the assignment problem.