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 ############################################################# 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 c518ea88..813848b0 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 = DEFAULT_MAX_REDISTRIBUTION_ITERATIONS 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 9e759f19..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) @@ -387,13 +393,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 +411,14 @@ function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus: break end it += 1 - if it > 1 + 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 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,13 +433,23 @@ 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 """ 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)), ) @@ -440,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] @@ -488,14 +506,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(; @@ -529,7 +543,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: diff --git a/test/test_dc_powerflow.jl b/test/test_dc_powerflow.jl index 5a487e9c..5c7d7c21 100644 --- a/test/test_dc_powerflow.jl +++ b/test/test_dc_powerflow.jl @@ -28,30 +28,42 @@ # 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, ) + @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 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, ) + @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 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, ) + @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 74633e4f..7f252e83 100644 --- a/test/test_multiperiod_dc_powerflow.jl +++ b/test/test_multiperiod_dc_powerflow.jl @@ -56,25 +56,25 @@ # 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 + 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) - 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 + 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) - 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 + 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