fix lpac for PSTs

fix lpac for PSTs
add unit tests for lpac
tag new version

Project.toml

@@ -1,7 +1,7 @@
 name = "CbaOPF"
 uuid = "b09c5755-86bc-4b69-9fb0-eea9d409a5ca"
 authors = ["Hakan Ergun <hakan.ergun@ieee.org>"]
-version = "0.1.2"
+version = "0.1.3"
 
 
 

src/core/base.jl

@@ -16,6 +16,10 @@ function ref_add_pst!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
             push!(bus_arcs_pst[i], (l,i,j))
         end
         nw_ref[:bus_arcs_pst] = bus_arcs_pst
+
+        if !haskey(nw_ref, :buspairs_pst)
+            nw_ref[:buspairs_pst] = _PM.calc_buspair_parameters(nw_ref[:bus], nw_ref[:pst], nw_ref[:conductor_ids], haskey(nw_ref, :conductors))
+        end
     end
 end
 

src/core/data.jl

@@ -22,6 +22,7 @@ function fix_data_single_network_pst!(data)
         pst["b_fr"] = 0
         pst["g_to"] = 0
         pst["b_to"] = 0
+        pst["tap"] = 1.0
     end
 end
 function to_pu_multi_network_pst!(data)
@@ -41,6 +42,7 @@ function fix_data_multi_network_pst!(data)
             pst["b_fr"] = 0
             pst["g_to"] = 0
             pst["b_to"] = 0
+            pst["tap"] = 1.0
         end
     end
 end

src/core/variable.jl

@@ -3,6 +3,7 @@ function variable_pst(pm; kwargs...)
     variable_active_pst_flow(pm, kwargs...)
     variable_reactive_pst_flow(pm, kwargs...)
     variable_pst_angle(pm, kwargs...)
+    variable_pst_cosine(pm, kwargs...)
 end
 
 "variable: `p[l,i,j]` for `(l,i,j)` in `arcs`"
