cutting_stock_column_generation.jl

#  Copyright 2019, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################

# Based on http://doi.org/10.5281/zenodo.3329388

using JuMP, GLPK, SparseArrays, Test

"""
    solve_pricing(dual_demand_satisfaction, maxwidth, widths, rollcost, demand, prices)
This function specifically implements the pricing problem for
`example_cutting_stock`. It takes, as input, the dual solution from the master
problem and the cutting stock instance. It outputs either a new cutting pattern,
or `nothing` if no pattern could improve the current cost.
"""
function solve_pricing(dual_demand_satisfaction, maxwidth, widths, rollcost, demand, prices)
    reduced_costs = dual_demand_satisfaction + prices
    n = length(reduced_costs)

    # The actual pricing model.
    submodel = Model(with_optimizer(GLPK.Optimizer))
    @variable(submodel, xs[1:n] >= 0, Int)
    @constraint(submodel, sum(xs .* widths) <= maxwidth)
    @objective(submodel, Max, sum(xs .* reduced_costs))

    optimize!(submodel)

    # If the net cost of this new pattern is nonnegative, no more patterns to add.
    new_pattern = round.(Int, value.(xs))
    net_cost = rollcost - sum(new_pattern .* (dual_demand_satisfaction .+ prices))

    if net_cost >= 0 # No new pattern to add.
        return nothing
    else
        return new_pattern
    end
end

