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.0"
 LinearAlgebra = "1"
 Logging = "1"
 NLsolve = "4"

src/PowerFlowData.jl

@@ -157,7 +157,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,
+    ::Union{NLSolveACPowerFlow, KLUACPowerFlow},
     sys::PSY.System;
     time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
@@ -440,7 +440,7 @@ 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::NLSolveACPowerFlow, sys::PSY.System; kwargs...) =
     PowerFlowData(pfem, sys; kwargs...)
 
 make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) =

src/PowerFlows.jl

@@ -4,7 +4,8 @@ export solve_powerflow
 export solve_ac_powerflow!
 export PowerFlowData
 export DCPowerFlow
-export ACPowerFlow
+export NLSolveACPowerFlow
+export KLUACPowerFlow
 export PTDFDCPowerFlow
 export vPTDFDCPowerFlow
 export PSSEExportPowerFlow
@@ -20,6 +21,7 @@ import PowerSystems
 import PowerSystems: System
 import LinearAlgebra
 import NLsolve
+import KLU
 import SparseArrays
 import InfrastructureSystems
 import PowerNetworkMatrices

src/nlsolve_ac_powerflow.jl

@@ -30,29 +30,33 @@ solve_ac_powerflow!(sys)
 solve_ac_powerflow!(sys, method=:newton)
 ```
 """
-function solve_ac_powerflow!(system::PSY.System; kwargs...)
+function solve_ac_powerflow!(
+    pf::Union{NLSolveACPowerFlow, KLUACPowerFlow},
+    system::PSY.System;
+    kwargs...,
+)
     #Save per-unit flag
     settings_unit_cache = deepcopy(system.units_settings.unit_system)
     #Work in System per unit
     PSY.set_units_base_system!(system, "SYSTEM_BASE")
     check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false)
     data = PowerFlowData(
-        ACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits),
+        NLSolveACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits),
         system;
         check_connectivity = get(kwargs, :check_connectivity, true),
     )
     max_iterations = DEFAULT_MAX_REDISTRIBUTION_ITERATIONS
-    res = _solve_powerflow!(data, check_reactive_power_limits; kwargs...)
-    if res.f_converged
-        write_powerflow_solution!(system, res.zero, max_iterations)
+    converged, x = _solve_powerflow!(pf, data, check_reactive_power_limits; kwargs...)
+    if converged
+        write_powerflow_solution!(system, x, max_iterations)
         @info("PowerFlow solve converged, the results have been stored in the system")
         #Restore original per unit base
         PSY.set_units_base_system!(system, settings_unit_cache)
-        return res.f_converged
+        return converged
     end
-    @error("The powerflow solver returned convergence = $(res.f_converged)")
+    @error("The powerflow solver returned convergence = $(converged)")
     PSY.set_units_base_system!(system, settings_unit_cache)
-    return res.f_converged
+    return converged
 end
 
 """
@@ -68,7 +72,7 @@ res = solve_powerflow(sys, method=:newton)
 ```
 """
 function solve_powerflow(
-    pf::ACPowerFlow,
+    pf::Union{NLSolveACPowerFlow, KLUACPowerFlow},
     system::PSY.System;
     kwargs...,
 )
@@ -82,18 +86,18 @@ function solve_powerflow(
         check_connectivity = get(kwargs, :check_connectivity, true),
     )
 
-    res = _solve_powerflow!(data, pf.check_reactive_power_limits; kwargs...)
+    converged, x = _solve_powerflow!(pf, data, pf.check_reactive_power_limits; kwargs...)
 
-    if res.f_converged
+    if converged
         @info("PowerFlow solve converged, the results are exported in DataFrames")
-        df_results = write_results(pf, system, data, res.zero)
+        df_results = write_results(pf, system, data, x)
         #Restore original per unit base
         PSY.set_units_base_system!(system, settings_unit_cache)
         return df_results
     end
-    @error("The powerflow solver returned convergence = $(res.f_converged)")
+    @error("The powerflow solver returned convergence = $(converged)")
     PSY.set_units_base_system!(system, settings_unit_cache)
-    return res.f_converged
+    return converged
 end
 
 function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64})
@@ -124,27 +128,32 @@ function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64})
 end
 
 function _solve_powerflow!(
+    pf::Union{NLSolveACPowerFlow, KLUACPowerFlow},
     data::ACPowerFlowData,
     check_reactive_power_limits;
     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)
-                    return res
+            converged, x = _nlsolve_powerflow(pf, data; nlsolve_kwargs...)
+            if converged
+                if _check_q_limit_bounds!(data, x)
+                    return converged, x
                 end
             else
-                return res
+                return converged, x
             end
         end
     else
-        return _nlsolve_powerflow(data; nlsolve_kwargs...)
+        return _nlsolve_powerflow(pf, data; nlsolve_kwargs...)
     end
 end
 
-function _nlsolve_powerflow(data::ACPowerFlowData; nlsolve_kwargs...)
+function _nlsolve_powerflow(
+    pf::NLSolveACPowerFlow,
+    data::ACPowerFlowData;
+    nlsolve_kwargs...,
+)
     pf = PolarPowerFlow(data)
     J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
 
@@ -153,5 +162,100 @@ function _nlsolve_powerflow(data::ACPowerFlowData; nlsolve_kwargs...)
     if !res.f_converged
         @error("The powerflow solver returned convergence = $(res.f_converged)")
     end
-    return res
+    return res.f_converged, res.zero
+end
+
+function _nlsolve_powerflow(pf::KLUACPowerFlow, data::ACPowerFlowData; nlsolve_kwargs...)
+    pf = PolarPowerFlow(data)
+    #J_function = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
+
+    maxIter = 30  # TODO
+    tol = 1e-6    # TODO
+    i = 0
+
+    Vm = data.bus_magnitude[:]
+    Va = data.bus_angles[:]
+    V = Vm .* exp.(1im * Va)
+
+    # 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)
+
+    Ybus = pf.data.power_network_matrix.data
+
+    Sbus =
+        data.bus_activepower_injection[:] - data.bus_activepower_withdrawals[:] +
+        1im * (data.bus_reactivepower_injection[:] - data.bus_reactivepower_withdrawals[:])
+
+    mis = V .* conj(Ybus * V) - Sbus
+    F = [real(mis[[pv; pq]]); imag(mis[pq])]
+
+    # nref = length(ref)
+    npv = length(pv)
+    npq = length(pq)
+
+    converged = (npv + npq) == 0  # if only ref buses present, we do not need to enter the loop
+
+    while i < maxIter && !converged
+        i += 1
+        diagV = LinearAlgebra.Diagonal(V)
+        diagIbus = LinearAlgebra.Diagonal(Ybus * V)
+        diagVnorm = LinearAlgebra.Diagonal(V ./ abs.(V))
+        dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm
+        dSbus_dVa = 1im * diagV * conj(diagIbus - Ybus * diagV)
+
+        j11 = real(dSbus_dVa[[pv; pq], [pv; pq]])
+        j12 = real(dSbus_dVm[[pv; pq], pq])
+        j21 = imag(dSbus_dVa[pq, [pv; pq]])
+        j22 = imag(dSbus_dVm[pq, pq])
+        J = [j11 j12; j21 j22]
+
+        factor_J = KLU.klu(J)
+        dx = -(factor_J \ F)
+
+        #J_function(J_function.Jv, x)
+        #dx = - (J_function.Jv \ F)
+
+        Va[pv] .+= dx[1:npv]
+        Va[pq] .+= dx[(npv + 1):(npv + npq)]
+        Vm[pq] .+= dx[(npv + npq + 1):(npv + 2 * npq)]
+
+        V = Vm .* exp.(1im * Va)
+
+        Vm = abs.(V)
+        Va = angle.(V)
+
+        mis = V .* conj(Ybus * V) - Sbus
+        F = [real(mis[[pv; pq]]); imag(mis[pq])]
+        converged = LinearAlgebra.norm(F, Inf) < tol
+    end
+
+    if !converged
+        @error("The powerflow solver with KLU did not converge after $i iterations")
+    else
+        @debug("The powerflow solver with KLU converged after $i iterations")
+    end
+
+    # mock the expected x format, where the values depend on the type of the bus:
+    n_buses = length(data.bus_type)
+    x = Float64[0.0 for _ in 1:(2 * n_buses)]
+    Sbus_result = V .* conj(Ybus * V)
+    for (ix, b) in enumerate(data.bus_type)
+        if b == 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
+            # 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
+            # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle
+            x[2 * ix - 1] = Vm[ix]
+            x[2 * ix] = Va[ix]
+        end
+    end
+
+    return converged, x
 end

src/post_processing.jl

@@ -608,7 +608,7 @@ dictionary will therefore feature just one key linked to one DataFrame.
         vector containing the reults for one single time-period.
 """
 function write_results(
-    ::ACPowerFlow,
+    ::Union{NLSolveACPowerFlow, KLUACPowerFlow},
     sys::PSY.System,
     data::ACPowerFlowData,
     result::Vector{Float64},

src/powerflow_types.jl

@@ -1,6 +1,10 @@
 abstract type PowerFlowEvaluationModel end
 
-Base.@kwdef struct ACPowerFlow <: PowerFlowEvaluationModel
+Base.@kwdef struct NLSolveACPowerFlow <: PowerFlowEvaluationModel
+    check_reactive_power_limits::Bool = false
+end
+
+Base.@kwdef struct KLUACPowerFlow <: PowerFlowEvaluationModel
     check_reactive_power_limits::Bool = false
 end