# 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_analysis`

— Method`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")
```

`JuliaSimControls.frequency_response_analysis`

— Method`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")
```