## Description

“The Armageddon Movie Problem”

This assignment provides you with hands-on practice for state estimation, with application to oil & gas

energy systems. The assignment is organized in a tutorial fashion, thereby allowing you to practice state

estimation on a relevant real-world energy system example.

Figure 1: Schematic view of a drilling system [1].

Movie-goers will recognize this problem from the

1998 science fiction thriller “Armageddon.” In this

movie, Bruce Willis, an experienced and gritty offshore oil-driller, must drill a hole into an asteroid

plummeting towards Earth to destroy it. If only he

knew about state estimation, then maybe he could

prevent the insufferable Ben Affleck from marrying

Liv Tyler… and save Earth.

## Reading

The following four articles on bCourses describe

state estimation in various energy applications: batteries [2], buildings [3], traffic [4], and geophysical

systems [5]. Please peruse for your own self-interest.

There are no questions regarding this reading.

## Background

The process of oil well drilling involves creating a

borehole several thousand meters deep in the ground

until an oil reservoir is reached. The drill string,

consisting of the assembly of drill pipes, drill collars,

and the rock-cutting tool referred to as drill bit (see

Figure 1) rotates around its vertical axis, penetrating through the rock. At the top of the drillstring,

the rotary table provides the necessary torque to put

the system into a rotational motion.

During this process, the downhole drill bit will stick and slip, causing oscillations in the drill string and

creating mechanical fatigue to the bit itself. For this reason, it is crucial to monitor the drill bit speed.

However, it’s nearly impossible to place sensors thousands of meters into the ground, especially in the

perilous environment near the bit which contains flows of rock cuttings, oil, gas, water, and mud. For this

Page 1 of 4

CE 295: Energy Systems and Control

reason, we seek to estimate the bit velocity by fusing a model of the drill string dynamics, and measurements

of the rotary table speed.

T(t), Torque

b·ωT (t) Viscous

Drag

b·ωB (t) Viscous

Drag

k·[θT(t)-θB (t)]

Rotational spring

ωT (t)

ωB (t)

JT Table/Top

Rotational Inertia

JB

Bit/Bottom

Rotational Inertia

Table/Top

Bit/

Bottom

Tf (t), Frictional Torque

Figure 2: Free-body diagram of drill string.

## Problem 1: Dynamic System Modeling

We model the drill-string as two rigid rotating bodies: the

table/top portion and the bit/bottom portion, connected by a

rotational spring, as shown in Fig. 2. The top and bottom

have rotational inertia JT and JB, respectively. An electric

motor provides the table torque T(t). The rock provides some

unmeasurable frictional torque Tf (t) on the drill bit at the

bottom. The top and bottom also experience torques from

viscous drag and the rotational spring.

(a). Define & write the modeling objective. What are the

controllable and uncontrollable inputs? What are the

measured and performance outputs? List the parameters.

(b). Use Newton’s second law in rotational coordinates to derive the equations of motion for the top/table and bottom/bit portions of the drill string. HINT: Fig 2 is a free

body diagram.

(c). Write all the dynamical equations into matrix state space

form. What are the A, B, C matrices? Hint: A ∈ R

4×4

.

## Problem 2: Observability Analysis

In this problem, we consider IF it’s possible to estimate the bit velocity from table velocity measurements.

Assume the following parameter values: JT = 100, JB = 25, k = 2, b = 5.

(a). Compute the Observability Matrix O described in Section 2 of Chapter 2, and include in your submitted

report. Are all the states observable from measurements of table velocity? Justify why or why not.

(b). Define a new state θ(t) = θT (t)−θB(t) as the relative angular displacement between the top and bottom.

Re-derive the state-space model replacing θT , θB with θ only. Provide the new A, B, C matrices.

(c). Compute the Observability Matrix O for this modified model, and include in your report. Is the bit

velocity observable from measurements of the table velocity? Justify why or why not.

## Problem 3: Measurement Data

Download HW2_Skeleton.m/HW2_Skeleton.ipynb and HW2_Data.csv from bCourses. Import HW2_Data.csv.

In one figure, with two subplots, plot the table torque T(t) and measured table velocity ωT (t) versus time.

Use legends, x- and y-labels, appropriate font sizes, and line widths. Provide the plot in your report.

## Problem 4: Luenberger Observer

In this problem we design a Luenberger observer to estimate bit speed from the 3-state model and measurements. You are encouraged to experiment with the eigenvalue placement to obtain a good balance between

estimation speed and robustness to noise.

Page 2 of 4

CE 295: Energy Systems and Control

(a). Compute the eigenvalues of the open loop system, A. Write down the Luenberger observer equations

in your report.

(b). Design the output error injection gain, L, as described in Chapter 2, Section 3. Select some desired

eigenvalues for the error system matrix (A − LC). Remember, the only requirements are that (i) the

real parts must be strictly negative (Re[λi(A − LC)] < 0, ∀ i), and (ii) complex numbers must be in

complex conjugate pairs. Provide your selected eigenvalues and observer gain L in the report. Justify

your choice of eigenvalues. HINT: You need the place and acker commands in Matlab and Python,

respectively. These require the Matlab Control Systems Toolbox or Python Control Systems module.

(c). Write down the state-space matrices for your Luenberger observer, denoted Alobs, Blobs, Clobs. Note, the

