换热网络优化中换热单元生成频次的影响分析及策略改进

苏戈曼，崔国民，肖媛，赵倩倩

（上海理工大学新能源科学与工程研究所，上海200093）

换 热 网 络（heat exchanger networks，HENs）是过程系统的重要环节，通过对冷热物流进行合理的换热匹配以实现热量回收，对提升节能水平和生产经济性有重要意义。HENs 的综合优化由Rudd[1]首次提出，而当前优化方法主要包括分布优化法与同步优化法[2]。分步优化法[3-5]通过将原问题分解成子问题以降低计算复杂性，但忽略了子问题间的制约关系，且不能兼顾优化整型变量和连续变量，优化结果较差。同步优化法[6-9]通过建立HENs的混合整数非线性规划模型（mixed-integer nonlinear programming，MINLP），对整型变量与连续变量进行同步优化，优化效果较好。

1990年Yee等[9]提出一种经典的同步优化方法，在该方法中，确定性和启发式优化算法被引入以求解分级超结构模型（stage-wise superstructure，SWS）。其中确定性方法包括外部逼近法[10]、分支定界法[11]、填充打洞函数法[12]等，借助精确的数学理论求解MINLP 问题，但受制于问题的非凸非线性，易在求解过程中陷入局部极小值。启发式方法包括模拟退火算法（simulated annealing，SA）[13-14]、遗传算法（genetic algorithm，GA）[15-16]、粒子群算法（particle swarm optimization，PSO）[17-18]等，这类方法对数学模型依赖性小，计算效率较高，全局搜索能力较强，更适用于HENs优化问题。在启发式方法的应用过程中，换热单元生成方式已有一定研究。Luo 等[19]在GA 基础上结合SA 并提出结构控制策略与局部优化策略，强化换热单元生成；Huo等[20]建立了GA和PSO的双层算法，利用GA生成换热单元，产生HENs网络结构；Zhang等[21]基于混沌蚁群算法提出一种结构进化策略，通过设定换热单元个数调控结构进化方向；Peng等[22]基于SA算法，通过蒙特卡罗随机抽样模型对换热单元进行随机生成与消去；此外，肖媛等[23]提出强制进化随机游走算 法 （random walk algorithm with compulsive evolution，RWCE）应用于SWS 模型，通过设置最小换热量和一定的概率模型来控制换热单元的生成和消去，实现了整型变量与连续变量的同步优化，优化效果较好。然而SWS 以“级”为单位，每级内换热匹配总数与顺序固定，使得换热单元生成受限，常需要通过增加级数来实现换热单元的生成。同时，变量维度随着级数增加呈指数级扩大，导致计算规模过大，优化效率低下。

2018 年，肖媛[24]提出节点非结构同步模型（node-wise non-structural superstructure， NWNSS），该模型取消“级”的概念，并以节点形式表征换热单元的匹配位置，通过冷热流体上节点的随机匹配生成新的换热单元，换热单元生成自由，增大了结构灵活发展的空间。该作者将RWCE 应用于NW-NSS 模型，与基于SWS 模型的RWCE 算法相比，前者充分发挥了算法随机灵活的特点，在表现形式上将换热单元的生成操作与连续变量优化进行分离处理，使得对优化路径的调节更为直观、方便。其中换热单元生成操作涉及生成概率与生成个数两个参数，分别控制优化过程中换热单元生成频率和结构进化时换热单元生成个数。本文将二者定义为换热单元生成频次，其对整型变量与连续变量影响重大，对优化过程的调控至关重要，但目前仍缺乏从该角度出发对整型变量与连续变量的影响进行考察的研究，也急需一种方法合理地调控整型变量与连续变量的优化频率。

为了探究换热单元生成频次对整型变量与连续变量优化的影响，本文采用基于NW-NSS 模型的RWCE算法，通过调整算法中的生成概率与生成个数以获得不同的生成频次，进而跟踪不同生成频次下整型变量与连续变量的优化状态及年综合费用。基于此，提出一种策略合理调整换热单元生成频次，从而提升算法优化性能。

1 NW-NSS模型

