Identification and Optimal Control for Surge and Swab Pressure Reduction While Performing Offshore Drilling Operations

In this paper, an unscented Kalman filter (UKF) coupled with a nonlinear model-predictive controller (NMPC) for a hydraulic wellbore model with multi-variable control and tracking is presented. In a wellbore, high drill string velocities in operational sequences such as tripping might result in surge and swab pressures in the annular section of the wellbore. To overcome these challenges, a controller incorporating safety and actuator limits should be used. A second-order model is used to predict axial drill string velocity downhole. With a NMPC specifying the block position trajectory, choke flow reference, desired backpressure pump flowrate and stand-pipe pressure, we can automatically supervise and control the pressure in the wellbore. To compensate for unmeasured states, an estimator is designed to predict the frictional pressure forces in the wellbore and filter noisy measurements. A stochastic approach for the hydraulic model is taken, including variance of the average fluctuations for the flow and pressure states. Comparing three NMPC configurations, the result of using an integration of the tracking error in the prediction model gave best offset-free tracking of the bottom-hole pressure. The controller compensates for the unknown fluctuations, and is shown to be robust towards model mismatch. Including the mechanical system in the NMPC prediction model, we can effectively constrain the predicted axial drill string velocity to reduce the pressure oscillations and achieve tracking of bottom hole pressure and choke differential pressure. The outcome is shown through extensive simulations to be an effective control strategy, reducing the pressure spikes while tripping.


Introduction
Drilling an offshore well is comprised of sequences to be executed in a safe and efficient manner to reduce pressure fluctuations in the wellbore. Tripping, either running-in or running-out sections of drill string (DS) pipe, are done to extend or shorten the DS assembly while drilling a well. In these operations, typically the cost due to time is larger than the production costs, requiring the highest rate-of-penetration possible. Furthermore, increasing tripping speed might lead to instability in the wellbore known as surge, and swab pressures occurring in the annular section (Rasmussen and Sangesland, 2007;Lyons et al., 2015). An offshore drilling process is illustrated in Figure 1.
Automatic pressure control is a measure of stabilizing the transient pressure in the wellbore. Implemented in the process control of the drillingrig, the conventional Proportional-Integral-Derivative (PID) controller is used to control the subsea or topside choke valve, creating a back-pressure in the well to effectively control the bottom hole pressure (BHP) (Gravdal et al., 2018). However, as pointed out in Gravdal et al. (2018), changing conditions in the wellbore due to well geometry, open-hole sections, formation and fluid properties, along with changing temperature profile along the well can limit the PID controller in terms of required re-tuning during the operations. These quantities are in many cases dependent on the movement of the DS, and it is emphasized that the control law should oppose limitations on the mechanical side of the drilling operation (Cayeux et al., 2014). Application of automatic control to obtain reference tracking of the bottom-hole pressure (BHP) with tight pressure margins is commonly referred to as Managed Pressure Drilling (MPD). On the other hand, for wider margins the fluid density is typically adjusted to maintain desired well pressure (Nygaard et al., 2007b).
The challenges of maintaining a stable wellbore are connected to the reliability of the measurement data and the physical model used for real-time wellbore state prediction (Cayeux et al., 2014). Sensors used for well control are mud-pulse telemetry, wired drillpipe transmitting pressure measurements to the surface, and the gyro in the bottom-hole assembly (BHA). Mud-pulse is restricted to cases where sufficient flow is maintained. Currently, drilling-rigs are equipped with higher fidelity sensor packages, such that more advanced control systems can be utilized (Cayeux et al., 2010). As such, a broad research field on soft sensors is available in terms of observer-based applications (e.g Stamnes et al. (2008)) and the use of Kalman filters (KF) (e.g. Nygaard et al. (2006); Gravdal et al. (2010)). Filtering techniques are methods to provide additional insight, parameter estimates of unobserved process variables with minimum variance and bring redundancy in the measurement data.
Control design in MPD applications traditionally restrict to the use of the choke valve (Meglio and Aarsnes, 2015). The choke valve is the variable restriction in the mud return flow from the annulus. A rotating circulation device is included at the top of the well to seal off the annulus between the DS and borehole wall (Downton, 2012). However, an advantage is to combine the choke with the mud circulation system and draw works to further increase the capabilities when tight pressure margins are present and longer reach wells are drilled (Godhavn, 2009). In terms of choke control, see for example Nygaard et al. (2007a); Stakvik et al. (2016); Zhou (2018).
In multi-variable control, a supervisory control system is included to set reference points to sub-level controllers (typically PID) directly actuating valves, pumps, etc. In MPD applications for controlling pressures at defined locations, model-based schemes are commonly applied. Typically, a first-order model com-  The variable v s denotes the tripping speed. Green dashed lines mark the boundary of each segment in the model. prised of ordinary differential equations is sufficient for capturing the transient pressure and flow effects (Kaasa et al., 2012;Gjerstad et al., 2013). Linear model-based control schemes have been studied extensively in Nygaard et al. (2007b);  ;Breyholtz et al. (2010) and Møgster et al. (2013). The latter utilizes the WeMod high fidelity well simulator with Equinor's SEPTIC Model Predicitive Control (MPC) software. In these studies, the DS velocity is manipulated directly.
In terms of nonlinear multi-variable control, the nonlinear MPC is using a nonlinear model to calculate the process inputs. A comparsion between a PI controller and a NMPC using a low-order model is given in , were the control design was additionally tested on a high-fidelity dynamic model used in offshore drilling. In Pedersen et al. (2018), choke pressure, pump flow and the separator are controlled to achieve multi-objective control with no DS dynamics. The work considers underbalanced drilling which is a variant of the MPD, allowing reservoir production during drilling. The work of Nandan and Imtiaz (2017) includes an NMPC for switching between tracking of BHP and kick-attenuation in the wellbore, showing promising results by using the choke valve for control when performing pipe connections.
In this paper, we combine both the mechanical and fluid domains by means of a hoisting model for the axial DS dynamics, along with a first-order lumped parameter model with first-principles mud flow dynamics from Gjerstad et al. (2013). The frictional pressure forces in the wellbore and inside DS are inherently coupled with the DS velocity in this model. To effectively compensate for pressure transients and achieving desired wellbore conditions, we develop a supervisory control system in terms of a NMPC, with state and parameter estimates supplied by an UKF. We assume sparse knowledge of the frictional pressure forces in the annulus, and show the performance of the suggested control strategy through extensive simulations.
The rest of the paper is organized as follows: Section 2 gives the preliminaries of the methods and theory used in this paper. Section 3 gives an overview of the hydraulic wellbore model which is extended with dynamics of the hoisting system. Section 4 includes a survey on the aspect of system identification for the wellbore when the flow is limited through the bit nozzle, the estimator design and a test case. Section 5 presents the control design, with nonlinear model-predictive control. Section 6 shows the result of the work through extensive simulation case studies with Section 6.3 giving a discussion on the results. Section 7 gives the concluding remarks on the work.
The nomenclature used in the paper is summarized in Table 1.

