時変システムダイナミクスを追跡するための ARX によるオンラインパラメーター推定

ライブスクリプトを開く

この例では、時変 ARX モデルのオンラインパラメーター推定を MATLAB® コマンドラインで実行する方法を説明します。モデルのパラメーターは、各タイムステップで新たな入力データによって更新されます。このモデルは線形プラントの時変ダイナミクスを再現しています。

プラント

プラントは次にように表すことができます。

$y (s) = G (s) u (s) + e (s)$

ここで、G は伝達関数、e はホワイトノイズです。プラントには 2 つの動作モードがあります。最初の動作モードでは伝達関数が次のようになります。

$G 1 (s) = \frac{4500}{(s + 5) (s^{2} + 1.2 s + 900)}$

G1(s) における減衰率の小さな極は、減衰が 0.02 で固有振動数が 30 rad/s です。2 番目の動作モードでは、これらの極の固有振動数が 60 rad/s になります。このモードでは、伝達関数は次のようになります。

$G 2 (s) = \frac{18000}{(s + 5) (s^{2} + 2.4 s + 3600)}$

プラントは t=10 秒まで最初のモードで動作した後、2 番目のモードに切り替わります。

以下は G1 と G2 のボード線図です。

wn = 30; % natural frequency of the lightly damped poles
ksi = 0.02; % damping of the poles
G1 = tf(1,conv([1/5 1],[1/wn^2 2*ksi/wn 1])); % plant in mode 1
wn = 60; % natural frequency in the second operating mode
G2 = tf(1,conv([1/5 1],[1/wn^2 2*ksi/wn 1])); % plant in mode 2
bode(G1,G2,{1,125});
legend('G1','G2','Location','Best');

Figure contains 2 axes objects. Axes object 1 with ylabel Magnitude (dB) contains 2 objects of type line. These objects represent G1, G2. Axes object 2 with ylabel Phase (deg) contains 2 objects of type line. These objects represent G1, G2.

ARX のオンラインパラメーター推定

ここではプラントの動作中にそのダイナミクスを推定するのが目的です。ARX はこの目的でよく使われるモデル構造です。ARX モデルは次の形式をもっています。

$y (t) = \frac{B (q)}{A (q)} u (t) + \frac{1}{A (q)} e (t)$

ここで、 $q^{- 1}$ は時間シフト演算子です。多項式の比率 B(q)/A(q) が (u(t) から y(t) への) 入出力モデルを表し、1/A(q) が (e(t) から y(t) への) ノイズモデルを表します。ここでは多項式 A(q) と B(q) の係数を推定します。e(t) はホワイトノイズです。

ARX モデルの構造は、線形モデルを推定する場合の最初の候補として適しています。関連する推定法は計算負荷が小さく数値的にロバストで、凸性を有しています。凸性によって、パラメーター推定が局所的最適値に留まるリスクが確実に排除されます。ただし、ARX モデルの構造ではノイズモデルに適した柔軟性が提供されません。

ノイズのモデル化に柔軟性が欠けていると、プラントの構造が ARX モデルの構造と一致しない場合やノイズがホワイトノイズではない場合に問題が生ずることがあります。この問題を解決するには 2 つの方法があります。

データのフィルター処理: アプリケーションにとってノイズモデルが重要でない場合、データのフィルター処理の手法を使用できます。詳細については、「データのフィルター処理」の節を参照してください。
異なるモデル構造: より大きな柔軟性をモデル構造にもたせるには、ARMAX、出力誤差、Box-Jenkins の各多項式モデルを使用します。

サンプル時間の選択

サンプル時間の選択は、離散時間モデルによって連続時間プラントの良好な近似を行ううえで重要です。一般的には、サンプリング周波数としてシステムの支配的なダイナミクスの 20 倍を選択します。プラントの最速の支配的ダイナミクスは 60 rad/s、すなわち約 10 Hz です。したがって、サンプリング周波数は 200 Hz になります。

Ts = 0.005; % [s], Sample time, Ts=1/200

システムの励起

