Linear analysis

Frequency-response analysis

No all dynamical systems are amenable to analytic linearization, and some that technically are, are better to linearize with finite methods. A commonly used technique is frequency-response analysis (FRA), where the systems response to various frequencies are directly measured. The input signal used during FRA can technically be anything, but it is important that the signal is sufficiently exciting in the frequency range of interest. We provide two input methods for FRA, one based on a logarithmic chirp signal, and one based on a set of sinusoidal inputs. We will illustrate their usages by means of two examples.

Sinusoidal input

To perform FRA using sinusoidal inputs, send a frequency vector as second argument to frequency_response_analysis. The system will be simulated once for each frequency.

using JuliaSimControls, Plots
import ModelingToolkitStandardLibrary.Blocks as lib
P0 = ss(tf(1, [1,1,1])) # An example system of second order
P  = lib.StateSpace(ssdata(P0)..., name=:P) # Create a ModelingToolkit StateSpace system
Ts = 0.001 # Sample rate

w = exp10.(LinRange(-1.2, 1, 12)) # Frequency vector
G = frequency_response_analysis(P, w, P.input.u[1], P.output.u[1]; Ts, amplitude=1, settling_time=20, limit_der=true, threads=true)

bodeplot(P0, w, lab="True response", l=(4, 0.5))
scatter!(G, sp=1, lab="Estimated response")
scatter!(w, rad2deg.(angle.(G.r)), sp=2, lab="Estimated response")

The result G is a frequency-response data object which represents the response of the system at each frequency in w. At this stage, you may obtain a rational transfer function like so:

using ControlSystemIdentification
Ghat = tfest(G, tf(2.0, [1,2,1])) # Provide initial guess that determines the order of the transfer function
bodeplot!(Ghat, G.w, lab="Rational estimate", l=:dash)

Chirp input

If you do not pass a frequency vector, but instead provide the keyword arguments f0 and f1, denoting the initial and final frequencies, a chirp input will be used. This method returns not only the estimated frequency response, but also an estimate of disturbances, $H(i\omega)$ (such as nonlinearities, stochasitcity or transient effects)

P0  = ss(tf(1, [1,1,1]))
P   = lib.StateSpace(ssdata(P0)..., name=:P)
res = frequency_response_analysis(P, P.input.u[1], P.output.u[1]; Ts, amplitude=1, f0=w[1]/2pi, f1 = w[end]/2pi, settling_periods=2)
G = res.G
bodeplot(P0, G.w,  lab="True system", l=(6, 0.5))
plot!(G,     sp=1, lab="Estimate", l=(3, 0.8))
plot!(res.H, sp=1, lab="Disturbance estimate", ylims=(1e-2, Inf))
# And the rational estimate
Ghat = tfest(G, tf(2.0, [1,2,1])) # Provide initial guess that determines the order of the transfer function
bodeplot!(Ghat, G.w, lab="Rational estimate", l=:dash, c=4)

You may also estimate a statespace system directly from the data using a subspace-based algorithm. This may be useful for MIMO systems. (Only single input, multiple output is supported by frequency_response_analysis at the moment, MIMO systems thus require several calls to the function.)

Gd = c2d(G, Ts) # Perform a bilinear transformation to discrete time frequency vector
Ph, _ = subspaceid(Gd, Ts, 2, zeroD=true) # (G, Ts, nx, ...)
bodeplot!(Ph, G.w, lab="Statespace estimate", l=:dashdot, c=5)

Note, the different lines are hard to distinguish in the plot since they appear on top of each other.

Handling multiple inputs

The analysis methods described above currently only support single-input systems. For multiple inputs, the analysis may be repeated several times. Parametric estimates, such as transfer function and statespace systems, may be appended in the horizontal direction in order to add inputs, e.g. (pseudocode),

G1f = frequency_response_analysis(..., input1)
G2f = frequency_response_analysis(..., input2)
G1 = subspaceid(G1f)
G2 = subspaceid(G2f)
G = [G1 G2] 
G,_ = balreal(G) # use baltrunc to simplify the model further

The last call to balreal is optional, it converts the model to a balanced realization, possibly removing unobservable and uncontrollable states in the process.

Estimating linearity

Before performing linear analysis, it may be useful to determine if a linear model is an accurate description of the dynamics of the system at the operating point and with the input considered. To this end, we may call any of the functions coherence or coherenceplot (we use the former here to be able to crop the data in the frequency domain before plotting).

