CN103279030B - 基于贝叶斯框架的动态软测量建模方法及装置 - Google Patents

Abstract

本发明公开一种基于贝叶斯框架的动态软测量建模方法及装置，该方法包括：步骤A：建立一阶脉冲响应模型以及支持向量机模型步骤B：测量并输入u_k(…,t)(k=1,…,m)，计算出状态变量x(t)并计算出y(t)；步骤C：输出y(t)。

Description

基于贝叶斯框架的动态软测量建模方法及装置

技术领域

本发明涉及计算机领域，尤其涉及工业生产过程中的基于贝叶斯框架的动态软测量建模方法及装置。

背景技术

在工业生产控制过程中，质量变量通常是由与产品的质量相关的指标构成，例如柴油的干点、聚丙烯的熔融指数，化学反应器中的pH值等。然而，质量变量通常无法进行在线测量，只能够通过长达几小时的实验室分析化验来得到，而实验室的化验数据延迟太大，无法用来作为控制回路的反馈信号，高品质的控制要求得不到满足；即便某些质量变量有可供选择的在线分析仪器，但这些仪器通常价格昂贵且后期维护费用高。基于上述原因，在化工过程、生物制药、钢铁锻造等领域，普遍采用软测量方法通过建立二次变量与质量变量之间的数学关系建立模型，实现对难测变量的实时预报输出。

目前的软测量方法，如主成分分析(PCA)、偏最小二乘(PLS)、人工神经网络(ANN)、支持向量机(SVM)等一大批基于统计学习理论的建模方法被应用到软测量建模中，并且取得了一定的成果。但是存在以下问题：

第一：通常只考虑与质量数据点时序匹配的过程数据点，然而实际的工业处理过程大多数都处在动态变化当中，当前时刻系统的某个状态必然是由之前一段时间的状态决定的，而不是某一个时间点的状态单独决定的，故测量结果不够精确；

第二：忽略了过程时延或即便考虑了过程时延的因素，单纯使用某个过程数据点进行建模的方法，由于缺乏对过程动态特性的描述，使得软测量模型的适用性仅仅局限在稳态工况下，一旦出现工作点变化或者扰动发生，模型的跟踪能力会大幅度降低。被利用的过程数据仅仅占全部数据的很小一部分，大部分含有丰富过程信息的数据被直接丢弃处理，故测量结果不够精确；

第三，由于静态软测量模型在二次变量发生剧烈变化时质量变量也会剧烈变化，因此从控制的角度出发，当软测量的输出作为控制回路的反馈信号时，控制性能会有较大幅度的恶化，最终影响测量的反馈结果。

发明内容

(一)发明目的

本发明旨在提供一种高精度的测量结果的基于贝叶斯框架的动态软测量建模方法及装置。

(二)技术方案

为达上述目的，本发明基于贝叶斯框架的动态软测量建模方法，包括以下步骤：

步骤A：建立一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{g}{({\mathrm{\α}}_1,{\mathrm{\τ}}_1)}^T{u}_1(\mathrm{...},t)\\ \stackrel{\→}{g}{({\mathrm{\α}}_2,{\mathrm{\τ}}_2)}^T{u}_2(\mathrm{...},t)\\ .\\ .\\ .\\ \stackrel{\→}{g}{({\mathrm{\α}}_m,{\mathrm{\τ}}_m)}^T{u}_m(\mathrm{...},t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)(k＝1,…,m)为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的向量，x_k(t)(k＝1,…,m)为状态变量的第k个分量在采样时间t的取值,模型为一阶脉冲响应模板；且

其中，α_k(k＝1,…,m)为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；ΔT为采样周期，τ_k/ΔTzeros表示的0的个数；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\σ}}^2}\}+b---\left(3\right)$

其中α_i,与b为支持向量机参数，σ为RBF核参数，为支持向量，l为支持向量的个数；

步骤B：测量u_k(…,t)(k＝1,…,m)并输入一阶脉冲响应模型，计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)；

步骤C：输出y(t)。

优选地，所述步骤A还包括以下具体步骤：

步骤A1：构建如公式(1)以及公式(3)所述的关于u_k(…,t)(k＝1,…,m)以及y(t)的系统结构关系；

步骤A2：获取与u_k(…,t)(k＝1,…,m)相对应的数据U_k(…,t)(k＝1,…,m)、以及与U_k(…,t)(k＝1,…,m)相映射的Y(t)构成的数据样本，并按采样时间形成递增的时间序列；

步骤A3：对数据样本中的数据进行标准化处理，以使数据在每一个维度上具有零均值和单位方差；

步骤A4：设定一阶脉冲响应模型以及支持向量机模型训练的初始参数α_k,τ_k、脉冲响应序列长度L、训练参数ε、C、RBF核参数σ以及迭代停止阈值μ：

步骤A5：将数据样本中的u_k(…,t)(k＝1,…,m)输入公式(1)，且将公式(1)的结果输入公式(3)，比较公式(3)的输出与Y(t)，并根据比较结果校正公式(1)以及公式(3)的参数。

为达上述目的，本发明基于贝叶斯框架的动态软测量建模装置，所述基于贝叶斯框架的动态软测量建模装置至少包括模型建立单元、数据采集计算单元以及数据输出单元；

所述模型建立单元，用以一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{g}{({\mathrm{\α}}_1,{\mathrm{\τ}}_1)}^T{u}_1(\mathrm{...},t)\\ \stackrel{\→}{g}{({\mathrm{\α}}_2,{\mathrm{\τ}}_2)}^T{u}_2(\mathrm{...},t)\\ .\\ .\\ .\\ \stackrel{\→}{g}{({\mathrm{\α}}_m,{\mathrm{\τ}}_m)}^T{u}_m(\mathrm{...},t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)(k＝1,…,m)为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的向量，x_k(t)(k＝1,…,m)为状态变量的第k个分量在采样时间t的取值,模型为一阶脉冲响应模板；且

其中，α_k(k＝1,…,m)为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；ΔT为采样周期，τ_k/ΔTzeros表示的0的个数；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\σ}}^2}\}+b---\left(3\right)$

其中α_i 与b为支持向量机参数，σ为RBF核参数，为支持向量,l为支持向量的个数；

所述数据采集计算单元，用以采集在线测量值u_k(…,t)(k＝1,…,m)，并输入到一阶脉冲响应模型计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)。

进一步地，所述模型建立单元包括：

结构建立模块，用以构建如公式(1)以及公式(3)所述的关于u_k(…,t)(k＝1,…,m)以及y(t)的系统结构关系；

样本数据获取模块，用以获取与u_k(…,t)(k＝1,…,m)相对应的数据U_k(…,t)(k＝1,…,m)、以及与U_k(…,t)(k＝1,…,m)相映射的Y(t)构成的数据样本，并按采样时间形成递增的时间序列；

数据样本标准化模块，用以对数据样本中的数据进行标准化处理，以使数据在每一个维度上具有零均值和单位方差；

初始参数设定模块，用以设定一阶脉冲响应模型以及支持向量机模型训练的初始参数α_k,τ_k、脉冲响应序列长度L、RBF核参数σ以及迭代停止阈值μ：

参数训练模块，用以将数据样本中的U_k(…,t)(k＝1,…,m)输入公式(1)，且将公式(1)的结果输入公式(3)，比较公式(3)的输出与Y(t)，并根据比较结果校正公式(1)以及公式(3)的参数。

(三)本发明基于贝叶斯框架的动态软测量建模方法及装置的有益效果

第一：本发明基于贝叶斯框架的动态软测量建模方法及装置，通过将一阶脉冲响应模型与支持向量机模型的结合，以一段连续采样周期的采用组成的向量的为输入，从而打破了传统方法以及装置仅以单独的时间点来演算难测的主导变量即质量变量y(t)的局限，从而能够对过程的动态信息进行有效描述，利用了生产过程中的动态信息，与传统的方法以及装置相比估计精度更高，在实施闭环反馈控制时也具有更加优良的控制性能；

第二：本发明所述的基于贝叶斯框架的动态软测量建模方法及装置，充分利用支持向量机模型自身对时延的动态调整以及向量的时间特性，从而能对输出时延、动态参数进行智能寻优，减少了现有算法对于过程动态先验信息的依赖，解决了现有部分动态软测量模型由于输出时延信息无法获取导致无法使用的问题；

第三：本发明所述的基于贝叶斯框架的动态软测量建模方法及装置，对一阶脉冲响应以及支持向量机模型的训练能够在不依赖传统测试数据集的基础上，根据客观准则选择模型训练中的各类参数，自动惩罚过拟合的模型，有效地改进由于非线性数据驱动模型过学习导致的模型泛化性能问题。

附图说明

图1为本发明实施例一所述的基于贝叶斯框架的动态软测量建模方法的原理图；

图2～图5为本发明实施例三所述基于贝叶斯框架的动态软测量建模方法的测量结果与测量结果的比较图；

图6是本发明第二实施例所述基于贝叶斯框架的动态软测量建模装置结构示意图。

具体实施方式

下面结合说明书附图以及实施例对本发明基于贝叶斯框架的动态软测量建模方法做进一步的说明。

实施例一：

本实施例基于贝叶斯框架的动态软测量建模方法包括以下步骤：

步骤A：建立一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{g}{({\mathrm{\α}}_1,{\mathrm{\τ}}_1)}^T{u}_1(\mathrm{...},t)\\ \stackrel{\→}{g}{({\mathrm{\α}}_2,{\mathrm{\τ}}_2)}^T{u}_2(\mathrm{...},t)\\ .\\ .\\ .\\ \stackrel{\→}{g}{({\mathrm{\α}}_m,{\mathrm{\τ}}_m)}^T{u}_m(\mathrm{...},t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)(k＝1,…,m)为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的向量，x_k(t)(k＝1,…,m)为状态变量的第k个分量在采样时间t的取值,模型为一阶脉冲响应模板；且

其中，α_k(k＝1,…,m)为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；ΔT为采样周期；所述的一阶脉冲响应模型是动态系统，用于输入从实际测量的系统中动态变量，且该变量通常是容易测量的过程变量，且以向量的形式输入，从而能表征连续一段时间内被测量的系统内的变化；再通过公式(1)得到中间的状态变量x(t)用以输入后续的支持向量机模型；τ_k/ΔTzeros代表的是大括号所在位置所包含的0的个数；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\σ}}^2}\}+b---\left(3\right)$

其中α_i,与b为支持向量机参数，σ为RBF核参数，为支持向量机的支持向量，l为支持向量的个数；支持向量机模型为非线性的静态模型，其输入是一个m维的向量，通过公式(3)的计算将得出难以测量的主导变量y(t)，y(t)通常为质量变量；

步骤B：输入在线测量值u_k(…,t)(k＝1,…,m)，计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)；

步骤C：输出y(t)。在具体的实践过程中，通常y(t)最终将写入到的是DCS分布式显示系统中，再有工作系统或工作人员参考输入，用于实际操作以及控制。

图1则表示的本实施例所述的基于一阶脉冲响应和支持向量机的动态软测模型的原理图。

本实施例所述的基于贝叶斯框架的动态软测量建模方法，将一阶脉冲响应模型与支持向量机模型相结合，求取y(t)时以被测系统中一段连续时间的测量的离散变量为输入，从而求得输出不在是仅单个的时间状态点相关联，从而更加符合实际生产过程，从而精确度更高，控制更加精确；此外充分利用了支持向量机模型自身考虑了时延的特性以及向量输入的时间特性，从而解决了现有技术当中的忽略过程时延或即便考虑了过程时延的因素，单纯使用某个过程数据点进行建模的方法导致的测量结果不够精确的缺点。

实施例二：

本实施例在实施例一的基础上，具体化了所述步骤A；

步骤A1：构建如公式(1)以及公式(3)所述的关于u_k(…,t)(k＝1,…,m)以及y(t)的系统结构关系；

步骤A2：获取与u_k(…,t)(k＝1,…,m)相对应的数据U_k(…,t)(k＝1,…,m)以及与U_k(…,t)(k＝1,…,m)相映射的Y(t)构成的数据样本，并按采样时间形成递增的时间序列；通常在采集所述的U_k(…,t)(k＝1,…,m)时遵照香农定理；

步骤A3：对数据样本中的数据进行标准化处理，以使数据在每一个维度上具有零均值和单位方差；

步骤A4：设定一阶脉冲响应模型以及支持向量机模型训练的初始参数α_k,τ_k、脉冲响应序列长度L、训练参数ε、C、RBF核参数σ以及迭代停止阈值μ：

步骤A5：将数据样本中的U_k(…,t)(k＝1,…,m)输入公式(1)，且将公式(1)的结果输入公式(3)，比较公式(3)的输出与Y(t)，并根据比较结果校正公式(1)以及公式(3)的参数。

通过训练调整公式(1)以及公式(3)中的参数，使得得出的计算模型能精确的用于难测的主导变量y(t)的计算。

实施例三：

本实施例在实施例二的基础上，进一步地详细了步骤A.

如图2所示，步骤S1：在生产过程机理分析或实际经验的基础上构建一个难测主导变量和多个易测辅助变量之间的系统结构关系，即

u_k(…,t)(k＝1,…,m)与y(t)之间的系统结构关系；并给出系统的参数约束范围。

系统结构关系为：系统输出y(t)代表难测主导变量，系统输入u_k(…,t)(k＝1,…,m)代表相关的易测辅助变量。系统由动态环节和非线性静态环节串联而成，而动态环节表现为一阶脉冲响应模型的形式：

其中，m为输入变量的维数，也为状态变量x(t)的维数；为输入变

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{g}{({\mathrm{\α}}_1,{\mathrm{\τ}}_1)}^T{u}_1(\mathrm{...},t)\\ \stackrel{\→}{g}{({\mathrm{\α}}_2,{\mathrm{\τ}}_2)}^T{u}_2(\mathrm{...},t)\\ .\\ .\\ .\\ \stackrel{\→}{g}{({\mathrm{\α}}_m,{\mathrm{\τ}}_m)}^T{u}_m(\mathrm{...},t)\end{array}\right]$

量的第k个分量在采样时间t之前L个采样周期内每采用周期一个采样数据组成的采样值组成的向量，x_k(t)(k＝1,…,m)为状态变量的第k个分量在采样时间t的取值。模型为一阶脉冲响应模板，即

其中，α_k(k＝1,…,m)为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；ΔT为过程的采样周期。要求根据先验知识，给出各过程子系统的纯时延常数τ_k的大致范围0≤τ_k≤d_k，d_k为第k个输入通道可能出现的最大时延。

非线性静态环节表现为支持向量机模型的形式：

$y\left(t\right)=f\left(x\right(t\left)\right)=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\σ}}^2}\}+b$

其中α_i,与b为支持向量机参数，σ为RBF核参数，为支持向量机的支持向量。

步骤S2：完成相关数据样本的收集。从在线采集的数据库以及离线化验数据库中搜集相关的易测辅助变量与难测主导变量数据构成一组数据样本，并按采样时间递增的顺序形成相应的时间序列。

难测主导变量采样时间序列表示为：[y(P₁),y(P₂),...,y(P_l)]

设采集的样本数目为l，P_i为难测主导变量第i个离线采样数据的

u_k(…,P_i)＝[u_k(P_i-L+1),u_k(P_i-L+2),...,u_k(P_i)],i＝1,…,l,k＝1,…,m采样时刻，y(P_i)(i＝1,…,l)为难测主导变量第i个离线采样数据，采样周期一般较长，并且通常具有不规则特征。

m个易测辅助变量采样时间序列表示为：

其中u_k(P_i)为第k个易测辅助变量P_i时刻的在线采样数据，采样周期应满足香农采样定理。

步骤S3：对输入样本数据进行标准化处理，使输入数据在每一个维度上具有零均值和单位方差。

步骤S4：设定模型训练的初始参数：所有动态子系统的α_k,τ_k均设为0，表征静态软测量。设置脉冲响应序列长度L，根据难测变量的变化范围以及精度要求设定合适支持向量机训练参数ε、C及RBF核参数σ；设置迭代停止阈值μ，初始化迭代次数ite＝0。

步骤S5：利用处理后的该组数据样本，基于脉冲响应序列计算状态变量x(t)，利用初始参数训练支持向量机。

步骤S6：令ite＝ite+1。根据如下式子更新最优的支持向量机惩罚参数C^MP。

