# Recitation 9: Benders decomposition for stochastic facility location


```julia
using DataFrames, CSV
using JuMP, Gurobi
using LinearAlgebra, Random, Printf
using Plots

const GRB_ENV = Gurobi.Env()
```

## 1. Problem setup

We want to solve the following stochastic facility location problem:

$$\begin{align}
\min\quad & \sum_{i=1}^Nc_ix_i+\sum_{s=1}^Sp_s\left(\sum_{i=1}^N\sum_{j=1}^Mt_{ij}y_{ij}^s + \sum_{j=1}^Mq_jz_{j}^s\right)\\
\text{s.t.}\quad & \sum_{i=1}^Ny_{ij}^s+z_j^s\ge d_j^s & \forall j\in [M], s\in [S]\\
& \sum_{j=1}^My_{ij}^s \le C_ix_i&\forall i\in[N], s\in[S]\\
& \mathbf{y,z}\ge 0, \mathbf{x}\in\{0,1\}^N
\end{align}$$

## 2. Data Setup

We need to define quite a bit of data, so we're going to need to think a little more carefully than usual about how to store the data and parameters. First, we can use a NamedTuple to keep track of parameters.


```julia
function create_parameters(N, S)
    parameters = (
        seed = 1234,
        num_facilities = N,
        num_customers = N,
        num_scenarios = S,
        square_size = 100, # km
        unit_cost_per_km = 0.1,
        # demand generation
        min_demand = 1000 / N, max_demand = 2000 / N,
        # cost of unmet demand
        min_unmet_cost = 10.0, max_unmet_cost = 20.0,
        # cost of facility
        min_facility_cost = 20.0, max_facility_cost = 30.0,
        # capacity at each facility
        min_capacity = 10.0, max_capacity = 30.0,
    )
    return parameters
end
```

Let's create a parameters object to play with!


```julia
params = create_parameters(100, 10)
```

Now we can define a data type representing the facility location data. It's nice to have a single container with all the data we need. We could also have used a data frame, but the advantage of this approach is it can handle various shapes (e.g. a distance matrix and a cost vector).


```julia
"All dimensions assumed to be self-consistent"
struct FacilityLocationData
    "List of facility coordinates as n_facilities x 2 matrix"
    facilities::Matrix
    "List of customer coordinates"
    customers::Matrix
    
    "Facility capacities"
    capacity::Vector
    "Cost of unmet demand for each customer"
    cost_unmet_demand::Vector
    "Cost of building each facility"
    facility_cost::Vector
    "Transportation_cost (n_facilites x n_customers)"
    transportation_cost::Matrix
    
    "Demand scenarios"
    demand::Matrix
    "Probability of each scenario"
    prob::Vector
end

"Build from parameters"
function FacilityLocationData(p::NamedTuple)
    Random.seed!(p.seed)
    facilities = rand(p.num_facilities, 2) * p.square_size
    customers = rand(p.num_customers, 2) * p.square_size
    t_cost = [norm(facilities[i, :] .- customers[j, :])
              for i = 1:p.num_facilities, j=1:p.num_customers] .* p.unit_cost_per_km
    f_cost = rand(p.num_facilities) .* (p.max_facility_cost - p.min_facility_cost) .+ p.min_facility_cost
    capacity = rand(p.num_facilities) .* (p.max_capacity - p.min_capacity) .+ p.min_capacity
    unmet_cost = rand(p.num_customers) .* (p.max_unmet_cost - p.min_unmet_cost) .+ p.min_unmet_cost
    demand = rand(p.num_customers, p.num_scenarios) .* (p.max_demand - p.min_demand) .+ p.min_demand
    prob = 1/p.num_scenarios * ones(p.num_scenarios)
    return FacilityLocationData(facilities, customers, capacity, unmet_cost, f_cost, t_cost, demand, prob)
end
```

OK, let's now create a facility location data object using our parameters, and plot some characteristics.


```julia
data = FacilityLocationData(params)
```

This is a really ugly and useless printout! We can actually make this a bit nicer by defining a better print method for the `FacilityLocationData` type. To do this, we override the default print method, which is called `Base.show`.


```julia
function Base.show(io::IO, data::FacilityLocationData)
    @printf(io, "Facility location data with %d facilities, %d customers, and %d demand scenarios",
            size(data.facilities, 1), size(data.customers, 1), size(data.demand, 2))
end
data
```

And now let's plot the facilities and customers!


```julia
scatter(data.facilities[:, 1], data.facilities[:, 2], label="Facilities")
scatter!(data.customers[:, 1], data.customers[:, 2], label="Customers")
```

## 2. Direct optimization


```julia
TIME_LIMIT = 180;
OPTIMALITY_GAP = 0.01;
```


```julia
"Solve the optimization problem directly"
function direct_solve(data::FacilityLocationData; verbose::Bool = true)
    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attributes(model, "TimeLimit" => TIME_LIMIT, "MIPGap" => OPTIMALITY_GAP,
                             "OutputFlag" => ifelse(verbose, 1, 0))
    N = size(data.facilities, 1); M = size(data.customers, 1); S = size(data.demand, 2)
    @variable(model, x[1:N], Bin)
    @variable(model, y[1:N, 1:M, 1:S] >= 0)
    @variable(model, z[1:M, 1:S] >= 0)
    @objective(model, Min,
               sum(data.facility_cost[i] * x[i] for i=1:N) +
               sum(data.prob[s] * data.transportation_cost[i, j] * y[i, j, s] for i=1:N, j in 1:M, s in 1:S) +
               sum(data.prob[s] * data.cost_unmet_demand[j] * z[j, s] for j in 1:M, s in 1:S))
    @constraint(model, [j in 1:M, s in 1:S],
                sum(y[i, j, s] for i = 1:N) + z[j, s] >= data.demand[j,s])
    @constraint(model, [i in 1:N, s in 1:S],
                sum(y[i, j, s] for j=1:M) <= data.capacity[i] * x[i])
    solvetime = @elapsed optimize!(model)
    opt = objective_value(model)
    bound = objective_bound(model)
    return opt, bound, solvetime
end
```


```julia
direct_solve(data, verbose=true)
```

## 3. Benders decomposition

### 3.1. Multi-cut approach

We want to solve the problem using Benders decomposition. Consider the following multi-cut approach.

Main problem:
$$\begin{align}
\min\quad &\sum_{i=1}^N c_i x_i + \sum_{s=1}^Sp_s\theta_s\\
\text{s.t.}\quad & \mathbf{x}\in\{0,1\}^N, \mathbf{\theta}\ge 0
\end{align}$$

Subproblem $s$:
$$\begin{align}
\min\quad & \sum_{i=1}^N\sum_{j=1}^Mt_{ij}y_{ij}^s + \sum_{j=1}^Mq_jz_{j}^s\\
\text{s.t.}\quad & \sum_{i=1}^Ny_{ij}^s+z_j^s\ge d_j^s & \forall j\in [M]\\
& \sum_{j=1}^My_{ij}^s \le C_ix_i&\forall i\in[N]\\
& \mathbf{y,z}\ge 0
\end{align}$$

Dual subproblem $s$:
$$\begin{align}
\max\quad &  \sum_{j=1}^M \mu_j d_j^s - \sum_{i=1}^N\lambda_iC_ix_i\\
\text{s.t.}\quad & \mu_j -\lambda_i \le t_{ij} & \forall i\in[N], j\in[M]\\
& \mu_j \le q_j &\forall j\in [M]\\
& \mathbf{\mu, \lambda}\ge 0
\end{align}$$

- If the dual subproblem is unbounded, we obtain an extreme ray $(\mu^*, \lambda^*)$ and add a feasibility cut:
$$\sum_{j=1}^M\mu^*_j d_j^s - \sum_{i=1}^N \lambda^*_i C_i x_i \le 0$$

- If the dual subproblem solves to optimality, we obtain an extreme point $(\mu^*, \lambda^*)$ and add an optimality cut:
$$θ_s \ge \sum_{j=1}^M\mu^*_j d_j^s - \sum_{i=1}^N\lambda^*_i C_i x_i$$

We can implement this decomposition as follows:


