﻿ 基于最小自由能法的理想气相反应平衡新算法
 西南石油大学学报(自然科学版)  2018, Vol. 40 Issue (1): 122-129

Novel Algorithm for Calculating Ideal Gas Phase Reactions Based on the Minimization of Free Energy
ZHANG Shouxin, ZHANG Zubin
College of Pipeline and Civil Engineering, China University of Petroleum, Qingdao, Shandong 266580, China
Abstract: In this work, we propose a gradient projected Lagrangian algorithm to address the inadequacies of current free-energy minimization-based solution method for calculating the composition of the ideal gas phase equilibrium. In this algorithm, a weak convergence criterion was used as the convergence criterion of the gradient projection method. Mathematical proofs were then provided for the rationality of the weak convergence conditions used in the gradient projection method, as well as the convergence of this algorithm. On this basis, the result calculated by the gradient projection method was used as the initial value in Newton's method with Lagrange multipliers, and improvements were also made on the iteration step used in Newton's method. We have thus resolved the difficulty of selecting an initial value for Newton's method, and improved the stability and computational speed of the algorithm. The results of example calculations demonstrated that the algorithm proposed in this work converges very quickly and exhibits a high level of computational accuracy.
Key words: free-energy minimization methods     ideal gas phase reactions     Claus reaction     Lagrange multiplier     gradient projection method

1 最小自由能法

 $G = \mathop \sum \limits_{i = 1}^K n_i \mu _i$ (1)

$G$ —Gibbs自由能，J；

$n_i$ —组分 $i$ 的物质的量，mol；

$\mu _i$ —组分 $i$ 的化学势，J·mol-1

 ${\mu _i} = \mu _i^\ominus + {\rm{R}}T\ln \frac{{{p_i}}}{{{p^\ominus }}}$ (2)

$\mu _i^ \ominus$ —组分 $i$ 标态化学势，J·mol-1

R—理想气体常数，R=8.314 J·(mol·K)-1

$T$ —体系温度，K；

$p_i$ —组分 $i$ 的分压，kPa；

$p^ \ominus$ —标准态的压力， $p^ \ominus$ =100 kPa。

 $\mu _i^ \ominus {\rm{ = }}\Delta _{\rm{f}} G_i^ \ominus$ (3)

$\Delta _{\rm{f}} G_i^ \ominus$ —单质生成化合物的标准摩尔Gibbs自由能的变化量，J·mol-1

 $\mu _i = \Delta _{\rm{f}} G_i^ \ominus + {\rm{R}}T\ln {\dfrac{{\gamma _i py_i }}{{p^ \ominus }}}$ (4)

$\gamma _i$ —组分 $i$ 的逸度系数；

$y_i$ —组分 $i$ 的摩尔分数；

