行驶动力学建模、仿真及主动悬架控制器设计

1、目录1. 计算机仿真系统模型的建立-.1-2. LOG控制器设计-.2-3. 计算实例-.3-4. MATLAB仿真过程-.4-5. 半车模型建模及仿真-.8-5.1随机线性最优控制-.9-5.2预瞄控制-11-5.3结果比较-.12-（6）（m3/cycle）；u为车辆前进速度（m/s）；Hz）。式中,，为系统状态矢量;W=（w（t）,为高斯白噪声输入矩以单轮车辆模型为例，介绍行驶动力学计算机建模、仿真分析以及利用线性二次最优控制理论进行主动悬架LQG控制器设计过程。1.计算机仿真系统模型的建立根据图7所示的主动悬架单轮车辆模型，运用牛顿运动定律，建立系统的运动方程，即:mbXbUa心区Xw

2、）Xg)(5)mwXWUaKs（xbXw）Kt（xw这里，采用一个滤波白噪声作为路面输入模型，即：xg（t）2f°Xg（t）2G0uw（t）式中，xg为路面垂向位移（m）;GO为路面不平度系数w为数字期望为零的高斯白噪声；f0为下截止频率（结合式（4）、式（5）和式（6），将系统运动方程和路面输入方程写成矩阵形式，即得出系统的空间状态方程：(7)XAXBUFW阵；U=（Ua（t）,为输入控制矩阵;00KsmbKs00Amw100Ks0mbKtKsKamwlwB000002f。1mb10mw00F0000；2,o；u2.LOG控制器设计车辆悬架设计中的主要指标包括：的车身垂向振动加速度

3、；影响车身姿态且与轮胎布置有关的悬架动行程。因此，制器设计中的性能指标J即为轮胎动位移、悬架动行程和车身垂向振动加速度的加权平方和在时域T内的积分值，其表达式为：TJ何1厲风Xg(t)2珏区Xw(t)2q3Xb2(t)dtT0代表轮胎接地性的轮胎动载荷;代表轮胎舒适性LQG控(8)式中，q1、q2和q3分别为轮胎动位移、悬架动行程和车身垂向振动加速度的加权系数。加权系数的选取决定了设计者对悬架性能的倾向，如对车身垂向振动加速度项选择较大的权值，则考虑更多的是提高车辆操纵稳定性。为方便起见，值q3=1。将性能指标J的表达式(8)改写成矩阵形式，即：T(XtQXUtRU2XTNU)dt这里取车身垂

4、向振动加速度的加权式中,JTiml(9)q2q20K22mb2mbq20Kf2mbqiq2KS22mbqiqii2mbKaKs0当车辆参数值和加权系数值确定后，最优控制反馈增益矩阵可有黎卡提(求出，其形式如下：Riccati)方程PAAtP(PBN)R1(BtPNt)Q0(10)最优反馈控制增益矩阵KBTPNT，由车辆参数和加权系数决定。根据任意时刻的反馈状态变量X（t）,就可得到t时刻作动器的最优控制力Ua,即：(11)Ua(t)KX(t)3. 计算实例这里，以某轿车的后悬架为例，给出一个完整的计算实例，包括车辆模型参数、仿真路面输入参数、控制器的设计参数以及计算结果。此例中车辆以20m/s

5、的速度在某典型路面上行驶，仿真时间T=50s。计算中输入的各参数及数值详见表2。表2单轮车辆模型仿真输入参数值车辆模型参数符号单位数值黄载质量mbKg320非簧载质量mwKg40悬架刚度KsN/m20000轮胎刚度KtN/m200000悬架工作空间SWScmm100仿真路面输入参数符号单位数值路面不平度系数G0m3/cycle5.0x10-6车速Um/s20下截止频率f0Hz0.1性能指标加权参数符号单位数值轮胎动位移q1一80000悬架动仃程q2一5车身加速度q3一1仿真计算中以式（6）所示的滤波白噪声作为路面输入模型。白噪声的生成可直接调用MATLAB函数WGN（M,N,P）（此函数需要安