1.1 问题描述

本文采用NW-NSS 模型，热流体数为NH，冷流体数为NC，每股流体上布置一定数量的节点，通过随机方式分别选定热流体、冷流体上的某节点，形成一个换热匹配，只要该匹配满足约束条件，则认为是有效匹配。以2×2 算例为例，其NW-NSS模型如图1所示。

其中，i、j 分别为热、冷流体编号；SUMH、SUMC分别为每股热、冷流体上的节点总数；NUMH,i、NUMC,j分别为每股热、冷流体上的节点编号；此外，设置PH为某一热节点位于所有热流体中的编号、PC为某一冷节点位于所有冷流体中的编号；设置NCi,NUMC表示第i股热流体上的第NUMi个节点对应的冷节点编号PC，NCj,NUMC表示第j 股冷流体上的第NUMj个节点对应的热节点编号PH。NW-NSS模型中热、冷节点自由匹配，且换热匹配节点编号不重合，因此可由NCj,NUMH,i、Qi,NUMH,i或NCj,NUMC,j、Qj,NUMC,j确定一个换热匹配和换热量大小。为便于观察，NW-NSS模型中的换热匹配可根据其在流股上的相对位置将其垂直处理，处理后的NW-NSS模型如图2所示。

1.2 目标函数

换热网络NW-NSS 模型以最小年综合费用（total annual cost，TAC）为目标，主要包括换热单元、冷公用工程和热公用工程的费用，其目标函数如式(1)所示。

式中，FFIX为换热单元固定投资费用；E 为判断换热单元存在与否的0-1逻辑变量；CA为换热单元面积费用系数；β 为换热单元面积费用指数；Ai,NUMH,i为换热单元的换热面积；ACU,i、AHU,j分别为冷、热公用工程换热面积；CCU、CHU分别为冷、热公用工程运行费用系数；QCU,i、QHU,j分别为冷、热公用工程换热量。

1.3 约束条件

（1）换热网络每股流体上的出口温度应满足目标温度，温度约束如式(2)、式(3)所示。式中，THout,i、TCout,j分别为热、冷流体出口温度；THtar,i、TCtar,j分别为热、冷流体的目标温度。

（2）换热网络结构中单股流体的总换热量等于该流体上所有换热单元换热量与流体末端公用工程换热量的总和，热平衡式如式(4)、式(5)所示。式中，W 为流体热容流率；THin,i、TCin,j分别为热、冷流体进口温度。

图1 NW-NSS模型

图2 垂直处理后的NW-NSS模型

（4）对于冷、热公用工程均有热平衡式如式(7)、式(8)所示，其中为热流体流经最后一个换热单元的出口温度，为冷流体流经第1 个换热单元的出口温度。

（5）此外，还包含换热面积、换热量非负约束［式(9)、式(10)］以及换热单元最小传热温差约束［式(11)、式(12)］，本文最小传热温差ΔTmin均取值为0℃。

2 强制进化随机游走算法

采用RWCE 算法的NW-NSS 模型在优化之初不存在任何换热匹配，通过连接随机选定的热、冷节点生成新的换热匹配，而后对该新结构进行下一步的选择、变异操作，最后判断是否满足终止条件。该过程中换热单元的进化和生成分开进行，RWCE算法具体流程如下。

2.1 种群初始化

产生包含N 个个体的初始种群，每个个体Qn,k（n=1,2,…,N；k=1,2,…,Nk）代表一个Nk维的换热网络结构（Nk为换热单元数目）。个体初始解不进行换热匹配，每股流体的末端均采用公用工程换热使其达到目标温度。

2.2 存在的换热匹配进化

对种群内已存在的换热匹配换热量（Qn,k,it＞0；n=1,2,…,N；k=1,2,…,Nk）进行随机游走，实现个体进化，同时设置换热量最小值以控制换热单元的消去。换热量随机游走公式如式(13)所示。

式中，Qn,k,it、Q'n,k,it+1为第n个个体第k个换热单元随机游走前后的换热量；α1、α2是(0,1)之间的随机数；ΔL 是最大游走步长；Zk为控制是否参与进化的0-1变量，如式(14)所示，δ为进化概率，r1为(0,1)之间的随机数。

