[WIP] Update initialization procedure

src/LITS.jl

@@ -17,6 +17,8 @@ export small_signal_analysis
 export get_state_series
 export get_voltagemag_series
 export print_init_states
+export print_device_states
+export get_dict_init_states
 
 ####################################### Package Imports ####################################
 import Logging
@@ -38,11 +40,11 @@ include("base/ports.jl")
 include("base/perturbations.jl")
 include("base/small_signal_results.jl")
 include("base/simulation_initialization.jl")
+include("base/file_system.jl")
 include("base/simulation.jl")
 
 #Common Models
 include("models/branch.jl")
-include("models/initialize_device.jl")
 include("models/device.jl")
 include("models/kirchoff_laws.jl")
 include("models/dynline_model.jl")
@@ -69,12 +71,22 @@ include("models/source_models.jl")
 
 #Initialization Parameters
 include("initialization/init_device.jl")
+
+#Initialization Generator
 include("initialization/generator_components/init_machine.jl")
 include("initialization/generator_components/init_shaft.jl")
 include("initialization/generator_components/init_avr.jl")
 include("initialization/generator_components/init_tg.jl")
 include("initialization/generator_components/init_pss.jl")
 
+#Initialization Inverter
+include("initialization/inverter_components/init_filter.jl")
+include("initialization/inverter_components/init_DCside.jl")
+include("initialization/inverter_components/init_converter.jl")
+include("initialization/inverter_components/init_frequency_estimator.jl")
+include("initialization/inverter_components/init_inner.jl")
+include("initialization/inverter_components/init_outer.jl")
+
 #System Model
 include("models/system.jl")
 

src/base/definitions.jl

@@ -66,6 +66,7 @@ Vi_cnv_var :: Voltage supplied from the converter in the I-component
     VI_inv_var = 11
     Vr_cnv_var = 12
     Vi_cnv_var = 13
+    P_ES_var = 14 #Power from energy source
 end
 
 Base.to_index(ix::inverter_inner_vars) = Int(ix)
@@ -100,8 +101,17 @@ const VOLTAGE_BUSES_NO = "voltage_buses_no"
 const CURRENT_BUSES_NO = "current_buses_no"
 const TOTAL_SHUNTS = "total_shunts"
 const AUX_ARRAYS = "aux_arrays"
+const LOOKUP = "lookup"
 
-const SIMULATION_ACCEPTED_KWARGS = [:initial_guess, :initialize_simulation]
+const SIMULATION_ACCEPTED_KWARGS = [:initial_guess, :initialize_simulation, :system_to_file]
 
 PSY.get_V_ref(value::PSY.AVRFixed) = value.Vf
-PSY.set_V_ref!(value::PSY.AVRFixed, val::Float64) = value.Vf = val
+PSY.set_V_ref!(value::PSY.AVRFixed, val::Float64) = value.Vf = val
+
+get_V_ref_control(value::PSY.DynamicGenerator) = PSY.get_V_ref(PSY.get_avr(value))
+get_V_ref_control(value::PSY.DynamicInverter) =
+    PSY.get_V_ref(PSY.get_reactive_power(PSY.get_outer_control(value)))
+
+get_P_ref_control(value::PSY.DynamicGenerator) = PSY.get_P_ref(PSY.get_prime_mover(value))
+get_P_ref_control(value::PSY.DynamicInverter) =
+    PSY.get_P_ref(PSY.get_active_power(PSY.get_outer_control(value)))

src/base/file_system.jl

@@ -0,0 +1,9 @@
+function check_folder(folder::String)
+    !isdir(folder) && throw(IS.ConflictingInputsError("Specified folder is not valid"))
+    try
+        mkdir(joinpath(folder, "fake"))
+        rm(joinpath(folder, "fake"))
+    catch e
+        throw(IS.ConflictingInputsError("Specified folder does not have write access [$e]"))
+    end
+end

src/base/perturbations.jl

@@ -1,23 +1,25 @@
 abstract type Perturbation end
 
-struct ThreePhaseFault <: Perturbation
+mutable struct ThreePhaseFault <: Perturbation
     time::Float64
     Ybus::SparseMatrixCSC{Complex{Float64}, Int64}
 end
 