6、装信号处理工具箱Communicationstoolbox），其中M为生成矩阵的行数，N为列数，P为白噪声的功率（单位为dB）。本例中取M=10001，N=1，P=20。这意味着仿真计算中去一条白噪声，共10001个采集点，噪声强度为20dB。设定采样时间为0.005s、车速为20m/s时，相当于仿真路面长度为1000m，仿真时间为50s。根据建立的系统状态方程式（7）及最优化性能指标函数式（9），利用已知的矩阵A、B、Q、R、N，调用MATLAB中的线性二次最优控制器设计函数K，S,E=LQR（A，B，Q,R，N），即可完成最优主动悬架控制器的设计。输出的结果中，K为最优控制反馈增益矩阵，S

7、为黎卡提方程的解，E为系统闭环特征根。根据表2给出的仿真输入参数，本例中求得的最优反馈增益矩阵K为:K=（711.88-1241.5-19284-2038.520864）同时，还得到了黎卡提方程的解：2.45590.02892.47458.66077.3090.02890.48860.02987.52627.23642.47450.02984.97448.67545.10338.66077.52628.67542710.12700.47.3097.23645.10332700.42693.7在相同的仿真条件下，可将所设计的主动悬架系统与一个被动系统进行对比分析。在被动悬架系统中，取悬架刚度Ks

8、=22000N/m，阻尼系数Cs=1000NS/m。除此之外，其他输入参数值均与主动悬架系统完全相同。4. MATLAB仿真过程1) 生成路面输入模型代码如下：a=wgn(10001,1,20);t=0:0.005:50;road_file(:,1)=t'road_file(:,2)=a;saveroad_fileroad_file2) 参数输入代码如下：loadroad_file.mat%载入路面数据模型Ks=22000;mb=320;Kt=200000;mw=40;f0=0.1;G0=0.000005;u=20;Kb=20000;Ks1=22000;Cs=1000;%输入仿真有关参

9、数A=0,0,-Ks/mb,Ks/mb,0;%建立主动悬架的状态矩阵0,0,Ks/mw,(-Kt-Ks)/mw,Kt/mw;1,0,0,0,0;0,1,0,0,0;0,0,0,0,-2*pi*f0;A1=-Cs/mb,Cs/mb,-Ks1/mb,Ks1/mb,0;%建立被动悬架的状态矩阵Cs/mw,-Cs/mw,Ks1/mw,(-Kt-Ks1)/mw,Kt/mw;1,0,0,0,0;0,1,0,0,0;0,0,0,0,-2*pi*f0;B=1/mb,0;-1/mw,0;0,0;0,0;0,2*pi*sqrt(G0*u);B1=0,0;0,0;0,0;0,0;0,2*pi*sqrt(G0*u);

10、C=1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,1,0;0,0,0,0,1;D=0,0;0,0;0,0;0,0;0,0;K=711.88,-1241.5,-19284,-2038.5,20864;K1=0,0,0,0,0;3）用Simulink创建仿真框图?状态变量XXb,Xw,Xb,Xw,Xg输入与系统模块，如下图：Gain输出模块，如下图:>du/dtDerivative>O悬架动行程>simoutToWorkspace轮胎动位移踣面输入d?车身加速度(xb)dt悬架动行程xbxw轮胎动位移xwxg路面输入xg整体程序框图如下:Gainroad