$C^{MP}=\frac{\underset{i=1}{\overset{l}{\Σ}}\underset{j=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)({\mathrm{\α}}_j^{*}-{\mathrm{\α}}_j)\exp \{-\frac{\left|\right|x\left(P_i\right)-x\left(P_j\right)|{|}^2}{2{\mathrm{\σ}}^2}\}}{\mathrm{\γ}}$

而γ的计算方法如下：

$\mathrm{\γ}=\underset{i=1}{\overset{l}{\Σ}}\frac{{\mathrm{\ρ}}_i}{1/C+{\mathrm{\ρ}}_i}$

其中ρ_i为矩阵的特征值，矩阵为l×l维矩阵，第i行第j列的元素的计算如下：

$r_i\exp \{-\frac{\left|\right|x\left(P_i\right),x\left(P_j\right)|{|}^2}{2{\mathrm{\σ}}^2}\}$

其中r_i＝r(y_i-f(x_i)-ε)+r(f(x_i)-y_i-ε)，r(u)＝{u[1-2s(u)]+2}·s(u)[1-s(u)]，s(u)为sigmoid函数1/(1+e^-u)。

步骤S7：利用更新的参数训练支持向量机。

步骤S8：根据下式计算最优的核参数σ^MP，并在之后利用新的参数C^MP再进行一次的训练：

