# Asymptotic Method of Solution for a Problem of Construction of Optimal Gas-Lift Process Modes

Authors: Fikrat A. Aliev; Mutallim M. Mutallimov; Idrak M. Askerov; Iouldouz S. Raguimov
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

Mathematical model in oil extraction by gas-lift method for the case when the reciprocal value of well's depth represents a small parameter is considered. Problem of optimal mode construction (i.e., construction of optimal program trajectories and controls) is reduced to the linear-quadratic optimal control problem with a small parameter. Analytic formulae for determining the solutions at the first-order approximation with respect to the small parameter are obtained. Comparison of the obtained results with known ones on a specific example is provided, which makes it, in particular, possible to use obtained results in realizations of oil extraction problems by gas-lift method.


---

## Body

## 1. Introduction

It is known [1–3] that the gas-lift technique of exploitation of oil wells is widely used when the gushing method does not work for the reason of insufficiency in pressure. The essence of the gas-lift method consists of the fact that by the mean of energy of injected underground gas it is possible to lift fluid to the surface.While the gas-lift method is widely used in oil extraction for a sufficiently long period of time, construction of an adequate mathematical model is rather an actual problem. Mathematical model describing the oil lifting process in pump-compressor tubes is described in [4].In [4, 5], using this model, an optimal control problem for gas-lift process is formulated, where the pressure and volume of the injected gas is used as a control parameter. In the same papers, the method of straight lines is used to reduce the optimal control problem to the linear-quadratic optimal control problem (LQOCP). It should be noted that, when the depth of a well is large, in order to obtain sufficiently accurate results it is necessary to divide this distance into large number of relatively small segments. The latter, in its turn increases the dimensions of the system and, consequently, results in higher volume of computations to solve such problems. Thus, generally speaking, solution of these problems yields higher approximation and calculation errors.Therefore, it seems quite rational to develop other solution methods with lower computational complicatedness [6–8]. One of the methods applicable to the considered problem for the case when the reciprocal value of well’s depth represents a small parameter is the asymptotic method described in [9].In the present paper, using the method of straight lines, the mathematical model of the gas-lift process is represented as a system of partial differential equations of hyperbolic type. In a particular case, we arrive to LQOCP for a system of ordinary differential equations containing a small parameter. Applying the algorithm introduced in [10, 11], the solution to LQOCP is obtained as a function of the small parameter. Consequently, as in [6, 9], at the first approximation, analytic formulas for the volume of injected gas (control) and production level (trajectories) are obtained. As it supposed, the presented approach will reduce substantially the amount of required computations.


## 2. Mathematical Formulation of the Problem