```julia
"Solve problem using multi-cut Benders decomposition"
function solve_benders_multi(data::FacilityLocationData; verbose::Bool=false)
    # define main problem
    MP = Model(() -> Gurobi.Optimizer(GRB_ENV));
    set_optimizer_attributes(MP, "TimeLimit" => 60, "MIPGap" => 1e-4, "OutputFlag" => 0)
    N = size(data.facilities, 1); M = size(data.customers, 1); S = size(data.demand, 2)
    @variable(MP, x[1:N], Bin)
    @variable(MP, θ[1:S] >= 0)
    @objective(MP, Min, sum(data.facility_cost[i] * x[i] for i=1:N) + sum(data.prob[s] * θ[s] for s in 1:S))

    lower_bound_all = []; upper_bound_all = []
    MP_time = []; SP_max_time = []; SP_time = []
    while true
        # solve master problem
        push!(MP_time, @elapsed optimize!(MP))
        lower_bound_new = objective_value(MP)
        push!(lower_bound_all, lower_bound_new)
        x_MP = value.(MP[:x])
        # solve S subproblems
        obj_SP = zeros(S)
        SP_time_all = zeros(S)
        for s = 1:S
            SP_dual = Model(() -> Gurobi.Optimizer(GRB_ENV))
            set_optimizer_attributes(SP_dual, "OutputFlag" => 0)
            @variable(SP_dual, λ[1:N] >= 0);
            @variable(SP_dual, μ[1:M] >= 0);
            @objective(SP_dual, Max,
                       sum(μ[j] * data.demand[j,s] for j in 1:M) -
                       sum(λ[i] * data.capacity[i] * x_MP[i] for i in 1:N))
            @constraint(SP_dual, [i in 1:N, j in 1:M], μ[j] - λ[i] <= data.transportation_cost[i,j])
            @constraint(SP_dual, [j in 1:M], μ[j] <= data.cost_unmet_demand[j])
            SP_time_all[s] = @elapsed optimize!(SP_dual)
            obj_SP_dual = objective_value(SP_dual)
            λ_val = value.(SP_dual[:λ])
            μ_val = value.(SP_dual[:μ])            
            if termination_status(SP_dual) == MOI.DUAL_INFEASIBLE # feasibility cut
                @constraint(MP, sum(μ_val[j] * data.demand[j, s] for j in 1:M) -
                            sum(λ_val[i] * data.capacity[i] * x[i] for i in 1:N) <= 0)
                obj_SP[s] = 999999999
            elseif termination_status(SP_dual) == MOI.OPTIMAL
                @constraint(MP, θ[s] >= sum(μ_val[j] * data.demand[j,s] for j in 1:M) -
                            sum(λ_val[i] * data.capacity[i] * x[i] for i in 1:N))
                obj_SP[s] = obj_SP_dual
            end
        end
        push!(SP_max_time, maximum(SP_time_all))
        push!(SP_time, sum(SP_time_all))
        upper_bound_new = sum(data.facility_cost[i] * x_MP[i] for i=1:N) + sum(data.prob[s] * obj_SP[s] for s in 1:S)
        push!(upper_bound_all, upper_bound_new)
        verbose && @printf("Sol: %.2f - Bound: %.2f\n", upper_bound_all[end], lower_bound_all[end])
        if sum(MP_time) + sum(SP_time) >= TIME_LIMIT ||
            (upper_bound_new-lower_bound_new)/lower_bound_new < OPTIMALITY_GAP
            break
        end
    end
    return upper_bound_all, lower_bound_all, MP_time, SP_time, SP_max_time
end
```


```julia
@time upper1, lower1, main_time1, subproblem_time1, subproblem_max_time1 = solve_benders_multi(data, verbose = true);
```

We can take a look at the convergence by iteration:


```julia
plot([upper1 lower1], label=["Upper bound" "Lower bound"], xlabel="Iteration")
```

And we can also look at it as a function of time:


```julia
plot(cumsum(main_time1 .+ subproblem_time1), [upper1 lower1], label=["Upper bound" "Lower bound"], xlabel="Time (s)")
```

In our implementation, we solved all the subproblems in series, but they are actually independent. We can estimate the runtime on a fully parallel system by using the maximum subproblem runtime instead of the sum of subproblem runtimes.


```julia
plot(cumsum(main_time1 .+ subproblem_max_time1),
    [upper1 lower1], label=["Upper bound" "Lower bound"], xlabel="Time in parallel (s)")
```

