KLU PF: use correct V0; return x in expected format; write results fu…

…nction supports KLU PF; implement tests for KLUACPowerFlow

Project.toml

@@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -22,6 +23,7 @@ DataStructures = "0.18"
 Dates = "1"
 InfrastructureSystems = "2"
 JSON3 = "1"
+KLU = "0.6.0"
 LinearAlgebra = "1"
 Logging = "1"
 NLsolve = "4"

src/PowerFlows.jl

@@ -5,6 +5,7 @@ export solve_ac_powerflow!
 export PowerFlowData
 export DCPowerFlow
 export NLSolveACPowerFlow
+export KLUACPowerFlow
 export PTDFDCPowerFlow
 export vPTDFDCPowerFlow
 export PSSEExportPowerFlow
@@ -20,6 +21,7 @@ import PowerSystems
 import PowerSystems: System
 import LinearAlgebra
 import NLsolve
+import KLU
 import SparseArrays
 import InfrastructureSystems
 import PowerNetworkMatrices

src/nlsolve_ac_powerflow.jl

@@ -165,25 +165,33 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k
     tol = 1e-6    # TODO
     i = 0
 
-    V = copy(pf.x0)
-    Vm = abs.(V)
-    Va = angle.(V)
+    Vm = data.bus_magnitude[:]
+    Va = data.bus_angles[:]
+    V = Vm .* exp.(1im * Va)
 
-    Ybus = pf.data.power_network_matrix
+    # Find indices for each bus type
+    ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type)
+    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type)
+    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type)
 
-    mis = V .* conj(Ybus * V) - Sbus  # TODO Sbus
+    Ybus = pf.data.power_network_matrix.data
+
+    Sbus = data.bus_activepower_injection[:] - data.bus_activepower_withdrawals[:] + 1im * (data.bus_reactivepower_injection[:] - data.bus_reactivepower_withdrawals[:])
+
+    mis = V .* conj(Ybus * V) - Sbus
     F = [real(mis[[pv; pq]]); imag(mis[pq])]
 
-    npv = length(pv) # TODO
-    npq = length(pq) # TODO
+    # nref = length(ref)
+    npv = length(pv)
+    npq = length(pq)
 
-    converged = false
+    converged = (npv + npq) == 0  # if only ref buses present, we do not need to enter the loop
 
     while i < maxIter && !converged
         i += 1
-        diagV = Diagonal(V)
-        diagIbus = Diagonal(Ybus * V)
-        diagVnorm = Diagonal(V ./ abs.(V))
+        diagV = LinearAlgebra.Diagonal(V)
+        diagIbus = LinearAlgebra.Diagonal(Ybus * V)
+        diagVnorm = LinearAlgebra.Diagonal(V ./ abs.(V))
         dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm
         dSbus_dVa = 1im * diagV * conj(diagIbus - Ybus * diagV)
 
@@ -193,7 +201,7 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k
         j22 = imag(dSbus_dVm[pq, pq])
         J = [j11 j12; j21 j22]
 
-        factor_J = klu(J)
+        factor_J = KLU.klu(J)
         dx = -(factor_J \ F)
 
         #J_function(J_function.Jv, x)
@@ -210,13 +218,35 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k
 
         mis = V .* conj(Ybus * V) - Sbus
         F = [real(mis[[pv; pq]]); imag(mis[pq])]
-        converged = norm(F, Inf) < tol
+        converged = LinearAlgebra.norm(F, Inf) < tol
     end
 
+
     if !converged
         @error("The powerflow solver with KLU did not converge after $i iterations")
     else
-        @info("The powerflow solver with KLU converged after $i iterations")
+        @debug("The powerflow solver with KLU converged after $i iterations")
+    end
+
+    # mock the expected x format, where the values depend on the type of the bus:
+    n_buses = length(data.bus_type)
+    x = Float64[0.0 for _ in 1:(2*n_buses)]
+    Sbus_result = V .* conj(Ybus * V)
+    for (ix, b) in enumerate(data.bus_type)
+        if b == PSY.ACBusTypes.REF
+            # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated
+            x[2 * ix - 1] = real(Sbus_result[ix]) + data.bus_activepower_withdrawals[ix]
+            x[2 * ix] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix]
+        elseif b == PSY.ACBusTypes.PV
+            # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle
+            x[2 * ix - 1] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix]
+            x[2 * ix] = Va[ix]
+        elseif b == PSY.ACBusTypes.PQ
+            # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle
+            x[2 * ix - 1] = Vm[ix]
+            x[2 * ix] = Va[ix]
+        end
     end
-    return converged, V  # TODO make sure the structure of V matches the structure of res.zero
+
+    return converged, x
 end

src/post_processing.jl