Preliminaries
In this section, we present the preliminaries on modelling of lumped fluid flow in a pipe, the general stochastic differential equation with continuous Wiener processes, and finally an overview of nonlinear estimation in terms of the unscented Kalman filter. This is intended to give the reader some familiarity with the topics presented and the developed material in this paper.
Conservation of mass in a control volume is defined by the continuity laws, which are given as (Egeland and Gravdahl, 2002) where V c is the volume, β is the bulk modulus, p is the pressure, q 2 is flow into the volume and q 1 is the Root-mean-square-of-error flow out of the volume. Using dρ = ρ β dp, and assuming equal density at the inlet and outlet of V c , we consider positive flow direction upwards in the vertically oriented control volume.
Equivalently, we can express the momentum balance in terms of the rate-of-change of flow rate as a function of the net fluid pressure in the control volume (assuming that pressure is uniform in the volume) (Kaasa et al., 2012) where M is a fluid constant, p f is the frictional pressure losses, p 1 is the upstream pressure in the control volume and p 2 is the pressure downstream of the control volume.

Stochastic differential equations
In general, many systems express a stochastic nature and therefore need to be treated with a deterministic and a stochastic part. A stochastic differential equation (SDE) in general form is given as where f is the deterministic part of the SDE, g is the diffusion term, θ is a system parameter vector, v(t) is a standard Wiener process (Brownian motion) depending on the time interval t ∈ [t 0 , t f ]. The difference v(t) − v(s) is normally distributed with zero mean and variance σ 2 = I(t − s). Then, defining h as the time step h = (t f −t 0 )/N and N is the number of increments from t 0 to t f , we have . . . , N (4) where j is the number of noise variables, dv k ∼ N (0, √ h) and t k = k h. 1 To integrate the second part, including the diffusion term, an approximation of the stochastic integral is used due to the nondifferentiable state of dv (the term is not related to any past or present values) commonly either the Itô or the Stratonovich form (Chirikjian, 2009).
The SDE in (3), is then approximated numerically by a Taylor series expansion of defined order (here order 1, Euler or Runge Kutta-Maruyama) to simulate the stochastic behaviour.
The discrete-time measurement is formulated as where y k is the sampled measurement values, h is the observation model, w k ∼ N (0, R) is the zero-mean discrete-time measurement noise and R is the noise covariance.

The unscented Kalman filter
Physical states are in many cases not observable directly from the measurements. The system identification tool to obtain the system states or parameters (or combined) are the estimator. Estimation schemes applicable to nonlinear dynamics include the Extended KF (EKF) and the UKF. The EKF involves linearization of the system and observation model, to propagate the state and error covariance in time. However, divergence of the filter estimates can occur if the linearized model poorly represents the actual model and when large steps out from the linearized point (x * , u * ) are taken (Brown and Hwang, 2012, Chapter 7). The linearization is performed as where φ is the linearized system transition matrix, H is the linearized output mapping matrix. To overcome the challenges, the UKF was introduced in Julier and Uhlmann (1997). Using a nonlinear transform, the UKF estimates the system probability density function through a deterministic, minimal set of sigma points (SP). The use of SPs enables better approximation of the true mean and covariance by using a 2nd-order approximation, unlike the 1st-order approximation of the extended KF.
In this section, we consider the discrete-time dynamics of x k and additive noise in the system. The notation x k|k−1 denotes the current sample given last sample time information about the mean. The time-update starts with a new draw of estimator SPs, X k , calculated according to the initial estimated mean ofx k−1 from last sample-time t k as where U = chol(P) is the Cholesky factorization of the state covariance matrix 2 , col i (U) is the ith column and n x is the number of system states. The weights, determining the impact of each SP state is given as where ω i µ is the mean weights, ω i P is the covariance weights, λ = α 2 (n x + κ) − n x , α determines the spread around the mean, β = 2 assuming Gaussian distribution ofx k , and κ = 3 − n x is the scaling factor (Brown and Hwang, 2012).
The predicted mean x k and covariance P k are computed based on the nonlinear transformed stochastic variable, expressed as wheref is the nonlinear discrete-time state transition function, u k is the discrete input, p = 2m x + 1 and Q is the UKF covariance matrix. Furthermore, the predicted observation and its covariance and the resulting cross-covariance are computed according to whereh is the (nonlinear) discrete-time observation model, y i k is the predicted sigma point measurement, y k is the predicted, weighted measurement, P yy,k is the observation covariance, P xy,k is the cross-covariance and R is the UKF measurement covariance. When a new measurement sampling from the process sensors is obtained, a filter measurement update is performed. The a posteriori estimates are given by the Kalman filter update equations (Haug, 2012) K k|k−1 = P xy,k P −1 yy,k (14) where K k is the Kalman gain at time t k ,P k+1 is the a posteriori covariance estimate of x k ,x k+1 is the a posteriori estimate of the system states, and y k is the measurement obtained at t k . During each measurement sample-time interval, the filter prediction is updated using the current input and last iteration predicted mean x k−1 .

