# Combinatorial Benders cuts

**URL:** <https://discourse.julialang.org/t/combinatorial-benders-cuts/103012>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [August 20, 2023, 7:04pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012 "2023-08-20T19:04:05Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [August 20, 2023, 7:04pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/1 "2023-08-20T19:04:05Z")

</div>

I have a big-M constraint in a stochastic programming problem that I want to replace with combinatorial Benders cuts. I have tried implementing combinatorial Bender decomposition, but I get unbounded problem error when it first tries to solve the master problem. I think I am nearly there, but unclear on how to rectify it. I have given the codes below for both benders and big-M method. May be someone will be able to help out?  
Here is my code for the big-M formulation:

```julia
using JuMP
using Cbc 
# Define the problem data
n = 10 # Number of optimization variables
S = 100 # Number of scenarios
b = rand(80:140, n) # Vector of constant values of X
α = 0.05 # Confidence level
M = 1000 # Big constant value
ω = rand(40:140,S,n) # Matrix of uncertain values of X. Each column represents a scenario.
P = ones(S) ./ S # Vector of probabilities for each scenario
model = Model(Cbc.Optimizer)
M=1000
# Define variables
@variable(model, 0 <= x[1:n] <= 1) # Optimization variables
@variable(model, γ) # Objective variable (free variable)
@variable(model, y[1:S], Bin) # Binary variables for each scenario

# Objective function
@objective(model, Min, γ)

# Constraints
for s in 1:S
    @constraint(model, M*y[s] >= dot((b - ω[s,:]), x) - γ)
end

@constraint(model, sum(y[s] * P[s] for s in 1:S) <= α)
@constraint(model,sum(x[i] for i =1:n) == 1)

# Solve the model
optimize!(model)

# Display the results
println("Objective Value: ", objective_value(model))
println("Optimal x: ", value.(x))

```

**Code for Benders**

```julia
using JuMP, Cbc

# Define the problem data
n = 10
S = 100
b = rand(80:140, n)
α = 0.05
ω = rand(40:140, S, n)
P = ones(S) ./ S

# Master problem
master = Model(Cbc.Optimizer)
@variable(master, γ)
@variable(master, y[1:S], Bin)
@objective(master, Min, γ)
@constraint(master, sum(y[s]*P[s] for s=1:S) <= α)

# Subproblem
subproblem = Model(Cbc.Optimizer)
@variable(subproblem, 0<=x[1:n] <= 1)
@constraint(subproblem, sum(x[i] for i=1:n) == 1)

max_iterations = 100
for iter = 1:max_iterations
    optimize!(master)

    worst_γ = -Inf
    T = Vector{Int}()
    for s in 1:S
        if value(y[s]) == 0
            @objective(subproblem, Min, dot(b - ω[s, :], x))
            optimize!(subproblem)
            current_γ = objective_value(subproblem)
            if current_γ > value(γ)
                worst_γ = current_γ
                push!(T, s)
            end
        end
    end

    if isempty(T) || worst_γ <= value(γ) # If there's no violating scenario
        break
    else
        @constraint(master, γ >= worst_γ * (sum([1 - y[s] for s in T]) - length(T) + 1))
    end
end

println("Objective Value: ", objective_value(master))
println("Optimal y: ", value.(y))

```

---

<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:** [August 20, 2023, 9:00pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/2 "2023-08-20T21:00:18Z")

</div>

In the first iteration, your problem is unbounded because `γ` can be arbitrarily negative. The common solution it to add a valid lower bound:

```julia
@variable(master, γ >= -1_000) # Or something

```

---

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [August 20, 2023, 9:42pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/3 "2023-08-20T21:42:36Z")

</div>

Thank you. I have tried your suggestion and added a lower bound for `γ` :  
`@variable(master, γ >= 0.0)`  
But I get a very different results for the objective and variables `x` than the big-M formulation. In the Benders version after all iterations, the objective value returned is the lower bound of`γ` and values of `x` are very different.  
Could there be any additional issues in my code?

---

