ドキュメンテーション センター

  • 評価版
  • 製品アップデート

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

非断熱連続攪拌タンク反応器:Simulink でのシミュレーションによる MATLAB ファイル モデリング

この例では、IDNLGREY モデルを Simulink 内に指定してシミュレーションを実行する方法を説明します。ここでは、モデル化の基礎として化学反応システムを使用します。この例の最初のモデル化と同定の部分は、Simulink を使用せずに実行できます。

非断熱連続攪拌タンク反応器のモデル化

連続攪拌タンク反応器 (CSTR) は、プロセス業界では比較的一般的に利用される化学システムです。ここでは、Bequette の書籍『Process Dynamics: Modeling, Analysis and Simulation』(Prentice-Hall、1998 年刊) に広範に記述されている被覆断熱 (非断熱) タンク反応器を調査します。容器は完全に混合されるものと仮定され、1 次不可逆発熱反応 A -> B が 1 回発生します。容器と周囲の冷却被覆の概略をプロット ウィンドウに示します。これはスケッチであり、実際には冷却液は通常は反応器の被覆全体の周囲を流れ、底部の周囲のみではありません。

図 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 というファイルに入力されます。

  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 initial states.
Ts            = 0;                                 % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, 'Name', ...
                'Non-adiabatic continuous stirred tank reactor model',  ...
                'TimeUnit', 'hours');

CSTR モデル構造の入力、状態、出力が、メソッド SET と SETINIT を使用して指定されます。また、既定の設定で初期状態を推定することを指定します。

nlgr = set(nlgr, 'InputName',  {'Concentration of A in inlet feed stream'   ...   % u(1).
                         'Inlet feed stream temperature'             ...   % u(2).
                         'Jacket coolant temperature'},              ...   % u(3).
          '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 = set(nlgr, 'OutputName', {'A Concentration' ...   % y(1); Concentration of A in reactor tank
                         'Reactor temp.'}, ...   % y(2).
          '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 inputs, 2 states, 2 outputs, and 4 free parameters (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   (est) in [-Inf, Inf]
    x(2)  Reactor temperature(t) [K]                          xinit@exp1   311.267   (est) in [-Inf, Inf]
 Outputs:                                                                                                
    y(1)  A Concentration(t) [kgmol/m^3]                                                                 
    y(2)  Reactor temp.(t) [K]                                                                           
 Parameters:                                                              value                          
    p1   Volumetric flow rate (volume/time) [m^3/h]                             1   (fix) in [0, Inf]    
    p2   Volume in reactor [m^3]                                                1   (fix) in [0, Inf]    
    p3   Pre-exponential nonthermal factor [1/h]                          3.5e+07   (est) in [0, Inf]    
    p4   Activation energy [kcal/kgmol]                                     11850   (est) in [0, Inf]    
    p5   Boltzmann's ideal gas constant [kcal/(kgmo..]                    1.98589   (fix) in [0, Inf]    
    p6   Heat of reaction [kcal/kgmol]                                      -5960   (fix) in [-Inf, 0]   
    p7   Heat capacity times density [kcal/(m^3*K)]                           480   (est) in [0, Inf]    
    p8   Overall heat transfer coefficient times tank area [kcal/(K*h)]       145   (est) in [0, Inf]    
Created:       09-Aug-2013 03:08:47                                                                      
Last modified: 09-Aug-2013 03:08:47                                                                      

入出力データ

多くの非線形システムのシステム同定実験の設計は、線形システムよりも通常はより複雑です。CSTR にもこれは当てはまります。制御可能な入力 u_3 がシステムを十分に励起するようにする一方で、「プラントに優しく」(化学プロセスを安定に維持し、テスト期間をできるだけ短縮して生産になるべく影響しないようにするなど) なるように選択しなければなりません。次の記事

M.W. Braun, 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", no. 10,
2002.

では、CSTR への入力信号の選択について議論しています。多重正弦入力 u_3 は多重疑似ランダム入力信号よりもいくつかの理由で優れていることが説明されています。以下の同定実験では、上記の資料の著者が提供した MATLAB 入力データ生成ツール (GUI) によって生成された入力信号を 2 つ使用します。1 つは推定用、もう 1 つは検証用です。

