# Electric-Parameter-Based Inversion of Dynamometer Card Using Hybrid Modeling for Beam Pumping System

Authors: Jingliang Wei; Xianwen Gao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Conventional fault diagnosis and production calculation of an oil well can be conducted with the surface dynamometer cards, which are obtained by load sensor installed on the horse head. This method to measure the dynamometer cards is limited by the sensor maintenance and calibration, battery replacement, and safety hazards for staff. As the basic parameter of the oil extraction industry, electric parameters have the advantages of low cost and high efficiency. So the inversions of dynamometer card with electric parameters are attracting more and more attention. In order to solve the problem of insufficient data and consider the real-time performance in the actual oil extraction process, this paper proposes a novel hybrid model which consists of two parts: the mechanism model of polished rod load and the suspension displacement calculated with the space vector equations of motor and a data-dependent kernel online sequential extreme learning machine (DDKOS-ELM) model proposed to correct the output error of the mechanism model, which improves the kernel function selection and makes it real-time. Thus, the highlights of this paper can be summed up in two points:(1) under the circumstance of the bottom dead point detection without sensors, the mechanism modeling of the polished rod load and suspension displacement has been carried out from the perspective of mathematical model of AC motor; (2) a novel data-driven model based on data-dependent kernel online sequential extreme learning machine (DDKOS-ELM) has been proposed to improve the kernel functions selection. The coefficients in the data-dependent kernel function are optimized with improved free search algorithm (IFSA). The proposed hybrid model has been applied to a normal working oil well and the prediction results show better accuracy than the pure data-driven model and mechanism model.


---

## Body

## 1. Introduction

The beam pumping unit plays a predominant part in oil extraction, which is the major artificial lift method and employed by virtually more than 80% oil wells all over the world [1]. As one of the key data to identify the various downhole operation conditions, dynamometer card puts forward an immense influence on the oil well fault diagnosis. It describes the relationship between displacement and load of the polished rod in each stroke and fully reflects the oil exploitation status of the pumping unit. Accordingly, the real-time measurement of dynamometer card has invariably been the major topic of discussion in oil well diagnosis. The commonly used methods of measuring the dynamometer card in actual production are realized by directly mounting load sensor. The existing detection methods of dynamometer card have the following defects:( 1 ) The installation of the measuring equipment of the dynamometer card requires the unloading of the suspension point, the operation is cumbersome, and there are safety risks.( 2 ) Generally, the overall dynamometer card measuring instrument uses battery power which requires large amount of daily maintenance. Because the production will be affected when the well is forced to stop working by the replacement of the battery and the installation of the load sensor, so the measured intervals are usually relatively long, and the downhole conditions cannot be reflected in real time. Therefore, the fault diagnosis of the pumping well may not be timely.( 3 ) The load of the suspension point is measured by the pressure strain gauge which must be calibrated periodically. Besides, the ground dynamometer card tester mainly reflects the underground working condition. It is difficult to fully reflect many important operation indicators such as the balance of pumping unit, the system efficiency, energy consumption, and motor running status. This leads to a single measurement function and low cost performance.As the significant reference variables which are the direct and real-time reflection of the working condition of the beam pumping system, the induction motor electronic parameters are often ignored by the diagnostic staff. There are some calculation methods for the transformation from motor parameters to dynamometer card in recent years. Li et al. [2] presented a method for drawing the dynamometer card based on the actual input power curve of the motor, but it needs to estimate the location of the bottom dead point and thus decreases the accuracy of the calculation. Qi [3] realized the indirect measurement of the dynamometer card of the pumping unit by combining the rotating speed method with the mechanical characteristic curve of the pumping unit, but the modified engineering formula is used to calculate the electromagnetic torque, which affects the accuracy of the method. Besides, the mechanical characteristic curve of the pumping unit is obtained from the load sensor installed on the polished rod. Zhang et al. [4] detected the dead point position of polished rod by installing two microswitches near the link of the pumping beam and bracket and then got the corresponding relation between the power of the motor and the displacement of the polished rod. Because of the mechanical friction between the beam and the microswitches, the accuracy will also be affected. Guo and Ma [5] established a mathematical model containing useful motor power, polished rod load, and polished rod displacement, which can figure out the relationship between the beam pumping unit output and the load. Luo et al. [6] developed a multiple factor nonlinear mathematical model to analyze the operating performance of the induction motor system and introduced the field-circuit coupling equation of induction motors and the mechanical model of the four-bar linkage system. The mechanism models above all have some uncertain environment factors that are hard to precisely describe because of the real-time changes with undulations of the pumping units’ operating conditions. Consequently, the data-driven methods have captured increasingly researchers’ attention. Li and Han [7] proposed a serial hybrid model of motor load torque by establishing a data-driven model of underground friction. Considering the torque factor is zero when the polished rod is at the top and bottom dead centers, Zhang et al. [8] proposed a power-dynamometer card inversion model based on big data technologies by collecting the measured power curves and corresponding dynamometer cards of typical wells under different working conditions within a period of time as the training data for deep belief network. However, in the case of insufficient data, the data-drive model is difficult to achieve excellent performance. After the above discussion, we adopt the hybrid modeling by combining the mechanism model with the data-driven model in order to solve the problem of dynamometer card inversion. Nowadays, hybrid modeling algorithms have been successfully put into use in a multitude of complex industry processes [9–22], and they mostly consist of serial structure and parallel structure. The serial hybrid model firstly performs data-driven modeling for the uncertain parameters, and then the outputs of data-driven model are taken as the input of the mechanism model, while the outputs of the parallel hybrid model are the sum of the outputs of data-driven model and mechanism model, and the data-driven model is built to compensate the error of the mechanism model. The latter is adopted as the method in this paper.In consideration and summary of above discussion, our main contributions in this paper can be summed up as follows. For the part of the mechanism model, the relationship among voltage, current, and torque in the stationary second rectangular coordinate system was derived based on the mathematical model of AC induction motor. The instantaneous torque of AC induction motor was identified according to the stator line voltage and the line current of the motor, and then the polished rod load was calculated from the motor torque in consideration of the crank weight. Finally, the corresponding relation between the polished rod load and suspension displacement is obtained from the signal of the crank rotation angle. For the part of the data-driven model, the inputs of it are the mechanism and electric motor parameters including stroke, the crank radius, stroke frequency, pulley diameter, unbalanced weight, crank quality, balance block weight, two line voltages and line currents of the stator, motor speed, and other technical parameters, and the outputs are the error estimation between the actual measured values and the outputs of the mechanism model.


