MATLABで学ぶ状態空間モデルの離散化｜差分近似とZOH（ゼロ次ホールド）の精度比較

はじめに：最先端の振動工学と「カルマンフィルタ」

近年、機械工学や振動工学のグローバルなトレンドとして、純粋な物理現象の深堀りに加え、デジタルツイン、状態監視（CBM）、1DCAE、MBD、メタマテリアルといったシステム横断的な技術が非常に注目されています。

これらの分野は、「振動工学 ＋ 制御工学 ＋ IoT・システム工学」といったように、これまでは別領域とみなされていた知識を掛け合わせて解を出す必要があります。これは、これからのエンジニアや研究者にとって、スキルを拡張する非常にエキサイティングなチャンスでもあります。

この最先端トレンドの核となる技術の一つが、制御工学でおなじみの「カルマンフィルタ（Kalman Filter）」を用いた振動推定技術です。
カルマンフィルタは、ノイズを含む観測データからシステムの状態を最適に推定する手法です。ここで「状態変数 ＝ 振動（変位や速度）」と定義することで、センサーを設置できない箇所の振動を計算上で高精度に推定（仮想センシング）できるようになり、デジタルツインや状態監視へダイレクトに応用できます。

制御工学の世界ではおなじみの理論ですが、近年では海外を中心に、これを振動工学が扱うような高周波帯域や多自由度の大規模モデルへ適用する研究が盛んに行われています。

本記事では、このカルマンフィルタを実装して振動推定を行うための「第一歩（前段階）」として、コンピュータで計算処理を行うために絶対に欠かせない「状態空間モデルの離散化」について分かりやすく解説します。

前回の「状態空間モデルの構築」に関する記事をまだお読みでない方は、ぜひ下記から先にご覧ください。

👉 MATLABで学ぶ状態空間モデル（N自由度モデル）

なお、この分野をより実践的に学ぶために、本記事では海外の優れたこちらの論文や、MATLABプログラムが豊富に掲載されている下記の専門書を参考に理論を構築しています。

状態空間モデルのおさらい（N自由度のバネ-マス-ダンパモデル）

前回の記事の要点だけを簡単におさらいします。図1のような $N$ 自由度の運動方程式は式(1)で表されます。

図1 バネマスモデル（N自由度）

$$[M][\ddot{y}(t)] + [C][\dot{y}(t)] + [K][y(t)] = [u(t)] \tag{1}$$

この式(1)を変形し、変位と速度を状態変数ベクトルとした「状態方程式」が式(2)、「出力方程式」が式(3)となります。

$$\begin{bmatrix} [\dot{x}_1] \\ [\dot{x}_2] \end{bmatrix} = \begin{bmatrix} [0] & [I] \\ -[M]^{-1}[K] & -[M]^{-1}[C] \end{bmatrix} \begin{bmatrix} [x_1] \\ [x_2] \end{bmatrix} + \begin{bmatrix} [0] \\ [M]^{-1} \end{bmatrix} [u(t)] \tag{2}$$

$$[y(t)] = \begin{bmatrix} [I] & [0] \end{bmatrix} \begin{bmatrix} [x_1] \\ [x_2] \end{bmatrix} \tag{3}$$

行列部分を一つの文字にまとめて一般化すると、連続時間系の状態空間モデルは式(4)のようにスッキリと表現できます。

$$\begin{cases} [\dot{x}(t)] = [A][x(t)] + [B][u(t)] \\ [y(t)] = [C][x(t)] \end{cases} \tag{4}$$

式(4)は時間が滑らかに進む「連続時間系（Continuous-time）」のモデルです。しかし、MATLABなどのコンピュータ上でシミュレーションやカルマンフィルタを実行するためには、時間を $k, k+1, \dots$ と一定間隔（サンプリング時間 $\Delta t$）で区切った「離散時間系（Discrete-time）」に変換する必要があります。

最終的なゴールは、以下の式(5)の形（離散化モデル）に変換し、係数行列 $[A_d]$ と $[B_d]$ を求めることです。

$$\begin{cases} [x(k+1)] = [A_d][x(k)] + [B_d][u(k)] \\ [y(k+1)] = [C_d][x(k+1)] \end{cases} \tag{5}$$

※出力行列 $[C]$ は時間を離散化しても変化しないため、そのまま $[C_d] = [C]$ として扱えます。

離散化手法1：差分方程式（オイラー法）による近似

もっとも直感的でわかりやすいのが、微分の定義を利用して差分近似する方法です。
式(4)の速度ベクトルの部分 $[\dot{x}(t)]$ は、時間変動率 $dx/dt$ のことですから、微小なサンプリング時間 $\Delta t$ を用いて式(6)のように近似できます（$k$ は現在のステップを表します）。

$$\frac{[x(k+1)] – [x(k)]}{\Delta t} = [A][x(k)] + [B][u(k)] \tag{6}$$

この式を「次のステップの状態 $[x(k+1)]$」について解き直すと、式(7)が得られます。

$$[x(k+1)] = ([I] + [A]\Delta t)[x(k)] + [B]\Delta t[u(k)] \tag{7}$$

式(7)と目標の式(5)を見比べると、$[A_d] = [I] + [A]\Delta t$、$[B_d] = [B]\Delta t$ とすれば良いことがわかります。

離散化手法2：ZOH（ゼロ次ホールド）による厳密な変換

先ほどの差分近似はわかりやすいですが、サンプリング時間 $\Delta t$ が大きいと誤差が蓄積するという弱点があります。
そこで、デジタル制御やシミュレーションで標準的に使われるのがZOH（Zero-Order Hold：ゼロ次ホールド）と呼ばれる手法です。

これは「サンプリング時間 $\Delta t$ の間、入力 $u(k)$ は一定値を保持する（ホールドされる）」という、実際のデジタル機器の動作に即した前提を置くことで、より厳密に解を求める手法です。
行列指数関数（マトリックス・エクススポネンシャル）を用いると、ZOHによる離散化行列 $[A_d], [B_d]$ は以下の式(8)の演算によって一括で求めることができます。（理論の詳細はWikipedia等をご参照ください）

$$\exp \left( \begin{bmatrix} [A] & [B] \\ [0] & [0] \end{bmatrix} \Delta t \right) = \begin{bmatrix} [A_d] & [B_d] \\ [0] & [I] \end{bmatrix} \tag{8}$$

MATLABでは expm 関数を使うことで、この行列指数関数を非常に簡単に計算できます。

MATLABで解く状態空間モデルの離散化検証

では、実際の例題を解いて、2つの離散化手法の精度を比較してみましょう。

【解析条件】

• 自由度： $N=10$
• 剛性： $k=10^5$
• 減衰： $c=1$
• 質量： $m=1$
• サンプリング周波数： $fs=1,048,576$（$dt = 1/fs$）※意図的に極めて細かく設定
• 初期状態： 質量1 ($m_1$) の変位のみ1、他は0（自由振動）
• 求めたい状態： 質量10 ($m_{10}$) の振動（変位）

まずはサンプリング周波数 $f_s$ を十分に大きく（時間刻みを細かく）設定して比較した結果が図2です。

図2 解法の比較（$f_s$ が十分大きい場合）

MATLABの高精度ソルバー「ode45」の計算結果を正解（真値）と仮定すると、”離散化1（差分近似）”と”離散化2（ZOH）”の両手法が見事に一致していることがわかります。

しかし、実用上は計算負荷を下げるため、もっとサンプリング周波数を小さく（時間刻みを粗く）設定したいケースが多くなります。そこで、サンプリング周波数 $f_s$ を徐々に小さくしていき、真値との誤差がどのように推移していくかを検証した結果が図3です。

図3 離散化手法の精度（誤差推移）検証

図3を見ると一目瞭然ですね！
差分近似（離散化1）はサンプリングを少しでも粗くすると急激に誤差が増大しますが、ZOH（離散化2）はサンプリング周波数を下げても圧倒的に高い精度を維持しています。
これが、カルマンフィルタや状態空間シミュレーションにおいて、ZOHによる離散化が強く推奨される理由です。

MATLABプログラム（実装コード）

ご自身の環境ですぐにシミュレーションを試せるよう、使用したMATLABコードを記載しておきます。ぜひコピペして動かし、未来の振動を予測する面白さを体感してみてください！

