diff --git a/content/zettels/system_identification.md b/content/zettels/system_identification.md index cd0584d..b75fa98 100644 --- a/content/zettels/system_identification.md +++ b/content/zettels/system_identification.md @@ -8,39 +8,174 @@ Tags : [Modal Analysis]({{< relref "modal_analysis.md" >}}) -## SISO Identification {#siso-identification} +## System Identification {#system-identification} -### Problem Description {#problem-description} +### Problem Definition {#problem-definition} - - -{{< figure src="/ox-hugo/siso_identification_schematic.png" caption="Figure 1: Block diagram of the SISO system identification" >}} +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 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="Figure 2: Simpler Block diagram of the SISO system identification" >}} +{{< figure src="/ox-hugo/siso_identification_schematic_simplier.png" caption="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\\). -If the open-loop system is unstable, first design a simple controller that stabilizes the system and then identify the closed-loop system. +There are different ways of obtaining a model \\(G\\) of the system summarized in Table . +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 | -### Design of the Excitation Signal {#design-of-the-excitation-signal} +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. -#### Introduction {#introduction} +### Open-loop identification {#open-loop-identification} + +In open-loop identification (Figure ), 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="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: + +```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 ). + + + +{{< figure src="/ox-hugo/system_identification_ol_comp_plant.png" caption="Figure 3: 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 . +At high frequency, the measurement noise dominates and the coherence is poor. + + + +{{< figure src="/ox-hugo/system_identification_ol_coh.png" caption="Figure 4: 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 ). + + + +{{< figure src="/ox-hugo/siso_identification_schematic_simplier_closed_loop.png" caption="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: + +```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 ). + + + +{{< figure src="/ox-hugo/mimo_identification_schematic_simplier_closed_loop.png" caption="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 {#design-of-the-excitation-signal} + + + + +### Introduction {#introduction} There are several choices for excitation signals: - Impulse, Steps - Sweep Sinus -- Random noise, Periodic signals +- Random noise, Periodic signals (PRBS) +- Multi-Sine -#### Random noise with specific ASD {#random-noise-with-specific-asd} +### Random noise with specific ASD {#random-noise-with-specific-asd} The ASD of the measured output is: @@ -82,7 +217,7 @@ u = lsim(G_u, u_norm, t); ``` -#### Choose Sampling Frequency and Duration of Excitation {#choose-sampling-frequency-and-duration-of-excitation} +### Choose Sampling Frequency and Duration of Excitation {#choose-sampling-frequency-and-duration-of-excitation}
@@ -134,39 +269,102 @@ T\_{\text{exc}} = \frac{10}{1} = 10\\,s
-### Computation of the Frequency Response Function {#computation-of-the-frequency-response-function} +### Multi-Sine {#multi-sine} -#### Windowing Function {#windowing-function} - - -#### Example {#example} - -`tfestimate` +### `generatemultisine` - Matlab Function {#generatemultisine-matlab-function} ```matlab -[G, f] = tfestimate(u, y, win, [], [], 1/Ts); -``` +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 -### Verification of the Identification Quality {#verification-of-the-identification-quality} + mag = (args.asd*2*sqrt(Fs)/sqrt(Ns)).^2; + end -`mscohere` + 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 -```matlab -[coh, f] = mscohere(u, y, win, [], [], 1/Ts); + %% 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 {#reference-books} -- (Pintelon and Schoukens 2012) -- (Schoukens, Pintelon, and Rolain 2012) +- +- ## Bibliography {#bibliography} -
-
Pintelon, Rik, and Johan Schoukens. 2012. System Identification : a Frequency Domain Approach. Hoboken, N.J. Piscataway, NJ: Wiley IEEE Press. doi:10.1002/9781118287422.
-
Schoukens, Johan, Rik Pintelon, and Yves Rolain. 2012. Mastering System Identification in 100 Exercises. John Wiley & Sons.
-
+<./biblio/references.bib> diff --git a/static/ox-hugo/identification_closed_loop.png b/static/ox-hugo/identification_closed_loop.png new file mode 100644 index 0000000..016c29f Binary files /dev/null and b/static/ox-hugo/identification_closed_loop.png differ diff --git a/static/ox-hugo/mimo_identification_schematic_simplier_closed_loop.png b/static/ox-hugo/mimo_identification_schematic_simplier_closed_loop.png new file mode 100644 index 0000000..cd77d7b Binary files /dev/null and b/static/ox-hugo/mimo_identification_schematic_simplier_closed_loop.png differ diff --git a/static/ox-hugo/siso_identification_schematic_simplier_closed_loop.png b/static/ox-hugo/siso_identification_schematic_simplier_closed_loop.png new file mode 100644 index 0000000..dffd119 Binary files /dev/null and b/static/ox-hugo/siso_identification_schematic_simplier_closed_loop.png differ diff --git a/static/ox-hugo/siso_identification_schematic_simplier_open_loop.png b/static/ox-hugo/siso_identification_schematic_simplier_open_loop.png new file mode 100644 index 0000000..d53a699 Binary files /dev/null and b/static/ox-hugo/siso_identification_schematic_simplier_open_loop.png differ diff --git a/static/ox-hugo/system_identification_ol_coh.png b/static/ox-hugo/system_identification_ol_coh.png new file mode 100644 index 0000000..42311fd Binary files /dev/null and b/static/ox-hugo/system_identification_ol_coh.png differ diff --git a/static/ox-hugo/system_identification_ol_comp_plant.png b/static/ox-hugo/system_identification_ol_comp_plant.png new file mode 100644 index 0000000..b77c313 Binary files /dev/null and b/static/ox-hugo/system_identification_ol_comp_plant.png differ