-get_affect(pert::ThreePhaseFault) =
+get_affect(::PSY.System, pert::ThreePhaseFault) =
     (integrator) -> PSY.get_ext(integrator.p)["Ybus"] = pert.Ybus
 
-struct ControlReferenceChange <: Perturbation
+mutable struct ControlReferenceChange <: Perturbation
     time::Float64
     device::PSY.DynamicInjection
     signal_index::Int64
     ref_value::Float64
 end
 
-function get_affect(pert::ControlReferenceChange)
+function get_affect(system::PSY.System, pert::ControlReferenceChange)
+    device = PSY.get_component(typeof(pert.device), system, PSY.get_name(pert.device))
+    pert.device = device
     return (integrator) -> begin
-        control_ref = PSY.get_ext(pert.device)[CONTROL_REFS]
+        control_ref = PSY.get_ext(device)[CONTROL_REFS]
         return control_ref[pert.signal_index] = pert.ref_value
     end
 end

src/base/simulation.jl

@@ -8,48 +8,57 @@ mutable struct Simulation
     tstops::Vector{Float64}
     callbacks::DiffEqBase.CallbackSet
     solution::Union{Nothing, DiffEqBase.DAESolution}
+    simulation_folder::String
     ext::Dict{String, Any}
 end
 
 function Simulation(
+    simulation_folder::String,
     system::PSY.System,
     tspan::NTuple{2, Float64},
     perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
-    initialize_simulation::Bool = true,
     kwargs...,
 )
+    check_folder(simulation_folder)
+    simulation_system = deepcopy(system)
     check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
     initialized = false
-    DAE_vector = _index_dynamic_system!(system)
-    var_count = get_variable_count(system)
+    DAE_vector = _index_dynamic_system!(simulation_system)
+    var_count = get_variable_count(simulation_system)
 
     flat_start = zeros(var_count)
-    bus_count = length(PSY.get_components(PSY.Bus, system))
+    bus_count = length(PSY.get_components(PSY.Bus, simulation_system))
     flat_start[1:bus_count] .= 1.0
     x0_init = get(kwargs, :initial_guess, flat_start)
 
+    initialize_simulation = get(kwargs, :initialize_simulation, true)
     if initialize_simulation
         @info("Initializing Simulation States")
-        _add_aux_arrays!(system, Real)
-        initialized = calculate_initial_conditions!(system, x0_init)
+        _add_aux_arrays!(simulation_system, Real)
+        initialized = calculate_initial_conditions!(simulation_system, x0_init)
     end
 
     dx0 = zeros(var_count)
-    callback_set, tstops = _build_perturbations(perturbations::Vector{<:Perturbation})
-
-    _add_aux_arrays!(system, Float64)
+    callback_set, tstops = _build_perturbations(simulation_system, perturbations)
+    _add_aux_arrays!(simulation_system, Float64)
     prob = DiffEqBase.DAEProblem(
         system!,
         dx0,
         x0_init,
         tspan,
-        system,
+        simulation_system,
         differential_vars = DAE_vector;
         kwargs...,
     )
 
