# Design of Telerobotic Drilling Control System with Haptic Feedback

Authors: Faraz Shah; Ilia G. Polushin
Journal: Journal of Control Science and Engineering
 (2013)

---

## Abstract

The paper deals with the design of control algorithms for virtual reality based telerobotic
system with haptic feedback that allows for the remote control of the vertical drilling operation. The human operator controls the vertical penetration velocity using a haptic device while simultaneously receiving the haptic feedback from the locally implemented virtual environment. The virtual environment is rendered as a virtual spring with stiffness updated based on the estimate of the stiffness of the rock currently being cut. Based on the existing mathematical models of drill string/drive systems and rock cutting/penetration process, a robust servo controller is designed which guarantees the tracking of the reference vertical penetration velocity of the drill bit. A scheme for on-line estimation of the rock intrinsic specific energy is implemented. Simulations of the proposed control and parameter estimation algorithms have been conducted; consequently, the overall telerobotic drilling system with a human operator controlling the process using PHANTOM Omni haptic device is tested experimentally, where the drilling process is simulated in real time in virtual environment.


---

## Body

## 1. Introduction

Drilling a borehole is a common method for extracting oil, gas, and natural resources from beneath the surface of the earth. Conventional oil well drilling has made significant progress over recent years, and currently is one of the most automated processes in the oil and gas industry. However, there are still some fundamental challenges associated with the drilling. One of the challenges is the choice of vertical penetration velocity of the drill bit. For efficient drilling operation, this velocity must depend upon the type of rock beds drilled. In particular, the velocity must be adjusted when mechanical characteristics of rock strata change. Often, it is difficult to estimate in real time the relative position of the drill bit with respect to different rock layers and, therefore, hard to predict the mechanical characteristics of the rock formations.The goal of this research is to design a telerobotic system with haptic feedback for control of the drilling process. Telerobotics for drilling well is a relatively novel idea, and it is substantial endeavor to automate one of the fundamental processes in the extraction of energy and resources. As telerobotics is integrated with drilling, it can greatly decrease the number of people working and monitoring operation on the site. This, in particular, can reduce the work site hazards. Also, telerobotics can bring actual analysis of in situ conditions (underground drilling environment) in real time to the human operator that works remotely, where (s)he will be able to monitor the current drilling conditions and, in particular, promptly enforce changes in the vertical speed of penetration of the drill bit in the oil well. Real-time control and optimization of the drilling speed are crucial for today's drilling industry, as it can reduce time and immense cost associated with the drilling an oil well. Introduction of haptic feedback would allow the human operator to feel the changes in mechanical characteristics of the rock and adjust the vertical velocity of penetration accordingly.In this paper, we address the problem of design of control algorithms for virtual reality based telerobotic system with haptic feedback that allows for the remote control of the vertical drilling operation. Based on a simplified mathematical model of the drilling process, control algorithms are designed which allow to achieve a desired rate of the vertical penetration, regardless of the mechanical properties of the rock. The control design includes an online parameter estimator of the intrinsic specific energy which is a parameter that describes the hardness of the rock. All these algorithms are consequently used in the design of a telerobotic drilling system with virtual environment-based haptic feedback that allows the human operator to feel the stiffness of the rock in contact with the drill bit. Simulations and semiexperimental results are performed which confirm the validity of the theoretical developments.The potential application domain of this research is not limited to onshore/offshore oil well drilling, but the same principles can be applied, in particular, to different types of mining robots [1], telerobotic systems for dredging and mining ocean [2–5], surgical drilling [6], and telerobotic systems for drilling the extraterrestrial terrain to discover and research the minerals and composition beneath [5, 7].The structure of the paper is as follows. In Section2, a mathematical model of the drilling process is derived which is subsequently used for the control design. Section 3 deals with the design of control algorithms for rotational and translational motion of the drilling systems, as well as the design of an online parameter estimator of the intrinsic specific energy of the rock. In Section 4, the structure of a telerobotic drilling system is described and the corresponding experimental results are presented. Finally, in Section 5, some conclusions are given and possible future directions are formulated.


## 2. Mathematical Model of Drilling System

In this section, mathematical models that describe the drilling system are presented. Specifically, the mathematical model of drill string and drive system is described in Section2.1, while the model of rock cutting and penetration is the subject of Section 2.2.


### 2.1. Mathematical Model of the Drill String and Drive System

The drill string is the assembly of rotating pipes which are responsible for transmitting rotation and weight to the bit and bridge up a connection between the bottom hole tools [8]. The components of a drill string along with drill pipes and the bottom hole assembly (BHA) are shown in Figure 1. A number of simplified mathematical models for drill string and drive systems were proposed in the literature, such as [9–12]. The model used in our work was developed in [9]. This model describes the drill string as a simple torsional pendulum, where the drill pipes are represented as torsional springs and the bottom hole assembly is described as a rigid body with inertia. The model is based on the following simplifying assumptions.(1)
The bottom hole assembly and the drill bits behave like rigid bodies.(2)
The moment of inertia of the drill pipe is considered to be small in comparison with the moments of inertia of the bottom hole assembly and the rotary table and, therefore, neglected.(3)
The nonzero time propagation of the torsional force disturbances along the drill string is neglected. The forces assume to propagate instantaneously along the drill string.Figure 1
Drill string components [8].Under the above described assumptions 1–3, the whole drill string and drive system with equivalent electro-mechanical components can be represented by its structural diagram shown in Figure2. This system is described by the following mathematical model [9]. First, the motion of the drill string is described by the following equation:
(1)J1ϕ¨1+c1ϕ˙1+k(ϕ1-ϕ2)-T=0.
Here, ϕ1 is the angular displacement of bit and drill collars (BHA), ϕ2 is the angular displacement of the rotary table, J1 is the equivalent moment of inertia of the collars (BHA) and the drill pipes, coefficient c1 represents equivalent viscous damping, k is the equivalent torsional stiffness of the drill pipes, and T is the torque-on-bit (TOB) generated during the rock cutting process (see Section 2.2 below). The dynamics of the rotary table and drive system is described by the following equation:
(2)J2ϕ¨2+c2ϕ˙2-k(ϕ1-ϕ2)-nTm=0,
where J2 is combined moment of inertia of the rotary table and of the rotor of the electric motor coupled together with a gearbox that has 1:n gear ratio, c2 is aggregated damping of all the components of the drive system, and Tm is the motor torque. Finally, the electric motor is described by the following equations:
(3)LI˙+RI+Vb-V=0,Vb=Kϕ˙3=Knϕ˙2,Tm=KI,
where I is the armature current, L is an equivalent armature inductance, R is an equivalent armature resistance, Vb is the back emf, V is the armature voltage, ϕ˙3 is the rotor angular velocity, and K is a constant that depends upon the motor characteristics.Figure 2
Representation of drill string/drive system with mechanical and electrical components [9].By combining all the above equations, the complete drill string/drive system can be written in the following state space form:(4)[ϕ˙1ω˙1ϕ˙2ω˙2I˙]=[01000-kJ1-c1J1kJ10000010kJ20-kJ2-c2J2KnJ2000-KnL-RL][ϕ1ω1ϕ2ω2I]+[0-TJ100VL].
Here, ω1 and ω2 are the angular velocities of the drill bit and the rotary table, respectively. Equation (4) is valid when the drill bit rotational velocity is greater than zero, that is, ω1>0. In order to reduce the number of equations, a variable ϕ is introduced as the difference of ϕ2 and ϕ1. In this case, the original system can be rewritten in the following reduced state space form:
(5)[ω˙1ϕ˙ω˙2I˙]=[-c1J1kJ100-10100-kJ2-c2J2KnJ200-KnL-RL][ω1ϕω2I]+[-TJ100VL].
Equation (5) defines the reduced order model of the drill string and drive system. The model (5) is used for the control design below.


### 2.2. Rock Cutting and Vertical Penetration Models

A standard drill bit usually exhibits two kinds of motions: rotational along its axis of rotation and vertical motion while penetrating through the rocks. As described in  [13], in the normal mode of operation of the drill bit, the bit rotational velocity ω is parallel to its axis of rotation, and the penetration velocity v is directed vertically straight through the rocks. Similarly, the weight-on-bit W acts in the vertical direction and the torque-on-bit T is applied in parallel to the direction of rotation of drill bit. The cutting components of the weight-on-bit and torque-on-bit depend on the radius of PDC drill bit a, intrinsic specific energy ϵ, a parameter ζ>0 which represents the ratio of the vertical force to the horizontal force between rock and cutter contact surfaces, and the depth of cut d. The depth of cut d plays significant role in the equations to follow that describe the cutting components of the torque-on-bit T and the weight-on-bit W. The equations for these two cutting components are as follows [13]:
(6)Tc=12a2ϵd,(7)Wc=aζϵd.In this work, the system is developed under simplifying assumption that the friction effects are negligible. In this case, both variablesT≈Tc and W≈Wc are proportional to the depth of cut d, according to (6) and (7). As illustrated in Figure 3, the depth of cut d is the thickness of rock ridge in front of the blade. It is assumed that the drill bit has n number of identical blades, and the difference of angular positions of these two successive blades is (2π/n). In this case, d is the combined depth of cut of all n blades in each revolution of drill bit, according to the formula
(8)d(t):=ndn(t),
where dn is the depth of cut of each blade. The depth of cut for each blade is in turn defined according to the formula
(9)dn(t):=U(t)-U(t-tn),
where U(t) and U(t-tn) are the vertical positions of the drill bit at current time instant t and a certain previous instant t-tn, respectively [10, 11]. The delay tn in the above formula is exactly the time that is required for the drill bit to rotate by an angle 2π/n to achieve its current angular position ϕ1(t); in other words, it also satisfies the following equation:
(10)ϕ(t)-ϕ(t-tn)=2πn.Figure 3
Section of the bottom hole profile located between two successive blades [10].Using (9) and (10) for calculating d(t) would significantly complicate the control design. In this work, we simplify this problem by assuming that both the vertical and angular velocities change slowly; specifically, it is assumed that both v(τ)≡U˙(τ) and ω1(τ)≡ϕ˙1(τ) are approximately constant during each period τ∈[t-tn,t]. Using this assumptions, (9) and (10) can be rewritten as follows:
(11)d(t)≈n·v(t)·tn,(12)ω1(t)·tn≈2πn.
Combining (11), (12), and assuming ω1(t)≠0, one gets the following approximate expression for d(t):
(13)d(t)≈2π·t(t)ω1(t).The above formula has a singularity atω1(t)=0. To remove this singularity, note that the drilling occurs when both ω1(t)>0 and v(t)>0. On the contrary, ω1(t)≤0, the drill bits do not cut the rock and therefore d(t)≡0 in this case. Based on the above considerations, one can approximately define the depth of cut according to the formula
(14)d(t)≈2π·v(t)max{ω1(t),ϵ0},
where ϵ0>0 is sufficiently small positive constant. The formula (14) does not have singularity at ω1(t)=0; it will be occasionally used for calculations of d(t) instead of (13) in the cases where avoiding singularity is important (in simulations, etc.).Finally, the vertical motion of the drill bit is described by the following equation [12]:
(15)Mdvdt=Ws-W-H0-Kfv.
Here, v is the vertical velocity of the drill bit, M is the combined mass of the drill string and BHA, H0 is the constant upward force applied from the top of drilling rig, and Ws is the submerged weight of the drill string and Bottom Hole Assembly (BHA). In this model, it is assumed that Ws and H0 to be constants and defined their difference with another constant W0 such that W0=Ws-H0. Also, W is the applied weight on bit from the interaction of rock defined by (7), and Kf>0 is the coefficient of viscous friction.


## 2.1. Mathematical Model of the Drill String and Drive System

The drill string is the assembly of rotating pipes which are responsible for transmitting rotation and weight to the bit and bridge up a connection between the bottom hole tools [8]. The components of a drill string along with drill pipes and the bottom hole assembly (BHA) are shown in Figure 1. A number of simplified mathematical models for drill string and drive systems were proposed in the literature, such as [9–12]. The model used in our work was developed in [9]. This model describes the drill string as a simple torsional pendulum, where the drill pipes are represented as torsional springs and the bottom hole assembly is described as a rigid body with inertia. The model is based on the following simplifying assumptions.(1)
The bottom hole assembly and the drill bits behave like rigid bodies.(2)
The moment of inertia of the drill pipe is considered to be small in comparison with the moments of inertia of the bottom hole assembly and the rotary table and, therefore, neglected.(3)
The nonzero time propagation of the torsional force disturbances along the drill string is neglected. The forces assume to propagate instantaneously along the drill string.Figure 1
Drill string components [8].Under the above described assumptions 1–3, the whole drill string and drive system with equivalent electro-mechanical components can be represented by its structural diagram shown in Figure2. This system is described by the following mathematical model [9]. First, the motion of the drill string is described by the following equation:
(1)J1ϕ¨1+c1ϕ˙1+k(ϕ1-ϕ2)-T=0.
Here, ϕ1 is the angular displacement of bit and drill collars (BHA), ϕ2 is the angular displacement of the rotary table, J1 is the equivalent moment of inertia of the collars (BHA) and the drill pipes, coefficient c1 represents equivalent viscous damping, k is the equivalent torsional stiffness of the drill pipes, and T is the torque-on-bit (TOB) generated during the rock cutting process (see Section 2.2 below). The dynamics of the rotary table and drive system is described by the following equation:
(2)J2ϕ¨2+c2ϕ˙2-k(ϕ1-ϕ2)-nTm=0,
where J2 is combined moment of inertia of the rotary table and of the rotor of the electric motor coupled together with a gearbox that has 1:n gear ratio, c2 is aggregated damping of all the components of the drive system, and Tm is the motor torque. Finally, the electric motor is described by the following equations:
(3)LI˙+RI+Vb-V=0,Vb=Kϕ˙3=Knϕ˙2,Tm=KI,
where I is the armature current, L is an equivalent armature inductance, R is an equivalent armature resistance, Vb is the back emf, V is the armature voltage, ϕ˙3 is the rotor angular velocity, and K is a constant that depends upon the motor characteristics.Figure 2
Representation of drill string/drive system with mechanical and electrical components [9].By combining all the above equations, the complete drill string/drive system can be written in the following state space form:(4)[ϕ˙1ω˙1ϕ˙2ω˙2I˙]=[01000-kJ1-c1J1kJ10000010kJ20-kJ2-c2J2KnJ2000-KnL-RL][ϕ1ω1ϕ2ω2I]+[0-TJ100VL].
Here, ω1 and ω2 are the angular velocities of the drill bit and the rotary table, respectively. Equation (4) is valid when the drill bit rotational velocity is greater than zero, that is, ω1>0. In order to reduce the number of equations, a variable ϕ is introduced as the difference of ϕ2 and ϕ1. In this case, the original system can be rewritten in the following reduced state space form:
(5)[ω˙1ϕ˙ω˙2I˙]=[-c1J1kJ100-10100-kJ2-c2J2KnJ200-KnL-RL][ω1ϕω2I]+[-TJ100VL].
Equation (5) defines the reduced order model of the drill string and drive system. The model (5) is used for the control design below.


## 2.2. Rock Cutting and Vertical Penetration Models

