From 198a4fe100683611838b2ece7cfab76f1cc0e8c9 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 5 Dec 2024 08:55:27 -0700 Subject: [PATCH 01/12] WIP: implement NR power flow function, rename ACPowerFlow to NLSolveACPowerFlow, introduce KLUACPowerFlow for the use in the new NR power flow function --- src/PowerFlowData.jl | 4 +- src/PowerFlows.jl | 2 +- src/nlsolve_ac_powerflow.jl | 99 ++++++++++++++++++++++++++-------- src/post_processing.jl | 2 +- src/powerflow_types.jl | 6 ++- test/test_nlsolve_powerflow.jl | 46 ++++++++-------- test/test_powerflow_data.jl | 12 ++--- test/test_psse_export.jl | 6 +-- 8 files changed, 118 insertions(+), 59 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index f064a02c..00c8c4a2 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -157,7 +157,7 @@ NOTE: use it for AC power flow computations. WARNING: functions for the evaluation of the multi-period AC PF still to be implemented. """ function PowerFlowData( - ::ACPowerFlow, + ::Union{NLSolveACPowerFlow, KLUACPowerFlow}, sys::PSY.System; time_steps::Int = 1, timestep_names::Vector{String} = String[], @@ -440,7 +440,7 @@ Create an appropriate `PowerFlowContainer` for the given `PowerFlowEvaluationMod """ function make_power_flow_container end -make_power_flow_container(pfem::ACPowerFlow, sys::PSY.System; kwargs...) = +make_power_flow_container(pfem::NLSolveACPowerFlow, sys::PSY.System; kwargs...) = PowerFlowData(pfem, sys; kwargs...) make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) = diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index 7d88f80e..feeb0f1d 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -4,7 +4,7 @@ export solve_powerflow export solve_ac_powerflow! export PowerFlowData export DCPowerFlow -export ACPowerFlow +export NLSolveACPowerFlow export PTDFDCPowerFlow export vPTDFDCPowerFlow export PSSEExportPowerFlow diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 813848b0..576b0cec 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -30,29 +30,29 @@ solve_ac_powerflow!(sys) solve_ac_powerflow!(sys, method=:newton) ``` """ -function solve_ac_powerflow!(system::PSY.System; kwargs...) +function solve_ac_powerflow!(pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, system::PSY.System; kwargs...) #Save per-unit flag settings_unit_cache = deepcopy(system.units_settings.unit_system) #Work in System per unit PSY.set_units_base_system!(system, "SYSTEM_BASE") check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false) data = PowerFlowData( - ACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits), + NLSolveACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits), system; check_connectivity = get(kwargs, :check_connectivity, true), ) max_iterations = DEFAULT_MAX_REDISTRIBUTION_ITERATIONS - res = _solve_powerflow!(data, check_reactive_power_limits; kwargs...) - if res.f_converged - write_powerflow_solution!(system, res.zero, max_iterations) + converged, x = _solve_powerflow!(pf, data, check_reactive_power_limits; kwargs...) + if converged + write_powerflow_solution!(system, x, max_iterations) @info("PowerFlow solve converged, the results have been stored in the system") #Restore original per unit base PSY.set_units_base_system!(system, settings_unit_cache) - return res.f_converged + return converged end - @error("The powerflow solver returned convergence = $(res.f_converged)") + @error("The powerflow solver returned convergence = $(converged)") PSY.set_units_base_system!(system, settings_unit_cache) - return res.f_converged + return converged end """ @@ -68,7 +68,7 @@ res = solve_powerflow(sys, method=:newton) ``` """ function solve_powerflow( - pf::ACPowerFlow, + pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, system::PSY.System; kwargs..., ) @@ -82,18 +82,18 @@ function solve_powerflow( check_connectivity = get(kwargs, :check_connectivity, true), ) - res = _solve_powerflow!(data, pf.check_reactive_power_limits; kwargs...) + converged, x = _solve_powerflow!(pf, data, pf.check_reactive_power_limits; kwargs...) - if res.f_converged + if converged @info("PowerFlow solve converged, the results are exported in DataFrames") - df_results = write_results(pf, system, data, res.zero) + df_results = write_results(pf, system, data, x) #Restore original per unit base PSY.set_units_base_system!(system, settings_unit_cache) return df_results end - @error("The powerflow solver returned convergence = $(res.f_converged)") + @error("The powerflow solver returned convergence = $(converged)") PSY.set_units_base_system!(system, settings_unit_cache) - return res.f_converged + return converged end function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) @@ -124,27 +124,28 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) end function _solve_powerflow!( + pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, data::ACPowerFlowData, check_reactive_power_limits; nlsolve_kwargs..., ) if check_reactive_power_limits for _ in 1:MAX_REACTIVE_POWER_ITERATIONS - res = _nlsolve_powerflow(data; nlsolve_kwargs...) - if res.f_converged - if _check_q_limit_bounds!(data, res.zero) - return res + converged, x = _nlsolve_powerflow(pf, data; nlsolve_kwargs...) + if converged + if _check_q_limit_bounds!(data, x) + return converged, x end else - return res + return converged, x end end else - return _nlsolve_powerflow(data; nlsolve_kwargs...) + return _nlsolve_powerflow(pf, data; nlsolve_kwargs...) end end -function _nlsolve_powerflow(data::ACPowerFlowData; nlsolve_kwargs...) +function _nlsolve_powerflow(pf::NLSolveACPowerFlow, data::ACPowerFlowData; nlsolve_kwargs...) pf = PolarPowerFlow(data) J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) @@ -153,5 +154,59 @@ function _nlsolve_powerflow(data::ACPowerFlowData; nlsolve_kwargs...) if !res.f_converged @error("The powerflow solver returned convergence = $(res.f_converged)") end - return res + return res.f_converged, res.zero +end + +function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_kwargs...) + pf = PolarPowerFlow(data) + J_function = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) + + maxIter = 30 + tol = 1e-6 + + V = pf.x0 + Vm = abs.(V) + Va = angle.(V) + + mis = V .* conj(Ybus * V) - Sbus + F = [real(mis[[pv;pq]]); imag(mis[pq])] + + npv = length(pv) + npq = length(pq) + + converged = false + + while i < maxIter && converged + i += 1 + #diagV = Diagonal(V) + #diagIbus = Diagonal(Ybus * V) + #diagVnorm = Diagonal(V./abs.(V)) + #dSbus_dVm = Diagonal(V) * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm + #dSbus_dVa = im * diagV * conj(diagIbus - Ybus * diagV) + #j11 = real(dSbus_dVa[[pv; pq], [pv; pq]]) + #j12 = real(dSbus_dVm[[pv; pq], pq]) + #j21 = imag(dSbus_dVa[pq, [pv; pq]]) + #j22 = imag(dSbus_dVm[pq, pq]) + #J = [j11 j12; j21 j22] + J_function(J_function.Jv, x) + dx = - (J_function.Jv \ F) + + Va[pv] += dx[1:npv] + Va[pq] += dx[npv+1:(npv+npq)] + Vm[pq] += dx[(npv+npq+1):end] + + V = Vm .* exp.(im*Va) + + Vm = abs.(V) + Va = angle.(V) + + mis = V .* conj(Ybus * V) - Sbus + F = [real(mis[[pv;pq]]); imag(mis[pq])] + converged = norm(F, Inf) < tol + end + + if !converged + @error("The powerflow solver with KLU did not converge after $maxIter iterations") + end + return converged, V end diff --git a/src/post_processing.jl b/src/post_processing.jl index 6c500a25..efa5fe6c 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( - ::ACPowerFlow, + ::NLSolveACPowerFlow, sys::PSY.System, data::ACPowerFlowData, result::Vector{Float64}, diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl index e1b9828a..65b6b241 100644 --- a/src/powerflow_types.jl +++ b/src/powerflow_types.jl @@ -1,6 +1,10 @@ abstract type PowerFlowEvaluationModel end -Base.@kwdef struct ACPowerFlow <: PowerFlowEvaluationModel +Base.@kwdef struct NLSolveACPowerFlow <: PowerFlowEvaluationModel + check_reactive_power_limits::Bool = false +end + +Base.@kwdef struct KLUACPowerFlow <: PowerFlowEvaluationModel check_reactive_power_limits::Bool = false end diff --git a/test/test_nlsolve_powerflow.jl b/test/test_nlsolve_powerflow.jl index eb2a2813..189cc5b3 100644 --- a/test/test_nlsolve_powerflow.jl +++ b/test/test_nlsolve_powerflow.jl @@ -30,35 +30,35 @@ 1.0213119628726421 -0.2803812119374241 ] - data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true) + data = PowerFlows.PowerFlowData(NLSolveACPowerFlow(), sys; check_connectivity = true) #Compare results between finite diff methods and Jacobian method - res1 = PowerFlows._solve_powerflow!(data, false) - @test LinearAlgebra.norm(result_14 - res1.zero) <= 1e-6 - @test solve_ac_powerflow!(sys; method = :newton) + converged1, x1 = PowerFlows._solve_powerflow!(NLSolveACPowerFlow(), data, false) + @test LinearAlgebra.norm(result_14 - x1) <= 1e-6 + @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys; method = :newton) # Test enforcing the reactive power Limits set_reactive_power!(get_component(PowerLoad, sys, "Bus4"), 0.0) - data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true) - res2 = PowerFlows._solve_powerflow!(data, true) - @test LinearAlgebra.norm(result_14 - res2.zero) >= 1e-6 - @test 1.08 <= res2.zero[15] <= 1.09 + data = PowerFlows.PowerFlowData(NLSolveACPowerFlow(), sys; check_connectivity = true) + converged2, x2 = PowerFlows._solve_powerflow!(NLSolveACPowerFlow(), data, true) + @test LinearAlgebra.norm(result_14 - x2) >= 1e-6 + @test 1.08 <= x2[15] <= 1.09 end @testset "NLsolve Power Flow 14-Bus Line Configurations" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) - base_res = solve_powerflow(ACPowerFlow(), sys) + base_res = solve_powerflow(NLSolveACPowerFlow(), sys) branch = first(PSY.get_components(Line, sys)) dyn_branch = DynamicBranch(branch) add_component!(sys, dyn_branch) - @test dyn_pf = solve_ac_powerflow!(sys) - dyn_pf = solve_powerflow(ACPowerFlow(), sys) + @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 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!(sys) + solve_ac_powerflow!(NLSolveACPowerFlow(), 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) @@ -66,7 +66,7 @@ end 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(ACPowerFlow(), sys) + res = solve_powerflow(NLSolveACPowerFlow(), sys) @test res["flow_results"].P_from_to[4] == 0.0 @test res["flow_results"].P_to_from[4] == 0.0 end @@ -78,7 +78,7 @@ end 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(ACPowerFlow(), sys_3bus) + df = solve_powerflow(NLSolveACPowerFlow(), 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 @@ -95,7 +95,7 @@ end @test_logs( (:error, "The powerflow solver returned convergence = false"), match_mode = :any, - @test !solve_ac_powerflow!(pf_sys5_re) + @test !solve_ac_powerflow!(NLSolveACPowerFlow(), pf_sys5_re) ) end @@ -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!(system) + pf = solve_ac_powerflow!(NLSolveACPowerFlow(), system) @test pf - pf_result_df = solve_powerflow(ACPowerFlow(), system) + pf_result_df = solve_powerflow(NLSolveACPowerFlow(), 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) @@ -166,13 +166,13 @@ end X_th = 1e-5, ) add_component!(sys, s2) - @test solve_ac_powerflow!(sys) + @test solve_ac_powerflow!(NLSolveACPowerFlow(), 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!(sys) + ) solve_ac_powerflow!(NLSolveACPowerFlow(), sys) end @testset "AC PowerFlow with Multiple sources at PV" begin @@ -245,11 +245,11 @@ end ) add_component!(sys, s3) - @test solve_ac_powerflow!(sys) + @test solve_ac_powerflow!(NLSolveACPowerFlow(), 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!( + @test_throws ErrorException("Sources do not match Q requirements for PV bus.") solve_ac_powerflow!(NLSolveACPowerFlow(), sys, ) end @@ -300,7 +300,7 @@ end ) add_component!(sys, g1) - @test solve_ac_powerflow!(sys) + @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys) @test isapprox( get_active_power(get_component(Source, sys, "source_1")), 0.5; @@ -394,7 +394,7 @@ end ) add_component!(sys, g1) - @test solve_ac_powerflow!(sys) + @test solve_ac_powerflow!(NLSolveACPowerFlow(), sys) @test isapprox( get_active_power(get_component(Source, sys, "source_2")), 0.5; diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl index 75964cb7..5f91cddc 100644 --- a/test/test_powerflow_data.jl +++ b/test/test_powerflow_data.jl @@ -1,6 +1,6 @@ @testset "PowerFlowData" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) - @test PowerFlowData(ACPowerFlow(), sys) isa PF.ACPowerFlowData + @test PowerFlowData(NLSolveACPowerFlow(), sys) isa PF.ACPowerFlowData @test PowerFlowData(DCPowerFlow(), sys) isa PF.ABAPowerFlowData @test PowerFlowData(PTDFDCPowerFlow(), sys) isa PF.PTDFPowerFlowData @test PowerFlowData(vPTDFDCPowerFlow(), sys) isa PF.vPTDFPowerFlowData @@ -22,18 +22,18 @@ end # TODO test that update_system! errors if the PowerFlowData doesn't correspond to the system sys_original = build_system(PSISystems, "RTS_GMLC_DA_sys") - data_original = PowerFlowData(ACPowerFlow(), sys_original) + data_original = PowerFlowData(NLSolveACPowerFlow(), sys_original) sys_modified = deepcopy(sys_original) modify_rts_system!(sys_modified) - data_modified = PowerFlowData(ACPowerFlow(), sys_original) + data_modified = PowerFlowData(NLSolveACPowerFlow(), sys_original) modify_rts_powerflow!(data_modified) # update_system! with unmodified PowerFlowData should result in system that yields unmodified PowerFlowData # (NOTE does NOT necessarily yield original system due to power redistribution) sys_null_updated = deepcopy(sys_original) PF.update_system!(sys_null_updated, data_original) - data_null_updated = PowerFlowData(ACPowerFlow(), sys_null_updated) + data_null_updated = PowerFlowData(NLSolveACPowerFlow(), sys_null_updated) @test IS.compare_values(powerflow_match_fn, data_null_updated, data_original) # Modified versions should not be the same as unmodified versions @@ -47,7 +47,7 @@ end # Constructing PowerFlowData from modified system should result in data_modified @test IS.compare_values( powerflow_match_fn, - PowerFlowData(ACPowerFlow(), sys_modified), + PowerFlowData(NLSolveACPowerFlow(), sys_modified), data_modified, ) @@ -56,6 +56,6 @@ end sys_modify_updated = deepcopy(sys_original) PF.update_system!(sys_modify_updated, data_modified) sys_mod_redist = deepcopy(sys_modified) - PF.update_system!(sys_mod_redist, PowerFlowData(ACPowerFlow(), sys_mod_redist)) + PF.update_system!(sys_mod_redist, PowerFlowData(NLSolveACPowerFlow(), sys_mod_redist)) @test IS.compare_values(powerflow_match_fn, sys_modify_updated, sys_mod_redist) end diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl index d6abd32d..df7d925d 100644 --- a/test/test_psse_export.jl +++ b/test/test_psse_export.jl @@ -223,8 +223,8 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; end function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow = false) - result1 = solve_powerflow(ACPowerFlow(), sys1) - result2 = solve_powerflow(ACPowerFlow(), sys2) + result1 = solve_powerflow(NLSolveACPowerFlow(), sys1) + result2 = solve_powerflow(NLSolveACPowerFlow(), sys2) reactive_power_tol = exclude_reactive_flow ? nothing : POWERFLOW_COMPARISON_TOLERANCE @test compare_df_within_tolerance("bus_results", result1["bus_results"], @@ -408,7 +408,7 @@ end # Updating with changed value should result in a different reimport (PowerFlowData version) exporter = PSSEExporter(sys, :v33, export_location) - pf2 = PowerFlowData(ACPowerFlow(), sys) + pf2 = PowerFlowData(NLSolveACPowerFlow(), sys) # This modifies the PowerFlowData in the same way that modify_rts_system! modifies the # system, so the reimport should be comparable to sys2 from above modify_rts_powerflow!(pf2) From d4a443acf0e96434c71290336f278a9221307a33 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 5 Dec 2024 12:25:16 -0700 Subject: [PATCH 02/12] WIP: some bugfixes in the PF function --- src/nlsolve_ac_powerflow.jl | 72 +++++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 31 deletions(-) diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 576b0cec..8f2c4fb6 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -159,54 +159,64 @@ end function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_kwargs...) pf = PolarPowerFlow(data) - J_function = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) + #J_function = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) maxIter = 30 tol = 1e-6 + i = 0 - V = pf.x0 + V = copy(pf.x0) Vm = abs.(V) Va = angle.(V) - mis = V .* conj(Ybus * V) - Sbus - F = [real(mis[[pv;pq]]); imag(mis[pq])] + Ybus = pf.data.power_network_matrix + + mis = V .* conj(Ybus * V) - Sbus # TODO Sbus + F = [real(mis[[pv; pq]]); imag(mis[pq])] + + npv = length(pv) # TODO + npq = length(pq) # TODO - npv = length(pv) - npq = length(pq) - converged = false - while i < maxIter && converged + while i < maxIter && !converged i += 1 - #diagV = Diagonal(V) - #diagIbus = Diagonal(Ybus * V) - #diagVnorm = Diagonal(V./abs.(V)) - #dSbus_dVm = Diagonal(V) * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm - #dSbus_dVa = im * diagV * conj(diagIbus - Ybus * diagV) - #j11 = real(dSbus_dVa[[pv; pq], [pv; pq]]) - #j12 = real(dSbus_dVm[[pv; pq], pq]) - #j21 = imag(dSbus_dVa[pq, [pv; pq]]) - #j22 = imag(dSbus_dVm[pq, pq]) - #J = [j11 j12; j21 j22] - J_function(J_function.Jv, x) - dx = - (J_function.Jv \ F) - - Va[pv] += dx[1:npv] - Va[pq] += dx[npv+1:(npv+npq)] - Vm[pq] += dx[(npv+npq+1):end] - - V = Vm .* exp.(im*Va) - + diagV = Diagonal(V) + diagIbus = Diagonal(Ybus * V) + diagVnorm = Diagonal(V ./ abs.(V)) + dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm + dSbus_dVa = 1im * diagV * conj(diagIbus - Ybus * diagV) + + j11 = real(dSbus_dVa[[pv; pq], [pv; pq]]) + j12 = real(dSbus_dVm[[pv; pq], pq]) + j21 = imag(dSbus_dVa[pq, [pv; pq]]) + j22 = imag(dSbus_dVm[pq, pq]) + J = [j11 j12; j21 j22] + + factor_J = klu(J) + dx = -(factor_J \ F) + + #J_function(J_function.Jv, x) + #dx = - (J_function.Jv \ F) + + Va[pv] .+= dx[1:npv] + Va[pq] .+= dx[(npv+1):(npv+npq)] + Vm[pq] .+= dx[(npv+npq+1):(npv+2*npq)] + + V = Vm .* exp.(1im * Va) + Vm = abs.(V) Va = angle.(V) mis = V .* conj(Ybus * V) - Sbus - F = [real(mis[[pv;pq]]); imag(mis[pq])] + F = [real(mis[[pv; pq]]); imag(mis[pq])] converged = norm(F, Inf) < tol end - + if !converged - @error("The powerflow solver with KLU did not converge after $maxIter iterations") + @error("The powerflow solver with KLU did not converge after $i iterations") + else + @info("The powerflow solver with KLU converged after $i iterations") end - return converged, V + return converged, V # TODO make sure the structure of V matches the structure of res.zero end From f2ad409b4553192d459aa4db9babc87682c88289 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 5 Dec 2024 12:27:25 -0700 Subject: [PATCH 03/12] small changes --- src/nlsolve_ac_powerflow.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 8f2c4fb6..4c32fa4e 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -161,8 +161,8 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k pf = PolarPowerFlow(data) #J_function = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) - maxIter = 30 - tol = 1e-6 + maxIter = 30 # TODO + tol = 1e-6 # TODO i = 0 V = copy(pf.x0) From 251f1b306d89e5deae5883fcf49d883705866728 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Fri, 6 Dec 2024 14:47:49 -0700 Subject: [PATCH 04/12] 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; From 210937c1fdc2f58460211ce8188645b04f787e5e Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Fri, 6 Dec 2024 15:01:14 -0700 Subject: [PATCH 05/12] reformat code --- src/PowerFlowData.jl | 58 +++--- src/nlsolve_ac_powerflow.jl | 26 +-- src/post_processing.jl | 148 +++++++------- test/test_nlsolve_powerflow.jl | 346 ++++++++++++++++----------------- test/test_powerflow_data.jl | 10 +- test/test_psse_export.jl | 68 +++---- 6 files changed, 328 insertions(+), 328 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 00c8c4a2..ac0f198f 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -66,11 +66,11 @@ flows and angles, as well as these ones. - `neighbors::Vector{Set{Int}}`: Vector with the sets of adjacent buses. """ struct PowerFlowData{ - M <: PNM.PowerNetworkMatrix, - N <: Union{PNM.PowerNetworkMatrix, Nothing}, + M<:PNM.PowerNetworkMatrix, + N<:Union{PNM.PowerNetworkMatrix,Nothing}, } <: PowerFlowContainer - bus_lookup::Dict{Int, Int} - branch_lookup::Dict{String, Int} + bus_lookup::Dict{Int,Int} + branch_lookup::Dict{String,Int} bus_activepower_injection::Matrix{Float64} bus_reactivepower_injection::Matrix{Float64} bus_activepower_withdrawals::Matrix{Float64} @@ -81,7 +81,7 @@ struct PowerFlowData{ bus_magnitude::Matrix{Float64} bus_angles::Matrix{Float64} branch_flow_values::Matrix{Float64} - timestep_map::Dict{Int, String} + timestep_map::Dict{Int,String} valid_ix::Vector{Int} power_network_matrix::M aux_network_matrix::N @@ -119,8 +119,8 @@ end # TODO -> MULTI PERIOD: AC Power Flow Data function _calculate_neighbors( Yb::PNM.Ybus{ - Tuple{Vector{Int64}, Vector{Int64}}, - Tuple{Dict{Int64, Int64}, Dict{Int64, Int64}}, + Tuple{Vector{Int64},Vector{Int64}}, + Tuple{Dict{Int64,Int64},Dict{Int64,Int64}}, }, ) I, J, V = SparseArrays.findnz(Yb.data) @@ -157,11 +157,11 @@ NOTE: use it for AC power flow computations. WARNING: functions for the evaluation of the multi-period AC PF still to be implemented. """ function PowerFlowData( - ::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + ::Union{NLSolveACPowerFlow,KLUACPowerFlow}, sys::PSY.System; - time_steps::Int = 1, - timestep_names::Vector{String} = String[], - check_connectivity::Bool = true) + time_steps::Int=1, + timestep_names::Vector{String}=String[], + check_connectivity::Bool=true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -174,7 +174,7 @@ function PowerFlowData( end # get data for calculations - power_network_matrix = PNM.Ybus(sys; check_connectivity = check_connectivity) + power_network_matrix = PNM.Ybus(sys; check_connectivity=check_connectivity) # get number of buses and branches n_buses = length(axes(power_network_matrix, 1)) @@ -185,8 +185,8 @@ function PowerFlowData( n_branches = length(branches) bus_lookup = power_network_matrix.lookup[2] - branch_lookup = Dict{String, Int}() - temp_bus_map = Dict{Int, String}( + branch_lookup = Dict{String,Int}() + temp_bus_map = Dict{Int,String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) branch_types = Vector{DataType}(undef, n_branches) @@ -250,9 +250,9 @@ NOTE: use it for DC power flow computations. function PowerFlowData( ::DCPowerFlow, sys::PSY.System; - time_steps::Int = 1, - timestep_names::Vector{String} = String[], - check_connectivity::Bool = true) + time_steps::Int=1, + timestep_names::Vector{String}=String[], + check_connectivity::Bool=true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -263,7 +263,7 @@ function PowerFlowData( end # get the network matrices - power_network_matrix = PNM.ABA_Matrix(sys; factorize = true) + power_network_matrix = PNM.ABA_Matrix(sys; factorize=true) aux_network_matrix = PNM.BA_Matrix(sys) # get number of buses and branches @@ -272,7 +272,7 @@ function PowerFlowData( bus_lookup = aux_network_matrix.lookup[1] branch_lookup = aux_network_matrix.lookup[2] - temp_bus_map = Dict{Int, String}( + temp_bus_map = Dict{Int,String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.ACBus, sys) ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) @@ -317,9 +317,9 @@ NOTE: use it for DC power flow computations. function PowerFlowData( ::PTDFDCPowerFlow, sys::PSY.System; - time_steps::Int = 1, - timestep_names::Vector{String} = String[], - check_connectivity::Bool = true) + time_steps::Int=1, + timestep_names::Vector{String}=String[], + check_connectivity::Bool=true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -333,7 +333,7 @@ function PowerFlowData( # get the network matrices power_network_matrix = PNM.PTDF(sys) - aux_network_matrix = PNM.ABA_Matrix(sys; factorize = true) + aux_network_matrix = PNM.ABA_Matrix(sys; factorize=true) # get number of buses and branches n_buses = length(axes(power_network_matrix, 1)) @@ -341,7 +341,7 @@ function PowerFlowData( bus_lookup = power_network_matrix.lookup[1] branch_lookup = power_network_matrix.lookup[2] - temp_bus_map = Dict{Int, String}( + temp_bus_map = Dict{Int,String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) @@ -385,9 +385,9 @@ NOTE: use it for DC power flow computations. function PowerFlowData( ::vPTDFDCPowerFlow, sys::PSY.System; - time_steps::Int = 1, - timestep_names::Vector{String} = String[], - check_connectivity::Bool = true) + time_steps::Int=1, + timestep_names::Vector{String}=String[], + check_connectivity::Bool=true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -401,7 +401,7 @@ function PowerFlowData( # get the network matrices power_network_matrix = PNM.VirtualPTDF(sys) # evaluates an empty virtual PTDF - aux_network_matrix = PNM.ABA_Matrix(sys; factorize = true) + aux_network_matrix = PNM.ABA_Matrix(sys; factorize=true) # get number of buses and branches n_buses = length(axes(power_network_matrix, 2)) @@ -409,7 +409,7 @@ function PowerFlowData( bus_lookup = power_network_matrix.lookup[2] branch_lookup = power_network_matrix.lookup[1] - temp_bus_map = Dict{Int, String}( + temp_bus_map = Dict{Int,String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index c05c7ecd..57468aca 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -30,16 +30,16 @@ solve_ac_powerflow!(sys) solve_ac_powerflow!(sys, method=:newton) ``` """ -function solve_ac_powerflow!(pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, system::PSY.System; kwargs...) +function solve_ac_powerflow!(pf::Union{NLSolveACPowerFlow,KLUACPowerFlow}, system::PSY.System; kwargs...) #Save per-unit flag settings_unit_cache = deepcopy(system.units_settings.unit_system) #Work in System per unit PSY.set_units_base_system!(system, "SYSTEM_BASE") check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false) data = PowerFlowData( - NLSolveACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits), + NLSolveACPowerFlow(; check_reactive_power_limits=check_reactive_power_limits), system; - check_connectivity = get(kwargs, :check_connectivity, true), + check_connectivity=get(kwargs, :check_connectivity, true), ) max_iterations = DEFAULT_MAX_REDISTRIBUTION_ITERATIONS converged, x = _solve_powerflow!(pf, data, check_reactive_power_limits; kwargs...) @@ -68,7 +68,7 @@ res = solve_powerflow(sys, method=:newton) ``` """ function solve_powerflow( - pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + pf::Union{NLSolveACPowerFlow,KLUACPowerFlow}, system::PSY.System; kwargs..., ) @@ -79,7 +79,7 @@ function solve_powerflow( data = PowerFlowData( pf, system; - check_connectivity = get(kwargs, :check_connectivity, true), + check_connectivity=get(kwargs, :check_connectivity, true), ) converged, x = _solve_powerflow!(pf, data, pf.check_reactive_power_limits; kwargs...) @@ -101,7 +101,7 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) within_limits = true for (ix, b) in enumerate(data.bus_type) if b == PSY.ACBusTypes.PV - Q_gen = zero[2 * ix - 1] + Q_gen = zero[2*ix-1] else continue end @@ -124,7 +124,7 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) end function _solve_powerflow!( - pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + pf::Union{NLSolveACPowerFlow,KLUACPowerFlow}, data::ACPowerFlowData, check_reactive_power_limits; nlsolve_kwargs..., @@ -235,16 +235,16 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k 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] + 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] + 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] + x[2*ix-1] = Vm[ix] + x[2*ix] = Va[ix] end end diff --git a/src/post_processing.jl b/src/post_processing.jl index 9b38e858..f8473b2f 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -147,7 +147,7 @@ function _get_fixed_admittance_power( if (l.bus == b) Vm_squared = if b.bustype == PSY.ACBusTypes.PQ - result[2 * ix - 1]^2 + result[2*ix-1]^2 else PSY.get_magnitude(b)^2 end @@ -163,11 +163,11 @@ function _get_limits_for_power_distribution(gen::PSY.StaticInjection) end function _get_limits_for_power_distribution(gen::PSY.RenewableDispatch) - return (min = 0.0, max = PSY.get_max_active_power(gen)) + return (min=0.0, max=PSY.get_max_active_power(gen)) end function _get_limits_for_power_distribution(gen::PSY.Storage) - return (min = 0.0, max = PSY.get_output_active_power_limits(gen).max) + return (min=0.0, max=PSY.get_output_active_power_limits(gen).max) end function _power_redistribution_ref( @@ -188,8 +188,8 @@ function _power_redistribution_ref( elseif length(sources) > 1 && length(non_source_devices) == 0 Psources = sum(PSY.get_active_power.(sources)) Qsources = sum(PSY.get_reactive_power.(sources)) - if isapprox(Psources, P_gen; atol = 0.001) && - isapprox(Qsources, Q_gen; atol = 0.001) + if isapprox(Psources, P_gen; atol=0.001) && + isapprox(Qsources, Q_gen; atol=0.001) @warn "Only sources found at reference bus --- no redistribution of active or reactive power will take place" return else @@ -204,7 +204,7 @@ function _power_redistribution_ref( return elseif length(devices_) > 1 devices = - sort(collect(devices_); by = x -> _get_limits_for_power_distribution(x).max) + sort(collect(devices_); by=x -> _get_limits_for_power_distribution(x).max) else error("No devices in bus $(PSY.get_name(bus))") end @@ -226,14 +226,14 @@ function _power_redistribution_ref( p_residual -= p_set_point end - if !isapprox(p_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + if !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) @debug "Ref Bus voltage residual $p_residual" removed_power = sum([ g.max for g in _get_limits_for_power_distribution.(devices[units_at_limit]) ]) reallocated_p = 0.0 it = 0 - while !isapprox(p_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + while !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) if length(devices) == length(units_at_limit) + 1 @warn "all devices at the active Power Limit" break @@ -255,7 +255,7 @@ function _power_redistribution_ref( reallocated_p += p_frac end p_residual -= reallocated_p - if isapprox(p_residual, 0; atol = ISAPPROX_ZERO_TOLERANCE) + if isapprox(p_residual, 0; atol=ISAPPROX_ZERO_TOLERANCE) break end it += 1 @@ -263,7 +263,7 @@ function _power_redistribution_ref( break end end - if !isapprox(p_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + if !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) remaining_unit_index = setdiff(1:length(devices), units_at_limit) @assert length(remaining_unit_index) == 1 remaining_unit_index device = devices[remaining_unit_index[1]] @@ -298,7 +298,7 @@ function _reactive_power_redistribution_pv( @warn "Found sources and non-source devices at the same bus. Reactive power re-distribution is not well defined for this case. Source reactive power will remain unchanged and remaining reactive power will be re-distributed among non-source devices." elseif length(sources) > 1 && length(non_source_devices) == 0 Qsources = sum(PSY.get_reactive_power.(sources)) - if isapprox(Qsources, Q_gen; atol = 0.001) + if isapprox(Qsources, Q_gen; atol=0.001) @warn "Only sources found at PV bus --- no redistribution of reactive power will take place" return else @@ -311,14 +311,14 @@ function _reactive_power_redistribution_pv( PSY.set_reactive_power!(first(devices_), Q_gen) return elseif length(devices_) > 1 - devices = sort(collect(devices_); by = x -> PSY.get_max_reactive_power(x)) + devices = sort(collect(devices_); by=x -> PSY.get_max_reactive_power(x)) else error("No devices in bus $(PSY.get_name(bus))") end total_active_power = sum(PSY.get_active_power.(devices)) - if isapprox(total_active_power, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + if isapprox(total_active_power, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) @debug "Total Active Power Output at the bus is $(total_active_power). Using Unit's Base Power" sum_basepower = sum(PSY.get_base_power.(devices)) for d in devices @@ -333,8 +333,8 @@ function _reactive_power_redistribution_pv( for (ix, d) in enumerate(devices) q_limits = PSY.get_reactive_power_limits(d) - if isapprox(q_limits.max, 0.0; atol = BOUNDS_TOLERANCE) && - isapprox(q_limits.min, 0.0; atol = BOUNDS_TOLERANCE) + if isapprox(q_limits.max, 0.0; atol=BOUNDS_TOLERANCE) && + isapprox(q_limits.min, 0.0; atol=BOUNDS_TOLERANCE) push!(units_at_limit, ix) @info "Unit $(PSY.get_name(d)) has no Q control capability. Q_max = $(q_limits.max) Q_min = $(q_limits.min)" continue @@ -361,14 +361,14 @@ function _reactive_power_redistribution_pv( PSY.set_reactive_power!(d, q_set_point) q_residual -= q_set_point - if isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + if isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) break end end - if !isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + if !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) it = 0 - while !isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + while !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) if length(devices) == length(units_at_limit) + 1 @debug "Only one device not at the limit in Bus" break @@ -407,7 +407,7 @@ function _reactive_power_redistribution_pv( PSY.set_reactive_power!(d, q_set_point) end q_residual -= reallocated_q - if isapprox(q_residual, 0; atol = ISAPPROX_ZERO_TOLERANCE) + if isapprox(q_residual, 0; atol=ISAPPROX_ZERO_TOLERANCE) break end it += 1 @@ -419,7 +419,7 @@ function _reactive_power_redistribution_pv( end # Last attempt to allocate reactive power - if !isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) + if !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) remaining_unit_index = setdiff(1:length(devices), units_at_limit) @assert length(remaining_unit_index) == 1 remaining_unit_index device = devices[remaining_unit_index[1]] @@ -436,7 +436,7 @@ function _reactive_power_redistribution_pv( @assert isapprox( sum(PSY.get_reactive_power.(devices)), Q_gen; - atol = ISAPPROX_ZERO_TOLERANCE, + atol=ISAPPROX_ZERO_TOLERANCE, ) return @@ -451,21 +451,21 @@ function write_powerflow_solution!( max_iterations::Int, ) buses = enumerate( - sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x)), + sort!(collect(PSY.get_components(PSY.Bus, sys)); by=x -> PSY.get_number(x)), ) for (ix, bus) in buses if bus.bustype == PSY.ACBusTypes.REF - P_gen = result[2 * ix - 1] - Q_gen = result[2 * ix] + P_gen = result[2*ix-1] + Q_gen = result[2*ix] _power_redistribution_ref(sys, P_gen, Q_gen, bus, max_iterations) elseif bus.bustype == PSY.ACBusTypes.PV - Q_gen = result[2 * ix - 1] - bus.angle = result[2 * ix] + Q_gen = result[2*ix-1] + bus.angle = result[2*ix] _reactive_power_redistribution_pv(sys, Q_gen, bus, max_iterations) elseif bus.bustype == PSY.ACBusTypes.PQ - Vm = result[2 * ix - 1] - θ = result[2 * ix] + Vm = result[2*ix-1] + θ = result[2*ix] PSY.set_magnitude!(bus, Vm) PSY.set_angle!(bus, θ) end @@ -481,7 +481,7 @@ function _get_branches_buses(data::ABAPowerFlowData) end # returns list of branches names and buses numbers: PTDF and virtual PTDF case -function _get_branches_buses(data::Union{PTDFPowerFlowData, vPTDFPowerFlowData}) +function _get_branches_buses(data::Union{PTDFPowerFlowData,vPTDFPowerFlowData}) PNM.get_branch_ax(data.power_network_matrix) PNM.get_bus_ax(data.power_network_matrix) return PNM.get_branch_ax(data.power_network_matrix), @@ -513,28 +513,28 @@ function _allocate_results_data( end bus_df = DataFrames.DataFrame(; - bus_number = buses, - Vm = bus_magnitude, - θ = bus_angles, - P_gen = P_gen_vect, - P_load = P_load_vect, - P_net = P_gen_vect - P_load_vect, - Q_gen = Q_gen_vect, - Q_load = Q_load_vect, - Q_net = Q_gen_vect - Q_load_vect, + bus_number=buses, + Vm=bus_magnitude, + θ=bus_angles, + P_gen=P_gen_vect, + P_load=P_load_vect, + P_net=P_gen_vect - P_load_vect, + Q_gen=Q_gen_vect, + Q_load=Q_load_vect, + Q_net=Q_gen_vect - Q_load_vect, ) DataFrames.sort!(bus_df, :bus_number) branch_df = DataFrames.DataFrame(; - line_name = branches, - bus_from = from_bus, - bus_to = to_bus, - P_from_to = P_from_to_vect, - Q_from_to = Q_from_to_vect, - P_to_from = P_to_from_vect, - Q_to_from = Q_to_from_vect, - P_losses = zeros(length(branches)), - Q_losses = zeros(length(branches)), + line_name=branches, + bus_from=from_bus, + bus_to=to_bus, + P_from_to=P_from_to_vect, + Q_from_to=Q_from_to_vect, + P_to_from=P_to_from_vect, + Q_to_from=Q_to_from_vect, + P_losses=zeros(length(branches)), + Q_losses=zeros(length(branches)), ) DataFrames.sort!(branch_df, [:bus_from, :bus_to]) @@ -553,7 +553,7 @@ results. container storing the systam information. """ function write_results( - data::Union{PTDFPowerFlowData, vPTDFPowerFlowData, ABAPowerFlowData}, + data::Union{PTDFPowerFlowData,vPTDFPowerFlowData,ABAPowerFlowData}, sys::PSY.System, ) @info("Voltages are exported in pu. Powers are exported in MW/MVAr.") @@ -572,7 +572,7 @@ function write_results( to_bus[i] = PSY.get_number(PSY.get_arc(br).to) end - result_dict = Dict{Union{String, Char}, Dict{String, DataFrames.DataFrame}}() + result_dict = Dict{Union{String,Char},Dict{String,DataFrames.DataFrame}}() for i in 1:length(data.timestep_map) temp_dict = _allocate_results_data( branches, @@ -614,7 +614,7 @@ function write_results( result::Vector{Float64}, ) @info("Voltages are exported in pu. Powers are exported in MW/MVAr.") - buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x)) + buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by=x -> PSY.get_number(x)) N_BUS = length(buses) bus_map = Dict(buses .=> 1:N_BUS) sys_basepower = PSY.get_base_power(sys) @@ -634,16 +634,16 @@ function write_results( if bustype == PSY.ACBusTypes.REF Vm_vect[ix] = data.bus_magnitude[ix] θ_vect[ix] = data.bus_angles[ix] - P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower - Q_gen_vect[ix] = result[2 * ix] * sys_basepower + P_gen_vect[ix] = result[2*ix-1] * sys_basepower + Q_gen_vect[ix] = result[2*ix] * sys_basepower elseif bustype == PSY.ACBusTypes.PV Vm_vect[ix] = data.bus_magnitude[ix] - θ_vect[ix] = result[2 * ix] + θ_vect[ix] = result[2*ix] P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower - Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower + Q_gen_vect[ix] = result[2*ix-1] * sys_basepower elseif bustype == PSY.ACBusTypes.PQ - Vm_vect[ix] = result[2 * ix - 1] - θ_vect[ix] = result[2 * ix] + Vm_vect[ix] = result[2*ix-1] + θ_vect[ix] = result[2*ix] P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower Q_gen_vect[ix] = data.bus_reactivepower_injection[ix] * sys_basepower end @@ -666,27 +666,27 @@ function write_results( end bus_df = DataFrames.DataFrame(; - bus_number = PSY.get_number.(buses), - Vm = Vm_vect, - θ = θ_vect, - P_gen = P_gen_vect, - P_load = P_load_vect, - P_net = P_gen_vect - P_load_vect, - Q_gen = Q_gen_vect, - Q_load = Q_load_vect, - Q_net = Q_gen_vect - Q_load_vect, + bus_number=PSY.get_number.(buses), + Vm=Vm_vect, + θ=θ_vect, + P_gen=P_gen_vect, + P_load=P_load_vect, + P_net=P_gen_vect - P_load_vect, + Q_gen=Q_gen_vect, + Q_load=Q_load_vect, + Q_net=Q_gen_vect - Q_load_vect, ) branch_df = DataFrames.DataFrame(; - line_name = PSY.get_name.(branches), - bus_from = PSY.get_number.(PSY.get_from.(PSY.get_arc.(branches))), - bus_to = PSY.get_number.(PSY.get_to.(PSY.get_arc.(branches))), - P_from_to = P_from_to_vect, - Q_from_to = Q_from_to_vect, - P_to_from = P_to_from_vect, - Q_to_from = Q_to_from_vect, - P_losses = P_from_to_vect + P_to_from_vect, - Q_losses = Q_from_to_vect + Q_to_from_vect, + line_name=PSY.get_name.(branches), + bus_from=PSY.get_number.(PSY.get_from.(PSY.get_arc.(branches))), + bus_to=PSY.get_number.(PSY.get_to.(PSY.get_arc.(branches))), + P_from_to=P_from_to_vect, + Q_from_to=Q_from_to_vect, + P_to_from=P_to_from_vect, + Q_to_from=Q_to_from_vect, + P_losses=P_from_to_vect + P_to_from_vect, + Q_losses=Q_from_to_vect + Q_to_from_vect, ) DataFrames.sort!(branch_df, [:bus_from, :bus_to]) diff --git a/test/test_nlsolve_powerflow.jl b/test/test_nlsolve_powerflow.jl index a11e7124..138b3058 100644 --- a/test/test_nlsolve_powerflow.jl +++ b/test/test_nlsolve_powerflow.jl @@ -30,23 +30,23 @@ -0.2803812119374241 ] - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) - data = PowerFlows.PowerFlowData(ACPowerFlow(), 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!(ACPowerFlow(), data, false) @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6 - @test solve_ac_powerflow!(ACPowerFlow(), sys; method = :newton) + @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(ACPowerFlow(), sys; check_connectivity = true) + 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 "AC Power Flow 14-Bus Line Configurations" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) base_res = solve_powerflow(ACPowerFlow(), sys) branch = first(PSY.get_components(Line, sys)) dyn_branch = DynamicBranch(branch) @@ -55,15 +55,15 @@ end 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) + 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!(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, rtol=0) + @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) + 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(ACPowerFlow(), sys) @@ -79,12 +79,12 @@ end fix_shunt = PSY.FixedAdmittance("FixAdmBus3", true, bus_103, 0.0 + 0.2im) add_component!(sys_3bus, fix_shunt) 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) + @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 "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) + 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") br = get_component(Line, pf_sys5_re, "6") @@ -107,8 +107,8 @@ end ) system = System( file; - bus_name_formatter = x -> strip(string(x["name"])) * "-" * string(x["index"]), - runchecks = false, + bus_name_formatter=x -> strip(string(x["name"])) * "-" * string(x["index"]), + runchecks=false, ) pf_bus_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_bus_results.csv") @@ -135,35 +135,35 @@ end @testset "AC Multiple sources at ref" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b = ACBus(; - number = 1, - name = "01", - bustype = ACBusTypes.REF, - angle = 0.0, - magnitude = 1.1, - voltage_limits = (0.0, 2.0), - base_voltage = 230, + number=1, + name="01", + bustype=ACBusTypes.REF, + angle=0.0, + magnitude=1.1, + voltage_limits=(0.0, 2.0), + base_voltage=230, ) add_component!(sys, b) #Test two sources with equal and opposite P and Q s1 = Source(; - name = "source_1", - available = true, - bus = b, - active_power = 0.5, - reactive_power = 0.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_1", + available=true, + bus=b, + active_power=0.5, + reactive_power=0.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s1) s2 = Source(; - name = "source_2", - available = true, - bus = b, - active_power = -0.5, - reactive_power = -0.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_2", + available=true, + bus=b, + active_power=-0.5, + reactive_power=-0.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s2) @test solve_ac_powerflow!(ACPowerFlow(), sys) @@ -176,70 +176,70 @@ end @testset "AC PowerFlow with Multiple sources at PV" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b1 = ACBus(; - number = 1, - name = "01", - bustype = ACBusTypes.REF, - angle = 0.0, - magnitude = 1.1, - voltage_limits = (0.0, 2.0), - base_voltage = 230, + number=1, + name="01", + bustype=ACBusTypes.REF, + angle=0.0, + magnitude=1.1, + voltage_limits=(0.0, 2.0), + base_voltage=230, ) add_component!(sys, b1) b2 = ACBus(; - number = 2, - name = "02", - bustype = ACBusTypes.PV, - angle = 0.0, - magnitude = 1.1, - voltage_limits = (0.0, 2.0), - base_voltage = 230, + number=2, + name="02", + bustype=ACBusTypes.PV, + angle=0.0, + magnitude=1.1, + voltage_limits=(0.0, 2.0), + base_voltage=230, ) add_component!(sys, b2) - a = Arc(; from = b1, to = b2) + a = Arc(; from=b1, to=b2) add_component!(sys, a) l = Line(; - name = "l1", - available = true, - active_power_flow = 0.0, - reactive_power_flow = 0.0, - arc = a, - r = 1e-3, - x = 1e-3, - b = (from = 0.0, to = 0.0), - rating = 0.0, - angle_limits = (min = -pi / 2, max = pi / 2), + name="l1", + available=true, + active_power_flow=0.0, + reactive_power_flow=0.0, + arc=a, + r=1e-3, + x=1e-3, + b=(from=0.0, to=0.0), + rating=0.0, + angle_limits=(min=-pi / 2, max=pi / 2), ) add_component!(sys, l) #Test two sources with equal and opposite P and Q s1 = Source(; - name = "source_1", - available = true, - bus = b1, - active_power = 0.5, - reactive_power = 0.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_1", + available=true, + bus=b1, + active_power=0.5, + reactive_power=0.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s1) s2 = Source(; - name = "source_2", - available = true, - bus = b2, - active_power = 0.5, - reactive_power = 1.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_2", + available=true, + bus=b2, + active_power=0.5, + reactive_power=1.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s2) s3 = Source(; - name = "source_3", - available = true, - bus = b2, - active_power = -0.5, - reactive_power = -1.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_3", + available=true, + bus=b2, + active_power=-0.5, + reactive_power=-1.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s3) @@ -253,46 +253,46 @@ end @testset "AC PowerFlow Source + non-source at Ref" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b = ACBus(; - number = 1, - name = "01", - bustype = ACBusTypes.REF, - angle = 0.0, - magnitude = 1.1, - voltage_limits = (0.0, 2.0), - base_voltage = 230, + number=1, + name="01", + bustype=ACBusTypes.REF, + angle=0.0, + magnitude=1.1, + voltage_limits=(0.0, 2.0), + base_voltage=230, ) add_component!(sys, b) #Test two sources with equal and opposite P and Q s1 = Source(; - name = "source_1", - available = true, - bus = b, - active_power = 0.5, - reactive_power = 0.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_1", + available=true, + bus=b, + active_power=0.5, + reactive_power=0.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s1) g1 = ThermalStandard(; - name = "init", - available = true, - status = false, - bus = b, - active_power = 0.1, - reactive_power = 0.1, - rating = 0.0, - active_power_limits = (min = 0.0, max = 0.0), - reactive_power_limits = nothing, - ramp_limits = nothing, - operation_cost = ThermalGenerationCost(nothing), - base_power = 100.0, - time_limits = nothing, - prime_mover_type = PrimeMovers.OT, - fuel = ThermalFuels.OTHER, - services = Device[], - dynamic_injector = nothing, - ext = Dict{String, Any}(), + name="init", + available=true, + status=false, + bus=b, + active_power=0.1, + reactive_power=0.1, + rating=0.0, + active_power_limits=(min=0.0, max=0.0), + reactive_power_limits=nothing, + ramp_limits=nothing, + operation_cost=ThermalGenerationCost(nothing), + base_power=100.0, + time_limits=nothing, + prime_mover_type=PrimeMovers.OT, + fuel=ThermalFuels.OTHER, + services=Device[], + dynamic_injector=nothing, + ext=Dict{String,Any}(), ) add_component!(sys, g1) @@ -300,93 +300,93 @@ end @test isapprox( get_active_power(get_component(Source, sys, "source_1")), 0.5; - atol = 0.001, + atol=0.001, ) @test isapprox( get_reactive_power(get_component(Source, sys, "source_1")), 0.1; - atol = 0.001, + atol=0.001, ) end @testset "AC PowerFlow Source + non-source at PV" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b1 = ACBus(; - number = 1, - name = "01", - bustype = ACBusTypes.REF, - angle = 0.0, - magnitude = 1.1, - voltage_limits = (0.0, 2.0), - base_voltage = 230, + number=1, + name="01", + bustype=ACBusTypes.REF, + angle=0.0, + magnitude=1.1, + voltage_limits=(0.0, 2.0), + base_voltage=230, ) add_component!(sys, b1) b2 = ACBus(; - number = 2, - name = "02", - bustype = ACBusTypes.PV, - angle = 0.0, - magnitude = 1.1, - voltage_limits = (0.0, 2.0), - base_voltage = 230, + number=2, + name="02", + bustype=ACBusTypes.PV, + angle=0.0, + magnitude=1.1, + voltage_limits=(0.0, 2.0), + base_voltage=230, ) add_component!(sys, b2) - a = Arc(; from = b1, to = b2) + a = Arc(; from=b1, to=b2) add_component!(sys, a) l = Line(; - name = "l1", - available = true, - active_power_flow = 0.0, - reactive_power_flow = 0.0, - arc = a, - r = 1e-3, - x = 1e-3, - b = (from = 0.0, to = 0.0), - rating = 0.0, - angle_limits = (min = -pi / 2, max = pi / 2), + name="l1", + available=true, + active_power_flow=0.0, + reactive_power_flow=0.0, + arc=a, + r=1e-3, + x=1e-3, + b=(from=0.0, to=0.0), + rating=0.0, + angle_limits=(min=-pi / 2, max=pi / 2), ) add_component!(sys, l) #Test two sources with equal and opposite P and Q s1 = Source(; - name = "source_1", - available = true, - bus = b1, - active_power = 0.5, - reactive_power = 0.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_1", + available=true, + bus=b1, + active_power=0.5, + reactive_power=0.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s1) s2 = Source(; - name = "source_2", - available = true, - bus = b2, - active_power = 0.5, - reactive_power = 1.1, - R_th = 1e-5, - X_th = 1e-5, + name="source_2", + available=true, + bus=b2, + active_power=0.5, + reactive_power=1.1, + R_th=1e-5, + X_th=1e-5, ) add_component!(sys, s2) g1 = ThermalStandard(; - name = "init", - available = true, - status = false, - bus = b2, - active_power = 0.1, - reactive_power = 0.1, - rating = 0.0, - active_power_limits = (min = 0.0, max = 0.0), - reactive_power_limits = nothing, - ramp_limits = nothing, - operation_cost = ThermalGenerationCost(nothing), - base_power = 100.0, - time_limits = nothing, - prime_mover_type = PrimeMovers.OT, - fuel = ThermalFuels.OTHER, - services = Device[], - dynamic_injector = nothing, - ext = Dict{String, Any}(), + name="init", + available=true, + status=false, + bus=b2, + active_power=0.1, + reactive_power=0.1, + rating=0.0, + active_power_limits=(min=0.0, max=0.0), + reactive_power_limits=nothing, + ramp_limits=nothing, + operation_cost=ThermalGenerationCost(nothing), + base_power=100.0, + time_limits=nothing, + prime_mover_type=PrimeMovers.OT, + fuel=ThermalFuels.OTHER, + services=Device[], + dynamic_injector=nothing, + ext=Dict{String,Any}(), ) add_component!(sys, g1) @@ -394,11 +394,11 @@ end @test isapprox( get_active_power(get_component(Source, sys, "source_2")), 0.5; - atol = 0.001, + atol=0.001, ) @test isapprox( get_reactive_power(get_component(Source, sys, "source_2")), 1.1; - atol = 0.001, + atol=0.001, ) end diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl index 5f91cddc..8ae0c305 100644 --- a/test/test_powerflow_data.jl +++ b/test/test_powerflow_data.jl @@ -1,5 +1,5 @@ @testset "PowerFlowData" begin - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) @test PowerFlowData(NLSolveACPowerFlow(), sys) isa PF.ACPowerFlowData @test PowerFlowData(DCPowerFlow(), sys) isa PF.ABAPowerFlowData @test PowerFlowData(PTDFDCPowerFlow(), sys) isa PF.PTDFPowerFlowData @@ -7,13 +7,13 @@ end @testset "PowerFlowData multiperiod" begin - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) time_steps = 24 # TODO: "multiperiod AC still to implement" - @test PowerFlowData(DCPowerFlow(), sys; time_steps = time_steps) isa PF.ABAPowerFlowData - @test PowerFlowData(PTDFDCPowerFlow(), sys; time_steps = time_steps) isa + @test PowerFlowData(DCPowerFlow(), sys; time_steps=time_steps) isa PF.ABAPowerFlowData + @test PowerFlowData(PTDFDCPowerFlow(), sys; time_steps=time_steps) isa PF.PTDFPowerFlowData - @test PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps = time_steps) isa + @test PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps=time_steps) isa PF.vPTDFPowerFlowData end diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl index df7d925d..3e30f111 100644 --- a/test/test_psse_export.jl +++ b/test/test_psse_export.jl @@ -1,5 +1,5 @@ test_psse_export_dir = joinpath(TEST_FILES_DIR, "test_exports") -isdir(test_psse_export_dir) && rm(test_psse_export_dir; recursive = true) +isdir(test_psse_export_dir) && rm(test_psse_export_dir; recursive=true) function _log_assert(result, msg, comparison_name) result || @@ -10,7 +10,7 @@ end If the expression is false, log an error; in any case, pass through the result of the expression. Optionally accepts a name to include in the error log. """ -macro log_assert(ex, comparison_name = nothing) +macro log_assert(ex, comparison_name=nothing) return :(_log_assert($(esc(ex)), $(string(ex)), $(esc(comparison_name)))) end @@ -24,7 +24,7 @@ function compare_df_within_tolerance( comparison_name::String, df1::DataFrame, df2::DataFrame, - default_tol = SYSTEM_REIMPORT_COMPARISON_TOLERANCE; + default_tol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE; kwargs..., ) result = true @@ -38,7 +38,7 @@ function compare_df_within_tolerance( my_tol = (Symbol(colname) in keys(kwargs)) ? kwargs[Symbol(colname)] : default_tol isnothing(my_tol) && continue inner_result = if (my_eltype <: AbstractFloat) - all(isapprox.(col1, col2; atol = my_tol)) + all(isapprox.(col1, col2; atol=my_tol)) else all(IS.isequivalent.(col1, col2)) end @@ -52,7 +52,7 @@ end compare_df_within_tolerance( df1::DataFrame, df2::DataFrame, - default_tol = SYSTEM_REIMPORT_COMPARISON_TOLERANCE; + default_tol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE; kwargs..., ) = compare_df_within_tolerance("unnamed", df1, df2, default_tol; kwargs...) @@ -64,13 +64,13 @@ function reverse_composite_name(name::String) end loose_system_match_fn(a::Float64, b::Float64) = - isapprox(a, b; atol = SYSTEM_REIMPORT_COMPARISON_TOLERANCE) || IS.isequivalent(a, b) + isapprox(a, b; atol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE) || IS.isequivalent(a, b) loose_system_match_fn(a, b) = IS.isequivalent(a, b) function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; - bus_name_mapping = Dict{String, String}(), + bus_name_mapping=Dict{String,String}(), # TODO when possible, also include: PSY.FixedAdmittance, PSY.Arc - include_types = [ + include_types=[ PSY.ACBus, PSY.Area, PSY.Line, @@ -83,14 +83,14 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; PSY.FixedAdmittance, ], # TODO when possible, don't exclude so many fields - exclude_fields = Set([ + exclude_fields=Set([ :ext, :ramp_limits, :time_limits, :services, :angle_limits, ]), - exclude_fields_for_type = Dict( + exclude_fields_for_type=Dict( PSY.ThermalStandard => Set([ :prime_mover_type, :rating, @@ -111,16 +111,16 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; :reactive_power_flow, ]), ), - generator_comparison_fns = [ # TODO rating + generator_comparison_fns=[ # TODO rating PSY.get_name, PSY.get_bus, PSY.get_active_power, PSY.get_reactive_power, PSY.get_base_power, ], - ignore_name_order = true, - ignore_extra_of_type = Union{PSY.ThermalStandard, PSY.StaticLoad}, - exclude_reactive_power = false) + ignore_name_order=true, + ignore_extra_of_type=Union{PSY.ThermalStandard,PSY.StaticLoad}, + exclude_reactive_power=false) result = true if exclude_reactive_power push!(exclude_fields, :reactive_power) @@ -129,7 +129,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; end # Compare everything about the systems except the actual components - result &= IS.compare_values(sys1, sys2; exclude = [:data]) + result &= IS.compare_values(sys1, sys2; exclude=[:data]) # Compare the components by concrete type for my_type in include_types @@ -178,7 +178,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; loose_system_match_fn, comp1, comp2; - exclude = my_excludes, + exclude=my_excludes, ) result &= comparison if !comparison @@ -189,7 +189,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; end # Extra checks for other types of generators - GenSource = Union{Generator, Source} + GenSource = Union{Generator,Source} gen1_names = sort(PSY.get_name.(PSY.get_components(GenSource, sys1))) gen2_names = sort(PSY.get_name.(PSY.get_components(GenSource, sys2))) if gen1_names != gen2_names @@ -205,13 +205,13 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; ) # Skip pairs we've already compared # e.g., if they're both ThermalStandards, we've already compared them - any(Union{typeof(gen1), typeof(gen2)} .<: include_types) && continue + any(Union{typeof(gen1),typeof(gen2)} .<: include_types) && continue for comp_fn in generator_comparison_fns comparison = IS.compare_values( loose_system_match_fn, comp_fn(gen1), comp_fn(gen2); - exclude = exclude_fields, + exclude=exclude_fields, ) result &= comparison if !comparison @@ -222,7 +222,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; return result end -function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow = false) +function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow=false) result1 = solve_powerflow(NLSolveACPowerFlow(), sys1) result2 = solve_powerflow(NLSolveACPowerFlow(), sys2) reactive_power_tol = @@ -232,8 +232,8 @@ function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow = fal @test compare_df_within_tolerance("flow_results", sort(result1["flow_results"], names(result1["flow_results"])[2:end]), sort(result2["flow_results"], names(result2["flow_results"])[2:end]), - POWERFLOW_COMPARISON_TOLERANCE; line_name = nothing, Q_to_from = reactive_power_tol, - Q_from_to = reactive_power_tol, Q_losses = reactive_power_tol) + POWERFLOW_COMPARISON_TOLERANCE; line_name=nothing, Q_to_from=reactive_power_tol, + Q_from_to=reactive_power_tol, Q_losses=reactive_power_tol) end function read_system_and_metadata(raw_path, metadata_path) @@ -250,8 +250,8 @@ function test_psse_round_trip( exporter::PSSEExporter, scenario_name::AbstractString, export_location::AbstractString; - do_power_flow_test = true, - exclude_reactive_flow = false, + do_power_flow_test=true, + exclude_reactive_flow=false, ) raw_path, metadata_path = get_psse_export_paths(joinpath(export_location, scenario_name)) @@ -265,7 +265,7 @@ function test_psse_round_trip( sys2, sys2_metadata = read_system_and_metadata(raw_path, metadata_path) @test compare_systems_loosely(sys, sys2) do_power_flow_test && - test_power_flow(sys, sys2; exclude_reactive_flow = exclude_reactive_flow) + test_power_flow(sys, sys2; exclude_reactive_flow=exclude_reactive_flow) end "Test that the two raw files are exactly identical and the two metadata files parse to identical JSON" @@ -274,7 +274,7 @@ function test_psse_export_strict_equality( metadata1, raw2, metadata2; - exclude_metadata_keys = ["case_name"], + exclude_metadata_keys=["case_name"], ) open(raw1, "r") do handle1 open(raw2, "r") do handle2 @@ -329,7 +329,7 @@ end export_location = joinpath(test_psse_export_dir, "v33", "system_240") exporter = PSSEExporter(sys, :v33, export_location) test_psse_round_trip(sys, exporter, "basic", export_location; - exclude_reactive_flow = true) # TODO why is reactive flow not matching? + exclude_reactive_flow=true) # TODO why is reactive flow not matching? # Exporting the exact same thing again should result in the exact same files write_export(exporter, "basic2") @@ -360,7 +360,7 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys2)) - test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? + test_power_flow(sys2, reread_sys2; exclude_reactive_flow=true) # TODO why is reactive flow not matching? end @testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" begin @@ -374,7 +374,7 @@ end export_location = joinpath(test_psse_export_dir, "v33", "rts_gmlc") exporter = PSSEExporter(sys, :v33, export_location) test_psse_round_trip(sys, exporter, "basic", export_location; - exclude_reactive_flow = true) # TODO why is reactive flow not matching? + exclude_reactive_flow=true) # TODO why is reactive flow not matching? # Exporting the exact same thing again should result in the exact same files write_export(exporter, "basic2") @@ -404,7 +404,7 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys2)) - test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? + test_power_flow(sys2, reread_sys2; exclude_reactive_flow=true) # TODO why is reactive flow not matching? # Updating with changed value should result in a different reimport (PowerFlowData version) exporter = PSSEExporter(sys, :v33, export_location) @@ -417,16 +417,16 @@ end reread_sys3, sys3_metadata = read_system_and_metadata(joinpath(export_location, "basic5")) @test compare_systems_loosely(sys2, reread_sys3; - exclude_reactive_power = true) # TODO why is reactive power not matching? + exclude_reactive_power=true) # TODO why is reactive power not matching? @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys3)) - test_power_flow(sys2, reread_sys3; exclude_reactive_flow = true) # TODO why is reactive flow not matching? + test_power_flow(sys2, reread_sys3; exclude_reactive_flow=true) # TODO why is reactive flow not matching? # Exporting with write_comments should be comparable to original system - exporter = PSSEExporter(sys, :v33, export_location; write_comments = true) + exporter = PSSEExporter(sys, :v33, export_location; write_comments=true) test_psse_round_trip(sys, exporter, "basic6", export_location; - exclude_reactive_flow = true) # TODO why is reactive flow not matching? + exclude_reactive_flow=true) # TODO why is reactive flow not matching? end @testset "Test exporter helper functions" begin From c6995565ec2b2b959d76fd9944b927b35641843d Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Fri, 6 Dec 2024 15:13:12 -0700 Subject: [PATCH 06/12] formatting --- src/PowerFlowData.jl | 58 ++--- src/nlsolve_ac_powerflow.jl | 47 ++-- src/post_processing.jl | 150 ++++++------- test/test_nlsolve_powerflow.jl | 387 +++++++++++++++++---------------- test/test_powerflow_data.jl | 10 +- test/test_psse_export.jl | 68 +++--- 6 files changed, 373 insertions(+), 347 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index ac0f198f..00c8c4a2 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -66,11 +66,11 @@ flows and angles, as well as these ones. - `neighbors::Vector{Set{Int}}`: Vector with the sets of adjacent buses. """ struct PowerFlowData{ - M<:PNM.PowerNetworkMatrix, - N<:Union{PNM.PowerNetworkMatrix,Nothing}, + M <: PNM.PowerNetworkMatrix, + N <: Union{PNM.PowerNetworkMatrix, Nothing}, } <: PowerFlowContainer - bus_lookup::Dict{Int,Int} - branch_lookup::Dict{String,Int} + bus_lookup::Dict{Int, Int} + branch_lookup::Dict{String, Int} bus_activepower_injection::Matrix{Float64} bus_reactivepower_injection::Matrix{Float64} bus_activepower_withdrawals::Matrix{Float64} @@ -81,7 +81,7 @@ struct PowerFlowData{ bus_magnitude::Matrix{Float64} bus_angles::Matrix{Float64} branch_flow_values::Matrix{Float64} - timestep_map::Dict{Int,String} + timestep_map::Dict{Int, String} valid_ix::Vector{Int} power_network_matrix::M aux_network_matrix::N @@ -119,8 +119,8 @@ end # TODO -> MULTI PERIOD: AC Power Flow Data function _calculate_neighbors( Yb::PNM.Ybus{ - Tuple{Vector{Int64},Vector{Int64}}, - Tuple{Dict{Int64,Int64},Dict{Int64,Int64}}, + Tuple{Vector{Int64}, Vector{Int64}}, + Tuple{Dict{Int64, Int64}, Dict{Int64, Int64}}, }, ) I, J, V = SparseArrays.findnz(Yb.data) @@ -157,11 +157,11 @@ NOTE: use it for AC power flow computations. WARNING: functions for the evaluation of the multi-period AC PF still to be implemented. """ function PowerFlowData( - ::Union{NLSolveACPowerFlow,KLUACPowerFlow}, + ::Union{NLSolveACPowerFlow, KLUACPowerFlow}, sys::PSY.System; - time_steps::Int=1, - timestep_names::Vector{String}=String[], - check_connectivity::Bool=true) + time_steps::Int = 1, + timestep_names::Vector{String} = String[], + check_connectivity::Bool = true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -174,7 +174,7 @@ function PowerFlowData( end # get data for calculations - power_network_matrix = PNM.Ybus(sys; check_connectivity=check_connectivity) + power_network_matrix = PNM.Ybus(sys; check_connectivity = check_connectivity) # get number of buses and branches n_buses = length(axes(power_network_matrix, 1)) @@ -185,8 +185,8 @@ function PowerFlowData( n_branches = length(branches) bus_lookup = power_network_matrix.lookup[2] - branch_lookup = Dict{String,Int}() - temp_bus_map = Dict{Int,String}( + branch_lookup = Dict{String, Int}() + temp_bus_map = Dict{Int, String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) branch_types = Vector{DataType}(undef, n_branches) @@ -250,9 +250,9 @@ NOTE: use it for DC power flow computations. function PowerFlowData( ::DCPowerFlow, sys::PSY.System; - time_steps::Int=1, - timestep_names::Vector{String}=String[], - check_connectivity::Bool=true) + time_steps::Int = 1, + timestep_names::Vector{String} = String[], + check_connectivity::Bool = true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -263,7 +263,7 @@ function PowerFlowData( end # get the network matrices - power_network_matrix = PNM.ABA_Matrix(sys; factorize=true) + power_network_matrix = PNM.ABA_Matrix(sys; factorize = true) aux_network_matrix = PNM.BA_Matrix(sys) # get number of buses and branches @@ -272,7 +272,7 @@ function PowerFlowData( bus_lookup = aux_network_matrix.lookup[1] branch_lookup = aux_network_matrix.lookup[2] - temp_bus_map = Dict{Int,String}( + temp_bus_map = Dict{Int, String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.ACBus, sys) ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) @@ -317,9 +317,9 @@ NOTE: use it for DC power flow computations. function PowerFlowData( ::PTDFDCPowerFlow, sys::PSY.System; - time_steps::Int=1, - timestep_names::Vector{String}=String[], - check_connectivity::Bool=true) + time_steps::Int = 1, + timestep_names::Vector{String} = String[], + check_connectivity::Bool = true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -333,7 +333,7 @@ function PowerFlowData( # get the network matrices power_network_matrix = PNM.PTDF(sys) - aux_network_matrix = PNM.ABA_Matrix(sys; factorize=true) + aux_network_matrix = PNM.ABA_Matrix(sys; factorize = true) # get number of buses and branches n_buses = length(axes(power_network_matrix, 1)) @@ -341,7 +341,7 @@ function PowerFlowData( bus_lookup = power_network_matrix.lookup[1] branch_lookup = power_network_matrix.lookup[2] - temp_bus_map = Dict{Int,String}( + temp_bus_map = Dict{Int, String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) @@ -385,9 +385,9 @@ NOTE: use it for DC power flow computations. function PowerFlowData( ::vPTDFDCPowerFlow, sys::PSY.System; - time_steps::Int=1, - timestep_names::Vector{String}=String[], - check_connectivity::Bool=true) + time_steps::Int = 1, + timestep_names::Vector{String} = String[], + check_connectivity::Bool = true) # assign timestep_names # timestep names are then allocated in a dictionary to map matrix columns @@ -401,7 +401,7 @@ function PowerFlowData( # get the network matrices power_network_matrix = PNM.VirtualPTDF(sys) # evaluates an empty virtual PTDF - aux_network_matrix = PNM.ABA_Matrix(sys; factorize=true) + aux_network_matrix = PNM.ABA_Matrix(sys; factorize = true) # get number of buses and branches n_buses = length(axes(power_network_matrix, 2)) @@ -409,7 +409,7 @@ function PowerFlowData( bus_lookup = power_network_matrix.lookup[2] branch_lookup = power_network_matrix.lookup[1] - temp_bus_map = Dict{Int,String}( + temp_bus_map = Dict{Int, String}( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 57468aca..6e08e67a 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -30,16 +30,20 @@ solve_ac_powerflow!(sys) solve_ac_powerflow!(sys, method=:newton) ``` """ -function solve_ac_powerflow!(pf::Union{NLSolveACPowerFlow,KLUACPowerFlow}, system::PSY.System; kwargs...) +function solve_ac_powerflow!( + pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + system::PSY.System; + kwargs..., +) #Save per-unit flag settings_unit_cache = deepcopy(system.units_settings.unit_system) #Work in System per unit PSY.set_units_base_system!(system, "SYSTEM_BASE") check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false) data = PowerFlowData( - NLSolveACPowerFlow(; check_reactive_power_limits=check_reactive_power_limits), + NLSolveACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits), system; - check_connectivity=get(kwargs, :check_connectivity, true), + check_connectivity = get(kwargs, :check_connectivity, true), ) max_iterations = DEFAULT_MAX_REDISTRIBUTION_ITERATIONS converged, x = _solve_powerflow!(pf, data, check_reactive_power_limits; kwargs...) @@ -68,7 +72,7 @@ res = solve_powerflow(sys, method=:newton) ``` """ function solve_powerflow( - pf::Union{NLSolveACPowerFlow,KLUACPowerFlow}, + pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, system::PSY.System; kwargs..., ) @@ -79,7 +83,7 @@ function solve_powerflow( data = PowerFlowData( pf, system; - check_connectivity=get(kwargs, :check_connectivity, true), + check_connectivity = get(kwargs, :check_connectivity, true), ) converged, x = _solve_powerflow!(pf, data, pf.check_reactive_power_limits; kwargs...) @@ -101,7 +105,7 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) within_limits = true for (ix, b) in enumerate(data.bus_type) if b == PSY.ACBusTypes.PV - Q_gen = zero[2*ix-1] + Q_gen = zero[2 * ix - 1] else continue end @@ -124,7 +128,7 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) end function _solve_powerflow!( - pf::Union{NLSolveACPowerFlow,KLUACPowerFlow}, + pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, data::ACPowerFlowData, check_reactive_power_limits; nlsolve_kwargs..., @@ -145,7 +149,11 @@ function _solve_powerflow!( end end -function _nlsolve_powerflow(pf::NLSolveACPowerFlow, data::ACPowerFlowData; nlsolve_kwargs...) +function _nlsolve_powerflow( + pf::NLSolveACPowerFlow, + data::ACPowerFlowData; + nlsolve_kwargs..., +) pf = PolarPowerFlow(data) J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) @@ -176,7 +184,9 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k 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[:]) + 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])] @@ -208,8 +218,8 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k #dx = - (J_function.Jv \ F) Va[pv] .+= dx[1:npv] - Va[pq] .+= dx[(npv+1):(npv+npq)] - Vm[pq] .+= dx[(npv+npq+1):(npv+2*npq)] + Va[pq] .+= dx[(npv + 1):(npv + npq)] + Vm[pq] .+= dx[(npv + npq + 1):(npv + 2 * npq)] V = Vm .* exp.(1im * Va) @@ -221,7 +231,6 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k converged = LinearAlgebra.norm(F, Inf) < tol end - if !converged @error("The powerflow solver with KLU did not converge after $i iterations") else @@ -230,21 +239,21 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k # 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)] + 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] + 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] + 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] + x[2 * ix - 1] = Vm[ix] + x[2 * ix] = Va[ix] end end diff --git a/src/post_processing.jl b/src/post_processing.jl index f8473b2f..8916a378 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -147,7 +147,7 @@ function _get_fixed_admittance_power( if (l.bus == b) Vm_squared = if b.bustype == PSY.ACBusTypes.PQ - result[2*ix-1]^2 + result[2 * ix - 1]^2 else PSY.get_magnitude(b)^2 end @@ -163,11 +163,11 @@ function _get_limits_for_power_distribution(gen::PSY.StaticInjection) end function _get_limits_for_power_distribution(gen::PSY.RenewableDispatch) - return (min=0.0, max=PSY.get_max_active_power(gen)) + return (min = 0.0, max = PSY.get_max_active_power(gen)) end function _get_limits_for_power_distribution(gen::PSY.Storage) - return (min=0.0, max=PSY.get_output_active_power_limits(gen).max) + return (min = 0.0, max = PSY.get_output_active_power_limits(gen).max) end function _power_redistribution_ref( @@ -188,8 +188,8 @@ function _power_redistribution_ref( elseif length(sources) > 1 && length(non_source_devices) == 0 Psources = sum(PSY.get_active_power.(sources)) Qsources = sum(PSY.get_reactive_power.(sources)) - if isapprox(Psources, P_gen; atol=0.001) && - isapprox(Qsources, Q_gen; atol=0.001) + if isapprox(Psources, P_gen; atol = 0.001) && + isapprox(Qsources, Q_gen; atol = 0.001) @warn "Only sources found at reference bus --- no redistribution of active or reactive power will take place" return else @@ -204,7 +204,7 @@ function _power_redistribution_ref( return elseif length(devices_) > 1 devices = - sort(collect(devices_); by=x -> _get_limits_for_power_distribution(x).max) + sort(collect(devices_); by = x -> _get_limits_for_power_distribution(x).max) else error("No devices in bus $(PSY.get_name(bus))") end @@ -226,14 +226,14 @@ function _power_redistribution_ref( p_residual -= p_set_point end - if !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + if !isapprox(p_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) @debug "Ref Bus voltage residual $p_residual" removed_power = sum([ g.max for g in _get_limits_for_power_distribution.(devices[units_at_limit]) ]) reallocated_p = 0.0 it = 0 - while !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + while !isapprox(p_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) if length(devices) == length(units_at_limit) + 1 @warn "all devices at the active Power Limit" break @@ -255,7 +255,7 @@ function _power_redistribution_ref( reallocated_p += p_frac end p_residual -= reallocated_p - if isapprox(p_residual, 0; atol=ISAPPROX_ZERO_TOLERANCE) + if isapprox(p_residual, 0; atol = ISAPPROX_ZERO_TOLERANCE) break end it += 1 @@ -263,7 +263,7 @@ function _power_redistribution_ref( break end end - if !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + if !isapprox(p_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) remaining_unit_index = setdiff(1:length(devices), units_at_limit) @assert length(remaining_unit_index) == 1 remaining_unit_index device = devices[remaining_unit_index[1]] @@ -298,7 +298,7 @@ function _reactive_power_redistribution_pv( @warn "Found sources and non-source devices at the same bus. Reactive power re-distribution is not well defined for this case. Source reactive power will remain unchanged and remaining reactive power will be re-distributed among non-source devices." elseif length(sources) > 1 && length(non_source_devices) == 0 Qsources = sum(PSY.get_reactive_power.(sources)) - if isapprox(Qsources, Q_gen; atol=0.001) + if isapprox(Qsources, Q_gen; atol = 0.001) @warn "Only sources found at PV bus --- no redistribution of reactive power will take place" return else @@ -311,14 +311,14 @@ function _reactive_power_redistribution_pv( PSY.set_reactive_power!(first(devices_), Q_gen) return elseif length(devices_) > 1 - devices = sort(collect(devices_); by=x -> PSY.get_max_reactive_power(x)) + devices = sort(collect(devices_); by = x -> PSY.get_max_reactive_power(x)) else error("No devices in bus $(PSY.get_name(bus))") end total_active_power = sum(PSY.get_active_power.(devices)) - if isapprox(total_active_power, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + if isapprox(total_active_power, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) @debug "Total Active Power Output at the bus is $(total_active_power). Using Unit's Base Power" sum_basepower = sum(PSY.get_base_power.(devices)) for d in devices @@ -333,8 +333,8 @@ function _reactive_power_redistribution_pv( for (ix, d) in enumerate(devices) q_limits = PSY.get_reactive_power_limits(d) - if isapprox(q_limits.max, 0.0; atol=BOUNDS_TOLERANCE) && - isapprox(q_limits.min, 0.0; atol=BOUNDS_TOLERANCE) + if isapprox(q_limits.max, 0.0; atol = BOUNDS_TOLERANCE) && + isapprox(q_limits.min, 0.0; atol = BOUNDS_TOLERANCE) push!(units_at_limit, ix) @info "Unit $(PSY.get_name(d)) has no Q control capability. Q_max = $(q_limits.max) Q_min = $(q_limits.min)" continue @@ -361,14 +361,14 @@ function _reactive_power_redistribution_pv( PSY.set_reactive_power!(d, q_set_point) q_residual -= q_set_point - if isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + if isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) break end end - if !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + if !isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) it = 0 - while !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + while !isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) if length(devices) == length(units_at_limit) + 1 @debug "Only one device not at the limit in Bus" break @@ -407,7 +407,7 @@ function _reactive_power_redistribution_pv( PSY.set_reactive_power!(d, q_set_point) end q_residual -= reallocated_q - if isapprox(q_residual, 0; atol=ISAPPROX_ZERO_TOLERANCE) + if isapprox(q_residual, 0; atol = ISAPPROX_ZERO_TOLERANCE) break end it += 1 @@ -419,7 +419,7 @@ function _reactive_power_redistribution_pv( end # Last attempt to allocate reactive power - if !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE) + if !isapprox(q_residual, 0.0; atol = ISAPPROX_ZERO_TOLERANCE) remaining_unit_index = setdiff(1:length(devices), units_at_limit) @assert length(remaining_unit_index) == 1 remaining_unit_index device = devices[remaining_unit_index[1]] @@ -436,7 +436,7 @@ function _reactive_power_redistribution_pv( @assert isapprox( sum(PSY.get_reactive_power.(devices)), Q_gen; - atol=ISAPPROX_ZERO_TOLERANCE, + atol = ISAPPROX_ZERO_TOLERANCE, ) return @@ -451,21 +451,21 @@ function write_powerflow_solution!( max_iterations::Int, ) buses = enumerate( - sort!(collect(PSY.get_components(PSY.Bus, sys)); by=x -> PSY.get_number(x)), + sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x)), ) for (ix, bus) in buses if bus.bustype == PSY.ACBusTypes.REF - P_gen = result[2*ix-1] - Q_gen = result[2*ix] + P_gen = result[2 * ix - 1] + Q_gen = result[2 * ix] _power_redistribution_ref(sys, P_gen, Q_gen, bus, max_iterations) elseif bus.bustype == PSY.ACBusTypes.PV - Q_gen = result[2*ix-1] - bus.angle = result[2*ix] + Q_gen = result[2 * ix - 1] + bus.angle = result[2 * ix] _reactive_power_redistribution_pv(sys, Q_gen, bus, max_iterations) elseif bus.bustype == PSY.ACBusTypes.PQ - Vm = result[2*ix-1] - θ = result[2*ix] + Vm = result[2 * ix - 1] + θ = result[2 * ix] PSY.set_magnitude!(bus, Vm) PSY.set_angle!(bus, θ) end @@ -481,7 +481,7 @@ function _get_branches_buses(data::ABAPowerFlowData) end # returns list of branches names and buses numbers: PTDF and virtual PTDF case -function _get_branches_buses(data::Union{PTDFPowerFlowData,vPTDFPowerFlowData}) +function _get_branches_buses(data::Union{PTDFPowerFlowData, vPTDFPowerFlowData}) PNM.get_branch_ax(data.power_network_matrix) PNM.get_bus_ax(data.power_network_matrix) return PNM.get_branch_ax(data.power_network_matrix), @@ -513,28 +513,28 @@ function _allocate_results_data( end bus_df = DataFrames.DataFrame(; - bus_number=buses, - Vm=bus_magnitude, - θ=bus_angles, - P_gen=P_gen_vect, - P_load=P_load_vect, - P_net=P_gen_vect - P_load_vect, - Q_gen=Q_gen_vect, - Q_load=Q_load_vect, - Q_net=Q_gen_vect - Q_load_vect, + bus_number = buses, + Vm = bus_magnitude, + θ = bus_angles, + P_gen = P_gen_vect, + P_load = P_load_vect, + P_net = P_gen_vect - P_load_vect, + Q_gen = Q_gen_vect, + Q_load = Q_load_vect, + Q_net = Q_gen_vect - Q_load_vect, ) DataFrames.sort!(bus_df, :bus_number) branch_df = DataFrames.DataFrame(; - line_name=branches, - bus_from=from_bus, - bus_to=to_bus, - P_from_to=P_from_to_vect, - Q_from_to=Q_from_to_vect, - P_to_from=P_to_from_vect, - Q_to_from=Q_to_from_vect, - P_losses=zeros(length(branches)), - Q_losses=zeros(length(branches)), + line_name = branches, + bus_from = from_bus, + bus_to = to_bus, + P_from_to = P_from_to_vect, + Q_from_to = Q_from_to_vect, + P_to_from = P_to_from_vect, + Q_to_from = Q_to_from_vect, + P_losses = zeros(length(branches)), + Q_losses = zeros(length(branches)), ) DataFrames.sort!(branch_df, [:bus_from, :bus_to]) @@ -553,7 +553,7 @@ results. container storing the systam information. """ function write_results( - data::Union{PTDFPowerFlowData,vPTDFPowerFlowData,ABAPowerFlowData}, + data::Union{PTDFPowerFlowData, vPTDFPowerFlowData, ABAPowerFlowData}, sys::PSY.System, ) @info("Voltages are exported in pu. Powers are exported in MW/MVAr.") @@ -572,7 +572,7 @@ function write_results( to_bus[i] = PSY.get_number(PSY.get_arc(br).to) end - result_dict = Dict{Union{String,Char},Dict{String,DataFrames.DataFrame}}() + result_dict = Dict{Union{String, Char}, Dict{String, DataFrames.DataFrame}}() for i in 1:length(data.timestep_map) temp_dict = _allocate_results_data( branches, @@ -608,13 +608,13 @@ dictionary will therefore feature just one key linked to one DataFrame. vector containing the reults for one single time-period. """ function write_results( - ::Union{NLSolveACPowerFlow,KLUACPowerFlow}, + ::Union{NLSolveACPowerFlow, KLUACPowerFlow}, sys::PSY.System, data::ACPowerFlowData, result::Vector{Float64}, ) @info("Voltages are exported in pu. Powers are exported in MW/MVAr.") - buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by=x -> PSY.get_number(x)) + buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x)) N_BUS = length(buses) bus_map = Dict(buses .=> 1:N_BUS) sys_basepower = PSY.get_base_power(sys) @@ -634,16 +634,16 @@ function write_results( if bustype == PSY.ACBusTypes.REF Vm_vect[ix] = data.bus_magnitude[ix] θ_vect[ix] = data.bus_angles[ix] - P_gen_vect[ix] = result[2*ix-1] * sys_basepower - Q_gen_vect[ix] = result[2*ix] * sys_basepower + P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower + Q_gen_vect[ix] = result[2 * ix] * sys_basepower elseif bustype == PSY.ACBusTypes.PV Vm_vect[ix] = data.bus_magnitude[ix] - θ_vect[ix] = result[2*ix] + θ_vect[ix] = result[2 * ix] P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower - Q_gen_vect[ix] = result[2*ix-1] * sys_basepower + Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower elseif bustype == PSY.ACBusTypes.PQ - Vm_vect[ix] = result[2*ix-1] - θ_vect[ix] = result[2*ix] + Vm_vect[ix] = result[2 * ix - 1] + θ_vect[ix] = result[2 * ix] P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower Q_gen_vect[ix] = data.bus_reactivepower_injection[ix] * sys_basepower end @@ -666,27 +666,27 @@ function write_results( end bus_df = DataFrames.DataFrame(; - bus_number=PSY.get_number.(buses), - Vm=Vm_vect, - θ=θ_vect, - P_gen=P_gen_vect, - P_load=P_load_vect, - P_net=P_gen_vect - P_load_vect, - Q_gen=Q_gen_vect, - Q_load=Q_load_vect, - Q_net=Q_gen_vect - Q_load_vect, + bus_number = PSY.get_number.(buses), + Vm = Vm_vect, + θ = θ_vect, + P_gen = P_gen_vect, + P_load = P_load_vect, + P_net = P_gen_vect - P_load_vect, + Q_gen = Q_gen_vect, + Q_load = Q_load_vect, + Q_net = Q_gen_vect - Q_load_vect, ) branch_df = DataFrames.DataFrame(; - line_name=PSY.get_name.(branches), - bus_from=PSY.get_number.(PSY.get_from.(PSY.get_arc.(branches))), - bus_to=PSY.get_number.(PSY.get_to.(PSY.get_arc.(branches))), - P_from_to=P_from_to_vect, - Q_from_to=Q_from_to_vect, - P_to_from=P_to_from_vect, - Q_to_from=Q_to_from_vect, - P_losses=P_from_to_vect + P_to_from_vect, - Q_losses=Q_from_to_vect + Q_to_from_vect, + line_name = PSY.get_name.(branches), + bus_from = PSY.get_number.(PSY.get_from.(PSY.get_arc.(branches))), + bus_to = PSY.get_number.(PSY.get_to.(PSY.get_arc.(branches))), + P_from_to = P_from_to_vect, + Q_from_to = Q_from_to_vect, + P_to_from = P_to_from_vect, + Q_to_from = Q_to_from_vect, + P_losses = P_from_to_vect + P_to_from_vect, + Q_losses = Q_from_to_vect + Q_to_from_vect, ) DataFrames.sort!(branch_df, [:bus_from, :bus_to]) diff --git a/test/test_nlsolve_powerflow.jl b/test/test_nlsolve_powerflow.jl index 138b3058..9d9d404f 100644 --- a/test/test_nlsolve_powerflow.jl +++ b/test/test_nlsolve_powerflow.jl @@ -1,4 +1,5 @@ -@testset "AC Power Flow 14-Bus testing" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "AC Power Flow 14-Bus testing" for ACPowerFlow in + (NLSolveACPowerFlow, KLUACPowerFlow) result_14 = [ 2.3255081760423684 -0.15529254415401786 @@ -30,40 +31,42 @@ -0.2803812119374241 ] - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) - data = PowerFlows.PowerFlowData(ACPowerFlow(), 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!(ACPowerFlow(), data, false) @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6 - @test solve_ac_powerflow!(ACPowerFlow(), sys; method=:newton) + @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(ACPowerFlow(), sys; check_connectivity=true) + 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 "AC Power Flow 14-Bus Line Configurations" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) +@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(ACPowerFlow(), sys) branch = first(PSY.get_components(Line, sys)) dyn_branch = DynamicBranch(branch) add_component!(sys, dyn_branch) @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 + @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) + 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!(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, rtol=0) + @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) + 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(ACPowerFlow(), sys) @@ -71,7 +74,10 @@ end @test res["flow_results"].P_to_from[4] == 0.0 end -@testset "AC Power Flow 3-Bus Fixed FixedAdmittance testing" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@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") @@ -79,12 +85,13 @@ end fix_shunt = PSY.FixedAdmittance("FixAdmBus3", true, bus_103, 0.0 + 0.2im) add_component!(sys_3bus, fix_shunt) 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) + @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 "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) +@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") br = get_component(Line, pf_sys5_re, "6") @@ -99,7 +106,8 @@ end ) end -@testset "AC Test 240 Case PSS/e results" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "AC Test 240 Case PSS/e results" for ACPowerFlow in + (NLSolveACPowerFlow, KLUACPowerFlow) file = joinpath( TEST_FILES_DIR, "test_data", @@ -107,8 +115,8 @@ end ) system = System( file; - bus_name_formatter=x -> strip(string(x["name"])) * "-" * string(x["index"]), - runchecks=false, + bus_name_formatter = x -> strip(string(x["name"])) * "-" * string(x["index"]), + runchecks = false, ) pf_bus_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_bus_results.csv") @@ -132,114 +140,118 @@ end @test norm(q_diff, 2) / length(q_diff) < DIFF_L2_TOLERANCE end -@testset "AC Multiple sources at ref" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "AC Multiple sources at ref" for ACPowerFlow in + (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b = ACBus(; - number=1, - name="01", - bustype=ACBusTypes.REF, - angle=0.0, - magnitude=1.1, - voltage_limits=(0.0, 2.0), - base_voltage=230, + number = 1, + name = "01", + bustype = ACBusTypes.REF, + angle = 0.0, + magnitude = 1.1, + voltage_limits = (0.0, 2.0), + base_voltage = 230, ) add_component!(sys, b) #Test two sources with equal and opposite P and Q s1 = Source(; - name="source_1", - available=true, - bus=b, - active_power=0.5, - reactive_power=0.1, - R_th=1e-5, - X_th=1e-5, + name = "source_1", + available = true, + bus = b, + active_power = 0.5, + reactive_power = 0.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s1) s2 = Source(; - name="source_2", - available=true, - bus=b, - active_power=-0.5, - reactive_power=-0.1, - R_th=1e-5, - X_th=1e-5, + name = "source_2", + available = true, + bus = b, + active_power = -0.5, + reactive_power = -0.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s2) @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!(ACPowerFlow(), 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" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "AC PowerFlow with Multiple sources at PV" for ACPowerFlow in + (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b1 = ACBus(; - number=1, - name="01", - bustype=ACBusTypes.REF, - angle=0.0, - magnitude=1.1, - voltage_limits=(0.0, 2.0), - base_voltage=230, + number = 1, + name = "01", + bustype = ACBusTypes.REF, + angle = 0.0, + magnitude = 1.1, + voltage_limits = (0.0, 2.0), + base_voltage = 230, ) add_component!(sys, b1) b2 = ACBus(; - number=2, - name="02", - bustype=ACBusTypes.PV, - angle=0.0, - magnitude=1.1, - voltage_limits=(0.0, 2.0), - base_voltage=230, + number = 2, + name = "02", + bustype = ACBusTypes.PV, + angle = 0.0, + magnitude = 1.1, + voltage_limits = (0.0, 2.0), + base_voltage = 230, ) add_component!(sys, b2) - a = Arc(; from=b1, to=b2) + a = Arc(; from = b1, to = b2) add_component!(sys, a) l = Line(; - name="l1", - available=true, - active_power_flow=0.0, - reactive_power_flow=0.0, - arc=a, - r=1e-3, - x=1e-3, - b=(from=0.0, to=0.0), - rating=0.0, - angle_limits=(min=-pi / 2, max=pi / 2), + name = "l1", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + arc = a, + r = 1e-3, + x = 1e-3, + b = (from = 0.0, to = 0.0), + rating = 0.0, + angle_limits = (min = -pi / 2, max = pi / 2), ) add_component!(sys, l) #Test two sources with equal and opposite P and Q s1 = Source(; - name="source_1", - available=true, - bus=b1, - active_power=0.5, - reactive_power=0.1, - R_th=1e-5, - X_th=1e-5, + name = "source_1", + available = true, + bus = b1, + active_power = 0.5, + reactive_power = 0.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s1) s2 = Source(; - name="source_2", - available=true, - bus=b2, - active_power=0.5, - reactive_power=1.1, - R_th=1e-5, - X_th=1e-5, + name = "source_2", + available = true, + bus = b2, + active_power = 0.5, + reactive_power = 1.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s2) s3 = Source(; - name="source_3", - available=true, - bus=b2, - active_power=-0.5, - reactive_power=-1.1, - R_th=1e-5, - X_th=1e-5, + name = "source_3", + available = true, + bus = b2, + active_power = -0.5, + reactive_power = -1.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s3) @@ -247,52 +259,56 @@ end #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!(ACPowerFlow(), 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" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "AC PowerFlow Source + non-source at Ref" for ACPowerFlow in + (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b = ACBus(; - number=1, - name="01", - bustype=ACBusTypes.REF, - angle=0.0, - magnitude=1.1, - voltage_limits=(0.0, 2.0), - base_voltage=230, + number = 1, + name = "01", + bustype = ACBusTypes.REF, + angle = 0.0, + magnitude = 1.1, + voltage_limits = (0.0, 2.0), + base_voltage = 230, ) add_component!(sys, b) #Test two sources with equal and opposite P and Q s1 = Source(; - name="source_1", - available=true, - bus=b, - active_power=0.5, - reactive_power=0.1, - R_th=1e-5, - X_th=1e-5, + name = "source_1", + available = true, + bus = b, + active_power = 0.5, + reactive_power = 0.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s1) g1 = ThermalStandard(; - name="init", - available=true, - status=false, - bus=b, - active_power=0.1, - reactive_power=0.1, - rating=0.0, - active_power_limits=(min=0.0, max=0.0), - reactive_power_limits=nothing, - ramp_limits=nothing, - operation_cost=ThermalGenerationCost(nothing), - base_power=100.0, - time_limits=nothing, - prime_mover_type=PrimeMovers.OT, - fuel=ThermalFuels.OTHER, - services=Device[], - dynamic_injector=nothing, - ext=Dict{String,Any}(), + name = "init", + available = true, + status = false, + bus = b, + active_power = 0.1, + reactive_power = 0.1, + rating = 0.0, + active_power_limits = (min = 0.0, max = 0.0), + reactive_power_limits = nothing, + ramp_limits = nothing, + operation_cost = ThermalGenerationCost(nothing), + base_power = 100.0, + time_limits = nothing, + prime_mover_type = PrimeMovers.OT, + fuel = ThermalFuels.OTHER, + services = Device[], + dynamic_injector = nothing, + ext = Dict{String, Any}(), ) add_component!(sys, g1) @@ -300,93 +316,94 @@ end @test isapprox( get_active_power(get_component(Source, sys, "source_1")), 0.5; - atol=0.001, + atol = 0.001, ) @test isapprox( get_reactive_power(get_component(Source, sys, "source_1")), 0.1; - atol=0.001, + atol = 0.001, ) end -@testset "AC PowerFlow Source + non-source at PV" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "AC PowerFlow Source + non-source at PV" for ACPowerFlow in + (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b1 = ACBus(; - number=1, - name="01", - bustype=ACBusTypes.REF, - angle=0.0, - magnitude=1.1, - voltage_limits=(0.0, 2.0), - base_voltage=230, + number = 1, + name = "01", + bustype = ACBusTypes.REF, + angle = 0.0, + magnitude = 1.1, + voltage_limits = (0.0, 2.0), + base_voltage = 230, ) add_component!(sys, b1) b2 = ACBus(; - number=2, - name="02", - bustype=ACBusTypes.PV, - angle=0.0, - magnitude=1.1, - voltage_limits=(0.0, 2.0), - base_voltage=230, + number = 2, + name = "02", + bustype = ACBusTypes.PV, + angle = 0.0, + magnitude = 1.1, + voltage_limits = (0.0, 2.0), + base_voltage = 230, ) add_component!(sys, b2) - a = Arc(; from=b1, to=b2) + a = Arc(; from = b1, to = b2) add_component!(sys, a) l = Line(; - name="l1", - available=true, - active_power_flow=0.0, - reactive_power_flow=0.0, - arc=a, - r=1e-3, - x=1e-3, - b=(from=0.0, to=0.0), - rating=0.0, - angle_limits=(min=-pi / 2, max=pi / 2), + name = "l1", + available = true, + active_power_flow = 0.0, + reactive_power_flow = 0.0, + arc = a, + r = 1e-3, + x = 1e-3, + b = (from = 0.0, to = 0.0), + rating = 0.0, + angle_limits = (min = -pi / 2, max = pi / 2), ) add_component!(sys, l) #Test two sources with equal and opposite P and Q s1 = Source(; - name="source_1", - available=true, - bus=b1, - active_power=0.5, - reactive_power=0.1, - R_th=1e-5, - X_th=1e-5, + name = "source_1", + available = true, + bus = b1, + active_power = 0.5, + reactive_power = 0.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s1) s2 = Source(; - name="source_2", - available=true, - bus=b2, - active_power=0.5, - reactive_power=1.1, - R_th=1e-5, - X_th=1e-5, + name = "source_2", + available = true, + bus = b2, + active_power = 0.5, + reactive_power = 1.1, + R_th = 1e-5, + X_th = 1e-5, ) add_component!(sys, s2) g1 = ThermalStandard(; - name="init", - available=true, - status=false, - bus=b2, - active_power=0.1, - reactive_power=0.1, - rating=0.0, - active_power_limits=(min=0.0, max=0.0), - reactive_power_limits=nothing, - ramp_limits=nothing, - operation_cost=ThermalGenerationCost(nothing), - base_power=100.0, - time_limits=nothing, - prime_mover_type=PrimeMovers.OT, - fuel=ThermalFuels.OTHER, - services=Device[], - dynamic_injector=nothing, - ext=Dict{String,Any}(), + name = "init", + available = true, + status = false, + bus = b2, + active_power = 0.1, + reactive_power = 0.1, + rating = 0.0, + active_power_limits = (min = 0.0, max = 0.0), + reactive_power_limits = nothing, + ramp_limits = nothing, + operation_cost = ThermalGenerationCost(nothing), + base_power = 100.0, + time_limits = nothing, + prime_mover_type = PrimeMovers.OT, + fuel = ThermalFuels.OTHER, + services = Device[], + dynamic_injector = nothing, + ext = Dict{String, Any}(), ) add_component!(sys, g1) @@ -394,11 +411,11 @@ end @test isapprox( get_active_power(get_component(Source, sys, "source_2")), 0.5; - atol=0.001, + atol = 0.001, ) @test isapprox( get_reactive_power(get_component(Source, sys, "source_2")), 1.1; - atol=0.001, + atol = 0.001, ) end diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl index 8ae0c305..5f91cddc 100644 --- a/test/test_powerflow_data.jl +++ b/test/test_powerflow_data.jl @@ -1,5 +1,5 @@ @testset "PowerFlowData" begin - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) @test PowerFlowData(NLSolveACPowerFlow(), sys) isa PF.ACPowerFlowData @test PowerFlowData(DCPowerFlow(), sys) isa PF.ABAPowerFlowData @test PowerFlowData(PTDFDCPowerFlow(), sys) isa PF.PTDFPowerFlowData @@ -7,13 +7,13 @@ end @testset "PowerFlowData multiperiod" begin - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false) + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) time_steps = 24 # TODO: "multiperiod AC still to implement" - @test PowerFlowData(DCPowerFlow(), sys; time_steps=time_steps) isa PF.ABAPowerFlowData - @test PowerFlowData(PTDFDCPowerFlow(), sys; time_steps=time_steps) isa + @test PowerFlowData(DCPowerFlow(), sys; time_steps = time_steps) isa PF.ABAPowerFlowData + @test PowerFlowData(PTDFDCPowerFlow(), sys; time_steps = time_steps) isa PF.PTDFPowerFlowData - @test PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps=time_steps) isa + @test PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps = time_steps) isa PF.vPTDFPowerFlowData end diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl index 3e30f111..df7d925d 100644 --- a/test/test_psse_export.jl +++ b/test/test_psse_export.jl @@ -1,5 +1,5 @@ test_psse_export_dir = joinpath(TEST_FILES_DIR, "test_exports") -isdir(test_psse_export_dir) && rm(test_psse_export_dir; recursive=true) +isdir(test_psse_export_dir) && rm(test_psse_export_dir; recursive = true) function _log_assert(result, msg, comparison_name) result || @@ -10,7 +10,7 @@ end If the expression is false, log an error; in any case, pass through the result of the expression. Optionally accepts a name to include in the error log. """ -macro log_assert(ex, comparison_name=nothing) +macro log_assert(ex, comparison_name = nothing) return :(_log_assert($(esc(ex)), $(string(ex)), $(esc(comparison_name)))) end @@ -24,7 +24,7 @@ function compare_df_within_tolerance( comparison_name::String, df1::DataFrame, df2::DataFrame, - default_tol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE; + default_tol = SYSTEM_REIMPORT_COMPARISON_TOLERANCE; kwargs..., ) result = true @@ -38,7 +38,7 @@ function compare_df_within_tolerance( my_tol = (Symbol(colname) in keys(kwargs)) ? kwargs[Symbol(colname)] : default_tol isnothing(my_tol) && continue inner_result = if (my_eltype <: AbstractFloat) - all(isapprox.(col1, col2; atol=my_tol)) + all(isapprox.(col1, col2; atol = my_tol)) else all(IS.isequivalent.(col1, col2)) end @@ -52,7 +52,7 @@ end compare_df_within_tolerance( df1::DataFrame, df2::DataFrame, - default_tol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE; + default_tol = SYSTEM_REIMPORT_COMPARISON_TOLERANCE; kwargs..., ) = compare_df_within_tolerance("unnamed", df1, df2, default_tol; kwargs...) @@ -64,13 +64,13 @@ function reverse_composite_name(name::String) end loose_system_match_fn(a::Float64, b::Float64) = - isapprox(a, b; atol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE) || IS.isequivalent(a, b) + isapprox(a, b; atol = SYSTEM_REIMPORT_COMPARISON_TOLERANCE) || IS.isequivalent(a, b) loose_system_match_fn(a, b) = IS.isequivalent(a, b) function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; - bus_name_mapping=Dict{String,String}(), + bus_name_mapping = Dict{String, String}(), # TODO when possible, also include: PSY.FixedAdmittance, PSY.Arc - include_types=[ + include_types = [ PSY.ACBus, PSY.Area, PSY.Line, @@ -83,14 +83,14 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; PSY.FixedAdmittance, ], # TODO when possible, don't exclude so many fields - exclude_fields=Set([ + exclude_fields = Set([ :ext, :ramp_limits, :time_limits, :services, :angle_limits, ]), - exclude_fields_for_type=Dict( + exclude_fields_for_type = Dict( PSY.ThermalStandard => Set([ :prime_mover_type, :rating, @@ -111,16 +111,16 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; :reactive_power_flow, ]), ), - generator_comparison_fns=[ # TODO rating + generator_comparison_fns = [ # TODO rating PSY.get_name, PSY.get_bus, PSY.get_active_power, PSY.get_reactive_power, PSY.get_base_power, ], - ignore_name_order=true, - ignore_extra_of_type=Union{PSY.ThermalStandard,PSY.StaticLoad}, - exclude_reactive_power=false) + ignore_name_order = true, + ignore_extra_of_type = Union{PSY.ThermalStandard, PSY.StaticLoad}, + exclude_reactive_power = false) result = true if exclude_reactive_power push!(exclude_fields, :reactive_power) @@ -129,7 +129,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; end # Compare everything about the systems except the actual components - result &= IS.compare_values(sys1, sys2; exclude=[:data]) + result &= IS.compare_values(sys1, sys2; exclude = [:data]) # Compare the components by concrete type for my_type in include_types @@ -178,7 +178,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; loose_system_match_fn, comp1, comp2; - exclude=my_excludes, + exclude = my_excludes, ) result &= comparison if !comparison @@ -189,7 +189,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; end # Extra checks for other types of generators - GenSource = Union{Generator,Source} + GenSource = Union{Generator, Source} gen1_names = sort(PSY.get_name.(PSY.get_components(GenSource, sys1))) gen2_names = sort(PSY.get_name.(PSY.get_components(GenSource, sys2))) if gen1_names != gen2_names @@ -205,13 +205,13 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; ) # Skip pairs we've already compared # e.g., if they're both ThermalStandards, we've already compared them - any(Union{typeof(gen1),typeof(gen2)} .<: include_types) && continue + any(Union{typeof(gen1), typeof(gen2)} .<: include_types) && continue for comp_fn in generator_comparison_fns comparison = IS.compare_values( loose_system_match_fn, comp_fn(gen1), comp_fn(gen2); - exclude=exclude_fields, + exclude = exclude_fields, ) result &= comparison if !comparison @@ -222,7 +222,7 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; return result end -function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow=false) +function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow = false) result1 = solve_powerflow(NLSolveACPowerFlow(), sys1) result2 = solve_powerflow(NLSolveACPowerFlow(), sys2) reactive_power_tol = @@ -232,8 +232,8 @@ function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow=false @test compare_df_within_tolerance("flow_results", sort(result1["flow_results"], names(result1["flow_results"])[2:end]), sort(result2["flow_results"], names(result2["flow_results"])[2:end]), - POWERFLOW_COMPARISON_TOLERANCE; line_name=nothing, Q_to_from=reactive_power_tol, - Q_from_to=reactive_power_tol, Q_losses=reactive_power_tol) + POWERFLOW_COMPARISON_TOLERANCE; line_name = nothing, Q_to_from = reactive_power_tol, + Q_from_to = reactive_power_tol, Q_losses = reactive_power_tol) end function read_system_and_metadata(raw_path, metadata_path) @@ -250,8 +250,8 @@ function test_psse_round_trip( exporter::PSSEExporter, scenario_name::AbstractString, export_location::AbstractString; - do_power_flow_test=true, - exclude_reactive_flow=false, + do_power_flow_test = true, + exclude_reactive_flow = false, ) raw_path, metadata_path = get_psse_export_paths(joinpath(export_location, scenario_name)) @@ -265,7 +265,7 @@ function test_psse_round_trip( sys2, sys2_metadata = read_system_and_metadata(raw_path, metadata_path) @test compare_systems_loosely(sys, sys2) do_power_flow_test && - test_power_flow(sys, sys2; exclude_reactive_flow=exclude_reactive_flow) + test_power_flow(sys, sys2; exclude_reactive_flow = exclude_reactive_flow) end "Test that the two raw files are exactly identical and the two metadata files parse to identical JSON" @@ -274,7 +274,7 @@ function test_psse_export_strict_equality( metadata1, raw2, metadata2; - exclude_metadata_keys=["case_name"], + exclude_metadata_keys = ["case_name"], ) open(raw1, "r") do handle1 open(raw2, "r") do handle2 @@ -329,7 +329,7 @@ end export_location = joinpath(test_psse_export_dir, "v33", "system_240") exporter = PSSEExporter(sys, :v33, export_location) test_psse_round_trip(sys, exporter, "basic", export_location; - exclude_reactive_flow=true) # TODO why is reactive flow not matching? + exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Exporting the exact same thing again should result in the exact same files write_export(exporter, "basic2") @@ -360,7 +360,7 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys2)) - test_power_flow(sys2, reread_sys2; exclude_reactive_flow=true) # TODO why is reactive flow not matching? + test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? end @testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" begin @@ -374,7 +374,7 @@ end export_location = joinpath(test_psse_export_dir, "v33", "rts_gmlc") exporter = PSSEExporter(sys, :v33, export_location) test_psse_round_trip(sys, exporter, "basic", export_location; - exclude_reactive_flow=true) # TODO why is reactive flow not matching? + exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Exporting the exact same thing again should result in the exact same files write_export(exporter, "basic2") @@ -404,7 +404,7 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys2)) - test_power_flow(sys2, reread_sys2; exclude_reactive_flow=true) # TODO why is reactive flow not matching? + test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Updating with changed value should result in a different reimport (PowerFlowData version) exporter = PSSEExporter(sys, :v33, export_location) @@ -417,16 +417,16 @@ end reread_sys3, sys3_metadata = read_system_and_metadata(joinpath(export_location, "basic5")) @test compare_systems_loosely(sys2, reread_sys3; - exclude_reactive_power=true) # TODO why is reactive power not matching? + exclude_reactive_power = true) # TODO why is reactive power not matching? @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys3)) - test_power_flow(sys2, reread_sys3; exclude_reactive_flow=true) # TODO why is reactive flow not matching? + test_power_flow(sys2, reread_sys3; exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Exporting with write_comments should be comparable to original system - exporter = PSSEExporter(sys, :v33, export_location; write_comments=true) + exporter = PSSEExporter(sys, :v33, export_location; write_comments = true) test_psse_round_trip(sys, exporter, "basic6", export_location; - exclude_reactive_flow=true) # TODO why is reactive flow not matching? + exclude_reactive_flow = true) # TODO why is reactive flow not matching? end @testset "Test exporter helper functions" begin From b0bc56522f9254cfcd9166c76002d13a0e9399a4 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Mon, 16 Dec 2024 14:45:09 -0700 Subject: [PATCH 07/12] Addressed review comments except for using functor --- Project.toml | 2 +- src/PowerFlowData.jl | 4 +- src/PowerFlows.jl | 4 +- src/definitions.jl | 3 + ...ac_powerflow.jl => newton_ac_powerflow.jl} | 45 ++++++----- src/post_processing.jl | 2 +- src/powerflow_types.jl | 12 ++- ...werflow.jl => test_newton_ac_powerflow.jl} | 74 +++++++++++-------- test/test_powerflow_data.jl | 15 ++-- test/test_psse_export.jl | 33 +++++---- 10 files changed, 111 insertions(+), 83 deletions(-) rename src/{nlsolve_ac_powerflow.jl => newton_ac_powerflow.jl} (89%) rename test/{test_nlsolve_powerflow.jl => test_newton_ac_powerflow.jl} (87%) diff --git a/Project.toml b/Project.toml index b339e3f3..708b8ecb 100644 --- a/Project.toml +++ b/Project.toml @@ -23,7 +23,7 @@ DataStructures = "0.18" Dates = "1" InfrastructureSystems = "2" JSON3 = "1" -KLU = "0.6.0" +KLU = "^0.6" LinearAlgebra = "1" Logging = "1" NLsolve = "4" diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 00c8c4a2..2fdfae9c 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -157,7 +157,7 @@ NOTE: use it for AC power flow computations. WARNING: functions for the evaluation of the multi-period AC PF still to be implemented. """ function PowerFlowData( - ::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + ::ACPowerFlow{<: ACPowerFlowSolverType}, sys::PSY.System; time_steps::Int = 1, timestep_names::Vector{String} = String[], @@ -440,7 +440,7 @@ Create an appropriate `PowerFlowContainer` for the given `PowerFlowEvaluationMod """ function make_power_flow_container end -make_power_flow_container(pfem::NLSolveACPowerFlow, sys::PSY.System; kwargs...) = +make_power_flow_container(pfem::ACPowerFlow{<: ACPowerFlowSolverType}, sys::PSY.System; kwargs...) = PowerFlowData(pfem, sys; kwargs...) make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) = diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index 499ac42c..e63802b2 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -6,6 +6,8 @@ export PowerFlowData export DCPowerFlow export NLSolveACPowerFlow export KLUACPowerFlow +export ACPowerFlow +export ACPowerFlowSolverType export PTDFDCPowerFlow export vPTDFDCPowerFlow export PSSEExportPowerFlow @@ -42,6 +44,6 @@ include("psse_export.jl") include("solve_dc_powerflow.jl") include("ac_power_flow.jl") include("ac_power_flow_jacobian.jl") -include("nlsolve_ac_powerflow.jl") +include("newton_ac_powerflow.jl") include("post_processing.jl") end diff --git a/src/definitions.jl b/src/definitions.jl index d8ecad1a..36e53b60 100644 --- a/src/definitions.jl +++ b/src/definitions.jl @@ -10,4 +10,7 @@ const DEFAULT_MAX_REDISTRIBUTION_ITERATIONS = 10 const ISAPPROX_ZERO_TOLERANCE = 1e-6 +const DEFAULT_NR_MAX_ITER::Int64 = 30 # default maxIter for the NR power flow +const DEFAULT_NR_TOL::Float64 = 1e-9 # default tolerance for the NR power flow + const AC_PF_KW = [] diff --git a/src/nlsolve_ac_powerflow.jl b/src/newton_ac_powerflow.jl similarity index 89% rename from src/nlsolve_ac_powerflow.jl rename to src/newton_ac_powerflow.jl index 6e08e67a..a909edb5 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -31,7 +31,7 @@ solve_ac_powerflow!(sys, method=:newton) ``` """ function solve_ac_powerflow!( - pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + pf::ACPowerFlow{<: ACPowerFlowSolverType}, system::PSY.System; kwargs..., ) @@ -41,7 +41,7 @@ function solve_ac_powerflow!( PSY.set_units_base_system!(system, "SYSTEM_BASE") check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false) data = PowerFlowData( - NLSolveACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits), + pf, system; check_connectivity = get(kwargs, :check_connectivity, true), ) @@ -72,7 +72,7 @@ res = solve_powerflow(sys, method=:newton) ``` """ function solve_powerflow( - pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + pf::ACPowerFlow{<: ACPowerFlowSolverType}, system::PSY.System; kwargs..., ) @@ -128,14 +128,14 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) end function _solve_powerflow!( - pf::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + pf::ACPowerFlow{<: ACPowerFlowSolverType}, data::ACPowerFlowData, check_reactive_power_limits; nlsolve_kwargs..., ) if check_reactive_power_limits for _ in 1:MAX_REACTIVE_POWER_ITERATIONS - converged, x = _nlsolve_powerflow(pf, data; nlsolve_kwargs...) + converged, x = _newton_powerflow(pf, data; nlsolve_kwargs...) if converged if _check_q_limit_bounds!(data, x) return converged, x @@ -145,12 +145,12 @@ function _solve_powerflow!( end end else - return _nlsolve_powerflow(pf, data; nlsolve_kwargs...) + return _newton_powerflow(pf, data; nlsolve_kwargs...) end end -function _nlsolve_powerflow( - pf::NLSolveACPowerFlow, +function _newton_powerflow( + pf::ACPowerFlow{NLSolveACPowerFlow}, data::ACPowerFlowData; nlsolve_kwargs..., ) @@ -160,29 +160,31 @@ function _nlsolve_powerflow( df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv) res = NLsolve.nlsolve(df, pf.x0; nlsolve_kwargs...) if !res.f_converged - @error("The powerflow solver returned convergence = $(res.f_converged)") + @error("The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))") end return res.f_converged, res.zero end -function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_kwargs...) - pf = PolarPowerFlow(data) - #J_function = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) - - maxIter = 30 # TODO - tol = 1e-6 # TODO +function _newton_powerflow(pf::ACPowerFlow{KLUACPowerFlow}, data::ACPowerFlowData; nlsolve_kwargs...) + # Fetch maxIter and tol from kwargs, or use defaults if not provided + maxIter = get(nlsolve_kwargs, :maxIter, DEFAULT_NR_MAX_ITER) + tol = get(nlsolve_kwargs, :tol, DEFAULT_NR_TOL) i = 0 - Vm = data.bus_magnitude[:] - Va = data.bus_angles[:] - V = Vm .* exp.(1im * Va) + pf = PolarPowerFlow(data) # 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) - Ybus = pf.data.power_network_matrix.data + Vm = data.bus_magnitude[:] + # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: + Vm[pq] = clamp.(Vm[pq], 0.9, 1.1) + Va = data.bus_angles[:] + V = Vm .* exp.(1im * Va) + + Ybus = data.power_network_matrix.data Sbus = data.bus_activepower_injection[:] - data.bus_activepower_withdrawals[:] + @@ -214,9 +216,6 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k factor_J = KLU.klu(J) dx = -(factor_J \ F) - #J_function(J_function.Jv, x) - #dx = - (J_function.Jv \ F) - Va[pv] .+= dx[1:npv] Va[pq] .+= dx[(npv + 1):(npv + npq)] Vm[pq] .+= dx[(npv + npq + 1):(npv + 2 * npq)] @@ -239,7 +238,7 @@ function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_k # 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)] + x = zeros(Float64, 2 * n_buses) Sbus_result = V .* conj(Ybus * V) for (ix, b) in enumerate(data.bus_type) if b == PSY.ACBusTypes.REF diff --git a/src/post_processing.jl b/src/post_processing.jl index 8916a378..8ae1d687 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( - ::Union{NLSolveACPowerFlow, KLUACPowerFlow}, + ::ACPowerFlow{<: ACPowerFlowSolverType}, sys::PSY.System, data::ACPowerFlowData, result::Vector{Float64}, diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl index 65b6b241..cc0c5c20 100644 --- a/src/powerflow_types.jl +++ b/src/powerflow_types.jl @@ -1,11 +1,17 @@ abstract type PowerFlowEvaluationModel end +abstract type ACPowerFlowSolverType end -Base.@kwdef struct NLSolveACPowerFlow <: PowerFlowEvaluationModel + +struct KLUACPowerFlow <: ACPowerFlowSolverType end +struct NLSolveACPowerFlow <: ACPowerFlowSolverType end + +Base.@kwdef struct ACPowerFlow{ACSolver <: ACPowerFlowSolverType} <: PowerFlowEvaluationModel check_reactive_power_limits::Bool = false end -Base.@kwdef struct KLUACPowerFlow <: PowerFlowEvaluationModel - check_reactive_power_limits::Bool = false +# Create a constructor that defaults to KLUACPowerFlow +function ACPowerFlow(; check_reactive_power_limits::Bool = false, ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow) + return ACPowerFlow{ACSolver}(check_reactive_power_limits) end struct DCPowerFlow <: PowerFlowEvaluationModel end diff --git a/test/test_nlsolve_powerflow.jl b/test/test_newton_ac_powerflow.jl similarity index 87% rename from test/test_nlsolve_powerflow.jl rename to test/test_newton_ac_powerflow.jl index 9d9d404f..7322b074 100644 --- a/test/test_nlsolve_powerflow.jl +++ b/test/test_newton_ac_powerflow.jl @@ -1,4 +1,4 @@ -@testset "AC Power Flow 14-Bus testing" for ACPowerFlow in +@testset "AC Power Flow 14-Bus testing" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) result_14 = [ 2.3255081760423684 @@ -32,36 +32,38 @@ ] sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) - data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true) + pf = ACPowerFlow{ACSolver}() + data = PowerFlows.PowerFlowData(pf, sys; check_connectivity = true) #Compare results between finite diff methods and Jacobian method - converged1, x1 = PowerFlows._solve_powerflow!(ACPowerFlow(), data, false) + converged1, x1 = PowerFlows._solve_powerflow!(pf, data, false) @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6 - @test solve_ac_powerflow!(ACPowerFlow(), sys; method = :newton) + @test solve_ac_powerflow!(pf, sys; method = :newton) # Test enforcing the reactive power Limits set_reactive_power!(get_component(PowerLoad, sys, "Bus4"), 0.0) - data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true) - converged2, x2 = PowerFlows._solve_powerflow!(ACPowerFlow(), data, true) + data = PowerFlows.PowerFlowData(pf, sys; check_connectivity = true) + converged2, x2 = PowerFlows._solve_powerflow!(pf, data, true) @test LinearAlgebra.norm(result_14 - x2, Inf) >= 1e-6 @test 1.08 <= x2[15] <= 1.09 end -@testset "AC Power Flow 14-Bus Line Configurations" for ACPowerFlow in +@testset "AC Power Flow 14-Bus Line Configurations" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) - base_res = solve_powerflow(ACPowerFlow(), sys) + pf = ACPowerFlow{ACSolver}() + base_res = solve_powerflow(pf, sys) branch = first(PSY.get_components(Line, sys)) dyn_branch = DynamicBranch(branch) add_component!(sys, dyn_branch) - @test dyn_pf = solve_ac_powerflow!(ACPowerFlow(), sys) - dyn_pf = solve_powerflow(ACPowerFlow(), sys) + @test dyn_pf = solve_ac_powerflow!(pf, sys) + dyn_pf = solve_powerflow(pf, 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!(ACPowerFlow(), sys) + solve_ac_powerflow!(pf, 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, rtol = 0) @@ -69,12 +71,12 @@ end 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(ACPowerFlow(), sys) + res = solve_powerflow(pf, sys) @test res["flow_results"].P_from_to[4] == 0.0 @test res["flow_results"].P_to_from[4] == 0.0 end -@testset "AC Power Flow 3-Bus Fixed FixedAdmittance testing" for ACPowerFlow in ( +@testset "AC Power Flow 3-Bus Fixed FixedAdmittance testing" for ACSolver in ( NLSolveACPowerFlow, KLUACPowerFlow, ) @@ -84,12 +86,13 @@ end 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(ACPowerFlow(), sys_3bus) + pf = ACPowerFlow{ACSolver}() + df = solve_powerflow(pf, 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 "AC Power Flow convergence fail testing" for ACPowerFlow in +@testset "AC Power Flow convergence fail testing" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) pf_sys5_re = PSB.build_system(PSB.PSITestSystems, "c_sys5_re"; add_forecasts = false) remove_component!(Line, pf_sys5_re, "1") @@ -98,15 +101,17 @@ end PSY.set_x!(br, 20.0) PSY.set_r!(br, 2.0) + pf = ACPowerFlow{ACSolver}() + # This is a negative test. The data passed for sys5_re is known to be infeasible. @test_logs( (:error, "The powerflow solver returned convergence = false"), match_mode = :any, - @test !solve_ac_powerflow!(ACPowerFlow(), pf_sys5_re) + @test !solve_ac_powerflow!(pf, pf_sys5_re) ) end -@testset "AC Test 240 Case PSS/e results" for ACPowerFlow in +@testset "AC Test 240 Case PSS/e results" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) file = joinpath( TEST_FILES_DIR, @@ -122,9 +127,11 @@ 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!(ACPowerFlow(), system) - @test pf - pf_result_df = solve_powerflow(ACPowerFlow(), system) + pf = ACPowerFlow{ACSolver}() + + pf1 = solve_ac_powerflow!(pf, system) + @test pf1 + pf_result_df = solve_powerflow(pf, 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) @@ -140,7 +147,7 @@ end @test norm(q_diff, 2) / length(q_diff) < DIFF_L2_TOLERANCE end -@testset "AC Multiple sources at ref" for ACPowerFlow in +@testset "AC Multiple sources at ref" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b = ACBus(; @@ -175,16 +182,17 @@ end X_th = 1e-5, ) add_component!(sys, s2) - @test solve_ac_powerflow!(ACPowerFlow(), sys) + pf = ACPowerFlow{ACSolver}() + @test solve_ac_powerflow!(pf, 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!(ACPowerFlow(), sys) + ) solve_ac_powerflow!(pf, sys) end -@testset "AC PowerFlow with Multiple sources at PV" for ACPowerFlow in +@testset "AC PowerFlow with Multiple sources at PV" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b1 = ACBus(; @@ -255,17 +263,19 @@ end ) add_component!(sys, s3) - @test solve_ac_powerflow!(ACPowerFlow(), sys) + pf = ACPowerFlow{ACSolver}() + + @test solve_ac_powerflow!(pf, 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!( - ACPowerFlow(), + pf, sys, ) end -@testset "AC PowerFlow Source + non-source at Ref" for ACPowerFlow in +@testset "AC PowerFlow Source + non-source at Ref" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b = ACBus(; @@ -312,7 +322,9 @@ end ) add_component!(sys, g1) - @test solve_ac_powerflow!(ACPowerFlow(), sys) + pf = ACPowerFlow{ACSolver}() + + @test solve_ac_powerflow!(pf, sys) @test isapprox( get_active_power(get_component(Source, sys, "source_1")), 0.5; @@ -325,7 +337,7 @@ end ) end -@testset "AC PowerFlow Source + non-source at PV" for ACPowerFlow in +@testset "AC PowerFlow Source + non-source at PV" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) sys = System(100.0) b1 = ACBus(; @@ -407,7 +419,9 @@ end ) add_component!(sys, g1) - @test solve_ac_powerflow!(ACPowerFlow(), sys) + pf = ACPowerFlow{ACSolver}() + + @test solve_ac_powerflow!(pf, sys) @test isapprox( get_active_power(get_component(Source, sys, "source_2")), 0.5; diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl index 5f91cddc..1f4152ad 100644 --- a/test/test_powerflow_data.jl +++ b/test/test_powerflow_data.jl @@ -1,6 +1,7 @@ @testset "PowerFlowData" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false) - @test PowerFlowData(NLSolveACPowerFlow(), sys) isa PF.ACPowerFlowData + @test PowerFlowData(ACPowerFlow{NLSolveACPowerFlow}(), sys) isa PF.ACPowerFlowData + @test PowerFlowData(ACPowerFlow{KLUACPowerFlow}(), sys) isa PF.ACPowerFlowData @test PowerFlowData(DCPowerFlow(), sys) isa PF.ABAPowerFlowData @test PowerFlowData(PTDFDCPowerFlow(), sys) isa PF.PTDFPowerFlowData @test PowerFlowData(vPTDFDCPowerFlow(), sys) isa PF.vPTDFPowerFlowData @@ -17,23 +18,23 @@ end PF.vPTDFPowerFlowData end -@testset "System <-> PowerFlowData round trip" begin +@testset "System <-> PowerFlowData round trip" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) # TODO currently only tested with ACPowerFlow # TODO test that update_system! errors if the PowerFlowData doesn't correspond to the system sys_original = build_system(PSISystems, "RTS_GMLC_DA_sys") - data_original = PowerFlowData(NLSolveACPowerFlow(), sys_original) + data_original = PowerFlowData(ACPowerFlow{ACSolver}(), sys_original) sys_modified = deepcopy(sys_original) modify_rts_system!(sys_modified) - data_modified = PowerFlowData(NLSolveACPowerFlow(), sys_original) + data_modified = PowerFlowData(ACPowerFlow{ACSolver}(), sys_original) modify_rts_powerflow!(data_modified) # update_system! with unmodified PowerFlowData should result in system that yields unmodified PowerFlowData # (NOTE does NOT necessarily yield original system due to power redistribution) sys_null_updated = deepcopy(sys_original) PF.update_system!(sys_null_updated, data_original) - data_null_updated = PowerFlowData(NLSolveACPowerFlow(), sys_null_updated) + data_null_updated = PowerFlowData(ACPowerFlow{ACSolver}(), sys_null_updated) @test IS.compare_values(powerflow_match_fn, data_null_updated, data_original) # Modified versions should not be the same as unmodified versions @@ -47,7 +48,7 @@ end # Constructing PowerFlowData from modified system should result in data_modified @test IS.compare_values( powerflow_match_fn, - PowerFlowData(NLSolveACPowerFlow(), sys_modified), + PowerFlowData(ACPowerFlow{ACSolver}(), sys_modified), data_modified, ) @@ -56,6 +57,6 @@ end sys_modify_updated = deepcopy(sys_original) PF.update_system!(sys_modify_updated, data_modified) sys_mod_redist = deepcopy(sys_modified) - PF.update_system!(sys_mod_redist, PowerFlowData(NLSolveACPowerFlow(), sys_mod_redist)) + PF.update_system!(sys_mod_redist, PowerFlowData(ACPowerFlow{ACSolver}(), sys_mod_redist)) @test IS.compare_values(powerflow_match_fn, sys_modify_updated, sys_mod_redist) end diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl index df7d925d..58fe743c 100644 --- a/test/test_psse_export.jl +++ b/test/test_psse_export.jl @@ -222,9 +222,9 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; return result end -function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow = false) - result1 = solve_powerflow(NLSolveACPowerFlow(), sys1) - result2 = solve_powerflow(NLSolveACPowerFlow(), sys2) +function test_power_flow(pf::ACPowerFlow{<:ACPowerFlowSolverType}, sys1::System, sys2::System; exclude_reactive_flow = false) + result1 = solve_powerflow(pf, sys1) + result2 = solve_powerflow(pf, sys2) reactive_power_tol = exclude_reactive_flow ? nothing : POWERFLOW_COMPARISON_TOLERANCE @test compare_df_within_tolerance("bus_results", result1["bus_results"], @@ -246,6 +246,7 @@ read_system_and_metadata(export_subdir) = read_system_and_metadata( get_psse_export_paths(export_subdir)...) function test_psse_round_trip( + pf::ACPowerFlow{<:ACPowerFlowSolverType}, sys::System, exporter::PSSEExporter, scenario_name::AbstractString, @@ -265,7 +266,7 @@ function test_psse_round_trip( sys2, sys2_metadata = read_system_and_metadata(raw_path, metadata_path) @test compare_systems_loosely(sys, sys2) do_power_flow_test && - test_power_flow(sys, sys2; exclude_reactive_flow = exclude_reactive_flow) + test_power_flow(pf, sys, sys2; exclude_reactive_flow = exclude_reactive_flow) end "Test that the two raw files are exactly identical and the two metadata files parse to identical JSON" @@ -318,17 +319,18 @@ end @test compare_systems_loosely(sys, deepcopy(sys)) end -@testset "PSSE Exporter with system_240[32].json, v33" begin +@testset "PSSE Exporter with system_240[32].json, v33" for (ACSolver, folder_name) in ((NLSolveACPowerFlow, "system_240_NLSolve"), (KLUACPowerFlow, "system_240_KLU")) sys = load_test_system() + pf = ACPowerFlow{ACSolver}() isnothing(sys) && return # PSS/E version must be one of the supported ones @test_throws ArgumentError PSSEExporter(sys, :vNonexistent, test_psse_export_dir) # Reimported export should be comparable to original system - export_location = joinpath(test_psse_export_dir, "v33", "system_240") + export_location = joinpath(test_psse_export_dir, "v33", folder_name) exporter = PSSEExporter(sys, :v33, export_location) - test_psse_round_trip(sys, exporter, "basic", export_location; + test_psse_round_trip(pf, sys, exporter, "basic", export_location; exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Exporting the exact same thing again should result in the exact same files @@ -360,20 +362,21 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys2)) - test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? + test_power_flow(pf, sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? end -@testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" begin +@testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" for (ACSolver, folder_name) in ((NLSolveACPowerFlow, "rts_gmlc_NLSolve"), (KLUACPowerFlow, "rts_gmlc_KLU")) sys = create_pf_friendly_rts_gmlc() + pf = ACPowerFlow{ACSolver}() set_units_base_system!(sys, UnitSystem.SYSTEM_BASE) # PSS/E version must be one of the supported ones @test_throws ArgumentError PSSEExporter(sys, :vNonexistent, test_psse_export_dir) # Reimported export should be comparable to original system - export_location = joinpath(test_psse_export_dir, "v33", "rts_gmlc") + export_location = joinpath(test_psse_export_dir, "v33", folder_name) exporter = PSSEExporter(sys, :v33, export_location) - test_psse_round_trip(sys, exporter, "basic", export_location; + test_psse_round_trip(pf, sys, exporter, "basic", export_location; exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Exporting the exact same thing again should result in the exact same files @@ -404,11 +407,11 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys2)) - test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? + test_power_flow(pf, sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Updating with changed value should result in a different reimport (PowerFlowData version) exporter = PSSEExporter(sys, :v33, export_location) - pf2 = PowerFlowData(NLSolveACPowerFlow(), sys) + pf2 = PowerFlowData(pf, sys) # This modifies the PowerFlowData in the same way that modify_rts_system! modifies the # system, so the reimport should be comparable to sys2 from above modify_rts_powerflow!(pf2) @@ -421,11 +424,11 @@ end @test_logs((:error, r"values do not match"), match_mode = :any, min_level = Logging.Error, compare_systems_loosely(sys, reread_sys3)) - test_power_flow(sys2, reread_sys3; exclude_reactive_flow = true) # TODO why is reactive flow not matching? + test_power_flow(pf, sys2, reread_sys3; exclude_reactive_flow = true) # TODO why is reactive flow not matching? # Exporting with write_comments should be comparable to original system exporter = PSSEExporter(sys, :v33, export_location; write_comments = true) - test_psse_round_trip(sys, exporter, "basic6", export_location; + test_psse_round_trip(pf, sys, exporter, "basic6", export_location; exclude_reactive_flow = true) # TODO why is reactive flow not matching? end From 42f02ab728c35fe8ab7a6e5f9a385999a82d5e96 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Mon, 16 Dec 2024 14:46:38 -0700 Subject: [PATCH 08/12] formatter --- src/PowerFlowData.jl | 8 ++++++-- src/newton_ac_powerflow.jl | 18 ++++++++++++------ src/post_processing.jl | 2 +- src/powerflow_types.jl | 9 ++++++--- test/test_powerflow_data.jl | 8 ++++++-- test/test_psse_export.jl | 17 ++++++++++++++--- 6 files changed, 45 insertions(+), 17 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 2fdfae9c..c22ca912 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -157,7 +157,7 @@ NOTE: use it for AC power flow computations. WARNING: functions for the evaluation of the multi-period AC PF still to be implemented. """ function PowerFlowData( - ::ACPowerFlow{<: ACPowerFlowSolverType}, + ::ACPowerFlow{<:ACPowerFlowSolverType}, sys::PSY.System; time_steps::Int = 1, timestep_names::Vector{String} = String[], @@ -440,7 +440,11 @@ Create an appropriate `PowerFlowContainer` for the given `PowerFlowEvaluationMod """ function make_power_flow_container end -make_power_flow_container(pfem::ACPowerFlow{<: ACPowerFlowSolverType}, sys::PSY.System; kwargs...) = +make_power_flow_container( + pfem::ACPowerFlow{<:ACPowerFlowSolverType}, + sys::PSY.System; + kwargs..., +) = PowerFlowData(pfem, sys; kwargs...) make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) = diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index a909edb5..e44d9aa6 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -31,7 +31,7 @@ solve_ac_powerflow!(sys, method=:newton) ``` """ function solve_ac_powerflow!( - pf::ACPowerFlow{<: ACPowerFlowSolverType}, + pf::ACPowerFlow{<:ACPowerFlowSolverType}, system::PSY.System; kwargs..., ) @@ -72,7 +72,7 @@ res = solve_powerflow(sys, method=:newton) ``` """ function solve_powerflow( - pf::ACPowerFlow{<: ACPowerFlowSolverType}, + pf::ACPowerFlow{<:ACPowerFlowSolverType}, system::PSY.System; kwargs..., ) @@ -128,7 +128,7 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) end function _solve_powerflow!( - pf::ACPowerFlow{<: ACPowerFlowSolverType}, + pf::ACPowerFlow{<:ACPowerFlowSolverType}, data::ACPowerFlowData, check_reactive_power_limits; nlsolve_kwargs..., @@ -160,12 +160,18 @@ function _newton_powerflow( df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv) res = NLsolve.nlsolve(df, pf.x0; nlsolve_kwargs...) if !res.f_converged - @error("The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))") + @error( + "The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))" + ) end return res.f_converged, res.zero end -function _newton_powerflow(pf::ACPowerFlow{KLUACPowerFlow}, data::ACPowerFlowData; nlsolve_kwargs...) +function _newton_powerflow( + pf::ACPowerFlow{KLUACPowerFlow}, + data::ACPowerFlowData; + nlsolve_kwargs..., +) # Fetch maxIter and tol from kwargs, or use defaults if not provided maxIter = get(nlsolve_kwargs, :maxIter, DEFAULT_NR_MAX_ITER) tol = get(nlsolve_kwargs, :tol, DEFAULT_NR_TOL) @@ -180,7 +186,7 @@ function _newton_powerflow(pf::ACPowerFlow{KLUACPowerFlow}, data::ACPowerFlowDat Vm = data.bus_magnitude[:] # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: - Vm[pq] = clamp.(Vm[pq], 0.9, 1.1) + Vm[pq] = clamp.(Vm[pq], 0.9, 1.1) Va = data.bus_angles[:] V = Vm .* exp.(1im * Va) diff --git a/src/post_processing.jl b/src/post_processing.jl index 8ae1d687..c8229201 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( - ::ACPowerFlow{<: ACPowerFlowSolverType}, + ::ACPowerFlow{<:ACPowerFlowSolverType}, sys::PSY.System, data::ACPowerFlowData, result::Vector{Float64}, diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl index cc0c5c20..fa9a1ae0 100644 --- a/src/powerflow_types.jl +++ b/src/powerflow_types.jl @@ -1,16 +1,19 @@ abstract type PowerFlowEvaluationModel end abstract type ACPowerFlowSolverType end - struct KLUACPowerFlow <: ACPowerFlowSolverType end struct NLSolveACPowerFlow <: ACPowerFlowSolverType end -Base.@kwdef struct ACPowerFlow{ACSolver <: ACPowerFlowSolverType} <: PowerFlowEvaluationModel +Base.@kwdef struct ACPowerFlow{ACSolver <: ACPowerFlowSolverType} <: + PowerFlowEvaluationModel check_reactive_power_limits::Bool = false end # Create a constructor that defaults to KLUACPowerFlow -function ACPowerFlow(; check_reactive_power_limits::Bool = false, ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow) +function ACPowerFlow(; + check_reactive_power_limits::Bool = false, + ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow, +) return ACPowerFlow{ACSolver}(check_reactive_power_limits) end diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl index 1f4152ad..2d45a148 100644 --- a/test/test_powerflow_data.jl +++ b/test/test_powerflow_data.jl @@ -18,7 +18,8 @@ end PF.vPTDFPowerFlowData end -@testset "System <-> PowerFlowData round trip" for ACSolver in (NLSolveACPowerFlow, KLUACPowerFlow) +@testset "System <-> PowerFlowData round trip" for ACSolver in + (NLSolveACPowerFlow, KLUACPowerFlow) # TODO currently only tested with ACPowerFlow # TODO test that update_system! errors if the PowerFlowData doesn't correspond to the system @@ -57,6 +58,9 @@ end sys_modify_updated = deepcopy(sys_original) PF.update_system!(sys_modify_updated, data_modified) sys_mod_redist = deepcopy(sys_modified) - PF.update_system!(sys_mod_redist, PowerFlowData(ACPowerFlow{ACSolver}(), sys_mod_redist)) + PF.update_system!( + sys_mod_redist, + PowerFlowData(ACPowerFlow{ACSolver}(), sys_mod_redist), + ) @test IS.compare_values(powerflow_match_fn, sys_modify_updated, sys_mod_redist) end diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl index 58fe743c..802eee3f 100644 --- a/test/test_psse_export.jl +++ b/test/test_psse_export.jl @@ -222,7 +222,12 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System; return result end -function test_power_flow(pf::ACPowerFlow{<:ACPowerFlowSolverType}, sys1::System, sys2::System; exclude_reactive_flow = false) +function test_power_flow( + pf::ACPowerFlow{<:ACPowerFlowSolverType}, + sys1::System, + sys2::System; + exclude_reactive_flow = false, +) result1 = solve_powerflow(pf, sys1) result2 = solve_powerflow(pf, sys2) reactive_power_tol = @@ -319,7 +324,10 @@ end @test compare_systems_loosely(sys, deepcopy(sys)) end -@testset "PSSE Exporter with system_240[32].json, v33" for (ACSolver, folder_name) in ((NLSolveACPowerFlow, "system_240_NLSolve"), (KLUACPowerFlow, "system_240_KLU")) +@testset "PSSE Exporter with system_240[32].json, v33" for (ACSolver, folder_name) in ( + (NLSolveACPowerFlow, "system_240_NLSolve"), + (KLUACPowerFlow, "system_240_KLU"), +) sys = load_test_system() pf = ACPowerFlow{ACSolver}() isnothing(sys) && return @@ -365,7 +373,10 @@ end test_power_flow(pf, sys2, reread_sys2; exclude_reactive_flow = true) # TODO why is reactive flow not matching? end -@testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" for (ACSolver, folder_name) in ((NLSolveACPowerFlow, "rts_gmlc_NLSolve"), (KLUACPowerFlow, "rts_gmlc_KLU")) +@testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" for (ACSolver, folder_name) in ( + (NLSolveACPowerFlow, "rts_gmlc_NLSolve"), + (KLUACPowerFlow, "rts_gmlc_KLU"), +) sys = create_pf_friendly_rts_gmlc() pf = ACPowerFlow{ACSolver}() set_units_base_system!(sys, UnitSystem.SYSTEM_BASE) From e6a0364de0ab375c5d9a535eab7777a1c5e79475 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Tue, 17 Dec 2024 14:14:58 -0700 Subject: [PATCH 09/12] small fix ACPowerFlow type; add test for results consistency for 2000 bus system --- src/powerflow_types.jl | 3 +-- test/test_newton_ac_powerflow.jl | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl index fa9a1ae0..e490e1f5 100644 --- a/src/powerflow_types.jl +++ b/src/powerflow_types.jl @@ -10,9 +10,8 @@ Base.@kwdef struct ACPowerFlow{ACSolver <: ACPowerFlowSolverType} <: end # Create a constructor that defaults to KLUACPowerFlow -function ACPowerFlow(; +function ACPowerFlow(ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow; check_reactive_power_limits::Bool = false, - ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow, ) return ACPowerFlow{ACSolver}(check_reactive_power_limits) end diff --git a/test/test_newton_ac_powerflow.jl b/test/test_newton_ac_powerflow.jl index 7322b074..7a0b31d8 100644 --- a/test/test_newton_ac_powerflow.jl +++ b/test/test_newton_ac_powerflow.jl @@ -433,3 +433,23 @@ end atol = 0.001, ) end + + +@testset "Compare larger grid results KLU vs NLSolve" begin + sys = build_system(MatpowerTestSystems, "matpower_ACTIVSg2000_sys") + + pf_default = ACPowerFlow() + pf_klu = ACPowerFlow(KLUACPowerFlow) + pf_nlsolve = ACPowerFlow(NLSolveACPowerFlow) + + res_default = solve_powerflow(pf_default, sys) # must be the same as KLU + res_klu = solve_powerflow(pf_klu, sys) + res_nlsolve = solve_powerflow(pf_nlsolve, sys) + + + @test all(isapprox.(res_klu["bus_results"][!, :Vm], res_default["bus_results"][!, :Vm], rtol=0, atol=1e-12)) + @test all(isapprox.(res_klu["bus_results"][!, :θ], res_default["bus_results"][!, :θ], rtol=0, atol=1e-12)) + + @test all(isapprox.(res_klu["bus_results"][!, :Vm], res_nlsolve["bus_results"][!, :Vm], rtol=0, atol=1e-8)) + @test all(isapprox.(res_klu["bus_results"][!, :θ], res_nlsolve["bus_results"][!, :θ], rtol=0, atol=1e-8)) +end \ No newline at end of file From da200a12b3af28a7cfa81f607841f78efcf66eaa Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Tue, 17 Dec 2024 14:18:59 -0700 Subject: [PATCH 10/12] small change --- src/newton_ac_powerflow.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index e44d9aa6..2d428553 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -177,8 +177,6 @@ function _newton_powerflow( tol = get(nlsolve_kwargs, :tol, DEFAULT_NR_TOL) i = 0 - pf = PolarPowerFlow(data) - # 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) From b385c6530d736a85e36e4f495507a2b78635df86 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 19 Dec 2024 10:51:38 -0700 Subject: [PATCH 11/12] optimize NR code to reduce memory allocations --- src/PowerFlows.jl | 2 +- src/newton_ac_powerflow.jl | 128 ++++++++++++++++++++++++++++--------- 2 files changed, 99 insertions(+), 31 deletions(-) diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index e63802b2..00907a8b 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -27,7 +27,7 @@ import KLU import SparseArrays import InfrastructureSystems import PowerNetworkMatrices -import SparseArrays: SparseMatrixCSC +import SparseArrays: SparseMatrixCSC, sparse import JSON3 import DataStructures: OrderedDict import Dates diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 2d428553..4ef55631 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -181,12 +181,31 @@ function _newton_powerflow( 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) + pvpq = [pv; pq] + + #nref = length(ref) + npv = length(pv) + npq = length(pq) + npvpq = npv + npq + n_buses = length(data.bus_type) Vm = data.bus_magnitude[:] # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: - Vm[pq] = clamp.(Vm[pq], 0.9, 1.1) + Vm[pq] .= clamp.(Vm[pq], 0.9, 1.1) Va = data.bus_angles[:] - V = Vm .* exp.(1im * Va) + V = zeros(Complex{Float64}, length(Vm)) + V .= Vm .* exp.(1im * Va) + + Va_pv = view(Va, pv) + Va_pq = view(Va, pq) + Vm_pq = view(Vm, pq) + + # pre-allocate dx + dx = zeros(Float64, npv + 2 * npq) + + dx_Va_pv = view(dx, 1:npv) + dx_Va_pq = view(dx, (npv + 1):(npv + npq)) + dx_Vm_pq = view(dx, (npv + npq + 1):(npv + 2 * npq)) Ybus = data.power_network_matrix.data @@ -194,43 +213,92 @@ function _newton_powerflow( 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])] - - # nref = length(ref) - npv = length(pv) - npq = length(pq) - - converged = (npv + npq) == 0 # if only ref buses present, we do not need to enter the loop + # Pre-allocate mis and F and create views for the respective real and imaginary sections of the arrays: + mis = zeros(Complex{Float64}, length(V)) + mis_pvpq = view(mis, pvpq) + mis_pq = view(mis, pq) + + F = zeros(Float64, npvpq + npq) + F_real = view(F, 1:npvpq) + F_imag = view(F, npvpq + 1:npvpq + npq) + + mis .= V .* conj(Ybus * V) .- Sbus + F_real .= real(mis_pvpq) # In-place assignment to the real part, using views + F_imag .= imag(mis_pq) # In-place assignment to the imaginary part, using views + + converged = npvpq == 0 # if only ref buses present, we do not need to enter the loop + + # preallocate Jacobian matrix + rows = vcat(1:npvpq, 1:npvpq, npvpq+1:npvpq+npq, npvpq+1:npvpq+npq) + cols = vcat(1:npvpq, npvpq+1:npvpq+npq, 1:npvpq, npvpq+1:npvpq+npq) + J = sparse(rows, cols, Float64(0)) + + # we need to define lookups for mappings of pv, pq buses onto the internal J indexing + pvpq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) + pvpq_lookup[pvpq] .= 1:npvpq + pq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1) + pq_lookup[pq] .= 1:npq + + # with the pre-allocated J and lookups, we can define views into the sub-matrices of the J matrix for updating the J matrix in the NR loop + j11 = view(J, pvpq_lookup[pvpq], pvpq_lookup[pvpq]) + j12 = view(J, pvpq_lookup[pvpq], npvpq .+ pq_lookup[pq]) + j21 = view(J, npvpq .+ pq_lookup[pq], pvpq_lookup[pvpq]) + j22 = view(J, npvpq .+ pq_lookup[pq], npvpq .+ pq_lookup[pq]) + + # pre-allocate the dSbus_dVm, dSbus_dVa to have the same structure as Ybus + # they will follow the structure of Ybus except maybe when Ybus has zero values in its diagonal, which we do not expect here + rows, cols, _ = SparseArrays.findnz(Ybus) + dSbus_dVm = sparse(rows, cols, Complex{Float64}(0)) + dSbus_dVa = sparse(rows, cols, Complex{Float64}(0)) + + # create views for the sub-arrays of Sbus_dVa, Sbus_dVm for updating the J: + Sbus_dVa_j11 = view(dSbus_dVa, pvpq, pvpq) + Sbus_dVm_j12 = view(dSbus_dVm, pvpq, pq) + Sbus_dVa_j21 = view(dSbus_dVa, pq, pvpq) + Sbus_dVm_j22 = view(dSbus_dVm, pq, pq) + + # we need views of the diagonals to avoid using LinearAlgebra.Diagonal: + diagV = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) + diag_idx = LinearAlgebra.diagind(diagV) + diagV_diag = view(diagV, diag_idx) + + diagIbus = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) + diagIbus_diag = view(diagIbus, diag_idx) + + diagVnorm = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) + diagVnorm_diag = view(diagVnorm, diag_idx) while i < maxIter && !converged i += 1 - 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) - - j11 = real(dSbus_dVa[[pv; pq], [pv; pq]]) - j12 = real(dSbus_dVm[[pv; pq], pq]) - j21 = imag(dSbus_dVa[pq, [pv; pq]]) - j22 = imag(dSbus_dVm[pq, pq]) - J = [j11 j12; j21 j22] + + ## use the new value of V to update dSbus_dVa, dSbus_dVm: + diagV_diag .= V + diagIbus_diag .= Ybus * V + diagVnorm_diag .= V ./ abs.(V) + dSbus_dVm .= diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm + dSbus_dVa .= 1im * diagV * conj(diagIbus - Ybus * diagV) + + # update the Jacobian by setting values through the pre-defined views for j11, j12, j21, j22 + j11 .= real(Sbus_dVa_j11) + j12 .= real(Sbus_dVm_j12) + j21 .= imag(Sbus_dVa_j21) + j22 .= imag(Sbus_dVm_j22) factor_J = KLU.klu(J) - dx = -(factor_J \ F) + dx .= -(factor_J \ F) - Va[pv] .+= dx[1:npv] - Va[pq] .+= dx[(npv + 1):(npv + npq)] - Vm[pq] .+= dx[(npv + npq + 1):(npv + 2 * npq)] + Va_pv .+= dx_Va_pv + Va_pq .+= dx_Va_pq + Vm_pq .+= dx_Vm_pq - V = Vm .* exp.(1im * Va) + V .= Vm .* exp.(1im * Va) - Vm = abs.(V) - Va = angle.(V) + Vm .= abs.(V) + Va .= angle.(V) - mis = V .* conj(Ybus * V) - Sbus - F = [real(mis[[pv; pq]]); imag(mis[pq])] + mis .= V .* conj(Ybus * V) .- Sbus + F_real .= real(mis_pvpq) # In-place assignment to the real part + F_imag .= imag(mis_pq) # In-place assignment to the imaginary part converged = LinearAlgebra.norm(F, Inf) < tol end From d7f7c88aaa20a516976c6f82e71a012b00dbbc39 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 19 Dec 2024 13:34:11 -0700 Subject: [PATCH 12/12] trying to reduce memory allocation --- src/newton_ac_powerflow.jl | 64 ++++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 30 deletions(-) diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 4ef55631..014384b3 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -191,10 +191,10 @@ function _newton_powerflow( Vm = data.bus_magnitude[:] # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: - Vm[pq] .= clamp.(Vm[pq], 0.9, 1.1) + @. Vm[pq] = clamp.(Vm[pq], 0.9, 1.1) Va = data.bus_angles[:] V = zeros(Complex{Float64}, length(Vm)) - V .= Vm .* exp.(1im * Va) + @. V = Vm .* exp.(1im * Va) Va_pv = view(Va, pv) Va_pq = view(Va, pq) @@ -212,7 +212,7 @@ function _newton_powerflow( Sbus = data.bus_activepower_injection[:] - data.bus_activepower_withdrawals[:] + 1im * (data.bus_reactivepower_injection[:] - data.bus_reactivepower_withdrawals[:]) - + # Pre-allocate mis and F and create views for the respective real and imaginary sections of the arrays: mis = zeros(Complex{Float64}, length(V)) mis_pvpq = view(mis, pvpq) @@ -222,9 +222,9 @@ function _newton_powerflow( F_real = view(F, 1:npvpq) F_imag = view(F, npvpq + 1:npvpq + npq) - mis .= V .* conj(Ybus * V) .- Sbus - F_real .= real(mis_pvpq) # In-place assignment to the real part, using views - F_imag .= imag(mis_pq) # In-place assignment to the imaginary part, using views + mis .= V .* conj.(Ybus * V) .- Sbus + @. F_real = real(mis_pvpq) # In-place assignment to the real part, using views + @. F_imag = imag(mis_pq) # In-place assignment to the imaginary part, using views converged = npvpq == 0 # if only ref buses present, we do not need to enter the loop @@ -244,18 +244,6 @@ function _newton_powerflow( j12 = view(J, pvpq_lookup[pvpq], npvpq .+ pq_lookup[pq]) j21 = view(J, npvpq .+ pq_lookup[pq], pvpq_lookup[pvpq]) j22 = view(J, npvpq .+ pq_lookup[pq], npvpq .+ pq_lookup[pq]) - - # pre-allocate the dSbus_dVm, dSbus_dVa to have the same structure as Ybus - # they will follow the structure of Ybus except maybe when Ybus has zero values in its diagonal, which we do not expect here - rows, cols, _ = SparseArrays.findnz(Ybus) - dSbus_dVm = sparse(rows, cols, Complex{Float64}(0)) - dSbus_dVa = sparse(rows, cols, Complex{Float64}(0)) - - # create views for the sub-arrays of Sbus_dVa, Sbus_dVm for updating the J: - Sbus_dVa_j11 = view(dSbus_dVa, pvpq, pvpq) - Sbus_dVm_j12 = view(dSbus_dVm, pvpq, pq) - Sbus_dVa_j21 = view(dSbus_dVa, pq, pvpq) - Sbus_dVm_j22 = view(dSbus_dVm, pq, pq) # we need views of the diagonals to avoid using LinearAlgebra.Diagonal: diagV = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) @@ -268,21 +256,37 @@ function _newton_powerflow( diagVnorm = sparse(1:n_buses, 1:n_buses, Complex{Float64}(1)) diagVnorm_diag = view(diagVnorm, diag_idx) + # pre-allocate the dSbus_dVm, dSbus_dVa to have the same structure as Ybus + # they will follow the structure of Ybus except maybe when Ybus has zero values in its diagonal, which we do not expect here + #rows, cols, _ = SparseArrays.findnz(Ybus) + #dSbus_dVm = sparse(rows, cols, Complex{Float64}(0)) + #dSbus_dVa = sparse(rows, cols, Complex{Float64}(0)) + + # preallocate dSbus_dVm, dSbus_dVa with correct structure: + dSbus_dVm = diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm + dSbus_dVa = 1im * diagV * conj.(diagIbus - Ybus * diagV) + + # create views for the sub-arrays of Sbus_dVa, Sbus_dVm for updating the J: + Sbus_dVa_j11 = view(dSbus_dVa, pvpq, pvpq) + Sbus_dVm_j12 = view(dSbus_dVm, pvpq, pq) + Sbus_dVa_j21 = view(dSbus_dVa, pq, pvpq) + Sbus_dVm_j22 = view(dSbus_dVm, pq, pq) + while i < maxIter && !converged i += 1 ## use the new value of V to update dSbus_dVa, dSbus_dVm: diagV_diag .= V diagIbus_diag .= Ybus * V - diagVnorm_diag .= V ./ abs.(V) - dSbus_dVm .= diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm - dSbus_dVa .= 1im * diagV * conj(diagIbus - Ybus * diagV) + @. diagVnorm_diag = V ./ abs.(V) + dSbus_dVm .= diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm + dSbus_dVa .= 1im * diagV * conj.(diagIbus - Ybus * diagV) # update the Jacobian by setting values through the pre-defined views for j11, j12, j21, j22 - j11 .= real(Sbus_dVa_j11) - j12 .= real(Sbus_dVm_j12) - j21 .= imag(Sbus_dVa_j21) - j22 .= imag(Sbus_dVm_j22) + @. j11 = real(Sbus_dVa_j11) + @. j12 = real(Sbus_dVm_j12) + @. j21 = imag(Sbus_dVa_j21) + @. j22 = imag(Sbus_dVm_j22) factor_J = KLU.klu(J) dx .= -(factor_J \ F) @@ -291,21 +295,21 @@ function _newton_powerflow( Va_pq .+= dx_Va_pq Vm_pq .+= dx_Vm_pq - V .= Vm .* exp.(1im * Va) + @. V = Vm .* exp.(1im * Va) Vm .= abs.(V) Va .= angle.(V) - mis .= V .* conj(Ybus * V) .- Sbus - F_real .= real(mis_pvpq) # In-place assignment to the real part - F_imag .= imag(mis_pq) # In-place assignment to the imaginary part + mis .= V .* conj.(Ybus * V) .- Sbus + @. F_real = real(mis_pvpq) # In-place assignment to the real part + @. F_imag = imag(mis_pq) # In-place assignment to the imaginary part converged = LinearAlgebra.norm(F, Inf) < tol end if !converged @error("The powerflow solver with KLU did not converge after $i iterations") else - @debug("The powerflow solver with KLU converged after $i iterations") + @info("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: