MATLABで学ぶモード同定法偏分反復法 - MATLABパイセンが教える振動・騒音・音響・機械工学

はじめに
4.4.1 偏分反復法（モード解析入門では 6.3.1項）
例題
MATLABプログラム
1. 関連

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。

再復習で学んだことをMATLABコードを示しながら解説したいと思います。今回は6.3.1項の「偏分反復法」P.362から復習します。ただし、モード解析入門の「偏分反復法」の解説は内容が非常に薄く、偏微分方程式も記載されていないので、今回は著者長松先生他の「音・振動のモード解析と制御」を参考にしたいと思います。参考にした章は東京都立大学の吉村先生の解説となります。ちなみに、吉村先生は長松先生の教え子で、日本ではじめてモード同定の研究をした（モード同定手法で博士論文を書いた）と言われています。

偏分反復法は実験モード解析およびモード同定法の基本となる考え方です。様々なモード同定法の手法がありますが、留数の同定に関しては偏分反復法の理論を採用しているものが多いです。（ただし、モード同定法は実験モード解析ソフトでブラックボックス化されているので、ソフトウェアの技術サポートに専門的な質問して回答してもらう必要があります。）

モード解析入門について他の章も解説しているので、過去の記事を読んでいない方、この章や節以外を復習したい方は下記記事を参照ください。

参考書解説一覧

MATLABとNVH(Noise, Vibration and Harshness)をこよなく愛する人生の先輩、それが私（MATLABパイセン）です。NVHを極めるにあたり、周辺の知識も網羅的に勉強することになり、その知識を共有すべく、本HPを運営しています。日本の製造業を応援すべく、機械エンジニアの「車輪の再開発」の防止と業務効率化の手助けをモットーに活動。専門である振動騒音工学や音響工学について、プログラムを示しながら解説。大学の授業よりもわかりやすく説明することを目指す。質問頂ければ、回答したいと思います。

私自身もカーブフィットソフトを作成しました。ご興味のある方は下記でダウンロードできますので、お試しください。
👉無料ダウンロード
👉マニュアルページ

4.4.1 偏分反復法（モード解析入門では 6.3.1項）

これまでの記事で、FRFや振動形状はモード特性（共振周波数＆振動モードなど）の重ね合わせで表現できるということを説明しましたね。この考え方がモード解析である、ということは理解頂けていると思います。そして、FRFからモード特性を抽出することをモード同定と呼びます。

モード円適合（カーブフィット）や半値幅法は、図1の赤色の周波数帯のように隣り合う共振周波数が十分離れており、1自由度系と見なせる場合にのみ適用可能な手法でした。

図1 FRF

一方、偏分反復法は多自由度系のモード同定手法です。つまり、図1の0～100Hzのデータのように複数の共振周波数を含んでいても、それぞれのモード特性を抽出することができる手法です。

前回の記事でプロニーの方法を説明しましたが、プロニーの方法は時間領域（インパルス応答データ）の手法でした。今回説明する偏分反復法は周波数領域の手法です。

ここで、点iと点jの間の周波数応答関数を考えます。周波数応答関数は式(4.6)のように表すことができます。

$$
G_{ij}(ω)=\sum_{r=1}^{n} \left( \frac{ a_r }{ j(ω-s_r)} + \frac{ a^*_r }{ j(ω-s^*_r)} \right) +\frac{Y}{ω^2}+Z \tag{4.6}
$$

arは留数(residu)、srは固有値を表します。srの固有値はモード減衰率σrと減衰固有振動数ωdrを用いて式(4.2)となります。

$$
s_r=-σ_r+jω_{dr} \tag{4.2}
$$

式(4.6)の留数arは複素数とすると、式(4.7)のように実部と虚部に分解することができます。

$$
a_r=U_r+jV_r \tag{4.7}
$$

すると、式(4.6)は式(4.38)に書き換えることができます。

$$
G_{ij}(ω)=\sum_{r=1}^{n} \left( \frac{ U_r+jV_r }{ σ_r + j(ω-ω_{dr})} + \frac{ U_r-jV_r }{ σ_r + j(ω+ω_{dr})} \right) +\frac{Y}{ω^2}+Z \tag{4.38}
$$

