CN105046326A - 基于功图主元分析的抽油机参数优化方法 - Google Patents

Abstract

本发明提供一种基于功图主元分析的抽油机参数优化方法，包括：1)确定抽油机的效率影响因素和性能变量、2)获得其样本数据、3)对载荷数据进行降维、4)由非载荷变量与载荷新主元构建网络输入变量、5)构建输入变量的样本值、6)归一化处理、7)构建前馈神经网络、8)利用无迹卡尔曼滤波对网络进行训练、9)构造父代和子代种群；10)对父代个体作遗传变异计算，以产生子代个体；11)对父代和子代个体求适应度函数；12)将父代和子代个体划分到层级不同的非支配集中；13)从这些非支配集中选择个体构成新的父代种群，循环10)-13)多次，得到优化后的效率影响因素值。优化后，可得到产液量最大时，耗电量最小。

Description

基于功图主元分析的抽油机参数优化方法

技术领域

本发明属于采油领域，具体涉及一种基于功图主元分析的抽油机参数优化方法。

背景技术

抽油机采油作为一种机械采油方式，主要由电动机、地面传动设备和井下抽油设备三部分组成，如图1所示。抽油机的整个采油过程主要分为上下两个冲程：上冲程，即驴头悬点向上运动，提起抽油杆柱和井下抽油设备，此过程中电动机需消耗大量的能量；下冲程，即驴头悬点向下运动，抽油机的抽油杆柱电动机做功。在抽油杆柱上下运动过程中，电动机的负载发生周期变化。抽油机的运行参数的选择对整个抽油机系统的能量消耗影响很大。为了使抽油机采油生产过程既能完成预定的产液量，又能使抽油机生产过程的耗电量最低，需要对抽油机运行参数进行节能优化。

发明内容

本发明是为了解决现有技术中存在的上述技术问题而做出，其目的在于提供一种基于功图主元分析的抽油机参数优化方法，以保证抽油机的生产状态最佳，从而达到减少能耗，提高系统效率的目的。

为了实现上述目的，本发明提供一种基于功图主元分析的抽油机参数优化方法，该方法包括的步骤如下：

1)确定抽油机采油过程生产效率影响因素构成效率观测变量集合其中α₁，α₂为决策变量，α₃～α₁₄₆载荷数据环境变量，为其他环境变量，选取抽油机系统的性能变量构成性能观测变量集合：{y₁,y₂,y₃,…y_l}；

2)获得所述生产效率影响因素和系统性能变量的样本数据，得到效率影响因素样本矩阵α和性能样本矩阵Y：

$\mathrm{\α}=\left[\begin{array}{cccc}{\mathrm{\α}}_{11}& {\mathrm{\α}}_{12}& ...& {\mathrm{\α}}_{1N}\\ {\mathrm{\α}}_{21}& {\mathrm{\α}}_{22}& ...& {\mathrm{\α}}_{2N}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{\stackrel{\‾}{M}1}& {\mathrm{\α}}_{\stackrel{\‾}{M}2}& ...& {\mathrm{\α}}_{\stackrel{\‾}{M}N}\end{array}\right]\≡\left[\begin{array}{c}{L}_1\\ L_2\\ ...\\ L_{\stackrel{\‾}{M}}\end{array}\right]$

$Y=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{22}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

其中为效率影响因素个数，N为样本个数，α_ik表示第i个效率影响因素变量的第k个观测值，i＝1,2,...,M；k＝1,2,...,N；

3)利用主元分析算法对载荷数据进行降维处理，从而构建新的载荷主元变量矩阵：

$\left[\begin{array}{c}{L}_{z1}\\ L_{z2}\\ L_{z3}\\ L_{\mathrm{zd}}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\α}}_{z11}& {\mathrm{\α}}_{z12}& ...& {\mathrm{\α}}_{z1N}\\ {\mathrm{\α}}_{z21}& {\mathrm{\α}}_{z22}& ...& {\mathrm{\α}}_{z2N}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{\mathrm{zd}1}& {\mathrm{\α}}_{\mathrm{zd}2}& ...& {\mathrm{\α}}_{\mathrm{zdN}}\end{array}\right]$

4)由影响因素观测变量集合中非载荷变量与载荷新主元观测变量集合{α_z1,α_z2,...,α_zd}构建网络输入变量集合：并令输入变量集合为：{x₁,x₂,x₃,...,x_M}，即，

5)构建输入变量集合{x₁,x₂,x₃,...,x_M}观测样本值：

$X=\left[\begin{array}{cccc}{X}_1& X_2& ...& X_N\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\α}}_{11}& {\mathrm{\α}}_{12}& ...& {\mathrm{\α}}_{1N}\\ {\mathrm{\α}}_{21}& {\mathrm{\α}}_{22}& ...& {\mathrm{\α}}_{2N}\\ {\mathrm{\α}}_{\mathrm{147,1}}& {\mathrm{\α}}_{\mathrm{147,2}}& ...& {\mathrm{\α}}_{147,N}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{\stackrel{\‾}{M}1}& {\mathrm{\α}}_{\stackrel{\‾}{M}2}& ...& {\mathrm{\α}}_{\stackrel{\‾}{M}N}\\ {\mathrm{\α}}_{z11}& {\mathrm{\α}}_{z12}& ...& {\mathrm{\α}}_{z1N}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{\mathrm{zd}1}& {\mathrm{\α}}_{\mathrm{zd}2}& ...& {\mathrm{\α}}_{\mathrm{zdN}}\end{array}\right]=\left[\begin{array}{cccc}{x}_{11}& x_{12}& ...& x_{1N}\\ x_{21}& x_{22}& ...& x_{2N}\\ x_{31}& x_{32}& ...& x_{3N}\\ ...& ...& ...& ...\\ x_{M1}& x_{M2}& ...& x_{\mathrm{MN}}\end{array}\right]$

$Y=\left[\begin{array}{cccc}{Y}_1& Y_2& ...& Y_N\end{array}\right]=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{22}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

其中，x₁～x₂为决策变量，x₃～x_M为新的环境变量；

6)对得到的训练输入样本X、输出样本Y进行归一化处理，得到新的训练输入矩阵输出矩阵

${\tilde{X}}_k={\left[\begin{array}{ccccc}{\tilde{x}}_{1k}& {\tilde{x}}_{2k}& {\tilde{x}}_{3k}& ...& {\tilde{x}}_{\mathrm{Mk}}\end{array}\right]}^T={\left[\begin{array}{ccccc}f\left(x_{1k}\right)& f\left(x_{2k}\right)& f\left(x_{3k}\right)& ...& f\left(x_{\mathrm{Mk}}\right)\end{array}\right]}^T\≡f\left(X_k\right)$

${\tilde{Y}}_k={\left[\begin{array}{cccc}{\tilde{y}}_{1k}& {\tilde{y}}_{2k}& ...& {\tilde{y}}_{\mathrm{lk}}\end{array}\right]}^T={\left[\begin{array}{cccc}g\left(y_{1k}\right)& g\left(y_{2k}\right)& ...& g\left(y_{\mathrm{lk}}\right)\end{array}\right]}^T\≡g\left(Y_k\right)$

$\tilde{X}=[{\tilde{X}}_1,{\tilde{X}}_2,...,{\tilde{X}}_N]=\left[\begin{array}{cccc}{\tilde{x}}_{11}& {\tilde{x}}_{12}& ...& {\tilde{x}}_{1N}\\ {\tilde{x}}_{21}& {\tilde{x}}_{22}& ...& {\tilde{x}}_{2N}\\ ...& ...& ...& ...\\ {\tilde{x}}_{M1}& {\tilde{x}}_{M2}& ...& {\tilde{x}}_{\mathrm{MN}}\end{array}\right]$

$\tilde{Y}=\left[\begin{array}{cccc}{\tilde{Y}}_1& {\tilde{Y}}_2& ...& {\tilde{Y}}_N\end{array}\right]=\left[\begin{array}{cccc}{\tilde{y}}_{11}& {\tilde{y}}_{12}& ...& {\tilde{y}}_{1N}\\ {\tilde{y}}_{21}& {\tilde{y}}_{22}& ...& {\tilde{y}}_{2N}\\ ...& ...& ...& ...\\ {\tilde{y}}_{l1}& {\tilde{y}}_{l2}& ...& {\tilde{y}}_{\mathrm{lN}}\end{array}\right]$

