2. 哈尔滨哈锅工业锅炉公司，黑龙江 哈尔滨 150046
2. Harbin HBC Industrial Branch Co., Harbin 150046, China
在钢铁、炼油、化工等高能耗产业中，换热网络(Heat Exchanger Network，HEN)是实现能量回收利用的重要手段，其优化对提高整个系统的能量利用率，降低工程中的能耗及综合费用具有重要意义。目前，换热网络综合技术主要分为夹点分析法和数学规划法两类[1]。夹点技术[2]最早由Linnhoff等提出，具有物理意义清晰，操作简单等优点，但不能权衡投资费用与运行费用，因此很难得到综合费用最优的网络设计。数学规划法是基于一些数值方法求解换热网络综合问题的常用方法，以同步综合分级超结构模型[3](Stage-wise Superstructure Model，SWS)最为经典，该数学模型以年综合费用最小为目标，是混合整数非线性规划(Mixed Integer Nonlinear Programming, MINLP)问题，即随着问题规模的增大，计算的复杂性将成倍增加。由于确定性算法对问题特征性信息依赖较强，因此在处理此类具有严重非凸、非线性的MINLP问题时极易陷入局部最优解，而启发式算法由于其随机特性具有操作方便、适应性强和计算效率高等优点，则越来越受到研究者的青睐。
常见的启发式算法，包括遗传算法(GA)[4~6]、模拟退火算法(SA)[7, 8]、粒子群算法(PSO)[9]、微分进化算(DE)[10, 11]等，已被广泛应用到换热网络最优化中。这些算法在处理整型变量或连续型变量方面各有优势，但对于处理较为复杂的换热网络问题时，往往受限于算法本身的进化机理，不能较好的处理换热网络中整型变量和连续型变量的同步优化问题，只能获得次优解。于是，一些学者提出双层混合优化算法优化换热网络，如Huo等[12]提出外层采用GA优化网络的结构，将得到的结构传到内层，采用PSO优化换热器的换热量分布；Khorasany和Fesanghary[13]提出外层使用和声搜索算法(HS)并且内层混合采用HS和序列二次规划法(SQP)的双层优化算法。这些双层混合优化算法虽然能够充分发挥各种启发式算法的优势，但操作性较差，略显复杂，并且对于大规模换热网络优化问题，计算效率不高。
布谷鸟搜索(Cuckoo Search, CS)算法[14]是一种新型元启发式算法，最初由学者Xin-she Yang和Suash Deb依据布谷鸟的“巢寄生”行为和莱维飞行(Lévy flights)机制提出于2009年。CS算法具有操作简单、易于实现、参数少、随机搜索路径优和寻优能力强等优点，并且在连续型变量的多峰目标函数优化方面相对于GA、PSO较为突出[14]，并已成功应用于实际工程优化等问题[15]，但目前暂未发现有关文献将其应用于换热网络综合领域。基于CS算法的优点及换热网络优化特点，本文提出一种改进的CS算法，这种算法能够实现换热网络中整型变量与连续型变量的同步优化，同时具有较强的全局搜索能力，通过具体算例验证改进CS算法的有效性，获得了优于其它优化算法取得的结果。
2 换热网络数学模型理论上分流比不分流的换热网络模型包含更为丰富的网络结构，得到最优换热结构和参数的可能性较大，但最优解对优化算法的性能要求很高。对于相同规模的换热网络综合问题，无分流模型相对于有分流模型具有较少的优化变量和约束条件，模型较为简单且不需要增加阀门、管路等新的费用项目。由于无分流模型并没有引入流股分流比优化变量，若单用CS算法则只需对流股匹配位置处的换热量(或换热面积)进行具有Lévy flights随机特性的增大或缩小优化，无需对流股分流比优化变量进行优化，因此在求解复杂换热网络综合问题时，CS算法的搜索时间成本和求解精度会得到一定控制和保证，能够在可接受的时间范围内得到较优的无分流网络结构设计，且部分设计的年综合费用甚至会低于文献中有分流的网络结构设计。基于此，本文采用无分流分级超结构模型[3]，以NH=2股热流体和NC=3股冷流体为例，设置级数NS=2，如图 1所示。
|
图 1 无分流分级超结构换热网络模型 Fig.1 Stage-wise superstructure model with no stream splits for heat exchanger networks |
分级超结构的级数NS可根据经验灵活设置，但最大级数不应超过max(NH, NC)，最大换热器个数为NH×NC×NS。若以年综合费用(Total Annual Cost，TAC)最小为优化目标，则目标函数表达式为：
| $\begin{array}{l} \min {\rm{TAC}} = {C_{{\rm{fix}}}}\sum\limits_{k \in NS} {\sum\limits_{j \in NC} {\sum\limits_{i \in NH} {{Z_{ijk}}} } } + {C_{{\rm{fix}}}}\sum\limits_{i \in NH} {{Z_{{\rm{cu}},i}}} + {C_{{\rm{fix}}}}\sum\limits_{j \in HC} {{Z_{{\rm{hu}},j}}{\rm{ + }}} {\kern 1pt} {\kern 1pt} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{{\rm{area}}}}\sum\limits_{k \in NS} {\sum\limits_{j \in NC} {\sum\limits_{i \in NH} {A_{ijk}^{\rm{e}}} } } + {C_{{\rm{area}}}}\sum\limits_{i \in NH} {A_{{\rm{cu}},i}^{\rm{e}}} + {C_{{\rm{area}}}}\sum\limits_{j \in NC} {A_{{\rm{hu}},j}^{\rm{e}}} {\rm{ + }}{\kern 1pt} {\kern 1pt} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_{{\rm{cu}}}}\sum\limits_{i \in NH} {{Q_{{\rm{cu}},i}}} + {C_{{\rm{hu}}}}\sum\limits_{j \in NC} {{Q_{{\rm{hu}},j}}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \end{array}$ | (1) |
目标函数由三部分组成，分别是所有换热器的固定投资费用，面积费用和所有冷、热公用工程的消耗费用。其中，目标函数前三项为固定投资费用，Z是0-1整型变量，代表换热器的无或有，其下标代表换热器的位置，Cfix为单个换热器的固定投资费用；第四项至第六项是换热器的面积费用，面积大小由A表示，Carea和指数e是面积费用的参数；最后两项分别指冷、热公用工程的消耗费用，Qcu, i和Qhu, j分别是冷、热公用工程的换热量，Ccu和Chu分别指冷、热公用工程的费用参数。
假设所有换热器采用逆流传热方式，若以换热器的换热量Q作为优化变量，则目标函数中所有换热器的换热面积可由以下公式推导：
| ${A_{ijk}} = {Q_{ijk}}/({K_{ij}}{\rm{LMT}}{{\rm{D}}_{ijk}}),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC,{\kern 1pt} {\kern 1pt} {\kern 1pt} k \in NS$ | (2) |
| ${A_{{\rm{cu}},i}} = {Q_{{\rm{cu}},i}}/({K_{{\rm{cu}},i}}{\rm{LMT}}{{\rm{D}}_{{\rm{cu}},i}}),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH$ | (3) |
| ${A_{{\rm{hu}},j}} = {Q_{{\rm{hu}},j}}/({K_{{\rm{hu}},j}}{\rm{LMT}}{{\rm{D}}_{{\rm{hu}},j}}),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC$ | (4) |
其中，K是冷、热两条流股换热的总传热系数，LMTD是换热器换热的对数平均温差。关于平均温差的计算，文献中一般采用近似计算法或对数平均温差法，在计算中采用对数平均温差法。注意，当冷、热两条流股的热容流率相同时，公式(2)~(4) 中的LMTD由算术平均温差计算。
换热网络综合是一个复杂的工程实际问题，将其转化为数学问题后，该数学模型具有大量的约束条件，其主要约束如下：
(1) 等式约束
单股流体热平衡：
| $FC{p_i}(T_i^{{\rm{IN}}} - T_i^{{\rm{OUT}}}) = \sum\limits_{k \in NS} {\sum\limits_{j \in NC} {{Q_{ijk}}} } + {Q_{{\rm{cu}},i}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH$ | (5) |
| $FC{p_j}(T_j^{{\rm{OUT}}} - T_j^{{\rm{IN}}}) = \sum\limits_{k \in NS} {\sum\limits_{i \in NH} {{Q_{ijk}}} } + {Q_{{\rm{hu}},j}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC$ | (6) |
能量回收换热器热平衡：
| ${Q_{ijk}} = FC{p_i}(T_{{\rm{h}},ijk}^{{\rm{in}}} - T_{{\rm{h}},ijk}^{{\rm{out}}}) = FC{p_j}(T_{{\rm{c}},ijk}^{{\rm{out}}} - T_{{\rm{c}},ijk}^{{\rm{in}}}),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC{\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k \in NS$ | (7) |
冷、热公用工程热平衡：
| ${Q_{{\rm{cu}},i}} = FC{p_i}(T_{i,NS}^{{\rm{out}}} - T_i^{{\rm{OUT}}}){\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH$ | (8) |
| ${Q_{{\rm{hu}},j}} = FC{p_j}(T_j^{{\rm{OUT}}} - T_{j,NS}^{{\rm{out}}}){\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC$ | (9) |
(2) 不等式约束
最小传热温差：
| $T_{{\rm{h}},ijk}^{{\rm{in}}} - T_{{\rm{c}},ijk}^{{\rm{out}}} \ge {\rm{\Delta }}{T_{\min }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in HC,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k \in NS$ | (10) |
| $T_{{\rm{h}},ijk}^{{\rm{out}}} - T_{{\rm{c}},ijk}^{{\rm{in}}} \ge {\rm{\Delta }}{T_{\min }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in HC,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k \in NS$ | (11) |
| $T_{i,NS}^{{\rm{in}}} - T_{{\rm{cu}}}^{{\rm{OUT}}} \ge {\rm{\Delta }}{T_{\min }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH$ | (12) |
| $T_{i,NS}^{{\rm{out}}} - T_{{\rm{cu}}}^{{\rm{IN}}} \ge {\rm{\Delta }}{T_{\min }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH$ | (13) |
| $T_{{\rm{hu}}}^{{\rm{IN}}} - T_{j,1}^{{\rm{out}}} \ge {\rm{\Delta }}{T_{\min }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC$ | (14) |
| $T_{{\rm{hu}}}^{{\rm{OUT}}} - T_{j,1}^{{\rm{in}}} \ge {\rm{\Delta }}{T_{\min }},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC$ | (15) |
式中，ΔTmin为最小传热温差。在程序中，ΔTmin会根据目标函数值得到优化，并允许有任何大于零的数值，但ΔTmin需要被预先设置，文中取0.5℃。
换热量和面积非负要求：
| ${A_{ijk}},{A_{{\rm{cu}},i}},{A_{{\rm{hu}},j}} \ge 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k \in NS$ | (16) |
| ${Q_{ijk}},{Q_{{\rm{cu}},i}},{Q_{{\rm{hu}},j}} \ge 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k \in NS$ | (17) |
目标温度误差要求：
| ${\rm{abs}}(T_i^{{\rm{OUT}}} - T_{i,NS}^{{\rm{out}}}) \le 0.1 \times {10^{ - 4}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH$ | (18) |
| ${\rm{abs}}(T_j^{{\rm{OUT}}} - T_{j,1}^{{\rm{out}}}) \le 0.1 \times {10^{ - 4}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \in NC$ | (19) |
在自然界中，布谷鸟是通过随机或拟随机的飞行方式寻找合适产卵的鸟窝位置，而布谷鸟产的卵又存在一定的概率被宿主鸟发现，为了便于描述布谷鸟搜索算法，对布谷鸟“巢寄生”行为做出三个理想的假设[14]：(1) 每只布谷鸟每次只产一个鸟蛋，并且将该鸟蛋放入随机选择的鸟巢中；(2) 存有高质量鸟蛋的最好鸟巢将会被保留到下一代；(3) 已知宿主鸟巢的数量是一定的，并且宿主鸟发现外来鸟蛋的概率为Pa∈[0, 1]。如果发现外来鸟蛋，则宿主鸟重新建立一个鸟巢。
基于以上三个假设，CS算法主要有两个循环过程，分别为循环体1即布谷鸟随机产卵过程和循环体2即宿主鸟发现外来鸟蛋并建立新巢过程。其中，循环体1中单个布谷鸟通过Lévy flights随机游走的位置和路径更新公式为：
| $\mathit{\boldsymbol{X}}_m^{g'} = \mathit{\boldsymbol{X}}_m^g + \alpha \oplus \mathit{\boldsymbol{Lévy}}(\beta ),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} m \in N$ | (20) |
式中，Xmg为第g代第m个个体，N为布谷鸟个数即种群规模；符号
| $Lévy(\beta ) \approx \frac{u}{{|v{|^{1/\beta }}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (0.3 ＜ \beta ＜ 1.99)$ | (21) |
其中，取β=1.5，u和v服从正态分布，即：
| $u \sim N\left( {0,\sigma _u^2} \right){\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} v \sim N\left( {0,\sigma _u^2} \right)$ | (22) |
其中，
| ${{\sigma }_{u}}={{\left\{ \frac{\mathit{\Gamma }\left( 1+\beta \right)\sin \left( {\pi \beta }/{2}\; \right)}{\mathit{\Gamma }\left[ {\left( 1+\beta \right)}/{2}\; \right]\beta {{2}^{{\left( \beta -1 \right)}/{2}\;}}} \right\}}^{{1}/{\beta }\;}},\quad {{\sigma }_{v}}=1$ | (23) |
其中，Γ是标准的Gamma函数。
公式(20) 本质上是一个随机游走方程，而一次随机游走就是一个马尔科夫链，它的下一步位置(Qmg')依赖于当前位置(Qmg)和转移概率(α⊕Lévy(β))。转移概率中的Lévy(β)项服从Lévy过程的重尾稳定分布，该分布具有无限的均值和无限的方差，并且可以实现经常性的短距离移动和偶然性的较长距离的跳跃交替，这可以保证在不确定的环境中极大限度的提高搜索效率。
通过循环体1位置更新后，进入循环体2中，每个个体以1-Pa的概率通过偏好随机游走组件产生新解，公式如下：
| $\mathit{\boldsymbol{X}}_m^{g''} = \mathit{\boldsymbol{X}}_m^g + \mathit{\boldsymbol{rand}} \oplus (\mathit{\boldsymbol{X}}_{m{\rm{1}}}^g - \mathit{\boldsymbol{X}}_{m{\rm{2}}}^g)m,m{\rm{1}},m{\rm{2}} \in N$ | (24) |
其中，Xm1g和Xm1g是从种群N中随机抽出的两个个体，rand是服从(0, 1) 均匀分布的随机数序列。
循环体1和循环体2更新结束后，保留适应值较好的一组解，完成一次迭代更新。
3.2 CS算法应用于换热网络综合在CS算法中，鸟巢中的鸟蛋相当于一个旧解，而布谷鸟产在鸟巢中的鸟蛋相当于一个新解，将CS算法应用于换热网络综合设计，则一个换热网络设计方案相当于一个鸟蛋，即一个解，但这个解较为复杂，包括换热网络中换热器的位置分布(整型变量)及每个换热器的换热量大小(连续型变量)。CS算法直接优化换热网络结果显示，大量小负荷换热器无法消去，TAC没有得到有效的下降，不能得到符合工业生产需求的换热网络设计。其主要原因是换热网络综合模型是典型的MINLP问题，而CS算法只能处理连续性优化问题，虽然具有Lévy flights随机特性跳出局部最优能力，但直接处理具有整型变量和连续型变量的复杂换热网络有极大难度，本文则基于无分流超结构模型，对CS算法进行针对换热网络模型具有组合特性的优化特点进行改进，使改进的布谷鸟搜索算法能够同时优化整型变量和连续型变量，高效求解换热网络综合问题。
3.2.1 整型变量优化策略整型变量的优化是一个寻找最优流股匹配的过程，即确定最优换热网络结构中换热器的位置及个数。对于小规模换热网络，虽然CS算法可以将不合理结构处存在的换热器的换热量优化到零，达到消去换热器的目的，从而促使结构不断进化，但是对于中、大型规模的换热网络，优化变量过多，这将极大降低CS算法搜索最优结构的效率，同时因部分换热器无法及时消掉易陷入局部最优解，因此在CS算法优化换热量这一连续变量的同时，换热器的高效生成与消去对整型变量的优化起着关键性的作用。针对以上问题，本小节提出以下三种策略来增强CS算法的结构进化能力。
策略1：改进步长控制量
针对复杂换热网络优化问题，虽然Lévy flights组件偶尔会较大步长飞行，具有跳出局部最优解的能力，但缺乏自适应性，严重影响全局最优化的整体效率。引用文献[15]中计算步长控制量α的方法，使所有个体不断向当代最优个体结构靠近，加快种群结构进化，则循环体1中的Lévy flights随机游走的位置更新公式为：
| $\mathit{\boldsymbol{Q}}_m^{g'} = \mathit{\boldsymbol{Q}}_m^g + {\alpha _1}(\mathit{\boldsymbol{Q}}_m^g - \mathit{\boldsymbol{Q}}_{best}^g) \oplus \frac{u}{{|v{|^{{\rm{1/}}\beta }}}} \oplus \mathit{\boldsymbol{randn}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} m \in N$ | (25) |
其中，αl是步长缩放因子，Q指换热网络中换热器的换热量向量，Qbestg是第g代种群中最优个体，randn是服从标准正态分布的随机数序列。注意，虽然公式(25) 能够加快种群结构的进化速度，但同时会在进化后期出现种群多样性缺失的现象，这个问题在改进循环体2部分得到一定解决。
策略2：设置最小热负荷
一般情况，较优的换热网络设计是不存在较小换热量的换热器，自动消去较小换热量的换热器不仅能够加快网络结构的进化，提高搜索效率，而且能够保证最优结构设计不存在极小的换热器，符合工业的设计要求。实现方法是设置一个换热器生成或消去的阀值，应用于循环体1和循环体2中，实现方法如下：
如果
策略3：一定概率接受差解
基于模拟退火算法的跳出机制，此处采用一种简化形式，即以一定的概率接受差解，此策略是给予具有潜在最优结构的一种保护，增加结构进化的可持续性，是一种跳出局部最优策略。该策略仅应用于循环体1中，方法如下：
如果目标函数值f(Qmg')≤f(Qmg)，则保留最优解Qmg'；如果目标函数值f(Qmg')＞f(Qmg)并且rand (0, 1)≤δ，则接受差解Qmg'。其中，δ是接受差解的概率，取值过大则无法保证算法搜索的贪婪性，导致结构进化低效，一般取δ=0.01。
3.2.2 CS算法循环体2的改进结合以上提出的整型变量优化策略，改进后的CS算法只能得到次优的网络设计方案。上文提到，换热网路综合属于MINLP问题，具有严重的非凸、非线性特性，则改进的CS算法应同时具有较强的整型变量和连续型变量处理能力才能获得最优解。虽然原始CS算法具有较强的连续处理能力，但改进公式(25) 使种群结构向目前最优结构进化，削弱了种群的多样性，因此对循环体2做如下改进。
如果rand (0, 1) >Pa，则产生一个随机数rand (0, 1)；如果rand (0, 1)＜Pb，则采用偏好随机游走组件更新解，更新公式为(24)。否则，采用Lévy flights随机游走组件更新解，更新公式为(20)；此处Pb为选择更新概率，取Pb=0.5，则循环体2过程以相等概率强化个体自身的进化(Lévy flights随机游走)和增进个体间信息的交流(偏好随机游走)，从而增强种群个体的多样性，给循环体1更多结构进化空间，避免在进化后期算法因种群多样性缺失而过早收敛。
3.2.3 辅助函数的建立当所得解违反约束条件时，基于罚函数法建立的辅助函数将被调用，目的是将不可行解通过惩罚的方式返回到可行搜索区间。在解的初始化阶段和优化过程中，辅助函数扮演一个重要角色。由于解的初始化是随机过程，并且解很可能会因为Lévy flights飞出边界，因此辅助函数可以保证解的可行性。此外，罚函数较为简单，易于实现，可以评估不可行解偏离可行区间的程度，并且能使不可行解从不同的方向返回到可行区间，从而增加搜索过程的随机性。辅助函数如下：
| $\min F(x,\sigma ) = f(x) + \sigma \{ \sum\limits_{i = 1}^p {{{[\max (0, - {g_i}(x))]}^{{\gamma _1}}} + \sum\limits_{j{\rm{ = }}1}^q {{{[{h_j}(x)]}^{{\gamma _2}}}} } \} $ | (26) |
其中，不等式约束g(x)≥0，等式约束h(x)=0；σ为罚因子，是一个非常大的整数，在此取107；γ1＞1，γ2＞1，均为给定常数，一般取γ1=γ2=2。
3.2.4 改进的CS算法优化换热网络步骤结合以上说明，将改进的布谷鸟搜索(Modified Cuckoo Search, MCS)算法应用于换热网络优化的步骤如下：
步骤1：(初始化)
首先，设置换热网络级数NS，最大迭代次数G，种群规模N，最小换热量阀值Qmin，发现概率Pa，步长缩放因子α1，循环体2中步长控制向量α2。然后，对每个种群个体分量初始化，公式如下：
| $Q_{m,ijk}^0 = \min ({Q_i},{Q_j}) \times rand(0,1),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i \in NH,{\kern 1pt} {\kern 1pt} j \in NC,{\kern 1pt} {\kern 1pt} k \in NS$ | (27) |
| $\left\{ \begin{align} & {{Q}_{i}}=FC{{p}_{i}}\left( T_{i}^{IN}-T_{i}^{OUT} \right),i\in NH \\ & {{Q}_{j}}=FC{{p}_{j}}\left( T_{j}^{OUT}-T_{j}^{IN} \right),j\in NC \\ \end{align} \right.$ | (28) |
其中，Qi和Qj为第i股热物流和第j股冷物流所具有的最大换热潜能。最后，计算每个个体的目标函数值TAC，并找出当前最优解Qbest0和最小TAC fbest0。
步骤2：(循环体1)
记录当代最优个体Qbestg，通过Lévy flights随机游走的位置更新公式(25) 更新第g代每个个体Qmg得新个体Qmg'，比较这两个个体的TAC f(Qmg)和f(Qmg')，如果f(Qmg')小于f(Qmg)，则Qmg'代替Qmg，否则产生一个随机数rand (0, 1)，如果随机数在差解接受概率δ之内，则同样Qmg'代替Qmg。最后，所有第g代个体完成一次更新，更新后的解仍用Qmg表示，找出当代最优个体，用Qbestg表示。
步骤3：(循环体2)
对循环体1更新后的每个个体Qmg(除了最优个体Qbestg)随机产生一个随机数rand(0, 1) 作为宿主鸟发现外来鸟蛋的概率，如果rand(0, 1)＜Pa，即外来鸟蛋没被宿主鸟发现，则Qmg不变；否则再产生一个随机数rand(0, 1)，如果rand(0, 1)＜0.5，则对发现个体进行公式(24) 偏好随机游走更新操作，否则进行公式(20) Lévy flights随机游走更新操作，获得新个体Qmg"，并比较旧解Qmg和新解Qmg"的TAC，仅将较优解代替Qmg。最后，完成循环体2的一次更新，同时也完成第g更新，更新后的解用Qmg+1表示，找出最优个体，用Qbestg表示。
步骤4：(终止判断)
经过循环体1和循环体2两次个体更新后，比较第g+1代最优解Qbestg+1与当前最优解Qbest的TAC，用TAC最小的解代替Qbest。最后，判断是否达到设定的迭代次数，如果是，则输出最优解Qbest及其TAC等相关结果，否则，返回执行步骤2进入下一次迭代操作，直到满足迭代次数要求为止。
改进的CS算法优化流程图如图 2所示：
|
图 2 改进的CS算法优代换热网络流程图 Fig.2 Physical parameters of case 1 |
MCS算法应用于换热网络综合的代码用Fortran95编写，在Compaq Visual Fortran Version 6.6软件平台编译、运行。电脑配置：Windows10，Intel(R) Core(TM) i5-2410M CUP @2.30GHz，4GB RAM。MCS算法中主要涉及参数有迭代次数G，种群规模N，最小换热量阀值Qmin，发现概率Pa，循环体1中步长放缩因子α1，循环体2中步长控制量α2，这些参数的选取均由多次实验确定。在以下两个算例计算中均取Pa=0.25，Qmin=200 kw，α2=10，其它参数根据不同换热网络规模多次实验确定。
4.1 算例1算例1是一个经典的9股流换热网络综合问题，包含4股热流体和5股冷流体，可用热公用工程为热油，冷公用工程为冷却水，流股及费用参数详见表 1。Linnhoff and Ahmad[17]最早对该问题进行研究，获得最优解为2.960×106$·a-1。随后，很多学者以此算例作为基准来检测自己提出的优化HEN综合问题的方法，GA[5]，DE[10, 18~20]，SA[8]，混合GA与PSO[12]以及混沌蚁群(chaotic ant swarm, CAS)[21]等启发式算法对该9股流问题进行研究，不断获得年综合费用较小的网络结构，以上文献优化结果及本文研究结果比较详见表 2。用此算例是为了说明MCS算法比原始的CS算法以及GA、DE等其它启发式方法在处理HEN综合问题方面具有优越性。
|
表 1 算例1物流及费用参数 Table 1 Physical parameters of case 1 |
设置网络级数NS=4，则解的维数是NH×NC×NS=80维，G=1.0×106，N=60，α1=0.02。按照以上参数模拟20次，得到年综合费用在2.924×106~2.930×106$·a-1，均值为2.927×106 $·a-1，稳定性在可接受范围，其中最小的年综合费用为2, 924, 117 $·a-1，耗时2, 211s，网络结构如图 3所示。从表 2比较结果可以看出，改进的CS算法比文献中采用GA，SA，DE，Powell & CAS等算法获得的网络结构费用都低，如比Peng等[8]通过双层SA算法获得的最优结果减少约11, 000 $·a-1，并且计算时间由65, 744 s缩减到2, 211 s；比陈上等[20]通过改进DE算法优化出的最小TAC减少约2, 600 $·a-1，比张春伟等[21]通过Powell & CAS获得的最优TAC减少3, 712 $·a-1。这说明改进的CS算法是有效的，并能够获得更优的换热网络设计。
|
图 3 算例1最优结构(TAC=2, 924, 117 $·a-1) Fig.3 Optimal HEN structure of case 1 (TAC=2, 924, 117 $·a-1) |
|
|
表 2 算例1最优结果比较 Table 2 Optimal results of case 1 |
图 4为CS算法改进前与改进后优化算例1得最小TAC进化曲线比较，其中，“CS”曲线表示采用策略1和策略2，“MCS1”曲线表示采用策略1、策略2及策略3；“MCS2”曲线表示采用了3种策略及循环体2的改进。图 4能够明显的展示CS算法的改进是有效性：原始CS算法(CS曲线)缺乏整型变量处理能力，在2.0×105代之前，进化缓慢，之后完全丧失进化的能力；引入整型变量优化策略后(MCS1曲线)，在2.0×105代之前，进化能力较强，之后由于种群多样性缺失，算法出现早熟收敛问题，获得最小TAC为2, 929, 547 $·a-1；在循环体2内引入自我进化机制后(MCS2)，种群多样性得到一定增强，虽然前期进化能力较MCS1算法弱，但持续进化能力较强，在5.0×105代之后，获得最优解为2, 924, 117 $·a-1，结果优于MCS1算法获得最优结果。
|
图 4 CS、MCS1和MCS2的最小TAC进化曲线 Fig.4 Evolution profiles of the minimum TAC by CS, MCS1, MCS2 |
算例2是一个工程实际问题，由6股热物流和10股冷物流组成，可用热公用工程有高温废气和高温蒸汽，冷公用工程是冷却水，详细参数见表 3。最早，Khorasany和Fesanghary[13]对该工程问题进行分析研究，考虑了网络的分流情况，采用混合HS与SQP获得最优结果为7, 435, 740 $·a-1；Huo等[12]采用GA & PSO优化出年综合费用为7, 385, 856 $·a-1的无分流的换热网络和7, 361, 190 $·a-1的有分流的换热网络；公用工程的消耗量、年综合费用等比较结果详见表 4。
|
|
表 3 算例2物流及费用参数 Table 3 Physical parameters of case 2 |
|
|
表 4 算例2最优结果比较 Table 4 Optimal results of case 2 |
设置NS=3，则解的维数是180维，G=2.0×105，N=30，α1=0.05。与算例1中的参数对比，减少NS是为了避免解的维数过高，降低问题的复杂性；减少N是为了缩减计算时间；增大α1是为了加快结构的收敛速度。按照以上参数设置，经过2.0×105次迭代，获得最优年综合费用为7, 243, 518 $·a-1，用时875 s，最优换热网络如图 5所示。从表 4可知，采用改进的CS算法优化出的最优网络结构费用比Huo等[12]给出的有分流最优结果便宜117, 672 $·a-1，并且计算时间得到极大地缩短，这是因为Huo等[12]采用改进的有分流超结构模型比无分流超结构模型复杂，使优化问题的复杂性增加，同时双层GA & PSO混合算法在优化中大型规模问题时，优化效率及精度体现的并不明显。而基于Lévy flights机制的改进的CS算法具有较强的全局搜索能力与局部搜索能力，可有效的跳出局部最优解，最终获得符合工业生产需求的TAC较小的换热网络结构。
|
图 5 算例 2 最优网络结构图 (TAC=7,243,518 $·a-1) Fig.5 Optimal HEN structure of case 2 (TAC=7,243,518 $·a-1) |
基于无分流分级超结构模型，本文将CS算法应用于换热网络综合领域，并通过两个经典算例对改进的CS算法进行性能验证。CS算法具有较强的连续型变量优化能力，但不能有效处理具有组合特性的换热网络综合问题，而改进后CS算法则能够实现整型变量和连续型变量的同步优化，并且具有较强的全局与局部搜索能力。算例1和算例2的优化结果表明，改进的CS算法能够独立解决中小型甚至大型换热网络综合问题，并且优化结果与GA，SA，DE，PSO及其它混合改进算法相比，在计算精度和效率方面具有一定优越性。同时改进的CS算法也存在一些参数调优方面的问题有待研究，如策略2中最小热负荷的选取对优化结果的灵敏性影响研究等。
符号说明：
| A | -换热器面积，m2 | α2 | -步长控制量 |
| C | -费用系数 | δ | -接受差解的概率 |
| FCp | -热容流率，kW·K-1 | 上标 | |
| G | -最大迭代次数 | e | -换热器面积费用指数 |
| N | -种群规模 | in | -换热器入口 |
| NH | -热物流股数 | out | -换热器出口 |
| NC | -冷物流股数 | g | -迭代次数 |
| NS | -换热网络级数 | IN | -物流入口 |
| Pa | -发现概率 | OUT | -物流出口 |
| Pb | -选择更新概率 | 下标 | |
| Q | -换热器的换热量，kW | area | -面积 |
| Qmin | -最小换热量生成与消去阀值，kW | c | -冷物流 |
| Q | -换热网络中换热器的换热向量 | cu | -冷公用工程 |
| rand | - (0, 1) 均匀分布随机数序列 | fix | -固定投资 |
| randn | -标准正态分布随机数序列 | h | -热物流 |
| rand(0, 1) | - (0, 1) 均匀分布随机数 | hu | -热公用工程 |
| T | -温度，℃ | i | -热物流编号 |
| ΔTmin | -最小传热温差，℃ | j | -冷物流编号 |
| Z | - 0-1整型变量 | k | -级数编号 |
| α1 | -步长缩放因子 | m | -种群中个体位置 |