Practical Stability Analysis of a Drilling Pipe under Friction with a PI-Controller

TL;DR

This paper analyzes the exponential stability of a drilling pipe with friction controlled by a PI controller, using advanced Lyapunov functionals for coupled ODE/PDE models, and demonstrates limitations in reducing stick-slip oscillations.

Contribution

It introduces a new Lyapunov functional approach for stability analysis of coupled ODE/PDE systems with PI control, addressing both linear and nonlinear torsional dynamics.

Findings

Exponential stability established for linear torsional model using LMI framework.

Nonlinear stick-slip oscillations are not mitigated by PI control.

Numerical simulations confirm the theoretical results and limitations.

Abstract

This paper deals with the exponential stability of a drilling pipe controlled by a PI controller. The model used leads to a coupled ODE / PDE and is consequently of infinite dimension. Using recent advances in time-delay systems, we derive a new Lyapunov functional based on an state extension made up of projections of the Riemann coordinates. Two cases will be considered. First, we will provide an exponential stability result expressed using the LMI framework. This result is dedicated to a linear version of the torsional dynamic. On a second hand, the nonlinear terms in the initial model, that generates the well-known stick-slip phenomenon is captured through a new stability theorem. Numerical simulations show the effectiveness of the method and that the stick-slip oscillations cannot be weaken using a PI controller.

Figures

Tables

Table 1. TABLE I: Parameters values and their meanings for the lumped parameter model taken from [ 13 , 32 , 41 ] .

Parameter meaning		Value
$I_r$	Rotary table and drive inertia	$2122$ kg.m2
$I_b$	Bit and drillstring inertia	$374$ kg.m2
$k$	Drillstring stiffness	$1111$ N.m.rad-1
${\lambda }_r$	Coupled damping at top	$425$ N.m.s.rad-1
${\lambda }_b$	Coupled damping at bottom	$23.2$ N.m.s.rad-1
$d_r$	Rotary table damping	$425$ N.m.s.rad-1
$d_b$	Bit damping	$50$ N.m.s.rad-1
${\gamma }_b$	Velocity decrease rate	$0.9$ s.rad-1
${\mu }_{cb}$	Coulomb friction coefficient	$0.5$
${\mu }_{sb}$	Static friction coefficient	$0.8$
$c_b$	Bottom damping constant	$0.03$ N.m.s.rad-1
$T_{sb}$	Static / Friction torque	$15145$ N.m

Table 2. TABLE II: Physical parameters, meanings and their values [ 1 , 40 ] .

Parameter meaning		Value
$L$	Pipe length	$2000$m
$G$	Shear modulus	$79.3\times {10}^9$ N.m-2
$E$	Young modulus	$200\times {10}^9$ N.m-2
$\Gamma$	Drillstring’s cross-section	$35\times {10}^{-4}$ m4
$J$	Second moment of inertia	$1.19\times {10}^{-5}$ m4
$I_B$	Bottom hole lumped inertia	$89$ kg.m2
$M_B$	Bottom hole mass	$40000$ kg
$\rho$	Density	$8000$ kg.m-3
$g$	Angular momentum	$2000$ N.m.s.rad-1
$h$	Viscous friction coefficient	$200$ kg.s-1
${\gamma }_a$	Distributed axial damping	$0.69$ s-1
${\gamma }_t$	Distributed angle damping	$0.27$ s-1
$\delta$	Weight on bit coefficient	$1$ m-1

Table 3. TABLE III: Normalized parameters.

Parameter	Expression	Value
$\varphi (x,t)$	$\Phi (xL,t)$	-
$\psi (x,t)$	$\Psi (xL,t)$	-
${\tilde{c}}_t$	$c_t{L}^{-1}$	$1.57$
${\tilde{c}}_a$	$c_a{L}^{-1}$	$2.5$
${\alpha }_1$	$\frac{GJ}{L{I}_B}$	$5.3$
${\beta }_1$	$\frac{E\Gamma }{L{M}_B}$	$8.75$
${\alpha }_2$	$I_B^{-1}$	$1.12\cdot {10}^{-2}$
${\beta }_2$	$\frac{\delta }{M_B}$	$2.5\cdot {10}^{-5}$
$\tilde{g}$	$\frac{g}{GJ}$	$2.1\cdot {10}^{-3}$
$\tilde{h}$	$\frac{h}{E\Gamma }$	$2.86\cdot {10}^{-7}$
${\tilde{u}}_1\left(t\right)$	$\frac{1}{GJ}{u}_1\left(t\right)$	-
${\tilde{u}}_2\left(t\right)$	$\frac{1}{E\Gamma }{u}_2\left(t\right)$	-

Table 4. TABLE IV: Estimated decay-rate as a function of the order N 𝑁 N used. Note that μ m ​ a ​ x = 7.73 ⋅ 10 − 3 subscript 𝜇 𝑚 𝑎 𝑥 ⋅ 7.73 superscript 10 3 \mu_{max}=7.73\cdot 10^{-3} is calculated using ( 23 ) for k p = 10 − 3 subscript 𝑘 𝑝 superscript 10 3 k_{p}=10^{-3} , k i = 10 subscript 𝑘 𝑖 10 k_{i}=10 .

Order	$N=0$​	$N=1$​	$N=2$​	$N=3$​	$N=6$​
$\mu$ $(\times {10}^{-3})$	$0.87$	$4.24$	$7.31$	$7.59$	$7.73$

Equations

