1 前言

相比于连续过程，间歇过程没有稳定的静态工作点，过程操作的基本方式是确保批次过程跟踪设定操作曲线。间歇过程的优化目标通常为缩短运行周期或提高终端产品产量和质量^[1]。在完成有效的动态控制基础上，间歇过程操作曲线优化成为提高经济效益的关键。

间歇过程的操作优化概括分为机理驱动型和数据驱动型^[2]。前者需要过程精确的机理模型^{[3, 4]}，而间歇过程因自身复杂的动态非线性特性而难以建模，使得机理驱动型优化的发展受到了限制。数据驱动型优化大体分为基于经验模型的优化，实验设计以及无模型优化。常用的经验模型有：人工神经网络^{[5, 6]}，偏最小二乘^{[7, 8]}，支持向量机^{[9, 10]}等。由于经验模型精度受限于样本数，一般而言模型外推性差；实验设计^[11~13]通过有目的的试验获得过程响应与输入间的数值关系，但实验设计需要利用实际运行装置，成本高且效率低。结合基于模型优化的高效性和实验设计的可行性，孔祥松等人提出了基于迭代点在线试验的无模型优化方法^[14~16]。该方法避免了模型失配，提高了优化效率，但它仅适用于快速且低成本的间歇过程，对于具有较长工作周期的间歇过程，文献[17]提出了基于批次数据的迭代优化方法。而有些间歇过程具有多阶段特征^{[18, 19]}。忽略间歇过程局部的阶段特性，仅关注过程变量对优化指标的影响，显然降低了优化性能。

针对上述问题，本文提出一种基于多阶段时段变量相关性分析的操作曲线优化方法。首先划分时间片矩阵并对其主元分析，对计算得到的加权时间片载荷矩阵聚类分析得出过程的阶段，然后在每一阶段利用时段变量对指标变量的作用强弱和方向实现操作曲线的修正，最后考虑阶段与指标的相关系数，计算出权重融合各阶段操作曲线。该方法被成功用于某化工产品的间歇结晶过程温度曲线优化中。

2 问题描述

间歇过程操作优化的目的是在一个批次过程结束时，产品的产量或质量得到提高、过程操作的周期能够缩短。在间歇过程直接控制层已经实现对给定曲线的动态跟踪基础上，优化层的目标是如何确定最优的操作曲线，其数学描述如下：

