A Primal-Dual Global Optimization approach on the Cutting Stock Problem

This doc reads a bit lengthy, because it uses complex data structures (maybe for visualization purposes).
Therein it uses economic interpretation, which is hard to understand for me.
Therefore, I explicitly wrote the Lagrangian dual of the primal MIP, which makes the column generation subprocedure more intelligible.

Comments are welcome.

import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.I as Identity

w, d, W, I = let # 📘 Data for the Cutting Stock Problem
    w = [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]; # length (in meters)
    d = [189, 33, 36, 41, 35, 37, 44, 49, 37, 36, 42, 33, 47, 35, 49, 42]; # demand (in quantities)
    W = 100; # width of a raw roll
    I = length(d); @assert length(w) == I
    w, d, W, I
end;

# 🟣 Part I
# The following is a primary large-scale deterministic MIP formulation of the Cutting Stock Problem.
# The chief advantage ✅ of this problem formulation is that it can bring a global optimal solution,
# as long as you have enough time for computing.
function deterministic_mip(J) # J is a prescribed (large) number,
    # if the resulting termination status is INFEASIBLE, enlarge J.
    det = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.@variable(det, b[1:J], Bin) # b[j] == 1 ⟺ the j-th roll is used
    JuMP.@variable(det, n[1:I, 1:J] ≥ 0, Int) # n[i, j] is the number of (length == w[i])-subrolls cut from the j-th roll
    JuMP.@constraint(det, sum(n; dims = 2) .≥ d) # meet the demand
    JuMP.@constraint(det, w'n .≤ b'W) # width limit of each raw roll
    JuMP.@objective(det, Min, sum(b)) # number of used raw rolls
    JuMP.optimize!(det); JuMP.assert_is_solved_and_feasible(det; allow_local = false)
    num = JuMP.objective_value(det)
    println("The global optimal value is $(num)")
end
function deterministic_mip(b0, n0) # This is merely a variant of the above formulation in which you specify initial solutions
    J = length(b0) # 🟠
    det = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.@variable(det, b[1:J], Bin)
    JuMP.@variable(det, n[1:I, 1:J] ≥ 0, Int)
    JuMP.@constraint(det, sum(n; dims = 2) .≥ d)
    JuMP.@constraint(det, w'n .≤ b'W)
    JuMP.@objective(det, Min, sum(b))
    JuMP.set_start_value.(b, b0) # 🟠
    JuMP.set_start_value.(n, n0) # 🟠
    JuMP.optimize!(det); JuMP.assert_is_solved_and_feasible(det; allow_local = false)
    num = JuMP.objective_value(det)
    println("The global optimal value is $(num)")
end
# Do a blind (the number 500 is a trial) and brute-force (entirely by the MIP solver) solve
deterministic_mip(500); # 🖥️: The global optimal value is 334.0

