2.4 Geometry in software
Software 2.4.1 (Enumerating vertices). For Example 2.3.1, the intersection of every pair of the five constraint lines (the last two being \(x_1 = 0\) and \(x_2 = 0\)), whether it is feasible, and the objective at the feasible ones; then the LP solution for comparison:
>>> import numpy as np >>> from itertools import combinations >>> from scipy.optimize import linprog >>> G = np.array([[2, 1], [2, 4], [-1, 1], [-1, 0], [0, -1]]) # G x <= h, >>> h = np.array([26, 56, 5, 0, 0]) # the last two rows are x1, x2 >= 0 >>> c = np.array([10, 15]) >>> for i, j in combinations(range(5), 2): # every pair of constraint lines ... M = G[[i, j]] ... if abs(np.linalg.det(M)) < 1e-12: continue # parallel lines ... v = np.linalg.solve(M, h[[i, j]]) + 0.0 ... ok = np.all(G @ v <= h + 1e-9) ... print((i + 1, j + 1), np.round(v, 3), 'vertex' if ok else '-', ... round(c @ v, 2) if ok else '') ... (1, 2) [ 8. 10.] vertex 230.0 (1, 3) [ 7. 12.] - (1, 4) [ 0. 26.] - (1, 5) [13. 0.] vertex 130.0 (2, 3) [ 6. 11.] vertex 225.0 (2, 4) [ 0. 14.] - (2, 5) [28. 0.] - (3, 4) [0. 5.] vertex 75.0 (3, 5) [-5. 0.] - (4, 5) [0. 0.] vertex 0.0 >>> r = linprog(-c, A_ub=G[:3], b_ub=h[:3], method='highs'); r.x, -r.fun (array([ 8., 10.]), 230.0)
Reading the output. Of the ten intersections, five satisfy all the constraints; they are the vertices \(B\), \(A\), \(C\), \(D\), \(O\) of Figure 2. The other five lie outside the region: \((7, 12)\) violates the second constraint, \((28, 0)\) the first. The largest objective value, \(230\) at \((8, 10)\), is what linprog returns.
MATLAB/Octave: nchoosek(1:5, 2) lists the pairs.
Software 2.4.2 (What can go wrong). linprog on an unbounded problem (with \(x_1\) free), an infeasible one, and Example 2.3.1 with the objective \(x_1 + 2x_2\), which has a tie:
>>> cases = { # linprog status: 0 optimal, 2 infeasible, 3 unbounded
... 'unbounded': dict(c=[-1, 4], A_ub=[[-3, 1], [1, 2], [0, 1]],
... b_ub=[6, 4, -3], bounds=[(None, None)]*2),
... 'infeasible': dict(c=[1, 1], A_ub=[[1, 1], [-1, -1]], b_ub=[2, -3]),
... 'tied': dict(c=[-1, -2], A_ub=G[:3], b_ub=h[:3])}
>>> for name, kw in cases.items():
... r = linprog(method='highs', **kw)
... print(f"{name:10s} status {r.status}", r.x if r.status == 0 else '',
... -r.fun if r.status == 0 else '')
...
unbounded status 3
infeasible status 2
tied status 0 [ 6. 11.] 28.0
Reading the output. The solver reports status \(3\) (unbounded) and \(2\) (infeasible) instead of a solution. For the tied objective it returns one optimal vertex, \(C(6, 11)\), with \(Z = 28\); the vertex \(B(8, 10)\) and every point between are equally good, and a solver gives no warning of this, which sensitivity analysis will detect.
MATLAB/Octave: the exitflag output, \(-2\) for infeasible and \(-3\) for unbounded.
Software 2.4.3 (Bases and vertices). Minimise \(20x_1 + 10x_2\) subject to \(x_1 + 2x_2 \leq 40\), \(3x_1 + x_2 \geq 30\), \(4x_1 + 3x_2 \geq 60\), \(x \geq 0\), in standard form with one slack and two surplus variables: all \(\binom 53 = 10\) choices of three basic columns, keeping the feasible ones:
>>> A = np.array([[1, 2, 1, 0, 0], # min 20 x1 + 10 x2, standard form: ... [3, 1, 0, -1, 0], # x1 + 2 x2 + s1 = 40 ... [4, 3, 0, 0, -1]]) # 3 x1 + x2 - s2 = 30, and so on >>> b = np.array([40, 30, 60]); cost = np.array([20, 10, 0, 0, 0]) >>> count = feasible = 0 >>> for B in combinations(range(5), 3): # every choice of three basic columns ... AB = A[:, B] ... if abs(np.linalg.det(AB)) < 1e-12: continue ... count += 1; xB = np.linalg.solve(AB, b) ... if np.all(xB >= -1e-9): ... feasible += 1; x = np.zeros(5); x[list(B)] = xB ... print(B, np.round(x[:2], 3), round(cost @ x, 2)) ... (0, 1, 2) [ 6. 12.] 240.0 (0, 1, 4) [ 4. 18.] 260.0 (0, 2, 3) [15. 0.] 300.0 (0, 3, 4) [40. 0.] 800.0 >>> count, feasible (10, 4)
Reading the output. All ten column triples are bases, and four give basic feasible solutions: the vertices \((6, 12)\), \((4, 18)\), \((15, 0)\) and \((40, 0)\), with costs \(240\), \(260\), \(300\) and \(800\). The minimum is \(240\) at \((6, 12)\), where the two \(\geq \) constraints are tight.
MATLAB/Octave: nchoosek(1:5, 3) and A(:, B) ∖ b.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.