GitHub

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

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. It is currently not possible to obtain symbolic Jacobians (ModelingToolkit.calculate_jacobian returns symbolic Jacobians for sufficiently simple systems).

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.

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.9999558726115829
-------------------------------------------------
1.0s^2 + 0.9999520094736114s + 0.9999495306487491

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)