なお、Y/ω^2は対称の共振周波数（モード同定しようとしている最低次数の共振周波数）よりも低い周波数を補正するためのパラメーターであり、これを剰余質量と呼びます。一方、Zは対称の共振周波数よりも高い周波数を補正するためのパラメーターであり、これを剰余剛性と呼びます。

式(4.38)を実部と虚部に分けます。

実部

$$
G^R_{ij}(ω)=\sum_{r=1}^{n} \left( \frac{ U_r σ_r+jV_r(ω-ω_{dr}) }{ σ_r^2 + (ω-ω_{dr})^2} + \frac{ U_r σ_r-V_r(ω+ω_{dr}) }{ σ_r^2 + (ω+ω_{dr})^2} \right) +\frac{Y}{ω^2}+Z \tag{4.39.1}
$$

虚部

$$
G^I_{ij}(ω)=\sum_{r=1}^{n} \left( \frac{ -U_r (ω-ω_{dr})+V_rσ_r }{ σ_r^2 + (ω-ω_{dr})^2} – \frac{ U_r (ω+ω_{dr})+V_rσ_r }{ σ_r^2 + (ω+ω_{dr})^2} \right) \tag{4.39.2}
$$

上記までは周波数応答関数の理論式の展開の話です。

便宜的に実測の周波数応答関数Eを下記のように表すとします。
（周波数応答関数は複素数で表せるので、実数と虚数で分けて考えられますよね、という確認です。）

$$
E(ω)=E^R(ω)+jE^I(ω)
$$

ここからが偏分反復法の考え方になります。実測の周波数応答関数Eを理論式の周波数応答関数Gijで表現できれば、モード特性（共振周波数、モード減衰率、留数）が求めることができますよね。そこで偏分反復法ではEとGijの誤差が最小になるような、GijのUr、Vr、ωdr、σrを最小二乗法で求めます。

つまり、評価関数λ（誤差の自乗和）は式(4.40.0)になります。
$$
λ=\sum_{ω=ω_{start}}^{ω_{end}} \left( (E^R(ω)-G^R_{ij}(ω))^2 + (E^I(ω)-G^I_{ij}(ω))^2 \right) \tag{4.40.0}
$$

ただし、式(4.40.0)ではうまくいきません。というのも、共振周波数以外の応答と比べて共振周波数での応答は数百倍～数千倍になることがあります。つまり、式(4.40.0)の評価関数λは共振周波数での値だけが一致するようなUr、Vr、ωdr、σrを同定すれば良いということになってしまいます。これではダメですよね。一般的には共振周波数以外の周波数帯域も実測と一致するような周波数応答関数を同定したいですよね。

そこで、評価関数λに重み関数W(ω)を導入し、式(4.40)のようにします。

$$
λ=\sum_{ω=ω_{start}}^{ω_{end}} \left( (E^R(ω)-G^R_{ij}(ω))^2 + (E^I(ω)-G^I_{ij}(ω))^2 \right)W(ω) \tag{4.40}
$$

ここからは行列形式で表すために、以下のように表記を変更します。

$$
E=\left(
\begin{array}{c}
E^R(ω_{start}) \\
\vdots \\
E^R(ω_{end}) \\
E^I(ω_{start}) \\
\vdots \\
E^I(ω_{end})
\end{array}
\right) ,
G=\left(
\begin{array}{c}
G^R_{ij}(ω_{start}) \\
\vdots \\
G^R_{ij}(ω_{end}) \\
G^I_{ij}(ω_{start}) \\
\vdots \\
G^I_{ij}(ω_{end})
\end{array}
\right) ,
W=\begin{pmatrix}
W(ω_{start}) & & & & & \\
& \ddots & & \Huge{0} & & \\
& & W(ω_{end}) & & & \\
& \Huge{0} & & W(ω_{start}) & & \\
& & & & \ddots & \\
& & & & & W(ω_{end})
\end{pmatrix} \tag{4.41}
$$

式(4.41)の重み関数行列を用いて、評価関数λを式(4.42)のように置き換えます。

$$
λ=[\bf{E-G}]^t \bf{W} [\bf{E-G}] \tag{4.42}
$$

さて、評価関数λを定義できました。繰り返しになりますが、評価関数λは実測値と理論値の誤差です。つまり評価関数λが最小になるようなUr、Vr、ωdr、σrを求めればよいということですね。評価関数λが最小となるようなUr、Vr、ωdr、σrは、λを各パラメータで偏微分して0となるようなパラメータを求めれば良いということになります。

