# Benders Subproblem Infeasibility

Hello everyone,

I’m currently working on an optimization problem where I’m using Benders decomposition to determine the placement of EV charging stations. The master problem determines the placement of charging stations, while the subproblem involves optimal power flow (using the PowerModels.jl datastructure). I’ve encountered an issue with subproblem infeasibility despite using feasibility cuts. I checked my formulation multiple times and I think the problem lies within my code. For testing purposes I’m not using the modern Benders implementation with lazy constraints yet.

I solve the master problem to determine integer values x_{i,t}, which denote the number of charging stations of type t at bus i. I use these when I build my linear approximated OPF (primal) subproblem.
When the subproblem is infeasible, I retrieve the dual extreme ray values with JuMP.shadow_price().
I multiply these with the RHS of the constraints of the subproblem in canonical form (each equality constraint is transformed to two inequality constraints).

Adding this particular feasibility cut doesn’t cut away the infeasible solution from the master problem (the master problem solution stays the same after multiple iterations). Even more, if I insert the solution from the master into the feasibility cut, this is not violated. This is the mathematical formulation of the cut:

\sum_{k\in G}\left( P_{k}^{gl} \overline{\gamma_{k}} -P_{k}^{gu} \overline{\psi_{k}} +Q_{k}^{gl} \overline{\epsilon_{k}} -Q_{k}^{gu} \overline{\zeta_{k}}\right) +\sum_{i\in N}\left( v_{i}^{l} \overline{ \eta_{i}} - v_{i}^{u} \overline{\tau_{i}}\right) \\ +\sum_{i\in N}\left(\left(\sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs} x_{i,t} \right)\overline{ \iota_{i}} -\left(\sum_{k\in L_{i}} P_{k}^{d} +\sum_{t\in T} P_{t}^{cs} x_{i,t}\right) \overline{\kappa_{i}}\right) \\+\sum_{i\in N}\left(\sum_{k\in L_{i}} Q_{k}^{d} \overline{\lambda_{i}} -\sum_{k\in L_{i}} Q_{k}^{d} \overline{\mu_{i}}\right) +\sum_{( i,j) \in E}\left( \theta_{ij}^{\Delta l} \overline{\omicron_{ij}} -\theta_{ij}^{\Delta u} \overline{\xi_{ij}}\right) \leq 0

I know there is a feasible solution as I can solve the instance by using the full formulation.
All the extreme ray values are 0, expect some values for \iota_i that are -1. But these values are incorrect as this is not even a feasible solution to the dual of the subproblem (dual values must be positive).

I’m wondering if I’m doing something wrong here in my code, as I’m new to Julia/JuMP.
Any feedback or help would be appreciated!

function solve(...)
# Preprocessing ...

model_data = benders(ps, evcs, transport, optimizer)
return model_data
end

function benders(ps, evcs, transport, optimizer)
master_model = build_master_relaxed(ps, evcs, transport, optimizer)
lower_bound = -Inf
upper_bound = Inf
gap = Inf

for k in 1:100
if gap < 0.0001
break
end

JuMP.optimize!(master_model.model)