適切な推定を行うには、プラントの入力がそのダイナミクスを持続的に励起しなければなりません。通常の場合、単一ステップ入力のような単純な入力では十分でありません。この例のプラントは振幅が 10、周期が 1 秒のパルスによって駆動されています。パルスの幅はその周期の 50% です。

プラントの入力信号と出力信号を次のように生成します。

t = 0:Ts:20; % Time vector
u = double(rem(t/1,1)-0.5 < 0); % pulse
y = zeros(size(u));
% Store random number generator's states for reproducible results.
sRNG = rng;
rng('default');
% Simulate the mode-switching plant with a zero-order hold.
G1d = c2d(G1,Ts,'zoh');
B1 = G1d.num{1}.'; 
A1 = G1d.den{1}.'; % B1 and A1 corresponds to G1.
G2d = c2d(G2,Ts,'zoh');
B2 = G2d.num{1}.';
A2 = G2d.den{1}.'; % B2 and A2 corresponds to G2.
idx = numel(B2):-1:1;
for ct=(1+numel(B2)):numel(t)
    idx = idx + 1;
    if t(ct)<10 % switch mode after t=10s
        y(ct) = u(idx)*B1-y(idx(2:end))*A1(2:end);
    else
        y(ct) = u(idx)*B2-y(idx(2:end))*A2(2:end);
    end
end
% Add measurement noise
y = y + 0.02*randn(size(y));
% Restore the random number generator's states.
rng(sRNG);

入出力データをプロットします。

figure();
subplot(2,1,1);
plot(t,u);
ylabel('Input (u)');
subplot(2,1,2);
plot(t,y);
ylim([-0.2 1.2])
ylabel('Output (y)');
xlabel('Time [s]');

Figure contains 2 axes objects. Axes object 1 with ylabel Input (u) contains an object of type line. Axes object 2 with xlabel Time [s], ylabel Output (y) contains an object of type line.

データのフィルター処理

プラントは次の形式をもちます。

$y (t) = G (q) u (t) + e (t)$

ここで、e(t) はホワイトノイズです。一方、ARX モデルの形式は次のようになります。

$y (t) = \frac{B (q)}{A (q)} u (t) + \frac{1}{A (q)} e (t)$

推定器は B(q) と A(q) を使用して G(q) を近似します。ただし、ノイズモデルの相違には注意してください。プラントのホワイトノイズ e(t) は y(t) に直接影響しますが、ARX モデルでは、1/A(q) でフィルター処理されたホワイトノイズ項が y(t) に影響すると仮定しています。この不一致は推定結果に悪影響を及ぼします。

ノイズモデルに関心がない場合、この不一致の影響を少なくする 1 つの方法はデータフィルターを使用することです。フィルター $F (q)$ を u(t) と y(t) の両方に適用して、 $u_{f} (t) = F (q) u (t)$ および $y_{f} (t) = F (q) y (t)$ を得ます。その後、プラントの入力 u(t) と出力 y(t) の代わりに、フィルター処理された信号 $u_{f} (t)$ および $y_{f} (t)$ を推定器で使用します。データフィルターを選択することで、e(t) が推定に及ぼす影響を小さくできます。

データフィルター F(q) は通常、アプリケーションにおける重要な周波数範囲および e(t) の特性に基づいたローパスフィルターまたはバンドパスフィルターです。ここでは、カットオフ周波数が 10 Hz の 4 次バタワースローパスフィルターが使用されています。これは、プラントにおける最速の支配的ダイナミクスの周波数 (60 rad/s) とほぼ同じです。ノイズ項に低周波数成分がないため、ここではローパスフィルターで十分です。

% Filter coefficients
Fa = [1 -3.1806 3.8612 -2.1122 0.4383]; % denominator
Fb = [4.1660e-04 1.6664e-03 2.4996e-03 1.6664e-03 4.1660e-04]; % numerator
% Filter the plant input for estimation
uf = filter(Fb,Fa,u);
% Filter the plant output
yf = filter(Fb,Fa,y);

推定の設定

