1 前 言

换热网络优化方法主要分为三大类^[1]：基于热力学的启发试探和夹点分析法^[2]，数学规划法^[3]和人工智能法^[4]。在数学规划方法中，换热网络合成被抽象成一个混合整数非线性规划问题^[5]，即在所有满足特定目标的超结构解集中寻找最优解。最优目标则通常设定为网络的最小年度费用。

采用数学规划方法解决换热网络综合问题的第一步是建立换热网络超结构。典型的超结构包括：超级结构^[6]、分级超结构^[5]、区间超结构^[7]和其他创新超结构^[8]。其中包含绝大多数结构可能性的分级超结构应用最为广泛，在换热网络的模型中较易表达且易获得可行解。理论上，数学规划法可以解决所有的换热网络优化问题，但在实际应用中，对换热网络的影响因素非常多，相互关系也十分复杂。建立考虑所有影响因素的换热流程优化数学模型难度极大，且针对建立的复杂数学模型，其全局最优求解问题也是一大难点。近年来，换热网络优化模型的求解策略研究也成为国内外研究热点。

大型的、非线性的、非凸型的数学规划模型面临求解困难或者求解效率低等难题。为得到较好的可行解，针对不同的混合整数非线性模型求解问题可采用如下两种方法：一是建立求解或分解策略简化模型^[9]。雷杨等针对延迟焦化换热网络建模问题，在混合整数非线性模型的基础上，利用枚举迭代法消除整数变量，实现混合整数非线性问题简化为非线性规划问题，并给出了在数学规划软件MATLAB上实现模型快速求解的算法；二是研究针对具体问题的特定算法^[10]。针对一般性的混合整数非线性模型，一般性求解器如DICOPT和LOGMIP即可给出较优的可行解^[10]。常用的求解算法包括：遗传算法^[11]、禁忌搜索^[12]、模拟退火^[13]、进化算法^[14]、调和搜寻算法^[15]以及其他混合算法^[16]。Zamora等^[17]提出一个求解混合整数非线性模型的有效方法，模型引入两个估计因子，并假设费用函数线性，采用算术平均温差取代对数平均温差，计算结果误差较小。Björk和Westerlund^[18]评价了等温混合假设在模型中的影响，建立了非等温混合假设的模型并给出了全局优化算法用于求解该模型。Huang等^[19]给出了考虑非等温混合情况下换热网络优化的混合整数非线性模型和相应的求解策略。胡向柏等^[20]为了克服换热网络全局最优化过程中极易陷入局部最优解陷阱的难题，将蒙特卡罗随机抽样技术应用到换热网络冷、热流体随机组合中，随机变化换热网络多维优化参数的优化顺序，实现全局最优化。Lei等^[21]以复杂分馏塔的取热参数为耦合变量构建分馏与换热网络同时优化的混合整数非线性模型，并调用多种商用求解器进行模型求解。

当约束条件与目标函数同时具有单调线性或具有凸函数性质时，此类数学规划问题可以得到有效解决^[22]。但是，当两者中至少有一个是非凸函数时，非线性规划的求解问题则变得十分复杂。换热网络优化的超结构模型通常涉及非凸表达和二元变量，现有的求解器很可能难以获得全局最优解，部分实例中甚至无法得到可行解。一方面考虑开发普适性的一般性算法；另一方面可针对应用背景特点从源头简化物理模型后开发相应的优化方法。

石化工业中典型的换热网络可采用超结构的物理模型，通常可行域为线性约束集，目标函数为非凸函数，是典型非线性规划问题。大多难以求出全局极小值，只能解得局部极小值。但是，对于非线性规划问题中的凸规划问题，其任意局部极小值一定是全局极小值，可利用一般求解非线性规划局部极小值技术求解凸函数规划全局最优值^[18]。针对石化工业典型换热网络综合模型的特点，本文提出基于凸化非凸项技术将非线性项转化为凸函数的全局优化解求解策略。

2 问题描述

石化装置换热网络作为热量回收的重要手段，具有流股数目众多、热负荷大、与公用工程系统热量交换频繁的特点。由于石化装置通常是将一股复杂石油馏分通过反应或者分离等工艺得到各种馏程要求不同的产品，如催化裂化、延迟焦化等。且在工程实践中，由于控制和空间的限制，工艺装置进料的分流数目通常受到限制，远低于相应的热物流数目。因此，其换热网络通常是由一股冷物流，多股热物流构成的^[23]，如图 1所示为延迟焦化装置的分馏与换热集成系统^[24]。

