「続々制御工学のこころ」を読む　～第8章システム同定 ④補足 : 部分空間同定法～

はじめに

データから対象システムの入出力関係をモデリングするための手法であるシステム同定について理解するために，足立修一著「続々制御工学のこころ」を手に取ることにした。数理統計学の知識に基づき，「数理統計学→システム同定」という順にシステム同定を理解しようとしている。

本記事は，「第8章システム同定」の内容から少し発展して，部分空間同定法における状態変数の推定方法について調べた内容である。

本書の紹介ページ

続々制御工学のこころ - 東京電機大学出版局科学技術と教育を出版からサポートする

本記事を読むことで得られること

本記事を読むことで得られることは，主に以下の内容である。

部分空間同定法における状態変数の推定方法

改訂履歴

「状態変数の推定」の項における「すなわち，解くべき最適化問題は以下のようになる。」以降の $\mathbfit{X}_f$ を $\mathbfit{Y}_f$ に修正した。tty273533165様，ありがとうございました(2025/10/31)。

モチベーション

本書の「8.6.2 部分空間同定法」では，状態変数 $\mathbfit{x}$ ，入力 $\mathbfit{u}$ ，出力 $\mathbfit{y}$ が求められれば，最小二乗法によって状態空間モデルの係数パラメータが求められるという説明があった。そこで状態変数 $\mathbfit{x}$ の推定方法について調べた。

今回参考にしたのは井手剛著「入門機械学習による異常検知 - Rによる実践ガイド -」の第7.4節である。
www.coronasha.co.jp

対象とするシステム

対象とするシステムは，線形状態空間モデルである。ただし，入力 $\mathbfit{u}$ は無いものとする。

$\begin{align} &\mathbfit{x}^{(t)} \approx \mathbfit{A} \mathbfit{x}^{(t-1)} \tag{1}\\ \\ &\mathbfit{y}^{(t)} \approx \mathbfit{C} \mathbfit{x}^{(t)} \tag{2} \\ \\ &\mathbfit{x}^{(t)} \in \mathbb{R}^m, \mathbfit{y}^{(t)} \in \mathbb{R}^M, \mathbfit{A} \in \mathbb{R}^{m \times m}, \mathbfit{C} \in \mathbb{R}^{M \times m} \\ \end{align}$

ここで「 $\approx$ 」は，「確率的なノイズは別にして等しい」ということを示す。

(1)式は状態方程式，(2)式は観測方程式を表す。多くの場合，ノイズは正規分布に従うとするので，状態方程式と観測方程式は，以下のような確率モデルになる。

$\begin{align} p(\mathbfit{x}^{(t)} \lvert \mathbfit{x}^{(t-1)}) &= \mathcal{N}(\mathbfit{x}^{(t)} \lvert \mathbfit{A} \mathbfit{x}^{(t-1)}, \mathbfit{Q}), \quad \mathbfit{Q} \in \mathbb{R}^{m \times m} \\ \\ p(\mathbfit{y}^{(t)} \lvert \mathbfit{x}^{(t)}) &= \mathcal{N}(\mathbfit{x}^{(t)} \lvert \mathbfit{C} \mathbfit{x}^{(t)}, \mathbfit{R}), \quad \mathbfit{R} \in \mathbb{R}^{M \times M} \\ \end{align}$

さらに，始点の $t=1$ に対応するために， $p(\mathbfit{x}^{(1)} ) = \mathcal{N}(\mathbfit{x}^{(1)} \lvert \mathbfit{x}_0, \mathbfit{Q}_0)$ とする。

問題設定

今回の問題設定における入力は，十分に長い観測データの系列 $\mathcal{D} = \{ \mathbfit{y}^{(1)}, \mathbfit{y}^{(2)}, \cdots \}$ である。
求めたいものは，状態変数の系列 $\{ \mathbfit{x}^{(1)}, \cdots, \mathbfit{x}^{(t)} \}$ である。

なお，状態変数の系列が求められれば，状態空間モデルの係数パラメータ $\mathbfit{A}, \mathbfit{C}, \mathbfit{Q}, \mathbfit{R}$ も求められる。

特異値分解・正準相関分析を用いた状態変数の推定