if JuMP.termination_status(master_model.model) != _MOI.OPTIMAL
println("Master problem is $(JuMP.termination_status(master_model.model))") break end lower_bound = JuMP.objective_value(master_model.model) x = JuMP.value.(master_model.model[:x]) subprob, con_ref = build_sub(ps, evcs, x, optimizer) JuMP.optimize!(subprob.model) status = JuMP.termination_status(subprob.model) if status == _MOI.OPTIMAL # optimimality cut ... (not included here) upper_bound = JuMP.objective_value(master_model.model) + JuMP.objective_value(subprob.model) gap = (upper_bound - lower_bound) / upper_bound else # feasiblity cut feasibility_cut = gen_feasibility_cut(ps, evcs, master_model.model, con_ref, x) JuMP.@constraint(master_model.model, feasibility_cut <= 0) end end end function build_master(ps, evcs, transport, optimizer) model = JuMP.Model(optimizer) model_data = ModelData(model, ps, evcs, transport) # Building master problem... return model_data end function build_sub(ps, evcs, x, optimizer) subproblem = JuMP.Model(optimizer) subproblem_data = ModelData(subproblem, ps, evcs, transport) con_ref = create_constraint_ref() model = subproblem_data.model for gen in keys(ps[:gen]) JuMP.@variable(model, pg[gen]) JuMP.@variable(model, qg[gen]) con_ref[:γ][gen] = JuMP.@constraint(model, pg[gen] >= ps[:gen][gen]["pmin"]) con_ref[:ψ][gen] = JuMP.@constraint(model, -pg[gen] >= -ps[:gen][gen]["pmax"]) con_ref[:ϵ][gen] = JuMP.@constraint(model, qg[gen] >= ps[:gen][gen]["qmin"]) con_ref[:ζ][gen] = JuMP.@constraint(model, -qg[gen] >= -ps[:gen][gen]["qmax"]) end JuMP.@variable(model, -ps[:branch][l]["rate_a"] <= p[(l, i, j) in ps[:arcs]] <= ps[:branch][l]["rate_a"]) JuMP.@variable(model, -ps[:branch][l]["rate_a"] <= q[(l, i, j) in ps[:arcs]] <= ps[:branch][l]["rate_a"]) JuMP.@variable(model, va[i in keys(ps[:bus])]) JuMP.@variable(model, vm[i in keys(ps[:bus])], start=1.0) # Kirchhoff constraints for bus in keys(ps[:bus]) con_ref[:η][bus] = JuMP.@constraint(model, vm[bus] >= ps[:bus][bus]["vmin"]) con_ref[:τ][bus] = JuMP.@constraint(model, -vm[bus] >= -ps[:bus][bus]["vmax"]) bus_loads = [ps[:load][l] for l in ps[:bus_loads][bus]] con_ref[:ι][bus] = JuMP.@constraint(model, sum(var(model, "pg[$g]") for g in ps[:bus_gens][bus]; init=0.0) -
sum(var(model, "p[$arc]") for arc in ps[:bus_arcs][bus]) >= sum(load["pd"] for load in bus_loads; init=0.0) + sum(evcs[:p_cs][j] * x[bus, j] for j in 1:evcs[:no_types_cs]) # charging station power output times number of charging stations of that type ) con_ref[:κ][bus] = JuMP.@constraint(model, -sum(var(model, "pg[$g]") for g in ps[:bus_gens][bus]; init=0.0) +
sum(var(model, "p[$arc]") for arc in ps[:bus_arcs][bus]) >= -sum(load["pd"] for load in bus_loads; init=0.0) - sum(evcs[:p_cs][j] * x[bus, j] for j in 1:evcs[:no_types_cs]) ) con_ref[:λ][bus] = JuMP.@constraint(model, sum(var(model, "qg[$g]") for g in ps[:bus_gens][bus]; init=0.0) -
sum(var(model, "q[$arc]") for arc in ps[:bus_arcs][bus]) >= sum(load["qd"] for load in bus_loads; init=0.0) ) con_ref[:μ][bus] = JuMP.@constraint(model, -sum(var(model, "qg[$g]") for g in ps[:bus_gens][bus]; init=0.0) +
sum(var(model, "q[$arc]") for arc in ps[:bus_arcs][bus]) >= -sum(load["qd"] for load in bus_loads; init=0.0) ) end for (i, branch) in ps[:branch] branch_f = branch["f_bus"] branch_t = branch["t_bus"] idx = (i, branch_f, branch_t) va_fr = var(model, "va[$branch_f]")
va_to = var(model, "va[$branch_t]") con_ref[:ξ][idx] = JuMP.@constraint(model, va_fr - va_to >= branch["angmin"]) con_ref[:𝚶][idx] = JuMP.@constraint(model, va_to - va_fr >= -branch["angmax"]) end # linearized ohm constraints for (i, branch) in ps[:branch] f_idx = (i, branch["f_bus"], branch["t_bus"]) t_idx = (i, branch["t_bus"], branch["f_bus"]) p_f = var(model, "p[$(f_idx)]")
p_t = var(model, "p[$(t_idx)]") q_f = var(model, "q[$(f_idx)]")
q_t = var(model, "q[$(t_idx)]") vm_f = var(model, "vm[$(branch["f_bus"])]")
vm_t = var(model, "vm[$(branch["t_bus"])]") va_f = var(model, "va[$(branch["f_bus"])]")
va_t = var(model, "va[$(branch["t_bus"])]") g, b = _PM.calc_branch_y(branch) # Active power con_ref[:ρ][f_idx] = JuMP.@constraint(model, p_f - (vm_f - vm_t) * g + (va_f - va_t) * b >= 0) con_ref[:σ][f_idx] = JuMP.@constraint(model, -p_f + (vm_f - vm_t) * g - (va_f - va_t) * b >= 0) con_ref[:ρ][t_idx] = JuMP.@constraint(model, p_t - (vm_t - vm_f) * g + (va_t - va_f) * b >= 0) con_ref[:σ][t_idx] = JuMP.@constraint(model, -p_t + (vm_t - vm_f) * g - (va_t - va_f) * b >= 0) # Reactive power con_ref[:ν][f_idx] = JuMP.@constraint(model, q_f - (vm_t - vm_f) * b + (va_f - va_t) * g >= 0) con_ref[:ω][f_idx] = JuMP.@constraint(model, -q_f + (vm_t - vm_f) * b - (va_f - va_t) * g >= 0) con_ref[:ν][t_idx] = JuMP.@constraint(model, q_t - (vm_f - vm_t) * b + (va_t - va_f) * g >= 0) con_ref[:ω][t_idx] = JuMP.@constraint(model, -q_t + (vm_f - vm_t) * b - (va_t - va_f) * g >= 0) end # Add opf objective (linear objective) expr_min_fuel_lin = objective_min_fuel_lin(subproblem_data) JuMP.@objective(subproblem_data.model, Min, expr_min_fuel_lin) return subproblem_data, con_ref end function create_constraint_ref() constraint_ref = Dict{Symbol, Any}() constraint_ref[:γ] = Dict{Any, Any}() constraint_ref[:ψ] = Dict{Any, Any}() constraint_ref[:ϵ] = Dict{Any, Any}() constraint_ref[:ζ] = Dict{Any, Any}() constraint_ref[:η] = Dict{Any, Any}() constraint_ref[:τ] = Dict{Any, Any}() constraint_ref[:ξ] = Dict{Any, Any}() constraint_ref[:𝚶] = Dict{Any, Any}() constraint_ref[:ι] = Dict{Any, Any}() constraint_ref[:κ] = Dict{Any, Any}() constraint_ref[:λ] = Dict{Any, Any}() constraint_ref[:μ] = Dict{Any, Any}() constraint_ref[:ρ] = Dict{Any, Any}() constraint_ref[:σ] = Dict{Any, Any}() constraint_ref[:ν] = Dict{Any, Any}() constraint_ref[:ω] = Dict{Any, Any}() return constraint_ref end function gen_feasibility_cut(ps, evcs, master_model, con_ref, x) current_violation = 0 expr_gen_bounds = sum( ps[:gen][i]["pmin"] * shadow(con_ref[:γ][i]) - ps[:gen][i]["pmax"] * shadow(con_ref[:ψ][i]) + ps[:gen][i]["qmin"] * shadow(con_ref[:ϵ][i]) - ps[:gen][i]["qmax"] * shadow(con_ref[:ζ][i]) for i in keys(ps[:gen]) ) expr_voltage_bounds = sum( ps[:bus][i]["vmin"] * shadow(con_ref[:η][i]) - ps[:bus][i]["vmax"] * shadow(con_ref[:τ][i]) for i in keys(ps[:bus]) ) expr_kirchhoff = 0 for (i, bus) in ps[:bus] bus_loads = [ps[:load][l] for l in ps[:bus_loads][i]] ∑_p_loads = sum(load["pd"] for load in bus_loads; init = 0.0) ∑_q_loads = sum(load["qd"] for load in bus_loads; init = 0.0) ∑_pcs = sum(evcs[:p_cs][t] * var(master_model, "x[$i,\$t]") for t in 1:evcs[:no_types_cs] if x[i, t] > 0.0; init=0.0)

