% 表題: RAND 法
%
% 履歴: 2004-05-29 杉山耕一朗

%%%
%%% Newton 法の探査ベクトルとの関連
%%%

\chapter{RAND 法}

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

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



\section{気相のみの場合}
\label{gas_phase}

適当な物質存在量を $m_{i}$ とする.  $m_{i}$ よりも熱化学平衡に近い物質存
在度を $n_{i} (= m_{i} + \delta m_{i})$ とする.  平衡時の熱力学ポテンシャ
ル $G$ を $m_{i}$ のまわりでテーラー展開し, $\delta m_{i}$ に関する 2 次
の微小量まで近似式 $Q(n_{i})$ を作る.
%
\begin{eqnarray}
Q(n_{i}) = G(m_{i}) + \sum_{i} \left( \DP{G}{m_{i}}\right) \delta m_{i}
+ \Dinv{2} \sum_{i} \sum_{j} \left( \frac{\partial^{2}G}{\partial m_{i}\partial m_{l}}\right) \delta m_{i} \delta m_{j}.  
\Deqlab{3.3.24}
\end{eqnarray}
%
ここで理想気体の化学ポテンシャルの定義式を用いると, \Deqref{3.3.24} 式の
右辺第 2 項のカッコ内は以下のように書ける.
%
\begin{eqnarray}
 \left(\DP{G}{m_{i}}\right) 
= \mu_{i}
= \mu_{i}^{\circ}(T) 
  + RT \ln{\left( \frac{m_{i}}{\sum m_{i}}\right)}
  + RT \ln\left(\frac{p}{p_{0}}\right).
\Deqlab{3.3.25}
\end{eqnarray}
%
但し $p_{0}$ は基準圧力で通常 $1.0 \times 10^{5}$ Pa である. 
\Deqref{3.3.25} 式を用いて \Deqref{3.3.24} 式の右辺第 3 項の括弧内は以
下のように書ける. 
%
\begin{eqnarray}
\left( \frac{\partial^{2}G}{\partial m_{i}\partial m_{l}}\right)
 = RT \left( \frac{\delta_{il}}{m_{i}} - \Dinv{\sum m_{i}}\right). 
\Deqlab{3.3.27}
\end{eqnarray}
%
\Deqref{3.3.25}, \Deqref{3.3.27} 式を 
\Deqref{3.3.24} 式に代入し, $RT$ で割ることによって, 
%
\begin{eqnarray}
\frac{Q(n_{i})}{RT} = \frac{G(m_{i})}{RT} 
 + \sum_{i}  \frac{\mu_{i}}{RT} \delta m_{i}
 + \Dinv{2} \left[ \sum_{i}\frac{(\delta m_{i})^{2}}{m_{i}}
 - \frac{\{\sum (\delta m_{i})\}^{2}}{\sum m_{i}}\right]
\Deqlab{3.3.28}
\end{eqnarray}
%
が導かれる. 

\Deqref{3.3.28} 式の $Q(n_{i})/RT$ と \Deqref{2} 式で気相の
み考えた場合の元素数の保存式
%
\begin{eqnarray}
 \sum_{i} a_{ie}n_{i} = B_{e}
  \Deqlab{mass-conservation_gas}
\end{eqnarray}
%
を用いて, ラグランジュの未定係数法を用いて極値の条件を与える.  $L$ を
%
\begin{eqnarray}
L \equiv
\frac{G(m_{i})}{RT} 
+ \sum_{i} \frac{\mu_{i}}{RT} \delta m_{i}
+ \Dinv{2} \left[ \sum_{i}\frac{(\delta m_{i})^{2}}{m_{i}} 
   - \frac{\{\sum (\delta m_{i})\}^{2}}{\sum m_{i}}\right]  
- \sum_{e} \chi_{e} 
 \left[\sum_{i} a_{ie}n_{i} - B_{e}\right] \nonumber 
\end{eqnarray}
%
と定義すると, \Deqref{3.3.28} 式が極値をとる条件は 
$\partial L / \partial n_{i} = 0$ なので, 
%
\begin{eqnarray}
 \frac{\mu_{i}}{RT}
+  \frac{n_{i}}{m_{i}} 
- \frac{\sum n_{i}}{\sum m_{i}}
- \sum_{e} \chi_{e} a_{ie} = 0
\Deqlab{3.2.30}
\end{eqnarray}
%
となる. \Deqref{3.2.30} 式と元素数保存の式 \Deqref{mass-conservation_gas} を連立することで, 
変数 $n_{i}$, $n$, $\chi_{e}$(ラグランジュの未定係数)に関する N + M 個の
線形方程式が得られ, 方程式系が閉じる.



