an integer linear programming problem from real life

planks.py

from scipy.optimize import milp, LinearConstraint
import numpy as np
from math import inf

costs = [1270, 1690, 2190]
lengths = [1500, 2000, 3000]
num_types = len(costs)
min_margin = 3

cut_amounts = [(710, 4), (1130, 4), (650, 4), (680, 2), (1190, 6)]
cuts = [length for length, amount in cut_amounts for _ in range(amount)]
num_cuts = len(cuts)

num_planks = 0
optimum_cost = None

while optimum_cost is None or optimum_cost > min(*costs) * num_planks:
    num_planks += 1

    type_vars = [list(range(i * num_types, (i+1) * num_types)) for i in range(num_planks)]
    z = num_planks * num_types
    asgn_vars = [list(range(z + i * num_planks, z + (i+1) * num_planks)) for i in range(num_cuts)]
    z = num_planks * num_types + num_cuts * num_planks

    A = np.zeros([0, z])
    low = np.zeros(0)
    upp = np.zeros(0)

    def row(d, l, u):
        global A, low, upp
        r = np.zeros(z)
        for k in d: r[k] = d[k]
        A = np.vstack([A, r])
        low = np.append(low, l)
        upp = np.append(upp, u)

    for i in range(num_planks): row({ ty: 1 for ty in type_vars[i] }, 1, 1)
    for i in range(num_cuts): row({ pl: 1 for pl in asgn_vars[i] }, 1, 1)
    for i in range(num_planks): row(
            { asgn_vars[j][i]: cuts[j]+min_margin for j in range(num_cuts) } | 
            { type_vars[i][j]: -lengths[j]-min_margin for j in range(num_types) },
            -inf,
            0)

    # objective function: the cost of the planks
    c = np.zeros(z)
    for i in range(num_planks):
        for j in range(num_types):
            c[type_vars[i][j]] = costs[j]

    res = milp(c=c, constraints=LinearConstraint(A, low, upp), integrality=np.ones(z))

    if res.success and (optimum_cost is None or res.fun < optimum_cost):
        print("===", res.fun)
        for i in range(num_planks):
            for j in range(num_types):
                if res.x[type_vars[i][j]] > 0.5:
                    print("purchase:", lengths[j])
                    margin = lengths[j]

            for j in range(num_cuts):
                if res.x[asgn_vars[j][i]] > 0.5:
                    print("\tcut:", cuts[j])
                    margin -= cuts[j]

            print("\tmargin:", margin)

        optimum_cost = res.fun