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_kirchh_p = (∑_p_loads + ∑_pcs) * shadow(con_ref[:ι][i]) - (∑_p_loads + ∑_pcs) * shadow(con_ref[:κ][i])
expr_kirchh_q = ∑_q_loads * shadow(con_ref[:λ][i]) - ∑_q_loads * shadow(con_ref[:μ][i])
expr_kirchhoff = expr_kirchhoff + expr_kirchh_p + expr_kirchh_q
current_violation = current_violation + (∑_p_loads + current_pcs) * shadow(con_ref[:ι][i]) - (∑_p_loads + current_pcs) * shadow(con_ref[:κ][i]) + 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"]
expr_voltage_angle = expr_voltage_angle + angmax * shadow(con_ref[:ξ][f_idx]) + angmin * shadow(con_ref[:𝚶][f_idx])
end
return expr_gen_bounds + expr_voltage_bounds + expr_kirchhoff + expr_voltage_angle
end
shadow(con_ref::JuMP.ConstraintRef) = JuMP.shadow_price(con_ref)