7)构建三层前馈神经网络，其输入变量集为输出变量集为隐含层神经元个数为s₁，输入层、隐含层、输出层通过权值、阈值进行连接，并且该神经网络的输入输出函数表达式为：

式中函数F(X)为S型函数；

8)利用无迹卡尔曼滤波对所述前馈神经网络进行训练，得到该神经网络的结构参数值，该训练过程包括：

①将所述神经网络中的所有权值和阈值组成状态变量I：

$I={[w_{11}^{1}L{w}_{{\mathrm{Ms}}_1}^1,b_1^{1}L{b}_{s1}^1,w_{11}^{2}L{w}_{s_{1}l}^2,b_1^{2}L{b}_l^2]}^T$

其中，M为输入层神经元数，s₁为隐层神经元数，l为输出层神经元数，输入层至隐层神经元的连接权值为阈值为隐层至输出层的连接权值为阈值为I中的元素个数为n；设定非线性方程：

其中，函数表达式参考步骤7)，为K时刻的神经网络归一化输入样本，令ω_k＝0，v_k＝0，为神经网络归一化输出样本；

②设定无迹卡尔曼计算过程中控制采样点的分布状态参数a、待选参数κ，以及非负权系数β；

③计算2n+1个采样点σ点以及σ点的相应权重，其中n为状态矩阵的I维度，λ＝a²(n+κ)-n，2n+1个采样点计算如下：

$\left\{\begin{array}{c}{I}^{\left(0\right)}=\stackrel{\‾}{I}\\ I^{\left(i\right)}=\stackrel{\‾}{I}+\sqrt{(n+\mathrm{\λ})P},i=1:n\\ I^{\left(i\right)}=\stackrel{\‾}{I}-\sqrt{(n+\mathrm{\λ})P},i=n+1:2n\end{array}\right.$

每个采样点的权值计算如下：

$\left\{\begin{array}{c}{W}_m^{\left(0\right)}=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}\\ w_0^{\left(0\right)}=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-a+\mathrm{\β})\\ W_m^{\left(i\right)}=W_c^{\left(i\right)}=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},i=1:2n\end{array}\right.$

④计算σ点的一步状态预测及状态变量协方差P_k+1|k；

$\begin{array}{c}{I}_{k+1|k}^{\left(i\right)}=I_{k|k}^{\left(i\right)}\\ {\hat{I}}_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_m^{\left(i\right)}\·I_{k+1|k}^{\left(i\right)}\\ P_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_c^{\left(i\right)}\·[I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}^{\left(i\right)}]{[I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}^{\left(i\right)}]}^T\end{array}$

⑤计算输出的一步提前预测以及协方差

$\begin{array}{c}{Y}_{k+1|k}^{\left(i\right)}=h(I_{k+1|k}^{\left(i\right)},X_k)\\ {\hat{Y}}_{k+1|k}=\stackrel{2n}{\underset{i=0}{\mathrm{\Σ}}{W}_m^{\left(i\right)}\·}h(I_{k+1|k}^{\left(i\right)},X_k)\\ P_{Y_{k+1}}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{w}_c^{\left(i\right)}(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k}){(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k})}^T\\ P_{X_{k+1}{Y}_{k+1}}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_c^{\left(i\right)}(I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}){(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k})}^T\end{array}$

⑥进行滤波更新获取新的状态矩阵、协方差矩阵、增益矩阵：

$\begin{array}{c}{R}_{k+1}=P_{X_{k+1}{Y}_{k+1}}/P_{Y_{k+1}}\\ {\hat{I}}_{k+1|k+1}={\hat{I}}_{k+1|k}+R_{k+1}(Y_{k+1}-{\hat{Y}}_{k+1|k})\\ P_{k+1|k+1}=P_{k+1|k}-R_{k+1}{\left(P_{X_{k+1}{Z}_{k+1}}\right)}^T\end{array}$

⑦对获取的新样本数据重新进行②～⑥步骤，直至所有样本对状态矩阵、协方差矩阵、增益矩阵进行了更新，从而得到适应于所有样本状态矩阵；

⑧对最后一组样本得到状态矩阵I，作为网络训练得到的权值和阈值；

⑨在得到网络参数各层权值、阈值后，确定所述前馈神经网络的函数模型为：

$\hat{Y}\left(X\right)=g^{-1}\left(\tilde{Y}\left(\tilde{X}\right)\right)=g^{-1}\left(\tilde{Y}\left(f\left(X\right)\right)\right)$

9)利用决策变量(x₁,x₂)构建多目标优化过程父代种群P_D，

$P_D=\left\{(x_{1m}^{P_D},x_{2m}^{P_D})\right|1\≤m\≤K\}$

其中，父代种群P_D中的个体的数量为K，并从x₁的取值范围x_1,min≤x₁≤x_1,max内随机取值赋予从x₂的取值范围x_2,min≤x₂≤x_2,max内随机取值赋予从而对父代种群P_D进行初始化；

10)从父代种群P_D中选出任意对个体，对于每对个体 进行遗传交叉计算或变异计算，并将计算结果赋予子代种群Q_D中相应的一对个体

11)将父代种群P_D与子代种群Q_D进行合并得到种群R＝P_D∪Q_D，即有 $R=\left\{(x_{1s}^R,x_{2s}^R)\right|1\≤s\≤2K\}=\{(x_{1m}^{P_D},x_{2m}^{P_D})|1\≤m\≤K\}\∪\left\{(x_{1n}^{Q_D},x_{2n}^{Q_D})\right|1\≤n\≤K\},$ 将种群R的每个个体与环境变量平均值i＝3,…,M合成输入样本 $X_s={\left[\begin{array}{ccccc}{x}_{1s}^R& x_{2s}^R& {\stackrel{\‾}{x}}_3& ...& {\stackrel{\‾}{x}}_M\end{array}\right]}^T,$ 并计算相应的适应度函数 $\mathrm{objFun}\left(X_s\right)={\left[\begin{array}{cc}h\left({\hat{y}}_1\left(X_s\right)\right)& {\hat{y}}_2\left(X_s\right)\end{array}\right]}^T,$ 其中，函数h(x)＝-x；

12)将种群R的所有个体所对应的适应度函数相互进行比较，将种群R的所有个体划分到具有不同层级的非支配集中，其中，对于层级较低的非支配集中的任一个体所对应的适应度函数objFun(X_s)和层级较高的非支配集中的任一个体所对应的适应度函数objFun(X_t)来说，均不存在且而对于同一层级的非支配集中的任两个个体来说，该两个不等式中至少有一个不成立；

13)按照层级从低到高的顺序从所述非支配集中选择K个个体，将选择出的K个个体的值赋予父代种群P_D中的个体，并执行步骤10)-步骤13)的过程GEN次，GEN为预先确定的循环次数，最终得到优化后的L组决策变量将优化后的决策变量以及所述环境变量的平均值构成优化后的输入样本 $X_m^{P_D^{\mathrm{GEN}}}={\left[\begin{array}{ccccc}{x}_{1m}^{P_D^{\mathrm{GEN}}}& x_{2m}^{P_D^{\mathrm{GEN}}}& {\stackrel{\‾}{x}}_3& ...& {\stackrel{\‾}{x}}_M\end{array}\right]}^T(1\≤m\≤K),$ 这K个样本保证了在产液量最大，耗电量最小。

本发明的有益效果是，利用功图主元分析方法进行降维，以简化计算过程，利用无迹卡尔曼滤波神经网络(UKFNN)建立油田机采过程实施动态演化高精度模型，并利用带精英策略的快速非支配排序遗传算法对建立的模型进行搜索，探寻抽油机生产过程中最佳工艺决策参数，给出面向节能降耗的抽油机生产最佳参数指导生产，从而达到节能降耗目的。通过选择优化后的运行参数，可以使抽油机在运行过程中保证在产液量基本固定的情况下，耗电量最小，从而可以降低油田生产成本并提高油田生产效率。