frictional torque is unknown and unmeasured and therefore not included in the Luenberger observer.

(d). In one figure, create two subplots. On top, plot the estimated bit speed ωˆB(t) versus the true value

omega_B_true from HW2_Data.csv. On bottom, plot the estimation error ωB(t) − ωˆB(t). Compute

the root mean square error (RMSE), given by:

RMSE =

s

1

T

X

t

[ωB(t) − ωˆB(t)]2

(1)

where T is the total simulation time, i.e. t ∈ [0, T]. Use these plots and the RMSE to tune the

eigenvalues of (A − LC). Provide this plot in your report, and document your eigenvalue placement

for (A − LC). HINT: There is no single best answer. Tune to your preference.

## Problem 5: Kalman Filter (KF) Design

In this problem we design a Kalman Filter, using the canonical model format

x˙(t) = Ax(t) + Bu(t) + Bww(t), (2)

ym(t) = Cx(t) + n(t), (3)

where x(t) = [θ(t), ωT (t), ωB(t)], u(t) = T(t), w(t) = Tf (t), and ym(t) = ωT (t) + n(t). The process noise

w(t) encapsulates unknown frictional torque on the drill bit. The measurement noise n(t) represents noise

in the table velocity sensor ωT (t).

(a). Write down the Kalman filter equations from Chapter 2, Section 4.

(b). Implement the Kalman filter equations in your code. Assume the measurement noise covariance is

N = 0.02, the mean initial condition is x0 = 03×1, and the initial condition covariance is Σ0 = I3×3.

Tune the process noise covariance W ∈ R

3×3

, and provide the value in your report. Explain how you

tuned covariance matrix W.

(c). In one figure, create two subplots. On top, plot the estimated bit speed ωˆB(t) versus the true value

omega_B_true

p

from HW2_Data.csv. Also plot the ± one-standard deviation bound, i.e. ωˆB(t) ±

Σ3,3(t). On bottom, plot the estimation error ωB(t) − ωˆB(t). Use these plots to tune the process

covariance W. Provide this plot in your report.

Page 3 of 4

CE 295: Energy Systems and Control

(d). Compute the eigenvalues of A − L(t)C, where t is the final time of 300 sec. How does this compare to

your pole placement in the Luenberger observer? FOR FUN (ungraded): Generate an animated plot

of the λi(A − L(t)C) on the real-complex plane. Watch the eigenvalues move.

## Problem 6: Extended Kalman Filter (EKF) Design

Now we consider a nonlinear rotational spring law. In Fig. 2 and Problem 1, we relate spring displacement

θ(t) and spring torque by Hooke’s law: kθ(t). This is a linear relationship. In reality, spring torque exhibits

a nonlinear relationship. In this problem, let us replace Hooke’s law with the more realistic nonlinear spring

torque relationship: k1θ(t) + k2θ

3

(t). Consider values of k1 = 2 and k2 = 0.25.

Design and implement an Extended Kalman Filter from CH2, Section 5 by copying, pasting, and modifying

your KF code. Use the covariance W you used from Problem 5. In your report:

• Derive expressions for the Jacobians: F(t) = ∂f

∂x (ˆx(t), u(t)), H(t) = ∂h

∂x (ˆx(t), u(t)).

• Write down the EKF equations.

• In one figure, create two subplots like Problem 5. On top, plot the estimated bit speed ωˆB(t) versus

the true value omega_B_true from HW2_Data.csv. Compute the RMSE. Also plot the ± one-standard

deviation bound, i.e. ωˆB(t) ±

p

Σ3,3(t). On bottom, plot the estimation error ωB(t) − ωˆB(t).

## Deliverables

Submit the following on bCourses. Be sure that the function files are named exactly as specified (including

spelling and case). Please do NOT zip files.

LASTNAME_FIRSTNAME_HW2.PDF

LASTNAME_FIRSTNAME_HW2.xyz which contains your respective Matlab or Python files, i.e. xyz ∈ {m, ipynb}.

## References

[1] C. Sagert, F. Di Meglio, M. Krstic, and P. Rouchon, “Backstepping and flatness approaches for stabilization of

the stick-slip phenomenon for drilling,” in Proceedings of the 19th IFAC World Congress: System, Structure and

Control, vol. 5, no. 1, 2013, pp. 779–784.

[2] S. J. Moura and H. Perez, “Extracting the Full Potential of Batteries: Electrochemistry and Controls,” ASME

Dynamic Systems and Control Magazine, vol. 2, no. 2, pp. S15–S21, July 2014.

[3] P. Radecki and B. Hencey, “Online building thermal parameter estimation via unscented kalman filtering,” in

2012 American Control Conference (ACC), 2012, pp. 3056–3062.

[4] D. Work, O.-P. Tossavainen, S. Blandin, A. Bayen, T. Iwuchukwu, and K. Tracton, “An ensemble kalman filtering

approach to highway traffic estimation using gps enabled mobile devices,” in 47th IEEE Conference on Decision

and Control, Dec 2008, pp. 5062–5068.

[5] J. Anderson, “Ensemble kalman filters for large geophysical applications,” IEEE Control Systems Magainze, vol. 29,

no. 3, pp. 66–82, June 2009.

Page 4 of 4