在进化操作后，通过最小换热量Qmin判断是否消去该换热匹配，具体操作如式(15)所示。

若某个换热匹配在优化过程中被消去，则对应热、冷节点编号也相应归零，而没有被消去的换热量则为随机游走后的有效换热量。

2.3 生成新的换热匹配

为保证结构进化能力，针对流体上未产生换热单元连接的位置，从冷、热流体上随机选择m对节点（m即为个体单次进化中换热单元生成个数）进行匹配，通过以一定概率随机生成换热单元实现整型变量与连续型变量的交替优化，并按式(16)赋予换热量

式中，Qmax为生成新生成换热单元换热量最大值；φ 为换热单元生成概率；r2为(0,1)之间的随机数。

2.4 选择与变异

式中，Qn,it为随机游走前个体n 的解；Qn,it+1为被选择进入下一代的解；φ 为接受差解概率；γ 为(0,1)内均匀分布的随机数。

2.5 终止条件

当个体游走的迭代步数达到ITmax时，结束迭代，否则返回个体随机游走继续优化。

3 换热单元生成频次的影响分析

基于NW-NSS 模型的RWCE 算法将换热单元生成与换热量优化分开处理，其中换热单元生成操作直接影响优化过程，并且该操作由生成概率与生成个数两个参数共同决定，因此本文将结合具体算例对这两个参数的作用进行分析。

3.1 换热单元生成概率的影响

以9SP[25]为例（算例参数见表1），观察并分析φ 对RWCE 优化过程的影响。设置优化参数N=1、SUMH=50、SUMC=50、ΔL=90、δ=0.5、φ=0.001、ITmax=1×108，m 均为2，φ 分别取0.01、0.05、0.1、0.2、0.3。

表1 9SP算例参数

换热网络优化以TAC 作为目标函数，TAC 的变化情况可反映算法优化质量与效率。鉴于此，本文通过观察TAC 下降情况来分析换热单元生成概率的影响。如表2所示，统计了因新生成换热单元导致最优TAC 更新的次数D1，未生成换热单元时最优TAC 更新的次数D2，因新生成换热单元导致当前TAC 下降的次数D3，未生成换热单元时当前TAC下降的次数D4，同时记录各生成概率下最终优化结果TAC。

表2 不同换热单元生成概率下TAC变化情况统计

从表2中可以看出，在整个优化过程中，D2与D4分别多于D1与D3，即TAC 的下降多数通过固定结构的连续变量优化得到。同时，随着φ 增大，D2、D4减少而D1、D3增加，表明整型变量相对优化频率增多，使结构变异能力增强、连续变量优化削弱。反之，φ 减小时，D2、D4增多而D1、D3减少，虽连续变量优化更为充分但整型变量优化相对优化频率减少，导致结构变异能力较弱。另外，φ过大或过小将分别造成连续变量优化不充分或整型变量优化不充分，且对应TAC 均较差，因此应根据算例特点选择适当的φ。

3.2 换热单元生成个数的影响

同样选取9SP 算例，采用RWCE 进行优化，φ取0.1，m分别取1、2、3、4、5，其他优化参数不变。通过统计在不同生成个数下的D1～D4（表3），分析换热单元生成个数对RWCE 优化过程的影响。

表3 不同换热单元生成个数下TAC变化情况统计

从表3中可以看出，当m为1时，D1、D3最多，生成换热单元后的结构变化较小，最易于被保留，结构更新较为频繁，但也易造成结构还未进行彻底的连续变量优化就被新结构取代；此后，随着m增多，D1、D3均减少，即单次生成的换热单元越多，新生成结构越难优于当前最优结构，甚至可能出现不可行解，在优化过程中可能被直接舍弃，导致结构更新次数减少。同时，当m 较小（3 个及以下）时，D2、D4随着m的增加而增大，即单次生成多个换热单元时结构更新较慢，使得当前的最优结构有足够的时间进行连续变量优化；然而，当m 较大（3个以上）时，D2、D4却随着m的增加而减小，即单次换热单元生成个数的增加造成结构更新次数减少，可能出现始终对同一固定结构进行连续变量优化，并逐渐趋于该固定结构下的局部最优解，因此难以再通过连续变量优化使得费用下降。为进一步验证m 过大对进化的影响，另设置m 为10、15、20，统计结果见表4。

