线性多步量化状态系统算法及在溢流阀上应用

2021-08-10傅东金樊志华沈汉武李志华

杭州电子科技大学学报(自然科学版) 2021年4期

傅东金，樊志华，沈汉武，李志华

(杭州电子科技大学机械工程学院，浙江杭州 310018)

0 引言

常微分方程(Ordinary Differential Equation, ODE)系统广泛应用于机械、电子、控制、航空航天等工程领域，是建模仿真、数值计算、高性能计算等领域的重要研究方向之一。近年来，国家十分重视工业软件的开发，其中求解算法是工业软件的底层核心[1]。实际工程系统中的常微分方程常表现出刚性问题，即系统中同时存在快速与慢速变化量，刚性问题一直是传统数值积分方法难以解决的问题[2]。溢流阀的仿真计算中，由于阀门容腔压力巨变、油液压缩和多变量耦合等问题，导致溢流阀的动力学模型十分复杂，呈现高度的非线性和刚性，仿真求解十分困难[3]，而高效精确的求解溢流阀动力学问题是其设计的基础[4]。几十年来，尽管提出了不少求解刚性问题的方法，如隐式龙格库塔法、Gear法及其他一些特色方法[5-7]，但这些都是基于时间离散的方法，其算法性能的提升受到一定的限制。量化状态系统(Quantized State System, QSS)[8]是一种在计算过程中认为状态变量是不连续的算法，其逻辑与传统时间离散的积分方法不同，在QSS中可以通过量化函数将状态变量以“量子”大小进行分割，并行计算每个状态变量状态改变所需要的时间，这一过程称之为“跃迁”，再以最小跃迁时间为统一步长，采取个体异步、整体同步的时间推进思想，从而达到推进积分的目的，是一种并行式事件驱动型算法[9]。在求解一般ODE方程时，由于在求解过程中不进行迭代，故而其求解速度得到大幅提升，并且因为其本身的求解特性导致它天然具有强稳定、高精度等优点。

在QSS基础上，又提出了后向量化状态系统(Back Ward Quantized State System, BQSS)、中心量化状态系统(Centered Quantized State System, CQSS)等算法[10-12]，在不同工程领域取得了一定效果，但这些算法依然存在一些问题。首先这些算法没有考虑多变量耦合情况下多状态变量并行计算跃迁条件与单变量计算跃迁条件的差异，导致积分难以推进；其次这些算法在每一积分步均以量子大小进行跃迁，求解精度不高；最后已有的QSS系列算法均是使用单步法，没有充分利用已有的信息来提高算法的求解性能。目前，国内学者对QSS研究还较少，比如文献[13-15]对QSS方法进行了初步应用。本文以1个典型刚性算例和1个工程实例为研究对象，借鉴QSS方法和Adams多步法，提出一种线性多步量化状态系统算法(Linear Multi-Step Quantized State System, LMSQSS)，以提高刚性方程数值求解的精度和效率。

1 量化状态系统

一组由常微分方程表示的状态方程系统如下：

(1)

式中，x(t)为状态向量，x(t)∈R，u(t)为输入向量，QSS方法将式(1)量化为：

(2)

式中，q(t)为状态变量的量化向量，每个量化变量qj(t)通过量化函数得到：

(3)

式中，ΔQj为量子，qj(t-)为前一次的量化变量值。

量子由初始设定，量子的大小与算法最终求解误差呈正相关，文献[16]给提出了选择量子的方法。

图1 状态变量xj(t)和量化变量qj(t)的轨迹图

QSS的状态变量xj(t)和量化变量qj(t)的轨迹如图1所示，从图1可以看出，qj(t)轨迹为分段常数，当xj(t)的变化超过了1个“量子”时，引起qj(t)的变化，这个行为称之为“跃迁”，这一特性保证了算法的稳定性和全局误差界限，因为qj(t)与xj(t)的差值不会超过1个量子。

从图1还可以看出，任意轨迹的状态变量与相对应的量化变量共存在9种状态，而QSS未考虑到多变量耦合情况下单独状态变量跃迁条件与非耦合情况下的不同，只考虑图1中的1，2，5，6，8这5种状态[17]，本文算法完善后的跃迁条件表示为如下9种状态：

2 线性多步量化状态系统

本文所提的LMSQSS算法通过显式形式表达，避免了迭代求解，提高了求解效率；并利用求解过程中已知信息进行修正，得到更加精确的结果；LMSQSS算法本身包含隐式求解思想，将量化变量值作为系统状态变量预测值，进一步提升了算法的稳定性和精确性。