$ \begin{array}{l} \mathop {_{}\max }\limits_{s(t)} \ \ \ {\rm{ }}J\left[ {x({t_{\rm f}})} \right]\\ _{}{\rm{s}}{\rm{.t}}{\rm{.}}\\ \left\{ \begin{array}{l} {\rm{0}} \le t \le {t_{\rm f}}\\ \dot x(t) = f\left[ {x(t), s(t)} \right]\\ y(t) = \phi \left[ {x(t), s(t)} \right]\\ s{(t)_{\min }} \le s(t) \le s{(t)_{\max }} \end{array} \right.{\rm{ }} \end{array} $

(1)

其中t_f为间歇过程周期，J是性能指标，s(t)为操作曲线设定值，x(t)是被控过程的状态变量，y(t)是输出变量，s(t)_min和s(t)_max分别是操作曲线的下界和上界。式(1)意义为在满足约束的条件下寻找一条合适的操作曲线，使得间歇过程结束时性能指标最优或次优。

上式的优化问题本质上是一个非线性优化问题，为了避免直接处理非线性，首先采用分段离散化方法^[20]将原操作曲线以及状态曲线用m - 1个分段点等时间间隔分成了m段(图 1)，即：

$ s(t) \buildrel \Delta \over = \left\{ {{s_l}\left| {l = 1, 2, } \right. \cdots , m} \right\} $

(2)

$ x(t) \buildrel \Delta \over = \left\{ {{x_l}\left| {l = 1, 2, } \right. \cdots , m} \right\} $

(3)

其中s_l为离散化后的操作曲线设定值，x_l为离散化后的状态变量，记为时段变量。由此非线性优化问题(1)转化为离散化如下形式：

$ \begin{array}{l} \mathop {\max }\limits_{\{ {s_l}\left| {l = 1, 2, \cdots , m} \right.\} } \ \ \ {\rm{ }}J\left[ {{\rm{(}}x({t_{\rm f}})} \right]\\ {\rm{s}}{\rm{.t}}{\rm{.}}\\ \left\{ \begin{array}{l} {x_{l + 1}} = f\left[ {{x_l}, {s_l}} \right]\\ {y_l} = \phi \left[ {{x_l}, {s_l}} \right]\\ {s_{l, \min }} \le {s_l} \le {s_{l, \max }} \end{array} \right.{\rm{ }} \end{array} $

(4)

其中s_{l, min}和s_{l, max}分别是操作曲线离散化后的下界和上界，y_l是离散化后的输出变量。优化的目标为找出合适的操作曲线设定值集合{s_l|l=1, 2, …, m}使性能指标最优。以下讨论中，作者将性能指标的终值J[x(t_f)]称为批次的指标变量。

3 多阶段优化策略及算法 3.1 间歇过程阶段的划分

多阶段是间歇过程具有的显著特点，时段变量对指标变量的影响具有时间区域特性。操作曲线的优化不仅要关注各时段对指标的作用，还要关注各阶段对指标的影响。文献[21]指出：间歇过程变量相关关系并非随时间一直变化，而是随过程操作进程或机理特性的变化呈现分段性。在不同的阶段中，变量相关性有着显著的差异，但是在同一个阶段中，不同时刻的过程变量相关关系却近似一致。

过程机理未知，而过程变量的相关性随着阶段的变化而变化，从过程测量数据求得的加权载荷矩阵的变化中可以得出过程的阶段，因而采用对时间片加权载荷矩阵聚类分析结合操作时间的方法划分过程阶段。基于上述分析，采用图 2所示的方法对过程数据进行阶段划分，其详细步骤如下：

1) 划分时间片矩阵；

对于一个包含b条变量曲线，过程周期为t_f的间歇过程，将每条曲线离散化为m个时段变量，则该间歇过程单批次的状态变量曲线可以用矩阵以X(b×m)表示。若考虑到过程的批次n，间歇过程的信息可用三维矩阵以X(n×b×m)表示。将该矩阵沿着时间方向进行垂直切割，得到m个二维矩阵X_l(n×b) (l = 1, 2, …, m)，记二维矩阵为时间片矩阵。

2) 求取加权时间片载荷矩阵；

对划分后的时间片矩阵X_l主元分解：

$ {\mathit{\boldsymbol{X}}_l} = {\mathit{\boldsymbol{U}}_l}\mathit{\boldsymbol{W}}_l^{\rm{T}} = {\mathit{\boldsymbol{u}}_{1, l}}\mathit{\boldsymbol{w}}_{1, l}^{\rm{T}} + {\mathit{\boldsymbol{u}}_{2, l}}\mathit{\boldsymbol{w}}_{2, l}^{\rm{T}} + \cdots + {\mathit{\boldsymbol{u}}_{b, l}}\mathit{\boldsymbol{w}}_{b, l}^{\rm{T}} $

(5)

其中U_l为时间片得分矩阵，W_l为时间片载荷矩阵，它表征了离散化后第l时段过程变量之间的相关性信息，考虑到每个主元的重要性，赋予每个主元不同的权重从而获得加权的载荷矩阵$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}_\mathit{l}}$。

$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}_l} = [{\mathit{\boldsymbol{w}}_{1, l}} \cdot {g_{1, l}}, {\mathit{\boldsymbol{w}}_{2, l}} \cdot {g_{2, l}}, \cdots , {\mathit{\boldsymbol{w}}_{b, l}} \cdot {g_{b, l}}] $

(6)

其中$ {g_{j, l}} = {{{\lambda _{j, l}}} / {\sum\limits_{j = 1}^b {{\lambda _{j, l}}} }}$，λ_{j, l}为l时段第j个主元的方差。$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}_\mathit{l}}$综合考虑了主元的重要性与变量间的相关性，充分利用了得分矩阵与载荷矩阵的信息，其变化代表了过程统计特征的变化。

3) 采用k均值聚类算法对加权时间片载荷矩阵聚类分析从而得到过程的阶段，具体步骤如下：

① 初始化聚类中心W_c^*；

② 将样本$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}_\mathit{l}}$按最小距离原则分配到邻近的聚类；

