基于粒子滤波算法的涡旋压缩机性能预测研究

李 超，刘忠良，金银霞，侯军军，魏 宁，尹贺龙

（兰州理工大学 石油化工学院，兰州 730050）

0 引言

涡旋压缩机在市场中的广泛应用使得其前景广阔，因为工作效率高、运行噪声小、可靠性高，在制冷空调、真空泵、余热回收等行业广受好评［1-2］。但涡旋压缩机性能状态的优良与否直接关系到工作的可靠性，因此对其进行性能预测和状态监测，有助于了解其运转状况，可以有效预防运行不良等问题。在涡旋压缩机性能分析和预测方面，国内外部分学者进行了相关研究工作。CHEN等［3］通过建立质量守恒和能量守恒方程，分析了几何参数对涡旋压缩机的泄漏量以及热力学性能的影响；DINIZ等［4］针对涡旋压缩机传热过程建立了集中参数的温度分布预测模型，该模型能够准确预测标定范围内运行工况下压缩机的温度分布；彭斌等［5］通过设计变截面涡旋压缩机数学模型，详细分析了压缩机工作过程中实际气体的温度、压力及质量随曲轴转角变化的规律；李超等［6］结合几何模型分析了涡旋压缩机热力过程，建立了热力计算数学方程，并分析了变转速条件下压缩机的热力参数变化，通过实验验证了数学模型的准确性；彭斌等［7］在考虑压缩腔内介质泄漏和换热的条件下建立了无油涡旋压缩机的热力学模型，分析了气体在整个工作过程中容积、温度、压力以及质量随主轴转角的变化关系，并进一步分析了不同传热和泄漏对热力性能的影响规律；李超等［8］基于质量守恒原理和热传导原理，建立了实际复杂工况下多变压缩过程指数的理论模型，分析计算得到压缩机工作过程中气体质量变化过程指数和传热过程指数的变化规律。

现有的涡旋压缩机性能分析与预测都是通过求解数学方程来进行计算研究，随着智能算法的兴起，为压缩机性能分析提供了更多研究方法。本文提出了一种基于系统状态空间模型智能算法的预测方法，即利用粒子滤波算法对涡旋压缩机进行性能分析和状态参数预测［9］。1999年，CARPENTER等［10］首次提出了粒子滤波算法，在近20年的发展中，粒子滤波已广泛应用于状态预测、性能估计、故障诊断、目标跟踪等各个领域。焦自权等［11］通过更新建议分布的方法改进了粒子滤波算法并建立了锂电池的剩余使用寿命预测模型；曹洁等［12］提出基于GPU平台的粒子群优化粒子滤波的并行算法，对风力机组变桨距系统故障进行诊断，验证了优化算法的有效性；李位星等［13］提出了一种粒子滤波目标跟踪算法，用该算法对复杂场景下目标跟踪算法存在的目标丢失和漂移问题行进实验分析，验证了该算法的有效性。

本文根据能量守恒定律和气体状态方程，综合考虑涡旋压缩机实际运行过程中泄漏和传热等影响因素，推导计算了涡旋压缩机的温度、压力和质量的数学模型及变化规律，然后建立了涡旋压缩机温度和压力变化的状态方程和观测方程，搭建标准的粒子滤波算法进行计算，最后通过结果对比分析，验证了本文提出的粒子滤波算法的实用性和有效性。

1 热力过程计算

涡旋压缩机运行过程中实际工况比较复杂，由于受到压缩机系统本身的制造误差以及周围环境的影响，使得在建立数学模型时，需要考虑的影响因素较多，在理论计算过程中，为了使数学模型易于求解，同时能够使得求解过程更加接近于压缩机实际工况，作如下假设：

（1）涡旋压缩机主轴旋转至任意角度下，介质状态是均匀的，即在工作腔内任意位置处的状态参数一致；

（2）假设工作腔内气体总的能量为焓值；

（3）主轴转速保持不变。

1.1 理论推导

取涡旋压缩机中任意单一腔体为控制体，结合能量守恒定律，得到dt时间内控制体中工质总能量变化如式（1）所示：

h ——工质比焓，kJ/kg；

v ——比体积，m3/kg。

控制体的质量平衡方程：

由于在实际工作过程中，进出压缩腔体的工质，重力势能及动能很小，可以在理论计算中忽略不计，则总能量Ecv可用内能Ucv表示,如式（3）所示：

根据热力学关系式h=u+Pv，式（3）可表示为：

式中 Cv——恒容比热，kJ/（kg·K）；

V ——瞬时容积，m3。

控制体热力学第一定律可表示为：

假定流体从控制容积中流出时为瞬时稳定流，即h=hout，涡旋压缩机工作腔内工质温度随主轴转角的变化关系如下式：

