1 前言

高密度聚乙烯(high-density polyethylene，HDPE)是5大通用合成树脂之一，具有优异的化学稳定性和耐低温性。我国是HPDE最大消费国，也是最大进口国^[1]，2019年消费量为15 478 kt，其中进口7 997 kt，产品自足率为49.39%。目前我国HDPE产品结构矛盾突出，中低端产能过剩，高端产品仍以进口为主。高密度聚乙烯实际上是混合物，由不同链长的高分子聚合物组成，组成成分不同会带来不同的材料特性，为对聚合物材料特性进行表征，质量分数的分布常常作为质量指标。质量分数的分布为所有不同链长分子质量占聚合物总质量的分数，它可以反映产品的微观性质^[2]。高密度聚乙烯模型的开发对提升产品的质量有极为重要的意义，而参数估计是实现基于质量分数分布的高密度聚乙烯建模的核心技术手段之一。

机理模型在反应过程的设计、优化和控制中起着重要的作用，参数估计是建立过程模型的关键组成部分。然而，由于参数可识别性和可估计性的问题，高质量参数值的估计并不容易^[3]。从另一个角度来看，参数估计问题往往是不定的，即不是所有的参数都可以唯一确定。常用的方法是通过将其他参数固定在合理的值上来估计参数的一个子集，而难点就在于确定可估计性高的参数子集。Yao等^[4]对参数子集选择方法进行了研究，其中正交化方法具有明显的几何含义并且易于实现，它将灵敏度矩阵S正交化排序，可以看作是每一步骤的D准则最大化的顺序方法，有效性被Chu等^[5]证明。正交化方法是参数子集选择中应用最广泛的技术^[6-9]，缺点是很难设置一个适当的投影向量范数阈值来分离可估参数和不可估参数。Machado等^[10]基于费雪信息矩阵(Fisher information matrix, FIM)准则进行参数子集选择，通过枚举法评估每个参数子集，此外还有行列式准则和一个修正的E准则之间的比率(a ratio between a determinant criterion and a modified E criterion，RDE)准则，但是随着参数数量上升，将会面临组合爆炸的问题。Brun等^[11]提出了相关性和共线性方法选择参数子集，通过灵敏度矩阵S的列中的线性关系来确定不可估计的参数。Jolliffe等^[12]提出主成分分析法(principal component analysis，PCA)，利用S^TS或者FIM的特征值来选择参数子集，缺点在于很难将特征值与特定的参数相匹配。此外Chen等^[13]提出用简约海森矩阵(reduced Hessian, H_r)代替FIM，并且证明了等效性，避免了灵敏度矩阵S的计算，并基于系统输出对参数进行排序，同时引入了关于待估参数的标准差与估计值比率(ratio of standard deviation to estimated value，RSDEV)准则用于可估和不可估参数的划分。

基于模型的实验设计(design of experiments，DOE)的目标是设计最佳的系统输入，并以经济的方式使收集到的输出中包含尽可能多的系统信息以便于后期的参数估计，采样时间的选择是DOE的一个重要组成部分，特别是在实验和样本分析成本高昂的情况下。有很多文献的最优采样使用了FIM的A准则、D准则，E准则^[14~16]，但是在参数之间存在高相关性的情况下，A准则是不可靠的，而D准则往往过于重视敏感参数^[17]。Kutalik等^[18]提出了一种次优采样时间选择方法，该方法是假设在选择采样时间之前，通过一些时间点进行参数估计，基于多重打靶法最大化FIM的行列式来选择采样时间。结合基于扩展灵敏度状态方程系统的优化采样设计和有限元的变化，Hoang等^[19]提出了一种先进的采样策略，优点是直接计算控制变量的参数灵敏度的精确高阶导数，并且使用联立求解。Chen等^[20]基于RSDEV准则，提出一种迭代增强的近似最优的时间序列采样策略。

HDPE模型十分复杂，使用自适应步长计算灵敏度矩阵S^[21]时，无法使用同一步长计算多个点的灵敏度信息，这使得灵敏度矩阵S的计算成本过于巨大，导致无法使用与S相关的方法去选择参数子集，同时也无法对模型中参数排序。因此本研究采用了Chen等^[13]提出的H_r矩阵代替FIM，并且结合RSDEV准则选出可估计参数子集，继而提出了一种基于拉丁超立方体采样的最优采样策略，结合RSDEV准则与H_r的关系建立了混合整数非线性规划(mixed integer nonlinear programming，MINLP)模型，并使用了一种新的迭代增强近似最优采样策略。数值结果表明，在RSDEV准则约束条件下，可以显著降低采样个数，同时得到可以接受的待估参数的相对误差。

2 高密度聚乙烯生产过程模型

HDPE由于生产成本低、机械和流变性能广泛，是应用最广泛的合成商品聚合物之一^[22]。多活性位Ziegler-Natta催化下的HDPE聚合反应是HDPE生产的重要工艺，聚合物在催化剂的活性位点生长，直到链转移发生，从而形成死聚合物链^[23]。

考虑模型基于在连续搅拌釜反应器(continuous stirred tank reactor，CSTR)中进行的HDPE浆工艺，由乙烯(单体)、氢气(链转移剂)、正己烷(稀释剂)和Ziegler-Natta催化剂系统构成，其中四氯化钛(TiCl₄)作为催化剂，三乙基铝(Al(C₂H₅)₃)作为助催化剂，连续注入CSTR，并且将浆液产物连续从反应器中去除。假设各个活性位吸附单体的能力相同，则链增长、链转移和链失活的动力学常数与链的长度无关^[24]。Ziegler-Natta催化剂的不同活性位点具有不同的传播速率常数，会产生很多不同链长的聚合物，从而导致一个非常宽广的质量分数的分布。建立一个由大量微分和代数方程(differential-algebraic equation，DAE)组成的联立方程模型，以描述该过程的动力学和质量分数的分布。