其中样本$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}_\mathit{l}}$与聚类中心W_c^*欧式距离定义如下：

$ {\rm{Dist}}(l, c) = \left\| {{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}}_l} - \mathit{\boldsymbol{W}}_c^*} \right\|{\rm{ }}(c = 1, 2, \cdots , C) $

(7)

其中c为过程的阶段数。

③ 使用每个聚类中的样本均值更新聚类中心：

$ \mathit{\boldsymbol{W}}_c^* = \frac{1}{{{m_{\rm{c}}}}}\sum\limits_{l = 1}^{{m_{\rm{c}}}} {{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}}_l}} $

(8)

其中m_c是隶属于c阶段的时段数(即第1, 2, …, m_c时段)。

④ 重复步骤②和③直到迭代中两次聚类中心的距离小于一个很小的阈值。

当不同阶段的过程变量有相似的变量相关性时，同一阶段中会包含不同阶段的加权载荷矩阵，此时阶段是不连续的。这种情况下可采用分类结果与操作时间相结合的方法共同识别过程的阶段，即时间上不连续的时段被聚类算法分到同一个阶段中时，结合采样时间仍将它们划分为两个阶段；另外，一个或连续几个时段与前后时段不属于同一个阶段产生“跳跃”现象，通常由测量误差或噪声引起，可当做聚类过程的劣点问题处理，仍将其归结到相邻的阶段。

3.2 各阶段操作曲线的修正

采用上述方法划分完间歇过程的阶段后，基于每一阶段的操作变量对过程进行降维。若分段后的阶段数据包含m_c个时段变量，则n个批次的数据构成了n×m_c维矩阵X_c，$ {\mathit{\boldsymbol{X}}_c} = [{\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}, \cdots , {\mathit{\boldsymbol{x}}_{{m_{\rm{c}}}}}]$；相应的n个批次指标变量，构成n×1维矩阵Q，并令Z_c = [X_cQ]。

当m_c趋向于无穷大时，分段后的轨线可以以任意的形态逼近原动态曲线。为保证离散化后优化命题与原优化命题的等价性，m_c应该取一个较大的值，但m_c过大时，会带来Z_c维数过高的问题。

考虑到Z_c为高维矩阵，为了用较少的变量更直观地反映阶段过程各时段操作条件对指标变量的影响，可使用主元分析对其进行降维处理，因此将$ {\mathit{\boldsymbol{Z}}_c} \in {\mathit{\boldsymbol{R}}^{n \times ({m_{\rm{c}}} + 1)}}$可分解为(m_c+ 1)个向量外积和：

$ {\mathit{\boldsymbol{Z}}_c} = {\mathit{\boldsymbol{D}}_c}{\mathit{\boldsymbol{P}}_c}^{\rm{T}} = {\mathit{\boldsymbol{d}}_1}\mathit{\boldsymbol{p}}_1^{\rm{T}} + {\mathit{\boldsymbol{d}}_2}\mathit{\boldsymbol{p}}_2^{\rm{T}} + \cdots + {\mathit{\boldsymbol{d}}_{({m_c} + 1)}}\mathit{\boldsymbol{p}}_{({m_c} + 1)}^{\rm{T}} $

(9)

其中D_c为得分矩阵，P_c为载荷矩阵，P_c^T为P_c的转置矩阵，$ {\mathit{\boldsymbol{d}}_i} \in {\mathit{\boldsymbol{R}}^n}$(i=1, 2, …, m_c+1)为得分向量(主元变量)，$ {\mathit{\boldsymbol{p}}_i} \in {R^{{m_{\rm{c}}} + 1}}$为载荷向量。各得分向量之间相互正交，各载荷向量之间也为正交关系，且每个载荷向量均为单位长度。

为得出各时段对指标的作用关系，必须将原时段变量表示成主元变量的线性组合，即对主元分解后的数据进行重构。考虑由前r(r≤m_c+1)个主元张成的r维平面，做原变量z_i关于d₁, d₂, …, d_r的线性回归，即求回归方程：

$ {\mathit{\boldsymbol{z}}_i} = {\beta _{1, i}}{\mathit{\boldsymbol{d}}_1} + {\beta _{2, i}}{\mathit{\boldsymbol{d}}_2} + \cdots + {\beta _{r, i}}{\mathit{\boldsymbol{d}}_r} + {\mathit{\boldsymbol{\varepsilon }}_i} $

