このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
非断熱連続攪拌タンク反応器: Simulink でのシミュレーションによる MATLAB ファイル モデリング
この例では、IDNLGREY モデルを Simulink® 内に指定してシミュレーションを実行する方法を説明します。ここでは、モデル化の基礎として化学反応システムを使用します。この例の最初のモデル化と同定の部分は、Simulink を使用せずに実行できます。
非断熱連続攪拌タンク反応器のモデル化
連続攪拌タンク反応器 (CSTR) は、プロセス業界では比較的一般的に利用される化学システムです。ここでは、被覆非断熱 (非断熱的) タンク反応器を調査します。この詳細については、Bequette の著書『Process Dynamics:Modeling, Analysis and Simulation』(Prentice-Hall, 1998) に記述されています。容器は完全に混合されるものと仮定され、単一の 1 次不可逆発熱反応 A -> B が発生します。容器と周囲の冷却被覆の概略をプロット ウィンドウに示します。これはスケッチであり、実際には冷却液は通常は反応器の被覆全体の周囲を流れ、底部の周囲のみではありません。
図 1: CSTR の概略図
CSTR のモデルは、より詳細な制御アプローチに必要です。反応物 A の入口供給流は一定の割合 F で流入します。攪拌後、最終製品が反応物 A と同じ割合で容器からタンクに流れます (反応器タンクの容積 V は一定に維持されます)。制御戦略では、被覆温度 u_3(t) を操作して、入口供給流濃度と温度 (それぞれ入力 u_1(t) と u_2(t)) から生じる外乱によらず、反応物 A の濃度 y_1(t) を目的のレベルに維持することが要求されます。タンクの温度 y_2(t) は反応器の動作中に大きく変化するため、このプロセス変数を合理的な範囲内に維持しておくことも推奨されます。
CSTR のモデル化
CSTR システムは基本的な計算およびエネルギー保存の法則を使用してモデル化されます。時間単位あたりの容器内の反応物 A の濃度の変化 d C_A(t)/dt (= d y_1(t)/dt) は、次のようにモデル化されます。
d C_A(t)
-------- = F/V*(C_Af(t)-C_A(t)) - r(t)
dtここで、最初の項は反応物 A の入口供給流と容器内との濃度の違いによる濃度の変化を示します。2 番目の項は、容器内の化学反応によって生じる濃度の変化 (反応率) を示します。単位容積あたりの反応率は Arrhenius の法則によって記述されます。
r(t) = k_0*exp(-E/(R*T(t)))*C_A(t)
化学反応の率は絶対温度とともに指数関数的に増加することを示しています。ここで、k_0 は不明な非熱定数で、E は活性化エネルギー、R は Boltzmann の理想気体定数、T(t) (= y_2(t)) は反応器内の温度です。
同様に、エネルギー釣り合い原理を使用して (反応器内の容積が一定と仮定)、反応器内の時間単位あたりの温度変化 d T(t)/dt は次のようにモデル化できます。
d T(t)
------ = F/V(T_f(t)-T(t)) - (H/c_p*rho)*r(t) - (U*A)/(c_p*rho*V)*(T(t)-T_j(t))
dt最初の項と 3 番目の項は、供給流温度 T_f(t) と被覆冷却液温度 T_j(t) が反応器の温度と異なるために発生する変化をそれぞれ示します。2 番目の項は、容器内の化学反応によって生じる反応器の温度への影響です。この方程式では、H は反応熱のパラメーター、c_p は熱容量項、rho は密度項、U は全体的な温度転移係数、A は熱交換面積 (冷却液/容器面積) です。
まとめると、CSTR には 3 つの入力信号があります。
u_1(t) = C_Af(t) Concentration of A in inlet feed stream [kgmol/m^3]. u_2(t) = T_f(t) Inlet feed stream temperature [K]. u_3(t) = T_j(t) Jacket coolant temperature [K].
および 2 つの出力信号があります。
y_1(t) = C_A(t) Concentration of A in reactor tank [kgmol/m^3]. y_2(t) = T(t) Reactor temperature [K].
これはまた、自然なモデル状態、y_1(t) = x_1(t) および y_2(t) = x_2(t) でもあります。
元のパラメーターのいくつかをまとめて、8 つの異なるモデル パラメーターを得ます。
F Volumetric flow rate (volume/time) [m^3/h]. Fixed.
V Volume in reactor [m^3]. Fixed.
k_0 Pre-exponential nonthermal factor [1/h]. Free.
E Activation energy [kcal/kgmol]. Free.
R Boltzmann's gas constant [kcal/(kgmol*K)]. Fixed.
H Heat of reaction [kcal/kgmol]. Fixed.
HD = c_p*rho Heat capacity times density [kcal/(m^3*K)]. Free.
HA = U*A Overall heat transfer coefficient times tank
area [kcal/(K*h)] Free.ここでのパラメーター 4 つは自由として指定されています (推定対象)。実際には、非熱頻度因子 k_0 と活性化エネルギー E はベンチ スケール実験から決定することもできます。これによって次の 2 つの不明パラメーターのみを考慮することになり、同定タスクが簡略化されます。熱容量 c_p と全体的な熱転移係数 U (rho と A は既知のため、HD と HA から推定されます)。
上記で使用した表記により、次の状態空間表現が CSTR から得られます。
d x_1(t)
-------- = F/V*(u_1(t)-x_1(t)) - k_0*exp(-E/(R*x_2(t)))*x_1(t)
dt d x_2(t)
-------- = F/V(u_2(t)-x_2(t)) - (H/HD)*k_0*exp(-E/(R*x_2(t)))*x_1(t)
dt - (HA/(HD*V))*(x_2(t)-u_3(t)) y_1(t) = x_1(t)
y_2(t) = x_2(t)非断熱連続攪拌タンク反応器 IDNLGREY オブジェクトの作成
この情報は、次の内容をもつ cstr_m.m という MATLAB® ファイルに入力されます。
function [dx, y] = cstr_m(t, x, u, F, V, k_0, E, R, H, HD, HA, varargin) %CSTR_M A non-adiabatic Continuous Stirred Tank Reactor (CSTR).
% Output equations.
y = [x(1); ... % Concentration of substance A in the reactor.
x(2) ... % Reactor temperature.
]; % State equations.
dx = [F/V*(u(1)-x(1))-k_0*exp(-E/(R*x(2)))*x(1); ...
F/V*(u(2)-x(2))-(H/HD)*k_0*exp(-E/(R*x(2)))*x(1)-(HA/(HD*V))*(x(2)-u(3)) ...
];モデル化の状況を反映する IDNLGREY オブジェクトが次に作成されます。
FileName = 'cstr_m'; % File describing the model structure. Order = [2 3 2]; % Model orders [ny nu nx]. Parameters = [1; 1; 35e6; 11850; ... % Initial parameters. 1.98589; -5960; 480; 145]; InitialStates = [8.5695; 311.267]; % Initial value of the initial states. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, 'Name', ... 'Stirred tank reactor', ... 'TimeUnit', 'hours');
CSTR モデル構造の入力、状態、出力が、メソッド SET と SETINIT を使用して指定されます。また、既定の設定で初期状態を推定することを指定します。
nlgr.InputName = {'Concentration of A in inlet feed stream' ... % u(1).
'Inlet feed stream temperature' ... % u(2).
'Jacket coolant temperature'}; % u(3).
nlgr.InputUnit = {'kgmol/m^3' 'K' 'K'};
nlgr = setinit(nlgr, 'Name', {'Concentration of A in reactor tank' ... % x(1).
'Reactor temperature'}); ... % x(2).
nlgr = setinit(nlgr, 'Unit', {'kgmol/m^3' 'K'});
nlgr = setinit(nlgr, 'Fixed', {false false});
nlgr.OutputName = {'A Concentration' ... % y(1); Concentration of A in reactor tank
'Reactor temp.'}; % y(2).
nglr.OutputUnit = {'kgmol/m^3' 'K'};
CSTR モデル構造のパラメーターが定義され、F、V、R、H が固定するように設定されます。
nlgr = setpar(nlgr, 'Name', {'Volumetric flow rate (volume/time)' ... % F. 'Volume in reactor' ... % V. 'Pre-exponential nonthermal factor' ... % k_0. 'Activation energy' ... % E. 'Boltzmann''s ideal gas constant' ... % R. 'Heat of reaction' ... % H. 'Heat capacity times density' ... % HD. 'Overall heat transfer coefficient times tank area'}); ... % HA. nlgr = setpar(nlgr, 'Unit', {'m^3/h' 'm^3' '1/h' 'kcal/kgmol' 'kcal/(kgmol*K)' ... 'kcal/kgmol' 'kcal/(m^3*K)' 'kcal/(K*h)'}); nlgr.Parameters(1).Fixed = true; % Fix F. nlgr.Parameters(2).Fixed = true; % Fix V. nlgr.Parameters(5).Fixed = true; % Fix R. nlgr.Parameters(6).Fixed = true; % Fix H.
物理的な理由によって、反応熱パラメーター (反応は発熱を伴うため、常に負) が正であることがわかっています。この知識 (多少粗削りですが) も CSTR モデル構造に取り入れます。
nlgr.Parameters(1).Minimum = 0; % F. nlgr.Parameters(2).Minimum = 0; % V. nlgr.Parameters(3).Minimum = 0; % k_0. nlgr.Parameters(4).Minimum = 0; % E. nlgr.Parameters(5).Minimum = 0; % R. nlgr.Parameters(6).Maximum = 0; % H. nlgr.Parameters(7).Minimum = 0; % HD. nlgr.Parameters(8).Minimum = 0; % HA.
入力した CSTR モデル構造の概要が次に、PRESENT コマンドから得られます。
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'cstr_m' (MATLAB file):
dx/dt = F(t, u(t), x(t), p1, ..., p8)
y(t) = H(t, u(t), x(t), p1, ..., p8) + e(t)
with 3 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 8).
Inputs:
u(1) Concentration of A in inlet feed stream(t) [kgmol/m^3]
u(2) Inlet feed stream temperature(t) [K]
u(3) Jacket coolant temperature(t) [K]
States: Initial value
x(1) Concentration of A in reactor tank(t) [kgmol/m^3] xinit@exp1 8.5695 (estimated) in [-Inf, Inf]
x(2) Reactor temperature(t) [K] xinit@exp1 311.267 (estimated) in [-Inf, Inf]
Outputs:
y(1) A Concentration(t)
y(2) Reactor temp.(t)
Parameters: Value
p1 Volumetric flow rate (volume/time) [m^3/h] 1 (fixed) in [0, Inf]
p2 Volume in reactor [m^3] 1 (fixed) in [0, Inf]
p3 Pre-exponential nonthermal factor [1/h] 3.5e+07 (estimated) in [0, Inf]
p4 Activation energy [kcal/kgmol] 11850 (estimated) in [0, Inf]
p5 Boltzmann's ideal gas constant [kcal/(kgmo..] 1.98589 (fixed) in [0, Inf]
p6 Heat of reaction [kcal/kgmol] -5960 (fixed) in [-Inf, 0]
p7 Heat capacity times density [kcal/(m^3*K)] 480 (estimated) in [0, Inf]
p8 Overall heat transfer coefficient times tank area [kcal/(K*h)] 145 (estimated) in [0, Inf]
Name: Stirred tank reactor
Status:
Created by direct construction or transformation. Not estimated.
More information in model's "Report" property.
<a href="matlab:if exist('nlgr','var'), if isa(nlgr,'idnlgrey'), disp(' '); disp(get(nlgr)), else disp(' '); disp('Unable to display properties for variable nlgr because it is no longer a/an idnlgrey object.');, end, else matlab.graphics.internal.getForDisplay('nlgr'), end">Model Properties</a>
入出力データ
多くの非線形システムのシステム同定実験の設計は、線形システムよりも通常はより複雑です。CSTR にもこれは当てはまります。制御可能な入力 u_3 がシステムを十分に励起するようにする一方で、「プラントに優しく」(化学プロセスを安定に維持し、テスト期間をできるだけ短縮して生産になるべく影響しないようにするなど) なるように選択しなければなりません。[1] の論文では、CSTR への入力信号の選択について議論しています。多重正弦入力 u_3 は多重疑似ランダム入力信号よりもいくつかの理由で優れていることが説明されています。次の同定実験では、上記の論文の著者が MATLAB 入力データ生成ツール (GUI) によって生成して提供した入力信号を 2 つ使用します。1 つは推定用、もう 1 つは検証用です。
この CSTR データを読み込んで 2 つの異なる IDDATA オブジェクト、推定用の ze と検証用の zv に取り込みます。
load cstrdata Ts = 0.1; % 10 samples per hour! ze = iddata(y1, u1, 0.1, 'Name', 'Estimation data'); ze.InputName = nlgr.InputName; ze.InputUnit = nlgr.InputUnit; ze.OutputName = nlgr.OutputName; ze.OutputUnit = nlgr.OutputUnit; ze.Tstart = 0; ze.TimeUnit = 'hour'; ze.ExperimentName = 'Estimation'; zv = iddata(y2, u2, 0.1, 'Name', 'Validation data'); zv.InputName = nlgr.InputName; zv.InputUnit = nlgr.InputUnit; zv.OutputName = nlgr.OutputName; zv.OutputUnit = nlgr.OutputUnit; zv.Tstart = 0; zv.TimeUnit = 'hour'; zv.ExperimentName = 'Validation';
推定データセット ze の入力と出力は、2 つのプロットに表示されます。
figure('Name', [ze.Name ': input data']); for i = 1:ze.Nu subplot(ze.Nu, 1, i); plot(ze.SamplingInstants, ze.InputData(:,i)); title(['Input #' num2str(i) ': ' ze.InputName{i}]); xlabel(''); axis tight; end xlabel([ze.Domain ' (' ze.TimeUnit ')']);
図 2: CSTR への推定データセットの入力
figure('Name', [ze.Name ': output data']); for i = 1:ze.Ny subplot(ze.Ny, 1, i); plot(ze.SamplingInstants, ze.OutputData(:,i)); title(['Output #' num2str(i) ': ' ze.OutputName{i}]); xlabel(''); axis tight; end xlabel([ze.Domain ' (' ze.TimeUnit ')']);
図 3: CSTR からの推定データセットの出力
同定実験に進む前に、生成された入力信号によって CSTR の出力から多数の反応器の非線形が示されることに注意してください (全体として約 80 ℃の温度変化と、反応器の「点火」現象の発生が明らかです)。これによって反応器が励起されますが (通常は同定の観点からは良好)、エンジニアが実世界の反応器を動作させる方法とは、特にこれほど発熱しない反応器について、おそらく異なります。ガイドライン (Braun et al.(2002) 参照) を使用し、実験を実際に行う前に再設計できます。この場合、実験の期間を短縮して、多重正弦入力信号を短いサイクル長で使用してみるのも興味深いでしょう。目的は、制御可能な入力信号の低周波数成分を除去して、反応器の出力での変動を削減することです。
初期 CSTR モデルの性能
初期 CSTR モデルはどの程度良好でしょうか。初期モデルを入力信号 ze と zv を使用してシミュレートして調査し、計算された出力を ze と zv に含まれる実出力 (その他のパラメーターを使用してノイズを付加して、上記の IDNLGREY モデルをシミュレートして得られた出力) と比較します。
clf compare(ze, nlgr); figure; compare(zv, nlgr);
図 4: 初期 CSTR モデルの実際の出力とシミュレーションの出力の比較
パラメーター推定
初期 CSTR モデルの実出力とシミュレーションの出力は、適度に一致しています。さらに改善するため、4 つの自由モデル パラメーターとモデルの初期状態ベクトルを、推定データセット ze を使用して推定します。反復情報を表示し、最大で 25 回の検索反復を実行するよう NLGREYEST に指示します。
opt = nlgreyestOptions('Display', 'on'); opt.SearchOptions.MaxIterations = 25; nlgr = nlgreyest(ze, nlgr, opt);
必要に応じて NLGREYEST (または PEM) をもう一度呼び出し、いつでも検索を続行できます。ここでは、反復情報を一切表示せず、最大で 5 回の反復のみを実行するよう NLGREYEST は指示を与えられます。
opt.Display = 'off';
opt.SearchOptions.MaxIterations = 5;
nlgr = nlgreyest(ze, nlgr, opt);
推定された CSTR モデルの性能
推定されたモデルの性能を評価するため、もう一度 COMPARE を使用します。
compare(ze, nlgr); figure; compare(zv, nlgr);
図 5: 推定された CSTR モデルの実際の出力とシミュレーションの出力の比較
目で見てすぐに、ze と zv の両方で、推定されたモデルの出力は実出力に近いことがわかります。検証データセットについては特に大きく改善され、モデルの適合は 2 つのモデル出力に対して負の値からそれぞれ約 70 % と約 99 % まで向上しています。
推定された CSTR モデルの詳細情報は、PRESENT から返されます。
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'cstr_m' (MATLAB file):
dx/dt = F(t, u(t), x(t), p1, ..., p8)
y(t) = H(t, u(t), x(t), p1, ..., p8) + e(t)
with 3 input(s), 2 state(s), 2 output(s), and 4 free parameter(s) (out of 8).
Inputs:
u(1) Concentration of A in inlet feed stream(t) [kgmol/m^3]
u(2) Inlet feed stream temperature(t) [K]
u(3) Jacket coolant temperature(t) [K]
States: Initial value
x(1) Concentration of A in reactor tank(t) [kgmol/m^3] xinit@exp1 8.62914 (estimated) in [-Inf, Inf]
x(2) Reactor temperature(t) [K] xinit@exp1 311.215 (estimated) in [-Inf, Inf]
Outputs:
y(1) A Concentration(t)
y(2) Reactor temp.(t)
Parameters: Value Standard Deviation
p1 Volumetric flow rate (volume/time) [m^3/h] 1 0 (fixed) in [0, Inf]
p2 Volume in reactor [m^3] 1 0 (fixed) in [0, Inf]
p3 Pre-exponential nonthermal factor [1/h] 3.55889e+07 17548.3 (estimated) in [0, Inf]
p4 Activation energy [kcal/kgmol] 11853.9 0.0703052 (estimated) in [0, Inf]
p5 Boltzmann's ideal gas constant [kcal/(kgmo..] 1.98589 0 (fixed) in [0, Inf]
p6 Heat of reaction [kcal/kgmol] -5960 0 (fixed) in [-Inf, 0]
p7 Heat capacity times density [kcal/(m^3*K)] 500.71 0.194139 (estimated) in [0, Inf]
p8 Overall heat transfer coefficient times tank area [kcal/(K*h)] 150.127 0.0697455 (estimated) in [0, Inf]
Name: Stirred tank reactor
Status:
Termination condition: Maximum number of iterations or number of function evaluations reached..
Number of iterations: 6, Number of function evaluations: 7
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Estimation data".
Fit to estimation data: [71.36;99.18]%
FPE: 0.02736, MSE: 0.6684
More information in model's "Report" property.
<a href="matlab:if exist('nlgr','var'), if isa(nlgr,'idnlgrey'), disp(' '); disp(get(nlgr)), else disp(' '); disp('Unable to display properties for variable nlgr because it is no longer a/an idnlgrey object.');, end, else matlab.graphics.internal.getForDisplay('nlgr'), end">Model Properties</a>
Simulink での推定された CSTR モデルのシミュレーション
IDNLGREY モデルをインポートして Simulink 内で使用できます。Simulink モデル "cstr_sim" は検証データセット zv をインポートして、データを Simulink IDNLGREY モデル ブロックに渡します。これをシミュレートすると、MATLAB ワークスペースに保存されている使用された入力信号とともに、出力を IDDATA オブジェクト zsim に格納します (以下のコードの最後の 5 行を使用して、適切なモデル入力が Simulink モデルに渡るようにします。これは、zv と nlgr へのアクセスが保証できない、idnlgreydemo9 が関数内で実行される場合にのみ必要です)。
open_system('cstr_sim'); if ~evalin('base', 'exist(''zv'', ''var'')') cstrws = get_param(bdroot, 'modelworkspace'); cstrws.assignin('zv', zv); cstrws.assignin('nlgr', nlgr); end
図 6: 推定された CSTR モデルを含む Simulink モデル
汎用 IDNLGREY Simulink ライブラリ ブロックが標準システム同定 Simulink ライブラリ内にあり、コピーして Simulink モデル内で使用できます。たとえば、CSTR の場合、閉ループ制御配置に使用できます。
シミュレートする前に、IDNLGREY ブロックを設定しなければなりません。IDNLGREY モデル (ここでは nlgr) を保持する MATLAB ワークスペース変数を入力するか、または idnlgrey コンストラクターを "IDNLGREY model" というパラメーター編集ボックスに使用して適切な IDNLGREY モデル オブジェクトを定義して実行されます。ここで、初期状態ベクトルを指定して使用することもできます (既定の設定では指定された IDNLGREY オブジェクトの内部保存された初期状態ベクトル)。
IDNLGREY モデル オブジェクトには、SIM、PREDICT、NLGREYEST などで使用される微分または差分方程式ソルバーの設定を指定するプロパティが保存されます (モデルのシミュレーションを制御するオプションについては nlgr.SimulationOptions を参照)。Simulink では、これらの設定は常にオーバーライドされ、Simulink が指定するソルバーのオプションが使用されます。IDNLGREY モデル オブジェクトが別のソルバーを指定して使用する場合、Simulink で得られるシミュレーション結果は、MATLAB の IDNLGREY/SIM で得られる結果とは異なる場合があります。これを示す例は、例 "idnlgreydemo10" で紹介されています。
ソルバー オプションを設定して、次に cstr_sim Simulink モデルのコマンド プロンプト シミュレーションを実行します (推定された CSTR モデルについて実施されます)(idnlgreydemo9 が関数内で実行されている場合は、zsim の取得に evalin の呼び出しが必要です)。
sim('cstr_sim'); if ~exist('zsim', 'var') zsim = evalin('base', 'zsim'); end zsim.InputName = nlgr.InputName; zsim.InputUnit = nlgr.InputUnit; zsim.OutputName = nlgr.OutputName; zsim.OutputUnit = nlgr.OutputUnit; zsim.TimeUnit = 'hour';
最後に、Simulink から得られたシミュレーション結果をプロットして例を終わります。
figure('Name', 'Simulink simulation result of estimated CSTR model'); for i = 1:zsim.Ny subplot(zsim.Ny, 1, i); plot(zsim.SamplingInstants, ze.OutputData(:,i)); title(['Output #' num2str(i) ': ' zsim.OutputName{i}]); xlabel(''); axis tight; end xlabel([zsim.Domain ' (' zsim.TimeUnit ')']);
図 7: Simulink で推定された CSTR モデルのシミュレーションを実行して得られた出力
まとめ
このチュートリアルでは、非断熱連続攪拌タンク反応器のモデル化と同定を説明しました。特に、IDNLGREY モデルを Simulink 内でインポートして使用する方法を示しました。
参考文献
[1] Braun, M.W., R. Ortiz-Mojica and D.E. Rivera, "Application of minimum crest factor multisinusoidal signals for 'plant-friendly' identification of nonlinear process systems," in Control Engineering Practice, Vol. 10, No. 3, 2002, pp. 301-313.
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)