**Discretizing stochastic dynamical systems** **using Lyapunov equations**

### Niklas Wahlström, Patrik Axelsson, Fredrik Gustafsson

### Division of Automatic Control Linköping University

### Linköping, Sweden

**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

**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}

**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}

^{Aτ}

^{T}

^{τ} *dτ*

^{τ}

AUTOMATIC CONTROL

**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}

^{Aτ}

^{T}

^{τ} *dτ* Problem

^{τ}

### How do we solve this integral in a numerically good manner

### for arbitrary A , B , S and T _{s} ?

**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

^{Aτ}

### e ^{Aτ} = I + *Aτ* + · · · + ^{1}

^{Aτ}

### p! A ^{p} ^{−} ^{1} *τ* ^{p} ^{−} ^{1} .

### Then, the integral has an analytical solution

### Q =

### Z _{T}

_{s}

### 0 e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ*

^{τ}

### =

### p − 1 i ∑ = 0

### p − 1 j ∑ = 0

### T ^{i} _{s} ^{+} ^{j} ^{+} ^{1} i!j! ( i + j + 1 ) ^{A}

### i S ( A ^{i} ) ^{T} .

### Re Im

### Figure : The poles of the system

AUTOMATIC CONTROL

**Theorem** 5(15)

### Theorem - Contribution 1 The solution to the integral

### Q =

### Z _{T}

_{s}

### 0 e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ* (1) satisfies the Lyapunov equation

^{τ}

### AQ + QA ^{T} = − BSB ^{T} + e ^{AT}

^{s}

### BSB ^{T} e ^{A}

^{T}

^{T}

^{s}

### (2)

**Theorem** 5(15)

### Theorem - Contribution 1 The solution to the integral

### Q =

### Z _{T}

_{s}

### 0 e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ* (1) satisfies the Lyapunov equation

^{τ}

### AQ + QA ^{T} = − BSB ^{T} + e ^{AT}

^{s}

### BSB ^{T} e ^{A}

^{T}

^{T}

^{s}

### (2)

### Idea

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

AUTOMATIC CONTROL

**Theorem** 5(15)

### Theorem - Contribution 1 The solution to the integral

### Q =

### Z _{T}

_{s}

### 0 e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ* (1) satisfies the Lyapunov equation

^{τ}

### AQ + QA ^{T} = − BSB ^{T} + e ^{AT}

^{s}

### BSB ^{T} e ^{A}

^{T}

^{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

**Theorem** 5(15)

### Theorem - Contribution 1 The solution to the integral

### Q =

### Z _{T}

_{s}

### 0 e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ* (1) satisfies the Lyapunov equation

^{τ}

### AQ + QA ^{T} = − BSB ^{T} + e ^{AT}

^{s}

### BSB ^{T} e ^{A}

^{T}

^{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

**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}

^{T}

^{s}

### | {z }

### − V

### ,

### which is identical to the solution of the integral

### Q =

### Z _{T}

_{s}

### 0 e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ.*

^{τ}

### Re Im

### Figure : The poles of the system

**Summary so far** 7(15)

### Re Im

### Q =

### p − 1 i ∑ = 0

### p − 1 j ∑ = 0

### T _{s} ^{i} ^{+} ^{j} ^{+} ^{1} i!j! ( i + j + 1 ) ^{A}

### i S ( A ^{i} ) ^{T}

### Re Im

### AQ + QA ^{T} = − V

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

AUTOMATIC CONTROL

**Systems with integrators + strictly stable poles** 8(15)

### Consider a system on the following block triangular form

### A =

### A _{11} A _{12} 0 A _{22}

### where all integrators are collected in A _{22} , i.e.

*λ* _{i} ( A _{11} ) 6= 0 ∀ i *λ* _{j} ( A _{22} ) = 0 ∀ j

### Re Im

### Figure : The poles of the system

**Systems with integrators + strictly stable poles** 9(15)

### The corresponding Lyapunov equation for this system reads

### A _{11} A _{12} 0 A _{22}

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### + ^{Q} ^{11} ^{Q} ^{12} Q ^{T} _{12} Q _{22}

### A ^{T} _{11} 0 A ^{T} _{12} A ^{T} _{22}

### =− ^{V} ^{11} ^{V} ^{12} V _{12} ^{T} V _{22}

### ,

### which gives

### A _{11} Q _{11} + Q _{11} A ^{T} _{11} = − V _{11} − A _{12} Q ^{T} _{12} − Q _{12} A ^{T} _{12} , A _{11} Q _{12} + Q _{12} A ^{T} _{22} = − V _{12} − A _{12} Q _{22} ,

### A _{22} Q ^{T} _{12} + Q ^{T} _{12} A ^{T} _{11} = − V ^{T} _{12} − Q _{22} A ^{T} _{12} , A _{22} Q _{22} + Q _{22} A ^{T} _{22} = − V _{22}

AUTOMATIC CONTROL

**Systems with integrators + strictly stable poles** 9(15)

### The corresponding Lyapunov equation for this system reads

### A _{11} A _{12} 0 A _{22}

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### + ^{Q} ^{11} ^{Q} ^{12} Q ^{T} _{12} Q _{22}

### A ^{T} _{11} 0 A ^{T} _{12} A ^{T} _{22}

### =− ^{V} ^{11} ^{V} ^{12} V _{12} ^{T} V _{22}

### ,

### which gives

### A _{11} Q _{11} + Q _{11} A ^{T} _{11} = − V _{11} − A _{12} Q ^{T} _{12} − Q _{12} A ^{T} _{12} , A _{11} Q _{12} + Q _{12} A ^{T} _{22} = − V _{12} − A _{12} Q _{22} ,

### A _{22} Q

^{T}

_{12} + Q

^{T}

_{12} A

^{T}

_{11} = − V _{12}

^{T}

### − Q _{22} A

^{T}

_{12} ,

### A _{22} Q _{22} + Q _{22} A ^{T} _{22} = − V _{22}

**Systems with integrators + strictly stable poles** 9(15)

### The corresponding Lyapunov equation for this system reads

### A _{11} A _{12} 0 A _{22}

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### + ^{Q} ^{11} ^{Q} ^{12} Q ^{T} _{12} Q _{22}

### A ^{T} _{11} 0 A ^{T} _{12} A ^{T} _{22}

### =− ^{V} ^{11} ^{V} ^{12} V _{12} ^{T} V _{22}

### ,

### which gives

### A _{11} Q _{11} + Q _{11} A ^{T} _{11} = − V _{11} − A _{12} Q ^{T} _{12} − Q _{12} A ^{T} _{12} , A _{11} Q _{12} + Q _{12} A ^{T} _{22} = − V _{12} − A _{12} Q _{22} ,

### A _{22} Q

^{T}

_{12} + Q

^{T}

_{12} A

^{T}

_{11} = − V _{12}

^{T}

### − Q _{22} A

^{T}

_{12} , A _{22} Q _{22} + Q _{22} A

^{T}

_{22} = − V _{22}

**Proposed solution:**

### - contribution 2

### 1. Find Q _{22} by using the analytical solution.

### 2. Find Q _{12} by solving a Sylvester equation. 3. Find Q _{11} by solving a Lyapunov equation.

AUTOMATIC CONTROL

**Systems with integrators + strictly stable poles** 9(15)

### The corresponding Lyapunov equation for this system reads

### A _{11} A _{12} 0 A _{22}

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### + ^{Q} ^{11} ^{Q} ^{12} Q ^{T} _{12} Q _{22}

### A ^{T} _{11} 0 A ^{T} _{12} A ^{T} _{22}

### =− ^{V} ^{11} ^{V} ^{12} V _{12} ^{T} V _{22}

### ,

### which gives

### A _{11} Q _{11} + Q _{11} A ^{T} _{11} = − V _{11} − A _{12} Q ^{T} _{12} − Q _{12} A ^{T} _{12} , A _{11} Q _{12} + Q _{12} A ^{T} _{22} = − V _{12} − A _{12} Q _{22} ,

### A _{22} Q

^{T}

_{12} + Q

^{T}

_{12} A

^{T}

_{11} = − V _{12}

^{T}

### − Q _{22} A

^{T}

_{12} , A _{22} Q _{22} + Q _{22} A

^{T}

_{22} = − V _{22}

**Proposed solution:**

### - contribution 2

### 1. Find Q _{22} by using the analytical solution.

### 2. Find Q _{12} by solving a Sylvester equation.

### 3. Find Q _{11} by solving a Lyapunov equation.

**Systems with integrators + strictly stable poles** 9(15)

### The corresponding Lyapunov equation for this system reads

### A _{11} A _{12} 0 A _{22}

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### + ^{Q} ^{11} ^{Q} ^{12} Q ^{T} _{12} Q _{22}

### A ^{T} _{11} 0 A ^{T} _{12} A ^{T} _{22}

### =− ^{V} ^{11} ^{V} ^{12} V _{12} ^{T} V _{22}

### ,

### which gives

### A _{11} Q _{11} + Q _{11} A ^{T} _{11} = − V _{11} − A _{12} Q ^{T} _{12} − Q _{12} A ^{T} _{12} , A _{11} Q _{12} + Q _{12} A ^{T} _{22} = − V _{12} − A _{12} Q _{22} ,

### A _{22} Q

^{T}

_{12} + Q

^{T}

_{12} A

^{T}

_{11} = − V _{12}

^{T}

### − Q _{22} A

^{T}

_{12} , A _{22} Q _{22} + Q _{22} A

^{T}

_{22} = − V _{22}

**Proposed solution:**

### - contribution 2

### 1. Find Q _{22} by using the analytical solution.

### 2. Find Q _{12} by solving a Sylvester equation.

### 3. Find Q _{11} by solving a Lyapunov equation.

AUTOMATIC CONTROL

**Systems with integrators + strictly stable poles** 9(15)

### The corresponding Lyapunov equation for this system reads

### A _{11} A _{12} 0 A _{22}

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### +

### Q _{11} Q _{12} Q ^{T} _{12} Q _{22}

### A ^{T} _{11} 0 A ^{T} _{12} A ^{T} _{22}

### =− ^{V} ^{11} ^{V} ^{12} V _{12} ^{T} V _{22}

### ,

### which gives

### A _{11} Q _{11} + Q _{11} A ^{T} _{11} = − V _{11} − A _{12} Q ^{T} _{12} − Q _{12} A ^{T} _{12} , A _{11} Q _{12} + Q _{12} A ^{T} _{22} = − V _{12} − A _{12} Q _{22} ,

### A _{22} Q

^{T}

_{12} + Q

^{T}

_{12} A

^{T}

_{11} = − V _{12}

^{T}

### − Q _{22} A

^{T}

_{12} , A _{22} Q _{22} + Q _{22} A

^{T}

_{22} = − V _{22}

**Proposed solution:**

### - contribution 2

### 1. Find Q _{22} by using the analytical solution.

### 2. Find Q _{12} by solving a Sylvester equation.

### 3. Find Q _{11} by solving a Lyapunov equation.

**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} ^{11} ^{M} ^{12} 0 M _{22}

