Dynamic.jl

assemble_f_dynamic(
    cm::ContGridModML.ContModel,
    dm::ContGridModML.DiscModel,
    fault_ix::Vector{<:Integer},
    dP::Union{Real, Vector{<:Real}},
    Af::SparseArrays.SparseMatrixCSC
) -> Matrix{Float64}

Assemble all the force vectors for the dynamical simulations.

assemble_matrices_dynamic(
    model::ContGridModML.ContModel
) -> NTuple{4, SparseArrays.SparseMatrixCSC{Float64, Int64}}

Create all the necessary finite element matrices from a given model.

The matrices returned are

• M_const the constant (non-parameter dependent) part of the mass matrix
• K_const the constant (non-parameter dependent) part of the stiffness matrix
• A the matrix needed to create the non-constant part of both the mass and stiffness matrices. For example the mass matrix can be calculated as M = M_const + A * m_quad * A'
• Af the matrix needed to create the force vector as f = Af * (p_quad + fault_quad)

cont_dyn(
    M::SparseArrays.SparseMatrixCSC,
    K::SparseArrays.SparseMatrixCSC,
    f::Vector{<:Real},
    u₀::Vector{<:Real},
    tf::Real;
    solve_kwargs
) -> SciMLBase.ODESolution

Run a dynamical simulation of the continuous model.

disc_dyn(
    dm::ContGridModML.DiscModel,
    fault_node::Integer,
    fault_size::Real,
    dt::Real,
    tf::Real;
    solve_kwargs
) -> SciMLBase.ODESolution

Run a dynamical simulation of the discrete model.

gen_idxs(
    cm::ContGridModML.ContModel,
    dm::ContGridModML.DiscModel,
    dP::Real,
    n_train::Integer,
    n_test::Integer;
    seed
) -> Tuple{Vector{<:Integer}, Vector{<:Integer}}

Randomly choose generators from a homoegenous distribution for the training and test sets.

generate_comp_idxs(
    cm::ContGridModML.ContModel,
    dm::ContGridModML.DiscModel,
    tri::Vector{<:Integer},
    tei::Vector{<:Integer},
    n::Int64
) -> Vector{<:Integer}

Generate a list of indices used for the comparison of the ground truth with the results from the continuous model.

An equally spaced grid is overlaid over the area. If the points are within the area, the closest bus in the discrete model is found and added to the list of indices. The generators of the test and training set are not eligible for comparison and are remove from the list of possible indices. The number of points can be controlled by n, which gives the number of points in the largest dimension.

grad(
    sol_cont::SciMLBase.ODESolution,
    sol_lambda::SciMLBase.ODESolution,
    g_proj::Vector{SparseArrays.SparseMatrixCSC}
) -> Vector{Float64}
grad(
    sol_cont::SciMLBase.ODESolution,
    sol_lambda::SciMLBase.ODESolution,
    g_proj::Vector{SparseArrays.SparseMatrixCSC},
    dt::Real
) -> Vector{Float64}

Calculate the gradient from the solution of the adjoint dynamics.

grad_proj(
    A::SparseArrays.SparseMatrixCSC,
    q_proj::SparseArrays.SparseMatrixCSC,
    evecs::Matrix{<:Real},
    n_coeffs::Integer
) -> Vector{SparseArrays.SparseMatrixCSC}

Calculate the projection matrix of the discrete solution onto the adjoint solution.

init_expansion(
    cm::ContGridModML.ContModel,
    eve::Matrix{<:Real},
    n_modes::Integer,
    n_coeffs::Integer,
    m::Vector{<:Real},
    d::Vector{<:Real}
) -> Tuple{Vector{T} where T<:Real, Matrix{<:Real}}

Create the initial values for the parameters.

The parameters are expanded in eigenvectors of the grid Laplacian. The first n_modes eigenvectors are chosen. The first n_coeffs coefficients are chosen by projecting the results of the heat equation diffusion onto the eigenvectors. The n_modes - n_coeffs are set to zero.

initial_conditions(
    cm::ContGridModML.ContModel
) -> Vector{T} where T<:Real

Create the initial conditions for the continuous simulations.

lambda_dyn(
    cont_sol::SciMLBase.ODESolution,
    disc_sol::SciMLBase.ODESolution,
    M::SparseArrays.SparseMatrixCSC,
    K::SparseArrays.SparseMatrixCSC,
    ω_proj::SparseArrays.SparseMatrixCSC,
    tf::Real,
    idxs::Vector{<:Integer};
    solve_kwargs
) -> SciMLBase.ODESolution

Solve the adjoint dynamics.

The continuous and discrete solutions are needed as well as the comparison indices to calculate the contributions from the loss function.

lap_eigenvectors(
    cm::ContGridModML.ContModel
) -> Matrix{Float64}

Calculate the eigenvectors of the unweighted Laplacian of the grid used for the continuous simulations.

learn_dynamical_parameters(
;
    dm_fn,
    cm_fn,
    dP,
    n_train,
    n_test,
    dt,
    tf,
    disc_solve_kwargs,
    cont_solve_kwargs,
    lambda_solve_kwargs,
    seed,
    n_points,
    m_init,
    d_init,
    n_coeffs,
    n_modes,
    n_epochs,
    max_function,
    train_ix,
    test_ix,
    n_batches
)

