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