CN105205224B - 基于模糊曲线分析的时间差高斯过程回归软测量建模方法 - Google Patents

Abstract

本发明公开了一种基于模糊曲线分析的时间差高斯过程回归软测量建模方法，适合应用于具有时滞特性的化工过程。本方法能够从过程历史数据库中提取稳定的时滞信息，引入与主导变量序列更加相关的建模数据。首先，基于模糊曲线分析(FCA)的方法直观判断输入序列对于主导序列的重要性，估计过程时滞参数，用离线条件的时滞参数集对建模数据重构；对于新输入数据，基于一定时刻之前的历史变量值，采用时间差‑高斯过程回归(TDGPR)模型对当前时刻主导变量值在线预测，该方法不存在模型更新的问题，可以很好地追踪输入输出漂移。本发明方法相比于稳态建模方法能够对关键变量进行更精准的预测，从而提高产品质量，降低生产成本。

Description

基于模糊曲线分析的时间差高斯过程回归软测量建模方法

技术领域

本发明涉及基于模糊曲线分析的时间差高斯过程回归(FCA-TDGPR)的软测量建模方法，属于复杂工业过程建模和软测量领域。

背景技术

在实际工业生产过程控制中，需要对一些关键变量进行严格的把关，以满足一系列质量控制要求，在现有的技术条件以及经济代价等问题的制约下直接获取关键变量十分困难。因此，软测量技术应运而生.它通过构造过程易测变量与难测变量之间的数学关系来推断和估计待测主导变量，用软件的方式进行“测量”。

传统的软测量建模方法大多考虑零时延特性，即考虑输入输出具有同样的采样间隔，且在数据库中按时刻一一对应。然而，通过各个传感器采集的输入数据以及通过实验室分析或者在线仪表获得的输出数据间往往存在着显著的时间滞后，如果继续采用稳态情况下的建模方法，建立的模型已经不能完全解释过程的特性，不符合实际过程的因果关系。为了确保软测量模型能够在较长时间内实现关键变量的精确预测，有必要采取措施引入过程的时滞动态信息。

本质上来说，现有的时滞估计方法大多数都在寻找与主导变量最密切相关的辅助变量用于建模，当投入实际应用中时，需要在算法复杂度与算法精度上取得折衷。针对过程时滞信息的估计问题，目前已有的方法包括互信息(Mutual Information，MI)方法、相关系数的方法等。本发明采用模糊曲线分析(Fuzzy Curve Analysis，FCA)的方法将变量时滞信息引入到软测量模型中，这种方法的特点是计算复杂度较低同时易于理解，能够直观有效地确定输入变量的重要性程度。

为了对过程动态实时有效地进行跟踪和控制，软测量模型性能可以通过周期性地重建来维护，主要方法包括滑动窗(Moving Window，MW)方法、迭代方法和实时学习(Just-in-time Learning，JITL)方法，但是这些方法往往需要频繁地对所建模型进行更新，而时间差(Time Difference，TD)模型不但能够处理随时间推移造成的模型性能下降，能取得更新模型般的追踪过程动态的效果，而且能够使得重建模型的可能性最小化。

近年来，基于数据驱动的软测量建模方法得到了越来越多的关注。一些常用的软方法如偏最小二乘(partial least squares，PLS)、主成分分析(principalcomponentanalysis，PCA)等能够很好地处理输入变量和输出变量之间的线性关系，人工神经网络(artificial neural networks，ANN)、支持向量机(support vectormachine，SVM)、最小二乘支持向量机(least support vector squares support vector machine，LS-SVM)能够有效地处理过程的非线性关系。近年来，高斯过程回归(Gaussian processregression，GPR)作为一种非参数概率模型，不仅可以给出预测值，还可以得到预测值的不确定程度。故本发明选择GPR模型作为基本的软测量建模模型，结合TD思想来有效地处理过程输入输出的漂移。

综上所述，建立一种考虑时滞问题的在线软测量模型，对于过程关键变量的严格控制有着重要的意义。

发明内容

对于时延过程，不考虑时滞的软测量模型不能够选取与主导变量序列最相关的辅助变量序列进行建模，采用这样的模型进行估计时，对主导变量的估计精度造成了很大影响。为了有效地提取出过程时滞信息，同时在不频繁更新模型的情况下建立在线的软测量模型，对关键变量的预测提供一种更为有效的在线策略。本发明提供一种基于FCA-TDGPR的软测量建模方法，该方法包括如下步骤：

