simulation flow upgrades

.github/workflows/formatter/Project.toml

@@ -6,5 +6,5 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [compat]
-JuliaFormatter = "v0.16.2"
-julia = "^1.2"
+JuliaFormatter = "v0.21"
+julia = "^1.7"

Project.toml

@@ -6,7 +6,6 @@ version = "0.8.5"
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
@@ -22,7 +21,6 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 DataStructures = "~0.18"
-DiffEqBase = "^6.8"
 DocStringExtensions = "~0.8.6"
 ForwardDiff = "~v0.10"
 InfrastructureSystems = "^1.14"

docs/src/execute.md

@@ -75,7 +75,7 @@ Any solver option available in `DifferentialEquations.jl` can be passed as keywo
 Most common solver options used are `dtmax` to control the maximum dt for adaptive timestepping. `abstol` and `reltol` are also commonly used to control the tolerance in the adaptive timestepping. `saveat` is also used to store the results at a specified time stamps. For example, the following code is valid to further specify your solver options:
 
 ```julia
-execute!(sim, IDA(), dtmax = 0.01, abstol = 1e-9, reltol = 1e-6, saveat = 0.01)
+execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.01, abstol = 1e-9, reltol = 1e-6, saveat = 0.01)
 ```
 
-In addition, the keyword argument `enable_progress_bar = false` can be used to disable the progress bar.
+In addition, the keyword argument `enable_progress_bar = false` can be used to disable the progress bar.

docs/src/quick_start_guide.md

@@ -64,7 +64,7 @@ This means that the state ω of the generator at bus 102, participates 50% in ei
 ## Execute the simulation
 
 ```@repl quick_start_guide
-    execute!(sim, IDA(), dtmax = 0.02, saveat = 0.02, enable_progress_bar = false)
+    execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02, saveat = 0.02, enable_progress_bar = false)
 ```
 
 ## Make a plot of the results

src/PowerSimulationsDynamics.jl

@@ -30,25 +30,25 @@ export LoadChange
 # export BusTrip
 
 # Export for routines
-export small_signal_analysis
-export get_state_series
-export get_voltage_magnitude_series
-export get_voltage_angle_series
-export show_states_initial_value
-export read_initial_conditions
-export get_real_current_series
-export get_imaginary_current_series
 export get_activepower_series
+export get_jacobian
+export get_imaginary_current_series
 export get_reactivepower_series
+export get_real_current_series
 export get_setpoints
+export get_state_series
+export get_voltage_angle_series
+export get_voltage_magnitude_series
+export read_initial_conditions
+export show_states_initial_value
+export small_signal_analysis
 
 ####################################### Package Imports ####################################
 import Logging
 import InfrastructureSystems
 import SciMLBase
 import DataStructures: OrderedDict
 import Random
-import DiffEqBase
 import ForwardDiff
 import SparseArrays
 import LinearAlgebra
@@ -79,6 +79,7 @@ include("base/branch_wrapper.jl")
 include("base/frequency_reference.jl")
 include("base/file_system.jl")
 include("base/simulation_model.jl")
+include("base/simulation_internal.jl")
 include("base/simulation_inputs.jl")
 include("base/perturbations.jl")
 include("base/caches.jl")

src/base/jacobian.jl

@@ -2,9 +2,12 @@
 # custom code since the structure of the system models is not compatible with the functionalities
 # in SparseDiffTools
 
-struct JacobianFunctionWrapper{F} <: Function
+struct JacobianFunctionWrapper{
+    F,
+    T <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}},
+} <: Function
     Jf::F
-    Jv::SparseArrays.SparseMatrixCSC{Float64, Int64}
+    Jv::T
     x::Vector{Float64}
 end
 
@@ -31,7 +34,7 @@ function JacobianFunctionWrapper(
     end
     Jv = SparseArrays.sparse(jac)
     Jf(Jv, x0)
-    return JacobianFunctionWrapper{typeof(Jf)}(Jf, Jv, x0)
+    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0)
 end
 
 function JacobianFunctionWrapper(
@@ -56,7 +59,7 @@ function JacobianFunctionWrapper(
     end
     Jv = SparseArrays.sparse(jac)
     Jf(Jv, x0)
-    return JacobianFunctionWrapper{typeof(Jf)}(Jf, Jv, x0)
+    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0)
 end
 
 function (J::JacobianFunctionWrapper)(x::AbstractVector{Float64})
@@ -65,34 +68,62 @@ function (J::JacobianFunctionWrapper)(x::AbstractVector{Float64})
 end
 
 function (J::JacobianFunctionWrapper)(
-    JM::SparseArrays.SparseMatrixCSC{Float64, Int64},
+    JM::U,
     x::AbstractVector{Float64},
-)
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
     J.x .= x
     J.Jf(JM, x)
     return
 end
 
 function (J::JacobianFunctionWrapper)(
-    JM::SparseArrays.SparseMatrixCSC{Float64, Int64},
+    JM::U,
     x::AbstractVector{Float64},
     p,
     t,
-)
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
     J.x .= x
     J.Jf(JM, x)
     return
 end
 
 function (J::JacobianFunctionWrapper)(
-    JM::SparseArrays.SparseMatrixCSC{Float64, Int64},
+    JM::U,
     dx::AbstractVector{Float64},
     x::AbstractVector{Float64},
     p,
     gamma,
     t,
-)
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
     J.x .= x
     JM .= gamma * LinearAlgebra.Diagonal(ones(length(x))) .- J.Jf(JM, x)
     return JM
 end
+
+function get_jacobian(
+    ::Type{T},
+    inputs::SimulationInputs,
+    x0_init::Vector{Float64},
+) where {T <: SimulationModel}
+    return JacobianFunctionWrapper(T(inputs, x0_init, JacobianCache), x0_init)
+end
+
+function _set_operating_point!(
+    x0_init::Vector{Float64},
+    inputs::SimulationInputs,
+    system::PSY.System,
+)
+    status = power_flow_solution!(x0_init, system, inputs)
+    status = initialize_static_injection!(inputs)
+    status = initialize_dynamic_injection!(x0_init, inputs, system)
+    status = initialize_dynamic_branches!(x0_init, inputs)
+    return status
+end
+function get_jacobian(T, system::PSY.System)
+    # Deepcopy avoid system modifications
+    simulation_system = deepcopy(system)
+    inputs = SimulationInputs(T, simulation_system, ReferenceBus)
+    x0_init = get_flat_start(inputs)
+    _set_operating_point!(x0_init, inputs, system)
+    return get_jacobian(T, inputs, x0_init)
+end

src/base/nlsolve_wrapper.jl

@@ -71,6 +71,7 @@ function _convergence_check(sys_solve::NLsolveWrapper, tol::Float64, solv::Symbo
 end
 
 function refine_initial_condition!(
+    initial_guess::Vector{Float64},
     sim::Simulation,
     model::SystemModel,
     jacobian::JacobianFunctionWrapper,
@@ -84,7 +85,6 @@ function refine_initial_condition!(
 
     @debug "Start NLSolve System Run"
     converged = false
-    initial_guess = get_initial_conditions(sim)
     inputs = get_simulation_inputs(sim)
     bus_range = get_bus_range(inputs)
     powerflow_solution = deepcopy(initial_guess[bus_range])