周波数領域同定: 周波数領域データを使用したモデルの推定
この例では、周波数領域データを使用してモデルを推定する方法を説明します。周波数領域データを使用したモデルの推定と検証は、時間領域データの場合と同じように機能します。これにより、時間領域、周波数領域およびスペクトル (FRF) データを使用してモデルを推定および解析するときに大きな柔軟性を得られます。両方の領域のデータを使用してモデルを同時に推定し、これらのモデルを比較および結合することができます。時間領域データを使用して推定したモデルは、スペクトル データを使用して検証できます (この反対も当てはまります)。
周波数領域データは、非線形モデルの推定や検証には使用できません。
はじめに
周波数領域実験データは、多くのアプリケーションで一般的です。データは、周波数アナライザーを使用して、周波数応答データ (周波数関数: FRF) としてプロセスから収集できます。また、たとえば、周期または帯域制限データを処理する場合など、入力および出力のフーリエ変換 (時間領域データの FFT) を使用することがより実用的です (帯域制限連続時間信号には、ナイキスト周波数を超える周波数要素はありません)。System Identification Toolbox™ では、周波数領域 I/O データが時間領域データと同じように表示されます (すなわち、iddata オブジェクトを使用)。オブジェクトの "Domain" プロパティは、"Frequency" に設定しなければなりません。周波数応答データは、周波数の関数で複素数ベクトルまたは振幅/位相ベクトルとして表されます。ツールボックスの IDFRD オブジェクトは、FRF をカプセル化するために使用します。このような IDDATA または IDFRD オブジェクト (および Control System Toolbox の FRD オブジェクト) は、いずれの推定ルーチン (procest、tfest など) でもシームレスに使用できます。
周波数領域データの調査
それではまず、周波数領域データの一部を読み込んでみましょう。
load demofr
この MAT ファイルには、周波数 W での周波数応答データと、振幅応答 AMP および位相応答 PHA が含まれます。まず、データをよく見てください。
subplot(211), loglog(W,AMP),title('Amplitude Response') subplot(212), semilogx(W,PHA),title('Phase Response')
この実験データは、IDFRD オブジェクトに格納されます。振幅と位相を複素数値の応答に変換します。
zfr = AMP.*exp(1i*PHA*pi/180); Ts = 0.1; gfr = idfrd(zfr,W,Ts);
Ts は基になるデータのサンプル時間です。データが連続時間に対応する場合、たとえば、入力が帯域制限である場合、Ts = 0 を使用します。
メモ: Control System Toolbox™ がある場合、IDFRD オブジェクトの代わりに FRD オブジェクトを使用できます。IDFRD には、FRD オブジェクトでは利用できない不確かさの測定や外乱スペクトルなど、より多くの情報オプションがあります。
IDFRD オブジェクト gfr には、データが含まれ、さまざまな方法でプロットおよび解析できます。データを表示するには、plot または bode を使用できます。
clf
bode(gfr), legend('gfr')
周波数応答 (FRF) データを使用したモデルの推定
モデルを推定するには、gfr をデータセットとして、ツールボックスのすべてのコマンドで透過的に使用することができます。唯一の制約は、ノイズ モデルを構築できないことです。これはつまり、多項式モデルには OE (出力誤差モデル) のみが適用され、状態空間モデルの場合、K = 0 を固定する必要があることを意味します。
m1 = oe(gfr,[2 2 1]) % Discrete-time Output error (transfer function) model ms = ssest(gfr) % Continuous-time state-space model with default choice of order mproc = procest(gfr,'P2UDZ') % 2nd-order, continuous-time model with underdamped poles compare(gfr,m1,ms,mproc) L = findobj(gcf,'type','legend'); L.Location = 'southwest'; % move legend to non-overlapping location
前述のように、スペクトル データを使用すると、さまざまな種類の線形モデルを連続および離散時間領域で推定できます。これらのモデルは、時間領域データを使用して検証できます。たとえば、時間領域 I/O データセット ztime は同じシステムから収集され、m1、ms および mproc の検証に使用できます。
compare(ztime,m1,ms,mproc) %validation in a different domain
また、検証データ ztime を使用してモデルの質を確認するため、残差も調査できます。ご覧のように、残差はほぼホワイトとなります。
resid(ztime,mproc) % Residuals plot
SPAFDR を使用したデータの圧縮
周波数応答データを使用する重要な理由として、ほとんど損失なく、情報を圧縮できることがあげられます。コマンド SPAFDR は、限られた周波数で円滑化された応答データを算出できます (例: 対数空間)。gfr データが 100 の対数空間周波数値に圧縮された例を示します。同じように、元の時間領域データを圧縮することもできます。
sgfr = spafdr(gfr) % spectral estimation with frequency-dependent resolution sz = spafdr(ztime); % spectral estimation using time-domain data clf bode(gfr,sgfr,sz) axis([pi/100 10*pi, -272 105]) legend('gfr (raw data)','sgfr','sz','location','southwest')
このボード線図は、円滑化されたデータの情報が正しく処理されていることを示します。100 ポイントを含むこれらのデータは、モデルの推定に正しく使用できます。以下に例を示します。
msm = oe(sgfr,[2 2 1]);
compare(ztime,msm,m1) % msm has the same accuracy as M1 (based on 1000 points)
周波数領域 I/O データを使用した推定
測定値は、入力および出力のフーリエ変換として取得できる場合があります。システムからのこのような周波数領域データは、信号 Y および U として取得されます。loglog プロットでは、似たようなデータに見えます。
Wfd = (0:500)'*10*pi/500; subplot(211),loglog(Wfd,abs(Y)),title('The amplitude of the output') subplot(212),loglog(Wfd,abs(U)),title('The amplitude of the input')
周波数応答データは基本的に、Y と U の割合となります。周波数領域データを IDDATA オブジェクトとして収集するには、以下を実行します。
ZFD = iddata(Y, U, 'Ts', 0.1, 'Frequency', Wfd)
ここでも、時間領域データおよび周波数応答データと同じように、周波数領域応答データセット ZFD をすべての推定ルーチンのデータとして使用できます。
mf = ssest(ZFD) % SSEST picks best order in 1:10 range when called this way mfr = ssregest(ZFD) % an alternative regularized reduction based state-space estimator clf compare(ztime,mf,mfr,m1)
データ表現の変換 (時間 - 周波数)
時間および周波数領域入力/出力データセットは、FFT および IFFT を使用して、いずれかの領域に変換できます。以下のコマンドは、IDDATA オブジェクトに対応します。
dataf = fft(ztime) datat = ifft(dataf)
時間および周波数領域入力/出力データセットは、SPAFDR、SPA および ETFE を使用して、周波数応答データに変換できます。
g1 = spafdr(ztime) g2 = spafdr(ZFD); clf; bode(g1,g2)
周波数応答データは、SPAFDR および SPA を使用して、より円滑なデータ (分解能とデータが少ない) に変換することもできます。
g3 = spafdr(gfr);
周波数応答データは、コマンド IDDATA を使用して、周波数領域入力/出力信号に変換できます。
gfd = iddata(g3) plot(gfd)
連続時間周波数領域データを使用した連続時間モデルの推定
時間領域データは、離散時間サンプル データとしてのみ格納し、処理できます。周波数領域データは、連続時間データを正しく表現できるという利点があります。たとえば、クイック サンプル データであるため、あるいは入力にナイキスト周波数以上の周波数要素がないため、連続時間信号にナイキスト周波数以上の周波数情報がないと想定し、またデータが定常状態の実験から収集されたものと想定します。ここで、データの離散フーリエ変換 (DFT) は、選択した周波数で、連続時間信号のフーリエ変換でもあります。したがって、これらのデータを使用して、連続時間モデルに直接適合できます。
これは、以下の例で示すことができます。
連続時間システムを評価します。
m0 = idpoly(1,1,1,1,[1 1 1],'ts',0)
周期的入力を使用して、このシステムの定常状態のシミュレーションから得られるデータを読み込みます。収集したデータは周波数領域に変換され、CTFDData.mat ファイルに保存されています。
load CTFDData.mat dataf % load continuous-time frequency-domain data.
このデータを見てみましょう。
plot(dataf)
set(gca,'XLim',[0.1 10])
dataf を推定に使用すると、既定の設定で連続時間モデルが作成されます。状態空間:
m4 = ssest(dataf,2); % Second order continuous-time model
nb = 2 の分子係数と nf = 2 の推定分母係数をもつ多項式モデルの場合、以下を使用します。
nb = 2; nf = 2; m5 = oe(dataf,[nb nf])
ステップ応答を、真のシステム m0 およびモデル m4 と m5 の不確かさと比較します。信頼区間はプロットにパッチで示されます。
clf h = stepplot(m0,m4,m5); showConfidence(h,1) legend('show','location','southeast')
この事例では必要ありませんが、連続時間データを使用して推定を行なう場合、通常は適合の範囲を制限周波数帯域 (データをローパス フィルターで処理する) に的を絞ることをお勧めします。システムの帯域幅は約 3 rad/sで、最大 6.2 rad/s まで正弦波で励起されています。適合において的を絞る適切な周波数範囲は、[0 7] rad/s となります。
m6 = ssest(dataf,2,ssestOptions('WeightingFilter',[0 7])) % state space model
m7 = oe(dataf,[1 2],oeOptions('WeightingFilter',[0 7])) % polynomial model of Output Error structure
opt = procestOptions('SearchMethod','lsqnonlin',... 'WeightingFilter',[0 7]); % Requires Optimization Toolbox(TM) m8 = procest(dataf,'P2UZ',opt) % process model with underdamped poles
opt = tfestOptions('SearchMethod','lsqnonlin',... 'WeightingFilter',[0 7]); % Requires Optimization Toolbox m9 = tfest(dataf,2,opt) % transfer function with 2 poles
h = stepplot(m0,m6,m7,m8,m9);
showConfidence(h,1)
legend('show')
まとめ
時間、周波数およびスペクトル データをシームレスに使用して、連続および離散時間領域で各種の線形モデルを推定する方法を説明してきました。モデルは、推定したときと異なる領域で検証および比較できます。データ形式 (時間、周波数、およびスペクトル) は、fft、ifft、spafdr、spa などのメソッドを使用して相互互換が可能です。さらに、tfest、ssest、procest の各推定ルーチンを使用して、直接連続時間推定を実行できます。推定および解析にいずれかの領域のデータをシームレスで使用できることが、System Identification Toolbox の重要な機能です。