## 2. Process Description

The typical structure diagram of beam pumping units is shown as Figure1, according to the working principle of oil extraction, the oil pump barrel is placed under the oil surface through tubing, with a standing valve installed on the bottom of it; piston is placed into the oil pump barrel through tubing from rod and with a travelling valve installed on the top of it. The sucker rod is connected to a polished rod suspended from the “horse head”. Under the influence of the crankshaft and four connecting rod of the pumping unit, the motor's rotational motion is turned into a reciprocating movement up and down, and then the pump is driven to extract oil. The load of the beam hanger is called polished rod load, which reflects the load of sucker rod column and the liquid column above the piston. In the upper stroke process, the sucker rod is pulled up to pull the plunger, the travelling valve is closed, the standing valve is opened, the oil is pumped into the pump, and the polished rod load is equal to the sum of the weight of sucker rod and the weight of liquid above the plunger, while in the lower stroke process, the sucker rod pushes down the plunger, the standing valve is closed first, then the travelling valve is pushed open under the pressure, and the pump discharges the oil.Figure 1
Typical structure diagram of the beam pumping units.In a stroke, the closed curve mapped out by the load and displacement of the suspension point is known as the dynamometer card, and the area of the curve is the size of the work done by the pump in a single stroke.When only considering the liquid column static load above the rod string and the plunger surface and the elastic deformation of sucker rod and tubing, regardless of the sand, wax, gas, and interference of other environmental factors, the curve is obtained as the theoretical dynamometer card, which is similar to parallelogram. As shown in Figure2, line AB is the loading section, which starts at the moment when the plunger reaches the bottom dead point after ascending, while the travelling valve is closed, the plunger pushes the liquid column above it move upward, and the weight is loaded on the polished rod through the sucker rod, elastic deformation of tubing and sucker rod occurs, the tubing is shortened, and the sucker rod elongates. During this period, although the polished rod is moving upward, the plunger relative to the pump barrel does not produce relative changes.Figure 2
Theoretical dynamometer card.When the load on the polished rod increases enough to pull the plunger, the elastic deformation of tubing and sucker rod ends, the high load BC section starts drawing, at this time the load of the polished rod stays the same, the plunger is pulled upward, standing valve is opened, and the pumped oil takes up the pump space until the top dead point.In unloading section CD, the plunger reaches the top dead point and starts to go down, the travelling valve is opened and the standing valve is closed, elastic deformation of tubing and sucker rod occurs again, the tubing is elongated, and the sucker rod is shortened. During this time, although the polished rod is moving downward, the plunger relative to the pump barrel does not produce relative changes either.When the load on the polished rod is reduced to the same as the weight of the plunger and the polished rod, the elastic deformation of tubing and sucker rod ends, the low load DA section starts drawing, at this time the load of the polished rod remains the same, the plunger and the polished rod move downward in a constant speed, travelling valve is opened, and the pumped oil flowed into the plunger from the pump barrel.It is not hard to conclude from the oil pumping process introduced above that different shapes of the dynamometer cards may cause corresponding changes because of the combined variation of gas, liquid, and machinery.


## 3. Mechanism Modeling of Beam Pumping Units

### 3.1. Pumping Motor

In the actual oil production field, the induction motor is employed as the main equipment for the oil pumping motor, and the space vector equation of ac motor can be described as follows [23]:(1)us=Rsis+dψsdtur=Rrir+dψrdt-jωrψrψs=Lsis+Lmirψr=Lmis+Lrirwhere us is the stator voltage space vector, ur is the rotor voltage space vector, is is the stator current space vector, ir is the rotor current space vector, ψs is the stator magnetic chain space vector, ψr is the rotor magnetic chain space vector, Rs is the resistance of motor stator, Rr is the resistance of motor rotor, Ls is the stator self-induction, Lr is the rotor self-induction, Lm is the mutual inductance of rotor to stator winding, and ωr is the angular velocity of rotor.Then theα axis of two phase static coordinates (α-β coordinate system) is overlapped to the A axis of three-phase static coordinate system (ABC coordinate system) and the β axis of two phase static coordinates lag α axis 90°, and then the electromagnetic torque, synchronous angular velocity, and electromagnetic power in the α-β coordinate system can be obtained by coordinate transformation as follows:(2)Te=pnψ→s×i=pnψsαisβ-ψsβisα(3)ωs=usβ-Rsisβψsα-usα-Rsisαψsβψs2(4)Pe=TeΩs=Teωspnwhere pn is the pole logarithm, ωs is the synchronous motor angular velocity, Ωs is the mechanical angular velocity, Te is the electromagnetic torque, and Pe is the electromagnetic power. The line voltage and line current component of the α-β coordinate system are obtained by the phase transformation of the stator line voltage and line current directly measured in the electrical parameter acquisition process.Considering the mechanical loss of the motor, the output torque and output power of the motor can be obtained from the following formulas:(5)ω=RrTepnψs2(6)Tm=1-sPe-p0ωs-ω/pnwhere ω is the angular velocity of the rotor slip, Tm is the output torque of motor, s is the slip ratio, and p0 is the mechanical loss of the motor.