### 3.2. Single-cut approach

An alternative approach is the following formulation:

$$\begin{align}
\min\quad &\sum_{i=1}^N c_i x_i + \theta\\
\text{s.t.}\quad & \mathbf{x}\in\{0,1\}^N, \mathbf{\theta}\ge 0
\end{align}$$

We keep the same $S$ subproblems, but this time we add cuts as follows:

- If **any** dual subproblem is unbounded, we obtain an extreme ray $(\mu^*, \lambda^*)$ and add a feasibility cut:
$$\sum_{j=1}^M\mu^*_j d_j^s - \sum_{i=1}^N \lambda^*_i C_i x_i \le 0$$

- If **all** dual subproblems solve to optimality, we obtain extreme points $(\mu^*_s, \lambda^*_s)$ and add an optimality cut:
$$θ \ge \sum_{s=1}^Sp_s\left(\sum_{j=1}^M\mu^*_j d_j^s - \sum_{i=1}^N\lambda^*_i C_i x_i\right)$$


```julia
"Solve problem using multi-cut Benders decomposition"
function solve_benders_single(data::FacilityLocationData; verbose::Bool=false)
    MP = Model(() -> Gurobi.Optimizer(GRB_ENV));
    set_optimizer_attributes(MP, "TimeLimit" => 60, "MIPGap" => 1e-4, "OutputFlag" => 0)
    N = size(data.facilities, 1); M = size(data.customers, 1); S = size(data.demand, 2)
    @variable(MP, x[1:N], Bin)
    @variable(MP, θ >= 0)
    @objective(MP, Min, sum(data.facility_cost[i] * x[i] for i=1:N) + θ)
    lower_bound_all = []; upper_bound_all = []
    MP_time = []; SP_time = []; SP_max_time = []
    while true
        push!(MP_time, @elapsed optimize!(MP))
        lower_bound_new = objective_value(MP)
        push!(lower_bound_all, lower_bound_new)
        x_MP = value.(MP[:x])
        # set up subproblem loop
        obj_SP = zeros(S)
        SP_time_all = zeros(S)
        λ_all = zeros(N, S)
        μ_all = zeros(M, S)
        optimality_cut = true
        for s in 1:S
            SP_dual = Model(() -> Gurobi.Optimizer(GRB_ENV));
            set_optimizer_attributes(SP_dual, "TimeLimit" => 60, "MIPGap" => 1e-4, "OutputFlag" => 0)
            @variable(SP_dual, λ[1:N] >= 0)
            @variable(SP_dual, μ[1:M] >= 0)
            @objective(SP_dual, Max, sum(μ[j] * data.demand[j,s] for j in 1:M)-
                       sum(λ[i] * data.capacity[i] * x_MP[i] for i in 1:N))
            @constraint(SP_dual, [i in 1:N, j in 1:M], μ[j] - λ[i] <= data.transportation_cost[i, j])
            @constraint(SP_dual, [j in 1:M], μ[j] <= data.cost_unmet_demand[j])
            SP_time_all[s] = @elapsed optimize!(SP_dual)
            obj_SP_dual = objective_value(SP_dual)
            λ_all[:, s] = value.(SP_dual[:λ])
            μ_all[:, s] = value.(SP_dual[:μ])
            if termination_status(SP_dual) == MOI.DUAL_INFEASIBLE # feasibility cut
                optimality_cut = false
                @constraint(MP, sum(μ_val[j] * data.demand[j, s] for j in 1:M) -
                            sum(λ_val[i] * data.capacity[i] * x[i] for i in 1:N) <= 0)
                obj_SP[s] = 999999999
            elseif termination_status(SP_dual) == MOI.OPTIMAL
                obj_SP[s] = obj_SP_dual
            end
        end
        if optimality_cut
            @constraint(MP, θ >= sum(- data.prob[s] * λ_all[i, s] * data.capacity[i] * x[i] for i=1:N, s=1:S) +
                        sum(data.prob[s] * μ_all[j, s] * data.demand[j, s] for j = 1:M, s = 1:S))
        end
        push!(SP_max_time, maximum(SP_time_all))
        push!(SP_time, sum(SP_time_all))
        upper_bound_new = sum(data.facility_cost[i] * x_MP[i] for i=1:N) + sum(data.prob[s] * obj_SP[s] for s in 1:S);
        push!(upper_bound_all, upper_bound_new)
        verbose && @printf("Sol: %.2f - Bound: %.2f\n", upper_bound_all[end], lower_bound_all[end])
        if sum(MP_time) + sum(SP_time) >= TIME_LIMIT ||
            (upper_bound_new - lower_bound_new) / lower_bound_new < OPTIMALITY_GAP
            break
        end
    end
    return upper_bound_all, lower_bound_all, MP_time, SP_time, SP_max_time
end
```

