Jd/smallsignal

Project.toml

@@ -4,6 +4,7 @@ version = "0.3.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -13,6 +14,7 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [compat]
 DiffEqBase = "6"
+ForwardDiff = "~v0.10"
 NLsolve = "4"
 PowerSystems = "~0.8"
 Sundials = "3"

src/LITS.jl

@@ -8,13 +8,19 @@ export Simulation
 export run_simulation!
 export ThreePhaseFault
 export ControlReferenceChange
+
+# Export for routines
+export small_signal_analysis
 export get_state_series
 export get_voltagemag_series
 export print_init_states
+
 ####################################### Package Imports ####################################
 import DiffEqBase
+import ForwardDiff
 import SparseArrays: SparseMatrixCSC
 import LinearAlgebra: BLAS
+import LinearAlgebra: eigen
 import Base.to_index
 import NLsolve
 import PowerSystems
@@ -23,7 +29,8 @@ const PSY = PowerSystems
 #Structs for General Devices and System
 include("base/definitions.jl")
 include("base/ports.jl")
-include("perturbations/perturbations.jl")
+include("base/perturbations.jl")
+include("base/small_signal_results.jl")
 include("base/simulation.jl")
 
 #Common Models
@@ -55,11 +62,6 @@ include("models/source_models.jl")
 #System Model
 include("models/system_model.jl")
 
-#Perturbations
-include("perturbations/common.jl")
-include("perturbations/ThreePhaseFault.jl")
-include("perturbations/PowerStepChange.jl")
-
 #Utils
 include("utils/plot_utils.jl")
 include("utils/print.jl")

src/perturbations/perturbations.jl → src/base/perturbations.jl

@@ -2,7 +2,7 @@ abstract type Perturbation end
 
 struct ThreePhaseFault <: Perturbation
     time::Float64
-    Ybus #::SparseMatrixCSC{Complex{Float64}, Int64}
+    Ybus::SparseMatrixCSC{Complex{Float64},Int64}
 end
 
 get_affect(pert::ThreePhaseFault) =
@@ -17,7 +17,7 @@ end
 
 function get_affect(pert::ControlReferenceChange)
     return (integrator) -> begin
-        return PSY.get_ext(pert.device)[CONTROL_REFS][pert.signal_index] =
-                pert.ref_value
+        control_ref = PSY.get_ext(pert.device)[CONTROL_REFS]
+        return control_ref[pert.signal_index] = pert.ref_value
     end
 end

src/base/simulation.jl

@@ -114,21 +114,6 @@ function _calculate_initial_conditions(sys::PSY.System, initial_guess::Vector{Fl
     return sys_solve.zero, NLsolve.converged(sys_solve)
 end
 
-function run_simulation!(sim::Simulation, solver; kwargs...)
-    if sim.reset
-        @error("Reset the simulation")
-    end
-
-    sim.solution = DiffEqBase.solve(
-        sim.problem,
-        solver;
-        callback = sim.callbacks,
-        tstops = sim.tstops,
-        kwargs...,
-    )
-    return
-end
-
 function _index_local_states!(
     component_state_index::Vector{Int64},
     local_states::Vector{Symbol},
@@ -145,13 +130,19 @@ function _attach_ports!(component::PSY.DynamicComponent)
     return
 end
 
-function _attach_inner_vars!(device::PSY.DynamicGenerator)
-    device.ext[INNER_VARS] = zeros(Real, 8)
+function _attach_inner_vars!(
+    device::PSY.DynamicGenerator,
+    ::Type{T} = Float64,
+) where {T<:Real}
+    device.ext[INNER_VARS] = zeros(T, 8)
     return
 end
 
-function _attach_inner_vars!(device::PSY.DynamicInverter)
-    device.ext[INNER_VARS] = zeros(Real, 13)
+function _attach_inner_vars!(
+    device::PSY.DynamicInverter,
+    ::Type{T} = Float64,
+) where {T<:Real}
+    device.ext[INNER_VARS] = zeros(T, 13)
     return
 end
 
@@ -318,3 +309,52 @@ function get_input_port_ix(
 ) where {T<:PSY.DynamicComponent}
     return _get_internal_mapping(device, INPUT_PORT_MAPPING, ty)
 end
+
+function run_simulation!(sim::Simulation, solver; kwargs...)
+    if sim.reset
+        @error("Reset the simulation")
+    end
+
+    sim.solution = DiffEqBase.solve(
+        sim.problem,
+        solver;
+        callback = sim.callbacks,
+        tstops = sim.tstops,
+        kwargs...,
+    )
+    return
+end
+
+function _change_vector_type(sys::PSY.System)
+    for d in PSY.get_components(PSY.DynamicInjection, sys)
+        _attach_inner_vars!(d, Real)
+    end
+end
+
+function _determine_stability(vals::Vector{Complex{Float64}})
+    for real_eig in real(vals)
+        real_eig >= 0.0 && return false
+    end
+    return true
+end
+
+function small_signal_analysis(sim::Simulation; kwargs...)
+    if sim.reset
+        @error("Reset the simulation")
+    end
+    _change_vector_type(sim.system)
+    var_count = LITS.get_variable_count(sim.system)
+    dx0 = zeros(var_count) #Define a vector of zeros for the derivative
+    sysf! = (out, x) -> LITS.system_model!(
+        out,            #output of the function
+        dx0,            #derivatives equal to zero
+        x,              #states
+        sim.system,     #Parameters
+        0.0,            #time equals to zero.
+    )
+    out = zeros(var_count) #Define a vector of zeros for the output
+    x_eval = get(kwargs, :operating_point, sim.x0_init)
+    res = ForwardDiff.jacobian(sysf!, out, x_eval)
+    vals, vect = eigen(res)
+    return SmallSignalOutput(res, vals, vect, _determine_stability(vals), x_eval)
+end

src/base/small_signal_results.jl

@@ -0,0 +1,7 @@
+struct SmallSignalOutput
+    jacobian::Matrix{Float64}
+    eigenvalues::Vector{Complex{Float64}}
+    eigenvectors::Matrix{Complex{Float64}}
+    stable::Bool
+    operating_point::Vector{Float64}
+end

src/utils/print.jl

@@ -1,3 +1,8 @@
-function Base.show(io::IO, simulation::Simulation)
+function Base.show(io::IO, ::Simulation)
     println(io, "Simulation()")
 end
+
+function Base.show(io::IO, smr::SmallSignalOutput)
+    val = smr.stable ? "is" : "is not"
+    println(io, "The system $(val) small signal stable")
+end