Untitled

avatar
unknown
plain_text
a year ago
15 kB
7
Indexable
clear all;
M=5;
m=1;
k=5400;
c1=10.8;
c2=6*10^-3;
L=0.1;
g=10;

%% question 3

y = -2*(M+m)*g/k:0.0001:0;
theta = -1.5*pi:0.01:1.5*pi;
[X, Y]=meshgrid(y,theta);
V = 0.5*k*X.^2+g.*X.*(M+m)-m*g*L*cos(Y);
NL = 30; % Number of level lines
HL = min(min(V));
BL = max(max(V));

contour(X,Y,V,NL,'LineWidth',1); hold all
contour(X,Y,V,[HL-0.001 HL+0.001],'k','LineWidth',5,'ShowText','on');
contour(X,Y,V,[BL-0.01 BL],'k','LineWidth',5,'ShowText','on');
h = surf(X,Y,V);

xlabel('y','interpreter','Latex')
ylabel('$\theta$','interpreter','Latex')
zlabel('V','interpreter','Latex')
set(h,'LineStyle','none')
colorbar
exportgraphics(gcf,'Figure_1.png','Resolution',300);
clear all;

%% question 4

M=5;
m=1;
k=5400;
c1=10.8;
c2=6*10^-3;
L=0.1;
g=10;
syms z1 z2 z3 z4 F TT
 
 
dzdtt(1,1) = z3;
dzdtt(2,1) = z4;
dzdtt(3,1) = (-1).*g+L.^(-1).*(m+2.*M+m.*cos(2.*z2)).^(-1).*((-2).*L.*((-1).*F+...
k.*z1+c1.*z3+L.*m.*z4.^2.*cos(z2))+2.*((-1).*TT+c2.*z4).*sin(z2));
 
dzdtt(4,1) = L.^(-2).*m.^(-1).*(m+2.*M+m.*cos(2.*z2)).^(-1).*((-2).*(m+M).*((...
-1).*TT+c2.*z4)+2.*L.*m.*((-1).*F+k.*z1+c1.*z3+L.*m.*z4.^2.*cos(...
z2)).*sin(z2));
 
 
y_of_t= [z1 z2];
 
A=jacobian(dzdtt,[z1 z2 z3 z4]);
B=jacobian(dzdtt,[F TT]);
C=jacobian(y_of_t,[z1 z2 z3 z4]);
D=jacobian(y_of_t,[F TT]);
 
A4=subs(A,[z1,z2,z3,z4,F,TT],[-1/90 0 0 0 0 0]) ;
B4=subs(B,[z1,z2,z3,z4,F,TT],[-1/90 0 0 0 0 0]);
C4=subs(C,[z1,z2,z3,z4,F,TT],[-1/90 0 0 0 0 0]);
D4=subs(D,[z1,z2,z3,z4,F,TT],[-1/90 0 0 0 0 0]);
 
A4=double(A4);
B4=double(B4);
C4=double(C4);
D4=double(D4);
 
 
%% question 6
clear dzdtt
syms z1 z2 TT
 
dzdtt(1,1) = z2;
dzdtt(2,1) = (TT-c2*z2-m*g*L*sin(z1))/(m*L^2);
 
y_of_t= z1;
 
A2=jacobian(dzdtt,[z1 z2]);
B2=jacobian(dzdtt,TT);
C2=jacobian(y_of_t,[z1 z2]);
D2=jacobian(y_of_t,TT);
 
A42=subs(A2,[z1,z2,TT],[0 0 0]) ;
B42=subs(B2,[z1,z2,TT],[0 0 0]);
C42=subs(C2,[z1,z2,TT],[0 0 0]);
D42=subs(D2,[z1,z2,TT],[0 0 0]);
 
A42=double(A42);
B42=double(B42);
C42=double(C42);
D42=double(D42);

%% Question 7.1

IC1 = [5*pi/180;0];
IC2 = [50*pi/180;0];
 
Linear_system = ss(A42,B42,C42,D42);
tt = linspace(0,30,1e3).';
Y1 = initial(Linear_system,IC1,tt);

figure(2);
plot(tt,Y1.*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('$\theta$(degre)','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_2.png','Resolution',300);
 
Y2 = initial(Linear_system,IC2,tt);
figure(3);
plot(tt,Y2.*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('$\theta$(degre)','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_3.png','Resolution',300);
 
%% Question 7.2

opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
tt = linspace(0,30,1e4);
 
[~,z_1] = ode45(@(t,z) EX_4_NLode_SDOF(t,z,m,c2,L,g), tt, IC1, opts);
 
[~,z_2] = ode45(@(t,z) EX_4_NLode_SDOF(t,z,m,c2,L,g), tt, IC2, opts);
 
figure(4);
plot(tt,z_1(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('$\theta$(degre)','interpreter','Latex')
title('Nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_4.png','Resolution',300);
 
figure(5);
plot(tt,z_2(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('$\theta$(degre)','interpreter','Latex')
title('Nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_5.png','Resolution',300);
 
%% question 8

ICL1 = [0;5*pi/180;0;0];
ICL2 = [-1/900;0;0;0];
ICL3 = [-1/900;5*pi/180;0;0];
Linear_system = ss(A4,B4,C4,D4);

tt = linspace(0,30,1e3).';
Y1 = initial(Linear_system,ICL1,tt);
Y2 = initial(Linear_system,ICL2,tt);
Y3 = initial(Linear_system,ICL3,tt);
 
figure(6);
yyaxis left
plot(tt,Y1(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,Y1(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_6.png','Resolution',300);
 
figure(7);
yyaxis left
plot(tt,Y2(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,Y2(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_7.png','Resolution',300);
 
figure(8);
yyaxis left
plot(tt,Y3(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,Y3(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_8.png','Resolution',300);
 
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
tt = linspace(0,30,1e4);
ICNL1 = [-1/90;5*pi/180;0;0];
ICNL2 = [-11/900;0;0;0];
ICNL3 = [-11/900;5*pi/180;0;0];
 
[~,z_1] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,0,0,0), tt, ICNL1, opts);
 
[~,z_2] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,0,0,0), tt, ICNL2, opts);
 
[~,z_3] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,0,0,0), tt, ICNL3, opts);
 
 
figure(9);
yyaxis left
plot(tt,z_1(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,z_1(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_9.png','Resolution',300);
 
figure(10);
yyaxis left
plot(tt,z_2(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,z_2(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_10.png','Resolution',300);
 
figure(11);
yyaxis left
plot(tt,z_3(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,z_3(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_11.png','Resolution',300);
 
%% Question 9

ICL = [0;5*pi/180;0;0];
ICNL = [-1/90;5*pi/180;0;0];
omega = [10,20,30];
 
tt6 = linspace(0,30,1e4);
 
U1 = [20*sin(omega(1)*tt6)', zeros(length(tt6),1)].';
U2 = [20*sin(omega(2)*tt6)', zeros(length(tt6),1)].';
U3 = [20*sin(omega(3)*tt6)', zeros(length(tt6),1)].';
 
Yf1= lsim(Linear_system,U1,tt6,ICL);
Yf2= lsim(Linear_system,U2,tt6,ICL);
Yf3= lsim(Linear_system,U3,tt6,ICL);
 
figure(12);
yyaxis left
plot(tt,Yf1(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,Yf1(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_12.png','Resolution',300);
 
figure(13);
yyaxis left
plot(tt,Yf2(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,Yf2(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_13.png','Resolution',300);
 
figure(14);
yyaxis left
plot(tt,Yf3(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,Yf3(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_14.png','Resolution',300);
 
 
[~,z_1] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,omega(1),20,0,0), tt, ICNL, opts);
 
[~,z_2] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,omega(2),20,0,0), tt, ICNL, opts);
 
[~,z_3] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,omega(3),20,0,0), tt, ICNL, opts);
 
figure(15);
yyaxis left
plot(tt,z_1(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,z_1(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_15.png','Resolution',300);
 
figure(16);
yyaxis left
plot(tt,z_2(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,z_2(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_16.png','Resolution',300);
 
figure(17);
yyaxis left
plot(tt,z_3(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt,z_3(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_17.png','Resolution',300);
 
 
%% Question 10.a
 
AF = [10 200];
AT = [0.2 0.8];
ICNL = [-1/90;0;0;0];
ICL = [0;0;0;0];
tt8=linspace(0,20,2e4);
[~,z_1] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,AF(1),0,pi/2), tt8, ICNL, opts);
[~,z_2] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,AF(2),0,pi/2), tt8, ICNL, opts);
Ys=step(Linear_system,tt8);
Yssqz= squeeze(Ys(:,:,1));
Ys1=Yssqz*AF(1);
Ys2=Yssqz*AF(2);
 
figure(18);
yyaxis left
plot(tt8,z_1(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,z_1(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_18.png','Resolution',300);
 
figure(19);
yyaxis left
plot(tt8,Ys1(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,Ys1(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_19.png','Resolution',300);
 
figure(20);
yyaxis left
plot(tt8,z_2(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,z_2(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_20.png','Resolution',300);
 
figure(21);
yyaxis left
plot(tt8,Ys2(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,Ys2(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_21.png','Resolution',300);
 
%% Question 10.b

[~,z_3] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,0,AT(1),pi/2), tt8, ICNL, opts);
[~,z_4] = ode45(@(t,z) EX_4_NLode(t,z,M,m,k,c1,c2,L,g,0,0,AT(2),pi/2), tt8, ICNL, opts);
Ys=step(Linear_system,tt8);
Yssqz= squeeze(Ys(:,:,2));
Ys3=Yssqz*AT(1);
Ys4=Yssqz*AT(2);
 
figure(22);
yyaxis left
plot(tt8,z_3(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,z_3(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_22.png','Resolution',300);
 
figure(23);
yyaxis left
plot(tt8,Ys3(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,Ys3(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_23.png','Resolution',300);
 
figure(24);
yyaxis left
plot(tt8,z_4(:,1),'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,z_4(:,2),'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('nonlinear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_24.png','Resolution',300);
 
figure(25);
yyaxis left
plot(tt8,Ys4(:,1).*180./pi,'LineWidth',2)
xlabel('$t$(s)','interpreter','Latex')
ylabel('y(m)','interpreter','Latex')
yyaxis right
plot(tt8,Ys4(:,2).*180./pi,'LineWidth',2)
ylabel('$\theta(degre)$','interpreter','Latex')
title('Linear system response to intial conditions','interpreter','Latex')
set(gca,'FontSize',12)
set(gcf,'color','w')
exportgraphics(gcf,'Figure_25.png','Resolution',300);



function dzdt = EX_4_NLode(t,z,M,m,k,c1,c2,L,g,omega,AF,AT,phi)
% computes the derivative of the states z
% z = [x; theta; dxdt; dthetadz]
F = AF*sin(omega*t+phi);
TT = AT;
 
dzdt(1,1) = z(3);
dzdt(2,1) = z(4);
dzdt(3,1) = (-1).*g+L.^(-1).*(m+2.*M+m.*cos(2.*z(2))).^(-1).*((-2).*L.*((-1).*F+ ...
k.*z(1)+c1.*z(3)+L.*m.*z(4).^2.*cos(z(2)))+2.*((-1).*TT+c2.*z(4)).*sin(z(2))); ...
dzdt(4,1) = L.^(-2).*m.^(-1).*(m+2.*M+m.*cos(2.*z(2))).^(-1).*((-2).*(m+M).*(( ...
-1).*TT+c2.*z(4))+2.*L.*m.*((-1).*F+k.*z(1)+c1.*z(3)+L.*m.*z(4).^2.*cos( ...
z(2))).*sin(z(2)));
end



function dzdt = EX_4_NLode_SDOF(t,z,m,c2,L,g)
% computes the derivative of the states z
% z = [theta; dthetadz]
TT = 0;
 
dzdt(1,1) = z(2);
dzdt(2,1) = (TT-c2*z(2)-m*g*L*sin(z(1)))/(m*L^2);
end
Editor is loading...
Leave a Comment