とりあえず、Ur、Vr、ωdr、σrを設計変数aとして下記を求めます。

$$
\frac{\partial λ}{\partial \bf{a}} = -2 [ \frac{\partial \bf{G}}{\partial \bf{a}} ]^t \bf{W} [\bf{E-G}] = \bf{0} \tag{4.44}
$$

先ほど、設計変数をまとめてaとしましたので、具体的に∂G/∂aとは？と疑問に思っている方もいるかと思います。例えば、aをUrに置き換えて、偏微分することで各パラメータの偏微分が求まります。（少し式が煩雑になるので、画像で貼らせていただきます。）

ここで、初期値の設計変数をaoとし、式(4.44)をaoにおいてテイラー展開し、2次の項以降を省略すると式(4.47)となります。

$$
\bf{G(a)}=\bf{G(a_0)}+[\frac{\partial \bf{G}}{\partial \bf{a}}]\bf{δa} \tag{4.47}
$$

式(4.47)を式(4.44)に代入すると式(4.49)が得られます。

$$
[ \frac{\partial \bf{G}}{\partial \bf{a}} ]^t \bf{W} [\bf{E-G(a_0)}] = [\frac{\partial \bf{G}}{\partial \bf{a}}]^t \bf{W} [\frac{\partial \bf{G}}{\partial \bf{a}}]\bf{δa} \tag{4.49}
$$

式(4.49)および式(4.47)は初期値aoでテイラー展開で線形近似して求めたので、初期値を設定する必要があります。またテイラー展開は初期値近傍において線形性が保たれるので、初期値aoから微小にa=ao+δaとして更新していくことになります。

さて、ここまでの数式とモードパラメータ（Ur、Vr、ωdr、σr）の同定方法をフローチャートにしましたのでご確認ください。図の線形項と非線形項は式(4.6)の分子（線形項）と分母（非線形項）に対応します。

理論の説明はここまでとして、例題を通して解説を続けます。

例題

4自由度のバネマス系でモード減衰比が4%の下図のような周波数応答関数を考えます。

上記を実測の周波数応答関数Eとします。ちなみに、真値の共振周波数ωdrとモード減衰比ζは下表です。これを偏分反復法で同定したい！ということになります。

	1次共振	2次共振	3次共振	4次共振
ωdr	34.73	100.00	153.21	187.94
ζ	4.00%	4.00%	4.00%	4.00%

偏分反復法で同定する場合、初期値を設定する必要があります。初期値のωdrはグラフから大まかに30Hz、95Hz、155Hz、190Hzとします。モード減衰はグラフから読み取れないので、とりあえず10%とします。

図のフローチャートに従ってモード特性を同定し、周波数応答関数を再計算した結果が下記です。

再計算した周波数応答関数は同定したいFRF（実測のFRF）と一致していますね。
同定したモードパラメーターを下表で示します。

	1次共振	2次共振	3次共振	4次共振
ωdr	34.70	99.91	153.03	187.77
σr	1.39	4.02	6.19	7.31
ζ	4.01%	4.03%	4.04%	3.89%
Ur	1.29E-07	3.03E-07	7.88E-06	8.28E-07
Vr	-7.49E-04	-1.67E-03	-1.41E-03	-4.79E-04

モード減衰比ζがほぼ4%になっており、共振周波数ωdrもほとんど誤差なく推定できていますね。このプログラムを使えば、高額な振動計測ソフトがなくてもモード同定できるようになります。

MATLABプログラム

実行プログラム

clear all;clc;close all

% freq=30:0.1:40; %計算周波数の定義
freq=0.1:0.1:40; %計算周波数の定義
df=freq(2) - freq(1);
w=2*pi*freq; %角振動数
m_vec=ones(1,4);
k_vec=ones(1,4)*10^4;
[M]=eval_Mmatrix(m_vec); %質量行列
[K]=eval_Kmatrix(k_vec); %剛性行列

[V,D]=eig(K,M); %固有値D、固有ベクトルV
for ii1=1:length(m_vec)
    V(:,ii1)=V(:,ii1)/V(1,ii1); % 加振点m1のモードベクトルの成分が1になるように正規化
