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

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.

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𝑎𝑄𝑡≥0,𝑥∈[0,2𝐿],(2.1) which for 𝑥=𝑧/2𝐿and 𝜀=1/2𝐿 can be written as𝜕𝑃𝜕𝑄𝜕𝑡=−𝑐𝐹²𝜕𝜕𝑄𝑧𝜀𝜕𝑡=−𝐹𝜕𝜕𝑃𝑧𝜀−2𝑎𝑄,𝑧∈(0,1),(2.2) with appropriate boundary and initial conditions𝑃𝑃((𝑧𝐿,0+)0=,𝑡𝑃)⁰=(𝑧𝑃),(𝐿𝑄−(0𝑧,,𝑡0))+=𝑃𝑄⁰_𝑝𝑙(𝑧),𝑃(0,𝑡)=𝑃₀(𝑡),𝑄(0,𝑡)=𝑄₀(𝑡),(𝑡),𝑄(𝐿+0,𝑡)=𝑄(𝐿−0,𝑡)+𝑄_𝑝𝑙(𝑡).(2.3) It is required to determine a control𝑈(𝑡)=𝑃₀(𝑡),𝑄₀(𝑡)^(2.4) minimizing the functional𝐽=12𝛼𝑄(2𝐿,𝑇)−𝑄^deb²+12^𝑇0𝑈^(𝑡)𝑅𝑈(𝑡)𝑑𝑡,(2.5) where 𝑃 is pressure, 𝑄 is gas-fluid mixture volume and𝑄_debis the desired yield.

Applying the so-called method of straight lines to (2.2) and taking 𝑙=1/𝑁, we obtain𝑑𝑃_𝑘𝑑𝑡=−𝑐²𝜀𝐹⋅𝑙𝑄^𝑘𝑑𝑄_𝑘−𝑄_𝑘−1,𝑑𝑡=−𝐹𝑙𝜀𝑃^𝑘−𝑃_𝑘−1−2𝑎𝑄_𝑘.(2.6) Note that for 𝑘=𝑁+1 (2.6) can be written as 𝑃̇_𝑁+1=−𝑐²²𝜀𝐹₂𝑙𝑄^𝑁+1+𝑐²²𝜀𝐹₂𝑙𝑄^𝑁+𝑐²𝜀𝐹₂𝑙𝑄^𝑝𝑙𝑄̇_𝑁+1,=−𝐹²𝜀𝑙𝑃^𝑁+1+𝐹²𝜀𝑙𝑃^𝑁−2𝑎₂𝑄_𝑁+1+𝐹²𝜀𝑙𝑃^𝑝𝑙,(2.7) where 𝑄_𝑝𝑙,𝑃_𝑝𝑙 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 𝑥,𝑢 satisfying the equation𝑥̇=𝐴₀+𝐴₁𝜀𝑥+𝐵𝜀𝑢+𝐶𝜀(2.8) with initial condition𝑥(0)=𝑥⁰(2.9) such that the value of the functional𝐽=12𝑥(𝑇)−𝑥^𝐾𝑥(𝑇)−𝑥+12^𝑇0𝑢^𝑅𝑢𝑑𝑡(2.10) is minimized.