Let's try this on our problem:


```julia
@time upper2, lower2, main_time2, subproblem_time2, subproblem_max_time2 = solve_benders_single(data, verbose = true);
```

And we can plot the progress, as before...


```julia
plot(cumsum(main_time2 .+ subproblem_time2), [upper2 lower2], label=["Upper bound" "Lower bound"], xlabel="Time (s)")
```

### 3.3 Pareto-optimal cuts

In lecture, we also discussed Pareto-optimal cuts, i.e. cuts that provide a tighter approximation to the recourse function. We mentioned that we can compute a Pareto-optimal cut by finding the best cut for a particular "core point" $x^0$. We first solve the dual subproblem:

$$\begin{align}
z^*=\max\quad &  \sum_{j=1}^M \mu_j d_j^s - \sum_{i=1}^N\lambda_iC_ix_i\\
\text{s.t.}\quad & \mu_j -\lambda_i \le t_{ij} & \forall i\in[N], j\in[M]\\
& \mu_j \le q_j& \forall j\in [M]\\
& \mathbf{\mu, \lambda}\ge 0
\end{align}$$

The, given the optimal objective, $x^*$, we solve the problem again with a new objective, i.e.

$$\begin{align}
\max\quad &  \sum_{j=1}^M \mu_j d_j^s - \sum_{i=1}^N\lambda_iC_ix^0_i\\
\text{s.t.}\quad & \mu_j -\lambda_i \le t_{ij} & \forall i\in[N], j\in[M]\\
& \mu_j \le q_j& \forall j\in [M]\\
& \sum_{j=1}^M \mu_j d_j^s - \sum_{i=1}^N\lambda_iC_ix_i=z^*\\
& \mathbf{\mu, \lambda}\ge 0
\end{align}$$

The core point we use is $x^0=\{0.5\}_{i=1}^N$.