end
wn=sqrt(diag(D)); %固有角振動数
fn=wn/(2*pi); %固有振動数

mr=diag(V.'*M*V); %モード質量行列
kr=diag(V.'*K*V); %モード剛性行列
c_cr=2*sqrt(mr.*kr);
modal_dampimg=1*0.04; %モード減衰比（とりあえず4%と設定）
cr=c_cr*modal_dampimg; %モード減衰係数

F=zeros(length(m_vec),1); F(1)=1; %外力ベクトル(全周波数帯を1で加振する)
Xj=zeros(length(m_vec),length(freq)); %変位ベクトル（計算前にベクトルを定義してメモリを確保）

% % % % モード法
ii=1; %加振点 （加振点が複数ある場合はfor ii=1:kength(F)のループを追加すればよい ）
for ii1=1:length(m_vec) %応答点
    for ii2=1:length(wn)
        Xj(ii1,:)=Xj(ii1,:)+V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
%       横ﾍﾞｸﾄﾙ　=横ﾍﾞｸﾄﾙ  +  ｽｶﾗｰ　　×ｽｶﾗｰ　　　./(　ｽｶﾗｰ×横ベクトル + ｽｶﾗｰ×横ベクトル +ｽｶﾗｰ）×ｽｶﾗｰ
    end
end

FRF_E=Xj(1,:);


figure
semilogy(w,abs(FRF_E),'k-','linewidth',2)
hold on
xlabel('角周波数 Hz')
ylabel('変位 [m]')
xlim([w(1) w(end)])
grid on


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ここからが偏分反復法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=100; %繰り返し回数
wdr=[30 95 155 190]; %初期値wdr
sigma_r=wdr*0.1; %初期値のσr,   σr=角共振周波数×モード減衰比

[Ur,Vr,Y,Z,wdr,sigma_r,modal_damp,G]=Differential_Iteration_Method(N,FRF_E,w,wdr,sigma_r);


semilogy(w,abs(G(1:length(w))+j*G(length(w)+1:end)),'r--','linewidth',2)
legend('同定したいFRF(実測のFRF)','100回探索後のFRF')

functionファイル

function [Ur,Vr,Y,Z,wdr,sigma_r,modal_damp,G]=Differential_Iteration_Method(N,FRF_E,w,wdr,sigma_r)

if N==1
    N=2;
end

lambda=zeros(1,N);
E=[real(FRF_E),imag(FRF_E)].'; % 実測のFRF行列E
W=diag(1./[abs(FRF_E) abs(FRF_E)]).^2; %重み関数行列W
Nmode=length(wdr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 線形項の計算
[dGdaLiner]=dGdaLiner_Matrix(wdr,sigma_r,w,Nmode); %式(4.45)
aLiner=inv(dGdaLiner.'*W*dGdaLiner)*dGdaLiner.'*W*E;%式(4.44)  および G=dGdaLiner*aLinerなので
G=dGdaLiner*aLiner; %周波数応答関数の再計算  
lambda(1,1)=(E-G).'*W*(E-G);

for ii1=2:N
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 非線形項の更新
    Ur=aLiner(1:Nmode); Vr=aLiner(Nmode+1:Nmode*2);
    [dGdaNONLiner]=dGdaNONLiner_Matrix(Ur,Vr,wdr,sigma_r,w,Nmode); %式(4.45)
    aNONLiner=[sigma_r wdr].';
    da=inv(dGdaNONLiner.'*W*dGdaNONLiner)*dGdaNONLiner.'*W*(E-G); %式(4.49)
    aNONLiner = aNONLiner + da; %式(4.46)
    modal_damp=aNONLiner(1:Nmode)./aNONLiner(Nmode+1:Nmode*2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 線形項の計算
sigma_r=aNONLiner(1:Nmode).';  wdr=aNONLiner(Nmode+1:Nmode*2).';
[dGdaLiner]=dGdaLiner_Matrix(wdr,sigma_r,w,Nmode); %式(4.45)
aLiner=inv(dGdaLiner.'*W*dGdaLiner)*dGdaLiner.'*W*E; %式(4.44)  および G=dGdaLiner*aLinerなので
G=dGdaLiner*aLiner; %周波数応答関数の再計算
lambda(1,ii1)=(E-G).'*W*(E-G);
end

G=dGdaLiner*aLiner;

% % % 最終結果
Ur=aLiner(1:Nmode);
Vr=aLiner(Nmode+1:Nmode*2);
Y=aLiner(Nmode*2+1);
Z=aLiner(Nmode*2+2);
sigma_r=aNONLiner(1:Nmode).';
wdr=aNONLiner(Nmode+1:Nmode*2).';
modal_damp=sigma_r./wdr;

end

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 [dGdaLiner]=dGdaLiner_Matrix(wdr,sigma_r,w,Nmode)
% % G=T*[Ur;Vr;Y;Z;Ur;Vr]となるようなTを求める

% w=2*pi*freq;
Re_T=zeros(length(w)*1,2*Nmode+2);
Im_T=zeros(length(w)*1,2*Nmode+2);

for ii1=1:length(w)
    for ii2=1:Nmode
        Ar=sigma_r(ii2)^2+(w(ii1)-wdr(ii2))^2;
        Br=sigma_r(ii2)^2+(w(ii1)+wdr(ii2))^2;
        % Ur
        Re_T(ii1,ii2)=sigma_r(ii2)*(Ar+Br)/(Ar*Br);
        % Vr
        Re_T(ii1,ii2+Nmode)=( -Ar*(w(ii1)+wdr(ii2)) + Br*(w(ii1)-wdr(ii2)) )/(Ar*Br);
        
        % Ur
        Im_T(ii1,ii2)=( -Ar*(w(ii1)+wdr(ii2)) - Br*(w(ii1)-wdr(ii2)) )/(Ar*Br);
        % Vr
        Im_T(ii1,ii2+Nmode)=sigma_r(ii2)*(-Ar+Br)/(Ar*Br);
    end
    Re_T(ii1,2*Nmode+1:end)=[1/w(ii1)^2 1];    
end
dGdaLiner=[Re_T
    Im_T];
end

function [dGdaNONLiner]=dGdaNONLiner_Matrix(Ur,Vr,wdr,sigma_r,w,Nmode)


% w=2*pi*freq;
Re_T=zeros(length(w),2*Nmode);
Im_T=zeros(length(w),2*Nmode);

for ii1=1:length(w)
    for ii2=1:Nmode
        Ar=sigma_r(ii2)^2+(w(ii1)-wdr(ii2))^2;
        Br=sigma_r(ii2)^2+(w(ii1)+wdr(ii2))^2;
        Cr=-sigma_r(ii2)^2+(w(ii1)-wdr(ii2))^2;
        Dr=-sigma_r(ii2)^2+(w(ii1)+wdr(ii2))^2;
        % σr
        Re_T(ii1,ii2)= ( Ur(ii2)*Cr-2*Vr(ii2)*(w(ii1)-wdr(ii2)) )*sigma_r(ii2)/Ar^2 +......
                       ( Ur(ii2)*Dr+2*Vr(ii2)*(w(ii1)+wdr(ii2)) )*sigma_r(ii2)/Br^2;
        % ωdr
        Re_T(ii1,ii2+Nmode)= ( Vr(ii2)*Cr+2*Ur(ii2)*(w(ii1)-wdr(ii2)) )*sigma_r(ii2)/Ar^2 +......
                             ( Vr(ii2)*Dr-2*Ur(ii2)*(w(ii1)+wdr(ii2)) )*sigma_r(ii2)/Br^2;
        % σr
        Im_T(ii1,ii2)= ( Vr(ii2)*Cr+2*Ur(ii2)*(w(ii1)-wdr(ii2)) )*sigma_r(ii2)/Ar^2 +......
                       ( -Vr(ii2)*Dr+2*Ur(ii2)*(w(ii1)+wdr(ii2)) )*sigma_r(ii2)/Br^2;
        % ωdr
        Im_T(ii1,ii2+Nmode)= ( -Ur(ii2)*Cr+2*Vr(ii2)*(w(ii1)-wdr(ii2)) )*sigma_r(ii2)/Ar^2 +......
                             ( Ur(ii2)*Dr+2*Vr(ii2)*(w(ii1)+wdr(ii2)) )*sigma_r(ii2)/Br^2;
    end
end
dGdaNONLiner=[Re_T
    Im_T];
end