附图说明

图1示出了抽油机的工作模型；

图2示出了本发明一个实施例所述的基于功图主元分析的抽油机参数优化方法的流程图；

图3示出了本发明一个实施例中的前馈神经网络的结构；

图4示出了测试样本产液量预测效果图；

图5示出了测试样本耗电量预测效果图；

图6示出了产液量偏好值与耗电量的pareto解集关系。

具体实施方式

在下面的描述中，出于说明的目的，为了提供对一个或多个实施例的全面理解，阐述了许多具体细节。然而，很明显，也可以在没有这些具体细节的情况下实现这些实施例。在其它例子中，为了便于描述一个或多个实施例，公知的结构和设备以方框图的形式示出。

图2是流程图，示出了本发明的一个实施例所述的基于功图主元分析的抽油机参数优化方法。如图2所示，本发明所述的基于功图主元分析的抽油机参数优化方法包括如下步骤：

步骤S1：确定抽油机采油过程生产效率影响因素构成效率观测变量集合其中α₁，α₂为决策变量，α₃～α₁₄₆载荷数据环境变量，为其他环境变量，选取抽油机系统的性能变量构成性能观测变量集合：{y₁,y₂,y₃,…y_l}。

在本发明的一个实施例中，选取决策变量α₁为冲次、决策变量α₂为有效冲程、α₃～α₁₄₆为载荷1至载荷144，其余环境变量包括：理论排量、功率因数、有功功率、无功功率、含水率中的一个或多个变量；选取抽油机生产过程性能变量y₁为产液量、y₂为耗电量。

步骤S2：获得所述生产效率影响因素和系统性能变量的样本数据，得到效率影响因素样本矩阵α和性能样本矩阵Y：

$\mathrm{\α}=\left[\begin{array}{cccc}{\mathrm{\α}}_{11}& {\mathrm{\α}}_{12}& ...& {\mathrm{\α}}_{1N}\\ {\mathrm{\α}}_{21}& {\mathrm{\α}}_{22}& ...& {\mathrm{\α}}_{2N}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{\stackrel{\‾}{M}1}& {\mathrm{\α}}_{\stackrel{\‾}{M}2}& ...& {\mathrm{\α}}_{\stackrel{\‾}{M}N}\end{array}\right]\≡\left[\begin{array}{c}{L}_1\\ L_2\\ ...\\ L_{\stackrel{\‾}{M}}\end{array}\right]$

$Y=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{22}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

其中为效率影响因素个数，N为样本个数，α_ik表示第i个效率影响因素变量的第k个观测值，i＝1,2,...,M；k＝1,2,...,N。

设所述决策变量、环境变量和性能变量的观测值采集周期的最大值为tmax，则这些变量中的任一变量的样本取为tmax时间内该变量的观测值的平均值。

步骤S3：利用主元分析算法对载荷数据进行降维处理，从而构建新的载荷主元变量。本发明中采用示功图描绘数据的144个载荷点做为部分环境变量进行建模。然而利用144维数据建模易产生维度灾难。故而利用主元分析算法对载荷数据进行降维处理。

在一个实施例中，利用主元分析算法对载荷数据进行降维处理的步骤可以包括：

①设置样本累计贡献率precent＝0.95；

②获取载荷数据每个L_k具有第k观测变量的N个观测数据,3≤k≤146；

③求出数据平均值并利用原始数据减去均值得到

④计算协方差矩阵

⑤计算协方差矩阵的特征值E₁,E₂,...,E₁₄₄与特征向量EV₁,EV₂,...,EV₁₄₄；

⑥由大到小依次排列特征值E′₁,E'₂,...,E'_M,对应特征向量为EV′₁,EV′₂,...,EV′₁₄₄，按特征值大小顺序取前d个特征值的特征向量构成矩阵[EV₁',EV₂',...,EV_d']，此时其中d<144；特征向量代表原数据的分布方向，其对应的特征值越大，则该向量越重要(即为主元)；其对应的特征值越小，则该向量越次要。

⑦由[EV′₁,EV′₂,...,EV′_d]与原始样本求取载荷新的主元，其新载荷主元观测变量构成集合:{α_z1,α_z2,...,α_zd}，其为d个新变量，且每个变量为N个观测值构成的新主元矩阵：

$\left[\begin{array}{c}{L}_{z1}\\ L_{z2}\\ L_{z3}\\ L_{\mathrm{zd}}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_{z11}& {\mathrm{\&alpha;}}_{z12}& ...& {\mathrm{\&alpha;}}_{z1N}\\ {\mathrm{\&alpha;}}_{z21}& {\mathrm{\&alpha;}}_{z22}& ...& {\mathrm{\&alpha;}}_{z2N}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{\mathrm{zd}1}& {\mathrm{\&alpha;}}_{\mathrm{zd}2}& ...& {\mathrm{\&alpha;}}_{\mathrm{zdN}}\end{array}\right]$

步骤S4：由影响因素观测变量集合中非载荷变量与载荷新主元观测变量集合{α_z1,α_z2,...,α_zd}构建网络输入变量集合：并令输入变量集合为：{x₁,x₂,x₃,...,x_M}，即，

步骤S5：构建输入变量集合{x₁,x₂,x₃,...,x_M}观测样本值：

其中，x₁～x₂为决策变量，x₃～x_M为新的环境变量。

步骤S6：对得到的训练输入样本X、输出样本Y进行归一化处理，得到新的训练输入矩阵输出矩阵

在一个实施例中，所述归一化处理的算法如下：

$\begin{array}{c}{\tilde{x}}_{\mathrm{ik}}=f\left(x_{\mathrm{ik}}\right)=({\tilde{x}}_{\max }-{\tilde{x}}_{\min })\&CenterDot;\frac{x_{\mathrm{ik}}-x_{i,\min }}{x_{i,\max }-x_{i,\min }}+{\tilde{x}}_{\min }\\ i=\mathrm{1,2},...,M;k=\mathrm{1,2},...,N\\ {\tilde{y}}_{\mathrm{jk}}=g\left(y_{\mathrm{jk}}\right)=({\tilde{y}}_{\max }-{\tilde{y}}_{\min })\&CenterDot;\frac{y_{\mathrm{jk}}-y_{j,\min }}{y_{j,\max }-y_{j,\min }}+{\tilde{y}}_{\min }\\ j=\mathrm{1,2},...,l;k=\mathrm{1,2},...,N\end{array}$

其中:为设定输入变量归一化后数据范围的最大值、最小值；

x_ik为归一化前的第i个输入变量第k个样本值；

为归一化后第i个输入变量第k个样本值；

x_i,min＝min{x_ik|1≤k≤N}

x_i,max＝max{x_ik|1≤k≤N}

为设定输出变量归一化后数据范围的最大值、最小值；

y_jk为归一化前第j个输出变量的第k个采集样本值；

为归一化后第j个输出变量的第k个值；

y_j,max＝max{y_jk|1≤k≤N}

y_j,min＝min{y_jk|1≤k≤N}

于是得到：

${\tilde{X}}_k={\left[\begin{array}{ccccc}{\tilde{x}}_{1k}& {\tilde{x}}_{2k}& {\tilde{x}}_{3k}& ...& {\tilde{x}}_{\mathrm{Mk}}\end{array}\right]}^T={\left[\begin{array}{ccccc}f\left(x_{1k}\right)& f\left(x_{2k}\right)& f\left(x_{3k}\right)& ...& f\left(x_{\mathrm{Mk}}\right)\end{array}\right]}^T\&equiv;f\left(X_k\right)$

$\tilde{X}=[{\tilde{X}}_1,{\tilde{X}}_2,...,{\tilde{X}}_N]=\left[\begin{array}{cccc}{\tilde{x}}_{11}& {\tilde{x}}_{12}& ...& {\tilde{x}}_{1N}\\ {\tilde{x}}_{21}& {\tilde{x}}_{22}& ...& {\tilde{x}}_{2N}\\ ...& ...& ...& ...\\ {\tilde{x}}_{M1}& {\tilde{x}}_{M2}& ...& {\tilde{x}}_{\mathrm{MN}}\end{array}\right]$

$\tilde{Y}=\left[\begin{array}{cccc}{\tilde{Y}}_1& {\tilde{Y}}_2& ...& {\tilde{Y}}_N\end{array}\right]=\left[\begin{array}{cccc}{\tilde{y}}_{11}& {\tilde{y}}_{12}& ...& {\tilde{y}}_{1N}\\ {\tilde{y}}_{21}& {\tilde{y}}_{22}& ...& {\tilde{y}}_{2N}\\ ...& ...& ...& ...\\ {\tilde{y}}_{l1}& {\tilde{y}}_{l2}& ...& {\tilde{y}}_{\mathrm{lN}}\end{array}\right]$

步骤S7：构建三层前馈神经网络，其输入变量集为输出变量集为隐含层神经元个数为s₁，输入层、隐含层、输出层通过权值、阈值进行连接，并且该神经网络的输入输出函数表达式为：

式中函数F(X)为S型函数。图3示出了本发明的一个实施例所使用的前馈神经网络的结构。

步骤S8：利用无迹卡尔曼滤波对所述前馈神经网络进行训练，得到该神经网络的结构参数值，具体说，该训练过程包括：

①将所述神经网络中的所有权值和阈值组成状态变量I：

$I={[w_{11}^{1}L{w}_{{\mathrm{Ms}}_1}^1,b_1^{1}L{b}_{s1}^1,w_{11}^{2}L{w}_{s_{1}l}^2,b_1^{2}L{b}_l^2]}^T$

其中，M为输入层神经元数，s₁为隐层神经元数，l为输出层神经元数，输入层至隐层神经元的连接权值为阈值为隐层至输出层的连接权值为阈值为I中的元素个数为n；设定非线性方程：

其中，函数表达式参考步骤S7，为K时刻的神经网络输入样本，令ω_k＝0，v_k＝0，为神经网络输出样本。

②设定无迹卡尔曼计算过程中控制采样点的分布状态参数a、待选参数κ，以及非负权系数β。

③计算2n+1个σ点(即采样点，一个采样点即为一组I值)以及σ点(采样点)的相应权重，其中n为状态矩阵的的I维度，λ＝a²(n+κ)-n，2n+1个采样点计算如下：

$\left\{\begin{array}{c}{I}^{\left(0\right)}=\stackrel{\&OverBar;}{I}\\ I^{\left(i\right)}=\stackrel{\&OverBar;}{I}+\sqrt{(n+\mathrm{\&lambda;})P},i=1:n\\ I^{\left(i\right)}=\stackrel{\&OverBar;}{I}-\sqrt{(n+\mathrm{\&lambda;})P},i=n+1:2n\end{array}\right.$

每个采样点的权值计算如下：

$\left\{\begin{array}{c}{W}_m^{\left(0\right)}=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}\\ w_0^{\left(0\right)}=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}+(1-a+\mathrm{\&beta;})\\ W_m^{\left(i\right)}=W_c^{\left(i\right)}=\frac{\mathrm{\&lambda;}}{2(n+\mathrm{\&lambda;})},i=1:2n\end{array}\right.$

④计算σ点的一步状态预测及状态变量协方差P_k+1|k；

$\begin{array}{c}{I}_{k+1|k}^{\left(i\right)}=I_{k|k}^{\left(i\right)}\\ {\hat{I}}_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_m^{\left(i\right)}\&CenterDot;I_{k+1|k}^{\left(i\right)}\\ P_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_c^{\left(i\right)}\&CenterDot;[I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}^{\left(i\right)}]{[I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}^{\left(i\right)}]}^T\end{array}$