「続々制御工学のこころ」を読む～第6章ビヨンド最小二乗法 ⑥補足：正準相関分析(実践編)～ - jiku logにおいて説明したように，正準相関分析では2つの時系列の相関が大きくなるように正準変数を求めると，その正準変数は2つの時系列に共通するような潜在変数となっている。今回の問題設定における状態変数 $\mathbfit{x}$ も，直接観測できない潜在的な変数である。そのため，正準相関分析と同様の考え方で，状態変数を推定する。

特徴抽出に向けた時系列データの変形

観測データから特徴的なパターンを抽出するために，「続々制御工学のこころ」を読む～第6章ビヨンド最小二乗法 ③補足：時系列データに対する特異値分解(理論編)～ - jiku logにならって，窓幅 $w$ の部分時系列を $N$ 個作成する。たとえば，時刻 $t$ を先頭とする部分時系列ベクトルは次のようになる。

$\begin{align} \mathbfit{Y}^{(t)} \equiv \begin{bmatrix} \mathbfit{y}^{(t)} \\ \mathbfit{y}^{(t+1)} \\ \vdots \\ \mathbfit{y}^{(t+w-1)} \\ \end{bmatrix} \in \mathbb{R}^{Mw} \\ \end{align}$

ただし $N \gg Mw$ である。

過去の領域と未来の領域のデータ作成

まず観測データから2種類の時系列データを作成する。時刻 $t$ を境に，過去の領域(領域"p"(past))と未来の領域(領域"f"(future))に分け，長さ $w$ の部分時系列のベクトルを列ベクトルとする行列 $\mathbfit{Y}_p$ と $\mathbfit{Y}_f$ を次のように構成する。

$\begin{align} \mathbfit{Y}_p &\equiv \begin{bmatrix} \mathbfit{Y}^{(t-N)} & \cdots & \mathbfit{Y}^{(t-2)} & \mathbfit{Y}^{(t-1)} \\ \end{bmatrix} \in \mathbb{R}^{Mw \times N} \\ \\ \mathbfit{Y}_f &\equiv \begin{bmatrix} \mathbfit{Y}^{(t)} & \mathbfit{Y}^{(t+1)} & \cdots & \mathbfit{Y}^{(t+N-1)} \\ \end{bmatrix} \in \mathbb{R}^{Mw \times N} \\ \\ \end{align}$

観測データと状態変数の関係式

過去・未来のそれぞれにおける観測データと状態変数の関係式

行列 $\mathbfit{Y}_p$ ， $\mathbfit{Y}_f$ と状態変数の関係を調べる。(2)式を用いると，

$\begin{align} \mathbfit{y}^{(1)} \approx \mathbfit{C} \mathbfit{x}^{(1)}, \quad \mathbfit{y}^{(2)} \approx \mathbfit{C} \mathbfit{A} \mathbfit{x}^{(1)}, \quad \mathbfit{y}^{(3)} \approx \mathbfit{C} \mathbfit{A}^2 \mathbfit{x}^{(1)}, ... \\ \end{align}$

ということが分かる。この関係式を用いると，

$\begin{align} \mathbfit{Y}^{(t)} &= \begin{bmatrix} \mathbfit{y}^{(t)} \\ \mathbfit{y}^{(t+1)} \\ \vdots \\ \mathbfit{x}^{(t+w-1)} \\ \end{bmatrix} \approx \begin{bmatrix} \mathbfit{C} \mathbfit{x}^{(t)} \\ \mathbfit{C} \mathbfit{A} \mathbfit{x}^{(t)} \\ \vdots \\ \mathbfit{C} \mathbfit{A}^{w-1} \mathbfit{x}^{(t)} \\ \end{bmatrix} = \begin{bmatrix} \mathbfit{C} \\ \mathbfit{C} \mathbfit{A} \\ \vdots \\ \mathbfit{C} \mathbfit{A}^{w-1} \\ \end{bmatrix} \mathbfit{x}^{(t)} \equiv \mathbfit{\Gamma} \mathbfit{x}^{(t)} \\ \end{align}$

のように， $\mathbfit{Y}^{(t)}$ を $\mathbfit{x}^{(t)}$ を使って表現できる。ただし $\mathbfit{\Gamma} \in \mathbb{R}^{Mw \times m}$ は拡大可観測行列である。
この関係を用いると，

