From 251f1b306d89e5deae5883fcf49d883705866728 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Fri, 6 Dec 2024 14:47:49 -0700 Subject: [PATCH] KLU PF: use correct V0; return x in expected format; write results function supports KLU PF; implement tests for KLUACPowerFlow --- Project.toml | 2 + src/PowerFlows.jl | 2 + src/nlsolve_ac_powerflow.jl | 60 ++++++++++++++++++++------- src/post_processing.jl | 2 +- test/test_nlsolve_powerflow.jl | 74 ++++++++++++++++------------------ 5 files changed, 85 insertions(+), 55 deletions(-) diff --git a/Project.toml b/Project.toml index b61464a0..b339e3f3 100644 --- a/Project.toml +++ b/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" diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index feeb0f1d..499ac42c 100644 --- a/src/PowerFlows.jl +++ b/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 diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 4c32fa4e..c05c7ecd 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/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 diff --git a/src/post_processing.jl b/src/post_processing.jl index efa5fe6c..9b38e858 100644 --- a/src/post_processing.jl +++ b/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}, diff --git a/test/test_nlsolve_powerflow.jl b/test/test_nlsolve_powerflow.jl index 189cc5b3..a11e7124 100644 --- a/test/test_nlsolve_powerflow.jl +++ b/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;