$p$ —体系压力，kPa。

 $\mu _i = \Delta _{\rm{f}} G_i^ \ominus + {\rm{R}}T\ln{\dfrac{{\gamma _i py_i }}{{p^ \ominus }}} = \Delta _{\rm{f}} G_i^ \ominus + {\rm{R}}T\ln {\dfrac{{py_i }}{{p^ \ominus }}} =\\\hspace{3em} \Delta _{\rm{f}} G_i^ \ominus + {\rm{R}}T\ln \left( {\dfrac{p}{{p^ \ominus }}\dfrac{{n_i }}{{\sum\limits_{i = 1}^K {n_i } }}} \right)$ (5)

 $\left\{ \begin{array}{l} \min {\rm{ }}G = \mathop \sum \limits_{i = 1}^K n_i \left[{\Delta _{\rm{f}} G_i^ \ominus + {\rm{R}}T\ln \left( {\dfrac{p}{{p^ \ominus }}\dfrac{{n_i }}{{\sum\limits_{i = 1}^K {n_i } }}} \right)_i } \right] \\ {\rm{s.t.}}~~~~~n_i \geqslant 0~~~~~i = 1, 2, 3, \cdots, K \\ {\rm{ }}\sum\limits_{i = 1}^K {\beta _{ji} n_i = \omega _j }~~~~~~j = 1, {\kern 1pt} 2{\kern 1pt}, {\kern 1pt} 3{\kern 1pt}, \cdots, {\kern 1pt} W \\ \end{array} \right.$ (6)

 $Q = {H_{\rm{P}}} - {H_{\rm{R}}} = \sum\limits_{i = 1}^K {{n_i}{H_i}\left( T \right)} - {H_{\rm{R}}}$ (7)

$Q$ —体系与外界的热交换值，J；

$H_{\rm{P}}$ —生成物体系的焓值，J；

$H_{\rm{R}}$ —反应物体系的焓值，J；

$H_i$ —组分 $i$ 的摩尔焓，J·mol-1

2 梯度投影-拉格朗日算法

2.1 梯度投影法

 $\left\{ {\begin{array}{*{20}c} {\min f\left( \boldsymbol{{x}} \right)} \\ {\rm{s.t.}~~~~~\boldsymbol{{A}}\boldsymbol{{x}} \leqslant \boldsymbol{{b}}} \\ {\boldsymbol{{E}}\boldsymbol{{x}} = \boldsymbol{{c}}} \\ \end{array}} \right.$ (8)

$\boldsymbol{{n}}$ 为组分摩尔数向量，对比问题(6)、(8)的形式可知，采用梯度投影法求解问题(6)时，有 $\boldsymbol{{A}}=-\boldsymbol{{I}}$ ，其中 $\boldsymbol{{I}}$ 为单位阵； $\boldsymbol{{b}}=0$ $\boldsymbol{{E}}=\mathit{\pmb{β}}$ $\boldsymbol{{c}}={{\textbf{ω}}}$ ，所以最终需要求解的问题为

 $\left\{ {\begin{array}{*{20}{c}} {\min G\left( \boldsymbol{{n}} \right)}\\ {\rm{s.t.}~~~~~ - \boldsymbol{{I}}\boldsymbol{{n}} \leqslant 0}\\ {\mathit{\pmb{β}} \boldsymbol{{n}} = {\textbf{ω}} } \end{array}} \right.$ (9)
2.2 拉格朗日乘数法

 ${S_i}(\boldsymbol{{X}}) = 0~~~~~~i = 1, 2, \cdots, {\kern 1pt} K + W + 2$ (10)

 $\left\{ \begin{array}{l} {n_{K + 1}}\sum\limits_{i = 1}^K {\left[{\left( {{{\rm{e}}^{{x_i}}}} \right){\beta _{ji}}} \right]} = {\omega _j}~~~~~j = 1, 2, 3, \cdots W\\ {\Delta _{\rm{f}}}G_i^ \ominus + {\rm{R}}T\ln {\dfrac{{{{\rm{e}}^{{x_i}}}p}}{{{p^ \ominus }}}} + \sum\limits_{j = 1}^W {{\lambda _j}} {\beta _{ji}} = 0{\rm{ }}~~~~~i = 1, 2, 3, \cdots {\kern 1pt}, K\\ \sum\limits_{i = 1}^K {{{\rm{e}}^{{x_i}}} = 1} \\ {n_{K + 1}}\sum\limits_{i = 1}^K {{{\rm{e}}^{{x_i}}}{H_i}\left( T \right)} - {H_{\rm{R}}} - Q = 0 \end{array} \right.{\rm{ }}$ (11)

${x_i} = \ln \dfrac{{{n_i}}}{{{n_{K + 1}}}}$ ，其中 $i = 1, 2, 3, \cdots, {\kern 1pt} K$

$\mathit{\pmb{λ}}$ —拉格朗日乘子。

2.3 算法改进 2.3.1 梯度投影法的改进

$\boldsymbol{{x}}$ 是式(8)的一个可行解，且使 ${\boldsymbol{{A}}_1}\boldsymbol{{x}} = {\boldsymbol{{b}}_1}$ ${\boldsymbol{{A}}_2}\boldsymbol{{x}} < {\boldsymbol{{b}}_2}$ ，而 ${\boldsymbol{{A}}^{\rm{T}}} = ({\boldsymbol{{A}}_1}^{\rm{T}}$ ${\boldsymbol{{A}}_2}^{\rm{T}})$ ${\boldsymbol{{b}}^{\rm{T}}} = ({\boldsymbol{{b}}_1}^{\rm{T}}$ ${\boldsymbol{{b}}_2}^{\rm{T}})$ ${\boldsymbol{{M}}^{\rm{T}}} = ({\boldsymbol{{A}}_1}^{\rm{T}}$ ${\boldsymbol{{E}}^{\rm{T}}})$ 满秩，令

 $\boldsymbol{{P}} = \boldsymbol{{I}} - {\boldsymbol{{M}}^{\rm{T}}}{(\boldsymbol{{M}}{\boldsymbol{{M}}^{\rm{T}}})^{ - 1}}\boldsymbol{{M}}$ (12)
 $\boldsymbol{{w}} = - {(\boldsymbol{{M}}{\boldsymbol{{M}}^{\rm{T}}})^{ - 1}}\boldsymbol{{M}}\nabla f\left( \boldsymbol{{x}} \right){\boldsymbol{{w}}^{\rm{T}}} = ({\boldsymbol{{u}}^{\rm{T}}}, {\boldsymbol{{v}}^{\rm{T}}})$ (13)

$\boldsymbol{{P}}\nabla f\left( \boldsymbol{{x}} \right) = 0$

(1) 若 $\forall {u_i} \geqslant 0$ ，记为 $\boldsymbol{{u}}\geqslant 0$ ，则 $\boldsymbol{{x}}$ 是一个K-T(Kuhn-Tucker)点，即 $\boldsymbol{{x}}$ 为约束极值点；

(2) 若 $\exists {u_i} < 0$ ，令 ${\hat {\boldsymbol{{M}}}^{\rm{T}}} = ({\hat{ \boldsymbol{{A}}}_1}^{\rm{T}}, {\boldsymbol{{E}}^{\rm{T}}})$ ，其中 ${{{\hat {\boldsymbol{{A}}}}}_1}$ 是由 ${\boldsymbol{{A}}_1}$ 中去掉第 $i$ 行后得到的矩阵，令

 $\hat {\boldsymbol{{P}}} = \boldsymbol{{I}} - \hat{ \boldsymbol{{M}}}^{\rm{T}}{(\hat {\boldsymbol{{M}}}\hat{ \boldsymbol{{M}}}^{\rm{T}})^{ - 1}}\hat {\boldsymbol{{M}}}$ (14)
 $\boldsymbol{{d}} = - \hat {\boldsymbol{{P}}}\nabla f\left( \boldsymbol{{x}} \right)$ (15)

$\boldsymbol{{d}}$ 是梯度投影法一个改进的可行方向[20]

$\boldsymbol{{P}}\nabla f\left( \boldsymbol{{x}} \right) = {\mathit{\pmb{ε}}_0}$ ，因为

 $\boldsymbol{{P}}\nabla f\left( \boldsymbol{{x}} \right) = [\boldsymbol{{I}}-{\boldsymbol{{M}}^{\rm{T}}}{(\boldsymbol{{M}}{\boldsymbol{{M}}^{\rm{T}}})^{-1}}\boldsymbol{{M}}]\nabla f\left( \boldsymbol{{x}} \right) = \nabla f\left( \boldsymbol{{x}} \right) +\\\hspace{2em} {\boldsymbol{{M}}^{\rm{T}}}\boldsymbol{{w}} =\nabla f\left( \boldsymbol{{x}} \right) + \boldsymbol{{A}}_1^{\rm{T}}\boldsymbol{{u}} + {\boldsymbol{{E}}^{\rm{T}}}\boldsymbol{{v}}$ (16)

 $\nabla f\left( \boldsymbol{{x}} \right) + \boldsymbol{{A}}_1^{\rm{T}}\boldsymbol{{u}} + {\boldsymbol{{E}}^{\rm{T}}}\boldsymbol{{v}} = {\mathit{\pmb{ε}} _0}$ (17)
 $\nabla [f\left( \boldsymbol{{x}} \right)-{\mathit{\pmb{ε}} _0}\boldsymbol{{x}}] + \boldsymbol{{A}}_1^{\rm{T}}\boldsymbol{{u}} + {\boldsymbol{{E}}^{\rm{T}}}\boldsymbol{{v}} = 0$ (18)

$\bar F(\boldsymbol{{x}}) = f\left( \boldsymbol{{x}} \right) - {\mathit{\pmb{ε}} _0}\boldsymbol{{x}}$ ，则有

 $\nabla \bar F\left( \boldsymbol{{x}} \right) + \boldsymbol{{A}}_1^{\rm{T}}\boldsymbol{{u}} + {\boldsymbol{{E}}^{\rm{T}}}\boldsymbol{{v}} = 0$ (19)
 $\boldsymbol{{P}}\nabla \bar F\left( \boldsymbol{{x}} \right) = 0$ (20)

 ${\boldsymbol{{w}}_0} = - {(\boldsymbol{{M}}{\boldsymbol{{M}}^{\rm{T}}})^{ - 1}}\boldsymbol{{M}}\nabla \bar F\left( \boldsymbol{{x}} \right)$ (21)
 ${\boldsymbol{{w}}_0}^{\rm{T}} = ({\boldsymbol{{u}}_0}^{\rm{T}}, {\boldsymbol{{v}}_0}^{\rm{T}})$ (22)

$\bar F(\boldsymbol{{x}}) = f\left( \boldsymbol{{x}} \right) - {\mathit{\pmb{ε}}_0}\boldsymbol{{x}}$ $\boldsymbol{{P}}\nabla f\left( \boldsymbol{{x}} \right) = {\mathit{\pmb{ε}}_0}$ 代入式(21)中可得 ${\boldsymbol{{w}}_0} = \boldsymbol{{w}}$ ，即 ${{\boldsymbol{{{u}}}}_0} = {\boldsymbol{{{u}}}}$ ${{\boldsymbol{{{v}}}}_0} = {\boldsymbol{{{v}}}}$

 $\left\{ {\begin{array}{*{20}{c}} {\min {\rm{ }}\bar F(\boldsymbol{{x}}) = f\left( \boldsymbol{{x}} \right) - {\mathit{\pmb{ε}}_0}\boldsymbol{{x}}}\\ {\rm{s.t.}~~~~~\boldsymbol{{A}}\boldsymbol{{x}} \leqslant \boldsymbol{{b}}{\rm{ }}}\\ {\boldsymbol{{E}}\boldsymbol{{x}} = \boldsymbol{{c}}} \end{array}} \right.$ (23)

 $f\left( {{\boldsymbol{{x}}_1}} \right) - f({\boldsymbol{{x}}_2}) > 0$ (24)

 $\bar F({\boldsymbol{{x}}_2}) - \bar F\left( {{\boldsymbol{{x}}_1}} \right) > 0$ (25)
 $f({\boldsymbol{{x}}_2}) - {\mathit{\pmb{ε}}_0}{\boldsymbol{{x}}_2} - f({\boldsymbol{{x}}_1}) + {\mathit{\pmb{ε}}_0}{\boldsymbol{{x}}_1} > 0$ (26)
 $f({\boldsymbol{{x}}_1}) - f({\boldsymbol{{x}}_2}) < {\mathit{\pmb{ε}}_0}({\boldsymbol{{x}}_1} - {\boldsymbol{{x}}_2})$ (27)

 $0 < f\left( {{\boldsymbol{{x}}_1}} \right) - f({\boldsymbol{{x}}_2}) < {\mathit{\pmb{ε}}_0}({\boldsymbol{{x}}_1} - {\boldsymbol{{x}}_2})$ (28)

2.3.2 牛顿法的改进

 ${\boldsymbol{{X}}^{\left( {k + 1} \right)}} = {\boldsymbol{{X}}^{\left( k \right)}} - {\delta ^{\left( {k + 1} \right)}}\nabla S{\left( {{\boldsymbol{{X}}^{\left( k \right)}}} \right)^{ - 1}}S\left( {{\boldsymbol{{X}}^{\left( k \right)}}} \right)$ (29)

2.4 算法描述

 $\boldsymbol{{A}} = \left[{\begin{array}{*{20}{c}} {\bar {\boldsymbol{{a}}}_1^{\rm{T}}}\\ {\begin{array}{*{20}{c}} {\bar {\boldsymbol{{a}}}_2^{\rm{T}}}\\ \vdots \end{array}}\\ {\bar {{\boldsymbol{{a}}}}_K^{\rm{T}}} \end{array}} \right], E = \left[{\begin{array}{*{20}{c}} {\bar {\boldsymbol{{a}}}_{K + 1}^{\rm{T}}}\\ {\begin{array}{*{20}{c}} {\bar {\boldsymbol{{a}}}_{K + 2}^{\rm{T}}}\\ \vdots \end{array}}\\ {\bar {\boldsymbol{{a}}}_{K + W}^{\rm{T}}} \end{array}} \right]$ (30)

 ${\bar {\boldsymbol{{a}}}_j} = {\left( {{{\bar {a}}_{j1}}{\kern 1pt} {{\bar {a}}_{j2}}{\kern 1pt} \cdots {\kern 1pt} {{\bar {a}}_{jK}}} \right)^{\rm{T}}}j = 1, 2, \cdots, {\kern 1pt} K + W$ (31)

 $\boldsymbol{{B}} = \left( {\begin{array}{*{20}{c}} \boldsymbol{{b}}\\ \boldsymbol{{c}} \end{array}} \right) = \left[\begin{array}{l} {B_1}\\ {B_2}\\ \vdots \\ {B_{K + W}} \end{array} \right]$
 ${\boldsymbol{{J}}_k} = J\left( {{x^{\left( k \right)}}} \right) = \left\{ {j{\rm{|}}\bar {\boldsymbol{{a}}}_j^{\rm{T}}{\boldsymbol{{x}}^{\left( k \right)}} = {B_j}} \right\}$ (32)

${\boldsymbol{{J}}_k}$ 为指标集。用 ${\boldsymbol{{M}}_k}$ 表示以 $\bar{\boldsymbol{{a}}}_j^{\rm{T}}j \in {\boldsymbol{{J}}_k}$ 为行所组成的矩阵。

(1) 给定体系的初始温度 $T_1$ ；给定计算精度 ${\varepsilon _1} > 0$ ，梯度投影弱收敛标准 ${\varepsilon _2} > 0$ ；给定参数 $h>0$ $0 < {\delta ^{(1)}} < 1$ $1 < \varphi < 2$ ；给定计算精度 ${\varepsilon _3} > 0$ ，松弛系数 $0 < \alpha \leqslant 1$ ，步长 $\Delta T > 0$ ，令 $i=1$ $k=1$ ；任意给定 ${\boldsymbol{{n}}^{\left( 0 \right)}} > 0$ 作为初值，采用梯度投影法求解

 $\left\{ \begin{array}{l} \min {\rm{ }}\sum\limits_j^W {{{\left[{\mathop \sum \limits_{i = 1}^K \left( {{\beta _{ji}}{n_i}-{\omega _j}} \right) } \right]}^2} } \\ {\rm{s.t.}}~~~~~{n_i} \geqslant 0~~~~~i = 1, 2, 3, {\kern 1pt} \ldots, {\kern 1pt} K \end{array} \right.$ (33)

(2) 令 $T=T_i$ ，计算 $\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right)$ ${\boldsymbol{{J}}_k} = \left\{ {j{\rm{|}}\bar{ \boldsymbol{{a}}}_j^{\rm{T}}{\boldsymbol{{x}}^{\left( k \right)}} = {B_j}} \right\}$ ，若 ${\left\| {\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right)} \right\|_\infty } < {\varepsilon _1}$ ，则 ${\boldsymbol{{n}}^{\left( k \right)}}$ 为近似的 ${\rm{K}} - {\rm{T}}$ 点，转步骤(5)；否则，若 ${\boldsymbol{{J}}_k} = \emptyset$ (空集)，令 ${\boldsymbol{P}} = \boldsymbol{{I}}$ ；若 ${\boldsymbol{{J}}_k} \ne \emptyset$ (空集)，令 $\boldsymbol{{P}} = \boldsymbol{{I}} - \boldsymbol{{M}}_k^{\rm{T}}{({\boldsymbol{{M}}_k}\boldsymbol{{M}}_k^{\rm{T}})^{ - 1}}{\boldsymbol{{M}}_k}$

(3) 若 ${\left\| {\boldsymbol{{P}}\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right)} \right\|_\infty } > {\varepsilon _2}$ ，令 ${\boldsymbol{{d}}_k} = - \boldsymbol{{P}}\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right)$ ，转步骤(4)；若 ${\left\| {\boldsymbol{{P}}\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right)} \right\|_\infty } \leqslant {\varepsilon _2}$ ，令

 $\boldsymbol{{w}} = - {({\boldsymbol{{M}}_k}\boldsymbol{{M}}_k^{\rm{T}})^{ - 1}}{\boldsymbol{{M}}_k}\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right) = \left( {\begin{array}{*{20}{c}} \boldsymbol{{u}}\\ \boldsymbol{{v}} \end{array}} \right)$ (34)