"""
    example_cutting_stock(maxwidth, widths, rollcost, demand, prices)
Solves the cutting stock problem (sometimes also called the cutting rod
problem) using a column-generation technique.
Intuitively, this problem is about cutting large rolls of paper into smaller
pieces. There is an exact demand of pieces to meet, and all rolls have the
same size. The goal is to meet the demand while maximising the profits (each
paper roll has a fixed cost, each sold piece allows earning some money),
which is roughly equivalent to using the smallest amount of rolls
to cut (or, equivalently, to minimise the amount of paper waste).
This function takes five parameters:
  * `maxwidth`: the maximum width of a roll (or length of a rod)
  * `widths`: an array of the requested widths
  * `rollcost`: the cost of a complete roll
  * `demand`: the demand, in number of pieces, for each width
  * `prices`: the selling price for each width
Mathematically, this problem might be formulated with two variables:
  * `x[i, j] ∈ ℕ`: the number of times the width `i` is cut out of the roll `j`
  * `y[j] ∈ 𝔹`: whether the roll `j` is used
Several constraints are needed:
  * the demand must be satisfied, for each width `i`:
    ∑j x[i, j] = demand[i]
  * the roll size cannot be exceed, for each roll `j` that is used:
    ∑i x[i, j] width[i] ≤ maxwidth y[j]
If you want to implement this naïve model, you will need an upper bound on the
number of rolls to use: the simplest one is to consider that each required
width is cut from its own roll, i.e. `j` varies from 1 to ∑i demand[i].
This example prefers a more advanced technique to solve this problem:
column generation. It considers a different set of variables: patterns of
width to cut a roll. The decisions then become the number of times each pattern
is used (i.e. the number of rolls that are cut following this pattern).
The intelligence comes from the way these patterns are chosen: not all of them
are considered, but only the "interesting" ones, within the master problem.
A "pricing" problem is used to decide whether a new pattern should be generated
or not (it is implemented in the function `solve_pricing`). "Interesting" means,
for a pattern, that the optimal solution may use this cutting pattern.
In more details, the solving process is the following. First, a series of
dumb patterns is generated (just one width per roll, repeated until the roll is
completely cut). Then, the master problem is solved with these first patterns
and its dual solution is passed on to the pricing problem. The latter decides
if there is a new pattern to include in the formulation or not; if so,
it returns it to the master problem. The master is solved again, the new dual
variables are given to the pricing problem, until there is no more pattern to
generate from the pricing problem: all "interesting" patterns have been
generated, and the master can take its optimal decision. In the implementation,
the variables deciding how many times a pattern is chosen are called `θ`.
For more information on column-generation techniques applied on the cutting
stock problem, you can see:
  * [Integer programming column generation strategies forthe cutting stock
    problem and its variants](https://tel.archives-ouvertes.fr/tel-00011657/document)
  * [Tackling the cutting stock problem](http://doi.org/10.5281/zenodo.3329388)
"""
function example_cutting_stock(; max_gen_cols::Int=5000)
    maxwidth = 100.0
    rollcost = 500.0
    prices = Float64[167.0, 197.0, 281.0, 212.0, 225.0, 111.0, 93.0, 129.0, 108.0, 106.0, 55.0, 85.0, 66.0, 44.0, 47.0, 15.0, 24.0, 13.0, 16.0, 14.0]
    widths = Float64[75.0, 75.0, 75.0, 75.0, 75.0, 53.8, 53.0, 51.0, 50.2, 32.2, 30.8, 29.8, 20.1, 16.2, 14.5, 11.0, 8.6, 8.2, 6.6, 5.1]
    demand = Int[38, 44, 30, 41, 36, 33, 36, 41, 35, 37, 44, 49, 37, 36, 42, 33, 47, 35, 49, 42]
    nwidths = length(prices)

    n = length(widths)
    ncols = length(widths)

    # Initial set of patterns (stored in a sparse matrix: a pattern won't include many different cuts).
    patterns = spzeros(UInt16, n, ncols)
    for i in 1:n
        patterns[i, i] = min(floor(Int, maxwidth / widths[i]), round(Int, demand[i]))
    end

    # Write the master problem with this "reduced" set of patterns.
    # Not yet integer variables: otherwise, the dual values may make no sense
    # (actually, GLPK will yell at you if you're trying to get duals for integer problems)
    m = Model(with_optimizer(GLPK.Optimizer))
    @variable(m, θ[1:ncols] >= 0)
    @objective(m, Min,
        sum(θ[p] * (rollcost - sum(patterns[j, p] * prices[j] for j in 1:n)) for p in 1:ncols)
    )
    @constraint(m, demand_satisfaction[j=1:n], sum(patterns[j, p] * θ[p] for p in 1:ncols) >= demand[j])

    # First solve of the master problem.
    optimize!(m)
    if termination_status(m) != MOI.OPTIMAL
        warn("Master not optimal ($ncols patterns so far)")
    end

    # Then, generate new patterns, based on the dual information.
    while ncols - n <= max_gen_cols # Generate at most max_gen_cols columns.
        if ! has_duals(m)
            break
        end

        new_pattern = solve_pricing(dual.(demand_satisfaction), maxwidth, widths, rollcost, demand, prices)

        # No new pattern to add to the formulation: done!
        if new_pattern === nothing
            break
        end

        # Otherwise, add the new pattern to the master problem, recompute the duals,
        # and go on waltzing one more time with the pricing problem.
        ncols += 1
        patterns = hcat(patterns, new_pattern)

        # One new variable.
        new_var = @variable(m, [ncols], base_name="θ", lower_bound=0)
        push!(θ, new_var[ncols])

        # Update the objective function.
        set_objective_coefficient(m, θ[ncols], rollcost - sum(patterns[j, ncols] * prices[j] for j=1:n))

        # Update the constraint number j if the new pattern impacts this production.
        for j in 1:n
            if new_pattern[j] > 0
                set_normalized_coefficient(demand_satisfaction[j], new_var[ncols], new_pattern[j])
            end
        end

        # Solve the new master problem to update the dual variables.
        optimize!(m)
        if termination_status(m) != MOI.OPTIMAL
            warn("Master not optimal ($ncols patterns so far)")
        end
    end

    # Finally, impose the master variables to be integer and resolve.
    # To be exact, at each node in the branch-and-bound tree, we would need to restart
    # the column generation process (just in case a new column would be interesting
    # to add). This way, we only get an upper bound (a feasible solution).
    set_integer.(θ)
    optimize!(m)

    if termination_status(m) != MOI.OPTIMAL
        warn("Final master not optimal ($ncols patterns)")
    end

    @test JuMP.objective_value(m) ≈ 78374.0 atol = 1e-3
end

example_cutting_stock()