問題の定義

伝達遅延が 2 秒から 1 秒へと変化する線形システムについて考えます。伝達遅延は、入力が測定出力に影響を与えるまでの所要時間です。この例では、オンライン推定の手法とデータ セグメント化の手法を用いて伝達遅延の変化を検出します。システムからの測定された入出力データは、データ ファイルpdmAbruptChangesData.matにあります。

データを読み込んでプロットします。

loadpdmAbruptChangesData.matz = iddata(z(:,1),z(:,2)); plot(z) gridon

伝達遅延の変化は 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 plotColors = {'r','g','b'};斧头= gca; cla(ax)fork = 3:-1:1 h(k) = animatedline('Color',Colors{k});% lines for a, b1 and b2 parametersendh(4) = animatedline('Marker','.','Color',[0 0 0]);% line for Llegend({'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]; gridonboxon%% Now perform recursive estimation and show resultsn0 = 6; L = NaN(N,nk); B_old = NaN(1,3);forct = 1:N [A,B] = step(Estimator,Output(ct),Input(ct));ifct>n0 L(ct) = norm(B-B_old); B_old = B;endaddpoints(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を使うことで、時変パラメーター推定値の平滑化に基づいた自動検出手法が扱いやすくなります。