実行メインファイル

clear all; clc; close all

m_vec = ones(1,10);
k_vec = ones(1,10)*10^5;
[M] = eval_Mmatrix(m_vec); % 質量行列
[K] = eval_Kmatrix(k_vec); % 剛性行列
C = K*10^-5;

[V,D] = eig(K,M); % 固有値D、固有ベクトルV
wn = sqrt(diag(D)); % 固有角振動数
fn = wn/(2*pi); % 固有振動数

% % 状態方程式（連続時間系）
% % dx/dt = Ac*x(t) + Bc*p(t)
Ac = [zeros(size(M)) eye(size(M));
      -inv(M)*K     -inv(M)*C];
Bc = [zeros(size(M));
      -inv(M)];

Cc = [-inv(M)*K -inv(M)*C];
p = zeros(size(M,1),1);

x0 = zeros(size(M,1)*2,1);
x0(1) = 1; % 初期変位
t0 = 0;
dt = 1/400;
tf = 10;
time = t0:dt:tf;

% % % % MATLABのode45を使った場合 (理論値・真値と仮定) % % % %
[T,Y] = ode45(@state_equation, time, x0); 
figure
plot(T, Y(:,10), 'k', 'linewidth', 2)
hold on
xlabel('Time [s]')
ylabel('Displacement [m]')

% % % % 状態方程式（離散時間系） % % % %

% --- 離散化1：差分方程式からの離散化 ---
fs = 2^20; % 高いサンプリング周波数で検証
t0_d = 0;
dt_d = 1/fs;
tf_d = 4;
time_d = t0_d:dt_d:tf_d;

Adis = eye(size(Ac)) + Ac*dt_d; % 離散化したA
Bdis = Bc*dt_d;                 % 離散化したB

xdis = zeros(size(Adis,1), length(time_d));
xdis(:,1) = x0;
for ii1 = 1:length(time_d)-1
    xdis(:,ii1+1) = Adis*xdis(:,ii1) + Bdis*p;
end
plot(time_d, xdis(10,:), 'b-', 'linewidth', 1)

% --- 離散化2：ZOH（Zero-Order Hold） ---
temp = [Ac, Bc;
        zeros(size(Bc,2), size(Ac,2)), zeros(size(Bc,2), size(Bc,2))];
temp2 = expm(temp*dt_d); % 行列指数関数
Adis2 = temp2(1:size(Ac,1), 1:size(Ac,2));
Bdis2 = temp2(1:size(Ac,1), size(Ac,2)+1:end);

xdis2 = zeros(size(Adis2,1), length(time_d));
xdis2(:,1) = x0;
for ii1 = 1:length(time_d)-1
    xdis2(:,ii1+1) = Adis2*xdis2(:,ii1) + Bdis2*p;
end
plot(time_d, xdis2(10,:), 'r--', 'linewidth', 1)

xlim([time_d(1) time_d(end)])
legend('時刻歴積分 ode45 (真値と仮定)', '離散化1：差分近似', '離散化2：ZOH')

外部関数：eval_Mmatrix.m

function [M]=eval_Mmatrix(m_vec)
    M=diag(m_vec);
end

外部関数：eval_Kmatrix.m

function [K]=eval_Kmatrix(k_vec)
    K=zeros(length(k_vec));
    for ii1=1:length(k_vec)
        if ii1==1
            K(ii1,ii1)=k_vec(ii1);
        else
            K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);
        end
    end
end

外部関数：state_equation.m

function dx=state_equation(t,x)
    m_vec=ones(1,10);
    k_vec=ones(1,10)*10^5;
    [M]=eval_Mmatrix(m_vec); %質量行列
    [K]=eval_Kmatrix(k_vec); %剛性行列
    C=K*10^-5;

    % 状態方程式
    % dx/dt = Ac*x(t) + Bc*p(t)
    Ac=[zeros(size(M)) eye(size(M));
        -inv(M)*K -inv(M)*C];
    Bc=[zeros(size(M));
        -inv(M)];

    p=zeros(size(M,1),1);
    dx=Ac*x+Bc*p;
end