引言

理想气相反应平衡计算有很广泛的用途，例如在炼厂、含硫天然气净化厂等石化企业中广泛应用克劳斯硫磺回收工艺中的克劳斯反应，便可以视为复杂均相的理想气相反应。最小自由能法是计算该问题的一种重要的方法，该方法将平衡问题转化为非线性约束优化问题^[1-2]，具有如下优点：反应产物组分的数量不受限制^[3]；计算结果与复杂的中间化学反应过程无关；结果优于平衡常数法^[4]。

自20世纪50年代起，国内外学者对最小自由法及其求解方法进行了研究。该问题的求解方法主要可以分为3类：(1)用线性规划或二次规划来逐次逼近非线性优化问题的方法，如分段线性化方法^[5]、序列二次规划法^[6-7]；(2)将约束优化问题转化为无约束优化问题的方法，如惩罚函数法^[8-10]、拉格朗日乘数法^{[3, 11-15]}、Powell法^[16]；(3)对约束优化问题直接进行处理的方法，如变尺度投影法^[17]、共轭梯度投影法^[18]、梯度投影法^[1]、广义简约梯度法^[19]。但是由于该优化问题的强非线性，目前的方法在求解该问题时均存在一定的不足，在初值选取、计算速度和计算精度方面仍存在改进空间。本文在对梯度投影法和拉格朗日乘数法分析的基础之上，将二者结合在一起，形成一种新算法：梯度投影-拉格朗日算法，该算法避免了人为选取计算初值，并能保证较快的计算速度和较好的计算精度。

1 最小自由能法

最小自由能法由White等人于1958年提出^[5]。热力学第二定律指出，体系在给定的温度和压力下，处于平衡状态时，其Gibbs自由能最小。

设体系的组分数为 $K$ ，元素种数为 $W$ ，则体系的Gibbs自由能为

$ 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。

组分 $i$ 的化学势的计算公式为

$ {\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$ 等于纯净物的摩尔Gibbs自由能 $G^ \ominus$ 。若规定单质在标态压力时的摩尔Gibbs自由能 $G^ \ominus=0$ ，则有

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

(3)

式中：

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

对于实际气体，用逸度 $f_i = \gamma _i py_i $ 替代组分 $i$ 的分压 $p_i$ ，则组分 $i$ 的化学势为

$ \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。

对于理想气体 $\gamma=1$ ，则有

$ \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)

当体系达到平衡时， $G$ 为最小值，并满足质量守恒和非负约束条件。设原子矩阵为 $\mathit{\pmb{β}}$ ， $\beta_{ji}$ 为组分 $i$ 中 $j$ 元素的原子数； ${{\textbf{ω}}}$ 为元素的物质的量总量向量， ${{\textbf{ω}}_j}$ 为第 $j$ 种元素的物质的量总量，s.t.表示优化问题的约束条件，则该化学平衡计算可转化为非线性最优化问题

$ \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)

要使上述最优化问题具有实际意义，还需要加入能量平衡方程^[14]

$ 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{{A}}$ 是 $K\times K$ 矩阵， $\boldsymbol{{b}} \in {\boldsymbol{{R}}^K}$ ， $\boldsymbol{{x}} \in {\boldsymbol{{R}}^K}$ ， $\boldsymbol{{E}}$ 是 $W\times K$ 矩阵， $\boldsymbol{{c}} \in {\boldsymbol{{R}}^W}$ ， $f:{\boldsymbol{{R}}^K} \to \boldsymbol{{R}}f \in {\boldsymbol{{C}}^1}$ 。

设 $\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 拉格朗日乘数法

通过拉格朗日乘数法可以将最小自由法对应的非线性约束优化问题(6)转化为如下非线性方程组的求解问题

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

(10)

该非线性方程组中的方程数为 $K+M+2$ ，未知数 $\boldsymbol{{X}}$ 为 $K$ 个 $x_i$ 、 $W$ 个 $\lambda_i$ 、 ${n_{K + 1}}$ 、 $T$ ，共有 $K+W+2$ 个，所以方程组有唯一解。该方程组的具体形式为

$ \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}$ 、 ${\lambda_i}$ 、 ${n_{K + 1}}$ 、 $T$ —待求解，其中 ${n_{K + 1}} = \sum\limits_{i = 1}^K {{n_i}} $ ；

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

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

该非线性方程组中的方程数和未知数均为 $K+W+2$ 个，因此有唯一解。

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{{u}}$ 对应于不等式约束， $\boldsymbol{{v}}$ 对应于等式约束。

设 $\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) \ne 0$ 且 $\boldsymbol{{u}}\geqslant 0$ 的情况进行以下分析。

设 $\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}}}}$ 。

在以 $\boldsymbol{{P}}\nabla f\left( \boldsymbol{{x}} \right) = {\mathit{\pmb{ε}}_0}$ 为收敛准则的条件下，梯度投影法结束计算时，有 $\boldsymbol{{P}}\nabla \bar F\left( \boldsymbol{{x}} \right) = 0$ 且 ${\boldsymbol{{u}}_0} = \boldsymbol{{u}} \geqslant 0$ ，因此由上述定理可知，此时得到的是

$ \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)

的最优解，设为 $\boldsymbol{{x}}_1$ ；设原问题的真实最优解为 $\boldsymbol{{x}}_2$ ，则

对于 $f\left( \boldsymbol{{x}} \right)$ ，有

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

(24)

对于 $\bar F(\boldsymbol{{x}})$ ，有

$ \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)

