Discretizing stochastic dynamical systems using Lyapunov equations

(1)

Discretizing stochastic dynamical systems using Lyapunov equations

Niklas Wahlström, Patrik Axelsson, Fredrik Gustafsson

Division of Automatic Control Linköping University

Linköping, Sweden

(2)

Motivation 2(15)

Stochastic dynamical systems are important in state estimation, system identification and control.

System models are often provided in continuous time, while a major part of the applied theory is developed for discrete-time systems.

Discretization of continuous-time models is hence fundamental.

AUTOMATIC CONTROL

(3)

Discretization 3(15)

Continuous-time model

˙x ( t ) = Ax ( t ) + Bw ( t ) E [ w ( t ) w ( τ ) ^T ] = Sδ ( t − τ )

Discrete-time model x _k + 1 = Fx _k + w _k

E [ w _k w ^T _l ] = Qδ _kl

(4)

Discretization 3(15)

Continuous-time model

˙x ( t ) = Ax ( t ) + Bw ( t ) E [ w ( t ) w ( τ ) ^T ] = Sδ ( t − τ )

Discrete-time model x _k + 1 = Fx _k + w _k

E [ w _k w ^T _l ] = Qδ _kl

This gives the relations

F = e ^AT

^s

, Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ

AUTOMATIC CONTROL

(5)

Discretization 3(15)

Continuous-time model

˙x ( t ) = Ax ( t ) + Bw ( t ) E [ w ( t ) w ( τ ) ^T ] = Sδ ( t − τ )

Discrete-time model x _k + 1 = Fx _k + w _k

E [ w _k w ^T _l ] = Qδ _kl

This gives the relations

F = e ^AT

^s

, Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ Problem

How do we solve this integral in a numerically good manner

for arbitrary A , B , S and T _s ?

(6)

Analytical solution - System has only integrators 4(15)

If the system has only integrators, A is nilpotent and e ^Aτ can be expressed exactly with a finite Taylor series

e ^Aτ = I + Aτ + · · · + ¹

p! A ^p ⁻ ¹ τ ^p ⁻ ¹ .

Then, the integral has an analytical solution

Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ

=

p − 1 i ∑ = 0

p − 1 j ∑ = 0

T ⁱ _s ⁺ ^j ⁺ ¹ i!j! ( i + j + 1 ) ^A

i S ( A ⁱ ) ^T .

Re Im

Figure : The poles of the system

AUTOMATIC CONTROL

(7)

Theorem 5(15)

Theorem - Contribution 1 The solution to the integral

Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ (1) satisfies the Lyapunov equation

AQ + QA ^T = − BSB ^T + e ^AT

^s

BSB ^T e ^A

^T

^s

(2)

(8)

Theorem 5(15)

Theorem - Contribution 1 The solution to the integral

Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ (1) satisfies the Lyapunov equation

AQ + QA ^T = − BSB ^T + e ^AT

^s

BSB ^T e ^A

^T

^s

(2)

Idea

Solve (2) to find solution for (1).

AUTOMATIC CONTROL

(9)

Theorem 5(15)

Theorem - Contribution 1 The solution to the integral

Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ (1) satisfies the Lyapunov equation

AQ + QA ^T = − BSB ^T + e ^AT

^s

BSB ^T e ^A

^T

^s

(2)

Idea

Solve (2) to find solution for (1).

Lemma Eq. (2) has a unique solution if and only if

λ _i ( A ) + λ _j ( A ) 6= 0 ∀ i, j

(10)

Theorem 5(15)

Theorem - Contribution 1 The solution to the integral

Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ (1) satisfies the Lyapunov equation

AQ + QA ^T = − BSB ^T + e ^AT

^s

BSB ^T e ^A

^T

^s

(2)

Idea

Solve (2) to find solution for (1).

Lemma Eq. (2) has a unique solution if and only if λ _i ( A ) + λ _j ( A ) 6= 0 ∀ i, j

Note: This is not fulfilled if the system has integrators!

AUTOMATIC CONTROL

(11)

Lyapunov equation - System with only strictly stable

poles 6(15)

If the system has only strictly stable poles, the Lyapunov equation has a unique solution

AQ + QA ^T = − BSB ^T + e ^AT

^s

BSB ^T e ^A

^T

^s

| {z }

− V

,

which is identical to the solution of the integral

Q =

Z _T

_s

0 e ^Aτ BSB ^T e ^A

^T

^τ dτ.

Re Im

Figure : The poles of the system

(12)

Summary so far 7(15)

Re Im

Q =

p − 1 i ∑ = 0

p − 1 j ∑ = 0

T _s ⁱ ⁺ ^j ⁺ ¹ i!j! ( i + j + 1 ) ^A

i S ( A ⁱ ) ^T

Re Im

AQ + QA ^T = − V

How do we find Q if we have both integrators and strictly stable poles?

AUTOMATIC CONTROL

(13)

Systems with integrators + strictly stable poles 8(15)

Consider a system on the following block triangular form

A =

A ₁₁ A ₁₂ 0 A ₂₂

where all integrators are collected in A ₂₂ , i.e.

λ _i ( A ₁₁ ) 6= 0 ∀ i λ _j ( A ₂₂ ) = 0 ∀ j

Re Im

Figure : The poles of the system

(14)

Systems with integrators + strictly stable poles 9(15)

The corresponding Lyapunov equation for this system reads

A ₁₁ A ₁₂ 0 A ₂₂

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

+ ^Q ¹¹ ^Q ¹² Q ^T ₁₂ Q ₂₂

A ^T ₁₁ 0 A ^T ₁₂ A ^T ₂₂

=− ^V ¹¹ ^V ¹² V ₁₂ ^T V ₂₂

,

which gives

A ₁₁ Q ₁₁ + Q ₁₁ A ^T ₁₁ = − V ₁₁ − A ₁₂ Q ^T ₁₂ − Q ₁₂ A ^T ₁₂ , A ₁₁ Q ₁₂ + Q ₁₂ A ^T ₂₂ = − V ₁₂ − A ₁₂ Q ₂₂ ,

A ₂₂ Q ^T ₁₂ + Q ^T ₁₂ A ^T ₁₁ = − V ^T ₁₂ − Q ₂₂ A ^T ₁₂ , A ₂₂ Q ₂₂ + Q ₂₂ A ^T ₂₂ = − V ₂₂

AUTOMATIC CONTROL

(15)

Systems with integrators + strictly stable poles 9(15)

The corresponding Lyapunov equation for this system reads

A ₁₁ A ₁₂ 0 A ₂₂

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

+ ^Q ¹¹ ^Q ¹² Q ^T ₁₂ Q ₂₂

A ^T ₁₁ 0 A ^T ₁₂ A ^T ₂₂

=− ^V ¹¹ ^V ¹² V ₁₂ ^T V ₂₂

,

which gives

A ₁₁ Q ₁₁ + Q ₁₁ A ^T ₁₁ = − V ₁₁ − A ₁₂ Q ^T ₁₂ − Q ₁₂ A ^T ₁₂ , A ₁₁ Q ₁₂ + Q ₁₂ A ^T ₂₂ = − V ₁₂ − A ₁₂ Q ₂₂ ,

A ₂₂ Q

^T

₁₂ + Q

^T

₁₂ A

^T

₁₁ = − V ₁₂

^T

− Q ₂₂ A

^T

₁₂ ,

A ₂₂ Q ₂₂ + Q ₂₂ A ^T ₂₂ = − V ₂₂

(16)

Systems with integrators + strictly stable poles 9(15)

The corresponding Lyapunov equation for this system reads

A ₁₁ A ₁₂ 0 A ₂₂

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

+ ^Q ¹¹ ^Q ¹² Q ^T ₁₂ Q ₂₂

A ^T ₁₁ 0 A ^T ₁₂ A ^T ₂₂

=− ^V ¹¹ ^V ¹² V ₁₂ ^T V ₂₂

,

which gives

A ₁₁ Q ₁₁ + Q ₁₁ A ^T ₁₁ = − V ₁₁ − A ₁₂ Q ^T ₁₂ − Q ₁₂ A ^T ₁₂ , A ₁₁ Q ₁₂ + Q ₁₂ A ^T ₂₂ = − V ₁₂ − A ₁₂ Q ₂₂ ,

A ₂₂ Q

^T

₁₂ + Q

^T

₁₂ A

^T

₁₁ = − V ₁₂

^T

− Q ₂₂ A

^T

₁₂ , A ₂₂ Q ₂₂ + Q ₂₂ A

^T

₂₂ = − V ₂₂

Proposed solution:

- contribution 2

1. Find Q ₂₂ by using the analytical solution.

2. Find Q ₁₂ by solving a Sylvester equation. 3. Find Q ₁₁ by solving a Lyapunov equation.

AUTOMATIC CONTROL

(17)

Systems with integrators + strictly stable poles 9(15)

The corresponding Lyapunov equation for this system reads

A ₁₁ A ₁₂ 0 A ₂₂

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

+ ^Q ¹¹ ^Q ¹² Q ^T ₁₂ Q ₂₂

A ^T ₁₁ 0 A ^T ₁₂ A ^T ₂₂

=− ^V ¹¹ ^V ¹² V ₁₂ ^T V ₂₂

,

which gives

A ₁₁ Q ₁₁ + Q ₁₁ A ^T ₁₁ = − V ₁₁ − A ₁₂ Q ^T ₁₂ − Q ₁₂ A ^T ₁₂ , A ₁₁ Q ₁₂ + Q ₁₂ A ^T ₂₂ = − V ₁₂ − A ₁₂ Q ₂₂ ,

A ₂₂ Q

^T

₁₂ + Q

^T

₁₂ A

^T

₁₁ = − V ₁₂

^T

− Q ₂₂ A

^T

₁₂ , A ₂₂ Q ₂₂ + Q ₂₂ A

^T

₂₂ = − V ₂₂

Proposed solution:

- contribution 2

1. Find Q ₂₂ by using the analytical solution.

2. Find Q ₁₂ by solving a Sylvester equation.

3. Find Q ₁₁ by solving a Lyapunov equation.

(18)

Systems with integrators + strictly stable poles 9(15)

The corresponding Lyapunov equation for this system reads

A ₁₁ A ₁₂ 0 A ₂₂

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

+ ^Q ¹¹ ^Q ¹² Q ^T ₁₂ Q ₂₂

A ^T ₁₁ 0 A ^T ₁₂ A ^T ₂₂

=− ^V ¹¹ ^V ¹² V ₁₂ ^T V ₂₂

,

which gives

A ₁₁ Q ₁₁ + Q ₁₁ A ^T ₁₁ = − V ₁₁ − A ₁₂ Q ^T ₁₂ − Q ₁₂ A ^T ₁₂ , A ₁₁ Q ₁₂ + Q ₁₂ A ^T ₂₂ = − V ₁₂ − A ₁₂ Q ₂₂ ,

A ₂₂ Q

^T

₁₂ + Q

^T

₁₂ A

^T

₁₁ = − V ₁₂

^T

− Q ₂₂ A

^T

₁₂ , A ₂₂ Q ₂₂ + Q ₂₂ A

^T

₂₂ = − V ₂₂

Proposed solution:

- contribution 2

1. Find Q ₂₂ by using the analytical solution.

2. Find Q ₁₂ by solving a Sylvester equation.

3. Find Q ₁₁ by solving a Lyapunov equation.

AUTOMATIC CONTROL

(19)

Systems with integrators + strictly stable poles 9(15)

The corresponding Lyapunov equation for this system reads

A ₁₁ A ₁₂ 0 A ₂₂

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

+

Q ₁₁ Q ₁₂ Q ^T ₁₂ Q ₂₂

A ^T ₁₁ 0 A ^T ₁₂ A ^T ₂₂

=− ^V ¹¹ ^V ¹² V ₁₂ ^T V ₂₂

,

which gives

A ₁₁ Q ₁₁ + Q ₁₁ A ^T ₁₁ = − V ₁₁ − A ₁₂ Q ^T ₁₂ − Q ₁₂ A ^T ₁₂ , A ₁₁ Q ₁₂ + Q ₁₂ A ^T ₂₂ = − V ₁₂ − A ₁₂ Q ₂₂ ,

A ₂₂ Q

^T

₁₂ + Q

^T

₁₂ A

^T

₁₁ = − V ₁₂

^T

− Q ₂₂ A

^T

₁₂ , A ₂₂ Q ₂₂ + Q ₂₂ A

^T

₂₂ = − V ₂₂

Proposed solution:

- contribution 2

1. Find Q ₂₂ by using the analytical solution.

2. Find Q ₁₂ by solving a Sylvester equation.

3. Find Q ₁₁ by solving a Lyapunov equation.

(20)

Van Loan’s method (standard method in literature) 10(15) 1. Form the 2n × 2n matrix

H = ^{A BSB}

T

0 − A ^T

.

2. Compute the matrix exponential

e ^HT

^s

= ^M ¹¹ ^M ¹² 0 M ₂₂

.

3. Then Q is given as Q = M ₁₂ M ^T ₁₁ .

Van Loan, C.F. (1978).Computing integrals involving the matrix exponential.

IEEE Transactions on Automatic Control, 23(3), 395-404.

Re Im

Figure : The poles of the system

AUTOMATIC CONTROL

(21)

Numerical Evaluation 11(15)

Marginally stable systems with 4 stable poles and 2

integrators are considered.

The 2D region

γ slow , γ fast ∈ [ 10 ⁻ ¹ , 10 ¹ ] is divided into 25 × 25 bins.

In total 100 systems are randomly generated for each bin.

Consider the sampling time T _s = 1 .

Re Im

γ slow

γ fast

Figure : The poles of the system

(22)

Numerical Evaluation - Proposed method 12(15)

10 ⁻¹ 10 ⁰ 10 ¹ 10 ⁻¹

10 ⁰ 10 ¹

γ fast

γ slo w

Proposed method

10 ⁻⁷ 10 ⁻⁶ 10 ⁻⁵

Numer ical error

If γ slow is small, the condition is closer to be violated

λ _i ( A ₁₁ ) + λ _j ( A ₁₁ ) 6= 0 ∀ i, j

Re Im

γ slow

γ fast

Figure : The poles of the system

AUTOMATIC CONTROL

(23)

Numerical Evaluation - Van Loan’s method 13(15)

10 ⁻¹ 10 ⁰ 10 ¹ 10 ⁻¹

10 ⁰ 10 ¹

γ _fast γ slo w

Van Loan’s method

10 ⁻⁷ 10 ⁻⁶ 10 ⁻⁵

Numer ical error

If γ fast is large, the matrix exponential is ill-conditioned

exp A BSB ^T 0 − A ^T

T s

Re Im

γ slow

γ fast

Figure : The poles of the system

(24)

Numerical Evaluation - Comparison 14(15)

10 ⁻¹ 10 ⁰ 10 ¹ 10 ⁻¹

10 ⁰ 10 ¹

γ fast

γ slo w

Proposed method

10 ⁻¹ 10 ⁰ 10 ¹ 10 ⁻¹

10 ⁰ 10 ¹

γ fast

γ slo w

Van Loan’s method

10 ⁻⁷ 10 ⁻⁶ 10 ⁻⁵

Numer ical error

The proposed method performs better if the slowest pole is fast.

Van Loan’s method performs better if the fastest pole is slow.

AUTOMATIC CONTROL

(25)

Conclusions 15(15)

An algorithm for computing an integral involving the matrix exponential common in optimal sampling has been proposed.

The algorithm is based on a Lyapunov equation and is justified with a novel theorem.

Numerical evaluations showed that the proposed algorithm has

advantageous numerical properties for slow sampling and

fast dynamics in comparison with a standard method in the

literature.

(26)

Discretizing stochastic dynamical systems using Lyapunov equations

Full text