(10)

其中ε_i是回归误差，$ {\mathit{\boldsymbol{Z}}_c} = [{\mathit{\boldsymbol{z}}_1}, {\mathit{\boldsymbol{z}}_2}, \cdots , {\mathit{\boldsymbol{z}}_{{m_{\rm{c}}} + 1}}]$，记回归系数$ {\mathit{\boldsymbol{\beta }}_i} = {[{\beta _{1, i}}, {\beta _{2, i}}, \cdots , {\beta _{r, i}}]^{\rm{T}}}$，记$ \mathit{\boldsymbol{d}} = \left[ {{\mathit{\boldsymbol{d}}_1}, {\mathit{\boldsymbol{d}}_2}, \cdots , {\mathit{\boldsymbol{d}}_r}} \right]$，则：

$ {\mathit{\boldsymbol{\beta}} _i} = {({\mathit{\boldsymbol{d}}^{\rm{T}}}\mathit{\boldsymbol{d}})^{ - 1}}{\mathit{\boldsymbol{d}}^{\rm{T}}}{\mathit{\boldsymbol{z}}_i} = {\left[ {\begin{array}{*{20}{c}} {\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{d}}_1^{\rm{T}}}\\ \vdots \\ {\mathit{\boldsymbol{d}}_r^{\rm{T}}} \end{array}} \right)}&{\left( {{\mathit{\boldsymbol{d}}_1}, \cdots , {\mathit{\boldsymbol{d}}_r}} \right)} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{d}}_1^{\rm{T}}}\\ \vdots \\ {\mathit{\boldsymbol{d}}_r^{\rm{T}}} \end{array}} \right]\\ {z_i} = \left[ {\begin{array}{*{20}{c}} {{1 / {{{\left\| {{\mathit{\boldsymbol{d}}_1}} \right\|}^2}}}}&{}&0\\ {}& \ddots &{}\\ 0&{}&{{1 / {{{\left\| {{\mathit{\boldsymbol{d}}_r}} \right\|}^2}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{d}}_1^{\rm{T}}{\mathit{\boldsymbol{z}}_i}}\\ \vdots \\ {\mathit{\boldsymbol{d}}_r^{\rm{T}}{\mathit{\boldsymbol{z}}_i}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{d}}_1^{\rm{T}}{\mathit{\boldsymbol{z}}_i}} / {{{\left\| {{\mathit{\boldsymbol{d}}_1}} \right\|}^2}}}}\\ \vdots \\ {{{\mathit{\boldsymbol{d}}_r^{\rm{T}}{\mathit{\boldsymbol{z}}_i}} / {{{\left\| {{\mathit{\boldsymbol{d}}_r}} \right\|}^2}}}} \end{array}} \right] $

(11)

所以：

$ {\beta _{h, i}} = \frac{{\mathit{\boldsymbol{d}}_h^{\rm{T}}{\mathit{\boldsymbol{z}}_i}}}{{{{\left\| {{\mathit{\boldsymbol{d}}_h}} \right\|}^2}}} = \frac{{\left\| {{\mathit{\boldsymbol{z}}_i}} \right\|}}{{\left\| {{\mathit{\boldsymbol{d}}_h}} \right\|}} \times \frac{{\mathit{\boldsymbol{d}}_h^{\rm{T}}{\mathit{\boldsymbol{z}}_i}}}{{\left\| {{\mathit{\boldsymbol{d}}_h}} \right\|\left\| {{\mathit{\boldsymbol{z}}_i}} \right\|}} = \frac{1}{{\sqrt {{\gamma _h}} }}\phi ({\mathit{\boldsymbol{d}}_h}, {\mathit{\boldsymbol{z}}_i}) = {p_{i, h}} $

(12)

式中d_h^T为d_h的转置矩阵，γ_h(h=1, 2, …, r)为主元d_h对应的特征值，ϕ(d_h, z_i)为主元d_h与原变量z_i的相关系数，$ {\mathit{\boldsymbol{p}}_i} = {[{p_{1, i}}, {p_{2, i}}, \cdots , {p_{{m_{\rm{c}}} + 1, i}}]^{\rm{T}}}$。因此回归方程为：