⑤计算输出的一步提前预测以及协方差

$\begin{array}{c}{Y}_{k+1|k}^{\left(i\right)}=h(I_{k+1|k}^{\left(i\right)},X_k)\\ {\hat{Y}}_{k+1|k}=\stackrel{2n}{\underset{i=0}{\mathrm{\&Sigma;}}{W}_m^{\left(i\right)}\&CenterDot;}h(I_{k+1|k}^{\left(i\right)},X_k)\\ P_{Y_{k+1}}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_c^{\left(i\right)}(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k}){(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k})}^T\\ P_{X_{k+1}{Y}_{k+1}}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_c^{\left(i\right)}(I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}){(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k})}^T\end{array}$

⑥进行滤波更新获取新的状态矩阵、协方差矩阵、增益矩阵：

$R_{k+1}=P_{X_{k+1}{Y}_{k+1}}/P_{Y_{k+1}}$ (n*s2的矩阵)

${\hat{I}}_{k+1|k+1}={\hat{I}}_{k+1|k}+R_{k+1}(Y_{k+1}-{\hat{Y}}_{k+1|k})$

$P_{k+1|k+1}=P_{k+1|k}-R_{k+1}{\left(P_{X_{k+1}{Z}_{k+1}}\right)}^T$ (n×n的矩阵)

⑦对获取的新样本数据重新进行②～⑥步骤，直至所有样本对状态矩阵、协方差矩阵、增益矩阵进行了更新，从而得到适应于所有样本状态矩阵。

⑧对最后一组样本得到状态矩阵I，作为网络训练得到的权值和阈值。

⑨在得到网络参数各层权值、阈值后，确定所述前馈神经网络的函数模型为：

$\hat{Y}\left(X\right)=g^{-1}\left(\tilde{Y}\left(\tilde{X}\right)\right)=g^{-1}\left(\tilde{Y}\left(f\left(X\right)\right)\right)$

步骤S9：利用决策变量(x₁,x₂)构建父代种群P_D，

$P_D=\left\{(x_{1m}^{P_D},x_{2m}^{P_D})\right|1\&le;m\&le;K\}$

其中，父代种群P_D中的个体的数量为K，并从x₁的取值范围x_1,min≤x₁≤x_1,max内随机取值赋予从x₂的取值范围x_2,min≤x₂≤x_2,max内随机取值赋予从而对父代种群P_D进行初始化。

步骤S10：从父代种群P_D中选出任意对个体，对于每对个体 进行遗传交叉计算或变异计算，并将计算结果赋予子代种群Q_D中相应的一对个体

在一个实施例中，每对个体之间的遗传交叉计算公式为：

$\begin{array}{c}{x}_{1m}^{Q_D}=0.5\&times;(x_{1m}^{P_D}(1+\mathrm{\&alpha;})+x_{1n}^{P_D}(1-\mathrm{\&alpha;}))\\ x_{2m}^{Q_D}=0.5\&times;(x_{2m}^{P_D}(1+\mathrm{\&alpha;})+x_{2n}^{P_D}(1-\mathrm{\&alpha;}))\\ x_{1n}^{Q_D}=0.5\&times;(x_{1m}^{P_D}(1-\mathrm{\&alpha;})+x_{1n}^{P_D}(1+\mathrm{\&alpha;}))\\ x_{2n}^{Q_D}=0.5\&times;(x_{2m}^{P_D}(1-\mathrm{\&alpha;})+x_{2n}^{P_D}(1+\mathrm{\&alpha;}))\end{array}$

随机数α∈[0,1]；

每个个体的变异计算公式为：

$\begin{array}{c}{x}_{1m}^{Q_D}=x_{1m}^{P_D}(1+\mathrm{\&beta;})\\ x_{2m}^{Q_D}=x_{2m}^{P_D}(1+\mathrm{\&beta;})\end{array}$

随机数β∈[0,1]。

步骤S11：将父代种群P_D与子代种群Q_D进行合并得到种群R＝P_D∪Q_D，即有 $R=\left\{(x_{1s}^R,x_{2s}^R)\right|1\&le;s\&le;2K\}=\{(x_{1m}^{P_D},x_{2m}^{P_D})|1\&le;m\&le;K\}\&cup;\{x_{1n}^{Q_D},x_{2n}^{Q_D}|1\&le;n\&le;K\},$ 将种群R的每个个体与环境变量平均值i＝3,…,M合成输入样本 $X_s={\left[\begin{array}{ccccc}{x}_{1s}^R& x_{2s}^R& {\stackrel{\&OverBar;}{x}}_3& ...& {\stackrel{\&OverBar;}{x}}_M\end{array}\right]}^T,$ 并计算相应的适应度函数 $\mathrm{objFun}\left(X_s\right)={\left[\begin{array}{cc}h\left({\hat{y}}_1\left(X_s\right)\right)& {\hat{y}}_2\left(X_s\right)\end{array}\right]}^T,$ 其中函数h(x)＝-x，这是因为由于优化计算过程中，得到的产液量值越大，耗电量越小,其性能变量越好。

