Implementation of the System Agnostic Localization of Oscillations (SALO) algorithm, presented in the article
The SALO algorithm identifies the source and frequency of a forced oscillation in a complex network of coupled dynamical agents, based on the position and velocity time series of each agent.
The code has been developped on Julia 1.6 and should work on any subsequent version, with up-to-date packages. The packages requires are:
- DelimitedFiles
- Distributed
- FFTW
- Ipopt
- JuMP
- LinearAlgebra
- PyPlot
- example_data: Contains the data necessary to run the two examples.
- plot_SALO.jl: Loads the scripts plotting the outcome of the SALO algorithm.
- run_examples.jl: Loads the functions running the examples for the SALO algorithms.
- SALO.jl: Loads the SALO and SALO-relaxed algorithms.
- SALO_parallel.jl Loads the parallelized versions of the SALO and SALO-relaxed algorithms.
- Lmax_SALO (Lmax_SALO_par)
- Lmax_SALOrelax (Lmax_SALOrelax_par)
- plot_SALO
- plot_SALOrelax_1
- plot_SALOrelax_2
- run_example_ntw20
- run_example_ntw3
- run_SALO (run_SALO_par)
- run_SALOrelax (run_SALOrelax_par)
./SALO.jl (./SALO_parallel.jl)
Lmax_SALO(x::Matrix{Float64}, Dx::Matrix{Float64}, xt::Array{Complex{Float64},2}, Dxt::Array{Complex{Float64,2}, l::Int64, k::Int64, A1h::Matrix{Float64}, a1h::Vector{Float64}, b::Float64=0., μ::Float64=1e-1, bp::Float64=1e-1)Lmax_SALO_par(id::String, x::Matrix{Float64}, Dx::Matrix{Float64}, xt::Array{Complex{Float64},2}, Dxt::Array{Complex{Float64,2}, l::Int64, k::Int64, A1h::Matrix{Float64}, a1h::Vector{Float64}, b::Float64=0., μ::Float64=1e-1, bp::Float64=1e-1)
Minimizes the quadratic error in the estimation of the forced trajectory, for a fixed frequency (k) and location (l) of the forcing. The optimization parameters are the dynamics matrix (A1), the damings (a2), and the forcing amplitude (γ).
INPUT:
(id: For the parallel version only. Name of the system under investigation, for data labelling purpose.)
x: Time series of the phase angles.
Dx: Time series of the phase frequencies.
xt: (Inverse) Fourier transform of x.
Dxt: (Inverse) Fourier transform of Dx.
l: Fixed location of the forcing.
k: Fixed index of the forcing frequency (f = k/T).
A1h: Warm start for A1.
a2h: Warm start for a2.
b: Regularization parameter to avoid overfitting.
μ: Initial value of the barrier parameter (in Ipopt).
bp: Initial value of the bound_push parameter (in Ipopt).
OUTPUT:
objective: Value of the optimized objective function.
A1: Best estimate of the dynamcis matrix.
a2: Best estimate of the dampings.
γ: Best estimate of the forcing amplitude.
./SALO.jl (./SALO_parallel.jl)
Lmax_SALOrelax(x::Matrix{Float64}, Dx::Matrix{Float64}, xt::Array{Complex{Float64},2}, Dxt::Array{Complex{Float64,2}, k::Int64, A1h::Matrix{Float64}, a2h::Vector{Float64}, b::Tuple{Float64,Float64}=(0.,0.), μ::Float64=1e-1, bp::Float64=1e-1)Lmax_SALOrelax_par(id::String, x::Matrix{Float64}, Dx::Matrix{Float64}, xt::Array{Complex{Float64},2}, Dxt::Array{Complex{Float64,2}, k::Int64, A1h::Matrix{Float64}, a2h::Vector{Float64}, b::Tuple{Float64,Float64}=(0.,0.), μ::Float64=1e-1, bp::Float64=1e-1)
Minimizes the quadratic error in the estimation of the forced trajectory, for a fixed frequency (k) of the forcing. The optimization parameters are the dynamics matrix (A1), the damings (a2), and the forcing amplitude (γ).
INPUT:
(id: For the parallel version only. Name of the system under investigation, for data labelling purpose.)
x: Time series of the phase angles.
Dx: Time series of the phase frequencies.
xt: (Inverse) Fourier transform of x.
Dxt: (Inverse) Fourier transform of Dx.
k: Fixed index of the forcing frequency (f = k/T).
A1h: Warm start for A1.
a2h: Warm starrt for a2.
b: Regularizaion parameter to avoid overfitting.
μ: Initial value of the barrier parameter (in Ipopt).
bp: Initial value of the bound_push parameter (in Ipopt).
OUTPUT:
objective: Value of the optimized objective function.
A1: Best estimate of the dynamcis matrix.
a2: Best estimate of the dampings.
γ: Best estimate of the forcing amplitude at each possible location.
./plot_SALO.jl
plot_SALO(Ls::Matrix{Float64}, L::Float64, l::Int64, k::Int64, ks::Vector{Int64}, T::Union{Float64,Int64}, n::Int64)
Plots the value of the objective vs. the estimated (discretized) forcing frequencies by the SALO approach.
INPUT:
Ls: Values of the objective for the various values of l (rows) and k (columns).
L: Minimal objective value.
l0: Corresponding value of l.
k0: Corresponding value of k.
ks: Values of k assessed.
T: Observation time.
n: Number of agents.
./plot_SALO.jl
plot_SALOrelax_1(Ls::Array{Float64,1}, L::Float64, γ::Array{Float64,1}, l0::Int64, k0::Int64, ks::Array{Int64,1}, T::Union{Float64,Int64})
Plots the value of the objective vs. the estimated (discretized) forcing frequencies by the SALO-relax approach.
INPUT:
Ls: Values of the objective for the various values of k.
L: Minimal objective value.
k0: Corresponding value of k.
ks: Values of assessed.
./plot_SALO.jl
plot_SALOrelax_2(Ls::Vector{Float64}, L::Float64, γ::Array{Float64,1}, l0::Int64)
Plots the value of the estimated forcing amplitude vs. the agent index, obtained by the SALO-relax approach.
INPUT:
Ls: Values of the objective for the various values of k.
L: Minimal objective value.
γ: Corresponding forcing amplitudes.
l0: Index with largest amplitude.
./run_examples.jl
run_example_ntw20()
Runs the SALO and SALO-relax algorithms on a 20-node system. All data supporting the examples can be found in the folder example_data.
./run_examples.jl
run_example_ntw3()
Runs the SALO and SALO-relax algorithms on a 3-node system. All data supporting the examples can be found in the folder example_data.
./SALO.jl (./SALO_parallel.jl)
run_SALO(Xs::Matrix{Float64}, τ::Float64, ls::Vector{Int64}, ks::Vector{Int64}, plot::Bool=false, b::Float64=0., μ::Float64=1e-1, bp::Float64=1e-1)run_SALO_par(id::String, Xs::Matrix{Float64}, τ::Float64, ls::Vector{Int64}, ks::Vector{Int64}, b::Float64=0., μ::Float64=1e-1, bp::Float64=1e-1)
Runs the SALO algorithm on the time series Xs, with candidate sources' indices in ls and forcing's candidate frequencies indices in ks. Returns the Maximum Likelihood for each possible pair (l,k) as well as the identified system parameters corresponding to the largest likelihood.
INPUT:
(id: For the parallel version only. Name of the system under investigation, for data labelling purpose.)
Xs: Time series of the phase angles (rows 1:n) and of the phase frequencies (rows n+1:2*n).
τ: Time step size.
ls: Values of l to be tried.
ks: Values of k to be tried.
plot: For the non-parallel version only. If true, generates the plots of the objective function vs. frequency.
b: Regularization parameter to avoid overfitting.
μ: Initial value of the barrier parameter (in Ipopt).
bp: Initial value of the bound_push parameter (in Ipopt).
OUTPUT (written in folder ./data/ for the parallel version):
Ls_l0: values of the objective function for the various values of k and l, for the SALO algorithm.
(L, A, d, γ, k, l): estimated system's paramters from SALO.
L: Minimal value of the objective.A: Estimate of the dynamics matrix.d: Estimate of the dampings.γ: Estimate of the forcing amplitude.k: Estimate frequency index (see theory).l: Estimate of the forcing location.
./SALO.jl (./SALO_parallel.jl)
run_SALOrelax(Xs::Matrix{Float64}, τ::Float64, ks::Vector{Int64}, plot::Bool=false, b::Tuple{Float64,Float64}=(0.,0.), μ::Float64=1e-1, bp::Float64=1e-1)run_SALOrelax_par(id::String, Xs::Matrix{Float64}, τ::Float64, ks::Vector{Int64}, b::Tuple{Float64,Float64}=(0.,0.), μ::Float64=1e-1, bp::Float64=1e-1)
Runs the SALO-relax algorithm on the time series Xs, with forcing's candidate frequencies indices in ks. Returns the Maximum Likelihood for each possible k as well as the identified system parameters corresponding to the largest likelihood.
INPUT:
(id: For parallel version only. Name of the system under investigation, for data labelling purpose.)
Xs: Time series of the phase angles (rows 1:n) and of the phase frequencies (rows n+1:2*n).
τ: Time step size.
ks: Values of k to try.
plot: For the non-parallel version only. If true, generates the plots of the objective function vs. k.
b: Regularization parameters to avoid overfitting.
μ: Initial value of the barrier parameter (in Ipopt).
bp: Initial value of the bound_push parameter (in Ipopt).
OUTPUT (written in folder ./data/ for the parallel version):
Ls_l1: values of the objective function for the various values of k, for the SALO-relax algorithm.
(L, A, d, γ, k, l): estimated system's paramters from SALO-relaxed.
L: Minimal value of the objective.A: Estimate of the dynamics matrix.d: Estimate of the dampings.γ: Estimate of the forcing amplitude.l: Estimate of the forcing location.k: Estimate frequency index (see theory).
This code is provided under a BSD license as part of the Multi-Infrastructure Control and Optimization Toolkit (MICOT) project, LA-CC-13-108.