Moving Horizon Estimation in Managed Pressure ... - Semantic Scholar
Moving Horizon Estimation in Managed Pressure Drilling using Distributed Models Agus Hasan and Lars Imsland Abstract— This paper presents a moving horizon estimation (MHE) approach for estimating flow and pressure inside the annulus during managed pressure drilling (MPD). MPD is an advanced pressure control method that is used to precisely control the annular pressure throughout the wellbore in oil well drilling. The MPD model, which is a hyperbolic-type partial differential equation (PDE), is a hydraulic model based on basic fluid dynamics. Using an early lumping approach, the model yields a high-dimensional system of an ordinary differential equation (ODE). We show the adjoint-based method combined with a line search algorithm can solve the MHE dynamic optimization problem efficiently. The estimator is tested using data from a field scale flow-loop test conducted in Stavanger, Norway by the oil company Statoil.
I. INTRODUCTION A. Process Description During oil well drilling, a carefully designed fluid (often called mud) is pumped down from the mud pit through the drill string, through the drill bit, up the annulus around the drill string, and back to the mud pit. The aim is not only to transport cuttings in the annulus, but also to manage the pressure in the well (especially the annular downhole pressure pdh ) so that the fracturing of the formation or collapse of the well can be avoided [1]. In conventional drilling, the aim is achieved by changing the mud density. This is a nontrivial task that may lead to delay and the pressure cannot be precisely controlled. In managed pressure drilling (MPD), conversely, a rotating control device seals the top of the annulus and the mud flow from the well is controlled by the rig pump and choke manifold (Fig. 1). Adjusting the choke opening is much easier and faster than changing the mud density. Hence, it allows for the precise control of pdh . Three pressure limits need to be considered during oil well drilling (order from the lowest to highest); collapse pressure, formation pressure, and fracture pressure. In both conventional drilling and MPD, the pressure reference at the drill bit is set above the formation pressure and below the fracture pressure. Therefore, fluid outflux into the formation known as lost circulation, occurs during MPD operations. If for some reasons, e.g., the rig pump loss of power, and if the formation is permeable, the reference pressure may fall below the pore pressure and the fluid influx from the formation known as kick may occur. Without proper maintenance, either kick or lost circulation can cause loss This work was supported by The Research Council of Norway and Statoil ASA. The authors are with Department of Cybernetics Engineering, Norwegian University of Science and Technology, 7034 Trondheim, Norway.
[email protected]
of well control (blowout), loss of assets, and even loss of life. The implication of such undesirable conditions are not only environmental, but it also concerns worker’s safety [2]. Therefore, monitoring and control are important issues in oil well drilling. Another variant in oil well drilling is by drilling the formation with pressure reference between the collapse pressure and formation pressure, known as under-balanced drilling (UBD). Therefore, the reservoir fluid flows into the wellbore and up to the surface. This makes the well cleaner and the productivity index higher. To ensure the operation runs safely, the fluid is controlled by adding additional equipment at the surface. Consequently, the number of workers involved, the operational time, and the operational costs also increased. This drawback makes UBD rarely used [3].
Fig. 1: Schematics of an automated MPD system [4].
B. Motivation and Objective of this Paper The monitoring system for all drilling techniques is usually relies only on the topside flow and pressure measurements, whilst the main interest is to know the pressure and the flow inside the wellbore especially at the bottom of the well. Hence, a monitoring and controlling system for the oil well drilling must be based on a model. Furthermore, the three pressure limits needs to be considered in the entire open hole; if the open hole is long, one may simultaneously be close to the fracture pressure at the casing shoe (top of open hole) and close to the formation pressure in the bottom. This motivates a distributed approach to pressure estimation.
The objective of this paper is to estimate the flow and pressure inside the wellbore in MPD using only topside measurements. To this end, a moving horizon estimation (MHE) approach for distributed models is employed where the dynamic optimization problem is solved efficiently by combining an adjoint-based method to compute the gradient and a line search algorithm to find the optimal solution. C. Literature Review Flow and pressure estimations during oil well drilling have been subject to many researches, e.g., [4], [5], [6]. Usually a model of the system is used together with inputs and outputs to calculate the estimates of the unmeasured flow, pressure and unknown or uncertain parameters like the frictional coefficient. The most known method is the Kalman filter (KF). The KF will provide an optimal estimate of the states in a linear system driven by white gaussian noise with known covariance. Some modified versions of the KF have been developed for nonlinear (Extended KF and Unscented KF) and large-scale systems (Ensemble KF). Although already widely used in industry, the KF-based method can fail to estimate the state even in the absence of plant-model mismatch. It has been shown MHE consistently provides improved state estimation and greater robustness to both poor guesses of the initial state and tuning parameters in comparison to the EKF [7]. The only drawback of MHE is the greater computational expense required to solve the MHE optimization [8]. More recent state and parameter estimations related to this paper were presented in [9], [11], [12]. However, only pressure at the bottom is estimated using low order lumped models. Thus, these methods cannot be used directly to determine if the pressure limits inside the annulus are satisfied. D. Contribution of this Paper The paper presents MHE approach for state and parameter estimations for a class of 2×2 linear hyperbolic systems with application to MPD. The optimization for MHE is solved efficiently using an adjoint-based method and the estimations converge to the actual value. Furthermore the design is tested using data from a field-scale flow loop test conducted in Stavanger, Norway by the oil company Statoil. E. Organization of this Paper The paper is organized as follows. The model for MPD based on the hydraulic model is presented in section II. In this section, the model is transformed into a hyperbolictype PDE and is simplified using an early lumping approach which yield a high-dimensional ODE system. Some basic concepts regarding MHE and how to solve the optimization problem is highlighted in section III. The experimental design is described in section IV, while the simulation results are presented in section V. The last section presents the conclusions and recommendations for future work.
II. MATHEMATICAL MODEL A. The Hydraulic Model The aim of the hydraulic model is to be able to estimate the down-hole pressure profile and provide a choke pressure set point to the MPD control system in real time. Another aim of the hydraulic model is control system design; however, this will not be discussed in this paper. To ensure robustness of the resulting control system, the model should not be more complex than the requirements on the control system dictate. Some simplification can be made for instance by lumping together parameters that are not possible to distinguish or calibrate independently from existing measurements [1]. The main assumption when deriving the hydraulic model is to consider the drilling fluid as a viscous fluid, i.e., the flow is completely described by fundamental equations such as equation of state, conservation of mass, conservation of momentum, and conservation of energy. The derivation of the hydraulic model can be seen in [13]. Here, we only give the result. The hydraulic model for pressure and flow inside the annulus can be written as follows β (1) pt (z, t) = − qz (z, t) Aa Fc (z, t) Aa + Aa g cos(η) (2) qt (z, t) = − pz (z, t) − ρ ρ q(0, t) = qp (t) (3) p(l, t)
= pc (t)
(4)
where the annular pressure p and volumetric flow rate q are the distributed states over (z, t) ∈ [0, l] × [0, ∞). The isothermal bulk modulus is denoted by β, and the frictional pressure drop is denoted by Fc . Furthermore, Aa , ρ, and g denote the area of the annulus, the fluid density, and the gravitational constant, respectively. η is defined as the angle between gravity and the positive flow direction, e.g., for a vertical well η = 180. The volumetric flow rate of the injected fluid from the main pump is denoted by qp . For simplicity, we assume that the choke is operated by a fast pressure controller so that we can treat pc rather than the choke opening as the control input. Furthermore, the model neglects dependence on temperature. Although temperature gradient does exist, the thermal liquid expansion coefficient is usually small. Therefore, the density changes due to temperature changes are negligible. An important issue in the hydraulic model is the frictional pressure drop Fc . The frictional pressure drop in the annulus is typically linear with respect to the flow, see e.g., [14]. Hence, the model for frictional pressure drop as a function of flow rate is given by Fc (z, t) = F × q(z, t)
(5)
where F is the frictional coefficient. In the flow-loop experiment, for a flow below 1000 liters per minute, the model agrees well with the experiment. The comparison is presented in Fig. 2. Since the frictional coefficient is an uncertain parameter of the annulus, there is an incentive to calibrate the hydraulic model with respect to this parameter. In section
The last step is to define
Pressure drop [bara]
2
u(x, t) = u ¯(xl, t)eax v(x, t) = v¯(xl, t)e−ax
1.5
1
where
Modeled Measured
Fl a= √ 2 βρ
0.5
0 650
700
750
800
850
900
950
1000
1050
Volumetric flow rate [lpm]
Fig. 2: Frictional pressure drop in the annulus.
III, we formulate the MHE for state and parameter estimation of the hydraulic model (1)-(4).
Plugging (21)-(22) into (17)-(20), yield s 1 β F ut (x, t) = − ux (x, t) − e2ax v(x, t) l ρ 2ρ s F 1 β vx (x, t) − e−2ax u(x, t) vt (x, t) = l ρ 2ρ u(0, t)
=
v(1, t)
=
B. A Canonical Form of the Hydraulic Model The hydraulic model (1)-(4) can be written into the following canonical form of a hyperbolic-type PDE ut (x, t)
= −1 (x)ux (x, t) + c1 (x)v(x, t)
(6)
vt (x, t)
= 2 (x)vx (x, t) + c2 (x)u(x, t)
(7)
u(0, t)
= qv(0, t) + U1 (t)
(8)
v(1, t)
= U2 (t)
(9)
First, the hydrostatic pressure due to gravity from the momentum equation can be eliminated by using the following change of variable p¯(z, t)
= p(z, t) − ρgl − ρgz cos(η)
β qz (z, t) Aa F Aa qt (z, t) = − p¯z (z, t) − q(z, t) ρ ρ q(0, t) = qp (t) p¯(l, t)
= −
= pc (t) − ρgl(1 + cos(η))
(24) (25)
−v(0, t) + qp (t) (26) −a A e a (pc (t) − ρgl(1 + cos η)) u(1, t)e−2a − √ βρ (27)
The above system resemble (6)-(9). The original state can be computed using these formulas z a z a q(z, t) = u (28) , t e− l z + v ,t elz l √ l a a βρ z z p(z, t) = , t e− l z − v ,t elz u Aa l l +ρgl + ρgz cos(η) (29) The distributed model (6)-(9) can be simplified into a lumped model using an early lumping approach as we describe in the following section.
(11)
C. Early Lumping Approach
(12) (13) (14)
Next, using the Riemann coordinate transformations [15] Aa 1 (15) q(z, t) + √ p¯(z, t) u ¯(z, t) = 2 βρ 1 Aa v¯(z, t) = q(z, t) − √ p¯(z, t) (16) 2 βρ (11)-(14) are transformed into s β F u ¯t (z, t) = − u ¯z (z, t) − (¯ u(z, t) + v¯(z, t)) (17) ρ 2ρ s β F v¯t (z, t) = v¯z (z, t) − (¯ u(z, t) + v¯(z, t)) (18) ρ 2ρ = −¯ v (0, t) + qp (t) (19) Aa v¯(l, t) = u ¯(l, t) − √ (pc (t) − ρgl(1 + cos(η))) (20) βρ
u ¯(0, t)
(23)
(10)
Plugging (10) into (1)-(4), yields p¯t (z, t)
(21) (22)
The idea of early lumping approach is to replace the infinite-dimensional PDE (6)-(9) with a finite set of ODE using the method of lines. Here, the spatial coordinate x is discretized using a non-uniform grid with N -nodes where the variable ui and vi for all i ∈ {1, · · · , N } denote the evaluation of u and v at the grid nodes. Hence, for i = 1, · · · , N , equations (6)-(7) can be written as dui dt dvi dt
= =
dui + c1 (xi )vi dxi dvi 2 (xi ) + c2 (xi )ui dxi
−1 (xi )
(30) (31)
The spatial derivative is discretized using an upwind finite difference scheme as follow dui dxi dvi dxi
= =
ui − ui−1 δx vi+1 − vi δx
(32) (33)
where δx = N1+1 . Remark that the discretization schemes for u and v are slightly different due to boundary conditions.
Plugging (32)-(33) to (30)-(31), yield du1 dt du2 dt duN dt dv1 dt dv2 dt
III. MOVING HORIZON ESTIMATION DESIGN
u1 − u0 + c1 (x1 )v1 δx u2 − u1 = −1 (x2 ) + c1 (x2 )v2 δx .. . = −1 (x1 )
(34) (35) (36)
uN − uN −1 = −1 (xN ) + c1 (xN )vN δx v2 − v1 = 2 (x1 ) + c2 (x1 )u1 δx v3 − v2 = 2 (x2 ) + c2 (x2 )u2 δx .. .
(37) (38)
ˆ ∆y = y − C x
(40)
x˙ = Ax + BU
(42)
[u1 · · · uN v1 · · · vN ]T
The state matrix is given by E1 A= C2
C1 E2
(43)
C2
E1
0
···
E2
=
tk −τ
ˆ k −τ,θ
ˆ˙ = Aˆ ˆ (tk − τ ) = x ˆ tk −τ (53) x x + Bu, x
−2 (x1 ) 2 (x1 ) · · · −2 (x2 ) · · · 1 0 .. .. .. δx . . . 0 0 ···
The input matrix is given by 1 1 (x1 ) 0 B= 0 −2 (xN ) δx while the control vector is given by U1 U= U2
with the positive semi-definite matrix Q ∈ Rp×p weighting the estimation error (51). Remark that the optimization ˆ tk −τ and variables for solving (52) are the initial value x the parameter θˆ (in this case, the frictional coefficient F ). B. Gradient Computation
1 (x1 ) (x1 ) −q 1δx 2q δx1 1 = diag(c1 ) + 0 0 = diag(c2 ) −1 (x1 ) 0 ··· (x ) − (x ) · ·· 1 2 1 1 2 = .. .. .. δx . . . 0
ˆt x
(44)
where C1
(51)
during the preceding time interval t ∈ [tk − τ, tk ] with the estimation horizon t > 0. Here, y contains p measurements and C ∈ Rp×2N denotes matrix which appropriately interˆ at the measurement locations. polates the estimator state x The dynamic optimization problem for estimating the current ˆ k can be formulated as follow state x Z tk J(ˆ xtk −τ,θˆ) = min ∆y | Q∆y dt (52) subject to
where the state vector is given by =
The idea of MHE is to estimate the current states by solving a least squares dynamic optimization problem, which penalizes the deviation between the measurements and the predicted outputs, with the initial value of the estimator system and/or the parameters of the system considered as the optimization variables [16]. Mathematically, we want to minimize the estimation error
(39)
U2 (t) − vN dvN = −2 (xN ) + c2 (xN )uN (41) dt δx The left boundary condition can be computed by interpolating v0 ≈ 2v1 − v2 ; hence, u0 = q(2v1 − v2 ) + U1 . The ODEs (34)-(41) can be written in a state-space representation as
x
A. Formulation of Dynamic Optimization
0 0
(45)
0 0 .. .
(46)
(47) −1 (xN ) 0 0 (48) .. . −2 (xN )
The dynamic optimization problem (52)-(53) can be solved using a gradient-based method. It has been known, e.g. see [17], [18], that the gradients of the cost function with respect to the initial condition and the parameter θˆ are given by ˆ ∂J(ˆ xtk −τ , θ) ˆ tk −τ ∂x ˆ ∂J(ˆ xt −τ , θ) k
∂ θˆ
= λ(tk − τ ) =
ˆ ∂H(ˆ xtk −τ , λ, θ) ∂ θˆ
(55)
where the adjoint state λ(tk − τ ) follows from solving the following costate equations ˆ˙ x
(49)
(54)
= Aˆ x + Bu ˆ (tk − τ ) = x ˆ tk −τ x ˆ ˆ x, λ, θ) ∂ H(ˆ λ˙ = − ˆ ∂x λ(tk ) = 0
(56) (57) (58) (59)
with the Hamiltonian is given by (50)
To enable online calibration of uncertainties in the friction characteristic, the frictional coefficient F is considered as a tuning parameter for the model. We highlighted the MHE method for state and parameter estimation in the following section.
ˆ = ∆y | Q∆y + λ| (Aˆ H(ˆ x, λ, θ) x + Bu)
(60)
(0) ˆ tk −ˆτ and θˆ(0) , one Starting with an initial estimate of x gradient iteration j consists of ˆ (j) , • integrating (56)-(57) forward in time to obtain x (j) • integrating (58)-(59) backward in time to obtain λ ,
Parameter qp ρ l Aa β
Value 1/60 1000 700 154.8 2×10−9
Unit m3 /s kg/m3 m mm Pa
TABLE I: Experimental data
•
(j)
g1
(j)
g2 •
(a) The main pump.
(b) The pipes.
(c) Downhole sensors.
(d) Topside sensors.
computing the gradient =
λ(j) (tk − τ )
(61)
=
ˆ ∂H(ˆ xtk −τ , λ, θ) ∂ θˆ
(62)
updating the estimated state (j+1)
¯ tk −τ x θˆ(j+1)
= =
(j)
(j)
ˆ tk −τ − α1 g1 x (j) θˆ(j) − α2 g 2
(63) (64)
where α1 and α2 are the line search parameters. The gradient ˆ iterations, which leads to the algorithm is stopped after N estimate of the current state ˆ
ˆk = x ˆ N (tk ) x
Fig. 4: Experiment setup
(65)
IV. EXPERIMENT DESIGN The pipes are placed nearly horizontally, with the elevation observed in Fig. 5. Therefore, the hydrostatic pressure has a relatively small effect, and the pressure difference is primarily due to the frictional pressure. 60 50
Vertical Depth [m]
The experiment is carried out in a field scale flow-loop test in Stavanger, Norway. The MPD system is configured as a U-tube (Fig. 3), which consists of the main pump (Fig. 4a), 700 meters of pipes (each for drill string and annulus) (Fig. 4b), the bottom hole assembly (Fig. 4c), and the top side sensors (Fig. 4d). The data used in this experiment are given in Table 1.
40 30 20 10 0 0
100
200
300
400
500
600
700
Well Path [m]
Fig. 5: Loop profile in XZ-coordinate. V. SIMULATION RESULTS
Fig. 3: Schematics of the flow-loop. In this experiment, water is injected by the main pump through the drill string and up the annulus. The water rate is changed so that the dynamic of the model can be observed. The aim is to estimate the current pressure and flow throughout the annulus using previous measurement from the top of the annulus.
The spatial discretization is realized with N = 10 grid nodes. The sampling time and prediction horizon are chosen as ∆t = 1 s and τ = t = 10 s, respectively. Since the measurement taken only at the top of the well, then p = 2 sensors are considered. Hence, the measurement vector is u|x=xN given as y = . The hydraulic model is solved v|x=xN using Matlab solver ode23, while the step sizes α1 and α2 are obtained from an adaptive line search algorithm [10]. Implementing the MHE algorithm described in section III, the estimation result for the pressure and flow inside the well can be found in Fig. 6, while the estimated frictional
coefficient can be seen in Fig. 7. It takes less than 1 s to run a one step MHE, which shows the real-time feasibility of the MHE scheme. 800
Estimated Measured
Volumetric flow rate [lpm]
Volumetric flow rate [lpm]
800
750
700
650
600 0
50
100 150 Time [s]
200
250
(a) Flow rate through bit.
750
700
650
R EFERENCES 600 0
25
20
50
100 150 Time [s]
200
250
(c) Pressure at bit.
100 150 Time [s]
30
Estimated Measured
15 0
50
200
250
(b) Flow rate through wellhead.
Pressure [bara]
Pressure [bara]
30
Estimated Measured
Estimated Measured
25
20
15 0
50
100 150 Time [s]
200
250
(d) Pressure at wellhead.
Fig. 6: Pressure and flow at the bottom (bit) and topside (wellhead) of the well.
93
Coefficient of friction
92.5 92 91.5 91 90.5 90 0
50
ACKNOWLEDGMENT Financial support from Statoil ASA and the Norwegian Research Council (NFR project 210432/E30 Intelligent Drilling) is gratefully acknowledged. The experimental data are results from ongoing internal research in the department Intelligent Drilling in Statoil RDI. The authors gratefully acknowledge the work done by Alexey Pavlov, Henrik Manum and Glenn-Ole Kaasa in planning and conducting experiments.
100
150
200
250
Time [s]
Fig. 7: Coefficient of friction as tuning parameter. We can see a good performance of MHE for flow, pressure, and friction characteristic estimations. VI. CONCLUSIONS AND RECOMMENDATIONS We have implemented MHE for estimating annular pressure and flow in MPD where the coefficient of friction serves as the tuning parameter for the model. The dynamic optimization problem is solved efficiently using a gradientbased method combined with a line search method. Because only one phase (fluid) is involved, the annular model is represented by a one-phase incompressible flow model. There is an interest to include the gas phase, e.g., when drilling a gas reservoir, leading to a two-phase flow model (refers to the drift flux model). This is a more challenging problem that should be addressed in the future research.
[1] J. M. Godhavn, Control Requirements for Automatic Managed Pressure Drilling System, SPE Drilling and Completion, vol. 25, no. 3, pp. 336–345, 2010. [2] M. Davoudi, J. R. Smith, J. E. Chirinos, and B. M. Patel, Evaluation of Alternative Initial Responses to Kicks Taken During ManagedPressure Drilling , SPE Drilling and Completion, vol. 26, no. 2, pp. 169–181, 2011. [3] J. E. Olsen, E. Vollen, and T. Tonnessen, Challenges in Implementing UBO Technology, SPE/IADC Underbalanced Technology Conference and Exhibition, 11-12 October, 2004, Houston, USA. [4] G. O. Kaasa, O. Y. Stamnes, O. M. Aamo, and L. Imsland, Simplified Hydraulics Model Used for Intelligent Estimation of Downhole Pressure for a Managed-Pressure-Drilling Control System, SPE Drilling and Completion, vol. 27, no. 1, pp. 127–138, 2012. [5] K. S. Bjorkevoll, R. Rommetveit, B. Aas, H. Gjeraldstveit, and A. Merlo, Transient Gel Breaking Model for Critical Wells Applications with Field Data Verification, SPE/IADC Drilling Conference, 19-21 February, 2003, Amsterdam, The Nedherland. [6] J. Petersen, R. Rommetveit, K. S. Bjorkevoll, and J. Froyen, A General Dynamic Model for Single and Multi-phase Flow Operations during Drilling, Completion, Well Control and Intervention, IADC/SPE Asia Pacific Drilling Technology Conference and Exhibition, 25-27 November, 2008, Jakarta, Indonesia. [7] E. L. Haseltine and J. B. Rawlings, Critical Evaluation of Extended Kalman Filtering and Moving-Horizon Estimation, Industrial and Engineering Chemistry Research, vol. 44, no. 8, pp. 2451–2460, 2005. [8] P. Kuhl, M. Diehl, T. Kraus, J. P. Schloder, and H. G. Bock, A Realtime Algorithm for Moving Horizon State and Parameter Estimation, Computer and Chemical Engineering, vol. 35, no. 1, pp. 71–83, 2011. [9] O. N. Stamnes, O. M. Kaasa, and G. O. Kaasa, Redesign of Adaptive Observers for Improved Parameter Identification in Nonlinear Systems, Automatica, vol. 47, no. 2, pp. 403–410, 2011. [10] J. Nocedal and S. Wright, Numerical Optimization, Springer, 2006. [11] M. Paasche, T. A. Johansen, and L. Imsland, Regularized and Adaptive Nonlinear Moving Horizon Estimation of Bottomhole Pressure during Oil Well Drilling, IFAC World Congress, 28 August-2 September, 2011, Milan, Italy. [12] D. Sui, R. Nybo, S. Hovland, and T. A. Johansen, A Moving Horizon Observer for Estimation of Bottomhole Pressure During Drilling, IFAC Workshop on Automatic Control in Offshore Oil and Gas Production, 31 May-1 June, 2012, Trondheim, Norway. [13] F. White, Fluid Mechanics, McGraw-Hill, New York, 2007. [14] I. S. Landet, A. Pavlov, O. M. Aamo, Modeling and Control of HeaveInduced Pressure Fluctuations in Managed Pressure Drilling, IEEE Transaction on Control System Technology, vol. 21, no. 4, pp. 1340– 1351, 2013. [15] G. Bastin and J. M. Coron, Further Result on Boundary Feedback Stabilization of 2×2 Hyperbolic System over a Bounded Interval, IFAC Symposium on Nonlinear Control System, 1-3 September, 2010, Bologna, Italy. [16] H. Michalska and D. Q. Mayne, Moving Horizon Observers and Observer-based Control, IEEE Transaction on Automatic Control, vol. 40, no. 6, pp. 995–1006, 1995. [17] Y. Cao, S. Li, L. Petzold, and R. Serban, Adjoint Sensitivity Analysis for Differential-Algebraic Equations: The Adjoint DAE System and its Numerical Solution, SIAM Journal on Scientific Computing, vol. 24, no. 3, pp. 1076–1089, 2003. [18] S. Rhein, T. Utz, and K. Graichen, Model Predictive Control and Moving Horizon Estimation of a Large-Scale Chemical Reactor Model, IFAC Workshop on Control of Systems Modeled by Partial Differential Equations, 25-27 September, 2013, Paris, France.
I. INTRODUCTION A. Process Description During oil well drilling, a carefully designed fluid (often called mud) is pumped down from the mud pit through the drill string, through the drill bit, up the annulus around the drill string, and back to the mud pit. The aim is not only to transport cuttings in the annulus, but also to manage the pressure in the well (especially the annular downhole pressure pdh ) so that the fracturing of the formation or collapse of the well can be avoided [1]. In conventional drilling, the aim is achieved by changing the mud density. This is a nontrivial task that may lead to delay and the pressure cannot be precisely controlled. In managed pressure drilling (MPD), conversely, a rotating control device seals the top of the annulus and the mud flow from the well is controlled by the rig pump and choke manifold (Fig. 1). Adjusting the choke opening is much easier and faster than changing the mud density. Hence, it allows for the precise control of pdh . Three pressure limits need to be considered during oil well drilling (order from the lowest to highest); collapse pressure, formation pressure, and fracture pressure. In both conventional drilling and MPD, the pressure reference at the drill bit is set above the formation pressure and below the fracture pressure. Therefore, fluid outflux into the formation known as lost circulation, occurs during MPD operations. If for some reasons, e.g., the rig pump loss of power, and if the formation is permeable, the reference pressure may fall below the pore pressure and the fluid influx from the formation known as kick may occur. Without proper maintenance, either kick or lost circulation can cause loss This work was supported by The Research Council of Norway and Statoil ASA. The authors are with Department of Cybernetics Engineering, Norwegian University of Science and Technology, 7034 Trondheim, Norway.
[email protected]
of well control (blowout), loss of assets, and even loss of life. The implication of such undesirable conditions are not only environmental, but it also concerns worker’s safety [2]. Therefore, monitoring and control are important issues in oil well drilling. Another variant in oil well drilling is by drilling the formation with pressure reference between the collapse pressure and formation pressure, known as under-balanced drilling (UBD). Therefore, the reservoir fluid flows into the wellbore and up to the surface. This makes the well cleaner and the productivity index higher. To ensure the operation runs safely, the fluid is controlled by adding additional equipment at the surface. Consequently, the number of workers involved, the operational time, and the operational costs also increased. This drawback makes UBD rarely used [3].
Fig. 1: Schematics of an automated MPD system [4].
B. Motivation and Objective of this Paper The monitoring system for all drilling techniques is usually relies only on the topside flow and pressure measurements, whilst the main interest is to know the pressure and the flow inside the wellbore especially at the bottom of the well. Hence, a monitoring and controlling system for the oil well drilling must be based on a model. Furthermore, the three pressure limits needs to be considered in the entire open hole; if the open hole is long, one may simultaneously be close to the fracture pressure at the casing shoe (top of open hole) and close to the formation pressure in the bottom. This motivates a distributed approach to pressure estimation.
The objective of this paper is to estimate the flow and pressure inside the wellbore in MPD using only topside measurements. To this end, a moving horizon estimation (MHE) approach for distributed models is employed where the dynamic optimization problem is solved efficiently by combining an adjoint-based method to compute the gradient and a line search algorithm to find the optimal solution. C. Literature Review Flow and pressure estimations during oil well drilling have been subject to many researches, e.g., [4], [5], [6]. Usually a model of the system is used together with inputs and outputs to calculate the estimates of the unmeasured flow, pressure and unknown or uncertain parameters like the frictional coefficient. The most known method is the Kalman filter (KF). The KF will provide an optimal estimate of the states in a linear system driven by white gaussian noise with known covariance. Some modified versions of the KF have been developed for nonlinear (Extended KF and Unscented KF) and large-scale systems (Ensemble KF). Although already widely used in industry, the KF-based method can fail to estimate the state even in the absence of plant-model mismatch. It has been shown MHE consistently provides improved state estimation and greater robustness to both poor guesses of the initial state and tuning parameters in comparison to the EKF [7]. The only drawback of MHE is the greater computational expense required to solve the MHE optimization [8]. More recent state and parameter estimations related to this paper were presented in [9], [11], [12]. However, only pressure at the bottom is estimated using low order lumped models. Thus, these methods cannot be used directly to determine if the pressure limits inside the annulus are satisfied. D. Contribution of this Paper The paper presents MHE approach for state and parameter estimations for a class of 2×2 linear hyperbolic systems with application to MPD. The optimization for MHE is solved efficiently using an adjoint-based method and the estimations converge to the actual value. Furthermore the design is tested using data from a field-scale flow loop test conducted in Stavanger, Norway by the oil company Statoil. E. Organization of this Paper The paper is organized as follows. The model for MPD based on the hydraulic model is presented in section II. In this section, the model is transformed into a hyperbolictype PDE and is simplified using an early lumping approach which yield a high-dimensional ODE system. Some basic concepts regarding MHE and how to solve the optimization problem is highlighted in section III. The experimental design is described in section IV, while the simulation results are presented in section V. The last section presents the conclusions and recommendations for future work.
II. MATHEMATICAL MODEL A. The Hydraulic Model The aim of the hydraulic model is to be able to estimate the down-hole pressure profile and provide a choke pressure set point to the MPD control system in real time. Another aim of the hydraulic model is control system design; however, this will not be discussed in this paper. To ensure robustness of the resulting control system, the model should not be more complex than the requirements on the control system dictate. Some simplification can be made for instance by lumping together parameters that are not possible to distinguish or calibrate independently from existing measurements [1]. The main assumption when deriving the hydraulic model is to consider the drilling fluid as a viscous fluid, i.e., the flow is completely described by fundamental equations such as equation of state, conservation of mass, conservation of momentum, and conservation of energy. The derivation of the hydraulic model can be seen in [13]. Here, we only give the result. The hydraulic model for pressure and flow inside the annulus can be written as follows β (1) pt (z, t) = − qz (z, t) Aa Fc (z, t) Aa + Aa g cos(η) (2) qt (z, t) = − pz (z, t) − ρ ρ q(0, t) = qp (t) (3) p(l, t)
= pc (t)
(4)
where the annular pressure p and volumetric flow rate q are the distributed states over (z, t) ∈ [0, l] × [0, ∞). The isothermal bulk modulus is denoted by β, and the frictional pressure drop is denoted by Fc . Furthermore, Aa , ρ, and g denote the area of the annulus, the fluid density, and the gravitational constant, respectively. η is defined as the angle between gravity and the positive flow direction, e.g., for a vertical well η = 180. The volumetric flow rate of the injected fluid from the main pump is denoted by qp . For simplicity, we assume that the choke is operated by a fast pressure controller so that we can treat pc rather than the choke opening as the control input. Furthermore, the model neglects dependence on temperature. Although temperature gradient does exist, the thermal liquid expansion coefficient is usually small. Therefore, the density changes due to temperature changes are negligible. An important issue in the hydraulic model is the frictional pressure drop Fc . The frictional pressure drop in the annulus is typically linear with respect to the flow, see e.g., [14]. Hence, the model for frictional pressure drop as a function of flow rate is given by Fc (z, t) = F × q(z, t)
(5)
where F is the frictional coefficient. In the flow-loop experiment, for a flow below 1000 liters per minute, the model agrees well with the experiment. The comparison is presented in Fig. 2. Since the frictional coefficient is an uncertain parameter of the annulus, there is an incentive to calibrate the hydraulic model with respect to this parameter. In section
The last step is to define
Pressure drop [bara]
2
u(x, t) = u ¯(xl, t)eax v(x, t) = v¯(xl, t)e−ax
1.5
1
where
Modeled Measured
Fl a= √ 2 βρ
0.5
0 650
700
750
800
850
900
950
1000
1050
Volumetric flow rate [lpm]
Fig. 2: Frictional pressure drop in the annulus.
III, we formulate the MHE for state and parameter estimation of the hydraulic model (1)-(4).
Plugging (21)-(22) into (17)-(20), yield s 1 β F ut (x, t) = − ux (x, t) − e2ax v(x, t) l ρ 2ρ s F 1 β vx (x, t) − e−2ax u(x, t) vt (x, t) = l ρ 2ρ u(0, t)
=
v(1, t)
=
B. A Canonical Form of the Hydraulic Model The hydraulic model (1)-(4) can be written into the following canonical form of a hyperbolic-type PDE ut (x, t)
= −1 (x)ux (x, t) + c1 (x)v(x, t)
(6)
vt (x, t)
= 2 (x)vx (x, t) + c2 (x)u(x, t)
(7)
u(0, t)
= qv(0, t) + U1 (t)
(8)
v(1, t)
= U2 (t)
(9)
First, the hydrostatic pressure due to gravity from the momentum equation can be eliminated by using the following change of variable p¯(z, t)
= p(z, t) − ρgl − ρgz cos(η)
β qz (z, t) Aa F Aa qt (z, t) = − p¯z (z, t) − q(z, t) ρ ρ q(0, t) = qp (t) p¯(l, t)
= −
= pc (t) − ρgl(1 + cos(η))
(24) (25)
−v(0, t) + qp (t) (26) −a A e a (pc (t) − ρgl(1 + cos η)) u(1, t)e−2a − √ βρ (27)
The above system resemble (6)-(9). The original state can be computed using these formulas z a z a q(z, t) = u (28) , t e− l z + v ,t elz l √ l a a βρ z z p(z, t) = , t e− l z − v ,t elz u Aa l l +ρgl + ρgz cos(η) (29) The distributed model (6)-(9) can be simplified into a lumped model using an early lumping approach as we describe in the following section.
(11)
C. Early Lumping Approach
(12) (13) (14)
Next, using the Riemann coordinate transformations [15] Aa 1 (15) q(z, t) + √ p¯(z, t) u ¯(z, t) = 2 βρ 1 Aa v¯(z, t) = q(z, t) − √ p¯(z, t) (16) 2 βρ (11)-(14) are transformed into s β F u ¯t (z, t) = − u ¯z (z, t) − (¯ u(z, t) + v¯(z, t)) (17) ρ 2ρ s β F v¯t (z, t) = v¯z (z, t) − (¯ u(z, t) + v¯(z, t)) (18) ρ 2ρ = −¯ v (0, t) + qp (t) (19) Aa v¯(l, t) = u ¯(l, t) − √ (pc (t) − ρgl(1 + cos(η))) (20) βρ
u ¯(0, t)
(23)
(10)
Plugging (10) into (1)-(4), yields p¯t (z, t)
(21) (22)
The idea of early lumping approach is to replace the infinite-dimensional PDE (6)-(9) with a finite set of ODE using the method of lines. Here, the spatial coordinate x is discretized using a non-uniform grid with N -nodes where the variable ui and vi for all i ∈ {1, · · · , N } denote the evaluation of u and v at the grid nodes. Hence, for i = 1, · · · , N , equations (6)-(7) can be written as dui dt dvi dt
= =
dui + c1 (xi )vi dxi dvi 2 (xi ) + c2 (xi )ui dxi
−1 (xi )
(30) (31)
The spatial derivative is discretized using an upwind finite difference scheme as follow dui dxi dvi dxi
= =
ui − ui−1 δx vi+1 − vi δx
(32) (33)
where δx = N1+1 . Remark that the discretization schemes for u and v are slightly different due to boundary conditions.
Plugging (32)-(33) to (30)-(31), yield du1 dt du2 dt duN dt dv1 dt dv2 dt
III. MOVING HORIZON ESTIMATION DESIGN
u1 − u0 + c1 (x1 )v1 δx u2 − u1 = −1 (x2 ) + c1 (x2 )v2 δx .. . = −1 (x1 )
(34) (35) (36)
uN − uN −1 = −1 (xN ) + c1 (xN )vN δx v2 − v1 = 2 (x1 ) + c2 (x1 )u1 δx v3 − v2 = 2 (x2 ) + c2 (x2 )u2 δx .. .
(37) (38)
ˆ ∆y = y − C x
(40)
x˙ = Ax + BU
(42)
[u1 · · · uN v1 · · · vN ]T
The state matrix is given by E1 A= C2
C1 E2
(43)
C2
E1
0
···
E2
=
tk −τ
ˆ k −τ,θ
ˆ˙ = Aˆ ˆ (tk − τ ) = x ˆ tk −τ (53) x x + Bu, x
−2 (x1 ) 2 (x1 ) · · · −2 (x2 ) · · · 1 0 .. .. .. δx . . . 0 0 ···
The input matrix is given by 1 1 (x1 ) 0 B= 0 −2 (xN ) δx while the control vector is given by U1 U= U2
with the positive semi-definite matrix Q ∈ Rp×p weighting the estimation error (51). Remark that the optimization ˆ tk −τ and variables for solving (52) are the initial value x the parameter θˆ (in this case, the frictional coefficient F ). B. Gradient Computation
1 (x1 ) (x1 ) −q 1δx 2q δx1 1 = diag(c1 ) + 0 0 = diag(c2 ) −1 (x1 ) 0 ··· (x ) − (x ) · ·· 1 2 1 1 2 = .. .. .. δx . . . 0
ˆt x
(44)
where C1
(51)
during the preceding time interval t ∈ [tk − τ, tk ] with the estimation horizon t > 0. Here, y contains p measurements and C ∈ Rp×2N denotes matrix which appropriately interˆ at the measurement locations. polates the estimator state x The dynamic optimization problem for estimating the current ˆ k can be formulated as follow state x Z tk J(ˆ xtk −τ,θˆ) = min ∆y | Q∆y dt (52) subject to
where the state vector is given by =
The idea of MHE is to estimate the current states by solving a least squares dynamic optimization problem, which penalizes the deviation between the measurements and the predicted outputs, with the initial value of the estimator system and/or the parameters of the system considered as the optimization variables [16]. Mathematically, we want to minimize the estimation error
(39)
U2 (t) − vN dvN = −2 (xN ) + c2 (xN )uN (41) dt δx The left boundary condition can be computed by interpolating v0 ≈ 2v1 − v2 ; hence, u0 = q(2v1 − v2 ) + U1 . The ODEs (34)-(41) can be written in a state-space representation as
x
A. Formulation of Dynamic Optimization
0 0
(45)
0 0 .. .
(46)
(47) −1 (xN ) 0 0 (48) .. . −2 (xN )
The dynamic optimization problem (52)-(53) can be solved using a gradient-based method. It has been known, e.g. see [17], [18], that the gradients of the cost function with respect to the initial condition and the parameter θˆ are given by ˆ ∂J(ˆ xtk −τ , θ) ˆ tk −τ ∂x ˆ ∂J(ˆ xt −τ , θ) k
∂ θˆ
= λ(tk − τ ) =
ˆ ∂H(ˆ xtk −τ , λ, θ) ∂ θˆ
(55)
where the adjoint state λ(tk − τ ) follows from solving the following costate equations ˆ˙ x
(49)
(54)
= Aˆ x + Bu ˆ (tk − τ ) = x ˆ tk −τ x ˆ ˆ x, λ, θ) ∂ H(ˆ λ˙ = − ˆ ∂x λ(tk ) = 0
(56) (57) (58) (59)
with the Hamiltonian is given by (50)
To enable online calibration of uncertainties in the friction characteristic, the frictional coefficient F is considered as a tuning parameter for the model. We highlighted the MHE method for state and parameter estimation in the following section.
ˆ = ∆y | Q∆y + λ| (Aˆ H(ˆ x, λ, θ) x + Bu)
(60)
(0) ˆ tk −ˆτ and θˆ(0) , one Starting with an initial estimate of x gradient iteration j consists of ˆ (j) , • integrating (56)-(57) forward in time to obtain x (j) • integrating (58)-(59) backward in time to obtain λ ,
Parameter qp ρ l Aa β
Value 1/60 1000 700 154.8 2×10−9
Unit m3 /s kg/m3 m mm Pa
TABLE I: Experimental data
•
(j)
g1
(j)
g2 •
(a) The main pump.
(b) The pipes.
(c) Downhole sensors.
(d) Topside sensors.
computing the gradient =
λ(j) (tk − τ )
(61)
=
ˆ ∂H(ˆ xtk −τ , λ, θ) ∂ θˆ
(62)
updating the estimated state (j+1)
¯ tk −τ x θˆ(j+1)
= =
(j)
(j)
ˆ tk −τ − α1 g1 x (j) θˆ(j) − α2 g 2
(63) (64)
where α1 and α2 are the line search parameters. The gradient ˆ iterations, which leads to the algorithm is stopped after N estimate of the current state ˆ
ˆk = x ˆ N (tk ) x
Fig. 4: Experiment setup
(65)
IV. EXPERIMENT DESIGN The pipes are placed nearly horizontally, with the elevation observed in Fig. 5. Therefore, the hydrostatic pressure has a relatively small effect, and the pressure difference is primarily due to the frictional pressure. 60 50
Vertical Depth [m]
The experiment is carried out in a field scale flow-loop test in Stavanger, Norway. The MPD system is configured as a U-tube (Fig. 3), which consists of the main pump (Fig. 4a), 700 meters of pipes (each for drill string and annulus) (Fig. 4b), the bottom hole assembly (Fig. 4c), and the top side sensors (Fig. 4d). The data used in this experiment are given in Table 1.
40 30 20 10 0 0
100
200
300
400
500
600
700
Well Path [m]
Fig. 5: Loop profile in XZ-coordinate. V. SIMULATION RESULTS
Fig. 3: Schematics of the flow-loop. In this experiment, water is injected by the main pump through the drill string and up the annulus. The water rate is changed so that the dynamic of the model can be observed. The aim is to estimate the current pressure and flow throughout the annulus using previous measurement from the top of the annulus.
The spatial discretization is realized with N = 10 grid nodes. The sampling time and prediction horizon are chosen as ∆t = 1 s and τ = t = 10 s, respectively. Since the measurement taken only at the top of the well, then p = 2 sensors are considered. Hence, the measurement vector is u|x=xN given as y = . The hydraulic model is solved v|x=xN using Matlab solver ode23, while the step sizes α1 and α2 are obtained from an adaptive line search algorithm [10]. Implementing the MHE algorithm described in section III, the estimation result for the pressure and flow inside the well can be found in Fig. 6, while the estimated frictional
coefficient can be seen in Fig. 7. It takes less than 1 s to run a one step MHE, which shows the real-time feasibility of the MHE scheme. 800
Estimated Measured
Volumetric flow rate [lpm]
Volumetric flow rate [lpm]
800
750
700
650
600 0
50
100 150 Time [s]
200
250
(a) Flow rate through bit.
750
700
650
R EFERENCES 600 0
25
20
50
100 150 Time [s]
200
250
(c) Pressure at bit.
100 150 Time [s]
30
Estimated Measured
15 0
50
200
250
(b) Flow rate through wellhead.
Pressure [bara]
Pressure [bara]
30
Estimated Measured
Estimated Measured
25
20
15 0
50
100 150 Time [s]
200
250
(d) Pressure at wellhead.
Fig. 6: Pressure and flow at the bottom (bit) and topside (wellhead) of the well.
93
Coefficient of friction
92.5 92 91.5 91 90.5 90 0
50
ACKNOWLEDGMENT Financial support from Statoil ASA and the Norwegian Research Council (NFR project 210432/E30 Intelligent Drilling) is gratefully acknowledged. The experimental data are results from ongoing internal research in the department Intelligent Drilling in Statoil RDI. The authors gratefully acknowledge the work done by Alexey Pavlov, Henrik Manum and Glenn-Ole Kaasa in planning and conducting experiments.
100
150
200
250
Time [s]
Fig. 7: Coefficient of friction as tuning parameter. We can see a good performance of MHE for flow, pressure, and friction characteristic estimations. VI. CONCLUSIONS AND RECOMMENDATIONS We have implemented MHE for estimating annular pressure and flow in MPD where the coefficient of friction serves as the tuning parameter for the model. The dynamic optimization problem is solved efficiently using a gradientbased method combined with a line search method. Because only one phase (fluid) is involved, the annular model is represented by a one-phase incompressible flow model. There is an interest to include the gas phase, e.g., when drilling a gas reservoir, leading to a two-phase flow model (refers to the drift flux model). This is a more challenging problem that should be addressed in the future research.
[1] J. M. Godhavn, Control Requirements for Automatic Managed Pressure Drilling System, SPE Drilling and Completion, vol. 25, no. 3, pp. 336–345, 2010. [2] M. Davoudi, J. R. Smith, J. E. Chirinos, and B. M. Patel, Evaluation of Alternative Initial Responses to Kicks Taken During ManagedPressure Drilling , SPE Drilling and Completion, vol. 26, no. 2, pp. 169–181, 2011. [3] J. E. Olsen, E. Vollen, and T. Tonnessen, Challenges in Implementing UBO Technology, SPE/IADC Underbalanced Technology Conference and Exhibition, 11-12 October, 2004, Houston, USA. [4] G. O. Kaasa, O. Y. Stamnes, O. M. Aamo, and L. Imsland, Simplified Hydraulics Model Used for Intelligent Estimation of Downhole Pressure for a Managed-Pressure-Drilling Control System, SPE Drilling and Completion, vol. 27, no. 1, pp. 127–138, 2012. [5] K. S. Bjorkevoll, R. Rommetveit, B. Aas, H. Gjeraldstveit, and A. Merlo, Transient Gel Breaking Model for Critical Wells Applications with Field Data Verification, SPE/IADC Drilling Conference, 19-21 February, 2003, Amsterdam, The Nedherland. [6] J. Petersen, R. Rommetveit, K. S. Bjorkevoll, and J. Froyen, A General Dynamic Model for Single and Multi-phase Flow Operations during Drilling, Completion, Well Control and Intervention, IADC/SPE Asia Pacific Drilling Technology Conference and Exhibition, 25-27 November, 2008, Jakarta, Indonesia. [7] E. L. Haseltine and J. B. Rawlings, Critical Evaluation of Extended Kalman Filtering and Moving-Horizon Estimation, Industrial and Engineering Chemistry Research, vol. 44, no. 8, pp. 2451–2460, 2005. [8] P. Kuhl, M. Diehl, T. Kraus, J. P. Schloder, and H. G. Bock, A Realtime Algorithm for Moving Horizon State and Parameter Estimation, Computer and Chemical Engineering, vol. 35, no. 1, pp. 71–83, 2011. [9] O. N. Stamnes, O. M. Kaasa, and G. O. Kaasa, Redesign of Adaptive Observers for Improved Parameter Identification in Nonlinear Systems, Automatica, vol. 47, no. 2, pp. 403–410, 2011. [10] J. Nocedal and S. Wright, Numerical Optimization, Springer, 2006. [11] M. Paasche, T. A. Johansen, and L. Imsland, Regularized and Adaptive Nonlinear Moving Horizon Estimation of Bottomhole Pressure during Oil Well Drilling, IFAC World Congress, 28 August-2 September, 2011, Milan, Italy. [12] D. Sui, R. Nybo, S. Hovland, and T. A. Johansen, A Moving Horizon Observer for Estimation of Bottomhole Pressure During Drilling, IFAC Workshop on Automatic Control in Offshore Oil and Gas Production, 31 May-1 June, 2012, Trondheim, Norway. [13] F. White, Fluid Mechanics, McGraw-Hill, New York, 2007. [14] I. S. Landet, A. Pavlov, O. M. Aamo, Modeling and Control of HeaveInduced Pressure Fluctuations in Managed Pressure Drilling, IEEE Transaction on Control System Technology, vol. 21, no. 4, pp. 1340– 1351, 2013. [15] G. Bastin and J. M. Coron, Further Result on Boundary Feedback Stabilization of 2×2 Hyperbolic System over a Bounded Interval, IFAC Symposium on Nonlinear Control System, 1-3 September, 2010, Bologna, Italy. [16] H. Michalska and D. Q. Mayne, Moving Horizon Observers and Observer-based Control, IEEE Transaction on Automatic Control, vol. 40, no. 6, pp. 995–1006, 1995. [17] Y. Cao, S. Li, L. Petzold, and R. Serban, Adjoint Sensitivity Analysis for Differential-Algebraic Equations: The Adjoint DAE System and its Numerical Solution, SIAM Journal on Scientific Computing, vol. 24, no. 3, pp. 1076–1089, 2003. [18] S. Rhein, T. Utz, and K. Graichen, Model Predictive Control and Moving Horizon Estimation of a Large-Scale Chemical Reactor Model, IFAC Workshop on Control of Systems Modeled by Partial Differential Equations, 25-27 September, 2013, Paris, France.