System modelling
Referring to Figure 1, the wellbore is discretized into a 2n − 1 degrees of freedom fluid dynamic model. The boundaries are q 1 = q c and the injection flow rate from the bit nozzle. The segments are explained as follows: • segments 1 to n − 1 consist of drill pipes, • the last wellbore segment is the BHA, consisting of drill collars, the cutter, and various logging tools, • positive flow direction is defined upwards in the annulus, • the lumped control volume pressure p i is uniform in each segment.
The wellbore representation is drawn from (Gjerstad et al., 2013), which derived a discretized hydraulic model where the frame of reference is fixed to the wellbore formation, such that the movement of the DS is assumed to alter the volumetric flow rate for the two lowermost segments. These consist of the largest geometrical changes of the DS assembly, causing larger flow variations when movement occurs.
The model was also derived considering two different diameters in each wellbore segment. This corresponds to two different flow rates and pressure states at these points in the annulus, which was included to approximate the pressure loss in terms of both laminar and turbulent flow in the annulus and DS. To reduce the number of states in the model, being derived from (1) and (2), the pressure forces for the main and secondary sections are lumped together. The main section (index 1) is either the drill pipe or collar, and the secondary section (index 2) is either the tool-joint or BHA. For annular flow, the secondary flow rate is the difference between the main section flow rate and the portion following the moving wall, yielding where A f 1 and A f 2 are the cross-sectional flow areas. The averaged flow velocity in each section of the annular segment is then given as The model in this paper includes the abovementioned properties, and we extended it to include the choke flow rate and the dynamics for the DS. We consider the case of no influx from the reservoir and that the wellbore is closed down-hole, i.e., the last control volume is closed. The normal forces of the fluid is then assumed to cancel the gravitational forces (also in the case of inclination), for each control volume. The length of the wellbore is assumed to be fixed during the time instant of tripping, such that we do not consider the extension due to drilling (i.e., the number of segments is fixed during operations).

Conservation of mass and momentum in wellbore
From Figure 1, the mass and momentum balances for each segment can be expressed as where i ∈ {2, . . . , n}, γ and are conversion factors (Pa and m 3 ·s −1 to bar and l·min −1 , p = γp, , q = q) 3 , V i is the segment volume, β i is the bulk modulus of the segment, p i segment pressure, q i is the flow rate out from the segment, q 1 = q c is the choke flow, q bpp is the input flow rate from the back pressure pump unit topside,V i is the volume change due to surge and swab effects of a moving DS, M i is a fluid constant, F f,i are the frictional pressure forces, A eq,i is the equivalent area for the acting pressure forces, and q v is the flow out from bit nozzle valves into the annular segment n, which will be further explained in Section 3.2. The flow rate out of segment 1 is equal to the choke flow through the valve and is then subject to the pressure downstream, i.e. the pressure in the mud tank (being atmospheric for open tanks). The valve flow is expressed as where u c ∈ [0, 1] is the valve control signal and C c is the lumped choke valve coefficient. The differential pressure over the choke is p 1 − p 0 , where p 0 is the atmospheric pressure in the mud tank. Drill string movement results in pressure fluctuations in the annulus. These pressure changes are largest in the region where the geometrical changes of the DS are largest, specifically when the DS geometry changes from pipe to collar and the BHA (Gjerstad et al., 2013). This is due to the fluid velocity being larger when the annulus volume in the collar/BHA section is small. To account for the increase/decrease in rate of change for the pressure in the control volume for the last two segments, the time derivative of V is expressed aṡ where v s is the axial DS velocity, A c is the collar crosssectional area and A p is the pipe cross-sectional area. The pressure forces in each segment are expressed with the equivalent flow area, defined as where d w,i is the wellbore diameter and d k,i , k ∈ {p, j, c, b} is the outer string diameter, at segment i. The constant M i is expressed as where ρ m dx = dm 1,i is the infinitesimal fluid density per length in each section of the segment integrated for the flow path dx in the well. Furthermore, we assume uniform density in the sections, and according to (Kaasa et al., 2012), the parameter M i is approximated in lumped hydraulic systems. Hence, in (26) In Gjerstad et al. (2013), the frictional pressure forces are derived with Herschel-Bulkley fluid properties. These type of fluids reflect closely the properties of mud flow, being an approximation with fluid yield point related to the Bingham plastic and the power law model (Whittaker and EXLOG Staff, 1985). We will not repeat the derivations of the wall-shear stresses, and hence, the reader is referred to the work in Gjerstad et al. (2013).
The frictional forces arise from the wall shear stress from the fluid and is given as where A u1,i is the boundary surface between mud and the surrounding borehole wall and DS in the main section, τ w1,i is the corresponding averaged shear stress value over the main section wall surface area, A u2,i is the boundary surface between the mud and the secondary section and τ w2,i is the averaged wall shear stress value for the secondary section, v e is the effective flow velocity in the annular segment being the sum of velocity for the moving component in the control volume and actual flow velocity. The two components τ w1,i and τ w2,i are approximated by the laminar and turbulent flow regime for the mud flow with moving wall. The transition between these to regimes is modelled as where ii= {1, 2}, and f tr is a transition function depending on the equivalent Reynolds number for Non-Newtonian fluid.

