1 引言
油气悬架融合了液压与气压传动技术和机械悬架技术，能够有效缓解外界冲击并快速衰减振动，其非线性刚度和阻尼特性在车辆行驶平顺性中起重要作用[1]。
国内外对悬架系统的结构参数优化方面做了大量研究，韩寿松等[2]通过AMESim、RecurDyn和SIMULINK联合仿真，实现了半主动油气悬架的阻尼参数优化设计；田文朋等[3]针对连通式油气悬架，应用改进多目标遗传算法进行优化求解，大幅提高了车辆平顺性；姚琳等[4]通过联合仿真，采用遗传算法对悬架系统进行优化，使导弹发射筒筒口位移均方根值大幅降低；Gomes等[5]针对1/4车模型，应用PSO算法对悬架的刚度与阻尼特性进行了优化研究；Qin等[6]针对多缸油气悬架系统，研究了阻尼直径与初始充气体积对车辆平顺性的影响，为悬架系统的优化提供了参考；Kwon等[7]基于HPS(hydraulic power steering)的整车模型，通过多目标优化设计，得到了兼顾舒适性与稳定性互相矛盾的帕累托前沿；Han等[8]采用改进的帕累托人工鱼群算法，得到最优悬架参数，提高了配备机械弹性轮车辆的乘坐舒适性。相较于遗传算法，PSO算法更容易实现且需要调整的参数较少，通过改进可以极大改善算法的收敛速度与解的多样性。
本文中提出了一种基于改进PSO算法与GWO算法结合的混合算法，以车身加速度、悬架动挠度、轮胎动载荷为优化目标，通过SIMULINK建模仿真，研究优化前后油气悬架性能的变化对车身姿态与舒适性的影响。该算法同时具有PSO算法快速收敛的特性，又具有GWO算法精度高，鲁棒性强的特点。
2 九自由度整车动力学模型
本文中的研究对象为某三轴运输车辆，为了更全面分析车辆的车身姿态、操纵稳定性以及乘坐舒适性，建立九自由度整车模型，如图1所示。
图1 九自由度整车模型示意图
Fig.1 9-DOF vehicle model
油气悬架的刚度与阻尼特性需要用非线性表达式表示，影响油气悬架刚度和阻尼的主要因素是位移与速度。结合试验数据，分别以位移、速度为自变量，产生的刚度力与阻尼力为因变量，考虑静态时方程的准确性，去除多项式的常数项并在Matlab中使用最小二乘法拟合，得到的结果如图2、图3所示。三阶与四阶多项式的拟合曲线与试验结果基本吻合，说明使用多项式模型能够准确表示油气悬架的刚度阻尼特性。
图2 Matlab刚度力拟合曲线
Fig.2 Matlab stiffness force fitting curve
图3 Matlab阻尼力拟合曲线
Fig.3 Matlab damping force fitting curve
为了简化模型，本文中选用三阶非线性表达式，表达式如下：
(1)
以簧下质量为研究对象，根据牛顿第二定律，以车身静平衡位置为原点建立数学模型，当i=1,2时， j=1；当i=3,4时， j=4；当i=5,6时， j=7。
(2)
以车身为研究对象，建立如下数学模型：
(3)
(4)
(5)
(6)
式中：S0i(i=1,2,…,6)为对应轮胎相对位移；S、S2、S3、S4、S5、S6分别为左前、右前、左中、右中、左后、右后油气悬架活塞相对位移；mi(i=1,2，…，6)为簧下质量；m为车身质量；Jx、Jy分别为绕质心的侧倾和俯仰转动惯量，对应的侧倾角和俯仰角为θ、 β；kti、cti(i=1,2，…，6)为轮胎的刚度和阻尼；kj(j=1,2，…，8,9)、cj(j=1,2，…，8,9)为油气悬架非线性刚度与阻尼特性表达式系数；zi(i=1,2，…，6)为各簧下质量的垂向位移；z为车身的垂向位移；l1、l2、l3、d为悬架支点到质心的距离；qi(i=1,2，…，6)为路面激励。
3 Simulink仿真模型
根据式(1)～式(5)在Simulink中搭建仿真模型，由于非线性项的存在且整车9个自由度相互耦合，模型较为复杂，不易区分，因此使用Subsystem模块并进行封装，得到图4所示的整车动力学模型，仿真参数如表1所示。
表1 车辆主要参数
Table 1 Main parameters of vehicle
车辆参数数值前桥簧下质量m1、m2/kg280中桥簧下质量m3、m4/kg475后桥簧下质量m5、m6/kg475车身质量m/kg5 200轮胎等效刚度kf/(N·m-1)1.86×106轮胎等效阻尼/(N·m-1)1 500前桥与质心距离l1/m2.8中桥与质心距离l3/m1.3后桥与质心距离l2/m2.6左右轮距2×d/m2.3侧倾转动惯量Jx/(kg·m2)6 910俯仰转动惯量Jβ/(kg·m2)9 420
图4 Simulink仿真模型示意图
Fig.4 Simulink simulation model
路面激励是最主要的输入信号，只有尽可能接近实际路况才能更准确地模拟车辆的行驶特性，本文中采用计算方便的白噪声法创建路面模型。为了使整车的路面输入具有相关性，在前后轮路面激励间加入延时信号，对于左右轮，建立如下状态方程[9]：
(7)
式中： f0为空间截止频率，取0.011 m-1；v0为行驶速度；n为空间频率；n0为参考空间频率；Gq(n0)为路面不平度系数；ω(t)为随机白噪声。
设置车速为30 km/h，得到C级路面的激励仿真结果曲线如图5所示。可以看出，前后轮激励曲线波形基本相同，只在时间上存在滞后，左右轮波形变化趋势相同。
图5 C级路面激励仿真结果曲线
Fig.5 Grade C pavement simulation result
4 优化算法的改进
4.1 PSO算法
PSO算法也称鸟群觅食算法，能够有效解决复杂的优化问题。PSO算法的基本流程是初始化一群随机粒子，即随机解，利用种群每个粒子间的信息交换，从每个粒子位置中寻找个体最优位置pbest，并通过pbest更新种群最优位置gbest。假设粒子群的搜索空间为dim维，其中粒子i的位置为Xi=(xi1,xi2,…,xidim)，粒子i在空间中的速度为Vi=(vi1,vi2,…,vidim),粒子在空间中的速度更新公式为
vidim(k+1)=ωvidim(k)+c1r1 [pbestidim(k)-xidim(k)+]
c2r2[gbestidim(k)-xidim(k)]
(8)
粒子在空间中的位置更新公式为
xidim(k+1)=xidim(k)+vidim(k+1)
(9)
式中：k为迭代次数；vidim(k)、xidim(k)分别为粒子i当前在第dim维空间的速度与位置；ω为惯性因子；c1、c2为加速常数；r1、r2为区间[0,1]内的随机数。
4.2 GWO算法
GWO(灰狼算法)[10]是近年Mirjalili受狼群启发提出的一种优化算法。狼群拥有非常严格的等级，如图6所示，α狼位于金字塔顶端；β狼次之，δ狼和ω狼位于金字塔的第二层和最底端。GWO算法中，每个灰狼的位置代表空间中的一组可行解，因此可以得到最优解、亚优解、季优解和一般解4个等级的解，分别对应着狼群的4个社会等级。
图6 狼群等级分层示意图
Fig.6 Wolf pack hierarchy
设dim维空间目标位置Xp=(xp1,xp2,…,xpdim)，当前灰狼位置Xi=(xi1,xi2,…,xidim)，当前灰狼向目标移动后的位置为：
Xi(k+1)=Xp(k)-A·|C·Xp(k)-Xi1(k)|
(10)
式中：A为区间[-a,a]内的均匀随机数；a为常数，初始值为2，并随着迭代次数由2线性递减至0；C为区间[0,2]内的随机数。
每次迭代都保留3个头狼的最优位置并计算头狼与狼群的距离，以此更新整个狼群位置。
(11)
(12)
式中：X1、X2、X3分别表示α狼、β狼、δ狼的位置。
4.3 PSO-GWO混合算法
PSO算法结构简单，收敛速度快，但精度低，鲁棒性差。而GWO算法正好相反，鲁棒性较强，但收敛速度慢，融合2种算法的优点得到一种混合算法[11]。PSO-GWO算法主要以PSO算法为主，完成一次PSO算法后，将得到的种群最优位置gbest引入狼群中的α狼，此时α狼的位置为
Xα(k+1)=gbest(k)
(13)
为了避免粒子群陷入局部收敛，对PSO算法进行以下改进[12]：
1) 将GWO算法得到的头狼位置也作为最优解加入速度的计算。
2) 加速常数c1、c2分别代表粒子个体与群体对该粒子速度的影响，迭代前期，为了使种群多样性，PSO算法的粒子主要受个体信息的影响，迭代后期，为了快速收敛，粒子主要受社会群体影响，采用非线性变化的加速常数。
3) 惯性因子较大时有利于全局搜索，较小时有利于局部搜索，采用非线性的惯性因子变化策略。
修改后的速度公式为
vidim(k+1)=ωvidim(k)+c1r1[pbestidim(k)-xidim(k)]+
c2r2[gbestidim(k)-xidim(k)]+
c3r3[Xidimα(k)-xidim(k)]
(14)
加速常数更新公式为
(15)
惯性因子更新公式为
(16)
式中：ωmin、ωmax分别为惯性因子的最小值与最大值；cis、cif分别为加速常数迭代的初始值与最终值；k、 kmax为当前迭代次数与最大迭代次数。
以PSO算法为外循环，GWO算法为内循环，设PSO算法最大迭代次数为k2次，GWO算法最大迭代次数为k1次，则从开始到结束GWO算法共迭代k2×k1次，算法流程如图7所示。
图7 PSO-GWO算法流程框图
Fig.7 Flow chart of PSO-GWO algorithm
5 油气悬架参数优化
5.1 优化算法函数的设置
1) 目标函数
一般来说，车身加速度、悬架动挠度以及轮胎动载荷分别影响车辆的舒适性、车身姿态与稳定性，选取三者作为油气悬架系统的评价指标。通常情况，三项指标相互矛盾，悬架系统既要满足车辆的行驶稳定性，又要尽可能提高车辆舒适性，同时悬架动挠度也要在允许的范围内。由于优化目标有3个，将各目标函数进行无量纲处理并加权求和，故目标函数如下：
(17)
式中：rms表示各项的均方根；a0、a为优化前后车身加速度； fd0、 fd分别为优化前后悬架平均动挠度； Fw0、 Fw分别为优化前后轮胎平均动载荷； ρ1、 ρ2、 ρ3为子目标函数的加权系数。
2) 设计变量
在仿真模型中改变油气悬架刚度和阻尼参数，以式(17)为目标，经过大量试算，选取对目标函数影响较大的参数为设计变量。结合图2、图3所示的刚度阻尼曲线，考虑车辆行驶稳定性与舒适性的要求，各优化变量及取值范围如表2所示，表中ki的单位为N·m-1，ci的单位为N·s-1·m-1。
表2 设计变量取值范围
Table 2 Value range of design variables
k1k2k45.6×106～3.3×1074.6×106～1.6×1073.2×1010～2.3×1011k5k7 k81.5×1010～9.2×10104.5×109～1.3×10103.9×109～9.2×109c1c4c79.1×104～4.5×1059.1×104～4.5×1059.1×104～4.5×105
3) 约束条件
悬架动挠度均方根与轮胎动载荷均方根必须限制在合适的范围，一般来说，悬架动挠度应小于限位行程，轮胎动载荷应小于路面静载[13]。对于车身的侧倾角与俯仰角，一般限制车身侧倾角不超过6°，车身俯仰角不超过3°。考虑到安全性，建立约束的数学模型为：
(18)
式中：[fd]为限位行程；[G]为路面静载。
4) 罚函数设计
式(18)的约束条件均为不等式约束，可以采用内罚函数法设计罚函数[14]。相对于外罚函数在不可行区域施加惩罚，内罚函数法可直接在可行域边界“筑起高墙”，将目标函数拦在可行域内，可构造如下增广目标函数
(19)
式中： μk为惩罚因子,μk+1=0.1μk；B1(x)为罚函数，由所有约束条件构成。
5.2 优化结果
PSO-GWO算法设定初始种群大小为100，内循环最大迭代次数为20次，外循环最大迭代次数为200次，加速常数c1s、c1f、c2s、c2f、c3s、c3f分别为1.5、3.5、1.5、3.5、0.5、1.5，惯性因子ωmin、ωmax分别为0.4、0.9，通过大量试算，取ρ1=4.9， ρ2=3.6， ρ3=3.0。为评估改进算法的性能，将本文中算法与标准PSO算法进行比较，得到的迭代曲线如图8所示，可以看出，相比于PSO算法，PSO-GWO算法的收敛速度更快，得到最优解的精度更高，当陷入局部最优时能迅速跳出，优化前后的变量取值如表3所示。
图8 迭代收敛曲线
Fig.8 Iterative convergence curve
表3 刚度阻尼参数取值
Table 3 Optimum value of stiffness and damping
参数k1k2k4优化前6.12×1065.25×1068.54×1010优化后2.88×1077.12×1069.59×1010参数k5k7k8优化前7.33×10108.83×1096.6×109优化后3.17×10101.20×10105.52×109
续表(表3)
参数c1c4c7优化前2.50×1053.77×1051.28×105优化后3.44×1054.07×1051.94×105
在前中后三桥各选取一侧悬架进行图示说明，将优化后的各变量值代入模型，得到优化前后的车身加速度、悬架动挠度及轮胎动载荷仿真结果如图9～图11所示。
图9 优化前后车身加速度曲线
Fig.9 Optimize the acceleration curve of front and rear body
图10 优化前后悬架动挠度曲线
Fig.10 Optimize the dynamic deflection curve of front and rear suspension
图11 优化前后轮胎动载荷曲线
Fig.11 Optimize the dynamic load curve of front and rear tires
为了量化对比优化前后的结果，计算得到优化前、使用PSO算法优化后、使用PSO-GWO算法优化后的车身加速度均方根值(见图12)、各悬架动挠度均方根值(见图13)、各轮胎动载荷均方根值(见图14)，如表4所示。
表4 悬架性能指标均方根值
Table 4 Comparison of root mean square of suspension performance index
C级路面优化前PSOPSO-GWO加速度0.970.860.73动挠度13.2012.9613.6411.8911.4411.718.388.759.3013.4213.1113.7111.8111.7311.658.368.849.10动载荷9 1139 0468 7169 7859 4969 2169 9509 7898 8369 2698 7278 6819 8349 5469 16310 7619 6389 127
图12 车身加速度均方根值曲线
Fig.12 Root mean square of body acceleration
图13 悬架动挠度均方根值曲线
Fig.13 Root mean square of suspension dynamic deflection
图14 轮胎动载荷均方根值曲线
Fig.14 Root mean square of tire dynamic load
由图表可以看出，在种不同的路面等级下，优化后车身加速度与悬架动挠度性能均有所提升，且路况越差，提升效果越明显。优化后轮胎动载荷均方根增大，但仍然远小于路面静载，不足以影响轮胎对地面的抓附能力。研究结果表明，使用本文中方法对油气悬架参数进行优化后，在基本不影响车辆操纵稳定性的前提下，车辆的乘坐舒适性与车身姿态得到了明显改善。
6 结论
1) 选用三阶非线性表达式表示油气悬架的非线性刚度与阻尼，根据九自由度整车模型建立Simulink仿真模型。采用白噪声法生成路面输入信号，前后轮之间加入延时信号，左右轮采用具有相关性的路面激励状态方程，保证路面激励接近实际路况。
2) 以PSO算法与GWO算法为基础进行研究，GWO算法作为内循环， PSO算法为外循环，改进后的算法能有效解决多维非线性的问题，得到的结果精度高，收敛速度快，大大减少计算时间。
3) 在不同路面等级下使用改进PSO-GWO算法优化后车身加速度、悬架动刚度有明显提升，轮胎动载荷有所恶化，但不影响车辆的稳定性。研究结果可为进一步对油气悬架阻尼孔直径、蓄能器充气压力等参数的优化提供参考。
[1] 周长城.车辆悬架设计及理论[M].北京：北京大学出版社，2011.
Zhou C C.Vehicle suspension design and theory [M].Beijing:Peking University Press,2011.
[2] 韩寿松,刘相波,晁智强,等.基于粒子群算法的半主动油气悬架阻尼参数优化[J].兵器装备工程学报,2020,41(07):53-57.
Han S S,Liu X B,Chao Z Q,et al.Optimization of damping parameters of semi-active hydro pneumatic suspension based on particle swarm optimization[J].Journal of Ordnance and Equipment Engineering.2020,41(07):53-57.
[3] 田文朋,刘钊,王伟.连通式油气悬架整车建模与平顺性优化[J].噪声与振动控制,2017,37(06):81-85.
Tian W P.Liu Z,Wang W.Vehicle modeling and ride comfort optimization of connected hydro pneumatic suspension[J].Noise and Vibration Control,2017,37(06):81-85.
[4] 姚琳,马大为,王玺,等.冷弹射系统动力学响应与油气悬架优化设计[J].兵工学报,2015,36(07):1228-1236.
Yao L,Ma D W,Wang X,et al.Dynamic response of cold ejection system and optimization design of hydro pneumatic suspension[J].Acta Armamentarii,2015,36(07):1228-1236.
[5] Herbert M G.Multi-objective optimization of quarter car passive suspension design in the frequency domain based on PSO[J].Engineering Computations,2016,33(05):1422-1434.
[6] Qin B N,Zeng R Y,Li X M,et al.Design and performance analysis of the hydropneumatic suspension system for a novel road-rail Vehicle[J].Applied Sciences,2021,11(05):77-92.
[7] Kwon K H,Seo M S,Kim H,et al.Multi-objective optimisation of hydro-pneumatic suspension with gas-oil emulsion for heavy-duty vehicles[J].Vehicle System Dynamics,2020,58(07):1-20.
[8] Xu H,Zhao Y Q,Ye C,et al.Integrated optimization for mechanical elastic wheel and suspension based on an improved artificial fish swarm algorithm[J].Advances in Engineering Software,2019,137(11):102722.1-102722.8.
[9] 刘涛.车身姿态补偿半主动控制[J].重庆理工大学学报(自然科学) ，2020，34(07):63-74.
Liu T.Semi-active control for body attitude compensation[J].Journal of Chongqing University of Technology(Natural Science)，2020，34(07):63-74.
[10] 张子豪,靳其兵.基于社会等级淘汰机制的GWO-PSO算法[J].南京理工大学学报,2021,45(02):164-170.
Zhang Z H,Jin Q B.Gwo_PSO Algorithm based on social hierarchy elimination mechanism[J].Journal of Nanjing University of Technology,2021,45(02):164-170.
[11] Gohil B N,Patel D R.A hybrid GWO-PSO algorithm for load balancing in cloud computing environment[C]//2018 Second International Conference on Green Computing and Internet of Things(ICGCIoT).New York:IEEE,2018:185-191.
[12] 张伟,谢源海,王亚刚.基于收敛性分析的偏差粒子群算法及PID仿真应用[J].控制工程,2021,28(07):1466-1473.
Zhang W,Xie Y H,Wang Y G.Deviation particle swarm optimization algorithm based on convergence analysis and PID simulation application[J].Control Engineering,2021,28(07):1466-1473.
[13] 李伟平,柳超,窦现东,等.油气悬架的不确定性多目标优化[J].湖南大学学报(自然科学版),2014,41(10):27-34.
Li W P,Liu C,Dou X D,et al.Multi-objective optimization based on uncertainty in hydro pneumatic suspension[J].Journal of Hunan University(Natural Science Edition),2014,41(10):27-34.
[14] 钱思尧.基于内罚函数法的无线供电系统磁芯布局优化设计[D].南京：南京邮电大学,2020.
Qian S Y.Magnetic core layout optimization design of wireless power supply system based on internal penalty function method[D].Nanjing：Nanjing University of Posts and Telecommunications,2020.