Linear analysis

Linear analysis refers to a set of techniques that operate on linear dynamical systems, e.g., linear statespace models on the form

\[\begin{aligned} \dot{x}(t) &= Ax(t) + Bu(t) \\ y(t) &= Cx(t) + Du(t) \end{aligned}\]

or transfer functions. Nonlinear dynamical systems are commonly linearized in an operating point of interest in order to obtain a linear model suitable for linear analysis.

Linear analysis is used for a wide range of applications, including

  • Frequency-response analysis or modal analysis
  • Stability analysis

This page details how to perform a number of common linear analysis tasks using JuliaSim Control.


Linearization refers to the process of approximating a nonlinear dynamical system with a linear dynamical system. Linearization is commonly performed in a stationary point, referred to as an operating point, or along a trajectory of the nonlinear system. A nonlinear dynamics model may be implemented directly as a function $\dot x = f(x, u)$, or as a ModelingToolkit model. The following two section detail how to linearize models in these two cases.

Differential equations

Linearization of differential equations encoded as Julia functions like $\dot x = f(x, u, p, t)$ can be performed by simply computing the Jacobians of $f$ with respect to $x$ and $u$, see Linearizing nonlinear dynamics for an example demonstrating how to do this.

ModelingToolkit models

ModelingToolkit models can be linearized using the function

lsys_matrices, ssys = linearize(sys::ODESystem, u::Vector{Num}, y::Vector{Num}; op::Dict)

where u and y denote the inputs and outputs, and op is Dict containing the operating point to linearize around. If op is not specified, or only specifies some of the variables n sys, default values are used for non-specified variables.

  • lsys_matrices is a NamedTuple of statespace matrices A,B,C,D that can be transformed to a ControlSystemsBase.StateSpace object using ss(lsys_matrices...) or a linear ODESystem using ModelingToolkitStandardLibrary.Blocks.StateSpace(lsys_matrices...).
  • ssys is a simplified version of the original sys that indicates the order of the state variables in the linearized statespace representation.

Internally, ForwardDiff.jl is used for linearization. To instead obtain symbolic Jacobians, use ModelingToolkit.linearize_symbolic.

More details on the linearization in ModelingToolkit is available in the documentation for ModelingToolkit.

Sometimes, numerical linearization fails, e.g., if the system to be linearized

  • contains discontinuities, in particular at the linearization point (Coulomb friction is a common example)
  • throws an error when ForwardDiff.jl is used

In these situations, it might be better to linearize the system using simulation-based methods, such as Frequency-response analysis.

A video tutorial on using linearization and analysis points is available below.

Batch linearization and gain scheduling

See the tutorial Batch Linearization and gain scheduling. In this tutorial, use of the functions batch_ss and trajectory_ss for linearization around several operating points, or along a trajectory, is demonstrated.

Analysis points

Analysis points provide an interface to give names to points of interest in a causal ModelingToolkit model, such as the model of a control system. This allows the user to linearize models and derive, e.g., sensitivity functions and loop-transfer functions with a simple interface. See Linear analysis with ModelingToolkitStandardLibrary and the video above for more details.

Frequency-response analysis

Not 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 JuliaSimControl, 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")
Example block output

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
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0s^2 + 0.999952009473611s + 0.9999495306487488

Continuous-time transfer function model
bodeplot!(Ghat, G.w, lab="Rational estimate", l=:dash)
Example block output

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)
Example block output

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)
Example block output

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))
Example block output

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

Example: Modal analysis of a series of masses and springs

This example will demonstrate how to make use of the function ModelingToolkit.linearize to perform modal analysis of a system with several masses (inertias) connected in series. This kind of system is common in engineering applications, it arises as, e.g.,

  • Finite approximation of continuous beams.
  • In drive trains, where several rotating inertias are connected through flexible shafts.

We start by creating a function mass_spring_damper_chain that allows us to connect several masses, springs and dampers in series.

using ModelingToolkit, ModelingToolkitStandardLibrary.Mechanical.Rotational, LinearAlgebra
import ModelingToolkitStandardLibrary.Blocks
connect = ModelingToolkit.connect

