Matlab Computation
Table of Contents
- 1. System Description and Analysis
- 2. Integral Force Feedback
- 3. Direct Velocity Feedback
- 4. Comparison
Notations:
- \(k\): Actuator’s Stiffness [N/m]
- \(c\): Actuator’s Damping [N/(m/s)]
- \(m\): Payload’s mass [kg]
- \(\omega_0 = \sqrt{\frac{k}{m}}\): Resonance of the (non-rotating) mass-spring system [rad/s]
- \(\Omega = \dot{\theta}\): rotation speed [rad/s]
Mathematical Notation | Matlab | Unit | |
---|---|---|---|
Actuator Stiffness | \(k\) | k |
N/m |
Actuator Damping | \(c\) | c |
N/(m/s) |
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: Figure caption
The control inputs are the forces applied by the actuators of the translation stage (\(F_u\) and \(F_v\)). As the translation stage is rotating around the Z axis due to the spindle, the forces are applied along \(\vec{i}_u\) and \(\vec{i}_v\).
The measurement is either the \(x-y\) displacement of the object located on top of the translation stage or the \(u-v\) displacement of the sample with respect to a fixed reference frame.
1.2 Equations
Based on the Figure 1.
\begin{equation} \begin{bmatrix} d_u \\ d_v \end{bmatrix} = \frac{\frac{1}{k}}{\left( \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \right)^2 + \left( 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \right)^2} \begin{bmatrix} \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} & 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ -2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\ \end{bmatrix} \begin{bmatrix} F_u \\ F_v \end{bmatrix} \end{equation}1.3 Numerical Values
k = 2; m = 1; c = 0.05; xi = c/(2*sqrt(k*m)); w0 = sqrt(k/m);
1.4 Campbell Diagram
Compute the poles
wrs = linspace(0, 2, 51); % [rad/s] polesvc = zeros(4, length(wrs)); for i = 1:length(wrs) wr = wrs(i); polei = pole(1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wrs(i)^2)/(w0^2))^2 + (2*wrs(i)*s/(w0^2))^2)); [~, i_sort] = sort(imag(polei)); polesvc(:, i) = polei(i_sort); end
1.5 Simscape Model
open('rotating_frame.slx');
k = 2; m = 1; c = 0.05; xi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); wr = 0.1;
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/Fu'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Fv'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 1, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 2, 'openoutput'); io_i = io_i + 1;
G = linearize(mdl, io, 0); %% Input/Output definition G.InputName = {'Fu', 'Fv'}; G.OutputName = {'du', 'dv'};
1.6 Comparison with the model
Ga = (1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ... [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2), 2*wr*s/(w0^2) ; ... -2*wr*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)];
figure; bode(G, 'k-', Ga, 'r--')
2 Integral Force Feedback
2.1 IFF with pure integrator
2.1.1 Numerical Values
k = 1; m = 1; xi = 0.01; c = 2*xi*sqrt(k*m); w0 = sqrt(k/m); wr = 0.1*w0;
2.1.2 Equations
The sensed forces are equal to:
\begin{equation} \begin{bmatrix} F_{um} \\ F_{vm} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} F_u \\ F_v \end{bmatrix} - (c s + k) \begin{bmatrix} d_u \\ d_v \end{bmatrix} \end{equation}Giffa = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ... [(s^2/w0^2 - wr^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)) + (2*wr*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*wr*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*wr*s/(w0^2), (s^2/w0^2 - wr^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))+ (2*wr*s/(w0^2))^2];
2.1.3 Poles and Zeros
syms wr w0 xi positive assumealso(w0 > wr) syms x
z = (x^2/w0^2 - wr^2/w0^2)*((x^2)/(w0^2) + 1 - (wr^2)/(w0^2)) + (2*wr*x/(w0^2))^2 == 0 p = ((x^2)/(w0^2) + 1 - (wr^2)/(w0^2))^2 + (2*wr*x/(w0^2))^2 == 0
solve(p, x)
solve(z, x)
The zeros are the roots of:
\begin{equation} \left( \frac{s^2}{{\omega_0}^2} - \frac{\Omega^2}{{\omega_0}^2} \right) \left( \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \right) + \left( 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \right)^2 = 0 \end{equation}Poles (without damping)
\begin{equation} \left(\begin{array}{c} -w_{0}\,1{}\mathrm{i}-\mathrm{wr}\,1{}\mathrm{i}\\ -w_{0}\,1{}\mathrm{i}+\mathrm{wr}\,1{}\mathrm{i}\\ w_{0}\,1{}\mathrm{i}-\mathrm{wr}\,1{}\mathrm{i}\\ w_{0}\,1{}\mathrm{i}+\mathrm{wr}\,1{}\mathrm{i} \end{array}\right) \end{equation}Zeros (without damping)
\begin{equation} \left(\begin{array}{c} -\sqrt{-\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2}\\ -\sqrt{\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2}\\ \sqrt{-\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2}\\ \sqrt{\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2} \end{array}\right) \end{equation}2.1.4 Simscape Model
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/Fu'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Fv'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 5, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 6, 'openoutput'); io_i = io_i + 1;
Giff = linearize(mdl, io, 0); %% Input/Output definition Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'Fmu', 'Fmv'};
figure; bode(Giff, Giffa)
2.1.5 IFF Plant
2.1.6 Loop Gain
Let’s take \(\Omega = \frac{\omega_0}{10}\).
ws = 0.1*w0; G_iff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))^2 + (2*ws*s/(w0^2))^2) * ... [(s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2)) + (2*ws*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*ws*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*ws*s/(w0^2), (s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))+ (2*ws*s/(w0^2))^2];
2.1.7 Root Locus
2.2 Modified IFF (pseudo integrator)
2.2.1 Control Law
Let’s take the integral feedback controller as a low pass filter (pseudo integrator):
\begin{equation} K_{\text{IFF}} = \begin{bmatrix} g\frac{\omega_i}{\omega_i + s} & 0 \\ 0 & g\frac{\omega_i}{\omega_i + s} \end{bmatrix} \end{equation}xi = 0.005; w0 = 1; ws = 0.1*w0; G_iff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))^2 + (2*ws*s/(w0^2))^2) * ... [(s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2)) + (2*ws*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*ws*s/(w0^2) ; ... (2*xi*s/w0 + 1)*2*ws*s/(w0^2), (s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))+ (2*ws*s/(w0^2))^2];
g = 100; wi = ws; K_iff = (g/(1+s/wi))*eye(2);
2.2.2 Loop Gain
2.2.3 Root Locus
2.2.4 Optimal Gain
The DC gain for Giff is (for \(\Omega < \omega_0\)):
\begin{equation} G_{\text{IFF}}(\omega = 0) = \frac{1}{1 - \frac{{\omega_0}^2}{\Omega^2}} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \end{equation}The maximum gain where is system is still stable is
\begin{equation} g_\text{max} = \frac{{\omega_0}^2}{\Omega^2} - 1 \end{equation}Let’s find the gain that maximize the simultaneous damping of the two modes.
K_opt = (opt_gain/(1+s/wi))*eye(2); G_cl = feedback(G_iff, K_opt);
2.3 Stiffness in parallel with the force sensor
2.3.1 Schematic
2.3.2 Equations
The equations should be the same as before by taking \(k = k^\prime + k_a\). I just have to determine the measured force by the sensor
2.3.3 Parameters
k = 1; m = 1; c = 0.05; xi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); wr = 0.1;
2.3.4 IFF Plant
open('rotating_frame.slx');
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/Fu'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Fv'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 5, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 6, 'openoutput'); io_i = io_i + 1;
kp = 0; cp = 0;
Giff = linearize(mdl, io, 0); %% Input/Output definition Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'Fmu', 'Fmv'};
kp = 0.5*m*wr^2; cp = 0.01;
Giffa = linearize(mdl, io, 0); %% Input/Output definition Giffa.InputName = {'Fu', 'Fv'}; Giffa.OutputName = {'Fmu', 'Fmv'};
kp = 1.5*m*wr^2; cp = 0.01;
Giffb = linearize(mdl, io, 0); %% Input/Output definition Giffb.InputName = {'Fu', 'Fv'}; Giffb.OutputName = {'Fmu', 'Fmv'};
Comparison with the model
figure; bode(Giff, 'k-', Giffa, 'b--', Giffb, 'r--')
2.3.5 Parallel Stiffness effect
Pure IFF controller can be used if:
\begin{equation} k_{p} > m \Omega^2 \end{equation}However, having large values of \(k_p\) may:
- decrease the actuator stroke
- decrease the attainable damping (section about optimal value)
2.3.6 Root locus
2.3.7 Optimal value of \(k_p\)
kps = [0, 0.5, 1, 2, 10]*m*wr^2; cp = 0.01;
2.3.8 Physical Explanation
Negative stiffness Zeros are for high loop gain: remove the force sensor Thus the stiffness in parallel should be higher than the virtual negative stiffness added by the gyroscopic effects
2.4 Comparison
Fix the parameters, and see how the root locus change and the final damped system
k = 1; m = 1; c = 0.05; xi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); wr = 0.1;
kp = 0; cp = 0;
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/Fu'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Fv'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 5, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 6, 'openoutput'); io_i = io_i + 1;
Giff = linearize(mdl, io, 0); %% Input/Output definition Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'Fmu', 'Fmv'};
Modified IFF
wi = w0;
With parallel Stiffness
kp = 2*m*wr^2; cp = 0.01;
Giff_kp = linearize(mdl, io, 0); %% Input/Output definition Giff_kp.InputName = {'Fu', 'Fv'}; Giff_kp.OutputName = {'Fmu', 'Fmv'};
With parallel stiffness => unconditional stability and more damping but have to add mechanical parts
3 Direct Velocity Feedback
3.1 Equations
The sensed relative velocity are equal to:
3.2 Numerical Values
k = 1; m = 1; c = 0.05; xi = c/(2*sqrt(k*m)); w0 = sqrt(k/m);
Ga = (s/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ... [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2), 2*wr*s/(w0^2) ; ... -2*wr*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)];
3.3 Simscape Model
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/Fu'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Fv'], 1, 'openinput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 3, 'openoutput'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Tuv'], 4, 'openoutput'); io_i = io_i + 1;
G = linearize(mdl, io, 0); %% Input/Output definition G.InputName = {'Fu', 'Fv'}; G.OutputName = {'Vu', 'Vv'};
3.4 DVF Plant
figure; bode(G, Ga)
3.5 Loop Gain
Let’s take \(\Omega = \frac{\omega_0}{10}\).
ws = 0.1*w0; G_dvf = (s/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ... [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2), 2*wr*s/(w0^2) ; ... -2*wr*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)];
3.6 Root Locus
4 Comparison
4.1 Parameters
k = 1; m = 1; c = 0.05; xi = c/(2*sqrt(k*m)); w0 = sqrt(k/m); wr = 0.1;
4.2 Root Locus
4.2.1 Pseudo Integrator IFF
kp = 0; cp = 0;
%% 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;
Kiff = tf(zeros(2)); Kdvf = tf(zeros(2));
Giff = linearize(mdl, io, 0); %% Input/Output definition Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'Fmu', 'Fmv'};
4.2.2 IFF With parallel Stiffness
kp = 2*m*wr^2; cp = 0.01;
%% 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;
Giff_kp = linearize(mdl, io, 0); %% Input/Output definition Giff_kp.InputName = {'Fu', 'Fv'}; Giff_kp.OutputName = {'Fmu', 'Fmv'};
4.2.3 DVF
%% 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;
Gdvf = linearize(mdl, io, 0); %% Input/Output definition Gdvf.InputName = {'Fu', 'Fv'}; Gdvf.OutputName = {'Vu', 'Vv'};
4.2.4 Root Locus
wi = w0;
4.3 Controllers - Optimal Gains
Estimate the controller gain that yields good damping in all cases.
Pseudo integrator IFF: Parallel Stiffness DVF:
opt_zeta_iff, opt_zeta_kp, opt_zeta_dvf
4.4 Transmissibility
4.4.1 Open Loop
wi = w0;
Kdvf = tf(zeros(2)); Kiff = tf(zeros(2));
kp = 0; cp = 0;
%% 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, '/Meas'], 1, 'output'); io_i = io_i + 1;
Tol = linearize(mdl, io, 0); %% Input/Output definition Tol.InputName = {'Dwx', 'Dwy'}; Tol.OutputName = {'Dx', 'Dy'};
4.4.2 Pseudo Integrator IFF
wi = w0;
Kdvf = tf(zeros(2)); Kiff = opt_gain_iff/(1 + s/wi)*tf(eye(2));
kp = 0; cp = 0;
%% 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, '/Meas'], 1, 'output'); io_i = io_i + 1;
Tiff = linearize(mdl, io, 0); %% Input/Output definition Tiff.InputName = {'Dwx', 'Dwy'}; Tiff.OutputName = {'Dx', 'Dy'};
4.4.3 IFF With parallel Stiffness
kp = 2*m*wr^2; cp = 0.01;
Kdvf = tf(zeros(2)); Kiff = opt_gain_kp/s*tf(eye(2));
%% 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, '/Meas'], 1, 'output'); io_i = io_i + 1;
Tiff_kp = linearize(mdl, io, 0); %% Input/Output definition Tiff_kp.InputName = {'Dwx', 'Dwy'}; Tiff_kp.OutputName = {'Dx', 'Dy'};
4.4.4 DVF
kp = 0; cp = 0;
Kdvf = opt_gain_kp*tf(eye(2)); Kiff = tf(zeros(2));
%% 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, '/Meas'], 1, 'output'); io_i = io_i + 1;
Tdvf = linearize(mdl, io, 0); %% Input/Output definition Tdvf.InputName = {'Dwx', 'Dwy'}; Tdvf.OutputName = {'Dx', 'Dy'};
4.4.5 Transmissibility
freqs = logspace(-2, 1, 1000); figure; hold on; plot(freqs, abs(squeeze(freqresp(Tiff(1,1), freqs))), ... 'DisplayName', 'IFF Pseudo int') plot(freqs, abs(squeeze(freqresp(Tiff_kp(1,1), freqs))), ... 'DisplayName', 'IFF Paral. stiff') plot(freqs, abs(squeeze(freqresp(Tdvf(1,1), freqs))), ... 'DisplayName', 'DVF') plot(freqs, abs(squeeze(freqresp(Tol(1,1), freqs))), 'k-', ... 'DisplayName', 'IFF Pseudo int') hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Frequency [rad/s]'); ylabel('Amplitude [m/m]'); legend('location', 'northwest');
4.5 Compliance
4.5.1 Open Loop
Kdvf = tf(zeros(2)); Kiff = tf(zeros(2));
kp = 0; cp = 0;
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; 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;
Col = linearize(mdl, io, 0); %% Input/Output definition Col.InputName = {'Fdx', 'Fdy'}; Col.OutputName = {'Dx', 'Dy'};
4.5.2 Pseudo Integrator IFF
wi = w0;
Kdvf = tf(zeros(2)); Kiff = opt_gain_iff/(1 + s/wi)*tf(eye(2));
kp = 0; cp = 0;
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; 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;
Ciff = linearize(mdl, io, 0); %% Input/Output definition Ciff.InputName = {'Fdx', 'Fdy'}; Ciff.OutputName = {'Dx', 'Dy'};
4.5.3 IFF With parallel Stiffness
kp = 2*m*wr^2; cp = 0.01;
Kdvf = tf(zeros(2)); Kiff = opt_gain_kp/s*tf(eye(2));
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; 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;
Ciff_kp = linearize(mdl, io, 0); %% Input/Output definition Ciff_kp.InputName = {'Fdx', 'Fdy'}; Ciff_kp.OutputName = {'Dx', 'Dy'};
4.5.4 DVF
kp = 0; cp = 0;
Kdvf = opt_gain_kp*tf(eye(2)); Kiff = tf(zeros(2));
%% Name of the Simulink File mdl = 'rotating_frame'; %% Input/Output definition clear io; 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;
Cdvf = linearize(mdl, io, 0); %% Input/Output definition Cdvf.InputName = {'Fdx', 'Fdy'}; Cdvf.OutputName = {'Dx', 'Dy'};
4.5.5 Compliance
freqs = logspace(-2, 1, 1000); figure; hold on; plot(freqs, abs(squeeze(freqresp(Ciff(1,1), freqs))), ... 'DisplayName', 'IFF Pseudo int') plot(freqs, abs(squeeze(freqresp(Ciff_kp(1,1), freqs))), ... 'DisplayName', 'IFF Paral. stiff') plot(freqs, abs(squeeze(freqresp(Cdvf(1,1), freqs))), ... 'DisplayName', 'DVF') plot(freqs, abs(squeeze(freqresp(Col(1,1), freqs))), 'k-', ... 'DisplayName', 'IFF Pseudo int') hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Frequency [rad/s]'); ylabel('Compliance [m/N]'); legend('location', 'northwest'); % ylim([0, 1e3]); % xlim([freqs(1), freqs(end)]);