从过程控制的角度来看，可用的操纵变量为催化剂、助催化剂、单体、氢、六烷的供给量和反应器温度，而被控变量为单体转化率和质量分数的分布。

2.1 乙烯聚合反应机理

乙烯聚合反应采用的催化剂是Ziegler-Natta催化剂，是一种复合催化剂。乙烯聚合反应机理包括催化剂活化、链引发反应、链增长反应、链转移反应以及链失活反应^[25]。

1) 催化剂活化。通过Al(C₂H₅)₃对钛进行烷基化形成活性位点，反应如下：

$ {{\text{C}}_\mathrm{P}}(j) + {\text{A}}\xrightarrow{{{k_{\mathrm{aA}}}(j)}}{{\text{P}}_0}(j) $

(1)

式中：C_P(j)为具有潜在活性位点的催化剂TiCl₄；A为助催化剂Al(C₂H₅)₃；k_aA(j)为催化剂活化反应的动力学常数，L⋅(mol⋅s)⁻¹；j为活性位编号；P₀(j)为具有还未聚合任何单位的空活性位点的催化剂。

2) 链引发反应。单体分子与活性中心相互反应形成链长为1的聚乙烯链，反应如下：

$ {{\text{P}}_0}(j) + {\text{M}}\xrightarrow{{{k_\mathrm{i}}(j)}}{{\text{P}}_1}(j) $

(2)

式中：M为单体；k_i(j)为链引发反应的动力学常数，L⋅(mol⋅s)⁻¹；P₁(j)为链长为1且仍具有活性位点的乙烯聚合物链。

3) 链增长反应。继续插入乙烯单元使链长不断增长，反应如下：

$ {{\text{P}}_n}(j) + {\text{M}}\xrightarrow{{{k_\mathrm{p}}(j)}}{{\text{P}}_{n + 1}}(j) $

(3)

式中：k_p(j)为链增长反应的动力学常数，L⋅(mol⋅s)⁻¹；P_n(j)为链长为n的活聚合物链。

4) 链转移反应。活聚合物链分别与单体、氢气、助催化剂反应及自反应形成死链，反应如下：

$ {{\text{P}}_n}(j) + {\text{M}}\xrightarrow{{{k_{\mathrm{tM}}}(j)}}{{\text{P}}_0}(j) + {{\text{D}}_n}(j) $

(4)

$ {{\text{P}}_n}(j) + {{\text{H}}_2}\xrightarrow{{{k_{\mathrm{tH}}}(j)}}{{\text{P}}_0}(j) + {{\text{D}}_n}(j) $

(5)

$ {{\text{P}}_n}(j) + {\text{A}}\xrightarrow{{{k_{\mathrm{tA}}}(j)}}{{\text{P}}_0}(j) + {{\text{D}}_n}(j) $

(6)

$ {{\text{P}}_n}(j)\xrightarrow{{{k_\mathrm{t}}(j)}}{{\text{P}}_0}(j) + {{\text{D}}_n}(j) $

(7)