function mass_spring_damper_chain(N = 4; input=false, random=false)
    @named r = Blocks.Step(offset=1, start_time=1, height=-1)
    @named t = Torque(use_support=false)
    @named sens = Rotational.AngleSensor()
    systems = input ? [t, sens, r] : [t, sens]
    local m, s, d, eqs
    for i in 1:N
        if random
            m2 = Inertia(J = 10rand(), name = Symbol("m$i"))
            s2 = Spring(c = 100 * rand(), name = Symbol("s$i"))
            d2 = Damper(d = rand(), name = Symbol("d$i"))
            m2 = Inertia(J = 10, name = Symbol("m$i"))
            s2 = Spring(c = 100, name = Symbol("s$i"))
            d2 = Damper(d = 1e-2, name = Symbol("d$i"))

        if i == 1
            eqs = if input
                [ModelingToolkit.connect(t.flange, m2.flange_a); ModelingToolkit.connect(r.output, :u, t.tau)]
                [ModelingToolkit.connect(t.flange, m2.flange_a)]

        # Connect to following spring-damper
        if i < N
            push!(eqs, ModelingToolkit.connect(m2.flange_b, s2.flange_a, d2.flange_a))

        # Connect mass to previous spring-damper
        if 1 < i
            push!(eqs, ModelingToolkit.connect(m2.flange_a, s.flange_b, d.flange_b))
        s = s2
        d = d2
        m = m2
        push!(systems, m)
        if i < N # Do not add any more spring-dampers at the last point
            push!(systems, d)
            push!(systems, s)
    push!(eqs, connect(m.flange_b, sens.flange))
    @named mass_spring_damper_chain = ODESystem(eqs, ModelingToolkit.get_iv(systems[1]);
mass_spring_damper_chain (generic function with 2 methods)

Next, we instantiate a model with N = 3 masses.

using JuliaSimControl, StaticArrays, Plots
N = 3
model = mass_spring_damper_chain(N; input=false) |> complete

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( t_{+}flange, m1_{+}flange_{a} \right) \\ \mathrm{connect}\left( m1_{+}flange_{b}, s1_{+}flange_{a}, d1_{+}flange_{a} \right) \\ \mathrm{connect}\left( m2_{+}flange_{b}, s2_{+}flange_{a}, d2_{+}flange_{a} \right) \\ \mathrm{connect}\left( m2_{+}flange_{a}, s1_{+}flange_{b}, d1_{+}flange_{b} \right) \\ \mathrm{connect}\left( m3_{+}flange_{a}, s2_{+}flange_{b}, d2_{+}flange_{b} \right) \\ \mathrm{connect}\left( m3_{+}flange_{b}, sens_{+}flange \right) \\ t_{+}phi_{support}\left( t \right) = 0 \\ t_{+}flange_{+}tau\left( t \right) = - t_{+}tau_{+}u\left( t \right) \\ sens_{+}phi_{+}u\left( t \right) = sens_{+}flange_{+}phi\left( t \right) \\ sens_{+}flange_{+}tau\left( t \right) = 0 \\ m1_{+}phi\left( t \right) = m1_{+}flange_{a_{+}phi}\left( t \right) \\ m1_{+}phi\left( t \right) = m1_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} m1_{+}phi\left( t \right)}{\mathrm{d}t} = m1_{+}w\left( t \right) \\ \frac{\mathrm{d} m1_{+}w\left( t \right)}{\mathrm{d}t} = m1_{+}a\left( t \right) \\ m1_{+}J m1_{+}a\left( t \right) = m1_{+}flange_{a_{+}tau}\left( t \right) + m1_{+}flange_{b_{+}tau}\left( t \right) \\ d1_{+}phi_{rel}\left( t \right) = - d1_{+}flange_{a_{+}phi}\left( t \right) + d1_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} d1_{+}phi_{rel}\left( t \right)}{\mathrm{d}t} = d1_{+}w_{rel}\left( t \right) \\ \frac{\mathrm{d} d1_{+}w_{rel}\left( t \right)}{\mathrm{d}t} = d1_{+}a_{rel}\left( t \right) \\ d1_{+}flange_{b_{+}tau}\left( t \right) = d1_{+}tau\left( t \right) \\ d1_{+}flange_{a_{+}tau}\left( t \right) = - d1_{+}tau\left( t \right) \\ d1_{+}tau\left( t \right) = d1_{+}d d1_{+}w_{rel}\left( t \right) \\ s1_{+}phi_{rel}\left( t \right) = - s1_{+}flange_{a_{+}phi}\left( t \right) + s1_{+}flange_{b_{+}phi}\left( t \right) \\ s1_{+}flange_{b_{+}tau}\left( t \right) = s1_{+}tau\left( t \right) \\ s1_{+}flange_{a_{+}tau}\left( t \right) = - s1_{+}tau\left( t \right) \\ s1_{+}tau\left( t \right) = s1_{+}c \left( - s1_{+}phi_{rel0} + s1_{+}phi_{rel}\left( t \right) \right) \\ m2_{+}phi\left( t \right) = m2_{+}flange_{a_{+}phi}\left( t \right) \\ m2_{+}phi\left( t \right) = m2_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} m2_{+}phi\left( t \right)}{\mathrm{d}t} = m2_{+}w\left( t \right) \\ \frac{\mathrm{d} m2_{+}w\left( t \right)}{\mathrm{d}t} = m2_{+}a\left( t \right) \\ m2_{+}J m2_{+}a\left( t \right) = m2_{+}flange_{a_{+}tau}\left( t \right) + m2_{+}flange_{b_{+}tau}\left( t \right) \\ d2_{+}phi_{rel}\left( t \right) = - d2_{+}flange_{a_{+}phi}\left( t \right) + d2_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} d2_{+}phi_{rel}\left( t \right)}{\mathrm{d}t} = d2_{+}w_{rel}\left( t \right) \\ \frac{\mathrm{d} d2_{+}w_{rel}\left( t \right)}{\mathrm{d}t} = d2_{+}a_{rel}\left( t \right) \\ d2_{+}flange_{b_{+}tau}\left( t \right) = d2_{+}tau\left( t \right) \\ d2_{+}flange_{a_{+}tau}\left( t \right) = - d2_{+}tau\left( t \right) \\ d2_{+}tau\left( t \right) = d2_{+}d d2_{+}w_{rel}\left( t \right) \\ s2_{+}phi_{rel}\left( t \right) = - s2_{+}flange_{a_{+}phi}\left( t \right) + s2_{+}flange_{b_{+}phi}\left( t \right) \\ s2_{+}flange_{b_{+}tau}\left( t \right) = s2_{+}tau\left( t \right) \\ s2_{+}flange_{a_{+}tau}\left( t \right) = - s2_{+}tau\left( t \right) \\ s2_{+}tau\left( t \right) = s2_{+}c \left( - s2_{+}phi_{rel0} + s2_{+}phi_{rel}\left( t \right) \right) \\ m3_{+}phi\left( t \right) = m3_{+}flange_{a_{+}phi}\left( t \right) \\ m3_{+}phi\left( t \right) = m3_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} m3_{+}phi\left( t \right)}{\mathrm{d}t} = m3_{+}w\left( t \right) \\ \frac{\mathrm{d} m3_{+}w\left( t \right)}{\mathrm{d}t} = m3_{+}a\left( t \right) \\ m3_{+}J m3_{+}a\left( t \right) = m3_{+}flange_{a_{+}tau}\left( t \right) + m3_{+}flange_{b_{+}tau}\left( t \right) \\ \end{array} \right] \end{equation} \]

