This repository has been archived on 2025-04-18. You can view files and clone it, but cannot push or open issues or pull requests.
phd-nass-uniaxial-model/matlab/uniaxial_5_active_damping.m

637 lines
31 KiB
Matlab

%% 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
%% Colors for the figures
colors = colororder;
%% Frequency Vector [Hz]
freqs = logspace(0, 3, 1000);
%% Load the PSD of disturbances
load('uniaxial_disturbance_psd.mat', 'f', 'psd_ft', 'psd_xf');
%% Load Plants Dynamics
load('uniaxial_plants.mat', 'G_vc_light', 'G_md_light', 'G_pz_light', ...
'G_vc_mid', 'G_md_mid', 'G_pz_mid', ...
'G_vc_heavy', 'G_md_heavy', 'G_pz_heavy');
%% Damped plants for three considered payload masses - Comparison of active damping techniques
% Integral Force Feedback
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
ax1 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('fm', 'f'), freqs, 'Hz'))), '-', 'color', colors(1,:), 'DisplayName', '$m_s = 1\,kg$');
plot(freqs, abs(squeeze(freqresp(G_vc_mid( 'fm', 'f'), freqs, 'Hz'))), '-.', 'color', colors(1,:), 'DisplayName', '$m_s = 25\,kg$');
plot(freqs, abs(squeeze(freqresp(G_vc_heavy('fm', 'f'), freqs, 'Hz'))), '--', 'color', colors(1,:), 'DisplayName', '$m_s = 50\,kg$');
plot(freqs, abs(squeeze(freqresp(G_md_light('fm', 'f'), freqs, 'Hz'))), '-', 'color', colors(2,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_md_mid( 'fm', 'f'), freqs, 'Hz'))), '-.', 'color', colors(2,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_md_heavy('fm', 'f'), freqs, 'Hz'))), '--', 'color', colors(2,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_pz_light('fm', 'f'), freqs, 'Hz'))), '-', 'color', colors(3,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_pz_mid( 'fm', 'f'), freqs, 'Hz'))), '-.', 'color', colors(3,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_pz_heavy('fm', 'f'), freqs, 'Hz'))), '--', 'color', colors(3,:), 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);
ldg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);
ldg.ItemTokenSize = [20, 1];
ax1b = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_light('fm', 'f'), freqs, 'Hz')))), '-', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_mid( 'fm', 'f'), freqs, 'Hz')))), '-.', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_heavy('fm', 'f'), freqs, 'Hz')))), '--', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_light('fm', 'f'), freqs, 'Hz')))), '-', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_mid( 'fm', 'f'), freqs, 'Hz')))), '-.', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_heavy('fm', 'f'), freqs, 'Hz')))), '--', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_light('fm', 'f'), freqs, 'Hz')))), '-', 'color', colors(3,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_mid( 'fm', 'f'), freqs, 'Hz')))), '-.', 'color', colors(3,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_heavy('fm', 'f'), freqs, 'Hz')))), '--', 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
xticks([1e0, 1e1, 1e2]);
yticks(-360:90:360);
ylim([-200, 20]);
linkaxes([ax1,ax1b],'x');
xlim([1, 1000]);
% Relative Motion Control
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
ax2 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('dL', 'f'), freqs, 'Hz'))), '-', 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_vc_mid( 'dL', 'f'), freqs, 'Hz'))), '-.', 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_vc_heavy('dL', 'f'), freqs, 'Hz'))), '--', 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_md_light('dL', 'f'), freqs, 'Hz'))), '-', 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_md_mid( 'dL', 'f'), freqs, 'Hz'))), '-.', 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_md_heavy('dL', 'f'), freqs, 'Hz'))), '--', 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_pz_light('dL', 'f'), freqs, 'Hz'))), '-', 'color', colors(3,:));
plot(freqs, abs(squeeze(freqresp(G_pz_mid( 'dL', 'f'), freqs, 'Hz'))), '-.', 'color', colors(3,:));
plot(freqs, abs(squeeze(freqresp(G_pz_heavy('dL', 'f'), freqs, 'Hz'))), '--', 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
ax2b = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_light('dL', 'f'), freqs, 'Hz')))), '-', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_mid( 'dL', 'f'), freqs, 'Hz')))), '-.', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_heavy('dL', 'f'), freqs, 'Hz')))), '--', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_light('dL', 'f'), freqs, 'Hz')))), '-', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_mid( 'dL', 'f'), freqs, 'Hz')))), '-.', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_heavy('dL', 'f'), freqs, 'Hz')))), '--', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_light('dL', 'f'), freqs, 'Hz')))), '-', 'color', colors(3,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_mid( 'dL', 'f'), freqs, 'Hz')))), '-.', 'color', colors(3,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_heavy('dL', 'f'), freqs, 'Hz')))), '--', 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
xticks([1e0, 1e1, 1e2]);
yticks(-360:90:360);
ylim([-200, 20]);
linkaxes([ax2,ax2b],'x');
xlim([1, 1000]);
% Direct Velocity Feedback
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
ax3 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('vn', 'f'), freqs, 'Hz'))), '-', 'color', colors(1,:), 'DisplayName', '$k_n = 0.01\,N/\mu m$');
plot(freqs, abs(squeeze(freqresp(G_vc_mid( 'vn', 'f'), freqs, 'Hz'))), '-.', 'color', colors(1,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_vc_heavy('vn', 'f'), freqs, 'Hz'))), '--', 'color', colors(1,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_md_light('vn', 'f'), freqs, 'Hz'))), '-', 'color', colors(2,:), 'DisplayName', '$k_n = 1\,N/\mu m$');
plot(freqs, abs(squeeze(freqresp(G_md_mid( 'vn', 'f'), freqs, 'Hz'))), '-.', 'color', colors(2,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_md_heavy('vn', 'f'), freqs, 'Hz'))), '--', 'color', colors(2,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_pz_light('vn', 'f'), freqs, 'Hz'))), '-', 'color', colors(3,:), 'DisplayName', '$k_n = 100\,N/\mu m$');
plot(freqs, abs(squeeze(freqresp(G_pz_mid( 'vn', 'f'), freqs, 'Hz'))), '-.', 'color', colors(3,:), 'HandleVisibility', 'off');
plot(freqs, abs(squeeze(freqresp(G_pz_heavy('vn', 'f'), freqs, 'Hz'))), '--', 'color', colors(3,:), 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [m/s/N]'); set(gca, 'XTickLabel',[]);
ldg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);
ldg.ItemTokenSize = [20, 1];
ax3b = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_light('vn', 'f'), freqs, 'Hz')))), '-', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_mid( 'vn', 'f'), freqs, 'Hz')))), '-.', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_heavy('vn', 'f'), freqs, 'Hz')))), '--', 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_light('vn', 'f'), freqs, 'Hz')))), '-', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_mid( 'vn', 'f'), freqs, 'Hz')))), '-.', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_heavy('vn', 'f'), freqs, 'Hz')))), '--', 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_light('vn', 'f'), freqs, 'Hz')))), '-', 'color', colors(3,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_mid( 'vn', 'f'), freqs, 'Hz')))), '-.', 'color', colors(3,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_heavy('vn', 'f'), freqs, 'Hz')))), '--', 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
xticks([1e0, 1e1, 1e2]);
yticks(-360:90:360);
ylim([-110, 110]);
linkaxes([ax3,ax3b],'x');
xlim([1, 1000]);
%% Design of Active Damping controllers to have reasonable damping
% IFF
K_iff_vc = 20/(s + 2*pi*0.01);
K_iff_vc.InputName = {'fm'};
K_iff_vc.OutputName = {'f'};
K_iff_md = 200/(s + 2*pi*0.01);
K_iff_md.InputName = {'fm'};
K_iff_md.OutputName = {'f'};
K_iff_pz = 4000/(s + 2*pi*0.01);
K_iff_pz.InputName = {'fm'};
K_iff_pz.OutputName = {'f'};
% RDC
K_rdc_vc = -1e3*s;
K_rdc_vc.InputName = {'dL'};
K_rdc_vc.OutputName = {'f'};
K_rdc_md = -1e4*s;
K_rdc_md.InputName = {'dL'};
K_rdc_md.OutputName = {'f'};
K_rdc_pz = -1e5*s;
K_rdc_pz.InputName = {'dL'};
K_rdc_pz.OutputName = {'f'};
% DVF
K_dvf_vc = -tf(1e3);
K_dvf_vc.InputName = {'vn'};
K_dvf_vc.OutputName = {'f'};
K_dvf_md = -tf(8e3);
K_dvf_md.InputName = {'vn'};
K_dvf_md.OutputName = {'f'};
K_dvf_pz = -tf(2e5);
K_dvf_pz.InputName = {'vn'};
K_dvf_pz.OutputName = {'f'};
%% Save Active Damping Controller
save('./mat/uniaxial_active_damping_controllers.mat', 'K_iff_vc', 'K_iff_md', 'K_iff_pz', ...
'K_rdc_vc', 'K_rdc_md', 'K_rdc_pz', ...
'K_dvf_vc', 'K_dvf_md', 'K_dvf_pz');
%% Compute Damped Plants
% IFF
G_iff_vc_light = feedback(G_vc_light, K_iff_vc, 'name', +1);
G_iff_vc_mid = feedback(G_vc_mid , K_iff_vc, 'name', +1);
G_iff_vc_heavy = feedback(G_vc_heavy, K_iff_vc, 'name', +1);
G_iff_md_light = feedback(G_md_light, K_iff_md, 'name', +1);
G_iff_md_mid = feedback(G_md_mid , K_iff_md, 'name', +1);
G_iff_md_heavy = feedback(G_md_heavy, K_iff_md, 'name', +1);
G_iff_pz_light = feedback(G_pz_light, K_iff_pz, 'name', +1);
G_iff_pz_mid = feedback(G_pz_mid , K_iff_pz, 'name', +1);
G_iff_pz_heavy = feedback(G_pz_heavy, K_iff_pz, 'name', +1);
% RDC
G_rdc_vc_light = feedback(G_vc_light, K_rdc_vc, 'name', +1);
G_rdc_vc_mid = feedback(G_vc_mid , K_rdc_vc, 'name', +1);
G_rdc_vc_heavy = feedback(G_vc_heavy, K_rdc_vc, 'name', +1);
G_rdc_md_light = feedback(G_md_light, K_rdc_md, 'name', +1);
G_rdc_md_mid = feedback(G_md_mid , K_rdc_md, 'name', +1);
G_rdc_md_heavy = feedback(G_md_heavy, K_rdc_md, 'name', +1);
G_rdc_pz_light = feedback(G_pz_light, K_rdc_pz, 'name', +1);
G_rdc_pz_mid = feedback(G_pz_mid , K_rdc_pz, 'name', +1);
G_rdc_pz_heavy = feedback(G_pz_heavy, K_rdc_pz, 'name', +1);
% DVF
G_dvf_vc_light = feedback(G_vc_light, K_dvf_vc, 'name', +1);
G_dvf_vc_mid = feedback(G_vc_mid , K_dvf_vc, 'name', +1);
G_dvf_vc_heavy = feedback(G_vc_heavy, K_dvf_vc, 'name', +1);
G_dvf_md_light = feedback(G_md_light, K_dvf_md, 'name', +1);
G_dvf_md_mid = feedback(G_md_mid , K_dvf_md, 'name', +1);
G_dvf_md_heavy = feedback(G_md_heavy, K_dvf_md, 'name', +1);
G_dvf_pz_light = feedback(G_pz_light, K_dvf_pz, 'name', +1);
G_dvf_pz_mid = feedback(G_pz_mid , K_dvf_pz, 'name', +1);
G_dvf_pz_heavy = feedback(G_pz_heavy, K_dvf_pz, 'name', +1);
%% Verify Stability
% IFF
if not(isstable(G_iff_vc_light) && isstable(G_iff_vc_mid) && isstable(G_iff_vc_heavy) && ...
isstable(G_iff_md_light) && isstable(G_iff_md_mid) && isstable(G_iff_md_heavy) && ...
isstable(G_iff_pz_light) && isstable(G_iff_pz_mid) && isstable(G_iff_pz_heavy))
warning("One of the damped plant with decentralized IFF is not stable.");
end
% RDC
if not(isstable(G_rdc_vc_light) && isstable(G_rdc_vc_mid) && isstable(G_rdc_vc_heavy) && ...
isstable(G_rdc_md_light) && isstable(G_rdc_md_mid) && isstable(G_rdc_md_heavy) && ...
isstable(G_rdc_pz_light) && isstable(G_rdc_pz_mid) && isstable(G_rdc_pz_heavy))
warning("One of the damped plant with decentralized RDC is not stable.");
end
% DVF
if not(isstable(G_dvf_vc_light) && isstable(G_dvf_vc_mid) && isstable(G_dvf_vc_heavy) && ...
isstable(G_dvf_md_light) && isstable(G_dvf_md_mid) && isstable(G_dvf_md_heavy) && ...
isstable(G_dvf_pz_light) && isstable(G_dvf_pz_mid) && isstable(G_dvf_pz_heavy))
warning("One of the damped plant with decentralized DVF is not stable.");
end
%% Save Damped Plants
save('./mat/uniaxial_damped_plants.mat', 'G_iff_vc_light', 'G_iff_md_light', 'G_iff_pz_light', ...
'G_rdc_vc_light', 'G_rdc_md_light', 'G_rdc_pz_light', ...
'G_dvf_vc_light', 'G_dvf_md_light', 'G_dvf_pz_light', ...
'G_iff_vc_mid', 'G_iff_md_mid', 'G_iff_pz_mid', ...
'G_rdc_vc_mid', 'G_rdc_md_mid', 'G_rdc_pz_mid', ...
'G_dvf_vc_mid', 'G_dvf_md_mid', 'G_dvf_pz_mid', ...
'G_iff_vc_heavy', 'G_iff_md_heavy', 'G_iff_pz_heavy', ...
'G_rdc_vc_heavy', 'G_rdc_md_heavy', 'G_rdc_pz_heavy', ...
'G_dvf_vc_heavy', 'G_dvf_md_heavy', 'G_dvf_pz_heavy');
%% Active Damping Robustness to change of sample's mass - Root Locus for all three damping techniques with 3 different sample's masses
% Soft Nano-Hexapod
figure;
hold on;
% IFF
plot(real(pole(G_vc_light('fm', 'f'))), imag(pole(G_vc_light('fm', 'f'))), 'x', 'color', colors(1,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_vc_light('fm', 'f'))), imag(zero(G_vc_light('fm', 'f'))), 'o', 'color', colors(1,:), ...
'DisplayName', 'IFF');
for g = logspace(0,2,400)
clpoles = pole(feedback(G_vc_light('fm', 'f'), g/s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(1,:), ...
'HandleVisibility', 'off');
end
% RDC
plot(real(pole(G_vc_light('dL', 'f'))), imag(pole(G_vc_light('dL', 'f'))), 'x', 'color', colors(2,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_vc_light('dL', 'f'))), imag(zero(G_vc_light('dL', 'f'))), 'o', 'color', colors(2,:), ...
'DisplayName', 'RDC');
for g = logspace(1,3,400)
clpoles = pole(feedback(G_vc_light('dL', 'f'), -g*s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(2,:), ...
'HandleVisibility', 'off');
end
% DVF
plot(real(pole(G_vc_light('vn', 'f'))), imag(pole(G_vc_light('vn', 'f'))), 'x', 'color', colors(3,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_vc_light('vn', 'f'))), imag(zero(G_vc_light('vn', 'f'))), 'o', 'color', colors(3,:), ...
'DisplayName', 'DVF');
for g = logspace(1,3,400)
clpoles = pole(feedback(G_vc_light('vn', 'f'), -g, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(3,:), ...
'HandleVisibility', 'off');
end
hold off;
axis square;
xlabel('Real Part'); ylabel('Imaginary Part');
ldg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1);
ldg.ItemTokenSize = [10, 1];
xlim([-30, 0]); ylim([0, 30]);
ytickangle(90)
% Medium-Stiff Nano-Hexapod
figure;
hold on;
% IFF
plot(real(pole(G_md_light('fm', 'f'))), imag(pole(G_md_light('fm', 'f'))), 'x', 'color', colors(1,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_md_light('fm', 'f'))), imag(zero(G_md_light('fm', 'f'))), 'o', 'color', colors(1,:), ...
'HandleVisibility', 'off');
for g = logspace(0,3,400)
clpoles = pole(feedback(G_md_light('fm', 'f'), g/s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(1,:), ...
'HandleVisibility', 'off');
end
% RDC
plot(real(pole(G_md_light('dL', 'f'))), imag(pole(G_md_light('dL', 'f'))), 'x', 'color', colors(2,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_md_light('dL', 'f'))), imag(zero(G_md_light('dL', 'f'))), 'o', 'color', colors(2,:), ...
'HandleVisibility', 'off');
for g = logspace(2,4,400)
clpoles = pole(feedback(G_md_light('dL', 'f'), -g*s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(2,:), ...
'HandleVisibility', 'off');
end
% DVF
plot(real(pole(G_md_light('vn', 'f'))), imag(pole(G_md_light('vn', 'f'))), 'x', 'color', colors(3,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_md_light('vn', 'f'))), imag(zero(G_md_light('vn', 'f'))), 'o', 'color', colors(3,:), ...
'HandleVisibility', 'off');
for g = logspace(2,4,400)
clpoles = pole(feedback(G_md_light('vn', 'f'), -g, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(3,:), ...
'HandleVisibility', 'off');
end
hold off;
axis square;
xlabel('Real Part'); ylabel('Imaginary Part');
xlim([-300, 0]); ylim([0, 300]);
ytickangle(90)
% Stiff Nano Hexapod
figure;
hold on;
% IFF
plot(real(pole(G_pz_light('fm', 'f'))), imag(pole(G_pz_light('fm', 'f'))), 'x', 'color', colors(1,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_pz_light('fm', 'f'))), imag(zero(G_pz_light('fm', 'f'))), 'o', 'color', colors(1,:), ...
'HandleVisibility', 'off');
for g = logspace(2,5,400)
clpoles = pole(feedback(G_pz_light('fm', 'f'), g/s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(1,:), ...
'HandleVisibility', 'off');
end
% RDC
plot(real(pole(G_pz_light('dL', 'f'))), imag(pole(G_pz_light('dL', 'f'))), 'x', 'color', colors(2,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_pz_light('dL', 'f'))), imag(zero(G_pz_light('dL', 'f'))), 'o', 'color', colors(2,:), ...
'HandleVisibility', 'off');
for g = logspace(3,6,400)
clpoles = pole(feedback(G_pz_light('dL', 'f'), -g*s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(2,:), ...
'HandleVisibility', 'off');
end
% DVF
plot(real(pole(G_pz_light('vn', 'f'))), imag(pole(G_pz_light('vn', 'f'))), 'x', 'color', colors(3,:), ...
'HandleVisibility', 'off');
plot(real(zero(G_pz_light('vn', 'f'))), imag(zero(G_pz_light('vn', 'f'))), 'o', 'color', colors(3,:), ...
'HandleVisibility', 'off');
for g = logspace(3,6,400)
clpoles = pole(feedback(G_pz_light('vn', 'f'), -g, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(3,:), ...
'HandleVisibility', 'off');
end
hold off;
axis square;
xlabel('Real Part'); ylabel('Imaginary Part');
xlim([-4000, 0]); ylim([0, 4000]);
ytickangle(90)
%% Root Locus for the three damping techniques
figure;
hold on;
% IFF
plot(real(pole(G_md_mid('fm', 'f'))), imag(pole(G_md_mid('fm', 'f'))), 'x', 'color', colors(1,:), ...
'DisplayName', 'IFF');
plot(real(zero(G_md_mid('fm', 'f'))), imag(zero(G_md_mid('fm', 'f'))), 'o', 'color', colors(1,:), ...
'HandleVisibility', 'off');
for g = logspace(1,4,500)
clpoles = pole(feedback(G_md_mid('fm', 'f'), g/s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(1,:), ...
'HandleVisibility', 'off');
end
% RDC
plot(real(pole(G_md_mid('dL', 'f'))), imag(pole(G_md_mid('dL', 'f'))), 'x', 'color', colors(2,:), ...
'DisplayName', 'RDC');
plot(real(zero(G_md_mid('dL', 'f'))), imag(zero(G_md_mid('dL', 'f'))), 'o', 'color', colors(2,:), ...
'HandleVisibility', 'off');
% Estimate the maximum damping added by RDC
gs = logspace(2,5,500);
phis = zeros(size(gs));
for i = 1:length(gs)
g = gs(i);
clpoles = pole(feedback(G_md_mid('dL', 'f'), -g*s, +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(2,:), ...
'HandleVisibility', 'off');
% Estimate damping of u-station mode
ustation_pole = clpoles(imag(clpoles)>1000);
phis(i) = atan2(abs(real(ustation_pole)), abs(imag(ustation_pole)));
end
[~, i_max] = max(phis);
plot([0, -5e3*sin(phis(i_max))], [0, 5e3*cos(phis(i_max))], 'k--', 'HandleVisibility', 'off');
clpoles_max = pole(feedback(G_md_mid('dL', 'f'), -gs(i_max)*s, +1));
ustation_pole = clpoles_max(imag(clpoles_max)>1000);
plot(real(ustation_pole), imag(ustation_pole), 'kx', ...
'HandleVisibility', 'off');
% Plot angle
plot(-8e2*sin(0:0.01:max(phis)), 8e2*cos(sin(0:0.01:max(phis))), 'k-', 'HandleVisibility', 'off')
text(-200, 850, '$\phi$', 'horizontalalignment', 'center');
text(real(ustation_pole)-100, imag(ustation_pole), '$\xi = \sin(\phi)$', 'horizontalalignment', 'right');
% DVF
plot(real(pole(G_md_mid('vn', 'f'))), imag(pole(G_md_mid('vn', 'f'))), 'x', 'color', colors(3,:), ...
'DisplayName', 'DVF');
plot(real(zero(G_md_mid('vn', 'f'))), imag(zero(G_md_mid('vn', 'f'))), 'o', 'color', colors(3,:), ...
'HandleVisibility', 'off');
for g = logspace(2,5,500)
clpoles = pole(feedback(G_md_mid('vn', 'f'), -tf(g), +1));
plot(real(clpoles), imag(clpoles), '.', 'color', colors(3,:), ...
'HandleVisibility', 'off');
end
hold off;
xlim([-2100, 0]); ylim([0, 2100]);
axis square;
xlabel('Real Part'); ylabel('Imaginary Part');
ldg = legend('location', 'northwest', 'FontSize', 8, 'NumColumns', 1);
ldg.ItemTokenSize = [10, 1];
%% Obtained damped transfer function from f to d for the three damping techniques
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
ax1 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_mid('d', 'f'), freqs, 'Hz'))), 'k-', 'DisplayName', 'OL');
plot(freqs, abs(squeeze(freqresp(G_iff_vc_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'IFF');
plot(freqs, abs(squeeze(freqresp(G_rdc_vc_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'RDC');
plot(freqs, abs(squeeze(freqresp(G_dvf_vc_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f$ [m/N]'); set(gca, 'XTickLabel',[]);
ax2 = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_mid('d', 'f'), freqs, 'Hz')))), 'k-');
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_vc_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_rdc_vc_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_dvf_vc_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
yticks(-360:90:360);
ylim([-270, 90]);
xticks([1e0, 1e1, 1e2]);
linkaxes([ax1,ax2],'x');
xlim([1, 500]);
%% Obtained damped transfer function from f to d for the three damping techniques
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
ax1 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid('d', 'f'), freqs, 'Hz'))), 'k-', 'DisplayName', 'OL');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'IFF');
plot(freqs, abs(squeeze(freqresp(G_rdc_md_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'RDC');
plot(freqs, abs(squeeze(freqresp(G_dvf_md_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f$ [m/N]'); set(gca, 'XTickLabel',[]);
ax2 = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_mid('d', 'f'), freqs, 'Hz')))), 'k-');
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_md_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_rdc_md_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_dvf_md_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
yticks(-360:90:360);
ylim([-270, 90]);
xticks([1e0, 1e1, 1e2]);
linkaxes([ax1,ax2],'x');
xlim([1, 500]);
%% Obtained damped transfer function from f to d for the three damping techniques
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
ax1 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_pz_mid('d', 'f'), freqs, 'Hz'))), 'k-', 'DisplayName', 'OL');
plot(freqs, abs(squeeze(freqresp(G_iff_pz_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'IFF');
plot(freqs, abs(squeeze(freqresp(G_rdc_pz_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'RDC');
plot(freqs, abs(squeeze(freqresp(G_dvf_pz_mid('d', 'f'), freqs, 'Hz'))), 'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f$ [m/N]'); set(gca, 'XTickLabel',[]);
leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
leg.ItemTokenSize(1) = 15;
ax2 = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_mid('d', 'f'), freqs, 'Hz')))), 'k-');
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_pz_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_rdc_pz_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_dvf_pz_mid('d', 'f'), freqs, 'Hz')))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
yticks(-360:90:360);
ylim([-270, 90]);
xticks([1e0, 1e1, 1e2]);
linkaxes([ax1,ax2],'x');
xlim([1, 500]);
%% Change of sensitivity to disturbance with all three active damping strategies
% FS
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid('d', 'fs'), freqs, 'Hz'))), 'k-');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'fs'), freqs, 'Hz'))), 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_rdc_md_mid('d', 'fs'), freqs, 'Hz'))), 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_dvf_md_mid('d', 'fs'), freqs, 'Hz'))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f_{s}$ [m/N]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);
xlim([1, 500]);
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid('d', 'ft'), freqs, 'Hz'))), 'k-');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'ft'), freqs, 'Hz'))), 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_rdc_md_mid('d', 'ft'), freqs, 'Hz'))), 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_dvf_md_mid('d', 'ft'), freqs, 'Hz'))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f_{t}$ [m/N]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);
xlim([1, 500]);
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid('d', 'xf'), freqs, 'Hz'))), 'k-', 'DisplayName', 'OL');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'xf'), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'IFF');
plot(freqs, abs(squeeze(freqresp(G_rdc_md_mid('d', 'xf'), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'RDC');
plot(freqs, abs(squeeze(freqresp(G_dvf_md_mid('d', 'xf'), freqs, 'Hz'))), 'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/x_{f}$ [m/m]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);
xlim([1, 500]);
%% Cumulative Amplitude Spectrum of the distance d with all three active damping techniques
figure;
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', 'black', 'DisplayName', 'OL');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_iff_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(1,:), 'DisplayName', 'IFF');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_rdc_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_rdc_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(2,:), 'DisplayName', 'RDC');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_dvf_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_dvf_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('CAS of $d$ [m]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);
xlim([1, 500]);
ylim([2e-10, 3e-6])
figure;
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', 'black', 'DisplayName', 'OL');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_iff_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(1,:), 'DisplayName', 'IFF');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_rdc_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_rdc_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(2,:), 'DisplayName', 'RDC');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_dvf_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_dvf_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
xticks([1e0, 1e1, 1e2]);
xlim([1, 500]);
ylim([2e-10, 3e-6])
figure;
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', 'black', 'DisplayName', 'OL');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_iff_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(1,:), 'DisplayName', 'IFF');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_rdc_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_rdc_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(2,:), 'DisplayName', 'RDC');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_dvf_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
psd_xf.*abs(squeeze(freqresp(G_dvf_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
'color', colors(3,:), 'DisplayName', 'DVF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
xticks([1e0, 1e1, 1e2]);
legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
xlim([1, 500]);
ylim([2e-10, 3e-6])