Main Content

識別したモデルのモード解析

この例では次を使用します。

ライブ スクリプトを開く

システムの状態空間モデルを識別します。このモデルを使用して周波数応答関数およびモーダル パラメーターを計算します。この例では、System Identification Toolbox™ のライセンスが必要です。

ハンマー励起

4 kHz でサンプリングされた 3 入力/3 出力のハンマー励起データを含むファイルを読み込みます。推定には最初の 104 サンプル、モデルの品質検証には 2×104 ～ 5×104 のサンプルを使用します。サンプル レートの逆数としてサンプル時間を指定します。データを @iddata オブジェクトとして保存します。

load modaldata XhammerMISO1 YhammerMISO1 fs

rest = 1:1e4;
rval = 2e4:5e4;
Ts = 1/fs;

Estimation = iddata(YhammerMISO1(rest,:),XhammerMISO1(rest,:),Ts);
Validation = iddata(YhammerMISO1(rval,:),XhammerMISO1(rval,:),Ts,'Tstart',rval(1)*Ts);

推定データと検証データをプロットします。

plot(Estimation,Validation)
legend(gca,'show')

関数 ssest を使用して、測定された出力とモデル出力の間のシミュレーション エラーを最小化するシステムの 7 次状態空間モデルを推定します。状態空間モデルが直達を持つように指定します。

Orders = 7;
opt = ssestOptions('Focus','simulation');

sys = ssest(Estimation,Orders,'Feedthrough',true,'Ts',Ts,opt);

(精度と複雑度の間のベストなトレードオフを与えるモデル次数を見つけるために、前のコードで Orders に 1:15 を設定しました。ssest は、次数を対話的に指定することができる特異値の対数プロットを出力します。この関数ではモデル次数として 7 も推奨されています)

検証データセットでモデルの品質を検証します。適合度の正規化した平方根平均二乗誤差 (NRMSE) 測定をプロットします。モデルは、検証データの出力信号を正確に表します。

compare(Validation,sys)

モデルの周波数応答関数を推定します。出力引数なしで modalfrf を使用して関数を表示します。

[frf,f] = modalfrf(sys);
modalfrf(sys)

システムは 3 つのモードを使用してよく表されていると仮定します。3 つのモードの固有振動数、減衰比、モード形状ベクトルを計算します。

Modes = 3;
[fn,dr,ms] = modalfit(sys,f,Modes)
fn = 3×1
103 ×

    0.3727
    0.8524
    1.3705


dr = 3×1

    0.0008
    0.0013
    0.0030


ms = 3×3 complex

   0.0036 - 0.0019i   0.0036 - 0.0003i   0.0021 + 0.0007i
   0.0043 - 0.0023i   0.0009 - 0.0001i  -0.0034 - 0.0010i
   0.0040 - 0.0021i  -0.0028 + 0.0003i   0.0011 + 0.0003i

再構成された周波数応答関数を計算して表示します。振幅をデシベル単位で表します。

[~,~,~,ofrf] = modalfit(sys,f,Modes);

clf
for ij = 1:3
    for ji = 1:3
        subplot(3,3,3*(ij-1)+ji)
        plot(f/1000,20*log10(abs(ofrf(:,ji,ij))))
        axis tight
        title(sprintf('In%d -> Out%d',ij,ji))
        if ij==3
            xlabel('Frequency (kHz)')
        end
    end
end

制御不安定処理

モーダル密度が高い周波数応答測定値を含むファイルを読み込みます。このデータは、フィードバック制御を使用して平衡に維持される不安定処理に対応します。識別用に idfrd オブジェクトとしてデータを格納します。ボード線図をプロットします。

load HighModalDensData FRF f

G = idfrd(permute(FRF,[2 3 1]),f,0,'FrequencyUnit','Hz');
figure
bodemag(G)
xlim([0.01,2e3])

32 個の極、32 個の零点をもつ伝達関数を識別します。

sys = tfest(G,32,32);

モデルの周波数応答と測定された応答を比較します。

bodemag(G,sys)
xlim([0.01,2e3])
legend(gca,'show')

最初の 10 の最小減衰振動モードの固有振動数と減衰比を抽出します。結果を表に保存します。

[fn,dr] = modalfit(sys,[],10);
T = table((1:10)',fn,dr,'VariableNames',{'Mode','Frequency','Damping'})
T=10×3 table
    Mode    Frequency     Damping 
    ____    _________    _________

      1      82.764       0.011304
      2      85.013       0.015632
      3      124.04       0.025252
      4      142.04       0.017687
      5      251.46      0.0062182
      6      332.79      0.0058266
      7      401.21      0.0043645
      8      625.14      0.0039247
      9      770.49       0.002795
     10      943.64      0.0019943

参考

| | (System Identification Toolbox)