From 2f90d829c49a13a37e023d61c82a03e6b0fa822b Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 16 Jan 2025 15:15:22 -0700 Subject: [PATCH 1/4] add field "converged" to PowerFlowData, return NaN when not converged --- src/PowerFlowData.jl | 11 +++++++- src/common.jl | 4 +++ src/newton_ac_powerflow.jl | 53 ++++++++++++++++++++---------------- src/nlsolve_powerflow.jl | 7 +++-- src/solve_dc_powerflow.jl | 3 ++ test/test_utils/legacy_pf.jl | 4 ++- 6 files changed, 54 insertions(+), 28 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index af1b5bc8..fd144b7b 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -102,6 +102,7 @@ struct PowerFlowData{ power_network_matrix::M aux_network_matrix::N neighbors::Vector{Set{Int}} + converged::Vector{Bool} end get_bus_lookup(pfd::PowerFlowData) = pfd.bus_lookup @@ -129,6 +130,7 @@ get_power_network_matrix(pfd::PowerFlowData) = pfd.power_network_matrix get_aux_network_matrix(pfd::PowerFlowData) = pfd.aux_network_matrix get_neighbor(pfd::PowerFlowData) = pfd.neighbors supports_multi_period(::PowerFlowData) = true +get_converged(pfd::PowerFlowData) = pfd.converged function clear_injection_data!(pfd::PowerFlowData) pfd.bus_activepower_injection .= 0.0 @@ -220,10 +222,10 @@ function PowerFlowData( branch_types[ix] = typeof(b) end - timestep_map = Dict(1 => "1") valid_ix = setdiff(1:n_buses, ref_bus_positions) neighbors = _calculate_neighbors(power_network_matrix) aux_network_matrix = nothing + converged = fill(false, time_steps) return make_powerflowdata( sys, @@ -239,6 +241,7 @@ function PowerFlowData( timestep_map, valid_ix, neighbors, + converged, ) end @@ -295,6 +298,7 @@ function PowerFlowData( 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) + converged = fill(false, time_steps) return make_dc_powerflowdata( sys, time_steps, @@ -307,6 +311,7 @@ function PowerFlowData( branch_lookup, temp_bus_map, valid_ix, + converged, ) end @@ -364,6 +369,7 @@ function PowerFlowData( 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) + converged = fill(false, time_steps) return make_dc_powerflowdata( sys, time_steps, @@ -376,6 +382,7 @@ function PowerFlowData( branch_lookup, temp_bus_map, valid_ix, + converged, ) end @@ -432,6 +439,7 @@ function PowerFlowData( 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) + converged = fill(false, time_steps) return make_dc_powerflowdata( sys, time_steps, @@ -444,6 +452,7 @@ function PowerFlowData( branch_lookup, temp_bus_map, valid_ix, + converged, ) end diff --git a/src/common.jl b/src/common.jl index cd2d1fb3..0eb32bf9 100644 --- a/src/common.jl +++ b/src/common.jl @@ -117,6 +117,7 @@ function make_dc_powerflowdata( branch_lookup, temp_bus_map, valid_ix, + converged, ) branch_type = Vector{DataType}(undef, length(branch_lookup)) for (ix, b) in enumerate(PNM.get_ac_branches(sys)) @@ -139,6 +140,7 @@ function make_dc_powerflowdata( timestep_map, valid_ix, neighbors, + converged, ) end @@ -156,6 +158,7 @@ function make_powerflowdata( timestep_map, valid_ix, neighbors, + converged, ) bus_type = Vector{PSY.ACBusTypes}(undef, n_buses) bus_angles = zeros(Float64, n_buses) @@ -249,5 +252,6 @@ function make_powerflowdata( power_network_matrix, aux_network_matrix, neighbors, + converged, ) end diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 570732f8..3b10675c 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -177,10 +177,13 @@ function solve_powerflow!( sorted_time_steps = sort(collect(keys(data.timestep_map))) # preallocate results - ts_converged = zeros(Bool, 1, length(sorted_time_steps)) + ts_converged = fill(false, length(sorted_time_steps)) ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) + # TODO If anything in the grid topology changes, + # e.g. tap positions of transformers or in service + # status of branches, Yft and Ytf must be updated! Yft = data.power_network_matrix.data_ft Ytf = data.power_network_matrix.data_tf @@ -190,25 +193,25 @@ function solve_powerflow!( for t in sorted_time_steps converged, V, Sbus_result = _ac_powereflow(data, pf; time_step = t, kwargs...) - ts_converged[1, t] = converged + ts_converged[t] = converged ts_V[:, t] .= V ts_S[:, t] .= Sbus_result - ref = findall( - x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, - data.bus_type[:, t], - ) - pv = findall( - x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, - data.bus_type[:, t], - ) - pq = findall( - x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, - data.bus_type[:, t], - ) - - # temporary implementation that will need to be improved: if converged + ref = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, + data.bus_type[:, t], + ) + pv = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, + data.bus_type[:, t], + ) + pq = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, + data.bus_type[:, t], + ) + + # temporary implementation that will need to be improved: # write results for REF data.bus_activepower_injection[ref, t] .= real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t] @@ -219,21 +222,22 @@ function solve_powerflow!( imag.(Sbus_result[pv]) .+ data.bus_reactivepower_withdrawals[pv, t] # results for PQ buses do not need to be updated -> already consistent with inputs - # write bus bus_types - # todo - # write voltage results data.bus_magnitude[pq, t] .= abs.(V[pq]) data.bus_angles[pq, t] .= angle.(V[pq]) data.bus_angles[pv, t] .= angle.(V[pv]) - else - # todo - 1 + 2 + data.bus_activepower_injection[:, t] .= NaN64 + data.bus_activepower_withdrawals[:, t] .= NaN64 + data.bus_reactivepower_injection[:, t] .= NaN64 + data.bus_reactivepower_withdrawals[:, t] .= NaN64 + data.bus_magnitude[:, t] .= NaN64 + data.bus_angles[:, t] .= NaN64 end end # write branch flows + # TODO if Yft, Ytf change between time steps, this must be moved inside the loop! Sft = ts_V[fb, :] .* conj.(Yft * ts_V) Stf = ts_V[tb, :] .* conj.(Ytf * ts_V) @@ -242,8 +246,7 @@ function solve_powerflow!( data.branch_activepower_flow_to_from .= real.(Stf) data.branch_reactivepower_flow_to_from .= imag.(Stf) - # todo: - # return df_results + data.converged .= ts_converged return end @@ -713,6 +716,8 @@ function _newton_powerflow( end if !converged + V .*= NaN64 + Sbus_result .*= NaN64 @error("The powerflow solver with KLU did not converge after $i iterations") else @info("The powerflow solver with KLU converged after $i iterations") diff --git a/src/nlsolve_powerflow.jl b/src/nlsolve_powerflow.jl index c04e5fca..7b765d13 100644 --- a/src/nlsolve_powerflow.jl +++ b/src/nlsolve_powerflow.jl @@ -23,12 +23,15 @@ function _newton_powerflow( df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv) res = NLsolve.nlsolve(df, pf.x0; nlsolve_solver_kwargs...) if !res.f_converged + V = fill(NaN64 + NaN64 * im, length(res.zero) ÷ 2) + Sbus_result = fill(NaN64 + NaN64 * im, length(res.zero) ÷ 2) @error( "The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))" ) + else + V = _calc_V(data, res.zero; time_step = time_step) + Sbus_result = V .* conj(data.power_network_matrix.data * V) end - V = _calc_V(data, res.zero; time_step = time_step) - Sbus_result = V .* conj(data.power_network_matrix.data * V) return (res.f_converged, V, Sbus_result) end diff --git a/src/solve_dc_powerflow.jl b/src/solve_dc_powerflow.jl index e911d329..b75ca9f2 100644 --- a/src/solve_dc_powerflow.jl +++ b/src/solve_dc_powerflow.jl @@ -53,6 +53,7 @@ function solve_powerflow!( # evaluate bus angles p_inj = power_injection[data.valid_ix, :] data.bus_angles[data.valid_ix, :] .= data.aux_network_matrix.K \ p_inj + data.converged .= true return end @@ -77,6 +78,7 @@ function solve_powerflow!( p_inj = power_injection[data.valid_ix, i] data.bus_angles[data.valid_ix, i] .= data.aux_network_matrix.K \ p_inj end + data.converged .= true return end @@ -100,6 +102,7 @@ function solve_powerflow!( data.power_network_matrix.K \ @view power_injection[data.valid_ix, :] data.branch_activepower_flow_from_to .= data.aux_network_matrix.data' * data.bus_angles data.branch_activepower_flow_to_from .= -data.branch_activepower_flow_from_to + data.converged .= true return end diff --git a/test/test_utils/legacy_pf.jl b/test/test_utils/legacy_pf.jl index 5436b334..5f1c0a50 100644 --- a/test/test_utils/legacy_pf.jl +++ b/test/test_utils/legacy_pf.jl @@ -111,10 +111,12 @@ function _newton_powerflow( end if !converged + V .*= NaN64 + Sbus_result = fill(NaN64 + NaN64 * im, length(V)) @error("The powerflow solver with KLU did not converge after $i iterations") else + Sbus_result = V .* conj(Ybus * V) @info("The powerflow solver with KLU converged after $i iterations") end - Sbus_result = V .* conj(Ybus * V) return (converged, V, Sbus_result) end From 685a6b0fa202f5383486fe20d82600e5c48bedd4 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 16 Jan 2025 17:09:56 -0700 Subject: [PATCH 2/4] wip: SolverData, penalty factors --- src/PowerFlowData.jl | 15 +++++++++++++++ src/common.jl | 11 +++++++++++ src/newton_ac_powerflow.jl | 7 ++++++- test/test_utils/legacy_pf.jl | 3 +++ 4 files changed, 35 insertions(+), 1 deletion(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index fd144b7b..021b99d8 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -1,3 +1,8 @@ +Base.@kwdef mutable struct SolverData + J::Union{SparseMatrixCSC{Float64, Int}, Nothing} = nothing + dSbus_dV_ref::Union{Vector{Float64}, Nothing} = nothing +end + abstract type PowerFlowContainer end """ @@ -103,6 +108,7 @@ struct PowerFlowData{ aux_network_matrix::N neighbors::Vector{Set{Int}} converged::Vector{Bool} + solver_data::Vector{SolverData} end get_bus_lookup(pfd::PowerFlowData) = pfd.bus_lookup @@ -131,6 +137,7 @@ get_aux_network_matrix(pfd::PowerFlowData) = pfd.aux_network_matrix get_neighbor(pfd::PowerFlowData) = pfd.neighbors supports_multi_period(::PowerFlowData) = true get_converged(pfd::PowerFlowData) = pfd.converged +get_solver_data(pfd::PowerFlowData) = pfd.solver_data function clear_injection_data!(pfd::PowerFlowData) pfd.bus_activepower_injection .= 0.0 @@ -226,6 +233,7 @@ function PowerFlowData( neighbors = _calculate_neighbors(power_network_matrix) aux_network_matrix = nothing converged = fill(false, time_steps) + solver_data = [SolverData() for _ in 1:time_steps] return make_powerflowdata( sys, @@ -242,6 +250,7 @@ function PowerFlowData( valid_ix, neighbors, converged, + solver_data, ) end @@ -299,6 +308,7 @@ function PowerFlowData( ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) converged = fill(false, time_steps) + solver_data = [SolverData() for _ in 1:time_steps] return make_dc_powerflowdata( sys, time_steps, @@ -312,6 +322,7 @@ function PowerFlowData( temp_bus_map, valid_ix, converged, + solver_data, ) end @@ -370,6 +381,7 @@ function PowerFlowData( ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) converged = fill(false, time_steps) + solver_data = [SolverData() for _ in 1:time_steps] return make_dc_powerflowdata( sys, time_steps, @@ -383,6 +395,7 @@ function PowerFlowData( temp_bus_map, valid_ix, converged, + solver_data, ) end @@ -440,6 +453,7 @@ function PowerFlowData( ) valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions) converged = fill(false, time_steps) + solver_data = [SolverData() for _ in 1:time_steps] return make_dc_powerflowdata( sys, time_steps, @@ -453,6 +467,7 @@ function PowerFlowData( temp_bus_map, valid_ix, converged, + solver_data, ) end diff --git a/src/common.jl b/src/common.jl index a05dc7a3..a0bae85c 100644 --- a/src/common.jl +++ b/src/common.jl @@ -150,6 +150,7 @@ function make_dc_powerflowdata( temp_bus_map, valid_ix, converged, + solver_data, ) branch_type = Vector{DataType}(undef, length(branch_lookup)) for (ix, b) in enumerate(PNM.get_ac_branches(sys)) @@ -173,6 +174,7 @@ function make_dc_powerflowdata( valid_ix, neighbors, converged, + solver_data, ) end @@ -191,6 +193,7 @@ function make_powerflowdata( valid_ix, neighbors, converged, + solver_data, ) bus_type = Vector{PSY.ACBusTypes}(undef, n_buses) bus_angles = zeros(Float64, n_buses) @@ -285,5 +288,13 @@ function make_powerflowdata( aux_network_matrix, neighbors, converged, + solver_data, ) end + +# work in progress - quick but not optimized function for POC +function penalty_factors(solver_data::SolverData) + J_t = transpose(solver_data.J) + f = J_t \ solver_data.dSbus_dV_ref + return f +end diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 3b10675c..49d968e2 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -175,7 +175,7 @@ function solve_powerflow!( ) pf = ACPowerFlow() # todo: somehow store in data which PF to use (see issue #50) - sorted_time_steps = sort(collect(keys(data.timestep_map))) + sorted_time_steps = get(kwargs, :time_steps, sort(collect(keys(data.timestep_map)))) # preallocate results ts_converged = fill(false, length(sorted_time_steps)) ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) @@ -540,6 +540,8 @@ function _newton_powerflow( tol = get(kwargs, :tol, DEFAULT_NR_TOL) i = 0 + solver_data = data.solver_data[time_step] + Ybus = data.power_network_matrix.data # Find indices for each bus type @@ -720,6 +722,9 @@ function _newton_powerflow( Sbus_result .*= NaN64 @error("The powerflow solver with KLU did not converge after $i iterations") else + solver_data.J = J + solver_data.dSbus_dV_ref = + [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pvpq]))] @info("The powerflow solver with KLU converged after $i iterations") end diff --git a/test/test_utils/legacy_pf.jl b/test/test_utils/legacy_pf.jl index 5f1c0a50..c7dabb4d 100644 --- a/test/test_utils/legacy_pf.jl +++ b/test/test_utils/legacy_pf.jl @@ -116,6 +116,9 @@ function _newton_powerflow( @error("The powerflow solver with KLU did not converge after $i iterations") else Sbus_result = V .* conj(Ybus * V) + solver_data.J = J + solver_data.dSbus_dV_ref = + [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pvpq]))] @info("The powerflow solver with KLU converged after $i iterations") end return (converged, V, Sbus_result) From d2fb7b61ab2bd4ef0964880ce9d3159ce664b0f1 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 16 Jan 2025 17:23:02 -0700 Subject: [PATCH 3/4] fix penalty factors function --- src/PowerFlows.jl | 1 + src/common.jl | 7 ------- src/post_processing.jl | 7 +++++++ 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl index 78cbeeb8..13672ffe 100644 --- a/src/PowerFlows.jl +++ b/src/PowerFlows.jl @@ -16,6 +16,7 @@ export PSSEExporter export update_exporter! export write_export export get_psse_export_paths +export penalty_factors import Logging import DataFrames diff --git a/src/common.jl b/src/common.jl index a0bae85c..0343c053 100644 --- a/src/common.jl +++ b/src/common.jl @@ -291,10 +291,3 @@ function make_powerflowdata( solver_data, ) end - -# work in progress - quick but not optimized function for POC -function penalty_factors(solver_data::SolverData) - J_t = transpose(solver_data.J) - f = J_t \ solver_data.dSbus_dV_ref - return f -end diff --git a/src/post_processing.jl b/src/post_processing.jl index 294d8df8..e3b82fa8 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -757,3 +757,10 @@ function set_system_time_step(sys::PSY.System, data::PowerFlowData; time_step = end end end + +# work in progress - quick but not optimized function for POC +function penalty_factors(solver_data::SolverData) + J_t = transpose(solver_data.J) + f = J_t \ solver_data.dSbus_dV_ref + return f +end From 2106c0cc4ce952c97f2f566cef95713770919f94 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 16 Jan 2025 17:25:09 -0700 Subject: [PATCH 4/4] bugfix shape penalty factors --- src/newton_ac_powerflow.jl | 2 +- test/test_utils/legacy_pf.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 49d968e2..cd41f717 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -724,7 +724,7 @@ function _newton_powerflow( else solver_data.J = J solver_data.dSbus_dV_ref = - [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pvpq]))] + [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pq]))] @info("The powerflow solver with KLU converged after $i iterations") end diff --git a/test/test_utils/legacy_pf.jl b/test/test_utils/legacy_pf.jl index c7dabb4d..84865b3d 100644 --- a/test/test_utils/legacy_pf.jl +++ b/test/test_utils/legacy_pf.jl @@ -118,7 +118,7 @@ function _newton_powerflow( Sbus_result = V .* conj(Ybus * V) solver_data.J = J solver_data.dSbus_dV_ref = - [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pvpq]))] + [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pq]))] @info("The powerflow solver with KLU converged after $i iterations") end return (converged, V, Sbus_result)