Kalman Filter là một mô hình Linear-Gaussian State Space Model thuộc nhóm thuật toán dự đoán chuỗi thời gian. Thuật toán được lấy tên theo Rudolf E. Kálmán, một nhà khoa học ảnh hưởng quan trọng trong quá trình phát triển thuật toán.

Nếu là một kỹ sư điều khiển hệ thống, bạn hiểu rằng điều khiển hệ thống không đơn thuần nghĩa là đo rồi điều khiển những gì mình đo được. Đo lường chính xác nhường nào, hệ thống điều khiển càng hoạt động ổn định nhường nấy.

Hãy quay ngược thời gian về những năm 1960 với những khó khăn của những kỹ sư thiết kế phi thuyền Apollo trong xử lý tín hiệu. Dữ liệu thô đo được từ máy tính, được đo từ các cảm biến như con quay hồi chuyển, gia tốc kế, dữ liệu rađa vốn dĩ đầy nhiễu, đầy lỗi ngẫu nhiên và đặc biệt lộn xộn không chính xác. Khi lao về mặt trăng ở tốc độ cao, ai mà chắc chắn rằng điều gì sẽ xảy ra với những gì bạn điều khiển, bạn sẽ không muốn Apollo lao thẳng vào nhà mình đâu!

Kalman Filter là một công cụ mạnh mẽ kết hợp thông tin không chắc chắn ở thời điểm hiện tại cùng với thông tin đầy nhiễu loạn của môi trường sang một dạng thông tin mới đáng tin cậy hơn để phục vụ dự đoán tương lai. Điểm mạnh của Kalman Filter là chạy rất nhanh và tính ổn định cao.

Một trong những ứng dụng đầu tiên của Kalman Filter là được áp dụng vào phi thuyền Apollo:

The Apollo computer used 2k of magnetic core RAM and 36k wire rope […]. Clock speed was under 100 kHz […]. The fact that the MIT engineers were able to pack such good software (one of the very first applications of the Kalman filter) into such a tiny computer is truly remarkable.
— Phỏng vấn Jack Crenshaw bởi Matthew Reed, TRS-80.org (2009)

1. Ý tưởng cơ bản thuật toán

Thuật toán có rất nhiều ứng dụng vào các bài toán khác nhau. Tuy nhiên để dể hình dung, bây giờ chúng ta hãy tưởng tượng mình là một nhà sáng chế phi thuyền Thetalo (một phiên bản lỗi từ Apollo).

Phi thuyền sẽ gửi về thông tin trạng thái sau mỗi khoảng thời gian nhất định thông qua bộ cảm biến, nhiệm vụ của chúng ta là dự đoán trạng thái trên phi thuyền để phục vụ tác vụ điều khiển hệ thống.

Gọi $\mathbf{x}_{t} = \left[ \begin{array}{*{20}{c}} \mathtt{position} \\ \mathtt{velocity} \end{array} \right] \in \mathbb{R}^{2}$ là trạng thái vị trí và vận tốc hiện tại của phi thuyền. Tại mỗi thời điểm $t$, bạn không hề có thông tin chính xác về $\mathbf{x}_{t}$.

Vì sao vậy? Hệ mà chúng ta điều khiển không nằm trong một môi trường lý tưởng, nhiều tác động từ môi trường xung quanh mà chúng ta không đoán trước ở sẽ ảnh hưởng những tính toán ban đầu.

Biết đâu được một cơn gió vô tình, một cơn mưa ngang qua, một trận sấm sét cuồng phong sẽ làm mọi tính toán sai lầm. Trực cảm cho ta thấy nó đang ở đâu đó trong không gian này, tại một số vùng nào đó khả năng xảy ra cao hơn, tại một số vùng nào đó khả năng xảy ra thấp hơn. Phân bố xác suất như một bản đồ về sự chắc chắn, một niềm tin về những sự kiện xảy ra.

Lựa chọn phân bố xác suất phù hợp là điều cần thiết như thể lựa chọn cho mình một niềm tin, trong mô hình Kalman Filter phân bố xác suất được chọn là phân bố chuẩn nhiều chiều. Sở dĩ chọn phân bố chuẩn nhiều chiều là vì đây là phân bố liên tục phù hợp bài toán đang xét, một vài tính chất đặc biệt của phân bố chuẩn nhiều chiều thuận lợi cho tác vụ ước lượng tham số.

