Objectives, Decision Variables, Constraints, and Feasibility › Implementation in SciPy and NumPy

# Decision variables: pallets of product 1 and product 2.
A = np.array(
    [
        [2.0, 1.0],  # machine-hours per pallet
        [1.0, 2.0],
    ]
)  # labour-hours per pallet
b = np.array([8.0, 8.0])  # available machine-hours, labour-hours
c = np.array([40.0, 30.0])  # profit per pallet (currency units)

enumerated = []
for point in itertools.product(range(5), repeat=2):
    x = np.array(point, dtype=float)
    slack = b - A @ x
    if np.all(slack >= -1e-12):
        enumerated.append(
            (int(x[0]), int(x[1]), float(c @ x), float(slack.min()))
        )

best_enum = max(enumerated, key=lambda row: row[2])
n_feasible_enum = len(enumerated)

lp = linprog(c=-c, A_ub=A, b_ub=b, bounds=[(0.0, None)] * 2, method="highs")
if not lp.success or lp.x is None or lp.fun is None:
    raise RuntimeError(f"LP failed: {lp.message}")
lp_profit = float(-lp.fun)

milp_fit = milp(
    c=-c,
    constraints=LinearConstraint(A, -np.inf, b),
    integrality=np.ones(2),
    bounds=Bounds(np.zeros(2), np.full(2, np.inf)),
    options={"mip_rel_gap": 0.0},
)

from book_plot_style import PALETTE, label_bars, polish_axes, use_book_style

use_book_style()
C = [PALETTE[key] for key in ("primary", "teal", "amber", "purple")]

fig = plt.figure(figsize=(6.0, 4.8))
ax = fig.add_subplot(111)
positions = [0, 1]
bars = ax.barh(
    positions,
    gap_values,
    color=[PALETTE["primary"], PALETTE["amber"]],
    height=0.5,
)
ax.set_yticks(positions, gap_labels)
ax.set_xlim(0.0, max(gap_values) * 1.18)
ax.set_ylim(-0.6, 1.8)
ax.set_xlabel("Profit (currency units)")
ax.set_title("Relaxation bound and integer incumbent")
label_bars(ax, bars, labels=[f"{value:.2f}" for value in gap_values])
ax.text(
    0.02,
    0.96,
    f"LP-relaxation gap: {gap_relative:.2%}\nFinal MILP gap: {milp_fit.mip_gap:.2%}",
    transform=ax.transAxes,
    va="top",
    color=PALETTE["ink"],
    bbox=dict(facecolor=PALETTE["paper"], alpha=0.92, edgecolor="none"),
)

polish_axes(ax)
plt.show()

def audit(x, tol=1e-9):
    """Independent check of a claimed solution: feasibility first, then value."""
    x = np.asarray(x, dtype=float)
    if x.shape != c.shape or not np.all(np.isfinite(x)):
        raise ValueError(
            "Expected a finite decision vector with one value per product"
        )
    slack = b - A @ x
    return {
        "x": x,
        "slack": slack,
        "overspend": float(np.sum(np.maximum(0.0, -slack))),
        "max_resource_violation": float(np.max(np.maximum(0.0, -slack))),
        "profit": float(c @ x),
        "integer": bool(np.all(np.abs(x - np.round(x)) <= tol)),
        "feasible": bool(np.all(x >= -tol) and np.all(slack >= -tol)),
    }

from book_plot_style import PALETTE, label_bars, polish_axes, use_book_style

use_book_style()
C = [PALETTE[key] for key in ("primary", "teal", "amber", "purple")]

fig = plt.figure(figsize=(6.0, 4.8))
ax = fig.add_subplot(111)
machine_bars = ax.bar(
    slack_positions - slack_width / 2.0,
    slack_matrix[:, 0],
    width=slack_width,
    color=PALETTE["primary"],
    label="Machine slack",
)
labour_bars = ax.bar(
    slack_positions + slack_width / 2.0,
    slack_matrix[:, 1],
    width=slack_width,
    color=PALETTE["amber"],
    label="Labour slack",
)
label_bars(
    ax, machine_bars, labels=[f"{value:.0f}" for value in slack_matrix[:, 0]]
)
label_bars(
    ax, labour_bars, labels=[f"{value:.0f}" for value in slack_matrix[:, 1]]
)
ax.set_ylim(-0.05, 1.18)
ax.set_xticks(slack_positions, audit_names)
ax.set_ylabel("Slack (hours)")
ax.set_title("Independent slack audit across solution methods")
ax.axhline(0.0, color=PALETTE["muted"], lw=1.0)
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.20), ncol=2)
polish_axes(ax)
plt.show()

def penalised(lam):
    """Negative profit plus a quadratic penalty on resource overuse."""

    def objective(x):
        over = np.maximum(0.0, A @ x - b)
        return float(-(c @ x) + lam * np.sum(over**2))

    return objective


penalty_rows = []
for lam in (0.5, 5.0, 50.0, 500.0):
    fit = minimize(
        penalised(lam),
        x0=np.array([1.0, 1.0]),
        method="L-BFGS-B",
        bounds=[(0.0, None)] * 2,
        options={"gtol": 1e-10, "maxiter": 10000},
    )
    if not fit.success or fit.x is None:
        raise RuntimeError(f"Penalty solve failed: {fit.message}")
    report = audit(fit.x)
    penalty_rows.append(
        (lam, fit.x[0], fit.x[1], report["overspend"], report["profit"])
    )