オンラインの ARMAX 多項式モデル推定

この例では次を使用します。

スクリプトを開く

この例では、オンラインの多項式モデル推定器を実装する方法を説明します。1 つの非線形化学反応プロセスについて、ARMAX モデルを 2 つ推定します。これらのモデルは、2 つの操作条件におけるプロセスの動作を再現します。このモデルの動作はオンラインで同定され、システムの操作時に適応 PI コントローラーのゲインを調整するために使用されます。

連続攪拌タンク反応器

連続攪拌タンク反応器 (CSTR) は、プロセス業界で一般的な化学システムです。CSTR システムの概略図は次のとおりです。

これは、被覆非断熱 (非断熱的) タンク反応器であり、この詳細については、Bequette の著書『Process Dynamics:Modeling, Analysis and Simulation』(Prentice-Hall, 1998) に記述されています。容器は完全に混合されるものと仮定され、単一の 1 次不可逆発熱反応 A -> B が発生します。反応物 A の入口供給流は一定流量で流入します。攪拌後、最終製品が反応物 A と同じ流量で容器から流出します (反応タンクの容積は一定)。この例で使用する CSTR の操作とその 2 状態非線形モデルの詳細については、非断熱連続攪拌タンク反応器: Simulink でのシミュレーションによる MATLAB ファイル モデリングの例で説明されています。

CSTR モデルの入力は次のとおりです。

$$ \begin{array} {ll}
u_1(t) = C_{Af}(t) \; & \textnormal{Concentration of A in inlet feed
stream} [kgmol/m^3] \\
u_2(t) = T_f(t) \; & \textnormal{Inlet feed stream temperature} [K] \\
u_3(t) = T_j(t) \; & \textnormal{Jacket coolant temperature} [K] \\
\end{array} $$