### 3.2. The Movement Rule of Pumping Units

The surface structure of the beam pumping unit consists of drive motor, reduction gearbox, and four-bar linkage. For the four-bar linkage mechanism, the fixed rod is the connection between the beam support point and the crankshaft center, and the three movable rods are crank, connecting rod, and backward beam. As shown in Figure3, R is the radius of crank, Lc is the length of the connection rod, Lb is the length of backward beam, Lf is the length of forward beam, K is the length of the fixed rod, Ω denotes the mechanical angular velocity, O and O1 are reference points, and H is the vertical height of the support point to the center of the gearbox. Setting the clockwise direction positive, the motor parameters are collected at 12 o’clock of crank, and θo denotes the angle of crank spinning clockwise.Figure 3
Four-bar linkage mechanism.The angle between the backward beam and connecting rod can be calculated by applying the cosine theorem to triangleΔLbLcJ and ΔKRJ:(7)β=cos-1⁡Lb2+Lc2-K2-R2+2KRcos⁡θ2LbLcThe angle between the crank and connecting rod can be calculated as follows:(8)α=360∘-β+φ+θThe angles ofφ when the horse head is at the top and bottom dead points are obtained as follows, respectively:(9)φbottom=cos-1⁡Lb2+K2-Lc+R22LbLc(10)φtop=cos-1⁡Lb2+K2-Lc-R22LbLcCounting from the bottom dead point (12 o’clock of crank), the motor parameters are collected every 40ms,φ can be calculated by multiplying time by the motor speed, and then the displacement is expressed as a percentage of the stroke:(11)S=φbottom-φφbottom-φtopSCwhereS is the displacement of suspension point and SC is the stroke of the sucker rod pumping system.When the beam pumping unit is working, the torque synthesized on the crankshaft by the suspension point load and the balance weight is in equilibrium with the output torque of motor.(12)P=Tm×ibelt×ireductiongearbox×ηbelt×ηreductiongearboxσ-Wcbr+Wcrcsin⁡θ0TF¯,σ=1Tm>0-1Tm<0where TF¯=Lf/Lbrsin⁡α/sin⁡β is the torque factor, ibelt and ηbelt are transmission ratio and transmission efficiency of the belt, ireductiongearbox and ηreductiongearbox are transmission ratio and transmission efficiency of the reduction gearbox, Wcb and Wc are the weight of the crank balance block and the crank, rc is the crank radius, and P is the load of the suspension point.Thus, the relation between the displacement and the crank angleθ0 is obtained from (11), and the load of the suspension point can be calculated by (12). Then the curve of dynamometer card can be obtained with variable substitution according to the displacement and load curve of the crank angle θ0.


## 3.1. Pumping Motor

In the actual oil production field, the induction motor is employed as the main equipment for the oil pumping motor, and the space vector equation of ac motor can be described as follows [23]:(1)us=Rsis+dψsdtur=Rrir+dψrdt-jωrψrψs=Lsis+Lmirψr=Lmis+Lrirwhere us is the stator voltage space vector, ur is the rotor voltage space vector, is is the stator current space vector, ir is the rotor current space vector, ψs is the stator magnetic chain space vector, ψr is the rotor magnetic chain space vector, Rs is the resistance of motor stator, Rr is the resistance of motor rotor, Ls is the stator self-induction, Lr is the rotor self-induction, Lm is the mutual inductance of rotor to stator winding, and ωr is the angular velocity of rotor.Then theα axis of two phase static coordinates (α-β coordinate system) is overlapped to the A axis of three-phase static coordinate system (ABC coordinate system) and the β axis of two phase static coordinates lag α axis 90°, and then the electromagnetic torque, synchronous angular velocity, and electromagnetic power in the α-β coordinate system can be obtained by coordinate transformation as follows:(2)Te=pnψ→s×i=pnψsαisβ-ψsβisα(3)ωs=usβ-Rsisβψsα-usα-Rsisαψsβψs2(4)Pe=TeΩs=Teωspnwhere pn is the pole logarithm, ωs is the synchronous motor angular velocity, Ωs is the mechanical angular velocity, Te is the electromagnetic torque, and Pe is the electromagnetic power. The line voltage and line current component of the α-β coordinate system are obtained by the phase transformation of the stator line voltage and line current directly measured in the electrical parameter acquisition process.Considering the mechanical loss of the motor, the output torque and output power of the motor can be obtained from the following formulas:(5)ω=RrTepnψs2(6)Tm=1-sPe-p0ωs-ω/pnwhere ω is the angular velocity of the rotor slip, Tm is the output torque of motor, s is the slip ratio, and p0 is the mechanical loss of the motor.


## 3.2. The Movement Rule of Pumping Units