Tại mỗi thời điểm $t$, giả định vector ngẫu nhiên $x_{t}$ có phân bố chuẩn nhiều chiều hay $x_{t} \sim \mathcal{N}(\hat{x}_{t}, \Sigma_{t} )$ với $\hat{x}_{t}$ là kỳ vọng của $x_{t}$, $\Sigma_{t}$ là ma trận hiệp phương sai với mỗi phần tử tại hàng $i$ cột $j$ thể hiện độ biến thiên cùng nhau của hai biến ngẫu nhiên thứ $i$ và $j$. Công thức tổng quát của hàm mật độ xác suất vector ngẫu nhiên $x_{t}$ với kỳ vọng $\hat{x}_{t}$, ma trận hiệp phương sai $\Sigma_{t}$ với số chiều $D$ (trong ví dụ đang xét, $D=2$) là:

$$p(x_t) = \frac{1}{\left( 2 \pi \right)^{D/2} \left| \Sigma_{t} \right|^{1 / 2} } \text{exp}\left[ -\frac{1}{2} (x_t - \hat{x}_{t})^{\intercal} \Sigma_{t}^{-1} (x_t - \hat{x}_{t}) \right]$$

Lập luận theo cách trên lại có một ưu thế nữa, trong nhiều bài toán khoa học kỹ thuật, đôi lúc bạn cần biết xác suất xảy ra trong một vùng nào đấy, chẳng hạn:

Đôi lúc chúng ta chỉ quan tâm đến một vùng nào đó của một đại lượng, chẳng hạn từ $95$ đến $110$ độ C thì CPU máy tính phải tắt vì nhiệt độ vượt ngưỡng an toàn, xác suất để nhiệt độ rơi vào vùng này là bao nhiêu nhỉ? Tương tự với một tập đại lượng, chúng ta quan tâm đến một vùng giá trị nào đó thật sự quan trọng trong bài toán đang xét.
Với hàm mật độ xác suất liên tục $p(x_t)$ chúng ta có thể dể dàng tính xác suất xảy ra trong vùng $E$ bằng tích phân $\int_{E} p(x_t) dx_t$, thật tiện lợi đúng không nào!

Giả định rằng tại thời điểm $t$ phi thuyền có:

$$\hat{x}_{t} = \left[ \begin{array}{*{20}{c}} -0.2 \\ 0.1 \end{array} \right]$$

$$ \Sigma_{t} = \left[ \begin{array}{*{20}{c}} 0.4 & 0.35 \\ 0.35 & 0.6 \end{array} \right]$$

Thông tin mà chúng ta có không chỉ có những gì chúng ta tiên lượng về nó (cách chúng ta đoán về $x_t$). Trên phi thuyền được gắn một cảm biến đặc biệt cho phép đo trạng thái hiện tại và gửi thông tin về trụ sở.

Tuy nhiên, mặc dù vậy, cảm biến này lại không chính xác, bộ cảm biến cho chúng ta biết trạng thái hiện tại là $y_{t} = \left[ \begin{array}{*{20}{c}} 2 \\ -2 \end{array} \right]$

Trong trường hợp này, chúng ta giả định rằng: $$y_{t} = G_{t} x_{t} + v_{t} $$ với $v_{t}$ là một vector ngẫu nhiên nhiễu có phân bố chuẩn nhiều chiều $v_{t} \sim \mathcal{N}(0,R_{t})$, chúng ta gọi $v_{t}$ là nhiễu mô hình quan sát (observation noise).

Giả định trên có thể hiểu $y_{t}$ được xem như là “nhiễu quan sát được từ $x_t$“, $y_t$ là kết quả của tổ hợp tuyến tính từng thành phần $x_{t}$ dựa trên ma trận $G_t$ và bị ảnh hưởng cộng thêm một lượng nhiễu môi trường với vector ngẫu nhiên $v_{t}$. Người ta gọi công thức trên là mô hình quan sát, nó thể hiện quá trình $x_t$ biến đổi thành $y_t$ đồng thời bị ảnh hưởng nhiễu từ môi trường.

Giờ đây, chúng ta có niềm tin của chúng ta về trạng thái phi thuyền qua hàm mật độ xác suất $p(x_t)$ và một sự kiện đã xảy ra là thông tin thu được từ bộ cảm biến - trạng thái $y_t$, liệu rằng chúng ta có thể kết hợp cả hai thông tin đang có thành một thông tin mới có ý nghĩa hơn, giúp chúng ta hiểu hơn về trạng thái hiện tại hệ thống hay không?

Biết $y_t$ sẽ giúp chúng ta cập nhật niềm tin về $x_t$ như thế nào? Nếu bạn đang tìm một giải pháp như trên, có lẽ định lý Bayes là một câu trả lời chúng ta đang tìm:

$$p(x_t | y_t) = \frac{p(y_t | x_t) p(x_t)}{p(y_t)}$$

$p(x_t)$ hàm mật độ xác suất tiên nghiệm, niềm tin của chúng ta không phụ thuộc vào sự kiện $y_t$ xảy ra.
$p(y_t | x_t)$ hàm mật độ xác suất khả dĩ, mật độ xác suất có điều kiện khi biết trạng thái $x_t$ xảy ra, vì chúng ta biết rằng $y_t$ là một vector ngẫu nhiên sao cho $y_{t} = G_{t} x_{t} + v_{t} $, khi biết $x_t$ xảy ra nghĩa là $G_t x_t$ là một vector hằng, vector ngẫu nhiên $y_t$ được biểu diễn bằng một vector hằng cộng thêm một vector ngẫu nhiên có phân bố chuẩn nhiều chiều $v_{t} \sim \mathcal{N}(0,R_{t})$, hay nói cách khác $y_t | x_t$ hiện tại là một phân bố chuẩn nhiều chiều. Phân bố $y_t | x_t \sim \mathcal{N}(G_t x_t ,R_{t})$.
$p(y_t)$ hàm mật độ xác suất biên không phụ thuộc vào $x_t$ đóng vai trò như một hằng số chuẩn hóa.

Phân bố $x_t | y_t$ được tìm như thế nào? Trên nền tảng Linear Gaussian Systems (LGS - mô hình phân bố chuẩn tuyến tính), lời giải của định lý Bayes cho bài toán này như sau:

Định lý Bayes cho Linear Gaussian Systems
Giả định có hai vector ngẫu nhiên $m \in \mathbb{R}^{D_m}$ và $n \in \mathbb{R}^{D_n}$ với $m$ là vectơ ẩn, $n$ là vector nhiễu quan sát được từ $m$ với: $ m \sim \mathcal{N}(\mu_m , \Sigma_m )$ $ n \| m \sim \mathcal{N}(C m + d, \Sigma_n )$ Phân bố hậu nghiệm $m \| n$ khi biết $n$ là $m \| n \sim \mathcal{N}(\mu_{m \| y}, \Sigma_{m \| n})$ với: $ \Sigma_{m \| n} = \left( \Sigma_m^{-1} + C^{\intercal} \Sigma_n^{-1} C \right)^{-1}$ $ \mu_{m \| n} = \Sigma_{m \| n} \left[ C^{\intercal}\Sigma_n^{-1}(n-d) + \Sigma_m^{-1} \mu_m \right] $ TL;DR: phần chứng minh định lý sẽ không được bàn ở đây, trong một dịp nào đó nếu có thể mình sẽ viết về các mô hình phân bố chuẩn nhiều chiều, bạn đọc quan tâm có thể đọc thêm tài liệu đính kèm. (Kevin P. Murphy phần 6 tài liệu [7], Chris Bracegirdle tài liệu [6] phần hệ luận Corollary 5)

Định lý Bayes cho Linear Gaussian Systems

Giả định có hai vector ngẫu nhiên $m \in \mathbb{R}^{D_m}$ và $n \in \mathbb{R}^{D_n}$ với $m$ là vectơ ẩn, $n$ là vector nhiễu quan sát được từ $m$ với:

$ m \sim \mathcal{N}(\mu_m , \Sigma_m )$
$ n | m \sim \mathcal{N}(C m + d, \Sigma_n )$

Phân bố hậu nghiệm $m | n$ khi biết $n$ là $m | n \sim \mathcal{N}(\mu_{m | y}, \Sigma_{m | n})$ với:

$ \Sigma_{m | n} = \left( \Sigma_m^{-1} + C^{\intercal} \Sigma_n^{-1} C \right)^{-1}$
$ \mu_{m | n} = \Sigma_{m | n} \left[ C^{\intercal}\Sigma_n^{-1}(n-d) + \Sigma_m^{-1} \mu_m \right] $

TL;DR: phần chứng minh định lý sẽ không được bàn ở đây, trong một dịp nào đó nếu có thể mình sẽ viết về các mô hình phân bố chuẩn nhiều chiều, bạn đọc quan tâm có thể đọc thêm tài liệu đính kèm. (Kevin P. Murphy phần 6 tài liệu [7], Chris Bracegirdle tài liệu [6] phần hệ luận Corollary 5)

Tuy nhiên nếu tính toán trực tiếp ma trận hiệp phương sai $\Sigma_{x_t | y_t}$ như trên thì có lẽ bạn vừa bỏ qua một vài “điều thú vị” trong đại số tuyến tính khiến giải pháp của chúng ta thú vị hơn!

Đồng nhất thức ma trận Woodbury (hay còn gọi matrix inversion lemma)
Cho $4$ ma trận $A$ $(n \times n)$, $U$ $(n \times k)$, $C$ $(k \times k)$, $V$ $(k \times n)$ Đồng nhất thức Woodbury nói rằng: $$ \left( A + UCV \right)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}$$Nghĩa là khi bạn có nhu cầu tính toán công thức bên trái, bạn có thể xử dụng công thức bên phải để thay thế! Nó có lợi điểm gì trong bài toán này nhỉ? TL;DR: phần chứng minh này bạn có thể xem thêm tài liệu [8].

Đồng nhất thức ma trận Woodbury (hay còn gọi matrix inversion lemma)

Cho $4$ ma trận $A$ $(n \times n)$, $U$ $(n \times k)$, $C$ $(k \times k)$, $V$ $(k \times n)$
Đồng nhất thức Woodbury nói rằng: $$ \left( A + UCV \right)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}$$Nghĩa là khi bạn có nhu cầu tính toán công thức bên trái, bạn có thể xử dụng công thức bên phải để thay thế! Nó có lợi điểm gì trong bài toán này nhỉ?

TL;DR: phần chứng minh này bạn có thể xem thêm tài liệu [8].

Trước khi giải thích ý nghĩa của đồng nhất thức ma trận Woodbury trong Kalman Filter, chúng ta hãy cùng nhau viết lại lời giải phân bố hậu nghiệm $x_t | y_t$ vừa tìm được, chúng ta gọi đây là “phân bố lọc” (Filtering Distribution):

$$ x^{\mathtt{F}}_t \sim x_t | y_t \sim \mathcal{N}(\hat{x}^{\mathtt{F}}, \Sigma^{\mathtt{F}})$$

Với tham số kỳ vọng và ma trận hiệp phương sai của phân bố chuẩn nhiều chiều:

$\color{WildStrawberry}{\Sigma^{\mathtt{F}}_t = \left( \Sigma_t^{-1} + G_t^{\intercal} R_t^{-1} G_t \right)^{-1} = \Sigma_t - \Sigma_t G_t^{\intercal} \left( R + G_t \Sigma_t G_t^{\intercal} \right)^{-1} G_t \Sigma_t}$
$\color{WildStrawberry}{\hat{x}^{\mathtt{F}}_t = \hat{x}_t + \Sigma_t G_t^{\intercal} \left( G_t \Sigma_t G_t^{\intercal} + R_t \right)^{-1} (y_t - G_t \hat{x}_t)}$

Ý nghĩa của đồng nhất thức ma trận Woodbury: Trong nhiều ứng dụng theo thời gian thực, thường thì chỉ có một số ít trạng bạn nhận được từ bộ cảm biến tại một thời điểm, hay nói cách khác số chiều của vector ngẫu nhiên $y_t$ rất nhỏ, lúc này chi phí tìm ma trận nghịch đảo của $\left( R_t + G_t \Sigma_t G_t^{\intercal} \right)^{-1}$ tương đối nhỏ, đồng thời các ma trận khác đã có rồi, việc tính toán còn lại chỉ là nhân ma trận sẽ nhanh hơn rất nhiều so với bài toán gốc nhân rồi tính một ma trận nghịch đảo rất lớn!

Giả định tiếp tục với ví dụ trên với: $$ R_{t} = \left[ \begin{array}{*{20}{c}} 0.26 & 0.2275 \\ 0.2275 & 0.39 \end{array} \right] $$

$$ G_t = \left[ \begin{array}{*{20}{c}} 1 & 0 \\ 0 & 1 \end{array} \right] $$

Cuối cùng chúng ta cũng kết hợp được hai thông tin thành một thông tin mới hữu ích hơn để hiểu rõ về trạng thái của phi thuyển!

Ngay bây giờ chúng ta muốn dự đoán tương lai, trạng thái của phi thuyền sẽ đi đâu về đâu?

Trạng thái mới có thể biểu diễn dưới dạng trạng thái dịch chuyển và bị thêm một ít nhiễu từ trạng thái cũ:

$$x_{t+1} = A_t x_t + B_t u_t + w_t$$

Với:

$A_t$: ma trận mô hình chuyển đổi trạng thái.
$B_t$: ma trận mô hình điều khiễn đầu vào (dùng để kết hợp với $u_t$ từ người dùng).
$u_t$: vector điều khiển (người dùng nhập để kết hợp với $B_t$).
$w_t$: vector ngẫu nhiên nhiễu hệ thống (system noise), có phân bố chuẩn nhiều chiều $w \sim \mathcal{N}(0, Q_t)$.

Thường thì $B_t$ và $u_t$ bị khuyết, $B_t u_t$ được xem là sửa lỗi hệ thống được thêm bởi người dùng, việc thêm vào cũng giống như việc thêm một đại lượng không đổi vào các công thức, việc biến đổi với các tính chất kỳ vọng và hiệp phương sai không khó. Nên để gọn phần trình bày dưới đây chúng ta sẽ xem như $B_t u_t$ bị khuyết.

Vì vậy chúng ta xem $x_{t+1}$ như là:

$$x_{t+1} = A_t x_t + w_t$$

Tuy nhiên chúng ta lại trực nhớ ra rằng, xác suất tiên nghiệm của chúng ta đã được thay thế bằng một xác suất hậu nghiệm nhìn có vẻ hợp lí hơn, do đó chúng ta biểu diễn $x_{t+1}$ dưới dạng:

$$x_{t+1} = A_t x_t^{\mathtt{F}} + w_t$$

Vì tổ hợp tuyến tính của các phân bố chuẩn là phân bố chuẩn, nên $x_{t+1}$ cũng là phân bố chuẩn, mà là phân bố chuẩn nhiều chiều thì chúng ta cần tìm hai tham số, kỳ vọng và ma trận hiệp phương sai.

May mắn thay, nhờ vào tính chất kỳ vọng và ma trận hiệp phương sai với vector ngẫu nhiên, ta có:

$$\hat{x}_{t+1} = \mathbb{E} \left[ A_t x_t^{\mathtt{F}} + w_t\right] = A_t \mathbb{E} x_t^{\mathtt{F}} + \mathbb{E} w_t = A_t \hat{x}_t^{\mathtt{F}} $$ $$\Sigma_{t+1} = \text{Var} \left[ A_t x_t^{\mathtt{F}} + w_t \right] = A_t \text{Var} \left[ x_t^{\mathtt{F}} \right] A_t^{\intercal} + Q_t = A_t \Sigma_t^{\mathtt{F}} A_t^{\intercal} + Q_t $$

Đặt $\color{WildStrawberry} K_t = \Sigma_t G_t^{\intercal} (G_t \Sigma_t G_t^{\intercal} + R_t)^{-1}$ hệ số này được gọi là Kalman Gain, chúng ta có thể viết phân bố mới ở dạng gọn hơn: $$ x_{t + 1} \sim \mathcal{N}(\hat{x}_{t+1}, \Sigma_{t+1})$$

Viết gọn lại:

$\color{WildStrawberry} \hat{x}_{t+1} = A_t \color{Violet} \hat{x}_t^{\mathtt{F}} \color{WildStrawberry} = A_t \color{Violet} \left( \hat{x}_t + K_t (y_t - G_t \hat{x}_t ) \right)$
$\color{WildStrawberry} \Sigma_{t+1} = A_t \color{Violet} \Sigma_t^{\mathtt{F}} \color{WildStrawberry} A_t^{\intercal} + Q_t = A_t \color{Violet} (\mathbf{I} - K_t G_t) \Sigma_t \color{WildStrawberry} A_t^{\intercal} + Q_t $

(với $\mathbf{I}$ là ma trận đơn vị)

Giả định rằng ta có:

$$ A_{t} = \left[ \begin{array}{*{20}{c}} 1.25 & 0 \\ 0 & -0.4 \end{array} \right] $$

$$ Q_t = \left[ \begin{array}{*{20}{c}} 0.12 & 0.105 \\ 0.105 & 0.18 \end{array} \right] $$

Chúng ta có phân bố dự đoán tương lai như sau:

2. Tổng kết và thủ tục đệ quy:

Kalman Filter sẽ thực hiện thủ tục đệ quy như gồm 2 thủ tục “đo lường” và “dự đoán”:

Thủ tục cập nhật đo lường (Measurement):
1. Tính Kalman Gain:
  $\color{WildStrawberry} K_t = \Sigma_t G_t^{\intercal} (G_t \Sigma_t G_t^{\intercal} + R_t)^{-1}$
2. Cập nhật trạng thái hiện tại thành phân bố lọc:
  $ \color{WildStrawberry} \hat{x}_t := \hat{x}_t^{\mathtt{F}} $ với $ \hat{x}_t^{\mathtt{F}} = \hat{x}_t + K_t (y_t - G_t \hat{x}_t ) $
  $\color{WildStrawberry} \Sigma_t := \Sigma_t^{\mathtt{F}} $ với $\Sigma_t^{\mathtt{F}} = (\mathbf{I} - K_t G_t) \Sigma_t $
  (với $\mathbf{I}$ là ma trận đơn vị)
Dự đoán tương lai (Prediction):
1. Cập nhật trạng thái tương lai:
  $\color{WildStrawberry} \hat{x}_{t+1} = A_t \hat{x}_t + B_t u_t$
  (trong các hệ điều khiển đơn giản $B_t u_t$ khuyết)
2. Cập nhập ma trận hiệp phương sai:
  $\color{WildStrawberry} \Sigma_{t+1} = A_t \Sigma_t A_t^{\intercal} + Q_t$

3. Cài đặt thuật toán và những vấn đề tính toán

Chi phí tính toán của thuật toán chủ yếu đến từ việc nhân và lấy ma trận nghịch đảo, như đã trình bày ở trên nếu như vector ngẫu nhiên quan sát được $y_t$ có số chiều nhỏ hơn nhiều so với $x_t$ chúng ta có thể tính toán bằng khai triển ma trận nghịch đảo như trên để tính toán hiệu quả.

Việc thiết lập ma trận $A_t$, $Q_t$, $G_t$, $R_t$ thường là do kỹ sư thiết lập dựa trên hiểu biết về dữ liệu của mình (các đại lượng vật lý, kinh tế,… thường có mối quan hệ mật thiết với nhau). Tuy nhiên có rất nhiều phương pháp học tham số cho mô hình Kalman Filter mà bạn đọc có thể tự tham khảo thêm (trong cộng đồng lý thuyết điều khiển thì vấn đề này được biết đến như systems identification, thường cách đơn giản nhất là bạn sẽ thu thập dữ liệu, xây dựng hàm mất mát và giải bài toán bình phương tối tiểu).

Nếu chú ý kĩ bạn sẽ thấy ma trận hiệp phương sai $\Sigma_{t+1}$ chỉ phụ thuộc vào $\Sigma_t$ không phụ thuộc vào sự kiện quan sát được trả về từ bộ cảm biến $y_t$, trong một số trường hợp đặc biệt, chúng ta có thể tính toán ma trận này trước khi $y_t$ xuất hiện, việc này đòi hỏi bạn cần hiểu về sự hội tụ của ma trận $\Sigma_t$ xảy ra như thế nào (bạn đọc quan tâm tìm hiểu thêm về Ricatti equations).