将梯度投影法的收敛标准 $\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$ ，使得式(28)成立，所以当 $\varepsilon$ 为适当小的非负实数时，由梯度投影法计算出的体系的状态可以近似为体系的平衡状态，所以可以作为牛顿法计算的初始状态。

2.3.2 牛顿法的改进

为了保证梯度投影法的计算速度，要求弱收敛准则 $\varepsilon$ 不能取的过小，所以此时以梯度投影法的计算结果为初值，采用牛顿法求解拉格朗日乘数法的过程中仍可能出现不收敛的情况。目前一般通过调整步长来降低牛顿法对初值的要求^[21]，但是由于调整后的步长始终小于1，牛顿法的收敛速度得不到保证。因为在弱收敛条件下梯度投影法计算出的平衡状态接近体系的真实平衡状态，所以可对牛顿法作如下改进

$ {\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)

式中， ${\delta ^{\left( {k + 1} \right)}} = \varphi {\delta ^{\left( k \right)}}$ ，且若 ${\delta ^{\left( {k + 1} \right)}} > 1$ ，则 ${\delta ^{\left( {k + 1} \right)}} = 1$ ； $0 < {\delta ^{\left( 1 \right)}} < 1$ ， $1 < \varphi < 2$ 。

2.4 算法描述

在式(8)中，设

$ \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)

得最优解 $\boldsymbol{{n}}$ ，令 ${\boldsymbol{{n}}^{\left( 1 \right)}} = \boldsymbol{{n}}$ 。

(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)

的解，令 ${\boldsymbol{{n}}^{\left( {k + 1} \right)}} = {n^{\left( k \right)}} + {\lambda _k}{\boldsymbol{{d}}_k}$ ， $k = k + 1$ ，返回(2)。

(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)

其中 $\boldsymbol{{X}}_j^{(k)}\!=\!{\left( {X_1^{(k)}, {\kern 1pt} X_2^{(k)}, {\kern 1pt} \cdots, {\kern 1pt} X_{j - 1}^{(k)}, {\kern 1pt} X_j^{(k)}\!+\!h{\kern 1pt} X_{j + 1}^{(k)}, {\kern 1pt} \cdots, {\kern 1pt} X_{K + W + 2}^{(k)}} \right)^{\rm{T}}}$ 。

(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 算例

本文以克劳斯反应为例对梯度投影-拉格朗日算法的有效性进行了验证，并将计算结果与梯度投影法和拉格朗日乘数法的计算结果进行了对比；同时分析了梯度投影-拉格朗日算法中关键参数和温度初值对计算时间的影响。

克劳斯反应酸气流量为100 kmol/h，空气流量207.30 kmol/h，入口气体的温度为40 ℃，系统总压151.2 kPa，不计热损失，设反应产物为SO₂、CO₂、S₂、COS、CS₂、CO、H₂，酸气组成如表 1所示。

梯度投影-拉格朗日算法中取 ${\varepsilon _1} = {10^{ - 5}}$ ， ${\varepsilon _2} = 10$ ， ${\varepsilon _3} = {10^{ - 5}}$ ， $\alpha=1$ ， $\Delta T = {10^{ - 3}}$ ， $h=0.01$ ， ${\delta ^{(1)}} = 0.5$ ， $\varphi=1$ ；梯度投影法和拉格朗日乘数法中相关参数与梯度投影-拉格朗日算法中选取的相同；拉格朗日乘数法采用拟牛顿法求解，经试算计算初值选为 ${x_i} = 1{\rm{ }}{\rm{ }}(i = 1, {\kern 1pt} 2, {\kern 1pt} \cdots, {\kern 1pt} K)$ ， ${\lambda _i} = 1({\rm{ }}{\rm{ }}i = 1, {\kern 1pt} 2, {\kern 1pt} \cdots, {\kern 1pt} W)$ ， ${n_{K + 1}}{\rm{ = }}1$ ， $T$ =1 428 K时，算法收敛，梯度投影-拉格朗日算法中温度初值亦取 $T_1$ =1 428 K，梯度投影-拉格朗日算法与梯度投影法计算初值由式(33)给出。

表 2为文献中及采用不同算法时克劳斯反应的计算结果，从表中可以看出，梯度投影-拉格朗日算法与拉格朗日乘数法的计算结果相同，但是其计算时间要明显小于后者，并且两者均远远小于梯度投影算法，说明梯度投影-拉格朗日算法是有效的，优于拉格朗日乘数法和梯度投影法。

梯度投影-拉格朗日算法中关键参数 $\varepsilon_2$ ， ${\delta ^{(1)}}$ ， $\varphi$ 及温度初值 $T$ 会影响计算时间，结果见图 1~图 4。

从图 1~图 3中可以看出，计算时间随着 $\varepsilon_2$ ， ${\delta ^{(1)}}$ ， $\varphi$ 值的增大而减少，但是为了保证算法的收敛性，应控制 $\varepsilon_2$ ， ${\delta ^{(1)}}$ ， $\varphi$ 的大小，根据计算经验推荐选取 $\varepsilon_2 \leqslant 10$ ， ${\delta ^{(1)}}=0.1$ ， $\varphi=1.1$ 。

从图 4可以看出，温度初值 $T$ 与平衡温度越接近，计算时间越少，因此，合理估算平衡温度有利于提高算法的收敛速度。从图 4还可以看出，在各参数大范围变化的情况下，算法是稳定的。

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$ ，并且合理估算平衡温度有利于提高算法的收敛速度。