A standard drill bit usually exhibits two kinds of motions: rotational along its axis of rotation and vertical motion while penetrating through the rocks. As described in  [13], in the normal mode of operation of the drill bit, the bit rotational velocity ω is parallel to its axis of rotation, and the penetration velocity v is directed vertically straight through the rocks. Similarly, the weight-on-bit W acts in the vertical direction and the torque-on-bit T is applied in parallel to the direction of rotation of drill bit. The cutting components of the weight-on-bit and torque-on-bit depend on the radius of PDC drill bit a, intrinsic specific energy ϵ, a parameter ζ>0 which represents the ratio of the vertical force to the horizontal force between rock and cutter contact surfaces, and the depth of cut d. The depth of cut d plays significant role in the equations to follow that describe the cutting components of the torque-on-bit T and the weight-on-bit W. The equations for these two cutting components are as follows [13]:
(6)Tc=12a2ϵd,(7)Wc=aζϵd.In this work, the system is developed under simplifying assumption that the friction effects are negligible. In this case, both variablesT≈Tc and W≈Wc are proportional to the depth of cut d, according to (6) and (7). As illustrated in Figure 3, the depth of cut d is the thickness of rock ridge in front of the blade. It is assumed that the drill bit has n number of identical blades, and the difference of angular positions of these two successive blades is (2π/n). In this case, d is the combined depth of cut of all n blades in each revolution of drill bit, according to the formula
(8)d(t):=ndn(t),
where dn is the depth of cut of each blade. The depth of cut for each blade is in turn defined according to the formula
(9)dn(t):=U(t)-U(t-tn),
where U(t) and U(t-tn) are the vertical positions of the drill bit at current time instant t and a certain previous instant t-tn, respectively [10, 11]. The delay tn in the above formula is exactly the time that is required for the drill bit to rotate by an angle 2π/n to achieve its current angular position ϕ1(t); in other words, it also satisfies the following equation:
(10)ϕ(t)-ϕ(t-tn)=2πn.Figure 3
Section of the bottom hole profile located between two successive blades [10].Using (9) and (10) for calculating d(t) would significantly complicate the control design. In this work, we simplify this problem by assuming that both the vertical and angular velocities change slowly; specifically, it is assumed that both v(τ)≡U˙(τ) and ω1(τ)≡ϕ˙1(τ) are approximately constant during each period τ∈[t-tn,t]. Using this assumptions, (9) and (10) can be rewritten as follows:
(11)d(t)≈n·v(t)·tn,(12)ω1(t)·tn≈2πn.
Combining (11), (12), and assuming ω1(t)≠0, one gets the following approximate expression for d(t):
(13)d(t)≈2π·t(t)ω1(t).The above formula has a singularity atω1(t)=0. To remove this singularity, note that the drilling occurs when both ω1(t)>0 and v(t)>0. On the contrary, ω1(t)≤0, the drill bits do not cut the rock and therefore d(t)≡0 in this case. Based on the above considerations, one can approximately define the depth of cut according to the formula
(14)d(t)≈2π·v(t)max{ω1(t),ϵ0},
where ϵ0>0 is sufficiently small positive constant. The formula (14) does not have singularity at ω1(t)=0; it will be occasionally used for calculations of d(t) instead of (13) in the cases where avoiding singularity is important (in simulations, etc.).Finally, the vertical motion of the drill bit is described by the following equation [12]:
(15)Mdvdt=Ws-W-H0-Kfv.
Here, v is the vertical velocity of the drill bit, M is the combined mass of the drill string and BHA, H0 is the constant upward force applied from the top of drilling rig, and Ws is the submerged weight of the drill string and Bottom Hole Assembly (BHA). In this model, it is assumed that Ws and H0 to be constants and defined their difference with another constant W0 such that W0=Ws-H0. Also, W is the applied weight on bit from the interaction of rock defined by (7), and Kf>0 is the coefficient of viscous friction.


## 3. Controller Design

The block diagram of the overall drilling system is shown in Figure4. As it can be seen from this figure, the block diagram has a complex structure and consists of several interconnected subsystems. Specifically, the vertical motion subsystem is described by (15); the output of this subsystem is the vertical velocity of penetration v(t). The subsystem that represents the rotational motion is described by (5); this subsystem has one control input which is the armature voltage V(t) and one output which is the angular velocity of the drill bits ω1(t). Both v(t) and ω1(t) are the inputs of the nonlinear static block that represents the cutting process; this subsystem generates the depth of cut d(t) according to (13). Both the torque-on-bit T and weight-on-bit W are proportional to d; they are fed back to rotational motion and vertical motion subsystems, respectively.Figure 4
The block diagram of the drilling system.Our goal is to design a control system that maintains a desired rate of drilling. Specifically, we are looking for the control algorithm for the armature voltageV that would guarantee that the velocity of the vertical penetration v(t) tends asymptotically to an arbitrary positive desired value vref>0. We start designing a control algorithm by considering the equation of vertical motion (15) in some detail.


### 3.1. Control of the Vertical Motion of a Drill Bit

The vertical motion of the drilling system is described by (15). For convenience, this equation is rewritten below in a slightly modified form, as follows:
(16)v˙=-KfMv-(Ws-H0)M-WM.
The idea of the controller developed in this work is to use the weight-on-bit W as the control input to the vertical motion subsystem (16). More specifically, combining formulas (7) and (13), one get the following expression for W:
(17)W=aζϵ2πω1v,
which essentially indicates that W is proportional to the vertical velocity v(t) and inversely proportional to the angular velocity of the rotational motion ω1(t). Substituting the last formula into (16), one gets
(18)v˙=Ws-H0M-1M(aζϵ2πω1+Kf)v.
Equation (18) is a linear differential equation with respect to v which, assuming ω1>0, has one stable equilibrium v=v0 defined by the formula
(19)Ws-H0M-1M(aζϵ2πω1+Kf)v0=0.
Solving the above equation with respect to v0, one gets
(20)v0=Ws-H0(aζϵ(2π/ω1)+Kf).
The above equation (20) indicates that the location of the stable equilibrium v=v0 of the vertical motion subsystem (16) can be controlled if one can control the rotational velocity ω1. Specifically, (20) defines one-to-one correspondence between ω1 from the range (0,+∞) and v0 from the range (0,(Ws-H0)/Kf). In particular, for any given vref∈(0,(Ws-H0)/Kf), there exists an unique ωref∈(0,+∞) such that if the angular velocity satisfies ω1(t)≡ωref, then vref is a globally exponentially stable equilibrium of the translational dynamics (16). For a given vref∈(0,(Ws-H0)/Kf), the corresponding ωref can be found using formula (20), as follows:
(21)ωref=2πaζϵ((Ws-H0)/vref)-Kf.
Therefore, the control goal of stabilization of the vertical penetration velocity v(t)→vref can be achieved by designing a controller for rotational motion that guarantees a sufficiently fast convergence of ω1(t)→ωref. The design of such a controlled is presented in the next section.


### 3.2. Stabilization of the Angular Velocity of the Drilling System

The rotational dynamics of the drilling system together with the electric drive are described by (5), which is repeated below for convenience,
(22)[ω˙1ϕ˙ω˙2I˙]=[-c1J1kJ100-10100-kJ2-c2J2KnJ200-KnL-RL][ω1ϕω2I]+[0001L]V+[-1J1000]T.
The above system has one control input which is the armature voltage of the electric drive V and one disturbance input which is the torque-on-bit T. Our objective in this section is to design a control law for V which would track the reference angular velocity of the drill ω1→ωref while rejecting the disturbance T.To solve the control problem formulated above, one can use the approach to feedforward robust servo control problem presented in [14, 15]. Below, the above approach is described in a simplified manner which, however, serves our purpose well. Consider a linear time invariant system of the form
(23)x˙=Ax+Bu+Dw,y=Cx+Fu+Hw,
where x∈ℝn is the state, u∈ℝm is the control input, y∈ℝp is the output, w∈ℝr are the disturbances, and A, B, C, D, F, and H are matrices of appropriate dimensions. Consider a control problem described as follows. Suppose the disturbances w(t) are measurable. Given a desired output signal yref(t), design a control algorithm that guarantees y(t)→yref(t) as t→+∞. This problem was addressed in [14, 15] in a very general setting. In this work, a simple case is addressed where both yref and w(t) are assumed to be constant signals, yref(t)≡yref and w(t)≡wm. In this case, the following two conditions are necessary and sufficient for the existence of a linear time-invariant controller that solves the above described problem.(i)
The pair(A,B) is stabilizable, which means that
(24)rank[B,AB,A2B,…,An-1B]=n;(ii)
Consider(25)rank[ABCF]=n+p.If the above two conditions hold (and only in this case), the linear time-invariant controller that solves the above described problem is given according to the formula(26)u=Kx+𝒢†yref+𝒢*wm,
where K∈ℝn×n is the feedback gain matrix which is to be chosen such that A-BK is stable and has the required dynamic properties
(27)𝒢=-C(A-BK)-1B,(28)𝒢*=𝒢†C(A-BK)-1D,
where 𝒢† is the Moore-Penrose pseudoinverse of the matrix 𝒢 in (27), defined by the formula
(29)𝒢†=𝒢T(𝒢𝒢T)-1.The above described control approach can be applied to the problem of stabilization of the angular velocity of drilling as follows. Equations (22), which describe the rotational dynamics of a drilling system, can be rewritten in the form (23), where x:=[ω1ϕω2I]T∈ℝ4, u:=V∈ℝ1, y:=ω1∈ℝ1, w:=T∈ℝ1, and the corresponding matrices are
(30)A:=[-c1J1kJ100-10100-kJ2-c2J2KnJ200-KnL-RL],B=[0001L],D=[-1J1000],(31)C=[1000],F=[0],H=[0].
Below, we consider the drilling system with specific values of the parameters that are listed in Table 1. With these values, the matrices A, B, and D become
(32)A:=[-0.11231.264700-10100-0.2231-0.20050.020400-8640-2],B=[000200],D=[-0.0027000],
while C, F, and H are given by (31).Table 1
Numerical values for drilling system parameters.
Parameter
Description
Value
Unit
J
1
BHA + drill string inertia
374
[kgm2]
J
2
Rotary table + drive inertia
2120
[kgm2]
c
1
BHA damping
42
[
Nms
/
rad
]
c
2
Rotary table damping
425
[
Nms
/
rad
]
k
Drill string stiffness
473
[
Nm
/
rad
]
R
Motor armature resistance
0.010
[
Ω
]
L
Motor armature inductance
0.005
[H]
K
Motor constant
6
[Vs]
n
Combined gear ratio for bevel and gear box
7.2
—
a
Drill bit radius
0.108
[m]
ζ
Ratio of drilling strength to drilling specific energy
0.7
—
M
Mass of drill string (28120 Kg) + BHA (25080 Kg)
53000
[kg]
W
s
-
H
0
Submerged weightWs  − applied  weight  from top of the Rig H0
100 or 1000
[N]
K
f
Viscous friction coefficient
20
[
Nm
/
rad
]For the above system, the necessary and sufficient conditions for stabilization (24), (25) are satisfied. Indeed, the stabilizability condition (24) is satisfied since
(33)rank[B,AB,A2B,…An-1B]=rank[0005.15427300.0000004.075472-8.96804.075-8.968-700.339200-400-34412146307]=4.
On the other hand, the rank condition (25) is also satisfied because
(34)rank[ABCF]=rank[-0.1122991.264706000-101000-0.223113-0.2004720.020377000-8640-220010000]=5.
Therefore, a controller of the form (26), (27), (28), and (29) guarantees that the angular velocity of the drill approach the reference angular velocity ω1→ωref as t→∞, while rejecting the disturbance Tb.The design of controller (26), (27), (28), and (29) begins by choosing the desired location of the closed-loop system's poles. For the purpose of simulations presented below, we consider two specific set of poles. The first set, denoted by P1, is chosen as follows:
(35)P1:=[-10-2+2i-2-2i-4].
The set P1 consists of two real poles and two complex conjugate poles. On the other hand, the set P2 contains only poles on the real axis as follows:
(36)P2=[-5.5-2-4.5-1].
The feedback gain matrix K1 such that the poles of A-BK1 are located according to P1 is
(37)K1=[32.2457.45-19.410.0784].
The coefficients 𝒢*, 𝒢† in (26) are calculated according to the formulas (27)–(29); the results are
(38)𝒢1*=0.123497,𝒢1†=60.0844.
On the other hand, the feedback matrix K2 such that the poles of A-BK2 are located according to P2 is
(39)K2=[-5.16716.943-30.620.0534].
The corresponding coefficients 𝒢2*, 𝒢2† are
(40)𝒢2*=0.037286,𝒢2†=9.603682.


### 3.3. Rock Stiffness Estimation

In the controller design presented above, it was assumed that the “hardness” of the rock, represented by the intrinsic specific energyϵ, is constant and exactly known. This knowledge of ϵ was used explicitly in the controller design, in particular, in formula (21). In practical geological drilling, however, the hardness of different layers of rock lying underneath the surface can be different and usually is not exactly known beforehand. More specifically, different characteristics of the rock, such as hardness, density and porosity, typically remain constant through each layer, but differs from layer to layer. On the other hand, control engineers frequently deal with the problem of designing a controller without a priori knowledge of the exact values of one or more parameters involved in the process. Often, the processes can be robustly controlled without the actual knowledge of some of the parameters. In other cases, the unknown parameters can be identified using specially designed estimators. Below, a simple online estimator of the rock intrinsic specific energy ϵ is designed following the methods described in [16], and the resulting estimate is then used in the controller for for drilling system.Specifically, during the cutting process, the torque-on-bitT is produced by bit rock interaction, according to the formula
(41)T=12a2ϵd,
where a is the radius of drill bit, d is the depth of cut, and ϵ>0 is the intrinsic specific energy. The intrinsic specific energy ϵ>0 depends on the properties of the media and typically unknown beforehand. However, since the torque on bit T(t) can typically be measured with advanced transducers located in the bottom hole assembly [17], a>0 is constant and known, and d(t) can be calculated according to formula (13), one can use the method described in the previous section to design an online estimation scheme for ϵ. In particular, considering (1/2)a2d(t) as the input and torque-on-bit T as the measured output, one can follow the procedure described in the previous section to design an estimator for an unknown parameter ϵ. The predicted torque-on-bit T^b is defined according to the formula
(42)T^(t):=12a2ϵ^d(t),
where ϵ^(t) is the current estimate of actual rock strength ϵ. The algorithm for online estimation of the intrinsic specific energy ϵ has a form
(43)ϵ^˙=γ0(T-T^)12a2d,
where γ0>0 is an arbitrary gain.A natural question regarding the algorithm (43) is if it guarantees the convergence of the parameter estimate to the true value of the parameter ϵ; mathematically, is ϵ^(t)→ϵ as t→+∞. It is known [16], that the convergence can be guaranteed if the “input” signal (1/2)a2d(t) is persistently exciting. A signal u(t) is said to be persistently exciting which is to say that there exist α0>0, T0>0 such that the inequality
(44)‍∫tt+T0u2(τ)dτ≥α0T0
holds for all t. In particular, u(t) is persistently exciting if u2(t)≥α0 for all t. Since d(t) is the depth of cut, we see that, during normal cutting process, d(t)≥d0>0, which results in persistent excitation of the input (1/2)a2d(t). The parameter convergence ϵ^(t)→ϵ, therefore, is guaranteed during normal cutting process. This is also confirmed by the simulation results presented below.The obtained estimate of the rock strengthϵ^ is then used in the control algorithm. Specifically, in the original formulation of the control algorithm, for a given reference vertical velocity vref, the reference rotational velocity ωref is calculated according to formula (21), which depends on the parameter ϵ. In case ϵ is unknown, it is substituted by its estimate ϵ^(t) obtained above. The new formula for ωref has the form
(45)ωref:=2πaζϵ^((Ws-H0)/vref)-Kf.
The obtained estimate of the rock stiffness ϵ^ will also be used to update the stiffness of the virtual spring in the haptic teleoperator drilling system described below.


### 3.4. Simulation Results

