はじめに
近年の海外大学のトレンドとして、純粋に振動工学や機械工学を深堀するような研究はされておらず、デジタルツインや状態監視、1DCAE、MBD、メタマテリアルといった研究やワードが目立つようになってきています。これらの分野は、振動工学＋制御工学＋IoT＋システムといったように、ひと昔前だと他分野と見なされた分野を横通しする必要があります。そのため、1人の教授が研究室運営をすることが多い日本の大学では先端の研究ができていないのが現状です。特に振動分野では上記研究テーマで優れた研究がなされていません。
今回紹介する「カルマンフィルタ」は状態変数を推定する手法で、状態変数＝振動とすると振動推定が可能になります。そのため、デジタルツインや状態監視などの分野で適用されています。
カルマンフィルタ自体は制御工学で非常に広く知られていて、制御工学では結構やりつくされた感がある理論です。ただし、制御工学では数自由度のモデルでの適用であり、振動工学で扱うような周波数帯域を扱えなかったのですが、海外大学では振動工学で扱うような周波数帯域への適用を目指しています。
上記のようなトレンドがありますが、振動工学を生業としている日本の方にはあまりなじみがないのではないでしょうか？そもそもトレンドを捉えていないのではないでしょうか？
日本の大学の文献はあまりいいのがないので、本記事では本リンクの論文を参考にしたいと思います。カルマンフィルタについては下記の著書を参考にしています。下記著書にはカルマンフィルタなどのMATLABプログラムが掲載されています。
なお、本記事ではカルマンフィルタを学習する前段階として、状態空間モデルを連続時間系で解く方法について解説したいと思います。前回の記事を読んでいない方は下記も参考にしてください。
状態空間モデル（N自由度のバネ-マス-ダンパモデル）
前回の記事と重複するので、簡単におさらい部分だけ記載します。図1のようなN自由度の運動方程式は式(1)でしたね。なお、細かい文字式についての説明は前回の記事を参照ください。
図1 バネマスモデル（N自由度）
$$
[M][\ddot{y} (t)] + [C][\dot{y} (t)] + [K][y(t)]=[u(t)] \tag{1}
$$
式(1)の状態方程式は式(3)になり、出力方程式は式(4)になりましたね。
$$
\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}
$$
$$
[y]=
\left( \begin{array}{cc}
[I] & [0]
\end{array} \right)
\left( \begin{array}{c}
[x_1] \\
[x_2]
\end{array} \right) \tag{4}
$$
一般化すると式(5)のようになりましたね。
\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}
状態空間モデル（状態方程式と出力方程式）は連続時間系なので、カルマンフィルタを適用するためには離散化する必要があります。ただし、離散化の話をしても、連続時間系を解く方法がわかっていないと離散化についての理解ができないと思うので、連続時間系の解法について説明します。
MATLABには連続時間系の微分方程式を解く「ode45」とオイラー法および4次のルンゲクッタ法を比較したいと思います。
オイラー法
オイラー法はΔt後の状態x(t+Δt)を下記のように近似計算する方法です。
$$
[x(t+Δt)]=[x(t)]+([A][x(t)]+[B][u(t)])Δt \tag{8}
$$
ルンゲクッタ法
ルンゲクッタ法はΔt後の状態x(t+Δt)をtのまわりでテイラー級数展開して近似計算する方法です。よく使われる方法としては、Δtの4次までの項を考慮して5次以降を無視した形が使われます。
$$
[x(t+Δt)]=[x(t)]+\frac{d_1+2d_2+2d_3+d_4}{6} \tag{9}
$$
ただし、
$$
\dot{x}(t)=f(x(t))
$$
\begin{eqnarray}
\left\{
\begin{array}{l}
d_1=f(x(t))Δt \\
d_2=f(x(t)+d_1/2)Δt \\
d_3=f(x(t)+d_2/2)Δt \\
d_4=f(x(t)+d_3)Δt
\end{array}
\right. \tag{10}
\end{eqnarray}
MATLABで解く状態空間モデル（N自由度モデル）
では、前回の例題を使って「ode45」と比較してみましょう。
- 自由度N=10
- k=10^5
- c=1
- m=1
- サンプリング周波数fs=400（dt=1/fs）
- 初期値のx：m1の変位を1、m2～m10の変位を0
- 求めたい状態変数：加速度
全自由度の変位と速度が求まりますが、ここではm10の変位の結果だけを紹介します。
図2 解法の比較
MATLABの「ode45」が正解だとすると、4次のルンゲクッタ法はそこそこ精度よく解析できていることがわかりますね。しかしオイラー法は振幅が全然一致しておらず、本モデルのような解析では使い物になりそうにないですね。
本記事では簡単に状態空間モデルを連続時間系のまま解く方法について説明しました。使ったプログラムは下記に添付いたします。
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;
% % % % % % % % % % % MATLABのode45を使った場合
[T,Y] = ode45(@state_equation, time, x0); %理論値
figure
% subplot(211)
plot(T,Y(:,10),'k','linewidth',2)
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]')
% % % % % % % % % % % オイラー法
time=T.'; dt=T(2)-T(1);
p=zeros(size(M,1),1);
x_Euler_Method=zeros(length(M)*2,length(time));
x=x0;
for ii1=1:1:length(time)
x_Euler_Method(:,ii1)=x;
dx=Ac*x+Bc*p;
x = x + dx*dt;
end
% subplot(211)
hold on
plot(time,x_Euler_Method(10,:),'--g','linewidth',1)
% % % % % % % % % % %
% % % % % % % % % % % ルンゲクッタ法
p=zeros(size(M,1),1);
x_RungeKutta_Method=zeros(length(M)*2,length(time));
x=x0;
d=zeros(length(M)*2,4);
for ii1=1:1:length(time)
x_RungeKutta_Method(:,ii1)=x;
% % % % % % % d1
xt=x;
f = Ac*x + Bc*p;
d(:,1) = f*dt;
% % % % % % % d2
x = xt + d(:,1)*0.5;
f = Ac*x + Bc*p;
d(:,2) = f*dt;
% % % % % % % d3
x = xt + d(:,2)*0.5;
f = Ac*x + Bc*p;
d(:,3) = f*dt;
% % % % % % % d3
x = xt + d(:,3);
f = Ac*x + Bc*p;
d(:,4) = f*dt;
x = xt + (d(:,1) + d(:,2)*2 +d(:,3)*2 +d(:,4) )/6;
end
% subplot(211)
hold on
plot(time,x_RungeKutta_Method(10,:),'-r','linewidth',2)
ylim([-1 1])
xlim([0 0.3])
legend('ode45','オイラー法','4次のルンゲクッタ法')
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
コメント