$ {\mathit{\boldsymbol{\hat z}}_i} = {p_{i, 1}}{\mathit{\boldsymbol{d}}_1} + {p_{i, 2}}{\mathit{\boldsymbol{d}}_2} + \cdots + {p_{i, r}}{\mathit{\boldsymbol{d}}_r} $

(13)

主元的个数由累计方差贡献率确定。操作曲线离散化后各时段间变量间时序相关性较大，因而前两个主元包含过程数据大部分信息；实际应用中，平面中可以更直观地得出各时段变量对指标变量的作用关系。因此可以考虑r=2时的特殊情形，此时式(13)简化为：

$ {\mathit{\boldsymbol{\hat z}}_i} = {p_{i, 1}}{\mathit{\boldsymbol{d}}_1} + {p_{i, 2}}{\mathit{\boldsymbol{d}}_2} $

(14)

其中$ {\mathit{\boldsymbol{\hat Z}}_c} = [{\mathit{\boldsymbol{\hat z}}_1}, {\mathit{\boldsymbol{\hat z}}_2}, \cdots , {\mathit{\boldsymbol{\hat z}}_{{m_{\rm{c}}} + 1}}] = [{\mathit{\boldsymbol{\hat x}}_1}, {\mathit{\boldsymbol{\hat x}}_2}, \cdots , {\mathit{\boldsymbol{\hat x}}_{{m_{\rm{c}}}}}, \hat Q]$。

图 3为前两个主元d₁和d₂构成的二维载荷平面，$ {\mathit{\boldsymbol{\hat x}}_{{l_{\rm{c}}}}}$是每阶段重构后的时段变量，其中l_c=1, 1, …, m_c，$ \mathit{\boldsymbol{\hat Q}} $是每阶段重构后的指标变量，θ_lc是重构后阶段时段变量与指标变量间的夹角，它直观的表达了每一阶段中各时段变量与指标变量间的相互关系，基于上述关系给出阶段操作曲线的优化方案^[17]：

$ s_{{l_{\rm c}}}^{{\rm{new}}}{ = _{}}s_{{l_{\rm c}}}^{{\rm{old}}}{{\rm{ + }}_{}}{\sigma _{{l_{\rm c}}}}\cos {\theta _{{l_{\rm c}}}} $

(15)

其中s_lc^new是优化后操作曲线设定值，s_lc^old是优化前操作曲线设定值，σ_lc是时段变量的标准差，表征间歇过程批次差异信息，cosθ_lc是各时段变量与指标变量的余弦相似度，度量时段变量与指标变量的相关性，它与载荷矩阵中的元素的关系如式(16)所示：

$ \cos {\theta _{{l_{\rm c}}}} = \frac{{{p_{{l_{\rm c}}, 1}}{p_{{m_{\rm c}} + 1, 1}} + {p_{{l_{\rm c}}, 2}}{p_{{m_{\rm c}} + 1, 2}}}}{{\sqrt {{{({p_{{l_{\rm c}}, 1}})}^2} + {{({p_{{l_{\rm c}}, 2}})}^2}} \sqrt {{{({p_{{m_{\rm c}} + 1, 1}})}^2} + {{({p_{{m_{\rm c}} + 1, 2}})}^2}} }} $

(16)

操作优化是在原曲线的基础上加上微小的矫正量，使其满足原来优化问题(4)的约束条件。

3.3 多阶段操作曲线的融合

考虑各阶段对指标变量的影响，选择合适的权系数将各阶段操作曲线进行融合。首先对各阶段的数据矩阵X_c主元分解：

$ {\mathit{\boldsymbol{X}}_c} = {\mathit{\boldsymbol{B}}_c}\mathit{\boldsymbol{V}}_c^{\rm{T}} = {\mathit{\boldsymbol{b}}_1}\mathit{\boldsymbol{v}}_1^{\rm{T}} + {\mathit{\boldsymbol{b}}_2}\mathit{\boldsymbol{v}}_2^{\rm{T}} + \cdots + {\mathit{\boldsymbol{b}}_{{m_{\rm{c}}}}}\mathit{\boldsymbol{v}}_{{m_{\rm{c}}}}^{\rm{T}} $

