Optimisation of a Gas-Lifted System with Nonlinear Model Predictive Control

Abstract

A gas-lifted system in a mature oil well can experience casing-heading instability, which reduces the mean oil production and it is not healthy for the downstream equipment. This instability was removed by implementing a terminal equality-constrained nonlinear model predictive control (NMPC) having input targets with control zones in the system. The input-dependent stability behaviour of the gas-lifted system was visualised through the bifurcation diagram and the step responses of the linearised model at various operating points. The controller was then presented. Then, the close-loop feasibility, as well as the convergence, were discussed. The controller stabilised the undisturbed gas-lifted system, improving production by 5.63% compared to the open-loop operation when the system was in casing-heading instability. For the two input case, the steady-state production, aided by the high-input target, reached 12.25 kg/s, which was far more than 9.57 kg/s for the one input case. This controller showed a 3.76% improvement over the PI controller for the same purpose.

1. Introduction

As the production of hydrocarbon from beneath the earth’s surface increases, the natural energy for lifting it to the downstream facilities at the desired rate declines. At a certain low production rate, it becomes necessary for an artificial approach to be used to lift the hydrocarbons [1]. Various artificial lift methods have been developed to maximise oil recovery while minimising the total cost of delivering crude oil to storage facilities. These artificial lift methods include sucker-rod pumps, electrical submersible pumps, gas-lifted systems among others [2]. In addition to these are water flooding and $C{O}_2$ flooding [3,4,5,6]. The most commonly used lift method in the field with a high gas-oil ratio, with densely populated wells, is the gas-lifted method [7]. This gas-lifted method accounts for over 70% of the artificial lift methods in Brazil [8]. Gas lift reports based on type (continuous and intermittent), application (offshore and onshore), companies (General Electric, Schlumbergger etc.) and regions (North America, Europe etc.) are detailed in [9].

The natural source for the lift energy for the gas-lifted system is provided by the reservoir pressure ($P_r$). When this energy is insufficient, oil production falls below the desired level. Lift gas is then delivered to the system through the lift gas valve and injected into the tubing through the injection valve. This mixes with the liquid in the tubing and reduces the density of the oil and the bottom hole pressure, hence the available energy can then lift the crude as desired. The lift gas is obtained from crude oil itself and compressed before delivering to the gas-lifted system through the gas-lift valve. Methane is usually preferred as it is the lightest of the crude oil constituents. Details of crude oil constituents are found in [10].

The optimal operation of a gas-lifted system is influenced by the system flow assurance. Flow assurance implies that oil from the reservoir flows reliably and consistently till it reaches the storage locations [11]. Flow assurance is influenced by many factors, such as hydrates, wax, asphatenes, scales, slugging flow, emulsion and corrosion [12]. In addition to these are casing-heading instability (a form of slugging flow) and faults in the system since [13] showed the relationship between system reliability and flow assurance. Casing-heading instability is the oscillatory behaviour of the gas-lifted system based on the input or inputs. The phenomenon of casing-heading instability is described briefly:

As gas is delivered to the annulus, the annulus pressure increases, leading to an increased pressure difference with the tubing. This makes the lift gas flow into the tubing to reduce the density of the tubing oil. Oil production increases, resulting in a further decrease in bottom-hole pressure. More gas then flows into the tubing. The increased flow reduces the annulus pressure till the pressure difference is insufficient to permit gas flow into the tubing. The oil then builds up in the tubing while gas builds up in the annulus until the annulus pressure becomes large enough to allow the flow of gas into the tubing. A new cycle then starts. This leads to the oscillatory behaviour of the gas-lifted states and variables. Further theoretical descriptions of casing-heading instability can be found in [14].

This oscillatory behaviour of the gas-lifted system reduces the mean oil production, and it is also not healthy for the downstream equipment. The behaviour could lead to a violation of the operational constraints of the gas-lifted system or environmental regulations. Hence it is not desired and must be removed. Am approach to remove this casing-heading instability involves: (a) increasing the pressure drop caused by friction and (b) control applications. Pressure drop due to friction can be increased by reducing the flowrate through production choke, increasing flow into the annulus and operating the injection valve at critical points [15]. Increasing the pressure drop due to friction can be limited by gas availability and can also be very conservative. The control approach is, therefore, the optimal approach to remove casing-heading instability and optimise the gas-lifted system.

The control approach to casing-heading instability removal has traditionally been implemented using a simple proportional-integral (PI) controller [15,16,17,18] or proportional-integral-derivative (PID) controller [19,20]. The controller uses the production choke valve as the manipulated variable in most cases. The controlled variables, however, can be annulus pressure, bottom-hole pressure, flowrate through the production choke, among others. This could result in the case where the output of the controller may be different from the plant input due to constraints in the controlled variable or manipulated variable. In cases where the measurements of the states or controlled variables are not available or reliable, they are estimated using observers and Kalman filters before control application [15,21,22].

The model-based approach, where an explicit model is used to predict the gas-lifted system states and variables for use in the control solution, has gained traction lately [12,23,24,25,26]. Model predictive control (MPC) uses the model of the system to predict the future state trajectory at a sample time over a given horizon. The predicted states are used to solve an online optimal control problem (OCP) to obtain the optimal input. The first element of the optimal input sequence is then applied to the system at the given time. At the next sample time, the horizon recedes, and the above steps are repeated [27]. This is performed to accommodate disturbances and model inaccuracies that make predicted states of being different from the actual states. Zone control and input target is implemented when there are more controlled variables than the manipulated variable, such as a gas-lifted system with one or two inputs and three states. The zones are selected based on process needs, and the control objective is the first to ensure the states are in the zones before enforcing the input targets [28].

