Eigensystem Realization Algorithm (ERA) によるシステム同定の理論と実装

概要

構造ヘルスモニタリングにおけるシステム同定の手法として用いられる、Eigensystem Realization Algorithm (ERA) の理論と実装例について整理する。

Eigensystem Realization Algorithm (ERA) とは

加速度センサ等で計測した構造物の振動データから、構造物のパラメータを同定するシステム同定の手法である。得られた固有振動数、モード減衰比、モード形は、構造ヘルスモニタリングなどにおいて構造特性の変化を分析する上で利用される。

ERAでは自由振動波形を利用する。自由振動を直接計測することが難しい場合も多いため、ランダム応答から自由振動波形を抽出する手法である Random Decrement (RD法)Natural Excitation Technique (NExT) と併せて利用されることが多い。

ERAの理論

離散時間線形システムとして、次のような状態空間モデルを対象とする。

$$ \begin{align} & x(k+1) = Ax(k)+Bu(k) \tag{1} \\ & y(k) = Cx(k) + Du(k) \tag{2} \end{align} $$

このシステムに対する各入力点のインパルス応答を考慮すると、マルコフパラメータとして次が得られる。

$$ Y(k) = \{D, CB, CAB, CA^{2}B, \ldots, CA^{k - 1}B\} \tag{3} $$

$Y(k)$ について、ハンケル行列は次のように表される。$m,n$ は任意であり、ハンケル行列は正方行列である必要はない。

$$ H(k) = \begin{bmatrix} Y(k) & Y(k+1) & \ldots & Y(k+n) \\ Y(k+1) & Y(k+2) & \ldots & Y(k+n+1) \\ \vdots & \vdots & \ddots & \vdots \\ Y(k+m) & Y(k+m+1) & \ldots & Y(k+m+n) \end{bmatrix} \tag{4} $$

$H(1)$ に対する特異値分解により、直行行列 $U,V$ および特異値を対角成分にもつ対角行列 $Σ$ が得られる。*1 $$ H(1) = U Σ V^{T} \tag{5} $$

$Σ$ の対角成分は左上から特異値の絶対値の大きいものが並ぶ。特異値の大きなモードを抽出するように任意のモード数 $r$ を定めることで、有意な成分(添字$r$)とノイズ成分(添字$n$)によるモードを分離することができ、次のように書き直すことができる。

$$ H(1) = \begin{bmatrix} U_r & U_n \end{bmatrix} \begin{bmatrix} Σ_r & 0 \\ 0 & Σ_n \end{bmatrix} \begin{bmatrix} V_r^{T} \\ V_n^{T} \end{bmatrix} \tag{6} $$

さらに、ノイズ成分を無視することで最小実現(モード数を低減した形)のひとつとして次のように得られる。

$$ H(1) = U_r Σ_r V_r^{T} \tag{7} $$

また、$H(1),H(2)$ は自由振動時において可観測行列 $W_O$ と可制御行列 $W_C$ を用いて次のように表される。

$$ H(1) = \begin{bmatrix} C \\ CA \\ CA^{2} \\ \vdots \\ CA^{m} \end{bmatrix} \begin{bmatrix} B & AB & A^{2}B & \ldots & A^{n}B \end{bmatrix} = W_O W_C \tag{8} $$

$$ H(2) = \begin{bmatrix} C \\ CA \\ CA^{2} \\ \vdots \\ CA^{m} \end{bmatrix} A \begin{bmatrix} B & AB & A^{2}B & \ldots & A^{n}B \end{bmatrix} = W_O A W_C \tag{9} $$

式(7)と式(8)から、次が得られる。

$$ \begin{align} & W_O = U_r Σ_r^{\frac{1}{2}} \tag{10} \\ & W_C = Σ_r^{\frac{1}{2}} V_r^{T} \tag{11} \end{align} $$

式(10)および式(11)を式(9)に代入し、$A$ について整理することで、最小実現の状態行列 $A_r$ が得られる。

$$ A_r = Σ_r^{-\frac{1}{2}} U_r^{T} H(2) V_r Σ_r^{-\frac{1}{2}} \tag{12} $$