焦炭塔顶逸出的温度为420℃的过热反应油气携带了大量的热量，通过分馏塔侧线循环取热、侧线产品、塔顶油气等进入换热系统，该换热系统的主要目的是加热原料油，从而提高塔底油进炉温度，节省热公用工程(燃料气)，工程中常考虑图 1所示的4条热物流作为热源。因此，该典型的石化工艺装置换热网络在分级换热超结构模型的基础上简化为如图 2所示的超结构。

为简化问题，提出以下假设：(1) 仅考虑逆流传热；(2) 恒传热系数；(3) 恒热容流率；(4) 工艺物流均无相变；(5) 每两股物流至多一次换热。

3 数学模型 3.1 基础模型

基于分级超结构构建以延迟焦化原料油预热流程为研究对象的基础模型，为简化模型，减少非线性项，使模型易于收敛，基础模型假设等温混合。模型同时优化换热面积和能耗费用，最终得到最小年度总费用和最优网络结构。基础模型所含方程中包含的目标函数、各物流总的热平衡、各换热级热平衡、进口温度定义、温度的单调变化约束、冷热公用工程方程，逻辑性约束(二元变量表示超结构内某个特定位置上是否存在换热匹配)、最小温差约束及必要的变量边界约束分别见参考文献^[21]中的方程(2)~(10)。

3.2 全局优化策略

基础模型中采用凹形面积费用公式和非线性对数平均温差近似公式，利用外逼近法和罚函数法的复合运算可以快速收敛得到局部最优值。在没有适合的全局优化求解策略的前提下，无法确定是否可以找到更优的结果，因为无法判断所得结果是局部优化还是全局最优解。针对焦化原料预热所建立的换热流程优化基础模型是典型的混合整数非线性规划问题，其中约束条件是线性函数，构成线性可行域，目标函数f(x)是唯一非线性项且为非凸函数。因此，考虑凸化模型中非凸项，将非线性规划问题转化为凸规划问题，是求解全局最优值的有效方法。该策略的特点在于以凸化非凸项技术为基础，代替常用的分支定界法全局求解策略。

基础模型的非线性项集中在面积费用公式和非线性对数平均温差近似公式。将非线性部分转移到约束条件中，利用A和A_cu将目标函数中面积费用公式线性化。同时，常用的三种对数平均温差近似公式中，考虑应用Paterson近似公式代替Chen近似公式。Paterson近似方程在应用中没有引入新的非凸项，此时，只有关于q的约束$q/{{A}^{1/\beta }}$为唯一的非凸项，其它项均为凸项[41]。

通过以上一系列的变形之后，$q/{{A}^{1/\beta }}$成为现有模型中的唯一非凸项。因此，模型的凸化目标即凸化$q/{{A}^{1/\beta }}$项。引入一个新的变量Q，令$Q = \ln q$，则q与Q存在以下转换关系。

$q=\exp (Q)$

(1)

非凸项$q/{{A}^{1/\beta }}$可以改写为如下凸函数。

$\frac{q}{{{A}^{1/\beta }}}=\frac{\exp (Q)}{{{A}^{1/\beta }}}$

(2)

根据凸函数定义可以验证，在非负数范围内为凸函数。但是，为实现项的凸化过程，模型同时引入了非线性等式约束，该函数为凹函数。因此，之前的变量取代过程实际是将非凸项转移的过程，非凸项仍存在于模型中，模型的凸化目标变成了。热量平衡约束式是非凸的，但凹函数可以近似为一系列线性方程，如图 3所示。通过引入逐段线性的方法将取代为线性方程可实现模型凸化。

逐段线性方程为近似方程，需要引入新的变量$\overset{\wedge }{\mathop{Q}}\,$区分准确和近似的变量。其中$\overset{\wedge }{\mathop{Q}}\,$总是小于真实值。

$\ln (q)\ge \overset{\wedge }{\mathop{Q}}\,=pl\left( \ln (q) \right)\Leftrightarrow \exp \overset{\wedge }{\mathop{(Q)}}\,\le q$

(3)

pl(f) 表示逐段线性方程f。小于真实值求解的过程会使原始模型产生一个松弛的可行域。根据$\overset{\wedge }{\mathop{Q}}\,$计算所得的目标函数值会稍小于根据所得优化结构计算得到的真实目标函数值，通过调整L值的大小使得所得结果与真实值的差异在所接受范围内，关于L值的调节机制如图 4所示，可确定目标函数值及其所代表的优化结构。

3.3 模型改进

应用全局优化策略对模型做如下改进：