Most of these model-based approaches are implemented with linear models, as this is usually sufficient enough for most control cases [18]. However, the investment in the gas-lifted system is prohibitively high, and regulations on environmental impact are tight. Hence there is the need to operate it as close to the optimal level as possible. With the advancement in nonlinear control and this tighter demand on the operation of the gas-lifted system, the nonlinear operation becomes necessary to meet these demands [29].

Nonlinear MPC (NMPC) was used in improving a gas-lifted system in [26], where a nonlinear predictive generalised minimum variance control was used in stabilising the casing heading. The control law was obtained from the minimisation of the cost function, which involves only the deviation from the outputs, and did not consider the input. In [30], an NMPC that considers input deviation and input targets as part of the cost function in addition to output deviation was used to improve oil production by considering the choke valve and gas lift valve as inputs. The approach, however, involves dynamic linearization, which makes the system more like a linear time-varying system whose stability is difficult to ascertain [30]. With dynamic linearisation, Ref. [8] used NMPC for gas lift optimsation using only injection gas flow as the input for the multi-well situation.

In this article, NMPCm having input target and control zones is used to stabilise the casing heading in a gas-lifted system. The NMPC that incorporates the system constraints in its formulation is compared with the PI controller that is commonly used in the literature. Unlike in [30], the states are considered in the cost function since we rely on its estimation with Extended Kalman filter (EKF) using annulus pressure ($P_a$), pressure of tubing top ($P_{wh}$) and mixture density (${\rho }_t$) [22].

The key contributions are the presentation of the end-constraint NMPC with simplified convergence and recursive feasibility details. Further, the use of an artificial reference for states that are computed online and not decided offline takes into account the optimal input at the given sample time in the NMPC. The use of zone control ensures that control effort focuses on the input target only when the states in the zones are being considered for the first time in a gas-lifted system here. One input case (production choke opening) and two input cases (production choke opening and gas lift valve flow) are used in the optimisation and disturbance attenuation by the controller. Finally, the case where the input target is unreachable, ensuring that oil production is increased as the nonlinear MPC drives the input towards the upper control boundary, is presented. A comparison is made with the open-loop simulation and the use of a PI controller.

This paper is organised as follows: Section 2 presents the analysis of the gas-lifted system visualising the phase portrait, bifurcation and limit cycle. Section 3 presents the terminal equality-constrained nonlinear model predictive control (NMPC) and discusses the feasibility and convergence of the controller. Section 4 presents the stabilisation of the casing-heading instability using the NMPC. Section 5 concludes the paper.

2. Gas-Lifted System: Models and Analysis

We present the gas-lifted system models and examine the input-dependent stability behaviour of the system using a bifurcation diagram and linearised models. The analysis is limited to the simulation of these behaviours. The model of the gas-lifted system used here is first presented.

2.1. Gas-Lifted System Models

Figure 1 shows a gas-lifted system with the states and variables. The states are the mass of gas in annulus ($x_1$), the mass of gas in tubing ($x_2$) and the mass of oil in tubing ($x_3$). The input is the percentage opening of the production choke. The range is from 0 to 1, which represents 0% and $100\%$, respectively. The gas-lifted system contains two fluid volumes: outer annulus, which holds $x_1$, and inner tubing, which holds $x_2$ and $x_3$.

The model of the gas-lifted system is a slight modification of the model presented in [22]. The model is a differential algebraic equation (DAE) of the form:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \dot{x}=f(x,z,u)\hfill \end{array}\end{array}$$

(1a)

$$\begin{array}{cc}\hfill & 0=g(x,z,u)\hfill \end{array}$$

(1b)

The differential equations are due to mass balances, while the algebraic equations are due to the mass flowrate, density and pressure equations.

The mass (differential) equations:

$$\begin{array}{cc}\hfill & \frac{d{x}_1}{dt}=w_{gl}-w_{iv}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & \frac{d{x}_2}{dt}=w_{iv}+w_{rg}-w_{pg}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & \frac{d{x}_3}{dt}=w_{ro}-w_{po}\hfill \end{array}$$

(4)

where $w_{gl}$ is the gas flowrate into the annulus, $w_{iv}$ is the gas flowrate from the annulus into the tubing, $w_{rg}$ is the gas flowrate from the reservoir into the tubing, $w_{ro}$ is the oil flowrate from the reservoir into the tubing, $w_{pg}$ is the gas flowrate through the choke and $w_{po}$ is the oil flowrate through the choke. All flowrate units are in kg/s, including the flowrate of the mixture through the choke ($w_{pc}$), while the states are in kg.

Flowrate algebraic equations:

$$\begin{array}{cc}\hfill & w_{iv}=C_{iv}\sqrt{max(0,{\rho }_a(P_a-P_w))}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & w_{pc}=C_{pc}\sqrt{max(0,{\rho }_t(P_{wh}-P_s))}f\left(u\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & w_{ro}=C_{iv}\sqrt{max(0,{\rho }_0(P_r-P_{bh}))}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & w_{pg}=\left(\frac{x_2}{x_2+x_3}\right)w_{pc}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & w_{po}=\left(\frac{x_3}{x_2+x_3}\right)w_{pc}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & w_{rg}=GOR{w}_{ro}\hfill \end{array}$$

(10)

where $P_a$ is the annulus pressure, $P_{wh}$ is the wellhead pressure, $P_w$ is the tubing pressure, $P_{bh}$ is the bottom-hole pressure, $P_r$ is the reservior pressure and $P_s$ is the separator pressure. All units are N/m. Similarly, ${\rho }_a$ is the annulus gas density, and ${\rho }_t$ is the tubing mixture density. All units are kg/m${}^3$. GOR is the ratio of gas to oil in the reservoir. It has no unit.

The pressure algebraic equation:

$$\begin{array}{cc}\hfill & P_a=\left(\frac{T_{a}R}{V_a{M}_w}+\frac{g{L}_a}{V_a}\right)x_1\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & P_{wh}=\frac{T_{t}R}{M_w}\left(\frac{x_2}{V_t-\frac{x_3}{{\rho }_0}}\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & P_w=P_{wh}+(x_2+x_3)\frac{g}{A_t}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & P_{bh}=P_w+{\rho }_{0}g{H}_{bh}\hfill \end{array}$$

(14)

Density algebraic equations:

$${\rho }_a=\frac{M_w}{T_{a}R}{P}_a$$

(15)

$${\rho }_t=\frac{x_2+x_3}{V_t}$$

(16)

$$f\left(u\right)={50}^{u-1}$$

(17)

The friction factor discussed in [31] is ignored here as we use MPC for our stabilisation. This is because the MPC performance is still acceptable, even with up to 20% model inaccuracy [32]. The parameters’ values, definitions and units are given in [22]. Equation (17) applies the input to the mixture flowrate to prevent the input from being zero. The max operator prevents the flowrate from becoming negative when the upstream pressure is less than the downstream pressure by returning a zero.

2.2. The Stable and Unstable Gas-Lifted System

The DAE system is solved using ODE15s in MATLAB using the data provided in the appendix of [33]. To simulate the existence of casing-heading instability and to ensure a consistent solution to the DAE for the selected initial conditions, the flowrate into the annulus and GOR were altered. The input is u = 0.65 (65% valve opening), while $w_{gl}$ (flowrate into the annulus) is fixed at 0.4 kg/s. Starting from the initial condition of $x_0$ = ${\left[\begin{array}{ccc}2100& 400& 5700\end{array}\right]}^T$ kg, and with GOR fixed at 0.01, the system is simulated over 5 h. For a GOR of 0.01 (shown later), the system remains in a stable mode till the inputs become u = 0.75.

Figure 2 shows the states of the gas-lifted system. The gas-lifted system is seen to be a positive system as all the states have positive values. All the algebraic variables also have positive values. In Figure 2, the system shows no oscillatory behaviour, as the states converged to a steady state value of $x_s$ = ${\left[\begin{array}{ccc}2131& 347& 6429\end{array}\right]}^T$ kg, showing the system is stable. Operating the gas-lifted system at $u=0.65$, therefore, prevents the casing-heading instability, but the system operates at a low production rate as the flowrate depends on the input.

For the unstable regime, we use a choke opening of 95% (0.95), which falls into the unstable region. Figure 3 shows the oscillating states of the system. The system oscillates with a period of about 53 min (0.88 h). All algebraic variables of the gas-lifted system also show oscillatory behaviour. This period can be as low as 100 s (1.7 min) on a laboratory scale to a few hours in some fields [22]. Further simulations show that any input such that $0.0\le u\le 0.75$ will produce a gas-lifted system that is completely stable. For $0.75\le u\le 0.80$, there is slight oscillation. However, the bifurcation point is at u = 0.80 after which the system goes into oscillation.

2.3. Bifurcation

The dynamics of the gas-lifted system are highly dependent on the parameters. One of these parameters is the GOR, and the parameter within our control is the input, which is the percentage valve opening when we consider one input or flowrate into the annulus as a second input. By parameterising the input, we can observe the behaviour of the gas-lifted system output as the input changes from u = 0% (0.00) to u = 100% (1.00). This reveals the critical point at which the system changes from a stable to an unstable oscillatory regime. This critical point is GOR-dependent, as shown in Figure 4 for GORs of 0.0001, 0.001, 0.01 and 0.1, respectively. The black dashed lines are the minimum and maximum flowrates of produced oil, while the red solid is the average flowrate.

In all four cases in Figure 4, during the stable region, the minimum, mean and maximum oil production are the same. This stable production starts from u = 0 (0%) to about u = 0.67 (67%) for both GOR of 0.0001 and 0.001, while it is u = 0.75 (75%) for a GOR of 0.01. Slight oscillation starts at u = 0.67 (67%), and the bifurcation point is at u = 0.73 (73%) for GORs of 0.0001 and 0.001. Slight oscillation starts at u = 0.75 (75%) for a GOR of 0.01, while its bifurcation point is at u = 0.8 (80%). Thereafter, full oscillation sets in. In the case of a GOR of 0.1, there is no oscillation for any input value. This is because high GOR implies more gas in the tubing and a low chance of casing-heading instability. Further, note that there exists a slight difference in the critical points for the GORs of 0.0001 and 0.001, implying that a GOR of 0.001 is almost the same as no gas in the reservoir. The critical point is the same when any of the variables of the system is selected in place of the oil production rate ($w_{po}$).

The behaviour of the equilibrium points around the bifurcation point changes from stable to unstable regime. The eigenvalues for the equilibrium points corresponding to u = 0.7, 0.8 and 0.9 for the system with GOR = 0.01 are [−0.0171, −0.0003 ± 0.00014], [−0.0201, −0.0001 ± 0.00017] and [−0.0261, 0.0001 ± 0.0020], respectively. These points correspond to points before bifurcation, at the bifurcation point and when the system is fully in oscillation, respectively. In the equilibrium points corresponding $u=0.7$ and $u=0.80$, the real parts of the eigenvalues do not cross the imaginary line. Beyond the bifurcation point, two conjugate eigenvalues crossed into the right-hand plane; hence the system becomes unstable. The eigenvalues and the result in Figure 4, which show a change in the behaviour of oil the production rate as the bifurcation point is crossed, indicate that the gas-lifted system with the given parameter values exhibit local bifurcation. The fact that there is a change in the stability of the system based on eigenvalues crossing the imaginary line during bifurcation favours nonlinear control application.

2.4. Linearised Gas-Lifted System

We examine the behaviour of the gas-lifted system around the equilibrium point to guide the comparison of the controller presented in this article with the PI controller. The equilibrium point of a nonlinear system is the point at which the system dynamics die off. Consequently, the behaviour of the system appears to be better studied around this point than any other point since it is more complicated around this region.

The gas-lifted system is linearised around the operating points corresponding to the input values from $u=0.65$ to $u=1.00$ in 0.05 increments. The result from calculating the linear matrices is compared to that obtained from the Simulink model using the Linear Analysis tool in MATLAB and is the same or very close. The resulting state space model for steady state corresponding to $u=0.70$ is given in Equation (18). These models are used to obtain step responses around the bifurcation point. Figure 5 shows these step responses around the bifurcation point.

$$\begin{array}{c}\hfill A=\left[\begin{array}{ccc}-0.003545& 0.007077& 0.0009127\\ 0.003545& -0.01176& -0.0009194\\ 0& -0.04282& -0.00243\end{array}\right],B=\left[\begin{array}{c}0\\ -2.77\\ -42.24\end{array}\right]\\ \\ \hfill C=\left[\begin{array}{ccc}0& 8748& 1128\\ \end{array}\right],D=\left[\begin{array}{c}0\end{array}\right]\end{array}$$

(18)

The step response for the linearised system in Figure 5 shows a gradual change in the oscillatory behaviour of the system. At $u=0.65$, the step response converged to a fixed value, while at $u=0.80$, where bifurcation occurs, the system experiences damped oscillation. Beyond this point, however, the system goes into oscillatory behaviour, as seen in Figure 5e, and becomes unstable at $u=0.90$ and beyond. This brings the issue of model uncertainty as input (which determines operating points) changes, making the linear control application insufficient, hence necessitating the use of nonlinear control.

The situation is made more obvious in Figure 6, which compares the nonlinear response with two linear responses. The initial conditions for the linear states correspond to $u=0.65$ and $u=0.75$, respectively. A step input of $u=0.05$ and $u=-0.05$ are applied to the linear models, respectively, while an input of $u=0.70$ is applied to the nonlinear model. It is expected that the three states converge at steady state to the same value since all three states operate at $u=0.70$. Only Figure 6b, however, shows convergence in the states from different initial conditions. The state that starts from a steady state corresponding to $u=0.65$ converges to a value below the nonlinear state, while the state that starts from $u=0.75$ converges above the nonlinear states in Figure 6a,c. This only creates further uncertainty as all three states are supposed to converge as they are all operated at $u=0.70$. This implies that operating the system using linear models will produce a different trajectory at a neighbourhood of the equilibrium point that is a bit far. To receive accurate state prediction and hence an optimised result, a nonlinear approach is to be used.

The gas-lifted system presents conflicting behaviour, as observed here. The bifurcation diagram shows a change in the stability of the system as input changes; the eigenvalues confirm this too. However, the limit cycle in the unstable region indicated that it is a stable limit cycle. Since the oscillatory behaviour we intend to remove (in the next section) is called casing-heading instability, this instability may be a Lyapunov instability. This is due to the fact that Lyapunov stability is assured for a system in which the energy, output or trajectory decays to zero when little or no perturbation is given to it. However, in the gas-lifted system, the system trajectory grows until it enters the limit cycle. Hence it is Lyapunov unstable. The instability may also be due to the effect it has on downstream equipment. The gas-lifted system is a DAE system, and one property of a DAE system is that the stability cannot easily be detected from a local stability analysis of the linearized system [34]. This could explain why the system exhibits casing-heading instability while the limit cycle is stable.

3. Terminal Equality-Constrained Nonlinear Model Predictive Control (NMPC) with Input Target and Control Zones

We present the controller used in this work, which is a modification of the finite horizon NMPC in [35] with end constraints.

Consider a nonlinear model of a system:

$$\begin{array}{cc}\hfill & x_{k+1}=f(x_k,u_k)\hfill \end{array}$$

(19a)

$$\begin{array}{cc}\hfill & y_k=h\left(x_k\right)\hfill \end{array}$$

(19b)

where x and y are states and outputs, respectively, while f and h are nonlinear functions of states and outputs, respectively. The infinite horizon NMPC (IHNMPC) problem is given in (20) and (21a–f):

$$\begin{array}{c}\hfill \underset{U_k,x_{s,k}}{\mathrm{minimize}}\\ \hfill & V_k=\sum _{j=0}^{\infty }\left|\right|{(x(k+j|k)-x_{s,k})}_{Q_x}^2\hfill \\ \hfill & +\left|\right|(u(k+m-1|k)-u_{des}){\left|\right|}_{Q_u}^2\hfill \\ \hfill & +\sum _{j=0}^{m-1}{\left|\right|\Delta u(k+j|k)\left|\right|}_R^2\hfill \end{array}$$

(20)

Subject to

$$\begin{array}{cc}\hfill & x_{min}\le x_{s,k}\le x_{max}\hfill \end{array}$$

(21a)

$$\begin{array}{cc}\hfill & u_{min}\le u(k+j|k)\le u_{max}\hfill \end{array}$$

(21b)

$$\begin{array}{cc}\hfill & \Delta u_{min}\le \Delta u(k+j|k)\le \Delta u_{max}\hfill \end{array}$$

(21c)

$$\begin{array}{cc}\hfill & x_{s,k}=f(x_{s,k},u(k+m-1|k))\hfill \end{array}$$

(21d)

$$\begin{array}{cc}\hfill & x_{k+1}=f(x_k,u_k)\hfill \end{array}$$

(21e)

$$\begin{array}{cc}\hfill & x\left(k\right)=x_k\hfill \end{array}$$

(21f)

where the first, second and third terms in (20) are cost functions relating to state deviation from artificial reference $\left(x_{s,k}\right)$, deviation of the last input in the control sequence from the desired input $\left(u_{des}\right)$ and input moves $(\Delta u)$, respectively. $Q_x,Q_u$ and R are penalties on states, input and input moves, respectively. m is the control horizon. Equation (21a) provides the boundary condition for the artificial reference, (21b) is the boundary condition for the absolute input and (21c) is the constraint on input moves. At every sampling time, (21d) computes the steady state using the last element in the control sequence. Constraint (21e) is the equality constraint that describes the state evolution of the nonlinear system, while (21f) shows that the predicted states start from the current measurement.

If the predicted state is forced to equal the reference state ($x_{s,k}$) from the end of prediction horizon P to infinity, implying a terminal constraint, then the control problems (20) and (21) becomes (22) and (23).

$$\begin{array}{c}\hfill \underset{U_k,x_{s,k}}{\mathrm{minimize}}\\ \hfill & V_k=\sum _{j=0}^{P-1}\left|\right|{(x(k+j|k)-x_{s,k})}_{Q_x}^2\hfill \\ \hfill & +\left|\right|(u(k+m-1|k)-u_{des}){\left|\right|}_{Q_u}^2\hfill \\ \hfill & +\sum _{j=0}^{m-1}{\left|\right|\Delta u(k+j|k)\left|\right|}_R^2\hfill \end{array}$$

(22)

Subject to

$$\begin{array}{cc}\hfill & x_{min}\le x_{s,k}\le x_{max}\hfill \\ \hfill & \end{array}$$

(23a)

$$\begin{array}{cc}\hfill & u_{min}\le u(k+j|k)\le u_{max}\hfill \end{array}$$

(23b)

$$\begin{array}{cc}\hfill & \Delta u_{min}\le \Delta u(k+j|k)\le \Delta u_{max}\hfill \end{array}$$

(23c)

$$\begin{array}{cc}\hfill & x_{s,k}=f(x_{s,k},u(k+m-1|k))\hfill \end{array}$$

(23d)

$$\begin{array}{cc}\hfill & x_{k+1|k}=f(x_k,u_k)\hfill \end{array}$$

(23e)

$$\begin{array}{cc}\hfill & x\left(k\right)=x_k\hfill \end{array}$$

(23f)

$$\begin{array}{cc}\hfill & x(k+p|k)-x_{s,k}=0\hfill \end{array}$$

(23g)

Equation (23g) is the terminal constraint that forces the state deviation from a desired steady state to be zero from P till infinity. The implication is that the first term in (22) is bounded, and the controller enjoys a stability advantage due to the infinite horizon cost function. The implementation is also realisable because the first term in (22) is finite due to the cost function becoming zero from P to infinity. The downside of this is the fear of instability resulting from the controller not being able to enforce constraint (23g). This can be minimised by making P very large.

Input $u_{des}$ is the desired or target input from the real-time optimisation (RTO) layer. If the corresponding steady state $x_{s,des}=f(x_{s,k},u_{des})$ is not within the bound $(x_{min},x_{max})$, it is expected that from (21a), the optimisation will be infeasible. Constraint (23g) will then be impossible to meet, and the first term in (20) will sum to infinity as time goes to infinity. The controller, therefore, should not be defined due to the infinite horizon cost function $V_k$ becoming unbounded. However, this is avoided as this controller finds a reference state inside and close to the state boundary to make sure constraint (21a) is met and the controller is defined.

The optimal control sequence for an m length of control horizon is given by (24), while the control sequence for state prediction is given in (25).

$$\begin{array}{cc}\hfill & U_k=[u{\left(k\right|k)}^T,u{(k+1|k)}^T,\cdots ,u{(k+m-2|k)}^T,\hfill \\ \hfill & u{(k+m-1|k)}^T]\hfill \\ \hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & U_k=[u{\left(k\right|k)}^T,u{(k+1|k)}^T,\cdots ,u{(k+m-2|k)}^T,\hfill \\ \hfill & u{(k+m-1|k)}^T,u{(k+m-1|k)}^T,\cdots ,u{(k+m-1|k)}^T,\cdots ]\hfill \end{array}$$

(25)

Equation (25) implies that the input remains constant from the end of the control horizon m till the prediction horizon P and the infinite time. The cost function for this controller contains a term that penalises deviation of the input at the end of the control horizon from the target input computed by the RTO layer. Hence after many sample times, when the desired input equals the optimal input from the RTO, the second term vanishes for an undisturbed, nominal system if the desired input is reachable.

The zones are usually selected based on the system’s operating needs that put into consideration the constraints in the states. If the weights are properly selected to give priority to the first term in (20), the deviation from the steady state values dominates the cost function. This ensures that when the system states are in their zones, the control objective starts forcing the input to track the input targets or get as close as possible.

3.1. Recursive Feasibility

Suppose there exists a solution to problems (20) and (21) at k, the optimal input sequence obtained is as given in (24). At $k+1$, a feasible solution to the problem is given in (26), which is the previous solution, (24), shifted forward for one sampling time.

$$\begin{array}{cc}\hfill & U_k=[u{(k+1|k)}^T,u{(k+2|k)}^T,\cdots ,u{(k+m-2|k)}^T,\hfill \\ \hfill & u{(k+m-1|k)}^T,u{(k+m-1|k)}^T]\hfill \end{array}$$

(26)

Consequently, (21b) and (21c) are met at $k+1$ since they are unchanged in (26). Then, constraint (21d) is met since $u(k+m-1|k)$ has not changed, implying that (21a) is also met. Constraints (21e) and (21f) are met, being the state evaluation function and constraints that states must start from the current measurement, respectively. If constraint (23g) is met at a previous time step, it will be met at the next time step. Hence recursive feasibility is assured.

3.2. Convergence

Equation (27) is the cost function corresponding to input sequence (26) referring to shifted optimal input at $k+1$.

$$\begin{array}{cc}\hfill & V_k=\sum _{j=1}^{\infty }\left|\right|{(x(k+j|k)-x_{s,k})}_{Q_x}^2\hfill \\ \hfill & +\left|\right|(u(k+m-1|k)-u_{des}){\left|\right|}_{Q_u}^2\hfill \\ \hfill & +\sum _{j=1}^{m-1}{\left|\right|\Delta u(k+j|k)\left|\right|}_R^2\hfill \end{array}$$

(27)

If $V_k^{opt}$ is the optimal cost function corresponding to input sequence (24), then (27) is equivalently:

$$\begin{array}{cc}\hfill & V_{k+1}=V_k^{opt}-{(x\left(k\right|k)-x_{s,k})}^T{Q}_x(x\left(k\right|k)-x_{s,k})\hfill \\ \hfill & -\Delta u{\left(k\right|k)}^{T}R\Delta u\left(k\right|k)\hfill \end{array}$$

(28)

Then:

$$\begin{array}{cc}\hfill & V_k^{opt}-V_{k+1}={(x\left(k\right|k)-x_{s,k})}^T{Q}_x(x\left(k\right|k)-x_{s,k})\hfill \\ \hfill & +\Delta u{\left(k\right|k)}^{T}R\Delta u\left(k\right|k)\hfill \end{array}$$

(29)

The right-hand side of (29) is positive or at least zero, being the summation of quadratic terms weighted by positive-definite matrices Q and R. This is true only when the second term on the left-hand side is smaller than the first term. The cost function is, therefore, decreasing or, at worst, non-increasing. It will decrease to zero if the input target is reachable. If the input target is unreachable, it decreases to a value that is not zero but minimum, showing the convergence of the controller.

If we consider the controller with terminal equality constraint, then:

$$\begin{array}{cc}\hfill & V_k^{opt}-V_{k+1}={(x\left(k\right|k)-x_{s,k})}^T{Q}_x(x\left(k\right|k)-x_{s,k})+\hfill \\ \hfill & \Delta u{\left(k\right|k)}^{T}R\Delta u\left(k\right|k)-\hfill \\ \hfill & {(x(k+p+1|k)-x_{s,k})}^T{Q}_x(x(k+p+1|k)-x_{s,k})\hfill \end{array}$$

(30)

Equation (30), however, is a sum and difference of quadratic terms, hence for the controller to converge, the right-hand term must be non-negative, which implies:

$$\begin{array}{cc}\hfill & {(x\left(k\right|k)-x_{s,k})}^T{Q}_x(x\left(k\right|k)-x_{s,k})+\hfill \\ \hfill & \Delta u{\left(k\right|k)}^{T}R\Delta u\left(k\right|k)\ge \hfill \\ \hfill & {(x(k+p+1|k)-x_{s,k})}^T{Q}_x(x(k+p+1|k)-x_{s,k})\hfill \end{array}$$

(31)

Equation (31) gives a condition for the cost function to be non-increasing for the terminally constrained control problem. Since the terminal state constraint forces $x(k+p+1)|k$ to equal $x_{s,k}$, the right-hand side of (31) is zero. Hence the cost function is non-increasing. The resulting optimal control problem (OCP) can now be converted to a nonlinear programming problem (NLP), and the solution for the OCP is found.

4. Stabilisation of Gas-Lifted System Using Terminal Inequality Constrained NMPC with Input Target and Control Zones

The pressure of gas in the annulus depends on the mass of the annulus gas. If the mass of gas in the annulus can be maintained above a certain value, pressure in the annulus can be kept above a certain value that favours the unidirectional and continuous flow of gas from the annulus to the tubing. The third-order model of the gas-lifted system has one of the states of the system being the mass of the annulus gas ($x_1$). Oscillatory behaviour can be minimised if the mass of the annulus gas can be kept at a permissible minimum value. An NMPC with terminal equality constraint whose first objective is to ensure the states are in their zones before enforcing the input target can set a zone for the annulus gas. This gives the controller the flexibility to maintain the mass of annulus gas $\left(x_1\right)$ at a value that favours flow from the annulus to the tubing while the controller solves the input that approaches the optimal input from the RTO as much as possible.

4.1. Undisturbed Gas-Lifted Well Stabilization Using NMPC with Input within Bound: One Input Case

We consider a case where the system is undisturbed, and the input is within the input bound. There is only one input, which is the percentage valve opening. The controller presented in (22) and (23) is used. The controller parameters are: m = 4, P = 100, $Q_x=diag(\left[111\right]$×${10}^{-1})$, $Q_u$ = $diag\left(\right[1]$×${10}^6)$ and $R_u=diag(\left[1\right]$×${10}^5)$. We show the close loop behaviour for the system with desired input $u_{des}=0.95$, which is in the unstable region. The estimated states from the EKF and the optimal states are compared. The parameters of the EKF are: the state covariance, $P_0=diag$([100 10 1000]), the state noise, $Q=diag$([1 0.1 10]) and the measurement noise, $R=diag$([0.7 0.7 0.7]). The measurements used for the EKF are the mixture density, annulus pressure and topside pressure [22], while the detailed EKF procedure is described in [36]. We assume a gaussian distribution for the noise. Further, the system is assumed to be disturbance-free for most simulation cases except when stated.

Zone control and input target are implemented for the close-loop system. The zones are selected based on the limitation imposed by the gas-lifted system volumes and the desire to increase the annulus pressure while reducing the well pressure. This implies a large lower bound of the mass of gas in annulus $x_1$ and a small upper bound of the mass of gas in tubing $x_2$ and mass of oil in tubing $x_3$. From several simulations around the steady state, the control zones are selected as $x_{min}$ = ${\left[\begin{array}{ccc}2000& 100& 1000\end{array}\right]}^T$ kg, $x_{max}$ = ${\left[\begin{array}{ccc}2600& 800& 7000\end{array}\right]}^T$ kg.

Figure 7 shows the estimated and optimal states. The estimated states are used in the controller for prediction and are obtained from measurements, as detailed in [36]. The optimal states are obtained by using the first element of the optimal input sequence on the gas-lifted system model. The blue dashed lines are the upper and lower limits for the states (control zone), while the black solid line is the system optimal states, and the red dash-dotted line is the estimated states.

Figure 7 shows that the control zones are respected. The controller can enforce zone control for narrower control zones but at the expense of a low-optimal input sequence, which results in a low oil production rate. Hence, for an adverse case where the controller is unable to stabilise the oscillatory behaviour, the lower bound of $x_1$ can be increased to a limit dictated by the annulus volume, while the upper bounds of $x_2$ and $x_3$ can be decreased. This will jeopardise feasibility as the control zones will be narrower. The estimated states (red long-dotted line) are seen to converge to the optimal states (dark continuous lines) obtained from the actual measurement. This shows the performance of the EKF.

Figure 8 is the oil production rate at 95% valve opening (u = 0.95). The optimal production from the controller (black solid line), unstable production (red dashed line) and mean production (blue dash-dotted) were obtained at GOR = 0.01, while the stable production (green dotted line) was obtained at GOR = 0.4. The mean production is the average production when the gas-lifted system is in the unstable region at u = 0.95. The steady-state optimal production rate from the controller is 9.38 kg/s, while the mean production is 8.88 kg/s, and the stable oil production is 14.95 kg/s. If the reservoir parameter’s GOR is high enough to keep the gas-lifted system in stable mode, oil production will be at its highest rate for an input u = 0.95. However, if the GOR is not sufficient enough to naturally remove the instability at the given input, NMPC is better. The NMPC increased steady-state oil production by 5.63% over mean oil production when the gas-lifted system is operated at u = 0.95.

4.2. Undisturbed Gas-Lifted Well Stabilisation Using Terminal Equality Constrained NMPC with Input out of Bound: One Input Case

When it is desired to approach the input limit as close as possible, the input from RTO is no longer selected as the desired input. Since this controller respects the control limits, a desired input outside this limit was selected to force the optimal input close to the control limit without violating the control limit. We start with $u_{des}$ = 0.65 and switch to $u_{des}$ = 1.25 (which is out of the control limit) after 90 min (1.5 h). Figure 9 shows the states of the system, while Figure 10 and Figure 11 show the input target and cost function, respectively.

The second part of Figure 9, where $u_{des}$ = 1.25, is out of the input bound. The control zones are respected except for the mass of oil in tubing that temporarily violated the zones after the switch from the stable input within bound to input out of bound. Figure 10 shows that the first part where u = 0.65 is both within bound and reachable. The optimal control input converged to the desired input, and the cost function corresponding to the portion in Figure 11 decays to zero monotonically. This is not evident in the diagram, but on zooming, it shows how the cost decreases to zero. The second part of Figure 10, however, is out of input bound and unreachable; the optimal input is unable to reach desired input as well as the input bound. The corresponding cost function in Figure 11 decreases to a point where it leaves an offset due to the desired input being unreachable. The steady-state oil production, however, increased to 9.57 kg/s compared to 9.38 kg/s when the desired input is kept at u = 0.95.

4.3. Disturbed Gas-Lifted Well Stabilisation Using NMPC with Input within Bound: One Input Case

A system where disturbance causes the state to go out of the zone is considered here. The new state bounds are $x_{min}$ =${\left[\begin{array}{ccc}2000& 100& 1000\end{array}\right]}^T$ kg, $x_{max}$ =${\left[\begin{array}{ccc}2600& 500& 7000\end{array}\right]}^T$ kg, which is a reduction of the upper limit of $x_2$. The initial condition is $x_0$ = ${\left[\begin{array}{ccc}2300& 450& 5800\end{array}\right]}^T$ kg, which is a reduction of the initial condition of $x_2$. The disturbance is a 5% decrease in $x_1$, 20% increase in $x_2$ and 20% decrease in $x_3$ occurring between the 60th and 65th minutes where the system is in a steady state. The desired input is u = 0.95, which is in the unstable region. Figure 12 shows the states of the system.

The disturbance caused the states to violate the control zones in Figure 12, but the NMPC restored the states to their zones. State $x_1$ was restored in about, 11 min while the disturbance effect was removed about 67 min after the disturbance was introduced. It took 2 min to restore the $x_2$ to the zones after the disturbance ends, while $x_3$ took 40 min. The disturbance effect was removed after 87 and 77 min for $x_2$ and $x_3$, respectively.

The input target ($u_{des}$) of 0.95 is not met as the optimal input sequence reached a steady state value of 0.901. The optimal production attained a steady state value of 8.78 kg/s, which is lower than 9.38 kg/s for the undisturbed system. On simulating using $u_{des}$ = 1.25, which is out of the input bound and unreachable, the optimal input from the controller attained a steady value of 0.92, while the oil production steadied at 9.09 kg/s, which again is less than the undisturbed system. This shows that the controller has been able to restore the states to their zones but could not attain the optimal production for the undisturbed system.

4.4. Disturbed Gas-Lifted System Stabilisation Using Terminal Equality Constrained NMPC Having Input within Bound: Two Input Case

Most control solutions to operate a gas-lifted system use only one input (the production choke) to control the flow through the system. The flow into the annulus ($w_{gl}$) is usually fixed by a regulatory controller (mostly PID), and the set point depends on the available gas from the compressor station. In a few cases, both valves are manipulated to control the flow in a multivariable approach for optimal operation, as shown in [30,37]. We examine the effect of the additional degree of freedom on disturbance attenuation and improvement of the oil production rate. All the values in the one input case are retained, but $w_{gl}$, which is fixed in the one input case, is here an input where the desired value is given as $u_{des}$ = 0.4 kg/s. Further, the simulation time is 400 min to give room for the inputs to converge to the desired values. Figure 13 shows the two inputs, and Figure 14 shows the corresponding cost function.

The introduction of an extra degree of freedom makes the optimal input for the valve opening in Figure 13 slightly closer to the input target than the one input case. For a longer simulation time, the input target will be met for both inputs. The cost function in Figure 14 declines to almost zero after the arrival of the disturbance. An increase in simulation time showed that this offset was removed after the disturbance. The production rate improved slightly to 8.87 kg/s over the one input case, which is 8.78 kg/s. This is caused by the optimal input for the valve opening reaching a final value of 0.93 kg/s compared to 0.901 kg/s for the one-input case. There is, therefore, just a very small advantage in using the two-input case.

The state plot is almost identical to Figure 12; hence it is not shown here. This is due to the controller being able to fix the second input at almost the constant value of u = 0.4 kg/s in Figure 14, which is the value of the $w_{gl}$ in the one-input case. The first optimal input, therefore, approaches the optimal input corresponding to Figure 12, hence the state trajectories are almost identical.

If, however, the unreachable input is considered, the advantage of the two-input case over the one-input case in terms of disturbance attenuation and improved flow of oil becomes obvious. In under 20 min, the unreachable input targets forced the two inputs to converge to the upper control limits, as shown in Figure 15. Following the arrival of the disturbance, the optimal inputs remained at the upper boundary of the control limit, hence maintaining states $x_2$ and $x_3$ in their zones despite the disturbance. State $x_1$, however, violated the state zones briefly and was quickly restored to the zones by the saturated inputs, as shown in Figure 16. The steady-state production reached 12.25 kg/s, which is far more than 9.57 kg/s for the one-input case. This is expected as the rate of production, and the stability of the gas-lifted system depends on the flowrate of the gas into the annulus, which is saturated at 1.0 kg/s.

4.5. Compare with PI Controller

In Section 2.4, it is observed that the gas-lifted system presented here is highly nonlinear, necessitating the use of a nonlinear control approach for its optimisation. We compare the performance of the end-constrained NMPC with the PI controller described in [22]. Balanced rather than aggressive tuning [38] of the PI controller is used here since the gas lift system is a very slow system. The controller parameters were selected based on simulation results that provided the best optimal oil production. The parameters selected are $K_p=10,K_i=0.1$. The system is linearised around an operating point corresponding to $u=0.70$, which is in the stable region based on Figure 5. The output is the downhole pressure whose measurement is available; hence there is no need for estimation. The setpoint for the downhole pressure is selected such that the steady-state input is $u=0.95$. Figure 17 shows the gas-lifted states following the application of the PI controller, while Figure 18 shows the PI controller output (gas-lifted system input).

The controller stabilised the gas lift states in (Figure 17) at a steady state after damping the oscillation. This is achieved by the PI controller gradually opening the valve further from the equilibrium value of $u=0.70$ to the steady state value of $u=0.95$, as shown in Figure 18a. The oil production increased slightly to 9.04 kg/s in (Figure 18b), which corresponds to a meagre 1.08% increase over the open-loop mean production rate. Choosing a lower setpoint for the downhole pressure can increase the oil production rate, but this might lead to computed controller output being different from input to the gas-lifted system since the valve cannot handle an input of around 1 (100%). Further, since the objective of casing-heading instability removal is to increase mean production and protect downstream equipment, using a PI controller is still better than open-loop operation despite the slight difference above. There is also a meagre 3.76% increase in the production rate when NMPC is used compared to the PI controller. Despite this value being small to compensate for the rigour of using the NMPC, the fact that NMPC incorporates constraints in its formulation makes its use better than the PI controller.

4.6. Comments on Validation

A key issue with the controller solution to gas lift optimisation is the experimental validation. While the experimental validation increases the fidelity of the control solution, it is limited by the available experimental setup. Most controller solutions, therefore, stop at the theoretical results of the designed controller on the gas-lifted system, as shown in [4,10,11,13,14,15,16,17,18,21], among others. A few other cases show the experimental validation in addition to the theoretical results of the controller, as seen in [9,12]. While the scope of this work is limited to the theoretical results, below we briefly compare the results of this zone-controlled NMPC for gas lift optimisation with existing results for increased outflow.

In [15], a highly oscillatory gas-lifted system was stabilised by a proportional integral (PI) controller. An improvement of 150% was obtained. In the one-input case in this article, an improvement of 5.63% in oil production was obtained when NMPC was used compared to the open loop. While this improvement is large compared to the result in this article, achieving the setpoint supplied to the PI controller will require aggressive valve opening compared to the NMPC. This is because the NMPC incorporates constraints on valve opening. In [22], where NMPC with two inputs was used to optimise the gas-lifted system, the oil production increased by 36%, which is in the region of 37% obtained in this article using the NMPC with two inputs. This improvement depends on the amplitude of oscillation of the system in casing-heading instability in addition to the control parameters selected and disturbance in the system.

5. Conclusions

End equality-constrained NMPC is used to stabilise a gas-lifted system. The stability and convergence of the controller as well as some behaviours of the system, were first discussed. The controller stabilised the gas-lifted system for one input (production valve) and two inputs (production valve and lift gas injection valve). The two-input case has advantages over the one-input case when the input targets are out of input bound as the controller forces the inputs to saturate at the upper boundary of the inputs. The saturated inputs helped attenuate the disturbance much quicker and also improved the flowrate of the produced oil. The NMPC performs slightly better than the PI controller both in terms of a slight increase in production rate and speed of damping of the oscillation. The key benefit of this approach compared to the use of linear MPC is that the prediction horizon is better optimised by both the use of a nonlinear predictive model and the use of artificial reference states. In addition, the zone control approach ensures the controller focuses on optimising the input when states are in their zones, unlike set point tracking applications, where control efforts are on ensuring that the given reference is tracked.

Future work can look into solving the gas-lifted problem optimisation in CasADi using ipopt solver, which uses multiple shooting techniques for the evaluation of the ODE. Validation of the gas-lift experimental setup will be looked at in the future.