add field "converged" to PowerFlowData, return NaN when not converged

src/PowerFlowData.jl

@@ -1,3 +1,8 @@
+Base.@kwdef mutable struct SolverData
+    J::Union{SparseMatrixCSC{Float64, Int}, Nothing} = nothing
+    dSbus_dV_ref::Union{Vector{Float64}, Nothing} = nothing
+end
+
 abstract type PowerFlowContainer end
 
 """
@@ -102,6 +107,8 @@ struct PowerFlowData{
     power_network_matrix::M
     aux_network_matrix::N
     neighbors::Vector{Set{Int}}
+    converged::Vector{Bool}
+    solver_data::Vector{SolverData}
 end
 
 get_bus_lookup(pfd::PowerFlowData) = pfd.bus_lookup
@@ -129,6 +136,8 @@ get_power_network_matrix(pfd::PowerFlowData) = pfd.power_network_matrix
 get_aux_network_matrix(pfd::PowerFlowData) = pfd.aux_network_matrix
 get_neighbor(pfd::PowerFlowData) = pfd.neighbors
 supports_multi_period(::PowerFlowData) = true
+get_converged(pfd::PowerFlowData) = pfd.converged
+get_solver_data(pfd::PowerFlowData) = pfd.solver_data
 
 function clear_injection_data!(pfd::PowerFlowData)
     pfd.bus_activepower_injection .= 0.0
@@ -220,10 +229,11 @@ function PowerFlowData(
         branch_types[ix] = typeof(b)
     end
 
-    timestep_map = Dict(1 => "1")
     valid_ix = setdiff(1:n_buses, ref_bus_positions)
     neighbors = _calculate_neighbors(power_network_matrix)
     aux_network_matrix = nothing
+    converged = fill(false, time_steps)
+    solver_data = [SolverData() for _ in 1:time_steps]
 
     return make_powerflowdata(
         sys,
@@ -239,6 +249,8 @@ function PowerFlowData(
         timestep_map,
         valid_ix,
         neighbors,
+        converged,
+        solver_data,
     )
 end
 
@@ -295,6 +307,8 @@ function PowerFlowData(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.ACBus, sys)
     )
     valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions)
+    converged = fill(false, time_steps)
+    solver_data = [SolverData() for _ in 1:time_steps]
     return make_dc_powerflowdata(
         sys,
         time_steps,
@@ -307,6 +321,8 @@ function PowerFlowData(
         branch_lookup,
         temp_bus_map,
         valid_ix,
+        converged,
+        solver_data,
     )
 end
 
@@ -364,6 +380,8 @@ function PowerFlowData(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
     )
     valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions)
+    converged = fill(false, time_steps)
+    solver_data = [SolverData() for _ in 1:time_steps]
     return make_dc_powerflowdata(
         sys,
         time_steps,
@@ -376,6 +394,8 @@ function PowerFlowData(
         branch_lookup,
         temp_bus_map,
         valid_ix,
+        converged,
+        solver_data,
     )
 end
 
@@ -432,6 +452,8 @@ function PowerFlowData(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
     )
     valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions)
+    converged = fill(false, time_steps)
+    solver_data = [SolverData() for _ in 1:time_steps]
     return make_dc_powerflowdata(
         sys,
         time_steps,
@@ -444,6 +466,8 @@ function PowerFlowData(
         branch_lookup,
         temp_bus_map,
         valid_ix,
+        converged,
+        solver_data,
     )
 end
 

src/PowerFlows.jl

@@ -16,6 +16,7 @@ export PSSEExporter
 export update_exporter!
 export write_export
 export get_psse_export_paths
+export penalty_factors
 
 import Logging
 import DataFrames

src/common.jl

@@ -149,6 +149,8 @@ function make_dc_powerflowdata(
     branch_lookup,
     temp_bus_map,
     valid_ix,
+    converged,
+    solver_data,
 )
     branch_type = Vector{DataType}(undef, length(branch_lookup))
     for (ix, b) in enumerate(PNM.get_ac_branches(sys))
@@ -171,6 +173,8 @@ function make_dc_powerflowdata(
         timestep_map,
         valid_ix,
         neighbors,
+        converged,
+        solver_data,
     )
 end
 
@@ -188,6 +192,8 @@ function make_powerflowdata(
     timestep_map,
     valid_ix,
     neighbors,
+    converged,
+    solver_data,
 )
     bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
@@ -281,5 +287,7 @@ function make_powerflowdata(
         power_network_matrix,
         aux_network_matrix,
         neighbors,
+        converged,
+        solver_data,
     )
 end

src/newton_ac_powerflow.jl

@@ -175,12 +175,15 @@ function solve_powerflow!(
 )
     pf = ACPowerFlow()  # todo: somehow store in data which PF to use (see issue #50)
 
-    sorted_time_steps = sort(collect(keys(data.timestep_map)))
+    sorted_time_steps = get(kwargs, :time_steps, sort(collect(keys(data.timestep_map))))
     # preallocate results
-    ts_converged = zeros(Bool, 1, length(sorted_time_steps))
+    ts_converged = fill(false, 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))
 
+    # TODO If anything in the grid topology changes, 
+    #  e.g. tap positions of transformers or in service 
+    #  status of branches, Yft and Ytf must be updated!
     Yft = data.power_network_matrix.data_ft
     Ytf = data.power_network_matrix.data_tf
 
@@ -190,25 +193,25 @@ function solve_powerflow!(
     for t in sorted_time_steps
         converged, V, Sbus_result =
             _ac_powereflow(data, pf; time_step = t, kwargs...)
-        ts_converged[1, t] = converged
+        ts_converged[t] = converged
         ts_V[:, t] .= V
         ts_S[:, t] .= Sbus_result
 
-        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
+            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:
             # write results for REF
             data.bus_activepower_injection[ref, t] .=
                 real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t]
@@ -219,21 +222,22 @@ function solve_powerflow!(
                 imag.(Sbus_result[pv]) .+ data.bus_reactivepower_withdrawals[pv, t]
             # results for PQ buses do not need to be updated -> already consistent with inputs
 
-            # write bus bus_types
-            # todo
-
             # write voltage results
             data.bus_magnitude[pq, t] .= abs.(V[pq])
             data.bus_angles[pq, t] .= angle.(V[pq])
             data.bus_angles[pv, t] .= angle.(V[pv])
-
         else
-            # todo
-            1 + 2
+            data.bus_activepower_injection[:, t] .= NaN64
+            data.bus_activepower_withdrawals[:, t] .= NaN64
+            data.bus_reactivepower_injection[:, t] .= NaN64
+            data.bus_reactivepower_withdrawals[:, t] .= NaN64
+            data.bus_magnitude[:, t] .= NaN64
+            data.bus_angles[:, t] .= NaN64
         end
     end
 
     # write branch flows
+    # TODO if Yft, Ytf change between time steps, this must be moved inside the loop!
     Sft = ts_V[fb, :] .* conj.(Yft * ts_V)
     Stf = ts_V[tb, :] .* conj.(Ytf * ts_V)
 
@@ -242,8 +246,7 @@ function solve_powerflow!(
     data.branch_activepower_flow_to_from .= real.(Stf)
     data.branch_reactivepower_flow_to_from .= imag.(Stf)
 
-    # todo:
-    # return df_results
+    data.converged .= ts_converged
 
     return
 end
@@ -537,6 +540,8 @@ function _newton_powerflow(
     tol = get(kwargs, :tol, DEFAULT_NR_TOL)
     i = 0
 
+    solver_data = data.solver_data[time_step]
+
     Ybus = data.power_network_matrix.data
 
     # Find indices for each bus type
@@ -713,8 +718,13 @@ function _newton_powerflow(
     end
 
     if !converged
+        V .*= NaN64
+        Sbus_result .*= NaN64
         @error("The powerflow solver with KLU did not converge after $i iterations")
     else
+        solver_data.J = J
+        solver_data.dSbus_dV_ref =
+            [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pq]))]
         @info("The powerflow solver with KLU converged after $i iterations")
     end
 

src/nlsolve_powerflow.jl

@@ -23,12 +23,15 @@ function _newton_powerflow(
     df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv)
     res = NLsolve.nlsolve(df, pf.x0; nlsolve_solver_kwargs...)
     if !res.f_converged
+        V = fill(NaN64 + NaN64 * im, length(res.zero) ÷ 2)
+        Sbus_result = fill(NaN64 + NaN64 * im, length(res.zero) ÷ 2)
         @error(
             "The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))"
         )
+    else
+        V = _calc_V(data, res.zero; time_step = time_step)
+        Sbus_result = V .* conj(data.power_network_matrix.data * V)
     end
-    V = _calc_V(data, res.zero; time_step = time_step)
-    Sbus_result = V .* conj(data.power_network_matrix.data * V)
     return (res.f_converged, V, Sbus_result)
 end
 

src/post_processing.jl

@@ -757,3 +757,10 @@ function set_system_time_step(sys::PSY.System, data::PowerFlowData; time_step =
         end
     end
 end
+
+# work in progress - quick but not optimized function for POC
+function penalty_factors(solver_data::SolverData)
+    J_t = transpose(solver_data.J)
+    f = J_t \ solver_data.dSbus_dV_ref
+    return f
+end

src/solve_dc_powerflow.jl

@@ -53,6 +53,7 @@ function solve_powerflow!(
     # evaluate bus angles
     p_inj = power_injection[data.valid_ix, :]
     data.bus_angles[data.valid_ix, :] .= data.aux_network_matrix.K \ p_inj
+    data.converged .= true
     return
 end
 
@@ -77,6 +78,7 @@ function solve_powerflow!(
         p_inj = power_injection[data.valid_ix, i]
         data.bus_angles[data.valid_ix, i] .= data.aux_network_matrix.K \ p_inj
     end
+    data.converged .= true
     return
 end
 
@@ -100,6 +102,7 @@ function solve_powerflow!(
         data.power_network_matrix.K \ @view power_injection[data.valid_ix, :]
     data.branch_activepower_flow_from_to .= data.aux_network_matrix.data' * data.bus_angles
     data.branch_activepower_flow_to_from .= -data.branch_activepower_flow_from_to
+    data.converged .= true
     return
 end
 

test/test_utils/legacy_pf.jl

@@ -111,10 +111,15 @@ function _newton_powerflow(
     end
 
     if !converged
+        V .*= NaN64
+        Sbus_result = fill(NaN64 + NaN64 * im, length(V))
         @error("The powerflow solver with KLU did not converge after $i iterations")
     else
+        Sbus_result = V .* conj(Ybus * V)
+        solver_data.J = J
+        solver_data.dSbus_dV_ref =
+            [vec(real.(dSbus_dVa[ref, :][:, pvpq])); vec(real.(dSbus_dVm[ref, :][:, pq]))]
         @info("The powerflow solver with KLU converged after $i iterations")
     end
-    Sbus_result = V .* conj(Ybus * V)
     return (converged, V, Sbus_result)
 end