diff --git a/figs/apa95ml_5kg_PI_coh.pdf b/figs/apa95ml_5kg_PI_coh.pdf index 4610359..13b7844 100644 Binary files a/figs/apa95ml_5kg_PI_coh.pdf and b/figs/apa95ml_5kg_PI_coh.pdf differ diff --git a/figs/apa95ml_5kg_PI_coh.png b/figs/apa95ml_5kg_PI_coh.png index 78401d7..758c2d0 100644 Binary files a/figs/apa95ml_5kg_PI_coh.png and b/figs/apa95ml_5kg_PI_coh.png differ diff --git a/figs/apa95ml_5kg_cedrat_coh.pdf b/figs/apa95ml_5kg_cedrat_coh.pdf new file mode 100644 index 0000000..13a757e Binary files /dev/null and b/figs/apa95ml_5kg_cedrat_coh.pdf differ diff --git a/figs/apa95ml_5kg_cedrat_coh.png b/figs/apa95ml_5kg_cedrat_coh.png new file mode 100644 index 0000000..34cc10a Binary files /dev/null and b/figs/apa95ml_5kg_cedrat_coh.png differ diff --git a/figs/apa95ml_5kg_pi_comp_fem.pdf b/figs/apa95ml_5kg_pi_comp_fem.pdf index adc911a..1464441 100644 Binary files a/figs/apa95ml_5kg_pi_comp_fem.pdf and b/figs/apa95ml_5kg_pi_comp_fem.pdf differ diff --git a/figs/apa95ml_5kg_pi_comp_fem.png b/figs/apa95ml_5kg_pi_comp_fem.png index 3822d06..80250cf 100644 Binary files a/figs/apa95ml_5kg_pi_comp_fem.png and b/figs/apa95ml_5kg_pi_comp_fem.png differ diff --git a/figs/bode_plot_force_sensor_voltage_comp_fem.pdf b/figs/bode_plot_force_sensor_voltage_comp_fem.pdf index cfee0dc..6b26237 100644 Binary files a/figs/bode_plot_force_sensor_voltage_comp_fem.pdf and b/figs/bode_plot_force_sensor_voltage_comp_fem.pdf differ diff --git a/figs/bode_plot_force_sensor_voltage_comp_fem.png b/figs/bode_plot_force_sensor_voltage_comp_fem.png index 06729ea..093b7ec 100644 Binary files a/figs/bode_plot_force_sensor_voltage_comp_fem.png and b/figs/bode_plot_force_sensor_voltage_comp_fem.png differ diff --git a/figs/dynamics_act_disp_comp_mass.pdf b/figs/dynamics_act_disp_comp_mass.pdf index cfad4c3..7f3db59 100644 Binary files a/figs/dynamics_act_disp_comp_mass.pdf and b/figs/dynamics_act_disp_comp_mass.pdf differ diff --git a/figs/dynamics_act_disp_comp_mass.png b/figs/dynamics_act_disp_comp_mass.png index 94202da..f73ae4f 100644 Binary files a/figs/dynamics_act_disp_comp_mass.png and b/figs/dynamics_act_disp_comp_mass.png differ diff --git a/figs/dynamics_force_force_sensor_comp_mass.pdf b/figs/dynamics_force_force_sensor_comp_mass.pdf index f0b9014..371f92f 100644 Binary files a/figs/dynamics_force_force_sensor_comp_mass.pdf and b/figs/dynamics_force_force_sensor_comp_mass.pdf differ diff --git a/figs/dynamics_force_force_sensor_comp_mass.png b/figs/dynamics_force_force_sensor_comp_mass.png index 9a60c15..cbc8dda 100644 Binary files a/figs/dynamics_force_force_sensor_comp_mass.png and b/figs/dynamics_force_force_sensor_comp_mass.png differ diff --git a/figs/iff_plant_identification_apa95ml.pdf b/figs/iff_plant_identification_apa95ml.pdf index b945999..501db0b 100644 Binary files a/figs/iff_plant_identification_apa95ml.pdf and b/figs/iff_plant_identification_apa95ml.pdf differ diff --git a/figs/iff_plant_identification_apa95ml.png b/figs/iff_plant_identification_apa95ml.png index 8d004f2..8296582 100644 Binary files a/figs/iff_plant_identification_apa95ml.png and b/figs/iff_plant_identification_apa95ml.png differ diff --git a/figs/test_bench_apa_identification.pdf b/figs/test_bench_apa_identification.pdf new file mode 100644 index 0000000..1ce4148 Binary files /dev/null and b/figs/test_bench_apa_identification.pdf differ diff --git a/figs/test_bench_apa_identification.png b/figs/test_bench_apa_identification.png new file mode 100644 index 0000000..d7f8fff Binary files /dev/null and b/figs/test_bench_apa_identification.png differ diff --git a/figs/test_bench_apa_identification.svg b/figs/test_bench_apa_identification.svg new file mode 100644 index 0000000..414213f --- /dev/null +++ b/figs/test_bench_apa_identification.svg @@ -0,0 +1,4099 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + PI - E505 + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/figs/test_bench_apa_schematic.pdf b/figs/test_bench_apa_schematic.pdf index 4ea5f7d..5215316 100644 Binary files a/figs/test_bench_apa_schematic.pdf and b/figs/test_bench_apa_schematic.pdf differ diff --git a/figs/test_bench_apa_schematic.png b/figs/test_bench_apa_schematic.png index d2740fd..80651fc 100644 Binary files a/figs/test_bench_apa_schematic.png and b/figs/test_bench_apa_schematic.png differ diff --git a/figs/test_bench_apa_schematic.svg b/figs/test_bench_apa_schematic.svg index e9f2942..46184a0 100644 --- a/figs/test_bench_apa_schematic.svg +++ b/figs/test_bench_apa_schematic.svg @@ -15,7 +15,7 @@ viewBox="0 0 445.25269 309.88137" sodipodi:docname="test_bench_apa_schematic.svg" inkscape:version="1.0.1 (3bc2e813f5, 2020-09-07)" - inkscape:export-filename="/home/thomas/Cloud/thesis/matlab/encoder-test-bench/figs/exp_setup_schematic.png" + inkscape:export-filename="/home/thomas/Cloud/thesis/matlab/test-bench-apa/figs/test_bench_apa_schematic.png" inkscape:export-xdpi="252" inkscape:export-ydpi="252"> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/figs/test_bench_apa_schematic_iff.pdf b/figs/test_bench_apa_schematic_iff.pdf index 441913e..c8cb812 100644 Binary files a/figs/test_bench_apa_schematic_iff.pdf and b/figs/test_bench_apa_schematic_iff.pdf differ diff --git a/figs/test_bench_apa_schematic_iff.png b/figs/test_bench_apa_schematic_iff.png index 542941b..2d87688 100644 Binary files a/figs/test_bench_apa_schematic_iff.png and b/figs/test_bench_apa_schematic_iff.png differ diff --git a/figs/test_bench_apa_schematic_iff.svg b/figs/test_bench_apa_schematic_iff.svg index f6cdfe9..3f30f41 100644 --- a/figs/test_bench_apa_schematic_iff.svg +++ b/figs/test_bench_apa_schematic_iff.svg @@ -26,7 +26,7 @@ image/svg+xml - + @@ -2706,18 +2706,88 @@ id="path5148" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/index.html b/index.html index 17de956..83c489d 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Test Bench - Amplified Piezoelectric Actuator @@ -30,75 +30,76 @@

Table of Contents

+

+This document is divided in the following sections: +

-
-

1 Experimental Setup

+
+

1 Experimental Setup

- +

-A schematic of the test-bench is shown in Figure 1. +A schematic of the test-bench is shown in Figure 1.

@@ -111,31 +112,31 @@ The APA95ML has three stacks that can be used as actuator or as sensors.

-Pictures of the test bench are shown in Figure 2 and 3. +Pictures of the test bench are shown in Figure 2 and 3.

-
+

test_bench_apa_schematic.png

Figure 1: Schematic of the Setup

-
+

setup_picture.png

Figure 2: Picture of the Setup

-
+

setup_zoom.png

Figure 3: Zoom on the APA

-
+

Here are the equipment used in the test bench:

@@ -151,11 +152,14 @@ Here are the equipment used in the test bench:
-
-

2 Estimation of electrical/mechanical relationships

+
+

2 Estimation of electrical/mechanical relationships

-In order to correctly model the piezoelectric actuator, we need to determine: + +

+

+In order to correctly model the piezoelectric actuator with Simscape, we need to determine:

  1. \(g_a\): the ratio of the generated force \(F_a\) to the supply voltage \(V_a\) across the piezoelectric stack
    2. @@ -166,24 +170,23 @@ In order to correctly model the piezoelectric actuator, we need to determine: We estimate \(g_a\) and \(g_s\) using different approaches:

      -
    1. Section 2.1: \(g_a\) is estimated from the datasheet of the piezoelectric stack
      2. -
    2. Section 2.2: \(g_a\) and \(g_s\) are estimated using the piezoelectric constants
      3. -
    3. Section 2.3: \(g_a\) and \(g_s\) are estimated experimentally
      4. +
    4. Section 2.1: \(g_a\) is estimated from the datasheet of the piezoelectric stack
      5. +
    5. Section 2: \(g_a\) and \(g_s\) are estimated using the piezoelectric constants
      6. +
    6. Section 2.3: \(g_a\) and \(g_s\) are estimated experimentally
