257 lines
6.3 KiB
Matlab
257 lines
6.3 KiB
Matlab
%% [BE] Digital control of cantilever beam
|
|
|
|
% Justin Bos
|
|
% Wissal Guarni
|
|
% Nolan Reynier Nomer
|
|
% Aleksander Taban
|
|
clear;
|
|
clc;
|
|
close all;
|
|
|
|
SAVE = 0;
|
|
|
|
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 ;
|
|
-1/10 1/10 1/30 -2/15] ;
|
|
|
|
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 ; % Longueur de la poutre
|
|
phi_L = phi(L) ;
|
|
|
|
A = [zeros(4), inv(E)*D ;
|
|
-inv(E)*(D'), zeros(4) ] ;
|
|
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 ;
|
|
|
|
sys = ss(A, B, C, 0) ;
|
|
|
|
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')) ;
|
|
eig_zoh1 = eig(c2d(sys, Ts1, 'zoh')) ;
|
|
eig_zoh2 = eig(c2d(sys, Ts2, 'zoh')) ;
|
|
eig_zoh3 = eig(c2d(sys, Ts3, 'zoh')) ;
|
|
%%
|
|
Question4 = figure ;
|
|
subplot(1,3,1) ;
|
|
plot(eig_cont, 'x');
|
|
axis square;
|
|
grid on;
|
|
xline(0, 'r--');
|
|
xlabel('Axe Réel');
|
|
ylabel('Axe Imaginaire');
|
|
title({'Valeurs propres en temps continu'});
|
|
|
|
subplot(1,3,2);
|
|
hold on;
|
|
plot(eig_tustin1,'x', 'LineWidth',1);
|
|
plot(eig_tustin2,'x', 'LineWidth',1);
|
|
plot(eig_tustin3,'x', 'LineWidth',1);
|
|
axis square;
|
|
theta = linspace(0, 2*pi, 100);
|
|
plot(cos(theta), sin(theta), 'k--', 'HandleVisibility','off','LineWidth',0.25);
|
|
xlim([-1.1, 1.1]);
|
|
ylim([-1.1, 1.1]);
|
|
grid on;
|
|
title({'Valeurs propres en temps discret',' avec la méthode de Tustin'});
|
|
legend("Ts1", "Ts2", "Ts3");
|
|
xlabel('Réel');
|
|
ylabel('Imaginaire');
|
|
hold off;
|
|
|
|
subplot(1,3,3);
|
|
hold on;
|
|
plot(eig_zoh1,'x', 'LineWidth',1);
|
|
plot(eig_zoh2,'x', 'LineWidth',1);
|
|
plot(eig_zoh3,'x', 'LineWidth',1);
|
|
axis square;
|
|
legend("Ts1", "Ts2", "Ts3");
|
|
theta = linspace(0, 2*pi, 100);
|
|
plot(cos(theta), sin(theta), 'k--', 'HandleVisibility','off','LineWidth',0.25);
|
|
xlim([-1.1, 1.1]);
|
|
ylim([-1.1, 1.1]);
|
|
grid on;
|
|
title({'Valeurs propres en temps discret',' avec un bloqueur d''ordre 0'});
|
|
xlabel('Réel');
|
|
ylabel('Imaginaire');
|
|
hold off;
|
|
|
|
if SAVE
|
|
exportgraphics(Question4, './latex/Illustrations/Question4.pdf')
|
|
end
|
|
%% Question 5
|
|
|
|
t = 0:0.01:10;
|
|
cmd = double (t >= 1) ;
|
|
[y, t] = lsim(sys,cmd,t) ;
|
|
|
|
Question5 = figure ;
|
|
plot(t, cmd, t, y);
|
|
xlabel("Temps (s)");
|
|
ylabel("Amplitude");
|
|
legend("$u(t)$", "$y(t)$", 'Interpreter', 'latex');
|
|
title("Simulation en boucle ouverte") ;
|
|
|
|
if SAVE
|
|
exportgraphics(Question5, './latex/Illustrations/Question5.pdf')
|
|
end
|
|
|
|
%% Question 7
|
|
|
|
Cw = @(zeta) [zeta^2*(2*zeta^3-5*zeta^2+10)/20 ...
|
|
-zeta^4*(2*zeta-5)/20 ...
|
|
zeta^3*(3*zeta^2-10*zeta+10)/60 ...
|
|
zeta^4*(3*zeta-5)/60 ...
|
|
0, 0, 0, 0];
|
|
Cw_L = Cw(L);
|
|
%%
|
|
|
|
TransientTime = zeros(size(0.01:0.01:10)) ;
|
|
for i = 1:size(TransientTime, 2)
|
|
TransientTime(i) = stepinfo(feedback(sys, i*0.01)).TransientTime ;
|
|
end
|
|
[minimum, min_index] = min(TransientTime)
|
|
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 ;
|
|
%%
|
|
Question8 = figure ;
|
|
plot(t, u, t, y, t, cmd, t, (Cw_L*x')');
|
|
xlabel("Temps (s)");
|
|
ylabel("Amplitude");
|
|
legend("$u(t)$", "$y(t)$", "$w_c(L,t)$", "$w(L,t)$", 'Interpreter', 'latex');
|
|
title("Simulation du retour de sortie avec pré-gain") ;
|
|
|
|
if SAVE
|
|
exportgraphics(Question8, './latex/Illustrations/Question8.pdf')
|
|
end
|
|
|
|
%% Question 9
|
|
|
|
sys_CL = feedback(sys,k);
|
|
[y, tOut] = lsim(sys_CL, -3.*cmd, t);
|
|
% plot(t, Cw(L)*y, t, u)
|
|
eig_cont_CL = eig(sys_CL);
|
|
figure;
|
|
subplot(1,2,1) ;
|
|
plot(eig_cont_CL, 'x');
|
|
hold off;
|
|
|
|
subplot(1,2,2);
|
|
eig_zoh1 = eig(c2d(sys_CL, Ts1, 'zoh')) ;
|
|
eig_zoh2 = eig(c2d(sys_CL, Ts2, 'zoh')) ;
|
|
eig_zoh3 = eig(c2d(sys_CL, Ts3, 'zoh')) ;
|
|
plot(eig_zoh1,'x');
|
|
hold on;
|
|
plot(eig_zoh3,'x');
|
|
hold off;
|
|
|
|
%%
|
|
sys_CL_Ts = c2d(sys_CL, 1, 'zoh');
|
|
figure
|
|
step(sys_CL_Ts)
|
|
|
|
%% Question 11
|
|
|
|
K_lqr = lqr(sys, eye(size(A)), 1);
|
|
H_lqr = inv(Cw_L*inv(-A+B*K_lqr)*B) ;
|
|
|
|
%% Question 12
|
|
|
|
sys_CL_lqr = ss(A-B*K_lqr, B*H_lqr, C, 0);
|
|
[y, t, x] = lsim(sys_CL_lqr, cmd, t);
|
|
u = H_lqr*cmd'-(K_lqr*x')' ;
|
|
Question12 = figure ;
|
|
plot(t, u, t, y, t, cmd, t, (Cw_L*x')');
|
|
xlabel("Temps (s)");
|
|
ylabel("Amplitude");
|
|
legend("$u(t)$", "$y(t)$", "$w_c(L,t)$", "$w(L,t)$", 'Interpreter', 'latex');
|
|
title("Simulation du retour d'états avec pré-gain") ;
|
|
|
|
if SAVE
|
|
exportgraphics(Question12, './latex/Illustrations/Question12.pdf')
|
|
end
|
|
|
|
%% Question 15
|
|
Bp = [0; 0; 0; 0; 1/2*L^4-L^3+L; L^3-1/2*L^4; 1/4*L^4-2/3*L^3+1/2*L^2; 1/4*L^4-1/3*L^3];
|
|
|
|
q0 = -5 * double (t' >= 4);
|
|
sys_perturbations = ss(A-B*K_lqr, [B*H_lqr Bp], C, [0 0]);
|
|
|
|
%%Question 16
|
|
[y, t, x] = lsim(sys_perturbations, [cmd ; q0], t');
|
|
u = H_lqr*cmd'-(K_lqr*x')' ;
|
|
Question16 = figure ;
|
|
plot(t, u, t, y, t, cmd, t, (Cw_L*x')');
|
|
xlabel("Temps (s)");
|
|
ylabel("Amplitude");
|
|
legend("$u(t)$", "$y(t)$", "$w_c(L,t)$", "$w(L,t)$", 'Interpreter', 'latex');
|
|
title("Simulation du retour d'états avec pré-gain et perturbations") ;
|
|
|
|
if SAVE
|
|
exportgraphics(Question16, './latex/Illustrations/Question16.pdf')
|
|
end
|
|
|
|
%% Question 17
|
|
eig_lqr = eig(sys_CL_lqr);
|
|
eig_4 = [eig_lqr ; -2];
|
|
A_aug = [A, zeros(8,1) ; -Cw_L, 0];
|
|
B_aug = [B; 0];
|
|
B_p_aug = [Bp; 0];
|
|
B_ref = [zeros(8,1); 1];
|
|
|
|
%%
|
|
K_aug = place(A_aug, B_aug, eig_4)
|
|
sys_integral = ss(A_aug - B_aug*K_aug, [B_ref B_p_aug], [C, 0], [0 0]);
|
|
|
|
[y, t, x] = lsim(sys_integral, [cmd ; q], t');
|
|
u = -(K_aug*x')' ;
|
|
Question19 = figure ;
|
|
plot(t, u, t, y, t, cmd, t, ([Cw_L 0]*x')');
|
|
xlabel("Temps (s)");
|
|
ylabel("Amplitude");
|
|
legend("$u(t)$", "$y(t)$", "$w_c(L,t)$", "$w(L,t)$", 'Interpreter', 'latex');
|
|
title("Simulation du retour d'états avec action intégrale et perturbations") ;
|
|
if SAVE
|
|
exportgraphics(Question19, './latex/Illustrations/Question19.pdf')
|
|
end
|
|
|
|
%%
|
|
sys_integral_d = c2d(sys_integral, 0.5, 'zoh');
|
|
t_d = 0:0.5:10; % Define the time vector for discrete simulation
|
|
cmd_d = double (t_d >= 1) ;
|
|
q0_d = -5 * double(t_d >= 4);
|
|
[y_d, t_d, x_d] = lsim(sys_integral_d, [cmd_d ; q0_d], t_d);
|
|
u_d = -(K_aug*x_d')' ;
|
|
Question20 = figure ;
|
|
plot(t_d, u_d, t_d, y_d, t_d, cmd_d, t_d, ([Cw_L 0]*x_d')');
|
|
xlabel("Temps (s)");
|
|
ylabel("Amplitude");
|
|
legend("$u(t)$", "$y(t)$", "$w_c(L,t)$", "$w(L,t)$", 'Interpreter', 'latex');
|
|
title("Simulation du retour d'états avec action intégrale et perturbations") ;
|