【MATLAB】状態空間モデルの解法を比較｜ode45・オイラー法・ルンゲクッタ法による振動解析

はじめに：振動工学と制御工学の融合が求められる時代

近年、海外の大学や最先端の産業界では、純粋な振動工学や機械工学の枠を超え、デジタルツイン、状態監視（CBM）、1DCAE、MBD、メタマテリアルといったシステム横断的な技術が大きなトレンドになっています。

これらの分野は、「振動工学 ＋ 制御工学 ＋ IoT ＋ システム工学」といったように、ひと昔前であれば「他分野」と見なされていた知識を融合（横通し）させる必要があります。裏を返せば、これからのエンジニアにとって、自分の専門分野に制御やデータ処理の知識を掛け合わせることで、市場価値を劇的に高められる大チャンスでもあります。

このトレンドの核となる技術の一つが、制御工学でおなじみの「カルマンフィルタ」です。カルマンフィルタを用いて「状態変数 ＝ 振動（変位や速度）」を推定することで、直接センサーを置けない場所の振動を予測でき、デジタルツインや状態監視にダイレクトに応用できます。

海外では、制御工学の枠にとどまらず、振動工学が扱うような高周波帯域や多自由度モデルへカルマンフィルタを適用する研究が盛んに進められています。本記事の理論構築にあたっては、海外のこちらの論文や、MATLABプログラムが豊富に掲載されている下記の専門書を参考にしています。

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

 

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

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

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

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

これを変形した「状態方程式」と「出力方程式」は以下のようになります。

$$ \begin{bmatrix} [\dot{x}_1(t)] \\ [\dot{x}_2(t)] \end{bmatrix} = \begin{bmatrix} [0] & [I] \\ -[M]^{-1}[K] & -[M]^{-1}[C] \end{bmatrix} \begin{bmatrix} [x_1(t)] \\ [x_2(t)] \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(t)] \\ [x_2(t)] \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)は「微分方程式」です。MATLABにはこれを高精度に解いてくれる非常に便利な関数 ode45 が用意されています。
しかし、将来的にカルマンフィルタをマイコンやC言語で組み込み実装する場合、便利なブラックボックス関数は使えません。そこで本記事では、微分の基礎である「オイラー法」と、実務で多用される「4次のルンゲクッタ法」を自作して ode45 と比較してみます。

 

オイラー法（Euler Method）

オイラー法は、現在の状態と傾きから、$\Delta t$ 後の状態 $[x(t+\Delta t)]$ を直線的に近似計算する最もシンプルな方法です。

$$ [x(t+\Delta t)] = [x(t)] + \left( [A][x(t)] + [B][u(t)] \right) \Delta t \tag{5} $$

計算は軽いですが、振動系のようなカーブを描く現象に対しては誤差が蓄積しやすいのが弱点です。

 

4次のルンゲクッタ法（Runge-Kutta Method）

ルンゲクッタ法は、区間 $\Delta t$ の間で傾きを4回計算（評価）し、それらを加重平均することで極めて高い精度で近似する手法です。シミュレーションの実務では定番中の定番です。

$$ [x(t+\Delta t)] = [x(t)] + \frac{d_1 + 2d_2 + 2d_3 + d_4}{6} \tag{6} $$

ここで、状態方程式を $\dot{x}(t) = f(x(t))$ と置いたとき、4つの傾き成分 $d_1 \sim d_4$ は次のように計算されます。

$$ \begin{cases} d_1 = f(x(t)) \Delta t \\ d_2 = f\left(x(t) + \frac{d_1}{2}\right) \Delta t \\ d_3 = f\left(x(t) + \frac{d_2}{2}\right) \Delta t \\ d_4 = f(x(t) + d_3) \Delta t \end{cases} \tag{7} $$

一見複雑に見えますが、プログラムに落とし込むと非常にシンプルで強力です。

 

MATLABで解く状態空間モデル（各手法の精度比較）

では、前回の例題を使って「ode45」「オイラー法」「ルンゲクッタ法」を比較してみましょう。

【解析条件】

• 自由度： $N=10$
• 剛性： $k=10^5$
• 減衰： $c=1$
• 質量： $m=1$
• サンプリング周波数： $f_s=400$（$\Delta t = 1/f_s$）
• 初期状態： 質量1 ($m_1$) の変位のみ1、他は0
• 求めたい状態変数： 各質点の変位

全自由度の変位と速度が求まりますが、ここでは一番端にある 質量10 ($m_{10}$) の変位 の結果だけをプロットして紹介します。

図2 解法の比較

MATLABの公式関数である「ode45」を正解（真値）とみなしてグラフを見ると、結果は一目瞭然です。

• 4次のルンゲクッタ法（赤線）： ode45の黒線と完全に重なっており、極めて高精度に解析できています。
• オイラー法（緑破線）： 振幅がどんどん発散（増大）してしまっており、今回のような振動系のシミュレーションには全く使い物にならないことがわかります。

このように、連続時間系のモデルを数値計算で解く場合、最低でも4次のルンゲクッタ法を用いるか、素直に ode45 などのソルバーを使う必要があるということが確認できました。

本記事では、状態空間モデルを連続時間系のまま解く手法について解説しました。この知識を持っておけば、将来C言語などで独自のシミュレータを作る際にも困りません。
使用した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];

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


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


% % % % % % % 2. オイラー法による数値積分 % % % % % % %
time_vec = T.'; 
dt_step = T(2)-T(1);
p = zeros(size(M,1), 1); % 入力(外力)はゼロ
x_Euler_Method = zeros(length(M)*2, length(time_vec));
x = x0;

for ii1 = 1:length(time_vec)
    x_Euler_Method(:,ii1) = x;
    dx = Ac*x + Bc*p;
    x = x + dx*dt_step; % 直線近似
end
plot(time_vec, x_Euler_Method(10,:), '--g', 'linewidth', 1)


% % % % % % % 3. 4次のルンゲクッタ法による数値積分 % % % % % % %
x_RungeKutta_Method = zeros(length(M)*2, length(time_vec));
x = x0;
d = zeros(length(M)*2, 4);

for ii1 = 1:length(time_vec)
    x_RungeKutta_Method(:,ii1) = x;
    
    xt = x;
    % d1
    f = Ac*x + Bc*p;
    d(:,1) = f * dt_step;
    % d2
    x = xt + d(:,1)*0.5;
    f = Ac*x + Bc*p;
    d(:,2) = f * dt_step;
    % d3
    x = xt + d(:,2)*0.5;
    f = Ac*x + Bc*p;
    d(:,3) = f * dt_step;
    % d4
    x = xt + d(:,3);
    f = Ac*x + Bc*p;
    d(:,4) = f * dt_step;
    
    % 次のステップの状態を計算
    x = xt + (d(:,1) + 2*d(:,2) + 2*d(:,3) + d(:,4)) / 6;
end

plot(time_vec, x_RungeKutta_Method(10,:), '-r', 'linewidth', 1.5)
ylim([-1 1])
xlim([0 0.3])
legend('ode45', 'オイラー法', '4次のルンゲクッタ法')

 

外部関数：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