Here, 𝐴₁=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐴¹¹𝐴000⋯0000⋯0000₁₂0𝐴𝐴₁₁₁₂00⋯0000⋯000000𝐴𝐴₁₁₁₂0⋯0000⋯0000⋯0⋯0⋯0𝐴⋯0₁₁⋯⋯⋯𝐴⋯0₁₁⋯0⋯0⋯0⋯⋯⋯0⋯0⋯0⋯00000⋯𝐴₁₂000⋯00000000⋯0𝐴𝐴₁₁₂₁00⋯00000000⋯00𝐴𝐴₂₂₂₁0⋯0000⋯0⋯0⋯0⋯0⋯⋯⋯0⋯0⋯0𝐴⋯0₂₂⋯⋯⋯𝐴⋯0₂₂⋯0⋯0⋯00000⋯0000⋯𝐴₂₁0000000⋯0000⋯0𝐴𝐴₂₂₂₁000000⋯0000⋯00𝐴𝐴₂₂₂₁0𝐴₂₂𝐴⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(2)₀=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣00−20𝑎¹00⋯0000⋯00000000000−200𝑎₁⋯⋯00000000⋯⋯00000000⋯00⋯00⋯00⋯00⋯⋯⋯⋯⋯000−⋯200𝑎₁⋯00⋯00⋯⋯⋯⋯00⋯00⋯00⋯0000000000⋯⋯0000000−200𝑎₂⋯⋯00000000⋯00⋯00⋯00⋯00⋯⋯⋯⋯00⋯00⋯00⋯00⋯⋯⋯⋯⋯000−⋯200𝑎₂⋯00⋯0000000000⋯⋯00000000⋯⋯0000000−200𝑎₂𝐴₁₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=⎡⎢⎢⎢⎣−0𝐹₁−𝐹𝑐²¹₁𝑙𝑙0⎤⎥⎥⎥⎦,𝐴¹²=⎡⎢⎢⎢⎣𝐹0₁𝐹𝑐²¹₁𝑙𝑙0⎤⎥⎥⎥⎦,𝐴²²=⎡⎢⎢⎢⎣−0𝐹₂−𝐹𝑐²²₂𝑙𝑙0⎤⎥⎥⎥⎦,𝐴²¹=⎡⎢⎢⎢⎣𝐹0₂𝐹𝑐²²₂𝑙𝑥=𝑃₁𝑙0⎤⎥⎥⎥⎦,𝑥₀,𝑄₁,𝑃₂,𝑄₂,…,𝑃_𝑁,𝑄_𝑁,…,𝑃_2𝑁,𝑄_2𝑁^,=𝑃⁰¹,𝑄⁰¹,𝑃⁰²,𝑄⁰²𝑢,=…,𝑃𝑃^0𝑁₀,𝑄^0𝑁,…,𝑃^02𝑁,𝑄^02𝑁^,𝐵=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣0𝑐₂₁𝑄₀,𝐹₁𝐹₁𝑙⋯⋯0000𝑙⋯⋯00000⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,𝐶=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⋯000⋯𝑐00²²𝐹₂𝐹₂𝑙⋯00𝑙⋯000⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦𝑃^𝑝𝑙𝑄_𝑝𝑙.(1) Let 𝑅=𝑅𝜀. Then the corresponding Euler-Lagrange control problem can be written as𝑥𝜆̇̇=𝐴⁰+0𝐴₁𝜀−−𝐴𝜀𝐵₁₀𝑅⁻¹𝐵′𝑥(0)=−𝑥𝐴¹¹₀𝜀𝑥𝜆+𝐶0𝜀,(2.12)𝜆(𝑇)=𝑁𝑥(𝑇),−𝑥.(2.13) 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 𝑥(0)and𝜆(𝑇). According to [12], 𝑥𝜆((𝑇𝑇))=𝑒^{𝐴0+𝜀𝐴1−𝜀𝐵𝑅−1𝐵′0−𝐴′0−𝜀𝐴′1}^𝑇+⎡⎢⎢⎣𝑒𝜆𝑥(⁰0)_{𝐴0+𝜀𝐴1−𝜀𝐵𝑅−1𝐵′0−𝐴′0−𝜀𝐴′1}^𝑇−𝐸0𝐸0⎤⎥⎥⎦𝐴⁰+𝜀𝐴₁−𝜀𝐵𝑅⁻¹𝐵^0−𝐴^0−𝜀𝐴^1⁻¹𝐶0𝜀.(3.1) Further, let us expand the expression𝑒^{𝐴0+𝜀𝐴1−𝜀𝐵𝑅−1𝐵′0−𝐴′0−𝜀𝐴′1}^𝑇=𝑒^{𝐴000−𝐴′0}^𝑇+𝜀^{𝐴1−𝐵𝑅−1𝐵′}_0−𝐴^′1^𝑇(3.2) from (3.1) with respect to 𝜀. If we denote 𝐻=₁^i.e𝐴⁰0−0𝐴_0,𝐻^2i.e=𝐴¹0−𝐵−𝑅𝐴⁻¹_1𝐵′,(3.3) then according to [12], the expansion of expression (3.2) with respect to 𝜀 can be represented as: 𝑒^{𝐻1𝑇+𝜀𝐻2𝑇}≈𝑒^𝐻1𝑇+𝜀¹⁰𝑒^{𝐻1𝑇(1−𝑠)}𝐻₂𝑇𝑒^{𝐻1𝑇𝑠}𝑑𝑠(3.4) Denote the integral in (3.4) by 𝐿₀. Then it is not difficult to show that matrix 𝐿₀ is a solution of the following Sylvester's equation:𝐻₁𝐿₀−𝐿₀𝐻₁=𝑒^𝐻1𝑇𝐻₂−𝐻₂𝑒^𝐻1𝑇.(3.5) Therefore, for expression (3.2) we obtain the expansion𝑒^{𝐻1𝑇+𝜀𝐻2𝑇}≈𝑒^𝐻1𝑇+𝜀𝐿₀.(3.6) Introducing notations 𝐿₀=𝐿¹𝐿₂𝐿₃𝐿₄,𝐴⁰+𝜀𝐴₁−𝜀𝐵𝑅⁻¹𝐵^0−𝐴^0−𝜀𝐴^1⁻¹=𝑆¹𝑆₂0𝑆₄,(3.7) expression (3.1) can be simplified to 𝑥𝜆((𝑇𝑇))=𝑒^𝐴0𝑇𝑒𝑥(0)_{−𝐴0𝑇}𝜆(0)+𝜀𝐿¹𝑥(0)+𝐿𝐿₂₃𝜆(0)+𝑒^𝐴0𝑇𝑆₁𝐶−𝑆₁𝐶𝑥(0)+𝐿₄𝜆(0).(3.8) Hence, adding boundary conditions from (2.12), we arrive to the following system of algebraic equations: 𝑥(𝑇)−𝜀𝐿𝜆₂𝜆(𝑇(0))−=𝑒𝑒^𝐴0𝑇𝑥_{−𝐴′0𝑇}(0)+𝜀𝑒^𝐴0𝑇𝑆₁𝐶−𝑆₁𝐶+𝐿₁𝑥(0),𝜆(𝑇)−+𝑁𝜀𝑥𝐿(₄𝑇)𝜆=(0−)𝑁=𝜀𝐿³𝑥(0),𝑥(3.9) which in the matrix form can be written as⎡⎢⎢⎣𝐸−𝜀𝐿²0−𝑒_{−𝐴′0𝑇}0𝑁0−𝜀𝐿₄−𝐸𝐸⎤⎥⎥⎦⎡⎢⎢⎣𝑥𝜆𝜆(((𝑇𝑇0)))⎤⎥⎥⎦=⎡⎢⎢⎣𝑒^𝐴0𝑇𝑥₀+𝜀𝑒^𝐴0𝑇𝜀𝐿𝑆₁₃𝐶−𝑆₁𝐶+𝐿₁𝑥₀𝑥₀𝑁𝑥⎤⎥⎥⎦.(3.10) Multiplying (3.10) by [^{𝐸000𝐸𝐸00𝐸}] on the left, we obtain⎡⎢⎢⎣𝐸−𝜀𝐿²𝑁−𝑒_{−𝐴′0𝑇}0𝑁0−𝜀𝐿₄−0𝐸⎤⎥⎥⎦⎡⎢⎢⎣𝑥𝜆𝜆(((𝑇𝑇0)))⎤⎥⎥⎦=⎡⎢⎢⎣𝑒^𝐴0𝑇𝑥₀+𝜀𝑁𝑒^𝐴0𝑇𝑆₁𝐶−𝑆₁𝐶+𝐿₁𝑥₀𝑥+𝜀𝐿₃𝑥₀𝑁𝑥⎤⎥⎥⎦.(3.11) Hence, if we denote the coefficient matrix in (3.11) by 𝑀, then it can be written as𝑀=⎡⎢⎢⎣𝐸−𝜀𝐿²𝑁−𝑒_{−𝐴′0𝑇}0𝑁0−𝜀𝐿₄−0𝐸⎤⎥⎥⎦=𝐾𝐾(𝜀)0₁−𝐸,(3.12) where𝐾(𝜀)=𝐸−𝜀𝐿²𝑁−𝑒^{−𝐴′0𝑇}−𝜀𝐿₄,𝐾¹=𝑁0.(3.13) It is not difficult to show that𝐾𝐾(𝜀)−0₁𝐸⁻¹=𝐾⁻¹𝐾(𝜀)0₁𝐾⁻¹(𝜀)−𝐸.(3.14) Hence, using the fact that𝐾⁻¹(𝜀)≈𝐾⁻¹(0)−𝐾⁻¹(0)𝐾̇(0)𝐾⁻¹(0)𝜀(3.15) we obtain the inverse matrix 𝑀⁻¹≈⎡⎢⎢⎢⎣𝐸+𝐿²𝑒𝑒^𝐴′0𝑇𝑁𝜀−𝐿₂𝑒^𝐴′0𝑇𝜀0_𝐴′0𝑇𝑁−𝑒^𝐴′0𝑇𝑁𝐿₄+𝑒𝑁^𝐴′0𝑇𝐿𝑁₂𝜀+𝑒^𝐴′0𝑇𝑁𝐿₂𝑒^𝐴′0𝑇𝑁𝜀𝑒^𝐴′0𝑇−𝑒^𝐴′0𝑇𝑁𝐿₂𝑒^𝐴′0𝑇𝜀+𝑒^𝐴′0𝑇𝐿₄𝜀0𝑒^𝐴′0𝑇𝑁𝜀−𝑁𝐿₂𝑒^𝐴′0𝑇𝜀−𝐸⎤⎥⎥⎥⎦(3.16) and, consequently, multiplying the both sides of (3.11) by 𝑀⁻¹ on the left, we obtain the following analytic formulae to determine values of 𝑥(𝑇),𝜆(0),𝜆(𝑇): ⎡⎢⎢⎣𝑥𝜆𝜆(((𝑇𝑇0)))⎤⎥⎥⎦=⎡⎢⎢⎢⎣𝑒^𝐴0𝑇𝑥₀𝑒^𝐴′0𝑇𝑁𝑁𝑒𝑒^𝐴0𝑇_𝐴0𝑇𝑥₀+𝑒^𝐴′0𝑇𝑁𝑥+𝜀⎡⎢⎢⎢⎢⎢⎢⎣𝑒𝑥_0𝐴0𝑇−𝑁𝑥⎤⎥⎥⎥⎦⎧⎪⎨⎪⎩𝑒_{−𝐴′0𝑇}𝑆₁𝐶−𝑆₁𝐶+𝐿₁𝑥₀+𝐿₂𝑒^𝐴′0𝑇𝑁𝑒^𝐴0𝑇𝑥₀−𝐿₂𝑒^𝐴′0𝑇𝑁𝑥𝑁𝑒^𝐴0𝑇𝑆₁𝐶−𝑒^𝐴′0𝑇𝑁𝑆₁𝐶+𝑒^𝐴′0𝑇𝑁𝐿₁𝑥₀+𝑒^𝐴′0𝑇𝐿₃𝑥₀−𝑒^𝐴′0𝑇𝐿₄𝑒^𝐴′0𝑇𝑁𝑒^𝐴0𝑇𝑥₀+𝑒^𝐴′0𝑇𝑁𝐿₂𝑒^𝐴′0𝑇𝑁𝑒^𝐴0𝑇𝑥₀𝑁𝐿₂−𝑒^𝐴′0𝑇𝑁𝐿₂𝑒^𝐴′0𝑇𝑁𝑥+𝑒^𝐴0𝑇𝐿₄𝑒^𝐴′0𝑇𝑁𝑥⎫⎪⎬⎪⎭𝑒^𝐴′0𝑇𝑁𝑒^𝐴0𝑇𝑥₀−𝑁𝐿₂𝑒^𝐴′0𝑇𝑁𝑥+𝑁𝐿₁𝑥₀+𝑁𝑒^𝐴0𝑇𝑆₁𝐶−𝑁𝑆₁𝐶⎤⎥⎥⎥⎥⎥⎥⎦.(3.17) Further, as in (3.1), using determined in (3.17) value 𝜆(0), it is possible to find 𝑥(𝑡_𝑖),𝜆(𝑡_𝑖) in the form𝑥𝑡^𝑖𝜆𝑡_𝑖=𝑒^{𝐴0+𝜀𝐴1−𝜀𝐵𝑅−1𝐵′0−𝐴′0−𝜀𝐴′1}^𝑡𝑖+⎡⎢⎢⎣𝑒𝜆𝑥(⁰0)_{𝐴0+𝜀𝐴1−𝜀𝐵𝑅−1𝐵′0−𝐴′0−𝜀𝐴′1}−𝐸0𝐸0⎤⎥⎥⎦𝐴₀+𝜀𝐴₁−𝜀𝐵𝑅⁻¹𝐵^0−𝐴₀−𝜀𝐴₁⁻¹𝐶0𝜀(3.18) for every 𝑡_𝑖∈[0𝑇].

