From 691a73a1697489d548206b751ad04c79be229c25 Mon Sep 17 00:00:00 2001 From: Roman Bolgaryn Date: Thu, 16 Jan 2025 09:42:27 -0700 Subject: [PATCH] small change to _update_system! --- src/ac_power_flow.jl | 6 ++- src/newton_ac_powerflow.jl | 74 +++++++++++++++++++++----------- src/nlsolve_powerflow.jl | 6 ++- src/post_processing.jl | 58 ++++++++++++++++++++----- test/test_newton_ac_powerflow.jl | 17 ++++++-- test/test_utils/legacy_pf.jl | 10 ++++- 6 files changed, 126 insertions(+), 45 deletions(-) diff --git a/src/ac_power_flow.jl b/src/ac_power_flow.jl index 6daf57c7..ee7b2f2a 100644 --- a/src/ac_power_flow.jl +++ b/src/ac_power_flow.jl @@ -61,9 +61,11 @@ function PolarPowerFlow(data::ACPowerFlowData; time_step::Int64 = 1) Q_net = zeros(n_buses) for ix in 1:n_buses P_net[ix] = - data.bus_activepower_injection[ix, time_step] - data.bus_activepower_withdrawals[ix, time_step] + data.bus_activepower_injection[ix, time_step] - + data.bus_activepower_withdrawals[ix, time_step] Q_net[ix] = - data.bus_reactivepower_injection[ix, time_step] - data.bus_reactivepower_withdrawals[ix, time_step] + data.bus_reactivepower_injection[ix, time_step] - + data.bus_reactivepower_withdrawals[ix, time_step] end x0 = _calculate_x0(1, data.bus_type[:, time_step], diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index f2eaefc5..570732f8 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -33,7 +33,7 @@ solve_ac_powerflow!(sys, method=:newton) function solve_powerflow!( pf::ACPowerFlow{<:ACPowerFlowSolverType}, system::PSY.System; - time_step::Int64=1, + time_step::Int64 = 1, kwargs..., ) #Save per-unit flag @@ -47,11 +47,16 @@ function solve_powerflow!( check_connectivity = get(kwargs, :check_connectivity, true), ) - converged, V, Sbus_result = _ac_powereflow(data, pf; time_step=time_step, kwargs...) + converged, V, Sbus_result = _ac_powereflow(data, pf; time_step = time_step, kwargs...) x = _calc_x(data, V, Sbus_result) if converged - write_powerflow_solution!(system, x, data, get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER)) + write_powerflow_solution!( + system, + x, + data, + get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER), + ) @info("PowerFlow solve converged, the results have been stored in the system") else @error("The powerflow solver returned convergence = $(converged)") @@ -163,15 +168,13 @@ function solve_powerflow( return results end - - # Multiperiod power flow - work in progress function solve_powerflow!( data::ACPowerFlowData; kwargs..., ) pf = ACPowerFlow() # todo: somehow store in data which PF to use (see issue #50) - + sorted_time_steps = sort(collect(keys(data.timestep_map))) # preallocate results ts_converged = zeros(Bool, 1, length(sorted_time_steps)) @@ -191,38 +194,48 @@ function solve_powerflow!( 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]) + 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 # write results for REF - data.bus_activepower_injection[ref, t] .= real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t] - data.bus_reactivepower_injection[ref, t] .= imag.(Sbus_result[ref]) .+ data.bus_reactivepower_withdrawals[ref, t] + data.bus_activepower_injection[ref, t] .= + real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t] + data.bus_reactivepower_injection[ref, t] .= + imag.(Sbus_result[ref]) .+ data.bus_reactivepower_withdrawals[ref, t] # write Q results for PV - data.bus_reactivepower_injection[pv, t] .= imag.(Sbus_result[pv]) .+ data.bus_reactivepower_withdrawals[pv, t] + data.bus_reactivepower_injection[pv, t] .= + 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 + 1 + 2 end end # write branch flows - Sft = ts_V[fb,:] .* conj.(Yft * ts_V) - Stf = ts_V[tb,:] .* conj.(Ytf * ts_V) + Sft = ts_V[fb, :] .* conj.(Yft * ts_V) + Stf = ts_V[tb, :] .* conj.(Ytf * ts_V) data.branch_activepower_flow_from_to .= real.(Sft) data.branch_reactivepower_flow_from_to .= imag.(Sft) @@ -235,7 +248,6 @@ function solve_powerflow!( return end - function _ac_powereflow( data::ACPowerFlowData, pf::ACPowerFlow{<:ACPowerFlowSolverType}; @@ -275,12 +287,14 @@ function _check_q_limit_bounds!( @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ" within_limits = false data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ - data.bus_reactivepower_injection[ix, time_step] = data.bus_reactivepower_bounds[ix, time_step][1] + data.bus_reactivepower_injection[ix, time_step] = + data.bus_reactivepower_bounds[ix, time_step][1] elseif Q_gen >= data.bus_reactivepower_bounds[ix, time_step][2] @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ" within_limits = false data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ - data.bus_reactivepower_injection[ix, time_step] = data.bus_reactivepower_bounds[ix, time_step][2] + data.bus_reactivepower_injection[ix, time_step] = + data.bus_reactivepower_bounds[ix, time_step][2] else @debug "Within Limits" end @@ -296,7 +310,8 @@ function _solve_powerflow!( kwargs..., ) for _ in 1:MAX_REACTIVE_POWER_ITERATIONS - converged, V, Sbus_result = _newton_powerflow(pf, data; time_step=time_step, kwargs...) + converged, V, Sbus_result = + _newton_powerflow(pf, data; time_step = time_step, kwargs...) if !converged || !check_reactive_power_limits || _check_q_limit_bounds!(data, Sbus_result, time_step) return converged, V, Sbus_result @@ -525,9 +540,18 @@ function _newton_powerflow( Ybus = data.power_network_matrix.data # Find indices for each bus type - ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, time_step]) - pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step]) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step]) + ref = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, + data.bus_type[:, time_step], + ) + pv = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, + data.bus_type[:, time_step], + ) + pq = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, + data.bus_type[:, time_step], + ) pvpq = [pv; pq] # nref = length(ref) diff --git a/src/nlsolve_powerflow.jl b/src/nlsolve_powerflow.jl index 77274ade..c04e5fca 100644 --- a/src/nlsolve_powerflow.jl +++ b/src/nlsolve_powerflow.jl @@ -1,4 +1,5 @@ -const _NLSOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity]) +const _NLSOLVE_AC_POWERFLOW_KWARGS = + Set([:check_reactive_power_limits, :check_connectivity]) function _newton_powerflow( pf::ACPowerFlow{NLSolveACPowerFlow}, data::ACPowerFlowData; @@ -13,7 +14,8 @@ function _newton_powerflow( ) end - nlsolve_solver_kwargs = filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs) + nlsolve_solver_kwargs = + filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs) pf = PolarPowerFlow(data) J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) diff --git a/src/post_processing.jl b/src/post_processing.jl index b771b062..9faa256c 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -709,25 +709,63 @@ if a new `PowerFlowData` is constructed from the resulting system it is the same See also `write_powerflow_solution!`. NOTE that this assumes that `data` was initialized from `sys` and then solved with no further modifications. """ -function update_system!(sys::PSY.System, data::PowerFlowData; time_step=1) +function update_system!(sys::PSY.System, data::PowerFlowData; time_step = 1) for bus in PSY.get_components(PSY.Bus, sys) - if bus.bustype == PSY.ACBusTypes.REF + bus_number = data.bus_lookup[PSY.get_number(bus)] + bus_type = data.bus_type[bus_number, time_step] # use this instead of bus.bustype to account for PV -> PQ + if bus_type == PSY.ACBusTypes.REF # For REF bus, voltage and angle are fixed; update active and reactive - P_gen = data.bus_activepower_injection[data.bus_lookup[PSY.get_number(bus)], time_step] - Q_gen = data.bus_reactivepower_injection[data.bus_lookup[PSY.get_number(bus)], time_step] + P_gen = data.bus_activepower_injection[bus_number, time_step] + Q_gen = data.bus_reactivepower_injection[bus_number, time_step] _power_redistribution_ref(sys, P_gen, Q_gen, bus, DEFAULT_MAX_REDISTRIBUTION_ITERATIONS) - elseif bus.bustype == PSY.ACBusTypes.PV + elseif bus_type == PSY.ACBusTypes.PV # For PV bus, active and voltage are fixed; update reactive and angle - Q_gen = data.bus_reactivepower_injection[data.bus_lookup[PSY.get_number(bus)], time_step] + Q_gen = data.bus_reactivepower_injection[bus_number, time_step] _reactive_power_redistribution_pv(sys, Q_gen, bus, DEFAULT_MAX_REDISTRIBUTION_ITERATIONS) - PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)], time_step]) - elseif bus.bustype == PSY.ACBusTypes.PQ + PSY.set_angle!(bus, data.bus_angles[bus_number, time_step]) + elseif bus_type == PSY.ACBusTypes.PQ # For PQ bus, active and reactive are fixed; update voltage and angle - Vm = data.bus_magnitude[data.bus_lookup[PSY.get_number(bus)], time_step] + Vm = data.bus_magnitude[bus_number, time_step] PSY.set_magnitude!(bus, Vm) - PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)], time_step]) + PSY.set_angle!(bus, data.bus_angles[bus_number, time_step]) + # if it used to be a PV bus, also set the Q value: + if bus.bustype == PSY.ACBusTypes.PV + Q_gen = data.bus_reactivepower_injection[bus_number, time_step] + _reactive_power_redistribution_pv(sys, Q_gen, bus, + DEFAULT_MAX_REDISTRIBUTION_ITERATIONS) + # now both the Q and the Vm, Va are correct for this kind of buses + end + end + end +end + +# This cannot work because we do not know the redistribution keys from the power values of devices +# to the power values of buses. Leaving this here anyways to document this issue. +function set_system_time_step(sys::PSY.System, data::PowerFlowData; time_step = 1) + @error("This function cannot work") + for bus in PSY.get_components(PSY.Bus, sys) + bus_number = data.bus_lookup[PSY.get_number(bus)] + bus_type = data.bus_type[bus_number, time_step] + sys_bus_type = sys.bustype # here we should be using the system bus type to know which inputs to modify + + # for REF bus, we don't need to set anything + if sys_bus_type == PSY.ACBusTypes.PV + # we need to modify the P value + # TODO if we also modify time-series data for voltage set-points, also set that value. + # However, we cannot do this easily because it can happen that the PV -> PQ transformation + # caused a different Vm value as the result, while the set-point Vm is still the same. + # We need a solution for this. + P_gen = data.bus_activepower_injection[bus_number, time_step] + Q_gen = 0.0 + # TODO write time-series value of the bus back to devices (not possible) + elseif sys_bus_type == PSY.ACBusTypes.PQ + P_gen = data.bus_activepower_injection[bus_number, time_step] + Q_gen = data.bus_reactivepower_injection[bus_number, time_step] + P_load = data.bus_activepower_withdrawals[bus_number, time_step] + Q_load = data.bus_reactivepower_withdrawals[bus_number, time_step] + # TODO write time-series value of the bus back to devices (not possible) end end end diff --git a/test/test_newton_ac_powerflow.jl b/test/test_newton_ac_powerflow.jl index d029a49b..beed6970 100644 --- a/test/test_newton_ac_powerflow.jl +++ b/test/test_newton_ac_powerflow.jl @@ -505,12 +505,21 @@ end pf_default, sys; check_connectivity = true) - + time_step = 1 - ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, time_step]) - pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step]) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step]) + ref = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, + data.bus_type[:, time_step], + ) + pv = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, + data.bus_type[:, time_step], + ) + pq = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, + data.bus_type[:, time_step], + ) pvpq = [pv; pq] npvpq = length(pvpq) diff --git a/test/test_utils/legacy_pf.jl b/test/test_utils/legacy_pf.jl index eafb8077..5436b334 100644 --- a/test/test_utils/legacy_pf.jl +++ b/test/test_utils/legacy_pf.jl @@ -47,8 +47,14 @@ function _newton_powerflow( # 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[:, time_step]) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step]) + pv = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, + data.bus_type[:, time_step], + ) + pq = findall( + x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, + data.bus_type[:, time_step], + ) pvpq = [pv; pq] #nref = length(ref)