引言

为满足市场需求的多样化, 聚合物生产过程中需频繁进行牌号切换操作，其过渡过程不可避免会减少有效生产时间，产生大量过渡废料，造成原材料的损耗以及生产效率的降低。因此聚合物牌号切换优化研究对聚合工业生产具有很高的实用价值。

目前国内外关于牌号切换的研究主要集中在聚合反应过程分析^[1-2]、牌号切换模型建立^[3-4]及优化目标函数构造^[5]等方面。求解聚合反应牌号切换的动态优化问题主要有控制变量参数化方法^[6]、正交配置法^[7]及迭代动态规划^[8]等。传统的优化算法在求解牌号切换问题时，易受问题限制，具有局限性。进化算法由于具有鲁棒性高、适应性强、不受问题特性制约且无需先验知识的优点，在牌号切换优化问题的求解上得到了越来越多的应用^[9]，但进化算法多采用加权法^[10]将多目标优化问题转换为单目标问题求解，若权值选取不当会导致优化效率降低、效果变差。

本课题组前期以实验室规模的连续搅拌釜式反应器(CSTR)中苯乙烯聚合反应过程的牌号切换操作为对象，建立了多目标动态优化模型^[11]。此优化模型可以使聚苯乙烯的质量指标在牌号切换操作结束后接近目标牌号值，但无法达到目标值，反应过程需要在目标牌号的操作条件下运行一段时间，才能使聚合产品的质量指标满足目标牌号的要求。因此，本文在文献[11]基础上将基于Pareto支配的多目标优化算法与约束处理相结合，提出了一种基于策略协同的约束多目标粒子群优化(CMOPSO)算法求解牌号切换多目标优化问题，使聚合物质量指标在切换操作结束时满足目标牌号要求的同时，兼顾牌号切换过渡时间及过渡用料的优化。通过采取能够提供多样性解集的多目标优化算法以及具有针对性的约束处理的策略协同方法，来获得更好的聚合反应牌号切换优化结果。

1 牌号切换过程动态优化模型 1.1 苯乙烯连续聚合反应过程动态模型

如图 1所示为实验室规模的CSTR中苯乙烯溶液聚合反应过程示意图。本文假设单体流量F_m以及反应釜内反应物的体积V不变。单体和引发剂分别与溶剂充分混合后进入反应釜中，反应釜的温度通过冷却水调节。

苯乙烯聚合反应是典型的自由基聚合反应，包括链引发、链增长、链转移和链终止等基元反应。用$ {C_{\rm{R}}} = \sum\limits_{j = 1}^\infty {{C_{{{\rm{R}}_j}}}} $表示活聚物的总浓度，$ {C_{\rm{P}}} = \sum\limits_{j = 2}^\infty {\;{C_{{{\rm{P}}_j}}}} $表示死聚物的浓度，其中C_Rj和C_Pj分别代表链长为j的两种聚合物浓度。各反应组分的浓度满足式(1) 所示的物料守恒方程。

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}{C_{\rm{I}}}}}{{{\rm{d}}t}}\left( {C_{\rm{I}}^0 - {C_{\rm{I}}}} \right)/\theta - {K_{\rm{d}}}{C_{\rm{I}}}\\ \begin{array}{*{20}{c}} {\frac{{{\rm{d}}{C_{\rm{M}}}}}{{{\rm{d}}t}} = \left( {C_{\rm{M}}^0 - {C_{\rm{M}}}} \right)/\theta - 2{K_{\rm{i}}}{C_{\rm{I}}} - }\\ {\left( {{K_{\rm{p}}} + {K_{{\rm{trm}}}}} \right){C_{\rm{M}}}{C_{\rm{R}}}} \end{array}\\ \frac{{{\rm{d}}{C_{\rm{R}}}}}{{{\rm{d}}t}} = - {C_{\rm{R}}}/\theta + 2{K_{\rm{i}}}{C_{\rm{I}}} - {K_{{\rm{tc}}}}C_{\rm{R}}^2\\ \frac{{{\rm{d}}{C_{\rm{P}}}}}{{{\rm{d}}t}} = - {C_{\rm{P}}}/\theta + {K_{{\rm{trm}}}}{C_{\rm{M}}}{C_{\rm{R}}} + \frac{{{K_{{\rm{tc}}}}}}{2}{C_{\rm{R}}} \end{array} \right. $

