基于巴氏系数的非高斯随机系统性能评估

黄国豆,张金芳,黄俊雄,翟宇卓

(华北电力大学控制与计算机工程学院,北京 102206)

1 引言

控制回路在自动化系统中起着最重要的作用。性能评估的意义在于实现、恢复和保持控制回路的最佳性能。因此,在工业运行中快速、准确地评估系统的性能具有重要的应用价值[1]。

实际工业过程中的噪声扰动大都不符合高斯分布,低阶统计量(均方误差)所包含的信息显然无法完全反映变量的高阶统计特性。对于非高斯随机系统,应该在均值和方差之外选择一个更合适的指标来反映变量的高阶统计特性[2]。信息熵作为一种不确定性度量,相较于均值或方差,对于任何随机变量都有更为一般的含义。所有高阶矩(不仅仅是二阶)的最小化都可以通过最小化熵来实现[3,4]。因此,近年来一种基于熵的控制系统性能评估方法发展迅速。

钟振芳[5]等通过计算最小熵控制下的概率密度函数(Probability Density Function,PDF),以此为基准构建了系统的性能评估指标。Zhang[6]使用最小Renyi熵作为评估指标,将闭环控制下的输出分解为只与扰动相关的反馈不变量和反馈变量。Zhou[4]提出了基于有理熵的系统性能评价指标,并给出了该指标理论基准值和估计基准值的计算方法。在文献[7,8]中,熵指标被用于评估典型串级控制系统的性能。然而,熵指标具有均值漂移的缺点,且熵基准计算过程十分复杂,所耗时间较长,因此缺乏实际工程意义。

针对以上问题,本文提出了一种基于巴氏系数的非高斯随机系统控制回路性能评估指标。在统计学中,巴氏距离(Bhattacharyya distance)多被用来测量两个离散或连续概率分布之间的相似性,常用于性能监测[9],图像处理[10],神经网络[11]等领域。由于该指标具有无界性,因此并不适合直接用作系统性能评价指标,但与其密切相关的巴氏系数却是对称且有界的,在反映两个分布之间的相似性方面并无太大差异,且计算更为简便。通过计算系统基准输出与实际输出之间的巴氏系数就可以得到系统的性能指标。

为了更快速准确地估计性能基准,结合摸石头过河算法(Wading Across Stream Algorithm,WSA)的优点,提出了一种改进的混合分布估计算法(Hybrid Estimation of Distribution Algorithm,HEDA)进行系统参数寻优和噪声PDF估计。通过对高斯和非高斯噪声扰动下的单回路反馈控制系统的仿真,验证了改进算法和新指标的有效性。

2 性能评估指标

2.1 巴氏系数

巴氏系数是衡量两个统计样本之间重叠程度的度量,与巴氏距离密切相关。同时,它还可以用来确定两个样本是否相近及可分离性[11]。

对于定义在同一域X上的概率分布P和Q,巴氏距离的定义为

DB(P,Q)=-ln(BC(P,Q))

(1)

其中,BC(P,Q)即为概率分布P和Q的巴氏系数。对于离散的概率分布,P和Q的巴氏系数定义如下

(2)

且满足条件0≤BC≤1,0≤D≤∞。

对于连续概率分布,P和Q的巴氏系数定义为

(3)

两个概率分布之间的巴氏系数值会随着两个样本相同部分的增多而变大。对于任意两个概率分布,巴氏系数越接近1,意味着两个概率分布越相似;反之,表示两个概率分布越不相似。

2.2 巴氏系数性能指标

为便于研究,考虑下图1所示的单输入单输出反馈控制回路。其中y(k),u(k)和v(k)分别代表过程的输出、输入和未知随机噪声扰动。

图1 SISO反馈控制回路结构框图

图2 点分布离散近似

假定通过寻优算法获取的基准数据yb为系统性能理想时的数据,实际输出的yt是待评估数据,那么基于巴氏系数的性能指标IBC

(4)

式中,p(x)和q(x)分别为yb与yt的概率分布。基于巴氏系数的性能指标具有如下性质:

1) 有界性:0≤IBC≤1。当基准与实际输出概率分布越相近时,性能指标越接近于1,系统越接近于理想性能;反之,性能指标越接近于0,系统性能越差,需要改进。

2) 对称性:IBC(p,q)=IBC(q,p)。

2.3 IBC性能指标计算

如果p,q同时满足高斯分布,p～N(μp,σp),q～N(μq,σq),那么其巴氏距离可以通过提取它们的均值和方差进行简便计算

根据式(1)可推出巴氏距离与巴氏系数的关系如下

BC(p,q)=e-DB(p,q)

(6)

此时的性能指标IBC即可由式(5)和(6)求得

IBC=e-D(p,q)

(7)

对于实际系统中更为常见的非高斯输出分布,仅依靠均值和方差显然不能完全描述该随机变量的统计特性。此时,将基准输出yb和实际输出yt进行离散化处理。

