基于高斯混合聚类的切换系统的辨识

柴秀俊 ,王宏伟,,王 林 ,嵇薪儒

(1.新疆大学电气工程学院,新疆乌鲁木齐 830047;2.大连理工大学控制科学与工程学院,辽宁大连 110024)

1 引言

混杂系统是一类复杂系统,包含相互影响并相互作用的连续动力学和离散动力学[1].切换系统是一类典型的混杂系统,它在通讯、电力系统和机器人控制等领域中有广泛的应用[2–3].目前,在研究切换系统的结构特性上已经获得了许多重要的研究成果,比如在切换系统的控制器设计[4–6]与稳定性分析[7–9]等方面.而在切换系统的辨识领域,还存在很多问题需要解决.

切换系统一般由若干个子系统和决定它们之间切换的规则组成,这就导致了从切换系统收集的输入输出数据来自于不同的子系统.因此对于切换系统的辨识,通常需要通过模式检测计算每个采样数据的运行模式号以确定切换系统的切换时间与子系统的驻留时间,然后再对采样数据进行参数辨识.

对于切换系统的辨识方法,一般分为子空间方法和基于输入输出(input and output,IO)模型的方法.文献[10]首先引入子空间方法用于分段线性系统的辨识,文中假设切换规则已知,主要研究局部模式的辨识.对于未知切换规则的切换系统辨识,文献[11]在子空间框架中使用变化检测技术对切换系统进行辨识,但是这种方法需要假设系统只在稳态下进行切换.其它的子空间方法还包括投影子空间方法[12],结构化子空间方法[13]等.对于子空间方法,最大的问题就是要求两个连续的系统切换之间存在一个最小停顿时间[10].基于IO模型的方法是另一种简单而有效的切换系统辨识方法,比如自组织模型法[14]、变分贝叶斯法[15]、聚类方法[16–17]、随机迭代法[18]和椭球法[19]等.基于IO模型的方法适用于具有任意切换规则的切换系统,很好的解决了子空间方法中要求最小停顿时间问题.

上述文献提出的辨识方法中大都假设子系统的数量是先验已知的,而一般情况下具有未知规则的切换系统,只有输入数据与输出数据是可以采样的,子系统的数量、切换时间、子系统的驻留时间往往是未知的.在此基础上,本文提出一种将高斯混合聚类与递推增广最小二乘法结合的两阶段辨识方法对切换系统进行辨识.本文将具有未知切换规则与未知子系统数量的切换系统作为研究对象,在模式检测阶段,首先建立高斯混合模型表示采样数据的分布,并使用轮盘法选择较为合适的初始模型参数.其次,分别计算采样数据属于各个子系统的后验概率,同时通过极大似然估计算法迭代更新高斯混合模型中的模型参数,使高斯混合模型最大化地拟合所有采样数据的分布.在此基础上,通过贝叶斯信息准则(Bayesian information criterion,BIC)确定切换系统中子系统的数量,同时根据最大后验概率准则计算采样数据的运行模式号,从而估计切换系统的切换规则.在参数辨识阶段,通过递推增广最小二乘法估计每个子系统的参数向量.最后,利用含有有色噪声的切换系统来验证所提辨识方法的有效性.

2 问题描述

在存在有色噪声的情况下,切换系统可以表示为

其中:u(k)和y(k)分别代表第k次均匀采样的切换系统的输入和输出;v(k)表示均值为零、方差为δ2的高斯白噪声;w(k)为与时间相关的分段切换函数;N代表子系统的数量;z−1代表单位延迟算子,z−1y(k)=y(k−1);Ai(z−1),Bi(z−1)和Ci(z−1)分别表示对应的系统多项式,具体形式为

在k次采样时,切换系统被切换至第i个子系统,式(1)可以转换成以下线性回归模型:

其中:

其中θi和φi分别表示第i个子系统的参数向量和信息向量.

由于切换系统不提供任何的切换信息,所以必须在辨识子系统参数之前对切换系统进行模式检测,以确认切换信息[20].因此切换系统的辨识可以分为两个阶段:1)在模式检测阶段,通过检测技术估计切换规则,并通过目标函数法确定切换模式总数;2)在参数辨识阶段,通过辨识算法获得每个子系统的参数估计值.本文将在第3节中讨论模式检测问题,在第4节中讨论参数辨识问题.

3 基于高斯混合聚类的模式检测

聚类分析的算法可以分为划分法、层次法、基于密度的方法、基于网格的方法以及基于模型的方法[21].其中,基于模型的方法就是给每一个聚类假定一个模型,然后去寻找能够很好的满足这个模型的数据集.高斯混合聚类就是借助高斯混合模型来表示数据分布的一种基于统计模型的聚类方法.本文将通过高斯混合聚类对切换系统进行模式检测.

由于只有输入数据与输出数据是可以采样的,所以将输入数据与输出数据定义为IO向量ψ(k):

假设对切换系统进行采样共得到m组IO向量,构成IO向量 集ψ={ψ(1),ψ(2),···,ψ(m)},则ψ的 分布可由高斯混合模型表示:

