同定手法を使用したシステムの急激な変化の検出
この例では、オンライン推定の手法と自動データ セグメント化の手法を使用して、システムの動作に見られる急激な変化を検出する方法を説明します。この例では、System Identification Toolbox™ の機能を使用し、Predictive Maintenance Toolbox™ は必要ありません。
問題の定義
伝達遅延が 2 秒から 1 秒へと変化する線形システムについて考えます。伝達遅延は、入力が測定出力に影響を与えるまでの所要時間です。この例では、オンライン推定の手法とデータ セグメント化の手法を用いて伝達遅延の変化を検出します。システムからの測定された入出力データは、データ ファイル pdmAbruptChangesData.mat
にあります。
データを読み込んでプロットします。
load pdmAbruptChangesData.mat z = iddata(z(:,1),z(:,2)); plot(z) grid on
伝達遅延の変化は 20 秒前後で起こりますが、これをプロットで確認するのは容易ではありません。
多項式 A
の係数を 1 つ、多項式 B
の係数を 2 つ、および遅延を 1 つもつ ARX 構造を使用してシステムをモデル化します。
ここで、A = [1 a]
および B = [0 b1 b2]
です。
このモデルには直達がないため、多項式 B
の最初の係数はゼロです。システム ダイナミクスの変化に合わせて、3 つの係数 a
、b1
、b2
の値は変化します。b1
がゼロに近くなると、多項式 B
で先頭にゼロが 2 つ並ぶため、実質的な伝達遅延は 2 サンプルになります。b1
がこれより大きい場合、実質的な伝達遅延は 1 サンプルになります。
したがって、伝達遅延の変化を検出するために、多項式 B
の係数における変化を監視することができます。
オンライン推定を使用した変化の検出
オンライン推定アルゴリズムは、新しいデータが入手可能になると、モデル パラメーターと状態推定値を再帰的に更新します。オンライン推定を実行するには、System Identification Toolbox ライブラリの Simulink ブロックを使用するか、コマンド ラインで recursiveARX
などの再帰同定ルーチンを使用します。オンライン推定は、機械の老化や天候の移り変わりといった時変ダイナミクスをモデル化したり、電気機械システムの故障を検出するために使用できます。
推定器がモデル パラメーターを更新するにつれ、システム ダイナミクスの変化 (遅延) が、パラメーター b1
および b2
の値における通常より大きい変化によって示されます。多項式 B
の係数の変化は、次を計算することで追跡されます。
ARX モデルのオンライン パラメーター推定には recursiveARX
オブジェクトを使用します。
na = 1; nb = 2; nk = 1; Estimator = recursiveARX([na nb nk]);
再帰推定アルゴリズムを NormalizedGradient
、適応ゲインを 0.9 として指定します。
Estimator.EstimationMethod = 'NormalizedGradient';
Estimator.AdaptationGain = .9;
iddata
オブジェクト z
から生のデータを抽出します。
Output = z.OutputData; Input = z.InputData; t = z.SamplingInstants; N = length(t);
アニメーション化されたラインを使用して、推定したパラメーター値と L
をプロットします。推定を行う前に、これらのアニメーション化されたラインを初期化します。ストリーミング データをシミュレートするには、推定器にデータを 1 サンプルずつ供給します。推定の前にモデル パラメーターを初期化し、そのうえでオンライン推定を実行します。
%% Initialize plot Colors = {'r','g','b'}; ax = gca; cla(ax) for k = 3:-1:1 h(k) = animatedline('Color',Colors{k}); % lines for a, b1 and b2 parameters end h(4) = animatedline('Marker','.','Color',[0 0 0]); % line for L legend({'a','b1','b2','Deviation'},'location','southeast') title('ARX Recursive Parameter Estimation') xlabel('Time (seconds)') ylabel('Parameter value') ax.XLim = [t(1),t(end)]; ax.YLim = [-2, 2]; grid on box on %% Now perform recursive estimation and show results n0 = 6; L = NaN(N,nk); B_old = NaN(1,3); for ct = 1:N [A,B] = step(Estimator,Output(ct),Input(ct)); if ct>n0 L(ct) = norm(B-B_old); B_old = B; end addpoints(h(1),t(ct),A(2)) addpoints(h(2),t(ct),B(2)) addpoints(h(3),t(ct),B(3)) addpoints(h(4),t(ct),L(ct)) pause(0.1) end
データの先頭にある n0
個 (6 個) のサンプルは、変化検出器 L
の計算には使用されません。この期間中は初期状態が不明であるためパラメーターの変化が大きくなります。
Signal Processing Toolbox の findpeaks
コマンドを使用して、L
のピークの位置をすべて検出します。
[v,Loc] = findpeaks(L); [~,I] = max(v); line(t(Loc(I)),L(Loc(I)),'parent',ax,'Marker','o','MarkerEdgeColor','r',... 'MarkerFaceColor','y','MarkerSize',12)
fprintf('Change in system delay detected at sample number %d.\n',Loc(I));
Change in system delay detected at sample number 21.
最大ピークの位置は多項式 B
の係数の最大の変化に対応し、したがってこれが伝達遅延の変化の位置となります。
オンライン推定の手法では推定法やモデル構造を選択するオプションがより多く提供されますが、データ セグメント化の手法は、急激な孤立した変化を自動的に検出するのに役立ちます。
データ セグメント化を使用した変化の検出
データ セグメント化アルゴリズムは、データを、動的動作の異なる領域に自動的にセグメント化します。これは、動作状態のエラーや変化に起因する急激な変化を捉える場合に役立ちます。segment
コマンドは、単出力データについてこの操作を容易にします。segment
は、システムの運用中に時変動作を捉える必要がない場合に使う、オンライン推定手法の代替方法です。
データ セグメント化の用途には、音声信号のセグメント化 (各セグメントが音素に対応)、障害検出 (セグメントが障害あり/なしの操作に対応)、システムの各種動作モードの推定などがあります。
segment
コマンドへの入力には、測定データ、モデルの次数、システムに影響するノイズの分散の推定 r2
などが含まれます。分散がまったく未知である場合、自動的に推定できます。オンライン推定で使ったものと同じ次数の ARX モデルを使用して、データのセグメント化を実行します。分散を 0.1 に設定します。
[seg,V,tvmod] = segment(z,[na nb nk],0.1);
セグメント化の手法は、AFMM (Adaptive Forgetting Through Multiple Models、複数モデルによる適応忘却) に基づいています。この手法の詳細については、Andersson, Int. J. Control Nov 1985 を参照してください。
時変システムの追跡には複数モデル法が使用されます。結果として得られる追跡モデルは複数のモデルの平均であり、segment
の 3 番目の出力引数 tvmod
として返されます。
追跡モデルのパラメーターをプロットします。
plot(tvmod) legend({'a','b_1','b_2'},'Location','best') xlabel('Samples'), ylabel('Parameter value') title('Time-varying estimates')
これらのパラメーターの軌跡が、recursiveARX
を使って推定されたものと類似している点に注意してください。
segment
は、tvmod
および q
(モデルが急激な変化を示す確率) を使用して、変化が発生した時点を判定します。これらの時点を使用して、追跡モデルに円滑化プロシージャを適用することにより、セグメント化されたモデルを構築します。
セグメント化されたモデルのパラメーター値は、segment
の最初の出力引数 seg
に返されます。後続の各行の値は、その対応する時点における、基礎となるセグメント化モデルのパラメーター値です。これらの値は、各行を通して一定しており、システム ダイナミクスが変化したと判定されたときにのみ変化します。したがって、seg
の値は区分的な定数となります。
a
、b1
、b2
の各パラメーターの推定値をプロットします。
plot(seg) title('Parameter value segments') legend({'a','b1','b2'},'Location','best') xlabel('Time (seconds)') ylabel('Parameter value')
サンプル番号 19 の周辺でパラメーター値に変化が見られます。b1
の値が、小さい (ゼロに近い) 値から大きい (1 に近い) 値に変わります。b2
の値は、その逆のパターンを示しています。B
のパラメーターの値におけるこうした変化は、伝達遅延の変化を示すものです。
segment
の 2 番目の出力引数 V
は、セグメント化モデルの損失関数 (すなわち、セグメント化モデルの推定予測誤差の分散) です。V
を使用してセグメント化モデルの品質を評価できます。
セグメント化アルゴリズムで最も重要な 2 つの入力は、r2
と、segment
の 4 番目の入力引数 q
であることに注意してください。この例では、既定値の 0.01 が適切であるため q
は指定されていません。r2
の値を小さくしたり q
の値を大きくすると、セグメント点の数が多くなります。適切な値を見つけるために r2
と q
を変化させて、最良の結果が得られる値を使用します。通常、セグメント化アルゴリズムは q
よりも r2
に対してより敏感です。
まとめ
オンライン推定の手法とデータ セグメント化の手法を使用した、システム ダイナミクスの急激な変化の検出について評価しました。オンライン推定手法では、推定プロセスにおいてより優れた柔軟性と制御性が提供されます。ただし、頻度の低い変化や急激な変化の場合には、segment
を使うことで、時変パラメーター推定値の平滑化に基づいた自動検出手法が扱いやすくなります。