1 前言

在多活性位Ziegler-Natta催化剂的作用下，由于每一个活性位性质不同使得生成的聚烯烃具有较宽的分子量分布^[1-2](molecular weight distribution，MWD)，而聚合物的性能，如韧性、硬度、弹性等与MWD有着密切的关系。工业生产中，需要调控聚合物的MWD来生产不同性能的聚合物。

凝胶渗透色谱法(gel permeation chromatography，GPC)可以准确分析聚合物的MWD、平均分子量和多分散指数(polymer dispersity index，PDI)。将GPC测得的MWD数据去卷积成多个Flory分布曲线，可以有效确定催化剂活性位个数、每个活性位的质量分数以及聚合物微观结构信息。目前，有大量方法用于MWD去卷积，如非线性最小二乘法^[3]、Levenberg-Marquardt方法^[4](L-M)、广义梯度下降法^[5](generalized reduced gradient，GRG2)以及一些商用软件(Peakfit^[6]、Microsoft Excel Solver^[7])等。但是这些传统非线性优化算法需要提供精确的初始值，若初始值偏离全局最优解，算法易陷入局部最优。随着智能优化算法的发展，将智能算法应用到MWD去卷积中可以克服传统优化算法的不足。NANTHAPOOLSUB等^[5]使用遗传算法将MWD与共聚组成分布(copolymer composition distribution，CCD)同时去卷积分析，对比GRG2算法，发现使用GRG2算法同时对MWD和CCD进行分析会导致估计参数与实际值存在较大误差；SINGH等^[8]使用进化算法对GPC数据进行去卷积分析估计了催化剂活性位个数；安许锋等^[9]等使用改进飞蛾火焰算法，以数均分子量和重均分子量为目标估计了乙烯聚合过程中的动力学参数。

免疫算法(immune algorithm，IA)是一种新型的仿生智能算法，其根据生物免疫系统原理使得算法具有较强的自组织、自调节和并行性等特点。在解决复杂优化问题时，IA算法具有精度高，收敛速度快，不易陷入局部最优解的优点，在多目标优化^[10]、输配电系统^[11]、车间调度^[12]以及复杂网络结构优化^[13]等领域获得广泛应用。

本文针对聚烯烃MWD去卷积分析过程中采用L-M算法需提供精准的初始值否则无法准确获得每个活性位动力学参数的问题，提出基于全局搜索的IA算法。利用IA算法随机产生初始值，根据模型计算值与GPC采样数据的误差平方和估计催化剂活性位的个数以及每个活性位占的比例。最后将IA算法与L-M算法结合，增强IA算法的局部寻优能力，提高去卷积分析的精度。

2 算法原理 2.1 分子量分布去卷积

由非均相Ziegler-Natta多活性位催化剂制得的聚合物一般具有较宽的分子量分布，其中每一个活性位生成的聚合物的分子链长分布(chain length distribution，CLD)服从Flory-Schulz最可几分布^[14]：

$ {f_r} = r{\tau ^2}{e^{ - r\tau }} $

(1)

式中：$\tau $为链转移速率和链增长速率之比，r为聚合度(即聚合物的链长)。

通常根据M_p=M_pm·r(M_pm为聚合单体的摩尔质量，g·mol^-1)将CLD分布转换为MWD分布，得：

${w_{{M_P}}} = {M_{\rm{p}}}{\hat \tau ^2}{e^{ - {M_{\rm{p}}}\hat \tau }}$

(2)

式中：$\hat \tau = \frac{1}{{{M_n}}} = \frac{2}{{{M_w}}}$，M_n为数均分子量，M_w为重均分子量，M_p为聚合物分子量，g·mol^-1。

多活性位催化剂生成的聚烯烃MWD为单活性中心结果的加权求和，在实验分析过中，一般对分子量M_p取对数形式，则最终聚合物MWD形式如下^[15]：

$\frac{{dW}}{{d\log {M_{\rm{p}}}}} = \sum\nolimits_i^{NS} {\ln 10 \cdot {m_i}} {(\frac{1}{{{M_{n, i}}}})^2}\exp [2\ln 10 \cdot \log ({M_{\rm{p}}}) - \frac{1}{{{M_{n, i}}}}\exp (\ln 10 \cdot \log {M_{\rm{p}}})]$