In this Section, some results of simulations of the drilling control system with intrinsic specific energy estimator are presented. The vertical motion of the drilling system is described by (15), and it is interconnected with the rotational dynamics (5) through nonlinear equation (13) that describes the depth of cut d(t). For a given reference velocity of the vertical penetration vref>0, the corresponding reference rotational velocity ωref is calculated according to formula (45). The controller (21), (26)–(29) has been implemented to guarantee that the angular velocity of the drill bits ω1(t) tracks ωref, which in turn stabilizes the vertical penetration velocity v(t) converges to vref. The algorithm (43) provides an estimate of the intrinsic specific energy parameter ϵ^ which is then used in the calculation of the reference angular velocity according to formula (45). Specific values of the parameters appearing in these equations are given in Table 1. The simulations are carried out using MATLAB, where the integration step for each simulation is equal to 0.005 s. The feedback gain matrix is chosen K=K1, where K1 is defined by (37).In the simulations described below, the performance of the system was evaluated for different values of actual intrinsic specific energyϵ, different gains γ0 and different values of the applied weight W0:=Ws-H0. Figures 5 and 6 show the response of the vertical penetration velocity v(t), the intrinsic specific energy estimate ϵ^(t), the torque-on-bit T(t), the predicted value of the torque-on-bit T^(t), and the rotational velocity ω1(t), all for the case where the applied weight on bit W0=5000 N, the intrinsic specific energy ϵ=20 MPa, and the desired vertical velocity vref is set to 0.005 m/s. The estimator gain is set to γ0=5·109. The plots show that v(t) converges to vref in less than 8 sec whereas the estimate ϵ^(t) converges to the actual value of ϵ in less than 4 sec. Figures 7 and 8 show the output responses of described parameters where W0=2500 N and the desired vertical velocity vref is set to 0.01 m/s. It can be clearly seen that the convergence becomes slower with reducing the applied weight on the drill string W0; specifically, both v(t) and ϵ^(t) approach their reference values in about 12 sec. Reducing W0 also results in that ω1ref increases, the steady-state value of Tb(t) drops to around 200 N, and the steady-state value of d(t) also drops to less than 2 mm. On the other hand, Figures 9 and 10 demonstrate the response of the system with the same parameters except the intrinsic specific energy ϵ is reduced to 5 MPa. This results in decreased convergence time for v(t) and ϵ^(t). The steady state value of rotational velocity ω1(t) is also decreased to under 10 rad/s, and steady state value of the depth of cut d(t) is increased to 6.5 mm. Figures 11 and 12 present the output response for the case where the estimator gain is decreased to γ0=5·108, while the rest of the parameters are the same as in the last simulation except the intrinsic specific energy is set to ϵ=10 MPa. The resulting response is predictably characterized by much slower convergence, which takes about 25 sec for v(t) and ϵ^(t) to approach their steady-state values. Finally, Figures 13 and 14 correspond to to the case where W0=5000 N, ϵ=20 MPa, and γ0=1·108.Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=5000 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) vresus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=5000 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=5 MPa, and γ0=5·109.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=5 MPa, and γ0=5·109.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=10 MPa, and γ0=5·108.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=10 MPa, and γ0=5·108.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=20 MPa, and γ0=108.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=20 MPa, and γ0=108.
(a)
(b)
(c)Overall, simulation results show that the control system with intrinsic specific energy estimation demonstrate good stability and performance characteristics for a wide range of the parameters. In particular, the vertical velocity converges to the desired value, and the estimate of the intrinsic specific energyϵ^(t) converges to an actual value of ϵ.


## 3.1. Control of the Vertical Motion of a Drill Bit

The vertical motion of the drilling system is described by (15). For convenience, this equation is rewritten below in a slightly modified form, as follows:
(16)v˙=-KfMv-(Ws-H0)M-WM.
The idea of the controller developed in this work is to use the weight-on-bit W as the control input to the vertical motion subsystem (16). More specifically, combining formulas (7) and (13), one get the following expression for W:
(17)W=aζϵ2πω1v,
which essentially indicates that W is proportional to the vertical velocity v(t) and inversely proportional to the angular velocity of the rotational motion ω1(t). Substituting the last formula into (16), one gets
(18)v˙=Ws-H0M-1M(aζϵ2πω1+Kf)v.
Equation (18) is a linear differential equation with respect to v which, assuming ω1>0, has one stable equilibrium v=v0 defined by the formula
(19)Ws-H0M-1M(aζϵ2πω1+Kf)v0=0.
Solving the above equation with respect to v0, one gets
(20)v0=Ws-H0(aζϵ(2π/ω1)+Kf).
The above equation (20) indicates that the location of the stable equilibrium v=v0 of the vertical motion subsystem (16) can be controlled if one can control the rotational velocity ω1. Specifically, (20) defines one-to-one correspondence between ω1 from the range (0,+∞) and v0 from the range (0,(Ws-H0)/Kf). In particular, for any given vref∈(0,(Ws-H0)/Kf), there exists an unique ωref∈(0,+∞) such that if the angular velocity satisfies ω1(t)≡ωref, then vref is a globally exponentially stable equilibrium of the translational dynamics (16). For a given vref∈(0,(Ws-H0)/Kf), the corresponding ωref can be found using formula (20), as follows:
(21)ωref=2πaζϵ((Ws-H0)/vref)-Kf.
Therefore, the control goal of stabilization of the vertical penetration velocity v(t)→vref can be achieved by designing a controller for rotational motion that guarantees a sufficiently fast convergence of ω1(t)→ωref. The design of such a controlled is presented in the next section.


## 3.2. Stabilization of the Angular Velocity of the Drilling System

The rotational dynamics of the drilling system together with the electric drive are described by (5), which is repeated below for convenience,
(22)[ω˙1ϕ˙ω˙2I˙]=[-c1J1kJ100-10100-kJ2-c2J2KnJ200-KnL-RL][ω1ϕω2I]+[0001L]V+[-1J1000]T.
The above system has one control input which is the armature voltage of the electric drive V and one disturbance input which is the torque-on-bit T. Our objective in this section is to design a control law for V which would track the reference angular velocity of the drill ω1→ωref while rejecting the disturbance T.To solve the control problem formulated above, one can use the approach to feedforward robust servo control problem presented in [14, 15]. Below, the above approach is described in a simplified manner which, however, serves our purpose well. Consider a linear time invariant system of the form
(23)x˙=Ax+Bu+Dw,y=Cx+Fu+Hw,
where x∈ℝn is the state, u∈ℝm is the control input, y∈ℝp is the output, w∈ℝr are the disturbances, and A, B, C, D, F, and H are matrices of appropriate dimensions. Consider a control problem described as follows. Suppose the disturbances w(t) are measurable. Given a desired output signal yref(t), design a control algorithm that guarantees y(t)→yref(t) as t→+∞. This problem was addressed in [14, 15] in a very general setting. In this work, a simple case is addressed where both yref and w(t) are assumed to be constant signals, yref(t)≡yref and w(t)≡wm. In this case, the following two conditions are necessary and sufficient for the existence of a linear time-invariant controller that solves the above described problem.(i)
The pair(A,B) is stabilizable, which means that
(24)rank[B,AB,A2B,…,An-1B]=n;(ii)
Consider(25)rank[ABCF]=n+p.If the above two conditions hold (and only in this case), the linear time-invariant controller that solves the above described problem is given according to the formula(26)u=Kx+𝒢†yref+𝒢*wm,
where K∈ℝn×n is the feedback gain matrix which is to be chosen such that A-BK is stable and has the required dynamic properties
(27)𝒢=-C(A-BK)-1B,(28)𝒢*=𝒢†C(A-BK)-1D,
where 𝒢† is the Moore-Penrose pseudoinverse of the matrix 𝒢 in (27), defined by the formula
(29)𝒢†=𝒢T(𝒢𝒢T)-1.The above described control approach can be applied to the problem of stabilization of the angular velocity of drilling as follows. Equations (22), which describe the rotational dynamics of a drilling system, can be rewritten in the form (23), where x:=[ω1ϕω2I]T∈ℝ4, u:=V∈ℝ1, y:=ω1∈ℝ1, w:=T∈ℝ1, and the corresponding matrices are
(30)A:=[-c1J1kJ100-10100-kJ2-c2J2KnJ200-KnL-RL],B=[0001L],D=[-1J1000],(31)C=[1000],F=[0],H=[0].
Below, we consider the drilling system with specific values of the parameters that are listed in Table 1. With these values, the matrices A, B, and D become
(32)A:=[-0.11231.264700-10100-0.2231-0.20050.020400-8640-2],B=[000200],D=[-0.0027000],
while C, F, and H are given by (31).Table 1
Numerical values for drilling system parameters.
Parameter
Description
Value
Unit
J
1
BHA + drill string inertia
374
[kgm2]
J
2
Rotary table + drive inertia
2120
[kgm2]
c
1
BHA damping
42
[
Nms
/
rad
]
c
2
Rotary table damping
425
[
Nms
/
rad
]
k
Drill string stiffness
473
[
Nm
/
rad
]
R
Motor armature resistance
0.010
[
Ω
]
L
Motor armature inductance
0.005
[H]
K
Motor constant
6
[Vs]
n
Combined gear ratio for bevel and gear box
7.2
—
a
Drill bit radius
0.108
[m]
ζ
Ratio of drilling strength to drilling specific energy
0.7
—
M
Mass of drill string (28120 Kg) + BHA (25080 Kg)
53000
[kg]
W
s
-
H
0
Submerged weightWs  − applied  weight  from top of the Rig H0
100 or 1000
[N]
K
f
Viscous friction coefficient
20
[
Nm
/
rad
]For the above system, the necessary and sufficient conditions for stabilization (24), (25) are satisfied. Indeed, the stabilizability condition (24) is satisfied since
(33)rank[B,AB,A2B,…An-1B]=rank[0005.15427300.0000004.075472-8.96804.075-8.968-700.339200-400-34412146307]=4.
On the other hand, the rank condition (25) is also satisfied because
(34)rank[ABCF]=rank[-0.1122991.264706000-101000-0.223113-0.2004720.020377000-8640-220010000]=5.
Therefore, a controller of the form (26), (27), (28), and (29) guarantees that the angular velocity of the drill approach the reference angular velocity ω1→ωref as t→∞, while rejecting the disturbance Tb.The design of controller (26), (27), (28), and (29) begins by choosing the desired location of the closed-loop system's poles. For the purpose of simulations presented below, we consider two specific set of poles. The first set, denoted by P1, is chosen as follows:
(35)P1:=[-10-2+2i-2-2i-4].
The set P1 consists of two real poles and two complex conjugate poles. On the other hand, the set P2 contains only poles on the real axis as follows:
(36)P2=[-5.5-2-4.5-1].
The feedback gain matrix K1 such that the poles of A-BK1 are located according to P1 is
(37)K1=[32.2457.45-19.410.0784].
The coefficients 𝒢*, 𝒢† in (26) are calculated according to the formulas (27)–(29); the results are
(38)𝒢1*=0.123497,𝒢1†=60.0844.
On the other hand, the feedback matrix K2 such that the poles of A-BK2 are located according to P2 is
(39)K2=[-5.16716.943-30.620.0534].
The corresponding coefficients 𝒢2*, 𝒢2† are
(40)𝒢2*=0.037286,𝒢2†=9.603682.