11、jle)C-y-Cif+Du錯出况二r呈FromWort(5p3ce>du.iirDeiVHlivt!>O1x'-Aa*0uy=SDuSte-Spa-2tJGan1dj'd:D出ivaciwd>O4)结果分析可以直接通过双击scope查看输出的波形图，为更好比较主动悬架与被动悬架的差别，下面通过输出到workspace的状态变量编程绘图并计算均方根值。代码如下：%绘制车身加速度曲线，并计算均方根值%ba-主动悬架车身加速度%ba1-被动悬架车身加速度ba=diff(X.data(:,1)./diff(X.time);ba1=diff(X1.data(:,1).

12、/diff(X1.time);subplot(2,1,1)plot(X.time(1:end-1),ba)subplot(2,1,2)plot(X1.time(1:end-1),ba1)BA=norm(ba,2)./(length(ba).A0.5);BA1=norm(ba1,2)./(length(ba1).A0.5);%绘制悬架动行程曲线，并计算其均方根值%sws-主动悬架动行程%sws1-被动悬架动行程figure。sws=X.data(:,3)-X.data(:,4);sws1=X1.data(:,3)-X1.data(:,4);subplot(2,1,1)plot(X.time,sw

13、s)subplot(2,1,2)plot(X.time,sws1)SWS=norm(1000*sws,2)./(length(sws).A0.5);SWS1=norm(1000*sws1,2)./(length(sws1).A0.5);%绘制轮胎动位移曲线，并计算其均方根值%dtd-主动悬架动位移%dtd1-被动悬架动位移figure。dtd=X.data(:,4)-X.data(:,5);dtd1=X1.data(:,4)-X1.data(:,5);subplot(2,1,1)plot(X.time,dtd)subplot(2,1,2)plot(X.time,dtd1)DTD=norm(10

14、00*dtd,2)./(length(dtd).A0.5);DTD仁norm(1000*dtd1,2)./(length(dtd1).A0.5);结果如下： 车身加速度曲线円度速加身车1050-5-10IJI1'11/Bl1r：-.佃MrjIE*i卜*|.Jl鬧"J.4和1.i'>p111*!卜T；1k丁/1rd'Vtf忙1ifj主动悬架车身加速度01234567时间被动悬架车身加速度8910卩度速加身车20100-10-20FiJ>lII，'山“MiHr-机丫21'flft卩1帀H'fl1l01234567时间8910 悬

15、架动行程曲线卩移位b1ft机1*J1丹fiJk.Af旨1L尸*e/j«J"f.飞10.10.050-0.05-0.112主动悬架动行程34567时间被动悬架动行程8910卩移位00.60.40.20-0.2ji;1ik,ar1Ji:,丄叽I,Ij.fl1-If''.'r4'-A1tIi123456时间789100轮胎动位移i/i'ili'T皿.1:l'ihKJ11'h'Ij.J,1卩11-叫N忖1*，卩Lj1-.:f1I1'0.040.020-0.02-0.04012主动悬架轮胎动位移34567

16、时间被动悬架轮胎动位移89100.10.050-0.05-0.13456时间78910012主动悬架与被动悬架性能指标均方根值比较性能指标主动悬架被动悬架车身加速度BA(m/s2)1.514.60悬架动行程SWS(mm)34.4267.10轮胎动位移DTD(mm)5.8734.575. 半车模型建模及仿真半车模型建模悬架的半车模型有四个自由度，可选取前后轮垂直位移、前后悬架与车身连接处垂直位移四个自由度，写出状态空间方程:Z4,Z3,Z2,Zi,Z4,Z3,Z2,Zi,Xio,X20作为系统状态变量,状态空间方程如下:X(t)AX(t)BU(t)FW(t)UafUar,为控制输入矩阵，即前后悬

17、架作动器力;WW1,为路面输入模型的高斯白噪声输入。w2LQG控制器选取如下二次型性能指标：TimTo(乙Z01)2q2(Z2Z1)2?21Z2q3(Z3Z02)2?22Z4dtq1,q3,前后悬架轮胎动位移加权系数;q2,q4,前后悬架动行程加权系数；1,2,车身加速度加权系数。写成矩阵形式:utru2XtNU)dt1Tlim(XtQXTT10由黎卡提方程求出反馈矩阵K，则UKX，因此状态空间方程可写成:?X(t)(ABK)X(t)FW(t)(ABK)-状态矩阵F输入矩阵由此形式可方便求出其时域与频域响应。5.1随机线性最优控制路面模型本例仿真车速为20m/s，轴距为2.8m，因此滞后时间为

18、0.14s。若取仿真时间T为20s，采样时间间隔为0.005s，因而仿真点数4028。程序代码如下：%生成路面模型a=wgn(4029,1,20);t=0:0.005:20.14;r=a,t'saveroad_file%车身模型参数输入clearclcmb=690;I=1222;mwf=40;mwr=45;Ksf=17000;Ksr=22000;Ktf=200000;Ktr=200000;a=1.3;b=1.5;%仿真路面参数输入G0=5e-6;u=20;f0=0.1;%性能指标加权系数q1=80000;q2=100;q3=80000;q4=100;%半车悬架模型建模，计算A、B、F矩

19、阵a1=1/mb+bA2/I;a2=1/mb-a*b/I;a3=1/mb+aA2/I;A=0000a1*Ksr-a1*Ksra2*Ksf-a2*Ksf00;0000-Ksr/mwr(Ksr-Ktr)/mwr00Ktr/mwr0;0000a2*Ksr-a2*Ksra3*Ksf-a3*Ksf00;000000-Ksf/mwf(Ksf-Ktf)/mwf0Ktf/mwf;1000000000;0100000000;0010000000;0001000000;00000000-2*pi*f00;000000000-2*pi*f0;B=a2a1;0-1/mwr;a3a2;-1/mwf0;zeros(6,2

20、);F=zeros(8,2);2*pi*sqrt(G0*u)0;02*pi*sqrt(G0*u);%LQG控制器设计，计算Q、R、N矩阵并求出反馈矩阵Kb1=a3A2+a2A2;b2=a2*a3+a1*a2;b3=a2A2+a1A2;sq=q4+KsrA2*b3-q4-KsrA2*b3Ksf*Ksr*b2-Ksf*Ksr*b200;-q4-KsrA2*b3q3+q4+KsrA2*b3-Ksf*Ksr*b2Ksf*Ksr*b2-q30;Ksf*Ksr*b2-Ksf*Ksr*b2q2+KsfA2*b1-q2-KsfA2*b100;-Ksf*Ksr*b2Ksf*Ksr*b2-q2-KsfA2*b1q

21、1+q2+KsfA2*b10-q1;0-q300q30;000-q10q1;Q=zeros(4,10);zeros(6,4)sq;R=b1b2;b2b3;N=zeros(4,2);Ksr*b2Ksr*b3;-Ksr*b2-Ksr*b3;Ksf*b1-Ksf*b2;zeros(2,2);K=lqr(A,B,Q,R,N);%求主动悬架的时域响应loadroad_file.matw1=road_file(:,2);w2=w1;t=road_file(:,1);C=eye(10);D=zeros(10,2);Y,X=lsim(A-B*K,F,C,D,w1'w2',t);输出矩阵C为10

22、X10单位矩阵，即将所有状态均输出，因此系统输出X是完全相同的。5.2预瞄控制预瞄控制中，前后轮的路面输入不再相同，后轮的输入比前轮输入滞后轮路面输入w1与后轮路面输入w2满足：W?sew采用PaDa近似后，可转化为时域方程：Ksf*b2;-Ksf*b1Y与系统状态变量个预瞄时间，及前(t)A(t)Bw/t)则预瞄控制系统空间状态方程可写为:AFDXBFEUw-i0A0B同样可由黎卡提方程求出其反馈矩阵程序代码如下:%预瞄控制tao=0.14;a0=12/taoA2;a1=6/tao;a2=1;An=01;-a0-a1;Bn=-2*a1;6*a0;Dn=00;10;En=1;1;A0=AF*D

23、n;zeros(2,10)An;B0=B;zeros(2,2);C0=C,zeros(10,2);zeros(2,12);%C0=eye(12);D0=zeros(12,1);F0=F*En;Bn;Q0=Q,zeros(10,2);zeros(2,12);N0=N;zeros(2,2);K0=lqr(A0,B0,Q0,R,N0);Y0,X0=lsim(A0-B0*K0,F0,C0,D0,w1,t);5.3结果比较运行结果中的变量X、X0分别存储无预瞄和有预瞄的状态变量，预瞄控制主要改善后悬架的性能，因此选择后悬架的加速度、后悬架动行程及后轮动位移作为比较项目。程序代码如下：%绘制车身加速度曲线

24、，并计算均方根值%ba-最优控制车身加速度%ba0-预瞄控制车身加速度ba=diff(X(:,1)/0.005;ba0=diff(X0(:,1)/0.005;subplot(2,1,1)plot(t(1:end-1),ba)subplot(2,1,2)plot(t(1:end-1),ba0)BA=norm(ba,2)./(length(ba)40.5);BAO=norm(ba0,2)./(length(ba0)90.5);%绘制悬架动行程曲线，并计算其均方根值%swsr-最优控制悬架动行程%swsr0-预瞄控制后悬架动行程figure()swsr=X(:,5)-X(:,6);swsr0=X0(:,5)-X0(:,6);subplot(2,1,1)plot(t,swsr)subplot(2,1,2)plot(t,swsr0)SWS=norm(1000*swsr,2)./(length(swsr)90.5);SWSO=norm(1000*swsr0,2)./(length(swsr0)90.5);%绘制轮胎