Feature/nr_pf_solver

Project.toml

@@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -22,6 +23,7 @@ DataStructures = "0.18"
 Dates = "1"
 InfrastructureSystems = "2"
 JSON3 = "1"
+KLU = "^0.6"
 LinearAlgebra = "1"
 Logging = "1"
 NLsolve = "4"

src/PowerFlowData.jl

@@ -180,7 +180,7 @@ NOTE: use it for AC power flow computations.
 WARNING: functions for the evaluation of the multi-period AC PF still to be implemented.
 """
 function PowerFlowData(
-    ::ACPowerFlow,
+    ::ACPowerFlow{<:ACPowerFlowSolverType},
     sys::PSY.System;
     time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
@@ -196,6 +196,8 @@ function PowerFlowData(
         end
     end
 
+    timestep_map = Dict(zip([i for i in 1:time_steps], timestep_names))
+
     # get data for calculations
     power_network_matrix = PNM.Ybus(sys; check_connectivity = check_connectivity)
 
@@ -457,7 +459,11 @@ Create an appropriate `PowerFlowContainer` for the given `PowerFlowEvaluationMod
 """
 function make_power_flow_container end
 
-make_power_flow_container(pfem::ACPowerFlow, sys::PSY.System; kwargs...) =
+make_power_flow_container(
+    pfem::ACPowerFlow{<:ACPowerFlowSolverType},
+    sys::PSY.System;
+    kwargs...,
+) =
     PowerFlowData(pfem, sys; kwargs...)
 
 make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) =

src/PowerFlows.jl

@@ -1,10 +1,13 @@
 module PowerFlows
 
 export solve_powerflow
-export solve_ac_powerflow!
+export solve_powerflow!
 export PowerFlowData
 export DCPowerFlow
+export NLSolveACPowerFlow
+export KLUACPowerFlow
 export ACPowerFlow
+export ACPowerFlowSolverType
 export PTDFDCPowerFlow
 export vPTDFDCPowerFlow
 export PSSEExportPowerFlow
@@ -20,10 +23,11 @@ import PowerSystems
 import PowerSystems: System, with_units_base
 import LinearAlgebra
 import NLsolve
+import KLU
 import SparseArrays
 import InfrastructureSystems
 import PowerNetworkMatrices
-import SparseArrays: SparseMatrixCSC
+import SparseArrays: SparseMatrixCSC, sparse
 import JSON3
 import DataStructures: OrderedDict
 import Dates
@@ -40,6 +44,7 @@ include("psse_export.jl")
 include("solve_dc_powerflow.jl")
 include("ac_power_flow.jl")
 include("ac_power_flow_jacobian.jl")
-include("nlsolve_ac_powerflow.jl")
+include("newton_ac_powerflow.jl")
+include("nlsolve_powerflow.jl")
 include("post_processing.jl")
 end

src/ac_power_flow.jl

@@ -55,16 +55,17 @@ 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

@@ -159,7 +159,7 @@ function make_powerflowdata(
 )
     bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
-    bus_magnitude = zeros(Float64, n_buses)
+    bus_magnitude = ones(Float64, n_buses)
 
     _initialize_bus_data!(
         bus_type,
@@ -190,16 +190,17 @@ function make_powerflowdata(
     )
 
     # Shapes to reuse
-    zeros_bus_time = () -> zeros(n_buses, time_steps)
-    zeros_branch_time = () -> zeros(n_branches, time_steps)
+    zeros_bus_time = () -> zeros(Float64, n_buses, time_steps)
+    ones_bus_time = () -> ones(Float64, n_buses, time_steps)
+    zeros_branch_time = () -> zeros(Float64, n_branches, time_steps)
 
     # 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_magnitude_1 = ones_bus_time()
     bus_angles_1 = zeros_bus_time()
 
     # Initial values related to first timestep allocated in the first column

src/definitions.jl

@@ -10,4 +10,7 @@ const DEFAULT_MAX_REDISTRIBUTION_ITERATIONS = 10
 
 const ISAPPROX_ZERO_TOLERANCE = 1e-6
 
+const DEFAULT_NR_MAX_ITER::Int64 = 30   # default maxIter for the NR power flow
+const DEFAULT_NR_TOL::Float64 = 1e-9 # default tolerance for the NR power flow
+
 const AC_PF_KW = []