## 3.3. Rock Stiffness Estimation

In the controller design presented above, it was assumed that the “hardness” of the rock, represented by the intrinsic specific energyϵ, is constant and exactly known. This knowledge of ϵ was used explicitly in the controller design, in particular, in formula (21). In practical geological drilling, however, the hardness of different layers of rock lying underneath the surface can be different and usually is not exactly known beforehand. More specifically, different characteristics of the rock, such as hardness, density and porosity, typically remain constant through each layer, but differs from layer to layer. On the other hand, control engineers frequently deal with the problem of designing a controller without a priori knowledge of the exact values of one or more parameters involved in the process. Often, the processes can be robustly controlled without the actual knowledge of some of the parameters. In other cases, the unknown parameters can be identified using specially designed estimators. Below, a simple online estimator of the rock intrinsic specific energy ϵ is designed following the methods described in [16], and the resulting estimate is then used in the controller for for drilling system.Specifically, during the cutting process, the torque-on-bitT is produced by bit rock interaction, according to the formula
(41)T=12a2ϵd,
where a is the radius of drill bit, d is the depth of cut, and ϵ>0 is the intrinsic specific energy. The intrinsic specific energy ϵ>0 depends on the properties of the media and typically unknown beforehand. However, since the torque on bit T(t) can typically be measured with advanced transducers located in the bottom hole assembly [17], a>0 is constant and known, and d(t) can be calculated according to formula (13), one can use the method described in the previous section to design an online estimation scheme for ϵ. In particular, considering (1/2)a2d(t) as the input and torque-on-bit T as the measured output, one can follow the procedure described in the previous section to design an estimator for an unknown parameter ϵ. The predicted torque-on-bit T^b is defined according to the formula
(42)T^(t):=12a2ϵ^d(t),
where ϵ^(t) is the current estimate of actual rock strength ϵ. The algorithm for online estimation of the intrinsic specific energy ϵ has a form
(43)ϵ^˙=γ0(T-T^)12a2d,
where γ0>0 is an arbitrary gain.A natural question regarding the algorithm (43) is if it guarantees the convergence of the parameter estimate to the true value of the parameter ϵ; mathematically, is ϵ^(t)→ϵ as t→+∞. It is known [16], that the convergence can be guaranteed if the “input” signal (1/2)a2d(t) is persistently exciting. A signal u(t) is said to be persistently exciting which is to say that there exist α0>0, T0>0 such that the inequality
(44)‍∫tt+T0u2(τ)dτ≥α0T0
holds for all t. In particular, u(t) is persistently exciting if u2(t)≥α0 for all t. Since d(t) is the depth of cut, we see that, during normal cutting process, d(t)≥d0>0, which results in persistent excitation of the input (1/2)a2d(t). The parameter convergence ϵ^(t)→ϵ, therefore, is guaranteed during normal cutting process. This is also confirmed by the simulation results presented below.The obtained estimate of the rock strengthϵ^ is then used in the control algorithm. Specifically, in the original formulation of the control algorithm, for a given reference vertical velocity vref, the reference rotational velocity ωref is calculated according to formula (21), which depends on the parameter ϵ. In case ϵ is unknown, it is substituted by its estimate ϵ^(t) obtained above. The new formula for ωref has the form
(45)ωref:=2πaζϵ^((Ws-H0)/vref)-Kf.
The obtained estimate of the rock stiffness ϵ^ will also be used to update the stiffness of the virtual spring in the haptic teleoperator drilling system described below.


## 3.4. Simulation Results

In this Section, some results of simulations of the drilling control system with intrinsic specific energy estimator are presented. The vertical motion of the drilling system is described by (15), and it is interconnected with the rotational dynamics (5) through nonlinear equation (13) that describes the depth of cut d(t). For a given reference velocity of the vertical penetration vref>0, the corresponding reference rotational velocity ωref is calculated according to formula (45). The controller (21), (26)–(29) has been implemented to guarantee that the angular velocity of the drill bits ω1(t) tracks ωref, which in turn stabilizes the vertical penetration velocity v(t) converges to vref. The algorithm (43) provides an estimate of the intrinsic specific energy parameter ϵ^ which is then used in the calculation of the reference angular velocity according to formula (45). Specific values of the parameters appearing in these equations are given in Table 1. The simulations are carried out using MATLAB, where the integration step for each simulation is equal to 0.005 s. The feedback gain matrix is chosen K=K1, where K1 is defined by (37).In the simulations described below, the performance of the system was evaluated for different values of actual intrinsic specific energyϵ, different gains γ0 and different values of the applied weight W0:=Ws-H0. Figures 5 and 6 show the response of the vertical penetration velocity v(t), the intrinsic specific energy estimate ϵ^(t), the torque-on-bit T(t), the predicted value of the torque-on-bit T^(t), and the rotational velocity ω1(t), all for the case where the applied weight on bit W0=5000 N, the intrinsic specific energy ϵ=20 MPa, and the desired vertical velocity vref is set to 0.005 m/s. The estimator gain is set to γ0=5·109. The plots show that v(t) converges to vref in less than 8 sec whereas the estimate ϵ^(t) converges to the actual value of ϵ in less than 4 sec. Figures 7 and 8 show the output responses of described parameters where W0=2500 N and the desired vertical velocity vref is set to 0.01 m/s. It can be clearly seen that the convergence becomes slower with reducing the applied weight on the drill string W0; specifically, both v(t) and ϵ^(t) approach their reference values in about 12 sec. Reducing W0 also results in that ω1ref increases, the steady-state value of Tb(t) drops to around 200 N, and the steady-state value of d(t) also drops to less than 2 mm. On the other hand, Figures 9 and 10 demonstrate the response of the system with the same parameters except the intrinsic specific energy ϵ is reduced to 5 MPa. This results in decreased convergence time for v(t) and ϵ^(t). The steady state value of rotational velocity ω1(t) is also decreased to under 10 rad/s, and steady state value of the depth of cut d(t) is increased to 6.5 mm. Figures 11 and 12 present the output response for the case where the estimator gain is decreased to γ0=5·108, while the rest of the parameters are the same as in the last simulation except the intrinsic specific energy is set to ϵ=10 MPa. The resulting response is predictably characterized by much slower convergence, which takes about 25 sec for v(t) and ϵ^(t) to approach their steady-state values. Finally, Figures 13 and 14 correspond to to the case where W0=5000 N, ϵ=20 MPa, and γ0=1·108.Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=5000 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) vresus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=5000 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=20 MPa, and γ0=5·109.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=5 MPa, and γ0=5·109.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=5 MPa, and γ0=5·109.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=10 MPa, and γ0=5·108.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=10 MPa, and γ0=5·108.
(a)
(b)
(c)Response of the vertical velocityv(t) (a) and the intrinsic specific energy estimate ϵ^(t) (b) for W0=2500 N, ϵ=20 MPa, and γ0=108.
(a)
(b)Response of rotational velocityω1(t) (a), torque-on-bit T(t) versus estimated torque-on-bit T^(t) (b), and the depth of cut d(t) (c) for W0=2500 N, ϵ=20 MPa, and γ0=108.
(a)
(b)
(c)Overall, simulation results show that the control system with intrinsic specific energy estimation demonstrate good stability and performance characteristics for a wide range of the parameters. In particular, the vertical velocity converges to the desired value, and the estimate of the intrinsic specific energyϵ^(t) converges to an actual value of ϵ.