所述环境参数的平均值的计算公式为：

${\stackrel{\&OverBar;}{x}}_i=\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ik}},i=3,...,M$

步骤S12：将种群R的所有个体所对应的适应度函数相互进行比较，将种群R的所有个体划分到具有不同层级的非支配集中，其中，对于层级较低的非支配集中的任一个体所对应的适应度函数objFun(X_s)和层级较高的非支配集中的任一个体所对应的适应度函数objFun(X_t)来说，均不存在且而对于同一层级的非支配集中的任两个个体来说，该两个不等式中至少有一个不成立。

步骤S13：按照层级从低到高的顺序从所述非支配集中选择K个个体，将选择出的K个个体的值赋予父代种群P_D中的个体，返回步骤S10。

在按照层级从低到高的顺序从所述非支配集中选择K个个体时，对于同一层级非支配集中的个体，选择个体拥挤度d_s较大的个体，所述个体拥挤度d_s的计算方法为：

对当前种群R中所有个体所对应的适应度函数值objFun(X_s)中的按从小到大的顺序排序，另外，对所有objFun(X_s)中的按从小到大的顺序排序，令每次排序的第一个和最后一个个体的拥挤距离为无穷大，种群个体的拥挤度d_s为

$d_s=\frac{h{\left({\hat{y}}_1\left(X_s\right)\right)}_{+}-h{\left({\hat{y}}_1\left(X_s\right)\right)}_{-}}{h{\left({\hat{y}}_1\right)}_{\max }-h{\left({\hat{y}}_1\right)}_{\min }}+\frac{{\hat{y}}_2{\left(X_s\right)}_{+}-{\hat{y}}_2{\left(X_s\right)}_{-}}{{\hat{y}}_{2,\max }-{\hat{y}}_{2,\min }}$

分别为在所述所排的序列中值的后一个值和前一个值；

分别为在所述所排的序列中值的后一个值和前一个值；

$\begin{array}{c}h{\left({\hat{y}}_1\right)}_{\max }=\max \left\{h\left({\hat{y}}_1\left(X_s\right)\right)\right|1\&le;s\&le;2K\}\\ h{\left({\hat{y}}_1\right)}_{\min }=\min \left\{h\left({\hat{y}}_1\left(X_s\right)\right)\right|1\&le;s\&le;2K\}\\ {\hat{y}}_{2,\max }=\max \left\{{\hat{y}}_2\left(X_s\right)\right|1\&le;s\&le;2K\}\\ {\hat{y}}_{2,\min }=\min \left\{{\hat{y}}_2\left(X_s\right)\right|1\&le;s\&le;2K\}.\end{array}$

步骤S14：循环执行步骤10)-步骤13)的过程GEN次，GEN为预先确定的循环次数，最终得到优化后的K组决策变量将优化后的决策变量以及所述环境变量的平均值构成优化后的输入样本 $X_m^{P_D^{\mathrm{GEN}}}={\left[\begin{array}{ccccc}{x}_{1m}^{P_D^{\mathrm{GEN}}}& x_{2m}^{P_D^{\mathrm{GEN}}}& {\stackrel{\&OverBar;}{x}}_3& ...& {\stackrel{\&OverBar;}{x}}_M\end{array}\right]}^T(1\&le;m\&le;K),$ 这K个样本保证了在产液量最大时，耗电量最小。

下面以大港油田港510-3抽油机为实验对象采用本发明的方法进行优化。

确定抽油机采油过程生产效率影响因素构成效率观测变量集合其中α₁，α₂为决策变量冲次、有效冲程，α₃～α₁₄₆功图载荷数据，α₁₄₇为环境变量平均功率因数，α₁₄₈为环境变量平均有功功率，α₁₄₉为环境变量平均无功功率，α₁₅₀为环境变量理论排量、α₁₅₁为环境变量计算泵效、α₁₅₂为环境变量含水率。选择抽油机生产过程的产液量y₁、耗电量y₂作为建立抽油机系统模型的性能变量{y₁,y₂}。

采集油井2013年12月26日至2014年12月23日数据。对所有变量按照24小时为采集时间间隔求取24小时数据的平均值，并作为该变量样本。采集数据部分地示于表1和表2。

表1

表2

利用PCA输入变量降维处理，得到新主元数据示于表3。

表3

构建新的采集建模样本数据[X,Y]。共获得可建模数据256组，示于表4中。

表4

将匹配后的数据与载荷新主元数据进行匹配，并进行归一化处理。可建模样本归一化后部分数据示于表5中。

表5

通过无迹卡尔曼滤波对前馈神经网络的权值、阈值进行估计，将神经网络权值、阈值作为无迹卡尔曼滤波的状态变量，神经网络的输出作为无迹卡尔曼滤波的测量变量，从而得到抽油机工艺系统的精确模型。设置网络隐含层节点数为s₁＝5。设置隐层神经元输出函数为F(x)＝1/(exp(-x)+1)，得到：

输入层到隐层的权值w1：

$w1={\left[\begin{array}{ccccc}110.7717& -199.056& 20.95446& -300363& -47389.5\\ 58.1012& -53.4158& 11.35771& -94834.5& -17495.6\\ 50.60542& -42.3788& 2.027481& -215749& -32139\\ 104.5038& -84.0086& 33.51972& -110822& -19312.9\\ -70.2498& 90.35373& -5.15677& 215379.2& 50189.49\\ 90.4108& -148.601& 15.23917& -188964& -31993\\ -72.727& 100.1093& -5.96039& -190040& -31293.2\\ -95.8122& 103.793& -13.775& 206328.1& 28616.68\\ 26.40261& 12.69371& -7.30595& -173993& -16516.4\\ 98.10831& -72.4617& 5.587048& -44474.1& 15335.2\end{array}\right]}^T$

隐层神经元阈值b1：

$b1=\left[\begin{array}{c}112.442\\ -132.338\\ 34.751\\ -1.88*{10}^5\\ -3.16*{10}^4\end{array}\right]$

隐层到输出层权值w2：

$w2=\left[\begin{array}{ccccc}0.516& -0.005& 0.598& 0.772& -0.376\\ 0.640& 0.029& 0.597& 1.472& -1.426\end{array}\right]$

输出层阈值b2：

$b2=\left[\begin{array}{c}-0.365\\ -0.484\end{array}\right]$

利用上述权值和阈值，构建数学模型，选择可建模样本中最后60组数据作为测试样本进行验证。图4示出了测试样本产液量预测效果图；图5示出了测试样本耗电量预测效果图。

由模型的相对误差可知，建模效果较好，随着样本的不断训练，模型精度越来越高，符合动态建模的特性。

利用决策变量(x₁,x₂)的K对数据构建初始种群P_D，设置种群大小50，即初始化的种群样本数量为K＝50；决策变量冲次(x₁)的优化范围2.0≤x₁≤4.5；决策变量有效冲程(x₂)的范围3.0≤x₂≤3.6。设置最大遗传代数GEN＝100。

计算环境参数平均值具体算法如下：

${\stackrel{\&OverBar;}{x}}_i=\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ik}},i=3,...,M$

具体的环境参数平均值示于表6中。

表6

计算得到的产液量偏好值与耗电量的pareto解集关系示于图6中。

Pareto前沿所对应的数值解集示于表7中。

表7

在优化过程采用向最小值优化方法策略。故而，产液量值取反优化，向最小值方向计算，从而得到pareto前沿。

将优化后的决策变量，以及环境变量的平均值带入建立工艺过程模型，计算优化后的决策变量的系统性能，与实际样本的系统性能平均值进行比较，在保证固定的产液量的前提下，耗电量降低，则上述方法有效。去产液量偏好值解集的平均值求取实际对应的产液量大小。其对比结果示于表8。