$B$ は可制御行列 $W_C$ の先頭の要素であるため、次により最小実現の入力行列 $B_r$ が得られる。 $I_{out}$ は行数および列数が出力点数 $n_{out}$ の単位行列である。

$$ B_r = Σ_r^{\frac{1}{2}} V_r^{T} \begin{bmatrix} I_{out} \\ 0 \\ \vdots \\ 0 \end{bmatrix} \tag{13} $$

$C$ は可観測行列 $W_O$ の先頭の要素であるため、次により最小実現の出力行列 $C_r$ が得られる。 $I_{in}$ は行数および列数が入力点数 $n_{in}$ の単位行列である。

$$ C_r = \begin{bmatrix} I_{in} & 0 & \ldots & 0 \end{bmatrix} U_r Σ_r^{\frac{1}{2}} \tag{14} $$

最小実現の直達行列 $D_r$ は、式(3)から次のように得られる。

$$ D_r = Y(0) \tag{15} $$

以上により、観測値として得られる $Y(k)$ を用いて、線形システムの各パラメータ $A_r,B_r,C_r,D_r$ を同定することができる。

ERAの実装

Pythonにて簡単な例を実装する。書籍 Data Driven Science & Engineering: Machine Learning, Dynamical Systems, and Control のデモとしてGitHubに公開されているコードを参考に実装した。

また、上記のリポジトリで公開されている状態空間モデルのサンプルを使用し、インパルス応答の算出結果を擬似的な観測値としている。実行するには下記より testSys_ABCD.mat ファイルをダウンロードし、カレントディレクトリに保存する必要がある。

環境

下記の外部ライブラリを使用する。

  • numpy
  • matplotlib
  • scipy
  • control

ソースコード

2入力、2出力、100自由度の線形システムを対象にインパルス応答を算出し、ERAによる低減後のモード数を10としてパラメータの同定を行った。

import numpy as np
import matplotlib.pyplot as plt
from scipy import io
from scipy.linalg import fractional_matrix_power
from control import matlab


def hankel(y):
    ntime, noutputs, ninputs = y.shape
    n = (ntime + 1) // 2
    H = np.zeros((noutputs * n, ninputs * n))
    for ioutput in range(noutputs):
        for iinput in range(ninputs):
            for i in range(n):
                for j in range(n):
                    H_row = noutputs * i + ioutput
                    H_col = ninputs * j + iinput
                    H[H_row, H_col] = y[i + j, ioutput, iinput]
    return H


def ERA(y, r):
    _, noutputs, ninputs = y.shape
    H1 = hankel(y[1:-1])
    H2 = hankel(y[2:])

    U, S, VT = np.linalg.svd(H1, full_matrices=False)
    V = VT.T
    Sigma = np.diag(S[:r])
    Ur = U[:, :r]
    Vr = V[:, :r]

    Ar = fractional_matrix_power(Sigma, -0.5) @ Ur.T @ H2 @ Vr @ fractional_matrix_power(Sigma, -0.5)
    Br = fractional_matrix_power(Sigma, -0.5) @ Ur.T @ H1[:, :ninputs]
    Cr = H1[:noutputs, :] @ Vr @ fractional_matrix_power(Sigma, -0.5)
    Dr = y[0, :, :]
    return Ar, Br, Cr, Dr