${\mathrm{\&theta;}}^{MP}=\left\{\begin{array}{cc}\sqrt{\mathrm{\&eta;}}\mathrm{\&theta;},& 0<\mathrm{\&eta;}<1\\ \mathrm{\&eta;}\mathrm{\&theta;},& \mathrm{\&eta;}\&GreaterEqual;1\end{array}\right.$

θ＝1/σ，，σ^MP＝1/θ^MP，

η的计算式如下

$\mathrm{\&eta;}=\left|\begin{array}{c}{\mathrm{\&lambda;}}^{MP}\underset{i,j=1}{\overset{l}{\&Sigma;}}(a_i-a_j)(a_i^{*}-a_j^{*}){\left(x\right(P_i)-x(P_j\left)\right)}^2\\ \frac{\&times;\exp (-{\mathrm{\&theta;}}^2{\left(x\right(P_i)-x(P_j\left)\right)}^2/2)}{{\mathrm{\&lambda;}}^{MP}/(l-{\mathrm{\&lambda;}}^{MP}trace({\stackrel{\&OverBar;}{A}}^{-1}))\&times;tace\left({\stackrel{\&OverBar;}{A}}^{-2}\right(\&part;\stackrel{\&OverBar;}{K}/\&part;\mathrm{\&theta;}))}\\ +trace\left({\stackrel{\&OverBar;}{A}}^{-1}\right(\&part;\stackrel{\&OverBar;}{K}/\&part;\mathrm{\&theta;}))/\mathrm{\&theta;}\end{array}\right|$

其中，

$\stackrel{\&OverBar;}{A}=\frac{1}{C}{I}_{l\&times;l}+\stackrel{\&OverBar;}{K}$

矩阵的定义见步骤S6。

步骤S10：根据如下方法分别求取m个子系统最优的动态参数 $\{{\mathrm{\&alpha;}}_k^{MP},{\mathrm{\&tau;}}_k^{MP}\},k=1,\mathrm{...},m:$

定义损失函数如下：

$F({\mathrm{\&alpha;}}_k,{\mathrm{\&tau;}}_k)=\underset{i=1}{\overset{l}{\&Sigma;}}{\&lsqb;y\left(P_i\right)-f\left(x\right(P_i))\&rsqb;}^2,0\&le;{\mathrm{\&alpha;}}_k<1,k=1,...,m$

$x\left(P_i\right)=\left[\begin{array}{c}{x}_1\left(P_i\right)\\ x_2\left(P_i\right)\\ .\\ .\\ .\\ x_m\left(P_i\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(\mathrm{...},P_i)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(\mathrm{...},P_i)\\ .\\ .\\ .\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(\mathrm{...},P_i)\end{array}\right]$

对于第j个子系统(j＝1,…,m)，其纯时延参数由均匀分布集合{0,1,…,d_j}确定，其中。对于每一个可能的纯时延利用黄金分割法在0≤α_k＜1范围内优化损失函数寻找最优的

而在得到后，最优的由下式确定：

${\mathrm{\&tau;}}_k^{MP}=\arg \underset{0\&le;{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k\&le;d_k}{\min }F\left({\mathrm{\&alpha;}}_k^{MP}\right({\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k),{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k)$

步骤S11：利用更新的参数训练支持向量机。

步骤S12：利用下式计算其evidence：

$evid\left(ite\right)=-\frac{1}{C}\underset{i=1}{\overset{l}{\&Sigma;}}\underset{j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)({\mathrm{\&alpha;}}_j^{*}-{\mathrm{\&alpha;}}_j)\exp \{-\frac{\left|\right|x\left(P_i\right)-x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}-\underset{i=1}{\overset{l}{\&Sigma;}}(1+{\mathrm{C\&rho;}}_i)$

其中ρ_i为矩阵的特征值，矩阵的定义见步骤S6。

步骤S13：当evid(ite)-evid(ite-1)＜μ或者evid(ite)＜evid(ite-1)时停止迭代，此时的模型即为所求的动态软测量模型；否则，返回步骤S6继续迭代。

在本实施例中，采用利用了贝叶斯框架，能够在不依赖传统测试数据集的基础上，根据客观准则选择模型训练中的各类参数，自动惩罚过拟合的模型，有效地改进由于非线性数据驱动模型过学习导致的模型泛化性能问题。

为了进一步说明本发明所述的基于贝叶斯框架的动态软测量建模方法的效果下面进行一个实例验证；

选定三输入-单输出生成过程系统作为本实例需要仿真的对象，且已经知道该三输入-单输出的系统的状态方程与输出方程，且分布为：

$\left\{\begin{array}{c}{x}_1\left(s\right)=\frac{3}{600s+1}{e}^{-5s}{u}_1\left(s\right)\\ x_2\left(s\right)=\frac{2}{125000s^2+750s+1}{u}_2\left(s\right)\\ x_3\left(s\right)=\frac{1}{60s+1}{u}_3\left(s\right)\end{array}\right.$

$y\left(t\right)=\sqrt{2x_1\left(t\right)+3x_2\left(t\right)+x_3\left(t\right)-9}$

用Matlab仿真软件模拟训练数据样本：三个易测辅助变量为连续采样输入，输入类型为随机均匀分布，模拟采样周期为一分钟，产生15000个点的一组序列(模拟250小时的采样数据)，作为模型输入u₁(k),u₂(k),u₃(k)其对应着本方法中所述的 $u_k\left(\begin{array}{c}...\\ t\end{array}\right)(k=1,\mathrm{...},m);$ 根据上述三组组输入以及仿真对象的状态方程、输出方程，得到15000个点的过程子系统输出序列y(k)其对应着本方法中所述的y(t)，并将对应序列点外加10％的随机噪声v(k)，对该组输出序列进行稀疏采样，采样周期为P＝60分钟或70分钟，产生211个点的一组序列，作为模型输出y(P_i)，用于模拟离线采样数据。其中，149个数据样本用来建模，剩余62个样本用来对所建模型进行精度估计与泛化性能测试。误差衡量指标为均方误差衡量

$RMSE=\sqrt{\frac{1}{62}\underset{i=1}{\overset{62}{\&Sigma;}}{\&lsqb;o\left(P_i\right)-y\left(P_i\right)\&rsqb;}^2}.$

规定三组实验中，脉冲响应长度L＝40,60,80。三个输入子系统的纯时延寻优范围均为{0,1，…，8}。

其中，图2～图5均为本实施例所述的方法的计算结果与测量结果的比较图，从图示中可看出，采用本法进行质量变量的监测精度高的优点。在本实例中当L＝80时，本方法的训练精度最高。基于贝叶斯框架的动态软测量建模模型(采用贝叶斯方法建模)的测试结果：RMSE＝0.0721；基于支持向量机的静态软测量模型(采用贝叶斯方法建模)的测试结果：RMSE＝0.1660。测试结果表明该模型与传统的静态建模方法相比，具有精度高的优点。

下面结合说明书附图以及实施例对本发明基于贝叶斯框架的动态软测量建模装置做进一步的说明。

第一实施例：

本实施例基于贝叶斯框架的动态软测量建模装置至少包括模型建立单元、数据采集计算单元以及数据输出单元；

所述模型建立单元，用以一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(\mathrm{...},t)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(\mathrm{...},t)\\ .\\ .\\ .\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(\mathrm{...},t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)(k＝1,…,m)为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的行向量，x_k(t)(k＝1,…,m)为状态变量的第k个分量在采样时间t的取值,模型为一阶脉冲响应模板；且

其中，为的转置，α_k(k＝1,…,m)为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；ΔT为过程的采样周期；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\&sigma;}}^2}\}+b---\left(3\right)$

其中α_i,与b为支持向量机参数，σ为RBF核参数，为支持向量机的支持向量；

所述数据采集计算单元，用以采集在线测量值u_k(…,t)(k＝1,…,m)，并输入到一阶脉冲响应模型计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)；

所述数据输出单元，用以输出y(t)。

本实施例所述的基于贝叶斯框架的动态软测量建模装置通常为在线预报器。

首先，本实施例基于贝叶斯框架的动态软测量建模装置，通过将一阶脉冲响应模型与支持向量机模型的结合，以一段连续采样周期的采用组成的向量的为输入，从而打破了传统方法以及装置仅以单独的时间点来演算难测的主导变量即质量变量y(t)的局限，从而能够对过程的动态信息进行有效描述，利用了生产过程中的动态信息，与传统的方法以及装置相比估计精度更高，在实施闭环反馈控制时也具有更加优良的控制性能；

其次，本实施例所述的基于贝叶斯框架的动态软测量建模装置，充分利用支持向量机模型自身的对时延的动态调整以及向量的时间特性，从而能对输出时延、动态参数进行智能寻优，减少了现有算法对于过程动态先验信息的依赖，解决了现有部分动态软测量模型由于输出时延信息无法获取导致无法使用的问题；

再次，本实施例所述的基于贝叶斯框架的动态软测量建模装置，对一阶脉冲响应以及支持向量机模型的训练能够在不依赖传统测试数据集的基础上，根据客观准则选择模型训练中的各类参数，自动惩罚过拟合的模型，有效地改进由于非线性数据驱动模型过学习导致的模型泛化性能问题。

第二实施例：

如图6所示，基于贝叶斯框架的动态软测量建模装置，其特征在于，所述基于贝叶斯框架的动态软测量建模装置至少包括模型建立单元、数据采集计算单元以及数据输出单元；

所述模型建立单元，用以一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ .\\ .\\ .\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(\mathrm{...},t)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(\mathrm{...},t)\\ .\\ .\\ .\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(\mathrm{...},t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)(k＝1,…,m)为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的向量，x_k(t)(k＝1,…,m)为状态变量的第k个分量在采样时间t的取值,模型为一阶脉冲响应模板；且

其中，α_k(k＝1,…,m)为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；ΔT为采样周期；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\&sigma;}}^2}\}+b---\left(3\right)$