l：逐段近似方程中的节点；L：预先设点好的线段节点数；：在节点l处的Q值；：在节点l处的q值；${{b}^{l}}=\left( {{\gamma }^{l+1}}-{{\gamma }^{l}} \right)/\left( {{v}^{l+1}}-{{v}^{l}} \right)$即为第l段线段的斜率值；${{a}^{l}}={{\gamma }^{l}}-{{b}^{l}}{{v}^{l}}$即当第l段线段斜率为零时，其线性方程的常数值。

非负变量：${{\overset{\wedge }{\mathop{Q}}\,}_{ijk}},\overset{\wedge }{\mathop{Q}}\,c{{u}_{i}}$为${{q}_{ijk}},qc{{u}_{i}}$的凸化变量；$q_{ijk}^{l}$为变量q在第l段线段中所占的大小；${{A}_{ijk}},Ac{{u}_{i}}$为换热器和公用工程所需的换热面积。

二元变量：$z_{ijk}^{l},zcu_{i}^{l}$表示q是否落在在第l段线段中。

运用Paterson近似公式后模型面积公式变形为：

$\begin{align} & \frac{{{q}_{ijk}}}{A_{i,j,k}^{1/{{\beta }_{i,j}}}}-\frac{2}{3}{{U}_{i,j}}dt_{i,j,k}^{1/2}dt_{i,j,k+1}^{1/2}-\frac{1}{6}{{U}_{i,j}}dt_{i,j,k}^{{}}-\frac{1}{6}{{U}_{i,j}}dt_{i,j,k+1}^{{}}\le 0, \\ & k\in ST,i\in HP,j\in CP \\ & \frac{qc{{u}_{_{i}}}}{Acu_{i}^{1/{{\beta }_{i,cu}}}}-\frac{2}{3}{{U}_{i,cu}}dtcu_{i}^{1/2}dtcup_{i}^{1/2}-\frac{1}{6}{{U}_{i,cu}}dtcu_{i}^{{}}-\frac{1}{6}{{U}_{i,cu}}dtcup_{i}^{{}}\le 0, \\ & i\in HP \\ \end{align}$

(4)

引入新的变量$\overset{\wedge }{\mathop{Q}}\,$后，面积公式表示为：

$\begin{align} & \frac{{{q}_{ijk}}}{A_{i,j,k}^{1/{{\beta }_{i,j}}}}-\frac{2}{3}{{U}_{i,j}}dt_{i,j,k}^{1/2}dt_{i,j,k+1}^{1/2}-\frac{1}{6}{{U}_{i,j}}dt_{i,j,k}^{{}}-\frac{1}{6}{{U}_{i,j}}dt_{i,j,k+1}^{{}}\le 0, \\ & k\in ST,i\in HP,j\in CP \\ & \frac{qc{{u}_{_{i}}}}{Acu_{i}^{1/{{\beta }_{i,cu}}}}-\frac{2}{3}{{U}_{i,cu}}dtcu_{i}^{1/2}dtcup_{i}^{1/2}-\frac{1}{6}{{U}_{i,cu}}dtcu_{i}^{{}}-\frac{1}{6}{{U}_{i,cu}}dtcup_{i}^{{}}\le 0, \\ & i\in HP \\ \end{align}$

(5)

与${{q}_{ijk}},qc{{u}_{i}}$的${{\overset{\wedge }{\mathop{Q}}\,}_{ijk}},\overset{\wedge }{\mathop{Q}}\,c{{u}_{i}}$线性关系表达式为：

$\begin{align} & {{Q}_{i,j,k}}=\sum\limits_{l=1}^{L-1}{\left( {{a}^{l}}z_{_{i,j,k}}^{l}+{{b}^{l}}q_{_{i,j,k}}^{l} \right)},k\in ST,i\in HP,j\in CP \\ & {{q}_{i,j,k}}=\sum\limits_{l=1}^{L-1}{q_{_{i,j,k}}^{l}},k\in ST,i\in HP,j\in CP \\ & {{\gamma }^{l}}z_{_{i,j,k}}^{l}\le q_{_{i,j,k}}^{l}\le {{\gamma }^{l+1}}z_{_{i,j,k}}^{l},k\in ST,i\in HP,j\in CP,l=1,2,....(L-1) \\ & \sum\limits_{l=1}^{L-1}{z_{_{i,j,k}}^{l}=1},k\in ST,i\in HP,j\in CP \\ \end{align}$

(6)

此时，目标方程简化为：

