Initial Commit

matlab/charge_voltage_estimation.m

@@ -0,0 +1,330 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+addpath('./mat/');
+
+% Steps
+
+load('force_sensor_steps.mat', 't', 'encoder', 'u', 'v');
+
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+nexttile;
+plot(t, v);
+xlabel('Time [s]'); ylabel('Measured voltage [V]');
+nexttile;
+plot(t, u);
+xlabel('Time [s]'); ylabel('Actuator Voltage [V]');
+
+
+
+% #+name: fig:force_sen_steps_time_domain
+% #+caption: Time domain signal during the 3 actuator voltage steps
+% #+RESULTS:
+% [[file:figs/force_sen_steps_time_domain.png]]
+
+% Three steps are performed at the following time intervals:
+
+t_s = [ 2.5, 23;
+       23.8, 35;
+       35.8, 50];
+
+
+
+% Fit function:
+
+f = @(b,x) b(1).*exp(b(2).*x) + b(3);
+
+
+
+% We are interested by the =b(2)= term, which is the time constant of the exponential.
+
+tau = zeros(size(t_s, 1),1);
+V0  = zeros(size(t_s, 1),1);
+
+for t_i = 1:size(t_s, 1)
+    t_cur = t(t_s(t_i, 1) < t & t < t_s(t_i, 2));
+    t_cur = t_cur - t_cur(1);
+    y_cur = v(t_s(t_i, 1) < t & t < t_s(t_i, 2));
+
+    nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
+    B0 = [0.5, -0.15, 2.2];                        % Choose Appropriate Initial Estimates
+    [B,rnrm] = fminsearch(nrmrsd, B0);     % Estimate Parameters ‘B’
+
+    tau(t_i) = 1/B(2);
+    V0(t_i)  = B(3);
+end
+
+
+
+% #+RESULTS:
+% | $tau$ [s] | $V_0$ [V] |
+% |-----------+-----------|
+% |      6.47 |      2.26 |
+% |      6.76 |      2.26 |
+% |      6.49 |      2.25 |
+
+% With the capacitance being $C = 4.4 \mu F$, the internal impedance of the Speedgoat ADC can be computed as follows:
+
+Cp = 4.4e-6; % [F]
+Rin = abs(mean(tau))/Cp;
+
+ans = Rin
+
+
+
+% #+RESULTS:
+% : 1494100.0
+
+% The input impedance of the Speedgoat's ADC should then be close to $1.5\,M\Omega$ (specified at $1\,M\Omega$).
+
+% #+begin_important
+%   How can we explain the voltage offset?
+% #+end_important
+
+% As shown in Figure [[fig:force_sensor_model_electronics]] (taken from cite:reza06_piezoel_trans_vibrat_contr_dampin), an input voltage offset is due to the input bias current $i_n$.
+
+% #+name: fig:force_sensor_model_electronics
+% #+caption: Model of a piezoelectric transducer (left) and instrumentation amplifier (right)
+% [[file:figs/force_sensor_model_electronics.png]]
+
+% The estimated input bias current is then:
+
+in = mean(V0)/Rin;
+
+ans = in
+
+
+
+% #+RESULTS:
+% : 1.5119e-06
+
+% An additional resistor in parallel with $R_{in}$ would have two effects:
+% - reduce the input voltage offset
+%   \[ V_{off} = \frac{R_a R_{in}}{R_a + R_{in}} i_n \]
+% - increase the high pass corner frequency $f_c$
+%   \[ C_p \frac{R_{in}R_a}{R_{in} + R_a} = \tau_c = \frac{1}{f_c} \]
+%   \[ R_a = \frac{R_i}{f_c C_p R_i - 1} \]
+
+
+% If we allow the high pass corner frequency to be equals to 3Hz:
+
+fc = 3;
+Ra = Rin/(fc*Cp*Rin - 1);
+
+ans = Ra
+
+
+
+% #+RESULTS:
+% : 79804
+
+% With this parallel resistance value, the voltage offset would be:
+
+V_offset = Ra*Rin/(Ra + Rin) * in;
+
+ans = V_offset
+
+% Add Parallel Resistor
+% A resistor $R_p \approx 100\,k\Omega$ is added in parallel with the force sensor as shown in Figure [[fig:force_sensor_model_electronics_without_R]].
+
+% #+name: fig:force_sensor_model_electronics_without_R
+% #+caption: Model of a piezoelectric transducer (left) and instrumentation amplifier (right) with added resistor $R_p$
+% [[file:figs/force_sensor_model_electronics_without_R.png]]
+
+
+load('force_sensor_steps_R_82k7.mat', 't', 'encoder', 'u', 'v');
+
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+nexttile;
+plot(t, v);
+xlabel('Time [s]'); ylabel('Measured voltage [V]');
+nexttile;
+plot(t, u);
+xlabel('Time [s]'); ylabel('Actuator Voltage [V]');
+
+
+
+% #+name: fig:force_sen_steps_time_domain_par_R
+% #+caption: Time domain signal during the actuator voltage steps
+% #+RESULTS:
+% [[file:figs/force_sen_steps_time_domain_par_R.png]]
+
+% Three steps are performed at the following time intervals:
+
+t_s = [1.9,  6;
+       8.5, 13;
+      15.5, 21;
+      22.6, 26;
+      30.0, 36;
+      37.5, 41;
+      46.2, 49.5]
+
+
+
+% Fit function:
+
+f = @(b,x) b(1).*exp(b(2).*x) + b(3);
+
+
+
+% We are interested by the =b(2)= term, which is the time constant of the exponential.
+
+
+tau = zeros(size(t_s, 1),1);
+V0  = zeros(size(t_s, 1),1);
+
+for t_i = 1:size(t_s, 1)
+    t_cur = t(t_s(t_i, 1) < t & t < t_s(t_i, 2));
+    t_cur = t_cur - t_cur(1);
+    y_cur = v(t_s(t_i, 1) < t & t < t_s(t_i, 2));
+
+    nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
+    B0 = [0.5, -0.2, 0.2];                        % Choose Appropriate Initial Estimates
+    [B,rnrm] = fminsearch(nrmrsd, B0);     % Estimate Parameters ‘B’
+
+    tau(t_i) = 1/B(2);
+    V0(t_i)  = B(3);
+end
+
+
+
+% #+RESULTS:
+% | $tau$ [s] | $V_0$ [V] |
+% |-----------+-----------|
+% |      0.43 |      0.15 |
+% |      0.45 |      0.16 |
+% |      0.43 |      0.15 |
+% |      0.43 |      0.15 |
+% |      0.45 |      0.15 |
+% |      0.46 |      0.16 |
+% |      0.48 |      0.16 |
+
+% Knowing the capacitance value, we can estimate the value of the added resistor (neglecting the input impedance of $\approx 1\,M\Omega$):
+
+
+Cp = 4.4e-6; % [F]
+Rin = abs(mean(tau))/Cp;
+
+ans = Rin
+
+
+
+% #+RESULTS:
+% : 101200.0
+
+% And we can verify that the bias current estimation stays the same:
+
+in = mean(V0)/Rin;
+
+ans = in
+
+% Sinus
+
+load('force_sensor_sin.mat', 't', 'encoder', 'u', 'v');
+
+u       = u(t>25);
+v       = v(t>25);
+encoder = encoder(t>25) - mean(encoder(t>25));
+t       = t(t>25);
+
+
+
+% The driving voltage is a sinus at 0.5Hz centered on 3V and with an amplitude of 3V (Figure [[fig:force_sensor_sin_u]]).
+
+
+figure;
+plot(t, u)
+xlabel('Time [s]'); ylabel('Control Voltage [V]');
+
+
+
+% #+name: fig:force_sensor_sin_u
+% #+caption: Driving Voltage
+% #+RESULTS:
+% [[file:figs/force_sensor_sin_u.png]]
+
+% The full stroke as measured by the encoder is:
+
+max(encoder)-min(encoder)
+
+
+
+% #+RESULTS:
+% : 5.005e-05
+
+% Its signal is shown in Figure [[fig:force_sensor_sin_encoder]].
+
+
+figure;
+plot(t, encoder)
+xlabel('Time [s]'); ylabel('Encoder [m]');
+
+
+
+% #+name: fig:force_sensor_sin_encoder
+% #+caption: Encoder measurement
+% #+RESULTS:
+% [[file:figs/force_sensor_sin_encoder.png]]
+
+% The generated voltage by the stack is shown in Figure
+
+
+figure;
+plot(t, v)
+xlabel('Time [s]'); ylabel('Force Sensor Output [V]');
+
+
+
+% #+name: fig:force_sensor_sin_stack
+% #+caption: Voltage measured on the stack used as a sensor
+% #+RESULTS:
+% [[file:figs/force_sensor_sin_stack.png]]
+
+% The capacitance of the stack is
+
+Cp = 4.4e-6; % [F]
+
+
+
+% The corresponding generated charge is then shown in Figure [[fig:force_sensor_sin_charge]].
+
+figure;
+plot(t, 1e6*Cp*(v-mean(v)))
+xlabel('Time [s]'); ylabel('Generated Charge [$\mu C$]');
+
+
+
+% #+name: fig:force_sensor_sin_charge
+% #+caption: Generated Charge
+% #+RESULTS:
+% [[file:figs/force_sensor_sin_charge.png]]
+
+
+% The relation between the generated voltage and the measured displacement is almost linear as shown in Figure [[fig:force_sensor_linear_relation]].
+
+
+b1 = encoder\(v-mean(v));
+
+figure;
+hold on;
+plot(encoder, v-mean(v), 'DisplayName', 'Measured Voltage');
+plot(encoder, encoder*b1, 'DisplayName', sprintf('Linear Fit: $U_s \\approx %.3f [V/\\mu m] \\cdot d$', 1e-6*abs(b1)));
+hold off;
+xlabel('Measured Displacement [m]'); ylabel('Generated Voltage [V]');
+legend();
+
+
+
+% #+name: fig:force_sensor_linear_relation
+% #+caption: Almost linear relation between the relative displacement and the generated voltage
+% #+RESULTS:
+% [[file:figs/force_sensor_linear_relation.png]]
+
+% With a 16bits ADC, the resolution will then be equals to (in [nm]):
+
+abs((20/2^16)/(b1/1e9))

matlab/open_closed_circuit.m

@@ -0,0 +1,67 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+addpath('./mat/');
+
+% Load Data
+
+oc = load('identification_open_circuit.mat', 't', 'encoder', 'u');
+sc = load('identification_short_circuit.mat', 't', 'encoder', 'u');
+
+% Transfer Functions
+
+Ts = 1e-4; % Sampling Time [s]
+win = hann(ceil(10/Ts));
+
+[tf_oc_est, f] = tfestimate(oc.u, oc.encoder, win, [], [], 1/Ts);
+[co_oc_est, ~] = mscohere(  oc.u, oc.encoder, win, [], [], 1/Ts);
+
+[tf_sc_est, ~] = tfestimate(sc.u, sc.encoder, win, [], [], 1/Ts);
+[co_sc_est, ~] = mscohere(  sc.u, sc.encoder, win, [], [], 1/Ts);
+
+figure;
+hold on;
+plot(f, co_oc_est, '-')
+plot(f, co_sc_est, '-')
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
+ylabel('Coherence'); xlabel('Frequency [Hz]');
+hold off;
+xlim([0.5, 5e3]);
+
+
+
+% #+name: fig:stiffness_force_sensor_coherence
+% #+caption:
+% #+RESULTS:
+% [[file:figs/stiffness_force_sensor_coherence.png]]
+
+
+
+figure;
+tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+hold on;
+plot(f, abs(tf_oc_est), '-', 'DisplayName', 'Open-Circuit')
+plot(f, abs(tf_sc_est), '-', 'DisplayName', 'Short-Circuit')
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
+hold off;
+ylim([1e-7, 3e-4]);
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(f, 180/pi*angle(tf_oc_est), '-')
+plot(f, 180/pi*angle(tf_sc_est), '-')
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
+ylabel('Phase'); xlabel('Frequency [Hz]');
+hold off;
+yticks(-360:90:360);
+axis padded 'auto x'
+
+linkaxes([ax1,ax2], 'x');
+xlim([0.5, 5e3]);

matlab/runtest.m

@@ -0,0 +1,17 @@
+tg = slrt;
+
+%%
+f = SimulinkRealTime.openFTP(tg);
+mget(f, 'int_enc.dat', 'data');
+close(f);
+
+%% Convert the Data
+data = SimulinkRealTime.utils.getFileScopeData('data/int_enc.dat').data;
+
+interferometer = data(:, 1);
+encoder = data(:, 2);
+u = data(:, 3);
+t = data(:, 4);
+
+save('./mat/int_enc_id_noise_bis.mat', 'interferometer', 'encoder', 'u' , 't');
+

matlab/setup.m

@@ -0,0 +1,23 @@
+%%
+s = tf('s');
+Ts = 1e-4; % [s]
+
+%% Pre-Filter
+G_pf = 1/(1 + s/2/pi/20);
+G_pf = c2d(G_pf, Ts, 'tustin');
+
+% %% Force Sensor Filter (HPF)
+% Gf_hpf = s/(s + 2*pi*2);
+% Gf_hpf = tf(1);
+% Gf_hpf = c2d(Gf_hpf, Ts, 'tustin');
+% 
+% %% IFF Controller
+% Kiff = 1/(s + 2*pi*2);
+% Kiff = c2d(Kiff, Ts, 'tustin');
+% 
+% %% Excitation Signal
+Tsim = 100; % Excitation time + Measurement time [s]
+
+t = 0:Ts:Tsim;
+u_exc = timeseries(chirp(t, 10, Tsim, 40, 'logarithmic'), t);
+% u_exc = timeseries(y_v, t);