一种基于密度峰值聚类的多模型时空建模方法与流程

1.本发明涉及工业分布参数系统建模及数据预测的技术领域，尤其涉及一种基于密度峰值聚类的多模型时空建模方法。


背景技术：

2.工业领域中的热过程、对流扩散反应过程、化学气相沉积过程等都是非线性分布参数系统，其物理模型通常使用偏微分方程来表示，然而大部分情况下，其精确的偏微分方程往往因为系统的复杂性而难以确定，导致难以对其进行精确的建模与预测。基于数据驱动的时空模型可以在没有先验知识的情况下描述分布参数系统的时空动态，因此，在非线性分布参数系统建模中得到了广泛的应用。
3.在现有的数据驱动的工业热过程时空分离建模方法中，通常建立的时空模型只能针对存在较小的工作区域或者单一工况的情况，然而现实中许多复杂的分布参数系统具有大范围、多工况的特性，这就导致传统的单一全局时空模型无法精确的模拟这类系统的全过程，如何建立此类分布参数系统精确的时空模型有一定的挑战。


技术实现要素：

4.有鉴于此，本发明提供一种基于密度峰值聚类的多模型时空建模方法。本发明首先利用dpc聚类算法将时空变量划分为若干局部子空间，其每个子空间代表原始系统局部时空特征；其次，采用kl方法学习相应的局部空间基函数；然后，通过elm建立所有局部时间系数模型，并利用获得的局部空间基函数和相应的时间系数模型来重构局部时空模型；最后，通过lasso回归计算相应子模型的权重，对各局部时空模型进行加权求和得到的集成模型来逼近原系统。相比其他全局空间基函数构造的时空模型的方法，基于dpc的多模型时空分离方法通过多个局部简化的时空模型来逼近原始系统，因此，建立的分布参数系统时空模型精度也较高。该方法适用于具有大范围、多工况的非线性分布参数系统，具有建模简单、适用性强、预测精度高等优点。
5.为实现上述目的，本发明提出了一种基于密度峰值聚类的多模型时空建模方法，该方法包括如下步骤：
6.s101：利用离散控制系统收集分布参数系统热过程中的历史时空数据作为数据集；所述数据集包括系统输入变量和温度分布数据其中，n为样本数，m为输入变量数，n为输出变量数，r为实数集；
7.s102：对所述数据集进行预处理，各变量标准化处理后得到新的数据集均值为0，方差为1，预处理后的新的数据集输入变量为u∈rn×m，输出变量为时空变量数据y∈rn×n；
8.s103：采用dpc聚类算法将预处理后的时空变量数据划分为多个不同的子空间{y1(si,t),
…
,yk(si,t)}，其中，k为子空间个数，所述子空间代表原系统的局部特征；所述预处理后的时空变量数据即为的预处理后的温度分布数据y，定义为{y(si,t)|i＝1,...,l；t＝1,2,...,n；s∈ω}，其中，l为预处理后的时空变量数据在空间方向的数据点个数，s表示时
空数据的位置，si表示第i个时空数据点的空间位置，ω为坐标空间；
9.s104：利用kl方法对dpc聚类算法划分的各个子空间时空数据学习相应的局部空间基函数
10.s105：将划分的各个子空间时空变量数据投影到通过步骤s104学习到所述局部空间基函数上，从而得到对应的子空间低维时序数据
11.s106：使用elm方法构建各个低维空间的未知时序动态特性，即各个子空间输入变量uk(t)与低维时序数据之间的关系模型；
12.s107：根据所述关系模型预测各个子空间低维时序数据输出将所述各个子空间的低维时序数据输出与空间基函数进行时空合成得到的各个子空间局部时空预测输出
13.s108：利用lasso回归计算相应局部时空模型的权重，对各局部时空模型进行加权求和得到的集成模型来逼近原系统，即全局时空预测输出
14.s109：将全局的时空预测输出进行反归一化处理，得到温度预测值并对时空模型的性能进行评价。
15.进一步地，步骤s102中，所述标准化具体为：其中，μu、μy为原始数据的均值，σu、σy为原始数据的标准差。
16.进一步地，步骤s103中采用dpc聚类算法将预处理后的时空变量数据划分为多个不同的子空间，具体过程为：
17.s201：数据矩阵为{y(si,t)|i＝1,...,l；t＝1,2,...,n；s∈ω}，按初始化参数截断距离dc表示为式(1)：
18.dc＝sda[round(n
×
(n-1)
×
p)]
ꢀꢀꢀꢀꢀ
(1)
[0019]
式(1)中，n为样本数，round表示四舍五入，p为调节参数，在1％～2％之间取值；
[0020]
s202：计算数据矩阵中任意两个数据点之间的距离，得到距离矩阵；
[0021]
s203：根据截断距离dc，通过式(2)或者式(3)中任意一个计算式计算任意数据点的局部密度ρi：
[0022][0023][0024]
式(2)、(3)中，d
ij
＝dist(xi,xj)为数据点xi和xj之间的欧式距离；
[0025]
s204：通过式(4)计算任意数据点的距离δi：
[0026][0027]
s205：以ρi为横轴，以δi为纵轴，画出ρ-δ决策图；
[0028]
s206：利用ρ-δ决策图，将ρi和δi都相对较高的点标记为簇中心，将ρi相对较低但是δi相对较高的点标记为噪声点；
[0029]
s207：将剩余点进行分配，分配时将每个剩余点分配到它最近邻且密度比其大的数据点所在的簇；
[0030]
s208：返回多模型子空间{y1(si,t),
…
yk(si,t)}，k为子空间个数。进一步地，步骤s104中利用kl方法对dpc聚类算法划分的各个子空间时空数据学习相应的局部空间基函数，具体过程为：
[0031]
s301：对于每一个子空间yk(si,t),k＝1,
…
,k解耦成空间基函数和相应的时间系数的内积形式，如式(5)：
[0032][0033]
式(5)中，为无限维空间基函数，为无限维时间系数；
[0034]
s302：将式(5)截断简化为如式(6)：
[0035][0036]
式(6)中，是yk(s,t)的n阶近似项；
[0037]
s303：根据空间基函数的正交性质，采用kl方法设计最小化目标函数，所述最小化目标函数如式(7)：
[0038][0039]
式(7)中，||f(s,t)||＝(f(s,t),f(s,t))
1/2
，
[0040]
s304：求解式(7)，将其转换为寻求式(8)特征值问题：
[0041][0042]
式(8)中，式(4)中rk(s,ζ)＝《tk(s,t),tk(ζ,t)》为协方差函数；为拉格朗日乘子；ζ表示与s不同的另外任意一个点；
[0043]
s305：将表示为时空输出的线性组合如式(9)：
[0044][0045]
s306：将式(9)代入式(8)，特征值问题式(8)转化为如式(10)：
[0046][0047]
式(10)中，为第i个特征向量，
[0048]
s307：计算出式(10)的特征向量，则子空间的空间基函数式(9)得到；
[0049]
由于c是对称半定矩阵，所得到的特征函数也是正交的，将特征值进行排序，使其满足阶数n根据以下选择准则进行计算：
[0050][0051]
式(11)中，参数n的选取在比例η超过99.9％时，用式(11)确定，代表n阶时空能够近似原始系统；
[0052]
s308：通过以上式(7)～式(11)得到子空间的局部空间基函数进一步地，步骤s105中，空间基函数是单位正交的满足式(12)：
[0053][0054]
式(12)中，是与的内积；
[0055]
得到对应的子空间低维时序数据具体计算如式(13)：
[0056][0057]
进一步地，步骤s106中，使用elm方法构建各个低维空间的未知时序动态特性，即各个子空间输入变量uk(t)与低维时序数据之间的关系模型，具体步骤为：
[0058]
s401：elm由输入层、隐藏层和输出层组成，输入层与隐藏层、隐藏层与输出层间由神经元连接；其中，l个隐藏层节点的单隐藏层神经网络可以表示为如下形式：
[0059][0060]
式(14)中，表示n个任意的样本其中的一个，表示n个任意的样本其中的一个，g(
·
)为激活函数，为输入权重，为输出权重，是第i个隐藏层单元的偏置，表示为w
ik
和的内积；
[0061]
s402：elm以零误差来逼近任意一个连续非线性的目标函数，即存在w
ik
、使得下式(15)成立：
[0062][0063]
s403：将式(15)中的n个方程转化为矩阵乘积形式如式(16)：
[0064]hk
·
βk＝tkꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0065]
式(16)中，hk是隐藏层节点的输出，βk为输出权重，tk为期望输出，其表示如式(17)：
[0066][0067]
s404：计算得到使得：
[0068][0069]
式(18)中，i＝1,
…
,l，等价于最小化损失函数：
[0070][0071]
式(19)的解通过式(20)求逆确定，即：
[0072][0073]
式(20)中，为矩阵hk的moore-penrose广义逆。
[0074]
进一步地，步骤s107中，将所述各个子空间的低维时序数据输出与空间基函数进行时空合成得到的各个子空间局部时空预测输出，具体如式(21)：
[0075][0076]
进一步地，步骤s108中，利用lasso回归计算相应子模型的权重，对各局部时空模型进行加权求和得到的集成模型来逼近原系统，得到全局的时空预测输出，具体步骤为：
[0077]
s501：各个子空间局部时空模型集成用数学形式描述如式(22)：
[0078]
[0079]
式(22)中，k为局部模型的个数，为第k个局部模型，wk＝[w1,w2,
…
,wk]
t
为第k个局部模型的权重；
[0080]
s502：采用l1正则化的线性回归方法——lasso算法来估计权值wk，具体计算如式(23)：
[0081][0082]
s503：全局时空预测模型输出通过求和形式进行组合，具体计算如式(24)：
[0083][0084]
式(24)中，为全局的时空预测输出。步骤s109中，模型性能评价指标有3个，包括：绝对误差ae、时间归一化绝对误差tnae、空间归一化绝对误差snae。
[0085]
本发明提供的有益效果是：适用于具有大范围、多工况的非线性分布参数系统，具有建模简单、适用性强、预测精度高等优点。
附图说明
[0086]
图1是本发明一种基于密度峰值聚类的多模型时空建模方法的结构示意图；
[0087]
图2是本发明一种基于密度峰值聚类的多模型时空建模方法的流程示意图；
[0088]
图3是本发明实施例中催化反应棒的过程简图；
[0089]
图4是本发明实施例采用的随机输入信号示意图；
[0090]
图5是本发明实施例采集到的输出数据示意图；
[0091]
图6是本发明实施例得到的ρ-δ决策图；
[0092]
图7是本发明实施例划分子空间示意图；
[0093]
图8是本发明预测输出示意图；
[0094]
图9是本发明绝对误差ae分布示意图；
[0095]
图10是各个模型每个时刻的总空间误差信息的snae对比情况；
[0096]
图11是各个模型每个传感器的总时间误差信息的tnae的对比情况。
具体实施方式
[0097]
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。
[0098]
本发明给出基于密度峰值聚类的多模型时空建模方法，基本思路请参考图1。关于图1解释如下：
[0099]
(1)在整个工业热过程中，输入量u和输出量y之间的精确的偏微分方程难以得到，本技术通过一种数据驱动的间接方式进行建立；
[0100]
(2)采用dpc聚类算法将预处理后的时空变量数据划分为多个不同的子空间{y1(si,t),
…
,yk(si,t)}；
[0101]
(3)将各个子空间的输出量y初步表示成无限空间基函数与无限时间系数的乘积
和形式
[0102]
(4)利用kl方法将各个子空间的无限空间基函数与无限时间系数降维至有限维；
[0103]
(5)采用elm方法建立各个子空间的输入量uk与有限维时间系数之间的关系模型；
[0104]
(6)通过elm建立的关系模型与空间基函数进行时空合成，可以得到各个子空间局部时空预测输出；
[0105]
(7)利用lasso回归计算相应子模型的权重，对各局部时空模型进行加权求和得到的集成模型来逼近原系统，即全局时空预测输出。
[0106]
该方法针对分布参数系统热过程的建模问题，首先采用dpc聚类算法将预处理后的时空变量数据划分为多个不同的子空间；然后，采用kl方法对时空数据进行时空分离，获得空间基函数和时间系数；之后，将elm算法用来近似低维空间的未知时序动态特性，即系统输入变量与低维时序数据之间的关系模型，并对各个子空间进行时空合成；最后，利用lasso回归计算相应子模型的权重，对各局部时空模型进行加权求和得到的集成模型来逼近原系统，并对建立预测模型的性能进行评价。
[0107]
本发明采用的技术方案的流程如图2所示，主要步骤如下：
[0108]
s101：利用离散控制系统收集分布参数系统热过程中的历史数据，主要包含统输入变量以及温度分布数据其中，n为样本数，m为输入变量数，n为输出变量数，r为实数集。
[0109]
s102：对数据集进行预处理，采用进行标准化，其中，μu、μy为原始数据的均值，σu、σy为原始数据的标准差。各变量标准化处理后得到新的数据集均值为0，方差为1，输入变量为u∈rn×m，输出变量为y∈rn×n。
[0110]
s103：采用dpc聚类算法将预处理后的时空变量数据划分为多个不同的子空间{y1(si,t),
…
,yk(si,t)}，其中，k为子空间个数，数据矩阵为{y(si,t)|i＝1,...,l；t＝1,2,...,n；s∈ω}，按初始化参数截断距离dc表示为式(1)：
[0111]
dc＝sda[round(n
×
(n-1)
×
p)]
ꢀꢀꢀꢀꢀꢀꢀ
(1)
[0112]
式(1)中，n为样本数，round表示四舍五入，p为调节参数，大小在1％～2％之间。
[0113]
计算任意两个数据点之间的距离，得到距离矩阵。
[0114]
根据截断距离dc，通过式(2)或者式(3)计算任意数据点的局部密度ρi：
[0115][0116][0117]
式(2)(3)中，d
ij
＝dist(xi,xj)为数据点xi和xj之间的欧式距离。
[0118]
通过式(4)计算任意数据点的距离δi：
[0119][0120]
以ρi为横轴，以δi为纵轴，画出ρ-δ决策图。
[0121]
利用ρ-δ决策图，将ρi和δi都相对较高的点标记为簇中心，将ρi相对较低但是δi相对较高的点标记为噪声点。
[0122]
将剩余点进行分配，分配时将每个剩余点分配到它最近邻且密度其大的数据点所在的簇。
[0123]
返回多模型子空间{y1(si,t),
…
yk(si,t)}，k为子空间个数。
[0124]
s104：利用kl方法对dpc聚类算法划分的各个子空间时空数据学习相应的局部空间基函数对于每一个子空间yk(si,t),k＝1,
…
,k可以解耦成空间基函数和相应的时间系数的内积形式，如式(5)：
[0125][0126]
式(5)中，为无限维空间基函数，为无限维时间系数。
[0127]
在实际应用中，式(5)可截断简化为如式(6)：
[0128][0129]
式(6)中，是yk(s,t)的n阶近似。
[0130]
根据空间基函数的正交性质，采用kl方法设计最小化目标函数，所述最小化目标函数如式(7)：
[0131][0132]
式(7)中，||f(s,t)||＝(f(s,t),f(s,t))
1/2
，
[0133]
求解式(7)可以转换为寻求式(8)特征值问题：
[0134][0135]
式(8)中，式(4)中rk(s,ζ)＝《tk(s,t),tk(ζ,t)》为协方差函数；为拉格朗日乘子；ζ表示与s不同的另外任意一个点。
[0136]
假设可以表示为时空输出的线性组合如式(9)：
[0137]
[0138]
进一步将式(9)代入式(8)，特征值问题式(8)可以转化为如式(10)：
[0139][0140]
式(10)中，为第i个特征向量，
[0141]
计算出特征向量，则子空间的空间基函数可以由式(9)得到。由于c是对称半定矩阵，所得到的特征函数也是正交的，假设特征值满足那么阶数n可根据以下选择准则进行计算：
[0142][0143]
式(11)中，参数n通常选取比例η超过99.9％时用式(11)确定，代表n阶时空能够近似原始系统。
[0144]
通过以上式(7)～式(11)便可以得到子空间的局部空间基函数
[0145]
s105：各个子空间的低阶时序数据可以通过投影到第四步得到空间基函数上来获得，如下：
[0146][0147]
式中，尖括号代表内积。
[0148]
s106：使用elm方法构建各个低维空间的未知时序动态特性，即各个子空间输入变量uk(t)与低维时序数据之间的关系模型。elm由输入层、隐藏层和输出层组成，输入层与隐藏层、隐藏层与输出层间由神经元连接。假设有n个任意的样本其中则对于一个有l个隐藏层节点的单隐藏层神经网络可以表示为如下形式：
[0149][0150]
式(14)中，g(
·
)为激活函数，为输入权重，为输出权重，是第i个隐藏层单元的偏置，表示为w
ik
和的内积。
[0151]
elm以零误差来逼近任意一个连续非线性的目标函数，即也就是存在w
ik
、使得下式(15)成立：
[0152][0153]
将式(15)中的n个方程可以转化为矩阵乘积形式如式(16)：
[0154]hk
·
βk＝tkꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0155]
式(16)中，hk是隐藏层节点的输出，βk为输出权重，tk为期望输出，其可以表示如式(17)：
[0156][0157]
为了能够训练单隐藏层神经网络，希望得到使得：
[0158][0159]
式(18)中，i＝1,
…
,l，这等价于最小化损失函数：
[0160][0161]
在elm的参数辨识过程中，输入权重w
ik
和隐藏层偏置都是随机产生的，与训练数据无关，一旦被随机确定，则输出矩阵hk就被唯一确定，只有输出权重βi是未知的。它的解可以通过式(20)求逆确定，即：
[0162][0163]
式(20)中，为矩阵hk的moore-penrose广义逆。
[0164]
s107：根据所述关系模型预测各个子空间低维时序数据输出将所述各个子空间的低维时序数据输出与空间基函数进行时空合成得到的各个子空间局部时空预测输出具体如式(21)：
[0165][0166]
s108：利用lasso回归计算相应子模型的权重，对各局部时空模型进行加权求和得到的集成模型来逼近原系统，即全局时空预测输出各个子空间局部时空模型集成用数学形式描述如式(22)：
[0167][0168]
式(22)中，k为局部模型的个数，为第k个局部模型，wk＝[w1,w2,
…
,wk]
t
为第k个局部模型的权重。
[0169]
采用l1正则化的线性回归方法——lasso算法来估计权值wk，该方法会使得部分学习到的特征权值为0，从而达到稀疏化和特征选择的目的，可以很好的解决每个局部时空模型之间可能相关性很大的问题，以得到一个稳健的权重，具体计算如式(23)：
[0170][0171]
全局时空预测模型输出可以通过求和形式进行组合，具体计算如式(24)：
[0172][0173]
式(24)中，为全局的时空预测输出。
[0174]
s109：将预测输出值进行反归一化处理，得到温度预测值并对模型的性能进行评价。
[0175]
以下结合一个具体的催化反应棒仿真过程来说明本发明的有效性。化反应棒是化学工业中的一个运输-扩散反应过程，它属于典型的分布参数系统热过程，催化反应棒反应过程简图如图3所示。假设催化反应棒的密度、比热容、热传导系数以及反应棒两端的温度都是恒定不变的，则它的偏微分方程的表示式如下：
[0176][0177]
式中，t(x,t)为反应棒的时空温度分布，βr为反应热量，γ为执行器数量，βu为热传递系统，u(t)为输入信号，b(t)
t
为输入信号分布。
[0178]
满足狄利克雷边界条件和初始条件：
[0179][0180]
过程参数通常设置如下：
[0181][0182]
本仿真实验中采用随机输入信号，其取值范围为[-3,5]，如图4所示。采用了20个传感器来采集系统输出信号，这20个传感器均匀分布于催化反应棒的水平方向上。采样间隔δt＝0.01s，采样时间t＝40s。采用有限差分方法，共采集到4000组仿真实验数据。采集到的输出数据t(s,t)，如图5所示，其中第6个和第11个传感器用于模型验证使用，而剩余的18个传感器用于建模训练使用。采集到的前3500组仿真数据用于训练，后500组仿真数据用
于验证效果。
[0183]
数据预处理后，首先利用dpc聚类算法将时空变量划分为若干局部子空间，其每个子空间代表原始系统局部时空特征，在本实例中通过dpc算法在平面坐标中画出ρ-δ决策图，如图6所示，从图中可以看出，第44号、第73号、第1463号这三个数据点由于同时具有较大的ρ值和δ值，是符合dpc确定聚类中心原则的，所以自动被确认为聚类中心，因此原始空间被划分为三个子空间，划分结果如图7所示。然后，采用kl方法学习相应的局部空间基函数，本实例中前三个主导基函数即可获得超过99％的系统主要动态特征，因此，共获取三个子空间且每个子空间有三个空间基函数。最后，在分布式参数系统热过程的降维子空间中，采用elm对降维子空间中时间系数进行建模，之后将空间基函数与elm时间数模型相结合完成局部时空模型的合成，使用lasso方法对局部时空模型进行局部加权求和得到最终集成模型，便可以得到预测的未来温度分布情况。
[0184]
为了测试模型的性能，500个测试输入信号被用来激励得到的集成多模型，本发明方法的预测输出如图8所示，绝对误差ae分布如图9所示，从图中可以看出，所建立的模型对测试数据的预测输出与实际输出基本相同，绝对误差很小，说明该多模型建模方法在经典的催化棒反应过程中具有良好的精度。
[0185]
为进一步验证模型性能，在相同的实验条件下，将其与两种其它的建模方法进行了比较。一种是传统的基于kl方法的单模型时空模型；另一种采用传统的k-means聚类算法将时空变量数据划分为不同的子空间，并对每个子空间构造基于kl方法的局部时空模型，由于k-means聚类算法需要提前指定聚类中心个数，此处设置的聚类中心个数与dpc确定的聚类中心个数一致，即聚类中心个数k＝3，另外k-means多模型集成的时候采用的是pcr回归方法。图10和图11分别展示了表示每个时刻的总空间误差信息的snae与表示每个传感器的总时间误差信息的tnae的对比情况。从图中可以看出，总体来说与其它两种方法相比，所提方法不管是snae还是tnae都更小，说明所提方法误差更小，建模精度得到了明显提升。
[0186]
本发明实施的有益效果是：该方法适用于具有大范围、多工况的非线性分布参数系统，具有建模简单、适用性强、预测精度高等优点。
[0187]
在不冲突的情况下，本发明中上述实施例及实施例中的特征可以相互结合。
[0188]
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。