其中α_i,与b为支持向量机参数，σ为RBF核参数，为支持向量机的支持向量；

所述数据采集计算单元，用以采集在线测量值，并输入到一阶脉冲响应模型计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)；

所述数据输出单元，用以输出y(t)。

所述模型建立单元包括：

结构建立模块，用以构建如公式(1)以及公式(3)所述的关于u_k(…,t)(k＝1,…,m)以及y(t)的系统结构关系；

数据样本获取模块，用以获取与u_k(…,t)(k＝1,…,m)相对应的数据U_k(…,t)(k＝1,…,m)、以及与U_k(…,t)(k＝1,…,m)相映射的Y(t)构成的数据样本，并按采样时间形成递增的时间序列；

数据样本标准化模块，用以对数据样本中的数据进行标准化处理，以使数据在每一个维度上具有零均值和单位方差；

初始参数设定模块，用以设定一阶脉冲响应模型以及支持向量机模型训练的初始参数α_k,τ_k、脉冲响应序列长度L、训练参数ε、C、RBF核参数σ以及迭代停止阈值μ：

参数训练模块，用以将数据样本中的U_k(…,t)(k＝1,…,m)输入公式(1)，且将公式(1)的结果输入公式(3)，比较公式(3)的输出与Y(t)，并根据比较结果校正公式(1)以及公式(3)的参数。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (2)

1.一种基于贝叶斯框架的动态软测量建模方法，其特征在于，所述方法包括以下步骤：

步骤A：建立一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(...,t)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(...,t)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(...,t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)，k＝1,…,m为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的向量，x_k(t)，k＝1,…,m为状态变量的第k个分量在采样时间t的取值,模型k＝1,…,m为一阶脉冲响应模板；且

其中，α_k，k＝1,…,m为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；τ与τ_k意义相同，也是纯时延参数；ΔT为采样周期，τ_k/ΔT表示的0的个数；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\&sigma;}}^2}\}+b---\left(3\right)$

其中与b为支持向量机参数，σ为RBF核参数，为支持向量，l为支持向量的个数；

步骤B：测量u_k(…,t)，k＝1,…,m并输入一阶脉冲响应模型，计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)；

步骤C：输出y(t)；

所述步骤A还包括以下具体步骤：

步骤A1：构建公式(1)以及公式(3)的关于u_k(…,t)，k＝1,…,m以及y(t)的系统结构关系；

步骤A2：获取与u_k(…,t)，k＝1,…,m相对应的数据U_k(…,t)，k＝1,…,m、以及与U_k(…,t)，k＝1,…,m相映射的Y(t)构成的数据样本，并按采样时间形成递增的时间序列；

步骤A3：对数据样本中的数据进行标准化处理，以使数据在每一个维度上具有零均值和单位方差；

步骤A4：设定一阶脉冲响应模型以及支持向量机模型训练的初始参数α_k,τ_k均为0、脉冲响应序列长度L、训练参数ε、C、RBF核参数σ以及迭代停止阈值μ，初始化迭代次数ite＝0；

步骤A5：利用处理后的所述数据样本，基于脉冲响应序列计算状态变量x(t)，利用初始参数训练支持向量机；

步骤A6：令ite＝ite+1，根据如下式子更新最优的支持向量机惩罚参数C^MP；

$C^{MP}=\frac{\underset{i=1}{\overset{l}{\&Sigma;}}\underset{j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)({\mathrm{\&alpha;}}_j^{*}-{\mathrm{\&alpha;}}_j)\exp \{-\frac{\left|\right|x\left(P_i\right)-x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}}{\mathrm{\&gamma;}}$

而γ的计算方法如下：

$\mathrm{\&gamma;}=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&rho;}}_i}{1/C+{\mathrm{\&rho;}}_i}$

其中ρ_i为矩阵的特征值，矩阵为l×l维矩阵，第i行第j列的元素的计算如下：

