Combinatorial Benders cuts

Specific Domains Optimization (Mathematical)
1

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:

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

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))
1 Like
2

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

@variable(master, γ >= -1_000)  # Or something
1 Like
3

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?

4

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:

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 for details

5

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

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!

6

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.

7

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.

8

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.