Enable sensitivity analysis

Project.toml

@@ -4,40 +4,56 @@ authors = ["Jose Daniel Lara, Rodrigo Henriquez"]
 version = "0.15.0"
 
 [deps]
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+ChainRulesCore = "1.23"
+ComponentArrays = "0.15"
 DataFrames = "1"
 DataStructures = "~0.18"
+DiffEqBase = "6.151"
 DocStringExtensions = "~0.9"
+Enzyme = "0.12"
 FastClosures = "^0.3"
 ForwardDiff = "~v0.10"
 InfrastructureSystems = "2"
 NLsolve = "4"
+NonlinearSolve = "3.11"
 PowerSystems = "4"
 PowerFlows = "^0.7"
 PowerNetworkMatrices = "^0.11"
 PrettyTables = "1, 2"
 SciMLBase = "2"
+SciMLSensitivity = "7.0"
+
 TimerOutputs = "~0.5"
 LinearAlgebra = "1"
 Logging = "1"
+
 Random = "1"
 SparseArrays = "1"
+Zygote = "0.6"
 julia = "^1.6"

docs/make.jl

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

docs/src/sensitivity_analysis.md

@@ -0,0 +1,9 @@
+# Sensitivity Analysis 
+
+PowerSimulationsDynamics has limited support for performing sensitivity analysis of system parameters with respect to a user defined loss function. See the tutorial for an example of sensitivity analysis. 
+
+* `ForwardDiffSensitivity()` is used as the method for differentiating through the solve. See `SciMLSensitivity.jl` for more details. 
+* The gradient of the function provided by PSID can be calculated using `Zygote.jl`.
+* The Jacobian is not passed to the ODE Problem during sensitivity analysis. 
+* Parameters for sensitivity analysis must not change the steady state operating condition. Parameters cannot be in the mass matrix. This limitation is expected to be relaxed in the future. 
+* Limited to mass matrix formulations and pure Julia solvers. Can check if a solver is compatible with `SciMLBase.isautodifferentiable()`

src/PowerSimulationsDynamics.jl

@@ -61,32 +61,46 @@ export is_valid
 export transform_load_to_constant_impedance
 export transform_load_to_constant_current
 export transform_load_to_constant_power
+export get_parameter_values
+export get_parameter_labels
+export get_sensitivity_functions
 
 ####################################### Package Imports ####################################
 import Logging
 import InfrastructureSystems
 import SciMLBase
+import SciMLSensitivity
 import DataStructures: OrderedDict
 import DataFrames: DataFrame
+import DiffEqBase
 import Random
 import ForwardDiff
 import SparseArrays
 import LinearAlgebra
 import Base.to_index
-import NLsolve
+import Base.length
+import NLsolve  #Using Nlsovle wrapped using NonlinearSolve interfaace
+import NonlinearSolve
 import PrettyTables
 import Base.ImmutableDict
 import PowerSystems
 import PowerFlows
 import PowerNetworkMatrices
 import TimerOutputs
 import FastClosures: @closure
+import Enzyme
+Enzyme.API.runtimeActivity!(true)  #Needed for "activity unstable" code: https://enzymead.github.io/Enzyme.jl/stable/faq/
+Enzyme.API.looseTypeAnalysis!(true)  #Required for using component arrays with Enzyme
+import ChainRulesCore
+import ComponentArrays
+import Zygote
 
 const PSY = PowerSystems
 const IS = InfrastructureSystems
 const PSID = PowerSimulationsDynamics
 const PF = PowerFlows
 const PNM = PowerNetworkMatrices
+const CRC = ChainRulesCore
 
 using DocStringExtensions
 
@@ -104,6 +118,7 @@ include("base/device_wrapper.jl")
 include("base/branch_wrapper.jl")
 include("base/frequency_reference.jl")
 include("base/simulation_model.jl")
+include("utils/parameters.jl")
 include("base/simulation_inputs.jl")
 include("base/perturbations.jl")
 include("base/caches.jl")
@@ -179,6 +194,10 @@ include("post_processing/post_proc_results.jl")
 include("post_processing/post_proc_loads.jl")
 include("post_processing/post_proc_source.jl")
 
+#Sensitivity Analysis
+include("base/sensitivity_rule.jl") #Custom rule to support duplicated kwargs for callbacks. See: https://github.com/EnzymeAD/Enzyme.jl/issues/1491
+include("base/sensitivity_analysis.jl")
+
 #Utils
 include("utils/psy_utils.jl")
 include("utils/immutable_dicts.jl")

src/base/branch_wrapper.jl

@@ -2,11 +2,11 @@
 Wraps DynamicBranch devices from PowerSystems to handle changes in controls and connection
 status, and allocate the required indexes of the state space.
 """
-struct BranchWrapper
+mutable struct BranchWrapper
     branch::PSY.DynamicBranch
     system_base_power::Float64
     system_base_frequency::Float64
-    connection_status::Base.RefValue{Float64}
+    connection_status::Float64
     bus_ix_from::Int
     bus_ix_to::Int
     ix_range::Vector{Int}
@@ -26,7 +26,7 @@ struct BranchWrapper
             branch,
             sys_base_power,
             sys_base_freq,
-            Base.Ref(1.0),
+            1.0,
             bus_ix_from,
             bus_ix_to,
             ix_range,

src/base/caches.jl

@@ -70,11 +70,11 @@ get_global_vars(jc::JacobianCache, ::Type{T}) where {T <: ForwardDiff.Dual} =
 struct SimCache{F} <: Cache
     f!::F
     bus_count::Int
-    ode_output::Vector{Float64}
-    branches_ode::Vector{Float64}
-    current_balance::Vector{Float64}
-    inner_vars::Vector{Float64}
-    global_vars::Vector{Float64}
+    ode_output::Vector{Real}
+    branches_ode::Vector{Real}
+    current_balance::Vector{Real}
+    inner_vars::Vector{Real}
+    global_vars::Vector{Real}
 end
 
 function SimCache(f!, inputs::SimulationInputs)
@@ -84,14 +84,19 @@ function SimCache(f!, inputs::SimulationInputs)
     bus_count = get_bus_count(inputs)
     inner_vars_count = get_inner_vars_count(inputs)
     n_global_vars = length(keys(get_global_vars_update_pointers(inputs)))
+    global_vars = setindex!(
+        zeros(Real, n_global_vars),
+        1.0,
+        GLOBAL_VAR_SYS_FREQ_INDEX,
+    )
     return SimCache{typeof(f!)}(
         f!,
         bus_count,
-        zeros(Float64, n_inj),
-        zeros(Float64, n_branches),
-        zeros(Float64, 2 * bus_count),
-        zeros(Float64, inner_vars_count),
-        setindex!(zeros(Float64, n_global_vars), 1.0, GLOBAL_VAR_SYS_FREQ_INDEX),
+        zeros(Real, n_inj),
+        zeros(Real, n_branches),
+        zeros(Real, 2 * bus_count),
+        zeros(Real, inner_vars_count),
+        global_vars,
     )
 end
 
@@ -103,11 +108,11 @@ function get_current_injections_i(sc::SimCache, ::Type{Float64})
     return view(sc.current_balance, ((sc.bus_count + 1):(sc.bus_count * 2)))
 end
 
-get_ode_output(sc::SimCache, ::Type{Float64}) = sc.ode_output
-get_branches_ode(sc::SimCache, ::Type{Float64}) = sc.branches_ode
-get_current_balance(sc::SimCache, ::Type{Float64}) = sc.current_balance
-get_inner_vars(sc::SimCache, ::Type{Float64}) = sc.inner_vars
-get_global_vars(sc::SimCache, ::Type{Float64}) = sc.global_vars
+get_ode_output(sc::SimCache, ::Type{T}) where {T <: Real} = sc.ode_output
+get_branches_ode(sc::SimCache, ::Type{T}) where {T <: Real} = sc.branches_ode
+get_current_balance(sc::SimCache, ::Type{T}) where {T <: Real} = sc.current_balance
+get_inner_vars(sc::SimCache, ::Type{T}) where {T <: Real} = sc.inner_vars
+get_global_vars(sc::SimCache, ::Type{T}) where {T <: Real} = sc.global_vars
 
 get_ω_sys(cache::Cache, T::Type{<:ACCEPTED_REAL_TYPES}) =
     get_global_vars(cache, T)[GLOBAL_VAR_SYS_FREQ_INDEX]

src/base/definitions.jl

@@ -126,7 +126,6 @@ get_vars_ix() = Dict{Int, Int}(GLOBAL_VAR_SYS_FREQ_INDEX => -1)
 
 const RELAXED_NLSOLVE_F_TOLERANCE = :1e-6
 const STRICT_NLSOLVE_F_TOLERANCE = :1e-9
-const NLSOLVE_X_TOLERANCE = :1e-9
 const MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE = :1e-4
 const MAX_INIT_RESIDUAL = 10.0
 const MAX_NLSOLVE_INTERATIONS = 10
@@ -212,7 +211,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
@@ -223,6 +222,16 @@ Defines the status of the simulation object
     SIMULATION_FAILED = 7
 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}
+const ACCEPTED_REAL_TYPES = Union{Float64, ForwardDiff.Dual, Real}