計測自動制御学会東北支部 40 周年記念学術講演会
2004 年 12 月 22 日，23 日 日本大学工学部

講演番号：C12

可変グリッド を用いた急激な変化を伴う分布定数系の状態推定に関する研究◯ 望月宏希, 吉田雅俊, 山下善之, 松本繁 （ 東北大学）

State estimation of distributed parameter system involving steep gradients using
adaptive multiresolution grid
∗Hiroki Mochizuki, Masatoshi Yoshida, Yoshiyuki Yamashita and Shigeru Matsumoto(Tohoku University)
Abstract– An efficient state estimation algorithm is developed by using adaptive-grid method. Since
the grids are dynamically adapted so that the computation error becomes small, the limited grids allow
high accurate online state estimation with extended Kalman filter. This method is particularly useful forthe distributed parameter systems involving steep gradients, which are frequently found in many chemicalprocesses. It is shown that the estimation of hot-spot in a fixed bed reactor can be done well even if themodel error exists.
Key Words: adaptive grid, distributed parameter system, state estimation

1

緒言V4

多くの化学プロセスは、その状態量が時間のみでなく位置によっても変化する分布定数系システムである。
また、状態量分布が急激に変化するようなプロセスも数多く存在している。このようなプロセスの数値解析や制御を行う場合には、状態量の急激な変化の取り扱いが重要となる。一般的な手法では、状態量の急激な変化を補償するために計算の際には非常に細かく離散化したモデルを用いて計算する必要があり、その結果計算量が膨大となってしまう。
例えば、化学プロセスにおいて頻繁に用いられる触媒充填層反応器では反応器入口付近に局所的な温度ピークが生じることがある。これはホットスポットと呼ばれ、触媒の劣化や製品の質に大きな影響をおよぼす。そのため、高品質の製品を得るためにはホットスポットを含めた反応器全体の濃度や温度分布を正確に把握しておく必要がある。しかし反応器内部の濃度情報を測定することは困難であるため、一般的には観測可能な情報から全体の濃度や温度分布などを推定する手法が用いられている。ところが 、ホットスポットのように急激な変化が生じるプロセスを従来法を用いて推定した場合、計算量が膨大となるためにオンラインで推定することが非常に困難である。
そこで本研究では、状態量分布が急激に変化するプロセスとしてホットスポットを対象に、状態量分布に合わせてグリッドを変化させる adaptive-gird 法を利用した拡張カルマンフィルタを用いることで高精度かつ効率的に状態量分布を推定する手法を提案する。

2V3V2

$%<=k=0

k=1:;



k=2

次に示すような 2 進格子点を考える。
noV j = zkj ∈ R : zkj = 2−j k, k ∈ N


"#

230



k=5BC

k=6

j=4,-!

j=3*+

k=7*

67

