Software 9.6.1 (Branch and bound). A depth-first branch and bound on Example 9.3.2, with each subproblem’s LP solved by linprog and printed with its bounds; then SciPy’s milp, which calls the HiGHS branch-and-cut solver:

>>> import numpy as np
>>> from scipy.optimize import linprog, milp, LinearConstraint, Bounds
>>> c = np.array([5., 4])                  # max 5 x1 + 4 x2
>>> A = np.array([[6., 4], [1, 2]]); b = np.array([24., 6])
>>> best = {'z': -np.inf, 'x': None}
>>> def branch(bounds, depth=0):
...     r = linprog(-c, A_ub=A, b_ub=b, bounds=bounds, method='highs')
...     tag = '  ' * depth + str([(lo, hi) for lo, hi in bounds])
...     if r.status != 0: print(tag, 'infeasible'); return
...     z, x = -r.fun, r.x
...     if z <= best['z'] + 1e-9: print(tag, f"bound {z:.2f}: pruned"); return
...     frac = [j for j in range(2) if abs(x[j] - round(x[j])) > 1e-9]
...     if not frac:
...         best.update(z=z, x=x + 0.0)
...         print(tag, f"integer {x + 0.0}, z = {z:.0f}"); return
...     j = frac[0]
...     print(tag, f"LP {np.round(x, 3)}, z = {z:.2f}, branch x{j+1}")
...     lo, hi = bounds[j]
...     down = list(bounds); down[j] = (lo, np.floor(x[j]))
...     up = list(bounds); up[j] = (np.ceil(x[j]), hi)
...     branch(down, depth + 1); branch(up, depth + 1)
...
>>> branch([(0, None), (0, None)])
[(0, None), (0, None)] LP [3.  1.5], z = 21.00, branch x2
  [(0, None), (0, 1.0)] LP [3.333 1.   ], z = 20.67, branch x1
    [(0, 3.0), (0, 1.0)] integer [3. 1.], z = 19
    [(4.0, None), (0, 1.0)] integer [4. 0.], z = 20
  [(0, None), (2.0, None)] bound 18.00: pruned
>>> best['z'], best['x']
(20.0, array([4., 0.]))
>>> r = milp(-c, constraints=LinearConstraint(A, ub=b), integrality=[1, 1])
>>> r.x, -r.fun                            # SciPy's branch and cut (HiGHS)
(array([4., 0.]), 20.0)

Reading the output. Five nodes reproduce the tree of Figure 9: the root, the node \(x_2 \leq 1\), its two integer children, and the pruned node \(x_2 \geq 2\). Both methods find \((4, 0)\) with \(z = 20\).

MATLAB/Octave: intlinprog(-c, [1 2], A, b).

Software 9.6.2 (A Gomory cut). The \(x_2\) row of the optimal tableau of Example 9.3.2, computed as row \(2\) of \(B^{-1}[A \mid I]\); its fractional parts on the slacks; the cut rewritten in \(x_1\), \(x_2\); and the LP with the cut added:

>>> from fractions import Fraction as F
>>> B = np.array([[6., 4], [1, 2]]); Binv = np.linalg.inv(B)
>>> row_x2 = Binv[1] @ np.hstack([A, np.eye(2)])   # x2 row: coefficients
>>> rhs = Binv[1] @ b
>>> [str(F(v).limit_denominator(100)) for v in row_x2], str(F(rhs))
(['0', '1', '-1/8', '3/4'], '3/2')
>>> frac = lambda v: v - np.floor(v)  # cut: sum frac(a_j) s_j >= frac(b)
>>> cut = frac(row_x2[2:]); cut, frac(rhs)  # on the slacks s1, s2
(array([0.875, 0.75 ]), 0.5)
>>> g = -(cut @ A); h = frac(rhs) - cut @ b       # g x >= h
>>> g, h
(array([-6., -5.]), -25.0)
>>> r = linprog(-c, A_ub=np.vstack([A, -g]), b_ub=np.r_[b, -h],
...             method='highs')
>>> r.x, -r.fun                            # the cut moves the LP optimum
(array([3.33333333, 1.        ]), 20.666666666666668)

Reading the output. The row is \(x_2 - \frac 18s_1 + \frac 34s_2 = \frac 32\); the cut \(0.875s_1 + 0.75s_2 \geq 0.5\) becomes \(-6x_1 - 5x_2 \geq -25\), that is \(6x_1 + 5x_2 \leq 25\). With it the LP optimum moves from \((3, 1.5)\) to \((3.33, 1)\), and the bound falls from \(21\) to \(20.67\); further cuts would be needed to reach \(20\).

MATLAB/Octave: B ∖ [A eye(2)] gives the tableau rows.

Software 9.6.3 (Inspectors and goals). The inspector problem of Example 1.2.2 as an IP; then a goal programme with variables \(x_1, x_2\) and deviations \(d_1^\mp \) (profit \(40x_1 + 30x_2\), goal at least \(1000\)) and \(d_2^\mp \) (labour \(3x_1 + 4x_2\), goal at most \(60\) hours), with capacity \(x_1 + x_2 \leq 30\), under two weightings:

>>> pieces = LinearConstraint([[320, 240]], lb=2000)    # pieces a day
>>> r = milp([68.8, 68.0], constraints=pieces,
...          integrality=[1, 1], bounds=Bounds([0, 0], [9, 11]))
>>> r.x, round(r.fun, 1)              # the inspectors of Section 1, as an IP
(array([4., 3.]), 479.2)
>>> Ae = [[40, 30, 1, -1, 0, 0],           # 40 x1 + 30 x2 + d1- - d1+ = 1000
...       [3, 4, 0, 0, 1, -1]]             # 3 x1 + 4 x2 + d2- - d2+ = 60
>>> for name, w in (('profit first ', [0, 0, 1000, 0, 0, 1]),
...                 ('labour x 20  ', [0, 0, 1, 0, 0, 20])):
...     r = linprog(w, A_eq=Ae, b_eq=[1000, 60], A_ub=[[1, 1, 0, 0, 0, 0]],
...                 b_ub=[30], method='highs')        # and x1 + x2 <= 30
...     print(name, np.round(r.x, 2))
...
profit first  [25.  0.  0.  0.  0. 15.]
labour x 20   [ 20.   0. 200.   0.   0.   0.]

Reading the output. The inspector IP gives \(4\) Grade 1 and \(3\) Grade 2 inspectors at Rs. 479.2, the answer found by enumeration in Section 1. In the goal programme the two goals conflict. With profit given priority, \(x_1 = 25\) meets the profit goal exactly and overruns labour by \(15\) hours (\(d_2^+ = 15\)). With an hour weighted as Rs. 20 of profit, \(x_1 = 20\) meets the labour goal and falls Rs. 200 short of the profit goal. The model does not choose between them; the weights, set by the decision maker, do.

MATLAB/Octave: intlinprog and linprog with the same matrices.