(3)

$\sum\nolimits_i^{NS} {{m_i}} = 1$

(4)

式中：W表示聚合物质量分率，NS表示催化剂活性位个数，${m_i}$为第i个活性位产生的聚合物所占比重，${M_{n, i}}$为第i个活性中心生成聚合物的数均分子量，g·mol^-1。

分子量分布去卷积就是将MWD曲线拟合成多个单活性中心的MWD，目标方程如下：

${x^2} = \sum\nolimits_k^n {(w_{GPC}^k} - \sum\nolimits_i^{NS} {w_{{\rm{model, }}i}^k} {)^2}$

(5)

式中：n为MWD的GPC采样数据个数，$w_{GPC}^k$为第k个采样点的值，$w_{{\rm{model, }}i}^k$为第i个活性位在第k采样点处的值。

根据SOARES等^[4]的观点，逐渐增加活性位个数，对式(5)进行拟合，若拟合的残差平方和无显著地减少，则可以认为当前的活性位个数已经足够描述聚合物MWD。通过去卷积分析，可以确定催化剂活性位个数、每个活性位所生成聚合物的质量分数以及数均分子量。

2.2 免疫算法

IA算法是受生物免疫系统启发而发展起来的一种新兴智能优化算法，通过对抗体-抗体和抗体-抗原之间的亲和度进行评价，保证了抗体种群的多样性。此外，采用克隆选择、免疫记忆以及对产生的抗体进行促进或者抑制，实现了算法的自我调节功能，克服了一般算法易陷于“早熟”的缺点。

由式(3)、(4)和(5)可知，MWD去卷积分析本质上属于非线性优化问题。现有研究多以最小二乘、梯度下降等传统算法对该问题进行求解，但此类算法对初始值依赖性较高，易陷入局部最优。IA算法在求解复杂非线性问题时，通过随机生成初始值，借助亲和度算子、克隆算子、变异算子和选择算子等一系列免疫操作算子求解优化问题，实现全局寻优，可有效解决传统算法对初始值依赖性较高的问题且求解精度显著提高。

免疫系统概念与IA算法的概念关系如表 1。因此，对MWD去卷积分析时，目标方程(5)为IA算法的抗原，参数X=[${m_i}$, ${M_{n, i}}$]为IA算法的抗体。通过对抗体X进行一系列的免疫操作，求得目标方程(5)的最优解。具体的算法操作步骤如下：

(1) 初始种群  随机生成N个初始抗体种群。

(2) 评价  抗体-抗原和抗体-抗体之间的亲和度评价、抗体浓度及抗体激励度。

(3) 克隆  选取激励度前Nc的个体进行克隆作为记忆库。

(4) 变异  对记忆库中的克隆抗体进行变异操作。

(5) 选择  亲和度高的变异抗体替代亲和度低的克隆抗体，同时保留当代种群中亲和度最大的抗体并更新记忆库。

(6) 种群刷新  随机产生N-Nc个抗体与记忆库中的Nc个抗体形成新的种群。

(7) 迭代终止  判断是否满足终止条件，满足则输出结果，否则转步骤2。

免疫算法流程如图 1所示，上述步骤中的免疫操作算子如下：

1) 种群初始化

初始产生的种群主要采用随机方式产生，定义如下：

$X = L + U + {\rm{rand}} \cdot (U - L)$

(6)

式中：L和U分别为可行解的下界和上界，rand为0~1的随机数。

2) 抗体-抗原亲和度算子

抗体与抗原之间的亲和度通常表示可行解对目标问题的识别程度，由于求解问题为最小化问题，因此抗体与抗原之间的亲和度定义为目标方程的倒数：

$A(X) = {\frac{{\rm{1}}}{{{x^2}{\rm{|}}}}_X}$

(7)

3) 抗体-抗体亲和度(相似度)

抗体与抗体之间的亲和度主要表示2个抗体之间的相似程度，计算公式如下：

$S({X_{\rm{i}}}, {X_j}) = \sqrt {{{({X_i} - {X_j})}^2}} $