Conservation of mass and momentum in the drill string
The inner volume of the DS is presented in Figure 2. The inner volume of the DS is lumped into a single control volume, assuming a uniform pressure inside the string (Kaasa et al., 2012). The averaged flow velocities relative to the moving DS are given as Similar to the dynamics of a control volume in the annulus, the inner volume of the DS is comprised of the flow rate and pressure state, given as where subscript I is for inner, the change in pressure from of DS heave is modelled withV I , is the main inner cross-sectional area and A eq,I = πd 2 p,i /4 is the equivalent area for the inner section with d p,i being the drill-pipe outer diameter for segment i. The inner frictional forces are defined as where A u[1−2],I are the boundary surfaces of the fluid flow for the main and secondary sections, and τ w,I1−2 are the averaged shear stress value. Furthermore, the transition functions defines the contribution for each wall shear stress section for the inner control volume, expressed as In the case of running mud through the DS, p in = p spp , where p spp is the stand-pipe pressure. When the stand pipe is disconnected, the pressure at the inlet of the DS is atmospheric.
The inner DS control volume is defined as a single section. The diameters for this segment are defined by d I,k , k ∈ {j, c, b} for the secondary section and d I,p is the main section diameter (1).
In rotary drilling for high deviation wells, mud-motor or BHA turbines are applied to achieve high rotation speeds at the BHA. This is done by high pressure pumps topside, forcing large amounts of drilling mud into the string, rotating the turbines in the BHA for increased bit rotation speed (Black et al., 1986). The total pressure drop over the BHA is defined by the pressure drop over the turbine in the BHA lumped with other restrictions. We can then express the pressure drop over the BHA as where K b is a constant. The nozzle contains check valves, such that the flow is directed only out from the DS to the wellbore. Thus, the pressure drop over the valve is then The flow out from the nozzles is then a conditional function, which can be expressed by the valve equation as where K v is the lumped nozzle valve constant.

Rig travelling block dynamics
To account for the interaction between the DS and travelling block, we assume that we can represent the dynamics of v s as a second-order mass spring damper model. The rig system is represented as a fixed element with an attached spring and mass, as seen in Figure 3. The sum of forces for the point mass m ds with sign convention equal to what is defined in the previous sections, the DS dynamics is derived as where x s is the perturbation in DS position, v s =ẋ s , m ds is the mass of the DS, k eq is the equivalent DS stiffness, u b is the position of the travelling block, c eq is the equivalent damping to the surrounding fluid and formation, g is the gravitational constant and F b is the buoyancy factor.
Using the transform ξ = x s − 1 keq m ds gF b and dividing by m ds , we can express (35) as where ξ is the deviation from the suspended equilibrium point, in tension due to the gravitational forces acting on the DS. We proceed to define the system s 1 = ξ, s 2 =ξ anḋ where ω n = k eq /m ds is the natural frequency of the system, and ζ = c eq /(2 m ds k eq ) is the relative damping factor. Increased damping due to increasing mud flow is not included in this analysis, as we assume a constant damping factor for the mechanical system. The reference trajectory controller for the travelling block can be derived to fully complete the rig dynamics. With a linear system representation for the rig dynamics, we choose common methods to define u b as where k 1 , k 3 and k 2 is the proportional, integral and derivative gains respectively constituting a common PID control law (Khalil, 1996). The gains are chosen such that the system P = A − BK is Hurwitz, and where A z = ∂f (s)/∂s| si,ss,u b,ss , and f (s) is given by (37). Then, following the Routh-Hurwitz Criterion for a stable characteristic polynomial det(Is − P) = 0 (Khalil, 1996, Chapter 12.4), the gains ensures the Hurwitz condition.
Inserting (38) into (37), the DS and travelling block model constitutes theṅ where z 1 = s 1 , z 2 = s 2 , z 3 = σ, and the supervisory control input is actuating the reference u r .

State-space formulation
The wellbore and DS model from Section 3 can be formulated as a 2n + 1 degrees of freedom state-space model, where the wellbore is modelled with n pressure states and n − 1 flow rate states, and the inner DS contains two states.
Let n x = 2n + 1, x 1 = [p 1 , q 2 , p 2 , . . . , q n , p n ] , x 2 = [q I , p I ] , u = [v s , q c , q bpp , p in ] be the state vectors and input vector, then the wellbore dynamics can be formulated in state-space form aṡ where the system matrices and vectors are derived as Similarly, for the inner DS we havė where the x 2 system matrices are defined as The rig dynamics is already in a state-space formulation, and the total system is then given by the vector

System identification
In this section, we address the system identification for the offshore drilling operations with special attention to running connections. This is due to having no circulation of mud flow when performing this operation, which implies that the nozzle flow through the bit is reduced or zero. We wish to obtain estimates of the system states to monitor the pressures in the discretized wellbore along with the wellbore flow. To design a state controller, we require full-state knowledge of the system.
The measured outputs when circulating mud are assumed to be y = {z 1 , ∆p, p n , q I }, and y = {z 1 , ∆p, p n } otherwise, where ∆p = p 1 − p 0 is the differential pressure over the choke. The observation model h during circulation is defined as where k is the discrete time index, and in the case of no circulation, p in = p 0 and h is reduced by removing the fourth row.
In terms of the proposed model for the drilling operations, linearization of (22) and (34) leads to a singularity. When the pressure difference is negative or zero in (34), the expression becomes invalid and requires special numerical treatment when evaluating the simulation model. Due to this, we design a state and parameter estimator with the UKF. The UKF uses the nonlinear model directly in time propagation of the mean at its covariance from (40) and (41).

Observability
We must ensure that the deterministic systems (39), (40) and (41) are observable with (42) at time t k when sampling is performed. Observability is a measure of the system property to reconstruct the state x k given an input-output map y k , u k . For a linear system, this implies that the observability matrix formed by the system and output matrices has full rank (see e.g., Chen (2013, Chapter 6.3)). For nonlinear systems this is not a straight forward procedure, and often reduces to the previously mentioned rank test for linear systems. The conditions for nonlinear observability is addressed in Hermann and Krener (1977); Kou et al. (1973) According to Hermann and Krener (1977), we can show locally weak observability by calculating the Lie derivatives 4 up to n x − 1, and check the rank of the resulting Jacobian.
The linear second-order system comprising the rig dynamics in (37) can be shown to be observable for the measurement z 1 . The hydraulic system has a potential observability problem, with too few linearly independent measurements compared to the number of estimated states in the x 1 , x 2 system. Therefore, before designing the estimator we need to ensure system observability with the limited measurements. Consider a wellbore DS system with n = 2 (two segments), annular frictional pressure force are assumed to be quadratic in flow rate f j = −k f j x 2 j , where j = {2, 4} and k f is a constant (Kaasa et al., 2012). Using u = [u z2 , (q c + q bpp ), p in ], we can arrange (40) and (41) 4 For a more detailed description of Lie derivatives, see e.g., Slotine and Weiping (1991). where , 0] , B 1 and B 2 are given as in (40) and (41) by removing the third column, e = [0,0, −B 3 , 0, −B I ] , and q v ∈ C 1 from (34) being a nonlinear function in x. The state x 3 then corresponds to p n and x 5 is the inner drill string pressure, p I .
The system in (43) is observable if the following input-output map exists and the Jacobian of H has full rank, i.e rank The system is nonlinear and due to this, the results only show local weak observability at the operating point x ss . Using the symbolic toolbox in MATLAB we see that the system in (43) is locally weakly observable for q v > 0 and that O does not have full rank for q v = 0.
In that case, the state x 4 = q I , x 5 = p I are nonobservable from the output y. This is reasonable, due to q v (x) connecting the inner DS system to the wellbore system. Hence, when circulation is off, the system is only locally weakly observable for cases x 5 > K b x 4 +x 3 .