$\begin{align} \mathbfit{Y}_p &= \begin{bmatrix} \mathbfit{Y}^{(t-N)} & \cdots & \mathbfit{Y}^{(t-2)} & \mathbfit{Y}^{(t-1)} \\ \end{bmatrix} \\ \\ &\approx \begin{bmatrix} \mathbfit{\Gamma} \mathbfit{x}^{(t-N)} & \cdots & \mathbfit{\Gamma} \mathbfit{x}^{(t-2)} & \mathbfit{\Gamma} \mathbfit{x}^{(t-1)} \\ \end{bmatrix} \equiv \mathbfit{\Gamma} \mathbfit{X}_p \\ \\ \mathbfit{Y}_f &= \begin{bmatrix} \mathbfit{Y}^{(t)} & \mathbfit{Y}^{(t+1)} & \cdots & \mathbfit{Y}^{(t+N-1)} \\ \end{bmatrix} \\ \\ &\approx \begin{bmatrix} \mathbfit{\Gamma} \mathbfit{x}^{(t)} & \mathbfit{\Gamma} \mathbfit{x}^{(t+1)} & \cdots & \mathbfit{\Gamma} \mathbfit{x}^{(t+N-1)} \\ \end{bmatrix} \equiv \mathbfit{\Gamma} \mathbfit{X}_f \\ \\ \end{align}$

ただし，

$\begin{align} \mathbfit{X}_p &\equiv \begin{bmatrix} \mathbfit{x}^{(t-N)} & \cdots & \mathbfit{x}^{(t-2)} & \mathbfit{x}^{(t-1)} \\ \end{bmatrix} \in \mathbb{R}^{m \times N} \\ \\ \mathbfit{X}_f &\equiv \begin{bmatrix} \mathbfit{x}^{(t)} & \mathbfit{x}^{(t+1)} & \cdots & \mathbfit{x}^{(t+N-1)} \\ \end{bmatrix} \in \mathbb{R}^{m \times N} \\ \\ \end{align}$

である。これらの関係を転置して書き直すと，

$\begin{align} \mathbfit{Y}_p^T \approx \mathbfit{X}_p^T \mathbfit{\Gamma}^T , \quad \mathbfit{Y}_f^T \approx \mathbfit{X}_f^T \mathbfit{\Gamma}^T \\ \end{align}$

となる。

過去と未来の状態変数の関係

過去の状態変数の行列 $\mathbfit{X}_p$ と未来の状態変数の行列 $\mathbfit{X}_f$ の関係について， $\mathbfit{X}_p$ の時刻を $N$ ステップ進めたものが $\mathbfit{X}_f$ なので， $\mathbfit{X}_f \approx \mathbfit{A}^N \mathbfit{X}_p$ ，すなわち

$\mathbfit{X}_f^T \approx \mathbfit{X}_p^T (\mathbfit{A}^N)^T \quad(\mathbfit{X}_f^T \in \mathbb{R}^{N \times m}, \mathbfit{A}^N \in \mathbb{R}^{m \times m}) \\$

となっている。

ここで，

$\mathbfit{X}_f^T$ の列ベクトルが張る列空間 $\mathrm{Span}\{ \mathbfit{X}_f^T \}$
$\mathbfit{X}_p^T$ の列ベクトルが張る列空間 $\mathrm{Span}\{ \mathbfit{X}_p^T \}$

の関係を考える。

$\mathbfit{X}_f^T$ は， $\mathbfit{X}_p^T$ の右から $(\mathbfit{A}^N)^T$ をかけたものであるが，これは $\mathbfit{X}_p^T$ の各行に $(\mathbfit{A}^N)^T$ をかけること，すなわち時間を $N$ ステップ進めることに相当している。 $(\mathbfit{A}^N)^T$ が可逆だとすると列空間が変化しないので， $\mathbfit{X}_p^T$ の列空間は，時間が経過しても平均的には変化しないと言える。

したがって，過去における状態変数と未来における状態変数ができる限り一致するように状態空間の行列を選ぶことで，状態変数についての最善の推定ができると言える。

状態変数の推定

煩雑さを避けるために，状態変数の次元数が $m=1$ であるとする。過去の観測データの行列と未来の観測データの行列を用いて，状態変数を以下のようにモデル化する。

$\begin{align} \mathbfit{X}_p^T = \mathbfit{Y}_p^T \mathbfit{\alpha} , \quad \mathbfit{X}_f^T = \mathbfit{Y}_f^T \mathbfit{\beta} \\ \end{align}$

ただし， $\mathbfit{\alpha}, \mathbfit{\beta} \in \mathbb{R}^{Mw}$ は，これから求める未知の係数ベクトルである。このベクトルは，過去における状態変数と未来における状態変数ができる限り一致するように選ぶので，正準相関分析と同様に，これらの相関係数を最大化するようにする。

