Skip to content

Commit

Permalink
adding debuging models with power balance slack variables (#55)
Browse files Browse the repository at this point in the history
  • Loading branch information
@ccoffrin
ccoffrin authored Jul 14, 2018
1 parent 9e57447 commit c81080e
Show file tree
Hide file tree
Showing 13 changed files with 396 additions and 3 deletions.
2 changes: 2 additions & 0 deletions src/ThreePhasePowerModels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ const LOGGER = getlogger(PowerModels)
include("core/constraint_template.jl")
include("core/ref.jl")
include("core/variable.jl")
include("core/objective.jl")

include("form/acp.jl")
include("form/dcp.jl")
Expand All @@ -31,5 +32,6 @@ include("prob/tp_opf.jl")
include("prob/tp_opf_bf.jl")
include("prob/tp_ots.jl")
include("prob/tp_pf.jl")
include("prob/tp_debug.jl")

end
26 changes: 26 additions & 0 deletions src/core/constraint_template.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,32 @@
# Three-phase specific constraints


""
function constraint_kcl_shunt_slack(pm::GenericPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
if !haskey(con(pm, nw, cnd), :kcl_p)
con(pm, nw, cnd)[:kcl_p] = Dict{Int,ConstraintRef}()
end
if !haskey(con(pm, nw, cnd), :kcl_q)
con(pm, nw, cnd)[:kcl_q] = Dict{Int,ConstraintRef}()
end

bus = ref(pm, nw, :bus, i)
bus_arcs = ref(pm, nw, :bus_arcs, i)
bus_arcs_dc = ref(pm, nw, :bus_arcs_dc, i)
bus_gens = ref(pm, nw, :bus_gens, i)
bus_loads = ref(pm, nw, :bus_loads, i)
bus_shunts = ref(pm, nw, :bus_shunts, i)

bus_pd = Dict(k => ref(pm, nw, :load, k, "pd", cnd) for k in bus_loads)
bus_qd = Dict(k => ref(pm, nw, :load, k, "qd", cnd) for k in bus_loads)

bus_gs = Dict(k => ref(pm, nw, :shunt, k, "gs", cnd) for k in bus_shunts)
bus_bs = Dict(k => ref(pm, nw, :shunt, k, "bs", cnd) for k in bus_shunts)

constraint_kcl_shunt_slack(pm, nw, cnd, i, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs)
end


""
function constraint_tp_voltage(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
constraint_tp_voltage(pm, nw, cnd)
Expand Down
10 changes: 10 additions & 0 deletions src/core/objective.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
"a quadratic penalty for bus power slack variables"
function objective_min_bus_power_slack(pm::GenericPowerModel)
return @objective(pm.model, Min,
sum(
sum(
sum( var(pm, n, c, :p_slack, i)^2 + var(pm, n, c, :q_slack, i)^2 for (i,bus) in nw_ref[:bus])
for c in conductor_ids(pm, n))
for (n, nw_ref) in nws(pm))
)
end
25 changes: 25 additions & 0 deletions src/core/variable.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,28 @@
"generates variables for both `active` and `reactive` slack at each bus"
function variable_bus_power_slack(pm::GenericPowerModel; kwargs...)
variable_active_bus_power_slack(pm; kwargs...)
variable_reactive_bus_power_slack(pm; kwargs...)
end

""
function variable_active_bus_power_slack(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
var(pm, nw, cnd)[:p_slack] = @variable(pm.model,
[i in ids(pm, nw, :bus)], basename="$(nw)_$(cnd)_p_slack",
start = PMs.getval(ref(pm, nw, :bus, i), "p_slack_start", cnd)
)
end

""
function variable_reactive_bus_power_slack(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
var(pm, nw, cnd)[:q_slack] = @variable(pm.model,
[i in ids(pm, nw, :bus)], basename="$(nw)_$(cnd)_q_slack",
start = PMs.getval(ref(pm, nw, :bus, i), "q_slack_start", cnd)
)
end




"variable: `w[i] >= 0` for `i` in `bus`es"
function variable_tp_voltage_magnitude_sqr(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded=true)
bus_cnd = [(i, c) for i in ids(pm, nw, :bus) for c in PMs.conductor_ids(pm)]
Expand Down
18 changes: 17 additions & 1 deletion src/form/acp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,28 @@ function variable_tp_voltage(pm::GenericPowerModel{T}; kwargs...) where T <: PMs
PMs.variable_voltage_magnitude(pm; kwargs...)
end


"do nothing, this model does not have complex voltage constraints"
function constraint_tp_voltage(pm::GenericPowerModel{T}; kwargs...) where T <: PMs.AbstractACPForm
end


""
function constraint_kcl_shunt_slack(pm::GenericPowerModel{T}, n::Int, c::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs) where T <: PMs.AbstractACPForm
vm = var(pm, n, c, :vm, i)
p_slack = var(pm, n, c, :p_slack, i)
q_slack = var(pm, n, c, :q_slack, i)
p = var(pm, n, c, :p)
q = var(pm, n, c, :q)
pg = var(pm, n, c, :pg)
qg = var(pm, n, c, :qg)
p_dc = var(pm, n, c, :p_dc)
q_dc = var(pm, n, c, :q_dc)

con(pm, n, c, :kcl_p)[i] = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs))*vm^2 + p_slack)
con(pm, n, c, :kcl_q)[i] = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs))*vm^2 + q_slack)
end


