はじめに
これまでに、モード円適合やプロニー法（ポリレファレンス法）、偏分反復法を解説しました。本記事では、日本語では直交多項式（有理分数多項式法）、英語ではRational Fraction Polynomial Method (RFP)と呼ばれるモード同定手法について解説します。
私の愛読している振動工学の書籍「モード解析入門」ではモード同定法についての内容が不十分です。そこで、RFP法が紹介されている和書と洋書を紹介いたします。
和書
著者 長松先生他の「音・振動のモード解析と制御」の4.4.3項の”直交多項式法”にRFP法が記載されています。ただし、この参考書の説明では私は理解できなかったです。参考にした章は東京都立大学の吉村先生の解説となります。ちなみに、吉村先生は長松先生の教え子で、日本ではじめてモード同定の研究をした（モード同定手法で博士論文を書いた）と言われています。
洋書
著書 Ewinの「Modal Tesiting」の4.4.3項の”Method II: Rational Fraction Polynomial Method (RFP)”でRFP法が記載されています。
モード解析入門について他の章も解説しているので、過去の記事を読んでいない方、この章や節以外を復習したい方は下記記事を参照ください。
下記記事が目次となっており、各章や節にアクセスできるようにしています。
4.4.3 直交多項式法（RFP法）
※本記事は「Modal Tesiting」と「音・振動のモード解析と制御」の数式と対応させています。「Modal Testing」から引用の数式は末尾にMT、「音・振動のモード解析と制御」はOTをつけています。文献をまたがって数式を引用しているので、ところどころ記号の整合性が取れていない箇所があるかもしれませんがご了承ください。
RFPの特徴は、IIRフィルタのように伝達関数Hを下式のように表現します。
$$
H(ω)= \frac{ b_0 + b_1(iω) + b_2(iω)^2 + \cdots + b_{2N-1}(iω)^{2N-1} }{ a_0 + a_1(iω) + a_2(iω)^2 + \cdots + a_{2N}(iω)^{2N} } \tag{4.38MT}
$$
伝達関数Hが式(4.38MT)のように表現できるのは良いとして、anとbnの係数をどのように実測の伝達関数から求めるのか？
ここからがモード同定手法の難しいところになります。
まず、式(4.38)の係数anとbnを暫定的に決定したとします。すると、角周波数ωkの実測の伝達関数の値をHkとすると、誤差ekは式(4.39aMT)のようになります。
$$
ek= \frac{ b_0 + b_1(iω_k) + b_2(iω_k)^2 + \cdots + b_{2N-1}(iω_k)^{2N-1} }{ a_0 + a_1(iω_k) + a_2(iω_k)^2 + \cdots + a_{2N}(iω_k)^{2N} }-H_k \tag{4.39aMT}
$$
今後の計算の都合上、便宜的に式(4.39aMT)を式(4.39bMT)に変形します。
$$
ek’= (b_0 + b_1(iω_k) + b_2(iω_k)^2 + \cdots + b_{2N-1}(iω_k)^{2N-1}) – H_k ( a_0 + a_1(iω_k) + a_2(iω_k)^2 + \cdots + a_{2N}(iω_k)^{2N} ) \tag{4.39bMT}
$$
式(4.39bMT)を行列で表現すると式(4.41aMT)となります。
$$
ek’=\left( \begin{array}{ccccc}
1 & (iω_k) & (iω_k)^2 & \cdots & (iω_k)^{2N-1}
\end{array} \right)
\left( \begin{array}{c}
b_0 \\
\vdots \\
b_{2N-1}
\end{array} \right) -H_k
\left( \begin{array}{ccccc}
1 & (iω_k) & (iω_k)^2 & \cdots & (iω_k)^{2N-1}
\end{array} \right)
\left( \begin{array}{c}
a_0 \\
\vdots \\
a_{2N-1}
\end{array} \right)-Hk(iω_k)^{2N} a_{2m} \tag{4.41aMT}
$$
これはωkのときの誤差ekなので、周波数帯域で適用できるように行列表記にすると、下式になります。なお、これ以降は1982年のIMACの論文から式を引用します。
$$
e=\left( \begin{array}{ccccc}
1 & (iω_1) & (iω_1)^2 & \cdots & (iω_1)^{2N-1} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
1 & (iω_L) & (iω_L)^2 & \cdots & (iω_L)^{2N-1}
\end{array} \right)
\left( \begin{array}{c}
b_0 \\
\vdots \\
b_{2N-1}
\end{array} \right) –
\left( \begin{array}{ccccc}
H_1 & H_1(iω_k) & H_1(iω_k)^2 & \cdots & H_1(iω_k)^{2N-1} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
H_L & H_L(iω_k) & H_L(iω_k)^2 & \cdots & H_L(iω_k)^{2N-1}
\end{array} \right)
\left( \begin{array}{c}
a_0 \\
\vdots \\
a_{2N-1}
\end{array} \right)-
\left( \begin{array}{c}
H_1(iω_k)^{2N} \\
\vdots \\
H_L(iω_k)^{2N}
\end{array} \right)
\tag {4.1000}
$$
式(4.1000)を行列表記すると式(4.1001)となる。
$$
[ε]=[P][B] -[T][A] – [W] \tag{4.1001}
$$
この誤差ベクトルεを最小にするようにパラメーターを同定することになるので、下記の評価関数λを最小化することになります。（εkをベクトル表記したものがεです。）
$$
λ=[ε]^H [ε] \tag{4.1002}
$$
評価関数λを最小化するような[A]と[B]をもとめたいということになりますので、下記の評価関数λを[A]と[B]で偏微分することになります。
$$
\frac{\partial λ}{\partial [A]}=0 , \frac{\partial λ}{\partial [B]}=0
$$
上式を解くと式(4.1003)となる。
$$
\left( \begin{array}{cc}
Re([P]^H [P]) & Re([P]^H [T]) \\
Re([P]^H [T])^t & Re([T]^H [T])
\end{array} \right)
\left( \begin{array}{c}
[B] \\
[A]
\end{array} \right) =
\left( \begin{array}{c}
Re([P]^H [W]) \\
-Re([T]^H [W])
\end{array} \right) \tag{4.40MT}
$$
式(4.40MT)は行列が悪条件となり解けないという問題が生じます。そこで、式(4.40MT)を書き換える必要があります。（ここから先は「音・振動のモード解析と制御」にも記載があります）
$$
H(ω_i)=\frac{\sum_{k=0}^{m} c_k φ_{i.k} } {\sum_{k=0}^{n} d_k θ_{i.k} } \tag{4.60OT}
$$
ここからが重要です。
式(4.60OT)のφkr、θkrはr次の直交多項式であり、式(4.61OT)と式(4.62OT)に示すようにモードの直交性のような特徴を持つとします。
\[
\sum_{k=1}^{m} φ_{i.k}^* φ_{i.j}=
\left\{
\begin{array}{l}
0 & (k≠j)\\
0.5 & (k=j)
\end{array}
\right. \tag{4.61OT}
\]
\[
\sum_{k=1}^{m} θ_{i.k}^* |E(ω_i)|^2 θ_{i.j}=
\left\{
\begin{array}{l}
0 & (k≠j)\\
0.5 & (k=j)
\end{array}
\right. \tag{4.62OT}
\]
式(1000)と同様に誤差ベクトルを示します。
$$
ε=
\left( \begin{array}{ccc}
φ_{ 1,0 } & \ldots & φ_{ 1,m } \\
\vdots & \ddots & \vdots \\
φ_{ L,0 } & \ldots & φ_{ L,m }
\end{array} \right)
\left( \begin{array}{c}
c_0 \\
\vdots \\
c_m
\end{array} \right) –
\left( \begin{array}{ccc}
E_1 θ_{ 1,0 } & \ldots & E_1 θ_{ 1,n-1 } \\
\vdots & \ddots & \vdots \\
E_m θ_{ L,0 } & \ldots & E_m θ_{ L,n-1 }
\end{array} \right)
\left( \begin{array}{c}
d_0 \\
\vdots \\
d_{n-1}
\end{array} \right)-
\left( \begin{array}{c}
E_1 θ_{ 1,n } \\
\vdots \\
E_L θ_{ L,n }
\end{array} \right) \tag{4.67OT}
$$
$$
ε=[P][C] -[T][D] – [W] \tag{4.66OT}
$$
評価関数λを最小化するようなcとdをもとめたいということになりますので、下記の評価関数λをcとdで偏微分することになります。
$$
\frac{\partial λ}{\partial c} , \frac{\partial λ}{\partial d}
$$
上式を解くと式(4.70OT)となる。
$$
\left( \begin{array}{cc}
I & -RE([P]^H [T]) \\
RE([P]^H [T])^t & I
\end{array} \right)
\left( \begin{array}{c}
C \\
D
\end{array} \right) =
\left( \begin{array}{c}
RE([P]^H [W]) \\
0
\end{array} \right) \tag{4.70OT}
$$
式(4.70OT)の方程式はつぎのように分けて解きます。
$$
(I – (-RE([P]^H [T])^t (-RE([P]^H [T]))[D]=-(-RE([P]^H [T])^t (RE([P]^H [W])) \tag{4.72OT}
$$
$$
[C]=[H]-(-RE([P]^H [T])[D] \tag{4.73OT}
$$
参考にしたIMACの論文では式(4.70OT)をForsytheの方法で解くことを提案しています。Forsytheの方法についての詳しい展開式は参考文献をご確認ください。ここでは、式展開だけ載せておきます。
$$
P_{i,k}=(j)^k R_{i,k} \tag{4.1003}
$$
$$
R_{i,-1}^+=0 \tag{4.1004}
$$
$$
R_{i,0}^+=1/(2 \sum_{i=1}^{L} q_i)^{1/2} \tag{4.1005}
$$
$$
S_{i,k}^+=ω_ivR_{i,k-1}^+ – V_{k-1} R_{i,k-2}^+ \tag{4.1006}
$$
MATLABプログラムでの検証
まずは、共振周波数が4つあるFRFから、モード特性を推定できるかを検証しました。使用するFRFは下記です。
ちなみに、真値の共振周波数ωdrとモード減衰比ζは下表です。これをRFP法で同定したい！ということになります。
| 1次共振 | 2次共振 | 3次共振 | 4次共振 | |
| ωdr | 34.73 | 100.00 | 153.21 | 187.94 |
| ζ | 4.00% | 4.00% | 4.00% | 4.00% |
同定結果が下記です。ちょっとだけ安定化ダイアグラムっぽく図示してみました。同定した共振周波数ωdrとモード減衰比ζは下表です。
| 1次共振 | 2次共振 | 3次共振 | 4次共振 | |
| ωdr | 34.73 | 100.00 | 153.21 | 187.94 |
| ζ | 4.00% | 4.00% | 4.00% | 4.00% |
誤差なく同定できていますね。ただ、上記結果はFRFの中に共振周波数が4つあると仮定して同定しています。実際には共振周波数の数は不明なので、仮に7としてもう一度計算してみます。
再計算したFRFの結果がいびつになっていますね。つまり、同定したい共振周波数の数によって、モード特性の結果が変わってしまいます。（詳しく説明しませんが、モード減衰比の結果もおかしくなります。）
そのため、同定したい共振周波数の数に「適正な値」を設定する必要があります。もしくは、市販のモード同定ソフトのように、極の安定性を計算し、安定しているモード特性だけを選択するとう手段があります。まあこの辺はおいおい説明していければと思っています。
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)])
ylim([min(abs(FRF_E)) max(abs(FRF_E))])
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%RFP法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=4;
[FRF_Re,modal_parameter]=Rational_Fraction_Polynomial_Method(FRF_E,freq,N);
semilogy(w,abs(FRF_Re),'r-','linewidth',2)
% % % % % 適当に安定化ダイアグラムっぽいものを書いてみました
Y=logspace(log10(min(abs(FRF_E))),log10(max(abs(FRF_E))),N);
for ii1=1:1:N
[FRF_Re,modal_parameter]=Rational_Fraction_Polynomial_Method(FRF_E,freq,ii1);
semilogy(modal_parameter(:,1),Y(ii1)*ones(ii1,1),'bo','linewidth',2)
semilogy([w(1) w(end)],Y(ii1)*[1 1],'k--','linewidth',1)
end
functionファイル
function [FRF_recalcu,modal_parameter]=Rational_Fraction_Polynomial_Method(FRFdata,freq,N)
% modal_parameter = Modal Parameters [freq,damp,Ci,Oi]:
% fn：共振周波数
% damp：モード減衰比
% Ci：モード振幅
% Oi：モード位相角
% % % % % % % % % % % % % % % % % 初期設定
w=2*pi*freq;
[r,c]=size(w);
if r<c
w=w.';
end
[r,c]=size(FRFdata);
if r<c
FRFdata=FRFdata.';
end
% % 角周波数の正規化
nom_w=max(w);
w=w./nom_w;
% % % % % % % % % % % % % % % % % % % % % % % % % %
m=2*N-1; %式(4.60OT)のm
n=2*N; %式(4.60OT)のn
%orthogonal function that calculates the orthogonal polynomials
[phi,coeff_A]=Complex_Orthogonal_Polynomials(FRFdata,w,1,m);
[theta,coeff_B]=Complex_Orthogonal_Polynomials(FRFdata,w,2,n);
[r,c]=size(phi);
P=phi(:,1:c); %式(4.66OT)のP
[r,c]=size(theta);
Theta=theta(:,1:c); %式(4.67OT)のθ
T=sparse(diag(FRFdata))*theta(:,1:c-1);
W=FRFdata.*theta(:,c);%式(4.66OT)のW
X=-2*real(P'*T); %式(4.70OT)の最左の行列(1,2)
G=2*real(P'*W); %式(4.70OT)の最右のベクトル(1,1)
D=-inv(eye(size(X))-X.'*X)*X.'*G; %式(4.72OT)
C=G-X*D; %式(4.73OT)
D=[D;1];
%calculation of FRF (alpha)
FRF_recalcu=zeros(length(w),1);
for n=1:length(w),
numer=sum(C.'.*P(n,:));
denom=sum(D.'.*Theta(n,:));
FRF_recalcu(n)=numer/denom;
end
A=coeff_A*C;
[r,c]=size(A);
A=A(r:-1:1).'; %{A} standard numerator polynomial coefficients
B=coeff_B*D;
[r,c]=size(B);
B=B(r:-1:1).'; %{B} standard denominator polynomial coefficients
%calculation of the poles and residues
[R,P,K]=residue(A,B);
[r,c]=size(R);
for n=1:(r/2),
Residuos(n,1)=R(2*n-1);
Polos(n,1)=P(2*n-1);
end
[r,c]=size(Residuos);
Residuos=Residuos(r:-1:1)*nom_w; %residues
Polos=Polos(r:-1:1)*nom_w; %poles
fn=abs(Polos); %Natural frequencies (rad/sec)
damp=-real(Polos)./abs(Polos); %Damping ratios
Ai=-2*(real(Residuos).*real(Polos)+imag(Residuos).*imag(Polos));
Bi=2*real(Residuos);
const_modal=complex(Ai,abs(Polos).*Bi);
Ci=abs(const_modal); %Magnitude modal constant
Oi=angle(const_modal).*(180/pi); %Phase modal constant (degrees)
modal_parameter=[fn, damp, Ci, Oi]; %Modal Parameters
end
% % % % % % % % % % % % % % % % % % % % %
function [P,coeff]=Complex_Orthogonal_Polynomials(FRFdata,w,phitheta,kmax)
if phitheta==1
q=ones(size(w)); %分子項（bの項）の重み関数
elseif phitheta==2
q=(abs(FRFdata)).^2; %分母項（aの項）の重み関数
end
R_m1=zeros(size(w)); %式(4.1004)
R_0=1/sqrt(2*sum(q)).*ones(size(w)); %式(4.1005)
Rminus=[R_m1,R_0]; %polynomials Half Function
coeff=zeros(kmax+1,kmax+2);
coeff(1,2)=1/sqrt(2*sum(q));
%generating orthogonal polynomials matrix and transform matrix
Rplus=[];
R=[Rminus Rplus];
for k=1:kmax,
Vkm1=2*sum(w.*R(:,k+1).*R(:,k).*q); %IMACの論文の(29)式
Sk=w.*R(:,k+1)-Vkm1*R(:,k); %式(4.1006)
Dk=sqrt(2*sum((Sk.^2).*q)); %IMACの論文の(31)式
Rplus=[Rplus Sk/Dk];
coeff(:,k+2)=-Vkm1*coeff(:,k);
coeff(2:k+1,k+2)=coeff(2:k+1,k+2)+coeff(1:k,k+1);
coeff(:,k+2)=coeff(:,k+2)/Dk;
R=[Rminus Rplus];
end
R=R(:,2:kmax+2); %orthogonal polynomials matrix
coeff=coeff(:,2:kmax+2); %transform matrix
P=R.*i.^(0:kmax); %式(4.1003)
jk=i.^(0:kmax);
coeff=(jk'*jk).*coeff;
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
コメント