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

**URL:** <https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, algorithm, mip\
**Created:** [June 6, 2025, 3:13am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686 "2025-06-06T03:13:57Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 6, 2025, 3:13am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/1 "2025-06-06T03:13:57Z")

</div>

> **[Column generation · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/algorithms/cutting_stock_column_generation/#Column-generation)**
>
> Documentation for JuMP.

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.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 6, 2025, 4:35am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/2 "2025-06-06T04:35:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 6, 2025, 7:45am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/3 "2025-06-06T07:45:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 6, 2025, 1:01pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/4 "2025-06-06T13:01:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![guimarqu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/guimarqu/32/211057_2.png) [@guimarqu](https://discourse.julialang.org/u/guimarqu)\
**Post date:** [June 6, 2025, 2:21pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/5 "2025-06-06T14:21:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 6, 2025, 11:57pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/6 "2025-06-06T23:57:03Z")

</div>

> [@guimarqu](#):
>
> At each iteration of column generation, you look at the column that has best reduced cost.

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.

> [@guimarqu](#):
>
> Assume y\_1 = 1 because a master solution containing 0 pattern is feasible.

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.

> [@guimarqu](#):
>
> and this is equivalent to

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

```julia
  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,

```julia
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](https://jump.dev/JuMP.jl/stable/tutorials/algorithms/cutting_stock_column_generation/#Choosing-the-initial-set-of-patterns)" with an Identity matrix, thus made the exposition shorter and smoother.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 7, 2025, 5:13am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/7 "2025-06-07T05:13:57Z")

</div>

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**
>
> ```julia
> 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**
>
> ```julia
> 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**
>
> ```julia
> 
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![guimarqu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/guimarqu/32/211057_2.png) [@guimarqu](https://discourse.julialang.org/u/guimarqu)\
**Post date:** [June 7, 2025, 11:55am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/8 "2025-06-07T11:55:00Z")

</div>

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](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 🙂

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 7, 2025, 2:15pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/9 "2025-06-07T14:15:06Z")

</div>

> [@guimarqu](#):
>
> how to dynamically add variables to a formulation

Yes, I’ve noticed this point

```julia
# 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.)

> [@guimarqu](#):
>
> the work of Vanderbeck & Wolsey

Thanks for sharing! I’ll have a read.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 8, 2025, 2:41am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/10 "2025-06-08T02:41:11Z")

</div>

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**
>
> ```julia
> 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

```julia
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. 😄

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [June 27, 2025, 2:55pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/11 "2025-06-27T14:55:47Z")

</div>

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](https://discourse.julialang.org/t/find-the-max-element-of-a-numeric-vector-iteratively-with-mask-and-early-break/130233))

My implementation is

> **code**
>
> ```julia
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 1, 2025, 2:03pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/12 "2025-07-01T14:03:45Z")

</div>

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**
>
> ```julia
> 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)")
> 
> ```

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 7, 2025, 1:49am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/13 "2025-07-07T01:49:15Z")

</div>

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

```julia
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

```julia
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**
>
> ```julia
> 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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 7, 2025, 3:49am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/14 "2025-07-07T03:49:03Z")

</div>

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](https://discourse.julialang.org/t/gurobi-reports-optimal-on-a-qp-having-a-huge-maxvio/129645))

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”.)

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 7, 2025, 4:07am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/15 "2025-07-07T04:07:07Z")

</div>

> [@odow](#):
>
> for Julia- (and JuMP-) related

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…

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [July 7, 2025, 5:05am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/16 "2025-07-07T05:05:12Z")

</div>

> [@WalterMadelim](#):
>
> I’m writing algorithmic code in JuMP…

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

> [@odow](#):
>
> 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.

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](https://discourse.julialang.org/t/simplicial-decomposition-equals-dual-cutting-plane-a-toy-example/129128/2), [here](https://discourse.julialang.org/t/a-code-example-on-primal-dual-algorithm-for-lp/128148/2), [here](https://discourse.julialang.org/t/approaching-a-linear-program-correctly/127354/3), and in the rest of those threads. Several of your (early) topics did not get a single reaction ([here](https://discourse.julialang.org/t/reformulation-perspectification-technique-an-example-with-a-nonconvex-feasible-region/126915), [here](https://discourse.julialang.org/t/reformulation-perspectification-technique-ideas-embodied-in-a-contrived-example/126804), [here](https://discourse.julialang.org/t/reformulation-perspectification-technique-a-code-example-using-jump/126729), [here](https://discourse.julialang.org/t/a-naive-explanation-of-infeasibility-and-unboundedness/127343)), 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.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 7, 2025, 11:10am UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/17 "2025-07-07T11:10:02Z")

</div>

> [@gdalle](#):
>
> the community may not be reacting to your longer-form posts the way you expect it to.

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

> [@gdalle](#):
>
> In the past, we have made the same remarks to you repeatedly

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`:

```julia
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**
>
> ```julia
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [July 7, 2025, 2:50pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/18 "2025-07-07T14:50:08Z")

</div>

@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](https://forem.julialang.org) and link it here.  
I’d be happy to read.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 7, 2025, 10:38pm UTC](https://discourse.julialang.org/t/a-primal-dual-global-optimization-approach-on-the-cutting-stock-problem/129686/19 "2025-07-07T22:38:49Z")

</div>

> [@RoyiAvital](#):
>
> After you achieved your goal

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.

> [@RoyiAvital](#):
>
> Could you summarize the work into a single post and delete the intermediate posts

I’ll consider. Thanks for your advice!

> [@RoyiAvital](#):
>
> I’d be happy to read.

Thanks for your interest on this topic.