质量平衡方程可进一步表示为：

气体状态方程如下式所示：

式中 R ——气体常数，kJ/（kg·K）。

利用全导数关系得到涡旋压缩机工作腔内工质压力随主轴转角的变化关系如下式：

1.2 质量流量

气体质量流量为：

ψ——流量系数；

A ——工质流动面积，mm2；

Ph——高压侧压力，kPa；

ρh——高压侧密度，kg/m3；

k ——绝热指数。

质量流量的临界压力比为：

1.3 泄漏模型

涡旋压缩机泄漏形式主要包括径向泄漏以及切向泄漏。动（静）涡旋盘齿顶与静（动）涡旋盘底部之间通过轴向气体力作用产生轴向间隙形成径向泄漏，如图所示1（a）所示；动、静涡旋盘啮合处因为径向气体力作用产生的径向间隙形成切向泄漏，如图1（b）所示。

图1 泄漏模型Fig.1 Leakage model

径向泄漏线的计算式为：

泄漏面积为：

切向泄漏面积可表示为涡旋齿高与径向间隙的乘积：

式中 δr——轴向间隙，mm；

φk，φk+1——动静涡旋盘啮合角，rad；

δf——径向间隙，mm。

1.4 传热模型

涡旋压缩机工作时，工质进入吸气腔后，工质在压缩腔内发生比较复杂的热传递过程，因为工质主要是通过涡旋齿的壁面发生热传递，由于实际传热过程复杂多变，理论计算时，简化传热模型，忽略热传导和热辐射的影响因素，只考虑热对流现象［3］。气体在涡旋腔体中的传热量为：

式中 hc——换热系数，kJ/（m2·K）。

涡旋压缩机的动盘与静盘的几何模型形似于螺旋板式换热器，所以传热模型可以选用矩形截面的螺旋板式换热器所采用的平均对流换热模型，那么涡旋压缩机的传热系数可以表示为：

其中

式中 Re ——雷诺数；

Pr ——普朗特数；

Def——当量直径，mm；

Raver——平均半径，mm。

2 粒子滤波算法

粒子滤波是一种基于蒙特卡洛采样思想的贝叶斯滤波算法，该算法通过在某系统中寻找一组带有权值的随机粒子样本，然后用该样本均值代替复杂的积分运算，从而获得近似系统随机变量的概率密度，得到系统状态的最小方差估计。

首先通过数学方程建立系统状态空间模型，如下式所示：

式中 xk，xk-1——系统在某一时刻的状态值；

f（*）——系统的非线性状态函数；

Qk，Rk—— 两组互不相干的噪声序列，其中Qk表示系统状态的过程噪声，Rk表示观测过程的测量噪声；

h（*）——系统的非线性观测函数。

根据一系列已有的测量信息y1:k={y1，y2，…，yk}递推计算当前k时刻的状态变量xk的状态估计值，从而获得系统的后验概率分布p（xk|y1:k），该过程主要包括预测和更新2个步骤，其中状态预测方程为：

状态更新方程为：

其中

由于p（xk|y1:k）是多变量、非标准的，所以无法直接获得后验概率p（xk|y1:k）并从目标分布中直接采样。因此需要引入一个已知且容易采样的重要性密度函数q（xk|y1:k），引入重要性采样方法可以提高采样的效率，这样就避免了直接从后验概率密度中抽取样本的困难。

粒子滤波算法就是通过蒙特卡洛模拟和序列重要性采样来近似表示后验概率密度。当粒子数足够多时，后验概率密度p（xk|y1:k）可由下式近似表示：

式中——第i个粒子对应的归一化权重；

δ（*）——狄克拉函数。

随着时间递推，每增加一次采样，p（xk|y1:k）都需要重新计算，计算量不断增加，为了解决这一问题，提出序列重要性采样。该采样方法在k+1时刻采样时，不改动过去时刻的状态序列样本集，而是采用递归的方式计算重要性权值。在系统观测条件独立的情况下，假设系统状态符合马尔科夫链过程，则重要性密度函数可改写为：

则可以得到粒子未归一化的重要性权值计算式为：

引入权重后，所有粒子不再是直接相加求平均，而是一种加权和的形式，不同粒子都有其相应的权重，如果粒子权重大，说明信任该粒子比较多。最后粒子的状态输出xk的状态估计值为：

由于序列重要性采样的过程中，存在粒子退化问题，为了解决粒子退化问题，提出重采样方法，如图2示出重采样的实现过程，重采样能够减小无效样本（权重小的粒子），增加有效样本（权重大的粒子）。但是引进重采样会增加算法时间，所以并非所有情况下的粒子都要进行重采样，通常利用特定的粒子退化程度来判断，这里引入判定准则为有效粒子数Neff，并且引入阈值Nth，Nth一般取2N/3，如果Neff＜Nth则进行重采样，否则不进行重采样。