@@ -608,7 +608,7 @@ dictionary will therefore feature just one key linked to one DataFrame.
         vector containing the reults for one single time-period.
 """
 function write_results(
-    ::NLSolveACPowerFlow,
+    ::Union{NLSolveACPowerFlow,KLUACPowerFlow},
     sys::PSY.System,
     data::ACPowerFlowData,
     result::Vector{Float64},

test/test_nlsolve_powerflow.jl

@@ -1,5 +1,4 @@
-@testset "NLsolve Power Flow 14-Bus testing" begin
-    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
+@testset "AC Power Flow 14-Bus testing" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     result_14 = [
         2.3255081760423684
         -0.15529254415401786
@@ -30,60 +29,61 @@
         1.0213119628726421
         -0.2803812119374241
     ]
-    data = PowerFlows.PowerFlowData(NLSolveACPowerFlow(), sys; check_connectivity = true)
+
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
+    data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true)
     #Compare results between finite diff methods and Jacobian method
-    converged1, x1 = PowerFlows._solve_powerflow!(NLSolveACPowerFlow(), data, false)
-    @test LinearAlgebra.norm(result_14 - x1) <= 1e-6
-    @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys; method = :newton)
+    converged1, x1 = PowerFlows._solve_powerflow!(ACPowerFlow(), data, false)
+    @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6
+    @test solve_ac_powerflow!(ACPowerFlow(), sys; method = :newton)
 
     # Test enforcing the reactive power Limits
     set_reactive_power!(get_component(PowerLoad, sys, "Bus4"), 0.0)
-    data = PowerFlows.PowerFlowData(NLSolveACPowerFlow(), sys; check_connectivity = true)
-    converged2, x2 = PowerFlows._solve_powerflow!(NLSolveACPowerFlow(), data, true)
-    @test LinearAlgebra.norm(result_14 - x2) >= 1e-6
+    data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true)
+    converged2, x2 = PowerFlows._solve_powerflow!(ACPowerFlow(), data, true)
+    @test LinearAlgebra.norm(result_14 - x2, Inf) >= 1e-6
     @test 1.08 <= x2[15] <= 1.09
 end
 
-@testset "NLsolve Power Flow 14-Bus Line Configurations" begin
+@testset "AC Power Flow 14-Bus Line Configurations" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
-    base_res = solve_powerflow(NLSolveACPowerFlow(), sys)
+    base_res = solve_powerflow(ACPowerFlow(), sys)
     branch = first(PSY.get_components(Line, sys))
     dyn_branch = DynamicBranch(branch)
     add_component!(sys, dyn_branch)
-    @test dyn_pf = solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
-    dyn_pf = solve_powerflow(NLSolveACPowerFlow(), sys)
-    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - base_res["bus_results"].Vm) <=
-          1e-6
+    @test dyn_pf = solve_ac_powerflow!(ACPowerFlow(), sys)
+    dyn_pf = solve_powerflow(ACPowerFlow(), sys)
+    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - base_res["bus_results"].Vm, Inf) <= 1e-6
 
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
     line = get_component(Line, sys, "Line4")
     PSY.set_available!(line, false)
-    solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
+    solve_ac_powerflow!(ACPowerFlow(), sys)
     @test PSY.get_active_power_flow(line) == 0.0
     test_bus = get_component(PSY.Bus, sys, "Bus 4")
-    @test isapprox(PSY.get_magnitude(test_bus), 1.002; atol = 1e-3)
+    @test isapprox(PSY.get_magnitude(test_bus), 1.002; atol = 1e-3, rtol=0)
 
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
     line = get_component(Line, sys, "Line4")
     PSY.set_available!(line, false)
-    res = solve_powerflow(NLSolveACPowerFlow(), sys)
+    res = solve_powerflow(ACPowerFlow(), sys)
     @test res["flow_results"].P_from_to[4] == 0.0
     @test res["flow_results"].P_to_from[4] == 0.0
 end
 
-@testset "NLsolve Power Flow 3-Bus Fixed FixedAdmittance testing" begin
+@testset "AC Power Flow 3-Bus Fixed FixedAdmittance testing" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     p_gen_matpower_3bus = [20.3512373930753, 100.0, 100.0]
     q_gen_matpower_3bus = [45.516916781567232, 10.453799727283879, -31.992561631394636]
     sys_3bus = PSB.build_system(PSB.PSYTestSystems, "psse_3bus_gen_cls_sys")
     bus_103 = get_component(PSY.Bus, sys_3bus, "BUS 3")
     fix_shunt = PSY.FixedAdmittance("FixAdmBus3", true, bus_103, 0.0 + 0.2im)
     add_component!(sys_3bus, fix_shunt)
-    df = solve_powerflow(NLSolveACPowerFlow(), sys_3bus)
+    df = solve_powerflow(ACPowerFlow(), sys_3bus)
     @test isapprox(df["bus_results"].P_gen, p_gen_matpower_3bus, atol = 1e-4)
     @test isapprox(df["bus_results"].Q_gen, q_gen_matpower_3bus, atol = 1e-4)
 end
 
-@testset "NLsolve Power Flow convergence fail testing" begin
+@testset "AC Power Flow convergence fail testing" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     pf_sys5_re = PSB.build_system(PSB.PSITestSystems, "c_sys5_re"; add_forecasts = false)
     remove_component!(Line, pf_sys5_re, "1")
     remove_component!(Line, pf_sys5_re, "2")
@@ -95,11 +95,11 @@ end
     @test_logs(
         (:error, "The powerflow solver returned convergence = false"),
         match_mode = :any,
-        @test !solve_ac_powerflow!(NLSolveACPowerFlow(), pf_sys5_re)
+        @test !solve_ac_powerflow!(ACPowerFlow(), pf_sys5_re)
     )
 end
 
-@testset "Test 240 Case PSS/e results" begin
+@testset "AC Test 240 Case PSS/e results" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     file = joinpath(
         TEST_FILES_DIR,
         "test_data",
@@ -114,9 +114,9 @@ end
     pf_bus_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_bus_results.csv")
     pf_gen_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_gen_results.csv")
 
-    pf = solve_ac_powerflow!(NLSolveACPowerFlow(), system)
+    pf = solve_ac_powerflow!(ACPowerFlow(), system)
     @test pf
-    pf_result_df = solve_powerflow(NLSolveACPowerFlow(), system)
+    pf_result_df = solve_powerflow(ACPowerFlow(), system)
 
     v_diff, angle_diff, number = psse_bus_results_compare(pf_bus_result_file, pf_result_df)
     p_diff, q_diff, names = psse_gen_results_compare(pf_gen_result_file, system)
@@ -132,7 +132,7 @@ end
     @test norm(q_diff, 2) / length(q_diff) < DIFF_L2_TOLERANCE
 end
 
-@testset "Multiple sources at ref" begin
+@testset "AC Multiple sources at ref" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     sys = System(100.0)
     b = ACBus(;
         number = 1,
@@ -166,16 +166,14 @@ end
         X_th = 1e-5,
     )
     add_component!(sys, s2)
-    @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
+    @test solve_ac_powerflow!(ACPowerFlow(), sys)
 
     #Create power mismatch, test for error
     set_active_power!(get_component(Source, sys, "source_1"), -0.4)
-    @test_throws ErrorException(
-        "Sources do not match P and/or Q requirements for reference bus.",
-    ) solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
+    @test_throws ErrorException("Sources do not match P and/or Q requirements for reference bus.") solve_ac_powerflow!(ACPowerFlow(), sys)
 end
 
-@testset "AC PowerFlow with Multiple sources at PV" begin
+@testset "AC PowerFlow with Multiple sources at PV" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     sys = System(100.0)
     b1 = ACBus(;
         number = 1,
@@ -245,16 +243,14 @@ end
     )
     add_component!(sys, s3)
 
-    @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
+    @test solve_ac_powerflow!(ACPowerFlow(), sys)
 
     #Create power mismatch, test for error
     set_reactive_power!(get_component(Source, sys, "source_3"), -0.5)
-    @test_throws ErrorException("Sources do not match Q requirements for PV bus.") solve_ac_powerflow!(NLSolveACPowerFlow(), 
-        sys,
-    )
+    @test_throws ErrorException("Sources do not match Q requirements for PV bus.") solve_ac_powerflow!(ACPowerFlow(), sys)
 end
 
-@testset "AC PowerFlow Source + non-source at Ref" begin
+@testset "AC PowerFlow Source + non-source at Ref" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     sys = System(100.0)
     b = ACBus(;
         number = 1,
@@ -300,7 +296,7 @@ end
     )
     add_component!(sys, g1)
 
-    @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
+    @test solve_ac_powerflow!(ACPowerFlow(), sys)
     @test isapprox(
         get_active_power(get_component(Source, sys, "source_1")),
         0.5;
@@ -313,7 +309,7 @@ end
     )
 end
 
-@testset "AC PowerFlow Source + non-source at PV" begin
+@testset "AC PowerFlow Source + non-source at PV" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
     sys = System(100.0)
     b1 = ACBus(;
         number = 1,
@@ -394,7 +390,7 @@ end
     )
     add_component!(sys, g1)
 
-    @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys)
+    @test solve_ac_powerflow!(ACPowerFlow(), sys)
     @test isapprox(
         get_active_power(get_component(Source, sys, "source_2")),
         0.5;