(17)

其中B_c为阶段得分矩阵，V_c为阶段载荷矩阵，V_c^T为V_c的转置矩阵。${\mathit{\boldsymbol{B}}_c} = [{\mathit{\boldsymbol{b}}_1}, {\mathit{\boldsymbol{b}}_2}, \cdots , {\mathit{\boldsymbol{b}}_{{m_{\rm{c}}}}}]$，$ {\mathit{\boldsymbol{V}}_c} = [{\mathit{\boldsymbol{v}}_1}, {\mathit{\boldsymbol{v}}_2}, \cdots , {\mathit{\boldsymbol{v}}_{{m_{\rm{c}}}}}]$。考虑各个主元变量与指标变量间的相关关系赋予不同的权重，得到加权后的得分矩阵$ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over B} _c}$。

$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over B} }}_c} = [{\mathit{\boldsymbol{b}}_1} \cdot {a_1}, {\mathit{\boldsymbol{b}}_2} \cdot {a_2}, \cdots , {\mathit{\boldsymbol{b}}_{{m_{\rm{c}}}}} \cdot {a_{{m_{\rm{c}}}}}] = [{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} }}_1}, {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} }}_2}, \cdots , {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} }}_{{m_{\rm{c}}}}}] $

(18)

其中a_mc是第m_c个主元变量与指标变量间的相关系数的绝对值，加权后的主元变量兼顾了主元本身的方差信息及其与指标变量之间的相关关系，基于上述加权主元变量，定义反应阶段整体信息的阶段代表性主元变量$ {\mathit{\boldsymbol{\tilde b}}_c} = \sum\limits_{k = 1}^{{m_{\rm{c}}}} {{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} }}}_k}} $。

从整个阶段的角度出发，对阶段代表性主元变量$ {\mathit{\boldsymbol{\tilde b}}_c}$与指标变量Q进行相关性分析，以此识别操作优化的关键阶段。采用相关系数的平方CP_c^2[22]衡量阶段过程信息与指标变量的相关度大小，即阶段过程运行行为对指标变量的解释能力：

$ CP_c^2({\mathit{\boldsymbol{\tilde b}}_c}, \mathit{\boldsymbol{Q}}) = {\left( {\frac{{{\mathop{\rm cov}} ({{\mathit{\boldsymbol{\tilde b}}}_c}, \mathit{\boldsymbol{Q}})}}{{\sqrt {D({{\mathit{\boldsymbol{\tilde b}}}_c})D(\mathit{\boldsymbol{Q}})} }}} \right)^2} = \frac{{{{{\mathop{\rm cov}} }^2}({{\mathit{\boldsymbol{\tilde b}}}_c}, \mathit{\boldsymbol{Q}})}}{{D({{\mathit{\boldsymbol{\tilde b}}}_c})D(\mathit{\boldsymbol{Q}})}} $

(19)

其中cov(·)是计算两个变量间协方差的函数，D(·)是求取向量方差的函数。

阶段代表性主元变量综合提取该阶段主要过程波动信息，而计算相关系数揭示该阶段与指标变量间的相关度，CP_c²指标兼顾了上面的两个因素，它是一个表征整个阶段波动信息的阶段综合变量。CP_c²是0到1之间的实数，它的值越大，对应阶段时段变量可以更好解释指标变量波动，该阶段对操作优化更为重要。

基于求取的CP_c²，定义每个阶段的权重系数ω_c，它是CP_c²经过简单的比率运算获得：

$ {\omega _c} = \frac{{CP_c^2}}{{\sum\limits_{c = 1}^C {CP_c^2} }} $

(20)

由式(15)(16)结合阶段的权重系数ω_c可以得出考虑阶段关键程度的操作曲线优化方案：

$ s_{{l_{\rm c}}}^{{\rm{new}}} = s_{{l_{\rm c}}}^{{\rm{old}}}{\rm{ + }}\frac{{C{\omega _{\rm c}}{\sigma _{{l_{\rm c}}}}({p_{{l_{\rm c}}, 1}}{p_{{m_{\rm c}} + 1, 1}} + {p_{{l_{\rm c}}, 2}}{p_{{m_{\rm c}} + 1, 2}})}}{{\sqrt {{{({p_{{l_{\rm c}}, 1}})}^2} + {{({p_{{l_{\rm c}}, 2}})}^2}} \sqrt {{{({p_{{m_{\rm c}} + 1, 1}})}^2} + {{({p_{{m_{\rm c}} + 1, 2}})}^2}} }} $