# 🟣 Part II
# We introduce a pattern perspective of the Cutting Stock Problem and write an iterative algorithm---Column Generation.
# Each pattern is encoded into a column of the matrix `P` below.
# We justify the action "recruiting a new column" through the lens of Dual Cutting Plane method.
P = Identity(I); # ⚠️ This initial columns should make the initial `outer_mip` feasible
COT = 1e-6; # cut-off tolerance---a hyper parameter
function outer_mip(P) # employed to find feasible/suboptimal solutions
    J = size(P, 2)
    outer = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(outer)
    JuMP.@variable(outer, n[1:J] ≥ 0, Int)
    JuMP.@constraint(outer, [i = 1:I], view(P, i, :)'n ≥ d[i]) # [π ≥ 0] 🟡
    JuMP.@objective(outer, Min, sum(n))
    JuMP.optimize!(outer); JuMP.assert_is_solved_and_feasible(outer; allow_local = false)
    return ub, n = JuMP.objective_value(outer), JuMP.value.(n)
end;
# establish in advance (part of) the pricing subproblem
subproblem = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(subproblem);
JuMP.@variable(subproblem, p[1:I] ≥ 0, Int); JuMP.@constraint(subproblem, w'p ≤ W); # a feasible pattern
# establish in advance (part of) Lagrangian dual problem of the primal `outer_mip` above
outerDual = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(outerDual);
JuMP.@variable(outerDual, π[1:I] ≥ 0); JuMP.@objective(outerDual, Max, d'π);
J = size(P, 2);
JuMP.@constraint(outerDual, [j = 1:J], 1 ≥ view(P, :, j)'π); # [⇐ n ≥ 0] 🟡
while true # the aim of this main loop is to recruit as many advantageous columns as possible,
    # in order to lower the primal objective value of the `outer_mip`.
    # This goal is in turn pursued by attempting to lower the Lagrangian dual bound of `outer_mip`.
    # Hence, a Dual Cutting Plane method (i.e. the `@constraint` below) is employed
    Π, lb = let
        JuMP.optimize!(outerDual); JuMP.assert_is_solved_and_feasible(outerDual; allow_local = false)
        JuMP.value.(π), JuMP.objective_value(outerDual)
    end
    p_new, vio_level = let
        JuMP.@objective(subproblem, Max, Π'p)
        JuMP.optimize!(subproblem); JuMP.assert_is_solved_and_feasible(subproblem; allow_local = false)
        JuMP.value.(p), JuMP.objective_value(subproblem) - 1
    end
    if vio_level > COT # ✅ Note that this `if` is a hallmark of the Cutting Plane method
        @info "At J = $J, (invalid) lb = $lb"
        JuMP.@constraint(outerDual, 1 ≥ p_new'π) # add a new cut to lower the lb of the outer_mip
        global P = [P p_new] # recruit a new advantageous column
        global J = size(P, 2) # update J
    else
        @info "At J = $J, (invalid) lb = $lb\nno more advantageous columns can be found, thus break"
        break
    end
end # during the execution of this loop, we observe that `lb` is indeed iteratively lowered.
# Note that `size(P, 2)` is enlarged compared to that initial Identity state.
# Next We generate a feasible solution (in the pattern-formulation) to the primal problem.
# This feasible solution cannot assure global optimality because we have no idea
# whether there is still an advantageous column unrecruited. And I'm unaware of
# an efficient algorithm to identify it without branching.
ub, n = outer_mip(P)

# 🟣 Part III
# The feasible solution found in Part II is suboptimal. But we can employ it as a starting
# solution to be delivered to the MIP solver, so that we end up with a global solution. 
# To this end, we need to convert it to (b0, n0).
j_vec = findall(n -> n > .5, n) # used columns indices
n, P = n[j_vec], P[:, j_vec] # used part
n = round.(Int, n)
P = round.(Int, P)
n0 = hcat([repeat(view(P, :, j), 1, n[j]) for j in eachindex(n)]...)
b0 = trues(size(n0, 2))
@assert all(sum(n0; dims = 2) .≥ d) # do a demand feasibility check
# you can compare this "prepared" solve with that blind one in Part I
deterministic_mip(b0, n0) # 🖥️: The global optimal value is 334.0
2 Likes

If you think part of the documentation of a package is too lengthy or can be improved, why not open an issue or PR?

1 Like

For 2 reasons, maybe there are some suggestions I can hear. And my study about this problem is not finished yet. I’ll try if I can discover more insights.
The cutting pattern formulation, although looks appealing, does not admit a deterministic formulation, thereby global optimality is not assured. In contrast, the first formulation, despite looks a bit “silly”, it is indeed a deterministic MIP that can be guaranteed to be solved to global optimality, provided enough time to compute.

The suboptimal solution derived from column generation can serve as an starting feasible solution that can be fed to the primary roll-wise large-scale deterministic formulation—added in my “Part III”. Therefore we obviate the need for branch-and-price, yet still end up with a global optimization. You can take a look when you have time @odow

Hey,

At each iteration of column generation, you look at the column that has best reduced cost.

In the original formulation of the tutorial, all rolls are identical. You therefore just need to solve one knapsack at each iteration.

Make sure the master can pick any number of columns (patterns)

Focusing on subproblem for roll j=1, given \pi_i the dual value of demand covering constraint of item I, to generate the column with most negative reduced cost, you optimize :

\min_{x,y} y_1 - \sum_i \pi_i x_{i1} s.t. knapsack constraint

(Penalize demand convering contraint violation by \pi on the original formulation and you’ll end up with same problem as above if you focus only on roll j=1).

Assume y_1 = 1 because a master solution containing 0 pattern is feasible.

It yields: \min_{x} 1 - \sum_i \pi_i x_{i1} s.t. knapsack constraint
If reduced cost is negative, you keep the pattern because it improves the LP restricted master.

and this is equivalent to

\max_{x} \sum_i \pi_i x_{i1} s.t. knapsack constraint.
If reduced cost is > 1, keep the pattern.

Hope it helps

1 Like

To be honest, I cannot fully get the name “reduced cost”, along with “shadow price”. I think it should be enough for people to fathom what a “dual variable” is. I can only imagine it being a quantity like “rate of change” (for a scalar constraint) and a “gradient” (for a vector constraint).
Another point in MIP is, Lagrangian relaxation is always tighter than natural LP relaxation.
The JuMP doc uses “do LP relaxation and get dual”. In contrast, I think my Lagrangian dual method is better and can be generalized to more complex problems.

why is a master solution containing 0 pattern is feasible?
In my code above, my initial pattern is simply an identity matrix (each type of demand is considered independently), which saves effort.

Whatever, if finally we reach the same formulation, then both points of view would be nice.


Several places that the JuMP’s doc can be improved, in my opinion, is

  i   w_i d_i
 ------------
  1  75.0  38
  2  75.0  44
  3  75.0  30
  4  75.0  41
  5  75.0  36
  6  53.8  33

The i in 1:5 can be merged to belong to the same i, otherwise it looks weird.

Second,

function solve_pricing(data::Data, π::Vector{Float64})
..............
    return nothing
end

The subproblem is enclosed in a function, meaning that each iteration we’ll rebuild this JuMP model. In my code, this do not happen. Therefore I saved this overhead.

Third, I completely replaced this subsection " Choosing the initial set of patterns" with an Identity matrix, thus made the exposition shorter and smoother.

Here is a code in which I adopted an alternative algorithm for the same problem.
Note that this post is independent of my initial post.
And it can also find the global optimal.
Besided, I generated a larger scale test instance.

A primary version
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.⋅ as ⋅
w, d, W, I = let # 📘 Data for the Cutting Stock Problem
    w = [59, 61, 62, 63, 65, 66, 67, 70, 72, 74, 77, 78, 79, 81, 83, 84, 86, 88, 89, 91, 92, 103, 104, 105, 106, 119, 120, 122, 125, 131, 133, 134, 136, 142, 144, 145, 148, 151, 154, 157, 159, 161, 166, 167, 168, 170, 174, 175, 176, 178, 180, 181]
    d = [77, 48, 62, 64, 66, 49, 25, 21, 56, 47, 53, 51, 23, 73, 71, 22, 23, 68, 71, 51, 66, 55, 58, 76, 71, 31, 53, 46, 57, 56, 63, 77, 53, 65, 28, 70, 72, 35, 46, 23, 34, 39, 36, 68, 31, 64, 40, 22, 36, 37, 39, 57]
    W = 223; # width of a raw roll
    I = length(d); @assert length(w) == I
    w, d, W, I
end;

J, LagIniBnd = 2000, 2100 # 📗 hyper parameter of the algorithm devised in this file

Mip = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
JuMP.@variable(Mip, 0 ≤ b[1:J] ≤ 1, Int);    # consider b[j] ∀ j
JuMP.@variable(Mip, 0 ≤ n[1:I, 1:J], Int);   # consider n[1:I, j] ∀ j
JuMP.@constraint(Mip, [j = 1:J], w'view(n, :, j) ≤ b[j]W); # ∀ j
JuMP.@objective(Mip, Min, sum(b[j] for j in 1:J)); # for each j
JuMP.@constraint(Mip, [i = 1:I], sum(n[i, j] for j in 1:J) ≥ d[i]); # [π ≥ 0] 🟡
# 💡: you can try a direct solve at this line, and see in which point it
# is inferior, compared with that final solve at the end of this code file

Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Min_j); # the j-th block
JuMP.@variable(Min_j, 0 ≤ b_j ≤ 1, Int);
JuMP.@variable(Min_j, 0 ≤ n_j[1:I], Int);
JuMP.@constraint(Min_j, w'n_j ≤ (b_j)W);

Lag = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Lag); # Lagrangian dual of `Mip`
JuMP.@variable(Lag, 0 ≤ π[1:I]);
JuMP.@variable(Lag, θ); # All J blocks are homogeneous, therefore only one subblock is needed
JuMP.@variable(Lag, Λ ≤ LagIniBnd); # ⚠️ add an large artificial bound
JuMP.@constraint(Lag, Λ ≤ d'π + (J)θ);
JuMP.@objective(Lag, Max, Λ);

function get_best_dual()
    for Lag_ite in 1:typemax(Int16)
        JuMP.optimize!(Lag); JuMP.assert_is_solved_and_feasible(Lag; allow_local = false)
        Π, ub = JuMP.value.(π), JuMP.objective_bound(Lag)
        JuMP.@objective(Min_j, Min, b_j - n_j'Π)
        JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
        bj, nj, Θ = JuMP.value(b_j), JuMP.value.(n_j), JuMP.objective_value(Min_j)
        lb = d'Π + (J)Θ
        @info "Lag_ite = $Lag_ite ▶ lb = $lb | $ub = ub"
        if lb + 1e-5 > ub
            BndDual = JuMP.dual(JuMP.UpperBoundRef(Λ))
            isapprox(BndDual, 0; atol=1e-7) || error("BndDual is $BndDual")
            return Π
        end
        JuMP.@constraint(Lag, θ ≤ bj - nj'π);
    end
end;

let # strengthen the original Mip
    Π, rhs = let
        Π = get_best_dual()
        JuMP.@objective(Min_j, Min, b_j - n_j'Π)
        JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
        rhs = JuMP.objective_value(Min_j)
        Π, rhs
    end
    JuMP.@constraint(Mip, [j = 1:J], b[j] - view(n, :, j)'Π ≥ rhs) # 🟣 add DWB cut
end

JuMP.optimize!(Mip); JuMP.solution_summary(Mip) # solve the DWB-equipped orginal Mip to global optimality

# 🗒️ Test Result: 
# 1. with Lagrangian bound:     Explored 1 nodes (282897 simplex iterations) in 108.98 seconds (104.16 work units), Find the OPTIMAL solution
# 2. without Lagrangian bound:  After Exploring 1245 nodes (1769432 simplex iterations) in 299.92 seconds (349.91 work units), ObjVal = 1.47000e+03, OBjBnd = 1.44800e+03

Edit on 2025/6/16

Norm-2 Regularization
J, LagIniBnd = 1531, 1530 # 📗 hyper parameter of the algorithm devised in this file

Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Min_j); # the j-th block
JuMP.@variable(Min_j, b_j, Bin);
JuMP.@variable(Min_j, 0 ≤ n_j[1:I], Int);
JuMP.@constraint(Min_j, w'n_j ≤ (b_j)W);

Lag = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Lag); # Lagrangian dual of `Mip`
JuMP.@variable(Lag, 0 ≤ π[i = 1:I]);
JuMP.@variable(Lag, θ); # All J blocks are homogeneous, therefore only one subblock is needed
JuMP.@variable(Lag, Λ ≤ LagIniBnd); # ⚠️ add an large artificial bound
JuMP.@expression(Lag, LagObjBias, d'π)
JuMP.@constraint(Lag, Λ ≤ LagObjBias + (J)θ);
JuMP.@variable(Lag, nm2_dist); # used in the regularization step

lb, ub = 0., LagIniBnd # valid finite initial bounds
Π = falses(I) # (trial) dual solution starts from zero
t0 = time()
for Lag_ite in 1:typemax(Int)
    # outer layer
    JuMP.@objective(Lag, Max, Λ); JuMP.optimize!(Lag); JuMP.assert_is_solved_and_feasible(Lag; allow_local = false)
    ub = JuMP.objective_bound(Lag)
    good_objc = JuMP.@constraint(Lag, .7 * ub + .3 * lb ≤ Λ)
    quadc = JuMP.@constraint(Lag, [nm2_dist; π - Π] ∈ JuMP.SecondOrderCone())
    JuMP.@objective(Lag, Min, nm2_dist); JuMP.optimize!(Lag); JuMP.assert_is_solved_and_feasible(Lag; allow_local = false)
    Π = JuMP.value.(π) # Currently we take this as our trial solution
    common_value = JuMP.value(LagObjBias)
    JuMP.delete(Lag, good_objc); JuMP.delete(Lag, quadc)
    ### subproblems
    c_b_j = 1.0
    c_n_j = [-Π[i] for i = 1:I]
    JuMP.@objective(Min_j, Min, +(c_b_j ⋅ b_j, c_n_j ⋅ n_j))
    JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
    # add cut to Lag
    ec_b_j = 1.0
    ec_n_j = [-π[i] for i = 1:I]
    JuMP.@constraint(Lag, θ ≤ +(ec_b_j ⋅ JuMP.value(b_j), ec_n_j ⋅ JuMP.value.(n_j)))
    ### lb | ub
    lb_new = common_value + (J)JuMP.objective_value(Min_j)
    lb = max(lb, lb_new)
    @info "Lag_ite = $Lag_ite ▶ lb = $lb | $ub = ub"
    if lb + 1e-5 > ub
        @info "LagDualOpt converges"
        break
    end
end
println("lb = $lb | $ub = ub") # lb = 1467.0833254250904 | 1467.0833333333337 = ub
println("QuadRegularization: $(time() - t0) seconds") # 0.795
Norm-1 Regularization

J, LagIniBnd = 1531, 1530 # 📗 hyper parameter of the algorithm devised in this file

Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Min_j); # the j-th block
JuMP.@variable(Min_j, b_j, Bin);
JuMP.@variable(Min_j, 0 ≤ n_j[1:I], Int);
JuMP.@constraint(Min_j, w'n_j ≤ (b_j)W);

Lag = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Lag); # Lagrangian dual of `Mip`
JuMP.@variable(Lag, 0 ≤ π[i = 1:I]);
JuMP.@variable(Lag, θ); # All J blocks are homogeneous, therefore only one subblock is needed
JuMP.@variable(Lag, Λ ≤ LagIniBnd); # ⚠️ add an large artificial bound
JuMP.@expression(Lag, LagObjBias, d'π)
JuMP.@constraint(Lag, Λ ≤ LagObjBias + (J)θ);
# used in the norm_1 regularization step
JuMP.@variable(Lag, π_m[i = 1:I])
JuMP.@variable(Lag, π_p[i = 1:I])
JuMP.@constraint(Lag, π_m + 2π == π_p)
JuMP.@expression(Lag, nm1_dist, sum(π_m + π_p))

lb, ub = 0., LagIniBnd # valid finite initial bounds
Π = falses(I) # (trial) dual solution starts from zero
t0 = time()
for Lag_ite in 1:typemax(Int)
    # outer layer
    JuMP.@objective(Lag, Max, Λ); JuMP.optimize!(Lag); JuMP.assert_is_solved_and_feasible(Lag; allow_local = false)
    ub = JuMP.objective_bound(Lag)
    good_objc = JuMP.@constraint(Lag, .7 * ub + .3 * lb ≤ Λ)
    JuMP.set_lower_bound.(π_m, -Π)
    JuMP.set_lower_bound.(π_p, Π)
    JuMP.@objective(Lag, Min, nm1_dist); JuMP.optimize!(Lag); JuMP.assert_is_solved_and_feasible(Lag; allow_local = false)
    Π = JuMP.value.(π) # Currently we take this as our trial solution
    common_value = JuMP.value(LagObjBias)
    JuMP.delete(Lag, good_objc)
    ### subproblems
    c_b_j = 1.0
    c_n_j = [-Π[i] for i = 1:I]
    JuMP.@objective(Min_j, Min, +(c_b_j ⋅ b_j, c_n_j ⋅ n_j))
    JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
    # add cut to Lag
    ec_b_j = 1.0
    ec_n_j = [-π[i] for i = 1:I]
    JuMP.@constraint(Lag, θ ≤ +(ec_b_j ⋅ JuMP.value(b_j), ec_n_j ⋅ JuMP.value.(n_j)))
    ### lb | ub
    lb_new = common_value + (J)JuMP.objective_value(Min_j)
    lb = max(lb, lb_new)
    @info "Lag_ite = $Lag_ite ▶ lb = $lb | $ub = ub"
    if lb + 1e-5 > ub
        @info "LagDualOpt converges"
        break
    end
end
println("lb = $lb | $ub = ub") # lb = 1467.0833287621804 | 1467.0833333333335 = ub
println("Norm_1 Regularization: $(time() - t0) seconds") # 0.73

I think the general theory should not appear in a JuMP tutorial. My opinion is that the goal of JuMP documentation is to show what you can do with this package. In the case of column generation, the tutorial should focus on how to dynamically add variables to a formulation.

Take a look at the Dantzig-Wolfe decomposition. You’ll see that it’s based on a Lagrangian relaxation of the original problem and you need column generation to optimize the reformulation (e.g. the work of Vanderbeck & Wolsey: https://www.math.u-bordeaux.fr/~fvanderb/papers/refDecompIPworkingPaper.pdf part 3),

If you reformulate the original problem of the tutorial with a Dantzig-Wolfe decomposition, get rid of the convex combination constraint in the master and solve only one subproblem, you’ll get exactly the same algorithm with the same master and subproblem as in the tutorial. Reduced costs, dual values and Lagrangian are linked.

Good luck, it’s a nice topic :slight_smile:

2 Likes

Yes, I’ve noticed this point

# Create a new column
push!(x, @variable(model, lower_bound = 0))
# Update the non-zeros in the coefficient matrix
set_normalized_coefficient(demand[i], x[end], count)

To be honest, I’ve been using JuMP for 2 years, yet never used these 2 methods. (neither this time in this topic). (But it also reflects that Cutting Plane is indeed simpler—purely from a coding perspective.)

Thanks for sharing! I’ll have a read.

Joining my code in #1 and #7 posts, I’m glad to declare that I discover a (presumably) the most efficient algorithm so far to solving the cutting stock problem. The main ideas include

  1. We need a deterministic MIP formulation to the Cutting Stock Problem being studied. Such that we can quest a global optimal solution with the help of Gurobi (or other MIP solver).
  2. From primal side, we employ the idea of Column Generation to find a good heuristic solution
  3. From dual side, we enforce the Lagrangian dual bound

My code is available here

Primal-Dual Algorithm for Cutting Stock
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.I as Identity

w, d, W, I = let # 📘 Data for the Cutting Stock Problem
    w = [59, 61, 62, 63, 65, 66, 67, 70, 72, 74, 77, 78, 79, 81, 83, 84, 86, 88, 89, 91, 92, 103, 104, 105, 106, 119, 120, 122, 125, 131, 133, 134, 136, 142, 144, 145, 148, 151, 154, 157, 159, 161, 166, 167, 168, 170, 174, 175, 176, 178, 180, 181]
    d = [77, 48, 62, 64, 66, 49, 25, 21, 56, 47, 53, 51, 23, 73, 71, 22, 23, 68, 71, 51, 66, 55, 58, 76, 71, 31, 53, 46, 57, 56, 63, 77, 53, 65, 28, 70, 72, 35, 46, 23, 34, 39, 36, 68, 31, 64, 40, 22, 36, 37, 39, 57]
    W = 223; # width of a raw roll
    I = length(d); @assert length(w) == I
    w, d, W, I
end;

# 🟣 §I: Primal Heuristic Solution via Column Generation
P = Identity(I); # ⚠️ This initial columns should make the initial `outer_mip` feasible
COT = 1e-6; # cut-off tolerance---a hyper parameter
function outer_mip(P) # employed to find feasible/suboptimal solutions
    J = size(P, 2)
    outer = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(outer)
    JuMP.@variable(outer, n[1:J] ≥ 0, Int)
    JuMP.@constraint(outer, [i = 1:I], view(P, i, :)'n ≥ d[i]) # [π ≥ 0] 🟡
    JuMP.@objective(outer, Min, sum(n))
    JuMP.optimize!(outer); JuMP.assert_is_solved_and_feasible(outer; allow_local = false)
    return ub, n = JuMP.objective_value(outer), JuMP.value.(n)
end;
subproblem = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(subproblem);
JuMP.@variable(subproblem, p[1:I] ≥ 0, Int); JuMP.@constraint(subproblem, w'p ≤ W); # a feasible pattern
outerDual = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(outerDual);
JuMP.@variable(outerDual, π[1:I] ≥ 0); JuMP.@objective(outerDual, Max, d'π);
J = size(P, 2);
JuMP.@constraint(outerDual, [j = 1:J], 1 ≥ view(P, :, j)'π); # [⇐ n ≥ 0] 🟡
while true
    Π, lb = let
        JuMP.optimize!(outerDual); JuMP.assert_is_solved_and_feasible(outerDual; allow_local = false)
        JuMP.value.(π), JuMP.objective_value(outerDual)
    end
    p_new, vio_level = let
        JuMP.@objective(subproblem, Max, Π'p)
        JuMP.optimize!(subproblem); JuMP.assert_is_solved_and_feasible(subproblem; allow_local = false)
        JuMP.value.(p), JuMP.objective_value(subproblem) - 1
    end
    if vio_level > COT # ✅ Note that this `if` is a hallmark of the Cutting Plane method
        @info "At J = $J, (invalid) lb = $lb"
        JuMP.@constraint(outerDual, 1 ≥ p_new'π) # add a new cut to lower the lb of the outer_mip
        global P = [P p_new] # recruit a new advantageous column
        global J = size(P, 2) # update J
    else
        @info "At J = $J, (invalid) lb = $lb\nno more advantageous columns can be found, thus break"
        break
    end
end
ub, n = outer_mip(P)
j_vec = findall(n -> n > .5, n) # used columns indices
n, P = round.(Int, n[j_vec]), round.(Int, P[:, j_vec])
n0 = hcat([repeat(view(P, :, j), 1, n[j]) for j in eachindex(n)]...) # 🟠
b0 = trues(size(n0, 2)) # 🟠
@assert all(sum(n0; dims = 2) .≥ d) # do a demand feasibility check
J = length(b0)
LagIniBnd = 1.1J

# 🟣 §II: Lagrangian Dual Bound 
Mip = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV));
JuMP.@variable(Mip, 0 ≤ b[1:J] ≤ 1, Int);
JuMP.@variable(Mip, 0 ≤ n[1:I, 1:J], Int);
JuMP.set_start_value.(b, b0) # 🟠
JuMP.set_start_value.(n, n0) # 🟠
JuMP.@constraint(Mip, [j = 1:J], w'view(n, :, j) ≤ b[j]W);
JuMP.@objective(Mip, Min, sum(b[j] for j in 1:J));
JuMP.@constraint(Mip, [i = 1:I], sum(n[i, j] for j in 1:J) ≥ d[i]); # [π ≥ 0] 🟡
Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Min_j); # the j-th block
JuMP.@variable(Min_j, 0 ≤ b_j ≤ 1, Int);
JuMP.@variable(Min_j, 0 ≤ n_j[1:I], Int);
JuMP.@constraint(Min_j, w'n_j ≤ (b_j)W);
Lag = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(Lag); # Lagrangian dual of `Mip`
JuMP.@variable(Lag, 0 ≤ π[1:I]);
JuMP.@variable(Lag, θ); # All J blocks are homogeneous, therefore only one subblock is needed
JuMP.@variable(Lag, Λ ≤ LagIniBnd); # ⚠️ add an large artificial bound
JuMP.@constraint(Lag, Λ ≤ d'π + (J)θ);
JuMP.@objective(Lag, Max, Λ);
function get_best_dual()
    for Lag_ite in 1:typemax(Int16)
        JuMP.optimize!(Lag); JuMP.assert_is_solved_and_feasible(Lag; allow_local = false)
        Π, ub = JuMP.value.(π), JuMP.objective_bound(Lag)
        JuMP.@objective(Min_j, Min, b_j - n_j'Π)
        JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
        bj, nj, Θ = JuMP.value(b_j), JuMP.value.(n_j), JuMP.objective_value(Min_j)
        lb = d'Π + (J)Θ
        @info "Lag_ite = $Lag_ite ▶ lb = $lb | $ub = ub"
        if lb + 1e-5 > ub
            BndDual = JuMP.dual(JuMP.UpperBoundRef(Λ))
            isapprox(BndDual, 0; atol=1e-7) || error("BndDual is $BndDual")
            return Π
        end
        JuMP.@constraint(Lag, θ ≤ bj - nj'π);
    end
end;
let # strengthen the original Mip
    Π, rhs = let
        Π = get_best_dual()
        JuMP.@objective(Min_j, Min, b_j - n_j'Π)
        JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
        rhs = JuMP.objective_value(Min_j)
        Π, rhs
    end
    JuMP.@constraint(Mip, [j = 1:J], b[j] - view(n, :, j)'Π ≥ rhs) # 🟣 add DWB cut
end

# 🟣 §III: Solve to Global Optimality with (Primal Solution @§I, Dual Bound @§II)
# This should be observably a lot easier than a direct crude solve
JuMP.optimize!(Mip); JuMP.solution_summary(Mip)

Let’s take a look at how eased Gurobi is

Optimize a model with 2988 rows, 77804 columns and 231944 nonzeros
Model fingerprint: 0xd4ea8bd8
Variable types: 0 continuous, 77804 integer (0 binary)
Coefficient statistics:
  Matrix range     [3e-01, 2e+02]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+01, 8e+01]

Loaded user MIP start with objective 1468

Presolve time: 0.33s
Presolved: 2988 rows, 77804 columns, 231944 nonzeros
Variable types: 0 continuous, 77804 integer (41104 binary)
Deterministic concurrent LP optimizer: primal and dual simplex
Showing primal log only...

Concurrent spin time: 0.00s

Solved with primal simplex

Root relaxation: cutoff, 53526 iterations, 2.29 seconds (2.11 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0     cutoff    0      1468.00000 1468.00000  0.00%     -    2s

Explored 1 nodes (53526 simplex iterations) in 2.79 seconds (2.59 work units)
Thread count was 8 (of 8 available processors)

Solution count 1: 1468

Optimal solution found (tolerance 1.00e-04)
Best objective 1.468000000000e+03, best bound 1.468000000000e+03, gap 0.0000%

I would be more than glad to see if anyone can raise a better algorithm on this problem that beats mine. :grinning_face_with_smiling_eyes:

2 Likes

The previous posts are based on the primary case where we minimize the number of raw rolls.
This is actually a degenerate case which leads to homogeneous blocks—meaning that only one subblock is needed when writting code.

Therefore, I extend the Cutting Stock Problem to a heterogeneous version—multiple different blocks—each being a (MIP) subproblem.

I’ve devised a satisfying (I think it’s pretty cool) algorithm which I’m going to call it a sandwiched MIP method (strengthening the original MIP formulation by sandwiching it).

Roughly speaking, Gurobi doesn’t perform well on the original compact formulation----say, converges to 1% after 2 hours—and become slower and slower due to huge size of branches.
While my algorithm spend nearly 1 hour to train a tight Lagrangian bound. Then after sandwiching, Gurobi returns me the OPTIMAL termination status with a 0.0098% gap.

Currently the main concern in my algorithm is a heuristic employed to identify some promising subproblem block—I’m unaware of if there is a good strategy. (x-ref Find the max element of a numeric vector iteratively with mask and early break)

My implementation is

code
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.⋅ as ⋅
import Random
Random.seed!(321); # ⚠️ when modify this, U should also modify hyperparam

I, J, w, d, W, c = let # 📄 load data
    δ = 1.7e-12 # used in test case generation
    I = 20; # [demand side] number of types of demand
    J = 400; # [resource side] number of raw rolls
    wl, wh = rand(13:δ:35), rand(77:δ:97)
    w = rand(wl:δ:wh, I) # width of each demand type
    dl, dh = rand(7:14), rand(57:77)
    d = rand(dl:dh, I) # quantity::Int of each demand type
    Wl, Wh = rand(70:δ:100), rand(400:δ:500)
    W = rand(Wl:δ:Wh, J) # width of each raw roll
    c = 0.01rand(0.99:δ:1.99, J) .* W # used cost of raw rolls (c)
    I, J, w, d, W, c
end;

if false # 📘 The original compact Mip formulation---as a benchmark
    Mip = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); 
    JuMP.@variable(Mip, 0 ≤ b[1:J] ≤ 1, Int);
    JuMP.@objective(Mip, Min, c ⋅ b);
    JuMP.@variable(Mip, 0 ≤ n[1:I, 1:J], Int);
    JuMP.@constraint(Mip, [j = 1:J], w ⋅ view(n, :, j) ≤ b[j]W[j]);
    # Note that the following complicating constraint 🟡 can also be a `==` one
    # We prefer write it as `≥` since the dual space will be halved
    JuMP.@constraint(Mip, [i = 1:I], sum(n[i, j] for j in 1:J) ≥ d[i]); # 🟡 [π ≥ 0]
    JuMP.unset_integer.(b)
    JuMP.unset_integer.(n)
    JuMP.optimize!(Mip) # 📘 [Natural LP Relaxation] Optimal objective 647.4510818
    JuMP.set_integer.(b)
    JuMP.set_integer.(n)
    JuMP.optimize!(Mip)
    # Gurobi's MIP logging
    # 6576  4630  647.68526   88  216  667.02107  647.47885  2.93%   6.3   10s
    # 1498325 1290411  650.72418 1504  203  654.50850  647.48232  1.07%  11.4 6482s
    # 📘 Conclusion: Gurobi spends 6482s to reach ub = 654.50850 | 647.48232 = lb 
end

# Dual Decomposition Assets
function Min_j_fixed_part(j)
    Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.@variable(Min_j, b_j, Bin);
    JuMP.@expression(Min_j, prim_obj_tbMin, c[j]b_j)
    JuMP.@variable(Min_j, 0 ≤ n_j[1:I], Int);
    JuMP.@constraint(Min_j, w ⋅ n_j ≤ W[j]b_j)
    return Min_j
end; function penalty(j, π, n_j::Array{JuMP.VariableRef})
    model = JuMP.owner_model(n_j[1])
    return JuMP.@expression(model, -π ⋅ n_j)
end; function penalty(j, π::Array{JuMP.VariableRef}, n_j)
    model = JuMP.owner_model(π[1])
    return JuMP.@expression(model, -π ⋅ n_j)
end; function draw_a_block!(bit_vec, Θ, Θ_old) # 🟠💡 This is an important heuristic.
    # I'm not aware of any better heuristic method to pick a promising block (that can
    # contribute a violating cut to the master)
    max_i, max_v = 0, -Inf
    for (i, undrawn) in enumerate(bit_vec)
        if undrawn
            v = Θ[i] - Θ_old[i] # 💡
            if max_v < v
                max_i, max_v = i, v
            end
        end
    end
    if max_i == 0 # cannot draw one due to all-fathomed
        return false, max_i
    else
        bit_vec[max_i] = false # indicating that max_i is drawn this time
        return true, max_i # can be drawn, a valid index
    end
end;

# build subproblems (MIPs)
Min_vec = [Min_j_fixed_part(j) for j = 1:J]; # The inner subproblems
# allocate these to collect columns identified for each block
n_cols = [[zeros(Int, I)] for j = 1:J]; pop!.(n_cols);
b_cols = [[false] for j = 1:J];         pop!.(b_cols);

# build master problem (LP)
an_UB = 685 # ⚠️ [hyperparam] a valid ub obtained from Gurobi's trial solve
out = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # The outer dual problem
JuMP.@variable(out, 0 ≤ π[1:I]); JuMP.@variable(out, θ[1:J]);
JuMP.@expression(out, common, d ⋅ π);
JuMP.@expression(out, outer_obj_tbMax, common + sum(θ));
JuMP.@constraint(out, outer_obj_tbMax ≤ an_UB);
JuMP.@objective(out, Max, outer_obj_tbMax);

# The following main loop generates will a best possible dual variable via CTPLNs
JuMP.set_silent.(Min_vec); JuMP.optimize!.(Min_vec);
JuMP.set_silent(out); JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
Θ_old = JuMP.value.(θ); # 💡 need this for my heuristic
cnt_outˈs_cut = 0; # count the total number of cuts added to the outer model
cnt_futile_draw = 0; # a bad block-drawing strategy 🟠💡 might make this index large
t0 = time();
for COT = 5 * (0.5 .^ (0:16)) # Cut-Off Tolerance ⚠️ CANNOT by too small
    while true # work with the present COT
        exists_j_at_COT = false
        JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
        ub = JuMP.objective_bound(out)
        Π, Θ = JuMP.value.(π), JuMP.value.(θ)
        Θˈs_bitvec = trues(J) # 💡 need this for my heuristic
        while true # see if there is any block who can bring us a COT-violating cut
            still_has_undrawn, j = draw_a_block!(Θˈs_bitvec, Θ, Θ_old)
            still_has_undrawn || break
            Min_j = Min_vec[j]; prim_obj_tbMin, n_j = Min_j[:prim_obj_tbMin], Min_j[:n_j]
            JuMP.@objective(Min_j, Min, prim_obj_tbMin + penalty(j, Π, n_j)) # 🟢
            JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
            ObjVal_j = JuMP.objective_value(Min_j)
            if ObjVal_j + COT < Θ[j] # the hallmark of CTPLN's progress
                # Since we implement a single-cut algorithm, we can compress Logging into one line
                @info "#cut = $cnt_outˈs_cut, ub ▶ $ub | cut found at $j with vio = $(Θ[j] - ObjVal_j)"
                Θ_old = Θ # save the RHS, as it is to be updated
                Prim_obj_tbMin, N_j = JuMP.value(prim_obj_tbMin), JuMP.value.(n_j)
                JuMP.@constraint(out, θ[j] ≤ Prim_obj_tbMin + penalty(j, π, N_j)) # 🟢
                push!(n_cols[j], round.(Int, N_j))
                push!(b_cols[j], round(Bool, JuMP.value(Min_j[:b_j])))
                cnt_outˈs_cut += 1; exists_j_at_COT = true
                break # [Early Break] go back to outer model to update Π
            else
                cnt_futile_draw += 1
            end
        end
        if ! exists_j_at_COT
            @info "COT = $COT, elapsed $(time() - t0), futile $cnt_futile_draw, ub $ub ◀◀◀"
            break
        end
    end
end
@info "▶▶▶ total time elapsed is $(time() - t0)"
# retrieve that best dual solution
JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false); Π = JuMP.value.(π)

# Use the columns you had collected to build a restricted version of the original Mip
# The aim is to derive a good primal-side feasible solution
rMip = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # restricted primal side Mip
JuMP.@variable(rMip, rb[1:J], Int); # `r` reads `restricted`
JuMP.@variable(rMip, rn[1:I, 1:J], Int);
let
    JuMP.@variable(rMip, lb[j = 1:J, v = 1:length(b_cols[j])] ≥ 0); # each v indexes a vertex
    JuMP.@variable(rMip, ln[j = 1:J, v = 1:length(n_cols[j])] ≥ 0); # `l` reads `lambda`
    JuMP.@constraint(rMip, [j = 1:J], sum(lb[j, :]) == 1);
    JuMP.@constraint(rMip, [j = 1:J], sum(ln[j, :]) == 1);
    JuMP.@constraint(rMip, [j = 1:J], rb[j] == sum(lb[j, v]b_cols[j][v]                 for v = 1:length(b_cols[j])));
    JuMP.@constraint(rMip, [i = 1:I, j = 1:J], rn[i, j] == sum(ln[j, v]n_cols[j][v][i]  for v = 1:length(n_cols[j])));
end;
JuMP.@constraint(rMip, [j = 1:J], w ⋅ view(rn, :, j) ≤ rb[j]W[j]);
JuMP.@constraint(rMip, [i = 1:I], sum(rn[i, j] for j in 1:J) ≥ d[i]); # 🟡
JuMP.@objective(rMip, Min, c ⋅ rb);
# This model may still take some time to solve to exact optimality, but it's a lot easier than the unrestricted original
# you may also interrupt early and be content with a good feasible (suboptimal) solution
JuMP.optimize!(rMip); # JuMP.assert_is_solved_and_feasible(rMip; allow_local = false)
b_start = JuMP.value.(rb);
n_start = JuMP.value.(rn);

# Finally build an sandwiched version of the original Mip such that
# it can carry on searching the global optimal solution to the original Mip
lbv = let
    lbv = zeros(J)
    for j = 1:J
        Min_j = Min_vec[j]
        b_j = Min_j[:b_j]
        n_j = Min_j[:n_j]
        JuMP.@objective(Min_j, Min, c[j]b_j - Π ⋅ n_j)
        JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
        lbv[j] = JuMP.objective_value(Min_j)
    end
    lbv
end;
swMip = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # sandwiched Mip
JuMP.@variable(swMip, b[1:J], Bin);
JuMP.@objective(swMip, Min, c ⋅ b);
JuMP.@variable(swMip, 0 ≤ n[1:I, 1:J], Int);
JuMP.@constraint(swMip, [j = 1:J], w ⋅ view(n, :, j) ≤ b[j]W[j]);
JuMP.@constraint(swMip, [i = 1:I], sum(n[i, j] for j in 1:J) ≥ d[i]); # 🟡
JuMP.@constraint(swMip, [j = 1:J], c[j]b[j] - Π ⋅ view(n, :, j) ≥ lbv[j]); # 🥪 The lower bread of the sandwich---DWB cuts
JuMP.set_start_value.(b, b_start); JuMP.set_start_value.(n, n_start); # 🥪 The upper bread of the sandwich
JuMP.optimize!(swMip) # 🖥️ U can monitor Gurobi's logging

The previous posts are restricted to only 1-scenario.
In this post I add an axis of scenarios so that the original formulation becomes a Two-stage SIP.
The method to decompose it into (j, s)-blocks and calculating the corresponding Lagrangian bound is

code
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.⋅ as ⋅
import Random
seed = -1214589140429765082 # may need to cherry-pick
Random.seed!(seed)

J = 9; # number of resources
W = rand(9:49, J); # the sole width type of each resource
o = rand(30:0.01:70, J); # occupied cost
u = 7 * rand(0.75:0.01:1.25, J)/30 .* W; # unit price
L = rand(6:17, J); # purchasing num limit
I = 10; # number of customer
w = rand(7:0.7:24, I); # width needed
S = 2; # number of scenes
d = rand(6:35, I, S); # number needed
c = rand(0:0.001:0.1, I, J); # transport cost

if false 
    Xtn = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # extensive formulation
    JuMP.@variable(Xtn, b_free[1:J]);
    JuMP.@variable(Xtn, m_free[1:J]);
    JuMP.@variable(Xtn, 0 ≤ b[1:J, 1:S] ≤ 1, Int); # use this site or not
    JuMP.@variable(Xtn, 0 ≤ m[1:J, 1:S], Int); # order an amount
    JuMP.@variable(Xtn, 0 ≤ n[1:I, 1:J, 1:S], Int); # number dispatched
    JuMP.@constraint(Xtn, [j = 1:J, s = 1:S], L[j]b[j, s] ≥ m[j, s]); # source num is limited
    JuMP.@constraint(Xtn, [j = 1:J, s = 1:S], m[j, s]W[j] ≥ w ⋅ view(n, :, j, s)); # source length is limited
    JuMP.@constraint(Xtn, [j = 1:J, s = 1:S], b_free[j] == b[j, s]); # [μ] 🟡
    JuMP.@constraint(Xtn, [j = 1:J, s = 1:S], m_free[j] == m[j, s]); # [π] 🟡
    JuMP.@constraint(Xtn, [i = 1:I, s = 1:S], sum(n[i, j, s] for j = 1:J) ≥ d[i, s]); # 🟡 [β ≥ 0] number of demand is met
    JuMP.@objective(Xtn, Min, sum(+(o ⋅ view(b, :, s)/S, u ⋅ view(m, :, s)/S, c ⋅ view(n, :, :, s)/S) for s = 1:S));
    JuMP.unset_integer.(b);
    JuMP.unset_integer.(m);
    JuMP.unset_integer.(n);
    # JuMP.set_silent(Xtn);
    JuMP.optimize!(Xtn); JuMP.assert_is_solved_and_feasible(Xtn)
    LP_val = JuMP.objective_value(Xtn);
    JuMP.set_integer.(b);
    JuMP.set_integer.(m);
    JuMP.set_integer.(n);
    JuMP.optimize!(Xtn); JuMP.assert_is_solved_and_feasible(Xtn)
    IP_val = JuMP.objective_value(Xtn);
    @info "LP $LP_val | $IP_val IP"
end

function penalty(j, s, β, π, μ, b, m, n::Array{JuMP.VariableRef}) # disaggregated to each js-block
    model = JuMP.owner_model(n[1])
    return JuMP.@expression(model, μ[j, s]b + π[j, s]m - sum(β[i, s] * n[i] for i = 1:I))
end; function penalty(j, s, β, π::Array{JuMP.VariableRef}, μ, b, m, n) # disaggregated to each js-block
    model = JuMP.owner_model(π[1])
    return JuMP.@expression(model, μ[j, s]b + π[j, s]m - sum(β[i, s] * n[i] for i = 1:I))
end; function Min_js_model(j, s)
    Min_js = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.@variable(Min_js, b, Bin)
    JuMP.@variable(Min_js, 0 ≤ m, Int)
    JuMP.@variable(Min_js, 0 ≤ n[1:I], Int)
    JuMP.@constraint(Min_js, L[j]b ≥ m)
    JuMP.@constraint(Min_js, W[j]m ≥ w ⋅ n);
    JuMP.@expression(Min_js, prim_obj_tbMin, +(o[j]b/S, u[j]m/S, sum(c[i, j]n[i]/S for i = 1:I)))
    return Min_js
end; function draw_a_block!(bitmat, Θ, Θ_old) # 💡 This is an important heuristic.
    max_j, max_s, max_v = 0, 0, -Inf
    J, S = size(bitmat)
    for j = Random.shuffle(1:J), s = Random.shuffle(1:S) # for j = 1:J, s = 1:S
        if bitmat[j, s] # undrawn yet
            v = Θ[j, s] - Θ_old[j, s] # 💡
            if max_v < v
                max_j, max_s, max_v = j, s, v
            end
        end
    end
    max_j == 0 && return false, max_j, max_s
    bitmat[max_j, max_s] = false # indicating it's drawn this time
    return true, max_j, max_s # can be drawn, a valid index
end;

an_UB = 1140.0 # ⚠️ [hyperparam] a valid ub obtained from Gurobi's trial solve
out = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # The outer dual problem
JuMP.@variable(out, 0 ≤ β[1:I, 1:S]);
JuMP.@variable(out, μ[1:J, 1:S]);
JuMP.@variable(out, π[1:J, 1:S]);
JuMP.@constraint(out, [j = 1:J], sum(view(μ, j, :)) == 0); # 🟡
JuMP.@constraint(out, [j = 1:J], sum(view(π, j, :)) == 0); # 🟡
JuMP.@expression(out, common, β ⋅ d);
JuMP.@variable(out, θ[1:J, 1:S]);
JuMP.@expression(out, outer_obj_tbMax, common + sum(θ));
JuMP.@constraint(out, outer_obj_tbMax ≤ an_UB);
JuMP.@objective(out, Max, outer_obj_tbMax);

Min_mat = [Min_js_model(j, s) for j = 1:J, s = 1:S];
JuMP.set_silent.(Min_mat); JuMP.optimize!.(Min_mat);
JuMP.set_silent.(out); JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
Θ = JuMP.value.(θ); Θ_old = copy(Θ);
bitmat = trues(J, S); lbm = similar(θ, Float64);
t0 = time();
t_disposable = 0.;
for COT = 5 * (0.5 .^ (0:17))
    while true # work with the present COT
        JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
        common_value = JuMP.value(common)
        ub = JuMP.objective_bound(out)
        Θ_old .= Θ; Θ = JuMP.value.(θ)
        Β = JuMP.value.(β)
        Μ = JuMP.value.(μ)
        Π = JuMP.value.(π)
        bitmat .= true; exists_j_at_COT = false
        while true # see if there is any block who can bring us a COT-violating cut
            td0 = time()
            still_has_undrawn, j, s = draw_a_block!(bitmat, Θ, Θ_old)
            still_has_undrawn || break # go back to the outermost loop to update COT
            Min_js = Min_mat[j, s]
            prim_obj_tbMin = Min_js[:prim_obj_tbMin]
            n, b, m = Min_js[:n], Min_js[:b], Min_js[:m]
            JuMP.@objective(Min_js, Min, prim_obj_tbMin + penalty(j, s, Β, Π, Μ, b, m, n)) # 🟢
            JuMP.optimize!(Min_js); JuMP.assert_is_solved_and_feasible(Min_js; allow_local = false)
            lbm[j, s] = ObjVal = JuMP.objective_value(Min_js)
            overestimate = Θ[j, s] - ObjVal
            if overestimate > COT
                @info "ub ▶ $ub | cut found at ($j, $s) with vio = $overestimate > $COT = COT"
                Prim_obj_tbMin = JuMP.value(prim_obj_tbMin)
                N = JuMP.value.(n)
                M = JuMP.value.(m)
                B = JuMP.value.(b)
                JuMP.@constraint(out, θ[j, s] ≤ Prim_obj_tbMin + penalty(j, s, β, π, μ, B, M, N)) # 🟢
                exists_j_at_COT = true
                break # [Early Break] go back to outer model to update Β, Μ, Π
            end
            t_disposable += time() - td0
        end
        if ! exists_j_at_COT
            lb = common_value + sum(lbm)
            @info "COT = $COT finished, lb = $lb | $ub = ub ◀◀◀"
            break
        end
    end
end
t_total = time() - t0;
t_thin = t_total - t_disposable;
thin_ratio = round(t_thin / t_total; digits = 2);
println("t_thin = $t_thin < t_total = $t_total (ratio = $thin_ratio)")

Add some comments to my previous post:

For an MIP problem, there are 2 solution approaches. One is solving it directly with a solver (default).
When the default approach performs unsatisfyingly, we may then consider writing decomposition algorithms ourselves. Nevertheless, there is no guarantee that the decomposition algorithm can necessarily outperform a default single solve. (just like other MIP techniques, e.g. whether to adopt certain types of cutting planes)

The idea of block decomposition is to reduce a monolithic MIP into an array of smaller blocks. An important concern is that: we should make sure that each block subproblem (might be an MIP) is easy enough to be solved in short time. e.g., we had best consider

JuMP.set_attribute(Subproblem, "MIPGap", 0);
JuMP.set_attribute(Subproblem, "MIPGapAbs", 0);
JuMP.set_attribute(Subproblem, "TimeLimit", 0.5); # some very short time

If these expectations are unmet, it might suggest that our decomposition manner is inappropriate. If we do not adjust MIPGap to zero, then we might need to consider the accumulation of errors—since we may have large number of blocks.

An Lagrangian dual bound occurs after we decide which set of constraints are complicating.
Therefore there are 4 scalars to be aware of, e.g. the following is an instance

Heuristic     3.71986e+05 # a heuristic suboptimal ObjVal---a byproduct of the `LagBound`
IP            3.719130000000e+05 # The optimal ObjVal of the original MIP
LagBound      3.71673708333333e+05 # The Lagrangian dual bound
LP            3.671042390e+05 # The natural LP relaxation's ObjVal

The easiest of the 4 is the last LP value.
The hardest is certainly IP (our aim).
My “sandwitch” method in my post above is precisely about the Heuristic and LagBound.

To derive the LagBound, apart from those that might happen inside the easy subproblems, we do not need to do any branching. This is why it sounds appealing to derive a Lagrangian bound, which might be tighter than LP. To derive the Heuristic and IP, we might need to do a lot of branchings, but the depth of the Heuristic might be shorter than that of IP.

Typically, we might expect that the rgap between Heuristic and LagBound is small, e.g. < 0.1%, or < 1% (which is very likely). Sometimes we can even expect that with this “sandwitch” strenghening techniques, the subsequent global optimality searching could be easier. (But this is not always). Chances are: the solver might perform better without these additional strengthening techniques. And here is an example.

Code
import JuMP.value as ı
import LinearAlgebra.⋅ as ⋅
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import Random
import SparseArrays.sprand as sprand
function draw_a_block!(bitvec, Θ, Θ_old) # This follows the Max convention
    max_j, max_v = 0, -Inf
    J = length(bitvec)
    for j = Random.shuffle(1:J)
        if bitvec[j] # undrawn yet
            v = Θ[j] - Θ_old[j] # 💡
            if max_v < v
                max_j, max_v = j, v
            end
        end
    end
    max_j == 0 && return false, max_j
    bitvec[max_j] = false # indicating it's drawn this time
    return true, max_j # can be drawn, a valid index
end; function out_model(s, an_UB)
    out = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # The outer dual problem
    JuMP.@variable(out, 0 ≤ β[1:I]); JuMP.@expression(out, common, β ⋅ view(d, :, s))
    JuMP.@variable(out, θ[1:J]); JuMP.@expression(out, out_obj_tbMax, common + sum(θ))
    JuMP.@constraint(out, out_obj_tbMax ≤ an_UB + 1) # ⚠️ fetch the primal feasible value from Gurobi's failed solve
    JuMP.@objective(out, Max, out_obj_tbMax)
    return out
end; function Min_j_model(j, Πb, Πm, Π0) 
    Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.set_silent(Min_j)
    JuMP.set_attribute(Min_j, "TimeLimit", 4) # This model is easy enough to be expected to be solved in 5 s
    JuMP.set_attribute(Min_j, "MIPGap", 0) # This model is easy enough to be expected to be solved in 5 s
    JuMP.set_attribute(Min_j, "MIPGapAbs", 0) # This model is easy enough to be expected to be solved in 5 s
    JuMP.@variable(Min_j, b, Bin)
    JuMP.@variable(Min_j, 0 ≤ m, Int); JuMP.@constraint(Min_j, L[j]b ≥ m)
    JuMP.@variable(Min_j, 0 ≤ n[1:I], Int); JuMP.@constraint(Min_j, W[j]m ≥ w ⋅ n)
    JuMP.@expression(Min_j, prim_obj_tbMin, Πb[j]b + Πm[j]m + Π0/S * view(c, :, j) ⋅ n) # penalty is "-β ⋅ n"
    return Min_j
end;
seed = 5115625092019099801
Random.seed!(seed)
J = 175; # number of resources
W = rand(9:67, J); # the sole width type of each resource
L = rand(2:6, J); # purchasing num limit
S = 2; # number of scenes
s = 2; # currently working at scene s
I = 10; # number of customer
w = rand(2:197, I); # width needed
d = rand(6:83, I, S); # number needed
# c = sprand(I, J, .85, n -> rand(0:9, n)); # transport cost
c = rand(0:39, I, J); # transport cost
Π0 = rand(1:4);
Πb = rand(500:2500, J);
Πm = rand(77:500, J);
an_UB = 1e6 - 1 # ⚠️ must be valid

function run_original(I, J, S, s, L, W, w, d, c, cb, cm, c0)
    model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.set_attribute(model, "MIPGap", 0)
    JuMP.set_attribute(model, "MIPGapAbs", 0)
    JuMP.@variable(model, 0 ≤ b[1:J] ≤ 1, Int)
    JuMP.@variable(model, 0 ≤ m[1:J], Int); JuMP.@constraint(model, [j = 1:J], L[j]b[j] ≥ m[j])
    JuMP.@variable(model, 0 ≤ n[1:I, 1:J], Int)
    JuMP.@objective(model, Min, +(cb ⋅ b, cm ⋅ m, (c0/S)c ⋅ n))
    JuMP.@constraint(model, [j = 1:J], m[j]W[j] ≥ w ⋅ view(n, :, j))
    JuMP.@constraint(model, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s]) # 🟡 [β ≥ 0]
    JuMP.unset_integer.([b; m; vec(n)])
    JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
    @info "LP_Relax_obj = $(JuMP.objective_value(model))"
    JuMP.set_integer.([b; m; vec(n)])
    JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
    @info "IP_obj = $(JuMP.objective_value(model))" occupy_cost = ı(cb ⋅ b) buy_cost = ı(cm ⋅ m) transport_cost = ı((c0/S)c ⋅ n)
end; run_original(I, J, S, s, L, W, w, d, c, Πb, Πm, Π0)

Min_vec = Vector{JuMP.Model}(undef, J)
for j = 1:J
    print("\rj = $j")
    Min_vec[j] = model = Min_j_model(j, Πb, Πm, Π0)
    JuMP.optimize!(model)
end
out = out_model(s, an_UB); JuMP.set_silent(out);
vset_vec = [Set{Vector{Int64}}() for j = 1:J];
Θ, Β, lbv = zeros(J), zeros(I), zeros(J);
out_COT = 1e4 # 💡 This need to be tuned for each specific application
while true # [The Entry]
    Θ_old = Θ; JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
    Θ, Β = ı.(out[:θ]), ı.(out[:β]) # 🔰
    bitvec = trues(J); exists_j_at_COT = false
    while true # see if there is any block who can bring us a COT-violating cut
        still_has_undrawn, j = draw_a_block!(bitvec, Θ, Θ_old)
        if ! still_has_undrawn
            ub_LagDualBnd = JuMP.objective_bound(out)
            lb_LagDualBnd = ı(out[:common]) + sum(lbv)
            @info "COT = $out_COT phase over, lb = $lb_LagDualBnd | $ub_LagDualBnd = ub"
            break # go back to update COT
        end
        Min_j = Min_vec[j]
        JuMP.@objective(Min_j, Min, Min_j[:prim_obj_tbMin] - Β ⋅ Min_j[:n]) # 🟢
        JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false) # ✅ This solve is easy so that exact
        lbv[j] = lbv_j = JuMP.objective_value(Min_j) # write this conveniently
        overestimate = Θ[j] - lbv_j
        if overestimate > out_COT
            B, M, N = ı(Min_j[:b]), ı(Min_j[:m]), ı.(Min_j[:n])
            push!( vset_vec[j], round.(Int, [B; M; N]) ) # primal side
            @info "cut off out's trial with vio = $overestimate from block $j"
            JuMP.@constraint(out, out[:θ][j] ≤ ı(Min_j[:prim_obj_tbMin]) - out[:β] ⋅ N) # 🟢 dual side
            exists_j_at_COT = true
            break # [Early Break]
        end
    end
    exists_j_at_COT && continue
    if out_COT == 1e-11
        break # [The Exit]
    elseif out_COT < 0.01 # 💡 This need to be tuned for each specific application
        out_COT = 1e-11 # The last resort
    else
        out_COT *= 0.2 # preciser
    end
end

vset_vec = collect.(vset_vec);
Rhs = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # convex-combination restricted model
JuMP.set_attribute(Rhs, "MIPGap", 0);
JuMP.set_attribute(Rhs, "MIPGapAbs", 0);
JuMP.@variable(Rhs, b[1:J], Int);
JuMP.@variable(Rhs, m[1:J], Int);
JuMP.@variable(Rhs, n[1:I, 1:J], Int);
JuMP.@variable(Rhs, 0 ≤ λ[j = 1:J, v = 1:length(vset_vec[j])]); JuMP.@constraint(Rhs, [j = 1:J], sum(λ[j, :]) == 1);
JuMP.@constraint(Rhs, [j = 1:J], [b[j]; m[j]; view(n, :, j)] == sum(vset_vec[j][v]λ[j, v] for v = 1:length(vset_vec[j])));
JuMP.@constraint(Rhs, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s]); # 🟡 [β ≥ 0]
JuMP.@objective(Rhs, Min, +(Πb ⋅ b, Πm ⋅ m, (Π0/S)c ⋅ n));
JuMP.optimize!(Rhs); JuMP.assert_is_solved_and_feasible(Rhs; allow_local = false)
JuMP.objective_value(Rhs) # primal side heuristic objective value
B, M, N = ı.(Rhs[:b]), ı.(Rhs[:m]), ı.(Rhs[:n]) # primal side heuristic solution

Rhs = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # strengthened original formulation
JuMP.set_attribute(Rhs, "MIPGap", 0);
JuMP.set_attribute(Rhs, "MIPGapAbs", 0);
JuMP.@variable(Rhs, 0 ≤ b[1:J] ≤ 1, Int)
JuMP.@variable(Rhs, 0 ≤ m[1:J], Int); JuMP.@constraint(Rhs, [j = 1:J], L[j]b[j] ≥ m[j])
JuMP.@variable(Rhs, 0 ≤ n[1:I, 1:J], Int)
JuMP.@objective(Rhs, Min, +(Πb ⋅ b, Πm ⋅ m, (Π0/S)c ⋅ n))
JuMP.@constraint(Rhs, [j = 1:J], m[j]W[j] ≥ w ⋅ view(n, :, j))
JuMP.@constraint(Rhs, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s]) # 🟡 [β ≥ 0]
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### 
JuMP.@constraint(Rhs, [j = 1:J], Πb[j]b[j] + Πm[j]m[j] + (Π0/S)view(c, :, j) ⋅ view(n, :, j) - Β ⋅ view(n, :, j) ≥ lbv[j]) # 🟣 DWB
JuMP.set_start_value.(b, B)
JuMP.set_start_value.(m, M)
JuMP.set_start_value.(n, N)
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### 
JuMP.optimize!(Rhs); JuMP.assert_is_solved_and_feasible(Rhs; allow_local = false)

Finally, one piece of my experience: an MIP with all integer data input will be much easier for Gurobi to solve (e.g. to exact optimality). For example, I have only 16GB RAM, yet Gurobi can solve a 100000-all-integer-variable model to exact optimality steadily.

And additionally: to perform dual decomposition with large number of blocks, one may be eager for parallel computing hardware to do experiments.

Hi @WalterMadelim, I’m going to ask again that you consider setting up a personal blog for these sorts of posts. This forum is intended for Julia- (and JuMP-) related support questions. It isn’t intended to be a venue for longer form posts about algorithms.

We do greatly appreciate your “is this a feature or a bug” posts, and the Gurobi support team knows who “Walter from JuMP” is because of your many interesting bug reports. (This one in particular generated some nice internal discussion Gurobi reports OPTIMAL on a QP having a huge MaxVio)

So take this message as both a “yes, please, keep up the high quality questions” and a “no, thanks, less of these posts which are more suited to a personal blog.”

an MIP with all integer data input will be much easier for Gurobi to solve (e.g. to exact optimality)

I’ll repeat again: there are no performance guarantees. This claim is not true in general. (It can be true in some very specific cases. You may want to look up “total unimodularity”.)

2 Likes

I’m writing algorithmic code in JuMP…

The motivation is that there is not existing algorithmic code so that I have to write them myself.
e.g. Benders’ decomposition or dual decomposition for general MIPs.

I don’t think my posts are off-topic, and I think they will be insightful…

You missed the important part of Oscar’s message (emphasis mine):

Take a look at the conversation above. You have received only two answers from someone who wasn’t Oscar or myself. Then it’s just you conversing with yourself over five consecutive posts. To me, this shows that 1) given your expertise, you don’t really need “support”, at least not on this topic and 2) the community may not be reacting to your longer-form posts the way you expect it to.
We are not denying that these algorithmic expositions are insightful. However, we are telling you that they are not a good fit for the present forum (unlike the posts where you uncover potential bugs, which are of very high value).

Oscar is the maintainer of JuMP.jl and I’m a moderator of this forum, as well as a researcher in combinatorial optimization. We do know what we are talking about, both in terms of the science and in terms of the relevance to the Julia community.
In the past, we have made the same remarks to you repeatedly: here, here, here, and in the rest of those threads. Several of your (early) topics did not get a single reaction (here, here, here, here), possibly because they read more like notes taken for your own future use than material designed to help others on a forum.
Your use of Discourse would be perfectly okay if forum users and forum moderators had infinite time and tolerance to noise, but such is not the case. Expository posts like these (“this is more a comment than a question”) draw attention and energy away from people who actually need help. Oscar in particular is an exceptional asset to our community, who spends a lot of time on Discourse already. Let’s make it easier for him to help people who need it.

6 Likes

Actually I’m not expecting anything. If anyone is interested, then maybe they will take a look. Otherwise they will be indifferent.

I think, if you just be indifferent, like other ordinary people, then we all have energy conserved, and devote them to those who need help. Would this be a better solution?

Any this is the literal meaning of discourse:

a long and serious treatment or discussion of a subject in speech or writing

In this topic I’m writing two decomposition algorithms: dual decomposition and Benders’ decomposition. I think they are highly relevant to most users. And it is because I haven’t seen anyone write them (e.g. with block structures, or for general integer cases). Therefore I’m writing them myself.

I think: if there is no one come out to say “hay, here is one error in your code, you should do this…”. Then it might suggest that my implementation is pretty state-of-the-art.

Or if anyone is too busy and be indifferent, then none of us lose anything, right?

All the posts under this topic are adopting the same instance as the one in my top post—it’s JuMP’s doc. Therefore I don’t think there is any irrelevance

The above posts are concerned with Dual Decomposition. And I’ve just finished a Bender Decomposition implementation. And this should be the end. I think this will be a good reference.

There are some novel ideas that is thought out of my head. I think they will be helpful to the numeric stability. If anyone else can think of better implementations, please inform me.

Scenario Decomposition---Benders method
import JuMP.value as ı
import LinearAlgebra.⋅ as ⋅
import LinearAlgebra.norm as norm
import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import Random

# seed = rand(Int)
seed = 5984114361341617715
Random.seed!(seed)
J = 3; # number of resources
W = rand(13:400, J); # the sole width type of each resource
L = rand(2:15, J); # purchasing num limit
S = 6; # number of scenes
I = 10; # number of customer
w = rand(2:77, I); # width needed
d = rand(2:13, I, S); # number needed
c = rand(0:0.001:39, I, J); # transport cost
cb = rand(500:0.001:2500, J);
cm = rand(77:0.001:500, J);

function run_original(I, J, S, L, W, w, d, c, cb, cm)
    model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.set_silent(model)
    JuMP.set_attribute(model, "MIPGap", 0)
    JuMP.set_attribute(model, "MIPGapAbs", 0)
    JuMP.@variable(model, 0 ≤ b[1:J] ≤ 1, Int)
    JuMP.@variable(model, 0 ≤ m[1:J], Int); JuMP.@constraint(model, [j = 1:J], L[j]b[j] ≥ m[j])
    JuMP.@variable(model, 0 ≤ n[1:I, 1:J, 1:S], Int)
    JuMP.@objective(model, Min, +(cb ⋅ b, cm ⋅ m, sum(c/S ⋅ view(n, :, :, s) for s = 1:S)))
    JuMP.@constraint(model, [j = 1:J, s = 1:S], m[j]W[j] ≥ w ⋅ view(n, :, j, s))
    JuMP.@constraint(model, [i = 1:I, s = 1:S], sum(view(n, i, :, s)) ≥ d[i, s])
    JuMP.unset_integer.([b; m; vec(n)])
    JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
    @info "LP_Relax_obj = $(JuMP.objective_value(model))"
    JuMP.set_integer.([b; m; vec(n)])
    JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
    @info "IP_obj = $(JuMP.objective_value(model))" occupy_cost = ı(cb ⋅ b) buy_cost = ı(cm ⋅ m) transport_cost = ı(sum((1/S)c ⋅ view(n, :, :, s) for s = 1:S))
end; run_original(I, J, S, L, W, w, d, c, cb, cm)

function draw_a_block!(bitvec, Θ, Θ_old) # This follows the Max convention
    max_j, max_v = 0, -Inf
    J = length(bitvec)
    for j = Random.shuffle(1:J)
        if bitvec[j] # undrawn yet
            v = Θ[j] - Θ_old[j] # 💡
            if max_v < v
                max_j, max_v = j, v
            end
        end
    end
    max_j == 0 && return false, max_j
    bitvec[max_j] = false # indicating it's drawn this time
    return true, max_j # can be drawn, a valid index
end; function value_fun(s, b, m)
    model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); JuMP.set_silent(model)
    JuMP.@variable(model, 0 ≤ n[1:I, 1:J], Int)
    JuMP.@constraint(model, [j = 1:J], m[j]W[j] ≥ w ⋅ view(n, :, j))
    JuMP.@constraint(model, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s]) 
    JuMP.@objective(model, Min, c/S ⋅ n)
    JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
    JuMP.objective_value(model)
end; function Rhs_model(s)
    Rhs = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.set_attribute(Rhs, "TimeLimit", 3)
    JuMP.set_attribute(Rhs, "MIPGap", 0)
    JuMP.set_attribute(Rhs, "MIPGapAbs", 0)
    JuMP.@variable(Rhs, 0 ≤ b[1:J] ≤ 1, Int)
    JuMP.@variable(Rhs, 0 ≤ m[1:J], Int); JuMP.@constraint(Rhs, [j = 1:J], L[j]b[j] ≥ m[j])
    JuMP.@objective(Rhs, Min, 0) # TODO +(Πb ⋅ b, Πm ⋅ m, (Π0/S)c ⋅ n) # TODO third term also being ((c ⋅ n)/S)Π0
    JuMP.@variable(Rhs, 0 ≤ n[1:I, 1:J], Int)
    JuMP.@constraint(Rhs, [j = 1:J], m[j]W[j] ≥ w ⋅ view(n, :, j))
    JuMP.@constraint(Rhs, [i = 1:I], sum(view(n, i, :)) ≥ d[i, s]) # 🟡 [β ≥ 0]
    return Rhs
end; function Max_s_model(s)
    Sur = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.@variable(Sur, πb[1:J])
    JuMP.@variable(Sur, πm[1:J])
    JuMP.@variable(Sur, 0 ≤ π0) # ⚠️ You must also double check the solution being positive
    JuMP.@constraint(Sur, [true; πb; πm; π0] ∈ JuMP.MOI.NormOneCone(2J+2))
    JuMP.@variable(Sur, θ) # ✅ also recruit cuts added to `Max_s_no_π0`
    JuMP.@objective(Sur, Max, θ) # TBD
    return Sur
end; function Mst_model()
    Mst = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # Benders' Master LP
    JuMP.@variable(Mst, 0 ≤ b[1:J] ≤ 1);
    JuMP.@variable(Mst, 0 ≤ m[1:J]); JuMP.@constraint(Mst, [j = 1:J], L[j]b[j] ≥ m[j]);
    JuMP.@variable(Mst, θ[1:S]); # initially θ[s] < c/S ⋅ view(n, :, :, s)
    JuMP.@expression(Mst, common, cb ⋅ b + cm ⋅ m)
    JuMP.@expression(Mst, Mst_obj_tbMin, common + sum(θ))
    JuMP.@constraint(Mst, Mst_obj_tbMin ≥ 0)
    JuMP.@objective(Mst, Min, Mst_obj_tbMin)
    return Mst
end;

Max_vec = [Max_s_model(s) for s = 1:S]; JuMP.set_silent.(Max_vec);
Rhs_vec = [Rhs_model(s) for s = 1:S]; JuMP.set_silent.(Rhs_vec);
let divisor = norm([cb; cm; 1.0], 1)
    for s = 1:S # add initial bounding cuts to `Sur` in `Max_vec`
        Sur, Rhs = Max_vec[s], Rhs_vec[s] # 🥈
        JuMP.@objective(Rhs, Min, +(cb/divisor ⋅ Rhs[:b], cm/divisor ⋅ Rhs[:m], c/divisor/S ⋅ Rhs[:n]))
        JuMP.optimize!(Rhs); JuMP.assert_is_solved_and_feasible(Rhs; allow_local = false)
        p_πb, p_πm, p_π0 = ı.(Rhs[:b]), ı.(Rhs[:m]), c/S ⋅ ı.(Rhs[:n]) # primal feasible value => inner valid
        JuMP.@constraint(Sur, Sur[:θ] ≤ p_πb ⋅ Sur[:πb] + p_πm ⋅ Sur[:πm] + Sur[:π0]p_π0) # 🟣 θ-Cuts
        JuMP.optimize!(Sur); JuMP.assert_is_solved_and_feasible(Sur; allow_local = false) # ⚫ BOUNDEDNESS
    end
end
Mst = Mst_model(); JuMP.set_silent(Mst); # 🥇
MstˈΘ = zeros(S);
COT = [100., 10.]; # these 2 cut-related hyperparameters are inter-s
COT_reduce_index = 1
cgs = falses(2); # cut generation status
B, M = zeros(J), zeros(J)
while true # [The Entry]
    MstˈΘ_old = MstˈΘ; JuMP.optimize!(Mst); JuMP.assert_is_solved_and_feasible(Mst; allow_local = false)
    B, M = ı.(Mst[:b]), ı.(Mst[:m]) # 🔰
    MstˈΘ = ı.(Mst[:θ]) # 🔰
    cgs .= false; Bensˈbitvec = trues(S)
    while true # inter-s
        still_has_undrawn, s = draw_a_block!(Bensˈbitvec, MstˈΘ_old, MstˈΘ) # here Min => swap
        still_has_undrawn || break
        Sur, Rhs = Max_vec[s], Rhs_vec[s] # 🥈
        JuMP.@objective(Sur, Max, Sur[:θ] - B ⋅ Sur[:πb] - M ⋅ Sur[:πm] - MstˈΘ[s]Sur[:π0]) # 🔰
        JuMP.optimize!(Sur); JuMP.assert_is_solved_and_feasible(Sur; allow_local = false)
        SurˈΘ, Πb, Πm, Π0 = ı.(Sur[:θ]), ı.(Sur[:πb]), ı.(Sur[:πm]), ı.(Sur[:π0]) # ✏️
        Π0 < 0 && error("Gurobi gives a Π0 < 0")
        JuMP.@objective(Rhs, Min, +(Πb ⋅ Rhs[:b], Πm ⋅ Rhs[:m], (Π0)c/S ⋅ Rhs[:n]))
        JuMP.optimize!(Rhs); JuMP.assert_is_solved_and_feasible(Rhs; allow_local = false)
        rhs_lb = JuMP.objective_bound(Rhs) # dual bound => outer valid
        p_πb, p_πm, p_π0 = ı.(Rhs[:b]), ı.(Rhs[:m]), c/S ⋅ ı.(Rhs[:n]) # primal feasible value => inner valid
        SurˈΘˈvdiff_at_Π = SurˈΘ - (p_πb ⋅ Πb + p_πm ⋅ Πm + (p_π0)Π0) # ✏️
        if SurˈΘˈvdiff_at_Π > COT[2]
            JuMP.@constraint(Sur, Sur[:θ] ≤ p_πb ⋅ Sur[:πb] + p_πm ⋅ Sur[:πm] + Sur[:π0]p_π0) # 🟣 θ-Cuts
            @info "cut off Sur's trial with vio $SurˈΘˈvdiff_at_Π at scene $s"
            cgs[2] = true
        end
        Sur_obj_lb = rhs_lb - B ⋅ Πb - M ⋅ Πm - MstˈΘ[s]Π0 # 🔰✏️
        if Sur_obj_lb > COT[1]
            JuMP.@constraint(Mst, Mst[:b] ⋅ Πb + Mst[:m] ⋅ Πm + Mst[:θ][s]Π0 ≥ rhs_lb) # ✅ Lag cut by block s
            @info "cut off Mst's trial with vio $Sur_obj_lb at scene $s"
            cgs[1] = true
            break # [early break]
        end
    end
    any(cgs) && continue
    if COT[COT_reduce_index] == 1e-11
        COT[3 - COT_reduce_index] == 1e-11 && break # Exit
    elseif COT[COT_reduce_index] < 9e-4
        COT[COT_reduce_index] = 1e-11
    else
        COT[COT_reduce_index] *= 0.2
    end
    COT_reduce_index = 3 - COT_reduce_index
end
@info "lb = $(JuMP.objective_bound(Mst)) | $(ı(Mst[:common]) + sum(value_fun(s, B, M) for s = 1:S)) = ub"
# TODO because I can't find a test case so that the resulting trial points are fractional
# Therefore omit the rest code
# The rest code are: define a Master's Integer Solution Generator
# And use it to generate integer solutions
# once cut off that, go back to the Master's (≡ LP) loop
# And repeat

@WalterMadelim , It would be great to have your 2 decomposition algorithms.
I think what the moderators try steer for is raising the SNR level.

I think some editing work would be great.
After you achieved your goal, Could you summarize the work into a single post and delete the intermediate posts?
You may also consider a Forem Post and link it here.
I’d be happy to read.

1 Like

My work is limited to single-thread environment. But the decomposition algorithm is supposed to be mainly helpful in a multi-thread parallel computing environment. (Otherwise it may be outperformed by a default single solve by a black-box solver). I’m going to study how to implement them in parallel then.

I’ll consider. Thanks for your advice!

Thanks for your interest on this topic.

2 Likes