- -
-

2.1 Estimation from Data-sheet

+
+

2.1 Estimation from Data-sheet

- +

-The stack parameters taken from the data-sheet are shown in Table 1. +The stack parameters taken from the data-sheet are shown in Table 1.

- +
@@ -269,12 +272,12 @@ Vmax = 170; % [V]
-
ka*Lmax/Vmax % [N/V]
+
ga = ka*Lmax/Vmax; % [N/V]
 

 
-27.647
+ga = 27.6 [N/V]
@@ -284,11 +287,11 @@ From the parameters of the stack, it seems not possible to estimate the relation -
-

2.2 Estimation from Piezoelectric parameters

+
+

2.2 Estimation from Piezoelectric parameters

- +

@@ -305,7 +308,7 @@ ka = 235e6; % Stack stiffness [N/m] The ratio of the developed force to applied voltage is:

\begin{equation} -\label{org8300de8} +\label{org6f1476c} F_a = g_a V_a, \quad g_a = d_{33} n k_a \end{equation}

@@ -337,7 +340,7 @@ ga = 5.6 [N/V] From (Fleming and Leang 2014) (page 123), the relation between relative displacement of the sensor stack and generated voltage is:

\begin{equation} -\label{org312083f} +\label{org43f15ff} V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h \end{equation}

@@ -374,11 +377,11 @@ gs = 35.4 [V/um]

-
-

2.3 Estimation from Experiment

+
+

2.3 Estimation from Experiment

- +

The idea here is to obtain the parameters \(g_a\) and \(g_s\) from the comparison of an experimental identification and the identification using Simscape. @@ -394,8 +397,8 @@ Similarly, it is fairly easy to experimentally obtain the gain from the stack di To link that to the strain of the sensor stack, the simscape model is used.

-
-

2.3.1 From actuator voltage \(V_a\) to actuator force \(F_a\)

+
+

2.3.1 From actuator voltage \(V_a\) to actuator force \(F_a\)

The data from the identification test is loaded. @@ -439,7 +442,7 @@ The gain from input voltage of the stack to the vertical displacement is determi

-
+

gain_Va_to_d.png

Figure 4: Transfer function from actuator stack voltage \(V_a\) to vertical displacement of the mass \(d\)

@@ -495,11 +498,11 @@ ga = 33.7 [N/V]

-The obtained comparison between the Simscape model and the identified dynamics is shown in Figure 5. +The obtained comparison between the Simscape model and the identified dynamics is shown in Figure 5.

-
+

compare_Gd_id_simscape.png

Figure 5: Comparison of the identified transfer function between actuator voltage \(V_a\) and vertical mass displacement \(d\)

@@ -507,8 +510,8 @@ The obtained comparison between the Simscape model and the identified dynamics i
-
-

2.3.2 From stack strain \(\Delta h\) to generated voltage \(V_s\)

+
+

2.3.2 From stack strain \(\Delta h\) to generated voltage \(V_s\)

Now, the gain from the stack strain \(\Delta h\) to the generated voltage \(V_s\) is estimated. @@ -543,7 +546,7 @@ Here, an amplifier with a gain of 20 is used.

-Then, the transfer function from \(V_a\) to \(V_s\) is identified and its DC gain is estimated (Figure 6). +Then, the transfer function from \(V_a\) to \(V_s\) is identified and its DC gain is estimated (Figure 6). 

Ts = t(end)/(length(t)-1);
@@ -561,7 +564,7 @@ win = hann(ceil(10/Ts));
-
+

gain_Va_to_Vs.png

Figure 6: Transfer function from actuator stack voltage \(V_a\) to sensor stack voltage \(V_s\)

@@ -617,7 +620,7 @@ gs = 3.5 [V/um] -
+

compare_Gf_id_simscape.png

Figure 7: Comparison of the identified transfer function between actuator voltage \(V_a\) and sensor stack voltage

@@ -626,8 +629,8 @@ gs = 3.5 [V/um]
-
-

2.4 Conclusion

+
+

2.4 Conclusion

The obtained parameters \(g_a\) and \(g_s\) are not consistent between the different methods. @@ -644,11 +647,11 @@ The one using the experimental data are saved and further used.

-
-

3 Simscape model of the test-bench

+
+

3 Simscape model of the test-bench

- +

The idea here is to model the test-bench using Simscape. @@ -659,19 +662,19 @@ Whereas the suspended mass and metrology frame can be considered as rigid bodies

-To model the APA, a Finite Element Model (FEM) is used (Figure 8) and imported into Simscape. +To model the APA, a Finite Element Model (FEM) is used (Figure 8) and imported into Simscape.

-
+

APA95ML_FEM.png

Figure 8: Finite Element Model of the APA95ML

-
-

3.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

+
+

3.1 Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates

We first extract the stiffness and mass matrices. @@ -997,7 +1000,6 @@ M = readmatrix('APA95ML_M.CSV');

Table 1: Stack Parameters
-

Then, we extract the coordinates of the interface nodes.

@@ -1007,7 +1009,7 @@ Then, we extract the coordinates of the interface nodes.

-The interface nodes are shown in Figure 9 and their coordinates are listed in Table 4. +The interface nodes are shown in Figure 9 and their coordinates are listed in Table 4.

@@ -1041,7 +1043,7 @@ The interface nodes are shown in Figure 9 and their co
- +
@@ -1124,7 +1126,7 @@ The interface nodes are shown in Figure 9 and their co
Table 4: Coordinates of the interface nodes
-
+

APA95ML_nodes_1.png

Figure 9: Interface Nodes chosen for the APA95ML

@@ -1136,8 +1138,8 @@ Using K, M and int_xyz, we can use the
-
-

3.2 Simscape Model

+
+

3.2 Simscape Model

The flexible element is imported using the Reduced Order Flexible Solid Simscape block. @@ -1152,17 +1154,22 @@ A Relative Motion Sensor block is added between the nodes 1 and 2 t One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame.

-
m = 5;
+
m = 5.5;
+

+
+ +
+
load('apa95ml_params.mat', 'ga', 'gs');
-
-

3.3 Dynamics from Actuator Voltage to Vertical Mass Displacement

+
+

3.3 Dynamics from Actuator Voltage to Vertical Mass Displacement

-The identified dynamics is shown in Figure 10. +The identified dynamics is shown in Figure 10.

@@ -1174,12 +1181,12 @@ clear io; io_i = 1; io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % Actuator Voltage [V] io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m] -Ghm = -linearize(mdl, io); +Ghm = linearize(mdl, io);
-
+

dynamics_act_disp_comp_mass.png

Figure 10: Dynamics from \(F\) to \(d\) without a payload and with a 5kg payload

@@ -1187,11 +1194,11 @@ Ghm = -linearize(mdl, io);
-
-

3.4 Dynamics from Actuator Voltage to Force Sensor Voltage

+
+

3.4 Dynamics from Actuator Voltage to Force Sensor Voltage

-The obtained dynamics is shown in Figure 11. +The obtained dynamics is shown in Figure 11.

@@ -1208,7 +1215,7 @@ Gfm = linearize(mdl, io);
-
+

dynamics_force_force_sensor_comp_mass.png

Figure 11: Dynamics from \(F\) to \(F_m\) for \(m=0\) and \(m = 10kg\)

@@ -1216,34 +1223,36 @@ Gfm = linearize(mdl, io);
-
-

3.5 Save Data for further use

+
+

3.5 Save Data for further use

save('matlab/mat/fem_simscape_models.mat', 'Ghm', 'Gfm')
- -
-
save('mat/fem_simscape_models.mat', 'Ghm', 'Gfm')
-
-
-
-

4 Huddle Test

+
+

4 Measurement of the ambient noise in the system

- + +

+

+This first measurement consist of measuring the displacement of the mass using the interferometer when no voltage is applied to the actuator. +

+ +

+This can help determining the actuator voltage necessary to generate a motion way above the measured noise and disturbances, and thus obtain a good identification.

-
-

4.1 Time Domain Data

+
+

4.1 Time Domain Data

-
+

huddle_test_time_domain.png

Figure 12: Measurement of the Mass displacement during Huddle Test

@@ -1251,8 +1260,8 @@ Gfm = linearize(mdl, io);
-
-

4.2 PSD of Measurement Noise

+
+

4.2 PSD of Measurement Noise

Ts = t(end)/(length(t)-1);
@@ -1268,7 +1277,7 @@ win = hanning(ceil(1*Fs));
-
+