$$\left\{\begin{array}[]{l}I_{r}\ddot{\Phi}_{r}+\lambda_{r}(\dot{\Phi}_{r}-\dot{\Phi}_{b})+k(\Phi_{r}-\Phi_{b})+d_{r}\dot{\Phi}_{r}=u_{1},\\ I_{b}\ddot{\Phi}_{b}+\lambda_{b}(\dot{\Phi}_{b}-\dot{\Phi}_{r})+k(\Phi_{b}-\Phi_{r})+d_{b}\dot{\Phi}_{b}=-T(\dot{\Phi}_{b}),\end{array}\right.$$

$$\begin{array}[]{l}T(\theta)=T_{l}(\theta)+T_{nl}(\theta),\\ T_{l}(\theta)=c_{b}\theta,\\ T_{nl}(\theta)=T_{sb}\left(\mu_{cb}+\left(\mu_{sb}-\mu_{cb}\right)e^{-\gamma_{b}|\theta|}\right)\operatorname{sign}(\theta).\end{array}$$

$${\Phi }_{tt}(x,t)=c_t^2{\Phi }_{xx}(x,t)-{\gamma }_t{\Phi }_t(x,t),$$

$${\Psi }_{tt}(x,t)=c_a^2{\Psi }_{xx}(x,t)-{\gamma }_a{\Psi }_t(x,t),$$

$$GJ{\Phi }_x(0,t)=g{\Phi }_t(0,t)-u_1(t),$$

$$E\Gamma {\Psi }_x(0,t)=h{\Psi }_t(0,t)-u_2(t).$$

$$GJ{\Phi }_x(L,t)=-I_B{\Phi }_{tt}(L,t)-T\left({\Phi }_t(L,t)\right),$$

$$E\Gamma {\Psi }_x(L,t)=-M_B{\Psi }_{tt}(L,t)-\delta T\left({\Phi }_t(L,t)\right),$$

$$\left\{\begin{array}[]{l}\phi_{tt}(x,t)=\tilde{c}_{t}^{2}\phi_{xx}(x,t)-\gamma_{t}\phi_{t}(x,t),\\ \phi_{x}(0,t)=\tilde{g}\phi_{t}(0,t)-\tilde{u}_{1}(t),\\ \phi_{t}(1,t)=z_{1}(t),\\ \dot{z}_{1}(t)=-\alpha_{1}\phi_{x}(1,t)-\alpha_{2}T(z_{1}(t)),\end{array}\right.$$

$$\left\{\begin{array}[]{l}\psi_{tt}(x,t)=\tilde{c}_{a}^{2}\psi_{xx}(x,t)-\gamma_{a}\psi_{t}(x,t),\\ \psi_{x}(0,t)=\tilde{h}\psi_{t}(0,t)-\tilde{u}_{2}(t),\\ \psi_{t}(1,t)=y_{1}(t),\\ \dot{y}_{1}(t)=-\beta_{1}\psi_{x}(1,t)-\beta_{2}T(z_{1}(t)),\end{array}\right.$$

$$\left\{\begin{array}[]{ll}\phi(x,0)=\phi^{0}(x),&\phi_{t}(x,0)=\phi^{1}(x),\\ \psi(x,0)=\psi^{0}(x),&\psi_{t}(x,0)=\psi^{1}(x),\\ z_{1}(0)=\phi^{1}(1),&y_{1}(0)=\psi^{1}(1).\end{array}\right.$$

$$\begin{array}[]{l}\mathbb{X}=H^{1}\times L^{2},\ \mathbb{X}_{1}=H^{2}\times H^{1},\\ (\phi,\phi_{t},\psi,\psi_{t},z_{1},y_{1})\in C^{0}(\mathbb{X}_{1}^{2}\times\mathbb{R}^{2}).\end{array}$$

$$T(\theta )=c_b\theta +T_{nl}({\Omega }_0)=c_b\theta +T_0.$$

$$\begin{array}[]{l}\tilde{u}_{1}(t)=-k_{p}\left(\phi_{t}(0,t)-\Omega_{0}\right)-k_{i}z_{2}(t),\\ \dot{z}_{2}(t)=\phi_{t}(0,t)-\Omega_{0},\end{array}$$

$$\left\{\begin{array}[]{l}\phi_{tt}(x,t)=c^{2}\phi_{xx}(x,t)-\gamma_{t}\phi_{t}(x,t),\\ \phi_{x}(0,t)=(\tilde{g}+k_{p})\phi_{t}(0,t)-k_{p}\Omega_{0}+k_{i}C_{2}Z(t),\\ \phi_{t}(1,t)=C_{1}Z(t),\\ \dot{Z}(t)=AZ(t)+B\left[\begin{smallmatrix}\phi_{t}(0,t)\\ \phi_{x}(1,t)\end{smallmatrix}\right]+B_{2}\left[\begin{smallmatrix}\Omega_{0}\\ T_{0}\end{smallmatrix}\right],\end{array}\right.$$

$$\begin{array}[]{lll}A\!=\!\left[\begin{matrix}-\frac{c_{b}}{I_{B}}&0\\ 0&0\end{matrix}\right],&B_{\hphantom{1}}\!=\!\left[\begin{matrix}0&-\alpha_{1}\\ 1&0\end{matrix}\right],&B_{2}\!=\!\left[\begin{matrix}0&-\alpha_{2}\\ -1&0\end{matrix}\right],\\ Z\!=\!\left[\begin{matrix}z_{1}&z_{2}\end{matrix}\right]^{\top}\!\!\!\!,&C_{1}=\left[\begin{matrix}1&0\end{matrix}\right],&C_{2}=\left[\begin{matrix}0&1\end{matrix}\right].\end{array}$$

$$X=({\varphi }_x,{\varphi }_t,z_1,z_2)\in C^1([0,\infty ),\mathcal{H})$$

$$\Vert X{\Vert }_{\mathcal{H}}^2=z_1^2+z_2^2+c^2\Vert {\varphi }_x(\cdot ){\Vert }^2+\Vert {\varphi }_t(\cdot ){\Vert }^2.$$

$$\forall t>0,\Vert \dot{X}(t){\Vert }_{\mathcal{H}}=0.$$

$$\forall t>0,\Vert X(t)-X_{\infty }{\Vert }_{\mathcal{H}}\le \gamma \Vert X_0-X_{\infty }{\Vert }_{\mathcal{H}}{e}^{-\mu t}$$

$$X_{\infty }=\left({\varphi }_x^{\infty },{\Omega }_0,{\Omega }_0,\frac{{\varphi }_x^{\infty }(0)-\tilde{g}{\Omega }_0}{k_i}\right),$$

$$\begin{array}[]{l}\Theta_{N}=c\hskip 2.84544pt\Theta_{1,N}+\Theta_{2,N}-Q_{N}\prec 0,\\ P_{N}+S_{N}\succ 0,\end{array}$$

$$\begin{array}[]{l}\Gamma_{0}=cR+\frac{\gamma_{t}}{2}U_{0}-Q\succeq 0,\\ \Gamma_{1}=cR+\frac{\gamma_{t}}{2}U_{1}-Q\succeq 0,\end{array}$$

$$\begin{array}[]{cl}\!\!\!\!\!\!\Theta_{1,N}\!\!\!\!&=H_{N}^{\top}\left[\begin{smallmatrix}S_{1}+R_{1}&0\\ 0&-S_{2}\end{smallmatrix}\right]H_{N}-G_{N}^{\top}\left[\begin{smallmatrix}S_{1}&0\\ 0&-S_{2}-R_{2}\end{smallmatrix}\right]G_{N},\\ \!\!\!\!\!\!\Theta_{2,N}\!\!\!\!&=\operatorname{He}\left(D_{N}^{\top}P_{N}F_{N}\right),\\ \end{array}$$

$$\begin{array}[]{cl}\!\!\!\!\!\!F_{N}\!\!\!\!&=\left[\begin{matrix}I_{2+2(N+1)}&0_{2+2(N+1),2}\end{matrix}\right],\\ \!\!\!\!\!\!D_{N}\!\!\!\!&=\left[\begin{matrix}J_{N}^{\top}&cM_{N}^{\top}\end{matrix}\right]^{\!\top}\!\!\!,\quad J_{N}=\left[\begin{matrix}A&0_{2,2(N+1)}&B\end{matrix}\right],\\ \!\!\!\!\!\!M_{N}\!\!\!\!&=\mathbb{1}_{N}H_{N}\!-\!\bar{\mathbb{1}}_{N}G_{N}\!-\!\left[\begin{matrix}0_{2(N+1),2}&\!\!\!\!L_{N}&\!\!\!\!0_{2(N+1),2}\end{matrix}\right],\\ \!\!\!\!\!\!U_{0}\!\!\!\!&=\left[\begin{smallmatrix}2S_{1}&S_{1}+S_{2}+R_{2}\\ S_{1}+S_{2}+R_{2}&2(S_{2}+R_{2})\end{smallmatrix}\right],\ U_{1}=\left[\begin{smallmatrix}2(S_{1}+R_{1})&S_{1}+S_{2}+R_{1}\\ S_{1}+S_{2}+R_{1}&2S_{2}\end{smallmatrix}\right],\\ \end{array}$$

$$\begin{array}[]{cl}\!\!\!\!\!\!G_{N}\!\!\!\!&=\left[\begin{matrix}\begin{smallmatrix}ck_{i}C_{2}\\ -ck_{i}C_{2}\end{smallmatrix}&0_{2,2(N+1)}&G\end{matrix}\right],\quad G=\left[\begin{smallmatrix}1+c(\tilde{g}+k_{p})&0\\ 1-c(\tilde{g}+k_{p})&0\end{smallmatrix}\right],\\ \!\!\!\!\!\!H_{N}\!\!\!\!&=\left[\begin{matrix}\begin{smallmatrix}C_{1}\\ C_{1}\end{smallmatrix}&0_{2,2(N+1)}&H\end{matrix}\right],\quad\ \ H=\left[\begin{smallmatrix}0&c\\ 0&-c\end{smallmatrix}\right],\\ \!\!\!\!\!\!Q_{N}\!\!\!\!&=\mathrm{diag}(0_{2},Q,3Q,\cdots,(2N+1)Q,0_{2}),\\ \!\!\!\!\!\!S_{N}\!\!\!\!&=\mathrm{diag}\left(0_{2},S,3S,\cdots,(2N+1)S\right),\\ \!\!\!\!\!\!L_{N}\!\!\!\!&=\displaystyle\left[\ell_{j,k}\Lambda\right]_{j,k\in[0,N]}-\frac{\gamma_{t}}{2}\mathrm{diag}\left(\left[\begin{smallmatrix}1&1\\ 1&1\end{smallmatrix}\right],\dots,\left[\begin{smallmatrix}1&1\\ 1&1\end{smallmatrix}\right]\right),\\ \!\!\!\!\!\!\mathbb{1}_{N}\!\!\!\!&=\left[\begin{smallmatrix}\Lambda\\ \vdots\\ \Lambda\end{smallmatrix}\right],\quad\bar{\mathbb{1}}_{N}=\left[\begin{smallmatrix}\Lambda\\ \vdots\\ (-1)^{N}\Lambda\end{smallmatrix}\right],\quad\Lambda=\left[\begin{matrix}c&0\\ 0&-c\end{matrix}\right],\end{array}$$

$$\ell_{k,j}=\left\{\begin{array}[]{ll}(2j+1)(1-(-1)^{j+k}),&\text{ if }j\leqslant k,\\ 0&\text{ otherwise},\end{array}\right.$$

$$\tilde{X}(t)=X(t)-X_{\infty }=\left({\tilde{\varphi }}_x(t),{\tilde{\varphi }}_t(t),{\tilde{z}}_1(t),{\tilde{z}}_2(t)\right),$$

$$\left\{\begin{array}[]{l}\varepsilon_{1}\|\tilde{X}\|^{2}_{\mathcal{H}}\leqslant V(\tilde{X})\leqslant\varepsilon_{2}\|\tilde{X}\|^{2}_{\mathcal{H}},\\ \dot{V}(\tilde{X})+2\mu V(\tilde{X})\leqslant-\varepsilon_{3}\|\tilde{X}\|^{2}_{\mathcal{H}},\end{array}\right.$$

$$\begin{array}[]{ll}\tilde{\chi}(x)=\left[\begin{matrix}\tilde{\phi}_{t}(x)+c\tilde{\phi}_{x}(x)\\ \tilde{\phi}_{t}(x)-c\tilde{\phi}_{x}(x)\end{matrix}\right],&\displaystyle\tilde{\mathfrak{X}}_{k}(t)=\int_{0}^{1}\tilde{\chi}(x,t)\mathcal{L}_{k}(x)dx,\end{array}$$

$$V_N(\tilde{X})={\tilde{Z}}_N^{\top }{P}_N{\tilde{Z}}_N+\mathcal{V}(\tilde{\chi })$$

$$\mathcal{V}(\tilde{\chi })={\int }_0^1{\tilde{\chi }}^{\top }(x)\left[\begin{array}{cc}{S}_1+x{R}_1& 0\\ 0& S_2+(1-x)R_2\end{array}\right]\tilde{\chi }(x)dx.$$

Full text

I Introduction

Studying the behavior of complex machineries is a real challenge since they usually present nonlinear and coupled behaviors [49]. A drilling mechanism is a very good example of this. Many nonlinear effects can occur on the drilling pipe such as bit bouncing, stick-slip or whirling [14]. These phenomena induce generally some vibrations, increasing the drillpipe fatigue and affecting therefore the life expectancy of the well. The first challenge is then to provide a dynamical model which reflects such behaviors.

Looking at the literature, many models exist from the simplest finite-dimensional ones presented in for instance [13, 41] to the more complex but more realistic infinite-dimensional systems. Finite-dimensional systems were an important first step since they showed which characteristics are responsible of vibrations in the well. Nevertheless, they are too far from the physical laws which are expressed in terms of partial differential equation. Then, a coupled finite/infinite dimensional model seems more natural in the context of a drilling pipe and it was proposed in [14, 20, 38] for instance.

The second challenge was then to design a controller to remove, or at least weaken, these undesirable effects. Many control techniques were applied on the finite-dimensional model from the simple PI controller investigated in [13, 15], to more advance controllers as sliding mode control [30] or $H_{\infty}$ [41]. Nevertheless, extending these controllers on the coupled finite/infinite dimensional system is not straightforward.

The last decade has seen many developments regarding the analysis of infinite-dimensional systems. The semigroup theory, investigated in [48] for instance, was a great tool to simplify the proof of existence and uniqueness of a solution to this kind of problems. That leads to an extension of the Lyapunov theory to some classes of Partial Differential Equations (PDE) [9, 16, 33]. These advances have given rise to the stability analysis of the linearized infinite-dimensional drilling pipe with a PI controller [44, 45]. Since this kind of controller provides only two degree of freedom, to better enhance the performances, slightly different controllers arose. One of the most famous is the modified PI controller in [47], but there is also a delayed PI or a flatness based control in [38]. More complex controllers, coming from the backstepping technique for PDE, originally developed in [25], were also applied in [10, 12, 35] for instance.

Nevertheless, these techniques almost always use a Lyapunov argument to conclude and they consequently suffer from the lack of an efficient Lyapunov functional for coupled systems. Recent advances in the domain of time-delay systems [42] have lead to a hierarchy of Lyapunov functionals which are very efficient for coupled Ordinary Differential Equation (ODE) / string equation [8]. Since it relies on a state extension, the stability analysis cannot be assessed manually but it translates into an optimization problem expressed using Linear Matrix Inequalities (LMIs) and consequently easily solvable.

This paper takes advantage of this enriched Lyapunov functional to revisit the stability analysis of a PI controlled infinite-dimensional model of a drilling pipe for the torsion only. The first contribution of this article results in Theorem 1, which provides a Linear Matrix Inequality (LMI) to ensure the asymptotic stability of the linear closed-loop system. The second theorem deals with the practical stability of the controlled nonlinear plant. It shows for example that if the linear system is stable, the nonlinear system is also stable. Moreover, it provides an accurate bound on the oscillations during the stick-slip.

The article is organized as follows. Section 2 discusses the different models presented in the literature and enlighten the importance of treating the infinite-dimensional problem. Section 3 is the problem statement. Section 4 is dedicated to the study of the linear system. This is a first step before dealing with the nonlinear system, which is the purpose of Section 5. Section 6 finally proposes simulations to demonstrate the effectiveness of this approach and conclude about the design of a PI controller.

Notations: For a multi-variable function $(x,t)\mapsto u(x,t)$, the notation $u_{t}$ stands for $\frac{\partial u}{\partial t}$ and $u_{x}=\frac{\partial u}{\partial x}$. We also use the notations $L^{2}=L^{2}((0,1);\mathbb{R})$ and for the Sobolov spaces: $H^{n}=\{z\in L^{2};\forall m\leqslant n,\frac{\partial^{m}z}{\partial x^{m}}\in L^{2}\}$. The norm in $L^{2}$ is $\|z\|^{2}=\int_{\Omega}|z(x)|^{2}dx=\left<z,z\right>$. For any square matrices $A$ and $B$, the following operations are defined: $\text{He}(A)=A+A^{\top}$ and $\text{diag}(A,B)=\left[\begin{smallmatrix}A&0\\ 0&B\end{smallmatrix}\right]$. The set of positive definite matrices of size $n$ is denoted by $\mathbb{S}^{n}_{+}$ and, for simplicity, a matrix $P$ belongs to this set if $P\succ 0$.

II Model Description

A drilling pipe is a mechanism used to pump oil deep under the surface thanks to a drilling pipe as illustrated in Figure 1. Throughout the thesis, $\Phi(\cdot,t)$ is the twisting angle along the pipe and then $\Phi(0,t)$ and $\Phi(L,t)$ are the angles at the top and at the bottom of the well respectively. The well is a long metal rode of around one kilometer and consequently, the rotational velocity applied at the top using the torque $u_{1}(t)$ is different from the one at the bottom. Moreover, the interaction of the bit with the rock at the bottom are modeled by torque $T$, which depends on $\Phi_{t}(L,t)$.

As the bit drills the rock, an axial compression of the rod occurs and is denoted $\Psi$. This compression arises because of the propagation along the rode of the vertical force $u_{2}$ applied at the top to push up and down the well.

This description leads naturally to two control objectives to prevent the mechanism from breaking. The first one is to maintain the rotational speed at the end of the pipe $\Phi(L,t)$ at a constant value denoted here $\Omega_{0}$, preventing any twisting of the pipe. The other one is to keep the penetration rate constant such that there is no compression along the rode.

Several models have been proposed in the literature to achieve these control objectives. They are of very different natures and lead to a large variety of analysis and control techniques. The book [38] (chap. 2) and the survey [40] provide overviews of these techniques, which are, basically, of four kinds. To better motivate the model used in the sequel, a brief overview of the existing modeling tools is proposed but the reader can refer to [40] and the original papers to get a better understanding of how the models are constructed.

II-A Lumped Parameter Models (LPM)

These models are the first obtained in the literature [13, 27, 41] and the full mechanism is described by a sequence of harmonic oscillators. They can be classified into two main categories:

The first kind assumes that the dynamics of the twisting angles $\Phi(0,t)=\Phi_{r}(t)$ (at the top) and $\Phi(L,t)=\Phi_{b}(t)$ (at the bottom) are described by two coupled harmonic oscillators. The torque $u_{1}$ driving the system is applied on the dynamic of $\Phi_{r}$ and the controlled angle is $\Phi_{b}$. The axial dynamic is not taken into account in this model. This model can be found in [13, 32, 41] for instance. 2. 2.

The other two degrees of freedom model is described in [27, 34] for example. There also are two coupled harmonic oscillators for $\Psi(L,t)$ and $\Phi(L,t)$ representing the axial and torsional dynamics. This model only considers the motions at the end of the pipe and forget about the physics occurring along the rode.

The first class of models can be described by the following set of equations:

[TABLE]

where the parameters are given in Table I. $T$ is a torque modeled by a nonlinear function of $\dot{\Phi}_{b}$, it describes the bit-rock interaction111See [38] (chap. 3) for a detailed description about various models for $T$.. A second-order LPM can be derived by only taking into account the two dominant poles of the previous model.

An example of on-field measurements, depicted in Figure 2, shows the effect of this torque $T$ on the angular speed. The periodic scheme which arises is called stick-slip. It emerges because of the difference between the static and Coulomb friction coefficients making an antidamping on the torque function $T$. Even though the surface angular velocity seems not to vary much, there is a cycle for the downhole one and the angular speed is periodically close to zero, meaning that the bit is stuck to the rock.

The stick-slip effect appears mostly when dealing with a low desired angular velocity $\Omega_{0}$ on a controlled drilling mechanism. Indeed, if the angular speed $\Phi_{t}(L,t)$ is small, the torque provided by the rotary table at $x=0$ increases the torsion along the pipe. This increase leads to a higher $\Phi_{t}(L,t)$ but the negative damping on the torque function implies a smaller $T$. Consequently $\Phi_{t}(L,t)$ increases, this phenomenon is called the slipping phase. Then, the control law reduces the torque in order to match $\Phi_{t}(L,t)$ to $\Omega_{0}$. Since the torque increases as well, that leads to a sticking phase where $\Phi(L,t)$ remains close to [math]. A stick-slip cycle then emerges. Notice that this is not the case for high values of $\Omega_{0}$ since torque $T$ does not vary much with respect to $\Phi_{t}(L,t)$ making the system easier to control. In Figure 2, one can see that the frequency of the oscillations is $0.17$Hz and its amplitude is between $10$ and $25$ radians per seconds.

Modeling this phenomenon is of great importance as friction effects are quite common when studying mechanical machinery. In [40], the authors compare some models for $T$ and conclude that they produce very similar results. The main characteristic is a decrease of $T$ as $\Phi_{t}(L,t)$ increases. One standard model refers to the preliminary work of Karnopp [23] and Armstrong-Helouvry [3, 4] with an exponential decaying friction term as described in [31] for instance. This law is written thereafter where $\theta=\Phi_{t}(L,\cdot)$ is expressed in rad.s*-1*:

[TABLE]

This model has been used in [32, 38] for instance.

Notice that an on-field description of this mechanism applied in the particular context of drilling systems is provided in [2] and concludes that these models are fair approximations of the nonlinear phenomena visible in similar structures.

At this stage, a lumped parameter model is interesting for its simplicity but does not take into account the infinite-dimensional nature of the problem and, as a consequence, is a good approximation in the case of small vibrations only [39]. A deeper modeling can be done considering continuum mechanics and leading to a distributed parameter system.

II-B Distributed parameter models (DPM)

To tackle the finite-dimensional approximation of the previous model, another one derived from mechanical equations leads to a set of PDEs as described in the works [14, 49]. This model has been enriched in [1, 2] where the system is presented from a control viewpoint and compared to on-field measurements. In the first papers, the model focuses on the propagation of the torsion only along the pipe. The axial propagation was introduced in the model by [1, 40]. The new model is made up of two one-dimensional wave equations representing each deformation for $x\in(0,L)$ and $t>0$:

[TABLE]

where again $\Phi$ is the twist angle, $\Psi$ is the axial movement, $c_{t}=\sqrt{G/\rho}$ is the propagation speed of the angle, $\gamma_{t}$ is the internal damping, $c_{a}=\sqrt{E/\rho}$ is the axial velocity and $\gamma_{a}$ is the axial distributed damping. A list of physical parameters and their values is given in Table II and Figure 1 helps giving a better understanding of the physical system. In other words, if $\Psi(\cdot,t)=0$, then there is no compression in the pipe, meaning that the bit is not bouncing; if $\Phi_{tt}(\cdot,t)=0$, then the angular speed along the pipe is the same, meaning that there is no increase or decrease of the torsion.

For the previous model to be well-posed, top and bottom boundary conditions (at $x=0$ and $x=L$) must be incorporated in (3). In this part, only the topside boundary condition is derived. There is a viscous damping at $x=0$, and consequently a mismatch between the applied torque at the top and the angular speed. The topside boundary condition for the axial part is built on the same scheme and the following conditions are obtained for $t>0$:

[TABLE]

The downside boundary condition ($x=L$) is more difficult to grasp and is consequently derived later when dealing with a more complex model.

II-C Neutral-type time-delay model

Studying an infinite-dimensional problem stated in terms of PDEs represents a relevant challenge. The equations obtained previously are damped wave equations, but, for the special case where $\gamma_{a}=\gamma_{t}=0$, the system can be converted into a neutral time-delay system as done in [38]. This new formulation enables to use other tools to analyze its stability as the Lyapunov-Krasovskii Theorem or a frequency domain approach making its stability analysis slightly easier.

Nevertheless, the main drawback of this formulation is the assumption that the damping occurs at the boundary and not all along the pipe. This useful simplification, even if it encountered in many articles [12, 38, 39], is known to change in a significant manner the behavior of the system [1]. Indeed, it appears that without internal damping, the wave equation rephrases easily as a system of transport equations. It is then directly possible to observe with a delay of $c^{-1}$ at the top of the pipe what happened at the bottom of the pipe. This makes the control easier.

II-D Coupled ODE/PDE model

To overcome the issue mentioned previously, a simpler model than the one derived in (3) is proposed in [39], where an harmonic oscillator is used to describe axial vibrations and the model results in a coupled ODE/PDE.

A second possibility, reported in [14, 38] for example is to propose a second order ODE as the bottom boundary condition ($x=L$) for $t>0$:

[TABLE]

where $T$ represents the torque applied on the drilling bit by the rocks, described in equation (2). Notice that equation (5a) is coming from the conservation of angular momentum where $GJ\Phi_{x}(L,t)$ is the torque coming from the top of the pipe. Equation (5b) is the direct application of Newton’s second law of motion where $E\Gamma\Psi_{x}(L,t)$ is the force transmitted from the top to the bit and $\delta T\left(\Phi_{t}(L,t)\right)$ is the weight on bit due to the rock interaction. Since (5) is a second order in time differential equation, note that (3) together with (5) indeed leads to a coupled ODE/PDE.

There exist other bottom boundary conditions leading to a more complex coupling between axial and torsional dynamics. They nevertheless introduce delays which requires to have a better knowledge of the drilling bit. To keep the content general, the boundary conditions (5) used throughout this paper is proposed accordingly with [14, 39, 45].

As a final remark, using some transformations based on (3a), (4a) and (5a), it is possible to derive a system for which backstepping controllers can be used [12, 37]. This is the main reason why this model is widely used today.

II-E Models comparison

We propose in this subsection to compare the coupled ODE/PDE model and the lumped parameter models for the torsion only. We consider here a linearization of the system for large $\Omega_{0}$ and consequently we neglect the stick-slip effect by setting $T=0$.

First, denote by $\mathcal{H}_{DPM}$ the transfer function from $u_{1}$ to $\Phi(L,\cdot)$ for the DPM and $\mathcal{H}_{LPM}$ from $u_{1}$ to $\Phi_{b}$ for the LPM. We also define by $\mathcal{H}_{LPM2}$ a truncation of $\mathcal{H}_{LPM}$ considering only the two dominant poles. The Bode diagrams of $\mathcal{H}_{DPM}$, $\mathcal{H}_{LPM}$ and $\mathcal{H}_{LPM2}$ are drawn in Figure 4.

Clearly, the LPMs catch the behavior of the DPM at steady states and low frequencies until the resonance, occurring around $\sqrt{k/I_{b}}$ rad.s*-1*. From a control viewpoint, the DPM has infinitely many harmonics as it can be seen on the plots but of lower magnitudes and damped as the frequency increases (around $-10$ dB at each decade). The magnitude plots is not sufficient to make a huge difference between the three models. Nevertheless, considering the phase, we see a clear difference. It appears that the DPM crosses the frequency $-180$ many times making the control margins quite difficult to assess. Moreover, that shows that the DPM is harder to control because of the huge difference of behavior after the resonance. They may consequently has a very different behavior when controlled. That is why we focus in this study on the the DPM, even if it is far more challenging to control than the LPMs.

III Problem statement

The coupled nonlinear ODE/PDE model derived in the previous subsection can be written as follows for $t>0$ and $x\in[0,1]$:

[TABLE]

[TABLE]

where the normalized parameters are given in Table III. Note that the range for the spacial variable is now $x\in[0,1]$ to ease the calculations. The initial conditions are as follows:

[TABLE]

Since the existence and uniqueness of a solution to the previous problem is not the main contribution of this paper, it is assumed in the sequel. This problem has been widely studied (see [10, 12, 37, 39, 44] among many others) and the solution belongs to the following space if the initial conditions $(\phi^{0},\phi^{1},\psi^{0},\psi^{1})$ satisfies the boundary conditions (see [9] for more details):

[TABLE]

Remark 1

One may note that $\theta\mapsto T_{nl}(\theta)$ is not well defined for $\theta=0$ because of the sign function. Nevertheless, since the nonlinearity acts directly on the variable $z$, it follows that there exists a unique solution to the ODE system in the sense of Filipov [19]. A more detailed discussion on this point is provided in [11]. $\square$

System (6)-(7) is a cascade of two subsystems:

System (6) is a coupled nonlinear ODE/string equation describing the torsion angle $\phi$. 2. 2.

System (7) is a coupled linear ODE/string equation subject to the external perturbation $T(z_{1})$ for $\Psi$.

It appears clearly that the perturbation on the second subsystem depends on the first subsystem in $\phi$. Since they are very similar, the same analysis than the one conducted in this paper applies for the second subproblem. That is why we only study the evolution of the torsion.

IV Exponential Stability of the Linear System

System (6) is a nonlinear system because of the friction term $T_{nl}$ introduced in (2). Nevertheless, for a high desired angular speed $\Omega_{0}$, $T_{nl}$ can be assumed constant as seen in Figure 3. Moreover, studying this linear system can be seen as a first step before studying the nonlinear system, which relies mostly on the stability theorem derived in this section.

The proposed linear model of $T$ around $\Omega_{0}\gg 1$ is

[TABLE]

For high values of $\Omega_{0}$, $T_{0}=T_{nl}(\Omega_{0})$ is close to $T_{smooth}(\Omega_{0})$ and at the limit when $\Omega_{0}$ tends to infinity, they are equal. Hence the nonlinear friction term for relatively large angular velocity does not influence much the system.

In our case, the proposed controller for this problem has been studied in a different setting in [13, 15, 45] and is a simple proportional/integral controller based on the single measurement of the angular velocity at the top of the drill (i.e. $\phi_{t}(0,t)$). The following variables are therefore introduced:

[TABLE]

where $k_{p}$ and $k_{i}$ are the gains of the PI controller. Combining equations (6), (8) and (9) leads to:

[TABLE]

where

[TABLE]

To ease the notations, $c=\tilde{c}_{t}$. We denote by

[TABLE]

with $\mathcal{H}=L^{2}\times L^{2}\times\mathbb{R}^{2}$ the infinite dimensional state of system (10). The control objective in the linear case is to achieve the exponential stabilization of an equilibrium point of (10) in angular speed, i.e. $\phi_{t}(1)$ is going exponentially to a given constant reference value $\Omega_{0}$. To this extent, the following norm on $\mathcal{H}$ is introduced:

[TABLE]

The definition of an equilibrium point of (10) and its exponential stability follows from the previous definitions.

Definition 1

$X_{\infty}\in\mathcal{H}$* is an equilibrium point of (10) if for the trajectory $X\in C^{1}([0,\infty),\mathcal{H})$ of (10) with initial condition $X_{\infty}$, the following holds:*

[TABLE]

Definition 2

Let $X_{\infty}$ be an equilibrium point of (10). $X_{\infty}$ is said to be $\mu$-exponentially stable if

[TABLE]

holds for $\gamma\geqslant 1$, $\mu>0$ and for any initial conditions $X_{0}\in\mathcal{H}$ satisfying the boundary conditions. Here, $X$ is the trajectory of (10) whose initial condition is $X_{0}$.

Hence, $X_{\infty}$ is said to be exponentially stable if there exists $\mu>0$ such that $X_{\infty}$ is $\mu$-exponentially stable.

Before stating the main result of this part, a lemma about the equilibrium point is proposed.

Lemma 1

Assume $k_{i}\neq 0$, then there exists a unique equilibrium point $X_{\infty}=(\phi^{\infty}_{x},\phi^{\infty}_{t},z_{1}^{\infty},z_{2}^{\infty})\in\mathcal{H}$ of (10) and it satisfies $\phi^{\infty}_{t}=\Omega_{0}$.

Proof:

An equilibrium point $X_{\infty}$ of (10) is such that: $(\phi^{\infty}_{xt},\phi^{\infty}_{tt},\frac{d}{dt}z^{\infty}_{1},\frac{d}{dt}z^{\infty}_{2})=0$.

Since $\frac{d}{dt}z^{\infty}_{2}=\phi_{t}^{\infty}(0,t)-\Omega_{0}=0$ and $\partial_{x}\phi_{t}^{\infty}=0$ and $\partial_{t}\phi_{t}^{\infty}=0$, $\phi^{\infty}_{t}=\Omega_{0}$ holds.

We also get from $\phi_{xx}^{\infty}=0$ and $\partial_{t}\phi_{x}^{\infty}=0$ that $\phi^{\infty}_{x}$ is a first order polynomial in $x$. Together with the boundary conditions, this system has a unique solution if $k_{i}\neq 0$ such as:

[TABLE]

where $\phi_{x}^{\infty}(x)=\frac{\gamma_{t}\Omega_{0}}{c^{2}}(x-1)-\frac{c_{b}}{\alpha_{1}I_{B}}\Omega_{0}-\frac{\alpha_{2}}{\alpha_{1}}T_{0}$ for $x\in[0,1]$. ∎

IV-A Exponential stability of the closed-loop system (10)

The main result of this part is then stated as follows.

Theorem 1

Let $N\in\mathbb{N}$. Assume there exists $P_{N}\in\mathbb{S}^{2+2(N+1)}$, $R=\mathrm{diag}(R_{1},R_{2})\succeq 0,S=\mathrm{diag}(S_{1},S_{2})\succ 0$, $Q\in\mathbb{S}^{2}_{+}$ such that the following LMIs hold:

[TABLE]

and

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

then the equilibrium point of system (10) is exponentially stable. As a consequence, $\lim_{t\to+\infty}|\phi_{t}(1,t)-\Omega_{0}|=0$.

The proof is given thereafter but some practical consequences and preliminary results are firstly derived.

Remark 2

The methodology used to derive the previous result has been firstly introduced in [42] to deal with time-delay systems. It has been proved in the former article that this theorem provides a hierarchy of LMI conditions. That means if the conditions of Theorem 1 are met for $N=N_{1}\geqslant 0$, then the conditions are also satisfied for all $N>N_{1}$. Then, the LMIs provide a sharper analysis as the order $N$ of the theorem increases. Nevertheless, the price to pay for a more precise analysis is the increase of the number of decision variables and consequently a higher computational burden. It has been noted in [42] for a time-delay system and in [8] for a coupled ODE/string equation system that very sharp results are obtained even for low orders $N$. $\square$

As we aim at showing that $X_{\infty}$ is an exponentially stable equilibrium point of system (10) in the sense of Definition 2, the following variable is consequently defined for $t\geq 0$:

[TABLE]

where $X$ is a trajectory of system (10).

The following lemma, given in [7], provides a way for estimating the exponential decay-rate of system (10) as soon as a Lyapunov functional $V$ is known.

Lemma 2

Let $V$ be a Lyapunov functional for system (10) and $\mu\geq 0$. Assume there exist $\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}>0$ such that the following holds:

[TABLE]

then the equilibrium point of system (10) is $\mu$-exponentially stable. If $\mu=0$, then it is exponentially stable.

Proof:

For simplicity, the following notations are used throughout the paper for $k\in\mathbb{N}$:

[TABLE]

where $\left\{\mathcal{L}_{k}\right\}_{k\in\mathbb{N}}$ is the orthogonal family of Legendre polynomials as defined in Appendix -A. $\tilde{\mathfrak{X}}_{k}$ is then the projection coefficient of $\tilde{\chi}$ along the Legendre polynomial of degree $k$. $\tilde{\chi}$ refers to the Riemann coordinates and presents useful properties as discussed in [8] for instance.

IV-A1 Choice of a Lyapunov functional candidate

The proposed Lyapunov functional is inspired from [7, 8, 36] and is divided into two parts as follows:

[TABLE]

where $\tilde{Z}_{N}=\left[\tilde{z}_{1}\ \tilde{z}_{2}\ \tilde{\mathfrak{X}}_{0}^{\top}\ \cdots\ \tilde{\mathfrak{X}}_{N}^{\top}\right]^{\top}\!\!$ and

[TABLE]

$\mathcal{V}$ is a traditional Lyapunov functional candidate for a transport system as shown in [16].

IV-A2 Exponential stability

Existence of $\varepsilon_{1}$: This part is inspired by [7]. The inequalities $P_{N}+S_{N}\succ 0,R\succeq 0$ imply the existence of $\varepsilon_{1}>0$ such that:

[TABLE]

with $\mathcal{I}_{N}=\mathrm{diag}\left\{(2k+1)I_{2}\right\}_{k\in[0,N]}$. This statement implies the following on $V_{N}$:

[TABLE]

where $\tilde{S}=S-\frac{\varepsilon_{1}}{2}I_{2}$ which ends the proof of existence of $\varepsilon_{1}$.

Existence of $\varepsilon_{2}$: Following the same line as previously, the following holds for a sufficiently large $\varepsilon_{2}>0$:

[TABLE]

Using these inequalities in (16) leads to:

[TABLE]

The last inequality is a direct application of Lemma 4.

Existence of $\varepsilon_{3}$: The time derivation of $\tilde{\chi}$ leads to:

[TABLE]

Noting that $\tilde{\phi}_{t}(x)=\frac{1}{2}\left[\begin{smallmatrix}1&1\end{smallmatrix}\right]\tilde{\chi}(x)$, we get:

[TABLE]

The derivative of $\mathcal{V}$ along the trajectories of (6) leads to:

[TABLE]

with

[TABLE]

An integration by part on $\mathcal{V}_{1}$ shows that:

[TABLE]

The previous calculations lead to the following derivative:

[TABLE]

By a convexity argument, if (12) is verified then $cR+U(x)\succeq Q\succ 0$ holds for $x\in[0,1]$. Consequently, noticing that $\tilde{\chi}(0)=G_{N}\tilde{\xi}_{N}$, $\tilde{\chi}(1)=H_{N}\tilde{\xi}_{N}$, we get:

[TABLE]

with

[TABLE]

Using Lemma 5 and the fact that $\dot{\tilde{Z}}_{N}=D_{N}\tilde{\xi}_{N}$ and $\tilde{Z}_{N}=F_{N}\tilde{\xi}_{N}$, we get the following:

[TABLE]

with $\Theta_{N}$ defined in (11). Since $\Theta_{N}\!\prec\!0$ and $Q\!\succ\!0$, we get:

[TABLE]

Then, $\dot{V}_{N}$ is upper bounded by:

[TABLE]

The last inequality comes from a direct application of Bessel’s inequality (31).

Conclusion : Using Lemma 2, we indeed get the exponential convergence of all trajectories of (10) to the desired equilibrium point. ∎

IV-B Exponential stability with a guaranteed decay-rate

It is possible to estimate the decay-rate $\mu$ of the exponential convergence with a slight modification of the LMIs as it is proposed in the following corollary.

Corollary 1

Let $N\in\mathbb{N}$, $\mu>0$ and $\gamma_{t}\geq 0$. If there exist $P_{N}\in\mathbb{S}^{2+2(N+1)}$, $R=\mathrm{diag}(R_{1},R_{2})\succeq 0,S=\mathrm{diag}(S_{1},S_{2})\succ 0,Q\in\mathbb{S}^{2}_{+}$ such that (12) and the following LMIs hold:

[TABLE]

with the parameters defined as in Theorem 1, then the equilibrium point of system (10) is $\mu$-exponentially stable.

Proof:

To prove the exponential stability with a decay-rate of at least $\mu>0$, we use Lemma 2. Similarly to the previous proof, we have the existence of $\varepsilon_{1},\varepsilon_{2}>0$. The existence of $\varepsilon_{3}$ is slightly different. First, note that for the Lyapunov functional candidate (16), we get:

[TABLE]

This inequality was obtained using equation (31). Using the notations of the previous theorem yields:

[TABLE]

Coming back to (20) and using (21) leads to:

[TABLE]

Injecting inequality (22) into the previous inequality leads to:

[TABLE]

Using the same techniques than for the previous proof leads to the existence of $\varepsilon_{3}>0$ such that:

[TABLE]

Lemma 2 concludes then on the $\mu$-exponential stability. ∎

Remark 3

Wave equations can sometimes be modeled as neutral time-delay systems [6, 38]. This kind of system is known to possess some necessary stability conditions as noted in [6, 9]. In the book [9], the following criterion is derived:

[TABLE]

This result implies that there exists a maximum decay-rate and if this maximum is negative, then the system is unstable. The LMI $\Theta_{N}^{\mu}\prec 0$ contains the same necessary condition, meaning that the neutral aspect of the system is well-captured. $\square$

IV-C Strong stability against small delay in the control

A practical consequence of the neutral aspect of system (10) is that it is very sensitive to delay in the control (9). Indeed, if the control is slightly delayed, a new necessary stability condition (equivalent to (23)) coming from frequency analysis can be derived (Cor 3.3 from [21]):

[TABLE]

It is more restrictive and taking $k_{p}\neq 0$ practically leads to a decrease of the robustness (even if some other performances might be enhanced). This phenomenon has been studied in many articles [22, 29, 44]. Hence, to be robust to delays in the loop, one then needs to ensure the following:

[TABLE]

This inequality on $k_{p}$ comes when considering the infinite dimensional problem and does not arise when dealing with any finite dimensional model. That point enlightens that it is more realistic to consider the infinite dimensional problem. For more information on that point, the interested reader can refer to [5].

V Practical Stability of system (6)-(9)

The experiments conducted previously shows for $\Omega_{0}$ large, the trajectory of the nonlinear system (6)-(9) goes exponentially to $X_{\infty}$, so does the linear system. In other words, the previous result is a local stability test for the nonlinear system in case of large desired angular velocity $\Omega_{0}$ and does not extend straightforward to a global analysis.

In general, for a nonlinear system, ensuring the global exponential stability of an equilibrium point could be complicated. In many engineering situations, global exponential stability is not the requirement. Indeed, considering uncertainties in the system or nonlinearities, it is far more reasonable and acceptable for engineers to ensure that the trajectory remains close to the equilibrium point. This property is called dissipativeness in the sense of Levinson [26] or practical stability (also called globally uniformly ultimately bounded in [24]).

Definition 3

System (6) is practically stable if there exists $X_{bound}\geq 0$ such that for $X$ the solution of (6) with initial condition $X_{0}\in\mathcal{H}$:

[TABLE]

Saying that system (6)-(9) is practically stable means that there exists $T_{\eta}>0$ such that for any $x\in[0,1]$, $\phi_{t}(x,t)$ stays close to $\Omega_{0}$ for $t\geq T_{\eta}$. This property has already been applied to a drilling system in [39] for instance. In the nonlinear case, the aim is then to design a control law reducing the amplitude of the stick-slip when it occurs.

The objective of this section is to derive an LMI test ensuring the practical stability of nonlinear system (6) together with the control defined in (9):

[TABLE]

Remark 4

To simplify the writing, $\phi$ can both means the solution of the linear or nonlinear system depending on the context. In this part, it refers to a solution of nonlinear system (25). $\square$

V-A Practical stability of system (25)

The idea behind practical stability is to embed the static nonlinearities by the use of suitable sector conditions, as it has been done for the saturation for instance [43]. Then, the use of robust tools will lead to some LMI tests ensuring the practical stability.

Lemma 3

For almost all $\tilde{z}_{1}\in\mathbb{R}$, the following holds:

[TABLE]

Proof:

These inequalities can be easily verified using (2). ∎

These new information are the basis of the following theorem.

Theorem 2

Let $N\in\mathbb{N}$ and $V_{max}>0$. If there exist $P_{N}\in\mathbb{S}^{2+2(N+1)}$, $R=\mathrm{diag}(R_{1},R_{2})\succeq 0,S=\mathrm{diag}(S_{1},S_{2})\succ 0,Q\in\mathbb{S}^{2}_{+},\tau_{0},\tau_{1},\tau_{2},\tau_{3}\geq 0$ such that (12) holds together with:

[TABLE]

where

[TABLE]

and all the parameters defined as in Theorem 1, then the equilibrium point $X_{\infty}$ of system (25) is practically stable. More precisely, equation (24) holds for $X_{bound}=\sqrt{V_{max}\varepsilon_{1}^{-1}}$ where $\varepsilon_{1}$ is defined in (17).

Proof:

First, let do the same change of variable as before:

[TABLE]

Since the nonlinearity affects only the ODE part of system (25), the difference with the previous part lies in the dynamic of $\tilde{Z}$:

[TABLE]

Using the same Lyapunov functional as in (16), the positivity is ensured in the exact same way. For the bound on the time derivative, following the same strategy as before, we easily get:

[TABLE]

We introduce a new extended state variable $\bar{\xi}_{N}=\left[\begin{matrix}\tilde{\xi}_{N}^{\top}&T_{nl}(\tilde{z}_{1}+\Omega_{0})&1\end{matrix}\right]^{\top}$. Using the notation of Theorem 2, equation (28) rewrites as:

[TABLE]

It is impossible to ensure $\bar{\Theta}_{N}\prec 0$ since its last $2$ by $2$ diagonal block is $0_{2}$.

We then use the definition of practical stability. We want to show that if there exists $\varepsilon_{3}>0$ such that $\dot{V}_{N}(\tilde{X})\leq-\varepsilon_{3}\|\tilde{X}\|_{\mathcal{H}}^{2}$ when $V_{N}(\tilde{X})\geq V_{max}$ then the system is practically stable.

Let $\mathcal{S}=\left\{\tilde{X}\in\mathcal{H}\ |\ V_{N}(\tilde{X})\leq V_{max}\right\}$, the previous assertion implies that this set is invariant and attractive. Using (17), we get that $\left\{\tilde{X}\in\mathcal{H}\ |\ \|\tilde{X}\|_{\mathcal{H}}\leq X_{bound}=V_{max}^{1/2}\varepsilon_{1}^{-1/2}\right\}\supseteq\mathcal{S}$, meaning that the system is practically stable.

A sufficient condition to be practically stable is that $V_{N}$ should be strictly decreasing when outside the ball of radius $V_{max}$. This condition rewrites as $V_{N}(\tilde{X})\geq V_{max}$ and the following holds:

[TABLE]

The previous inequality is obtained using Bessel inequality (31) on $\int_{0}^{1}\tilde{\chi}^{\top}(x)S\tilde{\chi}(x)dx$. Hence, $V_{N}(\tilde{X})\geq V_{max}$ if $\bar{\xi}_{N}^{\top}\Pi_{0}\bar{\xi}_{N}\leq 0$.

Noting that $\tilde{z}_{1}=\pi_{1}\bar{\xi}_{N},T_{nl}(\tilde{z}_{1}+\Omega_{0})=\pi_{2}\bar{\xi}_{N}$ and $1=\pi_{3}\bar{\xi}_{N}$, Lemma 3 rewrites as:

[TABLE]

Consequently, a sufficient condition to be practically stable is:

[TABLE]

The technique called S-variable, explained in [18] for instance, translates the previous inequalities into an LMI condition. Indeed, Theorem 1.1 from [18] shows that condition (29) is verified if there exists $\tau_{0},\tau_{1},\tau_{2},\tau_{3}>0$ such that

[TABLE]

Consequently, in a similar way than for Theorem 1, condition (27) implies that $\mathcal{S}$ is an invariant and attractive set for system (25). Then, the equilibrium point $X_{\infty}$ of system (25) is practically stable with $X_{bound}=\sqrt{V_{max}\varepsilon_{1}^{-1}}$. ∎

Note that if the torque function is not perfectly known, one can change the lower and upper bounds $T_{min}$ and $T_{max}$ to get a more conservative result but robust to uncertainties on $T_{nl}$. For (27) to be feasible, the constraint $\bar{\xi}_{N}^{\top}\Pi_{1}\bar{\xi}_{N}\leq 0$ must hold, meaning that an upper bound of $T_{nl}$ needs to be proposed.

V-B On the optimization of $X_{bound}$

The condition (27) is a Bilinear Matrix Inequality (BMI) since $\tau_{0}$, $\tau_{1}$, $\tau_{2}$ and $\tau_{3}$ are decision variables and it is therefore difficult to get its global optimum. Nevertheless, the following lemma gives a sufficient condition for the existence of a solution to this problem.

Corollary 2

There exists $V_{max}$ and $\tau_{0}>0$ such that Theorem 2 holds if and only if there exists $N>0$ such that LMIs (12) and (11) are satisfied.

In other words, the equilibrium point of system (25) is practically stable if and only if the linear system (10) is exponentially stable.

Proof:

Note first that expending (27) with $\tau_{3}=0$ leads to:

[TABLE]

where $\kappa_{P_{N}}\in\mathbb{R}^{(4+2(N+1))\times 2}$ depends only on $P_{N}$.

Proof of sufficiency: Assume there exists $N>0$ such that LMIs (12) and (11) are satisfied. Considering $\tau_{2}=0$ and using Schur complement on $\Xi_{N}$, $\Xi_{N}\prec 0$ is equivalent to:

[TABLE]

with $\tau_{0}V_{max}>\tau_{1}T_{max}^{2}$. Since $\Theta_{N}\prec 0$, considering $\tau_{0}$ small enough, $\tau_{1}$ large and $V_{max}>\tau_{1}\tau_{0}^{-1}T_{max}^{2}$ the previous condition is always satisfied and Theorem 2 applies.

Proof of necessity: Assume $\Xi_{N}\prec 0$ and (12) holds. Then its first diagonal block must be definite negative. Consequently $\Theta_{N}\prec 0$ and, according to Theorem 1, system (10) is exponentially stable. ∎

Remark 5

Note that (30) provides a necessary condition for Theorem 2 which is $\Theta_{N,\tau_{0}/2}\prec 0$. In other words, $\tau_{0}$ is related to the decay-rate of the linear system and we get the following condition: $\tau_{0}<2\mu_{max}$. $\square$

Thanks to Corollary 2, the following method should help solving the BMI if Theorem 1 is verified. Assuming that equations (11) and (12) are verified for a given $N\in\mathbb{N}$.

Fix $\tau_{0}=2\mu_{max}$ as defined in (23). 2. 2.

Check that equations (11), (12) and (27) are satisfied for $V_{max}$ a strictly positive decision variable. If this is not the case, then decrease $\tau_{0}$ and do this step again. 3. 3.

Thanks to Corollary 2, there exists a $\tau_{0}$ small enough for which equations (11), (12) and (27) hold. Freeze this value. 4. 4.

Since the problem is unbounded, it is possible to fix a variable without loss of generality, let $V_{max}=10^{4}$ for instance and solve the following optimization problem:

[TABLE] 5. 5.

Then compute $X_{bound}=\sqrt{V_{max}\varepsilon_{1}^{-1}}$ where $\varepsilon_{1}$ is defined in (17).

VI Examples and Simulations

The section is devoted to numerical simulations222Numerical simulations are done using a first order approximation with at least $80$ space-discretization points and $9949$ time-discretization points. Simulations are done using Yalmip [28] together with the SDP solver SDPT-3 [46]. The code is available at https://homepages.laas.fr/mbarreau. and draw some conclusions about the PI regulation. In the first subsection, we focus on the linear system and the second one is dedicated to the nonlinear case.

VI-A On the linear model

This subsection recaps the result on the linear model.

VI-A1 Estimation of the decay-rate

The main result is a direct application of Corollary 1 for $k_{p}=10^{-3}$ and $k_{i}=10$. Indeed, using a dichotomy-kind algorithm, Table IV is obtained. It shows the estimated decay-rate $\mu$ at a given order between [math] and $6$.

The first thing to note in Table IV is the hierarchy property, the decay-rate is an increasing function of the degree, as noted in Remark 2. Note also that the gap between order [math] and order $1$ is significant, showing that using projections indeed improves the results.

For orders higher than $2$, the estimated decay-rate increases slightly, and, up to a four digits precision, it reaches its maximum value at $N=6$. Since it is then the maximum allowable decay-rate obtained using equation (23), that tends to show that the Lyapunov functional used in this paper together with condition (21) are sharp and provide a good analysis. Figure 5 represents a simulation on the linear system and it confirms the same observation. Indeed, one can see that the energy of the system is well-bounded by an exponential curve and that the bound becomes more and more accurate as $N$ increases.

VI-A2 Stability of the closed loop system

We are interested now in estimating the stability area in terms of the gains $k_{p}$ and $k_{i}$ such that the decay-rate of the coupled system is $\mu_{max}$ for an order $N=5$. That leads to Figure 6 where it is easy to see that increasing the gain $k_{p}$ decreases the range of possible $k_{i}$ while it increases its speed (see equation (23)). It is quite natural to note that the larger $k_{i}$ is, the slower the system is, while increasing the proportional gain leads to a faster system. As a conclusion, for small values of $k_{p}$ and $k_{i}$, the system is stable and that was the conclusion of the two papers [44, 45] using a different Lyapunov functional. Note that with the previous papers, it was not possible to quantify the notion of “small enough gains $k_{p}$ and $k_{i}$” while it is possible to give an estimation with the method of this paper.

VI-B The stick-slip effect on the nonlinear model

The previous subsection shows that the linear model is globally asymptotically stable for some values of gains $k_{p}$ and $k_{i}$, so the nonlinear system (25) is locally asymptotically stable for a large desired angular velocity $\Omega_{0}$. This can been verified in Figure 7a for $k_{p}=10^{-3}$, $k_{i}=10$ and $\Omega_{0}=10$. We can see that the linear and the nonlinear systems behave similarly if their initial condition is close to the equilibrium.

The higher $\Omega_{0}$ is, the larger the basin of attraction is. Indeed, the regulation tries to bring the system into the “quasi”-linear zone of $T_{nl}$, where it is close to a constant $T_{min}$ as seen in Figure 3 and consequently, the stick-slip phenomenon may occur at the beginning but it is not effective for a long time as depicted in Figure 7b which is a numerical simulation on the nonlinear model with a zero initial condition.

The real challenge is then the case of low desired angular velocities $\Omega_{0}$. The result of a simulation on the nonlinear model for $\Omega_{0}=5$ rad.s*-1* ($k_{p}=10^{-3}$, $k_{i}=10$) is displayed in Figure 7c. First of all, note that the oscillations are with a frequency of $0.2$Hz and with an amplitude of roughly $15$ rad.s*-1*. This is very close to what has been estimated using Figure 2. Then, the model presented in Section II seems to be a valid approximation of the real behavior, at least concerning the stick-slip phenomenon.

VI-C Practical stability analysis

Now, one can evaluate the amplitude of the oscillations using Theorem 2 with $k_{p}=10^{-3}$ and $k_{i}=10$. The result for several values of $\tau_{0}$ and for an order between [math] and $7$ is depicted in Figure 8. Note that after order $8$, there are numerical errors in the optimization process and the result is not accurate. The maximum $\tau_{0}$ is $2\mu_{max}=0.0155$ and the higher $\tau_{0}$ is, the better is the optimization. It appears that for $k_{i}=10$, $k_{p}=10^{-3}$ and $\Omega_{0}=5$, the optimal $X_{bound}$ is around $28$.

Figure 10 shows the energy of the system as a function of time. One can see that the bound $X_{bound}$ is quite accurate since the error between the maximum of the auto-oscillations and $X_{bound}$ is around $53\%$. Moreover, note that $\max|z_{1}-\Omega_{0}|=11.7=0.4X_{bound}$, in other words, nearly half of the oscillations are concentrated in the variable $z_{1}$, which means the stick-slip mostly acts on the variable $z_{1}$ and does not affect much the rest of the system. Particularly, it seems very difficult to estimate the variation of $z_{1}$ knowing only $\phi_{t}(0,t)$.

The final observation is about the variation of $X_{bound}$ for different $\Omega_{0}$. This is depicted in Figure 9. Up to errors in the numerical optimization, it seems that there is not important variation of $X_{bound}$ when $\Omega_{0}$ increases. This is counter-intuitive and does not reflect the observations made with Figure 7. One explanation is that we didn’t state that $T_{nl}$ is a strictly decreasing function for positive $\theta$.

VI-D Design of a PI controller

Finally, the problem stated at the beginning of the paper was to find the best PI controller, meaning that it minimizes $X_{bound}$. The plot in Figure 11a shows that the value of the integral gain $k_{i}$ does impact the oscillations due to stick-slip since $X_{bound}$ increases from $25$ to $43$ for $k_{i}\in[0.5,16]$ and the LMIs become infeasible after this point (these values have been obtained with $\Omega_{0}=5$, $k_{p}=10^{-3}$ and $N=5$).

To stay robustly stable against small delay in the loop, as noted in Section IV-C, we should consider $0\leq k_{p}\leq 2.1\cdot 10^{-3}$ which does not offer a large set of choices. For $k_{i}=10$, we get Figure 11b. It appears that the minimum of $X_{bound}$ is obtained for $k_{p}$ close to [math]. Consequently, increasing the gain $k_{p}$ seems not to reduce the stick-slip effect. Consequently, even if a PI controller does not weaken the stick-slip effect, the equilibrium point of the controlled system is practically stable. Moreover, it does enable an oscillation around the desired equilibrium point and a local convergence to that point.

VII Conclusion

This paper focused on the analysis of the performance of a PI controller for a drilling pipe. First, a discussion between existing models in the literature allows us to conclude that the infinite-dimensional model is closer to the real drilling pipe and should then be used for simulations and analysis. Based on this last model, the exponential stability of the closed-loop system was ensured using a Lyapunov functional depending on projections of the infinite-dimensional state on a polynomials basis. This approach enables getting exponential stability of the linear system with an estimation of the decay-rate. This result was then extended to the nonlinear case considering practical stability. The example section shows that it is possible to get an estimate of the smallest attractive and invariant set and that this estimate is close to the minimal one. Further research would focus on similar analysis but using different controllers as the ones developed in [13, 30, 41] or in Chapter 10 of [38].

-A Legendre Polynomials and Bessel Inequality

The performance of the Lyapunov functional (16) highly depends on the projection methodology developed in [42]. To help the reader to better understand the proof of Theorem 1, a definition and some properties of Legendre polynomials is reminded. More information can be found in [17].

Definition 4

The orthonormal family of Legendre polynomials $\{\mathcal{L}_{k}\}_{k\in\mathbb{N}}$ on $L^{2}([0,1])$ embedded with the canonical inner product is defined as follows:

[TABLE]

where $\left(\begin{smallmatrix}k\\ l\end{smallmatrix}\right)=\frac{k!}{l!(k-l)!}$.

Their expended formulation is not useful in this paper but the two following lemmas are of main interest.

Lemma 4

For any function $\chi\in L^{2}$ and symmetric positive matrix $R\in\mathbb{S}^{2}_{+}$, the following Bessel-like integral inequality holds for all $N\in\mathbb{N}$:

[TABLE]

where $\mathfrak{X}_{k}$ is the projection coefficient of $\chi$ with respect to the Legendre polynomial $\mathcal{L}_{k}$ as defined in (15).

Lemma 5

For any function $\chi\in L^{2}$ satisfying (15), the following expression holds for any $N$ in $\mathbb{N}$:

[TABLE]

where $L_{N},\mathbb{1}_{N}$ and $\bar{\mathbb{1}}_{N}$ are defined in (13).

Proof:

This proof is highly inspired from [8]. Since $\chi\in L^{2}$ satisfies (15), equation (18) can be derived and the following holds:

[TABLE]

As noted in [17], the interesting properties of Legendre polynomials are stated below:

the boundary conditions on $\mathcal{L}_{k}$ ensures $\mathcal{L}_{k}(0)=(-1)^{k}$ and $\mathcal{L}_{k}(1)=1$; 2. 2.

the derivation rule for Legendre polynomials is $\frac{d}{dx}\mathcal{L}_{k}(x)=\sum_{j=0}^{k}\ell_{j,k}\mathcal{L}_{j}(x)$.

These two properties lead to the proposed result in equation (32). ∎