(1)

其中，θ=V/F为CSTR的平均停留时间；F=F_i+F_m为反应器进料的总流量，V为反应器体积；C_I、C_M分别为引发剂、单体的浓度；K_d、K_i、K_p、K_trm和K_tc是各基元反应的反应速率常数，满足Arrhenius方程K=Aexp(-E/R′T)，R′是气体常数。

反应器和夹套的温度满足热平衡方程

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}T}}{{{\rm{d}}t}} = \frac{F}{V}\left( {{T_{\rm{f}}} - T} \right) + \frac{{ - \Delta {H_{\rm{r}}}}}{{\rho Cp}}{K_{\rm{p}}}{C_{\rm{M}}}{C_{\rm{R}}} - \\ \;\;\;\;\;\;\;\;\;\frac{{{U_{\rm{r}}}{A_{\rm{r}}}}}{{\rho CpV}}\left( {T - {T_{{\rm{jt}}}}} \right)\\ \frac{{{\rm{d}}{T_{{\rm{jt}}}}}}{{{\rm{d}}t}} = \frac{{{F_{{\rm{jt}}}}}}{{{V_{{\rm{jt}}}}}}\left( {{T_{{\rm{jf}}}} - {T_{{\rm{jt}}}}} \right) + \frac{{{U_{\rm{r}}}{A_{\rm{r}}}}}{{{\rho _{{\rm{jt}}}}C{p_{{\rm{jt}}}}{V_{{\rm{jt}}}}}}\left( {T - {T_{{\rm{jt}}}}} \right) \end{array} \right. $

(2)

其中，ΔH_r、U_r和A_r分别表示反应热、传热系数和传热面积，Cp和Cp_jt分别表示反应混合物和冷却水的比热容，ρ和ρ_jt分别表示冷却水密度及反应釜内混合物的密度，T和T_f分别表示反应器和反应物进料的温度，T_jf和T_jt分别表示冷却剂进料和夹套的温度，V_jt是反应器夹套的体积，F_jt为冷却剂流量。

用矩来推算聚合物的分子量分布参数，包括数均分子量(M_n)、重均分子量(M_w)和分散指数(D_i)。死聚物分子量分布的m阶矩定义为

$ {Z_m} = \sum\limits_{j = 2}^{ + \infty } {{j^m}{P_j}} \;\;\;\;\left( {m = 0,1,2, \cdots } \right) $

(3)

分子量分布参数M_n、M_w和D_i分别由式(4)~(6) 计算

$ {M_{\rm{n}}} = {M_{\rm{m}}}{Z_1}/{Z_0} $

(4)

$ {M_{\rm{w}}} = {M_{\rm{m}}}{Z_2}/{Z_1} $

(5)

$ {D_{\rm{i}}} = {M_{\rm{w}}}/{M_{\rm{n}}} = \frac{{{Z_0}{Z_2}}}{{Z_1^2}} $

(6)

其中M_m是单体分子量。

1.2 基于策略协同的牌号切换优化

苯乙烯聚合反应牌号切换过程优化的目的是确保质量指标在过渡过程结束时稳定在目标牌号值的前提下，使过渡时间尽可能短以及原材料损失尽可能少。该问题具有多极小、多目标的特点，且牌号切换过程中质量指标与成本指标存在矛盾^[12]。如采用简单的加权法将多目标转换为单目标求解，很可能会因为权值选取不当而影响优化项的作用，甚至淹没相应的优化项，且仅取得单一的优化结果，具有局限性。

本文采用策略协同优化，通过多元集成、协同优化的方式提高优化质量及效率^[13]，构造了如式(7)~ (9) 所示的苯乙烯聚合反应过程的牌号切换多目标优化模型