可以看到，m过大时，几乎没有可以通过整型变量优化更新TAC 的结果（D1、D3极少），通过固定结构下的连续变量优化更新TAC 次数也大幅减少（与表3 中D2、D4相比），验证了单次生成过多换热单元使得结构性能变差且难以被保留，个体长期在原结构基础上进行连续变量优化，从而陷于该结构的局部极值，使得连续变量优化因难以突破局部极值而逐渐减少，最终优化结果也较差。因此，在优化过程中m不宜过大，取值一般不超过5。

表4 换热单元生成个数过大时TAC变化情况统计

4 生成参数动态调节策略

从上述分析得知，调节生成概率与生成个数可以协调连续变量优化和整型变量优化。个体在优化过程中的状态是变化的，若在过程中以相同的频次生成换热单元，则忽略了整型变量与连续变量优化的合理时机。因此提出生成参数动态调节策略，通过跟踪个体优化进程，调整各个体生成概率与生成个数，从而维持合理的换热单元生成频次。

4.1 生成参数动态调节策略

在RWCE 算法的优化过程中，整型变量频繁变化易缺乏充分的连续变量优化，反之连续变量优化时间过长则整型变量优化能力较弱，导致结果陷入局部极值。

本文提出生成参数动态调节策略，即在优化过程中跟踪观察个体TAC 变化情况，设置TAC 停滞的迭代步数阈值Cmax。当个体TAC停滞累积步数大于Cmax时视为个体长期停滞，此时给予一定幅度的结构变异扰动，调整生成概率与生成个数，当个体TAC 停滞累积步数小于Cmax时，仍以原生成概率与生成个数进行优化，以此改善换热单元生成频次。具体操作如式(19)～式(21)所示。

式中，φnew为改进后的生成概率；mnew为改进后的生成个数；Cit+1为个体TAC 停滞的迭代步数累积值。为防止因缺乏结构变异能力导致TAC 长期停滞，在TAC 长期停滞时，以较大φ'和较小m'优化，从而提高结构进化能力；为保证较优的结构有充分的连续变量优化，在TAC下降时，以较小φ和较大m优化，从而提高局部搜索能力。具有生成参数动态调节策略的改进RWCE 算法（improved RWCE，IRWCE）流程如图3所示。

图3 IRWCE算法流程

4.2 有效性分析

仍选用9SP算例，对比改进前后的RWCE算法优化曲线如图4所示。RWCE算法优化曲线如曲线1、曲线2所示，其中曲线1取φ=0.2、m=1，曲线2取φ=0.04、m=5，IRWCE算法优化曲线如曲线4所示，取Cmax=104、φ=0.04、m=5、φ'=0.2、m'=1，其他优化参数不变。

图4 9SP算例TAC曲线对比

改进前生成概率与生成个数均为定值，其优化路径（曲线1、曲线2）均在4000 万步左右停滞，最终优化结果费用较高。改进后生成概率与生成个数随个体TAC 进化情况动态变化，其优化路径（曲线3）以更快的速度得到了更优的解，并且在优化后期仍保持一定的下降趋势。通过对比可以看出，采用IRWCE 算法优化效率明显高于原始RWCE算法，在优化过程中既保证了充分的连续变量优化，又提高了结构变异能力。

5 算例验证

将IRWCE 算法应用于两个基准算例，并以年综合费用为评价优化质量的指标验证其优化性能。

5.1 算例1

算例1 取自文献[26]，包含10 股热流体与5 股冷流体，具体参数见表5。分别采用RWCE算法与IRWCE算法应用于该算例，优化参数见表6，改进前后除新增参数m'、φ'、Cmax外，其余优化参数均不变。