式中：k_tM(j）为向单体转移链反应的动力学常数，L⋅(mol⋅s)⁻¹；k_tH(j）为向氢气转移链反应的动力学常数，L^0.5⋅(mol^0.5⋅s⁻¹)；k_tA(j）为向助催化剂转移链反应的动力学常数，L⋅(mol⋅s)⁻¹；k_t(j）为自转移反应的动力学常数，s⁻¹；D_n(j)为链长为n的死聚合物链。

5) 链失活反应。活聚合物链与活性中心全部失活，反应如下：

$ {{\text{P}}_n}(j)\xrightarrow{{{k_\mathrm{d}}(j)}}{{\text{C}}_\mathrm{d}}(j) + {{\text{D}}_n}(j) $

(8)

$ {{\text{P}}_0}(j)\xrightarrow{{{k_\mathrm{d}}(j)}}{{\text{C}}_\mathrm{d}}(j) $

(9)

式中：k_d(j）为自失活反应的动力学常数，s⁻¹；C_d(j)为具有不能再聚合任何单体的死活性位点的催化剂。

在这个模型中，活性位点的数量被设为2，根据乙烯聚合反应机理，转移终止的伪动力学速率系数定义为

$ {k_{\mathrm{TD}}}(j) = {k_{\mathrm{tM}}}(j)c[{\text{M}}] + {k_{\mathrm{tH}}}(j){c^{0.5}}\left[ {{{\text{H}}_{\text{2}}}} \right] + {k_{\mathrm{tA}}}(j)c[{\text{A}}] + {k_\mathrm{t}}(j) + {k_\mathrm{d}}(j) $

(10)

式中：c [⋅]为括号中某成分的浓度，mol⋅L⁻¹。

2.2 物料衡算

为了简洁起见，模型中将在第j个活性位点的活和死聚合物的m阶矩分别定义为

$ {Y^m}(j) = \sum\limits_{n = 1}^\infty{{n^m}} c\left[ {{{\text{P}}_n}(j)} \right] $

(11)

$ {X^m}(j) = \sum\limits_{n = 1}^\infty{{n^m}} c\left[ {{{\text{D}}_n}(j)} \right] $

(12)

活聚合物链的零阶矩表示为

$ {Y^0}(j) = \sum\limits_{n = 1}^\infty{c\left[ {{{\text{P}}_n}(j)} \right]} $

(13)

对于长度为n=1的活链，反应净生产率为

$ {r_{{P_1}(j)}} = {k_\mathrm{i}}(j)c[{\text{M}}]c\left[ {{{\text{P}}_0}(j)} \right] + {k_{\mathrm{tM}}}(j)c[{\text{M}}]{Y^0}(j) - {k_\mathrm{p}}(j)c[M]c\left[ {{{\text{P}}_1}(j)} \right] - {k_{\mathrm{TD}}}(j)c\left[ {{{\text{P}}_1}(j)} \right] $

(14)

对于长度n≥2的活链来说，反应净生产率为

$ {r_{{P_n}(j)}} = {k_\mathrm{p}}(j)c[{\text{M}}]\left( {c\left[ {{{\text{P}}_{n - 1}}(j)} \right] - c\left[ {{{\text{P}}_n}(j)} \right]} \right) - {k_{\mathrm{TD}}}(j)c\left[ {{{\text{P}}_1}(j)} \right] $

(15)

同样地，长度为n=1和n≥2的死链的反应净生产率分别为

$ {r_{{D_1}(j)}} = {k_{\mathrm{TD}}}(j)c\left[ {{{\text{P}}_{\text{1}}}(j)} \right] - {k_{\mathrm{tM}}}(j)c[{\text{M}}]c\left[ {{{\text{P}}_1}(j)} \right] $

(16)

$ {r_{{D_n}(j)}} = {k_{\mathrm{TD}}}(j)c\left[ {{{\text{P}}_n}(j)} \right] $

(17)

从式(14)~(17)可得，长度n=1, 2, 3, …, n_max的活链和死链的粒数衡算方程分别为

$ \frac{{{\text{d}}\left( {c\left[ {{{\text{P}}_n}(j)} \right]V} \right)}}{{{\text{d}}t}} = {r_{{P_n}(j)}} \cdot V - c\left[ {{{\text{P}}_n}(j)} \right] \cdot {Q_{{\text{out }}}} $

(18)

$ \frac{{{\text{d}}\left( {c\left[ {{{\text{D}}_n}(j)} \right]V} \right)}}{{{\text{d}}t}} = {r_{{D_n}(j)}} \cdot V - c\left[ {{{\text{D}}_n}(j)} \right] \cdot {Q_{{\text{out }}}} $

(19)

式中：Q_out为体积流出率，%；V为反应堆的液体体积，L。类似的单体、氢气、己烷、催化剂潜在活性位点、催化剂无活性位点、活性点、共催化剂以及活聚合物链的零阶矩和一阶矩的物料衡算方程也可以写为

$ \frac{{{\text{d}}(c[{\text{M}}]V)}}{{{\text{d}}t}} = {F_{{\text{in}}}}\left( {\text{M}} \right) - \sum\limits_{j = 1}^{{N_\mathrm{s}}} {\left( {{k_\mathrm{i}}(j)c\left[ {{{\text{P}}_0}(j)} \right] + \left( {{k_\mathrm{p}}(j) + {k_{\mathrm{tM}}}(j)} \right){Y^0}(j)} \right)}\cdot c[{\text{M}}] \cdot V - c[{\text{M}}] \cdot {Q_{{\text{out }}}} $

(20)

$ \frac{{{\text{d}}\left( {c\left[ {{{\text{H}}_2}} \right]V} \right)}}{{{\text{d}}t}} = {F_{{\text{in}}}}\left( {{{\text{H}}_2}} \right) - \sum\limits_{j = 1}^{{N_\mathrm{s}}} {\left( {{k_{\mathrm{tH}}}(j){c^{^{0.5}}}\left[ {{{\text{H}}_2}} \right]{Y^0}(j)} \right)}\cdot V - c\left[ {{{\text{H}}_2}} \right] \cdot {Q_{{\text{out }}}} $

(21)

$ \frac{{{\text{d}}\left( {c\left[ {{{\text{C}}_{\text{6}}}{{\text{H}}_{{\text{14}}}}} \right]V} \right)}}{{{\text{d}}t}} = {F_{{\text{in}}}}\left( {{{\text{C}}_6}{{\text{H}}_{14}}} \right) - c\left[ {{{\text{C}}_{\text{6}}}{{\text{H}}_{{\text{14}}}}} \right] \cdot {Q_{{\text{out }}}} $

(22)

$ \frac{{{\text{d}}\left( {c\left[ {{{\text{C}}_\mathrm{p}}} \right]V} \right)}}{{{\text{d}}t}} = {F_{{\text{in}}}}\left( {{\mathrm{C}_\mathrm{P}}} \right) - \sum\limits_{j = 1}^{{N_\mathrm{s}}} {\left( {{k_{\mathrm{aA}}}(j)c[{\text{A}}]c\left[ {{{\text{C}}_\mathrm{p}}} \right]} \right)}\cdot V - c\left[ {{{\text{C}}_\mathrm{p}}} \right] \cdot {Q_{{\text{out}}}} $

(23)

$ \frac{{{\text{d}}\left( {c\left[ {{{\text{C}}_\mathrm{d}}} \right]V} \right)}}{{{\text{d}}t}} = \sum\limits_{j = 1}^{{N_\mathrm{s}}} {\left( {{k_\mathrm{d}}(j)\left( {c\left[ {{{\text{P}}_0}(j)} \right] + {Y^0}(j)} \right)} \right)}\cdot V - c\left[ {{{\text{C}}_\mathrm{d}}} \right] \cdot {Q_{{\text{out }}}} $

(24)

$ \begin{gathered} \frac{{{\text{d}}\left( {c\left[ {{{\text{P}}_0}(j)} \right]V} \right)}}{{{\text{d}}t}} = {k_{\mathrm{aA}}}(j)c[{\text{A}}]c\left[ {{{\text{C}}_\mathrm{P}}} \right] + \left( {{k_{\mathrm{tH}}}(j){c^{^{0.5}}}\left[ {{{\text{H}}_2}} \right] + {k_{\mathrm{tA}}}(j)c[{\text{A}}] + {k_\mathrm{t}}(j)} \right){Y^0}(j) -\\ \left. {{k_\mathrm{i}}(j)c[{\text{M}}]c\left[ {{{\text{P}}_0}(j)} \right] - {k_\mathrm{d}}(j)c\left[ {{{\text{P}}_0}(j)} \right]} \right) - c\left[ {{{\text{P}}_0}(j)} \right] \cdot {Q_{\mathrm{out}}} \\ \end{gathered} $

(25)

$ \frac{{{\text{d}}(c[{\text{A}}]V)}}{{{\text{d}}t}} = {F_{{\text{in}}}}\left( {\text{A}} \right) - \sum\limits_{j = 1}^{{N_\mathrm{s}}} {\left( {{k_{\mathrm{aA}}}(j)c\left[ {{{\text{C}}_\mathrm{P}}} \right] + {k_{\mathrm{tA}}}(j){Y^0}(j)} \right)}\cdot c[{\text{A}}] \cdot V - c[{\text{A}}] \cdot {Q_{{\text{out }}}} $

(26)

$ \frac{{{\text{d}}\left( {{Y^0}(j)V} \right)}}{{{\text{d}}t}} = \left( {{k_\mathrm{i}}(j)c[{\text{M}}]c\left[ {{{\text{P}}_0}(j)} \right] + {k_{\mathrm{tM}}}(j)c[{\text{M}}]{Y^0}(j) - {k_{\mathrm{TD}}}(j){Y^0}(j)} \right) \cdot V - {Y^0}(j) \cdot {Q_{{\text{out }}}} $

(27)

$ \frac{{{\text{d}}\left( {{Y^1}(j)V} \right)}}{{{\text{d}}t}} = \left( {{k_\mathrm{i}}(j)c[{\text{M}}]c\left[ {{{\text{P}}_0}(j)} \right] + {k_p}(j)c[{\text{M}}]{Y^0}(j) + {k_{\mathrm{tM}}}(j)c[{\text{M}}]{Y^0}(j) - {k_{\mathrm{TD}}}(j){Y^1}(j)} \right) \cdot V - {Y^1}(j) \cdot {Q_{out}} $

(28)

$ \frac{{{\text{d}}\left( {{X^1}(j)V} \right)}}{{{\text{d}}t}} = \left( {{k_{\mathrm{TD}}}(j){Y^1}(j) - {k_{\mathrm{tM}}}(j)c[{\text{M}}]c\left[ {{{\text{P}}_1}(j)} \right]} \right) \cdot V - {X^1}(j) \cdot {Q_{{\text{out }}}} $

(29)

式中：N_s为活性位点个数，N_s=2；F_in(M)、F_in(H₂)、F_in(C₆H₁₄)、F_in(C_P)、F_in(A)分别为M、H₂、C₆H₁₄、C_P、A的摩尔进料流率，mol⋅h⁻¹。

Ma等^[26]提出一种利用二维正交配置的质量分数分布重构方法，以捕捉质量分数分布的时间动态特征和链长度的分布特征。采用联立配置法离散化模型，将模型最终离散成具有86 116个变量，86 100条等式约束的大型非线性规划(nonlinear programming，NLP)问题。16个待估参数分别为k_p(1)、k_p(2)、k_tH(1)、k_tH(2)、k_tM(1)、k_tM(2)、k_tA(1)、k_tA(2)、k_aA(1)、k_aA(2)、k_i(1)、k_i(2)、k_t(1)、k_t(2)、k_d(1)、k_d(2)。模型采用以下目标函数进行参数估计：

$ {\min _{}}\;\;J = \frac{1}{{2{\sigma ^2}}}\sum\limits_{i = 1}^{{t_\mathrm{f}}} {{{\left( {{w_{\mathrm{cal}}} - {w_{\mathrm{rel}}}} \right)}^2}} $

(30)

式中：$ {\sigma ^2} $为测量噪声的方差，t_f为时间采样点，w_cal为模型质量分数的分布仿真值，w_rel为模型质量分数的分布实际值。

3 最优采样

由于HDPE采样后的分析成本过高，故尽可能减少采样数量，同时尽量获得多的信息，本研究将结合RSDEV准则和FIM的A准则来选择参数子集，然后基于拉丁超立方体采样，使用Chen等^[20]提出的一种基于H_r的近似最优的时间序列采样策略。

3.1 参数子集选择

在参数子集选择方法上，Chen等^[13]证明了H_r与FIM的等效性，故采用H_r代替FIM，并且将基于FIM准则选择参数子集。常用的FIM准则包括了A准则、D准则、修正后的E准则，此外还有RDE准则，如果参数数量太多，则会面临组合爆炸问题，为避免此问题，采用Chen等^[13]提出的RSDEV准则，即

$ \varepsilon= \frac{{{\delta _q}}}{{{\theta _q}}} $

(31)

式中：δ_q为第q个待估参数标准差；$ {\theta _q} $为第q个待估参数的估计值；ε为比率；δ_q=a_{q, q}^1/2，a_{q, q}为相应参数的协方差矩阵A的对角元素的值，即A矩阵的第q个对角元素，A=FIM⁻¹=H_r⁻¹。

将FIM的A准则和RSDEV准则结合在一起形成新的准则：

min (tr (FIM⁻¹))，且满足max (ε) < ε_exp。

式中：ε_exp为期望的比率阈值。

3.2 最优采样策略

将HDPE模型表达成微分代数方程组，即

$ \frac{{{\text{d}}{{\boldsymbol{x}}}{(}t{)}}}{{{\text{d}}t}} = f{(}{\boldsymbol{x}}{(}t{)},{\boldsymbol{y}}{(}t{)},{ \mathit{\pmb{ θ}} }{)},\;\;\;{\boldsymbol{x}}\left( {{t_0}} \right) = {{\boldsymbol{x}}_0} $

(32a)

$ g({\boldsymbol{x}}(t),{\boldsymbol{y}}(t),{\boldsymbol{ \mathit{\pmb{ θ}} }}) = 0 $

(32b)

$ {\boldsymbol{z}}{(}t{)} = \widehat h{(}{\boldsymbol{x}}{(}t{)},{\boldsymbol{y}}{(}t{)},{{ \mathit{\pmb{ θ}} }}{)} $

(32c)

式中：t为连续时间，s；x(t)为状态变量向量，x(t)∈R^nx，R^nx为nx维的实向量空间；y(t)为代数变量向量，y(t)∈R^ny，其中R^ny为ny维的实向量空间；z(t) 为模型预测向量，z(t)∈R^nz，其中R^nz为nz维的实向量空间；θ为参数向量，θ∈R^nq，其中R^nq为nq维的实向量空间；f表示从R^nx×R^ny×R^nq到R^nx的映射，g表示从R^nx×R^ny×R^nq到R^ny的映射，$ \widehat h $表示R^nx×R^ny×R^nq到R^nz的映射，t₀为初始时刻，x₀为初始时刻的状态变量向量。采用一种联立配置的方法，微分和代数状态变量用有限元上[t_i-1, t_i](i=1, …, N)的多项式族来逼近，且满足t₀ < t₁ < … < t_N=t_f。在每个有限元上使用高斯点或Radau点的拉格朗日插值多项式如下：

$ \begin{array}{*{20}{l}} {{{\boldsymbol{x}}^K}(t) = \sum\limits_{k = 0}^K {{l_k}} (\tau ){{\boldsymbol{x}}_{ik}},\;\;\;{{\boldsymbol{y}}^K}(t) = \sum\limits_{k = 0}^K {{{\bar l}_k}} (\tau ){{\boldsymbol{y}}_{ik}},\;\;\;t = {t_{i - 1}} + \left( {{t_i} - {t_{i - 1}}} \right)\tau ,\;\;\;\tau\in [0,1]} \\ {{l_k}(\tau ) = \prod\limits_{j = 0, \ne k}^K {\frac{{\left( {\tau- {\tau _j}} \right)}}{{\left( {{\tau _k} - {\tau _j}} \right)}}} ,\;\;\;{{\bar l}_k}(\tau ) = \prod\limits_{j = 1, \ne k}^K {\frac{{\left( {\tau- {\tau _j}} \right)}}{{\left( {{\tau _k} - {\tau _j}} \right)}}} } \end{array} $

(33)

式中：K为插值的阶次，也是配置点的数量，x_ik为在第i个有限元的第k个配置点的状态变量向量，y_ik为在第i个有限元的第k个配置点的代数变量向量，τ为有限元中t的相对值，$ {\tau _k} $、$ {\tau _j} $分别为有限元中第k、j个配置点t的相对值。在第l个采样时间点t_l处的预测可以用以下方程来近似：

$ {{\boldsymbol{z}}_{l,{\text{cal}}}} = \hat h\left( {{{\boldsymbol{x}}^K}\left( {{t_l}} \right),{{\boldsymbol{y}}^K}\left( {{t_l}} \right),{\boldsymbol{ \mathit{\pmb{ θ}} }}} \right) $

(34)

式中：$ {{\boldsymbol{z}}_{l,{\text{cal}}}} $为采样时间点$ {t_l} $处的模拟预测向量。

假设一个{1, …, ntp}的时间序列，它可以作为潜在的采样点，相应在t_l的测量是z_l，l=1, …, ntp。基于似然原理，假设z的每个分量的测量噪声呈正态分布，均值为0、方差为σ_p²，p=1, …, nz，通过求解问题(35)可以估计θ：

$ \begin{array}{*{20}{l}} {\min J = \frac{1}{2}\sum\limits_{l \in S} {{{\left[ {{{\boldsymbol{z}}_{l,{\text{cal}}}} - {{\boldsymbol{z}}_{l,{\text{rel}}}}} \right]}^{\text{T}}}} {\mathit{\pmb{\Sigma}} ^{ - 1}}\left[ {{{\boldsymbol{z}}_{l,{\text{cal}}}} - {{\boldsymbol{z}}_{l,{\text{rel}}}}} \right]} \\ {{\text{ s}}{\text{.t}}{\text{. }}\left\{ {\begin{array}{*{20}{l}} {{{(a)}_{}}\sum\limits_{j = 0}^K {{{l'}_j}} \left( {{\tau _k}} \right){{\boldsymbol{x}}_{ij}} - \left( {{t_i} - {t_{i - 1}}} \right)f\left( {{{\boldsymbol{x}}_{ik}},{{\boldsymbol{y}}_{ik}},{{ \mathit{\pmb{ θ}} }}} \right) = 0} \\ {{{(b)}_{}}{{\boldsymbol{x}}_{1,0}} = {{\boldsymbol{x}}_0},\;\;{{\boldsymbol{x}}_{i + 1,0}} = \sum\limits_{j = 0}^K {{l_j}} (1){{\boldsymbol{x}}_{ij}},\;\;i = 1, \cdots ,N - 1} \\ {{{{\text{(c)}}}_{}}g{\text{ }}\left( {{{\boldsymbol{x}}_{ik}},{{\boldsymbol{y}}_{ik}},{{ \mathit{\pmb{θ}} }}} \right) = 0} \\ {{{({\text{d}})}_{}}{{\boldsymbol{z}}_l} = \hat h\left( {{{\boldsymbol{x}}^K}\left( {{t_l}} \right),{{\boldsymbol{y}}^K}\left( {{t_l}} \right),{{ \mathit{\pmb{ θ}} }}} \right)} \end{array}} \right.} \end{array} $

(35)

式中：z_{l, rel}为采样时间点t_l处的实际预测向量，Σ为对角线方差矩阵，其中第p个元素为σ_p²；$ {l'_j} $($ {\tau _k} $)为导数。

设S={t_1, t_2, …, t_n}是采样点集合，是时间序列{1, 2, …, ntp}的一个子集。为了便于最优抽样分析，将S={1, 2, …, ntp}，并且将非线性规划问题(35)中的离散变量和方程被重写为

$ \begin{array}{*{20}{l}} {{\text{min}}\;\;J = \frac{1}{2}{{[{\boldsymbol{w}} - {{\boldsymbol{w}}^*}]}^\mathrm{T}}{\boldsymbol{V}}_{\boldsymbol{w}}^{ - 1}[{\boldsymbol{w}} - {{\boldsymbol{w}}^*}]} \\ {{\text{ s}}{\text{.t}}{\text{. }}\left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{c}}({\boldsymbol{ \mathit{\pmb{ γ}} }},{\boldsymbol{ \mathit{\pmb{ θ}} }}) = 0} \\ {{\boldsymbol{w}} = {\boldsymbol{h}}(\mathit{\pmb{\gamma}} ,{\boldsymbol{ \mathit{\pmb{ θ}} }})} \end{array}} \right.} \end{array} $

(36)

式中：w=[z_{1, cal}^T, z_{2, cal}^T, …, z_{ntp, cal}^T]^T，w*=[z_{1, rel}^T, z_{2, rel}^T, …, z_{ntp, rel}^T]^T，γ=[x₁₁^T, y₁₁^T, …, x_1K^T, y_1K^T, …x_NK^T, y_NK^T]^T，V_w是一个对角矩阵，对角元素由Σ的ntp个对角元素组成，c表示方程(35(a)-35(c))，h表示方程(35(d))。在Chen等 ^[13]的论文已经证明，估计参数的协方差矩阵A等于问题(36)的目标函数的H_r的逆，其中H_r计算为

$ {{\boldsymbol{H}}_{\bf r}} = {{\boldsymbol{Z}}^\mathrm{T}}\left[ {\begin{array}{*{20}{c}} 0&0&0 \\ 0&0&0 \\ 0&0&{{\boldsymbol{V}}_{\boldsymbol{w}}^{ - 1}} \end{array}} \right]{\boldsymbol{Z}},\;\;{\boldsymbol{Z}} = \left[ \begin{gathered} {\boldsymbol{I}} \\ - \left[ {\begin{array}{*{20}{c}} {\frac{{\partial \boldsymbol{c}}}{{\partial \mathit{\pmb{\gamma}} }}{{\boldsymbol{c}}^\mathrm{T}}}&0 \\ {\frac{{\partial \boldsymbol{h}}}{{\partial \mathit{\pmb{\gamma }}}}{{\boldsymbol{h}}^\mathrm{T}}}&{ - {\boldsymbol{I}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\frac{{\partial \boldsymbol{c}}}{{\partial \mathit{\pmb{\theta }}}}{{\boldsymbol{c}}^\mathrm{T}}} \\ {\frac{{\partial \boldsymbol{h}}}{{\partial \mathit{\pmb{\theta }}}}{{\boldsymbol{h}}^\mathrm{T}}} \end{array}} \right] \\ \end{gathered}\right] $

(37)

式中：I为单位矩阵。

确定δ_q/θ_q的值，即参数q的标准差除以其估计值，建立一个MINLP问题寻找输出测量时间序列的最小子集，从而可以找到一个满足上述比率的A矩阵。

首先定义了Z_r，它由Z的最后一个nz×ntp行组成：

$ {{\boldsymbol{Z}}_{\bf r}} = \left[ {\begin{array}{*{20}{c}} {z_{1,1}^\text{r}}& \cdots &{z_{1,\mathrm{nz}}^\text{r}}& \cdots &{z_{1,\mathrm{nz} \times (\mathrm{ntp} - 1) + 1}^\text{r}}& \cdots &{z_{1,\mathrm{nz} \times \mathrm{ntp}}^\text{r}} \\ {z_{2,1}^\text{r}}& \cdots &{z_{2,\mathrm{nz}}^\text{r}}& \cdots &{z_{2,\mathrm{nz} \times (\mathrm{ntp} - 1) + 1}^\text{r}}& \cdots &{z_{2,\mathrm{nz} \times \mathrm{ntp}}^\text{r}} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots\\ {z_{\mathrm{nq},1}^\text{r}}& \cdots &{z_{\mathrm{nq},\mathrm{nz}}^\text{r}}& \cdots &{z_{\mathrm{nq},\mathrm{nz} \times (\mathrm{ntp} - 1) + 1}^\text{r}}& \cdots &{z_{\mathrm{nq},\mathrm{nz} \times \mathrm{ntp}}^\text{r}} \end{array}} \right] $

(38)

用二元变量b=[b₁, b₂, …, b_ntp]^T来确定哪个采样时刻t_l(l=1, 2, …, ntp)从时间序列中被选择。因此，关于时间序列S的H_r可以定义为

$ {{\boldsymbol{H}}_{\bf r}}({\boldsymbol{b}}) = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{l = 1}^{\mathrm{ntp}} {{b_l}} \sum\limits_{p = 1}^{\mathrm{nz}} {\frac{1}{{\sigma _p^2}}} z_{1,p + \mathrm{nz} \times (l - 1)}^\text{r}z_{1,p + \mathrm{nz} \times (l - 1)}^\text{r}}& \cdots &{\sum\limits_{l = 1}^{\mathrm{ntp}} {{b_l}} \sum\limits_{p = 1}^{\mathrm{nz}} {\frac{1}{{\sigma _p^2}}} z_{1,p + \mathrm{nz} \times (l - 1)}^\text{r}z_{\mathrm{nq},p + \mathrm{nz} \times (l - 1)}^\text{r}} \\ \vdots & \ddots & \vdots\\ {\sum\limits_{l = 1}^{\mathrm{ntp}} {{b_l}} \sum\limits_{p = 1}^{\mathrm{nz}} {\frac{1}{{\sigma _p^2}}} z_{\mathrm{nq},p + \mathrm{nz} \times (l - 1)}^\text{r}z_{1,p + \mathrm{nz} \times (l - 1)}^\text{r}}& \cdots &{\sum\limits_{l = 1}^{\mathrm{ntp}} {{b_l}} \sum\limits_{p = 1}^{\mathrm{nz}} {\frac{1}{{\sigma _p^2}}} z_{\mathrm{nq},p + \mathrm{nz} \times (l - 1)}^\text{r}z_{\mathrm{nq},p + \mathrm{nz} \times (l - 1)}^\text{r}} \end{array}} \right] $

(39)

在当前参数θ下的协方差矩阵A=(H_r(b))⁻¹，通过求解以下MINLP问题来确定b，时间序列子集和合适的A矩阵元素a_{m, s}：

$ \begin{array}{*{20}{l}} {\min \sum\limits_{l = 1}^{\mathrm{ntp}} {{b_l}} } \\ {{\text{ s}}{\text{.t}}{\text{. }}\left\{ {\begin{array}{*{20}{l}} {\sum\limits_{m = 1}^{\mathrm{\mathrm{nq}}} {{a_{m,s}}} \left( {\sum\limits_{l = 1}^{\mathrm{ntp}} {{b_l}} \sum\limits_{p = 1}^{\mathrm{nz}} {\frac{1}{{\sigma _p^2}}} z_{q,p + \mathrm{nz} \times (l - 1)}^rz_{m,p + \mathrm{nz} \times (l - 1)}^r} \right) = \left\{ {\begin{array}{*{20}{l}} {0,\quad s \ne q} \\ {1,\quad s = q} \end{array}} \right.} \\ {0 \leqslant {a_{q,q}} \leqslant \varepsilon _{\exp }^2\theta _q^2,\quad q = 1, \cdots ,\mathrm{nq};\;\quad- \varepsilon _{\exp }^2{\theta _m}{\theta _s} \leqslant {a_{m,s}} \leqslant \varepsilon _{\exp }^2{\theta _m}{\theta _s},\quad m,s, = 1, \cdots ,\mathrm{nq},\quad m \ne s} \end{array}} \right.} \end{array} $

(40)

然而，使用问题式(40)的前提是，真实的估计参数值是先验已知的，但事实并非如此。相反，该策略通过在参数空间上进行随机采样来选择θ的候选值，并用每个候选值来解决问题(40)。从每个解中，记录从b中选择了哪些时间序列点，并观察从整个采样中选择时间序列点的频率，然后选择大频率的时间序列点进行参数估计。

并且Chen等^[20]证明了输出采样数量的增加会使得RSDEV准则ε的值减小或者保持不变。因此，MINLP(40)的解决方案提供了一个最小的时间序列点集，并且同时保证了一个适当的比率ε。

3.3 拉丁超立方体采样

采样策略多次在参数空间中随机采样选择θ的候选值，并用每个候选值来解决问题式(40)，然而随机采样使结果充满不确定性。建立精确的代理模型要求采样点集可以描述模型在输入范围上的整体趋势变化，因而需要采样点足够均匀，数量也足够代表模型的变化，为了达到这一点，本研究采用了拉丁超立方体采样(Latin hypercube sampling，LHS)取代随机采样。LHS是分层采样抽样技术中的一种，由Mckay等^[27]于1979年提出。分层抽样法将总体按一定特征或者规则划分为不同的层，然后从不同的层内进行独立、随机的抽样，相比于简单随机抽样，分层采样获得的样本代表性较好，抽样误差较小。

在LHS方法中在每个维度上都进行了均匀的划分，但是采样点的数目不受限制。假设采样空间的维度为n，想要获得m个采样点，那么要将n个维度中的每个维度上的取值空间划分为m个等份，然后在每个维度的每个等分内都随机产生一个数值，由此就获得了n组数量为m的数值集合，之后将这n组数值集合内的点随机组合，就获得了所需要的m个采样点^[28]。

在此策略中，维度由参数子集选择技术决定，设可估参数子集个数为nq，采样次数为nrs，将参数子集的参数空间分成nrs个等份，在每个等份中随机取一个值，将nq组数值集合内的点随机组合，获得输入采样。

3.4 算法步骤

步骤1：通过模型式(32)确定可估计参数子集，可估参数子集个数为nq。

步骤2：设定式(40)中的ε_exp，乘子η > 1，采样次数为nrs。

步骤3：设定可估计参数的参数空间，用拉丁超立方体采样在参数空间里采样nrs次生成候选参数，并解决nrs次候选参数的问题式(40)，然后根据所有采样的参数值计算时间序列点的频率。

步骤4：设最终采样集合S由nq×η个时间序列点组成，选择频率最高的nq×η个时间点组成S。

步骤5：用集合S中的采样点时刻的测量值来求解问题(36)，即参数估计；计算参数协方差矩阵A，并表示最大比率δ_q/θ_q，q=1, …, nq。

步骤6：如果δ_q/θ_q < ε_exp，转至步骤9。

步骤7：如果δ_q/θ_q≥ε_exp，设S_l={1, 2, …, ntp}−S，对S_l中的每个元素{v}，令S_t=S+{v}，并计算其H_r和协方差矩阵A，并保留其δ_q*/θ_q*的值。

步骤8：选择δ_q*/θ_q*的最大值对应的元素v_max，令S=S+{v_max}，返回步骤5。

步骤9：结束。

算法流程图如图(1)所示：

4 仿真实验

以质量分数的分布作为测量输出，通过加入方差1×10⁻⁴来模拟测量值，有限元个数为30个，每个有限元配置点K=3，共有90个潜在采样点。根据第3节的参数子集选择方法，设定ε_exp=0.05，被估计参数子集θ=[k_p(2), k_tH(1), k_aA(2)]，nq=3，其参数空间的范围分别为3 554.041≤k_p(2)≤3 871.933，L⋅(mol⋅s)⁻¹；182.361 ≤k_tH(1)≤216.283，L^0.5⋅mol^−0.5⋅s⁻¹；204.329≤k_aA(2) ≤242.987，L⋅(mol⋅s)⁻¹；设定式(40)中的ε_exp=0.05，乘子η=2，采样次数nrs为100。在最优采样策略的步骤3中，使用随机采样，得出所有时间点的频率，如图 2所示。

根据频率最高选择采样时间点S={9, 15, 43, 69, 77, 90}，采样时间点数量为6。使用选择的采样时间点估计参数θ=[k_p(2), k_tH(1), k_aA(2)]，其值分别为[3 423.207, 221.58, 244.6]，其δ_q/θ_q分别为[0.038 8, 0.045 5, 0.048 4]，误差分别为6.23%、6.27%、7.23%。

用拉丁超立方体采样代替步骤3中的随机采样，同样得到所有采样的参数值计算时间序列点的频率，如图 3所示。

根据频率最高选择采样时间点S={9, 12, 15, 43, 77, 90}，采样时间点数量为6。使用选择的采样时间点估计参数θ=[k_p(2), k_tH(1), k_aA(2)]，其值分别为[3 542.722, 213.84, 237.2]，其δ_q/θ_q分别为[0.033 3, 0.036 9, 0.040 5]，误差分别为2.96%、2.56%、3.97%。

为对比在步骤3中使用拉丁超立方体采样与使用随机采样对最终参数估计的影响，分别使用这2种采样方法仿真了10次，并且取其误差平均值作为对比。使用随机采样的参数估计的误差平均值为5.82%、6.12%、6.99%，其方差分别为0.000 325，0.000 436，0.000 387，使用拉丁超立方体采样的参数估计的误差平均值为3.63%、3.49%、4.11%，其方差分别为0.000 377、0.000 631、0.000 563。

5 结语

针对HDPE模型的样本分析成本过高且参数估计困难问题，提出使用简约Hessian矩阵代替FIM并结合RSDEV准则作为新的参数子集选择方法，其优势在于避免了HDPE模型灵敏度矩阵S计算成本过大和组合爆炸问题。此外，在保证参数估计精度的情况下，为了减少采样点个数，降低样本分析成本，提出一种基于拉丁超立方体采样的近似最优采样策略。通过拉丁超立方体采样获得参数空间候选值，并建立和求解最优采样子问题，构建了各潜在采样点的选择频率，使用一种启发式策略和迭代增强策略选择最终采样点，并实现了参数估计。为了验证采样策略结果的可靠性，使用了基于随机采样的最优采样策略进行对比。使用新的参数子集选择方法选出的待估参数为k_p(2)、k_tH(1)、k_aA(2)，同时仿真结果表明，最终在90个潜在采样点中仅采样其中6个，大大减少了采样点个数，降低了样本分析成本，使用随机采样的参数估计相对误差平均值分别为5.82%、6.12%、6.99%，使用拉丁超立方体采样的参数估计的相对误差平均值分别为3.63%、3.49%、4.11%。从仿真结果可以看出参数估计精度在可接受范围，而且基于拉丁超立方体采样比基于随机采样的参数估计结果更加精确、稳定。此外，该最优采样策略还适用于不能在线测量和对样本离线分析成本很高的情况。