integrate previous changes to powersystemdata

src/ac_power_flow.jl

@@ -55,16 +55,15 @@ function _calculate_x0(n::Int,
     return x0
 end
 
-function PolarPowerFlow(data::ACPowerFlowData)
-    time_step = 1  # TODO placeholder time_step
+function PolarPowerFlow(data::ACPowerFlowData; time_step::Int64 = 1)
     n_buses = first(size(data.bus_type))
     P_net = zeros(n_buses)
     Q_net = zeros(n_buses)
     for ix in 1:n_buses
         P_net[ix] =
-            data.bus_activepower_injection[ix] - data.bus_activepower_withdrawals[ix]
+            data.bus_activepower_injection[ix, time_step] - data.bus_activepower_withdrawals[ix, time_step]
         Q_net[ix] =
-            data.bus_reactivepower_injection[ix] - data.bus_reactivepower_withdrawals[ix]
+            data.bus_reactivepower_injection[ix, time_step] - data.bus_reactivepower_withdrawals[ix, time_step]
     end
     x0 = _calculate_x0(1,
         data.bus_type[:, time_step],

src/common.jl

@@ -199,7 +199,7 @@ function make_powerflowdata(
     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_stepsm)
+    bus_reactivepower_bounds_1 = Matrix{Vector{Float64}}(undef, n_buses, time_steps)
     bus_magnitude_1 = ones_bus_time()
     bus_angles_1 = zeros_bus_time()
 

src/newton_ac_powerflow.jl

@@ -1,4 +1,3 @@
-const _SOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity])
 """
 Solves a the power flow into the system and writes the solution into the relevant structs.
 Updates generators active and reactive power setpoints and branches active and reactive
@@ -34,6 +33,7 @@ solve_ac_powerflow!(sys, method=:newton)
 function solve_powerflow!(
     pf::ACPowerFlow{<:ACPowerFlowSolverType},
     system::PSY.System;
+    time_step::Int64=1,
     kwargs...,
 )
     #Save per-unit flag
@@ -47,8 +47,7 @@ function solve_powerflow!(
         check_connectivity = get(kwargs, :check_connectivity, true),
     )
 
-    solver_kwargs = filter(p -> !(p.first in _SOLVE_AC_POWERFLOW_KWARGS), kwargs)
-    converged, V, Sbus_result = _ac_powereflow(data, pf; solver_kwargs...)
+    converged, V, Sbus_result = _ac_powereflow(data, pf; time_step=time_step, kwargs...)
     x = _calc_x(data, V, Sbus_result)
 
     if converged
@@ -124,8 +123,8 @@ function solve_powerflow(
     sorted_time_steps = sort(collect(keys(data.timestep_map)))
     # preallocate results
     ts_converged = zeros(Bool, 1, length(sorted_time_steps))
-    ts_V = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps))
-    ts_S = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps))
+    ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
+    ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
 
     results = Dict()
 
@@ -176,8 +175,8 @@ function solve_powerflow!(
     sorted_time_steps = sort(collect(keys(data.timestep_map)))
     # preallocate results
     ts_converged = zeros(Bool, 1, length(sorted_time_steps))
-    ts_V = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps))
-    ts_S = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps))
+    ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
+    ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
 
     Yft = data.power_network_matrix.data_ft
     Ytf = data.power_network_matrix.data_tf
@@ -192,9 +191,9 @@ function solve_powerflow!(
         ts_V[:, t] .= V
         ts_S[:, t] .= Sbus_result
 
-        ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type)
-        pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type)
-        pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type)
+        ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, t])
+        pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, t])
+        pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, t])
 
         # temporary implementation that will need to be improved:
         if converged
@@ -222,7 +221,13 @@ function solve_powerflow!(
     end
 
     # write branch flows
-    data.branch_flow_values .= real.(ts_V[fb,:] .* conj.(Yft * ts_V))
+    Sft = ts_V[fb,:] .* conj.(Yft * ts_V)
+    Stf = ts_V[tb,:] .* conj.(Ytf * ts_V)
+
+    data.branch_activepower_flow_from_to .= real.(Sft)
+    data.branch_reactivepower_flow_from_to .= imag.(Sft)
+    data.branch_activepower_flow_to_from .= real.(Stf)
+    data.branch_reactivepower_flow_to_from .= imag.(Stf)
 
     # todo:
     # return df_results
@@ -243,7 +248,7 @@ function _ac_powereflow(
         converged, V, Sbus_result =
             _newton_powerflow(pf, data; time_step = time_step, kwargs...)
         if !converged || !check_reactive_power_limits ||
-           _check_q_limit_bounds!(data, Sbus_result)
+           _check_q_limit_bounds!(data, Sbus_result, time_step)
             return (converged, V, Sbus_result)
         end
     end
@@ -259,8 +264,8 @@ function _check_q_limit_bounds!(
 )
     bus_names = data.power_network_matrix.axes[1]
     within_limits = true
-    for (ix, b) in enumerate(data.bus_type)
-        if b == PSY.ACBusTypes.PV
+    for (ix, bt) in enumerate(data.bus_type[:, time_step])
+        if bt == PSY.ACBusTypes.PV
             Q_gen = imag(Sbus_result[ix])
         else
             continue
@@ -271,10 +276,10 @@ function _check_q_limit_bounds!(
             within_limits = false
             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][2]
+        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[ixm, time_step] = PSY.ACBusTypes.PQ
+            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"
@@ -469,23 +474,24 @@ end
 function _calc_x(
     data::ACPowerFlowData,
     V::Vector{Complex{Float64}},
-    Sbus_result::Vector{Complex{Float64}},
+    Sbus_result::Vector{Complex{Float64}};
+    time_step::Int64 = 1,
 )
     Vm = abs.(V)
     Va = angle.(V)
     n_buses = length(V)
     # mock the expected x format, where the values depend on the type of the bus:
     x = zeros(Float64, 2 * n_buses)
-    for (ix, b) in enumerate(data.bus_type)
-        if b == PSY.ACBusTypes.REF
+    for (ix, bt) in enumerate(data.bus_type[:, time_step])
+        if bt == PSY.ACBusTypes.REF
             # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated
             x[2 * ix - 1] = real(Sbus_result[ix]) + data.bus_activepower_withdrawals[ix]
             x[2 * ix] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix]
-        elseif b == PSY.ACBusTypes.PV
+        elseif bt == PSY.ACBusTypes.PV
             # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle
             x[2 * ix - 1] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix]
             x[2 * ix] = Va[ix]
-        elseif b == PSY.ACBusTypes.PQ
+        elseif bt == PSY.ACBusTypes.PQ
             # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle
             x[2 * ix - 1] = Vm[ix]
             x[2 * ix] = Va[ix]
@@ -519,16 +525,15 @@ function _newton_powerflow(
     Ybus = data.power_network_matrix.data
 
     # Find indices for each bus type
-    ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type)
-    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type)
-    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type)
+    ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, time_step])
+    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step])
+    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step])
     pvpq = [pv; pq]
 
     # nref = length(ref)
     npv = length(pv)
     npq = length(pq)
     npvpq = npv + npq
-    n_buses = length(data.bus_type)
 
     Vm = data.bus_magnitude[:, time_step]
     # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this:

src/nlsolve_powerflow.jl

@@ -1,4 +1,4 @@
-
+const _NLSOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity])
 function _newton_powerflow(
     pf::ACPowerFlow{NLSolveACPowerFlow},
     data::ACPowerFlowData;
@@ -13,11 +13,13 @@ function _newton_powerflow(
         )
     end
 
+    nlsolve_solver_kwargs = filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs)
+
     pf = PolarPowerFlow(data)
     J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
 
     df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv)
-    res = NLsolve.nlsolve(df, pf.x0; nlsolve_kwargs...)
+    res = NLsolve.nlsolve(df, pf.x0; nlsolve_solver_kwargs...)
     if !res.f_converged
         @error(
             "The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))"

src/post_processing.jl

@@ -461,6 +461,7 @@ function write_powerflow_solution!(
         system_bustype = PSY.get_bustype(bus)
         data_bustype = data.bus_type[ix]
         (system_bustype == data_bustype) && continue
+        println("sys bustype: $system_bustype, data bustype: $data_bustype")
         @assert system_bustype == PSY.ACBusTypes.PV
         @assert data_bustype == PSY.ACBusTypes.PQ
         @debug "Updating bus $(PSY.get_name(bus)) reactive power and type to PQ due to check_reactive_power_limits"
@@ -620,7 +621,8 @@ function write_results(
     ::ACPowerFlow{<:ACPowerFlowSolverType},
     sys::PSY.System,
     data::ACPowerFlowData,
-    result::Union{Vector{Float64}, Matrix{Float64}}, #todo
+    result::Union{Vector{Float64}, Matrix{Float64}};
+    time_step::Int64 = 1,
 )
     @info("Voltages are exported in pu. Powers are exported in MW/MVAr.")
     buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x))
@@ -634,27 +636,27 @@ function write_results(
     P_load_vect = fill(0.0, N_BUS)
     Q_load_vect = fill(0.0, N_BUS)
 
-    for (ix, bustype) in enumerate(data.bus_type)
-        P_load_vect[ix] = data.bus_activepower_withdrawals[ix] * sys_basepower
-        Q_load_vect[ix] = data.bus_reactivepower_withdrawals[ix] * sys_basepower
+    for (ix, bustype) in enumerate(data.bus_type[:, time_step])
+        P_load_vect[ix] = data.bus_activepower_withdrawals[ix, time_step] * sys_basepower
+        Q_load_vect[ix] = data.bus_reactivepower_withdrawals[ix, time_step] * sys_basepower
         P_admittance, Q_admittance = _get_fixed_admittance_power(sys, buses[ix], result, ix)
         P_load_vect[ix] += P_admittance * sys_basepower
         Q_load_vect[ix] += Q_admittance * sys_basepower
         if bustype == PSY.ACBusTypes.REF
-            Vm_vect[ix] = data.bus_magnitude[ix]
-            θ_vect[ix] = data.bus_angles[ix]
+            Vm_vect[ix] = data.bus_magnitude[ix, time_step]
+            θ_vect[ix] = data.bus_angles[ix, time_step]
             P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix] * sys_basepower
         elseif bustype == PSY.ACBusTypes.PV
-            Vm_vect[ix] = data.bus_magnitude[ix]
+            Vm_vect[ix] = data.bus_magnitude[ix, time_step]
             θ_vect[ix] = result[2 * ix]
-            P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower
+            P_gen_vect[ix] = data.bus_activepower_injection[ix, time_step] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
         elseif bustype == PSY.ACBusTypes.PQ
             Vm_vect[ix] = result[2 * ix - 1]
             θ_vect[ix] = result[2 * ix]
-            P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower
-            Q_gen_vect[ix] = data.bus_reactivepower_injection[ix] * sys_basepower
+            P_gen_vect[ix] = data.bus_activepower_injection[ix, time_step] * sys_basepower
+            Q_gen_vect[ix] = data.bus_reactivepower_injection[ix, time_step] * sys_basepower
         end
     end
 

test/test_newton_ac_powerflow.jl

@@ -45,25 +45,25 @@
     x1 = PowerFlows._calc_x(data, V1, S1)
     @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6
 
-    # Test that solve_ac_powerflow! succeeds
+    # Test that solve_powerflow! succeeds
     solved1 = deepcopy(sys)
     @test solve_powerflow!(pf, solved1)
 
     # Test that passing check_reactive_power_limits=false is the default and violates limits
     solved2 = deepcopy(sys)
-    @test solve_ac_powerflow!(solved2; check_reactive_power_limits = false)
+    @test solve_powerflow!(pf, solved2; check_reactive_power_limits = false)
     @test IS.compare_values(solved1, solved2)
     @test get_reactive_power(get_component(ThermalStandard, solved2, "Bus8")) >
           get_reactive_power_limits(get_component(ThermalStandard, solved2, "Bus8")).max
 
     # Test that passing check_reactive_power_limits=true fixes that
     solved3 = deepcopy(sys)
-    @test solve_ac_powerflow!(solved3; check_reactive_power_limits = true)
+    @test solve_powerflow!(pf, solved3; check_reactive_power_limits = true)
     @test get_reactive_power(get_component(ThermalStandard, solved3, "Bus8")) <=
           get_reactive_power_limits(get_component(ThermalStandard, solved3, "Bus8")).max
 
     # Test Newton method
-    @test solve_ac_powerflow!(deepcopy(sys); method = :newton)
+    @test solve_powerflow!(pf, deepcopy(sys); method = :newton)
 
     # Test enforcing the reactive power limits in closer detail
     set_reactive_power!(get_component(PowerLoad, sys, "Bus4"), 0.0)
@@ -505,10 +505,12 @@ end
         pf_default,
         sys;
         check_connectivity = true)
+
+    time_step = 1
 
-    ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type)
-    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type)
-    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type)
+    ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, time_step])
+    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step])
+    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step])
     pvpq = [pv; pq]
 
     npvpq = length(pvpq)

test/test_utils/legacy_pf.jl

@@ -47,8 +47,8 @@ function _newton_powerflow(
 
     # Find indices for each bus type
     #ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type)
-    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type)
-    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type)
+    pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step])
+    pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step])
     pvpq = [pv; pq]
 
     #nref = length(ref)