このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
加圧された変圧器からの電流信号のモデル化
この例では、測定された信号のモデル化を示します。400 kV の 3 相変圧器が加圧されたときの R 位相からの電流信号を解析します。測定はスウェーデンの Sydkraft AB 社によって実行されました。
ここでは電流信号のモデル化における関数 ar
の使用について説明します。最初に、信号のノンパラメトリック解析が実行されます。次に、妥当なモデル次数を選択するためのツールについて論じ、信号のモデル化における ar
の使用について説明します。選択した高調波の範囲だけにモデルを適合させる方法についても説明します。
はじめに
信号は自己回帰線形モデルのインパルス応答と見なせるため、ar
などのツールを使ってモデル化できます。
信号のデータは、オブジェクトの出力データを信号値に設定し、入力を空のままにしておくことで、iddata
オブジェクトにカプセル化することができます。たとえば x(t)
がモデル化される信号を表している場合、対応する iddata
オブジェクトは data = iddata(x,[],T);
として作成できます。ここで T
は x
のサンプル時間です。
n4sid
、ssest
、ar
、arx
などの標準的な同定ツールを使用して、"出力のみ" のデータの特性を推定できる場合があります。これらのモデルは、そのスペクトル推定能力と、過去の値の測定から信号の将来値を予測する能力を評価されます。
データの解析
まず電流信号のデータを変圧器から読み込むことから、このケース スタディを始めます。
load current_data.mat
次に、電流データ (i4r
) を iddata
オブジェクトにパッケージ化します。サンプル時間は 0.001 秒 (1 ms
) です。
i4r = iddata(i4r,[],0.001) % Second argument empty for no input
i4r = Time domain data set with 601 samples. Sample time: 0.001 seconds Outputs Unit (if specified) y1
では、このデータを解析してみましょう。まず、データの概要を見ます。
plot(i4r)
データのクローズアップ表示を以下に示します。
plot(i4r(201:250))
次に、信号の生のピリオドグラムを計算します。
ge = etfe(i4r) spectrum(ge)
ge = IDFRD model. Contains spectrum of signal in the "SpectrumData" property. Output channels: 'y1' Status: Estimated using ETFE on time domain data "i4r".
このピリオドグラムはいくつかの高調波を示していますが、あまり滑らかではありません。次のコマンドを実行して平滑化されたピリオドグラムを取得します。
ges = etfe(i4r,size(i4r,1)/4); spectrum(ge,ges); legend({'ge (no smoothing)','ges (with smoothing)'})
線形周波数スケールと Hz 単位を使用するようにプロットを構成します。
h = spectrumplot(ges); opt = getoptions(h); opt.FreqScale = 'linear'; opt.FreqUnits = 'Hz'; setoptions(h,opt); axis([0 500,-5 40]) grid on, legend('ges')
明らかに、50 Hz の周波数成分とその高調波が支配的です。
spa
を使用してデータのスペクトル解析を実行しましょう。(生のピリオドグラムを計算するだけの etfe
とは違って) spa はハン ウィンドウを使ってスペクトル振幅を計算します。標準的な推定は以下のとおりです (共振スペクトルに合わせて調整されていない既定のウィンドウ % サイズを使用した場合)。
gs = spa(i4r); hold on spectrumplot(gs); legend({'ges (using etfe)','gs (using spa)'}) hold off
信号のわずかな共振をすべて見るには、非常に大きなラグ ウィンドウが必要になることがわかります。標準的なスペクトル解析はあまり役に立ちません。パラメトリック自己回帰手法によって提供されるもののような、より高度なモデルが必要です。
電流信号のパラメトリック モデリング
ここで、パラメトリック AR 法によってスペクトルを計算してみましょう。2 次、4 次、および 8 次のモデルが次のようにして得られます。
t2 = ar(i4r,2); t4 = ar(i4r,4); t8 = ar(i4r,8);
これらのスペクトルを見てみましょう。
spectrumplot(t2,t4,t8,ges,opt); axis([0 500,-8 40]) legend({'t2 (2nd order AR)','t4 (4th order AR)','t8 (8th order AR)','ges (using spa)'});
パラメトリック スペクトルでは高調波を拾えないことがわかります。AR モデルがモデル化しにくい高い周波数に注目しすぎていることが原因です。(Ljung (1999) の例 8.5 を参照してください)。
高調波を検出するためには高次モデルを試さなければなりません。
有用な次数はいくつでしょうか。arxstruc
を使用すれば、それを突き止めることができます。
V = arxstruc(i4r(1:301),i4r(302:601),(1:30)'); % Checking all order up to 30
次のコマンドを実行して、インタラクティブに最良の次数を選択します。nn = selstruc(V,'log');
上の図が示すように、n=20
の場合に劇的に下がっています。したがって、以降の説明ではこの次数を採用することにします。
t20 = ar(i4r,20); spectrumplot(ges,t20,opt); axis([0 500 -25 80]) legend({'ges (using spa)','t20 (20th order AR)'});
これですべての高調波が検出されるようになりましたが、なぜレベルが落ちたのでしょうか。それは、t20
が非常に細く高いピークを含んでいるためです。t20
の中の粗い周波数点のグリッドでは、ピークの本当のレベルはわかりません。これは次のように表すことができます。
g20c = idfrd(t20,(551:650)/600*150*2*pi); % A frequency region around 150 Hz spectrumplot(ges,t20,g20c,opt) axis([0 500 -25 80]) legend({'ges (using spa)','t20 (20th order AR)','g20c (resp. around 150 Hz)'});
このプロットを見ると、モデル t20
が非常に正確であることがわかります。細かい周波数グリッド上にプロットすることで、信号の高調波が非常に正確に再現されます。
低次の高調波のみのモデル化
主な興味の対象が低次の高調波で、低次のモデルを使用する場合は、データのプレフィルタリングを適用しなければなりません。ここではカットオフ周波数が 155 Hz の 5 次バタワース フィルターを選択します。(これで 50、100、および 150 Hz モードがカバーされます)。
i4rf = idfilt(i4r,5,155/500); % 500 Hz is the Nyquist frequency
t8f = ar(i4rf,8);
ここで、フィルタリングされたデータ (8 次モデル) から取得したスペクトルをフィルタリングされていないデータ (8 次) のスペクトルおよびピリオドグラムと比較してみましょう。
spectrumplot(t8f,t8,ges,opt) axis([0 350 -60 80]) legend({'t8f (8th order AR, filtered data)',... 't8 (8th order AR, unfiltered data)','ges (using spa)'});
フィルタリングされたデータでは、スペクトル内の最初の 3 つのピークがはっきりとわかります。
共振の数値は以下のように計算できます。サンプリングされた周波数 om
の正弦波の根は、T
をサンプル時間として、単位円上 exp(i*om*T)
の位置にあります。したがって、以下のようにモデル化を進めます。
a = t8f.a % The AR-polynomial omT = angle(roots(a))' freqs = omT/0.001/2/pi'; % show only the positive frequencies for clarity: freqs1 = freqs(freqs>0) % In Hz
a = Columns 1 through 7 1.0000 -5.0312 12.7031 -20.6934 23.7632 -19.6987 11.5651 Columns 8 through 9 -4.4222 0.8619 omT = Columns 1 through 7 1.3591 -1.3591 0.9620 -0.9620 0.3146 -0.3146 0.6314 Column 8 -0.6314 freqs1 = 216.3063 153.1036 50.0665 100.4967
このように、最初の 3 つの高調波 (50、100、および 150 Hz) を確実に見つけられます。
さらに、モデル t8f
がたとえば 100 ms (100 ステップ) 前に信号をどの程度正しく予測できるかをテストし、サンプル 201 ~ 500 に対する適合を評価することもできます。
compare(i4rf,t8f,100,compareOptions('Samples',201:500));
見てわかるとおり、最初の 3 つの高調波のモデルは未来の出力値を 100 ステップ先でも非常に正確に予測します。
高次の高調波のみのモデル化
4 次と 5 次の高調波 (200 および 250 Hz 前後) だけに関心がある場合は、データをその帯域でフィルター処理してこの高い周波数範囲に限定します。
i4rff = idfilt(i4r,5,[185 275]/500); t8fhigh = ar(i4rff,8); spectrumplot(ges,t8fhigh,opt) axis([0 500 -60 40]) legend({'ges (using spa)','t8fhigh (8th order AR, filtered to high freq. range)'});
その結果、4 次および 5 次の高調波を的確に説明するモデル t8fhigh
が得られます。このように、適切なプレフィルタリングを適用することで、希望する周波数範囲内の信号を的確に説明する低次パラメトリック モデルを構築できます。
まとめ
どのモデルが最良でしょうか。一般に、高次のモデルの方が忠実度が高くなります。これを解析するために、20 次モデルの高調波を推定する能力について考えてみましょう。
a = t20.a % The AR-polynomial omT = angle(roots(a))' freqs = omT/0.001/2/pi'; % show only the positive frequencies for clarity: freqs1 = freqs(freqs>0) %In Hz
a = Columns 1 through 7 1.0000 0.0034 0.0132 0.0012 0.0252 0.0059 0.0095 Columns 8 through 14 0.0038 0.0166 0.0026 0.0197 -0.0013 0.0143 0.0145 Columns 15 through 21 0.0021 0.0241 -0.0119 0.0150 0.0246 -0.0221 -0.9663 omT = Columns 1 through 7 0 0.3146 -0.3146 0.6290 -0.6290 0.9425 -0.9425 Columns 8 through 14 1.2559 -1.2559 1.5726 -1.5726 1.8879 -1.8879 2.2027 Columns 15 through 20 -2.2027 2.5136 -2.5136 3.1416 2.8240 -2.8240 freqs1 = Columns 1 through 7 50.0639 100.1139 149.9964 199.8891 250.2858 300.4738 350.5739 Columns 8 through 10 400.0586 500.0000 449.4611
このモデルは高調波を非常によく検出しています。このモデルは 100 ステップ先を以下のように予測します。
compare(i4r,t20,100,compareOptions('Samples',201:500));
t8f
の適合が 80% であるのに対して t20
の適合は 93% です。
したがって、高調波のキャプチャの点でも予測能力の点でも、信号の完全なモデルとして t20
を選択するのが妥当です。ただし、特定の周波数範囲を対象とするモデルの場合は低次のモデルが良い結果を出していますが、データのプレフィルターが必要です。
その他の情報
System Identification Toolbox を使った動的システムの同定の詳細については、System Identification Toolbox 製品の情報ページを参照してください。