步骤1：收集过程的输入输出变量数据组成历史训练数据库，获取N组样本{X(t)，y(t)}，t＝1，2，…，N，对数据进行预处理，根据过程机理及经验来确定各辅助变量中存在的最大的时滞参数T_max；

步骤2：对于每个原变量x_i，i∈{1，2，…，m}，分别扩展为含时滞的输入变量集{x_i(t-λ)，λ＝0，1，…，T_max}，扩展方式为：

步骤3：根据模糊曲线分析(FCA)方法确定含时滞输入变量集中每一个变量的重要性程度，确定最优的时滞变量x_i(t-d_i)，确定过程为：

输入变量集{x_i，i＝1，2，…，m}及输出变量y，对于输入变量x_i，t时刻采集的样本值记作x_i(t)，对于(x_i(t)，y(t))，输入变量x_i的模糊隶属度函数定义为：

对于每一个x_i，{Φ_it，y(t)}提供了一条模糊规则，描述为{if x_i isΦ_it(x_i)，theny isy(t)}，Φ_it为变量x_i关于第t个数据点的输入变量模糊隶属度函数，式(2)选取的是高斯模糊隶属度函数，b取变量x_i值域范围的20％；故N个训练样本对应每个变量都有N条模糊规则，在模糊隶属度函数中，每个点对应的{x_i(t)，y(t)}处，有Φ_it＝1；

对于时延过程，通过引入时滞信息，原变量x_i变为T_max+1维，可表示为x_i(t-λ)，λ＝0，1，…，T_max，λ为引入的变量时延值；通过式(3)对扩展后的每个新变量质心去模糊化，可得到第i个变量时延值为λ条件下的模糊曲线C_i，λ；如式(4)所示，d_i为使模糊曲线C_i，λ覆盖范围最大的λ，C_i，λ(λ)_max和C_i，λ(λ)_min为模糊曲线上点值域的最大值和最小值；

若得到的C_i，λ(λ)范围越接近y的范围，那么输入变量x_i(t-λ)的重要程度越高，对C_i，λ(λ)覆盖范围进行排序，可以得到各自的重要性，由此得到最优时滞变量x_i(t-d_i)；

步骤4：利用上一步分析得到的x_i(t-d_i)构成时滞输入集X_d(t)＝[x₁(t-d₁)，x₂(t-d₂)，…，x_m(t-d_m)]^T，重建的软测量训练样本集为{X_d(t)，y(t)}，如果有新的输入样本X(t+1)到来，则基于历史数据库用同样参数进行重组，并转到步骤5，否则，等待新数据到来；

步骤5：对重组训练集、重组的新数据进行j次时间差分处理，其中j的大小可根据过程的主导变量获得周期和性质确定，差分方式为：

然后建立差分输入输出样本之间的高斯过程模型，高斯过程回归(GPR)算法为：

给定训练样本集X∈R^m×N和y∈R^N，m为输入变量维数，N为样本数目，输入和输出之间的关系满足：

y＝f(x)+ε (6)

其中f是未知的函数形式，ε是均值为0，方差为的高斯噪声；对于一个新输入样本，相应的概率预测输出也满足高斯分布，联合高斯分布为

K(X，X)为训练样本间的n维协方差方阵，是测试样本与训练样本的协方差向量，k(x_*，x_*)为测试样本的自协方差值，GPR可以选择不同的协方差函数描述样本分布特征，这里选择高斯协方差函数：

其中，v是控制协方差函数的度量；

高斯过程的超参数可通过极大似然估计得到：

首先将参数Θ_gp设置为一个合理范围内的随机值，然后用共轭梯度法得到优化的参数，获得最优参数后，对于测试样本x_*，可以用式(7)来估计GPR模型的输出值；

步骤6：当t+1时刻的新输入数据到来时，且样本经过时滞信息重组后，采用时间差高斯过程回归(TDGPR)方法，基于y_j(t+1-j)得到预测值y_j，pred(t+1)的计算方式为：

ΔX_d，j(t+1)＝X_d(t+1)-X_d(t+1-j)

Δy_j，pred(t+1)＝f_GPR(ΔX_d，j(t+1)) (10)

y_j，pred(t+1)＝y_j(t+1-j)+Δy_j，pred(t+1)

y_j，pred(t+1)为t+1时刻的预测值。

本发明是通过以下技术方案实现的：

