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')

Figure contains 6 axes objects. Axes object 1 with title y1 contains 2 objects of type line. These objects represent Estimation, Validation. Axes object 2 with title y2 contains 2 objects of type line. These objects represent Estimation, Validation. Axes object 3 with title y3 contains 2 objects of type line. These objects represent Estimation, Validation. Axes object 4 with title u1 contains 2 objects of type line. These objects represent Estimation, Validation. Axes object 5 with title u2 contains 2 objects of type line. These objects represent Estimation, Validation. Axes object 6 with title u3 contains 2 objects of type line. These objects represent Estimation, Validation.

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

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

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

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

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

compare(Validation,sys)

Figure contains 3 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Validation (y1), sys: 99.72%. Axes object 2 contains 2 objects of type line. These objects represent Validation (y2), sys: 99.92%. Axes object 3 contains 2 objects of type line. These objects represent Validation (y3), sys: 99.99%.

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

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

Figure contains 18 axes objects. Axes object 1 with title FRF11 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 with title FRF12 contains an object of type line. Axes object 4 contains an object of type line. Axes object 5 with title FRF13 contains an object of type line. Axes object 6 contains an object of type line. Axes object 7 with title FRF21 contains an object of type line. Axes object 8 contains an object of type line. Axes object 9 with title FRF22 contains an object of type line. Axes object 10 contains an object of type line. Axes object 11 with title FRF23 contains an object of type line. Axes object 12 contains an object of type line. Axes object 13 with title FRF31 contains an object of type line. Axes object 14 contains an object of type line. Axes object 15 with title FRF32 contains an object of type line. Axes object 16 contains an object of type line. Axes object 17 with title FRF33 contains an object of type line. Axes object 18 contains an object of type line.

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

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

    0.3727
    0.8525
    1.3706

dr = 3×1

    0.0008
    0.0018
    0.0029

ms = 3×3 complex

   0.0036 - 0.0019i   0.0039 - 0.0005i   0.0021 + 0.0006i
   0.0043 - 0.0023i   0.0010 - 0.0001i  -0.0033 - 0.0010i
   0.0040 - 0.0021i  -0.0031 + 0.0004i   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

Figure contains 9 axes objects. Axes object 1 with title In1 -> Out1 contains an object of type line. Axes object 2 with title In1 -> Out2 contains an object of type line. Axes object 3 with title In1 -> Out3 contains an object of type line. Axes object 4 with title In2 -> Out1 contains an object of type line. Axes object 5 with title In2 -> Out2 contains an object of type line. Axes object 6 with title In2 -> Out3 contains an object of type line. Axes object 7 with title In3 -> Out1 contains an object of type line. Axes object 8 with title In3 -> Out2 contains an object of type line. Axes object 9 with title In3 -> Out3 contains an object of type line.

制御不安定処理

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

load HighModalDensData FRF f

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

Figure contains an axes object. The axes object contains an object of type line. This object represents G.

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

sys = tfest(G,32,32);

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

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

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent G, sys.

最初の 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)