The surface structure of the beam pumping unit consists of drive motor, reduction gearbox, and four-bar linkage. For the four-bar linkage mechanism, the fixed rod is the connection between the beam support point and the crankshaft center, and the three movable rods are crank, connecting rod, and backward beam. As shown in Figure3, R is the radius of crank, Lc is the length of the connection rod, Lb is the length of backward beam, Lf is the length of forward beam, K is the length of the fixed rod, Ω denotes the mechanical angular velocity, O and O1 are reference points, and H is the vertical height of the support point to the center of the gearbox. Setting the clockwise direction positive, the motor parameters are collected at 12 o’clock of crank, and θo denotes the angle of crank spinning clockwise.Figure 3
Four-bar linkage mechanism.The angle between the backward beam and connecting rod can be calculated by applying the cosine theorem to triangleΔLbLcJ and ΔKRJ:(7)β=cos-1⁡Lb2+Lc2-K2-R2+2KRcos⁡θ2LbLcThe angle between the crank and connecting rod can be calculated as follows:(8)α=360∘-β+φ+θThe angles ofφ when the horse head is at the top and bottom dead points are obtained as follows, respectively:(9)φbottom=cos-1⁡Lb2+K2-Lc+R22LbLc(10)φtop=cos-1⁡Lb2+K2-Lc-R22LbLcCounting from the bottom dead point (12 o’clock of crank), the motor parameters are collected every 40ms,φ can be calculated by multiplying time by the motor speed, and then the displacement is expressed as a percentage of the stroke:(11)S=φbottom-φφbottom-φtopSCwhereS is the displacement of suspension point and SC is the stroke of the sucker rod pumping system.When the beam pumping unit is working, the torque synthesized on the crankshaft by the suspension point load and the balance weight is in equilibrium with the output torque of motor.(12)P=Tm×ibelt×ireductiongearbox×ηbelt×ηreductiongearboxσ-Wcbr+Wcrcsin⁡θ0TF¯,σ=1Tm>0-1Tm<0where TF¯=Lf/Lbrsin⁡α/sin⁡β is the torque factor, ibelt and ηbelt are transmission ratio and transmission efficiency of the belt, ireductiongearbox and ηreductiongearbox are transmission ratio and transmission efficiency of the reduction gearbox, Wcb and Wc are the weight of the crank balance block and the crank, rc is the crank radius, and P is the load of the suspension point.Thus, the relation between the displacement and the crank angleθ0 is obtained from (11), and the load of the suspension point can be calculated by (12). Then the curve of dynamometer card can be obtained with variable substitution according to the displacement and load curve of the crank angle θ0.


## 4. Data-Dependent Kernel Online Sequential Extreme Learning Machine

It can be observed from the description in Section3 that there are many factors affecting the accuracy of the mechanism model such as the flux estimation, system power loss, sliding, and viscous frictions. Considering the related data information, we employ the data based regression algorithm to compensate the mechanism model by training the deviation between the output of the mechanism model and the measured value.


### 4.1. Kernel Extreme Learning Machine

Extreme Learning Machine (ELM) [24] was initially proposed for single-hidden-layer feedforward neural networks (SLFNS) by Huang and extended to the generalized single-hidden-layer feedforward neural networks where the hidden layer needs not to be neuron alike. As a supervised algorithm, its parameters are randomly generated, and then the output weights can be analytically determined by the generalized inverse operation. For a set of training samples {(xi,ti)∣xi∈Rd,ti∈Rm,i=1,2,⋯,N}, there is a nonlinear relationship between the input xi and the output ti. This nonlinear projection relationship can be expressed as a linear projection in the feature space according to the model of ELM. The outputs of ELM with L hidden layer neuron nodes are mathematically formulated as(13)ti=fLxi=∑i=1LβjGaj,bj,xi,i=1,2,…,Nwhere aj∈Rn and bj∈R are the parameters of the hidden layer of SLFN and βj∈Rm is the weight which connects the jth hidden layer node to the output layer node. G(aj,bj,xi)=g(ajxi+bj) is the output of the jth hidden layer node with respect to the inputxi. The activation functions g(·) can be sigmoid as well as the radial biases, sine, cosine, wavelet, and many other piecewise continuous computational functions. The graphic architecture of the basic ELM model with single input and single output is shown as Figure 4.Figure 4
Structure diagram of basic ELM model with single input and output.Then the hidden layer output matrixH is written as(14)H=hx1⋮hxN=Ga1,b1,x1⋯GaL,bL,x1⋮⋱⋮Ga1,b1,xN⋯GaL,bL,xNThe training purpose of ELM model is to find the optimal output weight matrixβ with the following formula:(15)β=H+Twhere T=[t1,t2,⋯tN]T denotes the output matrix and H+ is the Moore-Penrose generalized inverse of matrix. In order to calculate the norm of the output weights, we minimize the training error as follows:(16)minβ∈Rm12β2+C2∑i=1Nεi2(17)s.t.εiT=tiT-hxiβ,i=1,…,NAfter solving the quadratic optimization problem above, the output of the ELM is obtained as follows:(18)fx=hxHTIC+HHT-1T,N<LhxIC+HTH-1HTT,N≥LIf the nonlinear mapping functionh(·) is unknown to users, the corresponding ELM kernel matrix is given to users as replacement. The kernel matrix for ELM is defined as follows:(19)KELMi,j=HHT=hxi·hxj=Kxi,xiThen the output function of kernel-based ELM can be expressed as follows:(20)y=fx=Kx,x1⋮Kx,xNTIC+KELM-1TApparently, the selection and construction of kernel functions greatly influence the performance of K-ELM. Besides, as the volume of datasets increases, the training time complexity ofO(N3) and the kernel matrix size of O(N2) become an inevitable problem to concern.


### 4.2. Online Sequential Extreme Learning Machine with Kernels

In the actual oil extraction process, the measured dynamometer cards data is arriving in a streaming fashion in the actual oil extraction process. Therefore, in order to ease the computation complexity of the inverse matrix in (18) and extend the kernel extreme learning machine to the online version, Simoneet al. [25] proposed a kernel online sequential extreme learning machine(KOS-ELM) by straightforwardly extending the kernel recursive least squares to the online sequential extreme learning machine(OS-ELM) framework.According to the OS-ELM, the dataset is divided into successive minibatchesXi,i=1,2,…,B, so that X=∪i=1BXi. Similar to OS-ELM, the first approximation of the weight from the first minibatch is computed as follows for KOS-ELM:(21)P1′=C-1+kx1,x1-1(22)β1=P1′y1In combination with recursive least squares (RLS), at each new minibatchk, the output weights are updated as(23)kk+1=kxk+1,x1,…,kxk+1,xkT(24)zk+1=Pkkk+1(25)rk+1=C-1+kxk+1,xk+1-zk+1Tkk+1(26)Pk+1=rk+1-1Pkrk+1+zk+1zk+1T-zk+1-zk+1T1(27)βk+1=βk-zk+1rk+1-1ek+1rk+1-1ek+1(28)ek+1=yk+1-kk+1TβkFor each sample (xi,yi), ei denotes the error of the filter before the update and -ri-1ei is the corresponding coefficient of the center xi. And all the previous coefficients are updated by a factor -ziri-1ei.


