Modif sur le .m avec W ppur wc

BE_Beam.m

@@ -1,16 +1,22 @@
 %% [BE] Digital control of cantilever beam
+
 % Justin Bos
-% 
 % Wissal Guarni
-% 
 % Nolan Reynier Nomer
-% 
 % Aleksander Taban
+clear;
+clc;
+close all;
+
+set(0, 'DefaultLineLineWidth', 1)
+
+%% Question 3
 
 E = [13/35 9/70 11/210 -13/420 ;
      9/70 13/35 13/420 -11/210 ;
      11/210 13/420 1/105 -1/140 ;
      -13/420 -11/210 -1/140 1/105] ;
+
 D = [-6/5 6/5 -1/10 -1/10 ;
      6/5 -6/5 1/10 11/10 ;
      -11/10 1/10 -2/15 1/30 ;
@@ -20,58 +26,91 @@ phi = @(zeta) [2*zeta^3-3*zeta^2+1 ;
             3*zeta^2-2*zeta^3 ; 
             zeta^3-2*zeta^2+zeta ; 
             zeta^3-zeta^2 ] ; 
-L=1;
+
-%%
+L = 1 ; % Longueur de la poutre
+phi_L = phi(L) ;
+
 A = [zeros(4), inv(E)*D ;
     -inv(E)*(D'), zeros(4) ] ;
-phi_L = phi(L) ;
 B = [zeros(4,1) ; -inv(E)*phi_L] ;
 C = [zeros(1,4) , -phi_L'] ;
-%%
+
+A, B, C
+
+%% Question 4
 Ts1 = 0.01 ;
 Ts2 = 0.02 ;
 Ts3 = 0.04 ;
 
-eig_cont = eig(A);
 sys = ss(A, B, C, 0) ;
-figure;
-subplot(1,3,1) ;
-plot(eig_cont, 'x');
-hold off;
 
-subplot(1,3,2);
+eig_cont = eig(sys);
 eig_tustin1 = eig(c2d(sys, Ts1, 'tustin')) ;
 eig_tustin2 = eig(c2d(sys, Ts2, 'tustin')) ;
 eig_tustin3 = eig(c2d(sys, Ts3, 'tustin')) ;
-plot(eig_tustin1,'x');
+eig_zoh1 = eig(c2d(sys, Ts1, 'zoh')) ;
+eig_zoh2 = eig(c2d(sys, Ts2, 'zoh')) ;
+eig_zoh3 = eig(c2d(sys, Ts3, 'zoh')) ;
+
+figure;
+subplot(1,3,1) ;
+plot(eig_cont, 'x');
+grid on;
+xline(0, 'r--');
+xlabel('Axe Réel');
+ylabel('Axe Imaginaire');
+title('Pôles du système');
+
+subplot(1,3,2);
 hold on;
-plot(eig_tustin2,'x');
+plot(eig_tustin1,'x', 'LineWidth',2);
-plot(eig_tustin3,'x');
+plot(eig_tustin2,'x', 'LineWidth',2);
+plot(eig_tustin3,'x', 'LineWidth',2);
+
+theta = linspace(0, 2*pi, 100);
+plot(cos(theta), sin(theta), 'k--', 'HandleVisibility','off');
+xlim([-1.1, 1.1]);
+ylim([-1.1, 1.1]);
+grid on;
+title('Pôles dans le plan en Z (Discret)');
+xlabel('Réel');
+ylabel('Imaginaire');
 hold off;
 
 subplot(1,3,3);
-eig_zoh1 = eig(c2d(sys, Ts1, 'zoh')) ;
-sys_zoh1 = c2d(sys, Ts1, 'zoh');
-eig_zoh2 = eig(c2d(sys, Ts2, 'zoh')) ;
-eig_zoh3 = eig(c2d(sys, Ts3, 'zoh')) ;
-plot(eig_zoh1,'x');
 hold on;
-pzmap(sys_zoh1);
+plot(eig_zoh1,'x', 'LineWidth',1);
-plot(eig_zoh3,'x');
+plot(eig_zoh2,'x', 'LineWidth',1);
+plot(eig_zoh3,'x', 'LineWidth',1);
+
+legend("Eigen Values Ts1", "Eigen Values Ts2", "Eigen Values Ts3");
+theta = linspace(0, 2*pi, 100);
+plot(cos(theta), sin(theta), 'k--', 'HandleVisibility','off');
+xlim([-1.1, 1.1]);
+ylim([-1.1, 1.1]);
+grid on;
+title('Valeurs propres avec un bloqueur d''ordre 0 (ZOH)');
+xlabel('Réel');
+ylabel('Imaginaire');
 hold off;
-%% 
+
 % Conlsuion : On en conclut qqc.
 
-figure;
+
+%% Question 5
+
 t = 0:0.01:10;
-u = double (t >= 1) ;
+cmd = double (t >= 1) ;
-lsim(sys,u, t);
+[y, t] = lsim(sys,cmd,t) ;
-title();
+
-legend();
+figure ;
-xlabel();
+plot(t, cmd, t, y);
-ylabel();
+xlabel("Temps (s)");
-%% 
+ylabel("Amplitude");
-% 
+legend("u(t)", "y(t)");
+title("Simulation en boucle ouverte") ;
+
+%% Question 7
 
 Cw = @(zeta) [zeta^2*(2*zeta^3-5*zeta^2+10)/20 ...
             -zeta^4*(2*zeta-5)/20 ...
@@ -85,12 +124,24 @@ for i = 0.01:0.01:10
     TransientTime(end+1) = stepinfo(feedback(sys, i)).TransientTime ;
 end
 [minimum, min_index] = min(TransientTime)
-k = min_index * 0.01;
+k = min_index * 0.01 
-lsim(feedback(sys,k), u, t) ;
+
-H = -3
+%% Question 8 
-%%
+
+H = inv(Cw_L*inv(-A+B*C*k)*B)
+[y, t, x] = lsim(feedback(sys,k), H*cmd, t) ;
+u = H*cmd'-k*y ;
+figure ;
+plot(t, u, t, y, t, cmd, t, (Cw_L*x')');
+xlabel("Temps (s)");
+ylabel("Amplitude");
+legend("u(t)", "y(t)", "wc(t)", "w(t)");
+title("Simulation du retour de sortie avec pré-gain") ;
+
+%% Question 9
+
 sys_CL = feedback(sys,k);
-[y, tOut] = lsim(sys_CL, -3.*u, t);
+[y, tOut] = lsim(sys_CL, -3.*cmd, t);
 % plot(t, Cw(L)*y, t, u)
 eig_cont_CL = eig(sys_CL);
 figure;
@@ -110,7 +161,21 @@ hold off;
 sys_CL_Ts = c2d(sys_CL, 10, 'zoh');
 figure
 step(sys_CL_Ts)
-%%
+
-K = lqr(sys, eye(size(A)), 1)
+%% Question 11
-sys_CL_lqr = ss(A-B*K, B*H, C, 0);
+
-lsim(sys_CL_lqr, u, t);
+K = lqr(sys, eye(size(A)), 1);
+H_lqr = inv(Cw_L*inv(-A+B*K)*B) ;
+
+%% Question 12
+
+sys_CL_lqr = ss(A-B*K, B*H_lqr, C, 0);
+[y, t, x] = lsim(sys_CL_lqr, cmd, t);
+u = H_lqr*cmd'-k*y ;
+figure ;
+plot(t, u, t, y, t, cmd, t, (Cw_L*x')');
+xlabel("Temps (s)");
+ylabel("Amplitude");
+legend("u(t)", "y(t)", "wc(t)", "w(t)");
+title("Simulation du retour d'états avec pré-gain") ;
+