def main():
    # 低減後のモード数
    reduced_model_order = 10

    # 波形の長さ
    ntime = 52

    # 時間間隔
    dt = 1.0

    # システムモデルの読み込み
    testSys_mat = io.loadmat("testSys_ABCD.mat")
    A = testSys_mat["A"]
    B = testSys_mat["B"]
    C = testSys_mat["C"]
    D = testSys_mat["D"]
    sys = matlab.ss(A, B, C, D, dt=dt)
    ninputs = sys.ninputs
    noutputs = sys.noutputs

    # 各点へのインパルス入力に対する応答の算出
    t = np.arange(ntime) * dt
    y = np.zeros((ntime, noutputs, ninputs))
    for iinput in range(ninputs):
        y[:, :, iinput : iinput + 1], _ = matlab.impulse(sys, T=t, input=iinput)

    # システムの同定
    Ar, Br, Cr, Dr = ERA(y, reduced_model_order)
    sysr = matlab.ss(Ar, Br, Cr, Dr, dt=dt)

    # 同定したシステムに対するインパルス応答の算出
    yr = np.zeros((ntime, noutputs, ninputs))
    for iinput in range(ninputs):
        yr[:, :, iinput : iinput + 1], _ = matlab.impulse(sysr, T=t, input=iinput)

    # 結果のプロット
    fig, axes = plt.subplots(2, 2, figsize=(8, 4))
    for ioutput in range(noutputs):
        for iinput in range(ninputs):
            axes[ioutput, iinput].plot(t, y[:, ioutput, iinput], linewidth=3, color="dimgray", label="Observation")
            axes[ioutput, iinput].plot(t, yr[:, ioutput, iinput], linewidth=1.5, color="tomato", label="ERA")
            axes[ioutput, iinput].set_xlim(0, None)
            axes[ioutput, iinput].set_ylim(-10, 10)
    axes[0, 0].set_title("Input point 0")
    axes[0, 1].set_title("Input point 1")
    axes[0, 0].set_ylabel("Output point 0")
    axes[1, 0].set_ylabel("Output point 1")
    axes[0, 0].legend()
    fig.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

出力結果

各点に対するインパルス入力とその出力をグラフにプロットした。元のシステムの自由度が100でERAでは10モードに低減されているが、元のシステム(Observation)とERAで同等のインパルス応答が得られており、支配的なモードの抽出ができている。

補足

多入力多出力システムにおけるハンケル行列の構成

多入力多出力システムの場合、ハンケル行列の各要素 $Y(k)$ は入出力に対する行列として記述される。出力点数$n_{out}$,入力点数 $n_{in}$ の場合、式(4)は次のような行数 $(m+1)n_{out}$,列数 $(n+1)n_{in}$ の行列となる。

$$ H(k) = \begin{bmatrix} y(k,0,0) & \ldots & y(k,0,n_{in} - 1) & y(k+1,0,0) & \ldots \\ \vdots & \ddots & \vdots & \vdots & \ddots \\ y(k,n_{out} - 1,0) & \ldots & y(k,n_{out} - 1,n_{in} - 1) & y(k+1,n_{out} - 1,0) & \ldots \\ y(k+1,0,0) & \ldots & y(k+1,0,n_{in} - 1) & y(k+2,0,0) & \ldots \\ \vdots & \ddots & \vdots & \vdots & \ddots \\ \end{bmatrix} $$

ハンケル行列の生成方法

ハンケル行列の生成については、scipy ライブラリにて scipy.linalg.hankel() として提供されているものもあるが、多入力多出力系に対応できないため、前述の実装例では hankel() 関数を下記のように実装している。

def hankel(y):
    ntime, noutputs, ninputs = y.shape
    n = (ntime + 1) // 2
    H = np.zeros((noutputs * n, ninputs * n))
    for ioutput in range(noutputs):
        for iinput in range(ninputs):
            for i in range(n):
                for j in range(n):
                    H_row = noutputs * i + ioutput
                    H_col = ninputs * j + iinput
                    H[H_row, H_col] = y[i + j, ioutput, iinput]
    return H

hankel() 関数は下記のようにすることで、可読性は低下するが計算効率は向上する。

def hankel(y):
    ntime, noutputs, ninputs = y.shape
    n = (ntime + 1) // 2
    ind = np.indices((n, n)).sum(axis=0)
    H = y[ind].transpose(0, 2, 1, 3).reshape(noutputs * n, ninputs * n)
    return H

参考文献

*1:インパルス応答の場合、$Y(0)$ はインパルス入力の直達成分である $D$ となるため、自由振動となる $Y(1)$ 以降を対象とする。そのため、特異値分解を行うハンケル行列は $H(0)$ ではなく $H(1)$ とする。