digital-brain/content/zettels/system_identification.md

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

Tags

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

System Identification

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 fig:siso_identification_schematic_simplier with one input \(u\) and one output \(y_m\) affected by some disturbances and noise \(d\).

{{< figure src="/ox-hugo/siso_identification_schematic_simplier.png" caption="<span class="figure-number">Figure 1: 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.

Table 1: Different ways to obtain a model

	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:

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

In open-loop identification (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.

{{< figure src="/ox-hugo/siso_identification_schematic_simplier_open_loop.png" caption="<span class="figure-number">Figure 2: 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:

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 fig:system_identification_ol_comp_plant).

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

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.

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

The result for the example is shown in Figure fig:system_identification_ol_coh. At high frequency, the measurement noise dominates and the coherence is poor.

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

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 fig:siso_identification_schematic_simplier_closed_loop).

{{< figure src="/ox-hugo/siso_identification_schematic_simplier_closed_loop.png" caption="<span class="figure-number">Figure 5: 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:

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

This can be generalized to a MIMO plant (Figure fig:siso_identification_schematic_simplier_closed_loop).

{{< figure src="/ox-hugo/mimo_identification_schematic_simplier_closed_loop.png" caption="<span class="figure-number">Figure 6: 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

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 <&pintelon12_system_ident> (chapter 5).

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

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:

%% 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.

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

Choose Sampling Frequency and Duration of Excitation

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}

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}

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}

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}

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.

%% 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 fig:system_identification_multi_sine_asd)

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

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

In the time domain, it is shown in Figure fig:system_identification_multi_sine_time.

{{< figure src="/ox-hugo/system_identification_multi_sine_time.png" caption="<span class="figure-number">Figure 8: 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.

% 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 fig:system_identification_multi_sine_frf. The quality of the obtained FRF is only good in the defined range.

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

generatemultisine - Matlab Function

The synthesis of multi-sine with minimal "crest factor" is taken from <&schroeder70_synth_low_peak_factor_signal>.

The Matlab code is adapted from here.

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

• <pintelon12_system_ident>
• <schoukens12_master>

Bibliography

<./biblio/references.bib>