Notebook

連続撹拌反応槽 (CSTR) の最適制御¶

JModelica.org Users Guide - version 2.10 6.5.2.1 を実習します。ソースコードは$JMODELICA_HOME/Python/pyjmi/examples/files/CSTR.mop にありますが、実習に必要な部分のみを切り取って解説します。本文書は Users Guid の例題を解説するものであり、モデルのソースコードやオペレーションで使用する Pythonのコードなどの著作権は Modelon AB にあります。

連続撹拌反応槽(Continuously Stirred Tank Reactor, CSTR)の概要¶

容積 $V$ の容器に、入口から濃度 $c_0$ 、温度 $T_0$ の物質が流量 $F_0$ で流れ込みます。容器内で連続的に撹拌されて化学反応した結果、出口で濃度 $c(T)$ 、温度 $T(t)$ となって流出します。化学反応の速度はアレニウスの式に従って濃度 $c(t)$ や温度 $T(t)$ に依存します。また、容器には冷却装置がついていて冷却温度 $T_c$ で冷却されています。

Derivation of Stirred Tank Reaction Optimal Control, Lorenz T.Biegler

のモデルが参考になります。

状態変数 (state variables)¶

$c(t)$ 濃度 (concentration)
$T(t)$ 温度 (temperature)

制御変数 (control variable)¶

$T_c(t)$ 冷却温度 (cooling temperature)

方程式 (equationns)¶

濃度 $c(t)$ の変化は次の方程式で表されます。右辺第１項は流れによる物質の出入り、右辺第２項は化学反応による濃度変化を表しています。

$\dot{c}(t) = \frac{F_0 (c_0 - c(t))}{V} - k_0 c(t) e^{-\frac{E_{div} R}{T(t)}}$

温度 $T(t)$ の変化は次式で表されます。右辺第１項は流れによる熱の出入り、第２項は化学反応による発熱、第3項は冷却部との熱伝達をあらわします。

$\dot{T}(t) = \frac{F_0(T_0-T(t))}{V} - \frac{dHk_0c(t)}{\rho C_p}e^{-\frac{E_{div}R}{T(t)}}+\frac{2U}{r \rho C_p} (T_c(t)-T(t))$

