From f5bccf5f61687c8e491ef1ec466e3075d1791eed Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Tue, 10 Sep 2019 18:18:28 +0200 Subject: [PATCH] Rename Data File, Export matlab scripts and data --- index.html | 78 ++++--- index.org | 26 ++- mat/{data_001.mat => data_ux.mat} | Bin mat/{data_002.mat => data_uy.mat} | Bin matlab/plant_identification.m | 341 ++++++++++++++++++++++++++++++ 5 files changed, 409 insertions(+), 36 deletions(-) rename mat/{data_001.mat => data_ux.mat} (100%) rename mat/{data_002.mat => data_uy.mat} (100%) create mode 100644 matlab/plant_identification.m diff --git a/index.html b/index.html index e15c45d..5cf4bd6 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Cercalo Test Bench @@ -276,28 +276,38 @@ for the JavaScript code in this tag.

Table of Contents

-
-

1 Identification

+
+

1 Identification

+

+ +

+
+

+All the files (data and Matlab scripts) are accessible here. +

+
-
-

1.1 Excitation Data

+
+ +
+

1.1 Excitation Data

fs = 1e4;
@@ -334,15 +344,15 @@ Discrete-time zero/pole/gain model.
-
-

1.2 Input / Output data

+
+

1.2 Input / Output data

The identification data is loaded

-
ux = load('mat/data_001.mat', 't', 'ux', 'yx', 'yy');
-uy = load('mat/data_002.mat', 't', 'uy', 'yx', 'yy');
+
ux = load('mat/data_ux.mat', 't', 'ux', 'yx', 'yy');
+uy = load('mat/data_uy.mat', 't', 'uy', 'yx', 'yy');
@@ -378,7 +388,7 @@ uy.yy = uy.yy-mean +

identification_ux.png

Figure 1: Identification signals when exciting the \(x\) axis (png, pdf)

@@ -386,7 +396,7 @@ uy.yy = uy.yy-mean +

identification_uy.png

Figure 2: Identification signals when exciting the \(y\) axis (png, pdf)

@@ -394,8 +404,8 @@ uy.yy = uy.yy-mean -

1.3 Estimation of the Frequency Response Function Matrix

+
+

1.3 Estimation of the Frequency Response Function Matrix

We compute an estimate of the transfer functions. @@ -409,7 +419,7 @@ We compute an estimate of the transfer functions.

-
+

frequency_response_matrix.png

Figure 3: Frequency Response Matrix (png, pdf)

@@ -417,8 +427,8 @@ We compute an estimate of the transfer functions.
-
-

1.4 Coherence

+
+

1.4 Coherence

[coh_ux_yx, f] = mscohere(ux.ux, ux.yx, hanning(ceil(1*fs)), [], [], fs);
@@ -429,7 +439,7 @@ We compute an estimate of the transfer functions.
-
+

identification_coherence.png

Figure 4: Coherence (png, pdf)

@@ -438,8 +448,8 @@ We compute an estimate of the transfer functions.
-
-

1.5 Extraction of a transfer function matrix

+
+

1.5 Extraction of a transfer function matrix

First we define the initial guess for the resonance frequencies and the weights associated. @@ -493,7 +503,7 @@ Ignore data above some frequency.

-
+

weights.png

Figure 5: Weights amplitude (png, pdf)

@@ -545,7 +555,7 @@ An we run the vectfit3 algorithm.
-
+

identification_matrix_fit.png

Figure 6: Transfer Function Extraction of the FRF matrix (png, pdf)

@@ -572,16 +582,16 @@ G = [G_ux_yx, G_uy_yx;
-
-

2 Plant Analysis

+
+

2 Plant Analysis

-
-

3 Control

+
+

3 Control

Author: Dehaeze Thomas

-

Created: 2019-09-10 mar. 18:15

+

Created: 2019-09-10 mar. 18:18

Validate

diff --git a/index.org b/index.org index 101a343..ece7dcd 100644 --- a/index.org +++ b/index.org @@ -26,6 +26,28 @@ * Identification +:PROPERTIES: +:header-args:matlab+: :tangle matlab/plant_identification.m +:header-args:matlab+: :comments org :mkdirp yes +:END: +<> + +** ZIP file containing the data and matlab files :ignore: +#+begin_src bash :exports none :results none + if [ matlab/plant_identification.m -nt data/plant_identification.zip ]; then + cp matlab/plant_identification.m plant_identification.m; + zip data/plant_identification \ + mat/data_ux.mat \ + mat/data_uy.mat \ + plant_identification.m + rm plant_identification.m; + fi +#+end_src + +#+begin_note + All the files (data and Matlab scripts) are accessible [[file:data/plant_identification.zip][here]]. +#+end_note + ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) <> @@ -68,8 +90,8 @@ Discrete-time zero/pole/gain model. ** Input / Output data The identification data is loaded #+begin_src matlab - ux = load('mat/data_001.mat', 't', 'ux', 'yx', 'yy'); - uy = load('mat/data_002.mat', 't', 'uy', 'yx', 'yy'); + ux = load('mat/data_ux.mat', 't', 'ux', 'yx', 'yy'); + uy = load('mat/data_uy.mat', 't', 'uy', 'yx', 'yy'); #+end_src We remove the first seconds where the Cercalo is turned on. diff --git a/mat/data_001.mat b/mat/data_ux.mat similarity index 100% rename from mat/data_001.mat rename to mat/data_ux.mat diff --git a/mat/data_002.mat b/mat/data_uy.mat similarity index 100% rename from mat/data_002.mat rename to mat/data_uy.mat diff --git a/matlab/plant_identification.m b/matlab/plant_identification.m new file mode 100644 index 0000000..743cc52 --- /dev/null +++ b/matlab/plant_identification.m @@ -0,0 +1,341 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +% Excitation Data + +fs = 1e4; +Ts = 1/fs; + + + +% We generate white noise with the "random number" simulink block, and we filter that noise. + + +Gi = (1)/(1+s/2/pi/100); + +c2d(Gi, Ts, 'tustin') + +% Input / Output data +% The identification data is loaded + +ux = load('mat/data_ux.mat', 't', 'ux', 'yx', 'yy'); +uy = load('mat/data_uy.mat', 't', 'uy', 'yx', 'yy'); + + + +% We remove the first seconds where the Cercalo is turned on. + +i0x = 20*fs; + +i0y = 10*fs; + +ux.t = ux.t( i0x:end) - ux.t(i0x); +ux.ux = ux.ux(i0x:end); +ux.yx = ux.yx(i0x:end); +ux.yy = ux.yy(i0x:end); + +uy.t = uy.t( i0y:end) - uy.t(i0x); +uy.uy = uy.uy(i0y:end); +uy.yx = uy.yx(i0y:end); +uy.yy = uy.yy(i0y:end); + +ux.ux = ux.ux-mean(ux.ux); +ux.yx = ux.yx-mean(ux.yx); +ux.yy = ux.yy-mean(ux.yy); + +uy.ux = uy.ux-mean(uy.ux); +uy.yx = uy.yx-mean(uy.yx); +uy.yy = uy.yy-mean(uy.yy); + +figure; +ax1 = subplot(1, 2, 1); +plot(ux.t, ux.ux); +xlabel('Time [s]'); +ylabel('Amplitude [V]'); +legend({'$u_x$'}); + +ax2 = subplot(1, 2, 2); +hold on; +plot(ux.t, ux.yx, 'DisplayName', '$y_x$'); +plot(ux.t, ux.yy, 'DisplayName', '$y_y$'); +hold off; +xlabel('Time [s]'); +ylabel('Amplitude [V]'); +legend() + +linkaxes([ax1,ax2],'x'); +xlim([ux.t(1), ux.t(end)]) + + + +% #+NAME: fig:identification_ux +% #+CAPTION: Identification signals when exciting the $x$ axis ([[./figs/identification_ux.png][png]], [[./figs/identification_ux.pdf][pdf]]) +% [[file:figs/identification_ux.png]] + + + +figure; +ax1 = subplot(1, 2, 1); +plot(uy.t, uy.uy); +xlabel('Time [s]'); +ylabel('Amplitude [V]'); +legend({'$u_y$'}); + +ax2 = subplot(1, 2, 2); +hold on; +plot(uy.t, uy.yy, 'DisplayName', '$y_y$'); +plot(uy.t, uy.yx, 'DisplayName', '$y_x$'); +hold off; +xlabel('Time [s]'); +ylabel('Amplitude [V]'); +legend() + +linkaxes([ax1,ax2],'x'); +xlim([uy.t(1), uy.t(end)]) + +% Estimation of the Frequency Response Function Matrix +% We compute an estimate of the transfer functions. + +[tf_ux_yx, f] = tfestimate(ux.ux, ux.yx, hanning(ceil(1*fs)), [], [], fs); +[tf_ux_yy, ~] = tfestimate(ux.ux, ux.yy, hanning(ceil(1*fs)), [], [], fs); +[tf_uy_yx, ~] = tfestimate(uy.uy, uy.yx, hanning(ceil(1*fs)), [], [], fs); +[tf_uy_yy, ~] = tfestimate(uy.uy, uy.yy, hanning(ceil(1*fs)), [], [], fs); + +figure; +ax11 = subplot(2, 2, 1); +hold on; +plot(f, abs(tf_ux_yx)) +title('Frequency Response Function $\frac{y_x}{u_x}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude') +hold off; + +ax12 = subplot(2, 2, 2); +hold on; +plot(f, abs(tf_ux_yy)) +title('Frequency Response Function $\frac{y_x}{u_y}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +hold off; + +ax21 = subplot(2, 2, 3); +hold on; +plot(f, abs(tf_uy_yx)) +title('Frequency Response Function $\frac{y_y}{u_x}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude') +xlabel('Frequency [Hz]') +hold off; + +ax22 = subplot(2, 2, 4); +hold on; +plot(f, abs(tf_uy_yy)) +title('Frequency Response Function $\frac{y_y}{u_y}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +xlabel('Frequency [Hz]') +hold off; + +linkaxes([ax11,ax12,ax21,ax22],'x'); +xlim([10, 1000]); +linkaxes([ax11,ax12,ax21,ax22],'y'); +ylim([1e-2, 1e3]) + +% Coherence + +[coh_ux_yx, f] = mscohere(ux.ux, ux.yx, hanning(ceil(1*fs)), [], [], fs); +[coh_ux_yy, ~] = mscohere(ux.ux, ux.yy, hanning(ceil(1*fs)), [], [], fs); +[coh_uy_yx, ~] = mscohere(uy.uy, uy.yx, hanning(ceil(1*fs)), [], [], fs); +[coh_uy_yy, ~] = mscohere(uy.uy, uy.yy, hanning(ceil(1*fs)), [], [], fs); + +figure; +ax11 = subplot(2, 2, 1); +hold on; +plot(f, coh_ux_yx) +set(gca, 'Xscale', 'log'); +title('Coherence $\frac{y_x}{u_x}$') +ylabel('Coherence') +hold off; + +ax12 = subplot(2, 2, 2); +hold on; +plot(f, coh_ux_yy) +set(gca, 'Xscale', 'log'); +title('Coherence $\frac{y_x}{u_y}$') +hold off; + +ax21 = subplot(2, 2, 3); +hold on; +plot(f, coh_uy_yx) +set(gca, 'Xscale', 'log'); +title('Coherence $\frac{y_y}{u_x}$') +ylabel('Coherence') +xlabel('Frequency [Hz]') +hold off; + +ax22 = subplot(2, 2, 4); +hold on; +plot(f, coh_uy_yy) +set(gca, 'Xscale', 'log'); +title('Coherence $\frac{y_y}{u_y}$') +xlabel('Frequency [Hz]') +hold off; + +linkaxes([ax11,ax12,ax21,ax22],'x'); +xlim([10, 1000]); +linkaxes([ax11,ax12,ax21,ax22],'y'); +ylim([0, 1]) + +% Extraction of a transfer function matrix +% First we define the initial guess for the resonance frequencies and the weights associated. + +freqs_res = [410, 250]; % [Hz] +freqs_res_weights = [10, 10]; % [Hz] + + + +% From the number of resonance frequency we want to fit, we define the order =N= of the system we want to obtain. + +N = 2*length(freqs_res); + + + +% We then make an initial guess on the complex values of the poles. + +xi = 0.001; % Approximate modal damping +poles = [2*pi*freqs_res*(xi + 1i), 2*pi*freqs_res*(xi - 1i)]; + + + +% We then define the weight that will be used for the fitting. +% Basically, we want more weight around the resonance and at low frequency (below the first resonance). +% Also, we want more importance where we have a better coherence. + +weight = ones(1, length(f)); +% weight = G_coh'; + +% alpha = 0.1; + +% for freq_i = 1:length(freqs_res) +% weight(f>(1-alpha)*freqs_res(freq_i) & omega<(1 + alpha)*2*pi*freqs_res(freq_i)) = freqs_res_weights(freq_i); +% end + + + +% Ignore data above some frequency. + +weight(f>1000) = 0; + +figure; +hold on; +plot(f, weight); +plot(freqs_res, ones(size(freqs_res)), 'rx'); +hold off; +xlabel('Frequency [Hz]'); +xlabel('Weight Amplitude'); +set(gca, 'xscale', 'log'); +xlim([f(1), f(end)]); + + + +% #+NAME: fig:weights +% #+CAPTION: Weights amplitude ([[./figs/weights.png][png]], [[./figs/weights.pdf][pdf]]) +% [[file:figs/weights.png]] + +% When we set some options for =vfit3=. + +opts = struct(); + +opts.stable = 1; % Enforce stable poles +opts.asymp = 1; % Force D matrix to be null +opts.relax = 1; % Use vector fitting with relaxed non-triviality constraint +opts.skip_pole = 0; % Do NOT skip pole identification +opts.skip_res = 0; % Do NOT skip identification of residues (C,D,E) +opts.cmplx_ss = 0; % Create real state space model with block diagonal A + +opts.spy1 = 0; % No plotting for first stage of vector fitting +opts.spy2 = 0; % Create magnitude plot for fitting of f(s) + + + +% We define the number of iteration. + +Niter = 5; + + + +% An we run the =vectfit3= algorithm. + +for iter = 1:Niter + [SER_ux_yx, poles, ~, fit_ux_yx] = vectfit3(tf_ux_yx.', 1i*2*pi*f, poles, weight, opts); +end +for iter = 1:Niter + [SER_uy_yx, poles, ~, fit_uy_yx] = vectfit3(tf_uy_yx.', 1i*2*pi*f, poles, weight, opts); +end +for iter = 1:Niter + [SER_ux_yy, poles, ~, fit_ux_yy] = vectfit3(tf_ux_yy.', 1i*2*pi*f, poles, weight, opts); +end +for iter = 1:Niter + [SER_uy_yy, poles, ~, fit_uy_yy] = vectfit3(tf_uy_yy.', 1i*2*pi*f, poles, weight, opts); +end + +figure; +ax11 = subplot(2, 2, 1); +hold on; +plot(f, abs(tf_ux_yx)) +plot(f, abs(fit_ux_yx)) +title('Frequency Response Function $\frac{y_x}{u_x}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude') +hold off; + +ax12 = subplot(2, 2, 2); +hold on; +plot(f, abs(tf_ux_yy)) +plot(f, abs(fit_ux_yy)) +title('Frequency Response Function $\frac{y_x}{u_y}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +hold off; + +ax21 = subplot(2, 2, 3); +hold on; +plot(f, abs(tf_uy_yx)) +plot(f, abs(fit_uy_yx)) +title('Frequency Response Function $\frac{y_y}{u_x}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude') +xlabel('Frequency [Hz]') +hold off; + +ax22 = subplot(2, 2, 4); +hold on; +plot(f, abs(tf_uy_yy)) +plot(f, abs(fit_uy_yy)) +title('Frequency Response Function $\frac{y_y}{u_y}$') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +xlabel('Frequency [Hz]') +hold off; + +linkaxes([ax11,ax12,ax21,ax22],'x'); +xlim([10, 1000]); +linkaxes([ax11,ax12,ax21,ax22],'y'); +ylim([1e-2, 1e3]) + + + +% #+NAME: fig:identification_matrix_fit +% #+CAPTION: Transfer Function Extraction of the FRF matrix ([[./figs/identification_matrix_fit.png][png]], [[./figs/identification_matrix_fit.pdf][pdf]]) +% [[file:figs/identification_matrix_fit.png]] + +% And finally, we create the identified state space model: + +G_ux_yx = minreal(ss(full(SER_ux_yx.A),SER_ux_yx.B,SER_ux_yx.C,SER_ux_yx.D)); +G_uy_yx = minreal(ss(full(SER_uy_yx.A),SER_uy_yx.B,SER_uy_yx.C,SER_uy_yx.D)); +G_ux_yy = minreal(ss(full(SER_ux_yy.A),SER_ux_yy.B,SER_ux_yy.C,SER_ux_yy.D)); +G_uy_yy = minreal(ss(full(SER_uy_yy.A),SER_uy_yy.B,SER_uy_yy.C,SER_uy_yy.D)); + +G = [G_ux_yx, G_uy_yx; + G_ux_yy, G_uy_yy]; + +save('mat/plant.mat', 'G');