Cutting Stock - A Column Generation Approach (CUTSTOCK)

cutstock.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_cutstock.html
## LICENSETYPE: Demo
## MODELTYPE: MIP, RMIP
## KEYWORDS: mixed integer linear programming, cutting stock, column generation, paper industry


Cutting Stock - A Column Generation Approach (CUTSTOCK)

The task is to cut out some paper products of different sizes from a
large raw paper roll, in order to meet a customer's order. The objective
is to minimize the required number of paper rolls.


P. C. Gilmore and R. E. Gomory, A linear programming approach to the
cutting stock problem, Part I, Operations Research 9 (1961), 849-859.

P. C. Gilmore and R. E. Gomory, A linear programming approach to the
cutting stock problem, Part II, Operations Research 11 (1963), 863-888.
"""

from __future__ import annotations

import os

import gamspy.math as gams_math
from gamspy import (
    Card,
    Container,
    Equation,
    Model,
    Options,
    Ord,
    Parameter,
    Sense,
    Set,
    Sum,
    Variable,
)


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
    )

    # Sets
    i = Set(
        m,
        "i",
        records=[f"w{idx}" for idx in range(1, 5)],
        description="widths",
    )
    p = Set(
        m,
        "p",
        records=[f"p{idx}" for idx in range(1, 1001)],
        description="possible patterns",
    )
    pp = Set(m, "pp", domain=p, description="dynamic subset of p")

    # Parameters
    r = Parameter(m, "r", records=100, description="raw width")
    w = Parameter(
        m,
        "w",
        domain=i,
        records=[["w1", 45], ["w2", 36], ["w3", 31], ["w4", 14]],
        description="width",
    )
    d = Parameter(
        m,
        "d",
        domain=i,
        records=[["w1", 97], ["w2", 610], ["w3", 395], ["w4", 211]],
        description="demand",
    )
    aip = Parameter(
        m,
        "aip",
        domain=[i, p],
        description="number of width i in pattern growing in p",
    )

    # Master model variables
    xp = Variable(
        m, "xp", domain=p, type="integer", description="patterns used"
    )
    z = Variable(m, "z", description="objective variable")
    xp.up[p] = Sum(i, d[i])

    # Master model equations
    numpat = Equation(
        m,
        "numpat",
        definition=z == Sum(pp, xp[pp]),
        description="number of patterns used",
    )
    demand = Equation(m, "demand", domain=i, description="meet demand")
    demand[i] = Sum(pp, aip[i, pp] * xp[pp]) >= d[i]

    master = Model(
        m,
        "master",
        equations=[numpat, demand],
        problem="rmip",
        sense=Sense.MIN,
        objective=z,
    )

    # Pricing model variables
    y = Variable(m, "y", domain=i, type="integer", description="new pattern")
    y.up[i] = gams_math.ceil(r / w[i])

    defobj = Equation(
        m, "defobj", definition=z == (1 - Sum(i, demand.m[i] * y[i]))
    )
    knapsack = Equation(
        m,
        "knapsack",
        description="knapsack constraint",
        definition=Sum(i, w[i] * y[i]) <= r,
    )

    pricing = Model(
        m,
        "pricing",
        equations=[defobj, knapsack],
        problem="mip",
        sense=Sense.MIN,
        objective=z,
    )

    pp[p] = Ord(p) <= Card(i)
    aip[i, pp[p]].where[Ord(i) == Ord(p)] = gams_math.floor(r / w[i])

    pi = Set(m, "pi", domain=p, description="set of the last pattern")
    pi[p] = Ord(p) == Card(pp) + 1

    while len(pp) < len(p):
        master.solve(options=Options(relative_optimality_gap=0))
        pricing.solve(options=Options(relative_optimality_gap=0))

        if z.records["level"].values[0] >= -0.001:
            break

        aip[i, pi] = gams_math.Round(y.l[i])
        pp[pi] = True
        pi[p] = pi[p.lag(1)]

    master.problem = "mip"
    master.solve(options=Options(relative_optimality_gap=0))

    import math

    assert math.isclose(master.objective_value, 453.0000, rel_tol=0.001)

    patrep = Parameter(
        m, "patrep", domain=["*", "*"], description="solution pattern report"
    )
    demrep = Parameter(
        m,
        "demrep",
        domain=["*", "*"],
        description="solution demand supply report",
    )

    patrep["# produced", p] = gams_math.Round(xp.l[p])
    patrep[i, p].where[patrep["# produced", p]] = aip[i, p]
    patrep[i, "total"] = Sum(p, patrep[i, p])
    patrep["# produced", "total"] = Sum(p, patrep["# produced", p])

    demrep[i, "produced"] = Sum(p, patrep[i, p] * patrep["# produced", p])
    demrep[i, "demand"] = d[i]
    demrep[i, "over"] = demrep[i, "produced"] - demrep[i, "demand"]

    print(patrep.records)
    print(demrep.records)


if __name__ == "__main__":
    main()