We call the function complete on the model in order to be able to refer to state names using the syntax model.m1.w below. With our ODESystem model in hand, we may call named_ss, which internally calls linearize, to get a NamedStateSpace object

lsys = named_ss(model, [model.t.tau.u], [model.m1.w,model.m2.w,model.m3.w])
mass_spring_damper_chain: NamedStateSpace{Continuous, Float64}
A = 
   0.0    0.0    0.0   1.0     0.0     0.0
   0.0    0.0    0.0   0.0     1.0     0.0
   0.0    0.0    0.0   0.0     0.0     1.0
 -10.0    0.0   10.0  -0.001   0.0     0.001
  -0.0  -10.0   10.0  -0.0    -0.001   0.001
  10.0   10.0  -20.0   0.001   0.001  -0.002
B = 
C = 
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  1.0  0.0
D = 

Continuous-time state-space model
With state  names: m1₊phi(t) m3₊phi(t) m2₊phi(t) m1₊w(t) m3₊w(t) m2₊w(t)
     input  names: t₊tau₊u(t)
     output names: m1₊w(t) m2₊w(t) m3₊w(t)

The system does not have a minimal realization, and we may remove one state using minreal, let's do that before we continue

lsys = minreal(lsys)

We may plot the Bode curve of the system to visualize its frequency response