<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:** [August 20, 2023, 10:04pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/4 "2023-08-20T22:04:55Z")

</div>

If they don’t obtain equivalent solutions, then your code for Benders is wrong. Are you sure the cut is the correct formulation?

However, why do you need a decomposition for this problem? The problem solves almost instantly. If your real problem is larger, you could try Gurobi, or an open-source solver like HiGHS.

Gurobi supports indicator constraints, so you could write:

```julia
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, 0 <= x[1:n] <= 1)
@variable(model, γ)          
@variable(model, y[1:S], Bin)
@objective(model, Min, γ)
@constraint(model, [s = 1:S], !y[s] --> {γ >= (b - ω[s,:])' * x})
@constraint(model, sum(P[s] * y[s] for s in 1:S) <= α)
@constraint(model, sum(x) == 1)
optimize!(model)
println("Objective Value: ", objective_value(model))
println("Optimal x: ", value.(x))

```

See [Constraints · JuMP](https://jump.dev/JuMP.jl/stable/manual/constraints/#Indicator-constraints) for details

---

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [August 20, 2023, 11:37pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/5 "2023-08-20T23:37:22Z")

</div>

I doubt that Gurobi will be able to solve a large instance of this problem. It may solve a small instance but if you were to set n = 10000 and S = 10000 it will struggle badly. I suspect performance will be worst with indicator constraints than big-M for large problems, but fast for small problems. This is the motivation behind using combinatorial benders for this problem.

I am trying to use the method in this paper [https://ideas.repec.org/a/inm/oropre/v54y2006i4p756-766.html](https://ideas.repec.org/a/inm/oropre/v54y2006i4p756-766.html)

There could be possibly something wrong in my bender cuts. This is what I’m trying to do: if the value of the objective in the subproblem is greater than the objective of the candidate master solution it adds a cut and T keeps a track of the scenario in which this happens. This effectively should result in the same outcome as the big-M constraint. But obviously something is incorrect!

---

<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:** [August 21, 2023, 2:53am UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/6 "2023-08-21T02:53:40Z")

</div>

I think you need to take another look at your Benders algorithm. It seems like you’re trying to mix combinatorial Benders with regular Benders.

---

<div class="post-metadata">

**Author:** ![papar](https://avatars.discourse-cdn.com/v4/letter/p/c5a1d2/32.png) [@papar](https://discourse.julialang.org/u/papar)\
**Post date:** [August 28, 2023, 7:17pm UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/7 "2023-08-28T19:17:11Z")

</div>

I am trying to apply the logic introduced by Matteo Fischett for combinatorial bender cuts to replace big-M constraints. The logic of bender cut here is quite basic:  
If ( \gamma\_T ) from the SP is larger than the MP’s ( \gamma ), it means our candidate solution has underestimated the necessary ( \gamma ) for the scenarios in ( T ). In this case, we add a Benders cut to the MP:  
\end{itemize}

[\gamma \geq \gamma\_{T} \left( \sum\_{s \in T} (1 - Y\_{s}) - |T| + 1 \right)]

This inequality essentially ensures that, for the same or similar future ( Y ) solutions, ( \gamma ) will need to be at least as large as ( \gamma\_T ) to accommodate the worst-case scenarios in ( T ).

Apologies, I don’t know how to correctly display latex on this forum.

---

<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:** [August 30, 2023, 2:56am UTC](https://discourse.julialang.org/t/combinatorial-benders-cuts/103012/8 "2023-08-30T02:56:17Z")

</div>

You can use LaTeX by putting it between dollar signs, like `$ \gamma $`.

For example: \gamma \geq \gamma\_{T} \left( \sum\_{s \in T} (1 - Y\_{s}) - |T| + 1 \right)

I don’t think this is the right cut for your problem. It will only work if \gamma\_{T} \ge 0.

Otherwise you might add a constraint like \gamma \geq -100 \left( \sum\_{s \in T} (1 - Y\_{s}) - |T| + 1 \right), and if two of the Y\_s are `1`, then the bit in brackets will be `-1`, and the constraint enforces \gamma \> 100.