Trong những ứng dụng thương mại lớn, Kalman Filter được thiết kế “tương đối phức tạp” vì tính ổn định tính toán số, thời gian chạy,… Nếu quan tâm đến những vấn đề này, bạn đọc có thể truy cập vào kho tài liệu đồ sộ về Kalman Filter ở đây: http://www.cs.unc.edu/~welch/kalman/

Hiện tại có rất nhiều thư viện cài đặt Kalman Filter sẳn như OpenCV (C/C++, Python), filterpy (Python), Matlab.

Dưới đây là một bài mô phỏng đơn giản trên phân bố chuẩn 1 chiều dựa trên ví dụ của Greg Welch và Gary Bishop, giả sử có một nguồn điện với hiệu điện thế không thay đổi theo thời gian (cho lý tưởng là vậy nhé, nghĩa là $x = 0.28$ volts mãi luôn). Bạn thực sự không biết giá trị hiệu điện thế của nguồn điện này đâu, bạn quyết định đi tìm giá trị này như một thú vui lúc rảnh rỗi. May mắn thay, bạn tìm được một cái vôn kế đời cũ không chính xác. Rất may là nhà sản xuất có để lại thông tin, vôn kế này bị nhiễu theo phân bố $\mathcal{N}(0, 0.01)$. Biết rằng mỗi lần đo kế tiếp, hiệu điện thế không đổi gì nhiều (thiết lập ma trận chuyển trạng thái $A$), mô hình quan sát của được đo từ vôn kế bằng giá trị thực tế cộng thêm nhiễu (giả định ma trận $G$). Không rõ độ lệch chuẩn của nguồn điện phóng ra khi sử dụng bình thường là bao nhiêu, bạn tạm đoán là $1$ volt (thiết lập ma trận $\Sigma$), tại thời điểm ban đầu chưa biết bạn đoán nguồn điện $0.5$ volt (thiết lập $\hat{x}$).

Trong ví dụ này, mình sẽ cài đặt phần mô phỏng nhìn hơi dỏm một tí (vì đây là phân bố chuẩn 1 chiều không cần phải cài đặt rườm rà như vậy, hơn nữa việc giả lập $y$ có thể tạo luôn một lần thay vì gọi từng lần). Tuy nhiên mình mong muốn mã nguồn này có thể phản ánh được nội dung trong bài viết, vì thế chúng ta cùng thử nhé:

#   _______________________________________________________________________
#  |[] ThetaLogShell                                                 |X]|!"|
#  |"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""|"|
#  | > ThetaLog.com                                                      | |
#  | > ENV: Python 3.7.0 Anaconda                                        | |
#  | > Kalman Filter (ví dụ đơn giản)                                    | |
#  |_____________________________________________________________________|/

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import style

style.use('bmh')

def simulate_volts(x, G, R):
    """
    Hàm lấy giá trị quan sát y từ mô hình quan sát
    y = Gx + v
    (với v là vector ngẫu nhiên v ~ N(0,R))

    :param x: numpy.array - vector thực tế x
    :param G: ma trận mô hình quan sát
    :param R: ma trận hiệp phương sai nhiễu quan sát

    :return y: numpy.array - vector ngẫu nhiên quan sát được
    """
    return G @ x + np.random.multivariate_normal(np.zeros(shape=R.shape[0]), R)

# Số lượng vòng lặp
n_iteration = 100

# Thật sự là -0.28 vols
x = np.array([[-0.28]])

# Ma trận chuyển đổi trạng thái
A = np.array([[1]])
# Ma trận hiệp phương sai nhiễu hệ thống
Q = np.array([[1e-6]]) 
# Ma trận mô hình quan sát
G = np.array([[1]])
# Ma trận hiệp phương sai nhiễu quan sát
R = np.array([[0.01]])

list_of_x = []
list_of_y = []

# Thiết lập phân bố ban đầu
x_hat = np.array([[0.5]])
sigma = np.array([[1.0]])