## 4. Telerobotic Drilling System with Haptic Feedback

In this section, a telerobotic drilling system with haptic feedback is designed and experimentally evaluated. Haptics can be defined as the physical or virtual interaction through touch sensation for the purpose of perception and manipulation of objects  [18, 19]. Haptic feedback provides the operator with kinaesthetic clues of the physical features of virtual or real remote environment. The structure of a telerobotic drilling system with haptic feedback is shown in Figure 15. In this system, the human operator controls the drilling process using a haptic device. Specifically, the position of an end-effector of the haptic device defines the reference vertical velocity of the drilling. The reference vertical velocity is then transmitted to the drilling control system, designed in above in Section 3, which stabilizes the actual vertical penetration velocity to the level equal to the reference vertical velocity. On the other hand, an estimate of the intrinsic specific energy ϵ(t), which is generated online by an estimator described in above in Section 3.3, is sent back to the haptic device. The end-effector of the haptic device interacts with a virtual spring of variable stiffness; the stiffness of this virtual spring is updated in real time proportionally to the current estimate of the intrinsic specific energy ϵ(t). Thus, the telerobotic drilling system provides haptic feedback to the human operator which creates an intuitive feeling of the hardness of the remotely drilled material.Figure 15
The structure of a telerobotic drilling system.


### 4.1. Experimental Setup

The above described telerobotic drilling system is implemented in a semiexperimental setup, as follows. The setup consists of a PC based on Intel Pentium 4 processor with operational frequency of 1 GHz and RAM of 1 GB, and a PHANTOM Omni Haptic device, a product from SensAble Technologies Inc. The PHANTOM Omni Haptic device is designed for kinematic interaction with the virtual or real environment while providing the kinesthetic feedback to the operator. The device is equipped with a pen-based stylus, and is able to provide three degrees-of-freedom force feedback. The human operator uses the haptic device to (i) generate a desired vertical velocityvref(t) which is used as an input to the drilling control system, and (ii) to haptically perceive the stiffness ϵ of the rock layers. The remaining parts of the above described telerobotic system, including the drill string and drive system, the drilling process, as well as the control and estimation algorithms, are simulated in real time in virtual environment which is implemented using the Open Haptics tool kit and Microsoft Visual C++.The human operator assigns the desired velocityvref(t) by controlling the position of the end-effector of the PHANTOM device along its vertical (Y) axis. More exactly, a specific range along the y-axis is assigned to each desired vertical velocity, as follows:vref(t)=0.001 m/s if the position of stylus yn(t) is ≥80 mm;vref(t)=0.003 m/s if the position of stylus yn(t) is ≤80 mm and ≥60 mm;vref(t)=0.005 m/s if the position of stylus yn(t) is ≤60 mm and ≥40 mm;vref(t)=0.008 m/s if the position of stylus yn(t) is ≤40 mm and ≥25 mm;vref(t)=0.01 m/s if the position of stylus yn(t) is ≤25 mm and ≥10 mm;vref(t)=0.015 m/s if the position of stylus yn(t) is ≤10 mm and ≥0 mm;vref(t)=0.018 m/s if the position of stylus yn(t) is ≤0 mm.Another function of the haptic device is to allow the human operator to feel the stiffness of the rocks. As explained above, this is achieved by updating the stiffness of the virtual spring proportionally to the current estimate of the rock stiffness (intrinsic specific energyϵ^(t)). The coefficient of proportionality between the estimate of the intrinsic specific energy (with units of Pascals) and the stiffness of the virtual spring (with units of is N/m) is set in our experiments equal to 10−7. The feedback force Fest(t) due to the virtual spring is therefore calculated according to the formula
(46)Fest(t)=10-7·ϵ^(t)·yn(t).


### 4.2. Experimental Results

In this Section, some experimental results are presented. In these experiments, we have attempted to simulate a real drilling case scenario where the composition, characteristics, and types (which all contribute to the intrinsic specific energy) of various rock strata vary at different depths during drilling. Specifically, in every experiment, several layers of rocks with different intrinsic specific energyϵ ranging from 4 MPa to 60 MPa are simulated. In all experiments presented here, the applied weight W0=5000 N; the rest of the parameters, if not explicitly mentioned, are same as in Section 3.4.In the experiment shown in Figures16, 17, and 18, three layers of rocks with different intrinsic specific energy ϵ are simulated. The top layer has the stiffness of 5 MPa and its thickness is 20 cm from the surface. The second layer has a stiffness value of 12 MPa and lies between 20 cm and 30 cm from the surface (total thickness is 10 cm). The third layer starts at the depth of 30 cm and continues downward. It has a stiffness value of 20 MPa. The experiment is performed with the estimator gain γ0=109. Figure 16 shows the actual intrinsic specific energy ϵ(t) and its estimate ϵ^(t) on the top graph, and the reflected force Fest(t) on the bottom graph. Due to high estimator gain, ϵ^(t) quickly tracks ϵ(t) for all three layers as the drill bit progressed cutting through these layers. Figure 17 shows the vertical velocity vout(t) and the reference vertical velocity vref(t) on the top graph, and the reference rotational velocity ω1d(t) and the actual drill bit rotational velocity ω1(t) at the bottom graph. Finally, Figure 18 shows the behaviour of the actual torque-on-bit T(t) and the estimated torque T^(t), along with depth of cut d(t). These plots show that the system is stable and demonstrates good performance; in particular, all the output variables track their desired (reference) trajectories.Experiment 1: Actual stiffnessϵ(t) versus the estimated stiffness ϵ^(t) (a); the reflected force Fest(t) (b).
(a)
(b)Experiment 1: Output vertical velocityvout(t) versus reference vertical velocity vref(t) (a); output rotational velocity of the drill bit ω1(t) versus reference rotational velocity ω1d(t) (b).
(a)
(b)Experiment 1: Torque-on-bitT(t) versus estimated torque-on-bit T^(t) (a); depth of cut dcut(t) (b).
(a)
(b)Another set of experimental results is presented in Figures19–21, where the estimator gain is increased to γ0=5·109, and the depth of the rock layers and their corresponding stiffness values have been altered. Specifically, the first rock layer has depth 20 cm and the intrinsic specific energy ϵ is set to 20 MPa for this layer. Second layer lies between 20 cm and 40cm with ϵ=40 MPa. The third layer lies below 40 cm, and its intrinsic specific energy ϵ=60 MPa. Figure 19 shows the corresponding plots of ϵ(t), ϵ^(t) and Fest(t). Figure 20 shows the response of vref(t) and vout(t) on the top graph, and the responses of ω1d(t) and ω1(t) on the bottom graph, respectively. The response of T^(t) and T(t) along with d(t) are shown in Figure 21. Overall, our experiments demonstrate stability and good performance of the designed telerobotic drilling system with haptic feedback, for a wide range of parameters and control gains.Experiment 2: Actual stiffnessϵ(t) versus estimated stiffness ϵ^(t) (a); reflected force Fest(t) (b).
(a)
(b)Experiment 2: Output vertical velocityvout(t) versus reference vertical velocity vref(t) (a); output rotational velocity of the drill bit ω1(t) versus reference rotational velocity ω1d(t) (b).
(a)
(b)Experiment 2: Torque-on-bitT(t) versus estimated torque-on-bit T^(t) (a); depth of cut dcut(t) (b).
(a)
(b)