基于FCA-TDGPR的软测量建模方法，所述方法包括以下过程：针对时延过程，首先，对于离线采集的辅助变量和主导变量历史时刻值，搜集足够多的样本数据，组成历史数据库。

在离线阶段通过FCA方法确定过程各辅助变量相对于主导变量序列的时间滞后参数，用于软测量建模数据重构；

然后，对于新的输入数据，基于一定时刻之前采集的历史变量值，采用时间差-高斯过程回归(TDGPR)模型对当前时刻主导变量值进行在线预测，从而实现对于关键变量的实时估计和控制，得到更加精确的结果，从而提高产量，降低生产成本。本发明提供的基于模糊曲线分析的时间差高斯过程回归软测量建模方法，以一种直观有效的方式，同时在计算复杂度低的情况下，可以从过程历史数据库中提取变量的时滞信息用于软测量建模数据重构，校正了输入输出间实际的因果对应关系，同时将在线软测量模型更新的可能性最小化，不存在更新模型的问题，采用GPR对输入输出数据的漂移进行追踪，能够很好地提取过程时滞信息，相对于传统的不考虑时滞的在线模型，能够得到更加精确的预测结果。

附图说明

图1为基于FCA-TDGPR的在线软测量方法建模流程图；

图2为脱丁烷塔过程示意图；

图3为TDGPR建模思路示意图；

图4为原变量与最优时滞变量的模糊曲线分布图；

图5为不同j值情况下的丁烷浓度预测结果散点图。

具体实施方式

下面结合图1所示的建模流程图，对本发明做进一步详述：

以实际化工过程为例，脱丁烷塔过程是石油炼制生产过程中脱硫和石脑油分离装置的重要组成部分，需要控制的主导变量之一为塔底的丁烷(C₄)浓度，过程示意图如图2所示。由于丁烷浓度值不能直接检测得到，分析与获得存在滞后问题，同时不同辅助变量有着不同程度的滞后。实验数据来源于实际工业过程，该数据包含样本2394组，共有7个辅助变量，如图2中标注所示：x₁为塔顶温度；x₂为塔顶压力；x₃为塔顶回流量；x₄为塔顶产品流出量；x₅为第6层塔板温度；x₆为塔底温度1；x₇为塔底温度2，1个主导变量是塔底丁烷浓度，将塔底丁烷浓度作为过程关键变量进行预测。

步骤1：收集历史输入输出数据组成包含N组连续样本的训练数据库，假设其按时刻表示为{X(t)，y(t)}，t＝1，2，…，N，对数据进行预处理，将2个塔底温度变量取平均后作为1个辅助变量，那么X(t)＝[x₁(t)，x₂(t)，x₃(t)，x₄(t)，x₅(t)，x₆(t)]^T，设置6个变量中存在的最大时滞T_max参数值为19；

步骤2：对于每个原始变量x_i，i∈{1，2，…，6}，按式(1)方式分别扩展为T_max+1维含时滞的输入变量{x_i(t-λ)，λ＝0，1，…，T_max}，共得到120维待分析的时滞变量集；

步骤3：根据FCA方法确定含时滞输入变量集中每一个变量的重要性程度：

对于(x_i(t)，y(t))，输入变量x_i的模糊隶属度函数定义为：

对于每一个x_i，{Φ_it，y(t)}提供了一条模糊规则，描述为{if x_i is Φ_it(x_i)，then y is y(t)}，Φ_it为变量x_i关于第t个数据点的输入变量模糊隶属度函数，式(12)选取的是高斯模糊隶属度函数，b取变量x_i值域范围的20％。故N个训练样本对应每个变量都有N条模糊规则。在模糊隶属度函数中，每个点对应的{x_i(t)，y(t)}处，有Φ_it＝1。

对于时延过程，通过引入时滞信息，原变量扩展后表示为x_i(t-λ)，λ＝0，1，…，T_max，λ为引入的变量时延值。通过式(13)对扩展后的每个新变量质心去模糊化，可得到第i个变量时延值为λ条件下的模糊曲线C_i，λ。如式(14)所示，d_i为使模糊曲线C_i，λ覆盖范围最大的λ，C_i，λ(λ)_max和C_i，λ(λ)_min为模糊曲线上点值域的最大值和最小值。