State and parameter estimator design
The deterministic systems (39)-(41) can be used to formulate a continuous-discrete SDE state transition and observation model as in (3) and (5). This is done to include a diffusion term g, representing the model errors and fluctuations in the state x capturing random effects of the wellbore system.
The frictional pressure drop in the annulus is commonly estimated by running tests, fitting an n-th order polynomial to the data and using this in the model (Kaasa et al., 2012;Landet et al., 2012). Due to the uncertain nature of the wellbore dynamics, we choose to represent the first-principles friction model in (40) by a quadratic estimate of the pressure forces arising from the friction in the annular volumes. The estima-tion model is then given as where dθ i−1 = dv θi−1 , i = 2, 3, . . . , n x − 3, f z is the deterministic part from (39), g is the diffusion terms determining the average fluctuations in x, v are Wiener processes satisfying v(t) − v(s) ∼ N (0, Idt),F f,i is the annular segment estimated pressure force set up by friction withθ being the estimated coefficient for each annular wellbore section,v i = q i /A f 1,i and K c is the clinging constant for the fluid attached to the moving DS. We assume to have prior knowledge of K c , which can be calculated for laminar flow as where α i = d 1,i /d w and d w is the wellbore diameter, and d 1,i is the main diameter for the segment. The term between the brackets in (46) is referred to as the effective velocity (Whittaker and EXLOG Staff, 1985). The number of parameters to be estimated should correlate with the number of measurements or be driven by a persistent excitation signal. Hence, we assume full knowledge of the fluid behavior in the DS volume. The parameters are assumed to be random walk processes with equal distribution as the Wiener process. The measured outputs are given by the criteria for operation (from Section 4), given as where w = [w 1 , w 2 , w 3 , w 4 ] , ∼ N (0, σ 2 w ), σ 2 w is the measurement noise variance.
The discrete-time estimation model in (45) is obtained by using the Runge-Kutta Maruyama order 4 integration method from Rößler (2006), with step length equal to the plant simulation model.
The estimator propagates the stochastic state vector X i through the nonlinear process model. Following the notation defined in the preliminaries of Section 2.2, the filtering algorithm is summarized in Alg. 1. The UKF algorithm runs at the same frequency as the plant model. However, the UKF measurement update is in this work set at a different frequency. The measurement sampling time in the UKF is denoted t m .
The algorithm assumes additive uncorrelated white noise (E[wv ] = E[wx 0 ] = 0 and E[vx 0 ] = 0). This Algorithm 1 Unscented Kalman filter 1: Initialize:P 0 ,x 0 , Q, R 2: while Run-time: do 3: if mod(time, tm) = 0 then 4: Time-update:x k−1 ,P k−1 5: Compute a new set of sigma points X i k from (6) 6: for i = 0 to p do 7: Predict x i k in (7) 8: end for 9: Predict new weighted mean x k (8) 10: Compute the predicted covariance P k (9) 11: Return: x k , P k 12: else 13: Measurement-update: y k is sampled 14: Compute K k ,x k andP k from (14), (16) and (15)  15: Return:x k ,P k 16: end if 17: end while allows us to simplify the computations by not calculating additional sigma points for v k , w k to be projected through the nonlinear transform and observation model.

State estimation with no mud circulation
Consider a wellbore system with n = 4 number of pressure states (discretizing the wellbore into four segments). To confirm our predictions on observability, we sample the plant model in (45) at 100 Hz, runningin-hole (RIH) with circulation off. The wellbore and DS data are given in Table 2, based on Lyons et al. (2015, Ch. 4.6.4), and the scenario is similar to what is seen in Figure 5 and 6 in (Gjerstad et al., 2013). The estimator data and initial conditions are given in Table 3. The measurement noise variance is calculated as σ 2 i,w = (W i y i,nom ) 2 where σ w is the standard deviation, y i,nom is a chosen nominal signal value for the individual measurement and W i is corresponding signal noise amplitude.
The pressure at the top of the DS is set to atmospheric pressure, and the choke position is set to 50%. The imperfect model estimates resulting from driving the block 28 meters downwards, simulating run-in of the drill string in the wellbore.
The estimated pressure state and true values are seen in Figure 4. The pressure estimatep I does converge, but is seen to have a constant bias. It is seen that the inner DS flow rate is excited by the DS velocity and is closer to its true value.
The back-pressure pump kept is ramped up at t = 2 min to 400 l/min. The DS velocity is negative when RIH is performed, and the velocity returns to zero as the new pipe section is put in place. The estimated states are seen to converge close to their true values.  The variance and measurement noise of the system generates a rapidly varying frictional pressure force for the inner-DS, as seen in Figure 5. Constant deviations for the estimated frictional pressure forces from the true value are seen, in the rightmost of Figure 5 forF f,i . However, the estimates converge but with a significant bias.
From this study, some of the uncertainties in the process of RIH or making pipe connection to process disturbances are seen. To properly address the issue of a nonobservable x 2 system, in the specific case of no mud-circulation, we will have to treat q v as a disturbance in the control system when circulation is off.