### 4.3. Data-Dependent Kernel

In the traditional kernel ELM for regression, the selection of kernel function mostly depends on experience and the optimization of kernel parameters generally depends on intelligent optimization algorithms such as genetic algorithm, clonal selection algorithm, and particle swarm optimization algorithm. To a certain extent, these methods improve the performance of kernel based learning algorithms. However, the impact of data is not taken into account, the distribution structure of data in the kernel mapping space does not change and the performance of kernel learning algorithm cannot be fundamentally raised to higher level because the kernel function is fixed in different training samples’ structures. Thus, a novel data-dependent kernel learning method is proposed to improve the performance by replacing the basic kernel with a data-dependent kernel function [26–28]. Through conformal transformation of the basic kernel, the data-dependent kernel is defined as(29)kx,x′=cxk0x,x′cx′where K0(x,x′) is a basic kernel function such as the sigmoid kernel or the Gaussian kernel and the positive value function of x, c(x) is defined as(30)cx=b0+∑i=1Nbik′x,aiwhere the whole samples are chosen as “empirical cores” {ai=xi∈Rd,i=1,2,…N}, N is the number of samples, bi represents the combination expansion coefficients of ai, k(x,ai)=e-δx-ai2, and δ is a free parameter. Thus, the geometrical structure of the data in the kernel mapping space is determined by the expansion coefficients bi(i=1,2,⋯,N). The kernel matrices corresponding to k and k0 are denoted as K and K0, respectively. Consequently, (29) can be rewritten as(31)K=CK0Cwhere the diagonal matrix C=diag{c(x1),c(x2),…c(xN)}. After defining c=[c(x1),c(x2),…c(xN)]T and b=[b1,b2,…bN]T, we have(32)c=1k′x1,a1⋯k′x1,aN1k′x2,a1⋯k′x2,aN⋮⋮⋱⋮1k′xN,a1⋯k′xN,aNb0b1⋮bN=KNb


### 4.4. Kernel Optimization

There are many methods to optimize the combination expansion coefficients and the Gaussian kernel parameterσ2. For two key uncertain parameters, the improved free search algorithm (IFSA) [29] is adopted to address the optimal solution, which is put forward to find the optimal objective function value with the following equations:(33)ftj=fxtji,fj=minftj(34)Phj=min⁡fjfjf denotes the position where an individual completes a search step; t denotes the search step; j denotes the number of individuals; and i denotes the search spatial dimension.As a major parameter of the algorithm, the search radiusRji is adjusted in real time. The initial value of it is set to be 1, and with the growth of the search step, the radius shrinks as (35). In addition, the sensitivity parameter determines the scope of the pheromone value and also confirms the initial search point of the next search round, so the modification of sensitivity parameters is shown as (36) and (37).(35)Rji=ρRjt-g,Rjt≥RminRmin,else(36)Sej=Semin+ΔSej(37)ΔSej=Semax-Semin×t-GGThus, the implementation steps of improved free search algorithm can be summarized as Algorithm1.Algorithm 1:Improved free search algorithm.
Input: The initial population size M, max algebra G, search step length T, the parameters σ2 and b; Output: The optimal parameter σ2 and b; Step1: Initialization 1.1 Normalize the input parameters σ2 and b with the range of [-1,1]; 1.2 The initial pheromone is obtained with the initial value, which is released at the same time; 1.3 The initial search results are obtained. Step2: Searching process 2.1 Calculate the sensitivity parameter with (36) and (37); 2.2 Confirm the new search point of the next search round using the sensitivity parameter; 2.3 Record the RMSE value of the error between the predicted value and the actual value as the fitness function. Fitness function is defined as follows; f j = 1 n ∑ i = 1 n y i - y ^ i 2 2.4 The new search results can be obtained by (33) and (34); 2.5 Preserve the best individuals; 2.6 Adjust the search radius according to (35). Step3: Termination condition judgement When the termination condition is satisfied, output the optimal parameterσ2 and b.By applying the data-dependent kernel to the kernel online sequential extreme learning machine for regression, the training framework of the proposed IFS-DDKOS-ELM error model is summarized as Algorithm2.Algorithm 2:IFS-DDKOS-ELM.
Input: Training samples {(xi,ti)∣xi∈Rd,ti∈Rm,i=1,2,⋯,N}, the Gaussian kernel function is selected as the basic kernel function K0(x,x′)=exp(-x-x′2/2σ2), the nonzero constant C, the number of hidden nodes L and the number of data blocks B in a data set. Output: The predictive value of the polished rod load and suspension point displacement ( 1 ) Randomly generate the input layer parameters (aj,bj)(j=1,2,…,L); ( 2 ) Construct the data-dependent kernel function Kx,x′ by modifying the initial basic kernel function according to (29), whose parameters are optimized by Algorithm 1; ( 3 )  for each k∈[1,B]  do ( 4 )  Update the output weight β by (27); ( 5 )  end for ( 6 ) Compute the predictive value according to (20).


## 4.1. Kernel Extreme Learning Machine