其中:N代表子系统的数量;αi代表第i个子系统对应簇类的混合系数,满足条件代表切换系统在第k次采样得到的第i个子系统的IO向量;pζ(ψ(k))代表根据ψ的分布建立的高斯混合模型.由于系统的输入信号为具有零均值,单位方差的不相关随机序列信号,输出信号与输入信号呈线性关系,所以ψi(k)的分布由概率密度函数p(ψi(k)|µi,Σi)表示:

3.1 高斯混合模型初始模型参数的选择

由式(3)–(4)可以看出高斯混合模型由混合系数,均值向量与协方差矩阵所决定.对于高斯混合聚类,其聚类的效果十分依靠高斯混合模型中初始均值向量的选择,本文将通过轮盘法寻找一组合适的初始均值向量.

从IO向量集ψ中随机选取一组IO向量作为第一个初始聚类中心,分别计算其它的IO向量与当前聚类中心的IO向量之间的距离:

其中:ψe为当前聚类中心的IO向量,ψf为其它的IO向量,λ为距离系数.

在聚类算法中,为了保证聚类的准确性,通常选择相互距离较大的一组初始聚类中心,但在ψ中存在异常IO向量的情况下,如果直接通过式(5)选择相互距离最大的一组初始聚类中心,将导致收敛速度慢或者局部收敛.为了解决这个问题,本文通过计算IO向量的选择概率来选择初始聚类中心,IO向量之间的距离越大,越有可能被选为下一个聚类中心.在式(5)的基础上,IO向量的选择概率为

在每一轮聚类中心的选择中随机生成一个大于0且逼近于0的数τ,按顺序分别令选择概率p(ψf)减去τ,直到满足条件p(ψf)−τ>0时选择该组IO向量作为下一个初始聚类中心.当存在两个或两个以上的聚类中心时,式(5)中选取其它的IO向量与当前所有聚类中心的IO向量之间最大的距离.之后重复计算式(5)–(6),直到找到N个聚类中心,并将其构建为高斯混合模型中的初始均值向量.

3.2 运行模式号的计算

通过上式构建的高斯混合模型,用于计算IO向量集ψ中m组IO向量的运行模式号.假设ψj(k)是ψ中的一组IO向量,其中j=1,2,···,N代表需要计算的运行模式号.由贝叶斯定理可知,ψj(k)属于第i个子系统的后验概率为

为了使高斯混合模型最大化的拟合样本集ψ,需要迭代更新模型参数αi,µi和Σi.由式(3)可得样本集ψ的极大似然函数:

由于模型参数在每一轮迭代更新过程中,需要重复使用后验概率pζ(i|ψj(k)),为了减小算法复杂度,对pζ(i|ψj(k))进行0–1离散化处理,若

则令pζ(i|ψj(k))=1,表示将第k个IO向量ψj(k)归为第i个子系统;否则令pζ(i|ψj(k))=0.至此,将样本集ψ化分为N个数据集C1,C2,···,CN.

对于极大似然函数L(ψ),分别令

可以得到µi与Σi的更新公式:

其中:ψij(k)表示属于第i类数据集Ci的IO向量;mi表示数据集Ci中IO向量的总数,满足条件m1+m2+···+mN=m.

对于αi的更新公式,由上述推论可知

在迭代更新模型参数的过程中,当L(ψ)的增长小于给定阈值γ(γ>0)时,模型参数αi,µi与Σi停止更新.此时根据式(7),由最大后验概率原则计算出每个IO向量ψj(k)的运行模式号:

由式(13)计算可得j=i,表示IO向量ψj(k)属于第i个子系统.

3.3 切换模式总数的确定

在上述的计算过程中,切换模式总数或者子系统的数量N是未知的,本文通过比较模型拟合采样数据的优良性估计切换模式的总数.

赤池信息准则(Akaike information criterion,AIC)是衡量统计模型拟合优良性的一种标准,它提供了权衡估计模型复杂度和拟合数据优良性的标准[22].通常情况下,AIC 准则定义为

其中:LN(ψ)表示子系统数量分别为N(N=1,2,···,Q,Q为设置的最大子系统数量)时IO向量集ψ的极大似然函数,m代表IO向量的总数.

BIC准则是在AIC准则的基础上,增大了惩罚项,在样本数量过多时,可有效的防止模型精度过高造成的模型复杂度过高[22].假设子系统的数量分别为N=1,2,···,Q,根据贝叶斯信息准则,可以得到IO向量集ψ的BIC值:

一般而言,式(15)中模型复杂度mln(N)与似然函数项LN(ψ)会随着子系统数量N的增加而增加,当N较小时,似然函数项的增速大于模型复杂度的增速,从而导致BIC变小.当N过大时,模型过于复杂会出现过拟合的现象,似然函数项LN(ψ)增速减缓,导致BIC增大.所以从N=1,2,···,Q中选择使BIC值最小的N,故有

通过式(16)可以计算出子系统的数量N.

为了便于理解,上述的建模及计算过程,可以总结为图1.

图1 模式检测流程图Fig.1 Flow chart of mode detection

