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

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

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}}
               0.9999558722776221
-------------------------------------------------
1.0s^2 + 0.9999520087503453s + 0.9999495302712531

Continuous-time transfer function model

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

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

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

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

        # Connect mass to previous spring-damper
        if 1 < i
            push!(eqs, ModelingToolkit.connect(m2.flange_a, s.flange_b, d.flange_b))
        end
        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)
        end
    end
    push!(eqs, connect(m.flange_b, sens.flange))
    @named mass_spring_damper_chain = ODESystem(eqs, ModelingToolkit.get_iv(systems[1]);
                                                systems)
end

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_{b_{+}tau}\left( t \right) + m1_{+}flange_{a_{+}tau}\left( t \right) \\ d1_{+}phi_{rel}\left( t \right) = d1_{+}flange_{b_{+}phi}\left( t \right) - d1_{+}flange_{a_{+}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_{b_{+}phi}\left( t \right) - s1_{+}flange_{a_{+}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_{b_{+}phi}\left( t \right) - s2_{+}flange_{a_{+}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_{b_{+}tau}\left( t \right) + m3_{+}flange_{a_{+}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 = 
 0.0
 0.0
 0.0
 0.1
 0.0
 0.0
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 = 
 0.0
 0.0
 0.0

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)

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

plot(impulse(lsys, 10))

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

dampreport(lsys)

|        Pole        |   Damping     |   Frequency   |   Frequency   | Time Constant |
|                    |    Ratio      |   (rad/sec)   |     (Hz)      |     (sec)     |
+--------------------+---------------+---------------+---------------+---------------+
| -3.48e-16          |  1            |  3.48e-16     |  5.54e-17     |  2.87e+15     |
| -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, topmargin=-10Plots.mm)

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

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

plot_msd(res1, "Mode 1")

plot_msd(res2, "Mode 2")

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)")]
plot(
    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
)

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)
ispassive(-K)

true

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))
plot(
    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),
)

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.

Index

• JuliaSimControl.frequency_response_analysis
• JuliaSimControl.frequency_response_analysis

Docstrings

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

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)

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