高次の最適化手法

臨界点付近での相挙動の評価など数値的に平衡状態を決定することが難しい場合には、フュガシティの一致を利用したSSIニュートン法ではなく、より高次の最適化手法を利用しなければならない状況が稀に生じる。そのような場合はフュガシティが一致する点を探す求根問題ではなく、ギブズエネルギーが最小になる点を探す最小化に置き換えて、二次の最適化アルゴリズムを利用する。

独立変数を液相のモル分率n_{iL}と仮定すると、ギブズエネルギー、ギブズエネルギーの勾配ベクトルgradient vector、ヘッセ行列Hessian matrixは次のように与えられる。

(106)   \begin{equation*} \overline{G}=\frac{G}{RT}=\sum^{n_c}_{i=1} \left( n_{iL}\ln{f_{iL}} +n_{iV}\ln{f_{iV}}\right) \end{equation*}

(107)   \begin{equation*} g_i=\frac{\partial\overline{G}}{\partial{n_{iL}}}\sum^{n_c}_{i=1} \left( n_{iL}\ln{f_{iL}} +n_{iV}\ln{f_{iV}}\right) \end{equation*}

(108)   \begin{equation*} H_{i,j}=\frac{\partial{g_i}}{\partial{n_{jL}}} =\frac{\partial\ln{f_{iL}}}{\partial{n_{jL}}}+\frac{\partial\ln{f_{iV}}}{\partial{n_{jV}}} =\frac{1}{f_{iL}}\frac{\partial{f_{iL}}}{\partial{n_{jL}}}+\frac{1}{f_{iV}}\frac{\partial{f_{iV}}}{\partial{n_{jV}}} \end{equation*}

ただし、勾配ベクトルとヘッセ行列については各成分のみを抜き出してある。なお、気相のモル分率n_{iV}を独立変数であると仮定した場合のヘッセ行列は次のようになる。

(109)   \begin{equation*} H_{i,j}=\frac{\partial{g_i}}{\partial{n_{jV}}} =-\frac{\partial{g_i}}{\partial{n_{jL}}} =-\left( \frac{1}{f_{iL}}\frac{\partial{f_{iL}}}{\partial{n_{jL}}}+\frac{1}{f_{iV}}\frac{\partial{f_{iV}}}{\partial{n_{jV}}} \right) \end{equation*}

さて、ギブズエネルギーの最小化問題を解くにはまず次のような線形方程式を解いてニュートン方向Newton directionなるものを決める必要がある。

(110)   \begin{equation*} H\Delta \overrightarrow{n}_L=\overrightarrow{g} \end{equation*}

ここでの\Delta\overrightarrow{n}_Lがニュートン方向であり更新ベクトルに相当する。ニュートン方向が求まればギブズエネルギーの最小値は次の漸化式によって逐次的に求めることができる。

(111)   \begin{equation*} \overrightarrow{n}_L^{n+1}=\overrightarrow{n}_L^{n}+\Delta \overrightarrow{n}_L \end{equation*}

収束の判定は例えば次のような条件を確認すればよい。

(112)   \begin{equation*} \left|f_{iL}-f_{iV} \right| <\epsilon \quad \forall{i} \end{equation*}

ただし、\epsilonは許容誤差である。