Control design for drilling operations
This section presents the controller design for the tripping operations in offshore drilling. As discussed in Section 4, in the case of no circulation (such as running pipe connections), an estimator control design must be built around the assumption of a disturbance flow from the bit nozzle. We assume that we have all the available measurements (mud circulation is on) and concern this work with maintaining BHP pressure at reference, while pulling out or running in the DS. The main focus is then to control the BHP, actuating the travelling block, z 1 , back-pressure pump q bpp , stand-pipe pressure p spp and choke valve opening u c , keeping the actuator and process constraints intact. The block position is used for pull-out-of-hole (POOH) or RIH, maintaining the constraints on DS velocity in the wellbore. Multivariable control, with both ∆p and p n is presented in the last simulation case study. A natural choice of control law with the requirements discussed above is the MPC. The MPC is a widely applied process control law for both linear and nonlinear applications. It is an advanced method of optimal control with receeding horizon control, meaning that the measured output at sample-time t s together with current and future constraints on x, y and u are taken into account. For a thorough discussion on the topic of nonlinear optimal control, see e.g., Findeisen and Allgoewer (2002); Grune and Pannek (2011).
The MPC has become an industry standard in process control, with special focus on advanced application where process constraints are important (Qin and Badgwell, 2003). The MPC control algorithm uses a form of the actual plant model (either true or estimated), to predict future behaviour of the process and optimal inputs. At each sampling instant, process states are obtained (given by a estimator/observer and sensors) and the MPC utilizes a state prediction to form the optimal control inputū on an extended process horizon (Findeisen and Allgoewer, 2002). A disadvantage is that the MPC can be influenced by modeling errors (Grune and Pannek, 2011). Vital ingredients in an MPC algorithm is process and actuator constraints, a cost function and a model.

Constrained BHP Control
To constrain the inputs and system, we choose to design the control law around the MPC framework. We can either linearize and use the conventional linear MPC, or use the nonlinear model in predicting future process behaviour. The latter seems to be the better choice due to the nonlinearity of the bit nozzle and choke valves. Also, a nonlinear version will be robust on a broader range of operating points than the linearized MPC. The disturbance from DS heave and eventually the flow through q v is then included when predicting future optimal inputs for the process.
The NMPC is based on using a plant model to predict the future control inputs, leading the measurement values to their respective references. An optimal control problem (OCP) is defined, based on a nonlinear prediction model, system equality and inequality constraints. This involves either a continuous or discrete nonlinear programming problem (NLP), which the NMPC solves to obtain the optimal, constrained input. An important real-time criterion for the NLP is that the solution must be obtained before the next sampling time occurs. For each future sample-time, the input is open-loop applied in the prediction model x j+1 = f (x j , u j , θ k ), in a finite prediction interval.
The NMPC uses the discretized augmented deterministic nonlinear model in (45) to predict future states and outputs. The discrete-time OCP for the NMPC is defined as (in the following, the tilde notation on x is removed for convenience) where j is sample interval index, N p = T p /t s and N c = T c /t s is the discrete prediction and control intervals respectively, T p is the prediction horizon, T c is the control horizon, t s is the MPC sampling time, W c is the set-point weight matrix, ξ j are the slack variables, T is the slack variable weight matrix, R c is the weighting matrix for the slew in control input, R r is the input reference trajectory weight matrix, and ∆u j = u j − u j−1 . The subscript c is used on the MPC to distinguish its variables from those of the UKF. The state, output and input constraints are defined as The hard constraints are in this case the upper and lower bound on the inputs. The rest are subject to slack variables. The prediction in f (x j , u j , θ k ) is given by (45). The NLP is then implemented as a sequential approach, using finite parametrization of the control inputs (Findeisen and Allgoewer, 2002).
The following additional constraints are defined: where the subscripts, L, H are the low and high band control limits, respectively. In addition to constraining the input change on the travelling block movement, the first constraint in (52a) denotes the maximum and minimum DS velocity, which is calculated from (39). Rapid changes in the travelling block might result in large pressure spikes. In (52b) the requirement is to keep the pressure in the DS higher than the turbine loss and BHP pressure, such that flow spikes into the wellbore are avoided while circulating. The slack variables ξ > 0 define soft constraints in (52a)-(52e). These are included to allow minor violation in short intervals of the constraints on x j , such that a feasible solution is found. However, it is important to enforce hard constraints on actuator limitations. The band control limits on the BHP and choke differential pressure can be set in (52d).
A PI-controller is used to adjust the choke valve opening. The NMPC sends the desired choke flow rate to the controller, adjusting the choke to comply with the requirement. To constrain the flow control algorithm, we implement a nonlinear constraint on the choke flow in the NMPC. The flow constraint given as The discrete form of (49) and (50) are approximated by repeated forward Euler integration on the interval [t k , t k + N p ]. An intermediate step-length for integrating the differential equation in (45) is implemented, such that the sampling interval from start to end is divided into τ s = t s /M . The scheme is formulated in Alg. 2. Similarly, a Runge Kutta method, or any other discretization scheme, can be used. Exact discretization is only applicable if linearization of (45) is done prior to discretization. The input u j is held constant on the intermediate integration interval.
Algorithm 2 Euler-discretization of f 1: Input: xa, u, θ, t k , ts, M 2: Set step length τs = ts/M , M : Number of iterations 3: x 1 0 = xa 4: for k = 0 to M − 1 do 5: Solving (49) gives the optimalū k for the process, based on the initial guess. The initial guess is effectively drawn from the previous sample time control sequence. The NMPC is implemented in MATLAB, using fmincon with SQP. The NMPC algorithm can be summarized as in Alg. 3.
In (50), we utilize the forward Euler integration for the prediction model as described in Alg. 2. We choose to distinguish between the prediction and control intervals, and set N c = T c /t s . This implies that we use the linear equality constraint, defined in (51).