"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
Expand Down
18 changes: 18 additions & 0 deletions src/form/shared.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,24 @@
# Three-phase specific constraints


""
function constraint_kcl_shunt_slack(pm::GenericPowerModel{T}, n::Int, c::Int, i, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs) where T <: PMs.AbstractWForms
w = var(pm, n, c, :w, i)
p_slack = var(pm, n, c, :p_slack, i)
q_slack = var(pm, n, c, :q_slack, i)
pg = var(pm, n, c, :pg)
qg = var(pm, n, c, :qg)
p = var(pm, n, c, :p)
q = var(pm, n, c, :q)
p_dc = var(pm, n, c, :p_dc)
q_dc = var(pm, n, c, :q_dc)

@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs))*w + p_slack)
@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs))*w + q_slack)
end



"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
"""
Expand Down
146 changes: 146 additions & 0 deletions src/prob/tp_debug.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
# These problem formulations are used to debug Three Phase datasets
# that do not converge using the standard formulations

export run_tp_opf_pbs, run_tp_pf_pbs

""
function run_tp_opf_pbs(data::Dict{String,Any}, model_constructor, solver; kwargs...)
return PMs.run_generic_model(data, model_constructor, solver, post_tp_opf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
end


""
function run_tp_opf_pbs(file::String, model_constructor, solver; kwargs...)
data = ThreePhasePowerModels.parse_file(file)
return PMs.run_generic_model(data, model_constructor, solver, post_tp_opf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
end


""
function post_tp_opf_pbs(pm::GenericPowerModel)
for c in PMs.conductor_ids(pm)
variable_tp_voltage(pm, cnd=c)
variable_bus_power_slack(pm, cnd=c)
PMs.variable_generation(pm, cnd=c)
PMs.variable_branch_flow(pm, cnd=c)
PMs.variable_dcline_flow(pm, cnd=c)
end

for c in PMs.conductor_ids(pm)
constraint_tp_voltage(pm, cnd=c)

for i in ids(pm, :ref_buses)
constraint_tp_theta_ref(pm, i, cnd=c)
end

for i in ids(pm, :bus)
constraint_kcl_shunt_slack(pm, i, cnd=c)
end

for i in ids(pm, :branch)
constraint_ohms_tp_yt_from(pm, i, cnd=c)
constraint_ohms_tp_yt_to(pm, i, cnd=c)

PMs.constraint_voltage_angle_difference(pm, i, cnd=c)

PMs.constraint_thermal_limit_from(pm, i, cnd=c)
PMs.constraint_thermal_limit_to(pm, i, cnd=c)
end

for i in ids(pm, :dcline)
PMs.constraint_dcline(pm, i, cnd=c)
end
end

objective_min_bus_power_slack(pm)
end



""
function run_tp_pf_pbs(data::Dict{String,Any}, model_constructor, solver; kwargs...)
return PMs.run_generic_model(data, model_constructor, solver, post_tp_pf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
end


""
function run_tp_pf_pbs(file::String, model_constructor, solver; kwargs...)
data = ThreePhasePowerModels.parse_file(file)
return PMs.run_generic_model(data, model_constructor, solver, post_tp_pf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
end

""
function post_tp_pf_pbs(pm::GenericPowerModel)
for c in PMs.conductor_ids(pm)
variable_tp_voltage(pm, bounded=false, cnd=c)
variable_bus_power_slack(pm, cnd=c)
PMs.variable_generation(pm, bounded=false, cnd=c)
PMs.variable_branch_flow(pm, bounded=false, cnd=c)
PMs.variable_dcline_flow(pm, bounded=false, cnd=c)
end

for c in PMs.conductor_ids(pm)
constraint_tp_voltage(pm, cnd=c)

for (i,bus) in ref(pm, :ref_buses)
@assert bus["bus_type"] == 3
constraint_tp_theta_ref(pm, i, cnd=c)
PMs.constraint_voltage_magnitude_setpoint(pm, i, cnd=c)
end

for (i,bus) in ref(pm, :bus)
constraint_kcl_shunt_slack(pm, i, cnd=c)

# PV Bus Constraints
if length(ref(pm, :bus_gens, i)) > 0 && !(i in ids(pm,:ref_buses))
# this assumes inactive generators are filtered out of bus_gens
@assert bus["bus_type"] == 2

PMs.constraint_voltage_magnitude_setpoint(pm, i, cnd=c)
for j in ref(pm, :bus_gens, i)
PMs.constraint_active_gen_setpoint(pm, j, cnd=c)
end
end
end

for i in ids(pm, :branch)
constraint_ohms_tp_yt_from(pm, i, cnd=c)
constraint_ohms_tp_yt_to(pm, i, cnd=c)
end

for (i,dcline) in ref(pm, :dcline)
PMs.constraint_active_dcline_setpoint(pm, i, cnd=c)

f_bus = ref(pm, :bus)[dcline["f_bus"]]
if f_bus["bus_type"] == 1
PMs.constraint_voltage_magnitude_setpoint(pm, f_bus["index"], cnd=c)
end

t_bus = ref(pm, :bus)[dcline["t_bus"]]
if t_bus["bus_type"] == 1
PMs.constraint_voltage_magnitude_setpoint(pm, t_bus["index"], cnd=c)
end
end

end

objective_min_bus_power_slack(pm)
end


""
function get_pbs_solution(pm::GenericPowerModel, sol::Dict{String,Any})
PMs.add_bus_voltage_setpoint(sol, pm)
PMs.add_generator_power_setpoint(sol, pm)
PMs.add_branch_flow_setpoint(sol, pm)
add_bus_slack_setpoint(sol, pm)
end


""
function add_bus_slack_setpoint(sol, pm::GenericPowerModel)
PMs.add_setpoint(sol, pm, "bus", "p_slack", :p_slack)
PMs.add_setpoint(sol, pm, "bus", "q_slack", :q_slack)
end


2 changes: 1 addition & 1 deletion test/data/matlab/case5_c_m_a.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
tppmc.gen = [
1 0.0 40.0 -30.0 30.0 0.0 40.0 -30.0 30.0 0.0 40.0 -30.0 30.0 40.0 30.0 40.0 30.0 40.0 30.0 1;
1 0.0 170.0 -127.5 127.5 0.0 170.0 -127.5 127.5 0.0 170.0 -127.5 127.5 170.0 127.5 170.0 127.5 170.0 127.5 1;
3 0.0 520.0 -390. 390.0 0.0 520.0 -390.0 390.0 0.0 520.0 -390.0 390.0 325.0 390.0 325.0 390.0 325.0 390.0 1;
3 0.0 520.0 -390.0 390.0 0.0 520.0 -390.0 390.0 0.0 520.0 -390.0 390.0 325.0 390.0 325.0 390.0 325.0 390.0 1;
4 0.0 200.0 -150.0 150.0 0.0 200.0 -150.0 150.0 0.0 200.0 -150.0 150.0 0.0 -10.0 0.0 -10.0 0.0 -10.0 1;
10 0.0 600.0 -450.0 450.0 0.0 600.0 -450.0 450.0 0.0 600.0 -450.0 450.0 470.0 -165.0 470.0 -165.0 470.0 -165.0 1;
];
Expand Down
62 changes: 62 additions & 0 deletions test/data/matlab/case5_c_m_b.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
% 3 Coupled meshed networks with off diagonal terms in branch Y
% tests that opf_pbs model converges on an otherwise infeasible model

function tppmc = case5_c_m_b
tppmc.version = '1'

tppmc.baseMVA = 100.0;
tppmc.baseKV = 230.0;

%% bus data
% bus_i type vmin_1 vmax_1 vmin_2 vmax_2 vmin_3 vmax_3 vm_1 va_1 vm_2 va_2 vm_3 va_3
tppmc.bus = [
1 2 0.90 1.10 0.90 1.10 0.90 1.10 1.00000 2.80377 1.00000 2.80377 1.00000 2.80377;
2 1 0.90 1.10 0.90 1.10 0.90 1.10 1.08407 -0.73465 1.08407 -0.73465 1.08407 -0.73465;
3 2 0.90 1.10 0.90 1.10 0.90 1.10 1.00000 -0.55972 1.00000 -0.55972 1.00000 -0.55972;
4 1 0.90 1.10 0.90 1.10 0.90 1.10 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000;
10 3 0.90 1.10 0.90 1.10 0.90 1.10 1.00000 3.59033 1.00000 3.59033 1.00000 3.59033;
];

%% load data
% load_bus pd_1 qd_1 pd_2 qd_2 pd_3 qd_3 status
tppmc.load = [
2 850.0 100.0 300.0 100.0 300.0 100.0 1;
3 300.0 100.0 300.0 300.0 300.0 500.0 1;
4 400.0 130.0 1150.0 130.0 400.0 130.0 1;
];

%% shunt data
% shunt_bus gs_1 bs_1 gs_2 bs_2 gs_3 bs_3 status
tppmc.shunt = [
4 2.0 50.0 2.0 50.0 2.0 50.0 1;
10 1.0 25.0 1.0 25.0 1.0 25.0 1;
];

%% generator data
% gen_bus pmin_1 pmax_1 qmin_1 qmax_1 pmin_2 pmax_2 qmin_2 qmax_2 pmin_3 pmax_3 qmin_3 qmax_3 pg_1 qg_1 pg_2 qg_2 pg_3 qg_3 status
tppmc.gen = [
1 0.0 40.0 -30.0 30.0 0.0 40.0 -30.0 30.0 0.0 40.0 -30.0 30.0 40.0 30.0 40.0 30.0 40.0 30.0 1;
1 0.0 170.0 -127.5 127.5 0.0 170.0 -127.5 127.5 0.0 170.0 -127.5 127.5 170.0 127.5 170.0 127.5 170.0 127.5 1;
3 0.0 520.0 -190.0 190.0 0.0 520.0 -190.0 190.0 0.0 520.0 -190.0 190.0 325.0 390.0 325.0 390.0 325.0 390.0 1;
10 0.0 600.0 -450.0 450.0 0.0 600.0 -450.0 450.0 0.0 600.0 -450.0 450.0 470.0 -165.0 470.0 -165.0 470.0 -165.0 1;
];

%% generator cost data
% 2 startup shutdown n c(n-1) ... c0
tppmc.gencost = [
2 0.0 0.0 3 0.0 14.0 0.0;
2 0.0 0.0 3 0.0 15.0 0.0;
2 0.0 0.0 3 0.0 30.0 0.0;
2 0.0 0.0 3 0.0 10.0 0.0;
];

%% branch data
% f_bus t_bus r_11 x_11 r_12 x_12 r_13 x_13 r_22 x_22 r_23 x_23 r_33 x_33 b_1 b_2 b_3 rate_a rate_b rate_c angmin angmax status
tppmc.branch = [
1 2 0.00281 0.0281 0.0001 0.0002 0.0001 0.0002 0.00581 0.0281 0.0001 0.0002 0.00281 0.0281 0.0 0.0 0.0 400 400 400 -30.0 30.0 1;
1 4 0.00304 0.0304 0.0001 0.0002 0.0001 0.0002 0.01004 0.0304 0.0001 0.0002 0.00304 0.0304 0.0 0.0 0.0 426 426 426 -30.0 30.0 1;
1 10 0.00064 0.0064 0.0001 0.0002 0.0001 0.0002 0.00164 0.0064 0.0001 0.0002 0.00064 0.0064 0.0 0.0 0.0 426 426 426 -30.0 30.0 1;
2 3 0.00108 0.0108 0.0001 0.0002 0.0001 0.0002 0.00208 0.0108 0.0001 0.0002 0.00108 0.0108 0.0 0.0 0.0 426 426 426 -30.0 30.0 1;
4 3 0.00297 0.0297 0.0001 0.0002 0.0001 0.0002 0.00597 0.0297 0.0001 0.0002 0.00297 0.0297 0.0 0.0 0.0 426 426 426 -30.0 30.0 1;
4 10 0.00297 0.0297 0.0001 0.0002 0.0001 0.0002 0.00497 0.0297 0.0001 0.0002 0.00297 0.0297 0.0 0.0 0.0 240 240 240 -30.0 30.0 1;
];
3 changes: 3 additions & 0 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,4 +40,7 @@ pajarito_solver = PajaritoSolver(mip_solver=cbc_solver, cont_solver=ipopt_solver
include("tp_pf.jl")

include("tp_opf_bf.jl")

include("tp_debug.jl")

end
Loading

0 comments on commit c81080e

Please sign in to comment.