4 基于递推增广最小二乘法的参数辨识

通过第3节中的高斯混合聚类方法,可以计算出IO向量的运行模式号,在此之后,本文使用递推增广最小二乘法[23–24]获得式(2)中参数向量的估计值.由于信息向量φi(k)中的v(k)不可测,所以通过其估计值ˆv(k)代替,即

则递推增广最小二乘法的公式如下:

步骤1初始化参数PPPi(0),,k=1,其中i=1,2,···,N;

步骤2采集数据y(k)和u(k),通过式(21)构造信息向量

步骤3依次计算式(18)–(19)和式(20),获得参数向量θi的估计值,并通过式(22)计算ˆv(k);

步骤4判断k>m? 若是则完成参数辨识过程,若否则进入下一步骤,其中m为信息向量的总数;

步骤5k=k+1,返回步骤2进行新一轮的迭代.

5 仿真结果

为了验证所提方法的有效性,本文以如下包含3个子系统的切换系统作为研究对象:

其中:

在仿真中,将具有零均值,单位方差的不相关随机序列信号作为输入u(k);将均值为零,方差为δ2=0.1的高斯白噪声作为噪声信号v(k).设置距离系数λ=2,最大子系统数量Q=10,似然函数增长阈值γ=10−5,初始协方差矩阵Σi=III(其中i=1,···,N,···,Q,III为n阶单位矩阵),初始混合系数

本文采用周期性随机切换信号,每50个样本为一个周期,对1000组样本进行了仿真.图2展示了在一次独立实验中得到的BIC曲线,由图可以看出当N=3时,BIC为最小值,表明切换系统由3个子系统构成.图3为模式检测过程中IO向量对应于每个子系统的后验概率与实际切换过程的比较图,为了更加直观的表示,图中令所有小于10−9的后验概率都以10−9表示,对于IO向量,使其后验概率最大的子系统即为此IO向量所属的子系统.图4为实际的切换信号与通过模式检测重建的估计切换信号的对比图,图中两者的匹配率为93.7%.

图2 不同子系统数量对应的BIC值Fig.2 The BIC value corresponds to the number of different subsystems

图3 IO向量对应于每个子系统的后验概率与实际切换过程Fig.3 The IO vector corresponds to the posterior probabilities of each subsystems and the actual switched process

图4 实际切换信号及其估计Fig.4 Actual switched signal and its estimation

图5为通过50次重复实验统计,将本文所提的方法与投影子空间法[12]和聚类方法[16]进行对比得到的分类错误率箱式图,从图中可以看出本文所提方法的分类错误率主要集中在7%到11%之间,低于其它两种算法.表1为通过递推增广最小二乘法的参数辨识的统计结果,其中第1列和第2列的项分别是平均值和方差.可以看出,递推增广最小二乘法可以较为准确地估计切换系统的参数向量.

图5 不同方法分类错误率的箱式图Fig.5 Box plot of classification error rate of different algorithms

图6展示了在噪声信号的方差δ2分别为{0.1,0.2,···,1}时,切换系统的噪信比和模式检测的平均分类错误率.从图中可以看出分类错误率会随着噪信比的增大而增大,当噪声信号的方差为0.1时,分类错误率大致为8.75%,当噪声信号的方差增大到为1时,分类错误率增大到26.7%.结果表明该方法在噪信比较低的情况下表现良好,随着噪信比的增加,分类错误率也会略微增加.

图6 噪声信号方差不同时的分类错误率与噪信比的统计图Fig.6 Statistical chart of classification error rate and noise signal ratio when the variance of noise signal is different

从上面的仿真中可以得出以下结论:1)对于具有未知切换规则的切换系统的辨识,本文所提的方法不仅能够确定切换系统中子系统的数量,还能较为准确地计算出采样数据的模式运行号,获得整个切换规则运行情况;2) 在对切换系统进行模式检测之后,利用递推增广最小二乘法能够准确地得到各个子系统与噪声项的估计参数,也反应出模式检测过程的有效性.

6 结论

针对具有未知切换规则与未知子系统数量的切换系统,提出一种基于高斯混合聚类与递推增广最小二乘法的二阶段辨识方法.本文借助高斯混合模型表示采样数据的分布,并根据最大后验概率原则计算采样数据的运行模式号,以确定切换系统的切换时间与子系统的驻留时间.由于在高斯混合聚类的过程中初始均值向量对聚类效果影响较大,所以利用轮盘法选择一组较为合适的聚类中心作为高斯混合模型的初始均值向量.同时借助采样数据的极大似然函数,根据BIC准则确定切换子系统的数量.最后通过递推增广最小二乘法获得了较为准确的参数估计值.

切换系统的辨识在工业炼铁中硅含量的预测[14],配水网络的故障检测[25]等场景具有广泛应用.然而,对于切换系统的辨识依然还存在很多问题,比如在切换系统中含有非稳定子系统,但可以通过平均驻留时间理论来保证全局稳定的情况下的参数辨识,以及非线性切换系统的辨识问题与收敛性分析.这些问题依然值得思考,接下来我们将进行更深入的研究.