Simulation study
The wellbore and DS system together with estimator and controller, were implemented in MATLAB. The primary goal of the study is to reduce surge and swab pressures while performing tripping, along with multitarget tracking of pressures in the annular wellbore. As seen in Figure 4, the surge pressures while RIH cause rapid fluctuations in the pressure states. The target is given as p n,L ≤ p n (t) ≤ p n,H where the upper (H) and lower (L) bounds can represent the bottom-hole fracture and pore pressures, respectively.
As described in Gravdal et al. (2018), the system under consideration is limited in the sense of sensor measurements, time-delay, and under constant change due to geo-physical effects in the wellbore. We consider an ideal setting where measurements are not time-delayed and are aware that this effect may give a considerable performance decrease in the system identification and control of the pressure and flow in the wellbore.
The designed estimator and control configuration performance is tested in extensive simulation, with four case studies. The initial conditions in each case study are equal. The common system parameters are shown in Table 4, where the auxillary system time constants (see Section 5.1) are T dw , T ch , T bpp and T spp , denoting the draw-works, choke valve, stand-pipe pressure, and back-pressure pump, respectively. All the case studies presented here are based on data from the test well in Table 2. The tuning for prediction, control horizon, sample-time and weighted setpoint tracking W c , R c are based on trial and error. However, we use scaling values, based on the nominal output values and the individual actuator range. Furthermore, tuning factors are used to adjust the weights given as where y i,span is the span of the output set-point values, u j,span is the span of the input (actuator incremental change) andw c,i ,r c,j are tuning factors for the diagonal elements of W c and R c , and R r = R c .
The simulation study is broken down in four case studies. The three first case studies compare the NMPC performance in three configurations, related to how model-errors and disturbances are handled. The first method is adopted from Kwakernaak and Sivan (1974, Chapter 3.7.2) for offset-free LQR, by augment-ing the state vector with the tracking error given as where e is the integral error and z ref , x 1,ref , x 2,ref are the desired state trajectory for the rig system, wellbore and drill-string states, respectively. The term e Q c e is included in (49), where Q c, [i,i] =q i W c, [i,i] , withq denoting the tuning factor. The second method includes a disturbance state estimator as proposed in Pannocchia et al. (2015), yieldinĝ where y k is the sampled measurement,d k is the estimated disturbance, and K x , K d are the Kalman gains for the estimate error. The estimator design is chosen depending on the application. The advantage of including integration of the disturbance is to achieve integral control by capturing the model-mismatch between the real and estimated plant dynamics (Borrelli and Morari, 2007). The disturbance state is kept constant over the prediction horizon in the NMPC. The NMPC for Case 3 uses a linear disturbance model as where A d and C d are chosen such that the states are observable, which applies according to restricting the dimensions dim(n d ) ≤ dim(n y ), where n d , n y are the number of disturbance states and available measurements, respectively. As in Section 4, we utilize the UKF, and at NMPC sample-time an estimate of the d k is given to the NMPC before solving the NLP. The NMPC-disturbance model setup is presented in Figure 7. Furthermore, we implement a reference trajectory specified by the user to predict smooth changes over the horizon in the NMPC. The trajectory is given as is the reference trajectory value from the last sampling instant, α j i ∈ [0, 1] is the smoothing coefficient, and i ∈ [1, 2, 3]. This allows for better model-mismatch compensation with a less aggressive controller (Qin and Badgwell, 2003).
A summary of the configurations are presented in Table 5. The three first case studies comprise of BHP reference tracking and ∆p band control while performing a POOH followed by a RIH sequence. Band control is initiated on ∆p by setting the second diagonal element in W c (and Q c , in Case 2) to zero. The NMPC and UKF parameters are presented in Table 6.
The choice of disturbance model for the NMPC is not obvious, and is left to the designer having knowledge over the process or measurement disturbances. The immediate disturbance for the hydraulic model is the fluctuations in g, and model-mismatch in the NMPC. We use an input disturbance model for the configuration in NMPC2, were the fluctuations in x 1,pn , i.e., the BHP state, is estimated. The choice is based on the disturbance from the DS heave on the BHP, which is counteracted by the flow rate from the drill-string in q v .
In Table 6, the NMPC prediction model is discretized with intermediate step length of 0.1 s. A less strict constraint is enforced on ∆p, which is allowed to vary between the lower bound of the BHP band, and ∆p L . The plots of Cases 2 and 3 are given in Figure 8.
For the simulations made in this paper, the biggest challenge of Case 1 was to compensate for the constant pressure change when tripping. From the top-most plot in Figure 8, NMPC1 is closer to achieve overall offsetfree tracking of the BHP. The extra term included in minimizing the cost-function is tuned to achieve faster disturbance rejection. The NMPC1 is seen to deviate from the two others in how the choke flow and DS pressure is controlled. Since tripping out pipe reduces the pressure in the wellbore, the NMPC1 requires Case 3: lower choke flow rate. To fulfill this, the PI flow controller reduces the choke opening in the first seconds of the sequence and the NMPC1 initiates an increase the stand-pipe pressure. The NMPC2 seems to require larger changes to the choke flow and DS pressure during tripping. All configurations are able to create a trajectory for the travelling block and meet the set-point. No noticeable requirements in flow rate were made to the backpressure pump, by the NMPC0-2, seen in the fifth plot.
NMPC1 including the integrated tracking error is quicker in the response to changes in set-points. A disadvantage of this can be larger BHP spikes. This is also seen from NMPC2 configuration. However, further tuning of Q c in the cost function can reduce this.
The estimated frictional pressure forces are presented in Figure 9. The estimated frictional pressure forces for the wellbore converges during steady-state and follows the trend of the true friction, seen for all the case studies. However, a noticeable bias is seen when rapid changes in the BHP is occuring.
Recall from (46) that we only assume to know the laminar klinging factor, and base the frictional law on a quadratic function. The true friction model is based on predicting the mud laminar and turbulent flow, with a transition rate. Hence, properties such as fluid yield stress is not accounted for in the estimated model. A result of that is that we get a bias, since we can get the correct magnitude and direction of the fluid effective velocity but must use θ i with large variance to achieve a value in acceptable range for allowing pressure increase in the NMPC prediction model. The root-mean-square values of the tracking error e =x − x ref (RMSE) for z 1 and p n , summarizes the closed-loop performance of the three NMPC configurations. The values are given in Table 7 for comparison.
A general observation from the values in Table 7, is that NMPC1 configuration performs better than the two others with the smallest overall tracking error.  However, the RMSE error is based on how quickly the controller meets the reference values, which disregards violation of critical constraints, overshoot limits, etc.

