数値的最適化法: RAND 法

杉山耕一朗 (北大理・地惑) sugiyama@gfd-dennou.org

2004/05/29

\underline{数値的最適化法: RAND 法 1}

目次

第 1 章大気の平衡状態の表現 2 第 2 章 RAND 法 3

2.1 気相のみの場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 多相系の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

第 3 章 Newton 法との比較 8 \underline{数値的最適化法: RAND 法第 1 章大気の平衡状態の表現 2}

第1章大気の平衡状態の表現

熱力学量として温度・圧力・組成を選択すると, 大気の平衡状態は大気中の元素数保存の下でギブス自由エネルギー G の最小化された状態として記述することができる.

N 個の化学種が p 個の相に分かれているものとする. ギブス自由エネルギー G と化学ポテンシャル μ は物質の理想性を仮定することで以下のように書ける:

G(T, p, n) = ∑_{i=1}^{N} ∑_{ϕ=1}^{p} n_{i}^{ϕ}μ_{i}^{ϕ}(T, p, n)  (1.1)
μ_{i}^{ϕ}(T, p, n) = μ_{i}^{ϕ}^{∘}(T) + RT ln x_{i}^{ϕ} + δ_{1ϕ}RT ln p/p_{0}  (1.2)

ここで添字 i は化学種を表し, 添字 ϕ は相を表す. T は温度, p は圧力, p_{0} は基準圧力, R は気体定数, n はモル数, μ^{∘} は標準化学ポテンシャル, x_{i}^{ϕ} = n_{i}^{ϕ}/ ∑_{k} n_{k}^{ϕ} はモル比である. また相 i = 1 を気相とした. 化学ポテンシャルは相, もしくは純物質/混合物に依らず同一な形式で書けるので, 系に存在する全ての化学種の標準化学ポテンシャル μ_{i}^{ϕ}^{∘} を予め把握しておけさえすれば, 混合気体や溶液, 純物質の化学ポテンシャルは一元的に計算することができる.

平衡状態を求めるためには, 温度・圧力を与え, 元素数保存の条件

∑_{ϕ} ∑_{i} a_{ie}^{ϕ}n_{i}^{ϕ} = B_{e}  (1.3)
n_{i}^{ϕ} ≥ 0  (1.4)

の下でギブス自由エネルギーを最小化するようなモル数 n を探索すれば良い. ただし e は元素を表し, a_{ie} は化学種 i に含まれる元素 e の個数を意味する.

oboro モデルでは, (1.3) 式の条件の下で (1.1) 式の極値を求めるという「制約つき問題の最適化法」として White et al. (1958) の RAND 法を用いる. この手法は現在では「Newton 法」として分類される最適化法である.

第2章 RAND 法

RAND 法は「制約つき最適化法」の「Newton 法」に分類される最適化法である. この RAND 法は元々, 気相のみの系での平衡状態を求めるために開発され, その後に多相系を扱えるよう拡張された. この方法はギブスの自由エネルギーを適当な組成のまわりで組成に関する 2 次の方程式に近似し, 反復的に極値を求めるというものである.

本章の議論は, 温度圧力一定条件で気相のみを扱った White et al. (1958) の研究をもとにしている. この方法は RAND 法として知られており, Boynton (1960), Kubert and Stephanou (1960), Oliver et al. (1962) によって凝縮相を含む場合にも拡張されてきた.

2.1 気相のみの場合

適当な物質存在量を m_{i} とする. m_{i} よりも熱化学平衡に近い物質存在度を n_{i}(= m_{i} + δm_{i}) とする. 平衡時の熱力学ポテンシャル G を m_{i} のまわりでテーラー展開し, δm_{i} に関する 2 次の微小量まで近似式 Q(n_{i}) を作る.

Q(n_{i}) = G(m_{i}) + ∑_{i} (\frac{∂G}{∂m_{i}} ) δm_{i} + \frac{1}{2} ∑_{i} ∑_{j} (\frac{∂^{2}G}{∂m_{i}∂m_{l}} ) δm_{i}δm_{j}.  (2.1)

ここで理想気体の化学ポテンシャルの定義式を用いると, (2.1) 式の右辺第 2 項のカッコ内は以下のように書ける.

(\frac{∂G}{∂m_{i}} ) = μ_{i} = μ_{i}^{∘}(T) + RT ln (\frac{m_{i}}{∑ m_{i}} ) + RT ln (\frac{p}{p_{0}} ) .  (2.2)

