Active Damping of Rotating Platforms using Integral Force Feedback - Matlab Computation
Table of ContentsClose
This report is also available as a pdf.
This document gathers the Matlab code used to for the conference paper (Dehaeze and Collette 2020) and the journal paper (Dehaeze and Collette 2021).
It is structured in several sections:
- Section 1: presents a simple model of a rotating suspended platform that will be used throughout this study.
- Section 2: explains how the unconditional stability of IFF is lost due to Gyroscopic effects induced by the rotation.
- Section 3: suggests a simple modification of the control law such that damping can be added to the suspension modes in a robust way.
- Section 4: proposes to add springs in parallel with the force sensors to regain the unconditional stability of IFF.
- Section 5: compares both proposed modifications to the classical IFF in terms of damping authority and closed-loop system behavior.
- Section 6: contains the notations used for both the Matlab code and the paper
The matlab code is accessible on Zonodo and Github (Dehaeze 2020). It can also be download as a .zip
file here.
To run the Matlab code, go in the matlab
directory and run the following Matlab files corresponding to each section.
Sections | Matlab File |
---|---|
Section 1 | s1_system_description.m |
Section 2 | s2_iff_pure_int.m |
Section 3 | s3_iff_hpf.m |
Section 4 | s4_iff_kp.m |
Section 5 | s5_act_damp_comparison.m |
1 System Description and Analysis
1.1 System description
The system consists of one 2 degree of freedom translation stage on top of a spindle (figure 1).
Figure 1: Schematic of the studied system
The control inputs are the forces applied by the actuators of the translation stage (
1.2 Equations
Based on the Figure 1, the equations of motions are:
Important
Where
With:
1.3 Numerical Values
Let’s define initial values for the model.
Copyk = 1; % Actuator Stiffness [N/m] c = 0.05; % Actuator Damping [N/(m/s)] m = 1; % Payload mass [kg]
Copyxi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); % [rad/s]
1.4 Campbell Diagram
The Campbell Diagram displays the evolution of the real and imaginary parts of the system as a function of the rotating speed.
It is shown in Figures 2 and 3, and one can see that the system becomes unstable for
Figure 2: Campbell Diagram - Real Part
Figure 3: Campbell Diagram - Imaginary Part
1.5 Simscape Model
In order to validate all the equations of motion, a Simscape model of the same system has been developed. The dynamics of the system can be identified from the Simscape model and compare with the analytical model.
The rotating speed for the Simscape Model is defined.
CopyW = 0.1; % Rotation Speed [rad/s]
Copyopen('rotating_frame.slx');
The transfer function from
Copy%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/K'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/G'], 2, 'openoutput'); io_i = io_i + 1;
CopyG = linearize(mdl, io, 0); %% Input/Output definition G.InputName = {'Fu', 'Fv'}; G.OutputName = {'du', 'dv'};
The same transfer function from
CopyGth = (1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2), 2*W*s/(w0^2) ; ... -2*W*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)];
Both transfer functions are compared in Figure 4 and are found to perfectly match.
Figure 4: Bode plot of the transfer function from
1.6 Effect of the rotation speed
The transfer functions from
CopyWs = [0, 0.2, 0.7, 1.1]*w0; % Rotating Speeds [rad/s]
CopyGs = {zeros(2, 2, length(Ws))}; for W_i = 1:length(Ws) W = Ws(W_i); Gs(:, :, W_i) = {(1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2), 2*W*s/(w0^2) ; ... -2*W*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)]}; end
They are compared in Figures 5 and 6.
Figure 5: Comparison of the transfer functions from
Figure 6: Comparison of the transfer functions from
2 Problem with pure Integral Force Feedback
Force sensors are added in series with the two actuators (Figure 7).
Two identical controllers
Figure 7: System with added Force Sensor in series with the actuators
2.1 Plant Parameters
Let’s define initial values for the model.
Copyk = 1; % Actuator Stiffness [N/m] c = 0.05; % Actuator Damping [N/(m/s)] m = 1; % Payload mass [kg]
Copyxi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); % [rad/s]
2.2 Equations
The sensed forces are equal to:
Which then gives:
Important
2.3 Comparison of the Analytical Model and the Simscape Model
The rotation speed is set to
CopyW = 0.1*w0; % [rad/s]
Copyopen('rotating_frame.slx');
And the transfer function from
Copy%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/K'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/G'], 1, 'openoutput'); io_i = io_i + 1;
CopyGiff = linearize(mdl, io, 0); %% Input/Output definition Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'fu', 'fv'};
The same transfer function from
CopyGiff_th = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
The two are compared in Figure 8 and found to perfectly match.
Figure 8: Comparison of the transfer functions from
2.4 Effect of the rotation speed
The transfer functions from
CopyWs = [0, 0.2, 0.7]*w0; % Rotating Speeds [rad/s]
CopyGsiff = {zeros(2, 2, length(Ws))}; for W_i = 1:length(Ws) W = Ws(W_i); Gsiff(:, :, W_i) = {1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2]}; end
The obtained transfer functions are shown in Figure 9.
Figure 9: Comparison of the transfer functions from
2.5 Decentralized Integral Force Feedback
The decentralized IFF controller consists of pure integrators:
The Root Locus (evolution of the poles of the closed loop system in the complex plane as a function of
Figure 10: Root Locus for the Decentralized Integral Force Feedback controller. Several rotating speed are shown.
3 Integral Force Feedback with an High Pass Filter
3.1 Plant Parameters
Let’s define initial values for the model.
Copyk = 1; % Actuator Stiffness [N/m] c = 0.05; % Actuator Damping [N/(m/s)] m = 1; % Payload mass [kg]
Copyxi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); % [rad/s]
3.2 Modified Integral Force Feedback Controller
Let’s modify the initial Integral Force Feedback Controller ; instead of using pure integrators, pseudo integrators (i.e. low pass filters) are used:
where
Let’s arbitrary choose the following control parameters:
Copyg = 2; wi = 0.1*w0;
And the following rotating speed.
CopyGiff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
The obtained Loop Gain is shown in Figure 11.
Figure 11: Loop Gain for the modified IFF controller
3.3 Root Locus
As shown in the Root Locus plot (Figure 12), for some value of the gain, the system remains stable.
Figure 12: Root Locus for the modified IFF controller
Figure 13: Root Locus for the modified IFF controller - Zoom
3.4 What is the optimal and ?
In order to visualize the effect of
Copywis = [0.01, 0.1, 0.5, 1]*w0; % [rad/s]
Figure 14: Root Locus for the modified IFF controller (zoomed plot on the left)
Figure 15: Root Locus for the modified IFF controller (zoomed plot on the left)
For the controller
The gain at which the system becomes unstable is
While it seems that small
There must be an optimum for
Copywis = logspace(-2, 1, 100)*w0; % [rad/s] opt_xi = zeros(1, length(wis)); % Optimal simultaneous damping opt_gain = zeros(1, length(wis)); % Corresponding optimal gain for wi_i = 1:length(wis) wi = wis(wi_i); Kiff = 1/(s + wi)*eye(2); fun = @(g)computeSimultaneousDamping(g, Giff, Kiff); [g_opt, xi_opt] = fminsearch(fun, 0.5*wi*((w0/W)^2 - 1)); opt_xi(wi_i) = 1/xi_opt; opt_gain(wi_i) = g_opt; end
Figure 16: Simultaneous attainable damping of the closed loop poles as a function of
4 IFF with a stiffness in parallel with the force sensor
4.1 Schematic
In this section additional springs in parallel with the force sensors are added to counteract the negative stiffness induced by the rotation.
Figure 17: Studied system with additional springs in parallel with the actuators and force sensors
In order to keep the overall stiffness
4.2 Equations
Important
With:
If we compare
4.3 Plant Parameters
Let’s define initial values for the model.
Copyk = 1; % Actuator Stiffness [N/m] c = 0.05; % Actuator Damping [N/(m/s)] m = 1; % Payload mass [kg]
Copyxi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); % [rad/s]
4.4 Comparison of the Analytical Model and the Simscape Model
The same transfer function from
CopyW = 0.1*w0; % [rad/s] kp = 1.5*m*W^2; cp = 0;
Copyopen('rotating_frame.slx');
Copy%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/K'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/G'], 1, 'openoutput'); io_i = io_i + 1; Giff = linearize(mdl, io, 0); %% Input/Output definition Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'fu', 'fv'};
Copyw0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); Giff_th = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2 ]; Giff_th.InputName = {'Fu', 'Fv'}; Giff_th.OutputName = {'fu', 'fv'};
Figure 18: Comparison of the transfer functions from
4.5 Effect of the parallel stiffness on the IFF plant
The rotation speed is set to
CopyW = 0.1*w0; % [rad/s]
And the IFF plant (transfer function from
- without parallel stiffness
- with a small parallel stiffness
- with a large parallel stiffness
The results are shown in Figure 19.
One can see that for
Copykp = 0; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
Copykp = 0.5*m*W^2; k = 1 - kp; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); Giff_s = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
Copykp = 1.5*m*W^2; k = 1 - kp; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); Giff_l = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
Figure 19: Transfer function from
4.6 IFF when adding a spring in parallel
In Figure 20 is displayed the Root Locus in the three considered cases with
One can see that for
Thus, decentralized IFF controller with pure integrators can be used if:
Figure 20: Root Locus
Figure 21: Root Locus
4.7 Effect of on the attainable damping
However, having large values of
To study the second point, Root Locus plots for the following values of
Copykps = [2, 20, 40]*m*W^2;
It is shown that large values of
Figure 22: Root Locus plot
Copyalphas = logspace(-2, 0, 100); opt_xi = zeros(1, length(alphas)); % Optimal simultaneous damping opt_gain = zeros(1, length(alphas)); % Corresponding optimal gain Kiff = 1/s*eye(2); for alpha_i = 1:length(alphas) kp = alphas(alpha_i); k = 1 - alphas(alpha_i); w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2]; fun = @(g)computeSimultaneousDamping(g, Giff, Kiff); [g_opt, xi_opt] = fminsearch(fun, 2); opt_xi(alpha_i) = 1/xi_opt; opt_gain(alpha_i) = g_opt; end
Figure 23: Attainable damping ratio and corresponding controller gain for different parameter
5 Comparison
Two modifications to adapt the IFF control strategy to rotating platforms have been proposed. These two methods are now compared in terms of added damping, closed-loop compliance and transmissibility.
5.1 Plant Parameters
Let’s define initial values for the model.
Copyk = 1; % Actuator Stiffness [N/m] c = 0.05; % Actuator Damping [N/(m/s)] m = 1; % Payload mass [kg]
Copyxi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); % [rad/s]
The rotating speed is set to
CopyW = 0.1*w0;
5.2 Root Locus
IFF with High Pass Filter
Copywi = 0.1*w0; % [rad/s] Giff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
IFF With parallel Stiffness
Copykp = 5*m*W^2; k = k - kp; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); Giff_kp = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2 ]; k = k + kp;
Figure 24: Root Locus plot - Comparison of IFF with additional high pass filter, IFF with additional parallel stiffness
5.3 Controllers - Optimal Gains
In order to compare to three considered Active Damping techniques, gains that yield maximum damping of all the modes are computed for each case.
The obtained damping ratio and control are shown below.
Obtained |
Control Gain | |
---|---|---|
Modified IFF | 0.83 | 1.99 |
IFF with |
0.83 | 2.02 |
5.4 Passive Damping - Critical Damping
Critical Damping corresponds to to
Copyc_opt = 2*sqrt(k*m);
5.5 Transmissibility And Compliance
Copyopen('rotating_frame.slx');
Copy%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1; io(io_i) = linio([mdl, '/fd'], 1, 'input'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;
CopyG_ol = linearize(mdl, io, 0); %% Input/Output definition G_ol.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; G_ol.OutputName = {'Dx', 'Dy'};
5.5.1 Passive Damping
Copykp = 0; cp = 0;
Copyc_old = c; c = c_opt;
CopyG_pas = linearize(mdl, io, 0); %% Input/Output definition G_pas.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; G_pas.OutputName = {'Dx', 'Dy'};
Copyc = c_old;
CopyKiff = opt_gain_iff/(wi + s)*tf(eye(2));
CopyG_iff = linearize(mdl, io, 0); %% Input/Output definition G_iff.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; G_iff.OutputName = {'Dx', 'Dy'};
Copykp = 5*m*W^2; cp = 0.01;
CopyKiff = opt_gain_kp/s*tf(eye(2));
CopyG_kp = linearize(mdl, io, 0); %% Input/Output definition G_kp.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; G_kp.OutputName = {'Dx', 'Dy'};
Figure 25: Comparison of the transmissibility
Figure 26: Comparison of the obtained Compliance
6 Notations
Mathematical Notation | Matlab | Unit | |
---|---|---|---|
Actuator Stiffness | k |
N/m | |
Actuator Damping | c |
N/(m/s) | |
Payload Mass | m |
kg | |
Damping Ratio | xi |
||
Actuator Force | F Fu Fv |
N | |
Force Sensor signal | f fu fv |
N | |
Relative Displacement | d du dv |
m | |
Resonance freq. when |
w0 |
rad/s | |
Rotation Speed | W |
rad/s | |
Low Pass Filter corner frequency | wi |
rad/s |
Mathematical Notation | Matlab | Unit | |
---|---|---|---|
Laplace variable | s |
||
Complex number | j |
||
Frequency | w |
[rad/s] |