CN106126808A - 一种基于遗传算法的纯水管路动态摩擦项参数辨识方法 - Google Patents

Abstract

本发明公开了一种基于遗传算法的纯水液压管路瞬态过程中动态摩擦项的参数辨识方法。搭建纯水管路瞬态数学模型，结合特征线法，确定动态摩擦项未知参数的可行解范围，同时应用遗传算法对管路动态摩擦项参数进行优化辨识。本发明提供的参数辨识方法实用可靠，同时结合特征线法对管路瞬态过程进行仿真研究，完善了纯水管路动态摩擦项数学模型，为进一步研究纯水液压系统的流量及压力脉动过程奠定基础。

Description

一种基于遗传算法的纯水管路动态摩擦项参数辨识方法

技术领域

本发明属于纯水液压管路压力脉动研究的技术领域，涉及纯水液压管路瞬态仿真，特别是遗传算法在纯水液压管路压力脉动动态摩擦项参数辨识的应用。

背景技术

在流体传动领域，当管道内稳态流动的液体由于阀门突然关闭或者泵的突然停转，会导致管道内液体压力突然升高，进而产生管路瞬态压力脉动过程。在管路瞬态压力脉动过程中，摩擦力不仅与瞬时的平均速度有关，还与在整个瞬态过程中随时间变化的平均速度有关。如果摩擦力按照恒定流动时摩擦损失处理，则与实际管路的流动情况不符，不能较好的反映管路的瞬态流动过程，即此时管路摩擦力计算公式Darcy-Weisbach公式不适用于压力脉动瞬态过程。因此，对管路动态摩擦阻力项数学模型的研究也是很重要的。但目前，管路摩擦阻力项数学模型的研究还很不完善，尤其是动态摩擦阻力项中加权系数的确定。

本发明应用参数辨识对管路动态摩擦阻力项加权系数进行研究，从而为揭示管路瞬态过程中摩擦项的作用机理提供途径，同时为进一步研究纯水液压系统的流量及压力脉动过程奠定基础。

发明内容

在纯水液压管路瞬态压力脉动研究中，动态摩擦项数学模型涉及多参数。本发明提供一种基于遗传算法的纯水液压管路动态摩擦项的多参数辨识方法，减小管路压力脉动仿真结果和试验结果之间的误差，实现了对管路压力脉动更准确的预测。

本发明是一种基于遗传算法，对纯水液压管路压力脉动过程中动态摩擦项的参数辨识方法；考虑管路边界条件和初始条件，这里研究对象为等径水平直管路，其管路上游为水箱，下游为阀门。初始条件为管路中初始流速恒定。其边界条件为：管路上游压力恒定为水箱压力；下游阀门关闭，因此阀门处流速为零；

因而本发明一种基于遗传算法的纯水管路动态摩擦项参数辨识方法：包括以下步骤：

步骤1.建立纯水液压水平管路一维瞬态数学模型；

根据质量守恒定律和动量守恒定律推导出连续性方程和运动方程，建立水平管路一维瞬态数学模型，如公式1所示:

$\left.\begin{array}{c}\frac{\∂p}{\∂t}+\frac{{\mathrm{\ρa}}^2}{S}\frac{\∂q}{\∂x}=0\\ \frac{\mathrm{\ρ}}{S}\frac{\∂q}{\∂t}+\frac{\∂p}{\∂x}+f\left(q\right)=0\end{array}\right\}---\left(1\right)$

式中：p为管路压强，q为管路流量，S为管路横截面积，ρ为液体密度，a为压力波传播速度，f(q)为摩擦阻力项，t为时间变量，x为管道中位置变量；

步骤2.利用特征线法将管路一维数学模型转化为特征线方程；

利用特征线法对管路瞬态模型偏微分方程进行全微分转化，得到左、右特征线常微分方程组：