w = exp10.(LinRange(-1, 1, 500)) # A frequency vector
bodeplot(lsys, w, xticks=exp10.(LinRange(-1, 1, 5)), plotphase=false, layout=1)
Example block output

We may also hit the system with a hammer and check its time response (impulse response)

plot(impulse(lsys, 10))
Example block output

Before we start with the modal analysis, we check the output of dampreport to get a quick overview of the systems vibration characteristics:

|        Pole        |   Damping     |   Frequency   |   Frequency   | Time Constant |
|                    |    Ratio      |   (rad/sec)   |     (Hz)      |     (sec)     |
| -4.1e-17           |  1            |  4.1e-17      |  6.52e-18     |  2.44e+16     |
| -0.0005 ±   3.16im |  0.000158     |  3.16         |  0.503        |  2e+03        |
| -0.0015 ±   5.48im |  0.000274     |  5.48         |  0.872        |  667          |

Next up we carry on with the modal analysis. We extract the system matrices, $A,B,C,D$ so that we may perform eigen analysis on the $A$ matrix and project the result down to the output space using $C$. The mode shapes are given by the eigen vectors to the $A$ matrix, and the frequencies are given by the eigenvalues. Since we are dealing with a real-valued system, each complex eigenvalue is paired with a complex conjugate, we thus access the frequency of the second mode as freqs[3] below.

A,B,C,D = ssdata(lsys)