## 4.1. Experimental Setup

The above described telerobotic drilling system is implemented in a semiexperimental setup, as follows. The setup consists of a PC based on Intel Pentium 4 processor with operational frequency of 1 GHz and RAM of 1 GB, and a PHANTOM Omni Haptic device, a product from SensAble Technologies Inc. The PHANTOM Omni Haptic device is designed for kinematic interaction with the virtual or real environment while providing the kinesthetic feedback to the operator. The device is equipped with a pen-based stylus, and is able to provide three degrees-of-freedom force feedback. The human operator uses the haptic device to (i) generate a desired vertical velocityvref(t) which is used as an input to the drilling control system, and (ii) to haptically perceive the stiffness ϵ of the rock layers. The remaining parts of the above described telerobotic system, including the drill string and drive system, the drilling process, as well as the control and estimation algorithms, are simulated in real time in virtual environment which is implemented using the Open Haptics tool kit and Microsoft Visual C++.The human operator assigns the desired velocityvref(t) by controlling the position of the end-effector of the PHANTOM device along its vertical (Y) axis. More exactly, a specific range along the y-axis is assigned to each desired vertical velocity, as follows:vref(t)=0.001 m/s if the position of stylus yn(t) is ≥80 mm;vref(t)=0.003 m/s if the position of stylus yn(t) is ≤80 mm and ≥60 mm;vref(t)=0.005 m/s if the position of stylus yn(t) is ≤60 mm and ≥40 mm;vref(t)=0.008 m/s if the position of stylus yn(t) is ≤40 mm and ≥25 mm;vref(t)=0.01 m/s if the position of stylus yn(t) is ≤25 mm and ≥10 mm;vref(t)=0.015 m/s if the position of stylus yn(t) is ≤10 mm and ≥0 mm;vref(t)=0.018 m/s if the position of stylus yn(t) is ≤0 mm.Another function of the haptic device is to allow the human operator to feel the stiffness of the rocks. As explained above, this is achieved by updating the stiffness of the virtual spring proportionally to the current estimate of the rock stiffness (intrinsic specific energyϵ^(t)). The coefficient of proportionality between the estimate of the intrinsic specific energy (with units of Pascals) and the stiffness of the virtual spring (with units of is N/m) is set in our experiments equal to 10−7. The feedback force Fest(t) due to the virtual spring is therefore calculated according to the formula
(46)Fest(t)=10-7·ϵ^(t)·yn(t).


## 4.2. Experimental Results

In this Section, some experimental results are presented. In these experiments, we have attempted to simulate a real drilling case scenario where the composition, characteristics, and types (which all contribute to the intrinsic specific energy) of various rock strata vary at different depths during drilling. Specifically, in every experiment, several layers of rocks with different intrinsic specific energyϵ ranging from 4 MPa to 60 MPa are simulated. In all experiments presented here, the applied weight W0=5000 N; the rest of the parameters, if not explicitly mentioned, are same as in Section 3.4.In the experiment shown in Figures16, 17, and 18, three layers of rocks with different intrinsic specific energy ϵ are simulated. The top layer has the stiffness of 5 MPa and its thickness is 20 cm from the surface. The second layer has a stiffness value of 12 MPa and lies between 20 cm and 30 cm from the surface (total thickness is 10 cm). The third layer starts at the depth of 30 cm and continues downward. It has a stiffness value of 20 MPa. The experiment is performed with the estimator gain γ0=109. Figure 16 shows the actual intrinsic specific energy ϵ(t) and its estimate ϵ^(t) on the top graph, and the reflected force Fest(t) on the bottom graph. Due to high estimator gain, ϵ^(t) quickly tracks ϵ(t) for all three layers as the drill bit progressed cutting through these layers. Figure 17 shows the vertical velocity vout(t) and the reference vertical velocity vref(t) on the top graph, and the reference rotational velocity ω1d(t) and the actual drill bit rotational velocity ω1(t) at the bottom graph. Finally, Figure 18 shows the behaviour of the actual torque-on-bit T(t) and the estimated torque T^(t), along with depth of cut d(t). These plots show that the system is stable and demonstrates good performance; in particular, all the output variables track their desired (reference) trajectories.Experiment 1: Actual stiffnessϵ(t) versus the estimated stiffness ϵ^(t) (a); the reflected force Fest(t) (b).
(a)
(b)Experiment 1: Output vertical velocityvout(t) versus reference vertical velocity vref(t) (a); output rotational velocity of the drill bit ω1(t) versus reference rotational velocity ω1d(t) (b).
(a)
(b)Experiment 1: Torque-on-bitT(t) versus estimated torque-on-bit T^(t) (a); depth of cut dcut(t) (b).
(a)
(b)Another set of experimental results is presented in Figures19–21, where the estimator gain is increased to γ0=5·109, and the depth of the rock layers and their corresponding stiffness values have been altered. Specifically, the first rock layer has depth 20 cm and the intrinsic specific energy ϵ is set to 20 MPa for this layer. Second layer lies between 20 cm and 40cm with ϵ=40 MPa. The third layer lies below 40 cm, and its intrinsic specific energy ϵ=60 MPa. Figure 19 shows the corresponding plots of ϵ(t), ϵ^(t) and Fest(t). Figure 20 shows the response of vref(t) and vout(t) on the top graph, and the responses of ω1d(t) and ω1(t) on the bottom graph, respectively. The response of T^(t) and T(t) along with d(t) are shown in Figure 21. Overall, our experiments demonstrate stability and good performance of the designed telerobotic drilling system with haptic feedback, for a wide range of parameters and control gains.Experiment 2: Actual stiffnessϵ(t) versus estimated stiffness ϵ^(t) (a); reflected force Fest(t) (b).
(a)
(b)Experiment 2: Output vertical velocityvout(t) versus reference vertical velocity vref(t) (a); output rotational velocity of the drill bit ω1(t) versus reference rotational velocity ω1d(t) (b).
(a)
(b)Experiment 2: Torque-on-bitT(t) versus estimated torque-on-bit T^(t) (a); depth of cut dcut(t) (b).
(a)
(b)


## 5. Conclusions

This paper deals with control design for a teleoperator system with haptic feedback for an oil well drilling process. A mathematical model of the drilling process was described, and the control algorithm was designed that guarantees the convergence of the vertical penetration velocity to an arbitrary reference value. The control algorithm has a cascaded structure, where the velocity of vertical penetration is controlled indirectly through stabilization of the rotational motion of the drill bit. In order to guarantee the convergence of the angular velocity to a desired value in the presence of disturbances in the form of torque-on-bit, a robust servo controller was designed. However, the design of such controller depends on the parameter of environment called the intrinsic specific energy, which is generally unknown beforehand. To solve this issue, an online parameter estimator was designed that provides an estimate of the intrinsic specific energy. This estimate is substituted for the actual value of the parameter in the control algorithm, and the corresponding adaptive control system is evaluated through simulations. Finally, a telerobotic drilling system with haptic feedback is designed and verified through semi-experiments. The haptic feedback for the human operator is provided by creating a virtual spring that interacts with the haptic device; the stiffness of the spring is adjusted in real time depending on the current estimate of the intrinsic specific energy. Semi-experiments are conducted using PHANTOM Omni Haptic device, where the drilling process model is implemented in C++ environment, and the haptic feedback is provided to the human operator.There exists a number of challenges associated with the real-life drilling operation that were not addressed in our paper. In particular, the frictional forces at the contact were neglected in our analysis, while in reality they may play significant role in the drilling process. Also, in real-life drilling systems, the rotational velocity and the penetration rate are typically measured at the surface while the torque-on-bit should ideally be measured above the bit; thus, there exists a problem of synchronizing the data obtained at the surface with those obtained at the bit. The issue of communication delay between the haptic device and the drilling process is also not addressed. These, as well as more detailed experimental evaluation of the designed system, are the topics for future research.



---

*Source: 901610-2013-10-05.xml*