若得到的C_i，λ(λ)范围越接近y的范围，那么输入变量x_i(t-λ)的重要程度越高，对C_i，λ(λ)覆盖范围进行排序，可以得到各自的重要性。基于该原理，本发明利用FCA方法提取过程各变量的最优时滞参数d_i并确定最优的时滞变量x_i(t-d_i)，进行软测量建模数据重构。

步骤4：利用上一步分析得到的时滞参数d₁，d₂，…，d_m重组用于在线建模的软测量训练输入样本集，X_d(t)＝[x₁(t-d₁)，x₂(t-d₂)，x₃(t-d₃)，x₄(t-d₄)，x₅(t-d₅)，x₆(t-d₆)]。如果有新的输入样本X(t+1)到来，则基于历史数据库用同样参数进行重组，并转到步骤5，否则，等待新数据到来。

步骤5：对重组后的建模计算时间差分，如式(15)所示，j为差分阶次(j的大小可根据过程的主导变量获得周期和性质确定)：

对于ΔX_d，j(t)与Δy_j(t)采用GPR的方法对其之间关系进行回归，满足Δy_j(t)＝f(ΔX_d，j(t))+e(t)，GPR方法能够通过给定的训练样本的输入输出数据得到映射关系，便可由新的输入数据得到相应的预测值和不确定程度，获得具有概率意义的输出结果。GPR算法为：

通常样本输出观测值y和噪声ε满足式(16)关系：

若确定均值函数和协方差函数，高斯过程的分布就能唯一确定。为了方便，通常将均值函数预处理为0。协方差函数能够把输出间的相关关系转化为输入数据之间的函数关系，由于相近的输入产生相近的输出，协方差函数的选择可以根据样本分布的特征选取，要符合距离相近的样本间相关性大，反之相关性小的特征。本文选择的协方差函数形式如式(17)所示：

式中，x_p，x_q∈R^m，v控制协方差函数的量度，π_d刻画了每个x^d的相对重要性。高斯过程的超参数的确定一般通过MLE方法对式(18)进行估计，超参数的优化，可以通过共轭梯度法实现。基于测试样本和训练数据，可以计算出测试数据x_*预测值的后验分布式服从式(19)的联合高斯分布，K(X，X)为训练样本间的n维协方差方阵，k(x_*，X)是测试样本与训练样本的协方差向量，k(x_*，x_*)为测试样本的自协方差值，f_gp为GPR预测值。

当t+1时刻的新数据到来时，采用TDGPR方法得到的预测值y_j，pred(t+1)计算方式为：

ΔX_j(t+1)＝X(t+1)-X(t+1-j)

Δy_j，pred(t+1)＝f_GPR(ΔX_j(t+1)) (20)

y_j，pred(t+1)＝y_j(t+1-j)+Δy_j，pred(t+1)

在实际工业过程中，会出现仪表损坏或实验室分析滞后的情况，主导变量的获得时间间隔很大且数量很少，以及数据库中已有主导变量分析值缺失等情况。这样一来，由图3所示，对于新到来的测试数据X(t+1)，以j时刻前数据库保存的y_j(t+1-j)为基准，可以得到t+1时刻的主导变量预测值，在线模型得到的预测输出y_j，pred(t+1)由式(20)计算，即可以得到塔底丁烷浓度的预测结果。

如图4所示，与未引入时滞的原变量相比，6个重构变量对于主导变量的贡献更大，为在线建模引入了更加相关的建模数据。同时，为了验证本发明用于在线估计的有效性，选取2394组数据中的前1519组重构1500组训练样本，后续875组用于在线测试样本，对丁烷浓度进行了软测量在线预测。

图5为经过时滞估计的FCA-TDGPR(本发明方法)与不经过时滞估计的t-TDGPR方法对于丁烷浓度的预测结果散点图。

由图5可知，当时间差j从1到10逐渐增加时，基于历史数据库进行预测的时间间隔逐步增大，预测精度不断下降，这是因为基于越近的分析值做出的预测能够越好地追踪当前时刻的过程动态。

尽管两种方法精度都在下降，但是本发明的预测结果相比于不考虑时滞的TDGPR方法，在时间差递增时也能较好地贴近丁烷浓度真值，说明了提取的时滞信息符合过程的实际因果关系，考虑了变量时滞估计的软测量模型预测精度更加精确。

而且验证了通过模糊曲线分析方法确定最优时滞参数后，重构的数据对主导变量贡献增强，因此在线模型的精度也显著增强，同时反映了GPR方法能够很好地解释过程的动态变化，基于TDGPR方法的在线软测量模型能够自适应地基于j时刻前的历史变量值估计当前时刻的丁烷浓度。