首先对随机变量进行固定频率采样,记录收集到的有效数据并确定其落入的区间,将区间均匀且不重叠地划分为多个子区间,统计落入每个子区间的数据个数,将频率记为每个子区间中心值的概率,即用不同子区间的中心值概率来近似逼近连续随机变量的概率分布。

对于概率p的计算,直接利用样本的直方图对其进行离散化处理,具体的实现步骤如下:

1) 以固定频率采样x,收集N个有效数据点,将所有数据都放在区间S=[E1,EM]内;

2) 将区间S平均分为M个子空间(S1∪S2∪…∪SM),每个子空间互不交叉,并将每个子空间的中心点记ei(i=1,2,3,…,M);

3) 统计每个子区间内采样数据的个数l作为频数,记为(l1,l2,…,lM);

4) 将频数与数据总量N(N≫M)相比得到频率,记为P：

(8)

BC=∑Mi=1i(x)i(x)

(9)

3 基准输出数据的获取

3.1 基准输出的定义

假定图1系统中的设定值为0,则系统输出为

(10)

Gv(z-1)=F(z-1)+z-τR(z-1)=(1+n1z-1+n2z-2+…+nτ-1z-(τ-1))+z-τR(z-1)

(11)

式中,F(z-1)是扰动传递函数Gv的脉冲响应系数,R(z-1)是满足恒等式(11)的剩余函数。

y(t)=Fv(t)+Lv(t-τ)=

(12)

对于线性非高斯随机系统,最小熵控制器的目标就是最小化输出变量的熵[12,13]。

H(yt)=H(Fvt+Lvt-τ)≥H(Fvt)

(13)

由式(13)可知,当且仅当L=0时,系统输出数据的熵达到最小值,此时控制系统的性能达到最优,相应的反馈不变量部分即为基准输出数据,记:

yb(t)=Fv(t)=(1+n1z-1+n2z-2+…+nvz-(τ-1))v(t)

(14)

实际输出数据yt可由Simulink仿真平台采集得到。

3.2 系统参数辨识和噪声PDF估计

实际工业过程中的控制系统相当于一个“黑箱”,内部的结构参数均是未知的,而根据式(14),计算基准输出数据必须知道模型参数[n1,n2,…,nv],因此就需要对给定系统进行参数辨识。

图1所示的随机线性系统离散时间模型可描述为

A(z-1)y(k)=z-τB(z-1)u(k)+C(z-1)v(k)

(15)

其中A,B,C可表示为

(16)

其中ana,bnb,cnc和τ是模型的结构参数;na,nb,nc分别是A(z-1),B(z-1)和C(z-1)的阶次,τ为系统延迟。

非高斯随机系统参数辨识问题本质上是高维参数空间优化的问题,因此系统参数可以通过优化算法获得。为了获得系统参数,首先要知道给定系统的阶次和时延。

对于迟延τ大小的估计,最简单常用方法是通过分析输入ut和输出yt信号之间的相关性,即

(17)

应用Akaike信息准则[14]得到系统的阶次。对于给定的CARMA模型,AIC准则为

(18)

3.2.1 改进的混合EDA算法

EDA算法是一种随机搜索启发式算法,用于创建解空间的概率模型,该模型基于模型采样的解的质量进行迭代更新[15]。在之前的研究中,EDA算法已经暴露出在优化过程中存在计算量大、局部搜索能力较差等缺点。本文针对以上缺点进行了改进。

1) 获取参数辨识空间

在参数未知的情况下,参数空间过大会影响算法收敛速度,特别是对于高维空间无法保证均匀采样的情况。因此,在初始化参数空间之前,需要对模型参数进行初步估计,以确保参数空间尽可能覆盖真实参数。对于非高斯系统,没有很好的参数初始化方法,因此采用高斯系统参数初始化方法进行粗略估计。

首先将CARMA模型写成最小二乘形式

y(k)=φT(k)θ+v(k)

(19)

采用递推增广最小二乘算法(RELS)对系统参数进行初步估计。目标是求目标函数J()达到极小值时的参数,J()定义如下

J()＾U3]2

(20)

2) WSA算法思想

为了提高算法的局部搜索能力,在传统EDA算法的中引入了WSA算法。WSA算法[16]是受“摸石头过河”思想的启发,首先在“岸边”慎重考察一下选择初始起点,向该起点附近邻域随机搜索若干解,找出其中最优解作为迭代结果;然后在以这个点为起点再向附近邻域随机搜索若干个解,选出最优解作为第三次迭代结果,以此类推直到满足终止条件[12]。

a)起点

在参数辨识空间已知的情况下,初始选择的解将会对整个算法产生较大的影响。找初始解的思路就是先在参数空间内均匀随机取值产生R个解,找出其中的最优解作为起点解。对于(20)中的待辨识参数向量要在区间[-3v,+3v]内均匀随机取值,产生R个解A=[1;2; …;R],并计算这R个解的目标函数值(误差熵值),根据目标函数值高低找出最优解*记为初始起点解,为方便描述,将起点解记*=[x*1,x*2,…x*n],n=na+nb+nc+1。