表8

由表8可知，产液量得到提高，耗电量得到下降，系统效率得到了提高。达到了节能增产的目的。说明该方法有效。

尽管已经结合详细示出并描述的优选实施例公开了本发明，但是本领域技术人员应当理解，对于上述本发明所提出的基于功图主元分析的抽油机参数优化方法，还可以在不脱离本发明内容的基础上做出各种改进。因此，本发明的保护范围应当由所附的权利要求书的内容确定。

Claims (8)

1.一种基于功图主元分析的抽油机参数优化方法，包括如下步骤：

1)确定抽油机采油过程生产效率影响因素构成效率观测变量集合其中α₁，α₂为决策变量，α₃～α₁₄₆载荷数据环境变量，α₁₄₇～α_M为其他环境变量，选取抽油机系统的性能变量构成性能观测变量集合：{y₁,y₂,y₃,…y_l}；

2)获得所述生产效率影响因素和系统性能变量的样本数据，得到效率影响因素样本矩阵α和性能样本矩阵Y：

$\mathrm{\&alpha;}=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_{11}& {\mathrm{\&alpha;}}_{12}& ...& {\mathrm{\&alpha;}}_{1N}\\ {\mathrm{\&alpha;}}_{21}& {\mathrm{\&alpha;}}_{22}& ...& {\mathrm{\&alpha;}}_{2N}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}1}& {\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}2}& ...& {\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}N}\end{array}\right]\&equiv;\left[\begin{array}{c}{L}_1\\ L_2\\ ...\\ L_{\stackrel{\&OverBar;}{M}}\end{array}\right]$

$Y=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{22}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

其中为效率影响因素个数，N为样本个数，α_ik表示第i个效率影响因素变量的第k个观测值，i＝1,2,...,M；k＝1,2,...,N；

3)利用主元分析算法对载荷数据进行降维处理，从而构建新的载荷主元变量矩阵：

$\left[\begin{array}{c}{L}_{z1}\\ L_{z2}\\ L_{z3}\\ L_{\mathrm{zd}}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_{z11}& {\mathrm{\&alpha;}}_{z12}& ...& {\mathrm{\&alpha;}}_{z1N}\\ {\mathrm{\&alpha;}}_{z21}& {\mathrm{\&alpha;}}_{z22}& ...& {\mathrm{\&alpha;}}_{z2N}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{\mathrm{zd}1}& {\mathrm{\&alpha;}}_{\mathrm{zd}2}& ...& {\mathrm{\&alpha;}}_{\mathrm{zdN}}\end{array}\right]$

4)由影响因素观测变量集合中非载荷变量与载荷新主元观测变量集合{α_z1,α_z2,...,α_zd}构建网络输入变量集合：并令输入变量集合为：{x₁,x₂,x₃,...,x_M}，即， $\{{\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,{\mathrm{\&alpha;}}_{147},...,{\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}},{\mathrm{\&alpha;}}_{z1},...,{\mathrm{\&alpha;}}_{\mathrm{zd}}\}=\{x_1,x_2,x_3,...,x_M\};$

5)构建输入变量集合{x₁,x₂,x₃,...,x_M}观测样本值：

$X=\left[\begin{array}{cccc}{X}_1& X_2& ...& X_N\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_{11}& {\mathrm{\&alpha;}}_{12}& ...& {\mathrm{\&alpha;}}_{1N}\\ {\mathrm{\&alpha;}}_{21}& {\mathrm{\&alpha;}}_{22}& ...& {\mathrm{\&alpha;}}_{2N}\\ {\mathrm{\&alpha;}}_{\mathrm{147,1}}& {\mathrm{\&alpha;}}_{\mathrm{147,2}}& ...& {\mathrm{\&alpha;}}_{147,N}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}1}& {\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}2}& ...& {\mathrm{\&alpha;}}_{\stackrel{\&OverBar;}{M}N}\\ {\mathrm{\&alpha;}}_{z11}& {\mathrm{\&alpha;}}_{z12}& ...& {\mathrm{\&alpha;}}_{z1N}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{\mathrm{zd}1}& {\mathrm{\&alpha;}}_{\mathrm{zd}2}& ...& {\mathrm{\&alpha;}}_{\mathrm{zdN}}\end{array}\right]\begin{array}{}\end{array}=\left[\begin{array}{cccc}{x}_{11}& x_{12}& ...& x_{1N}\\ x_{21}& x_{22}& ...& x_{2N}\\ x_{31}& x_{32}& ...& x_{3N}\\ ...& ...& ...& ...\\ x_{M1}& x_{M2}& ...& x_{\mathrm{MN}}\end{array}\right]$

$Y=\left[\begin{array}{cccc}{Y}_1& Y_2& ...& T_N\end{array}\right]=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{33}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

其中，x₁～x₂为决策变量，x₃～x_M为新的环境变量；

6)对得到的训练输入样本X、输出样本Y进行归一化处理，得到新的训练输入矩阵输出矩阵

${\tilde{X}}_k={\left[\begin{array}{ccccc}{\tilde{x}}_{1k}& {\tilde{x}}_{2k}& {\tilde{x}}_{3k}& ...& {\tilde{x}}_{\mathrm{Mk}}\end{array}\right]}^T={\left[\begin{array}{ccccc}f\left(x_{1k}\right)& f\left(x_{2k}\right)& f\left(x_{3k}\right)& ...& f\left(x_{\mathrm{Mk}}\right)\end{array}\right]}^T\&equiv;f\left(X_k\right)$

${\tilde{Y}}_k={\left[\begin{array}{cccc}{\tilde{y}}_{1k}& {\tilde{y}}_{2k}& ...& {\tilde{y}}_{\mathrm{lk}}\end{array}\right]}^T{\left[\begin{array}{cccc}g\left(y_{1k}\right)& g\left(y_{2k}\right)& ...& g\left(y_{\mathrm{lk}}\right)\end{array}\right]}^T\&equiv;g\left(Y_k\right)$

$\tilde{X}=[{\tilde{X}}_1,{\tilde{X}}_2,...,{\tilde{X}}_N]=\left[\begin{array}{cccc}{\tilde{x}}_{11}& {\tilde{x}}_{12}& ...& {\tilde{x}}_{1N}\\ {\tilde{x}}_{21}& {\tilde{x}}_{22}& ...& {\tilde{x}}_{2N}\\ ...& ...& ...& ...\\ {\tilde{x}}_{M1}& {\tilde{x}}_{M2}& ...& {\tilde{x}}_{\mathrm{MN}}\end{array}\right]$

$\tilde{Y}=\left[\begin{array}{cccc}{\tilde{Y}}_1& {\tilde{Y}}_2& ...& {\tilde{Y}}_N\end{array}\right]=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{22}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

7)构建三层前馈神经网络，其输入变量集为输出变量集为隐含层神经元个数为s₁，输入层、隐含层、输出层通过权值、阈值进行连接，并且该神经网络的输入输出函数表达式为：

式中函数F(X)为S型函数；

8)利用无迹卡尔曼滤波对所述前馈神经网络进行训练，得到该神经网络的结构参数值，该训练过程包括：

①将所述神经网络中的所有权值和阈值组成状态变量I：

$I={\left[\begin{array}{ccccc}{w}_{11}^{1}L& w_{{\mathrm{Ms}}_1}^1,b_1^{1}L& b_{s1}^1,w_{11}^{2}L& w_{s_{1}l}^2,b_1^{2}L& b_l^2\end{array}\right]}^T$

其中，M为输入层神经元数，s₁为隐层神经元数，l为输出层神经元数，输入层至隐层神经元的连接权值为 $w_{\mathrm{ik}}^1(i=\mathrm{0,1},...,\mathrm{Mkk}=\mathrm{1,2},...,s_1),$ 阈值为 $b_k^1(k=\mathrm{1,2},...,s_1),$ 隐层至输出层的连接权值为 $w_{\mathrm{kj}}^2(i=\mathrm{0,1},L,s_1;j=\mathrm{1,2},L,l),$ 阈值为 $b_j^2=(j=\mathrm{1,2},L,l),$ I中的元素个数为n；设定非线性方程：