$\boldsymbol{{u}}\geqslant 0$ ，则 ${\boldsymbol{{n}}^{\left( k \right)}}$ ${\rm{K}} - {\rm{T}}$ 点，转步骤(5)；若 $\boldsymbol{{u}}$ 有某个分量 ${u_j} < 0$ ，令 ${\hat {\boldsymbol{{M}}}_k}$ 是在 ${\boldsymbol{{M}}_k}$ 中去掉与 ${u_j}$ 对应的第 $j$ 行而得到的矩阵。令 $\hat {\boldsymbol{{P}}} = {\bf{\boldsymbol{{I}}}} - \hat {\boldsymbol{{M}}}_k^{\rm{T}}{({\hat { \boldsymbol{{M}}}_k}\hat {\boldsymbol{{M}}}_k^{\rm{T}})^{ - 1}}{\hat { \boldsymbol{{M}}}_k}$ ${\boldsymbol{{d}}_k} = - \hat{ \boldsymbol{{P}}}\nabla f\left( {{\boldsymbol{{n}}^{\left( k \right)}}} \right)$ ，转步骤(4)。

(4) 计算 ${\lambda _{\max }} = \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\mathop {\min }\limits_{} \left\{ {\frac{{{B_i} - \bar{ \boldsymbol{{a}}}_i^{\rm{T}}{\boldsymbol{{n}}^{\left( k \right)}}}}{{\bar{ \boldsymbol{{a}}}_i^{\rm{T}}{\boldsymbol{{d}}_k}}}{\rm{|}}i\bar {\in} {J_k}, \bar{ \boldsymbol{{a}}}_i^{\rm{T}}{\boldsymbol{{d}}_k} > 0} \right\}}\\ \end{array}\\ + \infty {\rm{}}(i, \bar {\boldsymbol{{a}}}_i^{\rm{T}}{\boldsymbol{{d}}_k} \leqslant 0) \end{array} \right.$

