[SOLVED] ECE5550 Applied Kalman Filtering STATE-SPACE DYNAMIC SYSTEMS Matlab

ECE5550: Applied Kalman Filtering

STATE-SPACE DYNAMIC SYSTEMS

2.1: Introduction to state-space systems

■ Representation of the dynamics of an nth-order system as a first-order differential equation in an n-vector called the state.

➠ n first-order equations.

■ Classic example: Second-order equation of motion.

■ Define a (non-unique) state vector (note that x˙(t) = dx(t)/dt, etc.)

■ We can write this as x˙(t) = Ax(t) + Bu(t), where A and B are constant matrices.

■ Complete the model by computing z(t) = Cx(t) + Du(t), where C and D are constant matrices.

■ Fundamental form. for deterministic, time-invariant, continuous-time linear state-space model:

x˙(t) = Ax(t) + Bu(t)

z(t) = Cx(t) + Du(t),

where u(t) is input, x(t) is the state, A, B, C, D are constant matrices.

• Systems with noise inputs are considered in notes chapter 3.

• Time-varying systems have A, B, C, D that change with time.

DEFINITION: The state of a system at time t0 is a minimum amount of information at t0 that, together with the input u(t), t ≥ t0, uniquely determines the behavior. of the system for all t ≥ t0.

■ State variables provide access to what is going on inside the system.

■ Convenient way to express equations of motion.

■ Matrix format great for computers.

■ Allows new analysis and synthesis tools.

SIMULATING IN SIMULINK: To investigate state-space systems, we can simulate them in Simulink. The block diagram below gives explicit access to the state and other internal signals. It is a direct implementation of the transfer function above, and the initial state may be set by setting the initial integrator values.

Example: The nearly constant position (NCP) model

■ Consider a relatively immobile object that we would like to track using a Kalman filter.

■ It gets bumped around by unknown forces.

■ We let our model state be

where ξ(t) is the x-coordinate and η(t) is the y-coordinate of position.

■ Our model’s state equation is then

x˙(t) = 0x(t) + w(t),

where w(t) is a random process-noise input (unlike known u(t)).

■ One possible output equation is

z(t) = x(t) + v(t),

where v(t) is a random sensor-noise input.

■ A possible Simulink implementation and output trajectory:

Example: The nearly constant velocity (NCV) model

■ Another model we might consider is that of an object with momentum.

■ The velocity is nearly constant, but gets perturbed by external forces.

■ We let our model state be

■ Our model’s state equation is then

■ One possible output equation is

■ A possible Simulink implementation and output trajectory:

Example: The coordinated turn model

■ A third model considers an object moving in a 2D plane with constant speed and angular rate Ω where Ω> 0 is counter-clockwise motion and Ω < 0 is clockwise motion.

ξ (¨ t) = −Ωη(˙ t)     and     η(¨ t) = Ωξ (˙ t),

■ We again let our model state be

■ Our model’s state equation is then

■ One possible output equation is

■ A possible Simulink implementation and output trajectory:

2.2: Time (dynamic) response

■ Develop more insight into the system response by looking at time-domain solution for x(t).

Homogeneous part

■ Start with x˙(t) = Ax(t) and some initial state x(0).

■ Take Laplace transform. X(s) = (s I − A)
−1
x(0).

■ So, we have: x(t) = L
−1
[(s I − A)
−1
]x(0). But,

so,

■ e
At : “Transition matrix” or “state-transition matrix.”

■ In MATLAB,

x = expm(A*t)*x0;

■ e
(A+B)t = e
Ate
Bt   iff   AB = B A. (i.e., not in general).

■ Will say more about e
At when we discuss the structure of A.

■ Computation of e
At = L
−1
[(s I − A)
−1
] straightforward for 2 × 2.

EXAMPLE: Find e
At when

■ Solve

■ This is the best way to find e
At if A 2 × 2.

Forced solution

■ Where did this come from?

■ Clearly, if z(t) = Cx(t) + Du(t),

More on the matrix exponential

■ Have seen the key role of e
At in the solution for x(t). Impacts the system response, but need more insight.

■ Consider what happens if the matrix A is diagonalizable, that is, there exists a matrix T such that T
−1AT = 3 =diagonal. Then,

and

■ Much simpler form. for the exponential, but how to find T, 3?

■ Write T
−1AT = 3 as T
−1A = 3 T
−1 with

wTi A = λiwi
T
, so wi
is a left eigenvector of A and note that wi
T
v j = δi,j
.

■ How does this help?

■ Very simple form, which can be used to develop intuition about dynamic response ≈ e
λi
t
.

■ Trajectory can be expressed as a linear combination of modes: vie
λi
t
.

■ Left eigenvectors decompose x(0) into modal coordinates wi
T
x(0).

■ e
λi
t
propagates mode forward in time. Stability?

■ vi corresponds to “relative phasing” of state’s part of the response.

EXAMPLE: Let’s consider a specific system

x˙(t) = Ax(t)

z(t) = Cx(t),

with x(t) ∈ R
16×1
, z(t) ∈ R (16-state, single output).

■ A lightly damped system.

■ Typical output to initial conditions are shown:

■ Waveform. is very complicated. Looks almost random.

■ However, the solution can be decomposed into much simpler modal components.

2.3: Discrete-time state-space systems

■ Computer monitoring of real-time systems requires analog-to-digital (A2D) and digital-to-analog (D2A) conversion.

■ Discrete-time systems can also be represented in state-space form.

xk+1 = Ad xk + Bduk

zk = Cd xk + Dduk.

■ The subscript. “d” is used here to emphasize that, in general, the “A”, “B”, “C” and “D” matrices are different for discrete-time and continuous-time systems, even if the underlying plant is the same.

■ I will usually drop the “d” and expect you to interpret the system from its context.

Time (dynamic) response

■ The full solution, found by induction from xk+1 = Axk + Buk, is

■ Clearly, if zk = Cxk + Duk,

Converting plant dynamics to discrete time.

■ Combine the dynamics of the zero-order hold and the plant.

■ The continuous-time dynamics of the plant are:

x˙(t) = Ax(t) + Bu(t)

z(t) = Cx(t) + Du(t).

■ Evaluate x(t) at discrete times. Recall

■ With malice aforethought, break up the integral into two pieces. The first piece will become Ad times x(kT ). The second part will become Bd times u(kT).

■ In the remaining integral, note that u(τ ) is constant from kT to (k + 1)T, and equal to u(kT ).

■ So, we let σ = (k + 1)T − τ ; τ = (k + 1)T − σ ; dτ = −dσ .

■ So, we have a discrete-time state-space representation from the continuous-time representation

xk+1 = Ad xk + Bduk

where

■ Similarly,

zk = Cxk + Duk.

• That is, Cd = C; Dd = D.

Calculating Ad, Bd, Cd and Dd

■ Cd and Dd require no calculation since Cd = C and Dd = D.

■ Ad is calculated via the matrix exponential Ad = e
AT . This is different from taking the exponential of each element in AT .

■ If MATLAB is handy, you can type in

Ad=expm(A*T)

■ If MATLAB is not handy, then you need to work a little harder. Recall from earlier that e
At = L
−1
[(s I − A)
−1
]. So,

which is probably the “easiest” way to work it out by hand.

■ Now we focus on computing Bd. Recall that

■ If A is invertible, this method works nicely; otherwise, we will need to perform. the integral.

■ Also, in MATLAB,

[Ad,Bd]=c2d(A,B,T)