MATLABで学ぶモード同定法　直交多項式法（RFP法）

はじめに

これまでに、モード円適合やプロニー法（ポリレファレンス法）、偏分反復法を解説しました。本記事では、日本語では直交多項式（有理分数多項式法）、英語ではRational Fraction Polynomial Method (RFP)と呼ばれるモード同定手法について解説します。

私の愛読している振動工学の書籍「モード解析入門」ではモード同定法についての内容が不十分です。そこで、RFP法が紹介されている和書と洋書を紹介いたします。

和書

著者 長松先生他の「音・振動のモード解析と制御」の4.4.3項の”直交多項式法”にRFP法が記載されています。ただし、この参考書の説明では私は理解できなかったです。参考にした章は東京都立大学の吉村先生の解説となります。ちなみに、吉村先生は長松先生の教え子で、日本ではじめてモード同定の研究をした（モード同定手法で博士論文を書いた）と言われています。

洋書

著書 Ewinの「Modal Tesiting」の4.4.3項の”Method II: Rational Fraction Polynomial Method (RFP)”でRFP法が記載されています。

 

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

下記記事が目次となっており、各章や節にアクセスできるようにしています。

長松昭男著「モード解析入門」のMATLABコード付き解説

長松昭男先生の「モード解析入門」をMATLABコード付きで解説しています。

mech-eng-uni.com

2022.01.08

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プログラム

実行ファイル

functionファイル