其中，函数表达式参考步骤7)，为K时刻的神经网络输入样本，令ω_k＝0，v_k＝0，为神经网络输出样本；

②设定无迹卡尔曼计算过程中控制采样点的分布状态参数a、待选参数κ，以及非负权系数β；

③计算2n+1个采样点(σ点)以及σ点的相应权重，其中n为状态矩阵的I维度，λ＝a²(n+κ)-n，2n+1个采样点计算如下：

$\left\{\begin{array}{cc}{I}^{\left(0\right)}=\stackrel{\&OverBar;}{I}& \\ I^{\left(i\right)}=\stackrel{\&OverBar;}{I}+\sqrt{(n+\mathrm{\&lambda;})P}& i=1:n\\ I^{\left(i\right)}=\stackrel{\&OverBar;}{I}-\sqrt{(n+\mathrm{\&lambda;})P}& i=n+1:2n\end{array}\right.$

每个采样点的权值计算如下：

$\left\{\begin{array}{c}{W}_m^{\left(0\right)}=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}\\ W_c^{\left(0\right)}=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}+(1-a+\mathrm{\&beta;})\\ W_m^{\left(i\right)}=W_c^{\left(i\right)}=\frac{\mathrm{\&lambda;}}{2(n+\mathrm{\&lambda;})}i=1:2n\end{array}\right.$

④计算σ点的一步状态预测及状态变量协方差P_k+1|k；

$I_{k+1|k}^{\left(i\right)}=I_{k|K}^{\left(i\right)}$

$\begin{array}{c}{\hat{I}}_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_m^{\left(i\right)}\&CenterDot;I_{k+1|k}^{\left(i\right)}\\ P_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_c^{\left(i\right)}\&CenterDot;[I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}^{\left(i\right)}]{[I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}^{\left(i\right)}]}^T\end{array}$

⑤计算输出的一步提前预测以及协方差

$Y_{k+1|k}^{\left(i\right)}=h(I_{k+1|k}^{\left(i\right)},X_k)$

${\hat{Y}}_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_m^{\left(i\right)}\&CenterDot;h(I_{k+1|k}^{\left(i\right)},X_k)$

$P_{Y_{k+1}}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{w}_c^{\left(i\right)}(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|j}){(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k})}^T$

$P_{X_{k+1}{Y}_{k+1}}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_c^{\left(i\right)}(I_{k+1|k}^{\left(i\right)}-{\hat{I}}_{k+1|k}){(Y_{k+1|k}^{\left(i\right)}-{\hat{Y}}_{k+1|k})}^T$

⑥进行滤波更新获取新的状态矩阵、协方差矩阵、增益矩阵：

$R_{k+1}=P_{X_{k+1}{Y}_{k+1}}/P_{Y_{k+1}}$

${\hat{I}}_{k+1|k+1}={\hat{I}}_{k+1|k}+R_{k+1}(Y_{k+1}-{\hat{Y}}_{k+1|k})$

$P_{k+1|k+1}=P_{k+1|k}-R_{k+1}{\left(P_{X_{k+1}{Z}_{k+1}}\right)}^T$

⑦对获取的新样本数据重新进行②～⑥步骤，直至所有样本对状态矩阵、协方差矩阵、增益矩阵进行了更新，从而得到适应于所有样本状态矩阵；

⑧对最后一组样本得到状态矩阵I，作为网络训练得到的权值和阈值；

⑨在得到网络参数各层权值、阈值后，确定所述前馈神经网络的函数模型为：

$\hat{Y}\left(X\right)=g^{-1}\left(\tilde{Y}\left(\tilde{X}\right)\right)=g^{-1}\left(\tilde{Y}\left(f\left(X\right)\right)\right)$

9)利用决策变量(x₁,x₂)构建多目标优化父代种群P_D

$P_D=\left\{(x_{1m}^{P_D},x_{2m}^{P_D})\right|1\&le;m\&le;K\}$

其中，父代种群P_D中的个体的数量为K，并从x₁的取值范围x_1,min≤x₁≤x_1,max内随机取值赋予从x₂的取值范围x_2,min≤x₂≤x_2,max内随机取值赋予从而对父代种群P_D进行初始化；

10)从父代种群P_D中选出任意对个体，对于每对个体 进行遗传交叉计算或变异计算，并将计算结果赋予子代种群Q_D中相应的一对个体

11)将父代种群P_D与子代种群Q_D进行合并得到种群R＝P_D∪Q_D，即有 $R=\left\{(x_{1s}^R,x_{2s}^R)\right|1\&le;s\&le;2K\}=\{(x_{1m}^{P_D},x_{2m}^{P_D})|1\&le;m\&le;K\}\&cup;\left\{(x_{1n}^{Q_D},x_{2n}^{Q_D})\right|1\&le;n\&le;K\},$ 将种群R的每个个体与环境变量平均值i＝3,…,M合成输入样本 $X_s={\left[\begin{array}{ccccc}{x}_{1s}^R& x_{2s}^R& {\stackrel{\&OverBar;}{x}}_3& ...& {\stackrel{\&OverBar;}{x}}_M\end{array}\right]}^T,$ 并计算相应的适应度函数 $\mathrm{objFun}\left(X_s\right)={\left[\begin{array}{cc}h\left({\hat{y}}_1\left(X_s\right)\right)& {\hat{y}}_2\left(X_s\right)\end{array}\right]}^T,$ 其中，函数h(x)＝-x；

12)将种群R的所有个体所对应的适应度函数相互进行比较，将种群R的所有个体划分到具有不同层级的非支配集中，其中，对于层级较低的非支配集中的任一个体所对应的适应度函数objFun(X_s)和层级较高的非支配集中的任一个体所对应的适应度函数objFun(X_t)来说，均不存在且而对于同一层级的非支配集中的任两个个体来说，该两个不等式中至少有一个不成立；

13)按照层级从低到高的顺序从所述非支配集中选择K个个体，将选择出的K个个体的值赋予父代种群P_D中的个体，并执行步骤10)-步骤13)的过程GEN次，GEN为预先确定的循环次数，最终得到优化后的K组决策变量将优化后的决策变量以及所述环境变量的平均值构成优化后的输入样本 $X_m^{P_D^{\mathrm{GEH}}}={\left[\begin{array}{ccccc}{x}_{1m}^{P_D^{\mathrm{GEN}}}& x_{2m}^{P_D^{\mathrm{GEN}}}& {\stackrel{\&OverBar;}{x}}_3& ...& {\stackrel{\&OverBar;}{x}}_M\end{array}\right]}^T,(1\&le;m\&le;K),$ 这K个样本保证了在产液量最大，耗电量最小。

2.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

所述决策变量α₁为冲次、决策变量α₂为有效冲程、α₃～α₁₄₆为载荷1至载荷144，其余环境变量包括：理论排量、功率因数、有功功率、无功功率、含水率中的一个或多个变量；所述抽油机生产过程性能变量y₁为产液量、y₂为耗电量。

3.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

设所述决策变量、环境变量和性能变量的观测值采集周期的最大值为tmax，则这些变量中的任一变量的样本取为tmax时间内该变量的观测值的平均值。

4.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

利用主元分析算法对载荷数据进行降维处理的步骤包括：

①设置样本累计贡献率precent＝0.95；

②获取载荷数据每个L_k具有第k观测变量的N个观测数据,3≤k≤146；

③求出数据平均值并利用原始数据减去均值得到

④计算协方差矩阵 $\mathrm{Cov}=\frac{1}{144}\underset{i=1}{\overset{144}{\mathrm{\&Sigma;}}}(L_k^{\&prime;T}\&CenterDot;L_k^{\&prime;});$

⑤计算协方差矩阵的特征值E₁,E₂,...,E₁₄₄与特征向量EV₁,EV₂,...,EV₁₄₄；

