From 30c345ec91564f25104f6c6bb6015551689a5e3a Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Tue, 6 Aug 2024 14:01:39 -0700 Subject: [PATCH 1/8] use function for initializing data --- src/PowerFlowData.jl | 80 +++++++++++++++++++------------------------- src/common.jl | 21 ++++++++++++ 2 files changed, 55 insertions(+), 46 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 669123ac..43470a12 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -95,7 +95,7 @@ NOTE: use it for AC power flow computations. # Arguments: - `::ACPowerFlow`: - use ACPowerFlow() to evaluate the AC OPF. + use ACPowerFlow() to evaluate the AC PF. - `sys::PSY.System`: container storing the system data to consider in the PowerFlowData structure. @@ -109,7 +109,7 @@ NOTE: use it for AC power flow computations. - `check_connectivity::Bool`: Perform connectivity check on the network matrix. Default value = true. -WARNING: functions for the evaluation of the multi-period AC OPF still to be implemented. +WARNING: functions for the evaluation of the multi-period AC PF still to be implemented. """ function PowerFlowData( ::ACPowerFlow, @@ -151,17 +151,14 @@ function PowerFlowData( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) - for (bus_no, ix) in bus_lookup - bus_name = temp_bus_map[bus_no] - bus = PSY.get_component(PSY.Bus, sys, bus_name) - bus_type[ix] = PSY.get_bustype(bus) - if bus_type[ix] == PSY.ACBusTypes.REF - bus_angles[ix] = 0.0 - else - bus_angles[ix] = PSY.get_angle(bus) - end - bus_magnitude[ix] = PSY.get_magnitude(bus) - end + _initialize_bus_data!( + bus_type, + bus_angles, + bus_magnitude, + temp_bus_map, + bus_lookup, + sys, + ) bus_activepower_injection = zeros(Float64, n_buses) bus_reactivepower_injection = zeros(Float64, n_buses) @@ -281,17 +278,14 @@ function PowerFlowData( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.ACBus, sys) ) - for (bus_no, ix) in bus_lookup - bus_name = temp_bus_map[bus_no] - bus = PSY.get_component(PSY.Bus, sys, bus_name) - bus_type[ix] = PSY.get_bustype(bus) - if bus_type[ix] == PSY.ACBusTypes.REF - bus_angles[ix] = 0.0 - else - bus_angles[ix] = PSY.get_angle(bus) - end - bus_magnitude[ix] = PSY.get_magnitude(bus) - end + _initialize_bus_data!( + bus_type, + bus_angles, + bus_magnitude, + temp_bus_map, + bus_lookup, + sys, + ) # define injection vectors related to the first timestep bus_activepower_injection = zeros(Float64, n_buses) @@ -410,17 +404,14 @@ function PowerFlowData( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) - for (bus_no, ix) in bus_lookup - bus_name = temp_bus_map[bus_no] - bus = PSY.get_component(PSY.Bus, sys, bus_name) - bus_type[ix] = PSY.get_bustype(bus) - if bus_type[ix] == PSY.ACBusTypes.REF - bus_angles[ix] = 0.0 - else - bus_angles[ix] = PSY.get_angle(bus) - end - bus_magnitude[ix] = PSY.get_magnitude(bus) - end + _initialize_bus_data!( + bus_type, + bus_angles, + bus_magnitude, + temp_bus_map, + bus_lookup, + sys, + ) # define injection vectors related to the first timestep bus_activepower_injection = zeros(Float64, n_buses) @@ -539,17 +530,14 @@ function PowerFlowData( PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys) ) - for (bus_no, ix) in bus_lookup - bus_name = temp_bus_map[bus_no] - bus = PSY.get_component(PSY.Bus, sys, bus_name) - bus_type[ix] = PSY.get_bustype(bus) - if bus_type[ix] == PSY.ACBusTypes.REF - bus_angles[ix] = 0.0 - else - bus_angles[ix] = PSY.get_angle(bus) - end - bus_magnitude[ix] = PSY.get_magnitude(bus) - end + _initialize_bus_data!( + bus_type, + bus_angles, + bus_magnitude, + temp_bus_map, + bus_lookup, + sys, + ) # define injection vectors related to the first timestep bus_activepower_injection = zeros(Float64, n_buses) diff --git a/src/common.jl b/src/common.jl index 6ef897b8..6f33f91c 100644 --- a/src/common.jl +++ b/src/common.jl @@ -81,6 +81,27 @@ function _get_reactive_power_bound!( end end end + +function _initialize_bus_data!( + bus_type::Vector{PSY.ACBusTypes}, + bus_angles::Vector{Float64}, + bus_magnitude::Vector{Float64}, + temp_bus_map::Dict{Int, String}, + bus_lookup::Dict{Int, Int}, + sys::PSY.System, +) + for (bus_no, ix) in bus_lookup + bus_name = temp_bus_map[bus_no] + bus = PSY.get_component(PSY.Bus, sys, bus_name) + bus_type[ix] = PSY.get_bustype(bus) + if bus_type[ix] == PSY.ACBusTypes.REF + bus_angles[ix] = 0.0 + else + bus_angles[ix] = PSY.get_angle(bus) + end + bus_magnitude[ix] = PSY.get_magnitude(bus) + end +end ############################################################################## # Matrix Methods ############################################################# From 518b14de7a71b33622c95e9dfcf3c3d9d92beee9 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Tue, 6 Aug 2024 14:01:44 -0700 Subject: [PATCH 2/8] fix redistribution of Q --- src/post_processing.jl | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/post_processing.jl b/src/post_processing.jl index 9e759f19..9d6bdcbb 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -387,13 +387,13 @@ function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus: else PSY.InfrastructureSystems.@assert_op fraction > 0 end - current_q = PSY.get_reactive_power(d) q_frac = q_residual * fraction q_set_point = clamp(q_frac + current_q, q_limits.min, q_limits.max) - reallocated_q += q_frac - if (q_frac >= q_limits.max + BOUNDS_TOLERANCE) || - (q_frac <= q_limits.min - BOUNDS_TOLERANCE) + # Assign new capacity based on the limits and the fraction + reallocated_q += q_set_point - current_q + if ((q_frac + current_q) >= q_limits.max + BOUNDS_TOLERANCE) || + ((q_frac + current_q) <= q_limits.min - BOUNDS_TOLERANCE) push!(units_at_limit, ix) @warn "Unit $(PSY.get_name(d)) set at the limit $(q_set_point). Q_max = $(q_limits.max) Q_min = $(q_limits.min)" end @@ -405,12 +405,14 @@ function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus: break end it += 1 - if it > 1 + if it > 5 + @warn "Maximum number of iterations for Q-redistribution reached. Number of devices at Q limit are: $(length(units_at_limit)) of $(length(devices)) available devices" break end end end + # Last attempt to allocate reactive power 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 @@ -425,6 +427,12 @@ function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus: end end + @assert isapprox( + sum(PSY.get_reactive_power.(devices)), + Q_gen; + atol = ISAPPROX_ZERO_TOLERANCE, + ) + return end @@ -529,7 +537,7 @@ end """ Returns a dictionary containing the DC power flow results. Each key conresponds -to the name of the considered time periods, storing a DataFrame with the OPF +to the name of the considered time periods, storing a DataFrame with the PF results. # Arguments: From 23d016247191bf96a1cdcffef875d6337027db6e Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Tue, 6 Aug 2024 14:04:02 -0700 Subject: [PATCH 3/8] add reciprocal for DC PowerFlow branches --- src/post_processing.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/post_processing.jl b/src/post_processing.jl index 9d6bdcbb..d4542de0 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -496,14 +496,10 @@ function _allocate_results_data( Q_from_to_vect = zeros(length(branches)) P_to_from_vect = zeros(length(branches)) Q_to_from_vect = zeros(length(branches)) + # branch_flow_values is always from_to direction for i in 1:length(branches) - if branch_flow_values[i] >= 0 - P_from_to_vect[i] = branch_flow_values[i] - P_to_from_vect[i] = 0 - else - P_from_to_vect[i] = 0 - P_to_from_vect[i] = branch_flow_values[i] - end + P_from_to_vect[i] = branch_flow_values[i] + P_to_from_vect[i] = -branch_flow_values[i] end bus_df = DataFrames.DataFrame(; From b520c2b10dbe3c3bdba3fd289f3d979508556ed2 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Wed, 7 Aug 2024 11:55:02 -0700 Subject: [PATCH 4/8] update tests --- test/test_dc_powerflow.jl | 9 +++------ test/test_multiperiod_dc_powerflow.jl | 12 +++--------- 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/test/test_dc_powerflow.jl b/test/test_dc_powerflow.jl index 5a487e9c..ca5a2f1e 100644 --- a/test/test_dc_powerflow.jl +++ b/test/test_dc_powerflow.jl @@ -28,8 +28,7 @@ # CASE 1: ABA and BA matrices solved_data_ABA = solve_powerflow(DCPowerFlow(), sys) @test isapprox( - solved_data_ABA["1"]["flow_results"].P_from_to + - solved_data_ABA["1"]["flow_results"].P_to_from, + solved_data_ABA["1"]["flow_results"].P_from_to, ref_flow_values, atol = 1e-6, ) @@ -38,8 +37,7 @@ # CASE 2: PTDF and ABA MATRICES solved_data_PTDF = solve_powerflow(PTDFDCPowerFlow(), sys) @test isapprox( - solved_data_PTDF["1"]["flow_results"].P_from_to + - solved_data_PTDF["1"]["flow_results"].P_to_from, + solved_data_PTDF["1"]["flow_results"].P_from_to, ref_flow_values, atol = 1e-6, ) @@ -48,8 +46,7 @@ # CASE 3: VirtualPTDF and ABA MATRICES solved_data_vPTDF = solve_powerflow(vPTDFDCPowerFlow(), sys) @test isapprox( - solved_data_vPTDF["1"]["flow_results"].P_from_to + - solved_data_vPTDF["1"]["flow_results"].P_to_from, + solved_data_vPTDF["1"]["flow_results"].P_from_to, ref_flow_values, atol = 1e-6, ) diff --git a/test/test_multiperiod_dc_powerflow.jl b/test/test_multiperiod_dc_powerflow.jl index 74633e4f..bce5ecd8 100644 --- a/test/test_multiperiod_dc_powerflow.jl +++ b/test/test_multiperiod_dc_powerflow.jl @@ -56,25 +56,19 @@ # case 1 for i in 1:length(data_1.timestep_map) - flow_from_to = results_1[data_1.timestep_map[i]]["flow_results"].P_from_to - flow_to_from = results_1[data_1.timestep_map[i]]["flow_results"].P_to_from - net_flow = flow_from_to + flow_to_from + net_flow = results_1[data_1.timestep_map[i]]["flow_results"].P_from_to @test isapprox(net_flow, flows[:, i], atol = 1e-5) end # case 2 for i in 1:length(data_1.timestep_map) - flow_from_to = results_1[data_2.timestep_map[i]]["flow_results"].P_from_to - flow_to_from = results_1[data_2.timestep_map[i]]["flow_results"].P_to_from - net_flow = flow_from_to + flow_to_from + net_flow = results_1[data_2.timestep_map[i]]["flow_results"].P_from_to @test isapprox(net_flow, flows[:, i], atol = 1e-5) end # case 3 for i in 1:length(data_1.timestep_map) - flow_from_to = results_1[data_3.timestep_map[i]]["flow_results"].P_from_to - flow_to_from = results_1[data_3.timestep_map[i]]["flow_results"].P_to_from - net_flow = flow_from_to + flow_to_from + net_flow = results_1[data_3.timestep_map[i]]["flow_results"].P_from_to @test isapprox(net_flow, flows[:, i], atol = 1e-5) end end From 8137ac6077d45f9830c228323269a9c1824e11ba Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Mon, 12 Aug 2024 09:45:25 -0700 Subject: [PATCH 5/8] update to_from dc testing --- test/test_dc_powerflow.jl | 15 +++++++++++++++ test/test_multiperiod_dc_powerflow.jl | 6 ++++++ 2 files changed, 21 insertions(+) diff --git a/test/test_dc_powerflow.jl b/test/test_dc_powerflow.jl index ca5a2f1e..5c7d7c21 100644 --- a/test/test_dc_powerflow.jl +++ b/test/test_dc_powerflow.jl @@ -32,6 +32,11 @@ ref_flow_values, atol = 1e-6, ) + @test isapprox( + solved_data_ABA["1"]["flow_results"].P_to_from, + -ref_flow_values, + atol = 1e-6, + ) @test isapprox(solved_data_ABA["1"]["bus_results"].θ, ref_bus_angles, atol = 1e-6) # CASE 2: PTDF and ABA MATRICES @@ -41,6 +46,11 @@ ref_flow_values, atol = 1e-6, ) + @test isapprox( + solved_data_PTDF["1"]["flow_results"].P_to_from, + -ref_flow_values, + atol = 1e-6, + ) @test isapprox(solved_data_PTDF["1"]["bus_results"].θ, ref_bus_angles, atol = 1e-6) # CASE 3: VirtualPTDF and ABA MATRICES @@ -50,5 +60,10 @@ ref_flow_values, atol = 1e-6, ) + @test isapprox( + solved_data_vPTDF["1"]["flow_results"].P_to_from, + -ref_flow_values, + atol = 1e-6, + ) @test isapprox(solved_data_vPTDF["1"]["bus_results"].θ, ref_bus_angles, atol = 1e-6) end diff --git a/test/test_multiperiod_dc_powerflow.jl b/test/test_multiperiod_dc_powerflow.jl index bce5ecd8..7f252e83 100644 --- a/test/test_multiperiod_dc_powerflow.jl +++ b/test/test_multiperiod_dc_powerflow.jl @@ -57,18 +57,24 @@ # case 1 for i in 1:length(data_1.timestep_map) net_flow = results_1[data_1.timestep_map[i]]["flow_results"].P_from_to + net_flow_tf = results_1[data_1.timestep_map[i]]["flow_results"].P_to_from @test isapprox(net_flow, flows[:, i], atol = 1e-5) + @test isapprox(net_flow_tf, -flows[:, i], atol = 1e-5) end # case 2 for i in 1:length(data_1.timestep_map) net_flow = results_1[data_2.timestep_map[i]]["flow_results"].P_from_to + net_flow_tf = results_1[data_2.timestep_map[i]]["flow_results"].P_to_from @test isapprox(net_flow, flows[:, i], atol = 1e-5) + @test isapprox(net_flow_tf, -flows[:, i], atol = 1e-5) end # case 3 for i in 1:length(data_1.timestep_map) net_flow = results_1[data_3.timestep_map[i]]["flow_results"].P_from_to + net_flow_tf = results_1[data_3.timestep_map[i]]["flow_results"].P_to_from @test isapprox(net_flow, flows[:, i], atol = 1e-5) + @test isapprox(net_flow_tf, -flows[:, i], atol = 1e-5) end end From 392da2f06368dd6964a7758e59aecbe58107104d Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Mon, 12 Aug 2024 09:45:39 -0700 Subject: [PATCH 6/8] add iterations kwarg to redistribution --- src/nlsolve_ac_powerflow.jl | 3 ++- src/post_processing.jl | 26 ++++++++++++++++++-------- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index c518ea88..bc9ccfc8 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -41,9 +41,10 @@ function solve_ac_powerflow!(system::PSY.System; kwargs...) system; check_connectivity = get(kwargs, :check_connectivity, true), ) + max_iterations = get(kwargs, :max_redistribution_iterations, 10) res = _solve_powerflow!(data, check_reactive_power_limits; kwargs...) if res.f_converged - write_powerflow_solution!(system, res.zero) + write_powerflow_solution!(system, res.zero, 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) diff --git a/src/post_processing.jl b/src/post_processing.jl index d4542de0..20c23810 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -175,6 +175,7 @@ function _power_redistribution_ref( P_gen::Float64, Q_gen::Float64, bus::PSY.Bus, + max_iterations::Int, ) devices_ = PSY.get_components(x -> _is_available_source(x, bus), PSY.StaticInjection, sys) @@ -199,7 +200,7 @@ function _power_redistribution_ref( if length(devices_) == 1 device = first(devices_) PSY.set_active_power!(device, P_gen) - _reactive_power_redistribution_pv(sys, Q_gen, bus) + _reactive_power_redistribution_pv(sys, Q_gen, bus, max_iterations) return elseif length(devices_) > 1 devices = @@ -258,7 +259,7 @@ function _power_redistribution_ref( break end it += 1 - if it > 10 + if it > max_iterations break end end @@ -276,11 +277,16 @@ function _power_redistribution_ref( end end end - _reactive_power_redistribution_pv(sys, Q_gen, bus) + _reactive_power_redistribution_pv(sys, Q_gen, bus, max_iterations) return end -function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus::PSY.Bus) +function _reactive_power_redistribution_pv( + sys::PSY.System, + Q_gen::Float64, + bus::PSY.Bus, + max_iterations::Int, +) @debug "Reactive Power Distribution $(PSY.get_name(bus))" devices_ = PSY.get_components(x -> _is_available_source(x, bus), PSY.StaticInjection, sys) @@ -405,7 +411,7 @@ function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus: break end it += 1 - if it > 5 + if it > max_iterations @warn "Maximum number of iterations for Q-redistribution reached. Number of devices at Q limit are: $(length(units_at_limit)) of $(length(devices)) available devices" break end @@ -439,7 +445,11 @@ end """ Updates system voltages and powers with power flow results """ -function write_powerflow_solution!(sys::PSY.System, result::Vector{Float64}) +function write_powerflow_solution!( + sys::PSY.System, + result::Vector{Float64}, + max_iterations::Int, +) buses = enumerate( sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x)), ) @@ -448,11 +458,11 @@ function write_powerflow_solution!(sys::PSY.System, result::Vector{Float64}) if bus.bustype == PSY.ACBusTypes.REF P_gen = result[2 * ix - 1] Q_gen = result[2 * ix] - _power_redistribution_ref(sys, P_gen, Q_gen, bus) + _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] - _reactive_power_redistribution_pv(sys, Q_gen, bus) + _reactive_power_redistribution_pv(sys, Q_gen, bus, max_iterations) elseif bus.bustype == PSY.ACBusTypes.PQ Vm = result[2 * ix - 1] θ = result[2 * ix] From 5b6589d930cd1052f17af2a55dc2bbbe564b26d6 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Mon, 12 Aug 2024 09:50:08 -0700 Subject: [PATCH 7/8] add default const --- src/definitions.jl | 1 + src/nlsolve_ac_powerflow.jl | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/definitions.jl b/src/definitions.jl index 03af2999..d8ecad1a 100644 --- a/src/definitions.jl +++ b/src/definitions.jl @@ -6,6 +6,7 @@ const MAX_INIT_RESIDUAL = 10.0 const MAX_NLSOLVE_INTERATIONS = 10 const BOUNDS_TOLERANCE = 1e-6 const MAX_REACTIVE_POWER_ITERATIONS = 10 +const DEFAULT_MAX_REDISTRIBUTION_ITERATIONS = 10 const ISAPPROX_ZERO_TOLERANCE = 1e-6 diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index bc9ccfc8..6eb29a2d 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -41,7 +41,8 @@ function solve_ac_powerflow!(system::PSY.System; kwargs...) system; check_connectivity = get(kwargs, :check_connectivity, true), ) - max_iterations = get(kwargs, :max_redistribution_iterations, 10) + max_iterations = + get(kwargs, :max_redistribution_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) From e04428c47dfc1f4c11302257679e8215a3f571e8 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Mon, 12 Aug 2024 10:06:35 -0700 Subject: [PATCH 8/8] remove kwarg --- src/nlsolve_ac_powerflow.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 6eb29a2d..813848b0 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -41,8 +41,7 @@ function solve_ac_powerflow!(system::PSY.System; kwargs...) system; check_connectivity = get(kwargs, :check_connectivity, true), ) - max_iterations = - get(kwargs, :max_redistribution_iterations, DEFAULT_MAX_REDISTRIBUTION_ITERATIONS) + 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)