\Deqref{3.2.30} 式は式変形を行うことで連立方程式の次数を減らすことができる. 
まず \Deqref{3.2.30} 式を $n_{i}$ に関する式に書き直す.
%
\begin{eqnarray}
 n_{i} = 
  m_{i} \left\{
	 - \frac{\mu_{i}}{RT} 
	 +  \frac{\sum n_{i}}{\sum m_{i}} 
	 + \sum_{e} \chi_{e} a_{ie}
  \right\}. 
\Deqlab{3.3.31}
\end{eqnarray}
%
この \Deqref{3.3.31} 式を $i$ に関して 1 から $N$ までの和をとり, 
元素数の保存式 \Deqref{mass-conservation_gas} 式を用いると, 
%
\begin{eqnarray}
 \sum_{e}B_{e}\chi_{e} = \sum_{i} m_{i} \frac{\mu_{i}}{RT}
\Deqlab{3.3.33}
\end{eqnarray}
%
となる. また \Deqref{3.3.31} 式の両辺に $a_{if}$ をかけて $i$ に関して和
を取り, 元素数の保存式 \Deqref{mass-conservation_gas} を用いると, 
%
\begin{eqnarray}
\sum_{e} \left\{ \sum_{i} a_{ie} a_{if} m_{i} \right\} \chi_{e}
 + B_{f} \left( \frac{\sum n_{i}}{\sum m_{i}} - 1 \right) 
 = \sum_{i} a_{if} m_{i} \frac{\mu_{i}}{RT}
\Deqlab{3.3.34-2}
\end{eqnarray}
%
となる.  $a_{if}$ の添字 $f$ ($0 \leq f \leq N$)はダミーサフィックスで, 
実質的には $e$ と同等である.


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

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


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

\begin{screen}
温度・圧力一定, 初期状態での化学種のモル数 $m_{i}$, 理想気体を仮定した場
合, 平衡状態での化学種のモル数 $n_{i}$ は以下の式から得られる.
%
\begin{eqnarray}
 n_{i} = 
  m_{i} \left\{
	 - \frac{\mu_{i}}{RT} 
	 + \frac{\sum n_{i}}{\sum m_{i}}
	 + \sum_{e} \chi_{e} a_{ie}
  \right\}. \nonumber
\end{eqnarray}
%
但し, 上式に含まれる変数 $\chi_{e}, \sum n_{i}$ は
以下の 2 式を連立させることで得られる. 
%
\begin{eqnarray}
&& \sum_{e}B_{e}\chi_{e} = \sum_{i} m_{i} \frac{\mu_{i}}{RT},
\nonumber \\
&& \sum_{e} \left\{ \sum_{i} a_{ie} a_{if} m_{i} \right\} \chi_{e}
 + B_{f} \left( \frac{\sum n_{i}}{\sum m_{i}} - 1 \right) 
 = \sum_{i} a_{if} m_{i} \frac{\mu_{i}}{RT}. \nonumber
\end{eqnarray}
%
\end{screen}

\Deqref{3.3.33}, \Deqref{3.3.34-2} 式を行列式の形式で書けば, 
%
\[
  \left[
   \begin{array}{ccccc}
    \sum_{i} a_{i1} a_{i1} m_{i} & \sum_{i} a_{i2} a_{i1} m_{i} 
     & \cdots & \sum_{i} a_{ie} a_{i1} m_{i} & B_{1} \\
    \sum_{i} a_{i1} a_{i2} m_{i} & \sum_{i} a_{i2} a_{i2} m_{i} 
     & \cdots & \sum_{i} a_{ie} a_{i2} m_{i} & B_{2} \\
    \vdots & \vdots & \ddots & \vdots & \vdots \\
    \sum_{i} a_{i1} a_{if} m_{i} & \sum_{i} a_{i2} a_{if} m_{i} 
     & \cdots & \sum_{i} a_{ie} a_{if} m_{i} & B_{f} \\
    B_{1} & B_{2} & \cdots & B_{e} & 0 
   \end{array}
 \right]
 \left[
   \begin{array}{c}
    \chi_{1} \\
    \chi_{2} \\
    \vdots \\
    \chi_{e} \\
    \left(\frac{\sum n_{i}}{\sum m_{i}} - 1 \right)
   \end{array}
 \right]
= 
 \left[
   \begin{array}{c}
    \sum a_{i1} m_{i} \frac{\mu_{i}}{RT} \\
    \sum a_{i2} m_{i} \frac{\mu_{i}}{RT} \\
    \vdots \\
    \sum a_{if} m_{i} \frac{\mu_{i}}{RT} \\
    \sum m_{i} \frac{\mu_{i}}{RT} \\
   \end{array}
 \right]
\]
%
と書くことができる. 
