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

目录1 .建立计算机仿真系统模型22. LOG控制器设计33 .算例44. MATLAB仿真过程55 .半电机模型的建模与仿真95.1随机线性最佳控制105.2预读取控制部125.3结果比较13I以单轮车辆模型为例，介绍了行驶动力学计算机建模、仿真分析以及利用线性二次最优控制理论的主动悬架LQG控制器的设计过程。1 .建立计算机仿真系统模型基于图7所示主动悬架单轮车辆模型，利用牛顿运动定律建立系统的运动方程(4)(5)在此，采用滤波白噪声作为路面输入模型(6)式中，xg是路面垂直位移(m) G0是路面凹凸系数(m3/cycle )，u是车辆的前进速度(m/s) w是希望数字为零的高斯白噪声，f0是截止频率(Hz )。图7单轮车辆模型组合式(4)、式(5)和式(6)，将系统的运动方程式和路面输入方程式作为矩阵形式，得到系统的空间状态方程式(7)其中W=(w(t ) )为系统状态向量或U=(Ua(t ) )为高斯白噪声输入矩阵是输入控制矩阵2. LOG控制器设计车辆悬架设计中的主要指标是轮胎接地性的轮胎动载荷轮胎舒适性的车身垂直振动加速度影响车身姿势，与轮胎配置相关的悬架行程。 因此，LQG控制器设计中的性能指标j是轮胎的动位移、悬挂行程、车身的垂直振动加速度的加权平方和在时域t中的积分值，其式子如下(8)式中，q1、q2、q3分别是轮胎动位移、悬挂行程、车体的垂直振动加速度的加权系数。 权重系数的选择决定了设计者对悬架性能的倾向，车体的垂直振动加速度项选择较大的权重，可以提高车辆的操纵稳定性。 为了方便，这里取车体的垂直振动加速度的加权值q3=1。将性能指标j的公式(8)改写为矩阵形式后，如下所示(9)式中在确定车辆参数值和加权系数值之后，使用反馈增益矩阵来获得最佳控制反馈增益矩阵。 其形式如下(10 )最佳反馈控制增益矩阵由车辆参数和权重系数确定。 通过任意时刻的反馈状态变量X(t )，得到t时刻致动器的最佳控制力Ua、即(11 ) .3 .算例在此，以某轿车后悬架为例，给出了包括车辆模型参数、仿真路面输入参数、控制器设计参数、计算结果在内的完整计算实例。 在该例子中，车辆在作为20m/s的速度的典型路面上行驶，模拟时间T=50s。 计算中输入的各参数和数值详见表2。表2单轮车辆模型仿真输入参数值车辆模型参数符号单位数值弹簧载荷质量非弹簧载荷质量悬架刚度轮胎刚性悬挂作业空间美国职棒大联盟宝马KsKtSWSc公斤公斤N/mN/mmm3204020000仿真路面输入参数符号单位数值路面不规则度系数车速下截止频率G0uf0m3/cyclem/sHz5.0x10-6200.1绩效指标加权参数符号单位数值轮胎运动位移悬挂行程车身加速度q1.q1q2q3.q38000051在仿真计算中，将式(6)所示的滤波白噪声作为路面输入模型。 白噪声的生成可以直接调用MATLAB函数WGN(M，n，p ) (该函数需要安装信号处理工具箱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是递归方程的解，e是系统闭环特征根。根据表2所示的仿真输入参数，本例中求出的最佳反馈增益矩阵k如下所示k=(711.88-1241.5-19284-2038.520864 )同时，得到了利卡蒂方程的解在相同的仿真条件下，可以对设计的主动悬架系统和被动系统进行比较分析。 在被动悬架系统中，悬架刚性Ks=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；save road_file road_file2 )参数输入代码如下所示load road_file.mat %加载路面数据模型Ks=22000； mb=320； Kt=； mw=40； f0=0.1； G0=0. 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；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；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 ) ；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 ；利用Simulink制作仿真框图l状态变量l输入l输入和系统模块，如下图所示l输出模块如下图所示整个l的程序框图如下所示4 )结果分析直接双击范围可以查看输出的波形图，为了比较主动悬架和被动悬架之间的差异，可以使用下一个要输出到工作区的状态变量对绘图进行编程，然后计算均方根值。代码如下所示绘制%车体加速度曲线，计算均方根值-主动悬架车体加速度1-被动悬架车身加速度ba=diff (x.data (:1 )./diff (x.time )ba1=diff (x1.data (:1 ) )./diff (x1.time )subplot (2，1，1 )plot (x.time (1:结束-1)，ba )subplot (2，1，2 )plot (x1.time (1:结束-1)，ba1 )BA=norm(ba，2)./(length(ba).0.5)ba1=norm (ba1，2 )/(length (ba1 ) 0.5 )绘制%悬挂笔划曲线并计算其均方根值%sws-主动悬架行程%sws1-被动悬架行程figure ()SWS=x.data (:3 )-x.data (:4 )SW S1=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 * SWS 1，2 )/(SWS1) 0. 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 * dtd 1，2 )/(dtd1) 0. 5；结果如下：1车身加速度曲线双悬架行程曲线3轮胎的动态位移l主动悬架与被动悬架性能指标均方根值的比较绩效指标主动悬架无源悬架车身加速度BA(m/s2)1.514.60悬架行程SWS(mm )34.4267.10轮胎动态位移DTD(mm )5.8734.575 .半电机模型的建模与仿真l半马达型号悬架的半挂车模型有4个自由度，可以选择前后轮的垂直位移、前后悬架与车身的连接部的垂直位移这4个自由度来写状态空间方程式作为系统状态变量，状态空间方程式如下lqg控制器选择以下辅助性能指标写矩阵形式：从反馈方程式求出反馈矩阵k后，状态空间方程式可以写成如下以这种形式，可以容易地获得时域和频域的响应。5.1随机线性最优控制l路面模型本例的模拟车速为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；save road_file r%车体模型参数输入clearclcmb=690；I=1222；mwf=40；mwr=45；Ksf=17000；Ksr=22000；Ktf=；Ktr=；a=1.3；b=1.5；%模拟路面参数输入G0=5e-6；u=20；f0=0.1；%绩效指标加权因子q1=80000；q2=100；q3=80000q4=100；建模%半挂车悬架模型，计算a、b、f矩阵a1=1/mb b2/I；a2=1/mb-a*b/I；a3=1/mb a2/I；a= 00000 a1* ksr-a1* ksr2* ksf-a2* ksf 0；000-ksr/mwr (ksr-ktr )/MW r0ktr/MW r 0；0000 a2* ksr-a2* ksr a3* ksf-a3* ksf 0；000000-ksf/mwf (ksf-KTF )/mwf0KTF/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 ) 2* pi * sqrt (g0* u ) 0； 0 2*pi*sqrt(G0*u)；计算% LQG控制器设计、q、r、n矩阵，求出反馈矩阵kb1=a32 a22；b2=a2*a3 a1*a2；b3=a22 a12；sq= q4ksr 2* b3- q4- ksr 2* b3ksf * ksr * b2- ksf * ksr * b 20；- q4- ksr 2* b3q4k Sr 2