Fix Bugs/Issues of PowerFlow (#45)

* use function for initializing data * fix redistribution of Q * add reciprocal for DC PowerFlow branches * update tests * update to_from dc testing * add iterations kwarg to redistribution * add default const * remove kwarg
NREL-Sienna · Aug 12, 2024 · 035a207 · 035a207
1 parent 558a337
commit 035a207
Show file tree

Hide file tree

Showing 7 changed files with 119 additions and 82 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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