@@ -56,6 +57,44 @@ function variable_pst_angle(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_defaul
     report && _PM.sol_component_value(pm, nw, :pst, :alpha, _PM.ids(pm, nw, :pst), alpha)
 end
 
+function variable_pst_cosine(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
+end
+
+function variable_pst_cosine(pm::_PM.LPACCPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
+    cs_pst = _PM.var(pm, nw)[:cs_pst] = JuMP.@variable(pm.model,
+        [bp in _PM.ids(pm, nw, :buspairs_pst)], base_name="$(nw)_cs",
+        start = _PM.comp_start_value(_PM.ref(pm, nw, :buspairs_pst, bp), "cs_start", 1.0)
+    )
+
+    if bounded
+        for (bp, buspair) in _PM.ref(pm, nw, :buspairs_pst)
+            angmin = buspair["angmin"]
+            angmax = buspair["angmax"]
+            if angmin >= 0
+                cos_max = cos(angmin)
+                cos_min = cos(angmax)
+            end
+            if angmax <= 0
+                cos_max = cos(angmax)
+                cos_min = cos(angmin)
+            end
+            if angmin < 0 && angmax > 0
+                cos_max = 1.0
+                cos_min = min(cos(angmin), cos(angmax))
+            end
+
+            JuMP.set_lower_bound(cs_pst[bp], cos_min)
+            JuMP.set_upper_bound(cs_pst[bp], cos_max)
+        end
+    end
+
+    report && _PM.sol_component_value_buspair(pm, nw, :buspairs_pst, :cs_pst, _PM.ids(pm, nw, :buspairs_pst), cs_pst)
+end
+
+
+
+
+
 function variable_flexible_demand(pm::_PM.AbstractPowerModel; kwargs...)
     variable_total_flex_demand(pm; kwargs...)
     variable_demand_reduction(pm; kwargs...)

src/form/lpac.jl

@@ -54,25 +54,25 @@ function constraint_ohms_y_from_pst(pm::_PM.AbstractLPACModel, n::Int, i::Int, f
     phi_to = _PM.var(pm, n, :phi, t_bus)
     va_fr = _PM.var(pm, n, :va, f_bus)
     va_to = _PM.var(pm, n, :va, t_bus)
-    cs = _PM.var(pm, n, :cs, f_bus)
+    cs = _PM.var(pm, n, :cs_pst, (f_bus, t_bus))
 
-    JuMP.@constraint(model, p_fr ==  g * (1.0 + 2*phi_fr) - g * (cs + phi_fr + phi_to) - b * (va_fr - va_to - alpha))
-    JuMP.@constraint(model, q_fr == -b * (1.0 + 2*phi_fr) + b * (cs + phi_fr + phi_to) - g * (va_fr - va_to - alpha))
+    JuMP.@constraint(pm.model, p_fr ==  g * (1.0 + 2*phi_fr) - g * (cs + phi_fr + phi_to) - b * (va_fr - va_to - alpha))
+    JuMP.@constraint(pm.model, q_fr == -b * (1.0 + 2*phi_fr) + b * (cs + phi_fr + phi_to) - g * (va_fr - va_to - alpha))
 end
 
 function constraint_ohms_y_to_pst(pm::_PM.AbstractLPACModel, n::Int, i::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr)
     alpha = _PM.var(pm, n,  :psta, i)
-    p_fr  = _PM.var(pm, n,  :ppst, f_idx)
-    q_fr  = _PM.var(pm, n,  :qpst, f_idx)
+    p_to  = _PM.var(pm, n,  :ppst, t_idx)
+    q_to  = _PM.var(pm, n,  :qpst, t_idx)
 
     phi_fr = _PM.var(pm, n, :phi, f_bus)
     phi_to = _PM.var(pm, n, :phi, t_bus)
     va_fr = _PM.var(pm, n, :va, f_bus)
     va_to = _PM.var(pm, n, :va, t_bus)
-    cs = _PM.var(pm, n, :cs, t_bus)
+    cs = _PM.var(pm, n, :cs_pst, (f_bus, t_bus))
 
-    JuMP.@constraint(model, p_to ==  g * (1.0 + 2 * phi_to) - g * (cs + phi_fr + phi_to) -b * (va_to - va_fr + alpha))
-    JuMP.@constraint(model, q_to == -b * (1.0 + 2 * phi_to) + b * (cs + phi_fr + phi_to) -g * (va_to - va_fr + alpha))
+    JuMP.@constraint(pm.model, p_to ==  g * (1.0 + 2 * phi_to) - g * (cs + phi_fr + phi_to) -b * (va_to - va_fr + alpha))
+    JuMP.@constraint(pm.model, q_to == -b * (1.0 + 2 * phi_to) + b * (cs + phi_fr + phi_to) -g * (va_to - va_fr + alpha))
 end
 
 

test/cbaopf.jl

@@ -88,12 +88,12 @@ resultOPF = CbaOPF.solve_cbaopf(data, _PM.SOCWRPowerModel, ipopt; setting = s)
     @test isapprox(resultOPF["solution"]["convdc"]["4"]["ptf_to"], -0.005477, atol = 1e-2)
 end
 
-# resultOPF = CbaOPF.solve_cbaopf(data_ac, _PM.LPACCPowerModel, highs; setting = s)
+resultOPF = CbaOPF.solve_cbaopf(data, _PM.LPACCPowerModel, ipopt; setting = s)
 
-# @testset  "Nodal LPACC - CBA OPF" begin
+@testset  "Nodal LPACC - CBA OPF" begin
 
-#     @test isapprox(resultOPF["objective"], 3113.2, atol = 1e-1)
-#     @test isapprox(resultOPF["solution"]["gen"]["3"]["pg"], 0, atol = 1e-2)
-#     @test isapprox(resultOPF["solution"]["branch"]["2"]["pt"], -0.0074, atol = 1e-2)
-#     @test isapprox(resultOPF["solution"]["convdc"]["4"]["ptf_to"], -0.005477, atol = 1e-2)
-# end
+    @test isapprox(resultOPF["objective"],  24519.2, atol = 1e-1)
+    @test isapprox(resultOPF["solution"]["gen"]["3"]["pg"], 1.15441, atol = 1e-2)
+    @test isapprox(resultOPF["solution"]["branch"]["2"]["pt"], -1.75847, atol = 1e-2)
+    @test isapprox(resultOPF["solution"]["convdc"]["4"]["ptf_to"], 0.831186, atol = 1e-2)
+end