### .

### 3. Then Q is given as Q = M _{12} M ^{T} _{11} .

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

**Numerical Evaluation** 11(15)

### Marginally stable systems with 4 stable poles and 2

### integrators are considered.

### The 2D region

*γ* slow *, γ* fast ∈ [ 10 ^{−} ^{1} , 10 ^{1} ] 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

**Numerical Evaluation - Proposed method** 12(15)

### 10 ^{−1} 10 ^{0} 10 ^{1} 10 ^{−1}

### 10 ^{0} 10 ^{1}

*γ* fast

*γ* slo w

### Proposed method

### 10 ^{−7} 10 ^{−6} 10 ^{−5}

### Numer ical error

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

*λ* _{i} ( A _{11} ) + *λ* _{j} ( A _{11} ) 6= 0 ∀ i, j

### Re Im

*γ* slow

*γ* fast

### Figure : The poles of the system

AUTOMATIC CONTROL

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

### 10 ^{−1} 10 ^{0} 10 ^{1} 10 ^{−1}

### 10 ^{0} 10 ^{1}

*γ* _{fast} *γ* slo w

### Van Loan’s method

### 10 ^{−7} 10 ^{−6} 10 ^{−5}

### 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

**Numerical Evaluation - Comparison** 14(15)

### 10 ^{−1} 10 ^{0} 10 ^{1} 10 ^{−1}

### 10 ^{0} 10 ^{1}

*γ* fast

*γ* slo w

### Proposed method

### 10 ^{−1} 10 ^{0} 10 ^{1} 10 ^{−1}

### 10 ^{0} 10 ^{1}

*γ* fast

*γ* slo w

### Van Loan’s method

### 10 ^{−7} 10 ^{−6} 10 ^{−5}

### 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

**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.

**Analytical solution - Example** 16(15) Consider a **double integrator (two zero eigenvalues of** A )

### ˙x ( t ) = ^{0} ^{1} 0 0

### | {z }

### A

### x ( t ) + ^{0} 1

### |{z}

### B

### w ( t ) , E [ w ( t ) w ( *τ* )] = *γδ* ( t − *τ* ) .

### Only zero eigenvalues ⇔ nilpotent, so

### e ^{Aτ} = I + *Aτ* = ^{1} ^{0} 0 1

^{Aτ}

### + ^{0} ^{1} 0 0

### *τ*

### Q can be computed analytically

### Q =

### Z _{T}

_{s}

### 0

### e ^{Aτ} BSB ^{T} e ^{A}

^{Aτ}

^{T}

^{τ} *dτ* = *γ*

^{τ}

### " _{T}

3
s
### 3 T

^{2}

_{s}

### T

^{2}

_{s}

### 2

### 2 T _{s}

### # .

AUTOMATIC CONTROL

**Discretization** 17(15)

### Continuous-time model

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

### ...or more strict

### dx ( t ) = Ax ( t ) + Bdw ( t ) E [ _{dw} ( _{t} ) _{dw} ( _{t} )

^{T}