huddle_test_pdf.png

Figure 13: Amplitude Spectral Density of the Displacement during Huddle Test

@@ -1277,24 +1286,30 @@ win = hanning(ceil(1*Fs));
-
-

5 Identification of the dynamics from actuator to displacement

+
+

5 Identification of the dynamics from actuator Voltage to displacement

- + +

+

+The setup used for the identification of the dynamics from \(V_a\) to \(d\) is shown in Figure 14.

-
    -
  • [ ] List of equipment
    • -
  • [ ] Schematic
    • -
  • [ ] Problem of matching between the models? (there is a factor 10)
    • -

-E505 with gain of 10. +A Voltage amplifier with a gain of \(10\) is used. +Two stacks are used as actuators while one stack is used as a force sensor.

+ + +
+

test_bench_apa_identification.png +

+

Figure 14: Test Bench used for the identification of the dynaimcs from \(V_a\) to \(d\)

-
-

5.1 Load Data

+
+
+

5.1 Load Data

The data from the “noise test” and the identification test are loaded. @@ -1322,14 +1337,13 @@ Now we add a factor 10 to take into account the gain of the voltage amplifier. 

um   = 10*um;
-ht.u = 10*ht.u;
-
-

5.2 Comparison of the PSD with Huddle Test

+
+

5.2 Comparison of the PSD with Huddle Test

Ts = t(end)/(length(t)-1);
@@ -1346,16 +1360,16 @@ win = hanning(ceil(1*Fs));
-
+

apa95ml_5kg_PI_pdf_comp_huddle.png

-

Figure 14: Comparison of the ASD for the identification test and the huddle test

+

Figure 15: Comparison of the ASD for the identification test and the huddle test

-
-

5.3 Compute TF estimate and Coherence

+
+

5.3 Compute TF estimate and Coherence

[tf_est, f] = tfestimate(um, -y, win, [], [], 1/Ts);
@@ -1364,10 +1378,10 @@ win = hanning(ceil(1*Fs));
-
+

apa95ml_5kg_PI_coh.png

-

Figure 15: Coherence

+

Figure 16: Coherence

@@ -1379,31 +1393,23 @@ Comparison with the FEM model

-
+

apa95ml_5kg_pi_comp_fem.png

-

Figure 16: Comparison of the identified transfer function and the one estimated from the FE model

+

Figure 17: Comparison of the identified transfer function and the one estimated from the FE model

-
-

6 Identification of the dynamics from actuator to force sensor

+
+

6 Identification of the dynamics from actuator Voltage to force sensor Voltage

- +

-Two measurements are performed: -

-
    -
  • Speedgoat DAC => Voltage Amplifier (x20) => 1 Piezo Stack => … => 2 Stacks as Force Sensor (parallel) => Speedgoat ADC
    • -
  • Speedgoat DAC => Voltage Amplifier (x20) => 2 Piezo Stacks (parallel) => … => 1 Stack as Force Sensor => Speedgoat ADC
    • -
- -

-The obtained dynamics from force actuator to force sensor are compare with the FEM model. +The same setup shown in Figure 14 is used.

The data are loaded: @@ -1413,76 +1419,88 @@ The data are loaded:

-
-
u = detrend(u, 0);
-v = detrend(v, 0);
-
-
- -
-
u = 20*u;
-
-

-Let’s use the amplifier gain to obtain the true voltage applied to the actuator stack(s) -

- -

-The parameters of the piezoelectric stacks are defined below: +Any offset is removed.

-
d33 = 3e-10; % Strain constant [m/V]
-n = 80; % Number of layers per stack
-eT = 1.6e-8; % Permittivity under constant stress [F/m]
-sD = 2e-11; % Elastic compliance under constant electric displacement [m2/N]
-ka = 235e6; % Stack stiffness [N/m]
+
u = detrend(u, 0); % Speedgoat DAC output Voltage [V]
+v = detrend(v, 0); % Speedgoat ADC input Voltage (sensor stack) [V]

-From the FEM, we construct the transfer function from DAC voltage to ADC voltage. +Here, the amplifier gain is 20.

-
Gfem_aa_s = exp(-s/1e4)*20*(2*d33*n*ka)*(G(3,1)+G(3,2))*d33/(eT*sD*n);
-Gfem_a_ss = exp(-s/1e4)*20*(  d33*n*ka)*(G(3,1)+G(2,1))*d33/(eT*sD*n);
-
-
- -
-
Gfem_aa_s = exp(-s/1e4)*20*(2*d33*n*ka)*Gfm*d33/(eT*sD*n);
-Gfem_a_ss = exp(-s/1e4)*20*(  d33*n*ka)*Gfm*d33/(eT*sD*n);
+
u = 20*u; % Actuator Stack Voltage [V]

-The transfer function from input voltage to output voltage are computed and shown in Figure 17. +The transfer function from the actuator voltage \(V_a\) to the force sensor stack voltage \(V_s\) is computed. 

Ts = t(end)/(length(t)-1);
 Fs = 1/Ts;
 
-win = hann(ceil(10/Ts));
+win = hann(ceil(5/Ts));
 
 [tf_est,  f] = tfestimate(u, v, win, [], [], 1/Ts);
 [coh,     ~] = mscohere(  u, v, win, [], [], 1/Ts);
+

+The coherence is shown in Figure 18. +

+ + +
+

apa95ml_5kg_cedrat_coh.png +

+

Figure 18: Coherence

+
+ +

+The Simscape model is loaded and compared with the identified dynamics in Figure 19. +The non-minimum phase zero is just a side effect of the not so great identification. +Taking longer measurements would results in a minimum phase zero. +

+ 
load('mat/fem_simscape_models.mat', 'Gfm');
-
+

bode_plot_force_sensor_voltage_comp_fem.png

-

Figure 17: Comparison of the identified dynamics from voltage output to voltage input and the FEM

+

Figure 19: Comparison of the identified dynamics from voltage output to voltage input and the FEM

-
-

6.1 System Identification

-
+
+ +
+

7 Integral Force Feedback

+
+

+ +

+ +
+

test_bench_apa_schematic_iff.png +

+

Figure 20: Schematic of the test bench using IFF

+
+
+
+

7.1 IFF Plant

+
+

+From the identified plant, a model of the transfer function from the actuator stack voltage to the force sensor generated voltage is developed. +

+ 
w_z = 2*pi*111; % Zeros frequency [rad/s]
 w_p = 2*pi*255; % Pole frequency [rad/s]
@@ -1490,50 +1508,63 @@ xi_z = 0.05;
 xi_p = 0.015;
 G_inf = 0.1;
 
-Gi = G_inf*(s^2 - 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2);
+Gi = G_inf*(s^2 + 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2);
+

+Its bode plot is shown in Figure 21. +

-
+ +

iff_plant_identification_apa95ml.png

-

Figure 18: Identification of the IFF plant

-
-
+

Figure 21: Bode plot of the IFF plant

+

+The controller used in the Integral Force Feedback Architecture is: +

+\begin{equation} + K_{\text{IFF}}(s) = \frac{g}{s + 2\cdot 2\pi} \cdot \frac{s}{s + 0.5 \cdot 2\pi} +\end{equation} +

+where \(g\) is a gain that can be tuned. +

-
-

6.2 Integral Force Feedback

-
+

+Above 2 Hz the controller is basically an integrator, whereas an high pass filter is added at 0.5Hz to further reduce the low frequency gain. +

-
+

+In the frequency band of interest, this controller should mostly act as a pure integrator. +

+ +

+The Root Locus corresponding to this controller is shown in Figure 22. +

+ + +

root_locus_iff_apa95ml_identification.png

-

Figure 19: Root Locus for IFF

-
+

Figure 22: Root Locus for IFF

-
-

7 Integral Force Feedback

-
+ +
+

7.2 First tests with few gains

+

- +The controller is now implemented in practice, and few controller gains are tested: \(g = 0\), \(g = 10\) and \(g = 100\).

-
-

test_bench_apa_schematic_iff.png +

+For each controller gain, the identification shown in Figure 20 is performed.

-

Figure 20: Schematic of the test bench using IFF

-
-
- -
-

7.1 First tests with few gains

-
iff_g10  = load('apa95ml_iff_g10_res.mat',  'u', 't', 'y', 'v');
 iff_g100 = load('apa95ml_iff_g100_res.mat', 'u', 't', 'y', 'v');
@@ -1556,26 +1587,37 @@ win = hann(ceil(10/Ts));
+

+The coherence between the excitation signal and the mass displacement as measured by the interferometer is shown in Figure 23. +

-
+ +

iff_first_test_coherence.png

-

Figure 21: Coherence

+

Figure 23: Coherence

+

+The obtained transfer functions are shown in Figure 24. +It is clear that the IFF architecture can actively damp the main resonance of the system. +

-
+

iff_first_test_bode_plot.png

-

Figure 22: Bode plot for different values of IFF gain

+

Figure 24: Bode plot for different values of IFF gain

