From cd6f9860c179c8225224fbf41617cec6418f045e Mon Sep 17 00:00:00 2001 From: Carleton Coffrin Date: Fri, 11 Aug 2023 14:27:05 -0600 Subject: [PATCH] add pst only opf problem spec --- CHANGELOG.md | 2 +- src/core/constraint_template.jl | 36 ++++++++++++ src/form/acp.jl | 30 ++++++++++ src/form/dcp.jl | 10 +--- src/prob/test.jl | 44 ++++++++++++++- test/opf-var.jl | 97 +++++++++++++++++++++------------ 6 files changed, 173 insertions(+), 46 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b0cf17e49..fbffbac20 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ PowerModels.jl Change Log ========================= ### Staged -- Add support for DCP formulation for OPF with PST variables (#874) +- Add support for ACP and DCP formulation in OPF with PST variables (#543,#874) - Fix implementation of `calc_theta_delta_bounds` when conductor parameter is used (#870) ### v0.19.9 diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 3072acccb..96c298e9c 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -481,6 +481,42 @@ function constraint_ohms_y_oltc_pst_to(pm::AbstractPowerModel, i::Int; nw::Int=n end + +"" +function constraint_ohms_y_pst_from(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default) + branch = ref(pm, nw, :branch, i) + f_bus = branch["f_bus"] + t_bus = branch["t_bus"] + f_idx = (i, f_bus, t_bus) + t_idx = (i, t_bus, f_bus) + + g, b = calc_branch_y(branch) + g_fr = branch["g_fr"] + b_fr = branch["b_fr"] + tm = branch["tap"] + + constraint_ohms_y_pst_from(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tm) +end + + +"" +function constraint_ohms_y_pst_to(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default) + branch = ref(pm, nw, :branch, i) + f_bus = branch["f_bus"] + t_bus = branch["t_bus"] + f_idx = (i, f_bus, t_bus) + t_idx = (i, t_bus, f_bus) + + g, b = calc_branch_y(branch) + g_to = branch["g_to"] + b_to = branch["b_to"] + tm = branch["tap"] + + constraint_ohms_y_pst_to(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tm) +end + + + "" function constraint_voltage_drop(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default) branch = ref(pm, nw, :branch, i) diff --git a/src/form/acp.jl b/src/form/acp.jl index ddddf85a2..2d9820e96 100644 --- a/src/form/acp.jl +++ b/src/form/acp.jl @@ -485,6 +485,36 @@ function constraint_ohms_y_oltc_pst_to(pm::AbstractACPModel, n::Int, f_bus, t_bu JuMP.@NLconstraint(pm.model, q_to == -(b+b_to)*vm_to^2 - -b/tm*(vm_to*vm_fr*cos(va_to-va_fr+ta)) + -g/tm*(vm_to*vm_fr*sin(va_to-va_fr+ta)) ) end + +"" +function constraint_ohms_y_pst_from(pm::AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tm) + p_fr = var(pm, n, :p, f_idx) + q_fr = var(pm, n, :q, f_idx) + vm_fr = var(pm, n, :vm, f_bus) + vm_to = var(pm, n, :vm, t_bus) + va_fr = var(pm, n, :va, f_bus) + va_to = var(pm, n, :va, t_bus) + ta = var(pm, n, :ta, f_idx[1]) + + JuMP.@NLconstraint(pm.model, p_fr == (g+g_fr)/tm^2*vm_fr^2 + (-g)/tm*(vm_fr*vm_to*cos(va_fr-va_to-ta)) + (-b)/tm*(vm_fr*vm_to*sin(va_fr-va_to-ta)) ) + JuMP.@NLconstraint(pm.model, q_fr == -(b+b_fr)/tm^2*vm_fr^2 - (-b)/tm*(vm_fr*vm_to*cos(va_fr-va_to-ta)) + (-g)/tm*(vm_fr*vm_to*sin(va_fr-va_to-ta)) ) +end + +"" +function constraint_ohms_y_pst_to(pm::AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tm) + p_to = var(pm, n, :p, t_idx) + q_to = var(pm, n, :q, t_idx) + vm_fr = var(pm, n, :vm, f_bus) + vm_to = var(pm, n, :vm, t_bus) + va_fr = var(pm, n, :va, f_bus) + va_to = var(pm, n, :va, t_bus) + ta = var(pm, n, :ta, f_idx[1]) + + JuMP.@NLconstraint(pm.model, p_to == (g+g_to)*vm_to^2 + -g/tm*(vm_to*vm_fr*cos(va_to-va_fr+ta)) + -b/tm*(vm_to*vm_fr*sin(va_to-va_fr+ta)) ) + JuMP.@NLconstraint(pm.model, q_to == -(b+b_to)*vm_to^2 - -b/tm*(vm_to*vm_fr*cos(va_to-va_fr+ta)) + -g/tm*(vm_to*vm_fr*sin(va_to-va_fr+ta)) ) +end + + "`angmin <= z_branch[i]*(t[f_bus] - t[t_bus]) <= angmax`" function constraint_voltage_angle_difference_on_off(pm::AbstractACPModel, n::Int, f_idx, angmin, angmax, vad_min, vad_max) i, f_bus, t_bus = f_idx diff --git a/src/form/dcp.jl b/src/form/dcp.jl index bc6eb00dd..6fd1520ce 100644 --- a/src/form/dcp.jl +++ b/src/form/dcp.jl @@ -14,12 +14,6 @@ function variable_bus_voltage_magnitude(pm::AbstractDCPModel; nw::Int=nw_id_defa report && sol_component_fixed(pm, nw, :bus, :vm, ids(pm, nw, :bus), 1.0) end -"" -function variable_branch_transform_magnitude(pm::AbstractDCPModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true) - tm = [branch["tap"] for (b,branch) in ref(pm, nw, :branch)] - report && sol_component_value(pm, nw, :branch, :tm, ids(pm, nw, :branch), tm) -end - "" function sol_data_model!(pm::AbstractDCPModel, solution::Dict) @@ -142,7 +136,7 @@ function constraint_ohms_yt_from(pm::AbstractDCMPPModel, n::Int, f_bus, t_bus, f end "" -function constraint_ohms_y_oltc_pst_from(pm::AbstractDCPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr) +function constraint_ohms_y_pst_from(pm::AbstractDCPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tm) p_fr = var(pm, n, :p, f_idx) va_fr = var(pm, n, :va, f_bus) va_to = var(pm, n, :va, t_bus) @@ -152,7 +146,7 @@ function constraint_ohms_y_oltc_pst_from(pm::AbstractDCPModel, n::Int, f_bus, t_ end "" -function constraint_ohms_y_oltc_pst_to(pm::AbstractDCPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to) +function constraint_ohms_y_pst_to(pm::AbstractDCPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tm) p_to = var(pm, n, :p, t_idx) va_fr = var(pm, n, :va, f_bus) va_to = var(pm, n, :va, t_bus) diff --git a/src/prob/test.jl b/src/prob/test.jl index d03ac0ccc..11aa5bbec 100644 --- a/src/prob/test.jl +++ b/src/prob/test.jl @@ -494,7 +494,49 @@ function _build_opf_strg_mi(pm::AbstractPowerModel) end -"opf with tap magnitude and angle as optimization variables" +"opf with transformer angles as optimization variables" +function _solve_opf_pst(file, model_type::Type, optimizer; kwargs...) + return solve_model(file, model_type, optimizer, _build_opf_pst; kwargs...) +end + +"" +function _build_opf_pst(pm::AbstractPowerModel) + variable_bus_voltage(pm) + variable_gen_power(pm) + + variable_branch_transform_angle(pm) + variable_branch_power(pm) + variable_dcline_power(pm) + + objective_min_fuel_and_flow_cost(pm) + + constraint_model_voltage(pm) + + for i in ids(pm, :ref_buses) + constraint_theta_ref(pm, i) + end + + for i in ids(pm, :bus) + constraint_power_balance(pm, i) + end + + for i in ids(pm, :branch) + constraint_ohms_y_pst_from(pm, i) + constraint_ohms_y_pst_to(pm, i) + + constraint_voltage_angle_difference(pm, i) + + constraint_thermal_limit_from(pm, i) + constraint_thermal_limit_to(pm, i) + end + + for i in ids(pm, :dcline) + constraint_dcline_power_losses(pm, i) + end +end + + +"opf with transformer tap magnitude and angle as optimization variables" function _solve_opf_oltc_pst(file, model_type::Type, optimizer; kwargs...) return solve_model(file, model_type, optimizer, _build_opf_oltc_pst; kwargs...) end diff --git a/test/opf-var.jl b/test/opf-var.jl index 39bf9d383..8db35bf9f 100644 --- a/test/opf-var.jl +++ b/test/opf-var.jl @@ -887,104 +887,129 @@ end end end -@testset "test opf with optimization of oltc and pst" begin + +@testset "test opf with optimization of pst" begin @testset "test ac polar opf" begin - @testset "3-bus case with fixed phase shift / tap" begin + + @testset "3-bus case with optimal phase shifting" begin file = "../test/data/matpower/case3_oltc_pst.m" data = PowerModels.parse_file(file) - result = PowerModels.run_opf(data, ACPPowerModel, nlp_solver) + result = PowerModels._solve_opf_pst(data, ACPPowerModel, nlp_solver) @test result["termination_status"] == LOCALLY_SOLVED - @test isapprox(result["objective"], 5820.1; atol = 1e0) + @test isapprox(result["objective"], 5755.3; atol = 1e0) + + @test haskey(result["solution"]["branch"]["1"], "ta") + + @test isapprox(result["solution"]["branch"]["1"]["ta"], 0.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["2"]["ta"], 0.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["3"]["ta"], 15.0/180*pi; atol = 1e-1) end - @testset "3-bus case with optimal phase shifting / tap changing" begin + @testset "3-bus case with optimal phase shifting with equal lb/ub" begin file = "../test/data/matpower/case3_oltc_pst.m" data = PowerModels.parse_file(file) - result = PowerModels._solve_opf_oltc_pst(data, ACPPowerModel, nlp_solver) + for (i, branch) in data["branch"] + branch["ta_min"] = branch["shift"] + branch["ta_max"] = branch["shift"] + end + result = PowerModels._solve_opf_pst(data, ACPPowerModel, nlp_solver) @test result["termination_status"] == LOCALLY_SOLVED - @test isapprox(result["objective"], 5738.6; atol = 1e0) + @test isapprox(result["objective"], 5820.1; atol = 1e0) + + @test isapprox(result["solution"]["branch"]["1"]["ta"], 0.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["2"]["ta"], 0.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["3"]["ta"], 5.0/180*pi; atol = 1e-1) + end + end + + + @testset "test dc polar opf" begin + @testset "3-bus case with optimal phase shifting" begin + file = "../test/data/matpower/case3_oltc_pst.m" + data = PowerModels.parse_file(file) + result = PowerModels._solve_opf_pst(data, DCPPowerModel, milp_solver) + + @test result["termination_status"] == OPTIMAL + @test isapprox(result["objective"], 5639.0; atol = 1e0) - @test haskey(result["solution"]["branch"]["1"], "tm") @test haskey(result["solution"]["branch"]["1"], "ta") - @test isapprox(result["solution"]["branch"]["1"]["tm"], 0.948; atol = 1e-2) @test isapprox(result["solution"]["branch"]["1"]["ta"], 0.000; atol = 1e-3) - @test isapprox(result["solution"]["branch"]["2"]["tm"], 1.100; atol = 1e-3) @test isapprox(result["solution"]["branch"]["2"]["ta"], 0.000; atol = 1e-3) - @test isapprox(result["solution"]["branch"]["3"]["tm"], 1.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["3"]["ta"], 15.0/180*pi; atol = 1e-1) end - - @testset "3-bus case with optimal phase shifting / tap changing with equal lb/ub" begin + @testset "3-bus case with optimal phase shifting with equal lb/ub" begin file = "../test/data/matpower/case3_oltc_pst.m" data = PowerModels.parse_file(file) for (i, branch) in data["branch"] branch["ta_min"] = branch["shift"] branch["ta_max"] = branch["shift"] - branch["tm_min"] = branch["tap"] - branch["tm_max"] = branch["tap"] end - result = PowerModels._solve_opf_oltc_pst(data, ACPPowerModel, nlp_solver) + result = PowerModels._solve_opf_pst(data, DCPPowerModel, milp_solver) - @test result["termination_status"] == LOCALLY_SOLVED - @test isapprox(result["objective"], 5820.1; atol = 1e0) + @test result["termination_status"] == OPTIMAL + @test isapprox(result["objective"], 5698.1; atol = 1e0) - @test isapprox(result["solution"]["branch"]["1"]["tm"], 1.00; atol = 1e-2) @test isapprox(result["solution"]["branch"]["1"]["ta"], 0.000; atol = 1e-3) - @test isapprox(result["solution"]["branch"]["2"]["tm"], 1.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["2"]["ta"], 0.000; atol = 1e-3) - @test isapprox(result["solution"]["branch"]["3"]["tm"], 1.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["3"]["ta"], 5.0/180*pi; atol = 1e-1) end end +end - @testset "test dc polar opf" begin + +@testset "test opf with optimization of oltc and pst" begin + + @testset "test ac polar opf" begin @testset "3-bus case with fixed phase shift / tap" begin file = "../test/data/matpower/case3_oltc_pst.m" data = PowerModels.parse_file(file) - result = PowerModels.run_opf(data, DCPPowerModel, milp_solver) + result = PowerModels.run_opf(data, ACPPowerModel, nlp_solver) - @test result["termination_status"] == OPTIMAL - @test isapprox(result["objective"], 5782.0; atol = 1e0) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["objective"], 5820.1; atol = 1e0) end - @testset "3-bus case with optimal phase shifting / tap changing" begin + @testset "3-bus case with optimal phase shifting and tap changing" begin file = "../test/data/matpower/case3_oltc_pst.m" data = PowerModels.parse_file(file) - result = PowerModels._solve_opf_oltc_pst(data, DCPPowerModel, milp_solver) + result = PowerModels._solve_opf_oltc_pst(data, ACPPowerModel, nlp_solver) - @test result["termination_status"] == OPTIMAL - @test isapprox(result["objective"], 5639.0; atol = 1e0) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["objective"], 5738.6; atol = 1e0) + @test haskey(result["solution"]["branch"]["1"], "tm") @test haskey(result["solution"]["branch"]["1"], "ta") - @test isapprox(result["solution"]["branch"]["1"]["tm"], 1.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["1"]["tm"], 0.948; atol = 1e-2) @test isapprox(result["solution"]["branch"]["1"]["ta"], 0.000; atol = 1e-3) - @test isapprox(result["solution"]["branch"]["2"]["tm"], 1.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["2"]["tm"], 1.100; atol = 1e-3) @test isapprox(result["solution"]["branch"]["2"]["ta"], 0.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["3"]["tm"], 1.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["3"]["ta"], 15.0/180*pi; atol = 1e-1) end - @testset "3-bus case with optimal phase shifting / tap changing with equal lb/ub" begin + @testset "3-bus case with optimal phase shifting and tap changing with equal lb/ub" begin file = "../test/data/matpower/case3_oltc_pst.m" data = PowerModels.parse_file(file) for (i, branch) in data["branch"] branch["ta_min"] = branch["shift"] branch["ta_max"] = branch["shift"] + branch["tm_min"] = branch["tap"] + branch["tm_max"] = branch["tap"] end - result = PowerModels._solve_opf_oltc_pst(data, DCPPowerModel, milp_solver) + result = PowerModels._solve_opf_oltc_pst(data, ACPPowerModel, nlp_solver) - @test result["termination_status"] == OPTIMAL - @test isapprox(result["objective"], 5698.1; atol = 1e0) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["objective"], 5820.1; atol = 1e0) - @test isapprox(result["solution"]["branch"]["1"]["tm"], 1.000; atol = 1e-3) + @test isapprox(result["solution"]["branch"]["1"]["tm"], 1.00; atol = 1e-2) @test isapprox(result["solution"]["branch"]["1"]["ta"], 0.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["2"]["tm"], 1.000; atol = 1e-3) @test isapprox(result["solution"]["branch"]["2"]["ta"], 0.000; atol = 1e-3)