$\begin{align} & TAC=\sum\limits_{i\in HP}{CCUqc{{u}_{i}}}+\sum\limits_{j\in CP}{CHUqh{{u}_{j}}}+\sum\limits_{i\in HP}{\sum\limits_{j\in CP}{\sum\limits_{k\in ST}{C{{F}_{i,j}}{{z}_{i,j,k}}}}} \\ & \begin{matrix} {} & {} \\ \end{matrix}+\sum\limits_{i\in HP}{C{{F}_{i,cu}}}{{z}_{cui}}+\sum\limits_{i\in HP}{\sum\limits_{j\in CP}{\sum\limits_{k\in ST}{{{C}_{i,j}}}}}{{A}_{i,j,k}}+\sum\limits_{i\in HP}{{{C}_{i,CU}}}Ac{{u}_{i}} \\ \end{align}$

(7)

应用逐段线性近似方程使原始问题的可行域得到松弛，造成所得计算结果小于真实值，可分别看作该问题的上界和下界，以两者的差异作为判据，若相差过大，则节点密度过小，需增加节点数量，反之则减小节点数量，如此反复调节直至差距在可接受范围内，此时认为已找到全局最优解。

4 案例研究 4.1 基础数据

在流程模拟软件PRO/II中建立复杂分馏塔的稳态模拟模型^[25]。通过模拟结果获得塔的操作条件和部分参数。由图 1辨识系统中的冷热物流，见表 1。表 2为物流基础。

回收期pt= 1 a；加热炉热效率η = 90%；年操作时间h = 8400 h；最小传热温差：ε = 10℃；冷公用工程冷却水进口温度TIN_CU = 20℃；出口温度TOUT_CU= 30℃；所有物流膜传热系数为100 W⋅m^-2⋅℃^-1。

换热器设备费用：

$COST\text{ }=\text{ }CF\text{ }+\text{ }C\cdot {{\left( Area \right)}^{B}}$

(8)

其中，CF为固定费用，5500 CNY；C为换热器费用系数，150 CNY⋅m^-2；B为面积费用因子，1。

4.2 计算结果与讨论

基础模型与运用全局优化策略计算的结果分别如图 5和图 6所示。采用GAMS23.0求解计算时间分别约为1和0.6秒(Windows XP，CPU G2020，2.89GHz，RAM 3.19 GB)。

基础模型与全局优化的结果对比见表 3和表 4。

从图 5和图 6可以看出，全局优化模型的网络结构发生了变化，冷物流在第四级发生了分流，冷热物流换热设备数增加，冷公用工程设备数减少。说明运用全局优化策略对优化结构产生了影响。从表 3可以发现，换热网络结构改变的同时，换热面积也发生了改变，全局优化模型总换热面积小于基础模型总换热面积。通过表 4可以看出全局优化模型年度总费用小于基础模型，这是由于全局优化模型总换热面积减小，导致换热器面积投资费用减少。全局优化策略的应用引起换热网络结构的变化，在总传热量不变的情况下传热温差增大，最终引起换热面积的减小，即在求解过程中实现了投资费用与操作参数的权衡。由于非线性规划问题中的凸规划问题，其任意局部极小值一定是全局极小值，利用一般非线性规划局部极小值求解技术，可以得到凸函数规划的全局最优值。实例研究结果表明，利用凸化非凸项技术，将非线性规划问题转化为凸规划问题，所得优化结构即为该问题的全局最优解。其中，凸化后的松弛目标值与相应结构真实值接近，即可认为该优化结构为所求最优结构。与基础模型比较，应用全局优化策略所得的年度总费用减少3.4×10⁶CNY⋅a^-1。

5 结 论

针对典型石化工艺过程换热流程优化陷入局部最优解或无解的困境，提出利用凸化非凸项技术，将非线性规划问题转化为凸规划问题，可求得混合整数非线性模型的全局最优解。除数学变形转移非凸项手段外，提出采用逐段线性化的方法，松弛原有可行区间，得到最优结构。并通过对上下界差异大小的判定，确定分段数值，最终得到目标函数的全局最优解。

通过延迟焦化装置的换热网络实例研究，对比基础模型与全局优化模型所得换热网络优化结果，发现全局模型总换热面积减小，导致换热器面积投资费用减少，最终年度总费用减少了3.4×10⁶ CNY⋅a^-1，凸化后的松弛目标值与相应结构真实值接近。结果证明了全局求解策略的有效性和实用性。

符号说明：

A	- 换热面积，m2	Q	- 换热量q的转换变量，kW
Acu	- 冷却器换热面积，m2	$\overset{\wedge }{\mathop{Q}}\,$	- 换热量Q的转换变量，kW
L	- 线性分段数	TAC	- 年度总费用，CNY⋅a-1
q	- 换热量，kW