但し p_{0} は基準圧力で通常 1.0 × 10^{5} Pa である. (2.2) 式を用いて (2.1) 式の右辺第 3 項の括弧内は以下のように書ける.

(\frac{∂^{2}G}{∂m_{i}∂m_{l}} ) = RT (\frac{δ_{il}}{m_{i}} − \frac{1}{∑ m_{i}} ) .  (2.3)

(2.2), (2.3) 式を (2.1) 式に代入し, RT で割ることによって,

\frac{Q(n_{i})}{RT} = \frac{G(m_{i})}{RT} + ∑_{i} \frac{μ_{i}}{RT} δm_{i} + \frac{1}{2} [∑_{i} \frac{(δm_{i})^{2}}{m_{i}} − \frac{\{∑(δm_{i})\}^{2}}{∑ m_{i}} ]  (2.4)

が導かれる.

(2.4) 式の Q(n_{i})/RT と (1.3) 式で気相のみ考えた場合の元素数の保存式

∑_{i} a_{ie}n_{i} = B_{e}  (2.5)

を用いて, ラグランジュの未定係数法を用いて極値の条件を与える. L を

L ≡ \frac{G(m_{i})}{RT} + ∑_{i} \frac{μ_{i}}{RT} δm_{i} + \frac{1}{2} [∑_{i} \frac{(δm_{i})^{2}}{m_{i}} − \frac{\{∑(δm_{i})\}^{2}}{∑ m_{i}} ] − ∑_{e} χ_{e} [∑_{i} a_{ie}n_{i} − B_{e}]

と定義すると, (2.4) 式が極値をとる条件は ∂L/∂n_{i} = 0 なので,

\frac{μ_{i}}{RT} + \frac{n_{i}}{m_{i}} − \frac{∑ n_{i}}{∑ m_{i}} − ∑_{e} χ_{e}a_{ie} = 0  (2.6)

となる. (2.6) 式と元素数保存の式 (2.5) を連立することで, 変数 n_{i}, n, χ_{e}(ラグランジュの未定係数) に関する N + M 個の線形方程式が得られ, 方程式系が閉じる.

(2.6) 式は式変形を行うことで連立方程式の次数を減らすことができる. まず (2.6) 式を n_{i} に関する式に書き直す.

n_{i} = m_{i} \{−\frac{μ_{i}}{RT} + \frac{∑ n_{i}}{∑ m_{i}} + ∑_{e} χ_{e}a_{ie}\} .  (2.7)

この (2.7) 式を i に関して 1 から N までの和をとり, 元素数の保存式 (2.5) 式を用いると,

∑_{e} B_{e}χ_{e} = ∑_{i} m_{i} \frac{μ_{i}}{RT}  (2.8)

となる. また (2.7) 式の両辺に a_{if} をかけて i に関して和を取り, 元素数の保存式 (2.5) を用いると,

∑_{e} \{∑_{i} a_{ie}a_{if}m_{i}\} χ_{e} + B_{f} (\frac{∑ n_{i}}{∑ m_{i}} − 1) = ∑_{i} a_{if}m_{i} \frac{μ_{i}}{RT}  (2.9)

となる. a_{if} の添字 f (0 ≤ f ≤ N) はダミーサフィックスで, 実質的には e と同等である.

式変形の結果得られた (2.8), (2.9) 式は M + 1 個の方程式から成り, それらを連立することで χ_{e} と (∑ n_{i}/ ∑ m_{i} − 1) の値が求まる. χ_{e} と (∑ n_{i}/ ∑ m_{i}) の値を (2.7) 式に代入することで各化学種のモル数が得られる. 得られた n_{i} を m_{i} とみなして反復的に計算を繰り返すことで, n_{i} を平衡組成へ収束させる.

この式変形の利点は, N + M + 1 個の方程式を連立すべきところを M + 1 個の方程式の連立で済ませられる点にある. この変形によって逆変換する行列の大きさが M + 1 となるので, 計算量や計算時間を短縮できるものと期待される.

最終的に解くべき式をまとめる.

[figure]

温度・圧力一定, 初期状態での化学種のモル数 m_{i}, 理想気体を仮定した場合, 平衡状態での化学種のモル数 n_{i} は以下の式から得られる.

n_{i} = m_{i} \{−\frac{μ_{i}}{RT} + \frac{∑ n_{i}}{∑ m_{i}} + ∑_{e} χ_{e}a_{ie}\} .