$ \left\{ \begin{array}{l} {f_1} = \min \int_0^{{t_{\rm{p}}}} {\left[ {{w_1}{{\left( {\frac{{{M_{\rm{n}}} - M_{\rm{n}}^ * }}{{M_{\rm{n}}^ * }}} \right)}^2} + } \right.} \\ \;\;\;\;\;\;\left. {{w_2}{{\left( {\frac{{{M_{\rm{i}}} - M_{\rm{i}}^ * }}{{M_{\rm{i}}^ * }}} \right)}^2}} \right]{\rm{d}}t\\ {f_2} = \min \int_0^{{t_{\rm{p}}}} {{F_{\rm{i}}}\left( t \right){\rm{d}}t} \\ {\rm{s}}{\rm{.t}}x\left( t \right) = f\left( {x\left( t \right),u\left( t \right)} \right),x\left( {{t_0}} \right) = {x_0}\\ y\left( t \right) = g\left( {x\left( t \right),u\left( t \right)} \right)\\ {x_{\min }} \le x\left( t \right) \le {x_{\max }}\\ {y_{\min }} \le y\left( t \right) \le {y_{\max }}\\ 0.001 \le {F_{\rm{i}}} \le 0.04\\ 0.001 \le {F_{{\rm{jt}}}} \le 0.02\\ 340 \le T \le 370 \end{array} \right. $

(7)

$ \frac{{\left| {{M_{\rm{n}}}\left( {{t_{\rm{p}}}} \right) - M_{\rm{n}}^ * } \right|}}{{M_{\rm{n}}^ * }} \le 1\% $

(8)

$ \frac{{\left| {{M_{\rm{w}}}\left( {{t_{\rm{p}}}} \right) - M_{\rm{w}}^ * } \right|}}{{M_{\rm{w}}^ * }} \le 1\% $

(9)

其中，t_p为优化计算时间，M_n^*、D_i^*、F_i^*分别为目标牌号时的数均分子量、分散性系数及引发剂流量，w₁、w₂为f₁的各项优化权值，优化控制变量为引发剂流量F_i以及冷却水流量F_jt。出于系统安全的考虑，将反应器的温度上限值设定为370 K；为避免单体转换率过低，反应器温度的下限值设定为340 K。f₁和f₂分别为反映质量和引发剂消耗量的性能指标。

在多目标处理策略上，直接采用基于Pareto支配的多目标优化算法求解此常微分方程组，在实现牌号切换优化操作的同时，得到一组牌号切换过程的优化操作变量值。同时本文引入约束策略，在目标函数中增加如式(8)、(9) 所示的终点约束，优先选择满足约束条件的粒子放入Pareto非支配解集，从而在维护粒子多样性的同时提高搜索效率。

2 CMOPSO算法 2.1 多目标粒子群算法及约束处理机制

多目标粒子群(MOPSO)算法作为一种基于群体智能的进化算法，将多目标优化问题看作一个整体，并根据Pareto支配概念择优，定义目标最小值如式(10) 所示

$ \left\{ \begin{array}{l} \forall i \in \left\{ {1, \cdots ,m} \right\},{f_i}\left( u \right) \le {f_i}\left( v \right)\\ \wedge \exists j \in \left\{ {1, \cdots ,m} \right\},{f_j}\left( u \right) < {f_j}\left( v \right) \end{array} \right. $

(10)

当且仅当式(10) 成立时，解u=(u₁, u₂, …, u_n)支配解v=(v₁, v₂, …, v_n)，记做u≺v。f(x)定义了由决策空间向目标空间映射的函数。对于多维决策变量，采用式(11) 进行种群初始化：

$ X_{ij}^0 = {L_j} + r\left( {{U_j} - {L_j}} \right)\;\;\;j = 1,2, \cdots ,D $

(11)

其中，解X_ij⁰的取值范围为[L_j, U_j]，D为决策变量的维度，r为均匀分布于[0, 1]区间的随机数。

计算初始种群的多个适应度值，根据支配关系确定初始个体最优解b_p，将每次迭代的非支配集存储到外部归档集(EA)中，并采用拥挤距离法对EA进行维护。设有两个子目标f₁和f₂，个体i的拥挤距离d_i是该个体相邻其他两个个体所围成的矩形的长宽之和。P_{i, f}为个体i在子目标f上的函数值。则d_i的表达式为

$ {d_i} = \left| {{P_{i + 1,{f_1}}} - {P_{i - 1,{f_1}}}} \right| + \left| {{P_{i + 1,{f_2}}} - {P_{i - 1,{f_2}}}} \right| $

(12)

EA中d_i小的粒子被d_i大的粒子代替，非支配解代替支配解。将EA中的粒子按其d_i值降序排列，从前10%粒子中随机选择一个粒子作为全局最优解b_g。

在粒子搜索范围内，粒子速度和位置分别按式(13)、(14) 更新

$ V_{ij}^{I + 1} = wV_{ij}^I + {c_1}{r_1}\left( {b_{{\rm{p}}\left( j \right)}^I - X_{ij}^I} \right) + {c_2}{r_2}\left( {b_{{\rm{g}}\left( j \right)}^I - X_{ij}^I} \right) $

(13)

$ X_{ij}^{I + 1} = X_{ij}^I + V_{ij}^{I + 1} $

(14)

其中，w为惯性权重，I为当前迭代次数，V_ij为粒子的速度，c₁和c₂表示粒子的个体学习能力和社会学习能力，r₁和r₂是均匀分布于[0, 1]区间的随机数。

为提高MOPSO的全局搜索能力，引入变异机制随机选择EA中的粒子并对其位置进行变异，变异概率p一般为[0, 1]区间内一个较小的数。对于粒子X，随机选择一维按式(15) 进行变异：

$ {{X'}_j} = {X_j} + \left( {{U_j} - {L_j}} \right)r $

(15)

其中，X′_j代表变异后的粒子，U_j、L_j表示粒子第j维上界和下界，r为均匀分布于[-1, 1]区间的随机数。若X′_j超出上界，则X′_j=U_j；若X′_j超出下界，则X′_j=L_j。

计算新粒子的多个适应度值，并进入下一次迭代，直至I满足最大迭代次数I_max。

在MOPSO算法中，根据Pareto支配概念及约束处理机制定义粒子X的过程约束违反程度如式(16) 所示^[14]。

$ Q\left( X \right) = \sum\limits_{i = 1}^m {\max \left( {0,{g_i}\left( X \right)} \right)} + \sum\limits_{i = 1}^m {\max \left( {0,\left| {{h_i}\left( X \right)} \right|} \right)} $

(16)

其中，g_i(X), i=1, 2, …m为m个不等式约束条件，h_i(X), i=1, 2, …n为n个等式约束条件。若粒子可行，则有Q(X)=0。

2.2 控制变量参数化方法及改进

由于MOPSO算法本身无法处理动态优化问题，采用控制变量参数化(CVP)方法将动态问题转化为非线性规划(NLP)问题，该方法在时域上将连续的控制变量用一系列分段基函数近似，从而将难以描述的动态问题转化为参数优化问题。选择不同形式的基函数可得到控制变量的多种逼近策略，本文采用基本的分段常量逼近策略，将各控制变量平均分成10段。采用龙格-库塔法求解其DAE方程，然后将求得的目标函数值送入MOPSO中进行优化。

由于状态变量没有直接包含在NLP中，CVP方法不能直接处理状态变量带约束问题。为此本文在MOPSO约束违反程度中增加一项对状态变量在切换终点处的约束违反惩罚项。对于式(8)~(9) 所示的终点约束，有如下的惩罚项：

$ {Q_{{\rm{vio}}}}\left( x \right) = \sum\limits_{i = 1}^L {\max \left( {0,{M_i}\left( X \right)} \right)} $

(17)

其中，M_i(x), i=1, 2, …L为L个定义的惩罚函数，本文中取L=2。当各状态变量均满足终点约束时，若各惩罚函数值为零，则惩罚项Q_vio(x)为零；否则Q_vio(x)为无穷大。

考虑终点约束的Pareto解集更新策略如下：

1) 若粒子满足终点约束且粒子可行，则根据粒子适应度值Pareto支配关系择优；

2) 若粒子满足终点约束但粒子不可行，则优先选择粒子过程约束违反程度小的粒子；

3) 若粒子不满足终点约束，由于粒子约束违反程度较大而被淘汰。

改进的CMOPSO流程图如图 2所示。

3 结果与讨论

本文以苯乙烯聚合反应过程中牌号A切换到牌号B的过渡过程为研究对象，采用改进的CMOPSO算法求解基于策略协同的约束多目标优化方案，并将文献[10]中单目标加权优化策略的结果作为方案一进行比较。牌号A、B的质量指标及稳态操作条件如表 1所示，目标牌号的稳态操作条件由稳态优化获得^[15]。

令t_P=300 min^[9]，各参数设置为：种群规模N=120，I_max=100，c₁、c₂=1；Pareto最优集规模M=20；w从0.95至0.4线性递减；p=0.1；w₁=0.4，w₂=0.6。

由于优化目标的矛盾性，过渡时间与过渡用料不可能同时达到最优，而CMOPSO搜索得到的Pareto解集包含多个非劣解，如图 3所示。

由图 3看出，多目标优化方案中的Pareto解集随着代表质量指标的f₁逐渐减小，代表经济指标的f₂逐渐增大，因此需要根据具体情况选择合适的切换方案。本文从解集中较优的质量指标和经济指标中各选取若干点分别作为方案二及方案三，根据切换操作策略进行仿真，并与方案一的仿真结果进行比较。

经过优化计算得到的操作变量F_i、F_jt的优化轨迹如图 4(a)、(b)所示。0≤t≤t_P时的牌号切换操作为动态优化阶段，操作变量的轨迹由PSO或CMOPSO计算求得；t≥t_p时过程操作转换到准稳态过渡阶段，以目标牌号的操作条件F_i^*、T^*为输入量。将所得操作变量优化轨迹进行仿真，得到状态变量M_n、M_w、D_i的优化轨迹如图 4(c)、(d)、(e)所示，以M_n满足|M_n^*-M_n|≤30作为牌号切换过渡阶段结束的标志。

为了对优化结果进行评价，将3种方案的优化结果与阶跃切换进行了比较。阶跃切换方案以牌号A的稳态为初始状态，使操作条件F_i和T阶跃切换至牌号B的稳态操作值，仿真结果如图 4中虚线所示。该阶跃方案是一种较为理想的情况，没有考虑到热量衡算问题。

由图 4(a)、(b)可知，3种方案在前300 min各有10段优化值，300 min后F_i和F_jt均为牌号B的稳态操作值。由图 4(c)、(d)、(e)可知，3种方案都能较好地完成聚合反应的牌号切换过程，且优化效果均明显优于阶跃切换。至动态优化阶段结束时，方案一的质量指标接近但略小于目标牌号稳态值，方案二的质量指标满足终点约束，与目标牌号质量指标之间的偏差不大于1%，显然优于方案一；方案三的质量指标同样满足终点约束，但与方案二相比偏向于更优的经济指标，且其质量指标的超调相对较大，准稳态过渡时间也相对较长；3种方案相比，方案二的牌号切换过渡过程也更加快速、平稳。

由于牌号切换过渡时间和过渡用料的矛盾关系，方案一的优化结果严重依赖于权值的选择，如权值选取不当可能导致优化过程中较优解的丢失。采用兼顾二者的CMOPSO算法产生一组非劣解，可从中选择符合需求的切换操作策略，所以最终证明方案二最优。

表 2为4种牌号切换方案下的仿真结果。显然本文提出的3种方案在过渡时间及原材料消耗量上均优于阶跃切换方案。方案一较之阶跃切换，牌号切换时间缩短了420 min，引发剂总消耗量和单体的消耗量分别减少了4.8 L和8.4 L，优化过程引发剂消耗量为3.12 L；方案二与方案一相比，切换时间缩短了105 min，引发剂总消耗量和单体消耗量分别减少1.3 L和2.1 L；方案三与阶跃切换相比，切换时间缩短了160 min，引发剂总消耗量和单体消耗量分别减少2.35 L和3.2 L。需要指出的是，方案三优化过程中引发剂消耗量为2.71 L，优于其他3种方案，但其过程的平稳性不及方案一及方案二，导致过程参数及质量指标同时到达目标稳态的总过渡时间较长，引发剂和单体消耗量也随之增多。

4 结束语

本文针对苯乙烯聚合反应牌号切换过渡过程，提出一种基于策略协同的牌号切换优化方法，引入多目标优化策略，有效避免了加权优化方法的局限性，缩短了过渡时间，减少了过渡用料。针对质量指标的要求，引入约束策略，并对常规的CVP方法进行了改进，将状态变量约束转化为约束违反惩罚项。利用CMOPSO寻优得到一组牌号切换方案，决策者可根据自身需求选择最适合的方案，以满足企业在实际生产过程中多样化的生产需求，提高企业的竞争力。