采用RWCE 算法，优化结果如图5所示，TAC为5467884USD/a。采用IRWCE 算法，得到结构如图6所示，TAC为5178068USD/a，较RWCE算法的结果节省了289816USD/a，结果对比见表7，较文献最优结果[29]减少了55219USD/a。

对比结构图5、图6 可以看出，采用生成参数动态调节策略后由于根据个体状态调整生成概率与生成个数，促进了个体结构进化，获得了性能优于原算法的结构。此外，这两个结构的C4 流体都未参与匹配，这是因为C4 流体温位较高，与热流体匹配时换热温差小甚至可能出现温度交叉，形成的结构性能差且无法被保留。

5.2 算例2

算例2取自文献[30]，包含7股热流体与3股冷流体，算例参数见表8。将RWCE 算法与IRWCE算法分别应用于算例2，优化参数见表9。采用RWCE 算法的结构如图7 所示， TAC 为8862550USD/a；采用IRWCE 算法的结构如图8 所示，TAC为8693753USD/a。

表5 15SP算例参数

表6 15SP算例优化参数

图5 采用RWCE算法得到的结构图（TAC为5467884USD·a-1）

图6 采用RWCE算法得到的结构图（TAC为5178068USD·a-1）

表7 算例1结果对比

因该算例固定投资费用为零，所以在采用IRWCE 算法的优化过程中，因增加结构扰动而新生成的换热单元更容易被保留下来，最终得到换热单元较多的优化结果（图8），更直观地体现了IRWCE 算法提升结构变异能力的作用。同时，由于该算例面积费用指数为1，换热单元面积费用正比于换热面积，因此流体匹配相同的相邻换热单元合并前后费用不变，合并图8中方框内的换热单元后结构如图9 所示。对比图9 与图7 结构，改进后费用下降了168797USD/a，提高了算法的优化质量。表10 记录了文献与本文结果的具体数据，采用IRWCE 算法取得了优于文献的结果，较文献最优结果[32]下降了14230USD·a-1。

表8 10SP算例参数

表9 10SP算例优化参数

图7 采用RWCE算法得到的结构图（TAC为886255USD·a-1）

图8 采用IRWCE算法得到的结构图（合并前，TAC为8693753USD·a-1）

图9 采用IRWCE算法得到的结构图（合并后，TAC为8693753USD·a-1）

表10 算例2结果对比

5.3 算例3

算例3为包含13股热流体和7股冷流体的大规模算例[33]。算例参数见表11，其中固定投资费用为4000USD/a、面积费用系数与指数分别为500 和0.83。将IRWCE 算法应用于该算例，优化参数见表12，优化后所得结构如图10 所示，TAC 为1403139USD/a。与现有文献TAC 对比见表13，采用IRWCE 算法优化后TAC 较文献[37]下降了9662USD/a。

6 结论

本文探究了HENs 优化中换热单元生成频次对优化的影响，并提出生成参数动态调节策略。得出以下结论。

表11 20SP算例参数

表12 20SP算例优化参数

（1）换热单元生成概率越大时，整型变量相对优化频率增加，结构变异能力较强，但不利于连续变量持续优化；生成概率越小时，固定结构下连续变量的相对优化频率增加，但结构变异能力较弱。

（2）单次进化中换热单元生成个数较少时，结构变化小易被保留，结构更新较快，但容易缺乏充分的连续变量优化；生成个数较多时，新结构变化较大易被舍弃，造成长期以固定结构优化，连续变量优化充分，但结构更新较慢。此外生成个数过多，加剧了结构更新迟缓，导致优化陷入固定结构下的局部极值，连续变量优化因难以突破局部极值而逐渐减少。

（3）生成参数动态调节策略根据个体TAC 变化情况调整换热单元生成频次，使得个体优化既保持适当的整型变量优化能力又满足充分的连续变量优化。改进后的算法优化15SP、10SP 和20SP 算例，均获得相对文献更优的结果。结果表明该策略使得换热单元生成频次更符合个体优化状态，有效地提高了算法的优化性能。

图10 采用IRWCE算法得到的结构图（TAC为1403139USD·a-1）

表13 算例3结果对比

符号说明