Untitled

clear all;
close all;

%% question 6
 
% G(s) = 1 / (4s^2 + 8 s + 10e3)
 
sys_tf = tf([1],[4 8 10e3]);
sys_ss = ss(sys_tf);
 
step(sys_tf)
 
[Y_tf_step,tt] = step(sys_tf,linspace(0,7,1e3).');
figure(2);clf
plot(tt,Y_tf_step)
xlabel t(s)
ylabel y(m)
 
TF_step_info = stepinfo(sys_tf,'SettlingTimeThreshold',0.02,'RiseTimeLimits',[0.0 1.0]);
% To extract the data from the structure:
tr = TF_step_info.RiseTime % rise time
ts = TF_step_info.SettlingTime % settling time
Mp = TF_step_info.Overshoot % Over shoot
 
 
sys_tf2 = tf(sys_ss)
% you can check what additional parameters are return by stepinfo by:
% > doc stepinfo, or explore it in the workspace.
 
%% question 7 ?????
%%steel
% Es=200e9;
% Rs=7800;
% As=0.2;
% Bs=1e-4;
 
%%aluminium
% Es=70e9;
% Rs=2700;
% As=3;
% Bs=2e-4;
 
%%Epoxy
Es=5e9;
Rs=1000;
As=2;
Bs=5e-4;
 
L = optimvar('L');
h = optimvar('h');
b= optimvar('b');
 
k=(b*h^3*Es)/(12*3*L^3);
m=Rs*b*h*L;
c=As*m+Bs*k;
omega0 = sqrt(k/m);
ksi=c/(2*(k*m)^0.5);
Mps=exp((-ksi*pi)/sqrt(1-ksi^2))*100;
Trs=(pi-atan(sqrt(1-ksi^2)/ksi))/(omega0*sqrt(1-ksi^2));
Tss= 4/(omega0*ksi);
prob = optimproblem;
prob.Objective = 0.2*Mps+100*Trs+10*Tss;
prob.Constraints.cons1 = L<= 0.5;
prob.Constraints.cons2 = L>= 0.05;
prob.Constraints.cons3 = h<= L/10;
prob.Constraints.cons4 = h >= L/30;
prob.Constraints.cons5 = b<= L/4;
prob.Constraints.cons6 = b >= 5*h;
x0.L = 0.05;
x0.h = 0.05/10;
x0.b=0.1;
sol = solve(prob,x0)
 
 
L=0.05;
h=0.0025;
b=0.0125;
k=(b*h^3*Es)/(12*3*L^3);
m=Rs*b*h*L;
c=As*m+Bs*k;
omega0 = sqrt(k/m);
ksi=c/(2*(k*m)^0.5);
Mps=exp(-ksi*pi/sqrt(1-ksi^2))*100;
Trs=(pi-atan(sqrt(1-ksi^2)/ksi))/(omega0*sqrt(1-ksi^2));
Tss=4/(omega0*ksi);
FINAL_score= 0.2*Mps+100*Trs+10*Tss
 
%% question 8 and 9
Ee=5e9;
Re=1000;
ae=2;
Be=5e-4;
kfinal=(Ee*b*h^3)/(36*L^3);
mfinal=Re*b*h*L;
cfinal=ae*m+Be*kfinal;
M=1;
cd1=5;
cd2=200;
cd3=1000;
a=5000;
T=0.01;
 
% positive impulse
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
tt = linspace(0,0.2,1e4);
IC1 = [0; 0; 0; 0];
 
[~, z_1] = ode45(@(t,z) ODE_EQUATION(t,z,M,mfinal,kfinal,cfinal ...
,cd1,cd2,cd3,a,T), tt, IC1, opts);
 
figure(3)
hold on
plot(tt,z_1(:,1),'linewidth',0.5);
plot(tt,z_1(:,2),'linewidth',0.5);
xlabel('$time$(secs)','interpreter','Latex')
ylabel('$x$(m)','interpreter','Latex')
leg=legend('Mass M', 'Mass m');
 
set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
hold off
 
 
%%relative error
figure(4)
hold on
plot(tt,((z_1(:,2)-z_1(:,1))./z_1(:,1))*100)
maximum_rel_error=abs(min(((z_1(:,2)-z_1(:,1))./z_1(:,1))*100))
xlabel('$time$(secs)','interpreter','Latex')
ylabel('Relative Error (%)')
set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
hold off
% negative impulse
[~, z_2] = ode45(@(t,z) ODE_EQUATION(t,z,M,mfinal,kfinal,cfinal ...
,cd1,cd2,cd3,-a,T), tt, IC1, opts);
 
figure(5)
hold on
plot(tt,z_2(:,1),'linewidth',0.5);
plot(tt,z_2(:,2),'linewidth',0.5);
xlabel('$time$(secs)','interpreter','Latex')
ylabel('$x$(m)','interpreter','Latex')
leg=legend('Mass M', 'Mass m');
 
set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
hold off
 
%relative error
 
figure(6)
hold on
plot(tt,((z_2(:,2)-z_2(:,1))./z_2(:,1))*100)
xlabel('$time$(secs)','interpreter','Latex')
ylabel('Relative Error (%)')
set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
hold off
 
 
 
function dzdt = ODE_EQUATION(t,z,M,m,k,c,cd1,cd2,cd3,a,T)
% computes the derivative of the states z
if(t<T/2)
A=a;
else
A=0;
end
 
P = A*sin(2*pi*t/T);
 
% z = [x; theta; dxdt; dthetadz]
 
dzdt(1,1) = z(3);
dzdt(2,1) = z(4);
dzdt(3,1) = (z(2)-z(1))*k/M +c*(z(4)-z(3))/M +P/M -(cd1*z(3)+cd2*z(3)^2 ...
+cd3*z(3)^3)/M;
dzdt(4,1) = (z(1)-z(2))*k/m+(z(3)-z(4))*c/m;
end