Description
Consider the mass-spring system given by the second order differential equation
x
00(t) + x(t) = 0 (1)
with the initial conditions
x(0) = x0, x0
(0) = v0 (2)
Introducing the velocity variable v(t) = x
0
(t), an equivalent formulation to (1-2) is
x
0
v
0
=
0 1
−1 0 x
v
,
x(0)
v(0)
=
x0
v0
(3)
The solution to the initial-value problem (3) is expressed as
x(t) = v0 sin(t) + x0 cos(t), v(t) = v0 cos(t) − x0 sin(t) (4)
The solution (4) is denoted in vector format as
x(t) =
x(t)
v(t)
(5)
and we refer to x(t) ∈ R
n
(here n=2) as the true state of the dynamical system. The initial condition is
specified as
x(0) =
x(0)
v(0)
=
1
1
(6)
Henceforth, it is assumed that we can obtain information on the true state x(t) through state measurements and a time discrete model. Observations are represented as
y = Hx + o ∈ R
m, (m = 1 or m = 2) (7)
where H ∈ R
m×n denotes the observation matrix and o ∈ R
m is a random vector of observational errors.
A time-discrete model to the true dynamics is obtained by applying the Euler’s method to the system
(3) with a time step ∆t
xˆk+1 = Mxˆk (8)
where
M =
1 ∆t
−∆t 1
(9)
The evolution of the true state xk = x(tk) in (5) is given by
xk+1 = Mxk + wk (10)
where wk ∈ R
n
is an unknown vector consisting of numerical discretization errors.
Our goal is to implement the Kalman filter to estimate and predict the true state xk based on the
model (8) and observations (7). The setup is as follows.
• A prior guess to the true (unknown) initial state at t = 0 is prescribed as
xˆ0 = x(0) + b
where the initial condition error in each component is b ∼ N(0, σ2
b
) with σb = 1.
• The numerical discretization time step is taken ∆t = 1.
• The observational operator is defined as
1. Case 1 (both state components are observed)
H =
1 0
0 1
(11)
2. Case 2 (only first state component is observed)
H = [1 0] (12)
• The observational error is specified as o ∼ N(0, σ2
o
) where σo = 0.1.
• Each component of the model error is specified as wk ∼ N(0, σ2
q
) where σq = ∆t = 1
• The time interval for the analysis it 0 ≤ t ≤ tN , that is k = 0 : N − 1 in (8).
The KF algorithm is formulated as follows:
—————————————————————
x
a
0 = x(0) + b % true initial state plus random error taken from N(0, 1)
Pa = I % initial analysis error covariance matrix is set to identity
for k = 0 : 999
% forecast step
x
f
k+1 = Mx
a
k
Pf = MPaMT + Q % forecast error covariance where Q = σ
2
q
In×n
% generate observation(s) at tk+1
y = Hx(tk+1)+o % generate observation(s) at tk+1 of true state plus random errors taken from N(0, σ2
o
)
% analysis step
x
a
k+1 = x
f
k+1 + K(y − Hx
f
k+1) % where K is the Kalman gain matrix
Pa = (I − KH)Pf % analysis error covariance at tk+1
end % end for loop
—————————————————————
Your job:
(10 points) Investigate the performance of the discrete model in pure forecast mode.
(30 points) Implement the KF algorithm, as described above.
(30 points) Provide a qualitative study of the analysis and forecast errors in the time interval 0 ≤ t ≤ 1000, that
is for N = 1000 time steps.
