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

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

generatemultisine - Matlab Function

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
    sum(rel_mag)
    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>