为了抑制仿真求解过程中的振荡，需要保证积分推进过程中xj的取值一直向着qj取值的方向靠拢，故qj的量化函数为：

(4)

(5)

根据不同跃迁速度确定跃迁时间为：

(6)

从而求出每推进一次仿真所需的时间为：

Δt=min{Δtj}

(7)

此时，状态变量的仿真推进时间为：

ti+1=ti+Δt

(8)

跃迁状态变量为：

(9)

非跃迁状态变量为：

(10)

式中，n表示多步法的阶数。

(11)

表1 各阶数对应的系数αk

LMSQSS算法具有不同的阶数，各阶系数通过Adams-Bashforth法得到，各阶数对应的系数如表1所示。

3 仿真算例实验

仿真实验在Windows7系统上的MATLAB2012a平台上进行，处理器为Intel i5-4260。分别采用LMSQSS算法、ODE45，ODE15s，ODE23s，QSS算法进行实验。ODE45是传统求解的首选算法，是显式变步长单步Runge_kutta法，当解的光滑性比较好时，其精度比较高；ODE15s和ODE23s都是针对刚性问题的优秀求解算法，ODE15s是隐式Gear变步长多步求解法，ODE23s是隐式Rosenbrock变步长单步求解法。为了研究不同阶数对LMSQSS算法性能的影响，实验选取2阶、3阶和4阶阶数进行仿真求解。对所有算法设置相同误差容限。

系统各状态变量相对误差计算公式如下：

(12)

式中，k=1,2,…,n，um[k]为各个算法求得的解，u[k]为DASSL求解器在很小误差设定(10-9)下求得的解，作为基准值用于计算相对误差[18]。

一个线性二阶ODE系统[19]如下：

(13)

图2 LMSQSS算法对算例的仿真结果

表2 不同算法对算例1仿真结果

从图2中可以看出，随着仿真时间的推进，变量x1和x2并未出现振荡或发散。

从表2可以看出，相比ODE45，2阶LMSQSS算法的效率提高了约435倍，精度提高22倍；相比传统刚性求解器中最好的ODE15s，效率提高了2倍，精度提高了69倍；相比ODE23s，效率提高了5倍，精度提高了82倍。在求解刚性二阶系统问题时，QSS算法的结果产生发散。

由表2可知，在相同量子情况下，随着阶数的增加，LMSQSS算法的求解精度明显增高，而CPU仿真时间变化并不大。

1—主阀左腔；2—主阀右腔弹簧；3—主阀右腔；4—过滤器；5—先导阀弹簧。图3 溢流阀结构简图

4 工程实例

4.1 工程实例模型

一个先导式溢流阀[20]的结构简图如图3所示。图3中，油液从主阀左腔流入，经主阀右腔和过滤器到达先导阀，在油液压力到达先导阀弹簧所设定压力值后，先导阀开启，油液流入，油液再经上部回油管路流回，致使主阀左右腔形成压力差，直到压力能够克服主阀弹簧弹力，使主阀芯移动，油液再经主阀口流出。采用功率键合图法对其进行建模[21-22]，得到数学模型的微分代数方程如下：

(16)

仿真初值：[x1,x2,x3,x4,x5,x6]=[V1,V2,P1,P2,Y1,Y2]=[0.1,6.37×10-3,0,0,0,0]，仿真时长为0.5 s。其中，V1为泵至主阀左腔管道容腔中用来补偿油液压缩量及管道受压变形量的液压油总容积，单位为cm3；V2为主阀右腔中用来补偿油液压缩量的液压油总容积，单位为cm3；P1和P2分别为主阀芯和先导阀阀芯的动量，单位为N·s/cm2；Y1和Y2分别为主阀芯和先导阀阀芯的位移，单位为cm。其余参数在表3中进行说明。

表3 溢流阀参量值

4.2 工程实例仿真结果分析

采用LMSQSS算法对溢流阀数学模型进行求解，得到图4阀芯位移轨迹图和图5入口压力图。从图4可以直观看到溢流阀的阀芯位移变化情况，其中主阀芯位移轨迹很平稳，并未因为系统耦合而发生振荡；先导阀轨迹出现短暂波动后(0.02 s)，趋于平稳。从图5中可以看出，阀的入口压力在油液进入后开始急剧上升，其动态超调量在22%左右，上升时间在0.002 0 s左右，峰值出现的时间在0.003 5 s左右。在先导阀开启后，压力迅速下降，短暂波动后，在0.006 0 s后趋于稳定。由此可见，LMSQSS算法具有较好的稳定性和收敛性，能够胜任溢流阀动力学的实际求解。