\underline{最適化法によるギブス自由エネルギー最小化: RAND 法第 1 章はじめに 1}

第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 の個数を意味する.

第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_{j}} ) δ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_{j}} ) = \frac{∂}{∂m_{j}} \{μ_{i}^{∘}(T) + RT ln (\frac{m_{i}}{∑ m_{i}} ) + RT ln (\frac{p}{p_{0}} )\}
= RT \frac{∂}{∂m_{j}} \{ln (\frac{m_{i}}{∑ m_{i}} )\}
= RT \frac{∑ m_{i}}{m_{i}} \frac{∂}{∂m_{j}} (\frac{m_{i}}{∑ m_{i}} )
= RT \frac{∑ m_{i}}{m_{i}} \{\frac{δ_{ij}}{∑ m_{i}} − \frac{m_{i}}{(∑ m_{i})^{2}} \}
= 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 となるので, 計算量や計算時間を短縮できるものと期待される.

(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} ]

と書くことができる.

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

[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.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)

が導かれる.

(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 法との比較 7

第3章 Newton 法との比較

RAND 法は「制約つき最適化法」の「Newton 法」に分類される最適化法であるが, 教科書的な Newton 法の表記とは多少異なる. 本章では RAND 法を Newton 法と比較し, その表記の差異が何に起因するか考察する.

3.1 Newton 法

g(x) = (g_{1}(x) = 0, ⋯, g_{r}(x))^{T} = 0 という等式制約のみの場合を考える. このとき, 制約想定下で最適解 x^{∗} に対して λ^{∗} が存在し, これらは Karush-Kuhn-Tucker 条件, すなわち変数 x, λ に連立非線形方程式,

∇_{x}L(x, λ) = ∇_{x}f(x) − (J(x))^{T} λ = 0  (3.1)
∇_{λ}L(x, λ) = g(x) = 0  (3.2)

を満たす. ただし λ は未定乗数, L(x, λ) ≡ f(x) − ∑ λ_{j}g_{j}(x) とし, Jacobi 行列は J(x) ≡ ∇g(x) とする. この Karush-Kuhn-Tucker 条件はラグランジュの未定乗数法の条件と同じものである. 制約つき最適化法は究極的にはこの方程式を解かねばならないと言える.

Newton 法を (3.1), (3.2) 式から成る連立方程式に適用することで次のアルゴリズムが得られる. 探索方向 d_{k} を連立一次方程式,

[ \begin{matrix}H(x_{k}, λ_{k}) & −J^{T} (x_{k}) \\ J(x_{k}) & 0 \end{matrix} ] [ \begin{matrix}d & k \\ λ & k+1 \end{matrix} ] = − [ \begin{matrix}∇_{x}f(x_{k}) \\ g(x_{k}) \end{matrix} ]

によって定め, ステップ幅を常に α = 1 にとして点列 {x_{k}}, {λ_{k}} を生成する. λ_{k+1} は連立方程式の解として直接得られることに注意を払う必要がある. ただし H(x, λ) は,

H(x, λ) = ∇_{x}^{2}L(x, λ) = ∇_{x}^{2}f(x) − ∑ λ_{i}∇_{x}^{2}g_{i}(x)  (3.3)

とする.

3.2 RAND 法と Newton 法の比較

RAND 法と前節の Newton 法を比較すると,

• RAND 法では目的関数をギブス自由エネルギー G 自体を x に関する 2 次式として展開するのに対し, Newton 法ではラグランジュ未定乗数法による関数 L を展開する.

最適化法によるギブス自由エネルギー最小化: RAND 法

第 3 章 Newton 法との比較 8

• 解くべく連立方程式が, 探索ベクトル d ではなく最適解 x に対する方程式系になっている.

という差異があることがわかる. ギブス自由エネルギーの最小化の場合には等式制約が x に関する 1 次関数になっているため, ギブス自由エネルギー G だけを x に関して展開すれば十分なのであり, 等式制約をそのまま使える x に関する方程式系として整理する方が都合が良いのである.

ギブス自由エネルギーを最小化する問題について, Newton 法の形式で書くことを試みる. (3.1) 式は,

H(x_{k})d_{k} − J(x_{k})λ = −∇G(x_{k})  (3.4)

となる. ただし J(x) = ∇g(x) = a_{ie} であり, 等式制約の x に関する 2 階微分はゼロなので, ヘッシアンは

H(x, λ) = H(x) = ∇^{2}G(x)  (3.5)

である. (3.2) 式は,

∇g(x_{k})d_{k} = −g(x_{k})  (3.6)

である. したがって前章の定式化を利用すると,

[ \begin{matrix}RT (\frac{δ_{il}}{m_{i}} − \frac{1}{∑ m_{i}} ) & −a \\ a & 0 \end{matrix} ] [ \begin{matrix}d & k \\ λ & k+1 \end{matrix} ] = − [ \begin{matrix}μ^{∘} + RT ln \begin{matrix}x & k \\ ∑ & x & k \end{matrix} + RT ln p/p_{0} \\ a & n & k \end{matrix} ]

となる.

\underline{最適化法によるギブス自由エネルギー最小化: RAND 法第 4 章参考文献 9}

第4章参考文献

Van Zeggerren, F., Storey, S.H., 1970:

The computation of chemical equilibria, Cambridge Univ. Press.

White, W.B., Johnson, S.M., Dantzig, G.B., 1958: Chemical equilibrium in complex mixture. J. Chem. Phys. 28, 751–755

藤田宏, 今野浩, 田辺國士, 著:

「最適化法」 , 岩波講座応用数学方法 7, 岩波書店, 1994.