Extreme Learning Machine (ELM) [24] was initially proposed for single-hidden-layer feedforward neural networks (SLFNS) by Huang and extended to the generalized single-hidden-layer feedforward neural networks where the hidden layer needs not to be neuron alike. As a supervised algorithm, its parameters are randomly generated, and then the output weights can be analytically determined by the generalized inverse operation. For a set of training samples {(xi,ti)∣xi∈Rd,ti∈Rm,i=1,2,⋯,N}, there is a nonlinear relationship between the input xi and the output ti. This nonlinear projection relationship can be expressed as a linear projection in the feature space according to the model of ELM. The outputs of ELM with L hidden layer neuron nodes are mathematically formulated as(13)ti=fLxi=∑i=1LβjGaj,bj,xi,i=1,2,…,Nwhere aj∈Rn and bj∈R are the parameters of the hidden layer of SLFN and βj∈Rm is the weight which connects the jth hidden layer node to the output layer node. G(aj,bj,xi)=g(ajxi+bj) is the output of the jth hidden layer node with respect to the inputxi. The activation functions g(·) can be sigmoid as well as the radial biases, sine, cosine, wavelet, and many other piecewise continuous computational functions. The graphic architecture of the basic ELM model with single input and single output is shown as Figure 4.Figure 4
Structure diagram of basic ELM model with single input and output.Then the hidden layer output matrixH is written as(14)H=hx1⋮hxN=Ga1,b1,x1⋯GaL,bL,x1⋮⋱⋮Ga1,b1,xN⋯GaL,bL,xNThe training purpose of ELM model is to find the optimal output weight matrixβ with the following formula:(15)β=H+Twhere T=[t1,t2,⋯tN]T denotes the output matrix and H+ is the Moore-Penrose generalized inverse of matrix. In order to calculate the norm of the output weights, we minimize the training error as follows:(16)minβ∈Rm12β2+C2∑i=1Nεi2(17)s.t.εiT=tiT-hxiβ,i=1,…,NAfter solving the quadratic optimization problem above, the output of the ELM is obtained as follows:(18)fx=hxHTIC+HHT-1T,N<LhxIC+HTH-1HTT,N≥LIf the nonlinear mapping functionh(·) is unknown to users, the corresponding ELM kernel matrix is given to users as replacement. The kernel matrix for ELM is defined as follows:(19)KELMi,j=HHT=hxi·hxj=Kxi,xiThen the output function of kernel-based ELM can be expressed as follows:(20)y=fx=Kx,x1⋮Kx,xNTIC+KELM-1TApparently, the selection and construction of kernel functions greatly influence the performance of K-ELM. Besides, as the volume of datasets increases, the training time complexity ofO(N3) and the kernel matrix size of O(N2) become an inevitable problem to concern.


## 4.2. Online Sequential Extreme Learning Machine with Kernels

In the actual oil extraction process, the measured dynamometer cards data is arriving in a streaming fashion in the actual oil extraction process. Therefore, in order to ease the computation complexity of the inverse matrix in (18) and extend the kernel extreme learning machine to the online version, Simoneet al. [25] proposed a kernel online sequential extreme learning machine(KOS-ELM) by straightforwardly extending the kernel recursive least squares to the online sequential extreme learning machine(OS-ELM) framework.According to the OS-ELM, the dataset is divided into successive minibatchesXi,i=1,2,…,B, so that X=∪i=1BXi. Similar to OS-ELM, the first approximation of the weight from the first minibatch is computed as follows for KOS-ELM:(21)P1′=C-1+kx1,x1-1(22)β1=P1′y1In combination with recursive least squares (RLS), at each new minibatchk, the output weights are updated as(23)kk+1=kxk+1,x1,…,kxk+1,xkT(24)zk+1=Pkkk+1(25)rk+1=C-1+kxk+1,xk+1-zk+1Tkk+1(26)Pk+1=rk+1-1Pkrk+1+zk+1zk+1T-zk+1-zk+1T1(27)βk+1=βk-zk+1rk+1-1ek+1rk+1-1ek+1(28)ek+1=yk+1-kk+1TβkFor each sample (xi,yi), ei denotes the error of the filter before the update and -ri-1ei is the corresponding coefficient of the center xi. And all the previous coefficients are updated by a factor -ziri-1ei.


## 4.3. Data-Dependent Kernel

In the traditional kernel ELM for regression, the selection of kernel function mostly depends on experience and the optimization of kernel parameters generally depends on intelligent optimization algorithms such as genetic algorithm, clonal selection algorithm, and particle swarm optimization algorithm. To a certain extent, these methods improve the performance of kernel based learning algorithms. However, the impact of data is not taken into account, the distribution structure of data in the kernel mapping space does not change and the performance of kernel learning algorithm cannot be fundamentally raised to higher level because the kernel function is fixed in different training samples’ structures. Thus, a novel data-dependent kernel learning method is proposed to improve the performance by replacing the basic kernel with a data-dependent kernel function [26–28]. Through conformal transformation of the basic kernel, the data-dependent kernel is defined as(29)kx,x′=cxk0x,x′cx′where K0(x,x′) is a basic kernel function such as the sigmoid kernel or the Gaussian kernel and the positive value function of x, c(x) is defined as(30)cx=b0+∑i=1Nbik′x,aiwhere the whole samples are chosen as “empirical cores” {ai=xi∈Rd,i=1,2,…N}, N is the number of samples, bi represents the combination expansion coefficients of ai, k(x,ai)=e-δx-ai2, and δ is a free parameter. Thus, the geometrical structure of the data in the kernel mapping space is determined by the expansion coefficients bi(i=1,2,⋯,N). The kernel matrices corresponding to k and k0 are denoted as K and K0, respectively. Consequently, (29) can be rewritten as(31)K=CK0Cwhere the diagonal matrix C=diag{c(x1),c(x2),…c(xN)}. After defining c=[c(x1),c(x2),…c(xN)]T and b=[b1,b2,…bN]T, we have(32)c=1k′x1,a1⋯k′x1,aN1k′x2,a1⋯k′x2,aN⋮⋮⋱⋮1k′xN,a1⋯k′xN,aNb0b1⋮bN=KNb


## 4.4. Kernel Optimization