${\lambda _k}$ 是线搜索问题

 $\left\{ {\begin{array}{*{20}{c}} {\mathop {\min }\limits_{} f\left( {{\boldsymbol{{n}}^{\left( k \right)}} + \lambda {\boldsymbol{{d}}_k}} \right)}\\ {\rm{s.t.}~~~~~{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0 \leqslant \lambda \leqslant {\lambda _{\rm{MAX}}}} \end{array}} \right.$ (35)

(5) 计算 $g({T_i})$ $g({T_i} + \Delta T)$ ，其中 $g(T) = {H_{\rm{P}}}(T) - {H_{\rm{R}}} - Q$ ，则

 ${T_{i + 1}} = {T_i} - \alpha \dfrac{{g({T_i})\Delta T}}{{g({T_i} + \Delta T) - g({T_i})}}$ (36)

(6) 如果 $\left| {{T_i}_{ + 1} - {T_i}} \right| < {\varepsilon _3}$ ，令 $T= T_i+1$ $x_j^{(1)} = \ln \dfrac{{n_j^{(k)}}}{{{n_{K + 1}}}}, j = 1, 2, \cdots, {\kern 1pt} K$ ，转到步骤(7)，否则，令 $T_i=T_i+1$ $i=i+1$ ，转到步骤(2)。

(7) 令 ${\boldsymbol{{X}}^{(1)}}\!=\!{\left( {x_1^{(1)}, {\kern 1pt} x_2^{(1)}, {\kern 1pt} \cdots, {\kern 1pt} x_K^{(1)}, {\lambda _1}, {\kern 1pt} {\lambda _2}, {\kern 1pt} \cdots, {\kern 1pt} {\lambda _W}, {n_{K + 1}}, T} \right)^{\rm{T}}}$ ，其中 ${\lambda _i} > 0，i = 1, 2, \cdots {\kern 1pt}, W$

(8) 令 $k$ =1，计算 ${\bar B_i} = {S_i}\left( {{\boldsymbol{{X}}^{(k)}}} \right)$ $i = 1, 2, \cdots, {\kern 1pt} K + W + 2$ ，若 $\mathop {\max }\limits_{1 \leqslant i \leqslant K + W + 2} \left| {{{\bar B}_i}} \right| < {\varepsilon _1}$ ，则 ${\boldsymbol{{X}}^{(k)}}$ 为方程组的一组解，计算结束， $\boldsymbol{{n}} = X_{K + W + 1}^{(k)}{\left( {{{\rm{e}}^{X_1^{(k)}}}, {\kern 1pt} {{\rm{e}}^{X_2^{(k)}}}, {\kern 1pt} \cdots, {\kern 1pt} {{\rm{e}}^{X_k^{(k)}}}} \right)^{\rm{T}}}$ $T = X_{K + W + 2}^{(k)}$ ；若 $\mathop {\max }\limits_{1 \leqslant i \leqslant K + W + 2} \left| {{{\bar B}_i}} \right| \geqslant {\varepsilon _1}$ ，计算

 $\bar {\boldsymbol{{A}}} = \left[{\begin{array}{*{20}{c}} {{S_1}\left( {\boldsymbol{{X}}_1^{(k)}} \right)}\\ {\begin{array}{*{20}{c}} {{S_2}\left( {\boldsymbol{{X}}_1^{(k)}} \right)}\\ \vdots \end{array}}\\ {{S_{K + W + 2}}\left( {\boldsymbol{{X}}_1^{(k)}} \right)} \end{array}\begin{array}{*{20}{c}} {{S_1}\left( {\boldsymbol{{X}}_2^{(k)}} \right)}\\ {\begin{array}{*{20}{c}} {{S_2}\left( {\boldsymbol{{X}}_2^{(k)}} \right)}\\ \vdots \end{array}}\\ {{S_{K + W + 2}}\left( {\boldsymbol{{X}}_2^{(k)}} \right)} \end{array}\begin{array}{*{20}{c}} \cdots \\ {\begin{array}{*{20}{c}} \cdots \\ \ddots \end{array}}\\ \ldots \end{array}\begin{array}{*{20}{c}} {{S_1}\left( {\boldsymbol{{X}}_{K + W + 2}^{(k)}} \right)}\\ {\begin{array}{*{20}{c}} {{S_2}\left( {\boldsymbol{{X}}_{K + W + 2}^{(k)}} \right)}\\ \vdots \end{array}}\\ {{S_{K + W + 2}}\left( {\boldsymbol{{X}}_{K + W + 2}^{(k)}} \right)} \end{array}} \right]$ (37)

(9) 解线性方程组 $\bar {\boldsymbol{{A}}}\boldsymbol{{Z}} = \bar {\boldsymbol{{B}}}$ ，其中 $\boldsymbol{{Z}} = {\left( {{z_1}, {\kern 1pt} {z_2}, {\kern 1pt}, \cdots, {\kern 1pt} {z_{K + W + 2}}} \right)^{\rm{T}}}$ ，并计算 $\eta = 1 - \mathop \sum \limits_{i = 1}^{K + W + 2} {z_i}$ ，若 $\eta=0$ ，令 ${\varepsilon _2} = \dfrac{1}{2}{\varepsilon _2}$ ，令 $i=1$ $k=1$ ，转(2)；若 $\eta \ne 0$ ，计算 $X_i^{(k + 1)} = X_i^{(k)} - {\delta ^{(k)}}h{z_i}/\eta {\rm{ }}{\rm{ }}~~~~i = 1, {\kern 1pt} 2, {\kern 1pt} \cdots, {\kern 1pt} K + W + 2$ ${\delta ^{(k + 1)}} = \varphi {\delta ^{(k)}}$ $k=k+1$ ，返回(8)。

3 算例

 图1 参数 $\varepsilon_2$ 对计算时间的影响 Fig. 1 Influence of parameters $\varepsilon_2$ on calculation time
 图2 参数 ${\delta ^{(1)}}$ 对计算时间的影响 Fig. 2 Influence of parameters ${\delta ^{(1)}}$ on calculation time
 图3 参数 $\varphi$ 对计算时间的影响 Fig. 3 Influence of parameters $\varphi$ on calculation time
 图4 温度初值 $T$ 对计算时间的影响 Fig. 4 Influence of initial temperature $T$ on calculation time

4 结论

(1) 将梯度投影法的收敛标准 $\boldsymbol{{P}}\nabla f\left( {{\boldsymbol{{x}}^{\left( k \right)}}} \right) = 0$ 松弛为 ${\left\| {\boldsymbol{{P}}\nabla f\left( {{\boldsymbol{{x}}^{\left( k \right)}}} \right)} \right\|_\infty } < \varepsilon {\rm{ }}$ 时， $\exists {\mathit{\pmb{ε}}_0}$ 满足 ${\left\| {{\mathit{\pmb{ε}}_0}} \right\|_\infty } < \varepsilon$ ，使得式(27)成立，即当 $\varepsilon$ 为适当小的非负实数时，由梯度投影法计算出的体系的状态可以近似为体系的平衡状态，所以可作为牛顿法的初始计算状态。

(2) 本文提出了梯度投影-拉格朗日算法，避免了人为选取计算初值，从算例的计算的结果中可以看出该算法的计算速度比拉格朗日乘数法和梯度投影法快，并且计算精度与拉格朗日乘数法相同，在各关键参数大范围变化的情况下，算法的收敛性稳定，说明该算法是一种计算基于最小自由能法的理想气相反应平衡组成的有效方法。

(3) 为了保证算法的收敛性，应控制关键参数 $\varepsilon_2$ ${\delta ^{(1)}}$ $\varphi$ 的大小，根据计算经验推荐选取 $\varepsilon_2\leqslant 10$ ${\delta ^{(1)}}$ =0.1， $\varphi=1.1$ ，并且合理估算平衡温度有利于提高算法的收敛速度。

 [1] 王武谦, 方晨昭, 韩方煜. 多相多组份化学平衡模拟的研究(自由能最小法)[J]. 计算机与应用化学, 1990, 7(4): 259–267. WANG Wuqian, FANG Chenzhao, HAN Fanyu. Study on the equlibrium simulation of multicomponent and mutiphase system (mimizating gibbs free energy)[J]. Computers and Applied Chemistry, 1990, 7(4): 259–267. doi: 10.16866/j.com.app.chem1990.04.004 [2] 李国栋. 基于Gibbs自由能最小法的反应过程优化设计的研究[D]. 青岛: 中国海洋大学, 2007. LI Goudong. Study on optimal design of reaction process based on Gibbs free energy minimization method[D]. Qingdao: Ocean University of China, 2007. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1070593 [3] 朱利凯, 鲍钧. 用最小自由能法计算克劳斯反应的平衡组成[J]. 石油学报(石油加工), 1990, 6(2): 75–82. ZHU Likia, BAO Jun. Calculation of equilibrium composition of CLAUS reactions by free energy minimization technique[J]. Acta Petrolei Sinca (Petroleum Processing Section), 1990, 6(2): 75–82. [4] 杜军驻, 仇汝臣, 冀刚, 等. 克劳斯硫磺回收工艺数学模型的研究[J]. 上海化工, 2013, 38(3): 7–12. DU Junzhu, QIU Ruchen, JI Gang, et al. Study on mathematical model of Claus sulfur recovery process[J]. Shanghai Chemical Industry, 2013, 38(3): 7–12. doi: 10.3969/j.-issn.1004-017X.2013.03.003 [5] WHITE W B, JOHNSON S M, DANTZIG G B. Chemical equilibrium in complex mixtures[J]. Journal of Chemical Physics, 1958, 28(5): 751–755. doi: 10.1063/1.1744264 [6] GAUTAM R, SEIDER W D. Computation of phase and chemical equilibrium:Part Ⅰ. Local and constrained minima in Gibbs free energy[J]. Aiche Journal, 2004, 25(6): 991–999. doi: 10.1002/aic.690250610 [7] LUCIA A, PADMANABHAN L, VENKATARAMAN S. Multiphase equilibrium flash calculations[J]. Computers and Chemical Engineering, 2000, 24(12): 2557–2569. doi: 10.1016/S0098-1354(00)00563-9 [8] 郭汉杰, 赵玉祥. 最小自由能原理的SUMT方法[J]. 北京科技大学学报, 1992, 14(5): 502–508. GUO Hanjie, ZHAO Yuxiang. Method of minimizing the global free energy with SUMT program[J]. Journal of University of Science and Technology Beijing, 1992, 14(5): 502–508. doi: 10.13374/j.issn1001-053x.1992.05.002 [9] 许小平, 张唯, 李欣. 惩罚函数法计算燃烧产物的平衡组分[J]. 宇航学报, 1994, 15(3): 90–95. XU Xiaoping, ZHANG Wei, LI Xin. Calculation of equilibrium composition of combustion products by penalty method[J]. Journal of Astronautics, 1994, 15(3): 90–95. [10] NÉRON A, LANTAGNE G, MARCOS B. Computation of complex and constrained equilibria by minimization of the Gibbs free energy[J]. Chemical Engineering Science, 2012, 82: 260–271. doi: 10.1016/j.ces.2012.07.041 [11] GORDON S. Computer program for calculation of complex chemical equilibrium compositions, rocket performance incident and reflected shocks and ChapmanJouguet detonation[C]. NASA SP-273, 1971. [12] 金汀. 用最小自由能法计算克劳斯工艺过程[J]. 天然气工业, 1992, 12(2): 66–70. JIN Ding. Calculating the CLAUS technology process by use of minimum free energy method[J]. Natural Gas Industry, 1992, 12(2): 66–70. [13] 徐金火, 汤渭龙, 沈复. 克劳斯硫磺回收流程的工艺模拟[J]. 石油大学学报(自然科学版), 1993, 17(5): 93–99. XU Jinhuo, TANG Weilong, SHEN Fu. Flowsheet simulation for CLAUS recovery process of sulfur[J]. Journal of the University of Petroleum, China, 1993, 17(5): 93–99. [14] 黄河, 李政, 倪维斗. Gibbs反应器模型及其计算机实现[J]. 动力工程, 2004, 24(6): 902–907. HUANG He, LI Zheng, NI Weidou. The Gibbs reactor model and its realization on the computer[J]. Power Engineering, 2004, 24(6): 902–907. doi: 10.3321/j.issn:1000-6761.2004.06.031 [15] 高磊, 翟悦, 聂颖, 等. 自由能最小化方法处理化学激光器燃烧产物组份的探讨[J]. 舰船防化, 2009(1): 43–47. GAO Lei, ZHAI Yue, NIE Ying, et al. Discussion on method of minimum Gibbs'free energy dealing with compositions of combustion products in chemical laser[J]. Chemical Defence on Ships, 2009(1): 43–47. [16] GEORGE B, BROWN L P, FARMER C H, et al. Computation of multicomponent, multiphase equilibrium[J]. Industrial & Engineering Chemistry Process Design & Development, 1976, 15(3): 163–169. doi: 10.-1021/i260059a003 [17] 周维彪, 许志宏. 多相多组元化学平衡相平衡计算(Ⅰ)——算法M-SVMP[J]. 化工学报, 1987, 38(1): 39–48. ZHOU Weibiao, XU Zhihong. Calculation of multicomponent multiphase equilibria (Ⅰ)-Modified algorithm MSVMP method[J]. Journal of Chemeical Industry and Engineering (China), 1987, 38(1): 39–48. [18] 周维彪, 许志宏. 多相多组元化学平衡相平衡计算(Ⅱ)——新算法GCG法[J]. 化工学报, 1987, 38(1): 49–55. ZHOU Weibiao, XU Zhihong. Calculation of multicomponent multiphase equilibria (Ⅱ)-A new general method GCG[J]. Journal of Chemeical Industry and Engineering (China), 1987, 38(1): 49–55. [19] SILVA A L D, MALFATTI C D F, MÜLLER I L. Thermodynamic analysis of ethanol steam reforming using Gibbs energy minimization method:A detailed study of the conditions of carbon deposition[J]. International Journal of Hydrogen Energy, 2009, 34(10): 4321–4330. doi: 10.1016/j.ijhydene.2009.03.029 [20] 唐焕文, 秦学志. 实用最优化算法[M]. 大连: 大连理工大学出版社, 2007: 183-184. [21] 刘停战, 刘伟, 何颖. 调整步长牛顿法[J]. 中国传媒大学学报(自然科学版), 2012, 19(1): 8–10. LIU Tingzhan, LIU Wei, HE Ying. Step-adjusting newton method[J]. Journal of Communication University of China Science and Technology, 2012, 19(1): 8–10. doi: 10.3969/-j.issn.1673-4793.2012.01.002