+++
title = "System Identification"
author = ["Dehaeze Thomas"]
draft = false
+++

Tags
: [Modal Analysis]({{< relref "modal_analysis.md" >}})


## System Identification {#system-identification}


### Problem Definition {#problem-definition}

The goal of a system identification is to extract a model (usually a LTI transfer function) from experimental data.
The system is represented in [Figure 1](#figure--fig:siso-identification-schematic-simplier) with one input \\(u\\) and one output \\(y\_m\\) affected by some disturbances and noise \\(d\\).

<a id="figure--fig:siso-identification-schematic-simplier"></a>

{{< figure src="/ox-hugo/siso_identification_schematic_simplier.png" caption="<span class=\"figure-number\">Figure 1: </span>Simpler Block diagram of the SISO system identification" >}}

The goal of system identification is to compute the transfer function \\(G\\) from known excitation signal \\(u\\) and from a measure of \\(y\_m\\).

There are different ways of obtaining a model \\(G\\) of the system summarized in Table <tab:system_identification_model_types>.
For mechatronics systems, _frequency domain identification_ is the norm.

<a id="table--tab:system-identification-model-types"></a>
<div class="table-caption">
  <span class="table-number"><a href="#table--tab:system-identification-model-types">Table 1</a>:</span>
  Different ways to obtain a model
</div>

|                                     | **Advantages**                                          | **Disadvantages**                               |
|-------------------------------------|---------------------------------------------------------|-------------------------------------------------|
| Physical Modeling                   | Understanding dynamics, no measurements                 | Low accuracy, time consuming                    |
| Time domain identification          | Parametric models, small data sets, non-linearity, MIMO | order/structure selection can be time consuming |
| **Frequency domain identification** | easy, intuitive, user-defined model fitting             | fitting MIMO parametric model can be difficult  |

As an example, a second order plant will be used:

```matlab
G = 1/(1 + 0.1*s/(2*pi*40) + s^2/(2*pi*40)^2); % Plant model
```

The sampling time of the recorded digital signal is 1ms.


### Open-loop identification {#open-loop-identification}

In open-loop identification ([Figure 2](#figure--fig:siso-identification-schematic-simplier-open-loop)), a test signal \\(u\\) is used to _excite_ the system in the frequency range of interest.
The signal \\(u\\) can typically be a swept sine, noise or multi-sine.

<a id="figure--fig:siso-identification-schematic-simplier-open-loop"></a>

{{< figure src="/ox-hugo/siso_identification_schematic_simplier_open_loop.png" caption="<span class=\"figure-number\">Figure 2: </span>Open Loop signal identification" >}}

The plant estimate \\(\hat{G}(j\omega)\\) is computed as follows:
\\[ \hat{G}(j\omega) = \frac{\Phi\_{yu}(\omega)}{\Phi\_{u}(\omega)} \\]
with \\(\phi\_{u}(\omega)\\) is the spectrum of \\(u\\) and \\(\Phi\_{yu}(\omega)\\) is the cross-spectrum of \\(u\\) and \\(y\\).

In Matlab, the transfer function can be extracted from the data using:

```matlab
Nfft = floor(1/Ts); % Size of the window [-]
win = hanning(Nfft); % Hanning window
Noverlap = floor(Nfft/2); % Overlap of 50%

% Identified Transfer Functions
[Gm, f] = tfestimate(data.du, data.y, win, Noverlap, Nfft, 1/Ts);
```

Then, the bode plot of the obtained transfer function is compared against the plant model including a 1.5 sample time delay ([Figure 3](#figure--fig:system-identification-ol-comp-plant)).

<a id="figure--fig:system-identification-ol-comp-plant"></a>

{{< figure src="/ox-hugo/system_identification_ol_comp_plant.png" caption="<span class=\"figure-number\">Figure 3: </span>Comparison of the plant transfer function with the identified FRF in Open-Loop" >}}


### Verification of the Identification Quality - Coherence {#verification-of-the-identification-quality-coherence}

In order to assess the quality of the obtained FRF, the _coherence_ can be computed using the `mscohere` Matlab function.

```matlab
[coh, f] = mscohere(data.du, data.y, win, Noverlap, Nfft, 1/Ts);
```

The result for the example is shown in [Figure 4](#figure--fig:system-identification-ol-coh).
At high frequency, the measurement noise dominates and the coherence is poor.

<a id="figure--fig:system-identification-ol-coh"></a>

{{< figure src="/ox-hugo/system_identification_ol_coh.png" caption="<span class=\"figure-number\">Figure 4: </span>Comparison of the plant transfer function with the identified FRF in Open-Loop" >}}


### Closed-Loop identification {#closed-loop-identification}

If the open-loop system is unstable, a first simple controller needs to be designed to stabilizes the system.
Then, the plant can be identified from closed-loop system identification ([Figure 5](#figure--fig:siso-identification-schematic-simplier-closed-loop)).

<a id="figure--fig:siso-identification-schematic-simplier-closed-loop"></a>

{{< figure src="/ox-hugo/siso_identification_schematic_simplier_closed_loop.png" caption="<span class=\"figure-number\">Figure 5: </span>Closed Loop signal identification" >}}

To perform plant identification in closed-loop:

-   Chose \\(d\_u\\) (sweep sine, multi-sine, noise, ...)
-   Measure \\(u\\) and \\(y\\)
-   Compute the process sensitivity:
    \\[ GS(j\omega) = \frac{\Phi\_{yd\_u}(\omega)}{\Phi\_{d\_u}(\omega)} \\]
-   Compute the Sensitivity:
    \\[ S(j\omega) = \frac{\Phi\_{ud\_u}(\omega)}{\Phi\_{d\_u}(\omega)} \\]
-   Then the plan estimate \\(\hat{G}(j\omega)\\) can be computed:
    \\[ \hat{G}(j\omega) = \frac{GS(j\omega)}{S(j\omega)} \\]

In Matlab, this can be done with the following code:

```matlab
Nfft = floor(1/Ts); % Size of the window [-]
win = hanning(Nfft); % Hanning window
Noverlap = floor(Nfft/2); % Overlap of 50%

[S,  f] = tfestimate(data.du, data.u, win, Noverlap, Nfft, 1/Ts);
[GS, ~] = tfestimate(data.du, data.y, win, Noverlap, Nfft, 1/Ts);

Gm = GS ./ S;
K = (1 ./ S - 1) ./ Gm;
T = 1 - S;
```


### Multi-Input Multi-Output Plant {#multi-input-multi-output-plant}

This can be generalized to a MIMO plant ([Figure 5](#figure--fig:siso-identification-schematic-simplier-closed-loop)).

<a id="figure--fig:siso-identification-schematic-simplier-closed-loop"></a>

{{< figure src="/ox-hugo/mimo_identification_schematic_simplier_closed_loop.png" caption="<span class=\"figure-number\">Figure 6: </span>Closed Loop signal identification (MIMO case)" >}}

Suppose a plan with \\(m\\) inputs and \\(n\\) outputs.
\\(\bm{G}\\) is therefore a \\(n \times m\\) plant, and the controller \\(\bm{K}\\) an \\(m \times n\\) system.
Input sensitivity is an \\(m \times m\\) system:
\\[ \bm{S}\_i = (\bm{I}\_{m} + \bm{K}\bm{G})^{-1} \\]
And process sensitivity is an \\(m \times m\\) system:
\\[ \bm{GS} = \bm{G} (I\_m + KG)^{-1} \\]

And the \\(n \times m\\) plant can be computed using:
\\[ GS \cdot S^{-1} = \bm{G} (I\_m + KG)^{-1} (\bm{I}\_{m} + \bm{K}\bm{G}) = G \\]

To estimate the full plant, \\(m\\) separate identification needs to be performed (one for each input):

-   Chose the excitation signal for the \\(i^{th}\\) input: \\(d\_{ui}\\) (sweep sine, multi-sine, noise, ...)
-   Measure \\(u\_i\\) and \\(\bm{y}\\)
-   Compute the \\(i^{th}\\) column of the process sensitivity (\\(j\\) from \\(1\\) to \\(n\\)):
    \\[ GS\_{ij}(j\omega) = \frac{\Phi\_{y\_jd\_{ui}}(\omega)}{\Phi\_{d\_{ui}}(\omega)} \\]
-   Compute the \\(i^{th}\\) column of the input sensitivity (\\(j\\) from \\(1\\) to \\(n\\)):
    \\[ S\_{ij}(j\omega) = \frac{\Phi\_{u\_jd\_{ui}}(\omega)}{\Phi\_{d\_{ui}}(\omega)} \\]

When the complete \\(GS\\) and \\(S\\) matrices are identified, the plan estimate \\(\hat{\bm{G}}(j\omega)\\) can be computed:
\\[ \hat{\bm{G}}(j\omega) = \bm{GS}(j\omega) \bm{S\_i}^{-1}(j\omega) \\]


## Design of the Excitation Signal {#design-of-the-excitation-signal}

<span class="org-target" id="org-target--sec-system-identification-excitation-signal"></span>


### Introduction {#introduction}

There are several choices for excitation signals:

-   Impulse, Steps
-   Sweep Sinus
-   Random noise, Periodic signals (PRBS)
-   Multi-Sine

A good review is given in (<a href="#citeproc_bib_item_1">Pintelon and Schoukens 2012</a>) (chapter 5).


### Random noise with specific ASD {#random-noise-with-specific-asd}

The ASD of the measured output is:

\begin{equation}
\Gamma\_{y\_m}(\omega) = \Gamma\_d(\omega) + \Gamma\_u(\omega) \cdot |G(j\omega)|
\end{equation}

And we want the effect of the excitation signal to be much higher than the effect of the exogenous signals (measurement noise, input noise, disturbances).

\begin{equation}
\Gamma\_u(\omega) \gg \Gamma\_d(\omega) \cdot |G(j\omega)|^{-1}
\end{equation}

Note that \\(\Gamma\_d(\omega)\\) can be estimated by measuring the system output in the absence of any excitation signal.
The plant magnitude \\(|G(j\omega)|\\) can be roughly estimated from a first identification with bad coherence.

In order to design a random excitation signal with specific spectral characteristics, first a signal with an ASD equal to one is generated (i.e. white noise with unity ASD):

```matlab
Ts = 1e-4; % Sampling Time [s]
t = 0:Ts:10; % Time Vector [s]

%% Signal with an ASD equal to one
u_norm = sqrt(1/2/Ts)*randn(length(t), 1);
```

Then, a transfer function whose magnitude \\(|G\_u(j\omega)|\\) has the same shape as the wanted excitation ASD \\(\Gamma\_u(\omega)\\) is designed:

```matlab
%% Transfer function representing the wanted ASD
G_u = tf([1], [1/2/pi/100 1]);
```

Finally `lsim` is used to compute the shaped excitation signal.

```matlab
%% Shape the ASD of the excitation signal
u = lsim(G_u, u_norm, t);
```


### Choose Sampling Frequency and Duration of Excitation {#choose-sampling-frequency-and-duration-of-excitation}

<div class="important">

The sampling frequency \\(F\_s\\) will determine the maximum frequency \\(F\_{\text{max}}\\) that can be estimated (see Nyquist theorem):

\begin{equation}
F\_{\text{max}} = \frac{1}{2} F\_s
\end{equation}

</div>

<div class="important">

The duration of excitation \\(T\_{\text{exc}}\\) will determine the minimum frequency \\(F\_{\text{min}}\\) that can be estimated:

\begin{equation}
F\_{\text{min}} = \frac{1}{T\_{\text{exc}}}
\end{equation}

It will also corresponds to the frequency resolution \\(\Delta f\\):

\begin{equation}
\Delta f = \frac{1}{T\_{\text{exc}}}
\end{equation}

</div>

In order to increase the estimation quality, averaging can be use with a longer excitation duration.
A factor 10 is usually good enough, therefore the excitation time can be taken as:

\begin{equation}
T\_{\text{exc}} \approx \frac{10}{F\_{\text{min}}}
\end{equation}

<div class="exampl">

Therefore, if the system has to be identified from 1Hz up to 500Hz, the sampling frequency should be:

\begin{equation}
F\_s = 2 F\_{\text{max}} = 1\\,\text{kHz}
\end{equation}

Then, the excitation duration should be (10 averaging):

\begin{equation}
T\_{\text{exc}} = \frac{10}{1} = 10\\,s
\end{equation}

</div>


### Multi-Sine {#multi-sine}

Multi-sine signal excitation has many advantages as compared to random noise:

-   the signal power at each frequency can be precisely chosen
-   the signal is periodic and therefore necessitate no windowing (therefore increasing the obtained FRF quality)

It can be generated with the following code.

```matlab
%% Generate multinsine signal
Fs = 1e3; % Sampling Frequency [Hz]
Ns = 1e3; % Signal length [-]

f = linspace(0, Fs/2, Ns/2); % Frequency Vector [Hz]

% Define the wanted ASD of the test signal [unit/sqrt(Hz)]
wanted_asd = 3*ones(1,Ns);
f_min = 10; % [Hz]
f_max = 300; % [Hz]
wanted_asd([1:round(Ns*f_min/Fs)]) = 0;
wanted_asd([round(Ns*f_max/Fs+1):end]) = 0;

% Generate the multi-sine signal
u = generate_multisine(Fs, Ns, ...
                       'asd', wanted_asd, ...
                       'type', 'schroeder');
```

The ASD of the generated signal is exactly as expected ([Figure 7](#figure--fig:system-identification-multi-sine-asd))

```matlab
[pxx, f] = pwelch(u, ones(Ns,1), [], Ns, Fs);
```

<a id="figure--fig:system-identification-multi-sine-asd"></a>

{{< figure src="/ox-hugo/system_identification_multi_sine_asd.png" caption="<span class=\"figure-number\">Figure 7: </span>Amplitude Spectral Density of the multi-sine signal" >}}

In the time domain, it is shown in [Figure 8](#figure--fig:system-identification-multi-sine-time).

<a id="figure--fig:system-identification-multi-sine-time"></a>

{{< figure src="/ox-hugo/system_identification_multi_sine_time.png" caption="<span class=\"figure-number\">Figure 8: </span>Generated Multi-Sine signal" >}}

Then, the open-loop identification is performed, and the FRF is computed using the following code (not that no window is being used!).
Only the first period (here of 1s) is discarded to remove transient effects.

```matlab
% Skip the first period (transient)
[Gm, f] = tfestimate(data.du(Ns:end), data.y(Ns:end), ones(Ns,1), [], Ns, Fs);
```

The obtained FRF is shown in [Figure 9](#figure--fig:system-identification-multi-sine-frf).
The quality of the obtained FRF is only good in the defined range.

<a id="figure--fig:system-identification-multi-sine-frf"></a>

{{< figure src="/ox-hugo/system_identification_multi_sine_frf.png" caption="<span class=\"figure-number\">Figure 9: </span>Obtained FRF using the multi-sine excitation signal" >}}


### `generatemultisine` - Matlab Function {#generatemultisine-matlab-function}

The synthesis of multi-sine with minimal "crest factor" is taken from (<a href="#citeproc_bib_item_3">Schroeder 1970</a>).

The Matlab code is adapted from [here](https://bholmesqub.github.io/thesis/chapters/identification-design/multi-sine/).

```matlab
function y = generate_multisine(Fs, Ns, args)
%% Input parsing
    arguments
        Fs % Sampling frequency of the generated test signal
        Ns % Length of the generated test signal
        args.asd   double {mustBeNumeric, mustBeNonnegative} = 0
        args.type  char   {mustBeMember(args.type,{'schroeder', 'normal'})} = 'schroeder'
    end

    %% Find amplitude response
    if args.asd == 0
        % If magnitude response is zero, set default to "unitary ASD"
        mag = (ones(1, Ns)*2*sqrt(Fs)/sqrt(Ns)).^2;
        mag(round(Ns/2+1):end) = 0;
    else
        if length(args.asd) ~= Ns
            error('ASD must be of length Ns');
        end

        mag = (args.asd*2*sqrt(Fs)/sqrt(Ns)).^2;
    end

    if any(mag(round(Ns/2+1):end))
        warning('Non-zero magnitude values present outside of frequency limits.');
        mag(round(Ns/2+1):end) = 0;
    end

    %% Find phase response
    switch args.type
      case 'schroeder'
        phase = schroederPhases(Ns, mag);
      case 'normal'
        phase = randn(1, Ns);
      otherwise
        error('type must be "Schroeder" or "normal"');
    end

    %% Generate multisine signal
    % Frequency domain representation
    Y = sqrt(mag/2).*exp(sqrt(-1).*phase);

    % IFFT to time domain
    y = ifft(forceFFTSymmetry(Y))*(Ns/2);

    %% Find zero crossing closest to zero
    % Heuristic method of finding minimum gradient zero crossing
    ySign = y>0;
    zeroInds = find((ySign(2:end) ~= ySign(1:end-1)));

    % Find the index with the smallest gradient around the zero crossing.
    zeroGrad = zeros(1,length(zeroInds));
    for nn=1:length(zeroInds)
        zeroGrad(nn) = abs(y(zeroInds(nn)) - y(zeroInds(nn)+1));
    end
    [~, minInd] = min(zeroGrad);

    yWrapped = [y(zeroInds(minInd):end) y(1:zeroInds(minInd)-1)];
    y = yWrapped;
end

%% Generate phases as defined in "Synthesis of Low-Peak-Factor Signals and
% Binary Sequences With Low Autocorrelation" by M. R. Schroeder
function phase = schroederPhases(Ns, mag)
    rel_mag = mag./sum(mag); % Normalize magnitude for Schroeder's algorithm
    phase = zeros(1, Ns);
    for nn=2:floor(Ns/2+1)
        ll=1:(nn-1);
        phase(nn) = -2*pi*sum((nn-ll).*rel_mag(ll));
    end
end

%% forceFFTSymmetry  A function to force conjugate symmetry on an FFT such that when an
% IFFT is performed the result is a real signal.
% The function has been written to replace MATLAB's ifft(X,'symmetric'), as this function
% is not compatible with MATLAB Coder.
function Y = forceFFTSymmetry(X)
    Y = X;
    XStartFlipped = fliplr(X(2:floor(end/2)));
    Y(ceil(end/2)+2:end) = real(XStartFlipped) - sqrt(complex(-1))*imag(XStartFlipped);
end
```


## Reference Books {#reference-books}

-   (<a href="#citeproc_bib_item_1">Pintelon and Schoukens 2012</a>)
-   (<a href="#citeproc_bib_item_2">Schoukens, Pintelon, and Rolain 2012</a>)


## Bibliography {#bibliography}

<style>.csl-entry{text-indent: -1.5em; margin-left: 1.5em;}</style><div class="csl-bib-body">
  <div class="csl-entry"><a id="citeproc_bib_item_1"></a>Pintelon, Rik, and Johan Schoukens. 2012. <i>System Identification : a Frequency Domain Approach</i>. Hoboken, N.J. Piscataway, NJ: Wiley IEEE Press. doi:<a href="https://doi.org/10.1002/9781118287422">10.1002/9781118287422</a>.</div>
  <div class="csl-entry"><a id="citeproc_bib_item_2"></a>Schoukens, Johan, Rik Pintelon, and Yves Rolain. 2012. <i>Mastering System Identification in 100 Exercises</i>. John Wiley &#38; Sons.</div>
  <div class="csl-entry"><a id="citeproc_bib_item_3"></a>Schroeder, Manfred. 1970. “Synthesis of Low-Peak-Factor Signals and Binary Sequences with Low Autocorrelation (Corresp.).” <i>IEEE Transactions on Information Theory</i> 16 (1). IEEE: 85–89.</div>
</div>