(8)

式中：$ {X}_{i} $为第i个抗体。

4) 抗体浓度

抗体的浓度就是指抗体种群中相似抗体所占据的比重，定义如下：

$D({X_i}) = \frac{1}{N}\sum\nolimits_{j = 0}^{N - 1} {S({X_i}, {X_j})} $

(9)

式中：N为抗体总数；，T为相似度阈值。

5) 抗体激励度

每个抗体激励度有抗体和抗原之间的亲和度与抗体浓度构成，定义如下：

$P = \alpha \frac{{A({X_{\rm{i}}})}}{{\sum\nolimits_{j = 0}^N {A({X_j})} }} + (1 - \alpha )\frac{{D({X_i})}}{{\sum\nolimits_{j = 0}^N {D({X_j})} }}$

(10)

式中：α为常数。显然，抗体与抗原的亲和度越高、抗体浓度越低则激励度越高，越容易受到激励而产生繁殖。这样鼓励亲和度高的抗体，同时抑制浓度高的抗体，可以保证抗体的多样性。

3 实验结果与分析

非线性最小二乘算法进行去卷积分析时，初始值的设定根据文献[4]的经验公式(11)。

当初始拟合的活性位个数NS=2时，需要给定的初始值为${m_1}$、${m_2}$、${M_{n, 1}}$和${M_{n, 2}}$。

