% 表題: 化学平衡計算の数値解法の歴史
%
% 履歴: 2004/03/18 杉山耕一朗


\chapter{数値計算手法の変遷}

化学平衡計算の数値的計算手法の変遷を簡単にまとめる. 計算手法は, 計算機の
能力の発展, および応用数学的な最適化問題の数値的取り扱い方法の確立ととも
に進化してきた. 

平衡問題の扱いが難しいのは, (1)方程式が非線形, (2)熱力学関数の極値が 1 
つとは限らないためである. 特に (2) が重要な問題で, しばしば局所的な最小
値(local minimum)に収束させてしまい, 大局的な最小値(global minimun)に収
束しないことが挙げられる. global minimum が見付かるという保証は全く無い.
そのため, 熱力学関数の最小値を見付けるという最適化手法の議論に留まらず, 
計算の初期値の与え方, 得られた解が local minimum か否かの判定にも努力が
払われてきた. 

以下, その変遷を簡単にまとめる. 


\section{60 年代まで}

60 年代までの化学平衡計算は, 例えば Van Zeggeren and Storey (1970) にま
とめられている. 当時の平衡計算は, 相平衡において化学ポテンシャルは一定, 
もしくは平衡状態ではギブスの自由エネルギーが最小化される, という条件を用
いている.

当時の化学平衡計算は, Van Zeggeren and Storey (1970) の言い方を借りると, 
以下に大別される. 計算機の能力の向上とともに, 最適化手法が現実的な計算方
法となったようだ.
%
\begin{description}
 \item[最適化手法(optimization technique)]~\\ 
%
	    ある温度・圧力下でギブス自由エネルギーが最小化される平衡組成
	    を最適化手法を用いて計算する方法である. 当時は最適化手法が様々
	    に開発されていた段階で, 多くの最適化手法が化学平衡計算に適用
	    されていた. 
%
 \item[非線形方程式の解(solution of non-linear equations)] ~\\
%
	    相平衡において各相の化学ポテンシャルは等しいという条件のもと
	    で平衡組成を計算する方法である. 具体的には, mass action
	    equation と mass balance equation の式を解く事で平衡組成を計
	    算する. 
%
\end{description}
%

当時の代表的な論文として以下が挙げられる. 
%
\begin{description}
 \item[Anthony and Himmelblau (1963)]~\\ 
%
	    Hooke and Jeeves (1961) の開発した最適化手法である「Direct
	    Search」を化学平衡問題に適用した. これは組成を順繰りに変化さ
	    せて熱力学関数が最小化される向きを調べる. 
%
 \item[White {\it et al}. (1958)]~\\ 
%
	    熱力学関数を適当な初期値のまわりで組成に関する 2 次関数に近
	    似し, その極小値を与える組成を求める. 以上の操作を反復的に繰
	    り返すことで平衡組成に収束させる. 彼らの手法は「RAND method」
	    と呼ばれる.
%
 \item[Kandiner and Brinkley (1950)]~\\
%
	    化学反応を想定し, 系のギブス自由エネルギーを化学量論係数
	    ($\nu_{i}$), 反応度($\varepsilon$), 化学ポテンシャル
	    ($\mu_{i}$) を用いて表現する. そのときの平衡条件は
	    $\partial G/ \partial \varepsilon = 0$ となる. 
	    その平衡条件を解く事で求まる圧平衡定数と化学ポテンシャルの関
	    係式(mass action equation)と, 物質量保存(mass balance
	    equation) を連立し, 平衡組成を反復的に求める. 

	    彼らは Hildebrand (1956) によって開発された非線形方程式の数
	    値的解法を化学平衡計算に応用した. 対象となる関数($\partial
	    G/\partial \varepsilon$)をテーラー展開して 1 次の微小量まで
	    考慮する.  そして関数の値が収束するまで反復計算を繰り返
	    す. 関数の傾きから関数が収束する方向を判定する方法といえる. 
	    彼らの手法は「Newton-Rapson method」と呼ばれる.
%	    
 \item[Huff {\it et al} (1951)]~\\ 
%
	    基本的に Kandiner and Brinkley (1950) と同じやり方.
	    Kandiner and Brinkley (1950) との違いは, (1) 気体元素を成分
	    (component)とみなす, (2)系の熱力学的状態を $p, V, T, H, S$
	    から 2 つ選ぶことで記述する, ことである. この方法は「NASA
	    method」と呼ばれる. 
%
\end{description}




\section{70 年代から 80 年代}

この時代は, 「古典的な(classic)」な平衡計算手法である, 相間での化学ポテ
ンシャルが等しいという条件の下での平衡計算が多数行われていた. 
%
\begin{description}
%
 \item[Fussel (1979), Fussel and Yanosik (1978)]~\\
%
	    Newton-Rapshon 法(NVNR)を用いて 2 相系を計算
%
 \item[Mehra, Heidemann, and Aziz (1982)]~\\
%
	    Newton iteration を用いた accelerated SSI 
%
 \item[Nghiem and Li (1984)]~\\
%
	    Newton-Raphson 法と準 SSI 法を併用
%
 \item[Rishes and Dalen (1984)]~\\
%
	    SSI method
%
\end{description}
%

そのような流れにあって, 相平衡計算のマイルストーンというべき論文は, 
Michelsen (1984a, b)であった. 

%
\begin{description}
%
 \item[Michelsen (1984a, b)]~\\
%
	    接平面距離(tangent plane distance; TPD)の安定点を見付ける方
	    法を確立し, 2 次の方法(second-order method)に応用した. 
	    相平衡計算のマイルストーンというべき論文. 
%
\end{description}
%

応用数学としての最適化問題が開発されたのは 70 年代後半である. 当時, メイ
ンフレームで解くべき重要な問題が最適化問題であった. 最適化手法として現在
もよく用いられている Newton 法, quesi-Newton 法の数値計算パッケージが整
備されたのもこの時代だと思われる. White {\it et al} (1958) による RAND 
法は基本的に Newton 法と同じものである. この頃から整備された数値計算パッ
ケージが使われるようになり, RAND 法の名前は論文中から消えている.


\section{80 年代後半から現在}

80 年代後半からは, 最適化手法そのものよりも, local minimum を回避する方
法の模索に力点が置かれるようになった. 最適化手法自体は, Newton 法,
quesi-Newton 法といった定番とも言える方法が確立したためである. 

この時代に新たに取り入れられたことは, 確率的な手法である. モンテカルロ法
が使われたり, また local minimum を見付けたあとに確率的に組成を振って別
のポテンシャルの谷を探す方法が取り入れられた.

\begin{description}
%
 \item[Trangenstein (1985)]~\\
%
	    Two-phase equilibria reservoir compositional simulation with constrained.
%
 \item[Vickery and Taylor (1986), DeGance (1993)]~\\
%
	    path continuation methods.
%
 \item[Sun and Seider (1995)]~\\
%
	    homotopy continuation method.
%
 \item[Nagarayan {\it et al}. (1991)]~\\
%
	    modified ``semi-local'' stability test and higher order
	    initialization procedure
%
 \item[Eubank {\it et al}. (1992)]~\\
%
	    integration of the Gibbs energy surface
%
 \item[Bullard and Biegler (1993), Han and Rangaiah (1996)]~\\
%
	    panalty-based approch. 線形計算手法を活用. 
%
 \item[Hua {\it et al}. (1996, 1998, 1999)]~\\
%
	    Interval Newton method
%
\end{description}


Nichita {\it et al}. (2001) の言葉を借りると, 最近の論文で特筆すべきもの
は以下の通りである. 

\begin{description}
%
 \item[Simulated annealing (SA) アルゴリズム]~\\
%
	    Corana {\it et al}. (1987) によって開発

	    Pan and Firoozabadi (1998) によって, 多くの相平衡計算に応用
	    された. 

	    Nichita {\it et al}. (2001) によって, Wax precipitation.

	    Zhu (2000) によって相安定性解析
%
 \item[McDonald and Floudas (1995 a, b, c, 1997)]~\\
%
	    決定論的な branch と bound method を用いる. 
	    この方法は「GLOPEQ コード」と呼ばれ, global minimum に収束す
	    ることを保証するが, 計算コストが膨大である. 
%
 \item[heuristic stochastic method (発見的な確立的方法)]~\\
%
	    based on space approch. この方法では傾きを評価する必要はない. 
	    収束する可能性がある一般的なアルゴリズムとして用いられている. 
	    収束するのに時間がかかる.
%
 \item[Tunneling gloal optimization method]
%
	    適当な初期値から計算を開始して quesi-Newton 法を用いて local
	    minimum を見付ける.  そのあとに確率的な手法でもって, さらに
	    深いポテンシャルの谷が存在しないか検討する.

	    Gometz and Romero (1994), Gomez {\it et al}. (1998, 2001a,
	    b), Nichita {\it et al.} (2002) がある. 
%
\end{description}