この CSTR データを読み込んで 2 つの異なる IDDATA オブジェクト、推定用の ze と検証用の zv に取り込みます。

load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'cstrdata'));
Ts = 0.1;   % 10 samples per hour!
ze = iddata(y1, u1, 0.1, 'Name', 'Non-adiabatic continuous stirred tank reactor estimation data');
set(ze, 'InputName', nlgr.InputName, 'InputUnit', nlgr.InputUnit);
set(ze, 'OutputName', nlgr.OutputName, 'OutputUnit', nlgr.OutputUnit);
set(ze, 'Tstart', 0, 'TimeUnit', 'hour', 'ExperimentName', 'Estimation');
zv = iddata(y2, u2, 0.1, 'Name', 'Non-adiabatic continuous stirred tank reactor validation data');
set(zv, 'InputName', nlgr.InputName, 'InputUnit', nlgr.InputUnit);
set(zv, 'OutputName', nlgr.OutputName, 'OutputUnit', nlgr.OutputUnit);
set(zv, 'Tstart', 0, 'TimeUnit', 'hour', 'ExperimentName', 'Validation');

推定データセット ze の入力と出力は、2 つのプロットに表示されます。

figure('Name', [ze.Name ': input data']);
for i = 1:ze.Nu
   subplot(ze.Nu, 1, i);
   plot(ze(:, [], 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(:, i, []));
   title(['Output #' num2str(i) ': ' ze.OutputName{i}]);
   xlabel('');
   axis('tight');
end
xlabel([ze.Domain ' (' ze.TimeUnit ')']);

図 3: CSTR からの推定データセットの出力

同定実験に進む前に、生成された入力信号によって CSTR の出力から多数の反応器の非線形が示されることに注意してください (全体として約 80 ℃の温度変化と、反応器の「点火」現象の発生が明らかです)。これによって反応器が励起されますが (通常は同定の観点からは良好)、エンジニアが実世界の反応器を動作させる方法とは、特にこれほど発熱しない反応器について、おそらく異なります。Braun 他著 (2002 年) に記述されているガイドラインを使用して、実験を再設計してから実際に実行できます。この場合、実験の期間を短縮して、多重正弦入力信号を短いサイクル長で使用してみるのも興味深いでしょう。目的は、制御可能な入力信号の低周波数成分を除去して、反応器の出力での変動を削減することです。

初期 CSTR モデルの性能

初期 CSTR モデルはどの程度良好でしょうか。初期モデルを入力信号 ze と zv を使用してシミュレートして調査し、計算された出力を ze と zv に含まれる実出力 (その他のパラメーターを使用してノイズを付加して、上記の IDNLGREY モデルをシミュレートして得られた出力) と比較します。COMPARE を 2 つの入力引数を指定して呼び出すと、各実験のモデルの初期状態ベクトル全体が推定されることに注意してください。

compare(merge(ze, zv), nlgr);

図 4: 初期 CSTR モデルの実際の出力とシミュレーションの出力の比較

パラメーター推定

初期 CSTR モデルの実出力とシミュレーションの出力は、適度に一致しています。さらに改善するため、4 つの自由モデル パラメーターとモデルの初期状態ベクトルを、推定データセット ze を使用して推定します。PEM で、反復情報を MATLAB コマンド ウィンドウに表示させ、最大 25 個の検索反復を実行します。

nlgr = pem(ze, nlgr, 'Display', 'On', 'MaxIter', 25);

必要に応じて、検索は PEM をもう一度呼び出して継続できます。このとき、PEM は反復情報を表示せず、最大 5 個の反復のみを実施します。

nlgr = pem(ze, nlgr, 'Display', 'Off', 'MaxIter', 5);

推定された CSTR モデルの性能

推定されたモデルの性能を評価するため、もう一度 COMPARE を使用します。

compare(merge(ze, 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 inputs, 2 states, 2 outputs, and 4 free parameters (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.60701   (est) in [-Inf, Inf]             
    x(2)  Reactor temperature(t) [K]                          xinit@exp1   311.188   (est) in [-Inf, Inf]             
 Outputs:                                                                                                             
    y(1)  A Concentration(t) [kgmol/m^3]                                                                              
    y(2)  Reactor temp.(t) [K]                                                                                        
 Parameters:                                                              value         standard dev                  
    p1   Volumetric flow rate (volume/time) [m^3/h]                                 1           0   (fix) in [0, Inf] 
    p2   Volume in reactor [m^3]                                                    1           0   (fix) in [0, Inf] 
    p3   Pre-exponential nonthermal factor [1/h]                          3.56314e+07     6827.28   (est) in [0, Inf] 
    p4   Activation energy [kcal/kgmol]                                       11854.9   0.0310496   (est) in [0, Inf] 
    p5   Boltzmann's ideal gas constant [kcal/(kgmo..]                        1.98589           0   (fix) in [0, Inf] 
    p6   Heat of reaction [kcal/kgmol]                                          -5960           0   (fix) in [-Inf, 0]
    p7   Heat capacity times density [kcal/(m^3*K)]                           500.479   0.0836654   (est) in [0, Inf] 
    p8   Overall heat transfer coefficient times tank area [kcal/(K*h)]        150.16   0.0332088   (est) in [0, Inf] 
                                                                                                                      
The model was estimated from the data set 'Non-adiabatic continuous stirred tank reactor estimation data', which      
contains 241 data samples.                                                                                            
Loss function 0.0035168 and Akaike's FPE 0.00363354                                                                   
Created:       09-Aug-2013 03:08:47                                                                                   
Last modified: 09-Aug-2013 03:09:09                                                                                   
if exist('simulink','builtin')~=5
    disp('The remainder of this tutorial requires Simulink.');
    return;
end

Simulink での推定された CSTR モデルのシミュレーション

IDNLGREY モデルをインポートして Simulink 内で使用できます。Simulink モデル "cstr_sim" は検証データセット zv をインポートして、データを Simulink IDNLGREY モデル ブロックに渡します。これをシミュレートすると、MATLAB ワークスペースに保存されている使用された入力信号とともに、出力を IDDATA オブジェクト zsim に格納します(以下のコードの最後の 5 行を使用して、適切なモデル入力が Simulink モデルに渡るようにします。これは、zv と nlgr へのアクセスが保証できない、idnlgreydemo9 が関数内で実行される場合にのみ必要です)。

open('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 Simlink ライブラリ ブロックが標準システム同定 Simulink ライブラリ内にあり、コピーして Simulink モデル内で使用できます。たとえば、CSTR の場合、閉ループ制御配置に使用できます。

シミュレートする前に、IDNLGREY ブロックを設定しなければなりません。IDNLGREY モデル (ここでは nlgr) を保持する MATLAB ワークスペース変数を入力するか、または idnlgrey コンストラクターを "IDNLGREY model" というパラメーター編集ボックスに使用して適切な IDNLGREY モデル オブジェクトを定義して実行されます。ここで、初期状態ベクトルを指定して使用することもできます (既定の設定では指定された IDNLGREY オブジェクトの内部保存された初期状態ベクトル)。

IDNLGREY モデル オブジェクトには、SIM、PREDICT、PEM などが使用する微分/差分方程式ソルバーの設定を指定するプロパティが保存されます (nlgr.Algorithm.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
set(zsim, 'InputName', nlgr.InputName, 'InputUnit', nlgr.InputUnit);
set(zsim, 'OutputName', nlgr.OutputName, 'OutputUnit', nlgr.OutputUnit);
set(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(:, i, []));
   title(['Output #' num2str(i) ': ' zsim.OutputName{i}]);
   xlabel('');
   axis('tight');
end
xlabel([zsim.Domain ' (' zsim.TimeUnit ')']);

図 7: Simulink で推定された CSTR モデルのシミュレーションを実行して得られた出力

まとめ

このチュートリアルでは、非断熱連続攪拌タンク反応器のモデル化と同定を説明しました。特に、IDNLGREY モデルを Simulink 内でインポートして使用する方法を示しました。

その他の情報

System Identification Toolbox™ を使った動的システムの同定についての詳細は、System Identification Toolbox 製品の情報ページを参照してください。

この情報は役に立ちましたか?