$\left\{ \begin{gathered} \frac{{{M_n}}}{{{M_{{\rm{pm}}}}}} = \frac{1}{{{m_1}{\tau _1} + {m_2}{\tau _2}}} \\ \frac{{{M_w}}}{{{M_{{\rm{pm}}}}}} = 2(\frac{{{m_1}}}{{{\tau _1}}} + \frac{{{m_2}}}{{{\tau _2}}}) \\ \frac{{{M_z}}}{{{M_{{\rm{pm}}}}}} = 3(\frac{{{m_1}}}{{\tau _1^2}} + \frac{{{m_2}}}{{\tau _2^2}})(\frac{{{m_1}}}{{{\tau _1}}} + \frac{{{m_2}}}{{{\tau _2}}}) \\ {m_1} + {m_2} = 1 \\ \end{gathered} \right.$

(11)

其中，${\tau _i} = \frac{1}{{{M_{n, i}}}}$，${M_n}$、${M_w}$、${M_z}$分别为数均分子量、重均分子量和质均分子量。

当NS > 2时，令${\tau _{NS + 1}} = \sum\nolimits_i^{NS} {{m_i}{\tau _i}} $，${m_{NS + 1}} = \frac{1}{{NS}}\sum\nolimits_i^{NS} {{m_i}} $，并对新的${m_i}$进行归一化，使得$\sum\nolimits_i^{NS + 1} {{m_i}} = 1$。

以某石化公司生产的聚乙烯为例，对GPC测试得到的n组聚乙烯MWD数据进行去卷积分析。下文图中仅列出代表性的2、4、5活性位的分峰曲线，分峰结果数据均列在表格中。首先，采用L-M算法对GPC曲线去卷积，初始值根据式(11)计算，结果如表 2所示。

图 2为采用L-M算法对不同活性位个数的去卷积结果。图 2(a)表明2个活性位不能准确描述MWD；NS增加到4，拟合效果明显变优(图 2(b))；当NS为5时，图 2(c)的拟合效果相较于4个活性位并无显著改善。去卷积分析的误差平方和x²列在表 3。由表 3可知，当活性位个数为4时，x²为0.0171，增加活性位个数到5和6，x²反而还增大，则说明使用L-M算法时，4个活性位即可描述聚合物的MWD。

分析表 2中L-M算法初始值发现，由于${\tau _{NS + 1}} = \sum\nolimits_i^{NS} {{m_i}{\tau _i}} $，随着活性位的个数NS逐渐增大，${\tau _{NS + 1}}$的变化却越来越小，即数均分子量${M_{n, i}}$(${M_{n, i}} = \frac{1}{{{\tau _{NS + 1}}}}$)基本上维持在5 923，而${M_{n, i}}$对应的活性位的质量分数${m_i}$一直在减小，导致初始值的设置不够准确，使得L-M算法陷入局部最优解。

针对初始值设置问题，提出了IA算法对MWD进行去卷积分析。IA算法通过在解空间产生随机解，根据启发式搜索方式，结合克隆、变异和选择操作算子，可以有效搜索到全局最优解。图 3为IA算法得到的不同活性位个数的拟合结果，GPC数据与模型的误差平方和x²也列在了表 3。

随着活性位个数增加，GPC曲线和去卷积分析结果也随之变化。由图 3可知，当活性位个数为4时，模型计算的曲线与GPC曲线基本吻合，与L-M算法一致；当活性位个数增加到5时，IA算法的拟合效果进一步提升，而L-M算法却出现下降，原因在于L-M算法的初始值设置不够准确，使得求解结果陷入局部最优。从表 3可知，活性位个数从4增加到5，GPC数据与模型计算的误差平方和x²从0.017 4减小到0.008 7，而当活性位个数为6时，其误差平方和反而增加。因此，采用IA算法，5个活性位足以描述聚合物的MWD。如图 3所示，按照活性位所生成聚合物分子量从小到大的顺序，以数字1~5标识。聚合反应过程中，标识为1的活性位所生成聚合物的数均分子量最小，表明该活性位的链增长速率与链转移速率之比最小。该活性中心具有倾向于生成短链聚合物、共聚能力强、易失活的特性。标识为2~5的活性位生成的聚合物分子量依次增大，共聚能力减弱，活性越来越稳定。实际聚合过程中，聚合物的分子量分布与催化剂、聚合工艺、均相关。本文的聚合物分子量分布曲线为特定Ziegler-Natta催化剂和聚合工艺条件下采集工业装置样品进行GPC测试获得的数据。本文目的是采用IA算法解决传统算法进行分子量分布去卷积分析时对参数初始值依赖性较高的问题，后续将在此文基础上，对聚合过程特性、工艺参数等对聚合物分子量分布、共聚组成分布以及序列组成分布的影响作进一步研究。

分别采用L-M和IA算法进行去卷积分析获得的各个活性位的质量分数和所生成聚合物的数均分子量如表 4所示。活性位个数分析结果不同，会导致各活性位的质量分数以及数均分子量变化，进而影响对聚合物微观结构的分析^[16-17]，如CCD、序列组成分布(sequence composition distribution，SCD)。此外，聚合动力学参数估计的前提是对催化剂活性位个数、所占比例等精准的计算。因此，为了提高IA算法的精度，结合L-M算法局部寻优的优点和IA算法无需初始值及全局搜索的优点，提出了IA-LM算法，以进一步提高IA算法对活性位质量分数和数均分子量的求解精度。图 4为IA-LM算法对MWD去卷积的分析结果，表 5为每个活性位的质量分数和数均分子量。图 4表明，IA-LM算法和IA算法的拟合效果基本上保持一致。由表 3可知，活性位个数从5增加到6时，GPC数据与模型的误差平方和保持不变，表明5个活性位足以描述聚合物的MWD，与IA算法的结论一致。对比两种算法的GPC数据与模型的误差平方和可知，相同活性位个数时，IA-LM算法的误差比IA算法更小，表明IA-LM算法求解精度更高。

4 结论

针对传统算法对MWD去卷积分析时需要准确初始值的问题，提出了基于IA算法的MWD去卷积分析方法。结论如下：

(1) 采用经典的L-M算法对MWD去卷积分析结果表明，该算法对初始值的依赖性较高，求解精度低。本文提出的IA算法可以在不用设定初始值的条件下，准确地分析出Ziegler-Natta催化剂活性位的个数以及每个活性位所生成聚合物的质量分数和数均分子量，且GPC数据与模型的最小误差平方和从0.171减小到0.008 7。

(2) 结合L-M算法的局部寻优和IA算法的全局搜索的优点，进一步提高了每个活性位生成聚合物的质量分数和数均分子量的去卷积分析精度，误差平方和由0.008 7减小到0.006 6。

研究结果为后续聚合物CCD、SCD去卷积分析以及聚合物精细链结构的预测奠定了基础。