但し, 上式に含まれる変数 χ_{e}, ∑ n_{i} は以下の 2 式を連立させることで得られる.

∑_{e} B_{e}χ_{e} = ∑_{i} m_{i} \frac{μ_{i}}{RT} ,

∑_{e} \{∑_{i} a_{ie}a_{if}m_{i}\} χ_{e} + B_{f} (\frac{∑ n_{i}}{∑ m_{i}} − 1) = ∑_{i} a_{if}m_{i} \frac{μ_{i}}{RT} .

[figure]

(2.8), (2.9) 式を行列式の形式で書けば,

[ \begin{matrix}∑_{i} a_{i1}a_{i1}m_{i} & ∑_{i} a_{i2}a_{i1}m_{i} & ⋯ & ∑_{i} a_{ie}a_{i1}m_{i} & B_{1} \\ ∑_{i} a_{i1}a_{i2}m_{i} & ∑_{i} a_{i2}a_{i2}m_{i} & ⋯ & ∑_{i} a_{ie}a_{i2}m_{i} & B_{2} \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ ∑_{i} a_{i1}a_{if}m_{i} & ∑_{i} a_{i2}a_{if}m_{i} & ⋯ & ∑_{i} a_{ie}a_{if}m_{i} & B_{f} \\ B_{1} & B_{2} & ⋯ & B_{e} & 0 \end{matrix} ] [ χ_{1} ; χ_{2} ; ⋮ ; χ_{e} ; (^{\frac{∑ n_{i}}{∑ m_{i}}} − 1) ] = [ ∑ a_{i1}m_{i} \frac{μ_{i}}{RT} ; ∑ a_{i2}m_{i} \frac{μ_{i}}{RT} ; ⋮ ; ∑ a_{if}m_{i} \frac{μ_{i}}{RT} ; ∑ m_{i} \frac{μ_{i}}{RT} ]

と書くことができる.

2.2 多相系の場合

前章と同様に, 適当な物質存在量 m_{i}^{ϕ} のまわりで 2 次のテーラー展開を行い, その極値を求める (m_{i}^{ϕ} においても質量保存を満足していると仮定する).

Q(n_{i}^{ϕ}) = G(m_{i}^{ϕ}) + ∑_{ϕ} ∑_{i} \underbrace{(\frac{∂G}{∂m_{i}^{ϕ}} )}_{A} δm_{i}^{ϕ} + \frac{1}{2} ∑_{ϕ} ∑_{i} ∑_{j} \underbrace{(\frac{∂^{2}G}{∂m_{i}^{ϕ}∂m_{l}^{ϕ}} )}_{B} δm_{i}^{ϕ}δm_{j}^{ϕ}.  (2.10)

ここで (2.10) 式の A 項は (1.2) 式を用いることによって,

(\frac{∂G}{∂m_{i}^{ϕ}} ) ≡ μ_{i}^{ϕ} = μ_{i}^{∘}^{ϕ}(T) + RT ln (\frac{m_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) + αRT ln (\frac{p}{p_{0}} )  (2.11)

と書ける. また (2.10) 式の B 項を (1.2) 式を用いて書き下すと,

(\frac{∂^{2}G}{∂m_{i}^{ϕ}∂m_{l}^{ϕ}} ) = RT (\frac{δ_{il}}{m_{i}^{ϕ}} − \frac{1}{∑_{i} m_{i}^{ϕ}} )  (2.12)

となる. (2.11), (2.12) 式を (2.10) 式に代入し, RT で割ることによって,

\frac{Q(n_{i}^{ϕ})}{RT} = ∑_{ϕ} \frac{G(m_{i}^{ϕ})}{RT} + ∑_{ϕ} ∑_{i} \frac{μ_{i}^{ϕ}}{RT} δm_{i} + \frac{1}{2} ∑_{ϕ} [∑_{i} \frac{(δm_{i}^{ϕ})^{2}}{m_{i}^{ϕ}} − \frac{\{∑_{i}(δm_{i}^{ϕ})\}^{2}}{∑_{i} m_{i}^{ϕ}} ]  (2.13)

\overline{main.tex 2004/05/29(杉山耕一朗)}

が導かれる.

(2.13) 式の Q(n_{i})/RT と元素数の保存式 (1.3) より, ラグランジュの未定係数法を用いて極値の条件を与える. L を

L ≡  ∑_{ϕ} \frac{G(m_{i}^{ϕ})}{RT} + ∑_{ϕ} ∑_{i} \frac{μ_{i}^{ϕ}}{RT} δm_{i} + \frac{1}{2} ∑_{ϕ} [∑_{i} \frac{(δm_{i}^{ϕ})^{2}}{m_{i}^{ϕ}} − \frac{\{∑_{i}(δm^{ϕ})\}^{2}}{∑_{i} m^{ϕ}} ]

− ∑_{e} χ_{e} ∑_{ϕ} [∑_{i} a_{ie}^{ϕ}n_{i}^{ϕ} − B_{e}^{ϕ}]

と定義すると, 極値をとる条件は ∂L/∂n_{i}^{ϕ} = 0 なので,

\frac{μ_{i}^{ϕ}}{RT} + \frac{n_{i}^{ϕ}}{m_{i}^{ϕ}} − \frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} − ∑_{e} χ_{e}a_{ie}^{ϕ} = 0  (2.14)