-
-

7.2 Second test with many Gains

-
+
+

7.3 Second test with many Gains

+
+

+Then, the same identification test is performed for many more gains. +

load('apa95ml_iff_test.mat', 'results');
@@ -1602,12 +1644,21 @@ g_iff = [0, 1, 5, 10, 50, 100];
+

+The obtained dynamics are shown in Figure 25. +

-
+ +

iff_results_bode_plots.png

+

Figure 25: Identified dynamics from excitation voltage to the mass displacement

+

+For each gain, the parameters of a second order resonant system that best fits the data are estimated and are compared with the data in Figure 26. +

+ 
G_id = {zeros(1,length(results))};
 
@@ -1625,15 +1676,21 @@ f_end = 500; % [Hz]
-
+

iff_results_bode_plots_identification.png

+

Figure 26: Comparison of the measured dynamic and the identified 2nd order resonant systems that best fits the data

+

+Finally, we can represent the position of the poles of the 2nd order systems on the Root Locus plot (Figure 27). +

-
+ +

iff_results_root_locus.png

+

Figure 27: Root Locus plot of the identified IFF plant as well as the identified poles of the damped system

Bibliography

@@ -1646,7 +1703,7 @@ f_end = 500; % [Hz]

Author: Dehaeze Thomas

-

Created: 2020-11-24 mar. 13:54

+

Created: 2020-11-24 mar. 18:24

