Fix Bugs/Issues of PowerFlow

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)

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 #############################################################
 

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
 
@@ -488,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(;
@@ -529,7 +533,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:

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,
     )

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