function main()
clc;
clear;
% Configuration et simulation initiale
[G, K] = configureSystem();
simulateResponses(G, K);
% Analyse supplémentaire pour la robustesse
analyzeSystemRobustness(G, K);
end
function [G, K] = configureSystem()
% Définition des matrices d'état du système
A = [-0.0226, -36.6000, -18.9000, -32.1000;
0, -1.9000, 0.9830, 0;
0.0123, -11.7000, -2.6300, 0;
0, 0, 1.0000, 0];
B = [0, 0, 0, 0;
-0.4140, 0, 0, 0;
-77.8000, 22.4000, 0, 0;
0, 0, 0, 0];
C = [0, 57.3000, 0, 0;
0, 0, 0, 57.3000;
0, 0, 0, 0;
0, 0, 0, 0];
D = zeros(4);
G = ss(A, B, C, D);
Q = C'*C;
R = eye(size(B,2));
K = lqr(A, B, Q, R);
return
end
function simulateResponses(G, K)
% Définition de la plage de fréquences
frequencies = logspace(-2, 6, 500);
fine_freq = logspace(-2, 6, 10000); % Plus de points pour interpolation
% Simulation des données pour M11 et M22
inf_M11 = 0.01 * ones(1, 500);
sup_M11 = logspace(log10(0.01), log10(0.25), 500);
inf_M22 = 0.1 * ones(1, 500);
sup_M22 = logspace(log10(0.1), log10(0.7), 500);
% Interpolation cubique pour obtenir des courbes lisses
sup_M11_interp = interp1(frequencies, sup_M11, fine_freq, 'pchip');
sup_M22_interp = interp1(frequencies, sup_M22, fine_freq, 'pchip');
% Tracé de M11
figure;
semilogx(fine_freq, 0.01 * ones(size(fine_freq)), 'b--', 'DisplayName',
'inf(mu) M11_H');
hold on;
semilogx(fine_freq, sup_M11_interp, 'r-', 'DisplayName', 'sup(mu) M11_H');
title('Performance Nominale du Système Bouclé de M11');
xlabel('Fréquence (rad/s)');
ylabel('M11');
legend('show');
grid on;
% Tracé de M22
figure;
semilogx(fine_freq, 0.1 * ones(size(fine_freq)), 'b--', 'DisplayName',
'inf(mu) M22_H');
hold on;
semilogx(fine_freq, sup_M22_interp, 'r-', 'DisplayName', 'sup(mu) M22_H');
title('Robustesse en Stabilité du Système Bouclé de M22');
xlabel('Fréquence (rad/s)');
ylabel('M22');
legend('show');
grid on;
end
function analyzeSystemRobustness(G, K)
% Définition de la plage de fréquences
omega = logspace(-2, 4, 500);
% Boucle fermée avec contrôleur LQR
CL = feedback(G*K, eye(size(G,1)));
SISO_Channel = CL(1,1); % Sélection du canal SISO pour l'analyse de marge
% Marge de gain et de phase pour le canal sélectionné
figure;
margin(SISO_Channel);
title('Marge de Gain et de Phase pour le Canal SISO Sélectionné');
% Valeurs singulières du système bouclé
sv = sigma(CL, omega);
sv = 20*log10(sv);
% Affichage des valeurs singulières
figure;
semilogx(omega, sv);
title('Valeurs Singulières du Système Bouclé');
xlabel('Fréquence (rad/s)');
ylabel('Magnitude (dB)');
grid on;
% Analyse des réponses en fréquence pour chaque entrée-sortie
figure;
bodemag(CL);
title('Diagramme de Bode du Système Bouclé');
% Réponse à un échelon unitaire pour évaluer la performance temporelle
figure;
step(CL);
title('Réponse à un Échelon du Système Bouclé');
end
% Lancement du script principal
main();
Amelioration code :
function main()
clc;
clear;
% Configuration et simulation initiale
[G, K, W1, W2] = configureSystem();
K_hinf = calculateHinfController(G, W1, W2);
simulateResponses(G, K);
analyzeSystemRobustness(G, K, K_hinf);
end
function [G, K, W1, W2] = configureSystem()
% Définition des matrices d'état du système
A = [-0.0226, -36.6000, -18.9000, -32.1000;
0, -1.9000, 0.9830, 0;
0.0123, -11.7000, -2.6300, 0;
0, 0, 1.0000, 0];
B = [0, 0, 0, 0;
-0.4140, 0, 0, 0;
-77.8000, 22.4000, 0, 0;
0, 0, 0, 0];
C = [0, 57.3000, 0, 0;
0, 0, 0, 57.3000;
0, 0, 0, 0;
0, 0, 0, 0];
D = zeros(4);
G = ss(A, B, C, D);
Q = C'*C;
R = eye(size(B,2));
K = lqr(A, B, Q, R);
% Fonctions de pondération pour la performance et la robustesse
W1 = tf([1 50], [1 10]);
W2 = tf([1 0.5], [1 50]);
W1 = ss(W1);
W2 = ss(W2);
return
end
function K_hinf = calculateHinfController(G, W1, W2)
% Recherche par grille simplifiée pour minimiser la norme H-infini
[A, B, ~, ~] = ssdata(G);
n = size(A,1);
m = size(B,2);
p = size(B,2); % Pour simplifier, considérons m = p
% Recherche par grille sur un ensemble de contrôleurs potentiels
K_test_range = -5:0.5:5;
min_hinf_norm = inf;
K_hinf = zeros(m, p);
for i = K_test_range
for j = K_test_range
K_test = i * eye(m, p) + j * ones(m, p);
CL = feedback(G, ss(K_test));
hinf_norm = norm(CL, Inf);
if hinf_norm < min_hinf_norm
min_hinf_norm = hinf_norm;
K_hinf = K_test;
end
end
end
disp('Optimized H-infinity Controller K:');
disp(K_hinf);
return
end
function simulateResponses(G, K)
% Simulation des réponses avec des données interpolées
frequencies = logspace(-2, 6, 500);
fine_freq = logspace(-2, 6, 10000); % Plus de points pour interpolation
inf_M11 = 0.01 * ones(1, 500);
sup_M11 = logspace(log10(0.01), log10(0.25), 500);
inf_M22 = 0.1 * ones(1, 500);
sup_M22 = logspace(log10(0.1), log10(0.7), 500);
sup_M11_interp = interp1(frequencies, sup_M11, fine_freq, 'pchip');
sup_M22_interp = interp1(frequencies, sup_M22, fine_freq, 'pchip');
figure;
semilogx(fine_freq, 0.01 * ones(size(fine_freq)), 'b--');
hold on;
semilogx(fine_freq, sup_M11_interp, 'r-');
title('Performance Nominale du Système Bouclé de M11');
xlabel('Fréquence (rad/s)');
ylabel('M11');
grid on;
figure;
semilogx(fine_freq, 0.1 * ones(size(fine_freq)), 'b--');
hold on;
semilogx(fine_freq, sup_M22_interp, 'r-');
title('Robustesse en Stabilité du Système Bouclé de M22');
xlabel('Fréquence (rad/s)');
ylabel('M22');
grid on;
end
function analyzeSystemRobustness(G, K, K_hinf)
omega = logspace(-2, 4, 500);
CL = feedback(G*K, eye(size(G,1)));
CL_hinf = feedback(G*K_hinf, eye(size(G,1)));
% Marge de gain et de phase pour le LQR
figure;
margin(CL(1,1));
title('Marge de Gain et de Phase pour le LQR');
% Marge de gain et de phase pour H-infini
figure;
margin(CL_hinf(1,1));
title('Marge de Gain et de Phase pour H-infini');
end
% Lancement du script principal
main();