lab = string.(permutedims(output_names(lsys)))
fig_modes = plot(real.(C*eigen(A).vectors)[:, 1:2:end]; c=(1:3)', layout=(1,3), plot_title="Mode shapes", ylims=(-1,1), lab=permutedims(["Mode $i" for i in 1:3]), m=:o)
hline!([0 0 0], l=(:black, ), primary=false,

freqs = abs.(eigvals(lsys.A))

u(x,t) = SA[cos(freqs[1]*t)]
res1 = lsim(lsys, u, 0:0.01:20)
fig_mode1 = plot(res1; ylims=(-0.5, 0.5), layout=1, sp=1, lab, ylabel="Mode 1")

u(x,t) = SA[cos(freqs[3]*t)]
res2 = lsim(lsys, u, 0:0.01:20)
fig_mode2 = plot(res2; ylims=(-0.5, 0.5), layout=1, sp=1, lab, ylabel="Mode 2")

plot(fig_modes, fig_mode1, fig_mode2, layout=(3,1), size=(800, 800))
Example block output

The top figure illustrates the mode shapes. In the first mode, the masses 1 and 3 vibrate in phase, and mass 2 is 180° out of phase. This corresponds to the time-domain plot in the middle, where we simulate the system with a sinusoidal input force at the frequency of the first mode. We also see that the amplitude of the second mass is slightly higher than the amplitudes of masses 1 and 3, which is also indicated by the mode shape. The second mode shape indicates that the first and third masses vibrate out of phase, while the second mass is almost still. This is verified in the last plot, where we see a time-domain simulation with the frequency of the second mode driving the system. The third mode is not visualized in simulation, it corresponds to an integrating mode and is not very interesting (all three masses move together as a single mass, a rigid-body mode).

For fun, let's also animate the movements of the inertias to illustrate the mode shapes. Visualizing the movement of rotating inertias is not very easy, we thus treat them as linear displacement masses instead, and displace the coordinates slightly to make the movement easier to track by eye.

function plot_msd(res, title)
    y = copy(res.y)
    y .+= 0:size(y, 1)-1 # Move coordinates to make a nice plot
    @gif for i in 1:10:size(y, 2)
        plot([y[1:2, i] y[2:3, i]], zeros(2,2); m=:square, c=:black, ms=7, legend=false, xlims=(-1, 3), ylims=(-0.1, 0.1), size=(400, 100), xticks=false, yticks=false, title, framestyle=:none)

plot_msd(res1, "Mode 1")
Example block output
plot_msd(res2, "Mode 2")
Example block output

Control design

Say that we are interested in designing a disturbance-rejection regulator with feedback from a sensor measuring the position of the first mass. Such a system is known to be passive, i.e., it is known to have a Nyquist curve that lies entirely in the right half plane. We can verify this with the following three plots

w = exp10.(LinRange(-1, 1, 5000))
sys1 = lsys[Symbol("m1₊w(t)"), Symbol("t₊tau₊u(t)")]
    bodeplot(sys1, w),
    nyquistplot(sys1, w, xlims=(-0.1, 0.1), ylims=(-3, 3)),
    passivityplot(sys1, w),
    layout = (1,3),
    size   = (800, 300),
    legend = false
Example block output

These plots all give us indications that the system is passive, in the Bode plot, we see that the phase curve never goes outside of the band ±90°, and in the Nyquist plot we see that the curve lies in the right half plane. In the last plot, the passivityplot, we see the passivity index as a function of frequency. This index is never above 1, indicating that the system is passive for all frequencies.

The feedback interconnection of two passive systems is known to be passive, we thus have the opportunity to design a passive controller, and will automatically have a robustly stable closed-loop system no matter how large model errors we have! To design a passive controller, we make use of the function spr_synthesize (SPR, Strictly Positive Real).

P0 = ExtendedStateSpace(sys1.sys, D21=[0 1], B1 = [sys1.B zeros(sys1.nx)]) # Add a measurement-noise model
# P0.A .-= 1e-6*I(P0.nx)
K, Gcl, sν = spr_synthesize(P0)

As we can see, the synthesized controller is passive. We can also check the closed-loop system from input disturbances to the output

Gzw = Gcl[1,1]
w = exp10.(LinRange(0, 1, 5000))
    bodeplot([sys1, Gzw], w, lab=["Open loop" "Closed loop"], plotphase=false, ylims=(1e-3, Inf)),
    nyquistplot(Gzw, w, lab="Closed loop"),
    plot_title = "Input disturbance to output",
    size   = (800, 350),
Example block output

Also this system is passive. We have thus designed a robust, passive controller using the knowledge that the system is passive.

For more details on passive controller synthesis, see $\mathcal{H}_2$ Synthesis of a passive controller.



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.


  • Ω: 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
  • offset: A bias/offset around which the signal is oscillating.
  • diff_order = 0: Order of differentiation to apply to the input signal before estimating the frequency response. This is useful to mitigate problems due to non-stationarity, such as trends and drift. The resulting frequency-response is automatically integrated in the frequency domain to account for the applied differentiation. Try setting diff_order = 1 if the result is poor for high frequencies when using offset.
  • threads = false: use threads to parallelize the analysis.
  • solver = Rodas4()
  • kwargs... are passed to solve.


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")
G, H, sol, d = frequency_response_analysis(sys::ODESystem, u, y; f0, f1, Tf = 5/f0, Ts = 0.1/f1, amplitude = 1, offset = 0)

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


  • 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).
  • offset: A bias/offset around which the chirp is oscillating.
  • settling_periods = 2: number of periods of the lowest frequency to run before the chirp starts.
  • solver = Rodas4()
  • kwargs... are passed to solve.


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")