There are many methods to optimize the combination expansion coefficients and the Gaussian kernel parameterσ2. For two key uncertain parameters, the improved free search algorithm (IFSA) [29] is adopted to address the optimal solution, which is put forward to find the optimal objective function value with the following equations:(33)ftj=fxtji,fj=minftj(34)Phj=min⁡fjfjf denotes the position where an individual completes a search step; t denotes the search step; j denotes the number of individuals; and i denotes the search spatial dimension.As a major parameter of the algorithm, the search radiusRji is adjusted in real time. The initial value of it is set to be 1, and with the growth of the search step, the radius shrinks as (35). In addition, the sensitivity parameter determines the scope of the pheromone value and also confirms the initial search point of the next search round, so the modification of sensitivity parameters is shown as (36) and (37).(35)Rji=ρRjt-g,Rjt≥RminRmin,else(36)Sej=Semin+ΔSej(37)ΔSej=Semax-Semin×t-GGThus, the implementation steps of improved free search algorithm can be summarized as Algorithm1.Algorithm 1:Improved free search algorithm.
Input: The initial population size M, max algebra G, search step length T, the parameters σ2 and b; Output: The optimal parameter σ2 and b; Step1: Initialization 1.1 Normalize the input parameters σ2 and b with the range of [-1,1]; 1.2 The initial pheromone is obtained with the initial value, which is released at the same time; 1.3 The initial search results are obtained. Step2: Searching process 2.1 Calculate the sensitivity parameter with (36) and (37); 2.2 Confirm the new search point of the next search round using the sensitivity parameter; 2.3 Record the RMSE value of the error between the predicted value and the actual value as the fitness function. Fitness function is defined as follows; f j = 1 n ∑ i = 1 n y i - y ^ i 2 2.4 The new search results can be obtained by (33) and (34); 2.5 Preserve the best individuals; 2.6 Adjust the search radius according to (35). Step3: Termination condition judgement When the termination condition is satisfied, output the optimal parameterσ2 and b.By applying the data-dependent kernel to the kernel online sequential extreme learning machine for regression, the training framework of the proposed IFS-DDKOS-ELM error model is summarized as Algorithm2.Algorithm 2:IFS-DDKOS-ELM.
Input: Training samples {(xi,ti)∣xi∈Rd,ti∈Rm,i=1,2,⋯,N}, the Gaussian kernel function is selected as the basic kernel function K0(x,x′)=exp(-x-x′2/2σ2), the nonzero constant C, the number of hidden nodes L and the number of data blocks B in a data set. Output: The predictive value of the polished rod load and suspension point displacement ( 1 ) Randomly generate the input layer parameters (aj,bj)(j=1,2,…,L); ( 2 ) Construct the data-dependent kernel function Kx,x′ by modifying the initial basic kernel function according to (29), whose parameters are optimized by Algorithm 1; ( 3 )  for each k∈[1,B]  do ( 4 )  Update the output weight β by (27); ( 5 )  end for ( 6 ) Compute the predictive value according to (20).


## 5. Hybrid Model Structure

Based on the above discussion, the parallel hybrid model is constructed in Figure5, which not only considers the processing technology but also makes full use of relevant data information.Figure 5
The prediction process of hybrid model.The inputs of the hybrid model are the mechanism and electric motor parameters including stroke, the crank radius, stroke frequency, pulley diameter, unbalanced weight, crank quality, balance block weight, three line voltages and line currents of the stator, motor speed, and other technical parameters. The outputs of the hybrid model are expressed as follows:(38)Po=P+ePwhere load P and displacement S are the outputs of the mechanism model and eP is the output of the IFS optimized data-dependent kernel extreme learning machine (IFS-DDKOS-ELM), which is the error estimation between the actual measured value and the output of the mechanism model.


## 6. Case Research and Experimental Results