Learn the inertia and damping distribution for the continuous model.

The frequency response of faults at multiple generators are compared on homogeneously spread buses across the grid. The gradient is calculated using an adjoint sensitivity method and the updates to the parameters are calculated using a constraint gradient descent method.

Arguments

• dm_fn::String = MODULE_FOLDER * "/data/dm.h5": File name of the discrete model
• cm_fn::String = MODULE_FOLDER * "/data/cm.h5": File name of the continuous model
• dP::Real = -9.0: Fault size to be simulated
• n_train::Integer = 12: Amount of faults to consider for training
• n_test::Integer = 4: Amount of faults to consider for testing
• dt::Real = 0.01: Step size at which the solutions of the ODEs are saved
• tf::Real = 25.0: Duration of the simulations
• disc_solve_kwargs::Dict{Symbol, <:Any} = Dict{Symbol, Any}(): Keyword arguments passed to the ODE solver for the discrete model
• cont_solve_kwargs::Dict{Symbol, <:Any} = Dict{Symbol, Any}(): Keyword arguments passed to the ODE solver for the continuous model
• lambda_solve_kwargs::Dict{Symbol, <:Any} = Dict{Symbol, Any}(:saveat => 0.1, :abstol => 1e-3, :reltol => 1e-2): Keyword arguments passed to the ODE solver of the adjoint equations
• seed::Union{Nothing, Integer} = 1709: Seed for the random number generator to be used to pick the training and test generators
• σ = 0.05: Standard deviation of the Gaussian used to distribute the fault
• n_points = 40: Number of points for comparison along the largest dimension
• n_coeffs = 1: Number of coefficients that are non-zero at the beginning of the training. They correspond to the n_coeffs lowest modes of the Laplacian.
• n_modes = 20: Number of modes the Laplacian that are used to expand the parameters.
• n_epochs = 8000: Number of epochs used for the training.
• max_function::Function = (x) -> 30 * 2^(x / 500): Function that changes the magnitude of the change vector of the parameters wrt the epoch.
• train_ix::Union{Nothing, Vector{<:Integer}} = nothing: Indices of the generators used for training if they are not supposed to be picked randomly.
• test_ix::Union{Nothing, Vector{<:Integer}} = nothing: Indices of the generators used for testing if they are not supposed to be picked randomly.

loss(
    sol_cont::SciMLBase.ODESolution,
    sol_disc::SciMLBase.ODESolution,
    ω_proj::SparseArrays.SparseMatrixCSC,
    idxs::Vector{<:Integer}
) -> Vector{Float64}

Calculate the value of the loss function.

projectors_dynamic(
    cm::ContGridModML.ContModel,
    dm::ContGridModML.DiscModel,
    ω_idxs::Vector{<:Integer}
) -> SparseArrays.SparseMatrixCSC{Float64, Int64}

Create projectors for nodal values onto quadrature points and onto comparison locations.

Returned projectors

• ω_proj Project the nodal values onto the given comparison points. This is used to calculate the loss function.

simul(
    disc_sol::SciMLBase.ODESolution,
    M_const::SparseArrays.SparseMatrixCSC,
    K_const::SparseArrays.SparseMatrixCSC,
    m::Array{T<:Real, 1},
    d::Array{T<:Real, 1},
    f::Array{T<:Real, 1},
    A::SparseArrays.SparseMatrixCSC,
    q_proj::SparseArrays.SparseMatrixCSC,
    ω_proj::SparseArrays.SparseMatrixCSC,
    g_proj::Vector{SparseArrays.SparseMatrixCSC},
    idxs::Vector{<:Integer},
    u₀::Array{T<:Real, 1},
    tf::Real;
    cont_kwargs,
    lambda_kwargs
) -> Tuple{Vector{T} where T<:Real, Real}

Do a full simulation step for one data set.

The continuous solution is calculated first and then used to obtain the adjoint solution. Afterwards, the gradient and value of the loss function are calculated and returned.

test_loss(
    disc_sols::Vector{SciMLBase.ODESolution},
    M_const::SparseArrays.SparseMatrixCSC,
    K_const::SparseArrays.SparseMatrixCSC,
    m::Vector{<:Real},
    d::Vector{<:Real},
    f_test::Matrix{<:Real},
    A::SparseArrays.SparseMatrixCSC,
    q_proj::SparseArrays.SparseMatrixCSC,
    ω_proj::SparseArrays.SparseMatrixCSC,
    idxs::Vector{<:Integer},
    u₀::Vector{<:Real},
    tf::Real;
    cont_kwargs
) -> Vector{Float64}

Calculate the loss values for the test data sets.

trapz(
    t₀::Real,
    tf::Real,
    dt::Real,
    int!::Function,
    N::Integer
) -> Vector{Float64}

Trapezoidal rule for the integration of a given function.

update(
    p::Array{T<:Real, 1},
    g::Array{T<:Real, 1},
    eve::Array{T<:Real, 2},
    i::Integer,
    f::Function
) -> Vector{T} where T<:Real

Update the parameters using a restricted gradient descent to ensure positiveness.