b)邻域搜索策略

根据“摸石头过河”的思想,当摸到一个“石头”后必然要以该“石头”为起点向周围摸索其它“石头”,以此类推进行搜索。因此当得到起点解后就在起点解的邻域半径Lk内随机搜索m个领域解组成B。

B=[′1;′2;…;′m]

(21)

3)交叉操作

传统寻优算法得到的最优解只包含一个解的信息,这在一定程度上忽略了其它优秀解中所包含的信息,改进方法就是交叉操作。根据适应度值对以上m个个体进行排序,选出前N*(N*≪m)个优秀解(适应度值最优的个体),D=[′1;′2;…;′N*]。计算D的中心点式(22)所示

(22)

*=a*+(1-a)θ

(23)

式中:a为[0,1]的随机数。

4)适应度值的选择

对于本文的输入输出模型,由式(19)可知,k时刻的误差可以定义为:

(24)

对于估计的误差序列e=[e1,e1,…,eL],使用二阶Renyi熵来估计过程参数,表示为:

(25)

通过最小化熵的误差序列,能够得到最优模型参数和输出与扰动之间关系。当误差熵达到最小值时对应的参数也达到最优。即:

opt=arg minθ∈ΩWH2(e)

(26)

式中,ΩW即为参数寻优空间。

5)算法总结

HEDA优化算法的程序步骤如下:

算法1寻优算法的程序步骤:

①用输入输出模型描述系统,通过分析u(t)和y(t)之间的相关性来估计系统延迟τ;确定模型阶次na,nb,nc;

②初步估计。由RELS算法得到初步的模型参数,并以此确定寻优空间ΩW为±3v;

③慎重选择初始起点。在参数空间ΩW内随机生成R个个体A=[1;2; …;R],计算每个个体的适应度值(误差熵),根据适应度值(对应的误差熵最小)选择最优解*0;

WhileI≤nmax

a)计算上述m个解的适应度值并按大小进行排序CI=[(I)*1;(I)*2;…;(I)*m];

b)从CI中筛选出最优个体*I和前N*个适应度值最优的个体DI=[(I)*1;(I)*2;…;N*(I)*];

c)修改最优个体*I。计算DI的中心点(均值)交叉操作并修改最优个体的位置为*I=a*I+(1-a)θI,a=rand;

d)If 满足终止条件

两次相邻迭代误差熵之差小于0.0001,结束循环;

Else

LI=αLI-1;

I=I+1;

End

End

3.2.2 基于参数灵敏度分析的算法验证

考虑如下ARMA系统,通过分析改进算法初始参数的敏感性,进一步比较HEDA和传统EDA。

(27)

上述仿真实例的待寻优参数向量可表示为θ=[-1.7,0.7,1.5,0.9]。大多数的算法初始参数是基于前人的工作[2-4,12-14]。一般情况下,EDA中初始种群的个体数为N=1000,每次迭代的最佳个体数m=200,最大迭代次数是nmax=120。在WSA中,交叉操作的优秀个体数为N*=30,R=800。本文重点讨论收缩系数α和初始邻域半径L0的初始值的敏感性。

表1 基于初始参数敏感性分析的模型参数辨识结果

综上可得,获取系统基准输出数据主要分为如下两步:

仿真验证

在Matlab/Simulink中搭建图1所示的典型SISO系统,考虑如下CARMA模型:

(28)

控制器的传递函数为:

(29)

表2 不同噪声扰动下的参数辨识结果

图3 实际和估计的高斯噪声PDF

图4 实际和估计的非高斯噪声PDF

图3,4可以看出,无论对于高斯噪声还是非高斯噪声,HEDA算法估计出的扰动分布与真实的系统扰动都更加接近;表2中的测试结果也进一步验证了改进混合算法在寻优速度和精度上的优越性。

当得到系统参数和噪声PDF的估计值后,可以很容易计算出性能指标,如表3所示。BC是通过优化算法求得的基于巴氏系数的性能指标;IBC为其理论值;IME是基于最小有理熵的性能指标;IMV为最小方差指标。根据表3,无论对于高斯还是非高斯噪声,新指标都能准确地评估系统性能,且与理论值的差异很小。对比最小熵指标,新指标在计算过程中更为简便,实时性更强,更重要的是有效避免了熵值平移不变性的缺点,具有更好的实际意义。

表3 不同噪声扰动下的性能指标

4 总结

针对最小熵性能指标的缺陷,克服巴氏距离无界性的问题,提出了一种基于巴氏系数的随机性能指标IBC。高斯扰动下,该指标可通过输出数据的均值和方差直接求得;非高斯情况下,利用离散化处理获得概率分布后求出指标近似值。为了更加快速准确地获得基准输出,对传统EDA算法进行了改进,通过引入初步估计和WSA算法思想,极大地提高了算法的效率。仿真结果也验证了HEDA算法在寻优速度和精度上的优越性。下一步研究可将新指标和寻优算法用于串级、非线性等系统的性能评估中。