In this section, the proposed hybrid model is firstly employed to predict the load and the displacement of the suspension point of the beam pumping unit. Secondly, the curves of polished rod load and suspension point displacement with respect to crank angleθ0 are drawn by the outputs of the hybrid model respectively. Finally, the dynamometer card is drawn according to the curves of the polished load and suspension point through variable substitution. Taking a normal oil well in Liaohe oil field of China for case research, the technical parameters of the beam pumping unit and the parameters of the induction motor are shown in Tables 1 and 2, respectively.Table 1
Technical parameters of the beam pumping unit.
Number of oil wells CYJ10-3-53HB Stroke number per minute (min−1) 5 Transmission ratio of belt-pulley gear box 150.82 Crank radius (m) 1.134 Length of the link rod (m) 3.60 Length of the beam forearm (m) 3.8 Length of the beam backward arm (m) 2.5 Vertical height H (m) 3.2 Distance between the reference points K (m) 4.53 Transmission efficiency 0.95 Unbalanced weight (kN) -1.12 Balance block weight (kN) 1.40Table 2
Electronic parameters of the induction motor.
Number of oil wells CYJ10-3-53HB Rated power (kW) 25 Rated voltage (V) 380 Rated rotor speed (rad⋅min-1) 585 Stator resistance (Ω) 0.36 Rotor resistance (Ω) 0.19 Stator inductance (mH) 33 Rotor inductance (mH) 35 Peak value of the mutual inductance between rotor and stator (mH) 41 Moment of inertia (kg⋅m2) 1.71 Rotational resistance coefficient 0.019 Pole-pairs number 3 Motor transmission efficiency 0.98The data to establish the hybrid model is collected with the indicator instrument, the electric parameters measurement instrument, and two proximity switches. The indicator instrument installed at the suspension point collects the surface dynamometer cards periodically when the beam pumping units are working. The collected dynamometer cards are taken as the reference data and the historical training data for the IFS-DDK-ELM error predictive model. The electric parameters measurement instrument is installed in the electric control cabinet connecting to the induction motor with the sampling interval of 47.8ms, and the collected current, voltage, and power of the induction motor are taken as the input of the hybrid model. One proximity switch is installed at 12 o’clock on the pedestal of the crank shaft to measure the cycle of one stroke, when the crank is rotated to the 12 o’clock position, triggering the proximity switch action, and the timer of the electric parameters measurement instrument starts ticking until the crank is rotated again to 12 o’clock. Therefore, a periodic displacement curve can be obtained. The other one is installed on the motor shaft to measure the motor speed.The data collection instruments are installed behind the electronic control cabinet (as shown in Figure6(a)), and it mainly includes the RTU (Remote Terminal Unit), the electrical parameter acquisition board with LCD screen, and the industrial PC (as shown in Figure 6(b)). The RTU is used to receive the surface dynamometer cards collected with the wireless integrated indicator instrument installed at the suspension point of the horse head. The function of the electrical parameter acquisition board is data collection and dynamometer cards inversion; the industrial PC is used to preserve the work data and generate work reports. And the communication mode among the three modules is TCP/IP.Figure 6
The data collection instruments installed behind the electronic control cabinet.
(a) (b)1800 groups of measurement data in total three months from Aug 9, 2017, to Nov 10, 2017, including the electric parameters and the dynamometer cards are collected to establish the data-driven model. Randomly select 1200 groups of the collected data as the training samples and the other 600 groups of data as the test samples. Figure7 shows the suspension point displacement curve in one stroke, Figure 8 is the torque factor curve, and Figure 9 is the predictive curve of the polished rod load obtained by the hybrid model. Consequently, the predictive result of dynamometer card is obtained as Figure 10. The test error of maximum load, minimum load, and the area of the dynamometer card obtained by the hybrid model are shown as Table 3. It can be concluded from Figure 10 and Table 3 that the predictive dynamometer card of the proposed hybrid model is basically consistent with the measured dynamometer card and can reflect the actual working operation vividly.Table 3
Comparison of performance indicators between predictive value and actual value.
Predictive value Actual value Test error (%) Maximum load (kN) 82.87 82.22 0.7905 Minimum load (kN) 44.49 44.32 0.3835 Area of dynamometer card 67.0445 67.2565 0.3152Figure 7
The suspension point displacement curve.Figure 8
The torque factor curve.Figure 9
Polished rod load curve.Figure 10
The prediction result of dynamometer card.In order to approve the effectiveness of the proposed hybrid prediction method, the mechanism predictive model without data-driven model and pure IFS-DDKOS-ELM model are chosen to compare the prediction performance, and the prediction performances are shown as Figure11(a). In addition, combination kernel function LSSVM and the KOS-ELM in [25] are selected as the comparison algorithms to verify the effectiveness of the IFS-DDKOS-ELM. The prediction effect of the KOS-ELM in [25] is proved to be better than ELM, K-ELM, and OS-ELM. So this paper directly chooses these two methods to compare the prediction performance of the IFS-DDKOS-ELM model. The comparison results are shown in Figure 11(b).Figure 11
The predictive contrast results of dynamometer card.
(a) (b)Root mean square error (RMSE) and mean absolute error (MAE) are used as the prediction evaluation indicators, which are defined as follows:(39)RMSE=1n∑i=1nyi-y^i2(40)MAE=1n∑i=1nyi-y^iwhere n denotes the number of test data, yi is the real value of test data, and y^i is the prediction value. After operating the algorithms introduced above 50 times, the average is selected as the experimental result, which are shown as Tables 4 and 5.Table 4
Comparison of the prediction results with the three models by RMSE and MAE.
Predictive methods RMSE MAE Hybrid model 0.2383 0.1731 Mechanism model 1.4829 1.1775 Pure IFS-DDKOS-ELM 1.4365 1.1780Table 5
Comparison of the prediction results with the three data-driven models by RMSE and MAE.
Predictive methods RMSE MAE IFS-DDKOS-ELM model 1.4365 1.1780 LS-SVM model 3.4873 2.9140 KOS-ELM model 2.5195 2.1492From the prediction performance of Figure11(b) and Table 5, it can be known that the IFS-DDKOS-ELM model of this paper has the better prediction effect for dynamometer card inversion than the other two data-driven models. Although the LS-SVM with combination kernel function and the KOS-ELM in [25] show good prediction effect, the prediction accuracy is further improved by employing the optimized data-dependent kernel function. From the prediction performance of Figure 11(a), it can be observed that the proposed hybrid model has better prediction accuracy than the pure IFS-DDKOS-ELM and mechanism models, which also can be concluded from the corresponding RMSE and MAE with different models shown in Table 4. Notably, because of the stable operating condition of the motor, the fluctuation of the prediction dynamometer card with pure mechanism model is relatively weak.In summary, the experiment results demonstrate that the proposed hybrid model can perform better dynamometer cards inversion effect than the pure mechanism model and the pure data-driven model. And it can produce precise dynamometer cards for fault diagnosis of the beam pumping unit.


## 7. Conclusion

In actual oil production, real-time and accurate collection of dynamometer cards has far-reaching influence on the oil well fault diagnosis. In view of the insufficient data in the actual oilfield, the data collection is not timely. This paper proposes a novel hybrid modeling for dynamometer cards inversion based on the electronic parameters collected from the motor. The simulation results can basically match with the practical measurement dynamometer cards. The advantages of the proposed hybrid model are summarized as follows.( 1 ) The mechanism model is carried out to indirectly measure the polished rod load and suspension displacement from the perspective of mathematical model of AC motor. Besides, in the process of data collection, the proximity switch is used to realize the cycle identification without sensor.( 2 ) In order to compensate the error of the mechanism model output, an online learning with a data-dependent kernel-based ELM (DDKOS-ELM) is firstly presented for regression. The comparison results show that the performance of the data-dependent kernel function is proved to be better than the combination kernel function.



---

*Source: 6730905-2018-06-27.xml*