图4和图5共同验证了基于模糊曲线分析的时间差高斯过程回归的软测量建模方法对于塔底丁烷浓度的在线预测具有良好的精度。

Claims (1)

1.基于模糊曲线分析的时间差高斯过程回归软测量建模方法，其特征在于，该方法步骤为：

步骤1：收集过程的输入输出变量数据组成历史训练数据库，获取N组样本{X(t),y(t)},t＝1,2,…,N,对数据进行预处理，根据过程机理及经验来确定各辅助变量中存在的最大的时滞参数T_max；

步骤2：对于每个原变量x_i,i∈{1,2,…,m}，分别扩展为含时滞的输入变量集{x_i(t-λ),λ＝0,1,…,T_max}，扩展方式为：

步骤3：根据模糊曲线分析(FCA)方法确定含时滞输入变量集中每一个变量的重要性程度，确定最优的时滞变量x_i(t-d_i)，确定过程为：

输入变量集{x_i,i＝1,2,…,m}及输出变量y，对于输入变量x_i,t时刻采集的样本值记作x_i(t),对于(x_i(t),y(t))，输入变量x_i的模糊隶属度函数定义为：

对于每一个x_i，{Φ_it,y(t)}提供了一条模糊规则,描述为{if x_i isΦ_it(x_i)，then y isy(t)}，Φ_it为变量x_i关于第t个数据点的输入变量模糊隶属度函数，式(2)选取的是高斯模糊隶属度函数，b取变量x_i值域范围的20％；故N个训练样本对应每个变量都有N条模糊规则，在模糊隶属度函数中，每个点对应的{x_i(t),y(t)}处，有Φ_it＝1；

对于时延过程，通过引入时滞信息，原变量x_i变为T_max+1维，可表示为x_i(t-λ),λ＝0,1,…,T_max，λ为引入的变量时延值；通过式(3)对扩展后的每个新变量质心去模糊化，可得到第i个变量时延值为λ条件下的模糊曲线C_i,λ；如式(4)所示，d_i为使模糊曲线C_i,λ覆盖范围最大的λ，C_i,λ(λ)_max和C_i,λ(λ)_min为模糊曲线上点值域的最大值和最小值；

若得到的C_i,λ(λ)范围越接近y的范围，那么输入变量x_i(t-λ)的重要程度越高，对C_i,λ(λ)覆盖范围进行排序，可以得到各自的重要性，由此得到最优时滞变量x_i(t-d_i)；

步骤4：利用上一步分析得到的x_i(t-d_i)构成时滞输入集X_d(t)＝[x₁(t-d₁),x₂(t-d₂),…,x_m(t-d_m)]^T，重建的软测量训练样本集为{X_d(t),y(t)}，如果有新的输入样本X(t+1)到来，则基于历史数据库用同样参数进行重组，并转到步骤5，否则，等待新数据到来；

步骤5：对重组训练集、重组的新数据进行j次时间差分处理，其中j的大小可根据过程的主导变量获得周期和性质确定，差分方式为：

然后建立差分输入输出样本之间的高斯过程模型，高斯过程回归GPR算法为：

给定训练样本集X∈R^m×N和y∈R^N，m为输入变量维数，N为样本数目，输入和输出之间的关系满足：

y＝f(x)+ε (6)

其中f是未知的函数形式，ε是均值为0，方差为的高斯噪声；对于一个新输入样本，相应的概率预测输出也满足高斯分布，联合高斯分布为

K(X,X)为训练样本间的n维协方差方阵，k(x_*,X)＝k(X,x_*)^T是测试样本与训练样本的协方差向量，k(x_*,x_*)为测试样本的自协方差值，GPR可以选择不同的协方差函数描述样本分布特征，这里选择高斯协方差函数：

其中，ν是控制协方差函数的度量；

高斯过程的超参数可通过极大似然估计得到：

首先将参数Θ_gp设置为一个合理范围内的随机值，然后用共轭梯度法得到优化的参数，获得最优参数后，对于测试样本x_*，可以用式(7)来估计GPR模型的输出值；

步骤6：当t+1时刻的新输入数据到来时，且样本经过时滞信息重组后，采用时间差高斯过程回归(TDGPR)方法，基于y_j(t+1-j)得到预测值y_j,pred(t+1)的计算方式为：

y_j,pred(t+1)为t+1时刻的预测值。