k=8(

()

1

k=1

k=2

k=3

k=0



890
1V1



@A

k=445

k=0./



2
3

k=3&

k=4&'



k=1

j=2
j=1

k=2

Fig. 1: Example of points in a dyadic grid.
fluidi−2UPi−1ifDi+1i+2f=i+1/2

Fig. 2: Definition of local variables.
jj
おける状態量 yU, yPj , yDを用いて解像度レベル j +1 における点 f での状態量 yfj+1 を補間することを考える。
jj
まず最初に (1)-(3) 式を定義する。これは、ẑU= 0, ẑD=jj
1, ŷU = 0, ŷD = 1 となるように規格化された変数である。

ŷPj =ẑPj =ẑfj =adaptive-grid 法>?

jyPj − yU

jj
yD− yUjzPj − zU

jjzD− zUjzfj − zU

jjzD− zU

(1)

(2)

(3)

さらに状態量 yfj+1 の補間値 ỹfj+1 はこれらの変数を用いて次式より得られる。

ここで、j, k はそれぞれ Fig. 1 に示すような解像度レベル及び空間位置を表し 、j ∈ N である。

jjjỹfj+1 = yU+ ŷfj+1 (yD− yU)

この 2 進格子点から解像度レベル j における点のみを取り出し 、Fig. 2 に示すような離散点 U, P, D の 3 点に

ここで ŷfj+1 は規格化された状態量であり、格子間隔が一定の場合には SMART と呼ばれる手法により次式か161(4)ら計算することができる [1]。
·µ¶¸
3j+1jj 3 jŷf = max ŷP , min 3ŷP , ŷP + , 1483

反応系と物理モデル

本研究では、触媒充填層反応器を用いて o-キシレンから無水フタル酸を得る反応系を対象とした [3]。この反応系は Fig. 3 に示すように逐次反応と並行反応が起こり、生成物である無水フタル酸の他に副生成物として二酸化炭素も生じる。

(5)

この式より ŷfj+1 が得られ、これを (4) 式に代入することで補間値 ỹfj+1 を求めることができる。

次に 、解像度レベル j からの補間値 ỹfj+1 とシミュOC

レーション値 yfj+1 との差として、次式の補間誤差係数

djf を定義する。

OC

¯R1

¯¯
¯djf = ¯yfj+1 − ỹfj+1 ¯

O2

(6)

O

phthalic anhydride

CH 3

djf が予め定められた閾値 δ よりも小さい場合には点zfj+1 の状態量 y j+1 には重要な情報は含まれていない
ため取り除くことができる。反対に djf が閾値 δ よりも大きい場合には点 zfj+1 は次の計算に用いるグリッドに加えられる。このようにして得られたグリッド を用いることで、状態量の情報を非常に効率的かつ精度を落とすこと無く圧縮することができる。adaptive-grid法の手順を以下に示す [2]。

R2O2CH 3R3CO 2O2H2Oo −xylene

Fig. 3: Kinetic scheme for the oxidation of o-xylene
半径方向の温度や濃度分布、また軸方向 z の拡散は無視できるものと仮定すると、このシステムの物理モ
デルは (7)-(9) 式で表現される。但し 、o-キシレンの濃度を CA 、無水フタル酸の濃度を CB 、反応器温度を Tとした。

1. 時間 t = tstart における状態量 ykJmin から (4) 式を用いて解像度 Jmin + 1 における補間値 ỹfJmin +1 を計算する。そして、補間値 ỹfJmin +1 とシミュレーション値 yfJmin +1 との差として補間誤差係数 dJf minを計算する。

∂CA∂CA= −v+ µ(1 − ²)RA
∂t
∂z
∂CB
∂CB²= −v+ µ(1 − ²)RB∂t∂z(1 − ²)ρs Cps + ²ρf Cpf ∂T∂T= −v+
ρf Cpf∂t∂z²

2. 解像度 Jmin + 1 から Jmax まで、同様の方法により補間誤差係数 djf を計算する。

µ(1 − ²)
2Uf cR̃ +(Tc − T )ρf cpfr0 ρf cpf

3. 係数 djf を閾値 δ と比較する。djf が δ よりも大きい場合には対応する格子点 zfj+1 を次ステップに用いるグリッドに加える。
j+14. こ うし て 得られ たグ リッド に 点 z2(f+i) , i =−NL , · · · , NR を加える。ここで NL , NR ∈ N であり、これらの点は次の時間積分によって動く可能性のある鋭いピークを補償するためのものである。

5. 点 zfj+1 よりも 1 つ上の解像度 j + 2 における点
j+2j+2z4f+1 , z4f −1 をグ リッド に追加する。これらの点はその領域において状態量がさらに急峻となった場合を補償するためのものである。

(7)
(8)

(9)

ここで、o-キシレンの生成速度 RA 、無水フタル酸の生成速度 RB 、反応により発生する熱量 R̃ は次式で示される。

RA = −R1 − R3(10)RB = R1 − R 2(11)R̃ =3X(−∆Hi )Ri(12)i=1

但し 、R1 , R2 , R3 は Fig. 3 に示されるそれぞれの反応速度である。
µ¶
E1
R1 = k1 exp −CA(13)Rg Tµ¶
E2R2 = k2 exp −CB(14)
Rg Tµ¶
E3
R3 = k3 exp −CA(15)Rg T

6. 格子点の作成時の補間誤差係数の計算に必要となる全ての点を加える。つまり、djf の計算に必要とjjなる点 zU, zPj , zDを加える。
7. 最小解像度 Jmin における全ての点を加える。これ
らは計算の間、常に存在している点である。
8. 以上の操作により新しいグリッドが作成され 、このグリッド を用いて時間刻み ∆t の間だけ微分方程式の時間積分を行う。本研究においては、時間積分は陰的な解法である逐次過緩和法を用いて計算した。

このシステムは分布定数系であるため、(7)-(9) 式に現れる対流項 (∂Ci /∂z, ∂T /∂z) を風上差分を用いて離散化した状態空間モデル ẋ=f (x), y=hx を用いて状態推定を行う。但し 、x は CA , CB , T を要素とする状
態量、y は出力であり、本研究においては o-キシレンと無水フタル酸の入口濃度及び反応器内数点の温度である。

9. Step1 から Step7 までを t = tend となるまで繰り返す。

162

5

状態推定

集中定数系に近似されたシステムにカルマンフィルタを適用し 、o-キシレン及び無水フタル酸の濃度分布と温度分布を推定する。状態量の推定値 x̂(k) は、モデルによる状態量の予測値 x̃(k) を用いて次式のフィルタ方程式より算出される。但し 、モデルによる予測値x̃(k) は前述の adaptive-grid 法を用いて導出している。

x̂(k) = x̃(k) + K(k) [y − hx̃(k)](16)x̃(k) = fk (x̂(k − 1))(17)

このとき、カルマンゲイン K(k) やリカッチ方程式は次式で示される。

K(k) = P (k)hT (k) ×{R(k) + h(k)P (k)hT (k)}−1(18)

S(k) = {I − K(k)h(k)}P (k)(19)

P (k + 1) = F (k)S(k)F T (k) + Q(k)(20)

ここで 、fk (x(k)) は f (x) の離散時間関数、F (k) はfk (x(k)) を状態量の推定値 x̂(k) の周りで線形化した行列である。
このカルマンゲイン K(k) や誤差共分散行列 P (k) の計算の際には、状態量の数が大きな問題となる。一般的なカルマンフィルタの場合には状態量の数が一定であり、カルマンゲインや誤差共分散行列の次元も一定で不変である。しかし本手法ではグリッド を状態量の分布に合わせて変化させているため、カルマンゲインや誤差共分散の次元もその都度変えてゆく必要がある。
特に、誤差共分散行列はステップ毎に更新して行かなければならず、これをどのように計算していくかが問題となる。そこで本手法では、次のような方法を用いることで誤差共分散行列を更新し 、カルマンゲインを算出している。

結果と考察

Fig. 4 に Jmax =8, Jmin =5 とした際の本手法によって得られた 1000 秒後の濃度分布及び温度分布を示す。
ここで z は軸方向長さであり、z=0 が入口、z=4 が出口を示している。この計算で用いられた格子点数はおよそ 100 点であった。そこで比較として、同じ 100 点で均等に離散化した格子点を用いて同様の推定を行った。その結果を Fig. 5 に示す。これらの結果を比較すると、濃度分布についてはど ちらの手法も比較的良好に推定できていると言える。一方温度分布について検討してみると、本手法ではおよそ 50 K という大きなホットスポットも非常に良好に推定できているが、100点で均等に離散化した一般的な手法ではホットスポッ
トの急激な変化を精度良く推定することができず、およそ 20 K もの誤差が生じている。また、この 2 つの手法による計算時間はほぼ同じ程度であり、これはオンラインで推定することが可能な程度であった。
次に格子点数を 2Jmax + 1 = 257 点として均等に配置した場合について計算を行った。その結果、推定の精度は本手法の場合とほぼ同程度であったが 、本手法を用いた場合の 10 倍以上という膨大な計算時間が必要であった。
続いて様々なパラメータにモデル誤差が存在する場
0.2
0.18

o−xyleneo−xylene(estimation)phthalic anhydride
phthalic anhydride(estimation)

(a)

0.16
0.14
C [mol/m3]

4

0.12
0.1
0.08
0.06
0.04

1. 最初に、状態量の数が最大となる場合の誤差共分
散行列の初期値 Pmax (0) を与える。状態量が o-キシレン及び無水フタル酸の濃度、反応器温度の 3つであり、また最大の格子点数は 2Jmax + 1 点なので、Pmax (0) は (3(2Jmax + 1) × 3(2Jmax + 1)) の正方行列となる。

0.02
000.5

1

1.5

2z [m]

2.5700

3.5

4reacotr temp.
estimated temp.
coolant temp.
sensor(b)

2. カルマンゲインを計算する際には、Pmax (k) から必要な要素、つまり adaptive-grid 法によって用いられたグリッドに対応した要素のみを抽出することで 誤差共分散行列 P (k) を構成する。adaptivegird 法によって作成された格子点数が n 点だとすると、この誤差共分散行列 P (k) は、(3n × 3n) の行列となる。

3690680

T [K]670660
650

3. 得られた誤差共分散行列 P (k) を用いてカルマンゲイン K(k) や次ステップの誤差共分散行列 P (k +1|k) を算出する。

640
630
6200

4. P (k) の作成の際に Pmax (k) から抽出した各要素の値を P (k + 1|k) に置き換えることで誤差共分散行列 Pmax (k) を更新し 、Pmax (k + 1) を得る。

0.5

1

1.5

2z [m]

2.5

3

3.5

4Fig. 4: Estimation results by the adaptive-grid
method : (a)Concentration profiles, (b)Temperatureprofiles.

5. 次ステップ以降においても Step2 から Step4 の操作を繰り返すことで、誤差共分散行列 Pmax (k) を更新して行く。

163

0.2
0.18
0.160.18

0.14

0.12

0.12

0.1
0.08

0.06

0.04

0.04

0.02

0.02

0.5

1

1.5

2
z [m]

710
700

2.5

3

3.5

004

1

1.5

2
z [m]

2.5

3

3.5

4

reactor temp.
estimated temp.
coolant temp.
sensor

(b)
690
680

680

670
T [K]

670
660

660
650

650

640

640

630

630
620
0

0.5

700

reactor temp.
estimated temp.
coolant temp.
sensor

(b)

690

T [K]

0.1
0.08

0.06

00

o−xyleneo−xylene(estimation)phthalic anhydride
phthalic anhydride(estimation)

(a)

0.16

0.14C [mol/m3]

C [mol/m3]0.2o−xyleneo−xylene(estimation)phthalic anhydride
phthalic anhydride(estimation)

(a)

0.5

11.52z [m]2.5

33.562004

0.5

11.52z [m]

2.5

3

3.5

4

Fig. 5: Estimation results by the uniform grid (100points) : (a)Concentration profiles, (b)Temperatureprofiles.

Fig. 6: Estimation results by the adaptive-grid
method with the model errors : (a)Concentration profiles, (b)Temperature profiles.

合について、同様の検討を行った。Fig. 6 に本反応系で進行する 3 つの反応の頻度因子全てに +20% のモデル誤差が存在する場合の濃度分布、温度分布を示す。その結果、濃度分布に関してはそれほど 良好な推定結果が得られていない。これは、頻度因子自体にモデル誤差を含んでおり、なおかつ入口以外の濃度情報は全く得られないとためであると考えられる。一方、温度分布に関しては内部数点の温度情報が得られるため、このような大きなモデル誤差が存在する場合であっても本手法を用いることで良好に推定を行うことが可能であった。

Ei : 反応 i の活性化エネルギー
P (k) : 誤差共分散行列
Q(k) : システムノイズ共分散行列
R(k) : 観測ノイズ共分散行列
Rg : 気体定数
T : 反応器温度
Tc : 冷媒温度
Uf c : 熱伝達係数
ki: 反応 i の頻度因子 r0 : 反応器内径v: フィード 流速z : 軸方向長さ−∆Hi : 反応 i の反応熱
² : 空隙率µ: 触媒活性ρf : 流体密度
ρs : 触媒密度6参考文献

結言[1] P. H. Gaskell and A. K. C. Lau, “CURVATURECOMPENSATED CONVECTIVE TRANSPORT:SMART, A NEW BOUNDEDNESS PRESERVING
TRANSPORT ALGORITHM”, International Journal for Numerical Methods in Fluids, 8, 617–641,
1988.
[2] P. Cruz, M. A. Alves, F. D. Magalhães and A.
Mendes, “Solution of hyperbolic PDEs using a stable adaptive multiresolution method”, Chemical Engineering Science, 58, 1777–1792, 2003.
[3] Wei Wu and Ming-Yuan Huang, “Nonlinear inferential control for an exothermic packed-bed reactor: geometric approaches”, Chemical Engineering Science,
58, 2023–2034, 2003.

本研究では、触媒充填層反応器を用いた際に生じるホットスポットを対象として adaptive-grid 法を利用した拡張カルマンフィルタを構成した。そして本手法と均等に離散化したグリッド を用いた場合について比較検討した結果、本手法を適用することで高精度かつ短時間に状態量を推定することが可能であることが示された。

7

記号

CA : o-キシレン濃度
Cpf : 流体熱容量

CB : 無水フタル酸濃度
Cps : 触媒熱容量

164