```julia
"Solve problem using multi-cut Benders decomposition"
function solve_benders_pareto(data::FacilityLocationData; verbose::Bool=false)
    MP = Model(() -> Gurobi.Optimizer(GRB_ENV));
    set_optimizer_attributes(MP, "TimeLimit" => 60, "MIPGap" => 1e-4, "OutputFlag" => 0)
    N = size(data.facilities, 1); M = size(data.customers, 1); S = size(data.demand, 2)
    @variable(MP, x[1:N], Bin)
    @variable(MP, θ >= 0)
    @objective(MP, Min, sum(data.facility_cost[i] * x[i] for i=1:N) + θ)
    lower_bound_all = []; upper_bound_all = []
    MP_time = []; SP_time = []; SP_max_time = []
    while true
        push!(MP_time, @elapsed optimize!(MP))
        lower_bound_new = objective_value(MP)
        push!(lower_bound_all, lower_bound_new)
        x_MP = value.(MP[:x])
        # set up subproblem loop
        obj_SP = zeros(S)
        SP_time_all = zeros(S)
        λ_all = zeros(N, S)
        μ_all = zeros(M, S)
        optimality_cut = true
        for s in 1:S
            SP_dual = Model(() -> Gurobi.Optimizer(GRB_ENV));
            set_optimizer_attributes(SP_dual, "TimeLimit" => 60, "MIPGap" => 1e-4, "OutputFlag" => 0)
            @variable(SP_dual, λ[1:N] >= 0)
            @variable(SP_dual, μ[1:M] >= 0)
            @objective(SP_dual, Max, sum(μ[j] * data.demand[j,s] for j in 1:M)-
                       sum(λ[i] * data.capacity[i] * x_MP[i] for i in 1:N))
            @constraint(SP_dual, [i in 1:N, j in 1:M], μ[j] - λ[i] <= data.transportation_cost[i, j])
            @constraint(SP_dual, [j in 1:M], μ[j] <= data.cost_unmet_demand[j])
            SP_time_all[s] = @elapsed optimize!(SP_dual)
            obj_SP_dual = objective_value(SP_dual)
            λ_all[:, s] = value.(SP_dual[:λ])
            μ_all[:, s] = value.(SP_dual[:μ])
            if termination_status(SP_dual) == MOI.DUAL_INFEASIBLE # feasibility cut
                optimality_cut = false
                @constraint(MP, sum(μ_val[j] * data.demand[j, s] for j in 1:M) -
                            sum(λ_val[i] * data.capacity[i] * x[i] for i in 1:N) <= 0)
                obj_SP[s] = 999999999
            elseif termination_status(SP_dual) == MOI.OPTIMAL
                obj_SP[s] = obj_SP_dual
                @constraint(SP_dual, sum(μ[j] * data.demand[j,s] for j in 1:M)-
                            sum(λ[i] * data.capacity[i] * x_MP[i] for i in 1:N) == obj_SP[s])
                @objective(SP_dual, Max, sum(μ[j] * data.demand[j,s] for j in 1:M)-
                           sum(λ[i] * data.capacity[i] * 0.5 for i in 1:N))
                SP_time_all[s] += @elapsed optimize!(SP_dual)
                λ_all[:, s] = value.(SP_dual[:λ])
                μ_all[:, s] = value.(SP_dual[:μ])
            end
        end
        if optimality_cut
            @constraint(MP, θ >= sum(- data.prob[s] * λ_all[i, s] * data.capacity[i] * x[i] for i=1:N, s=1:S) +
                        sum(data.prob[s] * μ_all[j, s] * data.demand[j, s] for j = 1:M, s = 1:S))
        end
        push!(SP_max_time, maximum(SP_time_all))
        push!(SP_time, sum(SP_time_all))
        upper_bound_new = sum(data.facility_cost[i] * x_MP[i] for i=1:N) + sum(data.prob[s] * obj_SP[s] for s in 1:S);
        push!(upper_bound_all, upper_bound_new)
        verbose && @printf("Sol: %.2f - Bound: %.2f\n", upper_bound_all[end], lower_bound_all[end])
        if sum(MP_time) + sum(SP_time) >= TIME_LIMIT ||
            (upper_bound_new - lower_bound_new) / lower_bound_new < OPTIMALITY_GAP
            break
        end
    end
    return upper_bound_all, lower_bound_all, MP_time, SP_time, SP_max_time
end
```

Let's solve the problem using Pareto cuts:


```julia
@time upper3, lower3, main_time3, subproblem_time3, subproblem_max_time3 = solve_benders_pareto(data, verbose = true);
```


```julia
plot(cumsum(main_time3 .+ subproblem_time3), [upper3 lower3], label=["Upper bound" "Lower bound"], xlabel="Time (s)")
```

Instead of solving the dual subproblem twice, another way to obtain Pareto-optimal cuts is to add the core point objective with a small weight $\varepsilon$, to the original dual objective.

### 3.4 Benders cuts with lazy constraints

The main drawback of the approaches above is that we have to restart solving the master problem from scratch every time we add the cuts. One way we can change that is by using lazy constraints.

Instead of providing the solver with constraints, we give it a "callback function", i.e., a function which takes in a solution and provides a constraint that it violates. The solver will regularly call this function to see if there are new constraints that need to be added.


```julia
"Solve problem using multi-cut Benders decomposition - verbose can be 0, 1 or 2"
function solve_benders_lazy(data::FacilityLocationData; verbose::Int=0)
    # define main problem
    MP = Model(() -> Gurobi.Optimizer(GRB_ENV));
    set_optimizer_attributes(MP, "TimeLimit" => TIME_LIMIT, "MIPGap" => OPTIMALITY_GAP,
                             "OutputFlag" => ifelse(verbose > 0, 1, 0), "LazyConstraints" => 1)
    N = size(data.facilities, 1); M = size(data.customers, 1); S = size(data.demand, 2)
    @variable(MP, x[1:N], Bin)
    @variable(MP, θ >= 0)
    @objective(MP, Min, sum(data.facility_cost[i] * x[i] for i=1:N) + θ)

    "Define the callback function"
    function add_benders_cuts(cb_data)
        status = callback_node_status(cb_data, MP)
        if status == MOI.CALLBACK_NODE_STATUS_INTEGER
            verbose == 2 && println("Adding lazy constraints...")
            # get value of current solution
            x_MP = [round(callback_value(cb_data, x[i])) for i = 1:N]
            λ_all = zeros(N, S)
            μ_all = zeros(M, S)
            optimality_cut = true
            # solve dual subproblems
            for s = 1:S
                SP_dual = Model(() -> Gurobi.Optimizer(GRB_ENV))
                set_optimizer_attributes(SP_dual, "OutputFlag" => 0)
                @variable(SP_dual, λ[1:N] >= 0);
                @variable(SP_dual, μ[1:M] >= 0);
                @objective(SP_dual, Max,
                           sum(μ[j] * data.demand[j,s] for j in 1:M) -
                           sum(λ[i] * data.capacity[i] * x_MP[i] for i in 1:N))
                @constraint(SP_dual, [i in 1:N, j in 1:M], μ[j] - λ[i] <= data.transportation_cost[i,j])
                @constraint(SP_dual, [j in 1:M], μ[j] <= data.cost_unmet_demand[j])
                optimize!(SP_dual)
                λ_all[:,s] = value.(SP_dual[:λ])
                μ_all[:,s] = value.(SP_dual[:μ])            
                if termination_status(SP_dual) == MOI.DUAL_INFEASIBLE # feasibility cut
                    optimality_cut = false
                    fea = @build_constraint(sum(μ_all[j, s] * data.demand[j, s] for j in 1:M) -
                                            sum(λ_all[i, s] * data.capacity[i] * x[i] for i in 1:N) <= 0)
                    MOI.submit(MP, MOI.LazyConstraint(cb_data), fea)
                elseif termination_status(SP_dual) == MOI.OPTIMAL
                    # do nothing
                end
            end
            opt = @build_constraint(θ >= sum(data.prob[s] * μ_all[j, s] * data.demand[j,s] for j = 1:M, s=1:S) -
                                    sum(data.prob[s] * λ_all[i, s] * data.capacity[i] * x[i] for i = 1:N, s=1:S))
            MOI.submit(MP, MOI.LazyConstraint(cb_data), opt)
        end
    end
    # set callback function and attach to model
    MOI.set(MP, MOI.LazyConstraintCallback(), add_benders_cuts)
    solvetime = @elapsed optimize!(MP)
    opt = objective_value(MP)
    bound = objective_bound(MP)
    return opt, bound, solvetime
end
```

We can solve the problem using lazy constraints:


```julia
solve_benders_lazy(data, verbose=0)
```

Note that here we are using JuMP's solver-independent solver callback syntax. That means we don't have as much control about what to do within the callback, and when to actually call it. For example, when we set `verbose=2` above, we see we generate a lot of lazy constraints even before the end of the presolve. Why might that happen?

If we want more fine-grained control over the solver, we can write a solver-dependent callback function, using the `Gurobi.jl` package directly - see [here](https://www.gurobi.com/documentation/9.1/refman/cb_codes.html) for more details!

Another drawback of using lazy constraints is that it is harder to parallelize, and we do not have full control over when the callback is run. There are many examples where using lazy constraints considerably speeds up the problem. A usual one is the traveling salesman problem, and related vehicle routing problems.

3137a60b7504696ceb78be42b7168ffe55f19a49
ipynb
Jupyter Notebook
Rec9_code/Recitation9-Complete.ipynb
leobix/mio2021
ec28e15b9f701745b3616ee9174d3573fbac0d14
MIT
Rec9_code/Recitation9-Complete.ipynb
leobix/mio2021
ec28e15b9f701745b3616ee9174d3573fbac0d14
MIT
Rec9_code/Recitation9-Complete.ipynb
leobix/mio2021
ec28e15b9f701745b3616ee9174d3573fbac0d14
MIT
Qwen/Qwen-72B
1. YES
2. YES