冷却部の熱伝達率を $U$ とし、伝熱面積を $A_c$ とすると、熱流量は $U A_c (T_c(t)-T(t)$ となり、容器内物質の熱容量は $\rho C_p V$ となります。温度変化率は、

$\frac{U A_c (T_c(t)-T(t))}{\rho C_p V}$ となり、これと上の温度の変化式の右辺第3項を比較すると

$r = 2V/A_c$ は長さの次元の量になります。

パラメータ¶

$V$ 容積 (volume), $r = 2V/A_c$ 長さの次元の量 (length)
$F_0$ 流量 (flow rate), $c_0$ 入口濃度 (inlet concentration), $T_0$ 入口温度 (inlet temperaturea)
$k_0$ 速度定数 (rate constant), $E_{div}$ 活性化エネルギー (activation energy), $dH$ 反応熱 (heat of reaction), $R$ 気体定数 (gas constant)
$\rho$ 密度 (density), $C_p$ 比熱 (specific heat), $U$ 熱伝達率 (heat transfer coefficient)

最適化制御問題 (Optimal Control Problem) の概要¶

微分代数方程式系(DAEs)の動的最適化問題 (Optimal Control Problem)では、初期化したモデル方程式の微分が全てゼロになるべきです。このような点を停留点(stationary point)といいます。

ここでは、低温で反応速度が小さい停留点 A から高温で反応速度が大きい停留点 B まで状態変数が移動し、

$\int_{t_s}^{t_e}{ q_c (c_{ref}-c)^2 + q_T (T_{ref}-T)^2 + q_{Tc} (Tc_{ref}-Tc)^2 } dt$ を最小にするように冷却温度 Tc を制御する問題を解きます。ここで、

$c_{ref}, T_{ref}, Tc_{ref}$ は、停留点 B の濃度、温度、冷却温度です。

$q_c, q_T, q_{Tc}$ は重み定数で、ここでは、

$q_c = q_T = q_{Tc} = 1$ とします。

次のような手順を行います。

CSTR モデルを作成し、動作を確認します。
CSTR_Init モデルを使って、2つの停留点 A, B を求めます。
CSTR_Init_Optimization モデルを使って、初期推定値の軌跡(Initial Guess Trajectory)を求めます。
CSTR_Opt2を使用し最適制御シミュレーションを行います。

使用するモデル¶

CSTR¶

連続撹拌反応槽のモデルです。微分代数方程式の初期値問題です。

初期条件として濃度 c_init と温度 T_init を設定します。
制御変数として冷却温度 Tc を設定します。
状態変数である濃度 c(t) と温度 T(t) が時間とともに変化します。

CSTR_Init¶

停留点(stationary point)を求めるためのモデルです。CSTRを継承して初期化方程式として

$\frac{dc}{dt} = 0, \ \frac{dT}{dt} = 0$ を加えています。冷却温度 Tc を設定してモデルを初期化すると、停留点の状態変数 c と T が得られます。

CSTR_Init_Optimization¶

direct collocation 法では、よい初期推定値(good initial guesses)を与えるとロバストに収束することがよくあります。問題が非凸(non-convex)である場合には初期値はより重要になります。初期推定値は全ての状態変数に対して最適化制御を行う時間に渡って必要となります。このモデルは、初期推定値の軌跡 (initial trajectories) を得るためのモデルです。変数としてCSTRモデルを持ち、次の方程式を含んでいます。

$\frac{d (cost)}{dt} = q_c (c_{ref}-c)^2 + q_T (T_{ref}-T)^2 + q_{Tc} (Tc_{ref}-Tc)^2$ これにより、最適化の目的関数

$cost = \int{ q_c (c_{ref}-c)^2 + q_T (T_{ref}-T)^2 + q_{Tc} (Tc_{ref}-Tc)^2 } dt$ がどのように変化するか確認できます。

このモデルは、オリジナルの例題から一部変更しました。コメントアウトしている部分

input u = T_ref
parameter Real T_ref = 320
cstr.Tc = T_ref

のようなっていましたが、u はCSTRモデルの制御変数、T_refは目的関数のパラメータとして分離して考えるようにしました。

input u
parameter Real T_ref = 320
cstr.Tc = u

Modelica で記述した部分 CSTR1.mo と、最適化の部分 CSTR2.mop （後述）を別々のファイルにしました。

In [1]:

%%writefile CSTR1.mo
package CSTR1

  model CSTR "A CSTR"

    // Parameters
    parameter Modelica.SIunits.VolumeFlowRate F0=100/1000/60 "Inflow";
    parameter Modelica.SIunits.Concentration c0=1000 "Concentration of inflow"; 
    // Modelica.Blocks.Interfaces.RealInput Tc "Cooling temperature";
    // parameter Modelica.SIunits.VolumeFlowRate F=100/1000/60 "Outflow"; 
    parameter Modelica.SIunits.Temp_K T0 = 350;
    parameter Modelica.SIunits.Length r = 0.219;
    parameter Real k0 = 7.2e10/60;
    parameter Real EdivR = 8750;
    parameter Real U = 915.6;
    parameter Real rho = 1000;
    parameter Real Cp = 0.239*1000;
    parameter Real dH = -5e4;
    parameter Modelica.SIunits.Volume V = 100 "Reactor Volume";

    // Control Variable
    Modelica.Blocks.Interfaces.RealInput Tc "Cooling temperature";    
    
    // Initial State
    parameter Modelica.SIunits.Concentration c_init = 1000;
    parameter Modelica.SIunits.Temp_K T_init = 350;
  
    // State 
    Real c(start=c_init,fixed=true,nominal=c0);
    Real T(start=T_init,fixed=true,nominal=T0);
    
  equation 
    der(c) = F0*(c0-c)/V-k0*c*exp(-EdivR/T);
    der(T) = F0*(T0-T)/V-dH/(rho*Cp)*k0*c*exp(-EdivR/T)+2*U/(r*rho*Cp)*(Tc-T);

  end CSTR;
  
  model CSTR_Init
    extends CSTR(c(fixed=false),T(fixed=false));
  initial equation
    der(c) = 0;
    der(T) = 0;
  end CSTR_Init;

  model CSTR_Init_Optimization
    CSTR cstr "CSTR component";
    Real cost(start=0,fixed=true);
    // input Real u = Tc_ref;
    input Real u;
    parameter Real c_ref = 500;
    parameter Real T_ref = 320;
    parameter Real Tc_ref = 350;
    parameter Real q_c = 1;
    parameter Real q_T = 1;
    parameter Real q_Tc = 1;

  equation
    // cstr.Tc = Tc_ref;
    cstr.Tc = u; 
    der(cost) = q_c*(c_ref-cstr.c)^2 + q_T*(T_ref-cstr.T)^2 + q_Tc*(Tc_ref-cstr.Tc)^2;
  end CSTR_Init_Optimization;

end CSTR1;

Overwriting CSTR1.mo

準備¶

必要なモジュールをインポートします。

In [2]:

%matplotlib notebook
import matplotlib.pyplot as plt
from pymodelica import compile_fmu
from pyfmi import load_fmu

CSTR モデルの動作を確認します。¶

デフォルトの状態で CSTR モデルのシミュレーションを実行します。制御変数として、冷却温度を $Tc=0$ とし、初期条件として濃度 $c_{init} = 1000$ 、温度 $T_{init} = 350$ とした場合の濃度 $c(t)$ 、温度 $T(t)$ 、冷却温度 $Tc(t)$ をプロットします。

In [3]:

fmu = compile_fmu("CSTR1.CSTR", "CSTR1.mo")
model = load_fmu("CSTR1_CSTR.fmu")
res = model.simulate(final_time=150)
t = res["time"]
c = res["c"]
T = res["T"]
Tc = res["Tc"]

plt.figure(1)
plt.subplot(311)
plt.plot(t, c)
plt.ylabel("Concentration")
plt.xlim(0,150)
plt.grid()

plt.subplot(312)
plt.plot(t,T)
plt.ylabel("Temperature")
plt.xlim(0,150)
plt.grid()

plt.subplot(313)
plt.plot(t,Tc)
plt.ylabel("Cooling Temperature")
plt.xlim(0,150)
plt.grid()

Final Run Statistics: --- 

 Number of steps                                 : 72
 Number of function evaluations                  : 92
 Number of Jacobian evaluations                  : 2
 Number of function eval. due to Jacobian eval.  : 4
 Number of error test failures                   : 0
 Number of nonlinear iterations                  : 88
 Number of nonlinear convergence failures        : 0

Solver options:

 Solver                   : CVode
 Linear multistep method  : BDF
 Nonlinear solver         : Newton
 Linear solver type       : DENSE
 Maximal order            : 5
 Tolerances (absolute)    : [ 0.001    0.00035]
 Tolerances (relative)    : 0.0001

Simulation interval    : 0.0 - 150.0 seconds.
Elapsed simulation time: 0.00426888465881 seconds.

CSTR_Init モデルを使って２つの停留点 (stationary points) を求めます。¶

まずモデルをコンパイルしてFMUを作成し、ロードします。

In [4]:

init_fmu = compile_fmu("CSTR1.CSTR_Init","CSTR1.mo")
init_model = load_fmu(init_fmu)

Stationary Point A を求める¶

冷却温度 Tc_0_A = 250 を設定して、CSTR_init モデルを初期化(initialize)します。初期化方程式(initial equations)の評価の結果、停留点濃度 c_0_A と停留点温度 T_0_A が得られます。

In [5]:

Tc_0_A = 250
init_model.set('Tc', Tc_0_A)
init_model.initialize()
[c_0_A,T_0_A] = init_model.get(['c','T'])
print('Stationary Point A')
print('Tc = %f' % Tc_0_A)
print('c = %f, T = %f' %(c_0_A, T_0_A) )

Stationary Point A
Tc = 250.000000
c = 956.271352, T = 250.051971

Stationary Point B を求める¶

同様に冷却温度 Tc_0_B = 280 を設定して、停留点濃度 c_0_B と停留点温度 T_0_B を求めます。

In [6]:

init_model.reset()
Tc_0_B = 280
init_model.set('Tc', Tc_0_B)
init_model.initialize()
[c_0_B,T_0_B] = init_model.get(['c','T'])
print('Stationary Point B')
print('Tc = %f' % Tc_0_B)
print('c = %f, T = %f' %(c_0_B, T_0_B) )

Stationary Point B
Tc = 280.000000
c = 338.775781, T = 280.099198

CSTR_Init_Optimization モデルを使って初期推定値の軌跡 (Initial Trajectories) を求めます。¶

停留点 A の状態からスタートし、冷却温度をステップ的に Tc_0_A から Tc_0_B に変えた場合の状態変数の変化を求めます。

冷却温度の制御信号を作成します。時刻0に Tc_0_A から Tc_0_B に変化するものとします。

In [7]:

def stepf(t):
  startTime = 0
  if t <= startTime:
    return Tc_0_A
  else:
    return Tc_0_B

CSTRモデルの初期状態として停留点Aの状態変数を設定し、目標関数のパラメータとして停留点Bの値を設定します。そして上で定義した制御変数の関数を設定して、シミュレーションを実行します。シミュレーション結果 init_res に含まれる初期推定値の軌跡 (initial trajectories) をプロットします。

In [8]:

# Compile the optimization initialization model
init_sim_fmu = compile_fmu("CSTR1.CSTR_Init_Optimization", "CSTR1.mo")
# Load the model
init_sim_model = load_fmu("CSTR1_CSTR_Init_Optimization.fmu")

# Set initial and reference values
init_sim_model.set('cstr.c_init', c_0_A)
init_sim_model.set('cstr.T_init', T_0_A)
init_sim_model.set('c_ref', c_0_B)
init_sim_model.set('T_ref', T_0_B)
init_sim_model.set('Tc_ref', Tc_0_B)

# Simulate with constant input Tc
init_res = init_sim_model.simulate(start_time=0., final_time=150, input=("u", stepf))

# Extract variable profiles
t_init_sim = init_res['time']
c_init_sim = init_res['cstr.c']
T_init_sim = init_res['cstr.T']
Tc_init_sim = init_res['cstr.Tc']

# Plot the initial guess trajectories
plt.figure(2)

plt.subplot(3, 1, 1)
plt.plot(t_init_sim, c_init_sim)
plt.grid()
plt.ylabel('Concentration')
plt.title('Initial guess obtained by simulation')
plt.xlim(-1,150)

plt.subplot(3, 1, 2)
plt.plot(t_init_sim, T_init_sim)
plt.grid()
plt.ylabel('Temperature')
plt.xlim(-1,150)

plt.subplot(3, 1, 3)
plt.plot(t_init_sim, Tc_init_sim)
plt.grid()
plt.ylabel('Cooling temperature')
plt.xlabel('time')
plt.xlim(-1,150)

Final Run Statistics: --- 

 Number of steps                                 : 30
 Number of function evaluations                  : 33
 Number of Jacobian evaluations                  : 1
 Number of function eval. due to Jacobian eval.  : 3
 Number of error test failures                   : 0
 Number of nonlinear iterations                  : 32
 Number of nonlinear convergence failures        : 0

Solver options:

 Solver                   : CVode
 Linear multistep method  : BDF
 Nonlinear solver         : Newton
 Linear solver type       : DENSE
 Maximal order            : 5
 Tolerances (absolute)    : [  1.00000000e-06   1.00000000e-03   3.50000000e-04]
 Tolerances (relative)    : 0.0001

Simulation interval    : 0.0 - 150.0 seconds.
Elapsed simulation time: 0.00321006774902 seconds.

Out[8]:

(-1, 150)

CSTR_Opt2 を使用して、最適化制御問題を解きます。¶

最適化モデル CSTR_Opt2 のソースコード¶

最小化する目的関数

$\int_{startTime}^{finalTime}{q_c (c_{ref}-c(t))^2 + q_T (T_{ref}-T(t))^2 + q_{Tc} (Tc_{ref}-Tc(t))^2}dt$ の被積分関数 (objectiveIntegrand)と、動的最適化の開始時間 (startTime)、終了時間 (finalTime) を設定します。被積分関数は

$10^{-4}$ を乗じてスケーリングしています。

制御変数 u と最適化するモデルのオブジェクト cstr、目的関数のパラメータのデフォルト値などを定義します。

In [9]:

%%writefile CSTR2.mop
optimization CSTR_Opt2(
  objectiveIntegrand=1e-4*(q_c*(c_ref-cstr.c)^2 + q_T*(T_ref-cstr.T)^2 + q_Tc*(Tc_ref-cstr.Tc)^2),
  startTime=0.0,
  finalTime=150)
 
  input Real u(start = 350,initialGuess=350,min=230,max=370)=cstr.Tc;
  CSTR1.CSTR cstr(c(initialGuess=300),T(initialGuess=300, max=350),Tc(initialGuess=350));

  parameter Real c_ref = 500;
  parameter Real T_ref = 320;
  parameter Real Tc_ref = 300;
  parameter Real q_c = 1;
  parameter Real q_T = 1;
  parameter Real q_Tc = 1;
end CSTR_Opt2;

Overwriting CSTR2.mop

最適化モデル CSTR_Opt2 をロードします。

In [10]:

from pyjmi import transfer_optimization_problem
op = transfer_optimization_problem("CSTR_Opt2", ["CSTR1.mo","CSTR2.mop"])

目的関数のパラメータとして停留点Bの変数を設定し、初期値として停留点Aの変数を設定します。

In [11]:

op.set('Tc_ref', Tc_0_B)
op.set('c_ref', float(c_0_B))
op.set('T_ref', float(T_0_B))
op.set('cstr.c_init', float(c_0_A))
op.set('cstr.T_init', float(T_0_A))

最適化モデルのオプションを設定し、シミュレーションを実行します。 init_trajとnominal_trajに初期推定値の軌跡(initial trajectories)を含む init_res を設定しています。

In [12]:

# Set options
opt_opts = op.optimize_options()
opt_opts['n_e'] = 19 # Number of elements
opt_opts['init_traj'] = init_res
opt_opts['nominal_traj'] = init_res
opt_opts['IPOPT_options']['tol'] = 1e-10
opt_opts['verbosity'] = 1
# Solve the optimal control problem
res_opt = op.optimize(options=opt_opts)

Initialization time: 0.11 seconds

Total time: 0.21 seconds
Pre-processing time: 0.00 seconds
Solution time: 0.09 seconds
Post-processing time: 0.01 seconds

最適化制御シミュレーションの結果をプロットします。

In [13]:

# Extract variable profiles
c_res=res_opt['cstr.c']
T_res=res_opt['cstr.T']
Tc_res=res_opt['cstr.Tc']
time_res = res_opt['time']
c_ref = res_opt['c_ref']
T_ref = res_opt['T_ref']
Tc_ref = res_opt['Tc_ref']

# Plot the results
plt.figure(3)

plt.subplot(3, 1, 1)
plt.plot(time_res, c_res)
plt.plot(time_res, c_ref, '--')
plt.grid()
plt.ylabel('Concentration')
plt.title('Optimized trajectories')
plt.xlim(0,150)

plt.subplot(3, 1, 2)
plt.plot(time_res, T_res)
plt.plot(time_res, T_ref, '--')
plt.grid()
plt.ylabel('Temperature')
plt.xlim(0,150)

plt.subplot(3, 1, 3)
plt.plot(time_res, Tc_res)
plt.plot(time_res, Tc_ref, '--')
plt.grid()
plt.ylabel('Cooling temperature')
plt.xlabel('time')
plt.xlim(0,150)
plt.show()

再び CSTR_Init_Optimization モデルを使って最適制御シミュレーションを検証する¶

最適化シミュレーションで得られた最適化した冷却温度を CSTR モデルに入力して、最適化の効果を検証します。

最適化シミュレーションの結果から入力信号を抽出して opt_input オブジェクトを生成します。

In [14]:

# Get optimized input
(_, opt_input) = res_opt.get_opt_input()

例題ではCSTRモデルを用いますが、ここでは目的関数を計算する機能をもつ CSTR_Init_Optimization モデルを使用します。モデルをロードします。

In [15]:

sim_model = load_fmu("CSTR1_CSTR_Init_Optimization.fmu")

初期値、目的関数を計算するためのパラメータ、ソルバーのオプションパラメータを設定してシミュレーションを実行します。

In [16]:

# Set initial and reference values
sim_model.set('cstr.c_init', c_0_A)
sim_model.set('cstr.T_init', T_0_A)
sim_model.set('c_ref', c_0_B)
sim_model.set('T_ref', T_0_B)
sim_model.set('Tc_ref', Tc_0_B)

sim_opts = sim_model.simulate_options()
sim_opts['CVode_options']['rtol'] = 1e-6
sim_opts['CVode_options']['atol'] = 1e-8

# Simulate with constant input Tc
res_sim = sim_model.simulate(start_time=0., final_time=150, input=("u", opt_input), options=sim_opts)

Final Run Statistics: --- 

 Number of steps                                 : 228
 Number of function evaluations                  : 343
 Number of Jacobian evaluations                  : 5
 Number of function eval. due to Jacobian eval.  : 15
 Number of error test failures                   : 32
 Number of nonlinear iterations                  : 342
 Number of nonlinear convergence failures        : 0

Solver options:

 Solver                   : CVode
 Linear multistep method  : BDF
 Nonlinear solver         : Newton
 Linear solver type       : DENSE
 Maximal order            : 5
 Tolerances (absolute)    : 1e-08
 Tolerances (relative)    : 1e-06

Simulation interval    : 0.0 - 150.0 seconds.
Elapsed simulation time: 0.0509309768677 seconds.

最適化シミュレーションの結果と検証用モデルによるシミュレーションの結果をプロットします。

In [17]:

# Extract variable profiles
c_sim=res_sim['cstr.c']
T_sim=res_sim['cstr.T']
Tc_sim=res_sim['cstr.Tc']
time_sim = res_sim['time']

# Plot the results
plt.figure(4)

plt.subplot(311)
plt.plot(time_res,c_res,'--')
plt.plot(time_sim,c_sim)
plt.legend(('optimized','simulated'))
plt.grid()
plt.ylabel('Concentration')
plt.xlim(0,150)

plt.subplot(312)
plt.plot(time_res,T_res,'--')
plt.plot(time_sim,T_sim)
plt.legend(('optimized','simulated'))
plt.grid()
plt.ylabel('Temperature')
plt.xlim(0,150)

plt.subplot(313)
plt.plot(time_res,Tc_res,'--')
plt.plot(time_sim,Tc_sim)
plt.legend(('optimized','simulated'))
plt.grid()
plt.ylabel('Cooling temperature')
plt.xlabel('time')
plt.xlim(0,150)
plt.show()

冷却温度 $Tc$ をステップ的に変化させた初期推定値の場合と最適化制御を行った場合の目的関数

$\int_{startTime}^{t}{q_c (c_{ref}-c(t))^2 + q_T (T_{ref}-T(t))^2 + q_{Tc} (Tc_{ref}-Tc(t))^2}dt$ をプロットします。

In [18]:

init_cost = init_res['cost']
sim_cost = res_sim['cost']
plt.figure(5)
plt.title("cost")
plt.plot(t_init_sim, init_cost, time_sim, sim_cost)
plt.legend(['Initial Guess', 'Optimiezed'])
plt.xlim(0,150)

Out[18]:

(0, 150)

In [ ]: