dehaeze26_closed_loop_shaping/matlab/index.org

Closed-Loop Shaping using Complementary Filters

• Complementary filter design
• Numerical Example
  • Plant
  • Requirements and choice of complementary filters
  • Controller analysis
  • Robustness and Performance analysis

Complementary filter design

Numerical Example

Plant

%% Test model
freqs = logspace(0, 3, 1000); % Frequency Vector [Hz]

m = 20;  % mass [kg]
k = 1e6; % stiffness [N/m]
c = 1e2; % damping [N/(m/s)]

% Plant dynamics
G = 1/(m*s^2 + c*s + k);

% Uncertainty weight
wI = generateWF('n', 2, 'w0', 2*pi*50, 'G0', 0.1, 'Ginf', 10, 'Gc', 1);

Requirements and choice of complementary filters

%% Analytical Complementary Filters
w0 = 2*pi*20;
alpha = 1;

Hh = (s/w0)^2*((s/w0)+1+alpha)/(((s/w0)+1)*((s/w0)^2 + alpha*(s/w0) + 1));
Hl = ((1+alpha)*(s/w0)+1)/(((s/w0)+1)*((s/w0)^2 + alpha*(s/w0) + 1));

Controller analysis

%% Obtained controller
omega = 2*pi*1000;

K = 1/(Hh*G) * 1/((1+s/omega+(s/omega)^2));
K = zpk(minreal(K));

Robustness and Performance analysis

num_delta_points = 50;
theta = linspace(0, 2*pi, num_delta_points);
delta_points = exp(1j * theta);

% Get frequency responses for all components
G_resp  = squeeze(freqresp(G,  freqs, 'Hz'));
K_resp  = squeeze(freqresp(K,  freqs, 'Hz'));
Hl_resp = squeeze(freqresp(Hl, freqs, 'Hz'));
wI_resp = squeeze(freqresp(wI, freqs, 'Hz'));

% Calculate nominal responses
nom_L = G_resp .* K_resp .* Hl_resp;
nom_S = 1 ./ (1 + nom_L);
nom_T = nom_L ./ (1 + nom_L);

% Store all the points in the complex plane that L can take
loop_region_points = zeros(length(freqs), num_delta_points);

% Initialize arrays to store magnitude bounds
S_mag_min = ones(length(freqs), 1) * inf;
S_mag_max = zeros(length(freqs), 1);
T_mag_min = ones(length(freqs), 1) * inf;
T_mag_max = zeros(length(freqs), 1);

% Calculate magnitude bounds for all delta values
for i = 1:num_delta_points
    % Perturbed loop gain
    loop_perturbed = nom_L .* (1 + wI_resp .* delta_points(i));
    loop_region_points(:,i) = loop_perturbed;

    % Perturbed sensitivity function
    S_perturbed = 1 ./ (1 + loop_perturbed);
    S_mag = abs(S_perturbed);

    % Update S magnitude bounds
    S_mag_min = min(S_mag_min, S_mag);
    S_mag_max = max(S_mag_max, S_mag);

    % Perturbed complementary sensitivity function
    T_perturbed = loop_perturbed ./ (1 + loop_perturbed);
    T_mag = abs(T_perturbed);

    % Update T magnitude bounds
    T_mag_min = min(T_mag_min, T_mag);
    T_mag_max = max(T_mag_max, T_mag);
end

% At frequencies where |wI| > 1, T min is zero
T_mag_min(abs(wI_resp)>1) = 1e-10;

%% Nyquist plot to check Robust Stability
figure;
hold on;
plot(real(squeeze(freqresp(G*K*Hl, freqs, 'Hz'))), imag(squeeze(freqresp(G*K*Hl, freqs, 'Hz'))), 'k', 'DisplayName', '$L(j\omega)$ - Nominal');
plot(alphaShape(real(loop_region_points(:)), imag(loop_region_points(:)), 0.1), 'FaceColor', [0, 0, 0], 'EdgeColor', 'none', 'FaceAlpha', 0.3, 'DisplayName', '$L(j\omega)$ - $\forall G \in \Pi_i$');
plot(-1, 0, 'k+', 'MarkerSize', 5, 'HandleVisibility', 'off');
hold off;
grid on;
axis equal
xlim([-1.4, 0.2]); ylim([-1.2, 0.4]);
xticks(-1.4:0.2:0.2); yticks(-1.2:0.2:0.4);
xlabel('Real Part'); ylabel('Imaginary Part');
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);
leg.ItemTokenSize(1) = 18;