Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

加圧された変圧器からの電流信号のモデル化

この例では、測定された信号のモデル化を示します。400 kV の 3 相変圧器が加圧されたときの R 位相からの電流信号を解析します。測定はスウェーデンの Sydkraft AB 社によって実行されました。

ここでは電流信号のモデル化における関数 ar の使用について説明します。最初に、信号のノンパラメトリック解析が実行されます。次に、妥当なモデル次数を選択するためのツールについて論じ、信号のモデル化における ar の使用について説明します。選択した高調波の範囲だけにモデルを適合させる方法についても説明します。

はじめに

信号は自己回帰線形モデルのインパルス応答と見なせるため、ar などのツールを使ってモデル化できます。

信号のデータは、オブジェクトの出力データを信号値に設定し、入力を空のままにしておくことで、iddata オブジェクトにカプセル化することができます。たとえば x(t) がモデル化される信号を表している場合、対応する iddata オブジェクトは data = iddata(x,[],T); として作成できます。ここで Tx のサンプル時間です。

n4sidssestararx などの標準的な同定ツールを使用して、"出力のみ" のデータの特性を推定できる場合があります。これらのモデルは、そのスペクトル推定能力と、過去の値の測定から信号の将来値を予測する能力を評価されます。

データの解析

まず電流信号のデータを変圧器から読み込むことから、このケース スタディを始めます。

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 製品の情報ページを参照してください。