金纵横，李世伟，尚群立

(浙江工业大学 信息工程学院，杭州 310023)

0 引言

摩擦存在于任何现实的机械设备中，其非线性特征是影响控制回路性能的主要因素之一[1-4]。深入研究气动薄膜执行机构模型对开发新型智能电气阀门定位器控制算法，提高生产效益有重要实用价值。

在先前的研究中，针对气动阀中薄膜执行机构的研究偏少且主要研究其摩擦行为[5-6]。文献[7]对气动执行机构和比例方向阀组成的气动系统进行机理建模。文献[8]采集气动执行机构输入和输出做基于数据驱动的系统辨识。对于新型摩擦模型不断被提出。文献[9]提出一种应用于切削过程的新型摩擦模型。文献[10]通过改进LuGre摩擦模型(LuGre-Mod)，并验证该模型可描述纯滑动域中的顺时针和逆时针磁滞回线。但新型摩擦模型并不能显著提高实验数据的预测精度，同时也大幅增加模型复杂度，提高工业生产成本和带来更大的计算工作量[11]。在实际应用中需要在摩擦预测精度和模型复杂度之间做进一步权衡。对于摩擦模型的选择，在工程应用中，仍普遍采用较简单的摩擦模型，例如：Choudhury[12]，Kano[13]，He[14]。文献[15]比较了应用于气动阀的共八种基于物理原理和基于实验数据的摩擦模型，在实际气动阀上按ANSI/ISA标准进行测试，并与理论仿真结果进行比较，结果显示Karnopp、LuGre、Kano三种摩擦模型具有最佳预期效果。对于摩擦模型参数辨识方法，文献[16]提出将拟牛顿法应用于气动阀摩擦模型辨识，进而得到该模型的Stribeck参数。文献[17]针对Karnopp模型参数(m，k，Fc，Fv，Fs)，设计多种辨识算法，但每种算法仅能同时辨识一到两个参数。为了降低测量噪声的影响，需要对采集到的数据做滤波处理。文献[18]对气室气压和阀杆位移用零相位二阶线性滤波器进行滤波。

上述研究从建模方法，提出新的摩擦模型，比较摩擦模型实验效果，对摩擦模型进行更精准的参数辨识等方面去研究气动执行机构等机械设备[19]。本文对Karnopp摩擦模型动静摩擦切换条件进行修改，并将一种基于多元线性回归和最小二乘法的参数辨识方法应用于估算气动薄膜执行机构中移动部件质量，弹簧刚度系数，粘滞摩擦系数，库伦摩擦，弹簧预紧力。通过仿真和现场实验验证参数辨识方法的有效性。结果表明，整合修改后摩擦模型的气动薄膜执行机构动力学模型能准确模拟其动态过程。

1 气动薄膜执行机构模型

工业生产中常见的气开式调节阀，如图1所示。

图1 气动阀示意图

膜片将气室气压转换为力，并与方向相反的弹簧弹力共同决定阀杆的位移，进而改变阀芯和阀座之间的开合程度，控制管路流量。对图1移动部件进行受力分析，得出气动薄膜执行机构的动力学方程[15]：

(1)

式中，m是气动薄膜执行机构移动部件的质量；Fe(t)=S·P(t)，其中Fe(t)是气室压力，S是膜片受力面积，P(t)是气室气压；Fk(t)=k·x(t)，其中Fk(t)是弹簧弹力，k是弹簧刚度系数，x(t)是阀杆位移；Ff(t)是摩擦力；Fl是管道内液体施加在阀杆上的力，在气动阀空载状态下可以忽略；Fi是弹簧预紧力，当气动薄膜执行机构不受气室压力时，用于压紧阀芯与阀座，关闭管路，该力不会随阀杆运动改变。

本次研究在气动薄膜执行机构空载状态下进行，省略后的动力学方程：

(2)

在先前的研究中，不同的作者对Ff(t)有不同的描述，从而提出不同的摩擦模型，但新型摩擦模型包含更多的参数，更复杂，且并不能显著提高仿真精度。行业中仍普遍使用较简单的摩擦模型。文献[10]实际测试八个摩擦模型中，Karnopp、LuGre、Kano三种摩擦模型具有最佳预期效果。结合实际需求和实验可提供的测量数据，本次研究选择Stribeck曲线作为阀杆相对运动速度不为零时所受摩擦力的大小，如图2所示。该曲线中摩擦力的数学表达式可以用式(3)表示：

图2 Stribeck曲线

(3)

(4)

(5)

当速度小于DV时，摩擦力：

Ff(Fr)=min(Fs，|Fr|)sign(Fr)

(6)

Fr=S·P(t)-t·x(t)-Fi

(7)

DV是临界速度，如果阀杆速度小于临界速度，则认为移动部件处于静止状态。式(5)的第一行为移动部件滑动时阀杆所受摩擦力，式(5)的第二行为移动部件静止和动静切换临界时阀杆所受摩擦力。临界速度的目的在于解决零速度检测问题，界定阀杆是处于静止状态还是滑动状态。文献[10]提及临界速度取值标准是观察实验和仿真的拟合度，这需多次调整临界速度值并仿真。