Therefore, expansions for 𝑥(𝑡_𝑖)and 𝑢(𝑡_𝑖) with respect to 𝜀 can be obtained in the form𝑥𝜆𝑡𝑡_𝑖𝑖=𝑒^{𝐴0𝑡𝑖}+𝜀𝑒^{𝐴0𝑡𝑖}𝑆₁𝐶−𝑆₁𝐶+𝐿₁₁𝑥₀−𝐿₂₂𝑒^{𝐴0𝑡𝑖}𝑁𝑒^{𝐴0𝑡𝑖}𝑥₀+𝐿₂₂𝑒^{𝐴0𝑡𝑖}𝑁𝑥,=+−𝑁𝜀𝑒𝑁𝐿^𝐴0𝑇𝑥₀₂₂+𝑁𝑥𝑢𝑡_𝑖𝑒^{𝐴0𝑡𝑖}𝑁𝑒^{𝐴0𝑡𝑖}𝑥₀−𝑁𝐿₂₂𝑒^{𝐴0𝑡𝑖}𝑁𝑥−𝑁𝐿₁₁𝑥₀−𝑁𝑒^{𝐴0𝑡𝑖}𝑆₁𝐶−𝑁𝑆₁𝐶,=−𝑅⁻¹𝐵′𝜆𝑡_𝑖.(3.19) Here 𝐿₁₁,𝐿₂₂ are block matrices of matrix-solutions 𝐿_𝑖 of Sylvester's equation [12, 13]𝐻₁𝐿_𝑖−𝐿_𝑖𝐻₁=𝑒^{𝐻1𝑡𝑖}𝐻₂−𝐻₂𝑒^{𝐻1𝑡𝑖}(3.20) for every 𝑡_𝑖.

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 𝐶=0,𝑁=2 is developed. It is clear that for 𝑁>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 when𝐶=0,𝑁=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 Figures 1, 2, 3, and 4 the graphs of the functions 𝑃₁,𝑄₁,𝑃₂ and𝑄₂ 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 functions𝑃₁,𝑄₁,𝑃₂, and 𝑄₂ for both algorithms coincide with a sufficiently high accuracy.

(a)

(b)

Figure 1 

Dependence of $P_1$ on t.

(a)

(b)

Figure 2 

Dependence of $P_2$ on t.

(a)

(b)

Figure 3 

Dependence of $Q_1$ on t.

(a)

(b)

Figure 4 

Dependence of $Q_2$ on $t$ .