(21)

为了解释本节所提优化算法，作出如图 4所示多阶段优化策略的流程图。

4 实例研究

以某化工产品间歇结晶过程为研究对象，其原理如图 5所示。间歇结晶过程分为预结晶、降温结晶、升温发汗、全熔融四个阶段。首先保持设定温度为定值，进料从结晶罐上方流入结晶器，然后降低设定温度，使得目标产品和部分杂质从进料溶液中以晶体的形式析出附着在结晶器换热管的外表面，升高温度使得杂质以液态形式流进储料罐并排出，最后再次升高温度使得目标产品达到熔融温度，以液态的形式流入储料罐。当设定温度不同时，最终得到的目标产品纯度差值也不同。

间歇结晶过程是一种多维的动态非线性过程，其操作周期为150 min左右，取1 min为时段间隔。时段变量的取值为每1 min内温度的平均值，指标变量为结晶产品出料与进料的纯度差值。截尾法将全部批次中最短长度数据作为标准，截去多余部分数据，这种方法操作简单，易于实现，是最常用的处理批次数据不等长问题的方法，但它要求变量轨迹在公共部分保持一致且过程的主要信息包含在公共时段内。采集结晶过程正常运行时的130批离散化后的测量温度数据，选择时长为138分钟的最短批次作为标准，其余批次采用截尾的方法与最短长度对齐。对于前50批结晶过程，额外采集罐内温度、进料量、顶部压力等5个变量离散化后的数据，对它们的加权时间片载荷矩阵聚类分析，得出过程的4个阶段，如图 6所示。

由分段结果知，第一阶段为第1至第32 min，第二阶段为第33至第81 min，第三阶段为第82至第99 min，第四阶段为第100至第138 min。对于剩余的80批数据，依据已有的阶段划分结果，在每一阶段中，根据如图 7所示的载荷图，计算温度与纯度差值的相似度从而利用式(15)对操作曲线实现阶段修正。

各阶段的指标CP_c²和权重ω_c取值如表 1所示，由表可知：第2, 3阶段对指标的影响较大，优化时应该给予较大的调整量。利用式(21)可以得到阶段融合后的操作曲线。

图 8为优化前后的操作曲线对比图，其中实线为原操作曲线，虚线为采用文献[17]中的方法优化得出的操作曲线，点虚线为采用文中方法优化得出的操作曲线。为验证优化后操作曲线的有效性，需建立目标产品纯度差值Q关于分段后温度X的偏最小二乘预测模型：

$ \mathit{\boldsymbol{Q}} = \mathit{\boldsymbol{X\eta }} $

(22)

其中η为回归系数矩阵。

表 2记录了预测纯度差值的变化。由表可知：文献[17]和本文所提方法优化出的温度曲线都在原曲线的基础上提升了产品预测纯度差值，而相比于文献[17]，本文所提方法得出的产品预测纯度差值提升量更大。温度影响结晶出料与进料的浓度差，文献[17]考虑各时段变量对指标的作用优化了操作曲线，而实际过程中，降温结晶和升温发汗对纯度差的提升影响更大，这两个阶段在优化时应给予更大的权重。本文在文献[17]的基础上考虑阶段以及阶段重要性的特性，因而纯度差获得了明显的提升。操作曲线多阶段融合优化考虑到过程的阶段特性，综合分析了各时段与整体阶段对指标的作用，更符合过程的实际运行特性，因而能够得到对应更好指标的操作曲线。

5 结论

针对多阶段间歇过程的操作优化问题，提出了一种基于主元相似度的数据驱动型多阶段曲线融合优化方法。在通过对加权时间片载荷矩阵聚类分析的基础上，划分出过程的阶段，再针对各阶段，获取时段变量对指标变量的作用强弱与方向，利用主元相似度对原操作曲线进行修正，最终考虑不同阶段对指标的作用，计算阶段权重并进行阶段曲线融合。仿真结果表明，本文所提方法相比较与已有方法，对优化性能有明显提高。