work in progress (add delays)

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} =
+    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

@@ -37,6 +37,17 @@ function (J::JacobianFunctionWrapper)(
     return
 end
 
+function (J::JacobianFunctionWrapper)(
+    JM::U,
+    x::AbstractVector{Float64},
+    h,
+    p,
+    t,
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
+    J(JM, x)
+    return
+end
+
 function (J::JacobianFunctionWrapper)(
     JM::U,
     dx::AbstractVector{Float64},
@@ -53,7 +64,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),
@@ -122,6 +133,43 @@ function JacobianFunctionWrapper(
     return JacobianFunctionWrapper{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)
+
+    h(p, t; idxs = nothing) = typeof(idxs) <: Number ? x0[idxs] : x0  #Possilby the problem? Should pass the h that is stored in the problem instead of redefining?
+    m_ = (residual, x) -> m!(residual, x, h, nothing, 0.0)
+    jconfig = ForwardDiff.JacobianConfig(m_, similar(x0), x0, ForwardDiff.Chunk(x0))
+    Jf = (Jv, x) -> begin
+        @debug "Evaluating Jacobian Function"
+        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
+    Jf(Jv, x0)
+    mass_matrix = get_mass_matrix(m!.inputs)
+    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0, mass_matrix)
+end
+
 function get_jacobian(
     ::Type{T},
     inputs::SimulationInputs,

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

@@ -329,7 +329,7 @@ end
 
 function _get_diffeq_problem(
     sim::Simulation,
-    model::SystemModel{ResidualModel},
+    model::SystemModel{ResidualModel, NoDelays},
     jacobian::JacobianFunctionWrapper,
 )
     x0 = get_initial_conditions(sim)
@@ -354,7 +354,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 +375,39 @@ function _get_diffeq_problem(
     return
 end
 
+function get_history_function(simulation_inputs)
+    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,              #Fails after time of delay if autodiff=true with non-zero delays
+            jac_prototype = jacobian.Jv,
+            # Necessary to avoid unnecessary calculations in Rosenbrock methods
+            #tgrad = (dT, u, h, p, t) -> dT .= false,      #Doesn't work with passing tgrad yet
+        ),
+        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

src/base/simulation_inputs.jl

@@ -18,6 +18,7 @@ struct SimulationInputs
     global_vars_update_pointers::Dict{Int, Int}
     global_state_map::MAPPING_DICT
     global_inner_var_map::Dict{String, Dict}
+    delays::Vector
 
     function SimulationInputs(
         sys::PSY.System,
@@ -39,6 +40,7 @@ struct SimulationInputs
 
         TimerOutputs.@timeit BUILD_TIMER "Wrap Dynamic Injectors" begin
             wrapped_injectors = _wrap_dynamic_injector_data(sys, lookup, injection_start)
+            delays = get_system_delays(sys)
             var_count = wrapped_injectors[end].ix_range[end]
         end
 
@@ -71,6 +73,11 @@ struct SimulationInputs
                 break
             end
         end
+
+        if !isempty(delays)
+            @info "System has delays. Use the correct solver for delay differential equations."
+        end
+
         new(
             wrapped_injectors,
             wrapped_static_injectors,
@@ -91,6 +98,7 @@ struct SimulationInputs
             global_vars,
             MAPPING_DICT(),
             Dict{String, Dict}(),
+            delays,
         )
     end
 end
@@ -173,6 +181,7 @@ function _wrap_dynamic_injector_data(sys::PSY.System, lookup, injection_start::I
         @debug "ix_range=$ix_range ode_range=$ode_range inner_vars_range= $inner_vars_range"
         dynamic_device = PSY.get_dynamic_injector(device)
         @assert dynamic_device !== nothing
+        #TODO - add bus_size to the dynamic injector -- allows you to index into the voltages.
         wrapped_injector[ix] = DynamicWrapper(
             device,
             dynamic_device,
@@ -191,6 +200,15 @@ function _wrap_dynamic_injector_data(sys::PSY.System, lookup, injection_start::I
     return wrapped_injector
 end
 
+function get_system_delays(sys::PSY.System)
+    delays = []
+    for injector in get_injectors_with_dynamics(sys)
+        device_delays = get_delays(PSY.get_dynamic_injector(injector))
+        !isnothing(device_delays) && push!(delays, device_delays)
+    end
+    return delays
+end
+
 function _wrap_dynamic_branches(sys::PSY.System, lookup::Dict{Int, Int})
     branches_start = 2 * get_n_buses(sys) + 1
     sys_base_power = PSY.get_base_power(sys)

src/base/simulation_model.jl

@@ -1,4 +1,9 @@
-abstract type SimulationModel end
+abstract type ModelVariations end
 
+abstract type SimulationModel <: ModelVariations end
 struct MassMatrixModel <: SimulationModel end
 struct ResidualModel <: SimulationModel end
+
+abstract type DelayModel <: ModelVariations end
+struct HasDelays <: DelayModel end
+struct NoDelays <: DelayModel end

src/base/system_model.jl

@@ -1,8 +1,11 @@
-struct SystemModel{T <: SimulationModel, C <: Cache}
+struct SystemModel{T <: SimulationModel, D <: DelayModel, C <: Cache}
     inputs::SimulationInputs
     cache::C
 end
 
-function SystemModel{T}(inputs, cache::U) where {T <: SimulationModel, U <: Cache}
-    return SystemModel{T, U}(inputs, cache)
+function SystemModel{T, D}(
+    inputs,
+    cache::U,
+) where {T <: SimulationModel, D <: DelayModel, U <: Cache}
+    return SystemModel{T, D, U}(inputs, cache)
 end

src/models/device.jl

@@ -25,6 +25,7 @@ function device!(
     global_vars::AbstractArray{T},
     inner_vars::AbstractArray{T},
     dynamic_device::DynamicWrapper{DynG},
+    h,
     t,
 ) where {DynG <: PSY.DynamicGenerator, T <: ACCEPTED_REAL_TYPES}
     if get_connection_status(dynamic_device) < 1.0
@@ -41,13 +42,13 @@ function device!(
     _update_inner_vars!(device_states, output_ode, sys_ω, inner_vars, dynamic_device)
 
     #Obtain ODEs and Mechanical Power for Turbine Governor
-    mdl_tg_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)
+    mdl_tg_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device, h, t)
 
     #Obtain ODEs for PSS
-    mdl_pss_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)
+    mdl_pss_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device, h, t)
 
     #Obtain ODEs for AVR
-    mdl_avr_ode!(device_states, output_ode, inner_vars, dynamic_device)
+    mdl_avr_ode!(device_states, output_ode, inner_vars, dynamic_device, h, t)
 
     #Obtain ODEs for Machine
     mdl_machine_ode!(
@@ -57,10 +58,12 @@ function device!(
         current_r,
         current_i,
         dynamic_device,
+        h,
+        t,
     )
 
     #Obtain ODEs for PSY.Shaft
-    mdl_shaft_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)
+    mdl_shaft_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device, h, t)
     return
 end
 
@@ -142,6 +145,7 @@ function device!(
     global_vars::AbstractArray{T},
     inner_vars::AbstractArray{T},
     dynamic_device::DynamicWrapper{DynI},
+    h,
     t,
 ) where {DynI <: PSY.DynamicInverter, T <: ACCEPTED_REAL_TYPES}
     if get_connection_status(dynamic_device) < 1.0
@@ -164,19 +168,27 @@ function device!(
     _update_inner_vars!(device_states, output_ode, sys_ω, inner_vars, dynamic_device)
 
     #Obtain ODES for DC side
-    mdl_DCside_ode!(device_states, output_ode, sys_ω, inner_vars, dynamic_device)
+    mdl_DCside_ode!(device_states, output_ode, sys_ω, inner_vars, dynamic_device, h, t)
 
     #Obtain ODEs for PLL
-    mdl_freq_estimator_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)
+    mdl_freq_estimator_ode!(
+        device_states,
+        output_ode,
+        inner_vars,
+        sys_ω,
+        dynamic_device,
+        h,
+        t,
+    )
 
     #Obtain ODEs for OuterLoop
-    mdl_outer_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)
+    mdl_outer_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device, h, t)
 
     #Obtain inner controller ODEs and modulation commands
-    mdl_inner_ode!(device_states, output_ode, inner_vars, dynamic_device)
+    mdl_inner_ode!(device_states, output_ode, inner_vars, dynamic_device, h, t)
 
     #Obtain converter relations
-    mdl_converter_ode!(device_states, output_ode, inner_vars, dynamic_device)
+    mdl_converter_ode!(device_states, output_ode, inner_vars, dynamic_device, h, t)
 
     #Obtain ODEs for output filter
     mdl_filter_ode!(
@@ -187,6 +199,8 @@ function device!(
         inner_vars,
         sys_ω,
         dynamic_device,
+        h,
+        t,
     )
 
     return
@@ -219,6 +233,7 @@ function device!(
     ::AbstractArray{T},
     ::AbstractArray{T},
     dynamic_device::DynamicWrapper{PSY.PeriodicVariableSource},
+    h,
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     ω_θ = PSY.get_internal_angle_frequencies(get_device(dynamic_device))
@@ -643,6 +658,7 @@ function device!(
     global_vars::AbstractArray{T},
     ::AbstractArray{T},
     dynamic_wrapper::DynamicWrapper{PSY.SingleCageInductionMachine},
+    h,
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     Sbase = get_system_base_power(dynamic_wrapper)
@@ -726,6 +742,7 @@ function device!(
     global_vars::AbstractArray{T},
     ::AbstractArray{T},
     dynamic_wrapper::DynamicWrapper{PSY.SimplifiedSingleCageInductionMachine},
+    h,
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     Sbase = get_system_base_power(dynamic_wrapper)
@@ -822,6 +839,7 @@ function device!(
     global_vars::AbstractArray{T},
     ::AbstractArray{T},
     dynamic_wrapper::DynamicWrapper{PSY.CSVGN1},
+    h,
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     Sbase = get_system_base_power(dynamic_wrapper)
@@ -930,6 +948,7 @@ function device!(
     global_vars::AbstractArray{T},
     ::AbstractArray{T},
     dynamic_wrapper::DynamicWrapper{PSY.ActiveConstantPowerLoad},
+    h,
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     Sbase = get_system_base_power(dynamic_wrapper)
@@ -1076,6 +1095,7 @@ function device!(
     global_vars::AbstractArray{T},
     inner_vars::AbstractArray{T},
     dynamic_device::DynamicWrapper{PSY.AggregateDistributedGenerationA},
+    h,
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     Freq_Flag = PSY.get_Freq_Flag(get_device(dynamic_device))