simulation flow upgrades (re-use inputs and initial conditions)

src/base/definitions.jl

@@ -220,7 +220,7 @@ const DIFFEQ_SOLVE_KWARGS = [
 """
 Defines the status of the simulation object
 """
-@enum BUILD_STATUS begin
+@enum STATUS begin
     BUILT = 0
     BUILD_IN_PROGRESS = 1
     BUILD_INCOMPLETE = 2
@@ -231,6 +231,26 @@ Defines the status of the simulation object
     SIMULATION_FAILED = 7
 end
 
+"""
+Defines the level of building simulation inputs
+"""
+@enum BUILD_INPUTS_LEVEL begin
+    BUILD_ONE = 1
+    BUILD_TWO = 2
+    BUILD_THREE = 3
+    BUILD_NONE = 4
+end
+
+"""
+Defines the level of initializing simulation
+"""
+@enum INITIALIZE_LEVEL begin
+    POWERFLOW_AND_DEVICES = 0
+    DEVICES_ONLY = 1
+    FLAT_START = 2
+    INITIALIZED = 3
+end
+
 const BUILD_TIMER = TimerOutputs.TimerOutput()
 
 const ACCEPTED_REAL_TYPES = Union{Float64, ForwardDiff.Dual}

src/base/jacobian.jl

@@ -245,9 +245,9 @@ function get_jacobian(
 ) where {T <: SimulationModel}
     # Deepcopy avoid system modifications
     simulation_system = deepcopy(system)
-    inputs = SimulationInputs(T, simulation_system, ReferenceBus())
+    inputs = SimulationInputs(T, simulation_system, ReferenceBus(), nothing, Val(BUILD_ONE))
     p = get_parameters(inputs)
-    x0_init = get_flat_start(inputs)
+    x0_init = _get_flat_start(inputs)
     set_operating_point!(x0_init, inputs, system)
     return get_jacobian(T, inputs, x0_init, p, sparse_retrieve_loop)
 end

src/base/nlsolve_wrapper.jl

@@ -119,7 +119,7 @@ function _check_residual(
     if sum_residual > tolerance
         state_map = make_global_state_map(inputs)
         for (k, val) in state_map
-            inputs.global_state_map[k] = val
+            get_global_state_map(inputs)[k] = val
         end
         gen_name = ""
         state = ""
@@ -136,6 +136,7 @@ function _check_residual(
                 Generator = $gen_name, state = $state.
                Residual error is too large to continue")
         else
+            bus_count = get_bus_count(inputs)
             bus_no = ix > bus_count ? ix - bus_count : ix
             component = ix > bus_count ? "imag" : "real"
             error("The initial residual in the state located at $ix has a value of $val.
@@ -150,15 +151,11 @@ function refine_initial_condition!(
     sim::Simulation,
     model::SystemModel,
     jacobian::JacobianFunctionWrapper,
+    ::Val{POWERFLOW_AND_DEVICES},
 )
     @assert sim.status != BUILD_INCOMPLETE
-
-    if sim.status == SIMULATION_INITIALIZED
-        @info "Simulation already initialized. Refinement not executed"
-        return
-    end
     converged = false
-    initial_guess = get_initial_conditions(sim)
+    initial_guess = get_x0(sim)
     inputs = get_simulation_inputs(sim)
     parameters = get_parameters(inputs)
     bus_range = get_bus_range(inputs)
@@ -205,6 +202,29 @@ function refine_initial_condition!(
     if maximum(pf_diff) > MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE
         @warn "The resulting voltages in the initial conditions differ from the power flow results"
     end
-
     return
 end
+
+function refine_initial_condition!(
+    sim::Simulation,
+    model::SystemModel,
+    jacobian::JacobianFunctionWrapper,
+    ::Val{DEVICES_ONLY},
+)
+    refine_initial_condition!(sim, model, jacobian, Val(POWERFLOW_AND_DEVICES))
+end
+
+function refine_initial_condition!(
+    sim::Simulation,
+    model::SystemModel,
+    jacobian::JacobianFunctionWrapper,
+    ::Val{FLAT_START},
+)
+end
+function refine_initial_condition!(
+    sim::Simulation,
+    model::SystemModel,
+    jacobian::JacobianFunctionWrapper,
+    ::Val{INITIALIZED},
+)
+end

src/base/perturbations.jl

@@ -328,8 +328,8 @@ function ybus_update!(
     return
 end
 
-function ybus_update!(integrator_params, branch::PSY.ACBranch, mult::Float64)
-    ybus_update!(integrator_params.ybus_rectangular, branch, integrator_params.lookup, mult)
+function ybus_update!(inputs, branch::PSY.ACBranch, mult::Float64)
+    ybus_update!(get_ybus(inputs), branch, get_lookup(inputs), mult)
     return
 end
 
@@ -372,7 +372,7 @@ function get_affect(inputs::SimulationInputs, ::PSY.System, pert::NetworkSwitch)
         # TODO: This code can be more performant using SparseMatrix methods
         for (i, v) in enumerate(pert.ybus_rectangular)
             @debug "Changing Ybus network"
-            inputs.ybus_rectangular[i] = v
+            get_ybus(inputs)[i] = v
         end
         return
     end

src/base/simulation.jl

@@ -1,15 +1,18 @@
 mutable struct Simulation{T <: SimulationModel}
-    status::BUILD_STATUS
+    status::STATUS
     problem::Union{Nothing, SciMLBase.DEProblem}
     tspan::NTuple{2, Float64}
     sys::PSY.System
     perturbations::Vector{<:Perturbation}
+    initialize_level::INITIALIZE_LEVEL
+    x0::Vector{Float64}
     x0_init::Vector{Float64}
-    initialized::Bool
     tstops::Vector{Float64}
     callbacks::Vector
     simulation_folder::String
+    build_inputs_level::BUILD_INPUTS_LEVEL
     inputs::Union{Nothing, SimulationInputs}
+    inputs_init::Union{Nothing, SimulationInputs}
     results::Union{Nothing, SimulationResults}
     console_level::Base.CoreLogging.LogLevel
     file_level::Base.CoreLogging.LogLevel
@@ -19,7 +22,8 @@ end
 
 get_system(sim::Simulation) = sim.sys
 get_simulation_inputs(sim::Simulation) = sim.inputs
-get_initial_conditions(sim::Simulation) = sim.x0_init
+get_x0(sim::Simulation) = sim.x0
+get_x0_init(sim::Simulation) = sim.x0_init
 get_tspan(sim::Simulation) = sim.tspan
 
 function Simulation(
@@ -35,18 +39,30 @@ function Simulation(
     frequency_reference,
 ) where {T <: SimulationModel}
     PSY.set_units_base_system!(sys, "DEVICE_BASE")
-
+    if initialize_simulation && !isempty(initial_conditions)
+        @warn "initial_conditions were provided with initialize_simulation. User's initial_conditions will be overwritten."
+        init_level = POWERFLOW_AND_DEVICES
+    elseif initialize_simulation
+        init_level = POWERFLOW_AND_DEVICES
+    elseif !initialize_simulation && isempty(initial_conditions)
+        init_level = FLAT_START
+    else
+        init_level = INITIALIZED
+    end
     return Simulation{T}(
         BUILD_INCOMPLETE,
         nothing,
         tspan,
         sys,
         perturbations,
+        init_level,
+        initial_conditions,
         initial_conditions,
-        !initialize_simulation,
         Vector{Float64}(),
         Vector{SciMLBase.AbstractDiscreteCallback}(),
         simulation_folder,
+        BUILD_ONE,
+        nothing,
         nothing,
         nothing,
         console_level,
@@ -194,9 +210,7 @@ end
 
 function reset!(sim::Simulation{T}) where {T <: SimulationModel}
     @info "Rebuilding the simulation after reset"
-    sim.inputs = SimulationInputs(T, get_system(sim), sim.frequency_reference)
     sim.status = BUILD_INCOMPLETE
-    sim.initialized = false
     build!(sim)
     @info "Simulation reset to status $(sim.status)"
     return
@@ -219,7 +233,14 @@ end
 
 function _build_inputs!(sim::Simulation{T}) where {T <: SimulationModel}
     simulation_system = get_system(sim)
-    sim.inputs = SimulationInputs(T, simulation_system, sim.frequency_reference)
+    sim.inputs = SimulationInputs(
+        T,
+        simulation_system,
+        sim.frequency_reference,
+        sim.inputs_init,
+        Val(sim.build_inputs_level),
+    )
+    sim.inputs_init = deepcopy(sim.inputs)
     @debug "Simulation Inputs Created"
     return
 end
@@ -232,53 +253,14 @@ function _get_flat_start(inputs::SimulationInputs)
     return initial_conditions
 end
 
-function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel}
-    simulation_inputs = get_simulation_inputs(sim)
-    x0_init = sim.x0_init
-    if isempty(x0_init) && sim.initialized
-        @warn "Initial Conditions set to flat start"
-        sim.x0_init = _get_flat_start(simulation_inputs)
-    elseif isempty(x0_init) && !sim.initialized
-        sim.x0_init = _get_flat_start(simulation_inputs)
-    elseif !isempty(x0_init) && sim.initialized
-        if length(sim.x0_init) != get_variable_count(simulation_inputs)
-            throw(
-                IS.ConflictingInputsError(
-                    "The size of the provided initial state space does not match the model's state space.",
-                ),
-            )
-        end
-    elseif !isempty(x0_init) && !sim.initialized
-        @warn "initial_conditions were provided with initialize_simulation. User's initial_conditions will be overwritten."
-        if length(sim.x0_init) != get_variable_count(simulation_inputs)
-            sim.x0_init = _get_flat_start(simulation_inputs)
-        end
-    else
-        @assert false
-    end
-    return
-end
-
 function _pre_initialize_simulation!(sim::Simulation)
-    _initialize_state_space(sim)
-    if sim.initialized != true
-        @info("Pre-Initializing Simulation States")
-        sim.initialized = precalculate_initial_conditions!(sim)
-        if !sim.initialized
-            error(
-                "The simulation failed to find an adequate initial guess for the initialization. Check the intialization routine.",
-            )
-        end
-    else
-        @warn("Using existing initial conditions value for simulation initialization")
-        sim.status = SIMULATION_INITIALIZED
-    end
+    _initialize_state_space(sim, Val(sim.initialize_level))
     return
 end
 
 function _get_jacobian(sim::Simulation{ResidualModel})
     inputs = get_simulation_inputs(sim)
-    x0_init = get_initial_conditions(sim)
+    x0_init = get_x0(sim)
     parameters = get_parameters(inputs)
     return JacobianFunctionWrapper(
         ResidualModel(inputs, x0_init, JacobianCache),
@@ -290,7 +272,7 @@ end
 
 function _get_jacobian(sim::Simulation{MassMatrixModel})
     inputs = get_simulation_inputs(sim)
-    x0_init = get_initial_conditions(sim)
+    x0_init = get_x0(sim)
     parameters = get_parameters(inputs)
     return JacobianFunctionWrapper(
         MassMatrixModel(inputs, x0_init, JacobianCache),
@@ -339,7 +321,8 @@ function _get_diffeq_problem(
     model::SystemModel{ResidualModel, NoDelays},
     jacobian::JacobianFunctionWrapper,
 )
-    x0 = get_initial_conditions(sim)
+    x0 = get_x0(sim)
+    sim.x0_init = deepcopy(x0)
     dx0 = zeros(length(x0))
     simulation_inputs = get_simulation_inputs(sim)
     p = get_parameters(simulation_inputs)
@@ -365,6 +348,8 @@ function _get_diffeq_problem(
     model::SystemModel{MassMatrixModel, NoDelays},
     jacobian::JacobianFunctionWrapper,
 )
+    x0 = get_x0(sim)
+    sim.x0_init = deepcopy(x0)
     simulation_inputs = get_simulation_inputs(sim)
     p = get_parameters(simulation_inputs)
     sim.problem = SciMLBase.ODEProblem(
@@ -376,7 +361,7 @@ function _get_diffeq_problem(
             # Necessary to avoid unnecessary calculations in Rosenbrock methods
             tgrad = (dT, u, p, t) -> dT .= false,
         ),
-        sim.x0_init,
+        x0,
         get_tspan(sim),
         p;
     )
@@ -385,7 +370,7 @@ function _get_diffeq_problem(
 end
 
 function get_history_function(simulation_inputs::Simulation{MassMatrixModel})
-    x0 = get_initial_conditions(simulation_inputs)
+    x0 = get_x0(simulation_inputs)
     h(p, t; idxs = nothing) = typeof(idxs) <: Number ? x0[idxs] : x0
     return h
 end
@@ -395,6 +380,8 @@ function _get_diffeq_problem(
     model::SystemModel{MassMatrixModel, HasDelays},
     jacobian::JacobianFunctionWrapper,
 )
+    x0 = get_x0(sim)
+    sim.x0_init = deepcopy(x0)
     simulation_inputs = get_simulation_inputs(sim)
     h = get_history_function(sim)
     p = get_parameters(simulation_inputs)
@@ -405,11 +392,11 @@ function _get_diffeq_problem(
             jac = jacobian,
             jac_prototype = jacobian.Jv,
         ),
-        sim.x0_init,
+        x0,
         h,
         get_tspan(sim),
         p;
-        constant_lags = filter!(x -> x != 0, unique(simulation_inputs.delays)),
+        constant_lags = filter!(x -> x != 0, unique(get_delays(simulation_inputs))),
     )
     sim.status = BUILT
 
@@ -455,10 +442,15 @@ function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
                     jacobian = _get_jacobian(sim)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Make Model Function" begin
-                    model = T(simulation_inputs, get_initial_conditions(sim), SimCache)
+                    model = T(simulation_inputs, get_x0(sim), SimCache)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Initial Condition NLsolve refinement" begin
-                    refine_initial_condition!(sim, model, jacobian)
+                    refine_initial_condition!(
+                        sim,
+                        model,
+                        jacobian,
+                        Val(sim.initialize_level),
+                    )
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Build Perturbations" begin
                     _build_perturbations!(sim)