expr_kirchhoff = expr_kirchhoff + expr_kirchh_p + expr_kirchh_q
end

expr_voltage_angle = 0
for (i, branch) in ps[:branch]
f_bus = branch["f_bus"]
t_bus = branch["t_bus"]
f_idx = (i, f_bus, t_bus)
angmin = branch["angmin"]
angmax = branch["angmax"]
end

return expr_gen_bounds + expr_voltage_bounds + expr_kirchhoff + expr_voltage_angle
end



Hi there, welcome to the forum.

For the dual ray, you need to use dual instead of shadow_price.

Here’s an example:

julia> using JuMP, HiGHS

julia> model = Model(HiGHS.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: HiGHS

julia> set_attribute(model, "presolve", "off")

julia> @variable(model, x[1:2] >= 0)
2-element Vector{VariableRef}:
x
x

julia> @objective(model, Min, -sum(x))
-x - x

julia> @constraint(model, c1, x <= 1)
c1 : x ≤ 1

julia> @constraint(model, c2, sum(x) >= 3)
c2 : x + x ≥ 3

julia> @constraint(model, c3, x <= 1)
c3 : x ≤ 1

julia> optimize!(model)
Running HiGHS 1.5.3 [date: 1970-01-01, git hash: 45a127b78]
Copyright (c) 2023 HiGHS under MIT licence terms
Solving LP without presolve or with basis
Using EKK dual simplex solver - serial
Iteration        Objective     Infeasibilities num(sum)
0    -1.9999961998e+00 Ph1: 3(3); Du: 2(2) 0s
3    -1.9999961998e+00 0s
Model   status      : Infeasible
Simplex   iterations: 3
Objective value     : -2.0000000000e+00
HiGHS run time      :          0.00

julia> dual_status(model)
INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4

julia> cons = [c1, c2, c3]
3-element Vector{ConstraintRef{Model, C, ScalarShape} where C}:
c1 : x ≤ 1
c2 : x + x ≥ 3
c3 : x ≤ 1

julia> ray = dual.(cons)
3-element Vector{Float64}:
-1.0
1.0
-1.0


shadow_price is really just a helper function for querying the dual solution of an LP at it is taught in some textbooks. See Constraints · JuMP.