%% Clear Workspace and Close figures clear; close all; clc; %% Intialize Laplace variable s = zpk('s'); %% Path for functions, data and scripts addpath('./mat/'); % Path for data addpath('./src/'); % Path for Functions %% Colors for the figures colors = colororder; %% Simscape model name mdl = 'rotating_model'; %% Load "Generic" system dynamics load('rotating_generic_plants.mat', 'Gs', 'Wzs'); % Decentralized Relative Damping Control % The transfer functions from $[F_u,\ F_v]$ to $[d_u,\ d_v]$ is identified and shown in Figure ref:fig:rotating_rdc_plant_effect_rot for several rotating velocities. %% Bode plot of the direct and coupling term for the "relative damping control" plant - Effect of rotation figure; freqs = logspace(-2, 1, 1000); figure; tiledlayout(3, 2, 'TileSpacing', 'Compact', 'Padding', 'None'); % Magnitude ax1 = nexttile([2, 1]); hold on; for i = 1:length(Wzs) plot(freqs, abs(squeeze(freqresp(Gs{i}('du', 'Fu'), freqs, 'rad/s'))), '-', 'color', colors(i,:)) end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]'); title('Direct terms: $d_u/F_u$, $d_v/F_v$'); ylim([1e-3, 1e2]); ax2 = nexttile([2, 1]); hold on; for i = 1:length(Wzs) plot(freqs, abs(squeeze(freqresp(Gs{i}('dv', 'Fu'), freqs, 'rad/s'))), '-', 'color', colors(i,:), ... 'DisplayName', sprintf('$\\Omega = %.1f \\omega_0$', Wzs(i))) end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); set(gca, 'XTickLabel',[]); set(gca, 'YTickLabel',[]); title('Coupling Terms: $d_u/F_v$, $d_v/F_u$'); ldg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1); ldg.ItemTokenSize = [10, 1]; ylim([1e-3, 1e2]); ax3 = nexttile; hold on; for i = 1:length(Wzs) plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}('du', 'Fu'), freqs, 'rad/s'))), '-', 'color', colors(i,:)) end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); xlabel('Frequency [rad/s]'); ylabel('Phase [deg]'); yticks(-180:90:180); ylim([-180 180]); xticks([1e-2,1e-1,1,1e1]) xticklabels({'$0.01 \omega_0$', '$0.1 \omega_0$', '$\omega_0$', '$10 \omega_0$'}) ax4 = nexttile; hold on; for i = 1:length(Wzs) plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}('dv', 'Fu'), freqs, 'rad/s'))), '-', 'color', colors(i,:)); end hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); xlabel('Frequency [rad/s]'); set(gca, 'YTickLabel',[]); yticks(-180:90:180); ylim([-180 180]); xticks([1e-2,1e-1,1,1e1]) xticklabels({'$0.01 \omega_0$', '$0.1 \omega_0$', '$\omega_0$', '$10 \omega_0$'}) linkaxes([ax1,ax2,ax3,ax4],'x'); xlim([freqs(1), freqs(end)]); linkaxes([ax1,ax2],'y'); % #+name: fig:rotating_rdc_plant_effect_rot % #+caption: Bode plot of the direct and coupling term for the "relative damping control" plant - Effect of rotation % #+RESULTS: % [[file:figs/rotating_rdc_plant_effect_rot.png]] % In order to see if large damping can be added with Relative Damping Control, the root locus is computed (Figure ref:fig:rotating_rdc_root_locus). % The closed-loop system is unconditionally stable and the poles can be damped as much as wanted. %% Root Locus for Relative Damping Control figure; Krdc = s*eye(2); % Relative damping controller gains = logspace(-2, 2, 300); % Tested gains Wz_i = [1,3,4]; hold on; for i = 1:length(Wz_i) plot(real(pole(Gs{Wz_i(i)}({'du', 'dv'}, {'Fu', 'Fv'})*Krdc)), imag(pole(Gs{Wz_i(i)}({'du', 'dv'}, {'Fu', 'Fv'})*Krdc)), 'x', 'color', colors(i,:), ... 'DisplayName', sprintf('$\\Omega = %.1f \\omega_0 $', Wzs(Wz_i(i))),'MarkerSize',8); plot(real(tzero(Gs{Wz_i(i)}({'du', 'dv'}, {'Fu', 'Fv'})*Krdc)), imag(tzero(Gs{Wz_i(i)}({'du', 'dv'}, {'Fu', 'Fv'})*Krdc)), 'o', 'color', colors(i,:), ... 'HandleVisibility', 'off','MarkerSize',8); for g = gains cl_poles = pole(feedback(Gs{Wz_i(i)}({'du', 'dv'}, {'Fu', 'Fv'}), g*Krdc, -1)); plot(real(cl_poles), imag(cl_poles), '.', 'color', colors(i,:), ... 'HandleVisibility', 'off','MarkerSize',4); end end hold off; axis square; xlim([-1.8, 0.2]); ylim([0, 2]); xticks([-1, 0]) xticklabels({'-$\omega_0$', '$0$'}) yticks([0, 1, 2]) yticklabels({'$0$', '$\omega_0$', '$2 \omega_0$'}) xlabel('Real Part'); ylabel('Imaginary Part'); leg = legend('location', 'northwest', 'FontSize', 8); leg.ItemTokenSize(1) = 8; % Damped Plant % Let's select a reasonable "Relative Damping Control" gain, and compute the closed-loop damped system. % The open-loop and damped plants are compared in Figure ref:fig:rotating_rdc_damped_plant. % The rotating aspect does not add any complexity for the use of Relative Damping Control. % It does not increase the low frequency coupling as compared to Integral Force Feedback. i = 2; %% Relative Damping Controller Krdc = 2*s*eye(2); Krdc.InputName = {'Du', 'Dv'}; Krdc.OutputName = {'Fu', 'Fv'}; %% Compute the damped plant G_cl_rdc = feedback(Gs{i}, Krdc, 'name'); %% Damped plant using Relative Damping Control figure; freqs = logspace(-3, 2, 1000); figure; tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); % Magnitude ax1 = nexttile([2, 1]); hold on; plot(freqs, abs(squeeze(freqresp(Gs{i}( 'Du', 'Fu'), freqs, 'rad/s'))), '-', 'color', zeros(1,3), ... 'DisplayName', 'OL - $G_d(1,1)$') plot(freqs, abs(squeeze(freqresp(G_cl_rdc('Du', 'Fu'), freqs, 'rad/s'))), '-', 'color', colors(1,:), ... 'DisplayName', 'RDC - $G_d(1,1)$') plot(freqs, abs(squeeze(freqresp(Gs{i}( 'Dv', 'Fu'), freqs, 'rad/s'))), '-', 'color', [zeros(1,3), 0.5], ... 'DisplayName', 'OL - $G_d(2,1)$') plot(freqs, abs(squeeze(freqresp(G_cl_rdc('Dv', 'Fu'), freqs, 'rad/s'))), '-', 'color', [colors(1,:), 0.5], ... 'DisplayName', 'RDC - $G_d(2,1)$') hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); set(gca, 'XTickLabel',[]); ylabel('Magnitude [m/N]'); ldg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 2); ldg.ItemTokenSize(1) = 20; ylim([1e-6, 1e2]); ax2 = nexttile; hold on; plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}('Du', 'Fu'), freqs, 'rad/s'))), '-', 'color', zeros(1,3)) plot(freqs, 180/pi*angle(squeeze(freqresp(G_cl_rdc('Du', 'Fu'), freqs, 'rad/s'))), '-', 'color', colors(1,:)) hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); xlabel('Frequency [rad/s]'); ylabel('Phase [deg]'); yticks(-180:90:180); ylim([-180 180]); linkaxes([ax1,ax2],'x'); xlim([freqs(1), freqs(end)]);