diff --git a/BE_Beam.m b/BE_Beam.m index 454c889..610fdba 100644 --- a/BE_Beam.m +++ b/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 ; @@ -19,59 +25,92 @@ D = [-6/5 6/5 -1/10 -1/10 ; 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; -%% + zeta^3-zeta^2 ] ; + +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_tustin3,'x'); +plot(eig_tustin1,'x', 'LineWidth',2); +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_zoh3,'x'); +plot(eig_zoh1,'x', 'LineWidth',1); +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) ; -lsim(sys,u, t); -title(); -legend(); -xlabel(); -ylabel(); -%% -% +cmd = double (t >= 1) ; +[y, t] = lsim(sys,cmd,t) ; + +figure ; +plot(t, cmd, t, y); +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; -lsim(feedback(sys,k), u, t) ; -H = -3 -%% +k = min_index * 0.01 + +%% 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) -sys_CL_lqr = ss(A-B*K, B*H, C, 0); -lsim(sys_CL_lqr, u, t); \ No newline at end of file + +%% Question 11 + +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") ; + diff --git a/BE_Beam.mlx b/BE_Beam.mlx index 8e6beb3..67cd24a 100644 Binary files a/BE_Beam.mlx and b/BE_Beam.mlx differ