⑥由大到小依次排列特征值E′₁,E'₂,...,E'_M,对应特征向量为EV′₁,EV′₂,...,EV′₁₄₄，按特征值大小顺序取前d个特征值的特征向量构成矩阵[EV₁',EV₂',...,EV_d']，此时 $\underset{i=1}{\overset{d}{\mathrm{\&Sigma;}}}{E}_i^{\&prime;}/\underset{i=1}{\overset{144}{\mathrm{\&Sigma;}}}{E}_i^{\&prime;}\&GreaterEqual;0.95,$ 其中d<144；

⑦由[EV′₁,EV′₂,...,EV′_d]与原始样本求取载荷新的主元，其新载荷主元观测变量构成集合:{α_z1,α_z2,...,α_zd}，其为_d个新变量，且每个变量为N个观测值构成的新主元矩阵：

$\left[\begin{array}{c}{L}_{z1}\\ L_{z2}\\ L_{z3}\\ L_{\mathrm{zd}}\end{array}\right]=\begin{array}{c}\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_{z11}& {\mathrm{\&alpha;}}_{z12}& ...& {\mathrm{\&alpha;}}_{z1N}\\ {\mathrm{\&alpha;}}_{z21}& {\mathrm{\&alpha;}}_{z22}& ...& {\mathrm{\&alpha;}}_{z2N}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{\mathrm{zd}1}& {\mathrm{\&alpha;}}_{\mathrm{zd}2}& ...& {\mathrm{\&alpha;}}_{\mathrm{zdN}}\end{array}\right]\end{array}$

5.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

所述归一化处理的算法如下：

${\tilde{x}}_{\mathrm{ik}}=f\left(x_{\mathrm{ik}}\right)=({\tilde{x}}_{\max }-{\tilde{x}}_{\min })\&CenterDot;\frac{x_{\mathrm{ik}}-x_{i,\min }}{x_{i,\max }-x_{i,\min }}+{\tilde{x}}_{\min }$

i＝1,2,...,M；k＝1,2,...,N

${\tilde{y}}_{\mathrm{jk}}=f\left(y_{\mathrm{jk}}\right)=({\tilde{y}}_{\max }-{\tilde{y}}_{\min })\&CenterDot;\frac{y_{\mathrm{jk}}-y_{j,\min }}{y_{j,\max }-y_{j,\min }}+{\tilde{y}}_{\min }$

j＝1,2,....,l；k＝1,2,...,N

其中:为设定输入变量归一化后数据范围的最大值、最小值；

x_ik为归一化前的第i个输入变量第k个样本值；

为归一化后第i个输入变量第k个样本值；

x_i,min＝min{x_ik|1≤k≤N}

x_i,max＝max{x_ik|1≤k≤N}

为设定输出变量归一化后数据范围的最大值、最小值；

y_jk为归一化前第j个输出变量的第k个采集样本值；

为归一化后第j个输出变量的第k个值；

y_j,max＝max{y_jk|1≤k≤N}

y_j,min＝min{y_jk|1≤k≤N}

于是得到：

${\tilde{X}}_k={\left[\begin{array}{ccccc}{\tilde{x}}_{1k}& {\tilde{x}}_{2k}& {\tilde{x}}_{3k}& ...& {\tilde{x}}_{\mathrm{Mk}}\end{array}\right]}^T={\left[\begin{array}{ccccc}f\left(x_{1k}\right)& f\left(x_{2k}\right)& f\left(x_{3k}\right)& ...& f\left(x_{\mathrm{Mk}}\right)\end{array}\right]}^T\&equiv;f\left(X_k\right)$

${\tilde{Y}}_k={\left[\begin{array}{cccc}{\tilde{y}}_{1k}& {\tilde{y}}_{2k}& ...& {\tilde{y}}_{\mathrm{lk}}\end{array}\right]}^T{\left[\begin{array}{cccc}g\left(y_{1k}\right)& g\left(y_{2k}\right)& ...& g\left(y_{\mathrm{lk}}\right)\end{array}\right]}^T\&equiv;g\left(Y_k\right)$

$\tilde{X}=[{\tilde{X}}_1,{\tilde{X}}_2,...,{\tilde{X}}_N]=\left[\begin{array}{cccc}{\tilde{x}}_{11}& {\tilde{x}}_{12}& ...& {\tilde{x}}_{1N}\\ {\tilde{x}}_{21}& {\tilde{x}}_{22}& ...& {\tilde{x}}_{2N}\\ ...& ...& ...& ...\\ {\tilde{x}}_{M1}& {\tilde{x}}_{M2}& ...& {\tilde{x}}_{\mathrm{MN}}\end{array}\right]$

$\tilde{Y}=\left[\begin{array}{cccc}{\tilde{Y}}_1& {\tilde{Y}}_2& ...& {\tilde{Y}}_N\end{array}\right]=\left[\begin{array}{cccc}{y}_{11}& y_{12}& ...& y_{1N}\\ y_{21}& y_{22}& ...& y_{2N}\\ ...& ...& ...& ...\\ y_{l1}& y_{l2}& ...& y_{\mathrm{lN}}\end{array}\right]$

6.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

多目标优化的每对个体之间的遗传交叉计算公式为：

$x_{1m}^{Q_D}=0.5\&times;(x_{1m}^{P_D}(1+\mathrm{\&alpha;})+x_{1m}^{P_D}(1-\mathrm{\&alpha;}))$

$x_{2m}^{Q_D}=0.5\&times;(x_{2m}^{P_D}(1+\mathrm{\&alpha;})+x_{2n}^{P_D}(1-\mathrm{\&alpha;}))$

$x_{1n}^{Q_D}=0.5\&times;(x_{1n}^{P_D}(1-\mathrm{\&alpha;})+x_{1n}^{P_D}(1+\mathrm{\&alpha;}))$

$x_{2n}^{Q_D}=0.5\&times;(x_{2m}^{P_D}(1-\mathrm{\&alpha;})+x_{2n}^{P_D}(1+\mathrm{\&alpha;}))$

随机数α∈[0,1]；

每个个体的变异计算公式为：

$x_{1m}^{Q_D}=x_{1m}^{P_D}(1+\mathrm{\&beta;})$

$x_{2m}^{Q_D}=x_{2m}^{P_D}(1+\mathrm{\&beta;})$

随机数β∈[0,1]。

7.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

所述环境参数的平均值的计算公式为：

${\stackrel{\&OverBar;}{x}}_i=\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ik}},i=3,...,M$

8.如权利要求1所述的基于功图主元分析的抽油机参数优化方法，其中，

在按照层级从低到高的顺序从所述非支配集中选择K个个体时，对于同一层级非支配集中的个体，选择个体拥挤度d_s较大的个体，所述个体拥挤度d_s的计算方法为：

对当前种群R中所有个体所对应的适应度函数值objFun(X_s)中的按从小到大的顺序排序，另外，对所有objFun(X_s)中的按从小到大的顺序排序，令每次排序的第一个和最后一个个体的拥挤距离为无穷大，种群个体 $R_s=(x_{1s}^R,x_{2s}^R)$ 的拥挤度d_s为

$d_s=\frac{h{\left({\hat{y}}_1\left(X_s\right)\right)}_{+}-h{\left({\hat{y}}_1\left(X_s\right)\right)}_{-}}{h{\left({\hat{y}}_1\right)}_{\max }-h{\left({\hat{y}}_1\right)}_{\min }}+\frac{{\hat{y}}_2{\left(X_s\right)}_{+}-{\hat{y}}_2{\left(X_s\right)}_{-}}{{\hat{y}}_{2,\max }-{\hat{y}}_{2,\min }}$

分别为在所述所排的序列中值的后一个值和前一个值；

分别为在所述所排的序列中值的后一个值和前一个值；

$h{\left({\hat{y}}_1\right)}_{\max }=\max \left\{h\left({\hat{y}}_1\left(X_s\right)\right)\right|1\&le;s\&le;2K\}$

$h{\left({\hat{y}}_1\right)}_{\min }=\min \left\{h\left({\hat{y}}_1\left(X_s\right)\right)\right|1\&le;s\&le;2K\}$

${\hat{y}}_{2,\max }=\max \left\{{\hat{y}}_2\left(X_s\right)\right|1\&le;s\&le;2K\}$

${\hat{y}}_{2,\min }=\min \left\{{\hat{y}}_2\left(X_s\right)\right|1\&le;s\&le;2K\}.$