すなわち，解くべき最適化問題は以下のようになる。

$\begin{align} &\max_{\mathbfit{\alpha}, \mathbfit{\beta}} \mathbfit{\alpha}^T \mathbfit{Y}_p \mathbfit{Y}_f^T \mathbfit{\beta} \\ \\ &\text{subject to } \mathbfit{\alpha}^T \mathbfit{Y}_p \mathbfit{Y}_p^T \mathbfit{\alpha} = \mathbfit{\beta}^T \mathbfit{Y}_f \mathbfit{Y}_f^T \mathbfit{\beta} = 1 \\ \end{align}$

「続々制御工学のこころ」を読む～第6章ビヨンド最小二乗法 ⑤補足：正準相関分析(理論編)～ - jiku log を参考に式変形を進めると，以下のような一般化固有値方程式を解くことになる。

$\begin{align} &\mathbfit{S}_{pf} \mathbfit{S}_{ff}^{-1} \mathbfit{S}_{fp} \mathbfit{\alpha} = \lambda^2 \mathbfit{S}_{pp} \mathbfit{\alpha} \\ \\ &\mathbfit{S}_{fp} \mathbfit{S}_{pp}^{-1} \mathbfit{S}_{pf} \mathbfit{\beta} = \lambda^2 \mathbfit{S}_{ff} \mathbfit{\beta} \\ \\ &\mathbfit{S}_{pf} \equiv \mathbfit{Y}_p \mathbfit{Y}_f^T, \mathbfit{S}_{pp} \equiv \mathbfit{Y}_p \mathbfit{Y}_p^T, \mathbfit{S}_{fp} \equiv \mathbfit{Y}_f \mathbfit{Y}_p^T, \mathbfit{S}_{ff} \equiv \mathbfit{Y}_f \mathbfit{Y}_f^T \\ \end{align}$

これらの一般化固有方程式は，行列

$\begin{align} \mathbfit{W} = \mathbfit{S}_{pp}^{-1/2} \mathbfit{S}_{pf} \mathbfit{S}_{ff}^{-1/2} \\ \end{align}$

の特異値分解に帰着される。

ここまでは $m=1$ としていたが，特異値分解により複数の特異値・特異ベクトルを求めることで， $m \gt 1$ の解も求められる。

$\mathbfit{W}$ の左特異ベクトルの特異値の大きい順に $\tilde{\mathbfit{\alpha}}^1,...,\tilde{\mathbfit{\alpha}}^m$ ，右特異ベクトルを $\tilde{\mathbfit{\beta}}^1,...,\tilde{\mathbfit{\beta}}^m$ とすると，

$\begin{align} \mathbfit{\alpha}^i = \mathbfit{S}_{pp}^{-1/2} \tilde{\mathbfit{\alpha}}^i, \quad \mathbfit{\beta}^i = \mathbfit{S}_{ff}^{-1/2} \tilde{\mathbfit{\beta}}^i \\ \end{align}$

のように，上記の最適化問題の解が求められる。

最終的に，求めたい状態変数の系列の推定値は，以下のように求められる。

$\begin{align} \hat{\mathbfit{X}}_p = \begin{bmatrix} \mathbfit{\alpha}^1 & \cdots & \mathbfit{\alpha}^m \\ \end{bmatrix} ^T \mathbfit{Y}_p ,\quad \hat{\mathbfit{X}}_f = \begin{bmatrix} \mathbfit{\beta}^1 & \cdots & \mathbfit{\beta}^m \\ \end{bmatrix} ^T \mathbfit{Y}_f \\ \end{align}$

まとめと感想

今回は，「第8章システム同定」の内容から少し発展して，部分空間同定法における状態変数の推定方法についてまとめた。

状態変数を推定するために，「状態変数は，過去の観測データ系列と未来の観測データ系列の共通的な潜在変数となる」という考え方はとても興味深かった。今回は観測値のデータ行列の積をもとに議論を進めていたが，正準相関分析では共分散行列や分散共分散行列を用いていた。ただ今回の議論において，期待値を取れば正準相関分析にかなり近くなる。

今回調べた方法では，基本的に特異値分解などの線形代数の計算のみで状態変数が求められるうえに，考え方も結構シンプルだったで便利な手法だと感じた。数値実験なども試してみたい。

本記事を最後まで読んでくださり，どうもありがとうございました。