Merge pull request #78 from NREL-Sienna/gks/update_pfd_for_ac_multipe…

…riod

Modify `PowerFlowData` fields for multi-period AC PF

src/PowerFlowData.jl

@@ -38,20 +38,33 @@ flows and angles, as well as these ones.
 - `bus_reactivepower_withdrawals::Matrix{Float64}`:
         "(b, t)" matrix containing the bus reactive power withdrawals. b:
         number of buses, t: number of time period.
-- `bus_reactivepower_bounds::Vector{Float64}`: Upper and Lower bounds for the reactive supply
-        at each bus.
-- `bus_type::Vector{PSY.ACBusTypes}`:
-        vector containing type of buses present in the system, ordered
-        according to "bus_lookup".
+- `bus_reactivepower_bounds::Matrix{Float64}`:
+        "(b, t)" matrix containing upper and lower bounds for the reactive supply at each
+        bus at each time period.
+- `bus_type::Matrix{PSY.ACBusTypes}`:
+        "(b, t)" matrix containing type of buses present in the system, ordered
+        according to "bus_lookup," at each time period.
 - `bus_magnitude::Matrix{Float64}`:
         "(b, t)" matrix containing the bus magnitudes, ordered according to
         "bus_lookup". b: number of buses, t: number of time period.
 - `bus_angles::Matrix{Float64}`:
         "(b, t)" matrix containing the bus angles, ordered according to
         "bus_lookup". b: number of buses, t: number of time period.
-- `branch_flow_values::Matrix{Float64}`:
-        "(br, t)" matrix containing the power flows, ordered according to
-        "branch_lookup". br: number of branches, t: number of time period.
+- `branch_activepower_flow_from_to::Matrix{Float64}`:
+        "(br, t)" matrix containing the active power flows measured at the `from` bus,
+        ordered according to "branch_lookup". br: number of branches, t: number of time
+        period.
+- `branch_reactivepower_flow_from_to::Matrix{Float64}`:
+        "(br, t)" matrix containing the reactive power flows measured at the `from` bus,
+        ordered according to "branch_lookup". br: number of branches, t: number of time
+        period.
+- `branch_activepower_flow_to_from::Matrix{Float64}`:
+        "(br, t)" matrix containing the active power flows measured at the `to` bus, ordered
+        according to "branch_lookup". br: number of branches, t: number of time period.
+- `branch_reactivepower_flow_to_from::Matrix{Float64}`:
+        "(br, t)" matrix containing the reactive power flows measured at the `to` bus,
+        ordered according to "branch_lookup". br: number of branches, t: number of time
+        period.
 - `timestep_map::Dict{Int, S}`:
         dictonary mapping the number of the time periods (corresponding to the
         column number of the previosly mentioned matrices) and their names.