针对临界速度的取值标准需要基于实验和仿真的拟合度这一问题，本文对简化后的Karnopp摩擦模型提出进一步的修改，并对修改后的摩擦模型做离散化处理。

Ff[n]=

(8)

对移动部件处于滑动阶段的气动薄膜执行机构动力学方程做离散化处理：

(9)

式中，

(10)

(11)

n是自然数，T是采样周期。

需要注意的是，在仿真模型中，需指定移动部件初始速度和初始加速度均为零。

2 参数辨识方法

2.1 基于多元线性回归的方法

定义参数向量：

α=[FimFcFvk]

(12)

当阀杆处于滑动阶段时，气动薄膜执行机构离散化动力学方程：

(13)

式(13)和参数向量α是线性的，所以定义回归向量β[n]：

(14)

参数向量α可以通过最小二乘法估计得到：

(15)

为了后续现场实验得出的参数辨识结果作对比，接下来将引入传统参数辨识方法。

2.2 传统方法

刚度系数k：气动薄膜执行机构做开环阶跃测试时，记录每个阶跃测试下阀杆静止时的平均位移和相应的平均气室气压。

(16)

式中，ΔP，Δx分别为两个阶跃测试对应气室气压差和阀杆位移差。

最大静摩擦Fs：气动薄膜执行机构做开环周期和幅值不变的三角或正弦测试时，记录每个周期中两个阀杆位移方向倒转处的位移和对应气室气压。

(17)

式中，ΔP，Δx分别为测试中每个周期两个阀杆位移方向倒转处对应气室气压差和阀杆位移差。

库仑摩擦Fc：气动薄膜执行机构做开环周期三角或梯形斜坡测试时，记录阀杆正向运动和反向运动期间，同一阀杆位移下的气室气压。

(18)

式中，ΔP为阀杆正向运动和反向运动期间，同一阀杆位移下的气室气压差。

粘滞摩擦系数Fv和预紧力Fi：气动薄膜执行机构做开环周期三角或梯形斜坡测试时，记录同一周期阀杆处于正向和反向运动状态下各多组阀杆位移和相应气室气压，并计算该周期正向和反向斜坡速度即为阀杆速度。

(19)

移动部件质量m通常由阀门供应商提供。

为保证结果正确性，所有参数辨识方法中每个测试需做多次，并取平均值。

3 参数辨识实验及分析

3.1 仿真模拟分析

由气动薄膜执行机构动力学方程和摩擦模型方程，通过MATLAB/Simulink搭建仿真模型，该模型如图3所示。

图3 气动薄膜执行机构Simulink模型

该仿真模型能够模拟气动薄膜执行机构的开环运行状态，模型中各参数设定值：m=2 kg，Fv=5 000 N·s/m，Fc=200 N，Fs=200 N，k=200 000 N/m，Fi=1 800 N。为验证摩擦模型修改前后的区别及参数辨识方法的可靠性，进行了如下仿真实验。

仿真实验一：

在开环无测量噪声情况下，原摩擦模型与修改后的摩擦模型在同一正弦激励下的阀杆位移过程的对比。

实验目的：观察原模型与修改后模型之间的区别，验证优化结果。

采样频率：10 000 Hz

图4中，实线表示模型修改后DV=0.000 5 m/s在的阀杆位移(m)随时间(s)变化曲线，虚线代表原模型在DV=0.000 5 m/s的阀杆位移(m)随时间(s)变化曲线。

图4 摩擦模型修改前后阀杆位移对比

图5 加入噪声的仿真结果

实验结果显示，摩擦模型修改前，当阀杆运动临近静止时，阀杆位移呈阶梯式运动。修改后的模型与原模型对比，修改后的模型阀杆位移更平滑。模型修改前后阀杆位移差异并不显著，原因在于在气动薄膜执行机构动力学模型中，相对于弹簧弹力，摩擦力占比并不显著。仿真实验二选择修改后的模型进行实验。

仿真实验二：

在开环无测量噪声和有测量噪声情况下，修改后的气动薄膜执行机构模型在同一激励下的参数辨识结果对比。

实验目的：模拟验证在实际工况下，参数辨识方法的可行性。

采样频率：10 000 Hz

实验具体数据如表1所示。

表1 仿真实验二参数辨识结果

结果表明，在有噪声的情况下，参数辨识方法得出的辨识值误差在可接受范围内，该方法可行。

3.2 现场实验

本章对前文所述的参数辨识方法在实验台架上进行验证，实验台架结构如图6所示。

实验装置：气动阀，I/P转换器，阀位变送器，压力变送器，其中气动阀阀杆行程16 mm，填料为聚四氟乙烯，膜片有效受力面积0.038 m2，阀门移动部件质量1.6 kg，上述气动阀参数均由阀门供应商提供。数据采集装置为美国国家仪器有限公司的NI-DAQ-9332机箱。

