1 行驶动力学计算机建模，仿真及主动悬架
控制器设计
1.1建立M文件
%车辆模型参数
mb=320; %簧载质量(kg) mw=40; %非簧载质量(kg) ks=20000; %悬架刚度（N/m） kt=200000; %轮胎刚度（N/m） swsc=+-100; %悬架工作空间（mm） %仿真路面输入参数
g0=5*10^-6; %路面不平度系数（m^3/cycle） u=20; %车速（m/s）
f0=0.1; %下截止频率（Hz） %性能指标加权系数
q1=80000; %轮胎动位移 q2=5; %悬架动行程 q3=1; %车身加速度 %系统输入
x0=wgn(10001,1,20); %生成高斯白噪声（采样点数 列数 功率） t=0:0.005:50; %仿真时间
Tw=[t' x0]; %系统输入序列[仿真时间序列 白噪声序列] %系统状态方程矩阵
A=[0 0 -ks/mb ks/mb 0;0 0 ks/mw -(ks+kt)/mw kt/mw;1 0 0 0 0;0 1 0 0 0;0 0 0 0 -2*pi*f0]; B=[1/mb;-1/mw;0;0;0];
F=[0;0;0;0;2*pi*(g0*u)^0.5]； %二次最优控制器设计
Q=[0 0 0 0 0;0 0 0 0 0;0 0 q2+ks^2/mb^2 -q2-ks^2/mb^2 0;0 0 -q2-ks^2/mb^2 q1+q2+ks^2/mb^2 -q1;0 0 0 -q1 q1]; R=1/mb^2;
N=[0;0;-ks/mb^2;ks/mb^2;0];
[K,S,E]=LQR(A,B,Q,R,N); %K为最优控制反馈增益矩阵 以上程序执行后可得到的结果为：
A=[0 0 -62.5 62.5 0;0 0 500 -5500 5000;1 0 0 0 0;0 1 0 0 0;0 0 0 0 -0.6]; B=[0.0031;-0.025;0;0;0]; F=[0;0;0;0;0.0628];
K =[711.88 -1241.5 -19284 -2038.5 208]; Tw[时间序列 白噪声幅值序列]
1.2 建立simulink模型
设置输入模块：（Tw序列）
设置状态空间方程参数：（其中B矩阵为M执行结果中B与F矩阵的和矩阵维数：4×2）
设置增益模块：（K矩阵）
1.3 执行结果
根据以上设置执行仿真，可得如下结果：
车身加速度（m/s^2）-时间（s）
悬架动行程(m)-时间（s）
轮胎动位移(m)-时间（s）
路面位移输入(m)-时间（s）
1.4 结果分析
根据上图可分别导出车身加速度BA，悬架动行程SWS，轮胎动位移DTD的数据。
三者都是6136×2的矩阵，第一列为幅值对应的时间，第二列为三者的幅值，对其第二列分别求均方根值，程序为：
Y=sqrt(sum(X(:,2).^2)/6136) 执行程序，可得主动悬架性能指标的均方根值如表1。
在相同仿真条件下，可对一个被动系统进行分析，取悬架刚度Ks=22000N/m，阻尼系数Cs=1000N*s/m。其他参数均与主动系统相同，计算可得被动系统性能指标如表1。
表1 主动悬架与被动悬架性能指标均方根值的比较 性能指标 车身加速度 悬架动行程 轮胎动位移 单位 m/s2 mm mm 均方根值 主动悬架 1.4288 44.3 5.9 被动悬架 1.7533 17.6 6.0 由表可见，在轮胎动位移基本相同的情况下，所设计的最优主动悬架显著地降低了车身的垂向振动加速度，与被动悬架系统相比，其均方根值减少了18.5%。而且，主动悬架系统的悬架动行程被很好的控制在设计的范围（±100mm）内，且比被动悬架的行程大了许多，意味着许用的悬架工作空间得到了更充分的利用。
2 操纵动力学计算机建模及分析
2.1操纵动力学模型建立
%模型参数
m=[2045 1008]; %质量
Izz=[5428 1031]; %横摆转动惯量 a=[1.488 1.234]; %前轴到质心距离 b=[1.712 1.022]; Caf=[77850 117440]; Car=[76510 144930]; uc=40;
%建立状态方程矩阵 %矩阵A
A11=-(Caf+Car)./(m*uc);
A12=-(a.*Caf-b.*Car)./(m*uc)-uc; A21=-(a.*Caf-b.*Car)./(Izz*uc);
A22=-(a.^2.*Caf+b.^2.*Car)./(Izz*uc); A1=[A11(1) A12(1);A21(1) A22(1)]; A2=[A11(2) A12(2);A21(2) A22(2)]; %矩阵B B11=Caf./m; B22=a.*Caf./Izz; B1=[B11(1);B22(1)]; B2=[B11(2);B22(2)]; %矩阵CD C=[0 1]; D=0;
执行结果为：
A1 =[-1.8870 -39.8149;0.06981 -1.8267] A2 =[-6.5072 -39.9207;0.0775 -8.0070] B1=[38.0685;21.3413] B2=[116.5079;140.5635] C=[0 1]; D=0;
%后轴到质心距离 %前轮侧偏刚度 %后轮侧偏刚度 %速度
2.2 时域分析
建立Simulink模型
阶跃输入参数设置：初始时刻（t=0s）输入为0rad，t>0s时输入为0.0058rad。
别克模型状态方程参数设置：(矩阵A1 B1 C D)
法拉利模型状态方程参数设置：（矩阵A2 B2 C D）
运行结果：
横摆角速度/[rad/s]
时间/s
角阶跃输入下的横摆角速度时域响应(紫色为法拉利，黄色为别克)
由图可知，在同样的转向盘转角输入下，法拉利跑车的瞬态响应比别克轿车的要好，主要体现在较短的响应时间，较小的超调量以及更好的阻尼特性等。
2.3频域分析
执行程序如下： figure()
sys1=ss(A1,B1,C,D); sys2=ss(A2,B2,C,D); bode(sys1,sys2); 结果为：
由图可知，对于幅频响应在同一行驶车速下，法拉利跑车的响应带宽大于别克轿车的响应带宽，从而也说明了前者具有更好的频率响应特性，对于相频特性，法拉利跑车的系统响应之后要比别克车的少，系统延迟小。
2.3稳定性分析
求解车辆在车速uc=15~60m/s的特征值 执行程序如下： figure() u=15:5:60;
for j=1:length(u);
A=[A11(1)*uc/u(j) (A12(1)+uc)*uc/u(j)-u(j);A21(1)*uc/u(j) A22(1)*uc/u(j)]; plot(real(eig(A)),imag(eig(A)),'go'); hold on; end
for j=1:length(u);
A=[A11(2)*uc/u(j) (A12(2)+uc)*uc/u(j)-u(j);A21(2)*uc/u(j) A22(2)*uc/u(j)]; plot(real(eig(A)),imag(eig(A)),'bx'); hold on; end
结果为：
由图可知，随着车速的增加，系统的特征值由复平面的左侧不断向虚轴靠近，因此系统越来越趋向于不稳定，对两种车型的比较可见，他们的特征跟随车速变化而变化的趋势也不同，由于法拉利跑车具有相对较大的稳定裕度，因而其特征根位置与别克轿车相比更远离虚轴，而别克轿车的转向特性的稳定裕量较小且容易失稳。
3 半车模型的推导分析
3.1模型建立
%参数设置
mhb=690; %二分之一车身质量（kg） Ihp=1222; %二分之一转动惯量(kg*m2) mwf=40.55; %前轮非簧载质量(kg) mwr=45.4; %后轮非簧载质量(kg) Ktf=192000; %前轮胎刚度(kN/m) Ktr=192000; %后轮胎刚度(kN/m) Ksf=17000; %前悬架刚度(kN/m) Ksr=22000; %后悬架刚度(kN/m)
Csf=1500; %前悬架阻尼系数(kN*s/m) Csr=1500; %后悬架阻尼系数(kN*s/m) a=1.25; %车身质心至前轴距离（m） b=1.51; %车身质心至后轴距离（m） L=2.76; %轴距（m） i=(-1)^0.5; %虚数
u=20; %车速(m/s)
t=L/u; %前后轮输入时间间隔（s）
%建立动力学模型 w=sym('2*pi*f') c1=1/mhb+a^2/Ihp; c2=1/mhb-a*b/Ihp; c3=1/mhb+b^2/Ihp;
A1=[Csf*i*w+Ksf+Ktf-mwf*w^2;c1*(Ksf+Csf*i*w);0;c2*(Ksf+Csf*i*w)]; A2=[-Csf*i*w-Ksf;w^2-c1*(Ksf+Csf*i*w);0;-c2*(Ksf+Csf*i*w)];
A3=[0;c2*(Ksr+Csr*i*w);Csr*i*w+Ksr+Ktr-mwr*w^2;c3*(Ksr+Csr*i*w)]; A4=[0;-c2*(Ksr+Csr*i*w);-Csr*i*w-Ksr;w^2-c3*(Ksr+Csr*i*w)]; B=[Ktf;0;Ktr/exp(i*w*t);0];
G1=det([B A2 A3 A4])/det([A1 A2 A3 A4]); G2=det([A1 B A3 A4])/det([A1 A2 A3 A4]); G3=det([A1 A2 B A4])/det([A1 A2 A3 A4]); G4=det([A1 A2 A3 B])/det([A1 A2 A3 A4]);
3.2仿真程序及结果
%求悬架动行程增益的频响函数 Gswsf=G1-G2;
Gswsr=G3-G4; fsf=inline(Gswsf); fsr=inline(Gswsr); f=0:0.01:15; figure();
plot(f,abs(fsf(f)),'b',f,abs(fsr(f)),'g'); 结果：
%求垂向加速度增益的频响函数 Gvaf=-w^2*G2; Gvar=-w^2*G4; fvf=inline(Gvaf); fvr=inline(Gvar); f=0:0.01:15; figure();
plot(f,abs(fvf(f)),'b',f,abs(fvr(f)),'g'); 结果：
%求轮胎动载荷增益的频响函数 Gdtlf=Ktf*(G1-1); Gdtlr=Ktr*(G3-1); fdf=inline(Gdtlf); fdr=inline(Gdtlr); f=0:0.01:15; figure();
plot(f,abs(fdf(f)),'b',f,abs(fdr(f)),'g'); 结果：
%求质心处垂向加速度
Gvac=( Gvar*a+ Gvaf*b)/(a+b) fvc= inline(Gvac); f=0:0.01:15; figure();
plot(f,abs(fvc(f))); 结果：
%求俯仰角增益
Gp=(Gvar-Gvaf)/(a+b)/-w^2; frl= inline(Gp); f=0:0.01:15; figure();
plot(f,abs(frl(f))); 结果：
结果分析：
与单轮车辆模型相比，半车模型不同的地方在于，后轮的输入将会比前轮延迟时间t=L/u。故前后两轮输入存在相位差，前后轮输入的共同作用决定了质心垂向加速度和俯仰角的频率响应特性，且影响规律并不是简单的线性叠加。
4整车模型的推导及分析
4.1模型建立
%设置参数
mb=1380; %整车质量（kg）
Ip=2440; %俯仰转动惯量(kg*m2) Ir=380; %侧倾转动惯量(kg*m2)
mwa=40.5; %A处的簧下质量质量（kg） mwb=40.5; %B处的簧下质量质量（kg） mwc=45.4; %C处的簧下质量质量（kg） mwd=45.4; %D处的簧下质量质量（kg） kta=192000; %A处的轮胎刚度（N/m） ktb=192000; %B处的轮胎刚度（N/m） ktc=192000; %C处的轮胎刚度（N/m） ktd=192000; %D处的轮胎刚度（N/m） ksa=17000; %A处的悬架刚度（N/m） ksb=17000; %B处的悬架刚度（N/m） ksc=22000; %C处的悬架刚度（N/m） ksd=22000; %D处的悬架刚度（N/m） Csa=1500; %A处的阻尼系数（N*s/m） Csb=1500; %B处的阻尼系数（N*s/m） Csc=1500; %C处的阻尼系数（N*s/m） Csd=1500; %D处的阻尼系数（N*s/m） a=1.25; %车身质心至前轴距离（m） b=1.51; %车身质心至后轴距离（m） L=2.76; %轴距（m） i=(-1)^0.5; %虚数
u=20; %车速(m/s)
t=L/u; %前后轮输入时间间隔（s） Bf=1.48; %前轴轮距 Br=1.48; %后轴轮距
%建立动力学模型 ea=Csa*i*w+ksa; eb=Csb*i*w+ksb; ec=Csc*i*w+ksc; ed=Csd*i*w+ksd;
A1=[ea+kta-mwa*w^2;0;0;0;ea;-a*ea;Bf/2*ea]; A2=[0;eb+ktb-mwb*w^2;0;0;eb;-a*eb;-Bf/2*eb]; A3=[0;0;ec+ktc-mwc*w^2;0;ec;b*ec;Br/2*ec]; A4=[0;0;0;ed+ktd-mwd*w^2;ed;b*ed;-Br/2*ed];
A5=[-ea;-eb;-ec;-ed;mb*w^2-(ea+eb+ec+ed);a*(ea+eb)-b*(ec+ed);Bf/2*(-ea+eb)+Br/2*(-ec+ed)]; A6=[a*ea;a*eb;-b*ec;-b*ed;a*(ea+eb)-b*(ec+ed);Ip*w^2-a^2*(ea+eb)+b^2*(ec+ed);a*Bf/2*(ea-eb)-b*Br/2*(ec-ed)];
A7=[-Bf/2*ea;Bf/2*eb;-Br/2*ec;Br/2*ed;Bf/2*(-ea+eb)+Br/2*(-ec+ed);a*Bf/2*(ea-eb)+b*Br/2*(-ec+ed);-Bf^2/4*(ea+eb)-Br^2/4*(ec+ed)];
B=[Kta;Ktb/exp(i*w*2*t);Ktc/exp(i*w*t);Ktd/exp(i*w*3*t);0;0;0]; G1=det([B A2 A3 A4 A5 A6 A7])/det([A1 A2 A3 A4 A5 A6 A7]); G2=det([A1 B A3 A4 A5 A6 A7])/det([A1 A2 A3 A4 A5 A6 A7]); G3=det([A1 A2 B A4 A5 A6 A7])/det([A1 A2 A3 A4 A5 A6 A7]); G4=det([A1 A2 A3 B A5 A6 A7])/det([A1 A2 A3 A4 A5 A6 A7]); G5=det([A1 A2 A3 A4 B A6 A7])/det([A1 A2 A3 A4 A5 A6 A7]); G6=det([A1 A2 A3 A4 A5 B A7])/det([A1 A2 A3 A4 A5 A6 A7]); G7=det([A1 A2 A3 A4 A5 A6 B])/det([A1 A2 A3 A4 A5 A6 A7]);
4.2仿真程序及结果
%侧倾角增益 Groll=G7;
frl= inline(Groll); f=0:0.01:15; figure();
plot(f,abs(frl(f))); 结果：
%俯仰角增益 Gp=G6;
fp= inline(Gp);
f=0:0.01:15; figure();
plot(f,abs(fp(f))); 结果：
%质心加速度增益 Gvc=-w^2*G5; fvc= inline(Gvc); f=0:0.01:15; figure();
plot(f,abs(fvc(f))); 结果：
由图可知，对于整车模型而言，当四轮输入均不相同时，车身某点的响应可以认为是四者同时作用的主模态各自频率响应的叠加。故其频响分析较为复杂，但当前轮输入相同，后轮输入相同，且前后轮间输入存在时间延迟时，即转化为半车模型，通过设置相应的输入，可以通过频响函数进行验证。