出力 (y(t)) (モデルの状態 (x(t) でもある) は次のとおりです。

$$ \begin{array} {ll}
y_1(t) = x_1(t) = C_{A}(t) \; & \textnormal{Concentration of A in reactor tank} [kgmol/m^3] \\
y_2(t) = x_2(t) = T(t) \; & \textnormal{Reactor temperature} [K] \\
\end{array} $$

制御の目標は、反応物 A の濃度 $C_{A}(t)$ を、時間とともに変化する目的のレベル $C_{Aref}(t)$ に維持することです。入口供給流の温度 $T_f(t)$ による外乱を排除するために、被覆温度 $T_j(t)$ が PI コントローラーにより操作されます。PI コントローラーの入力は、追従誤差信号 $C_{Aref}(t)-C_A(t)$ です。入口供給流の濃度 $C_{Af}(t)$ は一定であると仮定されます。Simulink® モデル iddemo_cstr は、CSTR プラントをブロック CSTR として実装します。

open_system('iddemo_cstr');

入力の供給流温度 $T_f(t)$ は、定数オフセットにホワイト ノイズの外乱を加えたものです。ノイズ パワーは 0.0075 [K^2] です。このノイズ レベルにより、目的の $C_{Aref}(t)$ から最大 2% の偏差が生じます。

この例の信号 $C_{Aref}(t)$ は、$t=400$ の時点で 1.5 [kgmol/m^3] から 2 [kgmol/m^3] にステップ変化します。このステップ変化に加え、t が [0,200) と [400,600) の範囲にある場合、$C_{Aref}(t)$ にはホワイト ノイズの摂動も含まれます。このホワイト ノイズ信号のパワーは 0.015 です。信号対ノイズ比が約 10 になるように、ノイズ パワーが経験的に調整されます。閉ループ同定で参照信号が十分に励起されない場合、一意のモデルを同定するために十分な情報が得られないことがあります。$C_{Aref}(t)$ の実装は iddemo_cstr/CA Reference ブロックにあります。

適応制御のオンライン推定

非線形モデルから、CSTR 出力 $C_{A}(t)$ は $C_{A}(t)$ のレベルが高いと制御入力 $T_j(t)$ の影響を受けやすくなることがわかっています。この感度の変化を検出するために、Recursive Polynomial Model Estimator ブロックを使用します。この情報は、$C_{A}(t)$ の変化に応じて PI コントローラーのゲインを調整するために使用されます。これは、不安定の要因となり得る高ゲイン制御ループを使用しないためです。

Recursive Polynomial Model Estimator ブロックを使用して、$T_j(t)$ から $C_{A}(t)$ への離散伝達関数をオンラインで推定します。適応制御アルゴリズムは、この伝達関数の DC ゲインを使用します。追従誤差 $C_{Aref}(t)-C_A(t)$ が、推定した伝達関数の正規化された DC ゲインで除算されます。この正則化は、PI コントローラーの設計用途である、初期操作点における追従誤差でゲイン 1 を得ることにより行われます。たとえば、DC ゲインが元の値の 2 倍である場合は、誤差信号を 2 で除算します。これは、PI コントローラーのゲインを 2 で除算することに対応します。この適応コントローラーは iddemo_cstr/Adaptive PI Controller で実装されています。

Recursive Polynomial Model Estimator ブロックの入力

'Recursive Polynomial Model Estimator' ブロックは、Simulink の System Identification Toolbox/Estimators ライブラリにあります。このブロックを使用して、ARMAX 構造をもつ線形モデルを推定します。ARMAX モデルの形式は次のとおりです。

$$
\begin{array} {l}A(q)\bar{y}(t)=B(q)\bar{u}(t)+C(q)\bar{e}(t) \\[0.1in]
A(q) = 1+a_1z^{-1}+a_2z^{-2}+\cdots+a_{na}z^{-na} \\
B(q) = (b_01+b_1z^{-1}+b_2z^{-2}+\cdots+a_{nb-1}z^{-nb+1})z^{-nk} \\
C(q) = 1+c_1z^{-1}+c_2z^{-2}+\cdots+c_{nc}z^{-nc} \\
\end{array}
$$

  • 再帰多項式モデル推定器ブロックの Inputs および Output 入力端子は、それぞれ $\bar{u}(t)$ と $\bar{y}(t)$ に対応します。CSTR モデルの $\bar{y}$ と $\bar{u}$ は、被覆温度および A の濃度平衡化操作点の偏差です。つまり、$\bar{y}=C_A(t)-\bar{C}_A(t)$、$\bar{u}=T_j(t)-\bar{T}_j(t)$ となります。ピーク振幅が 1 になるよう $\bar{u}$ と $\bar{y}$ をスケーリングして、推定問題の数値条件を改善することが推奨されます。平衡化操作点 $\bar{C}_A(t)$ および $\bar{T}_j(t)$ は、システムを操作するまで正確にはわかりません。これらの値は、1 次移動平均フィルターを使用して、測定信号から推定されて抽出されます。これらの前処理フィルターは、iddemo_cstr/Preprocess Tj ブロックと iddemo_cstr/Preprocess CA ブロックに実装されています。

open_system('iddemo_cstr/Preprocess Tj');

  • Recursive Polynomial Model Estimator ブロックのオプションの Enable 入力端子は、ブロックのパラメーター推定を制御します。パラメーター推定は、Enable 信号がゼロの場合は無効になります。Enable 信号の他のすべての値に対してパラメーター推定は有効です。この例では、$t\in[200,400)$ と $t\in[600,800)$ の時間区間で推定が無効になります。これらの区間において、測定された入力 $T_j(t)$ には閉ループ システム同定に十分な励起が含まれていません。

Recursive Polynomial Model Estimator ブロックの設定

2 次 ARMAX モデルを推定するようにブロック パラメーターを設定します。[モデル パラメーター] タブで、次のように指定します。

  • Model Structure: ARMAX。システムに作用する外乱の現在の値と過去の値 $T_f(t)$ が CSTR システムの出力 $C_A(t)$ に影響することが予想されるため、[ARMAX] を選択します。

  • Initial Estimate: None。既定では、すべての推定パラメーターについて値 0 が使用されます。

  • A(q) (na) 内のパラメーター数: 2。非線形モデルは 2 つの状態をとります。

  • B(q) (nb) 内のパラメーター数: 2。

  • C(q) (nc) 内のパラメーター数: 2。na、nb、nc の最大値はすべて 2 であるため、推定モデルは 2 次モデルに対応します。

  • 入力遅延 (nk): 1。多くの物理システムと同様に、CSTR システムには直達がありません。また、I/O 間に余分な時間遅延がありません。

  • パラメーター共分散行列: 1e4。初期推定値は非常に不確かであるため、共分散の値は高く指定してください。

  • サンプル時間: 0.1。CSTR モデルは約 0.25 Hz の帯域幅をもつことがわかっています。$T_s=0.1$ を選択するのは、1/0.1 がこの帯域幅の 20 倍 (5 Hz) より大きくなるからです。

[Algorithm and Block Options] をクリックして、次の推定オプションを設定します。

  • 推定法: Forgetting Factor。

  • 忘却係数: 1-5e-3。推定パラメーターが操作点とともに変化することが予測されるため、忘却係数を 1 未満の値に設定します。$\lambda = 1-5e-3$ を選択します。これは、記憶の時定数 $T_0=\frac{Ts}{1-\lambda}=100$ 秒に対応します。記憶の時間を 100 秒に設定することにより、各操作点で 200 秒の同定期間から大量のデータが取得され、同定に使用されます。

  • [Output estimation error] チェック ボックスをオンにします。このブロック出力を使用して推定を検証します。

  • [Add enable port] チェック ボックスをオンにします。基準端子に外部ノイズが加えられた場合にのみ、推定モデルのパラメーターを適応させます。外部ノイズが加えられなくなった場合、この端子経由でのパラメーター推定 n が無効になります。

  • 外部リセット: None。

Recursive Polynomial Model Estimator ブロックの出力

各タイム ステップで、再帰的多項式モデル推定器は $A(q)$、$B(q)$、$C(q)$ の推定および推定誤差 $\bar{e}$ を出力します。この多項式モデル推定器ブロックの Error 出力端子は $\bar{e}(t)$ を含み、これは 1 ステップ先の予測誤差とも呼ばれます。ブロックの Parameters 出力端子のバス信号は、多項式 A(q)、B(q)、C(q) の係数を含みます。選択した多項式の次数 ($na=2$、$nb=2$、$nc=2$、$nk=1$) を指定すると、Parameters バス要素には次が含められます。

$$ \begin{array} {l}
A(t) = [1 \; a_1(t) \; a_2(t)] \\
B(t) = [0 \; b_0(t) \; b_1(t)] \\
C(t) = [1 \; c_1(t) \; c_2(t)] \\
\end{array} $$

多項式 A(q)、B(q)、C(q) の推定パラメーターは、シミュレーション中に次のように変化します。

sim('iddemo_cstr');
close_system('iddemo_cstr/Preprocess Tj');
open_system('iddemo_cstr/ABC');

初期のパラメーター共分散行列に高い値が選択されているため、パラメーターの推定値は初期値 0 から即座に変化します。多項式 $A(q)$ および $B(q)$ のパラメーターは、$t=200$ での値に急速に近づきます。しかし、多項式 $C(q)$ のパラメーターにはある程度の変動が見られます。このような変動の理由の 1 つに、CSTR 出力 $C_A(t)$ に加える外乱 $T_f(t)$ が ARMAX 構造では十分にモデル化されないことが挙げられます。ここで考察している制御問題では、誤差モデル $C(q)$ は重要ではありません。これは、$T_j(t)$ から $C_A(t)$ への関係が伝達関数 $\frac{B(q)}{A(q)}$ で捕捉されているからです。したがって、この同定と制御の問題では、$C(q)$ における変動を気にする必要はありません。

$t\in[200,400)$ の間は推定器ブロックが無効になっていたため (Enable 入力端子に 0 の信号)、パラメーター推定値はこの期間一定に保たれます。パラメーター推定は、CSTR タンクが新規の操作点に切り替わる $t=400$ の時点で有効になります。$A(q)$ および $B(q)$ のパラメーターは、$t=600$ までに新しい値に収束し、以降は Enable 端子が 0 に設定されて一定に保持されます。この操作点での $A(q)$ と $B(q)$ の収束はゆっくりしたものになります。この遅い収束は、t=400 におけるパラメーター共分散行列の固有値が t=0 で設定された初期の 1e4 値と比較して小さいことが原因です。推定における信頼度の尺度であるパラメーター共分散は、タイム ステップごとに更新されます。t=0 における推定の信頼度が低かったため、アルゴリズムはパラメーター推定値をすばやく変更しました。改善されたパラメーター推定値によりシステムの動作が的確に把握され、その結果、パラメーター共分散行列の 1 ステップ先の予測誤差と固定値が小さくなります (信頼性が向上します)。システムの動作は t=400 で変化します。ただし、推定の信頼性が向上するため、ブロックでのパラメーター推定値の変更は遅くなります。Recursive Polynomial Model Estimator ブロックの [外部リセット] オプションを使用して、t=400 におけるパラメーター共分散の新しい値を指定できます。パラメーター共分散の値を表示するには、Recursive Polynomial Model Estimator ブロックの [Output parameter covariance matrix] チェック ボックスをオンにします。

推定モデルの検証

Recursive Polynomial Model Estimator ブロックの Error 出力は、推定モデルの 1 ステップ先の誤差です。

close_system('iddemo_cstr/ABC');
open_system('iddemo_cstr/1-Step Ahead Prediction Error');

システム同定の $T_j(t)$ チャネルに追加の摂動を与えない場合、1 ステップ先の誤差はより大きい値になります。誤差が大きくなる理由として、推定器ブロックが依存する入力チャネル $T_j(t)$ に十分な情報がないことが考えられます。ただし、この大きな誤差も $C_A(t)$ で測定した変動と比較すると制限された小さなものとなっています。これにより、推定したパラメーター値の信頼度が得られます。

推定モデルのより厳密な検査では、推定モデルのシミュレーションを実行し、実際のモデル出力と比較します。iddemo_cstr/Time-Varying ARMAX ブロックは、Online Polynomial Model Estimator ブロックで推定した時変 ARMAX モデルを実装しています。CSTR システムの出力と推定時変 ARMAX モデルの出力の誤差は次のとおりです。

close_system('iddemo_cstr/1-Step Ahead Prediction Error');
open_system('iddemo_cstr/Simulation Error');

ここでも、シミュレーション誤差は、$C_A(t)$ での変動と比較すると制限された小さなものとなっています。これはさらに、推定線形モデルが非線形 CSTR モデルの動作を予測できるという信頼度を与えます。

同定したモデルをさらに MATLAB で解析できます。操作点 $C_A=1.5 [kgmol/m^3]$ および $C_A=2 [kgmol/m^3]$ におけるモデル推定値はそれぞれ、$t = 200$ および $t = 600$ の時点で推定した多項式 A(q)、B(q)、C(q) を調べることにより取得できます。これらのモデルのボード線図を示します。

Ts = 0.1;
tidx = find(t>=200,1);
P200 = idpoly(AHat(:,:,tidx),BHat(:,:,tidx),CHat(:,:,tidx),1,1,[],Ts);
tidx = find(t>=600,1);
P600 = idpoly(AHat(:,:,tidx),BHat(:,:,tidx),CHat(:,:,tidx),1,1,[],Ts);
bodemag(P200,'b',P600,'r--',{10^-1,20});
legend('Estimated Model at C_A=1.5 [kgmol/m^3]', ...
       'Estimated Model at C_A=2.0 [kgmol/m^3]', ...
       'Location', 'Best');

推定モデルは濃度が高くなるほどゲインが高くなります。これは、非線形 CSTR プラントに関する事前知識と一致しています。$C_A(t)=2 [kgmol/m^3]$ における伝達関数は、低周波数においてゲインが $6dB$ 高くなります (振幅が 2 倍)。

まとめ

2 つの ARMAX モデルを推定して、2 つの操作条件における非線形 CSTR プランドの動作を再現しました。この推定は、適応コントローラーを使用して閉ループ動作時に行いました。推定結果を検証するために、2 つの信号を確認しました。それらは、1 ステップ先の予測誤差および CSTR プラント出力と推定モデルのシミュレーション間の誤差です。これらの誤差信号は、CSTR プラント出力と比較すると、一定の範囲内にあり、小さい値でした。これにより、推定された ARMAX モデル パラメーターに信頼性が提供されました。

bdclose('iddemo_cstr');