add pst only opf problem spec

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

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)

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

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)

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

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)