+    if get(kwargs, :system_to_file, false)
+        PSY.to_json(
+            simulation_system,
+            joinpath(simulation_folder, "initialized_system.json"),
+        )
+    end
     return Simulation(
-        system,
+        simulation_system,
         false,
         prob,
         perturbations,
@@ -58,18 +67,21 @@ function Simulation(
         tstops,
         callback_set,
         nothing,
+        simulation_folder,
         Dict{String, Any}(),
     )
 end
 
 function Simulation(
+    simulation_folder::String,
     system::PSY.System,
     tspan::NTuple{2, Float64},
     perturbation::Perturbation;
     initialize_simulation::Bool = true,
     kwargs...,
 )
     return Simulation(
+        simulation_folder,
         system,
         tspan,
         [perturbation];
@@ -86,20 +98,20 @@ function _add_aux_arrays!(system::PSY.System, T)
         3 => collect(zeros(T, get_n_injection_states(system))),  #injection_ode
         4 => collect(zeros(T, get_n_branches_states(system))),   #branches_ode
         5 => collect(zeros(Complex{T}, bus_count)),              #I_bus
-        6 => collect(zeros(T, 2 * bus_count)),                      #I_balance
+        6 => collect(zeros(T, 2 * bus_count)),                   #I_balance
     )
     system.internal.ext[AUX_ARRAYS] = aux_arrays
     return
 end
 
-function _build_perturbations(perturbations::Vector{<:Perturbation})
+function _build_perturbations(system::PSY.System, perturbations::Vector{<:Perturbation})
     isempty(perturbations) && return DiffEqBase.CallbackSet(), [0.0]
     perturbations_count = length(perturbations)
     callback_vector = Vector{DiffEqBase.DiscreteCallback}(undef, perturbations_count)
     tstops = Vector{Float64}(undef, perturbations_count)
     for (ix, pert) in enumerate(perturbations)
         condition = (x, t, integrator) -> t in [pert.time]
-        affect = get_affect(pert)
+        affect = get_affect(system, pert)
         callback_vector[ix] = DiffEqBase.DiscreteCallback(condition, affect)
         tstops[ix] = pert.time
     end
@@ -136,16 +148,16 @@ function _attach_inner_vars!(
     device::PSY.DynamicInverter,
     ::Type{T} = Real,
 ) where {T <: Real}
-    device.ext[INNER_VARS] = zeros(T, 13)
+    device.ext[INNER_VARS] = zeros(T, 14)
     return
 end
 
 function _attach_control_refs!(device::PSY.DynamicInjection)
     device_basepower = PSY.get_basepower(device)
     device.ext[CONTROL_REFS] = [
-        PSY.get_V_ref(PSY.get_avr(device)),
+        get_V_ref_control(device),
         PSY.get_ω_ref(device),
-        PSY.get_P_ref(PSY.get_prime_mover(device)),
+        get_P_ref_control(device),
         PSY.get_reactivepower(PSY.get_static_injector(device)),
     ]
     return
@@ -284,17 +296,20 @@ function _index_dynamic_system!(sys::PSY.System)
         total_states += device_n_states
         _add_states_to_global!(global_state_index, state_space_ix, d)
         btype != PSY.BusTypes.REF && continue
-        global_vars[:ω_sys_index] = global_state_index[d.name][:ω] #To define 0 if infinite source, bus_number otherwise,
+        global_vars[:ω_sys_index] = global_state_index[PSY.get_name(d)][:ω] #To define 0 if infinite source, bus_number otherwise,
         found_ref_bus = true
     end
     injection_n_states = state_space_ix[1] - branches_n_states - n_buses * 2
     @assert total_states == state_space_ix[1] - static_bus_var_count
     @debug total_states
     setdiff!(current_buses_no, voltage_buses_no)
     if !isempty(PSY.get_components(PSY.ACBranch, sys))
-        Ybus = PSY.Ybus(sys)[:, :]
+        Ybus_ = PSY.Ybus(sys)
+        Ybus = Ybus_[:, :]
+        lookup = Ybus_.lookup[1]
     else
         Ybus = SparseMatrixCSC{Complex{Float64}, Int64}(zeros(n_buses, n_buses))
+        lookup = Dict{Int.Int}()
     end
     sys_ext = Dict{String, Any}()
     counts = Base.ImmutableDict(
@@ -307,6 +322,7 @@ function _index_dynamic_system!(sys::PSY.System)
         :bus_count => n_buses,
     )
 
+    sys_ext[LOOKUP] = lookup
     sys_ext[LITS_COUNTS] = counts
     sys_ext[GLOBAL_INDEX] = global_state_index
     sys_ext[VOLTAGE_BUSES_NO] = voltage_buses_no
@@ -335,6 +351,7 @@ get_current_bus_no(sys::PSY.System) = PSY.get_ext(sys)[CURRENT_BUSES_NO]
 get_voltage_bus_no(sys::PSY.System) = PSY.get_ext(sys)[VOLTAGE_BUSES_NO]
 get_total_shunts(sys::PSY.System) = PSY.get_ext(sys)[TOTAL_SHUNTS]
 get_bus_count(sys::PSY.System) = PSY.get_ext(sys)[LITS_COUNTS][:bus_count]
+get_lookup(sys::PSY.System) = PSY.get_ext(sys)[LOOKUP]
 
 function _get_internal_mapping(
     device::PSY.DynamicInjection,
@@ -383,8 +400,9 @@ function _change_vector_type(sys::PSY.System)
 end
 
 function _determine_stability(vals::Vector{Complex{Float64}})
+    stable = true
     for real_eig in real(vals)
-        real_eig >= 0.0 && return false
+        real_eig > 0.0 && return false
     end
     return true
 end
@@ -423,6 +441,15 @@ function small_signal_analysis(sim::Simulation; kwargs...)
     # TODO: Make operation using BLAS!
     reduced_jacobian = fx - fy * inv(gy) * gx
     vals, vect = LinearAlgebra.eigen(reduced_jacobian)
+    sources = collect(PSY.get_components(PSY.Source, sim.system))
+    if isempty(sources)
+        @warn("No Infinite Bus found. Confirm stability directly checking eigenvalues.\nIf all eigenvalues are on the left-half plane and only one eigenvalue is zero, the system is small signal stable.")
+        info_evals = "Eigenvalues are:\n"
+        for i in vals
+            info_evals = info_evals * string(i) * "\n"
+        end
+        @info(info_evals)
+    end
     return SmallSignalOutput(
         reduced_jacobian,
         vals,

src/base/simulation_initialization.jl

@@ -15,37 +15,35 @@ function calculate_initial_conditions!(sys::PSY.System, initial_guess::Vector{Fl
     branches_start = get_branches_pointer(sys)
     branches_count = 1
 
+    for bus in PSY.get_components(PSY.Bus, sys)
+        #Write voltage initial guess
+        bus_n = PSY.get_number(bus)
+        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
+        initial_guess[bus_ix] = PSY.get_voltage(bus) * cos(PSY.get_angle(bus))
+        initial_guess[bus_ix + bus_size] = PSY.get_voltage(bus) * sin(PSY.get_angle(bus))
+    end
+
     for d in PSY.get_components(PSY.DynamicInjection, sys)
         bus = PSY.get_bus(d)
-        bus_n = PSY.get_number(PSY.get_bus(d)) # TODO: This requires that the bus numbers are indexed 1-N
+        bus_n = PSY.get_number(PSY.get_bus(d))
+        bus_ix = PSY.get_ext(sys)[LOOKUP][bus_n]
         n_states = PSY.get_n_states(d)
         ix_range = range(injection_start, length = n_states)
-        #injection_start = injection_start + n_states
-        #Write Buses Voltages Initial Guess
-        initial_guess[bus_n] = PSY.get_voltage(bus)*cos(PSY.get_angle(bus))
-        initial_guess[bus_n + bus_size] = PSY.get_voltage(bus)*sin(PSY.get_angle(bus))
-        initialize_device!(
-            initial_guess,
-            ix_range,
-            d,
-        )
+        injection_start = injection_start + n_states
+        x0_device = initialize_device(d)
+        @assert length(x0_device) == n_states
+        initial_guess[ix_range] = x0_device
     end
 
-
-
     dx0 = zeros(var_count) #Define a vector of zeros for the derivative
-
     inif! = (out, x) -> system!(
         out,    #output of the function
         dx0,    #derivatives equal to zero
         x,      #states
         sys,    #Parameters
         0.0,    #time equals to zero.
     )
-
-
-
-
+    #Refine initial solution
     sys_solve = NLsolve.nlsolve(
         inif!,
         initial_guess,
@@ -54,7 +52,19 @@ function calculate_initial_conditions!(sys::PSY.System, initial_guess::Vector{Fl
         method = :trust_region,
     ) #Solve using initial guess x0
     if !NLsolve.converged(sys_solve)
-        @warn("Initialization failed, initial conditions do not meet conditions for an stable equilibrium")
+        @warn("Initialization failed, initial conditions do not meet conditions for an stable equilibrium.\nTrying to solve again reducing numeric tolerance:")
+        sys_solve = NLsolve.nlsolve(
+            inif!,
+            initial_guess,
+            xtol = :1e-6,
+            ftol = :1e-6,
+            method = :trust_region,
+        ) #Solve using initial guess x0
+        if NLsolve.converged(sys_solve)
+            @info("Initialization succeeded with less tolerance. Saving solution")
+        else
+            @warn("Initialization failed again. Initial conditions do not meet conditions for an stable equilibrium\nSaving best result.")
+        end
     end
     initial_guess .= sys_solve.zero
     return NLsolve.converged(sys_solve)