Solving this matrix with the correct response using MATLAB and simulink

3 次查看（过去 30 天）
Jon
Jon 2024-3-29
this is the given problem.
this is the Simulink and here is the call back for the function command:
here is my code
clc
clear
%---------------EQUATIONS_FROM_FBD------------------------------------------------------------------------------------------------------------------------------------------------
%......................................................................SUM_OF_EXTERNAL_FORCES=M*X"................................................................................
% MfXf"=-Ksf(Xf-[Xb-(L1*theta)})-Bsf(Xf'-{Xb'-(L1*theta')})+(Kf*Xf);
% Vf'=(1/Mf)[(-Ksf(Xf-(L1*theta))-(Bsf(Vf-(L1*theta'))+(Kf*Xf)] ;
% MrXr"=-Ksr(Xr+[Xb+(L2*theta)])-Bsr(Xr'+[Xb'+(L2*theta'))+(Kr*Xr) ;
% Vr'=(1/Mr)[-Ksr(Xr+(L2*theta))-Bsr(Vr+(L2*theta'))+(Kr*Xr)] ;
% MbXb"=Ksf(Xf-[Xb-(L1*theta)])+Bsf(Xf'-[Xb'-(L1*theta')])+Ksr(Xr+[Xb+(L2*theta)])+Bsr(Xr'+[Xb'+(L2*theta'))+fa(t);
% Vb'=(1/Mb)[Ksf(Xf-(L1*theta))+Bsf(Vf-(L1*omega))+Ksr(Xr+(L2*theta))+Bsr(Vr+(L2*omega))+fa(t)];
syms Ksf Xf L1 theta Bsf Vf omega Ksr Bsr Kf Kr Vr Mr Ar Xr L2 Af Ab Mf Mb fa(t) Alpha Ic L3 Xb Vb
z=Mr*Ar==-Ksr*(Xr+(Xb+(L2*theta)))-Bsr*(Vr+(Vb+(L2*omega)))+(Kr*Xr);
expand(solve(z,Ar))
ans =
y=Mf*Af==-Ksf*(Xf-(Xb-(L1*theta)))-Bsf*(Vf-(Vb-(L1*omega)))+(Kf*Xf);
expand(solve(y,Af))
ans =
x=Mb*Ab==Ksf*(Xf-(Xb-(L1*theta)))+Bsf*(Vf-(Vb-(L1*omega)))+Ksr*(Xr+(Xb+(L2*theta)))+Bsr*(Vr+(Vb+(L2*omega)))+fa(t);
expand(solve(x,Ab))
ans =
%.....................................................................SUM_OF_EXTERNAL_TORQUE=J*THETA"...............................................................................................
% Ic*theta"={-[Ksf(Xf-(L1*theta))*L1]-[Bsf(Xf'-(L1*theta'))*L1]+[Ksr(Xr+(L2*theta))*L2]+[Bsr(Xr'+(L2*theta'))*L2]+[fa(t)*L3]};
% omega'=(1/Ic){[-Ksf*Xf + Ksf((L1)^2)*theta)] + [-Bsf*Vf*L1 + Bsf*Vf*L1 + Bsf((L1)^2)*omega)] + [Ksr*Xr*L2 + Ksr((L2)^2)*theta)] + [Bsr*Vr*L2 + Bsr((L2)^2)*omega)*L2] + [fa(t)*L3]};
w=Ic*Alpha==-Ksf*(Xf-(Xb-(L1*theta)))*L1-Bsf*(Vf-(Vb-(L1*omega)))*L1+Ksr*(Xr+(Xb+(L2*theta)))*L2+Bsr*(Vr+(Vb+(L2*omega)))*L2+fa(t)*L3;
expand(solve(w,Alpha))
ans =
%.....................................................................OUTPUT_EQUATIONS..................................................................................................................
syms f_Ksf f_Kf f_Bsf f_Ksr f_Kr f_Bsr
% f_Ksf=Ksf(Xf-(-L1*theta));
% f_Kf=Kf*Xf;
% f_Bsf=Bsf(Vf-(-L1*omega));
% f_Ksr=Ksr(Xr-(L2*theta));
% f_Kr=Kr*Xr;
% f_Bsr=Bsr(Vr-(L2*omega));
f_Ksf=Ksf*(Xf-(Xb-(-L1*theta)));
A=simplify(f_Ksf);
expand(A)
ans =
f_Bsf=Bsf*(Vf-(Vb-(-L1*omega)));
B=simplify(f_Bsf);
expand(B)
ans =
f_Ksr=Ksr*(Xr-(Xb+(L2*theta)));
C=simplify(f_Ksr);
expand(C)
ans =
f_Bsr=Bsr*(Vr-(Vb+(L2*omega)));
D=simplify(f_Bsr);
expand(D)
ans =
% Mb*Vb'-Ksf(Xf-(L1*theta))-Bsf(Xf'-(L1*theta'))-Ksr(Xr+(L2*theta))-Bsr(Xr'+(L2*theta'))=fa(t);
% Ic*omega'+[Ksf(Xf-(L1*theta))*L1]+[Bsf(Vf-(L1*omega))*L1]-[Ksr(Xr+(L2*theta))*L2]-[Bsr(Vr+(L2*omega))*L2]=[fa(t)*L3]
clc
clear
%-------------------------------------------------SYSTEM_PARAMETERS----------------------------------------------------------------------------------------------------------------------------------------------------------
Ic=1356; %kg-m^2
Mb=730; %kg
Mf=59; %kg
Mr=45; %kg
Kf=23000; %N/m
Ksf=18750; %N/m
Kr=16182; %N/m
Ksr=12574; %N/m
Bsf=100; %N*s/m
Bsr=100; %N*s/m
L1=1.45; %m
L2=1.39; %m
L3=0.67; %m
%-------------------------------------------------------------------------------------Matrix's-----------------------------------------------------------------------------------------------------------------------------------
%__________________________________________________INPUT__MATRIX_________________________________________________________________________________________________________________________________________________________________
%________________________________________________________________________________________________________________________________________________________________________________________________________________________________
%......Xf...................Vf......................Xr...................Vr..................Xb...........................Vb................................theta.............................omega..............................
%________________________________________________________________________________________________________________________________________________________________________________________________________________________________
a=[ 0 1 0 0 0 0 0 0; %Xf'
((-Ksf+Kf)/Mf) (Bsf/Mf) 0 0 (Ksf/Mf) (Bsf/Mf) ((Ksf*L1)/Mf) ((Bsf*L1)/Mf); %Vf'
0 0 0 1 0 0 0 0; %Xr'
0 0 ((-Ksr+Kr)/Mr) (-Bsr/Mr) (-Ksr/Mr) (-Bsr/Mr) ((-Ksr*L2)/Mr) ((-Bsr*L2)/Mr); %Vr'
0 0 0 0 0 1 0 0; %Xb'
(Ksf/Mb) (Bsf/Mb) (Ksr/Mb) (Bsr/Mb) ((Ksr-Ksf)/Mb) ((Bsr-Bsf)/Mb) (((-Ksf*L1)+(Ksr*L2))/Mb) (((-Bsf*L1)+(Bsr*L2))/Mb); %Vb'
0 0 0 0 0 0 0 1; %theta'
((-Ksf*L1)/Ic) ((-Bsf*L1)/Ic) ((Ksr*L2)/Ic) ((Bsr*L2)/Ic) (((Ksf*L1)+(Ksr*L2))/Ic) (((Bsf*L1)+(Bsr*L2))/Ic) (((Ksr*(L2^2))+(Ksf*(L1^2)))/Ic) (((Bsr*(L2^2))+(Bsf*(L1^2)))/Ic)]; %omega'
b=[0;
0;
0;
0;
0;
(1/Mb);
0;
((L3)/Ic);];
%________________________________________________OUTPUT__MATRIX_____________________________________________________________________________________________________________________________________________________________________
c=[ (Ksf) 0 0 0 (-Ksf) 0 (-Ksf*L1) 0; %f_Ksf
Kf 0 0 0 0 0 0 0; %f_Kf
0 (Bsf) 0 0 0 (-Bsf) 0 (Bsf*L1); %f_Bsf
0 0 (Ksr) 0 (-Ksr) 0 (-Ksr*L2) 0; %f_Ksr
0 0 Kr 0 0 0 0 0; %f_Kr
0 0 0 (Bsr) 0 (-Bsr) 0 (-Bsr*L2); %f_Bsr
(Ksf/Mb) (Bsf/Mb) (Ksr/Mb) (Bsr/Mb) 0 0 (((-Ksf*L1)+(Ksr*L2))/Mb) (((-Bsf*L1)+(Bsr*L2))/Mb); %Ab
0 0 0 0 (((Ksr*L2)-(Ksf*L1))/Ic) (((Bsr*L2)-(Bsf*L1))/Ic) (((Ksr*(L2^2))-(Ksf*(L1^2)))/Ic) (((Bsr*(L2^2))-(Bsf*(L1^2)))/Ic)]; %ALPHAm
d=[0;
0;
0;
0;
0;
0;
1/Mb;
(L3/Ic)];
CSTR=ss(a,b,c,d)
CSTR =
A =
x1 x2 x3 x4 x5 x6 x7 x8
x1 0 1 0 0 0 0 0 0
x2 72.03 1.695 0 0 317.8 1.695 460.8 2.458
x3 0 0 0 1 0 0 0 0
x4 0 0 80.18 -2.222 -279.4 -2.222 -388.4 -3.089
x5 0 0 0 0 0 1 0 0
x6 25.68 0.137 17.22 0.137 -8.46 0 -13.3 -0.008219
x7 0 0 0 0 0 0 0 1
x8 -20.05 -0.1069 12.89 0.1025 32.94 0.2094 46.99 0.2975
B =
u1
x1 0
x2 0
x3 0
x4 0
x5 0
x6 0.00137
x7 0
x8 0.0004941
C =
x1 x2 x3 x4 x5 x6 x7 x8
y1 1.875e+04 0 0 0 -1.875e+04 0 -2.719e+04 0
y2 2.3e+04 0 0 0 0 0 0 0
y3 0 100 0 0 0 -100 0 145
y4 0 0 1.257e+04 0 -1.257e+04 0 -1.748e+04 0
y5 0 0 1.618e+04 0 0 0 0 0
y6 0 0 0 100 0 -100 0 -139
y7 25.68 0.137 17.22 0.137 0 0 -13.3 -0.008219
y8 0 0 0 0 -7.161 -0.004425 -11.16 -0.01257
D =
u1
y1 0
y2 0
y3 0
y4 0
y5 0
y6 0
y7 0.00137
y8 0.0004941
Continuous-time state-space model.
Model Properties
sys_tf=tf(CSTR);%System transfer function
impulse(sys_tf)
sys_zpg=zpk(CSTR);%system zero pole gain
impulse(sys_zpg)
syms s real;
i=eye(8);
h=simplify(c*inv((s*i)-a)*b+d)
h =
h(1)
ans =
simplify(h(1))
ans =
%-------------------------SIMULATION-----------------------------------------------
sim('Pagels_Jon_Project_Matrix_simulation')
fat=10*(exp(-5*t));
Xf=q(:,1);
Vf=q(:,2);
Xr=q(:,3);
Vr=q(:,4);
Xb=q(:,5);
Vb=q(:,6);
theta=q(:,7);
omega=q(:,8);
%----------------------------FORCES------------------------------------------------
KSF=y(:,1);
KF=y(:,2);
BSF=y(:,3);
KSR=y(:,4);
KR=y(:,5);
BSR=y(:,6);
Ab=y(:,7);
Alpha_m=y(:,8);
%----------------------------FIGURES-------------------------------------------------
figure(1)
set(gcf,'color','white')
subplot(3,1,1);
plot(t,Xb,'-b','linewidth',0.5);grid
title('Displacement of Body Linear vs. Time')
xlabel ('Time (s)') ;
ylabel ('Displacement (m)');
subplot(3,1,2);
plot(t,theta,'-k','linewidth',0.5);grid
title('Displacement of Body Angular vs. Time ')
xlabel ('Time (s)') ;
ylabel ('Displacement (rad)');
subplot(3,1,3);
plot(t,fat,'-m','linewidth',0.5);grid
title('Impulse Force from Road vs. Time')
xlabel ('Time (s)') ;
ylabel ('Force (N/s)');
figure(2)
set(gcf,'color','white')
subplot(3,1,1);
plot(t,Xf,'-c','linewidth',0.5);grid
title('Displacement of Front Axel vs. Time')
xlabel ('Time (s)') ;
ylabel ('Displacement (m)');
subplot(3,1,2);
plot(t,KSF,'-m','linewidth',0.5);grid
hold on
subplot(3,1,2);
plot(t,KF,'-y','linewidth',0.5);grid
title('Front Spring Force vs. Time ')
xlabel ('Time (s)') ;
ylabel ('Spring Force-Front (N)');
hold off
subplot(3,1,3);
plot(t,BSF,'-g','linewidth',0.5);grid
title('Front Damper Force vs. Time')
xlabel ('Time (s)') ;
ylabel ('Damper Force-Front (N/s)');
figure(3)
set(gcf,'color','white')
subplot(3,1,1);
plot(t,Xr,'-c','linewidth',0.5);grid
title('Displacement of Rear Axel vs. Time')
xlabel ('Time (s)') ;
ylabel ('Displacement (m)');
subplot(3,1,2);
plot(t,KSR,'-m','linewidth',0.5);grid
hold on
subplot(3,1,2);
plot(t,KR,'-y','linewidth',0.5);grid
title('Rear Spring Force vs. Time ')
xlabel ('Time (s)') ;
ylabel ('Spring Force-Rear (N)');
hold off
subplot(3,1,3);
plot(t,BSR,'-g','linewidth',0.5);grid
title('Rear Damper Force vs. Time')
xlabel ('Time (s)') ;
ylabel ('Damper Force-Rear (N/s)');
I am not seeing the response in the springs and dampers everything plots as an exponential curve, I'd expect to see some oscilation. I'm not sure if I'm entering something incorrectly.

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programmatic Tuning 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by