$r_i\exp \{-\frac{\left|\right|x\left(P_i\right),x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}$

其中r_i＝r(y_i-f(x_i)-ε)+r(f(x_i)-y_i-ε)，r(u)＝{u[1-2s(u)]+2}·s(u)[1-s(u)]，s(u)为sigmoid函数1/(1+e^-u)；

步骤A7：利用更新的参数训练支持向量机；

步骤A8：根据下式计算最优的核参数σ^MP，并在之后利用新的参数C^MP再进行一次的训练：

${\mathrm{\&theta;}}^{MP}=\left\{\begin{array}{cc}\sqrt{\mathrm{\&eta;}}\mathrm{\&theta;},& 0<\mathrm{\&eta;}<1\\ \mathrm{\&eta;}\mathrm{\&theta;},& \mathrm{\&eta;}\&GreaterEqual;1\end{array}\right.$

θ＝1/σ，σ^MP＝1/θ^MP，

η的计算式如下:

$\mathrm{\&eta;}\left|\frac{\begin{array}{c}{\mathrm{\&lambda;}}^{MP}\underset{i,j=1}{\overset{l}{\mathrm{\&Sigma;}}}(a_i-a_j)(a_i^{*}-a_j^{*}){\left(x\right(P_i)-x(P_j\left)\right)}^2\\ \&times;\exp (-{\mathrm{\&theta;}}^2{\left(x\right(P_i)-x(P_j\left)\right)}^2/2)\end{array}}{\begin{array}{c}{\mathrm{\&lambda;}}^{MP}/(l-{\mathrm{\&lambda;}}^{MP}trace({\stackrel{\&OverBar;}{A}}^{-1}))\&times;trace\left({\stackrel{\&OverBar;}{A}}^{-2}\right(\&part;\stackrel{\&OverBar;}{K}/\&part;\mathrm{\&theta;})\\ +trace\left({\stackrel{\&OverBar;}{A}}^{-1}\right(\&part;\stackrel{\&OverBar;}{K}/\&part;\mathrm{\&theta;}))/\mathrm{\&theta;}\end{array}}\right|$

其中，

$\stackrel{\&OverBar;}{A}=\frac{1}{C}{I}_{l\&times;l}+\stackrel{\&OverBar;}{K}$

矩阵的定义见步骤A6；

步骤A9：根据如下方法分别求取m个子系统最优的动态参数 $\{{\mathrm{\&alpha;}}_k^{MP},{\mathrm{\&tau;}}_k^{MP}\},k=1,...,m:$

定义损失函数如下：

$F({\mathrm{\&alpha;}}_k,{\mathrm{\&tau;}}_k)=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}{\&lsqb;y\left(P_i\right)-f\left(x\right(P_i\left)\right)\&rsqb;}^2,0\&le;{\mathrm{\&alpha;}}_k<1,k=1,...,m$

$x\left(P_i\right)=\left[\begin{array}{c}{x}_1\left(P_i\right)\\ x_2\left(P_i\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ x_m\left(P_i\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(...,P_i)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(...,P_i)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(...,P_i)\end{array}\right]$

对于第j个子系统j＝1,…,m，其纯时延参数由均匀分布集合{0,1,…,d_j}确定，其中，对于每一个可能的纯时延利用黄金分割法在0≤α_k＜1范围内优化损失函数寻找最优的

而在得到后，最优的由下式确定：

${\mathrm{\&tau;}}_k^{MP}=\arg \underset{0\&le;{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k\&le;d_k}{\min }F\left({\mathrm{\&alpha;}}_k^{MP}\right({\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k),{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k)$

步骤A10：利用更新的参数训练支持向量机；

步骤A11：利用下式计算其evid(ite)：

$evid\left(ite\right)=-\frac{1}{C}\underset{i=1}{\overset{l}{\&Sigma;}}\underset{j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)({\mathrm{\&alpha;}}_j^{*}-{\mathrm{\&alpha;}}_j)\exp \{-\frac{\left|\right|x\left(P_i\right)-x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}-\underset{i=1}{\overset{l}{\&Sigma;}}(1+{\mathrm{C\&rho;}}_i)$

其中ρ_i为矩阵的特征值，矩阵的定义见步骤A6；

步骤A12：当evid(ite)-evid(ite-1)＜μ或者evid(ite)＜evid(ite-1)时停止迭代，此时的模型即为所求的动态软测量模型；否则，返回步骤A6继续迭代。

2.一种基于贝叶斯框架的动态软测量建模装置，其特征在于，所述基于贝叶斯框架的动态软测量建模装置至少包括模型建立单元、数据采集计算单元以及数据输出单元；

所述模型建立单元，用以一阶脉冲响应模型以及支持向量机模型；

所述的一阶脉冲响应模型为：

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ x_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(...,t)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(...,t)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(...,t)\end{array}\right]---\left(1\right)$

其中，x(t)为采样时刻t离散化的m维状态变量；u_k(…,t)，k＝1,…,m为输入变量的第k个分量在采样时间t之前L个采样周期内的每个周期一个采样值组成长度为L的向量，x_k(t)，k＝1,…,m为状态变量的第k个分量在采样时间t的取值,模型k＝1,…,m为一阶脉冲响应模板；且

其中，α_k，k＝1,…,m为一阶脉冲响应参数，其范围满足0≤α_k＜1；非负整数τ_k为一阶脉冲响应的纯时延参数；τ与τ_k意义相同，也是纯时延参数；ΔT为采样周期，τ_k/ΔT表示的0的个数；

所述的支持向量机模型为

$y\left(t\right)=\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)\exp \{-\frac{\left|\right|x\left(t\right)-x_i^{sv}|{|}^2}{2{\mathrm{\&sigma;}}^2}\}+b---\left(3\right)$

其中与b为支持向量机参数，σ为RBF核参数，为支持向量,l为支持向量的个数；

所述数据采集计算单元，用以采集在线测量值u_k(…,t)，k＝1,…,m并输入到一阶脉冲响应模型计算出状态变量x(t)并将x(t)输入到支持向量机模型中计算出y(t)；

所述数据输出单元，用以输出y(t)；

所述模型建立单元包括：

结构建立模块，用以构建公式(1)以及公式(3)的关于u_k(…,t)，k＝1,…,m以及y(t)的系统结构关系；

样本数据获取模块，用以获取与u_k(…,t)，k＝1,…,m相对应的数据U_k(…,t)，k＝1,…,m、以及与U_k(…,t)，k＝1,…,m相映射的Y(t)构成的数据样本，并按采样时间形成递增的时间序列；

数据样本标准化模块，用以对数据样本中的数据进行标准化处理，以使数据在每一个维度上具有零均值和单位方差；

初始参数设定模块，用以设定一阶脉冲响应模型以及支持向量机模型训练的初始参数α_k,τ_k均为0、脉冲响应序列长度L、RBF核参数σ以及迭代停止阈值μ，初始化迭代次数ite＝0；

参数训练模块，用以利用处理后的所述数据样本，基于脉冲响应序列计算状态变量x(t)，利用初始参数训练支持向量机；步骤如下：

步骤A6：令ite＝ite+1，根据如下式子更新最优的支持向量机惩罚参数C^MP；

$C^{MP}=\frac{\underset{i=1}{\overset{l}{\&Sigma;}}\underset{j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)({\mathrm{\&alpha;}}_j^{*}-{\mathrm{\&alpha;}}_j)\exp \{-\frac{\left|\right|x\left(P_i\right)-x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}}{\mathrm{\&gamma;}}$

而γ的计算方法如下：

$\mathrm{\&gamma;}=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&rho;}}_i}{1/C+{\mathrm{\&rho;}}_i}$

其中ρ_i为矩阵的特征值，矩阵为l×l维矩阵，第i行第j列的元素的计算如下：

$r_i\exp \{-\frac{\left|\right|x\left(P_i\right),x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}$

其中r_i＝r(y_i-f(x_i)-ε)+r(f(x_i)-y_i-ε)，r(u)＝{u[1-2s(u)]+2}·s(u)[1-s(u)]，s(u)为sigmoid函数1/(1+e^-u)；

步骤A7：利用更新的参数训练支持向量机；

步骤A8：根据下式计算最优的核参数σ^MP，并在之后利用新的参数C^MP再进行一次的训练：

${\mathrm{\&theta;}}^{MP}=\left\{\begin{array}{cc}\sqrt{\mathrm{\&eta;}}\mathrm{\&theta;},& 0<\mathrm{\&eta;}<1\\ \mathrm{\&eta;}\mathrm{\&theta;},& \mathrm{\&eta;}\&GreaterEqual;1\end{array}\right.$

θ＝1/σ，σ^MP＝1/θ^MP，

η的计算式如下:

$\mathrm{\&eta;}\left|\frac{\begin{array}{c}{\mathrm{\&lambda;}}^{MP}\underset{i,j=1}{\overset{l}{\mathrm{\&Sigma;}}}(a_i-a_j)(a_i^{*}-a_j^{*}){\left(x\right(P_i)-x(P_j\left)\right)}^2\\ \&times;\exp (-{\mathrm{\&theta;}}^2{\left(x\right(P_i)-x(P_j\left)\right)}^2/2)\end{array}}{\begin{array}{c}{\mathrm{\&lambda;}}^{MP}/(l-{\mathrm{\&lambda;}}^{MP}trace({\stackrel{\&OverBar;}{A}}^{-1}))\&times;trace\left({\stackrel{\&OverBar;}{A}}^{-2}\right(\&part;\stackrel{\&OverBar;}{K}/\&part;\mathrm{\&theta;})\\ +trace\left({\stackrel{\&OverBar;}{A}}^{-1}\right(\&part;\stackrel{\&OverBar;}{K}/\&part;\mathrm{\&theta;}))/\mathrm{\&theta;}\end{array}}\right|$

其中，

$\stackrel{\&OverBar;}{A}=\frac{1}{C}{I}_{l\&times;l}+\stackrel{\&OverBar;}{K}$

矩阵的定义见步骤A6；

步骤A9：根据如下方法分别求取m个子系统最优的动态参数 $\{{\mathrm{\&alpha;}}_k^{MP},{\mathrm{\&tau;}}_k^{MP}\},k=1,...,m:$

定义损失函数如下：

$F({\mathrm{\&alpha;}}_k,{\mathrm{\&tau;}}_k)=\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}{\&lsqb;y\left(P_i\right)-f\left(x\right(P_i\left)\right)\&rsqb;}^2,0\&le;{\mathrm{\&alpha;}}_k<1,k=1,...,m$

$x\left(P_i\right)=\left[\begin{array}{c}{x}_1\left(P_i\right)\\ x_2\left(P_i\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ x_m\left(P_i\right)\end{array}\right]=\left[\begin{array}{c}\stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_1,{\mathrm{\&tau;}}_1)}^T{u}_1(...,P_i)\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_2,{\mathrm{\&tau;}}_2)}^T{u}_2(...,P_i)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \stackrel{\&RightArrow;}{g}{({\mathrm{\&alpha;}}_m,{\mathrm{\&tau;}}_m)}^T{u}_m(...,P_i)\end{array}\right]$

对于第j个子系统j＝1,…,m其纯时延参数由均匀分布集合{0,1,…,d_j}确定，其中，对于每一个可能的纯时延利用黄金分割法在0≤α_k＜1范围内优化损失函数寻找最优的

而在得到后，最优的由下式确定：

${\mathrm{\&tau;}}_k^{MP}=\arg \underset{0\&le;{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k\&le;d_k}{\min }F\left({\mathrm{\&alpha;}}_k^{MP}\right({\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k),{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_k)$

步骤A10：利用更新的参数训练支持向量机；

步骤A11：利用下式计算其evid(ite)：

$evid\left(ite\right)=-\frac{1}{C}\underset{i=1}{\overset{l}{\&Sigma;}}\underset{j=1}{\overset{l}{\&Sigma;}}({\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&alpha;}}_i)({\mathrm{\&alpha;}}_j^{*}-{\mathrm{\&alpha;}}_j)\exp \{-\frac{\left|\right|x\left(P_i\right)-x\left(P_j\right)|{|}^2}{2{\mathrm{\&sigma;}}^2}\}-\underset{i=1}{\overset{l}{\&Sigma;}}(1+{\mathrm{C\&rho;}}_i)$

其中ρ_i为矩阵的特征值，矩阵的定义见步骤A6；

步骤A12：当evid(ite)-evid(ite-1)＜μ或者evid(ite)＜evid(ite-1)时停止迭代，此时的模型即为所求的动态软测量模型；否则，返回步骤A6继续迭代。