@@ -75,12 +88,15 @@ struct PowerFlowData{
     bus_reactivepower_injection::Matrix{Float64}
     bus_activepower_withdrawals::Matrix{Float64}
     bus_reactivepower_withdrawals::Matrix{Float64}
-    bus_reactivepower_bounds::Vector{Vector{Float64}}
-    bus_type::Vector{PSY.ACBusTypes}
+    bus_reactivepower_bounds::Matrix{Vector{Float64}}
+    bus_type::Matrix{PSY.ACBusTypes}
     branch_type::Vector{DataType}
     bus_magnitude::Matrix{Float64}
     bus_angles::Matrix{Float64}
-    branch_flow_values::Matrix{Float64}
+    branch_activepower_flow_from_to::Matrix{Float64}
+    branch_reactivepower_flow_from_to::Matrix{Float64}
+    branch_activepower_flow_to_from::Matrix{Float64}
+    branch_reactivepower_flow_to_from::Matrix{Float64}
     timestep_map::Dict{Int, String}
     valid_ix::Vector{Int}
     power_network_matrix::M
@@ -99,7 +115,14 @@ get_bus_type(pfd::PowerFlowData) = pfd.bus_type
 get_branch_type(pfd::PowerFlowData) = pfd.branch_type
 get_bus_magnitude(pfd::PowerFlowData) = pfd.bus_magnitude
 get_bus_angles(pfd::PowerFlowData) = pfd.bus_angles
-get_branch_flow_values(pfd::PowerFlowData) = pfd.branch_flow_values
+get_branch_activepower_flow_from_to(pfd::PowerFlowData) =
+    pfd.branch_activepower_flow_from_to
+get_branch_reactivepower_flow_from_to(pfd::PowerFlowData) =
+    pfd.branch_reactivepower_flow_from_to
+get_branch_activepower_flow_to_from(pfd::PowerFlowData) =
+    pfd.branch_activepower_flow_to_from
+get_branch_reactivepower_flow_to_from(pfd::PowerFlowData) =
+    pfd.branch_reactivepower_flow_to_from
 get_timestep_map(pfd::PowerFlowData) = pfd.timestep_map
 get_valid_ix(pfd::PowerFlowData) = pfd.valid_ix
 get_power_network_matrix(pfd::PowerFlowData) = pfd.power_network_matrix
@@ -195,11 +218,6 @@ function PowerFlowData(
         branch_types[ix] = typeof(b)
     end
 
-    bus_reactivepower_bounds = Vector{Vector{Float64}}(undef, n_buses)
-    for i in 1:n_buses
-        bus_reactivepower_bounds[i] = [0.0, 0.0]
-    end
-    _get_reactive_power_bound!(bus_reactivepower_bounds, bus_lookup, sys)
     timestep_map = Dict(1 => "1")
     valid_ix = setdiff(1:n_buses, ref_bus_positions)
     neighbors = _calculate_neighbors(power_network_matrix)
@@ -216,7 +234,6 @@ function PowerFlowData(
         branch_lookup,
         temp_bus_map,
         branch_types,
-        bus_reactivepower_bounds,
         timestep_map,
         valid_ix,
         neighbors,

src/ac_power_flow.jl

@@ -56,7 +56,8 @@ function _calculate_x0(n::Int,
 end
 
 function PolarPowerFlow(data::ACPowerFlowData)
-    n_buses = length(data.bus_type)
+    time_step = 1  # TODO placeholder time_step
+    n_buses = first(size(data.bus_type))
     P_net = zeros(n_buses)
     Q_net = zeros(n_buses)
     for ix in 1:n_buses
@@ -66,7 +67,7 @@ function PolarPowerFlow(data::ACPowerFlowData)
             data.bus_reactivepower_injection[ix] - data.bus_reactivepower_withdrawals[ix]
     end
     x0 = _calculate_x0(1,
-        data.bus_type,
+        data.bus_type[:, time_step],
         data.bus_angles,
         data.bus_magnitude,
         data.bus_activepower_injection,
@@ -112,9 +113,10 @@ function polar_pf!(
     Q_net::Vector{Float64},
     data::ACPowerFlowData,
 )
+    time_step = 1  # TODO placeholder time_step
     Yb = data.power_network_matrix.data
-    n_buses = length(data.bus_type)
-    for (ix, b) in enumerate(data.bus_type)
+    n_buses = first(size(data.bus_type))
+    for (ix, b) in enumerate(data.bus_type[:, time_step])
         if b == PSY.ACBusTypes.REF
             # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated
             P_net[ix] = X[2 * ix - 1] - data.bus_activepower_withdrawals[ix]

src/ac_power_flow_jacobian.jl

@@ -24,18 +24,19 @@ function PolarPowerFlowJacobian(data::ACPowerFlowData, x0::Vector{Float64})
 end
 
 function _create_jacobian_matrix_structure(data::ACPowerFlowData)
+    time_step = 1  # TODO placeholder time_step
     #Create Jacobian structure
     J0_I = Int32[]
     J0_J = Int32[]
     J0_V = Float64[]
 
-    for ix_f in eachindex(data.bus_type)
+    for ix_f in eachindex(data.bus_type[:, time_step])
         F_ix_f_r = 2 * ix_f - 1
         F_ix_f_i = 2 * ix_f
         for ix_t in data.neighbors[ix_f]
             X_ix_t_fst = 2 * ix_t - 1
             X_ix_t_snd = 2 * ix_t
-            nb = data.bus_type[ix_t]
+            nb = data.bus_type[ix_t, time_step]
             #Set to 0.0 only on connected buses
             if nb == PSY.ACBusTypes.REF
                 if ix_f == ix_t
@@ -92,17 +93,18 @@ function jsp!(
     ::Vector{Float64},
     data::ACPowerFlowData,
 )
+    time_step = 1  # TODO placeholder time_step
     Yb = data.power_network_matrix.data
     Vm = data.bus_magnitude
     θ = data.bus_angles
-    for (ix_f, b) in enumerate(data.bus_type)
+    for (ix_f, b) in enumerate(data.bus_type[:, time_step])
         F_ix_f_r = 2 * ix_f - 1
         F_ix_f_i = 2 * ix_f
 
         for ix_t in data.neighbors[ix_f]
             X_ix_t_fst = 2 * ix_t - 1
             X_ix_t_snd = 2 * ix_t
-            nb = data.bus_type[ix_t]
+            nb = data.bus_type[ix_t, time_step]
             if nb == PSY.ACBusTypes.REF
                 # State variables are Active and Reactive Power Generated
                 # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - x[2*i-1]

src/common.jl

@@ -118,11 +118,11 @@ function make_dc_powerflowdata(
     temp_bus_map,
     valid_ix,
 )
-    branch_types = Vector{DataType}(undef, length(branch_lookup))
+    branch_type = Vector{DataType}(undef, length(branch_lookup))
     for (ix, b) in enumerate(PNM.get_ac_branches(sys))
-        branch_types[ix] = typeof(b)
+        branch_type[ix] = typeof(b)
     end
-    bus_reactivepower_bounds = Vector{Vector{Float64}}()
+    bus_reactivepower_bounds = Vector{Vector{Float64}}(undef, n_buses)
     timestep_map = Dict(zip([i for i in 1:time_steps], timestep_names))
     neighbors = Vector{Set{Int}}()
     return make_powerflowdata(
@@ -135,8 +135,7 @@ function make_dc_powerflowdata(
         bus_lookup,
         branch_lookup,
         temp_bus_map,
-        branch_types,
-        bus_reactivepower_bounds,
+        branch_type,
         timestep_map,
         valid_ix,
         neighbors,
@@ -153,13 +152,11 @@ function make_powerflowdata(
     bus_lookup,
     branch_lookup,
     temp_bus_map,
-    branch_types,
-    bus_reactivepower_bounds,
+    branch_type,
     timestep_map,
     valid_ix,
     neighbors,
 )
-    # TODO: bus_type might need to also be a Matrix since the type can change for a particular scenario
     bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
     bus_magnitude = zeros(Float64, n_buses)
@@ -192,27 +189,43 @@ function make_powerflowdata(
         sys,
     )
 
-    # initialize data
-    init_1 = zeros(n_buses, time_steps)
-    init_2 = zeros(n_branches, time_steps)
+    # Shapes to reuse
+    zeros_bus_time = () -> zeros(n_buses, time_steps)
+    zeros_branch_time = () -> zeros(n_branches, time_steps)
 
-    # define fields as matrices whose number of columns is eqault to the number of time_steps
-    bus_activepower_injection_1 = deepcopy(init_1)
-    bus_reactivepower_injection_1 = deepcopy(init_1)
-    bus_activepower_withdrawals_1 = deepcopy(init_1)
-    bus_reactivepower_withdrawals_1 = deepcopy(init_1)
-    bus_magnitude_1 = deepcopy(init_1)
-    bus_angles_1 = deepcopy(init_1)
-    branch_flow_values_1 = deepcopy(init_2)
+    # Define fields as matrices whose number of columns is equal to the number of time_steps
+    bus_activepower_injection_1 = zeros_bus_time()
+    bus_reactivepower_injection_1 = zeros_bus_time()
+    bus_activepower_withdrawals_1 = zeros_bus_time()
+    bus_reactivepower_withdrawals_1 = zeros_bus_time()
+    bus_reactivepower_bounds_1 = Matrix{Vector{Float64}}(undef, n_buses, time_steps)
+    bus_magnitude_1 = zeros_bus_time()
+    bus_angles_1 = zeros_bus_time()
 
-    # initial values related to first timestep allocated in the first column
+    # Initial values related to first timestep allocated in the first column
     bus_activepower_injection_1[:, 1] .= bus_activepower_injection
     bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
     bus_activepower_withdrawals_1[:, 1] .= bus_activepower_withdrawals
     bus_reactivepower_withdrawals_1[:, 1] .= bus_reactivepower_withdrawals
     bus_magnitude_1[:, 1] .= bus_magnitude
     bus_angles_1[:, 1] .= bus_angles
-    branch_flow_values_1[:, 1] .= zeros(n_branches)
+
+    bus_reactivepower_bounds = Vector{Vector{Float64}}(undef, n_buses)
+    for i in 1:n_buses
+        bus_reactivepower_bounds[i] = [0.0, 0.0]
+    end
+    _get_reactive_power_bound!(bus_reactivepower_bounds, bus_lookup, sys)
+    bus_reactivepower_bounds_1[:, 1] .= bus_reactivepower_bounds
+
+    # Initial bus types are same for every time period
+    bus_type_1 = repeat(bus_type; outer = [1, time_steps])
+    @assert size(bus_type_1) == (n_buses, time_steps)
+
+    # Initial flows are all zero
+    branch_activepower_flow_from_to = zeros_branch_time()
+    branch_reactivepower_flow_from_to = zeros_branch_time()
+    branch_activepower_flow_to_from = zeros_branch_time()
+    branch_reactivepower_flow_to_from = zeros_branch_time()
 
     return PowerFlowData(
         bus_lookup,
@@ -221,12 +234,15 @@ function make_powerflowdata(
         bus_reactivepower_injection_1,
         bus_activepower_withdrawals_1,
         bus_reactivepower_withdrawals_1,
-        bus_reactivepower_bounds,
-        bus_type,
-        branch_types,
+        bus_reactivepower_bounds_1,
+        bus_type_1,
+        branch_type,
         bus_magnitude_1,
         bus_angles_1,
-        branch_flow_values_1,
+        branch_activepower_flow_from_to,
+        branch_reactivepower_flow_from_to,
+        branch_activepower_flow_to_from,
+        branch_reactivepower_flow_to_from,
         timestep_map,
         valid_ix,
         power_network_matrix,

src/nlsolve_ac_powerflow.jl

@@ -98,7 +98,11 @@ function solve_powerflow(
     return res.f_converged
 end
 
-function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64})
+function _check_q_limit_bounds!(
+    data::ACPowerFlowData,
+    zero::Vector{Float64},
+    time_step::Int,
+)
     bus_names = data.power_network_matrix.axes[1]
     within_limits = true
     for (ix, b) in enumerate(data.bus_type)
@@ -108,16 +112,18 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64})
             continue
         end
 
-        if Q_gen <= data.bus_reactivepower_bounds[ix][1]
+        if Q_gen <= data.bus_reactivepower_bounds[ix, time_step][1]
             @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
             within_limits = false
-            data.bus_type[ix] = PSY.ACBusTypes.PQ
-            data.bus_reactivepower_injection[ix] = data.bus_reactivepower_bounds[ix][1]
-        elseif Q_gen >= data.bus_reactivepower_bounds[ix][2]
+            data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ
+            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] = PSY.ACBusTypes.PQ
-            data.bus_reactivepower_injection[ix] = data.bus_reactivepower_bounds[ix][2]
+            data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ
+            data.bus_reactivepower_injection[ix, time_step] =
+                data.bus_reactivepower_bounds[ix, time_step][2]
         else
             @debug "Within Limits"
         end
@@ -128,13 +134,14 @@ end
 function _solve_powerflow!(
     data::ACPowerFlowData,
     check_reactive_power_limits;
+    time_step = 1,  # TODO placeholder time_step
     nlsolve_kwargs...,
 )
     if check_reactive_power_limits
         for _ in 1:MAX_REACTIVE_POWER_ITERATIONS
             res = _nlsolve_powerflow(data; nlsolve_kwargs...)
             if res.f_converged
-                if _check_q_limit_bounds!(data, res.zero)
+                if _check_q_limit_bounds!(data, res.zero, time_step)
                     return res
                 end
             else

src/post_processing.jl

@@ -514,19 +514,10 @@ function _allocate_results_data(
     Q_gen_vect::Vector{Float64},
     P_load_vect::Vector{Float64},
     Q_load_vect::Vector{Float64},
-    branch_flow_values::Vector{Float64})
-
-    # get flows on each line
-    P_from_to_vect = zeros(length(branches))
-    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)
-        P_from_to_vect[i] = branch_flow_values[i]
-        P_to_from_vect[i] = -branch_flow_values[i]
-    end
-
+    branch_activepower_flow_from_to::Vector{Float64},
+    branch_reactivepower_flow_from_to::Vector{Float64},
+    branch_activepower_flow_to_from::Vector{Float64},
+    branch_reactivepower_flow_to_from::Vector{Float64})
     bus_df = DataFrames.DataFrame(;
         bus_number = buses,
         Vm = bus_magnitude,
@@ -544,10 +535,10 @@ function _allocate_results_data(
         line_name = branches,
         bus_from = from_bus,
         bus_to = to_bus,
-        P_from_to = P_from_to_vect,
-        Q_from_to = Q_from_to_vect,
-        P_to_from = P_to_from_vect,
-        Q_to_from = Q_to_from_vect,
+        P_from_to = branch_activepower_flow_from_to,
+        Q_from_to = branch_reactivepower_flow_from_to,
+        P_to_from = branch_activepower_flow_to_from,
+        Q_to_from = branch_reactivepower_flow_to_from,
         P_losses = zeros(length(branches)),
         Q_losses = zeros(length(branches)),
     )
@@ -600,7 +591,10 @@ function write_results(
             data.bus_reactivepower_injection[:, i],
             data.bus_activepower_withdrawals[:, i],
             data.bus_reactivepower_withdrawals[:, i],
-            data.branch_flow_values[:, i],
+            data.branch_activepower_flow_from_to[:, i],
+            data.branch_reactivepower_flow_from_to[:, i],
+            data.branch_activepower_flow_to_from[:, i],
+            data.branch_reactivepower_flow_to_from[:, i],
         )
         result_dict[data.timestep_map[i]] = temp_dict
     end