Multi-variable tracking in Case 4
The last case study comprises keeping both BHP and ∆p in their band and at their respective set-points through a sequence of set-point adjustments, with the block reference at zero during the entire time. We choose the NMPC1 configuration, which showed the best results in terms of RMSE. The parameters for the NMPC in Case 4 are summarized in Table 8. The noise parameters are the same as in Table 3. We set a steady-state input goal for the back pressure pump at 0 l/min, which is implemented by setting the third diagonal element of R r nonzero in the cost function. This allows the NMPC to focus more on using the  main mud pump and choke to control the down hole pressures. The tracking result in Case 4 is presented in Figure 10. Seen in plot 1-3 in Figure 10, the NMPC1 manages tracking of the three controlled variables, retrieving the block at zero position after using it for manipulating the BHP and ∆p. The velocity of the DS is well below the constraint margins, seen in the fourth plot of Figure 10.
The estimated frictional pressure forces and the θ i parameter are presented in Figure 11. The same covariance properties are used in the UKF for Cases 1-3.
Stepping through both the desired BHP and choke differential pressure values, NMPC1 manages to keep both in their respective bands by manipulating the stand-pipe pressure and choke flow. The RMSE of the tracking error in Case 4 is RMSE = [0.00881, 0.322, 0.70039], which for e 3 is similar to the results in Case 1-3 (Table 7).
Since increasing the BHP has a lag effect on the upper volumes in the wellbore, the pressure wave must be compensated for in the choke while meeting the setpoint. As also shown in Møgster et al. (2013), control of multiple locations in the wellbore is an important feature of efficient MPD.

Estimator performance
The variance of the estimated parameter θ i and variance of the error for the calculated frictional pressure forces e i = F f,i −F f,i are given in Table 9. Table 9: Cases 1 to 4, mean and variance-of-error for estimator forθ i ,f 1,i , i = [2, 3, 4].

NMPC tuning and constraints
For the NMPC configurations shown in this paper, the initial tuning has been based on trial and error with process span variables. The estimator bias due to predicted frictional pressure forces is a challenge in terms of increased prediction error over the horizon, and must be accounted for when performing the tuning of T p and M . A more accurate prediction model would allow for a longer NMPC horizon.
Comparing the three different designs, NMPC0-2, the nominal design in NMPC0 yields less tuning variables, and tracking gives acceptable results in our case studies. However, with knowledge of disturbances in the system (such as sudden inflow from formation) we can achieve a better design with NMPC1 and NMPC2 by the advantage of the integrated error states or estimation of the disturbance state.
A measure of improving trajectory control is to include a full reference path for the hoisting (which can be a smooth 2nd-order polynomial for example) such that the NMPC follows this precisely. Then, the velocity constraint can be removed and the controller only requires three degrees of freedom for keeping BHP at its reference. The implementation in this work used position control for the travelling block. In practice, the driller adjusts the block velocity, and the block travelling limits are set by the operator in advance (Lyons et al., 2015).

Controller implementation robustness
Robustness measures in terms of the implementation has been done, such as iterating over a fixed horizon until either a limit is met or the MPC converges. This is done in an ad-hoc way, to at least ensure that more than one iteration is performed by the optimization algorithm.

Simulation performance
An averaged of 16.4 s for each NMPC NLP iteration was logged during the simulation study (Case 1 NMPC0). The sampling time was 1 s, as such, the goal of convergence in that time is not met. The corresponding real-time factor (RTF), given as RTF = t actual /t sim , was 17.44 for a total of 10 minutes simulation time.
Optimal conditions for real-time applications is when RTF < 1. The NMPC does not achieve real-time performance for the presented algorithm. We are aware of better methods of solving the optimization problem, but as a first implementation MATLABs fmincon gives an easy way for organizing the NMPC algorithm.

Conclusion
In this work, we have presented a multi-purpose control system used to minimize surge and swab pressures while tripping in an offshore well, along with tracking of choke differential pressure and BHP. The wellbore fluid frictional pressure forces was parametrized by a quadratic function in flow-rate and DS velocity. The unknown parameters in this function were estimated using an UKF, which also estimated the system states used by a NMPC. The NMPC manipulated the inputs to the back-pressure flow and stand pipe pressure directly, set the reference trajectory for a linear state-integrator controlling the block positioning system, and supplied a PI choke flow controller with desired reference flow. The actuator dynamics for the block, choke, back-pressure flow, and mud pump pressure were represented by first-order linear low-pass filters.
We have presented the results from comparing three NMPC configurations where two of these included measures to improve offset-free tracking. The best result in terms of RMSE values was achieved by NMPC1, with integrated tracking error implemented in the cost function. All three NMPC configurations showed good results in terms of constraining the BHP inside the desired envelope and at the set-point, while tripping.
For multi-target tracking, the NMPC was configured with a shorter prediction horizon and longer control horizon with a steady-state input goal for the backpressure pump. NMPC1 was used in this case study. Controlling both BHP and choke differential pressure was accomplished by the controller.
Since we used an SDE with RK4 Maruyama in time integration, average fluctuations were included. For the configurations and case-studies presented here, the NMPC has shown to be robust towards modelmismatch and variations in the states due to the nature of the SDE.
Including the dynamics of a hoisting system we can predict and constrain the estimated DS velocity to both manipulate the BHP and prevent excessive surge and swab pressures. The system presented in this paper is then automated such that the driller is only required to supply the desired reference points for BHP, choke differential pressure and travelling block position.