附录5 基于粒子群算法的EDA裂解反应精馏塔优化

年产 万吨醋酸乙烯酯副产 万吨PVA废水资源化利用项目 初步设计说明书基于粒子群算法的EDA裂解反应精馏塔优化长城能源年产10万吨醋酸乙烯项目 长城能源年产10万吨醋酸乙烯项目 目录第一章计算原理和方法31.1反应精馏过程稳态模拟方法31.1.1Aspen Plus软件31.2.1 粒子群算法理论41.2.2 粒子群算法的基本步骤61.2.3 粒子群算法的 MATLAB 编译61.3 软件集成平台的建立81.3.1 软件集成技术82.3.2 软件集成平台的建立8第二章 基于粒子群算法的反应精馏塔优化模型的建立102.1 反应精馏塔稳态模型的建立102.2 基于粒子群算法的反应精馏塔的优化设计12第三章 实例分析153.1 EDA裂解反应精馏稳态模型描述及设计方案153.2 决策变量的影响及其变化范围17附录 MATLAB源代码25长城能源年产10万吨醋酸乙烯项目 第一章计算原理和方法1.1反应精馏过程稳态模拟方法反应精馏可通过化工流程模拟软件来模拟。科研工作及工业生产过程中常用的化工流程模拟软件有 ChemCAD、Pro II、Aspen Plus 等软件，其中 Aspen Plus 软件在反应精馏过程模拟的研究中应用十分广泛。1.1.1Aspen Plus软件Aspen Plus 软件是目前最常用的化工过程模拟软件之一，由美国 AspenTech 公司开发。该软件集化工过程稳态模拟、动态模拟、热力学计算等功能于一体， 功能强大。其在反应精馏的过程模拟精确度方面也备受学界及工业届认可。它主要由以下三部分组成： （1） 物性数据库 Aspen Plus 中物性数据库众多，存储了大量的基础物性数据，包括纯组分、 电解质等。另外，对于缺少的物性，Aspen Plus 还可使用物性估算方法进行估计。 （2） 单元操作模块 单元操作模块中存在 50 多种单元操作模型，每个模型都有独特的形式。此 外，Aspen Plus 还拥有用 Fortran 语言编写的用户模型，它是指当 Aspen Plus 系统提供的模型不能满足用户的要求时，用户自行开发的模型。用户可根据自身需 求编写相应的 Fortran 程序来实现期望的功能。 （3） 模型分析工具 模型分析工具中模块众多，如使用较多的灵敏度分析模块、数据拟合模块等。 其中，灵敏度分析模块可用来研究因变量随自变量的变化趋势，从而判断自变量 对因变量的影响程度；此外，可以用此模块来进行简单的系统优化；数据拟合模 块可用来拟合试验数据，以确定物性模型参数，实现更精确的模拟；值得一提的是，Aspen Plus 软件不仅为反应精馏过程模拟提供了强大的物性数据库和单元操作模块，而且其在软件集成方面也具有很大的优势。 Aspen Plus软件与反应精馏塔的模拟由于 Aspen Plus 具有强大的物性数据库，其被广泛应用于反应精馏的模拟及优化上，近年来与反应精馏稳态或动态模拟有关的文献相继发表。在反应精馏的模拟过程中，一般采用 RadFrac 模型进行严格计算。相比于普通精馏过程，反应精馏则需要输入更多的参数，如反应段位置、动力学、反应段滞液量等，但基本的输入规则及运算步骤是类似的。图1-1 Aspen Plus用户界面1.2.1 粒子群算法理论 粒子群算法，也可称之为粒子群优化算法（Particle Swarm Optimization, PSO） 是美国的社会心理学家 James Kennedy 和电气工程师 Russell Eberhart 在对鸟类觅食过程中的迁徙和群聚行为的建模与仿真过程中受研究结果的启发而提出的一种基于群体智能的全局随机搜索算法。粒子群算法易理解、精度高、收敛快，因此一经发现即得到了学界的重视，目前已经被广泛应用于资源调度、管理决策、结构设计等多个领域的目标优化问题。粒子群算法有基本粒子群算法、标准粒子群算法、压缩粒子群算法、离散粒 子群算法等。本文主要介绍标准粒子群算法。 定义 N 个粒子组成一个群落，该群落处于一个 D 维的搜索空间内，其中每 一个粒子可由一个 D 维的向量表示：第 i 个粒子的搜索速度同样为一个 D 维的向量：定义当前第 i 个粒子已搜索到的最优值为个体极值 pbest，对应最优位置为个体最优位置Pb，即定义在整个粒子群内，当前已搜索到的最优值为全局极值 gbest，对应最优位置为全局最优位置Pg，即 粒子的速度及位置更新公式分别如式(1-1)、(1-2)所示 （1-1） （1-2）式中：c1和c2为学习因子，有时也被称作加速常数；r1和r2为0,1范围内的均匀随机数；n ij 代表粒子速度，其取值范围在nmax ,vmax 之内，n max为常数，是算法中的设定值； w 为惯性权重。当 w=1 时，标准粒子群算法与基本粒子群算法一致。w取 值较大时，算法全局收敛能力较强，而局部收敛能力则弱一些；反之，局部的收敛能力相对较强，相应的全局收敛能力较弱。w 取值的推荐范围为0.81.2。 惯性权重 w 的值可以在搜索过程中进行调节，以便权衡局部搜索能力和全 局搜索能力，使优化算法的效果达到最佳。目前，常用的 w 值调整策略为 “线性递减权值策略”，其数学表达式为： （1-3）式中：wmax-最大惯性权重 wmin-最小惯性权重 t-当前迭代次数 Tmax-最多迭代次数 在实际应用案例中，wmax=0.8，wmin=0.4 。1.2.2 粒子群算法的基本步骤 粒子群算法的本质为一种随机搜索算法，由于其在计算速度及全局搜索能力等方面的优势，使其成为是一种极具前景的优化算法。该算法的运算步骤如下所 述： 步骤 1 初始化粒子群，初始化每个粒子的速度Vi和位置Xi，并设定粒子的种群规模 N。其中，Xi为一个 D 维的向量，向量维度与目标函数变量数相等。 步骤 2 根据目标函数计算每一个粒子在搜索范围内所对应的适应度值fiti。 步骤 3 对比每个粒子的适度值 fiti与个体极值 pbest(i)，若适度值 fiti较小，则可替换掉 pbest(i)；反之，个体极值维持不变。 步骤 4 对比每个粒子的个体极值 pbest(i)与全局极值 gbest，若个体极值 pbest(i)较小，则可替换掉 gbest；反之，全局极值维持不变。 步骤 5 根据计算结果，对粒子的速度Vi和位置Xi进行迭代更新。步骤 6 对边界条件进行处理。 步骤 7 判断是否满足算法的终止条件；若满足条件，则终止运算过程并输 出优化结果；若条件未满足，则返回上述步骤 2。 1.2.3 粒子群算法的 MATLAB 编译 在科研工作和工程实践中，存在大量的数学运算，尤其是矩阵运算。目前常用的计算机语言，如 C 语言、C+、Fortran 等，在矩阵运算方面效率不高，市场对矩阵计算能力突出的软件需求量极大。1967 年，美国 Mathworks 公司推出了MATLAB（全称为：MATrix LABoratory）软件包，功能强大，编程及计算的效率极高，与 Mathematica、Maple 一起，被评为三大数学类科技应用软件。 MATLAB 软件包为用户提供了直观间接的开发环境，其优势具体体现在： (1) 编程效率高 MATLAB 的语言逻辑较清晰，与人类书面计算时的思维方式较为接近，因此编写程序简单高效，容易掌握。 (2) 开发环境友好 MATLAB 内含有大量工具箱，可直接调用，且其中大部分工具箱拥有MATLAB 的图形用户界面，操作简单易懂。编程环境简单，能及时报错并提供相应分析，且具有完备的调试系统。 (3) 扩充能力强 一方面，可以方便地调用 MATLAB 自带的库函数；另一方面，可以通过建立Mex文件，与其他计算机语言进行混合编程。 (4) 强大的图形处理功能 MATLAB 软件具有强大的图形处理功能，能够绘制如柱状图、饼状图、三维立体图等各类图形。而且最新的 MATLAB 软件在图形处理方面改进极大，增加了图形的光照处理、色度处理等功能。 (5) 强大的矩阵和数组运算功能 MATLAB 不仅包括各种数学运算需要用到的运算符，而且，最重要的是， MATLAB 语言中无需专门定义数组的维度 ，几乎常用的矩阵函数都存在于其本身自带的库函数中，用户不再需要自己编程，可直接操作，运算高效。这一特点也使其在数据处理、函数优化等方面优势巨大。 (6) 工具箱功能强大 MATLAB 软件内含有大量的工具箱，如控制系统、优化系统工具箱等，这 些工具箱都各领域专家编写，用户可直接调用，无需另行编写。MATLAB 各工 具箱及相应功能如图所示： 图 1-2 MATLAB 工具箱 本文中 MATLAB 软件将被用于反应精馏塔的目标优化，粒子群算法可以用MATLAB 软件进行编译。在计算过程中，MATLAB 软件可以通过 Aspen Plus 软件提供的外部程序接口调用和访问后者中的数据及变量，从而实现目标函数的优化。 1.3 软件集成平台的建立 1.3.1 软件集成技术 将多个现有软件通过一定的方法和策略集成于一个用户软件之中，这个过程称之为“软件集成”。经集成后的软件可以突破单一软件功能的限制，同时实现多个软件的功能，更能满足科学研究及工程应用的需求。 近二三十年，计算机学科在全球范围内发展迅速，随着而来地即是软件方面的突破性进展。许多功能强大的软件逐渐从学术领域伸展到工业实践中。如在化工流程模拟中应用广泛的 Aspen Plus、ProII 等软件；在流体力学计算方面功能强大的 Fluent、ANSYS CFX 等软件不足的是，这些软件功能强大的领域较少，只靠一个软件并不能很好的解决实际的工程问题。因此，软件集成技术应运而生。目前常用的软件集成技术可以分为三类：软件总线技术（Software Bus）、ActiveX技术以及 OLE 技术（Object Linking and Embedding，对象连接与嵌入）。 软件集成技术是软件开发领域的重大突破，该技术使许多专业性强的软件的应用领域得到扩展。通过软件集成技术，在较少的时间和投资成本下，复杂的工程问题得以顺利解决。本文基于 ActiveX 技术建立了软件集成平台，建立了一套数学优化模型，并采用粒子群优化算法对反应精馏塔的设计参数进行了优化。 2.3.2 软件集成平台的建立 在建立 Matlab 与 Aspen Plus 这两软件集成平台之前，首先在 Aspen Plus 中建立反应精馏塔的稳态模拟，并将模拟文件保存为 apwn 及 bkp 格式。然后在Matlab 中连接 Aspen 文件，具体步骤如下：（1） 在 Matlab 中建立一个 Aspen Plus COM 对象的句柄，语句如下： handles.aspen=actxserver(apwn.document); （2） 连接 Aspen Plus 相关文件，语句为： handles.filepathname=C:Filename.bkp; handles.aspen.InitFromFile2(handles.filepathname); handles.aspen.Visible=0; aspen=handles.aspen; 此时两软件集成平台已经建立完毕，可以在 Matlab 界面中访问及修改 Aspen Plus 中的数据。 Aspen Plus 采用树状结构来存储每个输入输出数据，用户可以通过 ToolsVariable Explorer中查看所需要操作的变量地址，如图所示，并在 Matlab 界面中 进行读取及赋值，语句如下： handles.aspen.Tree.FindNode(DateBlocksvariables).value 在对相关变量读取及赋值完成后，便可在 Matlab界面中运行，运行语句如 下： aspen.Run2 至此，便完成了两软件集成平台的建立。后续便可在 Matlab 界面编译优化 算法，通过软件集成平台来完成精馏塔的优化设计。 图 1-3 Aspen Plus 第二章 基于粒子群算法的反应精馏塔优化模型的建立 2.1 反应精馏塔稳态模型的建立 反应精馏塔包括精馏段、反应段和提馏段三部分。在 Aspen Plus 中，冷凝器 是第一块理论板，再沸器则作为最后一块理论板，传统的反应精馏塔如图 2-1 所示： 图 2-1 反应精馏塔通用结构图 可逆反应在第 II 部分的反应段进行，生成物由精馏段和提馏段进行分离，塔顶馏出轻组分,塔釜馏出重组分。 反应精馏塔的优化设计首先要保证产物纯度达到要求或满足其他约束条件， 然后再进行优化。优化模型选择以年操作费用(TAC)为目标，TAC 的计算方式如 下： （2-1）式中，操作费用 energy cost 为每年所用低压蒸汽、冷凝水的费用。设备费用capital cost 则包括反应精馏段填料、其余段填料、换热器、塔壳及其内件的费用。Payback period 指资本回收期，通常假设其数值为 3 年。一年工作8000小时计算。TAC 所涉及的相关参数的具体计算方式如下所示： 塔壳及其内件费用： （2-2） （2-3）式中，N-总板数； D-塔直径,m； H-塔高度,m； 反应精馏段填料费用： （2-4）式中，Nreact-反应段级数，0.5 表示填料等板高度 0.5m； 精馏段及提馏段填料费用： （2-5）换热器费用： （2-6）式中，Sr0.65 -再沸器的换热面积； Sc0.65 -冷凝器的换热面积； 设备费用： （2-7）操作费用： （2-8）式中，mLPStream -0.5MPa 低压蒸汽的量，吨/年； mColdwater-冷凝水的量，吨/年； 采用 Aspen Plus 进行稳态建模时，按照如下步骤进行： (1) 输入反应精馏体系中的所有组分； (2) 根据组分性质，选择适合的物性方法； (3) 输入进料流股温度、压力、流量及组成； (4) 设置反应精馏塔的压力、塔板数、回流比等参数； (5) 运行Aspen，将反应精馏塔塔顶冷凝器及塔釜再沸器用HeatX模块模拟， 进料流股分别为第二块板上升蒸汽及倒数第二块板下降液体，通过 EDR 模块，计算得到换热器的换热面积及公用工程消耗量； (6) 对反应精馏塔各操作参数进行单变量分析，得到初步优化结果及各操作 参数的变化范围。虽然对反应精馏塔进行单变量分析可以得到较优的操作参数及较低的 TAC， 但由于反应精馏塔的多操作参数及复杂的非线性，单变量分析得到的结果仍存在 较大的优化空间，因此仍需采用粒子群算法进行优化。 2.2 基于粒子群算法的反应精馏塔的优化设计 将目标函数的数学表达式表示如下： （2-8） （2-9）其中，Mi表示决策变量，在反应精馏塔的优化设计中，决策变量的选取需根据具体情况而定。一般需要对各个变量进行单变量的灵敏度分析，分析它们对 TAC 或产物纯度的影响程度，如果该变量对 TAC 或产物纯度的影响程度大，那么需要将其设定为决策变量，如果影响程度较小，则为了提高算法的效率，可以不将其作为决策变量。反应精馏塔优化的决策变量一般包括进料流量配比、总板数、回流比、进料位置、反应开始级、反应结束级、塔板持液量等参数，因为这些参数对目标函数的影响较大，然而最终决策变量的选取仍需依据单变量分析之后来决定。 式(3-10)则代表相应的约束条件，一般为产品纯度或转化率。由于在化工生 产过程中，需要规定其出料流股中某些物质的合格指标，以保证产品质量，因此在对反应精馏塔进行优化设计之前，必须使出料流股相应物质的纯度达到合格，反应转化率达到标准，再去寻求其他的优化。 因此，在适应度函数的选取方面，对于满足上述要求的 TAC 值，适应度函 数就取 TAC 的值，否则适应度函数取10000，即： 本文采用了粒子群算法来求解上述模型，综合与设计方法如图所示： 图 2-2 基于粒子群算法的反应精馏塔的综合与设计方法 在 MATLAB 环境下，粒子群算法的流程图如图所示：图 2-3 粒子群优化算法流程图 在对反应精馏塔进行优化设计的过程中，首先通过 Aspen Plus 建立稳态模型。再随机初始化一定规模大小的种群，种群中每个粒子都是决策变量的矢量。 通过 COM 技术实现对 Aspen Plus 的数据访问，将相关结构参数输送给 Aspen Plus 后进行仿真，仿真所得精馏塔设计结果返回 MATLAB。MATLAB 通过返回得来得数据对种群进行适应度即 TAC 的计算，由适应度的值来更新每个粒子的个体最优位置及全局最优位置，再对每个粒子的速度及位置进行进化，对边界条件处理后，产生新一代的个体，再次输送给 Aspen Plus 进行仿真，如此进行一系列的迭代，最终在粒子群算法的优化下，得出一最优解，即反应精馏塔的最优结构。 第三章 实例分析 3.1 EDA裂解反应精馏稳态模型描述及设计方案 碳四烃反应精馏稳态模型建立如图 3-1 所示：图 3-1 EDA裂解反应精馏稳态模型 采用 Aspen Plus 中的 RadFrac 模块来模拟EDA裂解反应精馏过程，物性方法选择NRTL。在该反应精馏体系中，发生的反应如下： 主反应： 反应（1）副反应：据文献，不同反应体系下的EDA选择性如下：表 3-1 不同体系下的EDA选择性选择与体系最相近的条件，反算得到副反应EDA水解的转化率为32.63%，采用R0401转化率反应器表示在反应精馏塔内发生反应（2）的程度。漆志文 漆志文,孙海军,沈才大,于遵宏.二醋酸亚乙酯裂解制醋酸乙烯动力学J.华东理工大学学报,1998(05):111-115.等人在管式反应器内进行动力学实验，在不同温度下，根据不同酸酐与EDA进料比进行实验（包含适宜的工业操作条件），由动力学实验数据，Runge-Kutta法积分，并用单纯行法对参数进行优化，得到了EDA裂解的动力学方程： 在此基础上进行化简可得： 该反应精馏塔初始化的模拟参数设置如表3-1 所示：表3-2 反应精馏塔物性参数和操作条件参数数值FEED流股流量802.5kmol/hFEED流股温度、压力135.9、1atmFEED流股组成（mol%）醋酸醋酐双醋酸亚乙酯9.83410-50.67740.3224操作压力1atm塔顶冷凝器全凝器总塔板数46塔釜采出量628.08回流比6FEED流股进料位置6反应段开始级6反应段结束级45反应段停留时间40sEDA转化率0.945塔釜醋酸浓度0.005EDA裂解反应精馏塔系统的目标函数表达式为：Min(TAC)=f(RR，BOT，N，FFEED，RStart，REnd) 决策变量初步选择为：回流比RR、塔釜采出量BOT、总塔板数N、FEED进料位置FFEED、反应段起始级RStart、反应段结束级REnd，共6个，其中前两个决策变量为连续变量，后四个离散变量。约束条件：塔釜醋酐醋酸摩尔分数YHOAC在0.01以下，EDA转化率大于0.94（以初始进料FEEF1为基准计算）。3.2 决策变量的影响及其变化范围 初步确定了EDA裂解反应精馏塔系统的决策变量后，还需要研究单一决策变量对系统的影响，由此判断某一决策变量是否对系统影响很大，如果影响不大，可不将其作为决策变量，以减少粒子维度，降低算法工作量。此外，对单一决策变量进行分析，可以确定决策变量的搜索范围，以缩小搜寻空间，加快收敛速度。以上述设置的稳态模型为初始模型，分别对上述6个决策变量进行单变量的分析，分析步骤为：将已优化变量更改后，再对下一变量进行分析及优选，以达到初步优化和确定最终决策变量及其范围的目的。 (1) 总板数的影响图 3-2 总板数对EDA转化率及塔釜醋酸浓度的影响 从图3-2可以看出，当塔板数增加时，塔釜醋酐浓度、EDA转化率都呈上升趋势。这是因为塔板数的多少决定了精馏系统的分离效果，总板数越多，分离效果越好，但是总板数达到一定数量后，EDA转化率和醋酐摩尔分数并不会显著增加，即再增加塔板数对精馏分离毫无意义，而且塔板数增多，会增加固定投资，导致 TAC 的值直线上升。因此在满足上述约束条件的前提下，优选塔板数为46。（2）反应段起始级图 3-3 反应段起始级对EDA转化率及塔釜醋酸浓度的影响 从图3-3可以看出，随着反应段起始级的增加，塔釜醋酐摩尔分数及EDA的转化率都呈缓慢下降趋势。这是在固定反应段结束级的前提下，反应段起始级增加，意味着反应段级数的减少，反应进行的越不充分，因此塔釜醋酐摩尔分数及EDA的转化率都降低，而且由于反应段级数减少，更贵的反应精馏填料被普通填料代替，设备费用会显著降低。因此初选的反应段起始级为8。（3）反应段结束级图 3-4 反应段结束级对EDA转化率及塔釜醋酸浓度的影响 从图3-4可以看出，随着反应段结束级的增加，EDA转化率及塔釜醋酐浓度都先上升后趋于稳定。反应段级数增加，反应进行的更充分，且更贵的反应精馏填料取代了普通填料所导致的，因此初选的反应段结束45。(4) FEED进料位置的影响 图3-5 FEED进料位置对EDA转化率及塔釜醋酸浓度的影响从图 3-5 可以看出，随着FEED进料位置的增加，EDA转化率先上升后降低，塔釜醋酸浓度基本满足先上升后降低，这是因为当FEED进料位置增加时，逐渐接近反应段，使得反应段中EDA的浓度增大，有利于反应，所以EDA转化率和会上升，醋酸摩尔分数会下降。而当FEED进料位置增加到一定程度后，再增加FEED进料位置，即远离了反应段，不利于原料中 EDA的充分反应，因此EDA的转化率会降低醋酸摩尔分数会升高。因此优选的进料位置为8。(5) 回流比的影响图3-6 回流比对EDA转化率及塔釜醋酸的影响 从图3-6可以看出，随着回流比的增大，塔釜醋酸摩尔分数减小，EDA转化率增大。这是因为回流比是影响产物纯度的主要因素，回流比越大，产品纯度越高，但是当回流比增大到一定程度后，回流比对醋酸浓度和EDA转化率的影响不大，可再沸器能耗会急剧上升，造成TAC值增大。因此，优选回流比为6。(6) 塔釜采出量的影响图3-7 塔釜采出量对EDA转化率及塔釜醋酸的影响 从图3-7可以看出，随着塔釜采出量的增大，塔釜醋酸摩尔分数增加，EDA转化率减小。当回流比增大到一定程度后，塔釜采出量对醋酸浓度和EDA转化率的影响较大。由于采出量过小会影响产量，因此，优选塔釜采出量为600kmol/h。因此经过初步的单变量分析及优化后，EDA裂解反应精馏塔的操作参数总结如3-3：表 3-3 初步优化后反应精馏塔操作参数 操作参数值总塔板数46反应段起始级8反应段结束级45FEED进料位置8回流比6塔釜采出量（kmol/h）600TAC（万元）1307.3但是上述分析只是单变量分析及优化，优化的结果并不是反应精馏塔的最优 设计结构，综合考虑这些变量，会表现出高度的非线性，也会存在较多局部最优 解，因此仍需要借助优化算法来确定其设计结构。 经过上述单变量分析，将上述 6 个决策变量的变化范围设定如表3-4所示：表 3-4 决策变量变化范围决策变量变化范围总塔板数40,60反应段起始级2,40反应段结束级20,60FEED进料位置1,30回流比5,10塔釜采出量500,8004.3 粒子群算法优化结果及对比 本例中，将粒子群算法的参数选择为：种群大小 N=100,最大进化代数 T=100，学习因子 c1=1.5,学习因子 c2=1.5，惯性权重最大值 Wmax=0.9，惯性权重最小值Wmin=0.4。在经过 100 代的进化之后，最优解与进化代数的关系如下图所示：图 3-8 最优解与进化代数关系 从图3-8中可以看出，随着进化代数的增加，最优解的TAC 值在进化初期波动比较剧烈，而在进化后期逐步稳定。这是因为进化后期，种群中大部分个体都接近全局最优值，得出的个体极值与全局极值相差不大，而且由于惯性权重随迭代次数的减少，粒子的速度也在逐渐降低，即越来越靠近全局最优解，因此进化后期各结果会趋于平稳。当进化代数为61代时，TAC 的值保持不变,为1065.2，可将得出的最优解近似为全局最优，此时最优解为5.1081，619.1197，51.0000，12.0000，9.0000，50.0000。算法优化后反应精馏塔的各项参数总结如表3-5所示： 表 3-5 算法优化后反应精馏塔操作参数操作参数值总塔板数51反应段起始级9反应段结束级50FEED进料位置12回流比5.1081塔釜采出量（kmol/h）619.1197将初步优化后的反应精馏塔操作参数和算法优化后的反应精馏塔的操作参 数在 Aspen Plus 软件中进行仿真运行，得出的结果进行对比如下表所示：表 3-6 初步优化结果与算法优化结果对比相关参数初步优化结果算法优化后结果总塔板数4651反应段板数3842再沸器热负荷（Gcal/h）21.312916.8043冷凝器热负荷（Gcal/h）-24.3684-19.8396EDA转化率0.95490.9471醋酸浓度0.003220.00548TAC（万元）1307.31065.2 从表中可以看出，本例子通过粒子群优化算法设计出来的反应精馏塔的总塔 板数及反应段板数虽然较多，但是冷凝器、再沸器的热负荷都比初步优化后的结果低，且在EDA转化率和醋酸摩尔分数达到要求的前提下，TAC降低了242.1万元，相当于初步优化结果的18.5%。附录 MATLAB源代码%粒子群算法求T0401烃反应精馏塔的经济最大值%粒子 x 的维度 D=6，分别为回流比、塔釜采出量、理论板数、进料位置、反应起始段位置、反应终止段位置%初始化clear; %清除所有变量close all; %清除所有图clc; %清屏tic%开始算法N=30; %种群粒子个数 ND=6; %每个粒子的维数即特征T=100; %最大迭代次数c1=1.5; %学习因子 1c2=1.5; %学习因子 2wmax=0.9; %惯性权重最大值wmin=0.4; %惯性权重最小值compute_time=zeros(T+1,2); %计算时间，算法+aspen%限制各属性位置范围% X1max=5; %回流比最大值X1min=10; %回流比最小值X2max=600; %采出量最大值X2min=800; %采出量最小值X3max=60; %理论板数最大值X3min=40; %理论板数最小值X4max=30; %废液进料位置最大值X4min=1; %废液进料位置最小值X5max=40; %反应起始段位置最大值X5min=2; %反应起始段位置最小值X6max=60; %反应终止段位置最大值X6min=20; %反应终止段位置最小值 %限制各属性速度范围%V1max=5; V1min=-5; V2max=5; V2min=-5; V3max=5; V3min=-5; V4max=5; V4min=-5; V5max=5; V5min=-5; V6max=5; V6min=-5; %构造位置最小值矩阵%Xmin=zeros(N,D); %N 行 D 列最小值矩阵Xmin(:,1)=X1min; %第 1 列属性的最小值Xmin(:,2)=X2min; %第 2 列属性的最小值Xmin(:,3)=X3min; %第 3 列属性的最小值Xmin(:,4)=X4min; %第 4 列属性的最小值Xmin(:,5)=X5min; %第 5 列属性的最小值Xmin(:,6)=X6min; %第 6 列属性的最小值%构造位置最大值向量%Xmax=X1max X2max X3max X4max X5max X6max;%构造各属性位置长度矩阵%detaX=ones(N,D);detaX(:,1)=X1max-X1min;detaX(:,2)=X2max-X2min;detaX(:,3)=X3max-X3min;detaX(:,4)=X4max-X4min;detaX(:,5)=X5max-X5min;detaX(:,6)=X6max-X6min;%构造速度最小值矩阵%Vmin=zeros(N,D); %N 行 D 列最小值矩阵Vmin(:,1)=V1min; %第 1 列属性的速度最小值Vmin(:,2)=V2min; %第 2 列属性的最小值Vmin(:,3)=V3min; %第 3 列属性的最小值Vmin(:,4)=V4min; %第 4 列属性的最小值Vmin(:,5)=V5min; %第 5 列属性的最小值Vmin(:,6)=V6min; %第 6 列属性的最小值%构造速度最大值向量%Vmax=V1max V2max V3max V4max V5max V6max;%构造各属性速度长度矩阵%detaV=ones(N,D);detaV(:,1)=V1max-V1min;detaV(:,2)=V2max-V2min;detaV(:,3)=V3max-V3min;detaV(:,4)=V4max-V4min;detaV(:,5)=V5max-V5min;detaV(:,6)=V6max-V6min;%初始化种群个体（限定位置和速度）%初始化位置%x=rand(N,D).*detaX+Xmin; %每一行代表粒子，列代表属性，x 控制在范围内x(:,3:6)=round(x(:,3:6); %对第 3-6 列进行四舍五入取整%初始化速度%v=rand(N,D).*detaV+Vmin; %速度初始化，控制在范围内v(:,3:6)=round(v(:,3:6); %对第 3-7 列进行四舍五入取整%判断粒子的设置是否合理%for i=1:N %判断总板数是不是最大，反应段起始位置是否小于终止位置while x(i,3)=max(x(i,3:6)|x(i,5)=x(i,6)|x(i,3)=x(i,6)x(i,:)=rand(1,D).*detaX(1,:)+Xmin(1,:); %更新不符合的值%对第 3列到第 6 列取整x(i,3:6)=round(x(i,3:6);end end%初始化个体最优位置和最优值%p=x; %初始个体最优位置%求初始个体最优值%pbest=ones(N,1);pbest=100000*pbest;Molefrac_HOAC=zeros(N,1);Molefrac_C4H6O3=zeros(N,1);Moleflow_topEDA=zeros(N,1);Moleflow_botEDA=zeros(N,1);Convert_EDA=zeros(N,1);suanfa_t0=toc; %初始化阶段，一部分算法部分计算时间tic%连接 aspen plushandles.aspen=actxserver(apwn.document);handles.filepathname=C:Users黄江辉DesktopT0401.bkp;handles.aspen.InitFromFile2(handles.filepathname);handles.aspen.Visible=0;aspen=handles.aspen;%求每个个体的 TAC 值for i=1:Naspen.Tree.FindNode(DataBlocksT0401InputBASIS_RR).value=x(i,1); %给摩尔回流比赋值aspen.Tree.FindNode(DataBlocksT0401InputBASIS_B).value=x(i,2); %给采出量赋值aspen.Tree.FindNode(DataBlocksT0401InputNSTAGE).value=x(i,3); %给总板数赋值aspen.Tree.FindNode(DataBlocksT0401InputFEED_STAGEFEED).value=x(i,4); %流股进料级%判断反应精馏塔的赋值方式，防止报错%读取 aspen 设置的值React_2=aspen.Tree.FindNode(DataBlocksT0401InputREAC_STAGE2#0).value; %反应终止段React_1=aspen.Tree.FindNode(DataBlocksT0401InputREAC_STAGE1#0).value; %反应起始段%设置 Aspen 的值if x(i,6)=React_1aspen.Tree.FindNode(DataBlocksT0401InputREAC_STAGE2#0).value=x(i,6); %反应终止段aspen.Tree.FindNode(DataBlocksT0401InputREAC_STAGE1#0).value=x(i,5); %反应起始段aspen.Tree.FindNode(DataBlocksT0401InputRES_STAGE2#0).value=x(i,6); %持液段终止aspen.Tree.FindNode(DataBlocksT0401InputRES_STAGE1#0).value=x(i,5); %持液段起始else aspen.Tree.FindNode(DataBlocksT0401InputREAC_STAGE1#0).value=x(i,5); %反应起始段aspen.Tree.FindNode(DataBlocksT0401InputREAC_STAGE2#0).value=x(i,6); %反应终止段aspen.Tree.FindNode(DataBlocksT0401InputRES_STAGE1#0).value=x(i,5); %持液段起始aspen.Tree.FindNode(DataBlocksT0401InputRES_STAGE2#0).value=x(i,6); %持液段终止end%运行 aspenaspen.Run2 %判断是否收敛if aspen.Tree.FindNode(DataBlocksT0401OutputPER_ERROR).value=0 %判断 aspen 是否收敛,收敛其值为 0pbest(i)=100000; %如果不收敛，则 TAC 值为 100000aspen.Close %关闭 aspen,再重新连接，防止由于上一次不收敛而影响下一个粒子handles.aspen=actxserver(apwn.document);handles.filepathname=C:Users黄江辉DesktopT0401.bkp;handles.aspen.InitFromFile2(handles.filepathname);handles.aspen.Visible=0;aspen=handles.aspen;elseDIS=aspen.Tree.FindNode(DataStreamsDISOutputMOLEFLMXMIXED).value;%DIS流率BOT=aspen.Tree.FindNode(DataStreamsDISOutputMOLEFLMXMIXED).value;%BOT流率Moleflow_topC4H6O3=aspen.Tree.FindNode(DataStreamsDISOutputMOLEFLOWMIXEDC4H6O3).value; %塔顶醋酐摩尔流率 Moleflow_botHOAC=aspen.Tree.FindNode(DataStreamsBOTOutputMOLEFLOWMIXEDHOAC).value; %塔釜醋酸摩尔流率 Molefrac_disC4H6O3(i)=Moleflow_topC4H6O3/DIS;%塔顶醋酐浓度Molefrac_botHOAC(i)=Moleflow_botHOAC/BOT;%塔釜醋酸浓度Moleflow_topEDA(i)=aspen.Tree.FindNode(DataStreamsDISOutputMOLEFLOWMIXEDEDA).value; %塔顶EDA摩尔流量Moleflow_botEDAH(i)=aspen.Tree.FindNode(DataStreamsBOTOutputMOLEFLOWMIXEDEDA).value; %塔釜EDA摩尔流量Moleflow_feedEDA(i)=aspen.Tree.FindNode(DataStreamsFEEDOutputMOLEFLOWMIXEDEDA).value; %进料EDA摩尔流量Convert_EDA(i)=1-(Moleflow_topEDA(i)+Moleflow_botEDAH(i)/258.76793124; %甲醇转化率if Convert_EDA(i)0.01pbest(i)=100000; %不达标，T