diff --git a/matlab/index.org b/matlab/index.org index 81f10de..b621d43 100644 --- a/matlab/index.org +++ b/matlab/index.org @@ -137,7 +137,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system end plot(Ws, zeros(size(Ws)), 'k--') hold off; - xlabel('Rotation Frequency [rad/s]'); ylabel('Real Part'); + xlabel('Rotational Speed [rad/s]'); ylabel('Real Part'); ax2 = subplot(1,2,2); hold on; @@ -145,7 +145,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system plot(Ws, imag(p_ws(p_i, :)), 'k-') end hold off; - xlabel('Rotation Frequency [rad/s]'); ylabel('Imaginary Part'); + xlabel('Rotational Speed [rad/s]'); ylabel('Imaginary Part'); #+end_src #+begin_src matlab :tangle no :exports results :results file replace @@ -172,7 +172,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system plot(Ws, real(p_ws(3, :)), '-', 'HandleVisibility', 'off') plot(Ws, zeros(size(Ws)), 'k--', 'HandleVisibility', 'off') hold off; - xlabel('Rotation Frequency $\Omega$'); ylabel('Real Part'); + xlabel('Rotational Speed $\Omega$'); ylabel('Real Part'); xlim([0, 2*w0]); xticks([0,w0/2,w0,3/2*w0,2*w0]) xticklabels({'$0$', '', '$\omega_0$', '', '$2 \omega_0$'}) @@ -193,7 +193,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system plot(Ws, imag(p_ws(3, :)), '-') plot(Ws, zeros(size(Ws)), 'k--') hold off; - xlabel('Rotation Frequency $\Omega$'); ylabel('Imaginary Part'); + xlabel('Rotational Speed $\Omega$'); ylabel('Imaginary Part'); xlim([0, 2*w0]); xticks([0,w0/2,w0,3/2*w0,2*w0]) xticklabels({'$0$', '', '$\omega_0$', '', '$2 \omega_0$'}) @@ -482,7 +482,7 @@ Which then gives: \end{align} #+end_important -** Simscape Model +** Comparison of the Analytical Model and the Simscape Model The rotation speed is set to $\Omega = 0.1 \omega_0$. #+begin_src matlab W = 0.1*w0; % [rad/s] @@ -515,7 +515,6 @@ And the transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is identified using Giff.OutputName = {'fu', 'fv'}; #+end_src -** Comparison of the Analytical Model and the Simscape Model The same transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is written down from the analytical model. #+begin_src matlab Giff_th = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... @@ -524,7 +523,6 @@ The same transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is written down fro #+end_src The two are compared in Figure [[fig:plant_iff_comp_simscape_analytical]] and found to perfectly match. - #+begin_src matlab :exports none freqs = logspace(-1, 1, 1000); @@ -1123,40 +1121,36 @@ While it seems that small $\omega_i$ do allow more damping to be added to the sy There must be an optimum for $\omega_i$. To find the optimum, the gain that maximize the simultaneous damping of the mode is identified for a wide range of $\omega_i$ (Figure [[fig:mod_iff_damping_wi]]). -#+begin_src matlab - wis = logspace(-2, 1, 31)*w0; % [rad/s] - opt_zeta = zeros(1, length(wis)); % Optimal simultaneous damping +#+begin_src matlab + wis = 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); - gains = linspace(0, (w0^2/W^2 - 1)*wi, 100); + Kiff = 1/(s + wi)*eye(2); - for g = gains - Kiff = (g/(wi+s))*eye(2); + fun = @(g)computeSimultaneousDamping(g, Giff, Kiff); - [w, zeta] = damp(minreal(feedback(Giff, Kiff))); - - if min(zeta) > opt_zeta(wi_i) && all(zeta > 0) - opt_zeta(wi_i) = min(zeta); - opt_gain(wi_i) = g; - end - end + [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 #+end_src #+begin_src matlab :exports none figure; yyaxis left - plot(wis, opt_zeta, '-o', 'DisplayName', '$\xi_{cl}$'); + plot(wis, opt_xi, '-', 'DisplayName', '$\xi_{cl}$'); set(gca, 'YScale', 'lin'); ylim([0,1]); ylabel('Attainable Damping Ratio $\xi$'); yyaxis right hold on; - plot(wis, opt_gain, '-x', 'DisplayName', '$g_{opt}$'); + plot(wis, opt_gain, '-', 'DisplayName', '$g_{opt}$'); plot(wis, wis*((w0/W)^2 - 1), '--', 'DisplayName', '$g_{max}$'); set(gca, 'YScale', 'lin'); ylim([0,10]); @@ -1792,7 +1786,63 @@ It is shown that large values of $k_p$ decreases the attainable damping. exportFig('figs-inkscape/root_locus_iff_kps.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true); #+end_src -** Optimal Gain +#+begin_src matlab + alphas = 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 +#+end_src + +#+begin_src matlab :exports none + figure; + yyaxis left + plot(alphas, opt_xi, '-', 'DisplayName', '$\xi_{cl}$'); + set(gca, 'YScale', 'lin'); + ylim([0,1]); + ylabel('Attainable Damping Ratio $\xi$'); + + yyaxis right + hold on; + plot(alphas, opt_gain, '-', 'DisplayName', '$g_{opt}$'); + set(gca, 'YScale', 'lin'); + ylim([0,2.5]); + ylabel('Controller gain $g$'); + + xlabel('$\alpha$'); + set(gca, 'XScale', 'log'); + legend('location', 'northeast'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/opt_damp_alpha.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:opt_damp_alpha +#+caption: +#+RESULTS: +[[file:figs/opt_damp_alpha.png]] + +** TODO Optimal Gain Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such that maximum damping are added to the poles of the closed loop system. #+begin_src matlab @@ -1808,7 +1858,7 @@ Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such t #+end_src #+begin_src matlab - opt_zeta = 0; + opt_xi = 0; opt_gain = 0; gains = logspace(-2, 4, 1000); @@ -1816,10 +1866,10 @@ Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such t for g = gains Kiff = (g/s)*eye(2); - [w, zeta] = damp(minreal(feedback(Giff, Kiff))); + [w, xi] = damp(minreal(feedback(Giff, Kiff))); - if min(zeta) > opt_zeta && all(zeta > 0) - opt_zeta = min(zeta); + if min(xi) > opt_xi && all(xi > 0) + opt_xi = min(xi); opt_gain = min(g); end end @@ -2061,54 +2111,41 @@ IFF With parallel Stiffness In order to compare to three considered Active Damping techniques, gains that yield maximum damping of all the modes are computed for each case. #+begin_src matlab :exports none - %% IFF with pseudo integrators - gains = linspace(0, (w0^2/W^2 - 1)*wi, 100); - opt_zeta_iff = 0; - opt_gain_iff = 0; + fun = @(g)computeSimultaneousDamping(g, Giff, (1/(wi+s))*eye(2)); - for g = gains - Kiff = (g/(wi+s))*eye(2); - - [w, zeta] = damp(minreal(feedback(Giff, Kiff))); - - if min(zeta) > opt_zeta_iff && all(zeta > 0) - opt_zeta_iff = min(zeta); - opt_gain_iff = g; - end - end + [opt_gain_iff, opt_xi_iff] = fminsearch(fun, 0.5*(w0^2/W^2 - 1)*wi); + opt_xi_iff = 1/opt_xi_iff; #+end_src #+begin_src matlab :exports none - %% IFF with Parallel Stiffness - gains = logspace(-2, 4, 100); - opt_zeta_kp = 0; - opt_gain_kp = 0; + fun = @(g)computeSimultaneousDamping(g, Giff_kp, 1/s*eye(2)); - for g = gains - Kiff = g/s*eye(2); - - [w, zeta] = damp(minreal(feedback(Giff_kp, Kiff))); - - if min(zeta) > opt_zeta_kp && all(zeta > 0) - opt_zeta_kp = min(zeta); - opt_gain_kp = g; - end - end + [opt_gain_kp, opt_xi_kp] = fminsearch(fun, 2); + opt_xi_kp = 1/opt_xi_kp; #+end_src The obtained damping ratio and control are shown below. #+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) - data2orgtable([opt_zeta_iff, opt_zeta_kp; opt_gain_iff, opt_gain_kp]', {'Modified IFF', 'IFF with $k_p$'}, {'Obtained $\xi$', 'Control Gain'}, ' %.2f '); + data2orgtable([opt_xi_iff, opt_xi_kp; opt_gain_iff, opt_gain_kp]', {'Modified IFF', 'IFF with $k_p$'}, {'Obtained $\xi$', 'Control Gain'}, ' %.2f '); #+end_src #+RESULTS: | | Obtained $\xi$ | Control Gain | |----------------+----------------+--------------| -| Modified IFF | 0.83 | 2.0 | -| IFF with $k_p$ | 0.83 | 2.01 | +| Modified IFF | 0.83 | 1.99 | +| IFF with $k_p$ | 0.83 | 2.02 | ** Passive Damping - Critical Damping +\begin{equation} + \xi = \frac{c}{2 \sqrt{km}} +\end{equation} + +Critical Damping corresponds to to $\xi = 1$, and thus: +\begin{equation} + c_{\text{crit}} = 2 \sqrt{km} +\end{equation} + #+begin_src matlab c_opt = 2*sqrt(k*m); #+end_src @@ -2226,12 +2263,13 @@ The obtained damping ratio and control are shown below. 'DisplayName', 'Open-Loop') hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylim([1e-2, 3e1]); xlabel('Frequency [rad/s]'); ylabel('Transmissibility [m/m]'); legend('location', 'southwest'); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/comp_transmissibility.pdf', 'width', 'wide', 'height', 'tall'); + exportFig('figs/comp_transmissibility.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:comp_transmissibility @@ -2240,7 +2278,7 @@ The obtained damping ratio and control are shown below. [[file:figs/comp_transmissibility.png]] #+begin_src matlab :tangle no :exports none :results none - exportFig('figs-inkscape/comp_transmissibility.pdf', 'width', 'half', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true); + exportFig('figs-inkscape/comp_transmissibility.pdf', 'width', 'half', 'height', 'normal', 'png', false, 'pdf', false, 'svg', true); #+end_src *** Compliance :ignore: @@ -2259,12 +2297,13 @@ The obtained damping ratio and control are shown below. 'DisplayName', 'Open-Loop') hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylim([1e-2, 3e1]); xlabel('Frequency [rad/s]'); ylabel('Compliance [m/N]'); legend('location', 'southwest'); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/comp_compliance.pdf', 'width', 'wide', 'height', 'tall'); + exportFig('figs/comp_compliance.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:comp_compliance @@ -2273,7 +2312,263 @@ The obtained damping ratio and control are shown below. [[file:figs/comp_compliance.png]] #+begin_src matlab :tangle no :exports none :results none - exportFig('figs-inkscape/comp_compliance.pdf', 'width', 'half', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true); + exportFig('figs-inkscape/comp_compliance.pdf', 'width', 'half', 'height', 'normal', 'png', false, 'pdf', false, 'svg', true); +#+end_src + +** DC Compliance +*** Pseudo Integrator IFF :ignore: +#+begin_src matlab :exports none + k = 1; + m = 1; + w0 = sqrt(k/m); +#+end_src + +#+begin_src matlab :exports none + Gwi = tf(zeros(4,4)); + Gwi.InputName = {'Fx', 'Fy', 'Fu', 'Fv'}; + Gwi.OutputName = {'dx', 'dy', 'fu', 'fv'}; + + Gp = ((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2; + + Gwi('dx', 'Fu') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp; + Gwi('dy', 'Fv') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp; + + Gwi('dx', 'Fx') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp; + Gwi('dy', 'Fy') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp; + + Gwi('dx', 'Fv') = (1/k)*(2*W*s/(w0^2))/Gp; + Gwi('dy', 'Fu') = -(1/k)*(2*W*s/(w0^2))/Gp; + + Gwi('dx', 'Fy') = (1/k)*(2*W*s/(w0^2))/Gp; + Gwi('dy', 'Fx') = -(1/k)*(2*W*s/(w0^2))/Gp; + + Gwi('fu', 'Fu') = ((s^2/w0^2 - W^2/w0^2)*(s^2/w0^2 + 2*xi*s/w0 + 1 - W^2/w0^2) + (2*(s/w0)*(W/w0))^2)/Gp; + Gwi('fv', 'Fv') = ((s^2/w0^2 - W^2/w0^2)*(s^2/w0^2 + 2*xi*s/w0 + 1 - W^2/w0^2) + (2*(s/w0)*(W/w0))^2)/Gp; + + Gwi('fu', 'Fv') = -(2*xi*s/w0 + 1)*(2*(s/w0)*(W/w0))/Gp; + Gwi('fv', 'Fu') = (2*xi*s/w0 + 1)*(2*(s/w0)*(W/w0))/Gp; + + Gwi('fu', 'Fx') = -(c*s + k)*(1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp; + Gwi('fv', 'Fy') = -(c*s + k)*(1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp; + + Gwi('fu', 'Fy') = -(c*s + k)*(1/k)*(2*W*s/(w0^2))/Gp; + Gwi('fv', 'Fx') = (c*s + k)*(1/k)*(2*W*s/(w0^2))/Gp; +#+end_src + +#+begin_src matlab + wis = logspace(-2, 1, 100)*w0; % [rad/s] + + opt_xi_wi = zeros(1, length(wis)); % Optimal simultaneous damping + opt_gain_wi = zeros(1, length(wis)); % Corresponding optimal gain + C_dc_gain_wi = zeros(1, length(wis)); % Compliance DC value + + for wi_i = 1:length(wis) + wi = wis(wi_i); + Kiff = 1/(s + wi)*eye(2); + + fun = @(g)computeSimultaneousDamping(g, Gwi({'fu', 'fv'}, {'Fu', 'Fv'}), Kiff); + + [g_opt, xi_opt] = fminsearch(fun, 0.5*wi*((w0/W)^2 - 1)); + opt_xi_wi(wi_i) = 1/xi_opt; + opt_gain_wi(wi_i) = g_opt; + + G_dc_gain = dcgain(lft(Gwi, -g_opt/(s + wi)*eye(2), 2, 2)); + C_dc_gain_wi(wi_i) = G_dc_gain(1,1); + end +#+end_src + +*** IFF With parallel Stiffness :ignore: +#+begin_src matlab + alphas = logspace(-2, 0, 100); + + opt_xi_kp = zeros(1, length(alphas)); % Optimal simultaneous damping + opt_gain_kp = zeros(1, length(alphas)); % Corresponding optimal gain + C_dc_gain_kp = zeros(1, length(alphas)); % DC gain of the compliance + + 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)); + + Gkp = tf(zeros(4,4)); + Gkp.InputName = {'Fx', 'Fy', 'Fu', 'Fv'}; + Gkp.OutputName = {'dx', 'dy', 'fu', 'fv'}; + + Gp = ((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))^2 + (2*W*s/(w0p^2))^2; + + Gkp('dx', 'Fu') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + Gkp('dy', 'Fv') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + + Gkp('dx', 'Fx') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + Gkp('dy', 'Fy') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + + Gkp('dx', 'Fv') = (1/(k + kp))*(2*W*s/(w0p^2))/Gp; + Gkp('dy', 'Fu') = -(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + + Gkp('dx', 'Fy') = (1/(k + kp))*(2*W*s/(w0p^2))/Gp; + Gkp('dy', 'Fx') = -(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + + Gkp('fu', 'Fu') = ((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)/Gp; + Gkp('fv', 'Fv') = ((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)/Gp; + + Gkp('fu', 'Fv') = -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p))/Gp; + Gkp('fv', 'Fu') = (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p))/Gp; + + Gkp('fu', 'Fx') = -(c*s + k)*(1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + Gkp('fv', 'Fy') = -(c*s + k)*(1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + + Gkp('fu', 'Fy') = -(c*s + k)*(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + Gkp('fv', 'Fx') = (c*s + k)*(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + + fun = @(g)computeSimultaneousDamping(g, Gkp({'fu', 'fv'}, {'Fu', 'Fv'}), Kiff); + + [g_opt, xi_opt] = fminsearch(fun, 2); + opt_xi_kp(alpha_i) = 1/xi_opt; + opt_gain_kp(alpha_i) = g_opt; + + G_dc_gain = dcgain(lft(Gkp, -g_opt/s*eye(2), 2, 2)); + C_dc_gain_kp(alpha_i) = G_dc_gain(1,1); + end +#+end_src + +*** Comparison :ignore: +#+begin_src matlab :exports none + figure; + yyaxis left + plot(alphas, C_dc_gain_kp, '-', 'DisplayName', '$|T(0)|$'); + set(gca, 'YScale', 'log'); + ylim([0, 1e3]); + ylabel('DC value of the Compliance'); + + yyaxis right + hold on; + plot(alphas, opt_xi_kp, '-', 'DisplayName', '$\xi_{opt}$'); + set(gca, 'YScale', 'lin'); + ylim([0,1]); + ylabel('Optimal Damping Ratio'); + + xlabel('$\alpha$'); + set(gca, 'XScale', 'log'); + legend('location', 'northeast'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace +exportFig('figs/compliance_dc_gain_wi.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:compliance_dc_gain_wi +#+caption: +#+RESULTS: +[[file:figs/compliance_dc_gain_wi.png]] + +#+begin_src matlab :exports none + figure; + yyaxis left + plot(wis, abs(C_dc_gain_wi), '-', 'DisplayName', '$|T(0)|$'); + set(gca, 'YScale', 'log'); + ylim([0, 1e3]); + ylabel('DC value of the Compliance'); + + yyaxis right + hold on; + plot(wis, opt_xi_wi, '-', 'DisplayName', '$\xi_{opt}$'); + set(gca, 'YScale', 'lin'); + ylim([0,1]); + ylabel('Optimal Damping Ratio'); + + xlabel('$\omega_i$'); + set(gca, 'XScale', 'log'); + legend('location', 'northeast'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace +exportFig('figs/compliance_dc_gain_kp.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:compliance_dc_gain_kp +#+caption: +#+RESULTS: +[[file:figs/compliance_dc_gain_kp.png]] + +#+begin_src matlab :exports none + figure; + hold on; + plot(opt_xi_wi, C_dc_gain_wi, '-', 'DisplayName', '$\omega_i$'); + plot(opt_xi_kp, C_dc_gain_kp, '-', 'DisplayName', '$k_p$'); + hold off + set(gca, 'YScale', 'log'); + ylim([0, 1e3]); + ylabel('DC value of the Compliance'); + xlabel('Optimal Damping Ratio'); + legend('location', 'northwest'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace +exportFig('figs/opt_damp_vs_dc_comp.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:opt_damp_vs_dc_comp +#+caption: +#+RESULTS: +[[file:figs/opt_damp_vs_dc_comp.png]] + +*** Comparison :ignore: +#+begin_src matlab :exports none + figure; + yyaxis left + plot(wis, opt_xi_wi, '-', 'DisplayName', '$\xi_{cl}$'); + set(gca, 'YScale', 'lin'); + ylim([0,1]); + ylabel('Attainable Damping Ratio $\xi$'); + + yyaxis right + hold on; + plot(wis, opt_gain_wi, '-', 'DisplayName', '$g_{opt}$'); + plot(wis, wis*((w0/W)^2 - 1), '--', 'DisplayName', '$g_{max}$'); + set(gca, 'YScale', 'lin'); + ylim([0,10]); + ylabel('Controller gain $g$'); + + xlabel('$\omega_i/\omega_0$'); + set(gca, 'XScale', 'log'); + legend('location', 'northeast'); +#+end_src + +#+begin_src matlab :exports none + figure; + yyaxis left + plot(alphas, opt_xi_kp, '-', 'DisplayName', '$\xi_{cl}$'); + set(gca, 'YScale', 'lin'); + ylim([0,1]); + ylabel('Attainable Damping Ratio $\xi$'); + + yyaxis right + hold on; + plot(alphas, opt_gain_kp, '-', 'DisplayName', '$g_{opt}$'); + set(gca, 'YScale', 'lin'); + ylim([0,2.5]); + ylabel('Controller gain $g$'); + + xlabel('$\alpha$'); + set(gca, 'XScale', 'log'); + legend('location', 'northeast'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/mod_iff_damping_kp.pdf', 'width', 'normal', 'height', 'normal'); +#+end_src + +#+name: fig:mod_iff_damping_kp +#+caption: +#+RESULTS: +[[file:figs/mod_iff_damping_kp.png]] + +#+begin_src matlab :tangle no :exports none :results none + exportFig('figs-inkscape/mod_iff_damping_kp.pdf', 'width', 'half', 'height', '650', 'png', false, 'pdf', false, 'svg', true); #+end_src * Notations @@ -2452,3 +2747,15 @@ This Matlab function is accessible [[file:src/rootLocusPolesSorted.m][here]]. #+begin_src matlab poles = poles.'; #+end_src + +** =computeSimultaneousDamping= +#+begin_src matlab :tangle src/computeSimultaneousDamping.m + function [xi_min] = computeSimultaneousDamping(g, G, K) + [w, xi] = damp(minreal(feedback(G, g*K))); + xi_min = 1/min(xi); + + if xi_min < 0 + xi_min = 1e8; + end + end +#+end_src diff --git a/matlab/matlab/rotating_frame.slx b/matlab/matlab/rotating_frame.slx index adebb6c..f995f8d 100644 Binary files a/matlab/matlab/rotating_frame.slx and b/matlab/matlab/rotating_frame.slx differ diff --git a/matlab/matlab/s4_iff_kp.m b/matlab/matlab/s4_iff_kp.m index 0f48a51..818835d 100644 --- a/matlab/matlab/s4_iff_kp.m +++ b/matlab/matlab/s4_iff_kp.m @@ -1,389 +1,79 @@ -%% Clear Workspace and Close figures -clear; close all; clc; +% Attainable Damping as a function of $k_p$ -%% Intialize Laplace variable -s = zpk('s'); +tic; +alphas = logspace(-2, 0, 10); +gains = linspace(0.5, 2.5, 100); -% Plant Parameters -% Let's define initial values for the model. +opt_zeta = zeros(1, length(alphas)); % Optimal simultaneous damping +opt_gain = zeros(1, length(alphas)); % Corresponding optimal gain -k = 1; % Actuator Stiffness [N/m] -c = 0.05; % Actuator Damping [N/(m/s)] -m = 1; % Payload mass [kg] - -xi = c/(2*sqrt(k*m)); -w0 = sqrt(k/m); % [rad/s] - -kp = 0; % [N/m] -cp = 0; % [N/(m/s)] - -% Comparison of the Analytical Model and the Simscape Model -% The same transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is written down from the analytical model. - -W = 0.1*w0; % [rad/s] - -kp = 1.5*m*W^2; -cp = 0; - -Kiff = tf(zeros(2)); -Kdvf = tf(zeros(2)); - -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, '/K'], 1, 'openinput'); io_i = io_i + 1; -io(io_i) = linio([mdl, '/G'], 2, 'openoutput'); io_i = io_i + 1; - -Giff = linearize(mdl, io, 0); - -%% Input/Output definition -Giff.InputName = {'Fu', 'Fv'}; -Giff.OutputName = {'fu', 'fv'}; - -w0p = 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'}; - -freqs = logspace(-1, 1, 1000); - -figure; -ax1 = subplot(2, 2, 1); -hold on; -plot(freqs, abs(squeeze(freqresp(Giff(1,1), freqs))), '-') -plot(freqs, abs(squeeze(freqresp(Giff_th(1,1), freqs))), '--') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]'); -title('$f_u/F_u$, $f_v/F_v$'); - -ax3 = subplot(2, 2, 3); -hold on; -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,1), freqs))), '-') -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_th(1,1), freqs))), '--') -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [rad/s]'); ylabel('Phase [deg]'); -yticks(-180:90:180); -ylim([-180 180]); -hold off; - -ax2 = subplot(2, 2, 2); -hold on; -plot(freqs, abs(squeeze(freqresp(Giff(1,2), freqs))), '-') -plot(freqs, abs(squeeze(freqresp(Giff_th(1,2), freqs))), '--') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]'); -title('$f_u/F_v$, $f_v/F_u$'); - -ax4 = subplot(2, 2, 4); -hold on; -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,2), freqs))), '-', ... - 'DisplayName', 'Simscape') -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_th(1,2), freqs))), '--', ... - 'DisplayName', 'Analytical') -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [rad/s]'); ylabel('Phase [deg]'); -yticks(-180:90:180); -ylim([-180 180]); -hold off; -legend('location', 'northeast'); - -linkaxes([ax1,ax2,ax3,ax4],'x'); -xlim([freqs(1), freqs(end)]); -linkaxes([ax1,ax2],'y'); - -% Effect of the parallel stiffness on the IFF plant -% The rotation speed is set to $\Omega = 0.1 \omega_0$. - -W = 0.1*w0; % [rad/s] - - - -% And the IFF plant (transfer function from $[F_u, F_v]$ to $[f_u, f_v]$) is identified in three different cases: -% - without parallel stiffness -% - with a small parallel stiffness $k_p < m \Omega^2$ -% - with a large parallel stiffness $k_p > m \Omega^2$ - -% The results are shown in Figure [[fig:plant_iff_kp]]. - -% One can see that for $k_p > m \Omega^2$, the systems shows alternating complex conjugate poles and zeros. - - -kp = 0; -cp = 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]; - -kp = 0.5*m*W^2; -cp = 0; - -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]; - -kp = 1.5*m*W^2; -cp = 0; - -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]; - -freqs = logspace(-2, 1, 1000); - -figure; - -ax1 = subplot(2, 1, 1); -hold on; -plot(freqs, abs(squeeze(freqresp(Giff(1,1), freqs))), 'k-') -plot(freqs, abs(squeeze(freqresp(Giff_s(1,1), freqs))), 'k--') -plot(freqs, abs(squeeze(freqresp(Giff_l(1,1), freqs))), 'k:') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]'); - -ax2 = subplot(2, 1, 2); -hold on; -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,1), freqs))), 'k-', ... - 'DisplayName', '$k_p = 0$') -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_s(1,1), freqs))), 'k--', ... - 'DisplayName', '$k_p < m\Omega^2$') -plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_l(1,1), freqs))), 'k:', ... - 'DisplayName', '$k_p > m\Omega^2$') -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [rad/s]'); ylabel('Phase [deg]'); -yticks(-180:90:180); -ylim([-180 180]); -hold off; -legend('location', 'southwest'); - -linkaxes([ax1,ax2],'x'); -xlim([freqs(1), freqs(end)]); - -% IFF when adding a spring in parallel -% In Figure [[fig:root_locus_iff_kp]] is displayed the Root Locus in the three considered cases with -% \begin{equation} -% K_{\text{IFF}} = \frac{g}{s} \begin{bmatrix} -% 1 & 0 \\ -% 0 & 1 -% \end{bmatrix} -% \end{equation} - -% One can see that for $k_p > m \Omega^2$, the root locus stays in the left half of the complex plane and thus the control system is unconditionally stable. - -% Thus, decentralized IFF controller with pure integrators can be used if: -% \begin{equation} -% k_{p} > m \Omega^2 -% \end{equation} - - -figure; - -gains = logspace(-2, 2, 100); - -subplot(1,2,1); -hold on; -set(gca,'ColorOrderIndex',1); -plot(real(pole(Giff)), imag(pole(Giff)), 'x', ... - 'DisplayName', '$k_p = 0$'); -set(gca,'ColorOrderIndex',1); -plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', ... - 'HandleVisibility', 'off'); -for g = gains - cl_poles = pole(feedback(Giff, (g/s)*eye(2))); - set(gca,'ColorOrderIndex',1); - plot(real(cl_poles), imag(cl_poles), '.', ... - 'HandleVisibility', 'off'); -end - -set(gca,'ColorOrderIndex',2); -plot(real(pole(Giff_s)), imag(pole(Giff_s)), 'x', ... - 'DisplayName', '$k_p < m\Omega^2$'); -set(gca,'ColorOrderIndex',2); -plot(real(tzero(Giff_s)), imag(tzero(Giff_s)), 'o', ... - 'HandleVisibility', 'off'); -for g = gains - cl_poles = pole(feedback(Giff_s, (g/s)*eye(2))); - set(gca,'ColorOrderIndex',2); - plot(real(cl_poles), imag(cl_poles), '.', ... - 'HandleVisibility', 'off'); -end - -set(gca,'ColorOrderIndex',3); -plot(real(pole(Giff_l)), imag(pole(Giff_l)), 'x', ... - 'DisplayName', '$k_p > m\Omega^2$'); -set(gca,'ColorOrderIndex',3); -plot(real(tzero(Giff_l)), imag(tzero(Giff_l)), 'o', ... - 'HandleVisibility', 'off'); -for g = gains - set(gca,'ColorOrderIndex',3); - cl_poles = pole(feedback(Giff_l, (g/s)*eye(2))); - plot(real(cl_poles), imag(cl_poles), '.', ... - 'HandleVisibility', 'off'); -end -hold off; -axis square; -xlim([-1, 0.2]); ylim([0, 1.2]); - -xlabel('Real Part'); ylabel('Imaginary Part'); -legend('location', 'northwest'); - -subplot(1,2,2); -hold on; -set(gca,'ColorOrderIndex',1); -plot(real(pole(Giff)), imag(pole(Giff)), 'x'); -set(gca,'ColorOrderIndex',1); -plot(real(tzero(Giff)), imag(tzero(Giff)), 'o'); -for g = gains - cl_poles = pole(feedback(Giff, (g/s)*eye(2))); - set(gca,'ColorOrderIndex',1); - plot(real(cl_poles), imag(cl_poles), '.'); -end - -set(gca,'ColorOrderIndex',2); -plot(real(pole(Giff_s)), imag(pole(Giff_s)), 'x'); -set(gca,'ColorOrderIndex',2); -plot(real(tzero(Giff_s)), imag(tzero(Giff_s)), 'o'); -for g = gains - cl_poles = pole(feedback(Giff_s, (g/s)*eye(2))); - set(gca,'ColorOrderIndex',2); - plot(real(cl_poles), imag(cl_poles), '.'); -end - -set(gca,'ColorOrderIndex',3); -plot(real(pole(Giff_l)), imag(pole(Giff_l)), 'x'); -set(gca,'ColorOrderIndex',3); -plot(real(tzero(Giff_l)), imag(tzero(Giff_l)), 'o'); -for g = gains - set(gca,'ColorOrderIndex',3); - cl_poles = pole(feedback(Giff_l, (g/s)*eye(2))); - plot(real(cl_poles), imag(cl_poles), '.'); -end -hold off; -axis square; -xlim([-0.04, 0.06]); ylim([0, 0.1]); - -xlabel('Real Part'); ylabel('Imaginary Part'); - -% Effect of $k_p$ on the attainable damping -% However, having large values of $k_p$ may: -% - decrease the actuator force authority -% - decrease the attainable damping - -% To study the second point, Root Locus plots for the following values of $k_p$ are shown in Figure [[fig:root_locus_iff_kps]]. - -kps = [2, 20, 40]*m*W^2; - - - -% It is shown that large values of $k_p$ decreases the attainable damping. - - -figure; - -gains = logspace(-2, 4, 500); - -hold on; -for kp_i = 1:length(kps) - kp = kps(kp_i); +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 ]; + (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]; - set(gca,'ColorOrderIndex',kp_i); - plot(real(pole(Giff)), imag(pole(Giff)), 'x', ... - 'DisplayName', sprintf('$k_p = %.1f m \\Omega^2$', kp/(m*W^2))); - set(gca,'ColorOrderIndex',kp_i); - plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', ... - 'HandleVisibility', 'off'); for g = gains - Kiffa = (g/s)*eye(2); - cl_poles = pole(feedback(Giff, Kiffa)); - set(gca,'ColorOrderIndex',kp_i); - plot(real(cl_poles), imag(cl_poles), '.', ... - 'HandleVisibility', 'off'); - end -end -hold off; -axis square; -xlim([-1.2, 0.2]); ylim([0, 1.4]); - -xlabel('Real Part'); ylabel('Imaginary Part'); -legend('location', 'northwest'); - -% Optimal Gain -% Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such that maximum damping are added to the poles of the closed loop system. - - -kp = 5*m*W^2; -cp = 0.01; - -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 ]; - -opt_zeta = 0; -opt_gain = 0; - -gains = logspace(-2, 4, 1000); - -for g = gains - Kiff = (g/s)*eye(2); - - [w, zeta] = damp(minreal(feedback(Giff, Kiff))); - - if min(zeta) > opt_zeta && all(zeta > 0) - opt_zeta = min(zeta); - opt_gain = min(g); + Kiff = g/s*eye(2); + + [w, zeta] = damp(minreal(feedback(Giff, Kiff))); + + if min(zeta) > opt_zeta(alpha_i) && all(zeta > 0) + opt_zeta(alpha_i) = min(zeta); + opt_gain(alpha_i) = g; + end end end +toc figure; +yyaxis left +plot(alphas, opt_zeta, '-o', 'DisplayName', '$\xi_{cl}$'); +set(gca, 'YScale', 'lin'); +ylim([0,1]); +ylabel('Attainable Damping Ratio $\xi$'); -gains = logspace(-2, 4, 1000); - +yyaxis right hold on; -plot(real(pole(Giff)), imag(pole(Giff)), 'kx'); -plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko'); -for g = gains - clpoles = pole(minreal(feedback(Giff, (g/s)*eye(2)))); - plot(real(clpoles), imag(clpoles), 'k.'); +plot(alphas, opt_gain, '-x', 'DisplayName', '$g_{opt}$'); +set(gca, 'YScale', 'lin'); +ylim([0,3]); +ylabel('Controller gain $g$'); + +xlabel('$\alpha$'); +set(gca, 'XScale', 'log'); +legend('location', 'northeast'); + + + +% Alternative using fminserach + +alphas = logspace(-2, 0, 100); + +opt_zeta = 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_zeta(alpha_i) = 1/xi_opt; + opt_gain(alpha_i) = g_opt; end -% Optimal Gain -clpoles = pole(minreal(feedback(Giff, (opt_gain/s)*eye(2)))); -set(gca,'ColorOrderIndex',1); -plot(real(clpoles), imag(clpoles), 'x'); -for clpole = clpoles' - set(gca,'ColorOrderIndex',1); - plot([0, real(clpole)], [0, imag(clpole)], '-', 'LineWidth', 1); -end -hold off; -axis square; -xlim([-1.2, 0.2]); ylim([0, 1.4]); -xlabel('Real Part'); ylabel('Imaginary Part');