オンラインパラメーター推定には recursiveARX コマンドを使用します。このコマンドは、ARX モデル構造のオンラインパラメーター推定を行うための System object™ を作成します。オブジェクトに次のプロパティを指定します。

モデル次数: [3 1 0]。プラントは 3 つの極をもつため na = 3 です。プラントには入力遅延がないため nk = 0 です。nb = 1 はシステムに零点がない状態に対応します。nb は、3 つの零点に対応する nb=4 から始めて数回の反復後に設定されたため、適切なモデルとなっています。推定するパラメーターの数はより少ない方が望ましく、nb=1 で十分な結果が得られます。
EstimationMethod: 'ForgettingFactor' (既定の設定)。この方法には ForgettingFactor というスカラーパラメーターが 1 つだけあり、このパラメーターには、パラメーター値についての限られた事前情報が必要です。
ForgettingFactor: 0.995。パラメーターが経時的に変化するため忘却係数 $λ$ は 1 より小さくなります。 $\frac{1}{1 - λ} = 200$ は、推定値に最も大きく影響する過去のデータサンプルの数です。

X = recursiveARX([3 1 0]); % [na nb nk]
X.ForgettingFactor = 0.995;

推定結果を保存する配列を作成します。これらはアルゴリズムの検証に役立ちます。

np = size(X.InitialParameterCovariance,1);
PHat = zeros(numel(u),np,np);
A = zeros(numel(u),numel(X.InitialA));
B = zeros(numel(u),numel(X.InitialB));
yHat = zeros(1,numel(u));

step コマンドを使用して、各タイムステップで 1 組の入出力データを使いパラメーター値を更新します。これは推定器のオンライン動作を示しています。

for ct=1:numel(t)
    % Use the filtered output and input signals in the estimator
    [A(ct,:),B(ct,:),yHat(ct)] = step(X,yf(ct),uf(ct));
    PHat(ct,:,:) = X.ParameterCovariance;
end

推定された伝達関数のボード線図を表示します。

G1Hat = idpoly(A(1000,:),B(1000,:),1,1,1,[],Ts); % Model snapshot at t=10s
G2Hat = idpoly(X); % Snapshot of the latest model, at t=20s
G2Hat.Ts = G1d.Ts; % Set the sample time of the snapshot
figure();
bode(G1,G1Hat);
xlim([0.5 120]);
legend('G1','Identified model at t=10s','Location','Best');

Figure contains 2 axes objects. Axes object 1 with ylabel Magnitude (dB) contains 2 objects of type line. These objects represent G1, Identified model at t=10s. Axes object 2 with ylabel Phase (deg) contains 2 objects of type line. These objects represent G1, Identified model at t=10s.

figure();
bode(G2,G2Hat);
xlim([0.5 120]);
legend('G2','Identified model at t=20s','Location','Best');

Figure contains 2 axes objects. Axes object 1 with ylabel Magnitude (dB) contains 2 objects of type line. These objects represent G2, Identified model at t=20s. Axes object 2 with ylabel Phase (deg) contains 2 objects of type line. These objects represent G2, Identified model at t=20s.

推定モデルの検証

次の手法を使用してパラメーター推定を検証します。

出力の推定値 yhat(t) の確認: step メソッドの 3 番目の出力引数は、出力 yhat(t) の 1 ステップ先の予測値です。これは現在のモデルパラメーターと、現在および過去の入出力の測定値を基にしています。y(t) と yhat(t) の間の相対誤差と絶対誤差が、適合度の目安となります。
パラメーター共分散の推定値 Phat(t) の確認: これは ForgettingFactor 推定法と KalmanFilter 推定法で使用できます。推定器の ParameterCovarianceMatrix プロパティに保存されています。Phat(t) の対角要素は、パラメーターの分散の推定値です。これは有界でなければならず、値が小さいほど優れています。
推定された時変モデルのシミュレーション: u(t) および推定されたパラメーターを使用してモデルをシミュレートし、シミュレートの出力 ysim(t) を取得します。そのうえで、y(t) と ysim(t) を比較します。ysim(t) はプラントの出力測定値なしで生成されるため、これは y(t) と yhat(t) の比較よりも厳密な検証となります。