$C^{+}:\left\{\begin{array}{c}\frac{\mathrm{\ρ}a}{S}\frac{dq}{dt}+\frac{dp}{dt}+af\left(q\right)=0\\ \frac{dx}{dt}=a\end{array}\right.---\left(2\right)$

$C^{-}:\left\{\begin{array}{c}\frac{\mathrm{\ρ}a}{S}\frac{dq}{dt}-\frac{dp}{dt}+af\left(q\right)=0\\ \frac{dx}{dt}=-a\end{array}\right.---\left(3\right)$

则C⁺，C^-分别表示为左、右特征线方程；

将一段长度为L的水平管道分成n段，每段节点从左到右依次表示为1,2,3……n+1，所以Δx＝L/n，Δt＝Δx/a；A点为t时刻x位置点，B点为t时刻x+2Δx，P为t+Δt时刻x+Δx位置点，F为t-Δt时刻x位置点，G为t-Δt时刻x+2Δx；分别对C⁺和C^-在[A,P]和[P,B]内进行分段积分，得到两个方程，如下

$C^{+}:\frac{\mathrm{\ρ}}{S}{q}_P+\frac{1}{a}{p}_P=\frac{\mathrm{\ρ}}{S}{q}_A+\frac{1}{a}{p}_A-\frac{\mathrm{\Δ}x}{a}f\left(q_A\right)---\left(4\right)$

$C^{-}:\frac{\mathrm{\ρ}}{S}{q}_P-\frac{1}{a}{p}_P=\frac{\mathrm{\ρ}}{S}{q}_B-\frac{1}{a}{p}_B-\frac{\mathrm{\Δ}x}{a}f\left(q_B\right)---\left(5\right)$

令

$C_L=\frac{\mathrm{\ρ}}{S}{q}_A+\frac{1}{a}{p}_A-\frac{\mathrm{\Δ}x}{a}f\left(q_A\right)---\left(6\right)$

$C_R=\frac{\mathrm{\ρ}}{S}{q}_B-\frac{1}{a}{p}_B-\frac{\mathrm{\Δ}x}{a}f\left(q_B\right)---\left(7\right)$

联立(4)、(5)、(6)和(7)四个公式得：

$p_P=\frac{a}{2}(C_L-C_R)---\left(8\right)$

$q_P=\frac{S}{2\mathrm{\ρ}}(C_L+C_R)---\left(9\right)$

式中q_A，q_B，q_P分别为A，B，P点流量，p_A，p_B，p_P分别为A，B，P点压强，f(q_A)，f(q_B)分别为A，B点摩擦项；

步骤3.确定动态摩擦项的未知参数；

摩擦力项f(q)可以表示为：

$f\left(q\right)=f_0+\frac{1}{2}\underset{i=1}{\overset{k}{\Σ}}{Y}_i---\left(10\right)$

式中f₀为稳态摩擦项，Y_i为动态摩擦项，k为加权函数的项数；

根据达西-威斯巴哈(Darcy-Weisbach)公式：

$f_0=\frac{f\mathrm{\ρ}q\left|q\right|}{4{\mathrm{\π}}^2{r}^5}---\left(11\right)$

式中：f为无量纲摩擦系数，r为管道半径；

动态摩擦项Y_i的表示方法为：

$\left\{\begin{array}{c}{Y}_i(t+\mathrm{\Δ}t)=Y_i\left(t\right)e^{-n_i\mathrm{\Δ}\mathrm{\τ}}+m_i{e}^{-n_i(\mathrm{\Δ}\mathrm{\τ}/2)}\[f_0{|}_{t+\mathrm{\Δ}t}-f_0{|}_t\]\\ Y_i\left(0\right)=0\end{array}\right.,(i=1,2,3...k)---\left(12\right)$

式中

其中n_i和m_i为动态摩擦项的加权系数，i为加权函数的项数,μ为粘度；

在式(6)和式(7)中，C_L和C_R的摩擦阻力项f(q_A)和f(q_B)可分别表示为：

$f\left(q_A\right)=\frac{{\mathrm{f\ρq}}_A\left|q_A\right|}{4{\mathrm{\π}}^2{r}^5}+\frac{1}{2}\underset{i=1}{\overset{k}{\Σ}}\[Y_i\left(t\right)e^{-n_i{\mathrm{\Δ\τ}}_A}+\frac{{\mathrm{fm}}_i}{4{\mathrm{\π}}^2{r}^5}{e}^{-n_i({\mathrm{\Δ\τ}}_A/2)}\left({\mathrm{\ρq}}_A\right|q_A|-{\mathrm{\ρq}}_F|q_F\left|\right)\]---\left(14\right)$

$f\left(q_B\right)=\frac{{\mathrm{f\ρq}}_B\left|q_B\right|}{4{\mathrm{\π}}^2{r}^5}+\frac{1}{2}\underset{i=1}{\overset{k}{\Σ}}\[Y_i\left(t\right)e^{-n_i{\mathrm{\Δ\τ}}_B}+\frac{{\mathrm{fm}}_i}{4{\mathrm{\π}}^2{r}^5}{e}^{-n_i({\mathrm{\Δ\τ}}_B/2)}\left({\mathrm{\ρq}}_B\right|q_B|-{\mathrm{\ρq}}_G|q_G\left|\right)\]---\left(15\right)$

式中

${\mathrm{\Δ\τ}}_B=\frac{\mathrm{\Δ}x}{r^{2}a}\frac{{\mathrm{\μ}}_B}{{\mathrm{\ρ}}_B}---\left(17\right)$

式中q_F，q_G分别为F点和G点的流量，其中系数n_i和m_i有如下关系：

m_i＝αm_i-1(i≥3)

n_i＝β²n_i-1(i≥2) (18)

式中：α和β为系数；

根据公式18得知对动态摩擦项需辨识的参数有：m₁，m₂，n₁，α和β；

步骤4.采用遗传算法对步骤3获得的需辨识参数进行优化。

进一步的，所述步骤4的遗传算法具体方法为：

步骤4.1：初始化m₁，m₂，n₁，α和β种群，在可解范围内随机生成n个个体作为初始种群；

步骤4.2：针对生成的初始种群，计算各需辨识参数对应种群中每个个体的适应度；

步骤4.3：选择适应度高的个体，交叉变异产生新一代的种群；

步骤4.4：判断该种群是否满足终止条件，该条件根据实际情况设定若满足终止条件，则从该种群中选出最优的个体即为优化后的需辨识参数；若不满足终止条件则返回步骤4.2。

进一步的，所述步骤4.2中个体适应度的计算方法为：

首先设定目标函数为基于最小二乘法，定义纯水管路瞬态过程的试验和仿真相应的峰谷值之间误差的平方和再开方后为最小；适应度函数为在目标函数的基础上取倒数转化为求全局最大值。

本发明实现了遗传算法对纯水管路瞬态过程中动态摩擦项的参数辨识，使得基于特征线法的纯水管路瞬态仿真研究能够很好地预测纯水管路瞬态过程。

附图说明

图1为纯水液压管路动态摩擦项参数辨识流程图；

图2为特征线法分网示意图；

图3为参数辨识后压力脉动随时间变化仿真曲线与试验曲线对比图。

具体实施方式

为使本发明的目的、技术方案和优点表达得更加清楚明白，下面结合附图及具体实施例对本发明再作进一步详细的说明。

对搭建的一维纯水液压管路瞬态数学模型进行特征线法转化，同时考虑管路两端的边界条件，在阀门突然关闭后，仿真得到阀门处压力脉动随时间变化曲线，如图2所示。管路试验参数如表1所列。

表1管路试验参数

管路内径r	2.1×10-2m
		管路长度L	72m
管壁厚度d	3×10-3m
		流体密度ρ	1000kg/m3
压力波传播速度a	1249m/s
		动力粘度μ	1.1404×10-3Pa·s
初始阀门处压力	4.998×105Pa
		初始流速	0.385m/s

在MATLAB中搭建纯水液压瞬态数学模型，基于所搭建的管路试验平台，根据本发明试验条件，采用四项指数式的动态摩擦项(加权函数的项数k＝4)，其动态摩擦项的表示方法为：

$\left\{\begin{array}{c}{Y}_i(t+\mathrm{\Δ}t)=Y_i\left(t\right)e^{-n_i\mathrm{\Δ}\mathrm{\τ}}+\frac{{\mathrm{fm}}_i}{4{\mathrm{\π}}^2{r}^5}{e}^{-n_i(\mathrm{\Δ}\mathrm{\τ}/2)}\[f_0{|}_{t+\mathrm{\Δ}t}-f_0{|}_t\]\\ Y_i\left(0\right)=0\end{array}\right.,(i=1,2,3...k)---\left(19\right)$

其中n_i和m_i为动态摩擦项加权系数。

系数n_i和m_i有如下关系：

m_i＝αm_i-1(i≥3,4) (20)

n_i＝β²n_i-1(i≥2,3,4)

式中α和β为系数。

因此对动态摩擦项需辨识的参数有：m₁，m₂，n₁，α和β。采用MATLAB中遗传算法工具箱，基于管路瞬态试验数据来验证本发明的可行性。

确定动态摩擦辨识参数m₁，m₂，n₁，α和β的可行解范围。根据实测管路脉动试验数据，可以确定m₁，m₂，n₁，α和β的可行解范围，分别为[1,10]、[1,10]、[10,100]、[1,10]、[1,10]。初始化基本遗传算法的相关控制参数，选择参数采用十进制编码和解码方式，然后，对初始种群、终止条件的设定、交叉概率和变异概率设定等；接着，设计适应度函数，遗传算法是利用个体对应的适应度函数值来评价个体的优劣，本发明定义适应度函数F：

$F=\frac{{10}^7}{\sqrt{H}}---\left(21\right)$

目标函数：

$H=\underset{i=1}{\overset{w}{\Σ}}{\left(p_{ex}\right(i)-p_{sim}(i\left)\right)}^2---\left(22\right)$

式中：w为实验数据中压力波峰数，p_ex(i)和p_sim(i)分别是试验曲线和仿真曲线第i次压力脉动极值。

这样，管路动态摩擦项的参数辨识问题转化为求目标函数H的最小值问题。根据经过编码的参数可行解随机产生N个初始解，以这N个初始解作为初始点开始迭代。将初始种群中的个体随机组队后，进行选择、交叉和变异等遗传操作，通过迭代计算得到优化个体并检测进化代数进行终止条件判断，当满足终止条件，则输出全局最优解；否则，继续进行迭代优化。经过反复判断，到终止代数，解码，输出最终辨识结果。

本发明中，设置初始种群数量为100，繁衍代数为80代。参数辨识后的结果：m₁＝2.0368，m₂＝7.1676，n₁＝50.5294，α＝2.0269,β＝3.3204。

根据上面的辨识结果和式(18)，可以得到m_i和n_i值，如表2所列。

表2辨识后的m_i和n_i值

图3为参数辨识后压力脉动随时间变化仿真曲线与试验曲线对比图。由图可知，辨识后压力脉动仿真结果与试验数据基本吻合，包括压力脉动峰值和脉动时间间隔，说明本发明方法可行实用。

Claims (3)

1.一种基于遗传算法的纯水管路动态摩擦项参数辨识方法：包括以下步骤：

步骤1.建立纯水液压水平管路一维瞬态数学模型；

根据质量守恒定律和动量守恒定律推导出连续性方程和运动方程，建立水平管路一维瞬态数学模型，如公式1所示:

$\left.\begin{array}{c}\frac{\∂p}{\∂t}+\frac{{\mathrm{\ρa}}^2}{S}\frac{\∂q}{\∂x}=0\\ \frac{\mathrm{\ρ}}{S}\frac{\∂q}{\∂t}+\frac{\∂p}{\∂x}+f\left(q\right)=0\end{array}\right\}---\left(1\right)$

式中：p为管路压强，q为管路流量，S为管路横截面积，ρ为液体密度，a为压力波传播速度，f(q)为摩擦阻力项，t为时间变量，x为管道中位置变量；

步骤2.利用特征线法将管路一维数学模型转化为特征线方程；

利用特征线法对管路瞬态模型偏微分方程进行全微分转化，得到左、右特征线常微分方程组：

$C^{+}:\left\{\begin{array}{c}\frac{\mathrm{\ρ}a}{S}\frac{dq}{dt}+\frac{dp}{dt}+af\left(q\right)=0\\ \frac{dx}{dt}=a\end{array}\right.---\left(2\right)$

$C^{-}:\left\{\begin{array}{c}\frac{\mathrm{\ρ}a}{S}\frac{dq}{dt}-\frac{dp}{dt}+af\left(q\right)=0\\ \frac{dx}{dt}=-a\end{array}\right.---\left(3\right)$

则C⁺，C^-分别表示为左、右特征线方程；

将一段长度为L的水平管道分成n段，每段节点从左到右依次表示为1,2,3……n+1，所以Δx＝L/n，Δt＝Δx/a；A点为t时刻x位置点，B点为t时刻x+2Δx，P为t+Δt时刻x+Δx位置点，F为t-Δt时刻x位置点，G为t-Δt时刻x+2Δx；分别对C⁺和C^-在[A,P]和[P,B]内进行分段积分，得到两个方程，如下

$C^{+}:\frac{\mathrm{\ρ}}{S}{q}_P+\frac{1}{a}{p}_P=\frac{\mathrm{\ρ}}{S}{q}_A+\frac{1}{a}{p}_A-\frac{\mathrm{\Δ}x}{a}f\left(q_A\right)---\left(4\right)$

$C^{-}:\frac{\mathrm{\ρ}}{S}{q}_P-\frac{1}{a}{p}_P=\frac{\mathrm{\ρ}}{S}{q}_B-\frac{1}{a}{p}_B-\frac{\mathrm{\Δ}x}{a}f\left(q_B\right)---\left(5\right)$

令

$C_L=\frac{\mathrm{\ρ}}{S}{q}_A+\frac{1}{a}{p}_A-\frac{\mathrm{\Δ}x}{a}f\left(q_A\right)---\left(6\right)$

$C_R=\frac{\mathrm{\ρ}}{S}{q}_B-\frac{1}{a}{p}_B-\frac{\mathrm{\Δ}x}{a}f\left(q_B\right)---\left(7\right)$

联立(4)、(5)、(6)和(7)四个公式得：

$p_P=\frac{a}{2}(C_L-C_R)---\left(8\right)$

$q_P=\frac{S}{2\mathrm{\ρ}}(C_L+C_R)---\left(9\right)$

式中q_A，q_B，q_P分别为A，B，P点流量，p_A，p_B，p_P分别为A，B，P点压强，f(q_A)，f(q_B)分别为A，B点摩擦项；

步骤3.确定动态摩擦项的未知参数；

摩擦力项f(q)可以表示为：

$f\left(q\right)=f_0+\frac{1}{2}\underset{i=1}{\overset{k}{\Σ}}{Y}_i---\left(10\right)$

式中f₀为稳态摩擦项，Y_i为动态摩擦项，k为加权函数的项数；

根据达西-威斯巴哈(Darcy-Weisbach)公式：

$f_0=\frac{f\mathrm{\ρ}q\left|q\right|}{4{\mathrm{\π}}^2{r}^5}---\left(11\right)$

式中：f为无量纲摩擦系数，r为管道半径；

动态摩擦项Y_i的表示方法为：

$\left\{\begin{array}{c}{Y}_i(t+\mathrm{\Δ}t)=Y_i\left(t\right)e^{-n_i\mathrm{\Δ}\mathrm{\τ}}+m_i{e}^{-n_i(\mathrm{\Δ}\mathrm{\τ}/2)}\[f_0{|}_{t+\mathrm{\Δ}t}-f_0{|}_t\]\\ Y_i\left(0\right)=0\end{array}\right.,(i=1,2,3...k)---\left(12\right)$

式中

其中n_i和m_i为动态摩擦项的加权系数，i为加权函数的项数,μ为粘度；

在式(6)和式(7)中，C_L和C_R的摩擦阻力项f(q_A)和f(q_B)可分别表示为：

$f\left(q_A\right)=\frac{{\mathrm{f\ρq}}_A\left|q_A\right|}{4{\mathrm{\π}}^2{r}^5}+\frac{1}{2}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\[Y_i\left(t\right)e^{-n_i{\mathrm{\Δ\τ}}_A}+\frac{{\mathrm{fm}}_i}{4{\mathrm{\π}}^2{r}^5}{e}^{-n_i\left({\mathrm{\Δ\τ}}_A/2\right)}\left({\mathrm{\ρq}}_A\left|q_A\right|-{\mathrm{\ρq}}_F\left|q_F\right|\right)\]---\left(14\right)$

$f\left(q_B\right)=\frac{{\mathrm{f\ρq}}_B\left|q_B\right|}{4{\mathrm{\π}}^2{r}^5}+\frac{1}{2}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\[Y_i\left(t\right)e^{-n_i{\mathrm{\Δ\τ}}_B}+\frac{{\mathrm{fm}}_i}{4{\mathrm{\π}}^2{r}^5}{e}^{-n_i\left({\mathrm{\Δ\τ}}_B/2\right)}\left({\mathrm{\ρq}}_B\left|q_B\right|-{\mathrm{\ρq}}_G\left|q_G\right|\right)\]---\left(15\right)$

式中

${\mathrm{\Δ\τ}}_B=\frac{\mathrm{\Δ}x}{r^{2}a}\frac{{\mathrm{\μ}}_B}{{\mathrm{\ρ}}_B}---\left(17\right)$

式中q_F，q_G分别为F点和G点的流量，其中系数n_i和m_i有如下关系：

m_i＝αm_i-1 (i≥3)

n_i＝β²n_i-1 (i≥2) (18)

式中：α和β为系数；

根据公式18得知对动态摩擦项需辨识的参数有：m₁，m₂，n₁，α和β；

步骤4.采用遗传算法对步骤3获得的需辨识参数进行优化。

2.如权利要求1所述的一种基于遗传算法的纯水管路动态摩擦项参数辨识方法，其特征在于所述步骤4的遗传算法具体方法为：

步骤4.1：初始化m₁，m₂，n₁，α和β种群，在可解范围内随机生成n个个体作为初始种群；

步骤4.2：针对生成的初始种群，计算各需辨识参数对应种群中每个个体的适应度；

步骤4.3：选择适应度高的个体，交叉变异产生新一代的种群；

步骤4.4：判断该种群是否满足终止条件，该条件根据实际情况设定若满足终止条件，则从该种群中选出最优的个体即为优化后的需辨识参数；若不满足终止条件则返回步骤4.2。

3.如权利要求2所述的一种基于遗传算法的纯水管路动态摩擦项参数辨识方法，其特征在于所所述步骤4.2中个体适应度的计算方法为：

首先设定目标函数为基于最小二乘法，定义纯水管路瞬态过程的试验和仿真相应的峰谷值之间误差的平方和再开方后为最小；适应度函数为在目标函数的基础上取倒数转化为求全局最大值。