Practical Stability Analysis of a Drilling Pipe under Friction with a PI-Controller
PPractical Stability Analysis of a Drilling Pipe under Friction with aPI-Controller
Matthieu Barreau, Fr´ed´eric Gouaisbaut and Alexandre Seuret
Abstract —This paper deals with the stability analysis of adrilling pipe controlled by a PI controller. The model is a coupledODE / PDE and is consequently of infinite dimension. Usingrecent advances in time-delay systems, we derive a new Lyapunovfunctional based on a state extension made up of projections ofthe Riemann coordinates. First, we will provide an exponentialstability result expressed using the LMI framework. This resultis dedicated to a linear version of the torsional dynamic. On asecond hand, the influence of the nonlinear friction force, whichmay generate the well-known stick-slip phenomenon, is analyzedthrough a new stability theorem. Numerical simulations showthe effectiveness of the method and that the stick-slip oscillationscannot be weaken using a PI controller.
I. I
NTRODUCTION
Studying the behavior of complex machineries is a realchallenge since they usually present nonlinear and coupledbehaviors [49]. A drilling mechanism is a very good exampleof this. Many nonlinear effects can occur on the drillingpipe such as bit bouncing, stick-slip or whirling [14]. Thesephenomena induce generally some vibrations, increasing thedrillpipe fatigue and affecting therefore the life expectancy ofthe well. The first challenge is then to provide a dynamicalmodel which reflects such behaviors.Looking at the literature, many models exist from thesimplest finite-dimensional ones presented in for instance[13], [41] to the more complex but more realistic infinite-dimensional systems. Finite-dimensional systems were an im-portant first step since they showed which characteristics areresponsible of vibrations in the well. Nevertheless, they aretoo far from the physical laws which are expressed in termsof partial differential equation. Then, a coupled finite/infinitedimensional model seems more natural in the context of adrilling pipe and it was proposed in [14], [20], [38] forinstance.The second challenge was then to design a controller toremove, or at least weaken, these undesirable effects. Manycontrol techniques were applied on the finite-dimensionalmodel from the simple PI controller investigated in [13], [15],to more advance controllers as sliding mode control [30] or H ∞ [41]. Nevertheless, extending these controllers on the cou-pled finite/infinite dimensional system is not straightforward.The last decade has seen many developments regardingthe analysis of infinite-dimensional systems. The semigrouptheory, investigated in [48] for instance, was a great tool tosimplify the proof of existence and uniqueness of a solutionto this kind of problems. That leads to an extension of the M. Barreau, A. Seuret, F. Gouaisbaut are with LAAS - CNRS, Universit´e deToulouse, CNRS, UPS, France e-mail: (mbarreau,aseuret,[email protected]).This work is supported by the ANR project SCIDiS contract number 15-CE23-0014.
Lyapunov theory to some classes of Partial Differential Equa-tions (PDE) [9], [16], [33]. These advances have given riseto the stability analysis of the linearized infinite-dimensionaldrilling pipe with a PI controller [44], [45]. Since this kindof controller provides only two degree of freedom, to betterenhance 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 backsteppingtechnique for PDE, originally developed in [25], were alsoapplied in [10], [12], [35] for instance.Nevertheless, these techniques almost always use a Lya-punov argument to conclude and they consequently sufferfrom the lack of an efficient Lyapunov functional for coupledsystems. Recent advances in the domain of time-delay systems[42] have lead to a hierarchy of Lyapunov functionals whichare 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 ittranslates into an optimization problem expressed using LinearMatrix Inequalities (LMIs) and consequently easily solvable.This paper takes advantage of this enriched Lyapunovfunctional to revisit the stability analysis of a PI controlledinfinite-dimensional model of a drilling pipe for the torsiononly. The first contribution of this article results in Theorem 1,which provides a Linear Matrix Inequality (LMI) to ensurethe asymptotic stability of the linear closed-loop system.The second theorem deals with the practical stability of thecontrolled nonlinear plant. It shows for example that if thelinear system is stable, the nonlinear system is also stable.Moreover, it provides an accurate bound on the oscillationsduring the stick-slip.The article is organized as follows. Section 2 discussesthe different models presented in the literature and enlightenthe importance of treating the infinite-dimensional problem.Section 3 is the problem statement. Section 4 is dedicated tothe study of the linear system. This is a first step before dealingwith the nonlinear system, which is the purpose of Section 5.Section 6 finally proposes simulations to demonstrate theeffectiveness of this approach and conclude about the designof a PI controller.
Notations:
For a multi-variable function ( x, t ) (cid:55)→ u ( x, t ) ,the notation u t stands for ∂u∂t and u x = ∂u∂x . We also usethe notations L = L ((0 , (cid:82) ) and for the Sobolov spaces: H n = { z ∈ L ; ∀ m (cid:54) n, ∂ m z∂x m ∈ L } . The norm in L is (cid:107) z (cid:107) = (cid:82) Ω | z ( x ) | dx = (cid:104) z, z (cid:105) . For any square matrices A and B , the following operations are defined: He ( A ) = A + A (cid:62) and diag ( A, B ) = [ A B ] . The set of positive definite matricesof size n is denoted by (cid:83) n + and, for simplicity, a matrix P belongs to this set if P (cid:31) . a r X i v : . [ m a t h . O C ] M a r otarytable KellyWinchDrill pipesDrill collarsBit Fig. 1: Schematic of a drilling mechanism originally takenfrom [39]. Data corresponding to physical values are given inTable II. II. M
ODEL D ESCRIPTION
A drilling pipe is a mechanism used to pump oil deep underthe surface thanks to a drilling pipe as illustrated in Figure 1.Throughout the thesis, Φ( · , t ) is the twisting angle along thepipe and then Φ(0 , t ) and Φ( L, t ) are the angles at the top andat the bottom of the well respectively. The well is a long metalrode of around one kilometer and consequently, the rotationalvelocity applied at the top using the torque u ( t ) is differentfrom the one at the bottom. Moreover, the interaction of the bitwith the rock at the bottom are modeled by torque T , whichdepends on Φ t ( L, t ) .As the bit drills the rock, an axial compression of the rodoccurs and is denoted Ψ . This compression arises because ofthe propagation along the rode of the vertical force u appliedat the top to push up and down the well.This description leads naturally to two control objectivesto prevent the mechanism from breaking. The first one is tomaintain the rotational speed at the end of the pipe Φ( L, t ) ata constant value denoted here Ω , preventing any twisting ofthe pipe. The other one is to keep the penetration rate constantsuch that there is no compression along the rode.Several models have been proposed in the literature toachieve these control objectives. They are of very differentnatures and lead to a large variety of analysis and controltechniques. The book [38] (chap. 2) and the survey [40]provide overviews of these techniques, which are, basically, offour kinds. To better motivate the model used in the sequel, abrief overview of the existing modeling tools is proposed butthe reader can refer to [40] and the original papers to get abetter understanding of how the models are constructed. 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 sequenceof harmonic oscillators. They can be classified into two maincategories:1) The first kind assumes that the dynamics of the twistingangles
Φ(0 , t ) = Φ r ( t ) (at the top) and Φ( L, t ) = Φ b ( t ) (at the bottom) are described by two coupled harmonicoscillators. The torque u driving the system is appliedon the dynamic of Φ r and the controlled angle is Φ b . Theaxial dynamic is not taken into account in this model.This model can be found in [13], [32], [41] for instance.2) The other two degrees of freedom model is describedin [27], [34] for example. There also are two coupledharmonic oscillators for Ψ( L, t ) and Φ( L, t ) represent-ing the axial and torsional dynamics. This model onlyconsiders the motions at the end of the pipe and forgetabout the physics occurring along the rode.The first class of models can be described by the followingset of equations: I r ¨Φ r + λ r ( ˙Φ r − ˙Φ b ) + k (Φ r − Φ b ) + d r ˙Φ r = u ,I b ¨Φ b + λ b ( ˙Φ b − ˙Φ r ) + k (Φ b − Φ r ) + d b ˙Φ b = − T ( ˙Φ b ) , (1)where the parameters are given in Table I. T is a torquemodeled by a nonlinear function of ˙Φ b , it describes the bit-rock interaction . A second-order LPM can be derived by onlytaking into account the two dominant poles of the previousmodel.An example of on-field measurements, depicted in Figure 2,shows the effect of this torque T on the angular speed. Theperiodic scheme which arises is called stick-slip . It emergesbecause of the difference between the static and Coulombfriction coefficients making an antidamping on the torquefunction T . Even though the surface angular velocity seemsnot to vary much, there is a cycle for the downhole one andthe angular speed is periodically close to zero, meaning thatthe bit is stuck to the rock.The stick-slip effect appears mostly when dealing witha low desired angular velocity Ω on a controlled drillingmechanism. Indeed, if the angular speed Φ t ( L, t ) is small,the torque provided by the rotary table at x = 0 increases thetorsion along the pipe. This increase leads to a higher Φ t ( L, t ) but the negative damping on the torque function implies asmaller T . Consequently Φ t ( L, t ) increases, this phenomenonis called the slipping phase. Then, the control law reducesthe torque in order to match Φ t ( L, t ) to Ω . Since the torqueincreases as well, that leads to a sticking phase where Φ( L, t ) remains close to . A stick-slip cycle then emerges. Notice thatthis is not the case for high values of Ω since torque T doesnot vary much with respect to Φ t ( L, t ) making the systemeasier to control. In Figure 2, one can see that the frequencyof the oscillations is . Hz and its amplitude is between and radians per seconds. See [38] (chap. 3) for a detailed description about various models for T .
10 20 30 40 50-5051015202530
Fig. 2: Nonlinear effect on the drilling mechanism due to thefriction torque at the bottom of the pipe. These measurementsare done on the field [41]. -5 0 5-1.5-1-0.5 Fig. 3: Nonlinear part of the torque. T nl is an approximationof the real friction torque coming from Karnopp’s work [23].Modeling this phenomenon is of great importance as frictioneffects are quite common when studying mechanical machin-ery. In [40], the authors compare some models for T andconclude that they produce very similar results. The maincharacteristic is a decrease of T as Φ t ( L, t ) increases. Onestandard model refers to the preliminary work of Karnopp [23]and Armstrong-Helouvry [3], [4] with an exponential decayingfriction term as described in [31] for instance. This law iswritten thereafter where θ = Φ t ( L, · ) is expressed in rad.s − : T ( θ ) = T l ( θ ) + T nl ( θ ) ,T l ( θ ) = c b θ,T nl ( θ ) = T sb (cid:0) µ cb + ( µ sb − µ cb ) e − γ b | θ | (cid:1) sign( θ ) . (2)This model has been used in [32], [38] for instance.Notice that an on-field description of this mechanism ap-plied in the particular context of drilling systems is providedin [2] and concludes that these models are fair approximationsof the nonlinear phenomena visible in similar structures.At this stage, a lumped parameter model is interesting forits simplicity but does not take into account the infinite-dimensional nature of the problem and, as a consequence, isa good approximation in the case of small vibrations only[39]. A deeper modeling can be done considering continuummechanics and leading to a distributed parameter system. Parameter meaning Value I r Rotary table and drive inertia kg.m I b Bit and drillstring inertia kg.m k Drillstring stiffness
N.m.rad − λ r Coupled damping at top
N.m.s.rad − λ b Coupled damping at bottom . N.m.s.rad − d r Rotary table damping
N.m.s.rad − d b Bit damping N.m.s.rad − γ b Velocity decrease rate . s.rad − µ cb Coulomb friction coefficient . µ sb Static friction coefficient . c b Bottom damping constant . N.m.s.rad − T sb Static / Friction torque
15 145
N.mTABLE I: Parameters values and their meanings for thelumped parameter model taken from [13], [32], [41].
B. Distributed parameter models (DPM)
To tackle the finite-dimensional approximation of the pre-vious model, another one derived from mechanical equationsleads to a set of PDEs as described in the works [14], [49].This model has been enriched in [1], [2] where the systemis presented from a control viewpoint and compared to on-field measurements. In the first papers, the model focuses onthe propagation of the torsion only along the pipe. The axialpropagation was introduced in the model by [1], [40]. Thenew model is made up of two one-dimensional wave equationsrepresenting each deformation for x ∈ (0 , L ) and t > : Φ tt ( x, t ) = c t Φ xx ( x, t ) − γ t Φ t ( x, t ) , (3a) Ψ tt ( x, t ) = c a Ψ xx ( x, t ) − γ a Ψ t ( x, t ) , (3b)where again Φ is the twist angle, Ψ is the axial movement, c t = (cid:112) G/ρ is the propagation speed of the angle, γ t is theinternal damping, c a = (cid:112) E/ρ is the axial velocity and γ a is the axial distributed damping. A list of physical parametersand their values is given in Table II and Figure 1 helps giving abetter understanding of the physical system. In other words, if Ψ( · , t ) = 0 , then there is no compression in the pipe, meaningthat the bit is not bouncing; if Φ tt ( · , t ) = 0 , then the angularspeed along the pipe is the same, meaning that there is noincrease or decrease of the torsion.For the previous model to be well-posed, top and bottomboundary conditions (at x = 0 and x = L ) must beincorporated in (3). In this part, only the topside boundarycondition is derived. There is a viscous damping at x = 0 ,and consequently a mismatch between the applied torque atthe top and the angular speed. The topside boundary conditionfor the axial part is built on the same scheme and the followingarameter meaning Value L Pipe length m G Shear modulus . × N.m − E Young modulus × N.m − Γ Drillstring’s cross-section × − m J Second moment of inertia . × − m I B Bottom hole lumped inertia kg.m M B Bottom hole mass
40 000 kg ρ Density kg.m − g Angular momentum
N.m.s.rad − h Viscous friction coefficient kg.s − γ a Distributed axial damping . s − γ t Distributed angle damping . s − δ Weight on bit coefficient m − TABLE II: Physical parameters, meanings and their values [1],[40].conditions are obtained for t > : GJ Φ x (0 , t ) = g Φ t (0 , t ) − u ( t ) , (4a) E ΓΨ x (0 , t ) = h Ψ t (0 , t ) − u ( t ) . (4b)The downside boundary condition ( x = L ) is more difficultto grasp and is consequently derived later when dealing witha more complex model. C. Neutral-type time-delay model
Studying an infinite-dimensional problem stated in terms ofPDEs represents a relevant challenge. The equations obtainedpreviously are damped wave equations, but, for the specialcase where γ a = γ t = 0 , the system can be convertedinto a neutral time-delay system as done in [38]. This newformulation enables to use other tools to analyze its stabilityas the Lyapunov-Krasovskii Theorem or a frequency domainapproach making its stability analysis slightly easier.Nevertheless, the main drawback of this formulation is theassumption that the damping occurs at the boundary andnot all along the pipe. This useful simplification, even if itencountered in many articles [12], [38], [39], is known tochange in a significant manner the behavior of the system[1]. Indeed, it appears that without internal damping, the waveequation rephrases easily as a system of transport equations.It is then directly possible to observe with a delay of c − atthe top of the pipe what happened at the bottom of the pipe.This makes the control easier. D. Coupled ODE/PDE model
To overcome the issue mentioned previously, a simplermodel than the one derived in (3) is proposed in [39], where an harmonic oscillator is used to describe axial vibrations andthe model results in a coupled ODE/PDE.A second possibility, reported in [14], [38] for example is topropose a second order ODE as the bottom boundary condition( x = L ) for t > : GJ Φ x ( L, t ) = − I B Φ tt ( L, t ) − T (Φ t ( L, t )) , (5a) E ΓΨ x ( L, t ) = − M B Ψ tt ( L, t ) − δT (Φ t ( L, t )) , (5b)where T represents the torque applied on the drilling bit bythe rocks, described in equation (2). Notice that equation (5a)is coming from the conservation of angular momentum where GJ Φ x ( L, t ) is the torque coming from the top of the pipe.Equation (5b) is the direct application of Newton’s secondlaw of motion where E ΓΨ x ( L, t ) is the force transmitted fromthe top to the bit and δT (Φ t ( L, t )) is the weight on bit dueto the rock interaction. Since (5) is a second order in timedifferential equation, note that (3) together with (5) indeedleads to a coupled ODE/PDE.There exist other bottom boundary conditions leading to amore complex coupling between axial and torsional dynamics.They nevertheless introduce delays which requires to havea better knowledge of the drilling bit. To keep the contentgeneral, the boundary conditions (5) used throughout thispaper 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 whichbackstepping controllers can be used [12], [37]. This is themain reason why this model is widely used today. E. Models comparison
We propose in this subsection to compare the coupledODE/PDE model and the lumped parameter models for thetorsion only. We consider here a linearization of the systemfor large Ω and consequently we neglect the stick-slip effectby setting T = 0 .First, denote by H DP M the transfer function from u to Φ( L, · ) for the DPM and H LP M from u to Φ b for the LPM.We also define by H LP M a truncation of H LP M consideringonly the two dominant poles. The Bode diagrams of H DP M , H LP M and H LP M are drawn in Figure 4.Clearly, the LPMs catch the behavior of the DPM at steadystates and low frequencies until the resonance, occurringaround (cid:112) k/I b rad.s − . From a control viewpoint, the DPMhas infinitely many harmonics as it can be seen on the plots butof lower magnitudes and damped as the frequency increases(around − dB at each decade). The magnitude plots isnot sufficient to make a huge difference between the threemodels. Nevertheless, considering the phase, we see a cleardifference. It appears that the DPM crosses the frequency − many times making the control margins quite difficult toassess. Moreover, that shows that the DPM is harder to controlbecause of the huge difference of behavior after the resonance.They may consequently has a very different behavior whencontrolled. That is why we focus in this study on the the DPM,even if it is far more challenging to control than the LPMs. -3 -2 -1 -200-150-100-50010 -3 -2 -1 -400-300-200-1000 Fig. 4: Bode diagram of H DP M , H LP M and H LP M .III. P ROBLEM STATEMENT
The coupled nonlinear ODE/PDE model derived in theprevious subsection can be written as follows for t > and x ∈ [0 , : φ tt ( x, t ) = ˜ c t φ xx ( x, t ) − γ t φ t ( x, t ) ,φ x (0 , t ) = ˜ gφ t (0 , t ) − ˜ u ( t ) ,φ t (1 , t ) = z ( t ) , ˙ z ( t ) = − α φ x (1 , t ) − α T ( z ( t )) , (6) ψ tt ( x, t ) = ˜ c a ψ xx ( x, t ) − γ a ψ t ( x, t ) ,ψ x (0 , t ) = ˜ hψ t (0 , t ) − ˜ u ( t ) ,ψ t (1 , t ) = y ( t ) , ˙ y ( t ) = − β ψ x (1 , t ) − β T ( z ( t )) , (7)where the normalized parameters are given in Table III. Notethat the range for the spacial variable is now x ∈ [0 , to easethe calculations. The initial conditions are as follows: φ ( x,
0) = φ ( x ) , φ t ( x,
0) = φ ( x ) ,ψ ( x,
0) = ψ ( x ) , ψ t ( x,
0) = ψ ( x ) ,z (0) = φ (1) , y (0) = ψ (1) . Since the existence and uniqueness of a solution to theprevious problem is not the main contribution of this paper,it is assumed in the sequel. This problem has been widelystudied (see [10], [12], [37], [39], [44] among many others)and the solution belongs to the following space if the initialconditions ( φ , φ , ψ , ψ ) satisfies the boundary conditions(see [9] for more details): (cid:88) = H × L , (cid:88) = H × H , ( φ, φ t , ψ, ψ t , z , y ) ∈ C ( (cid:88) × (cid:82) ) . Remark 1:
One may note that θ (cid:55)→ T nl ( θ ) is not well defined Parameter Expression Value φ ( x, t ) Φ( xL, t ) - ψ ( x, t ) Ψ( xL, t ) - ˜ c t c t L − . c a c a L − . α GJLI B . β E Γ LM B . α I − B . · − β δM B . · − ˜ g gGJ . · − ˜ h hE Γ 2 . · − ˜ u ( t ) 1 GJu ( t ) - ˜ u ( t ) 1 E Γ u ( t ) -TABLE III: Normalized parameters.for θ = 0 because of the sign function. Nevertheless, sincethe nonlinearity acts directly on the variable z , it follows thatthere exists a unique solution to the ODE system in the senseof Filipov [19]. A more detailed discussion on this point isprovided in [11]. (cid:3) System (6)-(7) is a cascade of two subsystems:1) System (6) is a coupled nonlinear ODE/string equationdescribing the torsion angle φ .2) System (7) is a coupled linear ODE/string equationsubject to the external perturbation T ( z ) for Ψ .It appears clearly that the perturbation on the second subsys-tem depends on the first subsystem in φ . Since they are verysimilar, the same analysis than the one conducted in this paperapplies for the second subproblem. That is why we only studythe evolution of the torsion.IV. E XPONENTIAL S TABILITY OF THE L INEAR S YSTEM
System (6) is a nonlinear system because of the friction term T nl introduced in (2). Nevertheless, for a high desired angularspeed Ω , T nl can be assumed constant as seen in Figure 3.Moreover, studying this linear system can be seen as a firststep before studying the nonlinear system, which relies mostlyon the stability theorem derived in this section.The proposed linear model of T around Ω (cid:29) is T ( θ ) = c b θ + T nl (Ω ) = c b θ + T . (8)For high values of Ω , T = T nl (Ω ) is close to T smooth (Ω ) and at the limit when Ω tends to infinity, they are equal.Hence the nonlinear friction term for relatively large angularvelocity does not influence much the system.n our case, the proposed controller for this problem hasbeen studied in a different setting in [13], [15], [45] and isa simple proportional/integral controller based on the singlemeasurement of the angular velocity at the top of the drill (i.e. φ t (0 , t ) ). The following variables are therefore introduced: ˜ u ( t ) = − k p ( φ t (0 , t ) − Ω ) − k i z ( t ) , ˙ z ( t ) = φ t (0 , t ) − Ω , (9)where k p and k i are the gains of the PI controller. Combiningequations (6), (8) and (9) leads to: φ tt ( x, t ) = c φ xx ( x, t ) − γ t φ t ( x, t ) ,φ x (0 , t ) = (˜ g + k p ) φ t (0 , t ) − k p Ω + k i C Z ( t ) ,φ t (1 , t ) = C Z ( t ) , ˙ Z ( t ) = AZ ( t ) + B (cid:104) φ t (0 ,t ) φ x (1 ,t ) (cid:105) + B (cid:2) Ω T (cid:3) , (10)where A = − c b I B
00 0 , B = − α , B = − α − ,Z = (cid:104) z z (cid:105) (cid:62) , C = (cid:104) (cid:105) , C = (cid:104) (cid:105) . To ease the notations, c = ˜ c t . We denote by X = ( φ x , φ t , z , z ) ∈ C ([0 , ∞ ) , H ) with H = L × L × (cid:82) the infinite dimensional state ofsystem (10). The control objective in the linear case is toachieve the exponential stabilization of an equilibrium pointof (10) in angular speed, i.e. φ t (1) is going exponentially to agiven constant reference value Ω . To this extent, the followingnorm on H is introduced: (cid:107) X (cid:107) H = z + z + c (cid:107) φ x ( · ) (cid:107) + (cid:107) φ t ( · ) (cid:107) . The definition of an equilibrium point of (10) and its expo-nential stability follows from the previous definitions.
Definition 1: X ∞ ∈ H is an equilibrium point of (10)if for the trajectory X ∈ C ([0 , ∞ ) , H ) of (10) with initialcondition X ∞ , the following holds: ∀ t > , (cid:107) ˙ X ( t ) (cid:107) H = 0 . Definition 2:
Let X ∞ be an equilibrium point of (10). X ∞ is said to be µ -exponentially stable if ∀ t > , (cid:107) X ( t ) − X ∞ (cid:107) H ≤ γ (cid:107) X − X ∞ (cid:107) H e − µt holds for γ (cid:62) , µ > and for any initial conditions X ∈ H satisfying the boundary conditions. Here, X is the trajectoryof (10) whose initial condition is X .Hence, X ∞ is said to be exponentially stable if there exists µ > such that X ∞ is µ -exponentially stable.Before stating the main result of this part, a lemma aboutthe equilibrium point is proposed. Lemma 1:
Assume k i (cid:54) = 0 , then there exists a uniqueequilibrium point X ∞ = ( φ ∞ x , φ ∞ t , z ∞ , z ∞ ) ∈ H of (10) andit satisfies φ ∞ t = Ω . Proof.
An equilibrium point X ∞ of (10) is such that: ( φ ∞ xt , φ ∞ tt , ddt z ∞ , ddt z ∞ ) = 0 .Since ddt z ∞ = φ ∞ t (0 , t ) − Ω = 0 and ∂ x φ ∞ t = 0 and ∂ t φ ∞ t = 0 , φ ∞ t = Ω holds.We also get from φ ∞ xx = 0 and ∂ t φ ∞ x = 0 that φ ∞ x is a firstorder polynomial in x . Together with the boundary conditions,this system has a unique solution if k i (cid:54) = 0 such as: X ∞ = (cid:18) φ ∞ x , Ω , Ω , φ ∞ x (0) − ˜ g Ω k i (cid:19) , where φ ∞ x ( x ) = γ t Ω c ( x − − c b α I B Ω − α α T for x ∈ [0 , . A. Exponential stability of the closed-loop system (10)The main result of this part is then stated as follows.
Theorem 1:
Let N ∈ (cid:78) . Assume there exists P N ∈ (cid:83) N +1) , R = diag( R , R ) (cid:23) , S = diag( S , S ) (cid:31) , Q ∈ (cid:83) such that the following LMIs hold: Θ N = c Θ ,N + Θ ,N − Q N ≺ ,P N + S N (cid:31) , (11)and Γ = cR + γ t U − Q (cid:23) , Γ = cR + γ t U − Q (cid:23) , (12)where Θ ,N = H (cid:62) N (cid:2) S + R − S (cid:3) H N − G (cid:62) N (cid:2) S − S − R (cid:3) G N , Θ ,N = He (cid:0) D (cid:62) N P N F N (cid:1) ,F N = (cid:104) I N +1) N +1) , (cid:105) ,D N = (cid:104) J (cid:62) N cM (cid:62) N (cid:105) (cid:62) , J N = (cid:104) A , N +1) B (cid:105) ,M N = (cid:49) N H N − ¯ (cid:49) N G N − (cid:104) N +1) , L N N +1) , (cid:105) ,U = (cid:104) S S + S + R S + S + R S + R ) (cid:105) , U = (cid:104) S + R ) S + S + R S + S + R S (cid:105) , (13) G N = (cid:104) ck i C − ck i C , N +1) G (cid:105) , G = (cid:104) c (˜ g + k p ) 01 − c (˜ g + k p ) 0 (cid:105) ,H N = (cid:104) C C , N +1) H (cid:105) , H = (cid:2) c − c (cid:3) ,Q N = diag(0 , Q, Q, · · · , (2 N + 1) Q, ) ,S N = diag (0 , S, S, · · · , (2 N + 1) S ) ,L N = [ (cid:96) j,k Λ] j,k ∈ [0 ,N ] − γ t ] , . . . , [ ]) , (cid:49) N = (cid:34) Λ ... Λ (cid:35) , ¯ (cid:49) N = (cid:34) Λ ... ( − N Λ (cid:35) , Λ = c − c , and (cid:96) k,j = (2 j + 1)(1 − ( − j + k ) , if j (cid:54) k, otherwise , then the equilibrium point of system (10) is exponentiallystable. As a consequence, lim t → + ∞ | φ t (1 , t ) − Ω | = 0 .he proof is given thereafter but some practical conse-quences and preliminary results are firstly derived. Remark 2:
The methodology used to derive the previousresult has been firstly introduced in [42] to deal with time-delay systems. It has been proved in the former article that thistheorem provides a hierarchy of LMI conditions. That meansif the conditions of Theorem 1 are met for N = N (cid:62) ,then the conditions are also satisfied for all N > N . Then,the LMIs provide a sharper analysis as the order N of thetheorem increases. Nevertheless, the price to pay for a moreprecise analysis is the increase of the number of decisionvariables and consequently a higher computational burden. Ithas been noted in [42] for a time-delay system and in [8] fora coupled ODE/string equation system that very sharp resultsare obtained even for low orders N . (cid:3) As we aim at showing that X ∞ is an exponentially stableequilibrium point of system (10) in the sense of Definition 2,the following variable is consequently defined for t ≥ : ˜ X ( t ) = X ( t ) − X ∞ = (cid:16) ˜ φ x ( t ) , ˜ φ t ( t ) , ˜ z ( t ) , ˜ z ( t ) (cid:17) , where X is a trajectory of system (10).The following lemma, given in [7], provides a way forestimating the exponential decay-rate of system (10) as soonas a Lyapunov functional V is known. Lemma 2:
Let V be a Lyapunov functional for system (10)and µ ≥ . Assume there exist ε , ε , ε > such that thefollowing holds: ε (cid:107) ˜ X (cid:107) H (cid:54) V ( ˜ X ) (cid:54) ε (cid:107) ˜ X (cid:107) H , ˙ V ( ˜ X ) + 2 µV ( ˜ X ) (cid:54) − ε (cid:107) ˜ X (cid:107) H , (14)then the equilibrium point of system (10) is µ -exponentiallystable. If µ = 0 , then it is exponentially stable. Proof of Theorem 1.
For simplicity, the following nota-tions are used throughout the paper for k ∈ (cid:78) : ˜ χ ( x ) = ˜ φ t ( x ) + c ˜ φ x ( x )˜ φ t ( x ) − c ˜ φ x ( x ) , ˜ X k ( t ) = (cid:90) ˜ χ ( x, t ) L k ( x ) dx, (15)where {L k } k ∈ (cid:78) is the orthogonal family of Legendre polyno-mials as defined in Appendix A. ˜ X k is then the projectioncoefficient of ˜ χ along the Legendre polynomial of degree k . ˜ χ refers to the Riemann coordinates and presents usefulproperties as discussed in [8] for instance.
1) Choice of a Lyapunov functional candidate:
The pro-posed Lyapunov functional is inspired from [7], [8], [36] andis divided into two parts as follows: V N ( ˜ X ) = ˜ Z (cid:62) N P N ˜ Z N + V ( ˜ χ ) (16)where ˜ Z N = (cid:104) ˜ z ˜ z ˜ X (cid:62) · · · ˜ X (cid:62) N (cid:105) (cid:62) and V ( ˜ χ ) = (cid:90) ˜ χ (cid:62) ( x ) (cid:104) S + xR S +(1 − x ) R (cid:105) ˜ χ ( x ) dx. V is a traditional Lyapunov functional candidate for a transportsystem as shown in [16].
2) Exponential stability: Existence of ε : This part is inspired by [7]. The inequalities P N + S N (cid:31) , R (cid:23) imply the existence of ε > such that: P N + S N (cid:23) ε diag (cid:0) I , I N (cid:1) ,S (cid:23) ε I , (17)with I N = diag { (2 k + 1) I } k ∈ [0 ,N ] . This statement impliesthe following on V N : V N ( ˜ X ) ≥ ˜ Z (cid:62) N ( P N + S N ) ˜ Z N − N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k S ˜ X k + (cid:90) ˜ χ ( x ) (cid:62) (cid:16) S − ε I (cid:17) ˜ χ ( x ) dx + ε (cid:107) ˜ χ (cid:107) ≥ ε (cid:0) ˜ z + ˜ z + (cid:107) ˜ χ (cid:107) (cid:1) − N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k ˜ S ˜ X k + (cid:90) ˜ χ (cid:62) ( x ) ˜ S ˜ χ ( x ) dx ≥ ε (cid:107) ˜ X (cid:107) H , where ˜ S = S − ε I which ends the proof of existence of ε . Existence of ε : Following the same line as previously, thefollowing holds for a sufficiently large ε > : P N (cid:22) ε diag (cid:0) I , I N (cid:1) , (cid:104) S + xR S +(1 − x ) R (cid:105) (cid:22) ε I . Using these inequalities in (16) leads to: V N ( ˜ X ) ≤ ε (cid:32) ˜ z + ˜ z + N (cid:88) k =0 k + 14 ˜ X (cid:62) k ˜ X k + 14 (cid:107) ˜ χ (cid:107) (cid:33) ≤ ε (cid:0) ˜ z + ˜ z + (cid:107) ˜ χ (cid:107) (cid:1) = ε (cid:107) ˜ X (cid:107) H . The last inequality is a direct application of Lemma 4.
Existence of ε : The time derivation of ˜ χ leads to: ˜ χ t ( x ) = Λ ˜ χ x ( x ) − γ t [ ] ˜ φ t ( x ) . Noting that ˜ φ t ( x ) = [ ] ˜ χ ( x ) , we get: ˜ χ t ( x ) = Λ ˜ χ x ( x ) − γ t ] ˜ χ ( x ) . (18)The derivative of V along the trajectories of (6) leads to: ˙ V ( ˜ χ ) = 2 c V ( ˜ χ ) − γ t V ( ˜ χ ) , with V ( ˜ χ ) = (cid:90) ˜ χ (cid:62) x ( x ) (cid:104) S + xR − S − (1 − x ) R (cid:105) ˜ χ ( x ) dx, V ( ˜ χ ) = (cid:90) ˜ χ (cid:62) ( x ) U ( x ) ˜ χ ( x ) dx,U ( x ) = (cid:104) S + xR ) S + S + xR +(1 − x ) R S + S + xR +(1 − x ) R S +(1 − x ) R ) (cid:105) . An integration by part on V shows that: V ( ˜ χ ) = ˜ χ (cid:62) (1) (cid:2) S + R − S (cid:3) ˜ χ (1) − ˜ χ (cid:62) (0) (cid:2) S − S − R (cid:3) ˜ χ (0) − (cid:90) ˜ χ (cid:62) ( x ) (cid:2) R R (cid:3) ˜ χ ( x ) dx. he previous calculations lead to the following derivative: ˙ V ( ˜ χ ) = c (cid:0) ˜ χ (cid:62) (1) (cid:2) S + R − S (cid:3) ˜ χ (1) − ˜ χ (cid:62) (0) (cid:2) S − S − R (cid:3) ˜ χ (0) (cid:17) − (cid:90) ˜ χ (cid:62) ( x ) (cid:16) cR + γ t U ( x ) (cid:17) ˜ χ ( x ) dx. (19)By a convexity argument, if (12) is verified then cR + U ( x ) (cid:23) Q (cid:31) holds for x ∈ [0 , . Consequently, noticing that ˜ χ (0) = G N ˜ ξ N , ˜ χ (1) = H N ˜ ξ N , we get: ˙ V ( ˜ χ ) ≤ c ˜ ξ (cid:62) N Θ ,N ˜ ξ N − (cid:90) ˜ χ (cid:62) ( x ) Q ˜ χ ( x ) dx, with ˜ ξ N = (cid:104) ˜ Z (cid:62) N ˜ φ t (0) ˜ φ x (1) (cid:105) (cid:62) . Using Lemma 5 and the fact that ˙˜ Z N = D N ˜ ξ N and ˜ Z N = F N ˜ ξ N , we get the following: ˙ V N ( ˜ X ) = He (cid:16) ˙˜ Z (cid:62) N P N ˜ Z N (cid:17) + ˙ V ( ˜ χ ) ≤ ˜ ξ (cid:62) N Θ N ˜ ξ N + N (cid:88) k =0 (2 k + 1) ˜ X k Q ˜ X k − (cid:90) ˜ χ (cid:62) ( x ) Q ˜ χ ( x ) dx, (20)with Θ N defined in (11). Since Θ N ≺ and Q (cid:31) , we get: Θ N (cid:22) − ε diag (cid:0) I , I , I , . . . , N +12 I , (cid:1) ,Q (cid:23) ε I . Then, ˙ V N is upper bounded by: ˙ V N ( ˜ X ) ≤ − ε (cid:0) ˜ z + ˜ z + (cid:107) ˜ χ (cid:107) (cid:1) + N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k (cid:16) Q − ε I (cid:17) ˜ X k − (cid:90) ˜ χ (cid:62) ( x ) (cid:16) Q − ε I (cid:17) ˜ χ ( x ) dx ≤ − ε (cid:107) ˜ X (cid:107) H . The last inequality comes from a direct application of Bessel’sinequality (31).
Conclusion :
Using Lemma 2, we indeed get the exponentialconvergence of all trajectories of (10) to the desired equilib-rium point.
B. Exponential stability with a guaranteed decay-rate
It is possible to estimate the decay-rate µ of the exponentialconvergence with a slight modification of the LMIs as it isproposed in the following corollary. Corollary 1:
Let N ∈ (cid:78) , µ > and γ t ≥ . If thereexist P N ∈ (cid:83) N +1) , R = diag( R , R ) (cid:23) , S =diag( S , S ) (cid:31) , Q ∈ (cid:83) such that (12) and the following LMIs hold: Θ N,µ = Θ N + 2 µF (cid:62) N ( P N + S N ) F N ≺ ,P N + S N (cid:31) , (21)with the parameters defined as in Theorem 1, then the equi-librium point of system (10) is µ -exponentially stable. Proof.
To prove the exponential stability with a decay-rateof at least µ > , we use Lemma 2. Similarly to the previousproof, we have the existence of ε , ε > . The existenceof ε is slightly different. First, note that for the Lyapunovfunctional candidate (16), we get: V N ( ˜ X φ ) ≥ ˜ Z (cid:62) N P N ˜ Z N + (cid:90) ˜ χ (cid:62) ( x ) S ˜ χ ( x ) dx ≥ ˜ Z (cid:62) N P N ˜ Z N + N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k S ˜ X k . This inequality was obtained using equation (31). Using thenotations of the previous theorem yields: V N ( ˜ X φ ) ≥ ˜ ξ (cid:62) N F (cid:62) N ( P N + S N ) F N ˜ ξ N . (22)Coming back to (20) and using (21) leads to: ˙ V N ( ˜ X φ ) ≤ ˜ ξ (cid:62) N Θ N,µ ˜ ξ N − µ ˜ ξ (cid:62) N (cid:0) F (cid:62) N P N F N + S N (cid:1) ˜ ξ N + N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k Q ˜ X k − (cid:90) ˜ χ (cid:62) ( x ) Q ˜ χ ( x ) dx. Injecting inequality (22) into the previous inequality leads to: ˙ V N ( ˜ X φ ) + 2 µV N ( ˜ X φ ) ≤ ˜ ξ (cid:62) N Θ N,µ ˜ ξ N + N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k Q ˜ X k − (cid:90) ˜ χ (cid:62) ( x ) Q ˜ χ ( x ) dx. Using the same techniques than for the previous proof leadsto the existence of ε > such that: ˙ V N ( ˜ X φ ) + 2 µV N ( ˜ X φ ) ≤ − ε (cid:107) ˜ X φ (cid:107) H . Lemma 2 concludes then on the µ -exponential stability. Remark 3:
Wave equations can sometimes be modeled asneutral time-delay systems [6], [38]. This kind of system isknown to possess some necessary stability conditions as notedin [6], [9]. In the book [9], the following criterion is derived: µ ≤ c (cid:12)(cid:12)(cid:12)(cid:12) c (˜ g + k p )1 − c (˜ g + k p ) (cid:12)(cid:12)(cid:12)(cid:12) = µ max . (23)This result implies that there exists a maximum decay-rate andif this maximum is negative, then the system is unstable. TheLMI Θ µN ≺ contains the same necessary condition, meaningthat the neutral aspect of the system is well-captured. (cid:3) 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 stabilitycondition (equivalent to (23)) coming from frequency analysisan be derived (Cor 3.3 from [21]): (cid:12)(cid:12)(cid:12)(cid:12) − c ˜ g c ˜ g (cid:12)(cid:12)(cid:12)(cid:12) + 2 (cid:12)(cid:12)(cid:12)(cid:12) ck p c ˜ g (cid:12)(cid:12)(cid:12)(cid:12) < . It is more restrictive and taking k p (cid:54) = 0 practically leads to adecrease of the robustness (even if some other performancesmight be enhanced). This phenomenon has been studied inmany articles [22], [29], [44]. Hence, to be robust to delaysin the loop, one then needs to ensure the following: ≤ k p ≤ c ( | c ˜ g | − | − c ˜ g | ) = ˜ g = 2 . · − . This inequality on k p comes when considering the infinitedimensional problem and does not arise when dealing withany finite dimensional model. That point enlightens that it ismore realistic to consider the infinite dimensional problem.For more information on that point, the interested reader canrefer to [5].V. P RACTICAL S TABILITY OF SYSTEM (6)-(9)The experiments conducted previously shows for Ω large,the trajectory of the nonlinear system (6)-(9) goes exponen-tially to X ∞ , so does the linear system. In other words, theprevious result is a local stability test for the nonlinear systemin case of large desired angular velocity Ω and does notextend straightforward to a global analysis.In general, for a nonlinear system, ensuring the global expo-nential stability of an equilibrium point could be complicated.In many engineering situations, global exponential stability isnot the requirement. Indeed, considering uncertainties in thesystem or nonlinearities, it is far more reasonable and accept-able for engineers to ensure that the trajectory remains close tothe equilibrium point. This property is called dissipativenessin the sense of Levinson [26] or practical stability (also calledglobally uniformly ultimately bounded in [24]). Definition 3:
System (6) is practically stable if there exists X bound ≥ such that for X the solution of (6) with initialcondition X ∈ H : ∀ η > , ∃ T η > , ∀ t ≥ T η , (cid:107) X ( t ) − X ∞ (cid:107) H ≤ X bound + η. (24)Saying that system (6)-(9) is practically stable means thatthere exists T η > such that for any x ∈ [0 , , φ t ( x, t ) staysclose to Ω for t ≥ T η . This property has already been appliedto a drilling system in [39] for instance. In the nonlinear case,the aim is then to design a control law reducing the amplitudeof the stick-slip when it occurs.The objective of this section is to derive an LMI testensuring the practical stability of nonlinear system (6) togetherwith the control defined in (9): φ tt ( x, t ) = c φ xx ( x, t ) − γ t φ t ( x, t ) ,φ x (0 , t ) = (˜ g + k p ) φ t (0 , t ) − k p Ω + k i C Z ( t ) ,φ t (1 , t ) = C Z ( t ) , ˙ Z ( t ) = AZ ( t ) + B (cid:104) φ t (0 ,t ) φ x (1 ,t ) (cid:105) + B (cid:104) Ω T nl ( z ( t )) (cid:105) , (25) Remark 4:
To simplify the writing, φ can both meansthe solution of the linear or nonlinear system depending on the context. In this part, it refers to a solution of nonlinearsystem (25). (cid:3) A. Practical stability of system (25)The idea behind practical stability is to embed the staticnonlinearities by the use of suitable sector conditions, as ithas been done for the saturation for instance [43]. Then, theuse of robust tools will lead to some LMI tests ensuring thepractical stability.
Lemma 3:
For almost all ˜ z ∈ (cid:82) , the following holds: T nl (˜ z + Ω ) ≤ T max , T min ≤ T nl (˜ z + Ω ) , − z + Ω ) T nl (˜ z + Ω ) ≤ . (26) Proof.
These inequalities can be easily verified using (2).These new information are the basis of the followingtheorem.
Theorem 2:
Let N ∈ (cid:78) and V max > . If there exist P N ∈ (cid:83) N +1) , R = diag( R , R ) (cid:23) , S = diag( S , S ) (cid:31) , Q ∈ (cid:83) , τ , τ , τ , τ ≥ such that (12) holds togetherwith: Ξ N = ¯Θ N − τ Π − τ Π − τ Π − τ Π ≺ ,P N + S N (cid:31) , (27)where ¯Θ N = diag(Θ N , ) − α He (cid:16) ( F m − T F m ) (cid:62) e (cid:62) P N ˜ F N (cid:17) , ˜ F N = (cid:2) F N N +1)+2 , (cid:3) , e = (cid:2) , N +1)+1 (cid:3) (cid:62) ,F m = (cid:2) , N +1)+4 (cid:3) , F m = (cid:2) , N +1)+4 (cid:3) , Π = V max F (cid:62) m F m − ˜ F (cid:62) N ( P N + S N ) ˜ F N , Π = π (cid:62) π − T max π (cid:62) π , Π = T min π (cid:62) π − π (cid:62) π , Π = − He(( π + Ω π ) (cid:62) π ) ,π = (cid:2) , , N +1)+3 , , (cid:3) , π = (cid:2) , N +1)+4 , , (cid:3) ,π = (cid:2) , N +1)+4 , , (cid:3) , and all the parameters defined as in Theorem 1, then theequilibrium point X ∞ of system (25) is practically stable.More precisely, equation (24) holds for X bound = (cid:113) V max ε − where ε is defined in (17). Proof.
First, let do the same change of variable as before: ˜ X = X − X ∞ . Since the nonlinearity affects only the ODE part of sys-tem (25), the difference with the previous part lies in thedynamic of ˜ Z : ˙˜ Z ( t ) = ddt (cid:16) Z ( t ) − (cid:104) z ∞ z ∞ (cid:105)(cid:17) = A ˜ Z ( t ) + B (cid:104) ˜ φ t (0 ,t )˜ φ x (1 ,t ) (cid:105) + B (cid:2) T nl (˜ z ( t )+Ω ) − T (cid:3) . Using the same Lyapunov functional as in (16), the pos-itivity is ensured in the exact same way. For the bound onhe time derivative, following the same strategy as before, weeasily get: ˙ V N ( ˜ X ) ≤ ˜ ξ (cid:62) N Θ N ˜ ξ N − α ( T nl (˜ z + Ω ) − T ) e (cid:62) P N ˜ Z N + N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k Q ˜ X k − (cid:90) ˜ χ (cid:62) ( x ) Q ˜ χ ( x ) dx. (28)We introduce a new extended state variable ¯ ξ N = (cid:104) ˜ ξ (cid:62) N T nl (˜ z + Ω ) 1 (cid:105) (cid:62) . Using the notation of Theorem 2,equation (28) rewrites as: ˙ V N ( ˜ X ) ≤ ¯ ξ (cid:62) N ¯Θ N ¯ ξ (cid:62) N + N (cid:88) k =0 (2 k + 1) ˜ X (cid:62) k Q ˜ X k − (cid:90) ˜ χ (cid:62) ( x ) Q ˜ χ ( x ) dx. It is impossible to ensure ¯Θ N ≺ since its last by diagonalblock is .We then use the definition of practical stability. We want toshow that if there exists ε > such that ˙ V N ( ˜ X ) ≤ − ε (cid:107) ˜ X (cid:107) H when V N ( ˜ X ) ≥ V max then the system is practically stable.Let S = (cid:110) ˜ X ∈ H | V N ( ˜ X ) ≤ V max (cid:111) , the previous assertionimplies that this set is invariant and attractive. Using (17),we get that (cid:110) ˜ X ∈ H | (cid:107) ˜ X (cid:107) H ≤ X bound = V / max ε − / (cid:111) ⊇ S ,meaning that the system is practically stable.A sufficient condition to be practically stable is that V N shouldbe strictly decreasing when outside the ball of radius V max .This condition rewrites as V N ( ˜ X ) ≥ V max and the followingholds: V N ( ˜ X ) − V max ≥ ˜ Z (cid:62) N P N ˜ Z N + (cid:90) ˜ χ (cid:62) ( x ) S ˜ χ ( x ) dx − V max ≥ − ¯ ξ (cid:62) N Π ¯ ξ N ≥ . The previous inequality is obtained using Bessel inequal-ity (31) on (cid:82) ˜ χ (cid:62) ( x ) S ˜ χ ( x ) dx . Hence, V N ( ˜ X ) ≥ V max if ¯ ξ (cid:62) N Π ¯ ξ N ≤ .Noting that ˜ z = π ¯ ξ N , T nl (˜ z + Ω ) = π ¯ ξ N and π ¯ ξ N , Lemma 3 rewrites as: ∀ i ∈ { , , } , ¯ ξ N Π i ¯ ξ N ≤ . Consequently, a sufficient condition to be practically stable is: ∀ ¯ ξ N (cid:54) = 0 s. t. ∀ i ∈ [0 , , ¯ ξ (cid:62) N Π i ¯ ξ N ≤ , ¯ ξ (cid:62) N ¯Θ N ¯ ξ N < . (29)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) isverified if there exists τ , τ , τ , τ > such that ¯Θ N − (cid:88) k =0 τ i Π i ≺ . Consequently, in a similar way than for Theorem 1, condi-tion (27) implies that S is an invariant and attractive set forsystem (25). Then, the equilibrium point X ∞ of system (25)is practically stable with X bound = (cid:113) V max ε − . Note that if the torque function is not perfectly known, onecan change the lower and upper bounds T min and T max to geta more conservative result but robust to uncertainties on T nl .For (27) to be feasible, the constraint ¯ ξ (cid:62) N Π ¯ ξ N ≤ must hold,meaning that an upper bound of T nl needs to be proposed. B. On the optimization of X bound The condition (27) is a Bilinear Matrix Inequality (BMI)since τ , τ , τ and τ are decision variables and it is thereforedifficult to get its global optimum. Nevertheless, the followinglemma gives a sufficient condition for the existence of asolution to this problem. Corollary 2:
There exists V max and τ > such thatTheorem 2 holds if and only if there exists N > such thatLMIs (12) and (11) are satisfied.In other words, the equilibrium point of system (25) ispractically stable if and only if the linear system (10) isexponentially stable. Proof.
Note first that expending (27) with τ = 0 leads to: Ξ N = Θ N,τ / κ P N κ (cid:62) P N τ − τ τ T max − τ T min − τ V max , (30)where κ P N ∈ (cid:82) (4+2( N +1)) × depends only on P N . Proof of sufficiency:
Assume there exists
N > such thatLMIs (12) and (11) are satisfied. Considering τ = 0 and usingSchur complement on Ξ N , Ξ N ≺ is equivalent to: Θ N,τ / − κ P N − τ τ T max − τ V max κ (cid:62) P N ≺ , with τ V max > τ T max . Since Θ N ≺ , considering τ small enough, τ large and V max > τ τ − T max the previouscondition is always satisfied and Theorem 2 applies. Proof of necessity:
Assume Ξ N ≺ and (12) holds. Thenits first diagonal block must be definite negative. Consequently Θ N ≺ and, according to Theorem 1, system (10) isexponentially stable. Remark 5:
Note that (30) provides a necessary conditionfor Theorem 2 which is Θ N,τ / ≺ . In other words, τ isrelated to the decay-rate of the linear system and we get thefollowing condition: τ < µ max . (cid:3) Thanks to Corollary 2, the following method should helpsolving the BMI if Theorem 1 is verified. Assuming thatequations (11) and (12) are verified for a given N ∈ (cid:78) .1) Fix τ = 2 µ max as defined in (23).2) Check that equations (11), (12) and (27) are satisfied for V max a strictly positive decision variable. If this is notthe case, then decrease τ and do this step again.3) Thanks to Corollary 2, there exists a τ small enoughfor which equations (11), (12) and (27) hold. Freeze thisvalue.rder N = 0 N = 1 N = 2 N = 3 N = 6 µ ( × − ) 0 .
87 4 .
24 7 .
31 7 .
59 7 . TABLE IV: Estimated decay-rate as a function of the order N used. Note that µ max = 7 . · − is calculated using (23)for k p = 10 − , k i = 10 .4) Since the problem is unbounded, it is possible to fix avariable without loss of generality, let V max = 10 forinstance and solve the following optimization problem: min P N ,S,R,Q,τ ,τ ,τ ,ε P − ε P subject to (11) , (12) and (27) ,P N + S N − diag( ε P , N +3 ) (cid:31) ,R (cid:23) , S (cid:31) , Q (cid:31) .
5) Then compute X bound = (cid:113) V max ε − where ε isdefined in (17).VI. E XAMPLES AND S IMULATIONS
The section is devoted to numerical simulations and drawsome conclusions about the PI regulation. In the first subsec-tion, we focus on the linear system and the second one isdedicated to the nonlinear case. A. On the linear model
This subsection recaps the result on the linear model.
1) Estimation of the decay-rate:
The main result is a directapplication of Corollary 1 for k p = 10 − and k i = 10 . Indeed,using a dichotomy-kind algorithm, Table IV is obtained. Itshows the estimated decay-rate µ at a given order between and .The first thing to note in Table IV is the hierarchy property,the decay-rate is an increasing function of the degree, as notedin Remark 2. Note also that the gap between order andorder is significant, showing that using projections indeedimproves the results.For orders higher than , the estimated decay-rate increasesslightly, and, up to a four digits precision, it reaches itsmaximum value at N = 6 . Since it is then the maximumallowable decay-rate obtained using equation (23), that tendsto show that the Lyapunov functional used in this papertogether with condition (21) are sharp and provide a goodanalysis. Figure 5 represents a simulation on the linear systemand it confirms the same observation. Indeed, one can see thatthe energy of the system is well-bounded by an exponentialcurve and that the bound becomes more and more accurate as N increases. Numerical simulations are done using a first order approximation withat least space-discretization points and 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 . Fig. 5: Energy of X with the linear system for k p =10 − , k i = 10 and Ω = 5 . The initial condition is φ ( x ) = 4 (cid:0)(cid:82) x φ ∞ x ( s ) ds + 0 . x ) (cid:1) , φ = 2Ω and Z (0) = 2 [ z ∞ z ∞ ] (cid:62) . -3 Fig. 6: Values of gains k p and k i leading to a stable systemwith the maximum decay-rate for N = 5 using Theorem 1.The black area is stable and the white area is said unstable upto an order .
2) Stability of the closed loop system:
We are interestednow 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 µ max for an order N = 5 . That leads to Figure 6 where it is easy tosee that increasing the gain k p decreases the range of possible k i while it increases its speed (see equation (23)). It is quitenatural 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 isstable and that was the conclusion of the two papers [44],[45] using a different Lyapunov functional. Note that with theprevious papers, it was not possible to quantify the notion of“small enough gains k p and k i ” while it is possible to give anestimation with the method of this paper. B. The stick-slip effect on the nonlinear model
The previous subsection shows that the linear model isglobally asymptotically stable for some values of gains k p and k i , so the nonlinear system (25) is locally asymptoticallystable for a large desired angular velocity Ω . This can beenverified in Figure 7a for k p = 10 − , k i = 10 and Ω = 10 .We can see that the linear and the nonlinear systems behavesimilarly if their initial condition is close to the equilibrium. (a) Systems (10) and (25): Ω = φ = 10 , φ ( x ) =(1 + 0 .
32 sin( x )) (cid:82) x φ ∞ x ( s ) ds and Z (0) = (cid:0) z ∞ z ∞ (cid:1) (cid:62) . (b) System (25): Ω = 20 , φ = 0 , φ = 0 and Z (0) = 0 , . (c) System (25): Ω = φ = 5 , φ ( x ) =(1 + 0 . x )) (cid:82) x φ ∞ x ( s ) ds and Z (0) = (cid:0) z ∞ z ∞ (cid:1) (cid:62) . Fig. 7: Numerical simulations for systems (10) and (25) with k p = 10 − , k i = 10 .The higher Ω is, the larger the basin of attraction is. Indeed,the regulation tries to bring the system into the “quasi”-linearzone of T nl , where it is close to a constant T min as seenin Figure 3 and consequently, the stick-slip phenomenon mayoccur at the beginning but it is not effective for a long timeas depicted in Figure 7b which is a numerical simulation onthe nonlinear model with a zero initial condition.The real challenge is then the case of low desired angularvelocities Ω . The result of a simulation on the nonlinearmodel for Ω = 5 rad.s − ( k p = 10 − , k i = 10 ) is displayed Fig. 8: Solution X bound of BMI (27) for k p = 10 − , k i = 10 and Ω = 5 and some values of τ . The limit is X bound = 28 .in Figure 7c. First of all, note that the oscillations are witha frequency of . Hz and with an amplitude of roughly rad.s − . This is very close to what has been estimated usingFigure 2. Then, the model presented in Section II seems to bea valid approximation of the real behavior, at least concerningthe stick-slip phenomenon. C. Practical stability analysis
Now, one can evaluate the amplitude of the oscillationsusing Theorem 2 with k p = 10 − and k i = 10 . The resultfor several values of τ and for an order between and is depicted in Figure 8. Note that after order , there arenumerical errors in the optimization process and the resultis not accurate. The maximum τ is µ max = 0 . andthe higher τ is, the better is the optimization. It appears thatfor k i = 10 , k p = 10 − and Ω = 5 , the optimal X bound isaround .Figure 10 shows the energy of the system as a function oftime. One can see that the bound X bound is quite accuratesince the error between the maximum of the auto-oscillationsand X bound is around . Moreover, note that max | z − Ω | = 11 . . X bound , in other words, nearly half of theoscillations are concentrated in the variable z , which meansthe stick-slip mostly acts on the variable z and does not affectmuch the rest of the system. Particularly, it seems very difficultto estimate the variation of z knowing only φ t (0 , t ) .The final observation is about the variation of X bound fordifferent Ω . This is depicted in Figure 9. Up to errors in thenumerical optimization, it seems that there is not importantvariation of X bound when Ω increases. This is counter-intuitive and does not reflect the observations made withFigure 7. One explanation is that we didn’t state that T nl is astrictly decreasing function for positive θ . D. Design of a PI controller
Finally, the problem stated at the beginning of the paperwas to find the best PI controller, meaning that it minimizes X bound . The plot in Figure 11a shows that the value of theintegral gain k i does impact the oscillations due to stick-slipsince X bound increases from to for k i ∈ [0 . , and Fig. 9: Solution X bound of BMI (27) for k p = 10 − , k i = 10 , τ = 0 . and N = 5 . Fig. 10: Energy of X with the nonlinear system for k p = 10 − , k i = 10 and Ω = 5 . The initial condition is the same than inFigure 5.the LMIs become infeasible after this point (these values havebeen obtained with Ω = 5 , k p = 10 − and N = 5 ).To stay robustly stable against small delay in the loop, asnoted in Section IV-C, we should consider ≤ k p ≤ . · − which does not offer a large set of choices. For k i = 10 , we getFigure 11b. It appears that the minimum of X bound is obtainedfor k p close to . Consequently, increasing the gain k p seemsnot to reduce the stick-slip effect. Consequently, even if a PIcontroller does not weaken the stick-slip effect, the equilibriumpoint of the controlled system is practically stable. Moreover,it does enable an oscillation around the desired equilibriumpoint and a local convergence to that point.VII. C ONCLUSION
This paper focused on the analysis of the performance ofa PI controller for a drilling pipe. First, a discussion betweenexisting models in the literature allows us to conclude that theinfinite-dimensional model is closer to the real drilling pipeand should then be used for simulations and analysis. Basedon this last model, the exponential stability of the closed-loopsystem was ensured using a Lyapunov functional depending onprojections of the infinite-dimensional state on a polynomialsbasis. This approach enables getting exponential stability ofthe linear system with an estimation of the decay-rate. Thisresult was then extended to the nonlinear case considering (a) X bound for k p = 10 − but with different values of k i -3 (b) X bound for k i = 10 but with different values of k p Fig. 11: Minimal X bound obtained using Theorem 2 with Ω =5 and N = 5 .practical stability. The example section shows that it is possibleto get an estimate of the smallest attractive and invariant setand that this estimate is close to the minimal one. Furtherresearch would focus on similar analysis but using differentcontrollers as the ones developed in [13], [30], [41] or inChapter 10 of [38]. A PPENDIX
A. Legendre Polynomials and Bessel Inequality
The performance of the Lyapunov functional (16) highlydepends on the projection methodology developed in [42]. Tohelp the reader to better understand the proof of Theorem 1,a definition and some properties of Legendre polynomials isreminded. More information can be found in [17].
Definition 4:
The orthonormal family of Legendre polyno-mials {L k } k ∈ (cid:78) on L ([0 , embedded with the canonicalinner product is defined as follows: L k ( x ) = ( − k k (cid:88) l =0 ( − l kl k + ll x l , where (cid:0) kl (cid:1) = k ! l !( k − l )! .Their expended formulation is not useful in this paper butthe two following lemmas are of main interest. emma 4: For any function χ ∈ L and symmetric positivematrix R ∈ (cid:83) , the following Bessel-like integral inequalityholds for all N ∈ (cid:78) : (cid:90) χ (cid:62) ( x ) Rχ ( x ) dx (cid:62) N (cid:88) k =0 (2 k + 1) X (cid:62) k R X k (31)where X k is the projection coefficient of χ with respect to theLegendre polynomial L k as defined in (15). Lemma 5:
For any function χ ∈ L satisfying (15), thefollowing expression holds for any N in (cid:78) : (cid:34) ˙ X ... ˙ X N (cid:35) = (cid:49) N χ (1) − ¯ (cid:49) N χ (0) − L N (cid:34) X ... X N (cid:35) , (32)where L N , (cid:49) N and ¯ (cid:49) N are defined in (13). Proof.
This proof is highly inspired from [8]. Since χ ∈ L satisfies (15), equation (18) can be derived and the followingholds: ˙ X k = Λ [ χ ( x ) L k ( x )] − Λ (cid:90) χ ( x ) L (cid:48) ( x ) dx − γ t X k . As noted in [17], the interesting properties of Legendre poly-nomials are stated below:1) the boundary conditions on L k ensures L k (0) = ( − k and L k (1) = 1 ;2) the derivation rule for Legendre polynomials is ddx L k ( x ) = (cid:80) kj =0 (cid:96) j,k L j ( x ) .These two properties lead to the proposed result in equa-tion (32). R EFERENCES[1] U. J. F. Aarsnes and O. M. Aamo. Linear stability analysis of self-excitedvibrations in drilling using an infinite dimensional model.
Journal ofSound and Vibration , 360:239 – 259, 2016.[2] U.J.F. Aarsnes and J.S. Roman. Torsional vibrations with bit offbottom: Modeling, characterization and field data validation.
Journalof Petroleum Science and Engineering , 163:712 – 721, 2018.[3] B. Armstrong-Helouvry. Stick-slip arising from stribeck friction. In
Pro-ceedings., IEEE International Conference on Robotics and Automation ,pages 1377–1382. IEEE, 1990.[4] B. Armstrong-H´elouvry, P. Dupont, and C. Canudas De Wit. A surveyof models, analysis tools and compensation methods for the control ofmachines with friction.
Automatica , 30(7):1083–1138, 1994.[5] M. Barreau.
Stability analysis of coupled ordinary differential systemswith a string equation - Application to a Drilling Mechanism . Universit´eF´ed´erale Toulouse Midi-Pyr´en´ees, 2019.[6] M. Barreau, F. Gouaisbaut, A. Seuret, and R. Sipahi. Input/outputstability of a damped string equation coupled with ordinary differen-tial system.
International Journal of Robust and Nonlinear Control ,28(18):6053–6069, 2018.[7] M. Barreau, A. Seuret, and F. Gouaisbaut. Exponential Lyapunovstability analysis of a drilling mechanism. In , pages 2952–2957, 2018.[8] M. Barreau, A. Seuret, F. Gouaisbaut, and L. Baudouin. Lyapunov sta-bility analysis of a string equation coupled with an ordinary differentialsystem.
IEEE Transactions on Automatic Control , 63(11):3850–3857,Nov 2018.[9] G. Bastin and J.-M. Coron.
Stability and boundary stabilization of 1-dhyperbolic systems , volume 88. Springer, 2016.[10] H.I. Basturk. Observer-based boundary control design for the sup-pression of stick–slip oscillations in drilling systems with only surfacemeasurements.
Journal of Dynamic Systems, Measurement, and Control ,139(10):104501, 2017. [11] A. Bisoffi, M. Da Lio, A. R. Teel, and L. Zaccarian. Global asymp-totic stability of a PID control system with coulomb friction.
IEEETransactions on Automatic Control , 63(8):2654–2661, 2017.[12] D. Bresch-Pietri and M. Krstic. Output-feedback adaptive control of awave PDE with boundary anti-damping.
Automatica , 50(5):1407–1415,2014.[13] C. Canudas de Wit, F. Rubio, and M. Corchero. DOSKIL: A New Mech-anism for Controlling Stick-Slip Oscillations in Oil Well Drillstrings.
IEEE Transactions on Control Systems Technology , 16(6):1177–1191,November 2008.[14] N. Challamel. Rock destruction effect on the stability of a drillingstructure.
Journal of Sound and Vibration , 233(2):235–254, 2000.[15] A. P. Christoforou and A. S. Yigit. Fully coupled vibrations of activelycontrolled drillstrings.
Journal of sound and vibration , 267(5):1029–1045, 2003.[16] J. M. Coron.
Control and nonlinearity . Number 136 in MathematicalSurveys and Monographs. American Mathematical Soc., 2007.[17] R. Courant and D. Hilbert.
Methods of mathematical physics . JohnWiley & Sons, Inc., 1989.[18] Y. Ebihara, D. Peaucelle, and D. Arzelier.
S-Variable Approach toLMI-Based Robust Control , volume 17 of
Communications and ControlEngineering . Springer, 2015.[19] A. F. Filippov. Classical solutions of differential equations with multi-valued right-hand side.
SIAM Journal on Control , 5(4):609–621, 1967.[20] E. Fridman, S. Mondi´e, and B. Saldivar. Bounds on the responseof a drilling pipe model.
IMA Journal of Mathematical Control andInformation , 27(4):513–526, 2010.[21] J. K. Hale and S. M. V. Lunel. Effects of small delays on stabilityand control. In
Operator theory and analysis , pages 275–301. Springer,2001.[22] A. Helmicki, C. A. Jacobson, and C. N. Nett. Ill-posed distributedparameter systems: A control viewpoint.
IEEE Trans. on AutomaticControl , 36(9):1053–1057, 1991.[23] D. Karnopp. Computer simulation of stick-slip friction in mechanicaldynamic systems.
Journal of dynamic systems, measurement, andcontrol , 107(1):100–103, 1985.[24] H.K. Khalil.
Nonlinear Systems . Pearson Education. Prentice Hall, 1996.[25] M. Krstic.
Delay compensation for nonlinear, adaptive, and PDEsystems . Springer, 2009.[26] N. Levinson. Transformation theory of non-linear differential equationsof the second order.
Annals of Mathematics , pages 723–737, 1944.[27] X. Liu, N. Vlajic, X. Long, G. Meng, and B. Balachandran. Coupledaxial-torsional dynamics in rotary drilling with state-dependent delay:stability and control.
Nonlinear Dynamics , 78(3):1891–1906, Nov 2014.[28] J. L¨ofberg. YALMIP: A toolbox for modeling and optimization inMATLAB. In
IEEE International Symposium on Computer AidedControl Systems Design , pages 284–289, 2005.[29] ¨O. Morg¨ul. On the stabilization and stability robustness against smalldelays of some damped wave equations.
IEEE Trans. on AutomaticControl , 40(9):1626–1630, 1995.[30] E. Navarro-L´opez. An alternative characterization of bit-sticking phe-nomena in a multi-degree-of-freedom controlled drillstring.
NonlinearAnalysis: Real World Applications , 10:3162–3174, 2009.[31] E. Navarro-L´opez and D. Cortes. Sliding-mode control of a multi-dofoilwell drillstring with stick-slip oscillations. In
Proceedings of theAmerican Control Conference , pages 3837 – 3842, 2007.[32] E. Navarro-L´opez and R. Suarez. Practical approach to modelling andcontrolling stick-slip oscillations in oilwell drillstrings. In
Proceedingsof the 2004 IEEE International Conference on Control Applications,2004. , volume 2, pages 1454–1460 Vol.2, 2004.[33] C. Prieur, S. Tarbouriech, and J. M. G. da Silva. Wave equationwith cone-bounded control laws.
IEEE Trans. on Automatic Control ,61(11):3452–3463, 2016.[34] T. Richard, C. Germay, and E. Detournay. A simplified model to explorethe root cause of stick-slip vibrations in drilling systems with drag bits.
Journal of Sound and Vibration , 305(3):432–456, 8 2007.[35] C. Roman, D. Bresch-Pietri, E. Cerpa, C. Prieur, and O. Sename. Back-stepping observer based-control for an anti-damped boundary wave PDEin presence of in-domain viscous damping. In , 2016.[36] M. Safi, L. Baudouin, and A. Seuret. Tractable sufficient stability con-ditions for a system coupling linear transport and differential equations.
Systems and Control Letters , 110:1–8, December 2017.[37] C. Sagert, F. Di Meglio, M. Krstic, and P. Rouchon. Backstepping andflatness approaches for stabilization of the stick-slip phenomenon fordrilling.
IFAC Proceedings Volumes , 46(2):779–784, 2013.38] B. Saldivar, I. Boussaada, H. Mounier, and S.-I. Niculescu.
Analysis andControl of Oilwell Drilling Vibrations: A Time-Delay Systems Approach .Springer, 2015.[39] B. Saldivar, S. Mondi´e, and J. C. ´Avila Vilchis. The control of drillingvibrations: A coupled PDE-ODE modeling approach.
InternationalJournal of Applied Mathematics and Computer Science , 2016.[40] B. Saldivar, S. Mondie, S.-I. Niculescu, H. Mounier, and I. Boussaada.A control oriented guided tour in oilwell drilling vibration modeling.
Annual Reviews in Control , 42:100–113, September 2016.[41] A. Serrarens, M.J.G. Molengraft, J.J. Kok, and L. van den Steen. H ∞ control for suppressing stick-slip in oil well drillstrings. IEEE ControlSystems , 18:19 – 30, 05 1998.[42] A. Seuret and F. Gouaisbaut. Hierarchy of LMI conditions for thestability analysis of time-delay systems.
Systems & Control Letters ,81:1–7, 2015.[43] S. Tarbouriech, G. Garcia, J. M. G. da Silva Jr, and I. Queinnec.
Stabilityand stabilization of linear systems with saturating actuators . SpringerScience & Business Media, 2011.[44] A. Terrand-Jeanne, V. Andrieu, M. Tayakout-Fayolle, and V. Dos San-tos Martins. Regulation of inhomogeneous drilling model with a P-Icontroller.
IEEE Transactions on Automatic Control , pages 1–1, 2020.[45] A. Terrand-Jeanne, V. Dos Santo Martins, and V. Andrieu. Regulation ofthe downside angular velocity of a drilling string with a P-I controller.
ECC 2018, Cyprus , pages 2647–2652, 2018.[46] K.-C. Toh, M. J. Todd, and R. H. T¨ut¨unc¨u. SDPT3 - a MATLABsoftware package for semidefinite programming.
Optimization methodsand software , 11(1-4):545–581, 1999.[47] W. R. Tucker and C. Wang. On the effective control of torsional vibra-tions in drilling systems.
Journal of Sound and Vibration , 224(1):101–122, 1999.[48] M. Tucsnak and G. Weiss.
Observation and control for operatorsemigroups . Springer, 2009.[49] W. Jr. Weaver, S. P. Timoshenko, and D. H. Young.
Vibration problemsin engineering . John Wiley & Sons, 1990.
M. Barreau got the Engineers degree in aeronauticalengineering from ISAE-ENSICA (Toulouse, France)and a double degree from KTH in Space Engineeringwith specialization in system engineering in 2015.He was a Ph.D. student in LAAS-CNRS ‘Labora-toire d’Analyse et d’Architecture des Systmes” inToulouse, France under the supervision of A. Seuretand F. Gouaisbaut. He received the “Diplme de Doc-torat” in System Theory in 2019 from Universit PaulSabatier, Toulouse. His research interests includestability analysis of infinite dimension systems witha particular focus on drilling systems.
F. Gouaisbaut was born in Rennes (France) inApril 26,1973. He received the Diplme dIngnieurin automatic control (Engineers degree) from EcoleCentrale de Lille, France, in September 1997 andthe Diplme dEtudes Approfondies (Masters Degree)in the same subject from the University of Scienceand Technology of Lille, France, in September 1997.From October 1998 to October 2001 he was aPh.D. student at the Laboratoire dAutomatique, GnieInformatique et Signal (LAGIS) in Lille, France. Hereceived the Diplme de Doctorat (Ph.D. degree) fromEcole Centrale de Lille and University of Science and Technology of Lille,France, in October 2001. Since October 2003, he is an associate professor atthe Paul Sabatier University (Toulouse). His research interests include timedelay systems, quantized systems and robust control.