絶対誤差 yf(t)-yhat(t) と相対誤差 (yf(t)-yhat(t))/yf(t) は次のようになります。

figure();
subplot(2,1,1);
plot(t,yf-yHat);
ylabel('Abs. Error');
subplot(2,1,2);
plot(t,(yf-yHat)./yf);
ylim([-0.05 0.05]);
ylabel('Rel. Error');
xlabel('Time [s]');

Figure contains 2 axes objects. Axes object 1 with ylabel Abs. Error contains an object of type line. Axes object 2 with xlabel Time [s], ylabel Rel. Error contains an object of type line.

絶対誤差は 1e-3 程度であり、これは測定された出力信号そのものに比べ小さな値です。下方にある相対誤差のプロットはこれを裏付けており、シミュレーションの始めを除いて誤差は 5% 未満となっています。

パラメーター共分散行列の対角要素は残差 y(t)-yhat(t) の分散でスケーリングされ、パラメーター推定値の分散を表現しています。分散の平方根はパラメーター推定値の標準偏差です。最初の 3 つの対角要素は、多項式 A(q) で推定された 3 つのパラメーターです。最後の要素は多項式 B(q) の単一のパラメーターです。A(q) の最初の推定されたパラメーターを見てみましょう。

noiseVariance = var(yf-yHat);
X.A(2) % The first estimated parameter. X.A(1) is fixed to 1

ans = -2.8635

sqrt(X.ParameterCovariance(1,1)*noiseVariance)

ans = 0.0175

標準偏差の 0.0175 は、パラメーター値の絶対値 2.86 と比較して小さな値です。これは、推定されたパラメーターの信頼度が良好であることを示します。

figure();
plot(t,sqrt(PHat(:,1,1)*noiseVariance));
ylabel('Standard deviation estimate for the parameter A(2)')
xlabel('Time [s]');

Figure contains an axes object. The axes object with xlabel Time [s], ylabel Standard deviation estimate for the parameter A(2) contains an object of type line.

推定全体にわたって不確かさは小さく、有界です。ただし、パラメーターの標準偏差も推定値である点に注意してください。こうした値は残差 y(t)-yhat(t) がホワイトであるという仮定に基づいています。これは、推定法、その関連パラメーター、推定されるモデルの構造および入力信号 u によって決まります。モデルの仮定された構造と実際の構造との差異や、持続的な入力励起の不足、あるいは推定法の非現実的な設定などが原因で、不確かさの推定値があまりに楽観的あるいは悲観的になる可能性があります。

最後に、推定されたパラメーターの保存された履歴を使用して、推定した ARX モデルをシミュレートします。このシミュレーションは、オンライン操作の間に検証用の推定ループと同時に実行することもできます。

ysim = zeros(size(y));
idx = numel(B2):-1:1;
for ct=(1+numel(B2)):numel(t)
    idx = idx + 1;
    ysim(ct) = u(idx(1))*B(idx(1),:)-ysim(idx(2:end))*A(ct,2:end)';
end
figure();
subplot(2,1,1);
plot(t,y,t,ysim);
ylabel('System Output');
legend('Measured','Estimated','Location','Best');
subplot(2,1,2);
plot(t,y-ysim);
ylim([-0.5 0.5]);
ylabel('Error, y(t)-ysim(t)');
xlabel('Time [s]');

Figure contains 2 axes objects. Axes object 1 with ylabel System Output contains 2 objects of type line. These objects represent Measured, Estimated. Axes object 2 with xlabel Time [s], ylabel Error, y(t)-ysim(t) contains an object of type line.

最初は誤差が大きくなりますが、最初の動作モードでは t=5 秒付近のより小さな値に落ち着きます。初期の大きな誤差を抑えるには、モデルパラメーターの初期推定と初期パラメーター共分散を推定器に指定します。プラントが 2 番目のモードに切り替わると、最初は誤差が大きくなりますが、その後やはり徐々に小さくなります。これにより、与えられた入力信号に対するモデル動作が、推定したモデルパラメーターによって適切に捉えられるという信頼が得られます。