diff --git a/index.org b/index.org index 27dca3c..f587ce9 100644 --- a/index.org +++ b/index.org @@ -28,11 +28,14 @@ * Introduction :ignore: -- Section [[sec:experimental_setup]]: -- Section [[sec:simscape_model]]: -- Section [[]]: -- Section [[]]: -- Section [[]]: +This document is divided in the following sections: +- Section [[sec:experimental_setup]]: the experimental setup is described +- Section [[sec:estimation_piezo_params]]: the parameters which are important for the Simscape model of the piezoelectric stack actuator and sensors are estimated +- Section [[sec:simscape_model]]: the Simscape model of the test bench is presented +- Section [[sec:huddle_test]]: as usual, a first measurement of the noise/disturbances present in the system is performed +- Section [[sec:motion_identification]]: the transfer function from the actuator voltage to the displacement of the mass is identified and compared with the model +- Section [[sec:force_sensor_identification]]: the tranfer function from the actuator voltage to the sensor stack voltage is identified and compare with the model +- Section [[sec:integral_force_feedback]]: the Integral Force Feedback control architecture is applied on the system using the force sensor stack in order to add damping to the suspension resonance * Experimental Setup <> @@ -68,6 +71,11 @@ Here are the equipment used in the test bench: #+end_note * Estimation of electrical/mechanical relationships +:PROPERTIES: +:header-args:matlab+: :tangle matlab/piezo_parameters.m +:header-args:matlab+: :comments org :mkdirp yes +:END: +<> ** Introduction :ignore: In order to correctly model the piezoelectric actuator with Simscape, we need to determine: 1. $g_a$: the ratio of the generated force $F_a$ to the supply voltage $V_a$ across the piezoelectric stack @@ -130,11 +138,15 @@ The relation between the applied voltage and the generated force can be estimate #+end_src #+begin_src matlab :results value replace - ka*Lmax/Vmax % [N/V] + ga = ka*Lmax/Vmax; % [N/V] +#+end_src + +#+begin_src matlab :results value replace :exports results :tangle no + sprintf('ga = %.1f [N/V]', ga) #+end_src #+RESULTS: -: 27.647 +: ga = 27.6 [N/V] From the parameters of the stack, it seems not possible to estimate the relation between the strain and the generated voltage. @@ -564,7 +576,6 @@ We first extract the stiffness and mass matrices. | 6e-07 | -4e-08 | 0.0003 | -3e-10 | -2e-09 | -2e-09 | -2e-06 | -9e-07 | 0.02 | 2e-08 | | -8e-09 | 5e-06 | 1e-09 | 3e-08 | 7e-11 | 2e-13 | 8e-09 | -9e-05 | 2e-08 | 1e-06 | - Then, we extract the coordinates of the interface nodes. #+begin_src matlab [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); @@ -613,7 +624,11 @@ A =Relative Motion Sensor= block is added between the nodes 1 and 2 to measure t One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame. #+begin_src matlab - m = 5; + m = 5.5; +#+end_src + +#+begin_src matlab + load('apa95ml_params.mat', 'ga', 'gs'); #+end_src ** Dynamics from Actuator Voltage to Vertical Mass Displacement @@ -628,7 +643,7 @@ The identified dynamics is shown in Figure [[fig:dynamics_act_disp_comp_mass]]. io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % Actuator Voltage [V] io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m] - Ghm = -linearize(mdl, io); + Ghm = linearize(mdl, io); #+end_src #+begin_src matlab :exports none @@ -642,7 +657,7 @@ The identified dynamics is shown in Figure [[fig:dynamics_act_disp_comp_mass]]. plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz')))); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); + ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); hold off; ax2 = nexttile; @@ -650,7 +665,6 @@ The identified dynamics is shown in Figure [[fig:dynamics_act_disp_comp_mass]]. plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz'))))); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); yticks(-360:90:360); - ylim([-360 0]); xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); hold off; linkaxes([ax1,ax2],'x'); @@ -692,15 +706,14 @@ The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_ plot(freqs, abs(squeeze(freqresp(Gfm, freqs, 'Hz')))); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); + ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); hold off; ax2 = nexttile; hold on; plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gfm, freqs, 'Hz'))))); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); - yticks(-360:90:360); - ylim([-390 30]); + yticks(-360:90:360); ylim([-360, 180]); xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); hold off; linkaxes([ax1,ax2],'x'); @@ -721,17 +734,21 @@ The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_ save('matlab/mat/fem_simscape_models.mat', 'Ghm', 'Gfm') #+end_src -#+begin_src matlab :eval no +#+begin_src matlab :exports none :eval no save('mat/fem_simscape_models.mat', 'Ghm', 'Gfm') #+end_src -* Huddle Test +* Measurement of the ambient noise in the system :PROPERTIES: :header-args:matlab+: :tangle matlab/huddle_test.m :header-args:matlab+: :comments org :mkdirp yes :END: <> ** Introduction :ignore: +This first measurement consist of measuring the displacement of the mass using the interferometer when no voltage is applied to the actuator. + +This can help determining the actuator voltage necessary to generate a motion way above the measured noise and disturbances, and thus obtain a good identification. + ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) <> @@ -803,7 +820,7 @@ The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_ #+RESULTS: [[file:figs/huddle_test_pdf.png]] -* TODO Identification of the dynamics from actuator to displacement +* Identification of the dynamics from actuator Voltage to displacement :PROPERTIES: :header-args:matlab+: :tangle matlab/motion_identification.m :header-args:matlab+: :comments org :mkdirp yes @@ -811,11 +828,14 @@ The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_ <> ** Introduction :ignore: -- [ ] List of equipment -- [ ] Schematic -- [ ] Problem of matching between the models? (there is a factor 10) +The setup used for the identification of the dynamics from $V_a$ to $d$ is shown in Figure [[fig:test_bench_apa_identification]]. -E505 with gain of 10. +A Voltage amplifier with a gain of $10$ is used. +Two stacks are used as actuators while one stack is used as a force sensor. + +#+name: fig:test_bench_apa_identification +#+caption: Test Bench used for the identification of the dynaimcs from $V_a$ to $d$ +[[file:figs/test_bench_apa_identification.png]] ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) @@ -837,23 +857,21 @@ E505 with gain of 10. ** Load Data The data from the "noise test" and the identification test are loaded. #+begin_src matlab - ht = load('huddle_test.mat', 't', 'u', 'y'); + ht = load('huddle_test.mat', 't', 'y'); load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y'); #+end_src Any offset value is removed: #+begin_src matlab - um = detrend(um, 0); % Input Voltage [V] + um = detrend(um, 0); % Generated DAC Voltage [V] y = detrend(y , 0); % Mass displacement [m] - ht.u = detrend(ht.u, 0); ht.y = detrend(ht.y, 0); #+end_src Now we add a factor 10 to take into account the gain of the voltage amplifier. #+begin_src matlab - um = 10*um; - ht.u = 10*ht.u; + Va = 10*um; #+end_src ** Comparison of the PSD with Huddle Test @@ -892,8 +910,8 @@ Now we add a factor 10 to take into account the gain of the voltage amplifier. ** Compute TF estimate and Coherence #+begin_src matlab - [tf_est, f] = tfestimate(um, -y, win, [], [], 1/Ts); - [co_est, ~] = mscohere( um, -y, win, [], [], 1/Ts); + [tf_est, f] = tfestimate(Va, -y, win, [], [], 1/Ts); + [co_est, ~] = mscohere( Va, -y, win, [], [], 1/Ts); #+end_src #+begin_src matlab :exports none @@ -951,7 +969,7 @@ Comparison with the FEM model #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/apa95ml_5kg_pi_comp_fem.pdf', 'width', 'wide', 'height', 'normal'); + exportFig('figs/apa95ml_5kg_pi_comp_fem.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:apa95ml_5kg_pi_comp_fem @@ -959,18 +977,14 @@ Comparison with the FEM model #+RESULTS: [[file:figs/apa95ml_5kg_pi_comp_fem.png]] -* Identification of the dynamics from actuator to force sensor +* Identification of the dynamics from actuator Voltage to force sensor Voltage :PROPERTIES: :header-args:matlab+: :tangle matlab/force_sensor_identification.m :header-args:matlab+: :comments org :mkdirp yes :END: <> ** Introduction :ignore: -Two measurements are performed: -- Speedgoat DAC => Voltage Amplifier (x20) => 1 Piezo Stack => ... => 2 Stacks as Force Sensor (parallel) => Speedgoat ADC -- Speedgoat DAC => Voltage Amplifier (x20) => 2 Piezo Stacks (parallel) => ... => 1 Stack as Force Sensor => Speedgoat ADC - -The obtained dynamics from force actuator to force sensor are compare with the FEM model. +The same setup shown in Figure [[fig:test_bench_apa_identification]] is used. ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) @@ -995,51 +1009,56 @@ The data are loaded: load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'v'); #+end_src +Any offset is removed. #+begin_src matlab - u = detrend(u, 0); - v = detrend(v, 0); + u = detrend(u, 0); % Speedgoat DAC output Voltage [V] + v = detrend(v, 0); % Speedgoat ADC input Voltage (sensor stack) [V] #+end_src +Here, the amplifier gain is 20. #+begin_src matlab - u = 20*u; -#+end_src - -** Adjust gain :ignore: -Let's use the amplifier gain to obtain the true voltage applied to the actuator stack(s) - -The parameters of the piezoelectric stacks are defined below: -#+begin_src matlab - d33 = 3e-10; % Strain constant [m/V] - n = 80; % Number of layers per stack - eT = 1.6e-8; % Permittivity under constant stress [F/m] - sD = 2e-11; % Elastic compliance under constant electric displacement [m2/N] - ka = 235e6; % Stack stiffness [N/m] -#+end_src - -From the FEM, we construct the transfer function from DAC voltage to ADC voltage. -#+begin_src matlab - Gfem_aa_s = exp(-s/1e4)*20*(2*d33*n*ka)*(G(3,1)+G(3,2))*d33/(eT*sD*n); - Gfem_a_ss = exp(-s/1e4)*20*( d33*n*ka)*(G(3,1)+G(2,1))*d33/(eT*sD*n); -#+end_src - -#+begin_src matlab - Gfem_aa_s = exp(-s/1e4)*20*(2*d33*n*ka)*Gfm*d33/(eT*sD*n); - Gfem_a_ss = exp(-s/1e4)*20*( d33*n*ka)*Gfm*d33/(eT*sD*n); + u = 20*u; % Actuator Stack Voltage [V] #+end_src ** Compute TF estimate and Coherence :ignore: -The transfer function from input voltage to output voltage are computed and shown in Figure [[fig:bode_plot_force_sensor_voltage_comp_fem]]. +The transfer function from the actuator voltage $V_a$ to the force sensor stack voltage $V_s$ is computed. #+begin_src matlab Ts = t(end)/(length(t)-1); Fs = 1/Ts; - win = hann(ceil(10/Ts)); + win = hann(ceil(5/Ts)); [tf_est, f] = tfestimate(u, v, win, [], [], 1/Ts); [coh, ~] = mscohere( u, v, win, [], [], 1/Ts); #+end_src +The coherence is shown in Figure [[fig:apa95ml_5kg_cedrat_coh]]. + +#+begin_src matlab :exports none + figure; + + hold on; + plot(f, coh, 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Coherence'); xlabel('Frequency [Hz]'); + hold off; + xlim([10, 5e3]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/apa95ml_5kg_cedrat_coh.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:apa95ml_5kg_cedrat_coh +#+caption: Coherence +#+RESULTS: +[[file:figs/apa95ml_5kg_cedrat_coh.png]] + +The Simscape model is loaded and compared with the identified dynamics in Figure [[fig:bode_plot_force_sensor_voltage_comp_fem]]. +The non-minimum phase zero is just a side effect of the not so great identification. +Taking longer measurements would results in a minimum phase zero. + #+begin_src matlab load('mat/fem_simscape_models.mat', 'Gfm'); #+end_src @@ -1055,7 +1074,7 @@ The transfer function from input voltage to output voltage are computed and show plot(f, abs(tf_est)) plot(freqs, abs(squeeze(freqresp(Gfm, freqs, 'Hz')))) set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); - ylabel('Amplitude'); set(gca, 'XTickLabel',[]); + ylabel('Amplitude [V/V]'); set(gca, 'XTickLabel',[]); hold off; ylim([1e-3, 1e1]); @@ -1085,7 +1104,40 @@ The transfer function from input voltage to output voltage are computed and show #+RESULTS: [[file:figs/bode_plot_force_sensor_voltage_comp_fem.png]] -** System Identification +* Integral Force Feedback +:PROPERTIES: +:header-args:matlab+: :tangle matlab/integral_force_feedback.m +:header-args:matlab+: :comments org :mkdirp yes +:END: +<> +** Introduction :ignore: + +The test bench used to try the Integral Force Feedback control architecture is shown in Figure [[fig:test_bench_apa_schematic_iff]]. + +#+name: fig:test_bench_apa_schematic_iff +#+caption: Schematic of the test bench using IFF +[[file:figs/test_bench_apa_schematic_iff.png]] + +** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) + <> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+end_src + +#+begin_src matlab :tangle no + addpath('./matlab/mat/'); +#+end_src + +#+begin_src matlab :eval no + addpath('./mat/'); +#+end_src + +** IFF Plant +From the identified plant, a model of the transfer function from the actuator stack voltage to the force sensor generated voltage is developed. + #+begin_src matlab w_z = 2*pi*111; % Zeros frequency [rad/s] w_p = 2*pi*255; % Pole frequency [rad/s] @@ -1093,9 +1145,11 @@ The transfer function from input voltage to output voltage are computed and show xi_p = 0.015; G_inf = 0.1; - Gi = G_inf*(s^2 - 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2); + Gi = G_inf*(s^2 + 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2); #+end_src +Its bode plot is shown in Figure [[fig:iff_plant_identification_apa95ml]]. + #+begin_src matlab :exports none freqs = logspace(1, 4, 1000); @@ -1104,8 +1158,7 @@ The transfer function from input voltage to output voltage are computed and show ax1 = nexttile([2,1]); hold on; - plot(f, abs(tf_est), '-') - plot(freqs, abs(squeeze(freqresp(Gi, freqs, 'Hz'))), '--') + plot(freqs, abs(squeeze(freqresp(Gi, freqs, 'Hz'))), 'k-') set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); ylabel('Amplitude'); set(gca, 'XTickLabel', []); hold off; @@ -1113,14 +1166,12 @@ The transfer function from input voltage to output voltage are computed and show ax2 = nexttile; hold on; - plot(f, 180/pi*angle(tf_est), '-', 'DisplayName', '2 Act - 1 Sen') - plot(freqs, 180/pi*angle(squeeze(freqresp(Gi, freqs, 'Hz'))), '--', 'DisplayName', '2 Act - 1 Sen, - FEM') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gi, freqs, 'Hz'))), 'k-') set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); ylabel('Phase'); xlabel('Frequency [Hz]'); hold off; ylim([-180, 180]); yticks(-180:90:180); - legend('location', 'northeast') linkaxes([ax1,ax2], 'x'); xlim([10, 5e3]); @@ -1131,12 +1182,22 @@ The transfer function from input voltage to output voltage are computed and show #+end_src #+name: fig:iff_plant_identification_apa95ml -#+caption: Identification of the IFF plant +#+caption: Bode plot of the IFF plant #+RESULTS: [[file:figs/iff_plant_identification_apa95ml.png]] +The controller used in the Integral Force Feedback Architecture is: +\begin{equation} + K_{\text{IFF}}(s) = \frac{g}{s + 2\cdot 2\pi} \cdot \frac{s}{s + 0.5 \cdot 2\pi} +\end{equation} +where $g$ is a gain that can be tuned. + +Above 2 Hz the controller is basically an integrator, whereas an high pass filter is added at 0.5Hz to further reduce the low frequency gain. + +In the frequency band of interest, this controller should mostly act as a pure integrator. + +The Root Locus corresponding to this controller is shown in Figure [[fig:root_locus_iff_apa95ml_identification]]. -** Integral Force Feedback #+begin_src matlab :exports none gains = logspace(0, 5, 1000); @@ -1164,37 +1225,11 @@ The transfer function from input voltage to output voltage are computed and show #+RESULTS: [[file:figs/root_locus_iff_apa95ml_identification.png]] -* Integral Force Feedback -:PROPERTIES: -:header-args:matlab+: :tangle matlab/integral_force_feedback.m -:header-args:matlab+: :comments org :mkdirp yes -:END: -<> - -** Introduction :ignore: - -#+name: fig:test_bench_apa_schematic_iff -#+caption: Schematic of the test bench using IFF -[[file:figs/test_bench_apa_schematic_iff.png]] - -** Matlab Init :noexport:ignore: -#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> -#+end_src - -#+begin_src matlab :exports none :results silent :noweb yes - <> -#+end_src - -#+begin_src matlab :tangle no - addpath('./matlab/mat/'); -#+end_src - -#+begin_src matlab :eval no - addpath('./mat/'); -#+end_src ** First tests with few gains +The controller is now implemented in practice, and few controller gains are tested: $g = 0$, $g = 10$ and $g = 100$. + +For each controller gain, the identification shown in Figure [[fig:test_bench_apa_schematic_iff]] is performed. #+begin_src matlab iff_g10 = load('apa95ml_iff_g10_res.mat', 'u', 't', 'y', 'v'); iff_g100 = load('apa95ml_iff_g100_res.mat', 'u', 't', 'y', 'v'); @@ -1206,15 +1241,17 @@ The transfer function from input voltage to output voltage are computed and show win = hann(ceil(10/Ts)); [tf_iff_g10, f] = tfestimate(iff_g10.u, iff_g10.y, win, [], [], 1/Ts); - [co_iff_g10, ~] = mscohere(iff_g10.u, iff_g10.y, win, [], [], 1/Ts); + [co_iff_g10, ~] = mscohere( iff_g10.u, iff_g10.y, win, [], [], 1/Ts); [tf_iff_g100, ~] = tfestimate(iff_g100.u, iff_g100.y, win, [], [], 1/Ts); - [co_iff_g100, ~] = mscohere(iff_g100.u, iff_g100.y, win, [], [], 1/Ts); + [co_iff_g100, ~] = mscohere( iff_g100.u, iff_g100.y, win, [], [], 1/Ts); [tf_iff_of, ~] = tfestimate(iff_of.u, iff_of.y, win, [], [], 1/Ts); - [co_iff_of, ~] = mscohere(iff_of.u, iff_of.y, win, [], [], 1/Ts); + [co_iff_of, ~] = mscohere( iff_of.u, iff_of.y, win, [], [], 1/Ts); #+end_src +The coherence between the excitation signal and the mass displacement as measured by the interferometer is shown in Figure [[fig:iff_first_test_coherence]]. + #+begin_src matlab :exports none figure; @@ -1238,6 +1275,8 @@ The transfer function from input voltage to output voltage are computed and show #+RESULTS: [[file:figs/iff_first_test_coherence.png]] +The obtained transfer functions are shown in Figure [[fig:iff_first_test_bode_plot]]. +It is clear that the IFF architecture can actively damp the main resonance of the system. #+begin_src matlab :exports none figure; @@ -1277,6 +1316,7 @@ The transfer function from input voltage to output voltage are computed and show [[file:figs/iff_first_test_bode_plot.png]] ** Second test with many Gains +Then, the same identification test is performed for many more gains. #+begin_src matlab load('apa95ml_iff_test.mat', 'results'); #+end_src @@ -1300,19 +1340,7 @@ The transfer function from input voltage to output voltage are computed and show end #+end_src -#+begin_src matlab :exports none - figure; - - hold on; - for i = 1:length(results) - plot(f, co_iff{i}, '-', 'DisplayName', sprintf('g = %0.f', g_iff(i))) - end - set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); - ylabel('Coherence'); xlabel('Frequency [Hz]'); - hold off; - legend(); - xlim([60, 600]) -#+end_src +The obtained dynamics are shown in Figure [[fig:iff_results_bode_plots]]. #+begin_src matlab :exports none figure; @@ -1347,10 +1375,12 @@ The transfer function from input voltage to output voltage are computed and show #+end_src #+name: fig:iff_results_bode_plots -#+caption: +#+caption: Identified dynamics from excitation voltage to the mass displacement #+RESULTS: [[file:figs/iff_results_bode_plots.png]] +For each gain, the parameters of a second order resonant system that best fits the data are estimated and are compared with the data in Figure [[fig:iff_results_bode_plots_identification]]. + #+begin_src matlab G_id = {zeros(1,length(results))}; @@ -1406,10 +1436,12 @@ The transfer function from input voltage to output voltage are computed and show #+end_src #+name: fig:iff_results_bode_plots_identification -#+caption: +#+caption: Comparison of the measured dynamic and the identified 2nd order resonant systems that best fits the data #+RESULTS: [[file:figs/iff_results_bode_plots_identification.png]] +Finally, we can represent the position of the poles of the 2nd order systems on the Root Locus plot (Figure [[fig:iff_results_root_locus]]). + #+begin_src matlab :exports none w_z = 2*pi*111; % Zeros frequency [rad/s] w_p = 2*pi*255; % Pole frequency [rad/s] @@ -1447,6 +1479,6 @@ The transfer function from input voltage to output voltage are computed and show #+end_src #+name: fig:iff_results_root_locus -#+caption: +#+caption: Root Locus plot of the identified IFF plant as well as the identified poles of the damped system #+RESULTS: [[file:figs/iff_results_root_locus.png]] diff --git a/matlab/force_sensor_identification.m b/matlab/force_sensor_identification.m new file mode 100644 index 0000000..e7f7d35 --- /dev/null +++ b/matlab/force_sensor_identification.m @@ -0,0 +1,95 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +addpath('./mat/'); + +% Load Data :ignore: +% The data are loaded: + +load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'v'); + + + +% Any offset is removed. + +u = detrend(u, 0); % Speedgoat DAC output Voltage [V] +v = detrend(v, 0); % Speedgoat ADC input Voltage (sensor stack) [V] + + + +% Here, the amplifier gain is 20. + +u = 20*u; % Actuator Stack Voltage [V] + +% Compute TF estimate and Coherence :ignore: +% The transfer function from the actuator voltage $V_a$ to the force sensor stack voltage $V_s$ is computed. + + +Ts = t(end)/(length(t)-1); +Fs = 1/Ts; + +win = hann(ceil(5/Ts)); + +[tf_est, f] = tfestimate(u, v, win, [], [], 1/Ts); +[coh, ~] = mscohere( u, v, win, [], [], 1/Ts); + + + +% The coherence is shown in Figure [[fig:apa95ml_5kg_cedrat_coh]]. + + +figure; + +hold on; +plot(f, coh, 'k-') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Coherence'); xlabel('Frequency [Hz]'); +hold off; +xlim([10, 5e3]); + + + +% #+name: fig:apa95ml_5kg_cedrat_coh +% #+caption: Coherence +% #+RESULTS: +% [[file:figs/apa95ml_5kg_cedrat_coh.png]] + +% The Simscape model is loaded and compared with the identified dynamics in Figure [[fig:bode_plot_force_sensor_voltage_comp_fem]]. +% The non-minimum phase zero is just a side effect of the not so great identification. +% Taking longer measurements would results in a minimum phase zero. + + +load('mat/fem_simscape_models.mat', 'Gfm'); + +freqs = logspace(1, 4, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(f, abs(tf_est)) +plot(freqs, abs(squeeze(freqresp(Gfm, freqs, 'Hz')))) +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude [V/V]'); set(gca, 'XTickLabel',[]); +hold off; +ylim([1e-3, 1e1]); + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(tf_est), ... + 'DisplayName', 'Identification') +plot(freqs, 180/pi*angle(squeeze(freqresp(Gfm, freqs, 'Hz'))), ... + 'DisplayName', 'FEM') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase'); xlabel('Frequency [Hz]'); +hold off; +ylim([-180, 180]); +yticks(-180:90:180); +legend('location', 'northeast') + +linkaxes([ax1,ax2], 'x'); +xlim([10, 5e3]); diff --git a/matlab/huddle_test.m b/matlab/huddle_test.m new file mode 100644 index 0000000..92e358a --- /dev/null +++ b/matlab/huddle_test.m @@ -0,0 +1,34 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +addpath('./mat/'); + +% Load Data :noexport: + +load('huddle_test.mat', 't', 'y'); + +y = y - mean(y(1000:end)); + +% Time Domain Data + +figure; +plot(t(1000:end), y(1000:end)) +ylabel('Output Displacement [m]'); xlabel('Time [s]'); + +% PSD of Measurement Noise + +Ts = t(end)/(length(t)-1); +Fs = 1/Ts; + +win = hanning(ceil(1*Fs)); + +[pxx, f] = pwelch(y(1000:end), win, [], [], Fs); + +figure; +plot(f, sqrt(pxx)); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]'); +xlim([1, Fs/2]); diff --git a/matlab/integral_force_feedback.m b/matlab/integral_force_feedback.m new file mode 100644 index 0000000..fe1f802 --- /dev/null +++ b/matlab/integral_force_feedback.m @@ -0,0 +1,307 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +addpath('./mat/'); + +% IFF Plant +% From the identified plant, a model of the transfer function from the actuator stack voltage to the force sensor generated voltage is developed. + + +w_z = 2*pi*111; % Zeros frequency [rad/s] +w_p = 2*pi*255; % Pole frequency [rad/s] +xi_z = 0.05; +xi_p = 0.015; +G_inf = 0.1; + +Gi = G_inf*(s^2 + 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2); + + + +% Its bode plot is shown in Figure [[fig:iff_plant_identification_apa95ml]]. + + +freqs = logspace(1, 4, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(freqs, abs(squeeze(freqresp(Gi, freqs, 'Hz'))), 'k-') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude'); set(gca, 'XTickLabel', []); +hold off; +ylim([1e-3, 1e1]); + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*angle(squeeze(freqresp(Gi, freqs, 'Hz'))), 'k-') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase'); xlabel('Frequency [Hz]'); +hold off; +ylim([-180, 180]); +yticks(-180:90:180); + +linkaxes([ax1,ax2], 'x'); +xlim([10, 5e3]); + + + +% #+name: fig:iff_plant_identification_apa95ml +% #+caption: Bode plot of the IFF plant +% #+RESULTS: +% [[file:figs/iff_plant_identification_apa95ml.png]] + +% The controller used in the Integral Force Feedback Architecture is: +% \begin{equation} +% K_{\text{IFF}}(s) = \frac{g}{s + 2\cdot 2\pi} \cdot \frac{s}{s + 0.5 \cdot 2\pi} +% \end{equation} +% where $g$ is a gain that can be tuned. + +% Above 2 Hz the controller is basically an integrator, whereas an high pass filter is added at 0.5Hz to further reduce the low frequency gain. + +% In the frequency band of interest, this controller should mostly act as a pure integrator. + +% The Root Locus corresponding to this controller is shown in Figure [[fig:root_locus_iff_apa95ml_identification]]. + + +gains = logspace(0, 5, 1000); + +figure; +hold on; +plot(real(pole(Gi)), imag(pole(Gi)), 'kx'); +plot(real(tzero(Gi)), imag(tzero(Gi)), 'ko'); +for i = 1:length(gains) + cl_poles = pole(feedback(Gi, (gains(i)/(s + 2*2*pi)*s/(s + 0.5*2*pi)))); + plot(real(cl_poles), imag(cl_poles), 'k.'); +end +ylim([0, 1800]); +xlim([-1600,200]); +xlabel('Real Part') +ylabel('Imaginary Part') +axis square + +% First tests with few gains +% The controller is now implemented in practice, and few controller gains are tested: $g = 0$, $g = 10$ and $g = 100$. + +% For each controller gain, the identification shown in Figure [[fig:test_bench_apa_schematic_iff]] is performed. + +iff_g10 = load('apa95ml_iff_g10_res.mat', 'u', 't', 'y', 'v'); +iff_g100 = load('apa95ml_iff_g100_res.mat', 'u', 't', 'y', 'v'); +iff_of = load('apa95ml_iff_off_res.mat', 'u', 't', 'y', 'v'); + +Ts = 1e-4; +win = hann(ceil(10/Ts)); + +[tf_iff_g10, f] = tfestimate(iff_g10.u, iff_g10.y, win, [], [], 1/Ts); +[co_iff_g10, ~] = mscohere( iff_g10.u, iff_g10.y, win, [], [], 1/Ts); + +[tf_iff_g100, ~] = tfestimate(iff_g100.u, iff_g100.y, win, [], [], 1/Ts); +[co_iff_g100, ~] = mscohere( iff_g100.u, iff_g100.y, win, [], [], 1/Ts); + +[tf_iff_of, ~] = tfestimate(iff_of.u, iff_of.y, win, [], [], 1/Ts); +[co_iff_of, ~] = mscohere( iff_of.u, iff_of.y, win, [], [], 1/Ts); + + + +% The coherence between the excitation signal and the mass displacement as measured by the interferometer is shown in Figure [[fig:iff_first_test_coherence]]. + + +figure; + +hold on; +plot(f, co_iff_of, '-', 'DisplayName', 'g=0') +plot(f, co_iff_g10, '-', 'DisplayName', 'g=10') +plot(f, co_iff_g100, '-', 'DisplayName', 'g=100') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Coherence'); xlabel('Frequency [Hz]'); +hold off; +legend(); +xlim([60, 600]) + + + +% #+name: fig:iff_first_test_coherence +% #+caption: Coherence +% #+RESULTS: +% [[file:figs/iff_first_test_coherence.png]] + +% The obtained transfer functions are shown in Figure [[fig:iff_first_test_bode_plot]]. +% It is clear that the IFF architecture can actively damp the main resonance of the system. + + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2, 1]); +hold on; +plot(f, abs(tf_iff_of), '-', 'DisplayName', 'g=0') +plot(f, abs(tf_iff_g10), '-', 'DisplayName', 'g=10') +plot(f, abs(tf_iff_g100), '-', 'DisplayName', 'g=100') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude'); set(gca, 'XTickLabel',[]); +hold off; +legend(); + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(-tf_iff_of), '-') +plot(f, 180/pi*angle(-tf_iff_g10), '-') +plot(f, 180/pi*angle(-tf_iff_g100), '-') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase'); xlabel('Frequency [Hz]'); +hold off; +yticks(-360:90:360); + +linkaxes([ax1,ax2], 'x'); +xlim([60, 600]); + +% Second test with many Gains +% Then, the same identification test is performed for many more gains. + +load('apa95ml_iff_test.mat', 'results'); + +Ts = 1e-4; +win = hann(ceil(10/Ts)); + +tf_iff = {zeros(1, length(results))}; +co_iff = {zeros(1, length(results))}; +g_iff = [0, 1, 5, 10, 50, 100]; + +for i=1:length(results) + [tf_est, f] = tfestimate(results{i}.u, results{i}.y, win, [], [], 1/Ts); + [co_est, ~] = mscohere(results{i}.u, results{i}.y, win, [], [], 1/Ts); + + tf_iff(i) = {tf_est}; + co_iff(i) = {co_est}; +end + + + +% The obtained dynamics are shown in Figure [[fig:iff_results_bode_plots]]. + + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2, 1]); +hold on; +for i = 1:length(results) + plot(f, abs(tf_iff{i}), '-', 'DisplayName', sprintf('g = %0.f', g_iff(i))) +end +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude'); set(gca, 'XTickLabel',[]); +hold off; +legend(); + +ax2 = nexttile; +hold on; +for i = 1:length(results) + plot(f, 180/pi*angle(-tf_iff{i}), '-') +end +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase'); xlabel('Frequency [Hz]'); +yticks(-360:90:360); +axis padded 'auto x' + +linkaxes([ax1,ax2], 'x'); +xlim([60, 600]); + + + +% #+name: fig:iff_results_bode_plots +% #+caption: Identified dynamics from excitation voltage to the mass displacement +% #+RESULTS: +% [[file:figs/iff_results_bode_plots.png]] + +% For each gain, the parameters of a second order resonant system that best fits the data are estimated and are compared with the data in Figure [[fig:iff_results_bode_plots_identification]]. + + +G_id = {zeros(1,length(results))}; + +f_start = 70; % [Hz] +f_end = 500; % [Hz] + +for i = 1:length(results) + tf_id = tf_iff{i}(sum(ff_end)); + f_id = f(sum(ff_end)); + + gfr = idfrd(tf_id, 2*pi*f_id, Ts); + G_id(i) = {procest(gfr,'P2UDZ')}; +end + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2, 1]); +hold on; +for i = 1:length(results) + set(gca,'ColorOrderIndex',i) + plot(f, abs(tf_iff{i}), '-', 'DisplayName', sprintf('g = %0.f', g_iff(i))) + set(gca,'ColorOrderIndex',i) + plot(f, abs(squeeze(freqresp(G_id{i}, f, 'Hz'))), '--', 'HandleVisibility', 'off') +end +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude'); set(gca, 'XTickLabel',[]); +hold off; +legend(); + +ax2 = nexttile; +hold on; +for i = 1:length(results) + set(gca,'ColorOrderIndex',i) + plot(f, 180/pi*angle(tf_iff{i}), '-') + set(gca,'ColorOrderIndex',i) + plot(f, 180/pi*angle(squeeze(freqresp(G_id{i}, f, 'Hz'))), '--') +end +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([60, 600]); + + + +% #+name: fig:iff_results_bode_plots_identification +% #+caption: Comparison of the measured dynamic and the identified 2nd order resonant systems that best fits the data +% #+RESULTS: +% [[file:figs/iff_results_bode_plots_identification.png]] + +% Finally, we can represent the position of the poles of the 2nd order systems on the Root Locus plot (Figure [[fig:iff_results_root_locus]]). + + +w_z = 2*pi*111; % Zeros frequency [rad/s] +w_p = 2*pi*255; % Pole frequency [rad/s] +xi_z = 0.05; +xi_p = 0.015; +G_inf = 2; + +Gi = G_inf*(s^2 - 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2); + + +gains = logspace(0, 5, 1000); + +figure; +hold on; +plot(real(pole(Gi)), imag(pole(Gi)), 'kx', 'HandleVisibility', 'off'); +plot(real(tzero(Gi)), imag(tzero(Gi)), 'ko', 'HandleVisibility', 'off'); +for i = 1:length(results) + set(gca,'ColorOrderIndex',i) + plot(real(pole(G_id{i})), imag(pole(G_id{i})), 'o', 'DisplayName', sprintf('g = %0.f', g_iff(i))); +end +for i = 1:length(gains) + cl_poles = pole(feedback(Gi, (gains(i)/(s + 2*pi*2)))); + plot(real(cl_poles), imag(cl_poles), 'k.', 'HandleVisibility', 'off'); +end +ylim([0, 1800]); +xlim([-1600,200]); +xlabel('Real Part') +ylabel('Imaginary Part') +axis square +legend('location', 'northwest'); diff --git a/matlab/mat/fem_simscape_models.mat b/matlab/mat/fem_simscape_models.mat index 47fd28a..e5badb3 100644 Binary files a/matlab/mat/fem_simscape_models.mat and b/matlab/mat/fem_simscape_models.mat differ diff --git a/matlab/motion_identification.m b/matlab/motion_identification.m new file mode 100644 index 0000000..3d89dd4 --- /dev/null +++ b/matlab/motion_identification.m @@ -0,0 +1,100 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +addpath('./mat/'); + +% Load Data +% The data from the "noise test" and the identification test are loaded. + +ht = load('huddle_test.mat', 't', 'y'); +load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y'); + + + +% Any offset value is removed: + +um = detrend(um, 0); % Generated DAC Voltage [V] +y = detrend(y , 0); % Mass displacement [m] + +ht.y = detrend(ht.y, 0); + + + +% Now we add a factor 10 to take into account the gain of the voltage amplifier. + +Va = 10*um; + +% Comparison of the PSD with Huddle Test + +Ts = t(end)/(length(t)-1); +Fs = 1/Ts; + +win = hanning(ceil(1*Fs)); + +[pxx, f] = pwelch(y, win, [], [], Fs); +[pht, ~] = pwelch(ht.y, win, [], [], Fs); + +figure; +hold on; +plot(f, sqrt(pxx), 'DisplayName', '5kg'); +plot(f, sqrt(pht), 'DisplayName', 'Huddle Test'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]'); +legend('location', 'northeast'); +xlim([1, Fs/2]); + +% Compute TF estimate and Coherence + +[tf_est, f] = tfestimate(Va, -y, win, [], [], 1/Ts); +[co_est, ~] = mscohere( Va, -y, win, [], [], 1/Ts); + +figure; + +hold on; +plot(f, co_est, 'k-') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Coherence'); xlabel('Frequency [Hz]'); +hold off; +xlim([10, 5e3]); + + + +% #+name: fig:apa95ml_5kg_PI_coh +% #+caption: Coherence +% #+RESULTS: +% [[file:figs/apa95ml_5kg_PI_coh.png]] + +% Comparison with the FEM model + +load('mat/fem_simscape_models.mat', 'Ghm'); + +freqs = logspace(0, 4, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(f, abs(tf_est), 'DisplayName', 'Identification') +plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), 'DisplayName', 'FEM') +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); +ylabel('Amplitude [m/V]'); set(gca, 'XTickLabel', []); +legend('location', 'northeast') +hold off; + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(tf_est)) +plot(freqs, 180/pi*angle(squeeze(freqresp(Ghm, freqs, 'Hz')))) +set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); +ylabel('Phase'); xlabel('Frequency [Hz]'); +hold off; +ylim([-180, 180]); +yticks(-180:90:180); + +linkaxes([ax1,ax2], 'x'); +xlim([10, 5e3]); diff --git a/matlab/piezo_amplified_3d.slx b/matlab/piezo_amplified_3d.slx index e1e0709..15f066f 100644 Binary files a/matlab/piezo_amplified_3d.slx and b/matlab/piezo_amplified_3d.slx differ diff --git a/matlab/piezo_parameters.m b/matlab/piezo_parameters.m index 348c4e4..0c19cf8 100644 --- a/matlab/piezo_parameters.m +++ b/matlab/piezo_parameters.m @@ -38,7 +38,7 @@ ka = 235e6; % [N/m] Lmax = 20e-6; % [m] Vmax = 170; % [V] -ka*Lmax/Vmax % [N/V] +ga = ka*Lmax/Vmax; % [N/V] % Estimation from Piezoelectric parameters % <> diff --git a/matlab/simscape_model.m b/matlab/simscape_model.m new file mode 100644 index 0000000..a90b8ec --- /dev/null +++ b/matlab/simscape_model.m @@ -0,0 +1,120 @@ +%% Clear Workspace and Close figures +clear; close all; clc; + +%% Intialize Laplace variable +s = zpk('s'); + +addpath('./mat/'); + +% Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates +% We first extract the stiffness and mass matrices. + +K = readmatrix('APA95ML_K.CSV'); +M = readmatrix('APA95ML_M.CSV'); + + + +% #+caption: First 10x10 elements of the Mass matrix +% #+RESULTS: +% | 0.03 | 7e-08 | 2e-06 | -3e-09 | -0.0002 | -6e-08 | -0.001 | 8e-07 | 6e-07 | -8e-09 | +% | 7e-08 | 0.02 | -1e-06 | 9e-05 | -3e-09 | -4e-09 | -1e-06 | -0.0006 | -4e-08 | 5e-06 | +% | 2e-06 | -1e-06 | 0.02 | -3e-08 | -4e-08 | 1e-08 | 1e-07 | -2e-07 | 0.0003 | 1e-09 | +% | -3e-09 | 9e-05 | -3e-08 | 1e-06 | -3e-11 | -3e-13 | -7e-09 | -5e-06 | -3e-10 | 3e-08 | +% | -0.0002 | -3e-09 | -4e-08 | -3e-11 | 2e-06 | 6e-10 | 2e-06 | -7e-09 | -2e-09 | 7e-11 | +% | -6e-08 | -4e-09 | 1e-08 | -3e-13 | 6e-10 | 1e-06 | 1e-08 | 3e-09 | -2e-09 | 2e-13 | +% | -0.001 | -1e-06 | 1e-07 | -7e-09 | 2e-06 | 1e-08 | 0.03 | 4e-08 | -2e-06 | 8e-09 | +% | 8e-07 | -0.0006 | -2e-07 | -5e-06 | -7e-09 | 3e-09 | 4e-08 | 0.02 | -9e-07 | -9e-05 | +% | 6e-07 | -4e-08 | 0.0003 | -3e-10 | -2e-09 | -2e-09 | -2e-06 | -9e-07 | 0.02 | 2e-08 | +% | -8e-09 | 5e-06 | 1e-09 | 3e-08 | 7e-11 | 2e-13 | 8e-09 | -9e-05 | 2e-08 | 1e-06 | + +% Then, we extract the coordinates of the interface nodes. + +[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA95ML_out_nodes_3D.txt'); + +% Simscape Model +% The flexible element is imported using the =Reduced Order Flexible Solid= Simscape block. + +% To model the actuator, an =Internal Force= block is added between the nodes 3 and 12. +% A =Relative Motion Sensor= block is added between the nodes 1 and 2 to measure the displacement and the amplified piezo. + +% One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame. + +m = 5.5; + +load('apa95ml_params.mat', 'ga', 'gs'); + +% Dynamics from Actuator Voltage to Vertical Mass Displacement +% The identified dynamics is shown in Figure [[fig:dynamics_act_disp_comp_mass]]. + + +%% Name of the Simulink File +mdl = 'piezo_amplified_3d'; + +%% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % Actuator Voltage [V] +io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m] + +Ghm = linearize(mdl, io); + +freqs = logspace(0, 4, 5000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz')))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); +hold off; + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz'))))); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +yticks(-360:90:360); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); + +% Dynamics from Actuator Voltage to Force Sensor Voltage +% The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_mass]]. + + +%% Name of the Simulink File +mdl = 'piezo_amplified_3d'; + +%% Input/Output definition +clear io; io_i = 1; +io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % Voltage Actuator [V] +io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; % Sensor Voltage [V] + +Gfm = linearize(mdl, io); + +freqs = logspace(1, 5, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(freqs, abs(squeeze(freqresp(Gfm, freqs, 'Hz')))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); +hold off; + +ax2 = nexttile; +hold on; +plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gfm, freqs, 'Hz'))))); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +yticks(-360:90:360); ylim([-360, 180]); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); + +save('mat/fem_simscape_models.mat', 'Ghm', 'Gfm')