Python lab: model, solve, verify

1. Set up the notebook

Run cells in order with Shift+Enter. Python 3 with NumPy, SciPy ≥1.9, and Matplotlib is required. If an import fails, install these packages in your selected environment, then restart the kernel.

import numpy as np
import scipy
from scipy.optimize import milp, Bounds, LinearConstraint
print("SciPy", scipy.__version__)

SciPy 1.18.1

2. Define the objective and variable order

Use v=[a,b,yA,yB]. Money is in $1,000/week. Why do the opening costs have positive signs in c?

c = np.array([-8., -6., 10., 8.])
integrality = [0, 0, 1, 1]
bounds = Bounds([0,0,0,0], [6,8,1,1])
print("Minimize negative net profit; only yA and yB are binary.")

Minimize negative net profit; only yA and yB are binary.

3. Build each row of the matrix

Read every row aloud: labor, sales upper bound, A linking, B linking. Sales is an upper limit, not a minimum demand.

A = np.array([[2,1,0,0],
              [1,1,0,0],
              [1,0,-6,0],
              [0,1,0,-8]], dtype=float)
upper = np.array([12,10,0,0], dtype=float)
constraints = LinearConstraint(A, -np.inf, upper)
print(A)

[[ 2.  1.  0.  0.]
 [ 1.  1.  0.  0.]
 [ 1.  0. -6.  0.]
 [ 0.  1.  0. -8.]]

4. Solve and interpret status

Predict which lines open before running. A returned vector alone is not an optimality certificate.

result = milp(c, integrality=integrality, bounds=bounds,
              constraints=constraints, options={"mip_rel_gap":0})
print("Status:", result.status, result.message)
assert result.success
print("[a,b,yA,yB] =", result.x)
print("Net profit =", -result.fun)
print("Profit upper bound =", -result.mip_dual_bound)
print("MIP gap =", result.mip_gap)
assert np.isclose(-result.fun,46)

Status: 0 Optimization terminated successfully. (HiGHS Status 7: Optimal)
[a,b,yA,yB] = [2. 8. 1. 1.]
Net profit = 46.0
Profit upper bound = 46.0
MIP gap = 0.0

5. Audit the result independently

Check the four rows, bounds, binary decisions, and original profit expression. Keep tolerances; never round continuous production automatically.

v = result.x
tol = 1e-7
print("Row use:", A @ v)
print("Row slack:", upper - A @ v)
assert np.all(A @ v <= upper + tol)
assert np.all(v >= bounds.lb-tol) and np.all(v <= bounds.ub+tol)
assert np.allclose(v[2:], np.round(v[2:]), atol=tol)
profit = 8*v[0]+6*v[1]-10*v[2]-8*v[3]
assert np.isclose(profit,-result.fun)
print("All checks passed; net profit =",profit)

Row use: [12. 10. -4.  0.]
Row slack: [0. 0. 4. 0.]
All checks passed; net profit = 46.0

6. Enumerate opening choices

Verify the answer with four continuous LPs, one for each fixed pair of opening decisions.

from scipy.optimize import linprog
from itertools import product
records = []
for ya,yb in product([0,1], repeat=2):
    r = linprog([-8,-6], A_ub=[[2,1],[1,1]], b_ub=[12,10],
                bounds=[(0,6*ya),(0,8*yb)], method="highs")
    assert r.success
    z = -r.fun-10*ya-8*yb
    records.append(z)
    print((ya,yb), "production",r.x,"net profit",z)
assert np.isclose(max(records),46)

(0, 0) production [0. 0.] net profit -0.0
(0, 1) production [0. 8.] net profit 40.0
(1, 0) production [6. 0.] net profit 38.0
(1, 1) production [2. 8.] net profit 46.0

7. Relax integrality

Why can a relaxed opening fraction pay only part of a fixed cost? Explain why this value is an upper bound for maximization.

relaxed = milp(c, integrality=[0,0,0,0], bounds=bounds,
               constraints=constraints)
assert relaxed.success
print("Relaxed decisions:",relaxed.x)
print("LP upper bound:",-relaxed.fun)
assert np.isclose(-relaxed.fun,158/3)
assert -relaxed.fun >= -result.fun - tol

Relaxed decisions: [2.         8.         0.33333333 1.        ]
LP upper bound: 52.666666666666664

8. Change one assumption at a time

Reset to the baseline before each scenario. Try sales cap 8, B opening cost 20, or minimum A batch 3. Predict first.

def solve_case(sales=10, fixed_b=8, minimum_a=0, one_line=False):
    rows=A.copy(); rhs=upper.copy(); rhs[1]=sales
    if minimum_a:
        rows=np.vstack([rows,[-1,0,minimum_a,0]])
        rhs=np.r_[rhs,0]
    if one_line:
        rows=np.vstack([rows,[0,0,1,1]])
        rhs=np.r_[rhs,1]
    r=milp([-8,-6,10,fixed_b],integrality=integrality,bounds=bounds,
           constraints=LinearConstraint(rows,-np.inf,rhs),
           options={"mip_rel_gap":0})
    assert r.success
    return r
for label,kwargs,expected in [("Sales cap 8",{"sales":8},40),
                              ("B cost 20",{"fixed_b":20},38),
                              ("A minimum 3",{"minimum_a":3},42)]:
    r=solve_case(**kwargs)
    print(label,r.x,-r.fun)
    assert np.isclose(-r.fun,expected)

Sales cap 8 [0. 8. 0. 1.] 40.0
B cost 20 [6. 0. 1. 0.] 38.0
A minimum 3 [3. 6. 1. 1.] 42.0

9. Add a purchasable capacity module

Add binary k. It supplies 4 labor units and costs p per week. Rewrite the labor row as 2a+b-4k ≤12. Do not add capacity for free.

def module_case(price=7, force=None):
    rows=np.column_stack([A,np.zeros(4)])
    rows[0,-1]=-4
    lb=[0,0,0,0,0]; ub=[6,8,1,1,1]
    if force is not None: lb[-1]=ub[-1]=force
    r=milp([-8,-6,10,8,price],integrality=[0,0,1,1,1],
           bounds=Bounds(lb,ub),constraints=LinearConstraint(rows,-np.inf,upper),
           options={"mip_rel_gap":0})
    assert r.success
    return r
for price in [7,8,9,10]:
    r=module_case(price)
    no=module_case(price,0); buy=module_case(price,1)
    print("Price",price,"best",-r.fun,"no buy",-no.fun,"buy",-buy.fun)
    assert np.isclose(-r.fun,max(-no.fun,-buy.fun))
assert np.isclose(-module_case(7).fun,47)
assert np.isclose(-module_case(9).fun,46)

Price 7 best 47.0 no buy 46.0 buy 47.0
Price 8 best 46.0 no buy 46.0 buy 46.0
Price 9 best 46.0 no buy 46.0 buy 45.0
Price 10 best 46.0 no buy 46.0 buy 44.0

10. Complete independent exercises

Use the separate exercise sheet first. Record your model change, plan, objective, and feasibility checks. Then submit the Exit Ticket before reading its solutions.

print("Next: exercises.html and the online Exit Ticket.")

Next: exercises.html and the online Exit Ticket.

Python code walkthrough

A narrated, step-by-step code walkthrough with the verified expected results beside each cell. Pause and run the full cells in the guided notebook. Approximately 5 minutes, using the same narration as chapter 4.

Download the code walkthrough MP4 · Open the complete runnable lab