最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

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

プラント

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

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

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

G1(s)=4500(s+5)(s2+1.2s+900)

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

G2(s)=18000(s+5)(s2+2.4s+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');

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

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

y(t)=B(q)A(q)u(t)+1A(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 つの方法があります。

  1. データのフィルター処理: アプリケーションにとってノイズ モデルが重要でない場合、データのフィルター処理の手法を使用できます。詳細については、「データのフィルター処理」の節を参照してください。

  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]');

データのフィルター処理

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

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

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

y(t)=B(q)A(q)u(t)+1A(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) の両方に適用して、uf(t)=F(q)u(t) および yf(t)=F(q)y(t) を得ます。その後、プラントの入力 u(t) と出力 y(t) の代わりに、フィルター処理された信号 uf(t) および yf(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 より小さくなります。11-λ=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();
bode(G2,G2Hat);
xlim([0.5 120]);
legend('G2','Identified model at t=20s','Location','Best');

推定モデルの検証

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

  1. 出力の推定値 yhat(t) の確認: step メソッドの 3 番目の出力引数は、出力 yhat(t) の 1 ステップ先の予測値です。これは現在のモデル パラメーターと、現在および過去の入出力の測定値を基にしています。y(t) と yhat(t) の間の相対誤差と絶対誤差が、適合度の目安となります。

  2. パラメーター共分散の推定値 Phat(t) の確認: これは ForgettingFactor 推定法と KalmanFilter 推定法で使用できます。推定器の ParameterCovarianceMatrix プロパティに保存されています。Phat(t) の対角要素は、パラメーターの分散の推定値です。これは有界でなければならず、値が小さいほど優れています。

  3. 推定された時変モデルのシミュレーション: 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]');

絶対誤差は 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]');

推定全体にわたって不確かさは小さく、有界です。ただし、パラメーターの標準偏差も推定値である点に注意してください。こうした値は残差 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]');

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

まとめ

ここでは ARX モデルのオンライン パラメーター推定を行いました。このモデルはモードの切り替わるプラントのダイナミクスを再現しました。システム出力の推定値、シミュレート値、測定値の間にある誤差と、パラメーター共分散の推定値を検討することで、推定されたモデルを検証しました。