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

2022.08.29

はじめに

近年の海外大学のトレンドとして、純粋に振動工学や機械工学を深堀するような研究はされておらず、デジタルツインや状態監視、1DCAE、MBD、メタマテリアルといった研究やワードが目立つようになってきています。これらの分野は、振動工学＋制御工学＋IoT＋システムといったように、ひと昔前だと他分野と見なされた分野を横通しする必要があります。そのため、1人の教授が研究室運営をすることが多い日本の大学では先端の研究ができていないのが現状です。特に振動分野では上記研究テーマで優れた研究がなされていません。

今回紹介する「カルマンフィルタ」は状態変数を推定する手法で、状態変数＝振動とすると振動推定が可能になります。そのため、デジタルツインや状態監視などの分野で適用されています。

カルマンフィルタ自体は制御工学で非常に広く知られていて、制御工学では結構やりつくされた感がある理論です。ただし、制御工学では数自由度のモデルでの適用であり、振動工学で扱うような周波数帯域を扱えなかったのですが、海外大学では振動工学で扱うような周波数帯域への適用を目指しています。

上記のようなトレンドがありますが、振動工学を生業としている日本の方にはあまりなじみがないのではないでしょうか？そもそもトレンドを捉えていないのではないでしょうか？

日本の大学の文献はあまりいいのがないので、本記事では本リンクの論文を参考にしたいと思います。カルマンフィルタについては下記の著書を参考にしています。下記著書にはカルマンフィルタなどのMATLABプログラムが掲載されています。

 

なお、本記事ではカルマンフィルタを学習する前段階として、状態空間モデルの解き方を解説したいと思います。

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

振動工学や実験モード解析を中心に勉強・研究している方は「状態空間モデル」についてなじみがないと思いますので、ここで簡単に説明いたします。状態空間モデルとは、時系列分析に用いられるモデルの一種で、状態方程式と出力方程式とで表されます。本記事では図1のN自由度のバネ-マス-ダンパモデルを用いて式展開をしていきます。

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

図1のバネ-マス-ダンパモデルの運動方程式は入力u(t)、変位y(t)とすると式(1)のように表すことができます。なお、[]は行列およびベクトルを表します。

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

式(1)の運動方程式を式(2)のように変位x1と速度x2という2変数に分けます。（ここだけ、時刻(t)を省略していませんが、特に意味はありません。）

\begin{eqnarray}
\left\{
\begin{array}{l}
[x_1(t)] = [y(t)] \\
[x_2(t)] = [ \dot{y} (t)]
\end{array}
\right.
\tag{2}
\end{eqnarray}

式(2)を用いて状態方程式は式(3)のように表せます。

$$
\left( \begin{array}{c}
[ \dot{x_1}] \\
[ \dot{x_2}]
\end{array} \right) =
\left( \begin{array}{cc}
[0]  & [I] \\
-[M]^{-1}[K] & -[M]^{-1}[C]
\end{array} \right)
\left( \begin{array}{c}
[x_1] \\
[x_2] [u]
\end{array} \right) +
\left( \begin{array}{c}
[0] \\
[M]^{-1}
\end{array} \right)  \tag{3}
$$

出力方程式は式(4)になります。振動関係の出力方程式では式(4)を扱うことが多いですが、たまたま出力方程式が変位x1になったということを覚えておいてください。（今後扱う問題では、例えばm1だけの加速度を出力するような出力方程式扱うことになります。）

$$
[y]=
\left( \begin{array}{cc}
[I] & [0]
\end{array} \right)
\left( \begin{array}{c}
[x_1] \\
[x_2]
\end{array} \right) \tag{4}
$$

式(3)と式(4)をまとめると式(5)のようになります。ここで、式(5)を文章で整理すると、下記がわかります。

  • 状態方程式：入力uと状態ベクトルxとの関係
  • 出力方程式：状態ベクトルxと出力yとの関係

つまり、出力方程式は状態ベクトルから求められる形であればどのような式形態でもよいということになり、「式(4)だけが運動方程式の出力方程式ではない」ということがわかるかと思います。

\begin{eqnarray}
\left\{
\begin{array}{l}
[\dot{x} (t)] = [A][x(t)]+[B][u(t)] \\
[y(t)] = [C] [x(t)]
\end{array}
\right.
\tag{5}
\end{eqnarray}

 

では、出力方程式を変位ではなく加速度にするにはどのようになるでしょうか？

答えは式(6)です。

$$
[y]=
\left( \begin{array}{cc}
-[M]^{-1}[K] & -[M]^{-1}[C]
\end{array} \right)
\left( \begin{array}{c}
[x_1] \\
[x_2]
\end{array} \right) \tag{6}
$$

本当は、状態空間モデル（状態方程式と出力方程式）は連続時間系なので、数値解析で解けるように離散化する必要がありますが、MATLABには連続時間系の微分方程式を解く「ode45」というfunctionがあるので、本日の理論解説はここまでとします。

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

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

  • 自由度N=10
  • k=10^5
  • c=1
  • m=1
  • サンプリング周波数fs=400（dt=1/fs）
  • 初期値のx：m1の変位を1、m2～m10の変位を0
  • 求めたい状態変数：加速度

全自由度の加速度が求まりますが、ここではm10の結果だけを紹介します。

図2　m10の加速度の結果

 

本記事では簡単に状態空間モデルの解法について説明いたしました。使ったプログラムは下記に添付いたします。

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;
t0 = 0;
dt=1/400;
tf = 10;
time=t0:dt:tf;

[T,Y] = ode45(@state_equation, time, x0); %理論値
figure
subplot(211)
plot(T,Y(:,10))
xlabel('Time [s]')
ylabel('Displacement [m]')
subplot(212)
y_k=Cc*Y.';
plot(T,y_k(10,:))
xlabel('Time [s]')
ylabel('Acc. [m/s^2]')

 

functionファイル

function [M]=eval_Mmatrix(m_vec)

M=diag(m_vec);

end

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

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;

% [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=[-Sa*inv(M)*K -Sa*inv(M)*C];
% Dc=Sa*inv(M);
p=zeros(size(M,1),1);
% size(p)
% size(Bc)
% size(Bc*p)
% 
% size(Ac)
% size(Ac*x)
dx=Ac*x+Bc*p;
end

 

 

コメント