である. (2.14) 式と (1.3) 式を連立することで, 変数 n_{i}^{ϕ}, χ_{e} (ラグランジュの未定係数) に関する ΦN + M 個の線形方程式が得られ, 方程式系が閉じる.

気相のみの場合と同様に (2.14) 式は式変形を行うことで連立方程式の次数を減らすことができる. まず (2.14) 式を n_{i}^{ϕ} に関する式に書き直す.

n_{i}^{ϕ} = m_{i}^{ϕ} \{−\frac{μ_{i}^{ϕ}}{RT} + \frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} + ∑_{e} χ_{e}a_{ie}^{ϕ}\} .  (2.15)

この (2.15) 式を i に関して 1 から N までの和をとり, 元素数の保存式 (1.3) を用いると,

∑_{i} ∑_{e} a_{ie}^{ϕ}m_{i}^{ϕ}χ_{e} = ∑_{i} m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT}  (2.16)

となる. また (2.15) 式の両辺に a_{if}^{ϕ} をかけて i と ϕ に関して和を取り, 元素数の保存式 (1.3) を用いると,

∑_{e} \{∑_{i} ∑_{ϕ} a_{ie}^{ϕ}a_{if}^{ϕ}m_{i}^{ϕ}\} χ_{e} + ∑_{ϕ} ∑_{i} a_{if}^{ϕ}m_{i}^{ϕ} (\frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) = ∑_{i} ∑_{ϕ} a_{if}^{ϕ}m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT} + B_{f} (2.17)

である. a_{if} の添字 f (0 ≤ f ≤ M) はダミーサフィックスで, 実質的には e と同等である.

(2.16), (2.17) 式は M + Φ の方程式であり, また変数の数は χ_{e} (M 個) と ∑_{i} n_{i}^{ϕ} (Φ 個) なので, この 2 式を連立することによって方程式が閉じる. 得られた n_{i}^{ϕ} を m_{i}^{ϕ} とみなして反復的に計算を繰り返すことで, n_{i}^{ϕ} を平衡組成へ収束させる.

最後に解くべき式は以下のようになる.

[figure]

温度・圧力一定, 初期状態での化学種のモル数 m_{i}^{ϕ}, 理想気体を仮定した場合, 平衡状態での化学種のモル数は以下の式から得られる.

n_{i}^{ϕ} = m_{i}^{ϕ} \{−\frac{μ_{i}^{ϕ}}{RT} + \frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} + ∑_{e} χ_{e}a_{ie}^{ϕ}\} .

但し, 上式に含まれる変数 χ_{e}, n^{ϕ} は以下の 2 式を連立させることで得られる.

∑_{i} ∑_{e} a_{ie}^{ϕ}m_{i}^{ϕ}χ_{e} = ∑_{i} m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT}

∑_{e} \{∑_{i} ∑_{ϕ} a_{ie}^{ϕ}a_{if}^{ϕ}m_{i}^{ϕ}\} χ_{e} + ∑_{ϕ} ∑_{i} a_{if}^{ϕ}m_{i}^{ϕ} (\frac{∑_{i} n_{i}^{ϕ}}{∑_{i} m_{i}^{ϕ}} ) = ∑_{i} ∑_{ϕ} a_{if}^{ϕ}m_{i}^{ϕ} \frac{μ_{i}^{ϕ}}{RT} + B_{f}

[figure]

数値的最適化法: RAND 法第 3 章 Newton 法との比較  8

第3章 Newton 法との比較

\overline{main.tex 2004/05/29(杉山耕一朗)}