MATLABで解くN自由度振動系の状態空間モデル｜変位から加速度まで自由自在に抽出する方法

はじめに：振動工学と制御工学が融合するエキサイティングな時代

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

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

このトレンドの核となる技術の一つが、制御工学でおなじみの「カルマンフィルタ」です。カルマンフィルタを用いて「状態変数 ＝ 振動（変位や速度）」を推定することで、直接センサーを置けない場所の振動を予測でき、デジタルツインや状態監視にダイレクトに応用できます。
海外では、制御工学の枠にとどまらず、振動工学が扱うような高周波帯域や多自由度の大規模モデルへカルマンフィルタを適用する研究が盛んに進められています。

カルマンフィルタの理論構築にあたっては、海外のこちらの論文や、MATLABプログラムが豊富に掲載されている下記の専門書が非常に参考になります。

本記事では、このカルマンフィルタを自在に操るための「前段階（必須知識）」として、振動工学の数式を制御工学のフォーマットに変換する「状態空間モデル」の作り方と、連続時間系での解き方について分かりやすく解説します。

状態空間モデルへの変換（N自由度のバネ-マス-ダンパモデル）

振動工学や実験モード解析を中心に勉強してきた方にとって、「状態空間モデル」という言葉は少し馴染みがないかもしれません。状態空間モデルとは、システム（構造物）の内部状態を表す「状態方程式」と、観測したいデータを抽出する「出力方程式」の2つの式で構成されるモデルです。

ここでは、図1に示す $N$ 自由度のバネ-マス-ダンパモデルを用いて式展開をしていきます。

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

このモデルの運動方程式は、変位を $y(t)$、入力を $u(t)$ とすると、式(1)のように表すことができます。（$[ ]$ は行列およびベクトルを表します）

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

式(1)の運動方程式を、変位 $x_1$ と速度 $x_2$ という2つの変数に分割します。

$$ \begin{cases} [x_1(t)] = [y(t)] \\ [x_2(t)] = [\dot{y}(t)] \end{cases} \tag{2} $$

この関係式を用いて、式(1)を1階の微分方程式に書き換えると、以下の「状態方程式」（式(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] \tag{3} $$

次に、このシステムから「変位データだけを取り出したい」とします。そのためのレンズ（フィルター）の役割を果たすのが「出力方程式」（式(4)）です。

$$ [y] = \begin{bmatrix} [I] & [0] \end{bmatrix} \begin{bmatrix} [x_1] \\ [x_2] \end{bmatrix} \tag{4} $$

式(3)と式(4)の行列部分をひとまとめにして一般化すると、よく見る式(5)の形になります。これが状態空間モデルの基本形です。

$$ \begin{cases} [\dot{x}(t)] = [A][x(t)] + [B][u(t)] \quad \text{（状態方程式）} \\ [y(t)] = [C][x(t)] \quad \text{（出力方程式）} \end{cases} \tag{5} $$

ここで、このモデルの素晴らしい汎用性に注目してください。

• 状態方程式： システムの物理的性質（質量・剛性・減衰）と入力の関係を表す（不変）。
• 出力方程式： 内部状態（変位と速度）から、自分が欲しいデータだけを自由に抽出する（自由自在に変更可能）。

例えば、センサーで「加速度」を計測している場合、出力方程式で加速度を取り出さなければなりません。
加速度 $[\ddot{y}]$ は、式(1)を変形すると $[\ddot{y}] = [\dot{x}_2] = -[M]^{-1}[K][x_1] – [M]^{-1}[C][x_2]$ となります（外力がゼロの場合）。
したがって、出力方程式の $[C]$ 行列を式(6)のように書き換えるだけで、一瞬にして加速度を出力するモデルに早変わりします。

$$ [y] = \begin{bmatrix} -[M]^{-1}[K] & -[M]^{-1}[C] \end{bmatrix} \begin{bmatrix} [x_1] \\ [x_2] \end{bmatrix} \tag{6} $$

これが状態空間モデルの最大の強みです。

本来、カルマンフィルタなどのデジタル制御に適用するためには、この連続時間系のモデルを「離散化」する必要があります。しかし、まずは連続時間系のまま微分方程式を解くイメージを掴むため、MATLABの強力な常微分方程式ソルバー ode45 を使ってシミュレーションを行ってみましょう。

MATLABで解く状態空間モデル（N自由度モデル）

では、例題として図1のモデルを用いて、以下のパラメータで計算してみます。

【解析条件】

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

全自由度の変位と加速度が計算されますが、ここでは一番端にある 質量10 ($m_{10}$) の変位と加速度の結果をグラフ化して紹介します。

図2 m10の変位と加速度の結果

ode45 に状態方程式の関数を渡すだけで、非常に簡単かつ高精度に振動の時刻歴応答（変位と加速度）が計算できていることがわかります。

本記事では、カルマンフィルタを学ぶための準備として、状態空間モデルの考え方と連続時間系での解法について解説しました。次回以降のステップとして、このモデルの「離散化」へと進んでいきます。
今回使用した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)];

% 出力方程式の行列C（加速度を出力する場合）
Cc = [-inv(M)*K -inv(M)*C]; 

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

% MATLABの ode45 で連続時間系モデルを解く
[T,Y] = ode45(@state_equation, time, x0); 

figure
% 変位のプロット
subplot(2,1,1)
plot(T, Y(:,10), 'k', 'linewidth', 1.5)
xlabel('Time [s]')
ylabel('Displacement [m]')
title('質量10 (m_{10}) の変位応答')

% 加速度のプロット（出力方程式から算出）
subplot(2,1,2)
y_k = Cc * Y.'; % 行列演算で全時刻の加速度を一括計算
plot(T, y_k(10,:), 'r', 'linewidth', 1.5)
xlabel('Time [s]')
ylabel('Acc. [m/s^2]')
title('質量10 (m_{10}) の加速度応答')

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