行驶动力学建模仿真及主动悬架控制器设计DOC文档格式.docx

1、加权系数的选取决定了设计者对悬架性能的倾向，如对车身垂向振动加速度项选择较大的权值，则考虑更多的是提高车辆操纵稳定性。为方便起见，这里取车身垂向振动加速度的加权值q3=1。 将性能指标J的表达式（8）改写成矩阵形式，即： （9） 当车辆参数值和加权系数值确定后，最优控制反馈增益矩阵可有黎卡提（Riccati）方程求出，其形式如下： （10） 最优反馈控制增益矩阵，由车辆参数和加权系数决定。根据任意时刻的反馈状态变量X（t），就可得到t时刻作动器的最优控制力Ua，即： （11）3. 计算实例 这里，以某轿车的后悬架为例，给出一个完整的计算实例，包括车辆模型参数、仿真路面输入参数、控制器的设计参数

2、以及计算结果。此例中车辆以20m/s的速度在某典型路面上行驶，仿真时间T=50s。计算中输入的各参数及数值详见表2。表2 单轮车辆模型仿真输入参数值车辆模型参数符号单位数值簧载质量非簧载质量悬架刚度轮胎刚度悬架工作空间mbmwKsKtSWScKgN/mmm3204020000200000仿真路面输入参数路面不平度系数车速下截止频率G0Uf0m3/cyclem/sHz5.0x10-6200.1性能指标加权参数轮胎动位移悬架动行程车身加速度q1q2q38000051 仿真计算中以式（6）所示的滤波白噪声作为路面输入模型。白噪声的生成可直接调用MATLAB函数WGN（M，N，P）（此函数需要安装信号

3、处理工具箱Communications 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为最优控制反馈增益矩阵，S为

4、黎卡提方程的解，E为系统闭环特征根。 根据表2给出的仿真输入参数，本例中求得的最优反馈增益矩阵K为： K=（711.88 -1241.5 -19284 -2038.5 20864） 同时，还得到了黎卡提方程的解：在相同的仿真条件下，可将所设计的主动悬架系统与一个被动系统进行对比分析。在被动悬架系统中，取悬架刚度Ks=22000N/m，阻尼系数Cs=1000NS/m。除此之外，其他输入参数值均与主动悬架系统完全相同。4. MATLAB仿真过程1）生成路面输入模型代码如下：a=wgn（10001,1,20）;t=0:0.005:50;road_file（:,1）=t;,2）=a;save road

5、_file road_file2）参数输入load road_file.mat %载入路面数据模型Ks=22000;mb=320;Kt=200000;mw=40;f0=0.1;G0=0.000005;u=20; Kb=20000;Ks1=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,-

6、Cs/mw,Ks1/mw,（-Kt-Ks1）/mw,Kt/mw;B=1/mb,0; -1/mw,0; 0,0; 0,2*pi*sqrt（G0*u）;B1=0,0;C=1,0,0,0,0; 0,0,1,0,0; 0,0,0,1,0; 0,0,0,0,1;D=0,0; 0,0;K=711.88,-1241.5,-19284,-2038.5,20864;K1=0,0,0,0,0;3）用Simulink创建仿真框图 状态变量 输入与系统模块，如下图： 输出模块，如下图： 整体程序框图如下：4）结果分析可以直接通过双击scope查看输出的波形图，为更好比较主动悬架与被动悬架的差别，下面通过输出到work

7、space的状态变量编程绘图并计算均方根值。% 绘制车身加速度曲线，并计算均方根值%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.5）;BA1=norm（ba1,2）./（length（ba1）.0.5）;% 绘制悬架动行程曲线，并计算其

8、均方根值%sws-主动悬架动行程%sws1-被动悬架动行程figure（）sws=X.data（:,3）-X.data（:,4）;sws1=X1.data（:,3）-X1.data（:plot（X.time,sws）plot（X.time,sws1）SWS=norm（1000*sws,2）./（length（sws）.0.5）;SWS1=norm（1000*sws1,2）./（length（sws1）.0.5）;% 绘制轮胎动位移曲线，并计算其均方根值%dtd-主动悬架动位移%dtd1-被动悬架动位移dtd=X.data（:,4）-X.data（:,5）;dtd1=X1.data（:,4）-X