for i in range(n_iteration):
    # Bước dự đoán (Prediction)
    x_hat = A @ x_hat
    sigma = A @ sigma @ A.T + Q
    print("Sigma: %r " % sigma)
    list_of_x.append(float(x_hat))
    
    # Bước đo lường (Measurement)
    y = simulate_volts(x, G, R)
    kalman_gain = sigma @ G.T @ np.linalg.inv(G @ sigma @ G.T + R)
    x_hat = x_hat + kalman_gain @ (y - G @ x_hat)
    sigma = (np.eye(len(sigma)) - kalman_gain @ G) @ sigma
    list_of_y.append(float(y))

# Vẽ nào
plt.plot(list_of_y,'k+',label='dữ liệu quan sát')
plt.axhline(x, color='b',label='giá trị thực tế')
plt.plot(list_of_x,'r-',label='giá trị dự đoán')
plt.legend()

plt.xlabel('Thời gian $t$')
plt.ylabel('Hiệu điện thế (volts)')
plt.title('Kalman Filter')

plt.show()

Kết quả (quan sát $\Sigma_t$ từ Terminal bạn có thấy gì đặc biệt không nào?):

4. Phần mở rộng và ứng dụng:

Kalman Filter phiên bản cơ bản bị giới hạn bởi giả định mô hình tuyến tính, dưới đây là hai mô hình phi tuyến quan trọng của Kalman Filter:

Extended Kalman filter
Unscented Kalman filter (bản nâng cấp từ EKF, tốt hơn)

Ứng dụng Kalman Filter rất rộng, thuật toán 60 năm tuổi nhưng vẫn tiếp tục ảnh hưởng đến cuộc sống của chúng ta hàng ngày:

Ứng dụng vào các bộ GPS thương mại. IEEE
Xe tự láy (Advanced Driver Assistance Systems) Why should car drivers love Kalman filtering?
Thực tế ảo Human Motion - MIT Press
Robotic SLAM (simultaneous localization and mapping)
Dự báo thời tiết 2016 Weather forecasting
Theo dõi đối tượng trong thị giác máy tính Object Tracking - MathWorks
Kinh tế lượng Advances in Econometrics - Cambridge

5. Tổng kết

Kalman Filter là một thuật toán dự đoán chuỗi thời gian. Ý tưởng cơ bản thuật toán dựa trên định lý Bayes cho phân bố chuẩn nhiều chiều. Thuật toán có nhiều ứng dụng quan trọng trong xử lý tín hiệu, kinh tế lượng,…

Bài viết gói gọn trong Kalman Filter cơ bản, các phần mở rộng EKF, UKF được ứng dụng rộng rãi hơn trong các ứng dụng phức tạp.

Chúc bạn có nhiều trải nghiệm thú vị với Kalman Filter!

Tham khảo

Greg Welch, Gary Bishop. An Introduction to the Kalman Filter. https://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdf
Thomas J. Sargent John Stachurski. A First Look at the Kalman Filter. Quantitative Economics. https://lectures.quantecon.org/jl/kalman.html
Bzarg. How a Kalman filter works, in pictures. https://www.bzarg.com/p/how-a-kalman-filter-works-in-pictures/
Matlab. Understanding Kalman Filters, Part 1: Why Use Kalman Filters?. https://www.youtube.com/watch?v=mwn8xhgNpFY
Wikipedia contributors. “Kalman filter.” Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 1 Nov. 2018. Web. 5 Nov. 2018.
Chris Bracegirdle. Bayes’ Theorem for Gaussians. September 2010. http://web4.cs.ucl.ac.uk/staff/C.Bracegirdle/bayesTheoremForGaussians.pdf
Kevin P. Murphy. Multivariate Gaussians. Last updated September 28, 2007. https://www.cs.ubc.ca/~murphyk/Teaching/CS340-Fall07/reading/gauss.pdf
Wikipedia contributors. “Woodbury matrix identity.” Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 4 Nov. 2018. Web. 13 Nov. 2018.
Bilgin Esme. Kalman Filter For Dummies. http://bilgin.esme.org/BitsAndBytes/KalmanFilterforDummies
Statlect. Covariance matrix. https://www.statlect.com/fundamentals-of-probability/covariance-matrix
Greg Welch, Gary Bishop. The Kalman Filter. http://www.cs.unc.edu/~welch/kalman/
Simon D. Levy. The Extended Kalman Filter: An Interactive Tutorial for Non-Experts. https://home.wlu.edu/~levys/kalman_tutorial/