为了获取参数辨识所需数据，确保回归向量β[n]中所有系数都能影响气动薄膜执行机构模型动力学方程，通过上位机给定正弦激励信号，同时为了保证数值推算的准确性，设定采样频率为1 000 Hz。实验采集的气室气压，阀杆位移，如图7所示。由于I/P转换器的动态特性，采集的气室压力与给定的正弦激励信号在波形上有些许差异。

图7 实验结果

为了获取气动薄膜执行机构模型参数的最优值，下面将引入三个统计学参数作为评判标准。

式(20)为实际气室压力和估计气室压力的标准误差。

(20)

式(21)为实际气室压力和估计气室压力的相关系数。

(21)

最后一个评判标准是参数向量α中每个参数对应的p-value。当p-value小于0.05时，对应参数符合大于95%置信区间标准。

3.3 实验结果分析

对实验平台得出的实验数据进行分析处理。

其中在不同Δυ下的各项参数辨识值变化曲线如图8所示。

图8 不同Δυ下各参数辨识值

估计外力与实际外力的标准误差和相关系数如图9所示。

图9 估计外力与实际外力的标准误差和相关系数

不同Δυ下的各项参数p-value变化曲线如图10所示。

图10 不同Δυ下各参数p-value

由第二节参数辨识方法，结合图8～10，在Δυ介于[0.004，0.005 8]区间内，图8中参数向量α中各参数辨识值变化趋于平稳，图9中标准误差达到最小范围，相关系数接近于1，图10中各参数p-value趋于0，说明该区间内模型中各参数是显著的。在Δυ=5.72×10-3m/s时，标准误差趋于最小值，相关系数达到最大值，各参数对应p-value均小于0.05。图9所示，当Δυ<4×10-3m/s时，标准误差和相关系数均显著高于区间[0.004，0.005 8]，且Fc出现负数的情况，原因在于阀杆相对运动速度较小的时间段多集中在阀杆位移方向反转的过度区间，该区间气室气压抖动剧烈，如图7气室气压曲线所示，这极大影响辨识的精度。当Δυ>5.8×10-3m/s时，m，Fv，Fc，k，Fi辨识结果值均随着Δυ的递增波动加剧，原因在于随着Δυ的递增，符合参数辨识的样本越少。综上所述，当Δυ=5.72×10-3m/s时，各参数的辨识结果值最接近真实值。

其中基于多元线性回归方法和传统方法进行参数辨识的结果数据对比如表2所示。

表2 现场实验参数辨识结果

表2中粘滞摩擦系数的辨识结果差异较大，原因有以下两点：

1)在气动薄膜执行机构的动力学方程中，相对于占比更高的弹簧弹力，粘滞摩擦力占比极小，对气动薄膜执行机构运动的贡献并不显著。

2)式(19)在估算粘滞摩擦系数和预紧力时，还有弹簧刚度系数和库伦摩擦参与计算，两个变量的估算本身就存在误差。

4 气动薄膜执行机构模型验证

为了验证气动薄膜执行机构模型的可靠性，需要比对仿真和现场实验结果。仿真模型输入信号均为现场采集后经滤波的气室气压。

测试一：开环阶跃输入响应测试，验证模型能否重现死区和评估气动薄膜执行机构动态响应。图11为该激励信号作用下仿真与实验的阀杆位移对照结果。

图11 阶跃输入响应测试

测试二：输入信号为随机信号。选择该信号的原因在于当阀杆处于滑动状态时候，阀杆位移，速度和加速度都随时间不断变化，此时，气动阀模型中所有参数都会影响移动部件的力平衡。图12为该激励信号下仿真与实验的阀杆位置对照结果。

图12 正弦扫频测试

需要注意的是，由于I/P转换器的动态特性，在阀杆位移变小，I/P转换器对外排气时，会出现图11和图12中排气过量导致位移出现尖峰的情况。

图11和图12中仿真和实验阀杆位移没有完全重合的原因有以下两点：

1)实际生产中，气动薄膜执行机构阀杆不同部位与填料接触部分是不均匀的，摩擦系数可变，而人为量化的摩擦系数均为定值。

2)I/P转换器在给气室充排气瞬间，气源供气压力不稳定且气路和气室存在一定程度漏气。

3)实验使用的设备存在不同程度误差。

5 结束语

本研究对应用于气动薄膜执行机构的Karnopp摩擦模型动静摩擦切换条件进行修改。在传统方法的基础是，将多元线性回归和最小二乘法对气动薄膜执行机构参数进行辨识，通过仿真和现场实验验证模型修改前后区别和参数辨识方法的可行性，得出如下结论：

1)应用于气动薄膜执行机构的传统Karnopp摩擦模型存在的实际问题进行分析，并对其动静摩擦切换条件进行修改，解决了DV取值问题。

2)通过带测量噪声的仿真和现场采集气动薄膜执行机构输入输出的真实实验数据，验证基于多元线性回归和最小二乘法的参数辨识方法是有效的。该方法相对于传统方法，极大减少建模工作量。

3)修改后的模型和参数辨识方法并不局限于本研究，经适当修改可应用于其它领域，提高其实用性和普适性。