Pull Jd/psy4 update

Project.toml

@@ -1,7 +1,7 @@
 name = "PowerSimulationsDynamics"
 uuid = "398b2ede-47ed-4edc-b52e-69e4a48b4336"
 authors = ["Jose Daniel Lara, Rodrigo Henriquez"]
-version = "0.14.0"
+version = "0.14.1"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

docs/make.jl

@@ -22,6 +22,7 @@ makedocs(;
         "Small Signal" => "small.md",
         "Reference Frames" => "reference_frames.md",
         "Perturbations" => "perturbations.md",
+        "Time Delays" => "time_delays.md",
         "Industrial Renewable Models" => "generic.md",
         "Generator Component Library" => Any[
             "Machine" => "component_models/machines.md",

docs/src/component_models/loads.md

@@ -160,4 +160,25 @@ Finally, the withdrawed current from the bus is:
 I_r = \left(\frac{S_\text{motor}}{S_\text{base}}\right) (i_{ds} - v_{qs} B_{sh}) \\
 I_i = \left(\frac{S_\text{motor}}{S_\text{base}}\right) (i_{qs} + v_{ds} B_{sh}) 
 \end{align*}
-```
+```
+
+### Active Constant Power Load Model
+
+The following 12-state model Active Load model  that measures the AC side using a Phase-Lock-Loop (PLL) and regulates a DC voltage to supply a resistor $r_L$. This model induces a constant power load-like behavior as it tries to maintain a fixed DC voltage to supply ``P = v_\text{DC}^2 / r_L``. The model is based on [the following reference](https://www.sciencedirect.com/science/article/pii/S0142061516000740). 
+
+The complete model is given by:
+```math
+\begin{align}
+    \dot{\theta} &= \Omega_b (\omega_\text{pll} - \omega_s) \tag{4a} \\
+    \dot{\epsilon} &= v_\text{o}^q \tag{4b}\\
+    \omega_\text{pll} &= \omega^\star + k^p_\text{pll} v_\text{o}^q + k_\text{pll}^i \epsilon \tag{4c}\\
+    \dot{\zeta} &= v_\text{DC}^\star - v_\text{DC} \tag{4d} \\
+    i_\text{cv}^{d,\star} &= k_\text{DC}^p ( v_\text{DC}^\star - v_\text{DC}) + k_\text{DC}^i \zeta \tag{4e}  \\
+    \frac{c_\text{DC}}{\Omega_b} \dot{v}_\text{DC} &= \frac{p_\text{cv}}{v_\text{DC}} - \frac{v_\text{DC}}{r_L} \tag{4f} \\
+    \dot{\gamma}_d &= i_\text{cv}^d - i_\text{cv}^{d,\star} \tag{4g}\\
+    \dot{\gamma}_q &= i_\text{cv}^q - i_\text{cv}^{q,\star} \tag{4h} \\
+    v_\text{cv}^{d,\star} &= k_\text{pc}( i_\text{cv}^d - i_\text{cv}^{d,\star}) + k_\text{ic} \gamma_d + \omega_\text{pll} l_f i_\text{cv}^q \tag{4i}\\
+    v_\text{cv}^{q,\star} &= k_\text{pc}( i_\text{cv}^q - i_\text{cv}^{q,\star}) + k_\text{ic} \gamma_q - \omega_\text{pll} l_f i_\text{cv}^d \tag{4j}
+\end{align}
+```
+Equations (4a)--(4c)  describes the PLL dynamics to lock the active load to the grid. Equations (4d)-(4e)  describes the DC Voltage Controller to steer the DC voltage to ``v_\text{DC}^\star``, while equation (4f) describes the DC voltage dynamics at the capacitor assuming an ideal converter. Finally, equations (4g)--(4j) describes the dynamics of the AC Current Controller. Additionally six states are defined for the LCL filter in a similar fashion of GFM inverters.

docs/src/component_models/machines.md

@@ -219,12 +219,12 @@ The Sauer Pai model defines 6 differential equations as follows:
 \begin{align}
 \dot{\psi}_d &= \Omega_b(r_ai_d + \omega \psi_q + v_d) \tag{9a} \\
 \dot{\psi}_q &= \Omega_b(r_ai_q - \omega \psi_d + v_q) \tag{9b} \\
-\dot{e}_q' &= \frac{1}{T_{d0}'} \left[(-e_q' - (x_d - x_d')(i_d + \gamma_d2 * \dot{\psi}_d'') + v_f)] \tag{9c}\\
-\dot{e}_q' &= \frac{1}{T_{q0}'} \left[(-e_d' + (x_q - x_q')(i_q + \gamma_q2 * \dot{\psi}_q''))] \tag{9d}\\
-\dot{\psi}_d'' &= \frac{1}{T_{d0}''} \left[(-\psi_d'' + e_q' - (x_d' - x_l)*i_d)] \tag{9e} \\
-\dot{\psi}_q'' &= \frac{1}{T_{q0}''} \left[(-\psi_q'' - e_d' - (x_q' - x_l)*i_q)]  \tag{9f} \\
-i_d &= \frac{1}{x_d''} (\gamma_d1 * e_q' - \psi_d + (1 - \gamma_d1) * \psi_d'') \tag{9g} \\
-i_q &= \frac{1}{x_q''} ((-\gamma_q1 * e_d' - \psi_q + (1 - \gamma_q1) * \psi_q'') \tag{9h} \\
+\dot{e}_q' &= \frac{1}{T_{d0}'} \left[(-e_q' - (x_d - x_d')(i_d + \gamma_{d2} \cdot  \dot{\psi}_d'') + v_f)\right] \tag{9c}\\
+\dot{e}_q' &= \frac{1}{T_{q0}'} \left[(-e_d' + (x_q - x_q')(i_q + \gamma_{q2} \cdot \dot{\psi}_q''))\right] \tag{9d}\\
+\dot{\psi}_d'' &= \frac{1}{T_{d0}''} \left[(-\psi_d'' + e_q' - (x_d' - x_l)\cdot i_d)\right] \tag{9e} \\
+\dot{\psi}_q'' &= \frac{1}{T_{q0}''} \left[(-\psi_q'' - e_d' - (x_q' - x_l)\cdot i_q)\right]  \tag{9f} \\
+i_d &= \frac{1}{x_d''} (\gamma_{d1} \cdot  e_q' - \psi_d + (1 - \gamma_{d1}) * \psi_d'') \tag{9g} \\
+i_q &= \frac{1}{x_q''} ((-\gamma_{q1} \cdot  e_d' - \psi_q + (1 - \gamma_{q1}) \cdot  \psi_q'') \tag{9h} \\
 \tau_e &= \psi_d i_q - \psi_q i_d \tag{9i}
 \end{align}
 ```
@@ -233,9 +233,9 @@ with
 
 ```math
 \begin{align*}
-  \gamma_d1 &= \frac{x_d'' - x_l}{x_d' - x_l} \\
-  \gamma_q1 &= \frac{x_q'' - x_l}{x_q' - x_l} \\
-  \gamma_d2 &= \frac{1 - \gamma_d1}{x_d' - x_l} \\
-  \gamma_q2 &= \frac{1 - \gamma_q1}{x_q' - x_l}
+  \gamma_{d1} &= \frac{x_d'' - x_l}{x_d' - x_l} \\
+  \gamma_{q1} &= \frac{x_q'' - x_l}{x_q' - x_l} \\
+  \gamma_{d2} &= \frac{1 - \gamma_d1}{x_d' - x_l} \\
+  \gamma_{q2} &= \frac{1 - \gamma_q1}{x_q' - x_l}
 \end{align*}
 ```

docs/src/component_models/outer_control.md

@@ -37,7 +37,7 @@ is the point of common coupling.
 ## Active Power Droop (P-droop) and Q-droop ```[OuterControl]```
 
 The following model represent a ``P\text{-}f`` droop model to represent how active
-power is going to be deployed. The constructor is ```OuterControl{ActivePowerControl, ReactivePowerDroop}```.
+power is going to be deployed. The constructor is ```OuterControl{ActivePowerDroop, ReactivePowerDroop}```.
 It defines a new SRF denoted as ``\theta_{\text{olc}}`` for the active power controller and
 uses a simple voltage droop for dispatching reactive power. Both active and reactive power are measured via a low-pass filter:
 

docs/src/models.md

@@ -65,7 +65,7 @@ The implementation of Synchronous generators as components uses the following st
 share values across components.
 
 ```@raw html
-<img src="https://github.com/nrel-sienna/PowerSystems.jl/blob/master/docs/src/assets/gen_metamodel.png?raw=true" width="75%">
+<img src="https://raw.githubusercontent.com/nrel-sienna/PowerSystems.jl/main/docs/src/assets/gen_metamodel.png" width="75%">
 ```
 
 ## Inverter Models
@@ -84,7 +84,7 @@ components in a way that resembles current practices for synchronoues machine mo
 The following figure summarizes the components of a inverter and which variables they share:
 
 ```@raw html
-<img src="https://github.com/nrel-sienna/PowerSystems.jl/blob/master/docs/src/assets/inv_metamodel.png?raw=true" width="75%">
+<img src="https://raw.githubusercontent.com/nrel-sienna/PowerSystems.jl/main/docs/src/assets/inv_metamodel.png" width="75%">
 ```
 
 Contrary to the generator, there are many control structures that can be used to model

docs/src/time_delays.md

@@ -0,0 +1,22 @@
+# Delays
+
+PowerSimulationsDynamics supports models with constant delays in a mass matrix formulation:
+
+```math
+\begin{align}
+M\frac{dx(t)}{dt} = f(x(t), x(t-\tau_1), ... , x(t-\tau_N))   
+\end{align}
+```
+
+For more information on solving such models, refer to the documentation for [DelayDiffEq.jl](https://github.com/SciML/DelayDiffEq.jl) package.
+
+The following models include time delays:
+
+* `DEGOV`
+
+There is currently limited support for including models with time delays. The following limitations apply:
+
+* Only constant delays are supported (state dependent delays are not).
+* System models with delays must use `MassMatrixModel`  formulation (`ResidualModel` is not currently compatible).
+* System models with delays are not compatible with small signal analysis tools.
+* The system formulation with delays is not compatible with automatic differentiation for calculating the gradient with respect to time. The setting `autodiff=false` should be set when passing the solver (e.g. `MethodofSteps(Rodas5(autodiff=false))`).

src/base/device_wrapper.jl

@@ -20,6 +20,12 @@ index(::Type{<:PSY.InnerControl}) = 4
 index(::Type{<:PSY.Converter}) = 5
 index(::Type{<:PSY.Filter}) = 6
 
+get_delays(::PSY.DynamicInjection) = nothing
+get_delays(
+    dynamic_injector::PSY.DynamicGenerator{M, S, A, PSY.DEGOV, P},
+) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS} =
+    Float64[PSY.get_Td(PSY.get_prime_mover(dynamic_injector))]
+
 """
 Wraps DynamicInjection devices from PowerSystems to handle changes in controls and connection
 status, and allocate the required indexes of the state space.

src/base/jacobian.jl

@@ -3,6 +3,7 @@
 # in SparseDiffTools
 
 struct JacobianFunctionWrapper{
+    D <: DelayModel,
     F,
     T <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}},
 } <: Function
@@ -13,12 +14,12 @@ struct JacobianFunctionWrapper{
 end
 
 # This version of the function is type unstable should only be used for non-critial ops
-function (J::JacobianFunctionWrapper)(x::AbstractVector{Float64})
+function (J::JacobianFunctionWrapper{NoDelays})(x::AbstractVector{Float64})
     J.x .= x
     return J.Jf(J.Jv, x)
 end
 
-function (J::JacobianFunctionWrapper)(
+function (J::JacobianFunctionWrapper{NoDelays})(
     JM::U,
     x::AbstractVector{Float64},
 ) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
@@ -27,7 +28,17 @@ function (J::JacobianFunctionWrapper)(
     return
 end
 
-function (J::JacobianFunctionWrapper)(
+function (J::JacobianFunctionWrapper{HasDelays})(
+    JM::U,
+    x::AbstractVector{Float64},
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
+    h(p, t; idxs = nothing) = typeof(idxs) <: Number ? x[idxs] : x
+    J.x .= x
+    J.Jf(JM, x, h, 0.0)
+    return
+end
+
+function (J::JacobianFunctionWrapper{NoDelays})(
     JM::U,
     x::AbstractVector{Float64},
     p,
@@ -37,7 +48,29 @@ function (J::JacobianFunctionWrapper)(
     return
 end
 
-function (J::JacobianFunctionWrapper)(
+function (J::JacobianFunctionWrapper{HasDelays})(
+    JM::U,
+    x::AbstractVector{Float64},
+    h,
+    p,
+    t,
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
+    J.Jf(JM, x, h, t)
+    return
+end
+
+function (J::JacobianFunctionWrapper{NoDelays})(
+    JM::U,
+    x::AbstractVector{Float64},
+    h,
+    p,
+    t,
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
+    J.Jf(JM, x, h, t)
+    return
+end
+
+function (J::JacobianFunctionWrapper{NoDelays})(
     JM::U,
     dx::AbstractVector{Float64},
     x::AbstractVector{Float64},
@@ -53,7 +86,7 @@ function (J::JacobianFunctionWrapper)(
 end
 
 function JacobianFunctionWrapper(
-    m!::SystemModel{MassMatrixModel},
+    m!::SystemModel{MassMatrixModel, NoDelays},
     x0_guess::Vector{Float64};
     # Improve the heuristic to do sparsity detection
     sparse_retrieve_loop::Int = 0, #max(3, length(x0_guess) ÷ 100),
@@ -84,7 +117,12 @@ function JacobianFunctionWrapper(
     end
     Jf(Jv, x0)
     mass_matrix = get_mass_matrix(m!.inputs)
-    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0, mass_matrix)
+    return JacobianFunctionWrapper{NoDelays, typeof(Jf), typeof(Jv)}(
+        Jf,
+        Jv,
+        x0,
+        mass_matrix,
+    )
 end
 
 function JacobianFunctionWrapper(
@@ -119,7 +157,53 @@ function JacobianFunctionWrapper(
     end
     Jf(Jv, x0)
     mass_matrix = get_mass_matrix(m!.inputs)
-    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0, mass_matrix)
+    return JacobianFunctionWrapper{NoDelays, typeof(Jf), typeof(Jv)}(
+        Jf,
+        Jv,
+        x0,
+        mass_matrix,
+    )
+end
+
+function JacobianFunctionWrapper(
+    m!::SystemModel{MassMatrixModel, HasDelays},
+    x0_guess::Vector{Float64};
+    # Improve the heuristic to do sparsity detection
+    sparse_retrieve_loop::Int = 0, #max(3, length(x0_guess) ÷ 100),
+)
+    x0 = deepcopy(x0_guess)
+    n = length(x0)
+    Jf =
+        (Jv, x, h, t) -> begin
+            @debug "Evaluating Jacobian Function"
+            m_ = (residual, x) -> m!(residual, x, h, nothing, t)
+            jconfig =
+                ForwardDiff.JacobianConfig(m_, similar(x0), x0, ForwardDiff.Chunk(x0))
+            ForwardDiff.jacobian!(Jv, m_, zeros(n), x, jconfig)
+            return
+        end
+    jac = zeros(n, n)
+    if sparse_retrieve_loop > 0
+        for _ in 1:sparse_retrieve_loop
+            temp = zeros(n, n)
+            rng_state = copy(Random.default_rng())
+            Jf(temp, x0 + Random.randn(n))
+            copy!(Random.default_rng(), rng_state)
+            jac .+= abs.(temp)
+        end
+        Jv = SparseArrays.sparse(jac)
+    elseif sparse_retrieve_loop == 0
+        Jv = jac
+    else
+        throw(IS.ConflictingInputsError("negative sparse_retrieve_loop not valid"))
+    end
+    mass_matrix = get_mass_matrix(m!.inputs)
+    return JacobianFunctionWrapper{HasDelays, typeof(Jf), typeof(Jv)}(
+        Jf,
+        Jv,
+        x0,
+        mass_matrix,
+    )
 end
 
 function get_jacobian(

src/base/nlsolve_wrapper.jl

@@ -8,11 +8,25 @@ NLsolveWrapper() = NLsolveWrapper(Vector{Float64}(), false, true)
 converged(sol::NLsolveWrapper) = sol.converged
 failed(sol::NLsolveWrapper) = sol.failed
 
-function _get_model_closure(model::SystemModel{MassMatrixModel}, ::Vector{Float64})
+function _get_model_closure(
+    model::SystemModel{MassMatrixModel, NoDelays},
+    ::Vector{Float64},
+)
     return (residual, x) -> model(residual, x, nothing, 0.0)
 end
 
-function _get_model_closure(model::SystemModel{ResidualModel}, x0::Vector{Float64})
+function _get_model_closure(
+    model::SystemModel{MassMatrixModel, HasDelays},
+    x0::Vector{Float64},
+)
+    h(p, t; idxs = nothing) = typeof(idxs) <: Number ? x0[idxs] : x0
+    return (residual, x) -> model(residual, x, h, nothing, 0.0)
+end
+
+function _get_model_closure(
+    model::SystemModel{ResidualModel, NoDelays},
+    x0::Vector{Float64},
+)
     dx0 = zeros(length(x0))
     return (residual, x) -> model(residual, dx0, x, nothing, 0.0)
 end

src/base/simulation.jl

@@ -289,7 +289,10 @@ end
 function _get_jacobian(sim::Simulation{MassMatrixModel})
     inputs = get_simulation_inputs(sim)
     x0_init = get_initial_conditions(sim)
-    return JacobianFunctionWrapper(MassMatrixModel(inputs, x0_init, JacobianCache), x0_init)
+    return JacobianFunctionWrapper(
+        MassMatrixModel(inputs, x0_init, JacobianCache),
+        x0_init,
+    )
 end
 
 function _build_perturbations!(sim::Simulation)
@@ -329,7 +332,7 @@ end
 
 function _get_diffeq_problem(
     sim::Simulation,
-    model::SystemModel{ResidualModel},
+    model::SystemModel{ResidualModel, NoDelays},
     jacobian::JacobianFunctionWrapper,
 )
     x0 = get_initial_conditions(sim)
@@ -354,7 +357,7 @@ end
 
 function _get_diffeq_problem(
     sim::Simulation,
-    model::SystemModel{MassMatrixModel},
+    model::SystemModel{MassMatrixModel, NoDelays},
     jacobian::JacobianFunctionWrapper,
 )
     simulation_inputs = get_simulation_inputs(sim)
@@ -375,6 +378,37 @@ function _get_diffeq_problem(
     return
 end
 
+function get_history_function(simulation_inputs::Simulation{MassMatrixModel})
+    x0 = get_initial_conditions(simulation_inputs)
+    h(p, t; idxs = nothing) = typeof(idxs) <: Number ? x0[idxs] : x0
+    return h
+end
+
+function _get_diffeq_problem(
+    sim::Simulation,
+    model::SystemModel{MassMatrixModel, HasDelays},
+    jacobian::JacobianFunctionWrapper,
+)
+    simulation_inputs = get_simulation_inputs(sim)
+    h = get_history_function(sim)
+    sim.problem = SciMLBase.DDEProblem(
+        SciMLBase.DDEFunction{true}(
+            model;
+            mass_matrix = get_mass_matrix(simulation_inputs),
+            jac = jacobian,
+            jac_prototype = jacobian.Jv,
+        ),
+        sim.x0_init,
+        h,
+        get_tspan(sim),
+        simulation_inputs;
+        constant_lags = filter!(x -> x != 0, unique(simulation_inputs.delays)),
+    )
+    sim.status = BUILT
+
+    return
+end
+
 function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
     check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
     # Branches are a super set of Lines. Passing both kwargs will