ch = coherence(res.d)[0rad:w[end]*rad]
plot(ch, title="Magnitude-squared coherence", yscale=:identity, ylims=(0, 1.01))

A coherence close to 1 is a sign that the data describes a linear input-output system. Note, some decrease in coherence towards the edges of the frequency span is expected, it can be reduced by extending the frequency span for the analysis and cropping the resulting estimate. The FRD objects can be indexed by frequency like

G[1Hz : 3Hz]   # Hz
G[1rad : 3rad] # rad/s

Index

    Docstrings

    JuliaSimControls.frequency_response_analysisMethod
    G(iω) = frequency_response_analysis(G::ODESystem, Ω::AbstractVector, u, y; kwargs...)

    Frequency-response analysis of G u->y. Returns the frequency-response $G(iω)$ as a FRD object that contains the comple response G.r and the frequency vector G.w.

    Note: this is a time-consuming process.

    Arguments

    • Ω: A vector of frequencies
    • u,y: input and output of G to perform the analysis between. Currently, only SIMO analysis supported.
    • Ts: Sample rate
    • settling_time = 2: In seconds, rounded up to closest integer periods
    • nbr_of_periods = 5: to collect data from.
    • amplitude = 1: Scalar or vector of same length as Ω. Very low freqs might require smaller amplitude
    • threads = false: use threads to parallelize the analysis.
    • solver = Rodas4()
    • kwargs... are passed to solve.

    Example:

    using ModelingToolkitStandardLibrary.Blocks, Plots
    P0 = ss(tf(1, [1,1,1]))
    P = Blocks.StateSpace(ssdata(P0)..., name=:P)
    
    w = exp10.(LinRange(-1.2, 1, 12)) # Frequency vector
    G = frequency_response_analysis(P, w, P.u[1], P.y[1], Ts=0.001, amplitude=1, settling_time=20, limit_der=true, threads=true)
    
    bodeplot(P0, w, lab="True response")
    plot!(G, sp=1, lab="Estimated response")
    plot!(w, rad2deg.(angle.(G.r)), sp=2)
    
    # At this stage, you may obtain a rational transfer function like so:
    using ControlSystemIdentification
    Ghat = tfest(G, tf(2.0,[1,2,1])) # Provide initial guess that determines the order of the transfer function
    bodeplot!(Ghat, w, lab="Rational estimate")
    source
    JuliaSimControls.frequency_response_analysisMethod
    G, H, sol, d = frequency_response_analysis(sys::ODESystem, u, y; f0, f1, Tf = 5/f0, Ts = 0.1/f1, amplitude = 1)

    Linearize sys through simulation. Internally, the system is simulated with an exponential chirp input that sweeps the frequencies from f0 to f1.

    The returned system G is of type FRD and contains the comple response G.r and the frequency vector G.w. H::FRD is an estimate of the error as a function of frequency. If the error is large, try increasing Tf or change the amplitude.

    The default chirp duration Tf is 5 periods of the lowest frequency.

    The system identification is performed by sampling the output of the system with a frequency 10f1

    Arguments:

    • sys: System to be linearized
    • u,y: input and output of G to perform the analysis between. Currently, only SIMO analysis supported.
    • f0: Start frequency
    • f1: End frequency
    • Tf: Duration of the chirp experiment
    • Ts: The sample time of the identified model.
    • amplitude: The amplitude of the chirp signal. May be a number of a symbolic expression of t = ModelingToolkit.get_iv(sys).
    • settling_periods = 2: number of periods of the lowest frequency to run before the chirp starts.
    • solver = Rodas4()
    • kwargs... are passed to solve.

    Example

    using ModelingToolkitStandardLibrary.Blocks, Plots
    P0 = ss(tf(1, [1,1,1]))
    P = Blocks.StateSpace(ssdata(P0)..., name=:P)
    res = frequency_response_analysis(P, P.u[1], P.y[1], Ts = 0.001, amplitude=1, f0=w[1]/2pi, f1 = w[end]/2pi, settling_periods=2)
    G = res.G
    bodeplot(P0, G.w)
    plot!(G, sp=1)
    
    # At this stage, you may obtain a rational transfer function like so:
    using ControlSystemIdentification
    Ghat = tfest(G, tf(2.0,[1,2,1])) # Provide initial guess that determines the order of the transfer function
    bodeplot!(Ghat, G.w, lab="Rational estimate")
    source