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

1、目录1. 计算机仿真系统模型的建立.- 1 -2. LOG 控制器设计.-2 -3. 计算实例.-3 -4. MATLAB 仿真过程.-4 -5. 半车模型建模及仿真.-8 -5.1 随机线性最优控制.-9 -5.2 预瞄控制 .-11-5.3 结果比较.-12 -Word 文档以单轮车辆模型为例， 介绍行驶动力学计算机建模、 仿真分析以及利用线性二次最优控制理论进行主动悬架 LQG 控制器设计过程。1. 计算机仿真系统模型的建立根据图 7 所示的主动悬架单轮车辆模型，运用牛顿运动定律，建立系统的运动方程，即：&U aK s ( xbxw)mb xb（ 4）&U aK s ( xbxw ) K

2、t (xwX g )mw xw（ 5）这里，采用一个滤波白噪声作为路面输入模型，即：&2f0 xg (t )2 G0uw(t )xg (t )（ 6）式中，xg 为路面垂向位移 （ m）；G0 为路面不平度系数 （ m3/cycle ）；u 为车辆前进速度 （ m/s）； w 为数字期望为零的高斯白噪声； f0 为下截止频率（ Hz）。图 7单轮车辆模型结合式（ 4）、式（ 5）和式（ 6），将系统运动方程和路面输入方程写成矩阵形式，即得出系统的空间状态方程：&AXBU FW（ 7）X&xbxwxg )TW= （ w（t），为高斯白噪声输入矩式中， X( xbxw，为系统状态矢量；阵； U=（

3、 Ua （t），为输入控制矩阵；Word 文档00K sK s01mbmbmb000K sKt K sK a10mwmwmwBmwA100000F0010000000002 f0；0；2 G0u2. LOG 控制器设计车辆悬架设计中的主要指标包括：代表轮胎接地性的轮胎动载荷； 代表轮胎舒适性的车身垂向振动加速度；影响车身姿态且与轮胎布置有关的悬架动行程。因此，LQG控制器设计中的性能指标J 即为轮胎动位移、 悬架动行程和车身垂向振动加速度的加权平方和在时域 T 的积分值，其表达式为：T&212q2 xb (t) xw (t )2(t) dtJ lim q1 xw (t) xg (t)q3 xb

4、TT 0（ 8）式中， q1、q2 和 q3 分别为轮胎动位移、悬架动行程和车身垂向振动加速度的加权系数。加权系数的选取决定了设计者对悬架性能的倾向，如对车身垂向振动加速度项选择较大的权值，则考虑更多的是提高车辆操纵稳定性。为方便起见， 这里取车身垂向振动加速度的加权值 q3=1。将性能指标J 的表达式（ 8）改写成矩阵形式，即：J lim1 T (XTQXUTRU2 X T NU )dtTT 0（ 9）式中，000000000000K s2q2K s200Qq2mb20mb200q2K s2K s2q1NK amb2q1 q21K smb2R000q1q1mb20；当车辆参数值和加权系数值确

5、定后，最优控制反馈增益矩阵可有黎卡提（Riccati ）方程求出，其形式如下：PA ATP (PB N)R 1(BTP NT) Q 0（ 10）最优反馈控制增益矩阵KBT PN T ，由车辆参数和加权系数决定。根据任意时刻的Word 文档反 馈 状 态 变 量 X （ t ）， 就 可 得 到 t时 刻 作 动 器 的 最 优 控 制 力 Ua ， 即 ：U a (t )KX (t)（ 11）3. 计算实例这里， 以某轿车的后悬架为例，给出一个完整的计算实例，包括车辆模型参数、仿真路面输入参数、 控制器的设计参数以及计算结果。此例中车辆以20m/s 的速度在某典型路面上行驶，仿真时间T=50s

6、。计算中输入的各参数及数值详见表2。表 2 单轮车辆模型仿真输入参数值车辆模型参数符号单位数值簧载质量mbKg320非簧载质量mwKg40悬架刚度KsN/m20000轮胎刚度KtN/m200000悬架工作空间SWScmm100仿真路面输入参数符号单位数值路面不平度系数G0m3/cycle5.0x10-6车速Um/s20下截止频率f0Hz0.1性能指标加权参数符号单位数值轮胎动位移q180000悬架动行程q25车身加速度q31仿真计算中以式（6）所示的滤波白噪声作为路面输入模型。白噪声的生成可直接调用MATLAB函数 WGN（ M ，N ，P）（此函数需要安装信号处理工具箱Communicati

7、ons toolbox ），其中 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

8、为最优控制反馈增益矩阵， S 为黎卡提方程的解，E 为系统闭环特征根。根据表 2 给出的仿真输入参数，本例中求得的最优反馈增益矩阵K 为：K= （ 711.88-1241.5-19284-2038.520864）同时，还得到了黎卡提方程的解：2.45590.02892.47458.66077.3090.02890.48860.02987.52627.2364S2.47450.02984.97448.67545.10338.66077.52628.67542710.12700.47.3097.23645.10332700.42693.7在相同的仿真条件下，可将所设计的主动悬架系统与一个被动系统进

9、行对比分析。在被动悬架系统中，取悬架刚度Ks=22000N/m ，阻尼系数Cs=1000NS/m 。除此之外，其他输入参数Word 文档值均与主动悬架系统完全相同。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;K

10、s1=22000;Cs=1000;%输入仿真有关参数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

11、;0,0;0,0;0,2*pi*sqrt(G0*u);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;Word 文档0,0;0,0;0,0;K=711.88,-1241.5,-19284,-2038.5,20864;K1=0,0,0,0,0;3）用 Simulink 创建仿真框图?状态变量 X xb , xw , xb , xw , xg 输入与系统模块，如下图：输出模块，如下图：车身加速度悬架动行程轮胎动位移路面输入d?( xb )dtxbxwxwxgxgWord 文档整体程序框图如下：4）结果分析可以直接通过双击

12、scope 查看输出的波形图， 为更好比较主动悬架与被动悬架的差别， 下面通过输出到 workspace 的状态变量编程绘图并计算均方根值。代码如下：% 绘制车身加速度曲线，并计算均方根值%ba- 主动悬架车身加速度%ba1- 被动悬架车身加速度ba=diff(X.data(:,1)./diff(X.time);ba1=diff(X1.data(:,1)./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).0

13、.5);BA1=norm(ba1,2)./(length(ba1).0.5);% 绘制悬架动行程曲线，并计算其均方根值%sws- 主动悬架动行程%sws1- 被动悬架动行程figure()sws=X.data(:,3)-X.data(:,4);Word 文档sws1=X1.data(:,3)-X1.data(:,4);subplot(2,1,1)plot(X.time,sws)subplot(2,1,2)plot(X.time,sws1)SWS=norm(1000*sws,2)./(length(sws).0.5);SWS1=norm(1000*sws1,2)./(length(sws1).0

14、.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(1000*dtd,2)./(length(dtd).0.5);DTD1=norm(1000*dtd1,2)./(length(dtd1).0.5);结果如下：车身加速度曲线主动悬架车身加速度10)2s/5m(度0速加身-5车-1

15、0012345678910时 间(s)被动悬架车身加速度20)2s/10m(度0速加身 -10车-20012345678910时 间 (s)悬架动行程曲线Word 文档)m(主动悬架动行程0.10.05移0位-0.05-0.1)m(012345678910时 间(s)被动悬架动 行 程0.60.4移 0.2位0-0.2012345678910时 间 (s)轮胎动位移主动悬架轮胎动位移0.040.02)m(移0位-0.02-0.04012345678910时 间 (s)被动悬架轮胎动位移0.10.05)m(移0位-0.05-0.1012345678910时 间 (s)主动悬架与被动悬架性能指标均

16、方根值比较性能指标主动悬架被动悬架车身加速度 BA(m/s 2)1.514.60悬架动行程 SWS(mm)34.4267.10轮胎动位移 DTD(mm)5.8734.575. 半车模型建模及仿真半车模型建模Word 文档悬架的半车模型有四个自由度， 可选取前后轮垂直位移、 前后悬架与车身连接处垂直位移四个自由度，写出状态空间方程：?以 X z4 , z3 , z2 , z1, z4 , z3 , z2 , z1, x10, x20 作为系统状态变量，状态空间方程如下：?X (t )AX (t )BU (t )FW (t)UU af ,为控制输入矩阵，即前 后悬架作动器力；U arW w1 ,

17、为路面输入模型的高斯 白噪声输入。 w2LQG 控制器选取如下二次型性能指标：J lim 1T? 2? 2 q1( z1 z01 ) 2q2 (z2 z1) 21 z2q3 ( z3 z02 )22 z4 dtTT 0q1 ,q3 , 前后悬架轮胎动位移加权系数；q2 ,q4, 前后悬架动行程加权系数；1 ,2 ,车身加速度加权系数。写成矩阵形式：TJlim 1( X T QX U T RU 2X T NU ) dtTT 0由黎卡提方程求出反馈矩阵K，则 UKX ，因此状态空间方程可写成：?X (t ) ( ABK ) X (t)FW (t )( ABK ) - 状态矩阵F 输入矩阵由此形式可

18、方便求出其时域与频域响应。5.1 随机线性最优控制路面模型本例仿真车速为20m/s，轴距为2.8m，因此滞后时间为0.14s。若取仿真时间T 为 20s，采样时间间隔为0.005s，因而仿真点数4028。程序代码如下：% 生成路面模型a=wgn(4029,1,20);t=0:0.005:20.14;r=a,t;saveroad_filerWord 文档% 车身模型参数输入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;

19、f0=0.1;% 性能指标加权系数q1=80000;q2=100;q3=80000;q4=100;% 半车悬架模型建模，计算A、 B、 F矩阵a1=1/mb+b2/I;a2=1/mb-a*b/I;a3=1/mb+a2/I;A=0 0 0 0 a1*Ksr -a1*Ksr a2*Ksf -a2*Ksf 0 0;0 0 0 0 -Ksr/mwr (Ksr-Ktr)/mwr 0 0 Ktr/mwr 0;0 0 0 0 a2*Ksr -a2*Ksr a3*Ksf -a3*Ksf 0 0;0 0 0 0 0 0 -Ksf/mwf (Ksf-Ktf)/mwf 0 Ktf/mwf;1000000000;01

20、00000000;0010000000;0001000000;0 0 0 0 0 0 0 0 -2*pi*f0 0;0 0 0 0 0 0 0 0 0 -2*pi*f0;B=a2 a1;0 -1/mwr;a3 a2;-1/mwf 0;zeros(6,2);F=zeros(8,2);2*pi*sqrt(G0*u) 0;0 2*pi*sqrt(G0*u);% LQG控制器设计，计算 Q、 R、 N矩阵并求出反馈矩阵 K b1=a32+a22;b2=a2*a3+a1*a2;b3=a22+a12;sq=q4+Ksr2*b3 -q4-Ksr2*b3 Ksf*Ksr*b2 -Ksf*Ksr*b2 0 0;

21、 -q4-Ksr2*b3 q3+q4+Ksr2*b3 -Ksf*Ksr*b2 Ksf*Ksr*b2 -q3 0;Word 文档Ksf*Ksr*b2 -Ksf*Ksr*b2 q2+Ksf2*b1 -q2-Ksf2*b1 0 0;-Ksf*Ksr*b2 Ksf*Ksr*b2 -q2-Ksf2*b1 q1+q2+Ksf2*b1 0 -q1;0 -q3 0 0 q3 0;0 0 0 -q1 0 q1;Q=zeros(4,10);zeros(6,4) sq;R=b1 b2;b2 b3;N=zeros(4,2);Ksr*b2 Ksr*b3;-Ksr*b2 -Ksr*b3;Ksf*b1 Ksf*b2;-Ks

22、f*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 为 10X10 单位矩阵，即将所有状态均输出，因此系统输出Y 与系统状态变量X 是完全相同的。5.2 预瞄控制预瞄控制中， 前后轮的路面输入不再相同，后轮的输入比前轮输入滞后一个预瞄时间，及前轮路面输入w1 与后轮路面输入w2 满足：w2se采用 PaDa

23、近似后，可转化为时域方程：?(t )A(t )B w1(t )则预瞄控制系统空间状态方程可写为：?AFD XFEXB?0AUw10B同样可由黎卡提方程求出其反馈矩阵K p程序代码如下：% 预瞄控制tao=0.14;a0=12/tao2;a1=6/tao;a2=1;Word 文档An=0 1;-a0 -a1;Bn=-2*a1;6*a0;Dn=0 0;1 0;En=1;1;A0=A F*Dn;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,z

24、eros(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 分别存储无预瞄和有预瞄的状态变量，预瞄控制主要改善后悬架的性能，因此选择后悬架的加速度、后悬架动行程及后轮动位移作为比较项目。程序代码如下：% 绘制车身加速度曲线，并计算均方根值%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).0.5);BA0=norm(ba0,2)./(length(ba0).0.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(100