As in [4], mathematical model of gas-fluid mixture flow in pipes is described by the system of hyperbolic-type partial differential equations:
(2.1)∂P∂t=-c2F̅∂Q∂x∂Q∂t=-F̅∂P∂x-2aQt≥0,x∈[0,2L],
which for x=z/2Land ε=1/2L can be written as(2.2)∂P∂t=-c2F∂Q∂zε∂Q∂t=-F∂P∂zε-2aQ,z∈(0,1),
with appropriate boundary and initial conditions(2.3)P(z,0)=P0(z),Q(z,0)=Q0(z),P(0,t)=P0(t),Q(0,t)=Q0(t),P(L+0,t)=P(L-0,t)+Ppl(t),Q(L+0,t)=Q(L-0,t)+Qpl(t).
It is required to determine a control(2.4)U(t)=[P0(t),Q0(t)]′
minimizing the functional(2.5)J=12α(Q(2L,T)-Qdeb)2+12∫0TU′(t)RU(t)dt,
where P is pressure, Q is gas-fluid mixture volume andQdebis the desired yield.Applying the so-called method of straight lines to (2.2) and taking l=1/N,  we obtain(2.6)dPkdt=-c2εF·l(Qk-Qk-1),dQkdt=-Fεl(Pk-Pk-1)-2aQk.
Note that for k=N+1 (2.6) can be written as (2.7)P.N+1=-c22εF2lQN+1+c22εF2lQN+c2εF2lQpl,Q.N+1=-F2εlPN+1+F2εlPN-2a2QN+1+F2εlPpl,
where Qpl,Ppl  denote gas-fluid outlay (yield) and pressure at the bottom of a well, respectively. As in [10], a linear-quadratic optimal control problem is formulated for this system. It is required to find x,u satisfying the equation(2.8)ẋ=(A0+A1ε)x+Bεu+Cε
with initial condition(2.9)x(0)=x0
such that the value of the functional(2.10)J=12(x(T)-x̅)′K(x(T)-x̅)+12∫0Tu′Rudt
is minimized.Here,(2)A1=[A11000⋯0000⋯0000A12A1100⋯0000⋯00000A12A110⋯0000⋯000000A12A11⋯0000⋯0000⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯0000⋯A11000⋯00000000⋯A12A1100⋯00000000⋯0A21A220⋯00000000⋯00A21A22⋯0000⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯0000⋯0000⋯A220000000⋯0000⋯A21A22000000⋯0000⋯0A21A2200000⋯0000⋯00A21A22],A0=[0000⋯0000⋯00000-2a100⋯0000⋯00000000⋯0000⋯0000000-2a1⋯0000⋯0000⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯0000⋯0000⋯00000000⋯0-2a100⋯00000000⋯0000⋯00000000⋯000-2a2⋯0000⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯0000⋯0000⋯00000000⋯0000⋯0-2a2000000⋯0000⋯00000000⋯0000⋯000-2a2],A11=[0-c12F1l-F1l0],A12=[0c12F1lF1l0],A22=[0-c22F2l-F2l0],A21=[0c22F2lF2l0],x=[P1,Q1,P2,Q2,…,PN,QN,…,P2N,Q2N]′,x0=[P10,Q10,P20,Q20,…,PN0,QN0,…,P2N0,Q2N0]′,u=[P0Q0],B=[0c12F1lF1l000⋯⋯⋯⋯000000],C=[00⋯⋯000c22F2lF2l000⋯⋯00][PplQpl].
Let R=R̅ε. Then the corresponding Euler-Lagrange control problem can be written as(2.12)[ẋλ̇]=[A0+A1ε-εBR̅-1B'0-A01-A11ε][xλ]+[Cε0],(2.13)x(0)=x0,λ(T)=N(x(T)-x̅).
Thus, initial problem (2.2)–(2.5) is reduced to finding the solution of problem (2.12), where ε is a small parameter. Therefore, the asymptotic method (see [9]) allowing to expand the solution of system (2.12) with respect to small parameterε can be applied.


## 3. Application of Asymptotic Method

Let us apply the asymptotic method to the Euler-Lagrange equation (2.12) with boundary conditions for x(0)andλ(T).  According to [12], (3.1)[x(T)λ(T)]=e[A0+εA1-εBR̅-1B′0-A0′-εA1′]T[x0λ(0)]+[e[A0+εA1-εBR̅-1B′0-A0′-εA1′]T-[E00E]][A0+εA1-εBR̅-1B′0-A0′-εA1′]-1[Cε0].
Further, let us expand the expression(3.2)e[A0+εA1-εBR̅-1B′0-A0′-εA1′]T=e[A000-A0′]T+ε[A1-BR̅-1B′0-A1′]T
from (3.1) with respect to ε. If we denote  (3.3)H1=i.e[A000-A0′],H2=i.e[A1-BR̅-1B'0-A1′],
then according to [12], the expansion of expression (3.2) with respect to ε can be represented as:
(3.4)eH1T+εH2T≈eH1T+ε∫01eH1T(1-s)H2TeH1Tsds
Denote the integral in (3.4) by L0. Then it is not difficult to show that matrix L0 is a solution of the following Sylvester's equation:(3.5)H1L0-L0H1=eH1TH2-H2eH1T.
Therefore, for expression (3.2) we obtain the expansion(3.6)eH1T+εH2T≈eH1T+εL0.
Introducing notations (3.7)L0=[L1L2L3L4],[A0+εA1-εBR̅-1B′0-A0′-εA1′]-1=[S1S20S4],
expression (3.1) can be simplified to (3.8)[x(T)λ(T)]=[eA0Tx(0)e-A0'Tλ(0)]+ε[L1x(0)+L2λ(0)+eA0TS1C-S1CL3x(0)+L4λ(0)].
Hence, adding boundary conditions from (2.12), we arrive to the following system of algebraic equations: (3.9)x(T)-εL2λ(0)=eA0Tx(0)+ε(eA0TS1C-S1C+L1x(0)),λ(T)-(e-A0′T+εL4)λ(0)=εL3x(0),λ(T)-Nx(T)=-Nx̅
which in the matrix form can be written as(3.10)[E-εL200-e-A0′T-εL4EN0-E][x(T)λ(0)λ(T)]=[eA0Tx0+ε(eA0TS1C-S1C+L1x0)εL3x0Nx̅].
Multiplying (3.10) by [E000EE00E]  on the left, we obtain(3.11)[E-εL20N-e-A0′T-εL40N0-E][x(T)λ(0)λ(T)]=[eA0Tx0+ε(eA0TS1C-S1C+L1x0)Nx̅+εL3x0Nx̅].
Hence, if we denote the coefficient matrix in  (3.11)  by M, then it can be written as(3.12)M=[E-εL20N-e-A0′T-εL40N0-E]=[K(ε)0K1-E],
where(3.13)K(ε)=[E-εL2N-e-A0′T-εL4],K1=[N0].
It is not difficult to show that(3.14)[K(ε)0K1-E]-1=[K-1(ε)0K1K-1(ε)-E].
Hence, using the fact that(3.15)K-1(ε)≈K-1(0)-K-1(0)K.(0)K-1(0)ε
we obtain the inverse matrix (3.16)M-1≈[E+L2eA0′TNε-L2eA0′Tε0eA0′TN-eA0′TL4eA0′TNε+eA0′TNL2eA0′TNεeA0′T-eA0′TNL2eA0′Tε+eA0′TL4ε0N+NL2eA0′TNε-NL2eA0′Tε-E]
and, consequently, multiplying the both sides of (3.11)  by M-1 on the left, we obtain the following analytic formulae to determine values of x(T),λ(0),λ(T): (3.17)[x(T)λ(0)λ(T)]=[eA0Tx0eA0′TNeA0Tx0+eA0′TNx̅NeA0Tx0-Nx̅]+ε[eA0TS1C-S1C+L1x0+L2eA0′TNeA0Tx0-L2eA0′TNx̅{e-A0′TNeA0TS1C-eA0′TNS1C+eA0′TNL1x0+eA0′TL3x0-eA0′TL4eA0′TNeA0Tx0+eA0′TNL2eA0′TNeA0Tx0-eA0′TNL2eA0′TNx̅+eA0'TL4eA0′TNx̅}NL2eA0′TNeA0Tx0-NL2eA0′TNx̅+NL1x0+NeA0TS1C-NS1C].
Further, as in (3.1), using determined in (3.17) value λ(0), it is possible to find  x(ti),λ(ti) in the form(3.18)[x(ti)λ(ti)]=e[A0+εA1-εBR̅-1B′0-A0′-εA1′]ti[x0λ(0)]+[e[A0+εA1-εBR̅-1B′0-A0′-εA1′]-[E00E]][A0+εA1-εBR̅-1B′0-A0-εA1]-1[Cε0]
for every ti∈[0T].Therefore, expansions forx(ti)and u(ti) with respect to ε can be obtained in the form(3.19)x(ti)=eA0ti+ε(eA0tiS1C-S1C+L11x0-L22eA0'tiNeA0tix0+L22eA0'tiNx̅),λ(ti)=-NeA0Tx0+Nx̅+ε(NL22eA0'tiNeA0tix0-NL22eA0'tiNx̅-NL11x0-NeA0tiS1C-NS1C),u(ti)=-R-1B'λ(ti).
Here L11,L22 are block matrices of matrix-solutions Li of Sylvester's equation [12, 13](3.20)H1Li-LiH1=eH1tiH2-H2eH1ti
for every ti.


## 4. Computational Experiments

On the basis of the obtained formulae, we have developed an algorithm and computer program using MATLAB system. In [5] an algorithm for solving problem (2.8)–(2.10) for the case when C=0,N=2 is developed. It is clear that for N>2 the algorithm will yield more accurate solutions. However, for computational comparison of the asymptotic method with the method described in [5], authors have considered only the case whenC=0,N=2. It is possible to prove that in this case Sylvester's equation has infinitely many solutions. While using a special program in MATLAB one particular solution is determined.In Figures1, 2, 3, and 4 the graphs of the functions P1,Q1,P2 andQ2 for both algorithms are given. Namely, the graphs for the functions obtained using the algorithm introduced in [5, 14, 15] are shown in solid lines and by the algorithm elaborated using the asymptotic method are given in dashed lines. A comparative analysis shows that the results for the functionsP1,Q1,P2, and Q2 for both algorithms coincide with a sufficiently high accuracy.Dependence ofP1 on  t.
(a)(b)Dependence ofP2 on t.
(a)(b)Dependence ofQ1 on t.
(a)(b)Dependence ofQ2 on  t.
(a)(b)



---

*Source: 191053-2010-07-27.xml*