9、1.data（:plot（X.time,dtd）plot（X.time,dtd1）DTD=norm（1000*dtd,2）./（length（dtd）.0.5）;DTD1=norm（1000*dtd1,2）./（length（dtd1）.0.5）;结果如下：1车身加速度曲线2悬架动行程曲线3轮胎动位移 主动悬架与被动悬架性能指标均方根值比较性能指标主动悬架被动悬架车身加速度BA（m/s2）1.514.60悬架动行程SWS（mm） 34.4267.10轮胎动位移DTD（mm） 5.8734.575. 半车模型建模及仿真 半车模型建模悬架的半车模型有四个自由度，可选取前后轮垂直位移、前后悬架与车身

10、连接处垂直位移四个自由度，写出状态空间方程：以作为系统状态变量，状态空间方程如下： LQG控制器选取如下二次型性能指标：写成矩阵形式：由黎卡提方程求出反馈矩阵K，则，因此状态空间方程可写成：由此形式可方便求出其时域与频域响应。5.1随机线性最优控制 路面模型本例仿真车速为20m/s，轴距为2.8m，因此滞后时间为0.14s。若取仿真时间T为20s，采样时间间隔为0.005s，因而仿真点数4028。程序代码如下：% 生成路面模型a=wgn（4029,1,20）;20.14;r=a,t;save road_file r% 车身模型参数输入clearclcmb=690;I=1222;mwf=40;m

11、wr=45;Ksf=17000;Ksr=22000;Ktf=200000;Ktr=200000;a=1.3;b=1.5;% 仿真路面参数输入G0=5e-6;% 性能指标加权系数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

12、 -a3*Ksf 0 0; 0 0 0 0 0 0 -Ksf/mwf （Ksf-Ktf）/mwf 0 Ktf/mwf; 1 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0; 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控制器设

13、计，计算Q、R、N矩阵并求出反馈矩阵Kb1=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; -q4-Ksr2*b3 q3+q4+Ksr2*b3 -Ksf*Ksr*b2 Ksf*Ksr*b2 -q3 0; 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

14、 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;-Ksf*b1 -Ksf*b2;zeros（2,2）;K=lqr（A,B,Q,R,N）;% 求主动悬架的时域响应load road_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,w1w2,t）;输出矩阵C为10X10单位矩阵，即将所有状态均输出，因此系

15、统输出Y与系统状态变量X是完全相同的。5.2预瞄控制预瞄控制中，前后轮的路面输入不再相同，后轮的输入比前轮输入滞后一个预瞄时间，及前轮路面输入w1与后轮路面输入w2满足：采用PaDa近似后，可转化为时域方程：则预瞄控制系统空间状态方程可写为：同样可由黎卡提方程求出其反馈矩阵% 预瞄控制tao=0.14;a0=12/tao2;a1=6/tao;a2=1;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;C0=C,zeros（10,2）;zeros（2,12）;%C0=eye（12）;D0=ze

16、ros（12,1）;F0=F*En;Bn;Q0=Q,zeros（10,2）;N0=N;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（:plot（t（1:end-1）,ba0）BA0=norm（ba0,2）./（length（ba0）.0.5）;%s

17、wsr-最优控制悬架动行程%swsr0-预瞄控制后悬架动行程swsr=X（:,5）-X（:,6）;swsr0=X0（:,5）-X0（:plot（t,swsr）plot（t,swsr0）SWS=norm（1000*swsr,2）./（length（swsr）.0.5）;SWS0=norm（1000*swsr0,2）./（length（swsr0）.0.5）;%dtdr-最优控制动轮胎位移%dtdr0-预瞄控制后悬架轮胎动位移dtdr=X（:,6）-X（:,10）;dtdr0=X0（:,6）-X0（:plot（t,dtdr）plot（t,dtdr0）DTD=norm（1000*dtdr,2）./（length（dtdr）.0.5）;DTD0=norm（1000*dtdr0,2）./（length（dtdr0）.0.5）;、后悬架加速度曲线、后悬架动行程曲线、后轮动位移比较、性能指标均方根值比较后悬架加速度/（m/s2）后悬架动行程（mm）后轮动位移（mm）无预瞄1.525416.81656.0459有预瞄16.81646.0427可以看出，并没有得到理想的仿真结果，说明仿真过程有问题。但能力有限，未找出问题所在。