From f041b63ebfdbfd9c0b2e26623c5540fffc0f7096 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 15:51:29 +0330 Subject: [PATCH 01/22] =?UTF-8?q?Create=20Base=20model=20for=20robust=20H?= =?UTF-8?q?=E2=88=9E=20control?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...se model for robust H\342\210\236 control" | 165 ++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 "robust_hinfinity_control/Base model for robust H\342\210\236 control" diff --git "a/robust_hinfinity_control/Base model for robust H\342\210\236 control" "b/robust_hinfinity_control/Base model for robust H\342\210\236 control" new file mode 100644 index 0000000..119f0b5 --- /dev/null +++ "b/robust_hinfinity_control/Base model for robust H\342\210\236 control" @@ -0,0 +1,165 @@ +%% ========================================================================= +% final.m - Base model for robust H∞ control +% ========================================================================= +% This script serves as the fundamental implementation of H∞ control design. +% It includes: +% - Plant modeling +% - Controller synthesis using `hinfsyn` +% - Robust performance analysis using `mu()` +% Sensitivity weighting functions are applied to improve robustness. +% +% System Equations: +% 10 +% Plant = ------------------ * pade +% (s - 1)(s + 1) +% +% 2.4 s^2 +% Wt = ------------ +% ( s + 1 )^2 +% ========================================================================= + +%% 1. Define system parameters +tau = 1.5 / 2; % Delay parameter +w = logspace(-2, 2, 200); % Frequency range + +%% 2. Plant modeling +p0 = 10 ./ (-w.^2 - 1); +pm = p0 .* (-tau * 1j * w + 1) ./ (tau * 1j * w + 1); +p1 = exp(-2 * tau * 1j * w) .* p0; +W1 = abs(p1 ./ p0 - 1); +W2 = abs(p1 ./ pm - 1); + +%% 3. Sensitivity weighting functions +n1 = 2.4 * tau * [1 0]; +d1 = [tau 1]; +[m1, ~] = bode(n1, d1, w); + +tauf = 0.52; +n2 = 0.32 * [1 0 0 0]; +d2 = conv([tauf^2 2 * tauf * 0.7 1], [tauf 1]); +[m2, ~] = bode(n2, d2, w); + +% Uncomment to visualize frequency response +% semilogx(w, 20 * log10(W1), ':', w, 20 * log10(W2), '--', ... +% w, 20 * log10(m1), 'm', w, 20 * log10(m2), 'k--'), grid +% axis([0.01 100 -40 20]) + +%%%% ========================================================================= +% Define plant models +%%%% ========================================================================= +numerator_plant = 10 * [-tau 1]; +denominator_plant = conv([1 0 -1], [tau 1]); + +numerator_wt = n2; +denominator_wt = d2; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% Define interconnection structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic_st'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic_st); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic_st = pck(A, B, C, D); + +%%%% ========================================================================= +% H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic_st, 1, 1, 0.1, 5, 1e-4, 2); + +%%%% ========================================================================= +% Mixed Sensitivity Design +%%%% ========================================================================= +numerator_ws = [1 0.5]; +denominator_ws = 100 * [1 0.005]; + +ws = nd2sys(numerator_ws, denominator_ws); +wu = nd2sys(1e-6, 1); + +systemnames = 'plant ws wt wu'; +inputvar = '[yd; u]'; +outputvar = '[wt; ws; wu ; yd-plant]'; +input_to_plant = '[ u ]'; +input_to_ws = '[ yd-plant ]'; +input_to_wt = '[ plant ]'; +input_to_wu = '[ u ]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +% Nominal performance configuration +systemnames = 'plant ws wu'; +inputvar = '[yd; u]'; +outputvar = '[ ws; wu ; yd-plant]'; +input_to_plant = '[ u ]'; +input_to_ws = '[ yd-plant ]'; +input_to_wu = '[ u ]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 10, 0.001, 2); + +%%%% ========================================================================= +% Robust Performance Analysis using μ (Mu) Analysis +%%%% ========================================================================= +systemnames = 'plant wt ws'; +inputvar = '[pout; yd; u]'; +outputvar = '[ u; ws; plant-yd ]'; +input_to_plant = '[ u + wt ]'; +input_to_wt = '[ pout ]'; +input_to_ws = '[ yd - plant ]'; +sysoutname = 'himat_ic'; +cleanupsysic = 'yes'; +sysic + +[k1, clp1] = hinfsyn(himat_ic, 1, 1, 0.8, 50, 1e-1, 2); +omega = logspace(-3, 3, 60); +clp1g = frsp(clp1, omega); +clp1gs = vsvd(clp1g); % Singular value calculation + +% Plot singular values +vplot('liv,m', clp1gs), grid +title('Singular Value Plot for the Closed-Loop System') +xlabel('Frequency (rad/sec)') + +%%%% ========================================================================= +% D-K Iteration for Robust Controller Optimization +%%%% ========================================================================= +deltaset = [1 1; 1 1]; +clp_g1 = clp1g; +[bnds1, pv1] = mu(clp_g1, deltaset); + +vplot('liv,m', vnorm(clp_g1), bnds1); grid +title('Maximum Singular Value and Mu Plot'); +xlabel('Frequency (rad/sec)') + +% Perform D-K iteration +DK_DEF_NAME = 'himat_def'; % Define DK iteration parameters +dkit + +% Controller order reduction +[k_dk3bal, hsv] = sysbal(k_dk3); +[k_dk3red] = strunc(k_dk3bal, 13); +clpred_13 = starp(himat_ic, k_dk3red); +clpred_13g = frsp(clpred_13, OMEGA_DK); +[bnds] = mu(clpred_13g, [1 1; 1 1], 'c'); + +pkvnorm(sel(bnds, 1, 1)) From 9b1240f7f9372e563b911cbe713b3fccfbbe5daa Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:12:50 +0330 Subject: [PATCH 02/22] =?UTF-8?q?Rename=20Base=20model=20for=20robust=20H?= =?UTF-8?q?=E2=88=9E=20control=20to=20Base=20model=20for=20robust=20H?= =?UTF-8?q?=E2=88=9E=20control?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Base model for robust H\342\210\236 control" | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename "robust_hinfinity_control/Base model for robust H\342\210\236 control" => "ROBUST H\342\210\236 CONTROL DESIGN/Base model for robust H\342\210\236 control" (100%) diff --git "a/robust_hinfinity_control/Base model for robust H\342\210\236 control" "b/ROBUST H\342\210\236 CONTROL DESIGN/Base model for robust H\342\210\236 control" similarity index 100% rename from "robust_hinfinity_control/Base model for robust H\342\210\236 control" rename to "ROBUST H\342\210\236 CONTROL DESIGN/Base model for robust H\342\210\236 control" From 0d50448e37c71f0051293749aeb404cd3e24cc83 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:16:33 +0330 Subject: [PATCH 03/22] =?UTF-8?q?Rename=20Base=20model=20for=20robust=20H?= =?UTF-8?q?=E2=88=9E=20control=20to=20Base=20model=20for=20robust=20H?= =?UTF-8?q?=E2=88=9E=20control?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Base model for robust H\342\210\236 control" | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename "ROBUST H\342\210\236 CONTROL DESIGN/Base model for robust H\342\210\236 control" => "Robust H\342\210\236 Control Design/Base model for robust H\342\210\236 control" (100%) diff --git "a/ROBUST H\342\210\236 CONTROL DESIGN/Base model for robust H\342\210\236 control" "b/Robust H\342\210\236 Control Design/Base model for robust H\342\210\236 control" similarity index 100% rename from "ROBUST H\342\210\236 CONTROL DESIGN/Base model for robust H\342\210\236 control" rename to "Robust H\342\210\236 Control Design/Base model for robust H\342\210\236 control" From 26f3beac72eb1347a67a202c9423a44b269c7107 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:20:27 +0330 Subject: [PATCH 04/22] Create Modified Sensitivity Weighting Parameters --- ...Modified Sensitivity Weighting Parameters" | 149 ++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 "Robust H\342\210\236 Control Design/Modified Sensitivity Weighting Parameters" diff --git "a/Robust H\342\210\236 Control Design/Modified Sensitivity Weighting Parameters" "b/Robust H\342\210\236 Control Design/Modified Sensitivity Weighting Parameters" new file mode 100644 index 0000000..87962fc --- /dev/null +++ "b/Robust H\342\210\236 Control Design/Modified Sensitivity Weighting Parameters" @@ -0,0 +1,149 @@ +%% ========================================================================= +% final_80_1.m - Modified Sensitivity Weighting Parameters +% ========================================================================= +% This script is similar to `final.m`, but it modifies the weighting functions +% W_s and W_t to optimize performance under different frequency ranges. +% Robust performance is evaluated using `mu()` analysis. +% +% System Equations: +% 10 +% Plant = ------------------ * pade +% (s - 1)(s + 1) +% +% 2.4 s^2 +% Wt = ------------ +% ( s + 1 )^2 +% ========================================================================= + +%% 1. Define System Parameters +clear all +tau = 1 / 2; % Delay parameter +w = logspace(-2, 2, 200); % Frequency range + +%% 2. Plant Modeling +p0 = 10 ./ (-w.^2 - 1); +pm = p0 .* (-tau * 1j * w + 1) ./ (tau * 1j * w + 1); +p1 = exp(-2 * tau * 1j * w) .* p0; +W1 = abs(p1 ./ p0 - 1); +W2 = abs(p1 ./ pm - 1); + +%% 3. Sensitivity Weighting Functions +n1 = 2.4 * tau * [1 0]; +d1 = [tau 1]; +[m1, ~] = bode(n1, d1, w); + +tauf = 0.34; +n2 = 0.09 * [1 0 0 0]; +d2 = conv([tauf^2 2 * tauf * 0.7 1], [tauf 1]); +[m2, ~] = bode(n2, d2, w); + +%%%% ========================================================================= +% Define Plant Models +%%%% ========================================================================= +numerator_plant = 10 * [-tau 1]; +denominator_plant = conv([1 0 -1], [tau 1]); + +numerator_wt = n2; +denominator_wt = d2; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic_st'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic_st); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic_st = pck(A, B, C, D); + +%%%% ========================================================================= +% H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic_st, 1, 1, 0.1, 5, 1e-4, 2); +disp('---------->>>>>>>>>>> Robust Stability'); + +%%%% ========================================================================= +% Mixed Sensitivity Design +%%%% ========================================================================= +numerator_ws = [1 0.8]; +denominator_ws = 25 * [1 0.008]; + +ws = nd2sys(numerator_ws, denominator_ws); +wu = nd2sys(1e-6, 1); + +systemnames = 'plant ws wt wu'; +inputvar = '[yd; u]'; +outputvar = '[wt; ws; wu ; yd-plant]'; +input_to_plant = '[ u ]'; +input_to_ws = '[ yd-plant ]'; +input_to_wt = '[ plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 10, 0.001, 2); +disp('---------->>>>>>>>>>> Mixed Sensitivity'); + +%%%% ========================================================================= +% Robust Performance Analysis using μ (Mu) Analysis +%%%% ========================================================================= +systemnames = 'plant wt ws'; +inputvar = '[pout; yd; u]'; +outputvar = '[ u; ws; plant-yd ]'; +input_to_plant = '[ u + wt ]'; +input_to_wt = '[ pout ]'; +input_to_ws = '[ yd - plant]'; +sysoutname = 'himat_ic'; +cleanupsysic = 'yes'; +sysic + +[k1, clp1] = hinfsyn(himat_ic, 1, 1, 0.8, 50, 1e-1, 2); +omega = logspace(-3, 3, 60); +clp1g = frsp(clp1, omega); +clp1gs = vsvd(clp1g); % Singular value calculation + +% Plot singular values +vplot('liv,m', clp1gs), grid +title('Singular Value Plot for the Closed-Loop System') +xlabel('Frequency (rad/sec)') + +%%%% ========================================================================= +% D-K Iteration for Robust Controller Optimization +%%%% ========================================================================= +deltaset = [1 1; 1 1]; +clp_g1 = clp1g; +[bnds1, pv1] = mu(clp_g1, deltaset); + +vplot('liv,m', vnorm(clp_g1), bnds1); grid +title('Maximum Singular Value and Mu Plot'); +xlabel('Frequency (rad/sec)') + +% Perform D-K iteration +DK_DEF_NAME = 'himat_def'; % Define DK iteration parameters +dkit + +% Controller Order Reduction +[k_dk5bal, hsv] = sysbal(k_dk5); +[k_dk5red] = strunc(k_dk5bal, 7); +clpred_7 = starp(himat_ic, k_dk5red); +clpred_7g = frsp(clpred_7, OMEGA_DK); +[bnds] = mu(clpred_7g, [1 1; 1 1], 'c'); + +pkvnorm(sel(bnds, 1, 1)) From 7be1ab4381db3c8fee0e04c68f39d9a7a6cfe6fe Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:21:44 +0330 Subject: [PATCH 05/22] Create Alternative Model with Parameter Variation --- ...lternative Model with Parameter Variation" | 149 ++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 "Robust H\342\210\236 Control Design/Alternative Model with Parameter Variation" diff --git "a/Robust H\342\210\236 Control Design/Alternative Model with Parameter Variation" "b/Robust H\342\210\236 Control Design/Alternative Model with Parameter Variation" new file mode 100644 index 0000000..fcab40f --- /dev/null +++ "b/Robust H\342\210\236 Control Design/Alternative Model with Parameter Variation" @@ -0,0 +1,149 @@ +%% ========================================================================= +% final_80_2.m - Alternative Model with Parameter Variation +% ========================================================================= +% This script investigates how changes in the time delay (τ) affect system performance. +% It tests the stability and robustness of the H∞ controller under different scenarios. +% The core control design remains the same, but different plant dynamics are introduced. +% +% System Equations: +% 10 +% Plant = ------------------ * pade +% (s - 1)(s + 1) +% +% 2.4 s^2 +% Wt = ------------ +% ( s + 1 )^2 +% ========================================================================= + +%% 1. Define System Parameters +clear all +tau = 1 / 2; % Delay parameter +w = logspace(-2, 2, 200); % Frequency range + +%% 2. Plant Modeling +p0 = 10 ./ (-w.^2 - 1); +pm = p0 .* (-tau * 1j * w + 1) ./ (tau * 1j * w + 1); +p1 = exp(-2 * tau * 1j * w) .* p0; +W1 = abs(p1 ./ p0 - 1); +W2 = abs(p1 ./ pm - 1); + +%% 3. Sensitivity Weighting Functions +n1 = 2.4 * tau * [1 0]; +d1 = [tau 1]; +[m1, ~] = bode(n1, d1, w); + +tauf = 0.34; +n2 = 0.09 * [1 0 0 0]; +d2 = conv([tauf^2 2 * tauf * 0.7 1], [tauf 1]); +[m2, ~] = bode(n2, d2, w); + +%%%% ========================================================================= +% Define Plant Models +%%%% ========================================================================= +numerator_plant = 10 * [-tau 1]; +denominator_plant = conv([1 0 -1], [tau 1]); + +numerator_wt = n1; +denominator_wt = d1; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic_st'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic_st); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic_st = pck(A, B, C, D); + +%%%% ========================================================================= +% H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic_st, 1, 1, 0.1, 5, 1e-4, 2); +disp('<<<<<<<<<<<<<<<<<<<< Robust Stability >>>>>>>>>>>>>>>>>>>>>'); + +%%%% ========================================================================= +% Mixed Sensitivity Design +%%%% ========================================================================= +numerator_ws = [1 0.5]; +denominator_ws = 30 * [1 0.005]; + +ws = nd2sys(numerator_ws, denominator_ws); +wu = nd2sys(1e-6, 1); + +systemnames = 'plant ws wt wu'; +inputvar = '[yd; u]'; +outputvar = '[wt; ws; wu ; yd-plant]'; +input_to_plant = '[ u ]'; +input_to_ws = '[ yd-plant ]'; +input_to_wt = '[ plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 50, 0.001, 2); +disp('<<<<<<<<<<<<<<<<<<<< Mixed Sensitivity >>>>>>>>>>>>>>>>>>>>>'); + +%%%% ========================================================================= +% Robust Performance Analysis using μ (Mu) Analysis +%%%% ========================================================================= +systemnames = 'plant wt ws'; +inputvar = '[pout; yd; u]'; +outputvar = '[ u; ws; plant-yd ]'; +input_to_plant = '[ u + wt ]'; +input_to_wt = '[ pout ]'; +input_to_ws = '[ yd - plant]'; +sysoutname = 'himat_ic'; +cleanupsysic = 'yes'; +sysic + +[k1, clp1] = hinfsyn(himat_ic, 1, 1, 0.8, 5, 1e-1, 2); +omega = logspace(-3, 3, 60); +clp1g = frsp(clp1, omega); +clp1gs = vsvd(clp1g); % Singular value calculation + +% Plot singular values +vplot('liv,m', clp1gs), grid +title('Singular Value Plot for the Closed-Loop System') +xlabel('Frequency (rad/sec)') + +%%%% ========================================================================= +% D-K Iteration for Robust Controller Optimization +%%%% ========================================================================= +deltaset = [1 1; 1 1]; +clp_g1 = clp1g; +[bnds1, pv1] = mu(clp_g1, deltaset); + +vplot('liv,m', vnorm(clp_g1), bnds1); grid +title('Maximum Singular Value and Mu Plot'); +xlabel('Frequency (rad/sec)') + +% Perform D-K iteration +DK_DEF_NAME = 'himat_def'; % Define DK iteration parameters +dkit + +% Controller Order Reduction +[k_dk5bal, hsv] = sysbal(k_dk5); +[k_dk5red] = strunc(k_dk5bal, 7); +clpred_7 = starp(himat_ic, k_dk5red); +clpred_7g = frsp(clpred_7, OMEGA_DK); +[bnds] = mu(clpred_7g, [1 1; 1 1], 'c'); + +pkvnorm(sel(bnds, 1, 1)) From 9ab536e40345f711a12763f9f4716a8030727735 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:29:01 +0330 Subject: [PATCH 06/22] =?UTF-8?q?Create=20Alternative=20Weighting=20for=20?= =?UTF-8?q?Robust=20H=E2=88=9E=20Control?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...eighting for Robust H\342\210\236 Control" | 130 ++++++++++++++++++ 1 file changed, 130 insertions(+) create mode 100644 "Robust H\342\210\236 Control Design/Alternative Weighting for Robust H\342\210\236 Control" diff --git "a/Robust H\342\210\236 Control Design/Alternative Weighting for Robust H\342\210\236 Control" "b/Robust H\342\210\236 Control Design/Alternative Weighting for Robust H\342\210\236 Control" new file mode 100644 index 0000000..c8a63bd --- /dev/null +++ "b/Robust H\342\210\236 Control Design/Alternative Weighting for Robust H\342\210\236 Control" @@ -0,0 +1,130 @@ +%% ========================================================================= +% finex6.m - Alternative Weighting for Robust H∞ Control +% ========================================================================= +% This script uses a slightly different system structure while maintaining +% the robust H∞ control design. The plant transfer function is defined differently, +% resulting in a variation of control performance. +% The `hinfsyn` and `mu()` analysis remain the same for robust performance evaluation. +% +% System Equations: +% 10 +% Plant = ------------------ * pade +% (s - 1)(s + 1) +% +% 2.4 s +% Wt = ------------ +% ( s + 4 ) +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +% Time delay parameter +tau = 1.5; +tau = tau / 2; + +% Frequency range for analysis +w = logspace(-2, 2, 200); + +%%%% ========================================================================= +% 2. Define Plant Transfer Function and Weighting Functions +%%%% ========================================================================= + +% Plant Transfer Function: +% 10 (1.5s - 1) +% P(s) = ----------------- +% (s - 1)(s + 1)(0.25s + 1) +numerator_plant = 10 * [1.5 -1]; +denominator_plant = conv([1 0 -1], [0.25 1]); + +% Weighting function W_t: +% 2.4 s +% W_t = ------------ +% (s + 4) +numerator_wt = [2.4 0]; +denominator_wt = [1 4]; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% 3. Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +%%%% ========================================================================= +% 4. System Information and Minimal Realization +%%%% ========================================================================= +minfo(plant_ic) +spoles(plant_ic) + +pause % Press any key to continue +clc + +% Ensure minimal realization to avoid uncontrollable hidden modes +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 5. H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 8, 1e-4, 2); + +pause % Press any key to continue + +%%%% ========================================================================= +% 6. Bode Plot of Controller +%%%% ========================================================================= +clf +omega = logspace(-2, 6, 300); +k1g = frsp(k1, omega); +title('Bode Plot of Controller, k1') +xlabel('Frequency (rad/sec)') +vplot('bode', k1g) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 7. Frequency Response and Peak Norm Analysis +%%%% ========================================================================= +omega_2 = logspace(0, 5, 50); +g1g = frsp(g1, omega_2); + +% Find peak norm of the frequency response +g1pk = pkvnorm(g1g); +fprintf('\n Final Gamma: %g, \n Max Singular Value: %g \n\n', gf1, g1pk) + +%%%% ========================================================================= +% 8. Controller Transfer Function +%%%% ========================================================================= +g1gs = vsvd(g1g); +seesys(k1) + +% Convert state-space representation to transfer function +[ak, bk, ck, dk] = unpck(k1); +[numk, denk] = ss2tf(ak, bk, ck, dk); +[zk, pk, kk] = tf2zp(numk, denk); + +% Print transfer function +printsys(numk, denk, 's') + +%%%% ========================================================================= +% 9. End of Script +%%%% ========================================================================= +disp('This concludes the example.') From 8a2bf1ca6f5e457b36abcb7d4b4eff2094071f66 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:31:07 +0330 Subject: [PATCH 07/22] Create Comprehensive Robustness Testing --- .../Comprehensive Robustness Testing" | 158 ++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 "Robust H\342\210\236 Control Design/Comprehensive Robustness Testing" diff --git "a/Robust H\342\210\236 Control Design/Comprehensive Robustness Testing" "b/Robust H\342\210\236 Control Design/Comprehensive Robustness Testing" new file mode 100644 index 0000000..5c773e8 --- /dev/null +++ "b/Robust H\342\210\236 Control Design/Comprehensive Robustness Testing" @@ -0,0 +1,158 @@ +%% ========================================================================= +% test.m - Comprehensive Robustness Testing +% ========================================================================= +% This script serves as a test framework to compare the robustness of different +% controller configurations. +% - It validates the designed controller by simulating multiple system variations. +% - It performs `D-K Iteration` to reduce controller order while maintaining robustness. +% +% System Equations: +% 10 +% Plant = ------------------ * pade +% (s - 1)(s + 1) +% +% 2.4 s^2 +% Wt = ------------ +% ( s + 1 )^2 +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +tau = 1 / 2; % Delay parameter +w = logspace(-2, 2, 200); % Frequency range + +%%%% ========================================================================= +% 2. Plant Transfer Function and Weighting Functions +%%%% ========================================================================= +p0 = 10 ./ (-w.^2 - 1); +pm = p0 .* (-tau * 1j * w + 1) ./ (tau * 1j * w + 1); +p1 = exp(-2 * tau * 1j * w) .* p0; +W1 = abs(p1 ./ p0 - 1); +W2 = abs(p1 ./ pm - 1); + +% Sensitivity weighting functions +n1 = 2.4 * tau * [1 0]; +d1 = [tau 1]; + +tauf = 0.34; +n2 = 0.09 * [1 0 0 0]; +d2 = conv([tauf^2 2 * tauf * 0.7 1], [tauf 1]); + +%%%% ========================================================================= +% 3. Define Plant Models +%%%% ========================================================================= +numerator_plant = 10 * [-tau 1]; +denominator_plant = conv([1 0 -1], [tau 1]); + +numerator_wt = n2; +denominator_wt = d2; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% 4. Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic_st'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic_st); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic_st = pck(A, B, C, D); + +%%%% ========================================================================= +% 5. H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic_st, 1, 1, 0.1, 5, 1e-4, 2); +disp('---------->>>>>>>>>>> Robust Stability'); + +%%%% ========================================================================= +% 6. Mixed Sensitivity Design +%%%% ========================================================================= +numerator_ws = [1 5]; +denominator_ws = 50 * [1 0.05]; + +ws = nd2sys(numerator_ws, denominator_ws); +wu = nd2sys(1e-6, 1); + +systemnames = 'plant ws wt wu'; +inputvar = '[yd; u]'; +outputvar = '[wt; ws; wu ; yd-plant]'; +input_to_plant = '[ u ]'; +input_to_ws = '[ yd-plant ]'; +input_to_wt = '[ plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 10, 0.001, 2); +disp('---------->>>>>>>>>>> Mixed Sensitivity'); + +%%%% ========================================================================= +% 7. Robust Performance Analysis using μ (Mu) Analysis +%%%% ========================================================================= +systemnames = 'plant wt ws'; +inputvar = '[pout; yd; u]'; +outputvar = '[ u; ws; plant-yd ]'; +input_to_plant = '[ u + wt ]'; +input_to_wt = '[ pout ]'; +input_to_ws = '[ yd - plant]'; +sysoutname = 'himat_ic'; +cleanupsysic = 'yes'; +sysic + +% Compute closed-loop system +clp1 = starp(himat_ic, k1); + +omega = logspace(-3, 3, 60); +clp1g = frsp(clp1, omega); +clp1gs = vsvd(clp1g); % Singular value calculation + +% Plot singular values +vplot('liv,m', clp1gs), grid +title('Singular Value Plot for the Closed-Loop System') +xlabel('Frequency (rad/sec)') + +%%%% ========================================================================= +% 8. Robustness Performance Analysis +%%%% ========================================================================= +deltaset = [1 1; 1 1]; +clp_g1 = clp1g; +[bnds1, pv1] = mu(clp_g1, deltaset); + +vplot('liv,m', vnorm(clp_g1), bnds1); grid +title('Maximum Singular Value and Mu Plot'); +xlabel('Frequency (rad/sec)') + +%%%% ========================================================================= +% 9. D-K Iteration for Robust Controller Optimization +%%%% ========================================================================= +DK_DEF_NAME = 'himat_def1'; % Define DK iteration parameters +dkit + +%%%% ========================================================================= +% 10. Controller Order Reduction +%%%% ========================================================================= +[k_dk3bal, hsv] = sysbal(k_dk3); +[k_dk3red] = strunc(k_dk3bal, 13); +clpred_13 = starp(himat_ic, k_dk3red); +clpred_13g = frsp(clpred_13, OMEGA_DK); +[bnds] = mu(clpred_13g, [1 1; 1 1], 'c'); + +pkvnorm(sel(bnds, 1, 1)) From 6869379c05968f4bdb6ea5c59a53bc51ec4d4299 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:45:24 +0330 Subject: [PATCH 08/22] Create Fundamental Delay Analysis --- .../Fundamental Delay Analysis | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 Delay Analysis Systems/Fundamental Delay Analysis diff --git a/Delay Analysis Systems/Fundamental Delay Analysis b/Delay Analysis Systems/Fundamental Delay Analysis new file mode 100644 index 0000000..c1f1148 --- /dev/null +++ b/Delay Analysis Systems/Fundamental Delay Analysis @@ -0,0 +1,45 @@ +%% ========================================================================= +% delay.m - Fundamental Analysis of Time Delay Effects +% ========================================================================= +% This script analyzes the influence of different time delays on the +% frequency response of a dynamical system. The effect of varying delay +% values is visualized using Bode plots and logarithmic magnitude plots. +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +tau1 = 0.1; +tau2 = 0.2; +w = logspace(-1, 2, 200); % Frequency range + +%%%% ========================================================================= +% 2. Define Plant Transfer Function with Delays +%%%% ========================================================================= +p0 = 1 ./ (1j * w) ./ (1j * w + 1); +p1 = exp(-tau1 * 1j * w) .* p0; +W1 = abs(p0 ./ p1 - 1); +p2 = exp(-tau2 * 1j * w) .* p0; +W2 = abs(p0 ./ p2 - 1); + +%%%% ========================================================================= +% 3. Define Sensitivity Weighting Function +%%%% ========================================================================= +tau = 1 / 10; +numerator = 2.4 * tau * [1 0]; +denominator = [tau 1]; + +% Compute Bode Magnitude Response +[mag, ph] = bode(numerator, denominator, w); + +%%%% ========================================================================= +% 4. Plot Frequency Response +%%%% ========================================================================= +semilogx(w, 20 * log10(W1), ':', ... + w, 20 * log10(W2), '--', ... + w, 20 * log10(mag)), grid +title('Effect of Time Delay on Frequency Response') +xlabel('Frequency (rad/sec)') +ylabel('Magnitude (dB)') +legend('W1 (tau1=0.1)', 'W2 (tau2=0.2)', 'Bode Magnitude Response') From d2320ba65a5134d68eef08201511c0cb7f7cdc5f Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:48:27 +0330 Subject: [PATCH 09/22] Create Frequency Response Analysis with Delay Variation --- ...ncy Response Analysis with Delay Variation | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 Delay Analysis Systems/Frequency Response Analysis with Delay Variation diff --git a/Delay Analysis Systems/Frequency Response Analysis with Delay Variation b/Delay Analysis Systems/Frequency Response Analysis with Delay Variation new file mode 100644 index 0000000..0353a93 --- /dev/null +++ b/Delay Analysis Systems/Frequency Response Analysis with Delay Variation @@ -0,0 +1,45 @@ +%% ========================================================================= +% finex1.m - Frequency Response Analysis with Delay Variation +% ========================================================================= +% This script extends delay analysis by testing multiple delay values +% (τ₁, τ₂) and observing their impact on the system's frequency +% response. It uses logarithmic plots to highlight differences. +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +tau1 = 3; % First delay value +tau2 = 3; % Second delay value +w = logspace(-2, 1, 200); % Frequency range + +%%%% ========================================================================= +% 2. Define Plant Transfer Function with Delays +%%%% ========================================================================= +p0 = 10 ./ (-w.^2 - 1); +p1 = exp(-tau1 * 1j * w) .* p0; +W1 = abs(p0 ./ p1 - 1); +p2 = exp(-tau2 * 1j * w) .* p0; +W2 = abs(p0 ./ p2 - 1); + +%%%% ========================================================================= +% 3. Define Sensitivity Weighting Function +%%%% ========================================================================= +tau = 1.5; +numerator = 2.4 * tau * [1 0]; +denominator = [tau 1]; + +% Compute Bode Magnitude Response +[mag, ph] = bode(numerator, denominator, w); + +%%%% ========================================================================= +% 4. Plot Frequency Response +%%%% ========================================================================= +semilogx(w, 20 * log10(W1), ':', ... + w, 20 * log10(W2), '--', ... + w, 20 * log10(mag)), grid +title('Effect of Delay Variation on Frequency Response') +xlabel('Frequency (rad/sec)') +ylabel('Magnitude (dB)') +legend('W1 (τ₁=3)', 'W2 (τ₂=3)', 'Bode Magnitude Response') From 563e224b5766fb874ffdf3755114dce58768b9c9 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:51:58 +0330 Subject: [PATCH 10/22] Create Expanded Delay Influence with New Weighting --- ...xpanded Delay Influence with New Weighting | 57 +++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 Delay Analysis Systems/Expanded Delay Influence with New Weighting diff --git a/Delay Analysis Systems/Expanded Delay Influence with New Weighting b/Delay Analysis Systems/Expanded Delay Influence with New Weighting new file mode 100644 index 0000000..faa899b --- /dev/null +++ b/Delay Analysis Systems/Expanded Delay Influence with New Weighting @@ -0,0 +1,57 @@ +%% ========================================================================= +% finex3.m - Expanded Delay Influence with Alternative Weighting +% ========================================================================= +% This script builds on `finex1.m` by incorporating alternative weightings +% to assess how different weighting factors affect the impact of delay on +% system stability and performance. +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +tau = 1.8; % Delay parameter +w = logspace(-2, 2, 200); % Frequency range + +%%%% ========================================================================= +% 2. Define Plant Transfer Function with Delays +%%%% ========================================================================= +p0 = 10 ./ (-w.^2 - 1); +pm = p0 .* (-tau / 2 * 1j * w + 1) ./ (tau / 2 * 1j * w + 1); +p1 = exp(-tau * 1j * w) .* p0; +W1 = abs(p1 ./ p0 - 1); +W2 = abs(p1 ./ pm - 1); + +%%%% ========================================================================= +% 3. Define Alternative Weighting Functions +%%%% ========================================================================= +tau = 0.9; +numerator_1 = 2.4 * tau * [1 0]; +denominator_1 = [tau 1]; + +% Compute Bode Magnitude Response for first weighting function +[m1, p1] = bode(numerator_1, denominator_1, w); + +% Alternative weighting +t1 = 0.75; +t2 = 0.6; + +numerator_2 = 0.65 * [1 0 0 0]; +denominator_2 = [t2^3 3*t2^2 3*t2 1]; + +% Compute Bode Magnitude Response for second weighting function +[m2, p2] = bode(numerator_2, denominator_2, w); + +%%%% ========================================================================= +% 4. Plot Frequency Response +%%%% ========================================================================= +semilogx(w, 20 * log10(W1), ':', ... + w, 20 * log10(W2), '--', ... + w, 20 * log10(m1), 'm', ... + w, 20 * log10(m2), 'k--'), grid + +title('Effect of Delay with Alternative Weighting') +xlabel('Frequency (rad/sec)') +ylabel('Magnitude (dB)') +legend('W1 (Delay effect)', 'W2 (Modified response)', ... + 'Bode Magnitude (Weight 1)', 'Bode Magnitude (Weight 2)') From efb4b6357f9b0461053b2cbb366b83ef4771fe83 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:54:44 +0330 Subject: [PATCH 11/22] Create Sensitivity Analysis of Delay Impact --- .../Sensitivity Analysis of Delay Impact | 46 +++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 Delay Analysis Systems/Sensitivity Analysis of Delay Impact diff --git a/Delay Analysis Systems/Sensitivity Analysis of Delay Impact b/Delay Analysis Systems/Sensitivity Analysis of Delay Impact new file mode 100644 index 0000000..1c5449e --- /dev/null +++ b/Delay Analysis Systems/Sensitivity Analysis of Delay Impact @@ -0,0 +1,46 @@ +%% ========================================================================= +% finex5.m - Sensitivity Analysis of Delay Impact +% ========================================================================= +% This script studies how minor variations in time delay influence the +% magnitude response of a system. It provides insights into the stability +% sensitivity of delayed systems. +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +tau1 = 1; % First delay value +tau2 = 1; % Second delay value +w = logspace(-1, 2, 200); % Frequency range + +%%%% ========================================================================= +% 2. Define Plant Transfer Function with Delays +%%%% ========================================================================= +p0 = 10 ./ (-w.^2 - 1); +p1 = exp(-tau1 * 1j * w) .* p0; +W1 = abs(p0 ./ p1 - 1); +p2 = exp(-tau2 * 1j * w) .* p0; +W2 = abs(p0 ./ p2 - 1); + +%%%% ========================================================================= +% 3. Define Sensitivity Weighting Function +%%%% ========================================================================= +tau = 0.25; +numerator = 2.4 * tau * [1 0]; +denominator = [tau 1]; + +% Compute Bode Magnitude Response +[mag, ph] = bode(numerator, denominator, w); + +%%%% ========================================================================= +% 4. Plot Frequency Response +%%%% ========================================================================= +semilogx(w, 20 * log10(W1), ':', ... + w, 20 * log10(W2), '--', ... + w, 20 * log10(mag)), grid + +title('Sensitivity Analysis of Delay Impact') +xlabel('Frequency (rad/sec)') +ylabel('Magnitude (dB)') +legend('W1 (τ₁=1)', 'W2 (τ₂=1)', 'Bode Magnitude Response') From 81eeb2f79e4ae3216c1904a76a7bd17ad9c76d11 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 17:41:57 +0330 Subject: [PATCH 12/22] =?UTF-8?q?Create=20H=E2=88=9E=20Control=20with=20Al?= =?UTF-8?q?ternative=20Weighting=20and=20Plant=20Structure?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...Alternative Weighting and Plant Structure" | 122 ++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100644 "Robust H\342\210\236 Control Design/H\342\210\236 Control with Alternative Weighting and Plant Structure" diff --git "a/Robust H\342\210\236 Control Design/H\342\210\236 Control with Alternative Weighting and Plant Structure" "b/Robust H\342\210\236 Control Design/H\342\210\236 Control with Alternative Weighting and Plant Structure" new file mode 100644 index 0000000..7426405 --- /dev/null +++ "b/Robust H\342\210\236 Control Design/H\342\210\236 Control with Alternative Weighting and Plant Structure" @@ -0,0 +1,122 @@ +%% ========================================================================= +% finex2.m - H∞ Control Design for Robust Stability +% ========================================================================= +% This script implements H∞ control synthesis for a robust stability example. +% - Defines the plant model and weighting functions. +% - Designs an H∞ controller using `hinfsyn`. +% - Evaluates performance via Bode plots and singular value analysis. +% +% System Equations: +% 10 +% Plant = ------------------ +% (s - 1)(s + 1) +% +% 3.6 s +% Wt = ------------ +% 1.5 s + 1 +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +%%%% ========================================================================= +% 2. Define Plant Transfer Function and Weighting Functions +%%%% ========================================================================= +% Plant Transfer Function: +% 10 +% P(s) = ----------------- +% (s + 1)(s - 1) +numerator_plant = [10]; +denominator_plant = conv([1 1], [1 -1]); + +% Weighting function W_t: +% 3.6 s +% W_t = ------------ +% 1.5 s + 1 +numerator_wt = 3 * [1 0.25]; +denominator_wt = [1 1.5]; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% 3. Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +%%%% ========================================================================= +% 4. System Information and Minimal Realization +%%%% ========================================================================= +minfo(plant_ic) +spoles(plant_ic) + +pause % Press any key to continue +clc + +% Ensure minimal realization to avoid uncontrollable hidden modes +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 5. H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 8, 1e-4, 2); + +pause % Press any key to continue + +%%%% ========================================================================= +% 6. Bode Plot of Controller +%%%% ========================================================================= +clf +omega = logspace(-2, 6, 300); +k1g = frsp(k1, omega); +title('Bode Plot of Controller, k1') +xlabel('Frequency (rad/sec)') +vplot('bode', k1g) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 7. Frequency Response and Peak Norm Analysis +%%%% ========================================================================= +omega_2 = logspace(0, 5, 50); +g1g = frsp(g1, omega_2); + +% Find peak norm of the frequency response +g1pk = pkvnorm(g1g); +fprintf('\n Final Gamma: %g, \n Max Singular Value: %g \n\n', gf1, g1pk) + +%%%% ========================================================================= +% 8. Controller Transfer Function +%%%% ========================================================================= +g1gs = vsvd(g1g); +seesys(k1) + +% Convert state-space representation to transfer function +[ak, bk, ck, dk] = unpck(k1); +[numk, denk] = ss2tf(ak, bk, ck, dk); +[zk, pk, kk] = tf2zp(numk, denk); + +% Print transfer function +printsys(numk, denk, 's') + +%%%% ========================================================================= +% 9. End of Script +%%%% ========================================================================= +disp('This concludes the example.') From 8503beeb7295f62269c2caba78c5a840a032298a Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Mon, 17 Mar 2025 17:48:06 +0330 Subject: [PATCH 13/22] Create Robust Stability Evaluation --- .../Robust Stability Evaluation | 171 ++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 Robustness Analysis Optimization/Robust Stability Evaluation diff --git a/Robustness Analysis Optimization/Robust Stability Evaluation b/Robustness Analysis Optimization/Robust Stability Evaluation new file mode 100644 index 0000000..479d520 --- /dev/null +++ b/Robustness Analysis Optimization/Robust Stability Evaluation @@ -0,0 +1,171 @@ +%% ========================================================================= +% robust_st.m - Robust Stability Analysis with H∞ Control +% ========================================================================= +% This script evaluates the robust stability of a system using H∞ control. +% - Constructs an augmented plant model. +% - Performs frequency response analysis. +% - Designs an H∞ controller with `hinfsyn` and evaluates stability. +% +% System Equations: +% 1 +% Plant = ------------------ +% (s - 1)(s - 2) +% +% (s + 1)^2 +% Wt = ------------ +% 25 +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +%%%% ========================================================================= +% 2. Define Plant Transfer Function and Weighting Functions +%%%% ========================================================================= +% Plant Transfer Function: +% 1 +% P(s) = ----------------- +% (s - 1)(s - 2) +numerator_plant = [1]; +denominator_plant = conv([1 -1], [1 -2]); + +% Weighting function W_t: +% (s + 1)^2 +% W_t = ------------ +% 25 +numerator_wt = [1 2 1]; +denominator_wt = 25; + +% Augmented plant transfer function: P(s) * Wt(s) +numerator_aug = numerator_wt; +denominator_aug = denominator_wt * denominator_plant; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +aug_pl = nd2sys(numerator_aug, denominator_aug); + +%%%% ========================================================================= +% 3. Frequency Response Analysis +%%%% ========================================================================= +omega = logspace(-1, 2, 200); % Frequency range +wplantg = frsp(plant, omega); % Frequency response of the plant + +% Plot the plant frequency response +clf; +subplot(111) +vplot('liv,lm', wplantg) +title('Plant Frequency Response') +xlabel('Frequency (rad/s)') +ylabel('Magnitude') + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 4. Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'aug_pl wt'; +inputvar = '[yref; u]'; +outputvar = '[aug_pl; yref - wt]'; +input_to_aug_pl = '[ u ]'; +input_to_wt = '[ aug_pl]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +%%%% ========================================================================= +% 5. System Information and Minimal Realization +%%%% ========================================================================= +minfo(plant_ic) +spoles(plant_ic) + +pause % Press any key to continue +clc + +% Ensure minimal realization to avoid uncontrollable hidden modes +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 6. H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 2, 1e-4, 2); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 7. Controller Analysis +%%%% ========================================================================= +% Check poles of the designed controller +rifd(spoles(k1)) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 8. Bode Plot of Controller +%%%% ========================================================================= +clf +omega_3 = logspace(-2, 6, 300); +k1g = frsp(k1, omega_3); +title('Bode Plot of Controller, k1') +xlabel('Frequency (rad/sec)') +vplot('bode', k1g) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 9. Closed-Loop Stability Check +%%%% ========================================================================= +rifd(spoles(g1)) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 10. Closed-Loop Frequency Response +%%%% ========================================================================= +omega_2 = logspace(0, 5, 50); +g1g = frsp(g1, omega_2); + +% Find peak norm of the frequency response +g1pk = pkvnorm(g1g); +fprintf('\n Final Gamma: %g, \n Max Singular Value: %g \n\n', gf1, g1pk) + +% Singular value plot of the closed-loop system +g1gs = vsvd(g1g); +clf +vplot('liv,m', g1g, g1gs) +title('Singular Value and Sigma Plot of the H∞ Closed Loop') +xlabel('Frequency (rad/s)') +ylabel('Magnitude') +grid; + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 11. Controller Transfer Function +%%%% ========================================================================= +seesys(k1) + +% Convert state-space representation to transfer function +[ak, bk, ck, dk] = unpck(k1); +[numk, denk] = ss2tf(ak, bk, ck, dk); +[zk, pk, kk] = tf2zp(numk, denk); + +% Print transfer function +printsys(numk, denk, 's') + +%%%% ========================================================================= +% 12. End of Script +%%%% ========================================================================= +disp('This concludes the example.') From 55987125afc174f8f3df7d2976978f9d81cb40e2 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 04:57:01 +0330 Subject: [PATCH 14/22] Create Alternative Robust Stability Model --- .../Alternative Robust Stability Model | 122 ++++++++++++++++++ 1 file changed, 122 insertions(+) create mode 100644 Robustness Analysis Optimization/Alternative Robust Stability Model diff --git a/Robustness Analysis Optimization/Alternative Robust Stability Model b/Robustness Analysis Optimization/Alternative Robust Stability Model new file mode 100644 index 0000000..db5513d --- /dev/null +++ b/Robustness Analysis Optimization/Alternative Robust Stability Model @@ -0,0 +1,122 @@ +%% ========================================================================= +% robust_st2.m - Alternative Robust Stability Model +% ========================================================================= +% This script explores an alternative system configuration for robust stability. +% - Uses a different plant model while maintaining the H∞ control structure. +% - Evaluates system robustness using frequency response analysis. +% - Compares results with `robust_st.m`. +% +% System Equations: +% ( s - 1 ) +% Plant = ------------------ +% (s + 1)(s - 0.2) +% +% (s + 0.1) +% Wt = ------------ +% (s + 1) +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +%%%% ========================================================================= +% 2. Define Plant Transfer Function and Weighting Functions +%%%% ========================================================================= +% Plant Transfer Function: +% (s - 1) +% P(s) = ----------------- +% (s + 1)(s - 0.2) +numerator_plant = [1 -1]; +denominator_plant = conv([1 1], [1 -0.2]); + +% Weighting function W_t: +% (s + 0.1) +% W_t = ------------ +% (s + 1) +numerator_wt = [1 0.1]; +denominator_wt = [1 1]; + +% Convert systems to Mu-synthesis format +plant = nd2sys(numerator_plant, denominator_plant); +wt = nd2sys(numerator_wt, denominator_wt); +wu = nd2sys(1e-6, 1); + +%%%% ========================================================================= +% 3. Define Interconnection Structure +%%%% ========================================================================= +systemnames = 'plant wt wu'; +inputvar = '[yref; u]'; +outputvar = '[wt; wu; yref-plant]'; +input_to_plant = '[ u ]'; +input_to_wt = '[plant]'; +input_to_wu = '[u]'; +sysoutname = 'plant_ic'; +cleanupsysic = 'yes'; +sysic + +%%%% ========================================================================= +% 4. System Information and Minimal Realization +%%%% ========================================================================= +minfo(plant_ic) +spoles(plant_ic) + +pause % Press any key to continue +clc + +% Ensure minimal realization to avoid uncontrollable hidden modes +[A, B, C, D] = unpck(plant_ic); +[A, B, C, D] = minreal(A, B, C, D); +plant_ic = pck(A, B, C, D); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 5. H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 8, 1e-4, 2); + +pause % Press any key to continue + +%%%% ========================================================================= +% 6. Bode Plot of Controller +%%%% ========================================================================= +clf +omega = logspace(-2, 6, 300); +k1g = frsp(k1, omega); +title('Bode Plot of Controller, k1') +xlabel('Frequency (rad/sec)') +vplot('bode', k1g) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 7. Frequency Response and Peak Norm Analysis +%%%% ========================================================================= +omega_2 = logspace(0, 5, 50); +g1g = frsp(g1, omega_2); + +% Find peak norm of the frequency response +g1pk = pkvnorm(g1g); +fprintf('\n Final Gamma: %g, \n Max Singular Value: %g \n\n', gf1, g1pk) + +%%%% ========================================================================= +% 8. Controller Transfer Function +%%%% ========================================================================= +g1gs = vsvd(g1g); +seesys(k1) + +% Convert state-space representation to transfer function +[ak, bk, ck, dk] = unpck(k1); +[numk, denk] = ss2tf(ak, bk, ck, dk); +[zk, pk, kk] = tf2zp(numk, denk); + +% Print transfer function +printsys(numk, denk, 's') + +%%%% ========================================================================= +% 9. End of Script +%%%% ========================================================================= +disp('This concludes the example.') From 65e2e2eef9d6f93affaad33561000ea009bca1ea Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:01:51 +0330 Subject: [PATCH 15/22] Update Robust Stability Evaluation --- .../Robust Stability Evaluation | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Robustness Analysis Optimization/Robust Stability Evaluation b/Robustness Analysis Optimization/Robust Stability Evaluation index 479d520..3db1a6c 100644 --- a/Robustness Analysis Optimization/Robust Stability Evaluation +++ b/Robustness Analysis Optimization/Robust Stability Evaluation @@ -24,16 +24,10 @@ clc % 2. Define Plant Transfer Function and Weighting Functions %%%% ========================================================================= % Plant Transfer Function: -% 1 -% P(s) = ----------------- -% (s - 1)(s - 2) numerator_plant = [1]; denominator_plant = conv([1 -1], [1 -2]); % Weighting function W_t: -% (s + 1)^2 -% W_t = ------------ -% 25 numerator_wt = [1 2 1]; denominator_wt = 25; From d492c1db0e7760bea4e98e904bbcce6831f3ae28 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:02:29 +0330 Subject: [PATCH 16/22] Update Alternative Robust Stability Model --- .../Alternative Robust Stability Model | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Robustness Analysis Optimization/Alternative Robust Stability Model b/Robustness Analysis Optimization/Alternative Robust Stability Model index db5513d..cea8156 100644 --- a/Robustness Analysis Optimization/Alternative Robust Stability Model +++ b/Robustness Analysis Optimization/Alternative Robust Stability Model @@ -24,16 +24,10 @@ clc % 2. Define Plant Transfer Function and Weighting Functions %%%% ========================================================================= % Plant Transfer Function: -% (s - 1) -% P(s) = ----------------- -% (s + 1)(s - 0.2) numerator_plant = [1 -1]; denominator_plant = conv([1 1], [1 -0.2]); % Weighting function W_t: -% (s + 0.1) -% W_t = ------------ -% (s + 1) numerator_wt = [1 0.1]; denominator_wt = [1 1]; From e84c8ad63692053db80e9eeace50f387105f42b8 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:06:24 +0330 Subject: [PATCH 17/22] Create Robustness Evaluation with Mu-Analysis --- .../Robustness Evaluation with Mu-Analysis | 108 ++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 Robustness Analysis Optimization/Robustness Evaluation with Mu-Analysis diff --git a/Robustness Analysis Optimization/Robustness Evaluation with Mu-Analysis b/Robustness Analysis Optimization/Robustness Evaluation with Mu-Analysis new file mode 100644 index 0000000..b577aa4 --- /dev/null +++ b/Robustness Analysis Optimization/Robustness Evaluation with Mu-Analysis @@ -0,0 +1,108 @@ +%% ========================================================================= +% mu_analyse.m - Robustness Evaluation with Mu-Analysis +% ========================================================================= +% This script evaluates the robustness of a control system using µ-analysis. +% - Constructs system interconnections. +% - Uses `mu()` to compute robustness bounds. +% - Compares performance under uncertain conditions. +% +% System Equations: +% +% 1 +% G(s) = ----- +% (s - 1) +% +% Uncertainty Weighting Function: +% 0.5s + 1 +% Wu(s) = -------- +% 0.03125s + 1 +% +% Performance Weighting Function: +% 0.25s + 0.6 +% Wp(s) = ------------ +% s + 0.006 +% +% Controller 1: +% 0.9s + 1 +% K1(s) = -------- +% s +% +% Controller 2: +% 2.8s + 1 +% K2(s) = -------- +% s +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +%%%% ========================================================================= +% 2. Define Plant and Weighting Functions +%%%% ========================================================================= +G = nd2sys(1, [1 -1]); % Plant transfer function + +Wu = nd2sys([0.5 1], [0.03125 1], 0.25); % Uncertainty weighting function + +Wp = nd2sys([0.25 0.6], [1 0.006]); % Performance weighting function + +%%%% ========================================================================= +% 3. Define System Interconnections +%%%% ========================================================================= +systemnames = 'G Wu Wp'; +sysoutname = 'P'; +inputvar = '[ z; d; u ]'; +input_to_G = '[ z + u ]'; +input_to_Wu = '[ u ]'; +input_to_Wp = '[ d + G ]'; +outputvar = '[ Wu; Wp; G + d ]'; +cleanupsysic = 'yes'; +sysic + +%%%% ========================================================================= +% 4. Define Controllers +%%%% ========================================================================= +K1 = nd2sys([0.9 1], [1 0], -10); % Controller 1 + +K2 = nd2sys([2.8 1], [1 0], -1); % Controller 2 + +%%%% ========================================================================= +% 5. Compute Closed-Loop Systems +%%%% ========================================================================= +omega = logspace(-2, 2, 80); +M1 = starp(P, K1); +M2 = starp(P, K2); +M1g = frsp(M1, omega); +M2g = frsp(M2, omega); + +%%%% ========================================================================= +% 6. Perform Nominal Performance Analysis +%%%% ========================================================================= +figure; +subplot(1,2,1) +vplot('liv,m', vnorm(sel(M1g,2,2)), '-', ... + vnorm(sel(M2g,2,2)), '--', ... + vnorm(sel(M1g,1,1)), '.', 1, '-.') +axis([0.1 1000 0 2]) +xlabel('Frequency (rad/sec)') +title('Nominal Performance') +grid on; + +pause(1) + +%%%% ========================================================================= +% 7. Perform Robust Performance Analysis +%%%% ========================================================================= +uncblk = [1 1]; +fictblk = [1 1]; +deltaset = [uncblk; fictblk]; + +bnds1 = mu(M1g, deltaset); +bnds2 = mu(M2g, deltaset); + +subplot(1,2,2) +vplot('liv,m', bnds1, '-', bnds2, '--', 1, '-.') +axis([0.1 1000 0 2]) +xlabel('Frequency (rad/sec)') +title('Robust Performance') +grid on; From 2800dab659bce79c694f4b7d8dd315ecf614df4f Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:08:49 +0330 Subject: [PATCH 18/22] Create Mixed Sensitivity Optimization for Robust Control --- ...ensitivity Optimization for Robust Control | 121 ++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100644 Robustness Analysis Optimization/Mixed Sensitivity Optimization for Robust Control diff --git a/Robustness Analysis Optimization/Mixed Sensitivity Optimization for Robust Control b/Robustness Analysis Optimization/Mixed Sensitivity Optimization for Robust Control new file mode 100644 index 0000000..84a1ba7 --- /dev/null +++ b/Robustness Analysis Optimization/Mixed Sensitivity Optimization for Robust Control @@ -0,0 +1,121 @@ +%% ========================================================================= +% mixed_sen1.m - Mixed Sensitivity Optimization for Robust Control +% ========================================================================= +% This script designs a robust control system using a combination of: +% - H∞ control to maximize robustness. +% - H2 control to minimize control effort. +% - Uses `hinfsyn` and `h2syn` for controller design. +% +% System Equations: +% +% 1 +% Plant = --------------- +% s^2 +0.2s + 1 +% +% 10(s + 1) +% ws = ----------------- +% (s+0.01)(s + 10) +% +% 2*(s+10)^2 +% wt = ------------ +% 1000 +% +% wu = 1 +% ========================================================================= + +%% 1. Define System Parameters +clear all +clc + +%%%% ========================================================================= +% 2. Define Plant State-Space Model +%%%% ========================================================================= +An = [0 1 0 0; + -1 -.2 0 0; + 0 0 0 1; + -1 0 -0.1 -10.01]; + +Bn = [0 0; + 0 1; + 0 0; + 1 0]; + +Cn = [ 0 0 0.5 0.5; + 0 0 0 0; + .198 .0369 0 0; + -1 0 0 0]; + +Dn = [0 0; + 0 1; + 0 2e-3; + 1 0]; + +% Convert state-space representation to system format +plant_ic = pck(An, Bn, Cn, Dn); + +%%%% ========================================================================= +% 3. System Information and Minimal Realization +%%%% ========================================================================= +minfo(plant_ic) +spoles(plant_ic) + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 4. H∞ Controller Synthesis +%%%% ========================================================================= +[k1, g1, gf1] = hinfsyn(plant_ic, 1, 1, 0.1, 2, 1e-4, 2); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 5. H2 Controller Synthesis +%%%% ========================================================================= +% Syntax: [k,g,norms] = h2syn(plant, nmeas, ncon, ricmethod, quiet) +[k2, g2, n2] = h2syn(plant_ic, 1, 1, 2, 1); +norm2 = norm(n2, 2); + +disp('The 2-norm of objective function is:'); +disp(norm2); + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 6. Bode Plot of Controllers +%%%% ========================================================================= +clf +omega = logspace(-2, 6, 200); +k1_g = frsp(k1, omega); +k2_g = frsp(k2, omega); +vplot('bode', k1_g, k2_g), grid + +title('Magnitude Plot of Controllers: H∞ and H2') +xlabel('Frequency (rad/sec)') + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 7. Controller Transfer Functions +%%%% ========================================================================= +% H∞ Controller +[ak, bk, ck, dk] = unpck(k1); +[numk1, denk1] = ss2tf(ak, bk, ck, dk); +[z1, p1, k1] = ss2zp(ak, bk, ck, dk); + +disp('H∞ Controller Transfer Function:') +printsys(numk1, denk1, 's') + +% H2 Controller +[ak2, bk2, ck2, dk2] = unpck(k2); +[numk2, denk2] = ss2tf(ak2, bk2, ck2, dk2); +[z2, p2, k2] = ss2zp(ak2, bk2, ck2, dk2); + +disp('H2 Controller Transfer Function:') +printsys(numk2, denk2, 's') + +pause % Press any key to continue +clc From 9407214eebf0353c53381b24b8443802f29d7abd Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:18:42 +0330 Subject: [PATCH 19/22] =?UTF-8?q?Create=20H=E2=88=9E=20Control=20for=20Spa?= =?UTF-8?q?ce=20Shuttle=20Flight=20Dynamics?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...Control for Space Shuttle Flight Dynamics" | 171 ++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 "Aerospace Robust Control/H\342\210\236 Control for Space Shuttle Flight Dynamics" diff --git "a/Aerospace Robust Control/H\342\210\236 Control for Space Shuttle Flight Dynamics" "b/Aerospace Robust Control/H\342\210\236 Control for Space Shuttle Flight Dynamics" new file mode 100644 index 0000000..b15cdad --- /dev/null +++ "b/Aerospace Robust Control/H\342\210\236 Control for Space Shuttle Flight Dynamics" @@ -0,0 +1,171 @@ +%% ========================================================================= +% space_shuttle.m - H∞ Control for Space Shuttle Flight Dynamics +% ========================================================================= +% This script models flight dynamics and control for a space shuttle. +% - Simulates actuators, sensor noise, and pilot input. +% - Designs an H∞ controller for flight stability. +% - Uses µ-analysis to evaluate robustness against disturbances. +% +% System Description: +% +% Aircraft Nominal Model: +% - 4 States, 7 Outputs, 6 Inputs +% +% Actuator Models: +% - Rudder and Elevon models with delay dynamics +% +% Weighting Functions: +% - Wind Gust +% - Sensor Noise +% - Pilot Command +% - Error and Actuator Weightings +% ========================================================================= + +%% 1. Define Aircraft Model (SSLAFCS) +clear all +clc + +% State-space matrices for the aircraft dynamics +A = [-0.0946 0.141 -0.99 0.0364; + -3.59 -0.428 0.281 0; + 0.395 -0.0126 -0.0814 0; + 0 1 -0.141 0]; + +B = [-0.0124 0.0102 -0.000000109; + 6.57 1.26 -0.00000413; + 0.378 -0.256 0.000000453; + 0 0 0]; + +C = [ 0 1 0 0; + 0 0 1 0; + -68 -1.74 -4.06 -3.72e-5; + 0 0 0 1]; + +% Additional matrices for uncertainty and performance evaluation +Bp = [ 0.0128 0 0; + 0 -0.0311 -3.12; + 0 -0.19 -0.0644; + 0 0 0]; + +Cp = [1 0 0 0; 0 0 0 0; 0 0 0 0]; +D12 = [0 0 0.00000115; 1 0 0; 0 1 0]; +D21 = [zeros(2,3); 11.1 -11.1 -11.1; 0 0 0]; +D22 = [zeros(2,3); 26.7 -2.95 -0.0000781; 0 0 0]; + +% Construct system using `pck` +acnom = pck(A, [Bp B], [Cp; C], [zeros(3,3) D12; D21 D22]); + +%%%% ========================================================================= +% 2. Define Actuator Models +%%%% ========================================================================= +% Define rudder and elevon actuator dynamics +wrud = 21.0; zetarud = 0.75; +wele = 14.0; zetaele = 0.72; +wdel = 173.0; zetadel = 0.866; + +% Delay approximation transfer function +delaytf = nd2sys([1 -2*zetadel*wdel wdel*wdel], [1 2*zetadel*wdel wdel*wdel]); + +% Rudder actuator system +actrud = mmult(nd2sys([wrud^2 0], [1 2*zetarud*wrud wrud^2]), delaytf); + +% Elevon actuator system +actele = mmult(nd2sys([wele^2 0], [1 2*zetaele*wele wele^2]), delaytf); + +%%%% ========================================================================= +% 3. Define Weighting Functions +%%%% ========================================================================= +% Wind disturbance +wgust = nd2sys([.5 1], [1 1], 30); + +% Sensor noise +wr = nd2sys([100 1], [2 1], .0003); +wp = wr; +wphi = nd2sys([100 1], [.5 1], 0.0007); +wny = nd2sys([20 1], [.1 1], .25); +wnoise = daug(wp, wr, wny, wphi); + +% Pilot command weighting +wphicmd = nd2sys([.5 1], [2 1], .5); + +% Actuator weighting +wact = diag([4 1 .005 2 .2 .009]); + +% Ideal response model +wmod = 1.2; zmod = .7; +idmod = nd2sys(wmod^2, [1 2*zmod*wmod wmod^2]); + +% Performance variables +nyerr = nd2sys([1 1], [10 1], .8); +cterr = nd2sys([1 1], [100 1], 500); +baerr = nd2sys([1 1], [100 1], 250); +fix = [0 0 1 0 0; 0 1 0 -0.037 0; 0 0 0 1 -1]; +wperf = mmult(daug(nyerr, cterr, baerr), fix); + +%%%% ========================================================================= +% 4. System Interconnection +%%%% ========================================================================= +systemnames = 'acnom wr wl wnoise wact wgust wphicmd wperf idmod actele actrud'; +inputvar = '[pertin{9}; noise{4}; gust; comd; elec; rudc]'; +outputvar = '[wr; wperf; wact; wphicmd; acnom(4:7) + wnoise]'; +sysoutname = 'olic'; +sysic + +disp('The Interconnection Structure is in "olic"') +minfo(olic) + +%%%% ========================================================================= +% 5. H∞ Controller Design +%%%% ========================================================================= +olic_h = sel(olic, [10:23], [10:17]); +k_h = hinfsyn(olic_h, 5, 2, 0, 1, 0.01, 2); + +%%%% ========================================================================= +% 6. Load and Compare Controllers +%%%% ========================================================================= +load shutcont % Load µ-controller and another trade-off controller +minfo(k_x) +minfo(k_mu) + +%%%% ========================================================================= +% 7. Compute Frequency Response of Controllers +%%%% ========================================================================= +omega = logspace(-2, 3, 50); +clp_h = starp(olic, k_h, 5, 2); +clp_x = starp(olic, k_x, 5, 2); +clp_mu = starp(olic, k_mu, 5, 2); + +clp_hg = frsp(clp_h, omega); +clp_xg = frsp(clp_x, omega); +clp_mug = frsp(clp_mu, omega); + +% Nominal performance plot +vplot('liv,m', vnorm(sel(clp_hg, [10:18], [10:15])), '-.', ... + vnorm(sel(clp_xg, [10:18], [10:15])), '--', ... + vnorm(sel(clp_mug, [10:18], [10:15])), '-'); +grid; title('Nominal Performance: All Controllers') +legend('H∞', 'H_x', 'µ') + +pause + +%%%% ========================================================================= +% 8. Robust Stability Analysis +%%%% ========================================================================= +delsetrs = ones(9,2); +clp_hgrs = sel(clp_hg, 1:9, 1:9); +bnds_h = mu(clp_hgrs, delsetrs); + +% Robust performance analysis +delsetrp = [delsetrs; 6 9]; +bnds_h = mu(clp_hg, delsetrp); + +vplot('liv,d', bnds_h, '-.r'); +grid; title('Robust Stability of Closed-Loop System') +legend('H∞') + +pause + +%%%% ========================================================================= +% 9. Simulation Preparation +%%%% ========================================================================= +disp('Run mk_olsim to continue with GUI-based simulations.') From 7c3b4efecf81c370838daf1a424d6309cae349e2 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:22:06 +0330 Subject: [PATCH 20/22] =?UTF-8?q?Create=20Robust=20H=E2=88=9E=20Control=20?= =?UTF-8?q?for=20the=20HIMAT=20Aircraft?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...42\210\236 Control for the HIMAT Aircraft" | 108 ++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 "Aerospace Robust Control/Robust H\342\210\236 Control for the HIMAT Aircraft" diff --git "a/Aerospace Robust Control/Robust H\342\210\236 Control for the HIMAT Aircraft" "b/Aerospace Robust Control/Robust H\342\210\236 Control for the HIMAT Aircraft" new file mode 100644 index 0000000..132c4ee --- /dev/null +++ "b/Aerospace Robust Control/Robust H\342\210\236 Control for the HIMAT Aircraft" @@ -0,0 +1,108 @@ +%% ========================================================================= +% ex_himat.m - Robust H∞ Control for the HIMAT Aircraft +% ========================================================================= +% This script designs a robust H∞ controller for the HIMAT aircraft. +% - Defines the aircraft's state-space model. +% - Uses `hinfsyn` for controller synthesis. +% - Performs µ-analysis and D-K iteration to optimize robustness. +% +% System Description: +% - HIMAT Aircraft Model (Highly Maneuverable Aircraft Technology) +% - Actuator and disturbance weighting functions +% - Robust H∞ synthesis and performance evaluation +% ========================================================================= + +%% 1. Define Aircraft Model (HIMAT) +clear all +clc + +% HIMAT state-space model +himat = [ ... + -0.0226 -36.6000 -18.9000 -32.1000 0 0 4.0000; + 0 -1.9000 0.9830 0 -0.4140 0 0; + 0.0123 -11.7000 -2.6300 0 -77.8000 22.4000 0; + 0 0 1.0000 0 0 0 0; + 0 57.3000 0 0 0 0 0; + 0 0 0 57.3000 0 0 0; + 0 0 0 0 0 0 -Inf]; + +%%%% ========================================================================= +% 2. Define Weighting Functions +%%%% ========================================================================= +% Uncertainty weight +WDEL = nd2sys([50,5000], [1,10000]); +WDEL = daug(WDEL, WDEL); + +% Performance weighting +WP = nd2sys([0.5,1.5], [1,0.03]); +WP = daug(WP, WP); + +%%%% ========================================================================= +% 3. Define System Interconnections +%%%% ========================================================================= +systemnames = 'himat WDEL WP'; +inputvar = '[pert(2); dist(2); control(2)]'; +outputvar = '[ WDEL ; WP; himat + dist ]'; +input_to_himat = '[ control + pert ]'; +input_to_WDEL = '[ control ]'; +input_to_WP = '[ himat + dist ]'; +sysoutname = 'himat_ic'; +cleanupsysic = 'yes'; +sysic + +%%%% ========================================================================= +% 4. H∞ Controller Synthesis +%%%% ========================================================================= +minfo(himat_ic); +[k1, clp1] = hinfsyn(himat_ic, 2, 2, 0.8, 2, 1e-1, 2); + +%%%% ========================================================================= +% 5. Frequency Response Analysis +%%%% ========================================================================= +minfo(k1) +omega = logspace(-3, 3, 60); +clp1g = frsp(clp1, omega); +clp1gs = vsvd(clp1g); % Singular value calculation + +figure; +vplot('liv,m', clp1gs), grid +title('Singular Value Plot for the Closed-Loop System') +xlabel('Frequency (rad/sec)') + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 6. Robust Performance Analysis (µ-Analysis) +%%%% ========================================================================= +deltaset = [2 2; 2 2]; +clp_g1 = clp1g; +[bnds1, pv1] = mu(clp_g1, deltaset); + +figure; +vplot('liv,m', vnorm(clp_g1), bnds1); +grid +title('Maximum Singular Value and µ Plot') +xlabel('Frequency (rad/sec)') + +pause % Press any key to continue +clc + +%%%% ========================================================================= +% 7. Perform D-K Iteration for Robust Optimization +%%%% ========================================================================= +% Set DK parameter +DK_DEF_NAME = 'himat_def'; % 3 iteration is set +dkit + +clc + +%%%% ========================================================================= +% 8. Reduce Controller Order using Balanced Truncation +%%%% ========================================================================= +[k_dk3bal, hsv] = sysbal(k_dk3); +[k_dk3red] = strunc(k_dk3bal, 13); +clpred_13 = starp(himat_ic, k_dk3red); +clpred_13g = frsp(clpred_13, OMEGA_DK); +[bnds] = mu(clpred_13g, [2 2; 2 2], 'c'); +pkvnorm(sel(bnds, 1, 1)); From 7edf5c56994ef7128c84409019b164623589b5ee Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:34:58 +0330 Subject: [PATCH 21/22] Create Default Configuration for D-K Iteration --- .../Default Configuration for D-K Iteration | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 Aerospace Robust Control/Default Configuration for D-K Iteration diff --git a/Aerospace Robust Control/Default Configuration for D-K Iteration b/Aerospace Robust Control/Default Configuration for D-K Iteration new file mode 100644 index 0000000..4b63d11 --- /dev/null +++ b/Aerospace Robust Control/Default Configuration for D-K Iteration @@ -0,0 +1,37 @@ +%% ========================================================================= +% himat_def.m - Default Configuration for D-K Iteration +% ========================================================================= +% This script provides user-defined parameters for D-K iteration in the +% HIMAT aircraft's robust control design. +% - Defines uncertainty block structure for µ-analysis. +% - Specifies the number of measurements and control inputs. +% - Sets frequency range and iteration limits for optimization. +% ========================================================================= + +%% 1. Define Nominal Plant Interconnection Structure +NOMINAL_DK = himat_ic; % Nominal interconnection for D-K iteration + +%% 2. Define Number of Measurements and Control Inputs +NMEAS_DK = 1; % Number of measured outputs +NCONT_DK = 1; % Number of control inputs + +%% 3. Define Block Structure for µ-Analysis +BLK_DK = [1 1; 1 1]; % Uncertainty block structure + +%% 4. Define Frequency Range for D-K Iteration +OMEGA_DK = logspace(-3, 3, 60); % Logarithmic frequency range + +%% 5. Define Automatic Iteration Parameters +% AUTOINFO_DK: [Start Iteration, End Iteration, Visual Flag, Scaling Orders] +% - Start Iteration: 1 +% - End Iteration: 5 +% - Visual Flag: 1 (Show results) +% - Scaling Order: 5 (Maximum order of dynamic D-scalings) +AUTOINFO_DK = [1 5 1 5 5]; + +%% 6. Define Naming Convention for Exported Variables +NAME_DK = ''; % No suffix added to exported variables + +%% ========================================================================= +% End of himat_def.m +% ========================================================================= From 499c0b6f2b62683808fb54f3e57d91a06625d173 Mon Sep 17 00:00:00 2001 From: Kosar aminjafari <91948388+kosarAmin@users.noreply.github.com> Date: Tue, 18 Mar 2025 05:35:34 +0330 Subject: [PATCH 22/22] Create Alternative D-K Iteration Configuration --- .../Alternative D-K Iteration Configuration | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 Aerospace Robust Control/Alternative D-K Iteration Configuration diff --git a/Aerospace Robust Control/Alternative D-K Iteration Configuration b/Aerospace Robust Control/Alternative D-K Iteration Configuration new file mode 100644 index 0000000..6ac96e1 --- /dev/null +++ b/Aerospace Robust Control/Alternative D-K Iteration Configuration @@ -0,0 +1,36 @@ +%% ========================================================================= +% himat_def1.m - Alternative D-K Iteration Configuration +% ========================================================================= +% This script is a variant of `himat_def.m`, modifying D-K iteration settings +% for the HIMAT aircraft's robust control tuning. +% - Adjusts iteration starting values. +% - Uses alternative frequency response constraints. +% ========================================================================= + +%% 1. Define Nominal Plant Interconnection Structure +NOMINAL_DK = himat_ic; % Nominal interconnection for D-K iteration + +%% 2. Define Number of Measurements and Control Inputs +NMEAS_DK = 1; % Number of measured outputs +NCONT_DK = 1; % Number of control inputs + +%% 3. Define Block Structure for µ-Analysis +BLK_DK = [1 1; 1 1]; % Uncertainty block structure + +%% 4. Define Frequency Range for D-K Iteration +OMEGA_DK = logspace(-3, 3, 60); % Logarithmic frequency range + +%% 5. Define Automatic Iteration Parameters +% AUTOINFO_DK: [Start Iteration, End Iteration, Visual Flag, Scaling Orders] +% - Start Iteration: 2 +% - End Iteration: 5 +% - Visual Flag: 1 (Show results) +% - Scaling Order: 5 (Maximum order of dynamic D-scalings) +AUTOINFO_DK = [2 5 1 5 5]; + +%% 6. Define Naming Convention for Exported Variables +NAME_DK = ''; % No suffix added to exported variables + +%% ========================================================================= +% End of himat_def1.m +% =========================================================================