图2 粒子重采样示意Fig.2 Schematic diagram of particle resampling

有效粒子数Neff计算式为：

粒子滤波算法步骤：

（2）对于k=1，2，…，循环以下步骤。

步骤1：重要性采样。对于i=1，2，…，N，从重要概率密度采样q（xk|xk-1，yk），生成N个粒

步骤2：权值计算。按式（26）计算粒子权值，并按式（27）归一化。

步骤4：状态估计。按式（28）计算k时刻的状态估计。

3 模型预测与结果分析

利用粒子滤波算法对压缩机状态参数进行预测时，关键是确定预测模型的状态变量及状态方程，本文所研究的涡旋压缩机性能预测是以压缩气体温度和压力的变化作为涡旋压缩机性能状态的直接表征量，输入变量包括压缩机传热和泄漏等，图3示出了粒子滤波算法预测流程。

图3 粒子滤波算法流程示意Fig.3 Flow chart of particle filtering algorithm

3.1 状态方程建立

为准确估计模型参数，第1.1节推导了涡旋压缩机几何参数影响下压缩气体温度及压力的变化关系式，通过欧拉法来求解模型中相关参数的变化，在模拟计算涡旋压缩机的工作过程中，其性能参数的变化关系式都是关于主轴转角θ的一阶微分方程，因此可采用欧拉法列出温度和压力变化的状态方程。

欧拉基本公式为：

x（t）在t=ti处的一阶泰勒展开式为：

舍去二阶导数之后的余项，并令xi+1≈x（ti+1），xi≈x（ti），得到：

得到：

式中 Δt ——步长。

因此得到温度和压力在任意曲柄转角θj的计算式为：

则涡旋压缩机温度和压力的状态空间模型如式（38）（39）所示。

式中 Qj，Rj——零均值的高斯噪声。

3.2 状态估计精度评价

为了评价算法的估计精度，通过改变粒子数目，采用均方根误差对算法进行精度分析，均方根误差是预测值与真实值偏差的平方和与观测次数比值的平方根，因为该评价模型对一组测量中的特大或特小误差反映非常敏感，所以能够很好地反映出预测精度，其值越小表示预测精度越高［14-17］。均方根误差计算式如下式所示：

式中 T ——迭代步数；

xt，t——预测模型的真实值和估计值。

分别计算3组温度与压力的均方根误差值，结果见表1，2。

表1 温度均方根误差Tab.1 Root mean square error of temperature

表2 压力均方根误差Tab.2 Root mean square error of pressure

3.3 预测结果与分析

本文算法运行环境为MATLAB r2017b，计算机安装Windows10旗舰版64位操作系统［18］。温度仿真过程迭代200步，粒子数N=300，观测噪声Q=0.5，预测噪声R=100；压力仿真过程迭代200步，粒子数N=300，观测噪声Q=5，预测噪声R=300。

标准粒子滤波算法对压缩腔内气体温度变化预测如图4所示，可以看出预测曲线与真实值相近，尤其在吸气阶段和压缩阶段吻合度较高，在排气阶段出现了最大标准偏差值为1.2。对压缩腔中气体压力变化预测如图5所示，可以看出整体预测效果良好，吸气、压缩、排气3个阶段最大标准偏差值分别为2.0，2.4和1.6，从表1，2中可以看出，随着粒子数增加，均方根误差值减小，说明增加粒子数目，使算法稳定性改善，预测精度增加。

图4 温度预测模型Fig.4 Diagram of temperature prediction model

图5 压力预测模型Fig.5 Diagram of pressure prediction model

温度和压力的相对误差如图6所示，可以看出，真实值与预测值相对误差基本在10%左右，说明预测整体性能较好，但是由于算法不稳定，导致误差波动较大，后期可进一步通过优化算法，增加算法稳定性，提升预测精度。

图6 相对误差分析Fig.6 Relative error analysis

4 结论

（1）温度预测过程吸气、压缩、排气3个阶段最大标准偏差分别为0.5，0.4，1.2 ℃，压力预测过程吸气、压缩、排气3个阶段最大标准偏差值分别为2，2.4，1.6 kPa，说明标准粒子滤波算法对压缩气体的温度和压力变化预测曲线与真实值吻合度较高。

（2）预测过程中，随着粒子数从200增加到400，温度变化均方根误差从最大6.204 7减小至3.865 3，压力变化均方根误差从最大14.807 1减小至11.584 3，说明增加粒子数能够有效提高算法预测精度；

（3）从误差分析结果可以看出，真实值与预测值相对误差基本在10%左右，且在气体不断压缩过程中，误差波动逐渐变小，稳定性上升，预测精度提高。