现代控制论文-系统辨识-贝叶斯

基于贝叶斯方法的系统辨识摘要：建立动态系统在状态空间的模型，使其能正确反映系统输入、输出之间的基本关系，是对系统进行分析和控制的出发点。由于系统比较复杂，往往不能通过解析的方法直接建模，而主要是在系统输入输出的试验数据或运行数据的基础上，从一类给定的模型中确定一个与被研究系统本质特征等价的模型。如果模型的结构已经确定，只需确定其参数，就是参数估计问题。本文首先从原理上推导出贝叶斯参数估计的方法，然后建立一个方阵模型，最后用贝叶斯方法识别出该模型的参数。关键字：状态空间、参数识别、贝叶斯、最大后验概率估计一、贝叶斯的起源和基本原理1.1贝叶斯的起源贝叶斯ThomasBayes，英国数学家.1702年出生于伦敦，做过神甫。1742年成为英国皇家学会会员。1763年4月7日逝世。贝叶斯在数学方面主要研究概率论。他首先将归纳推理法用于概率论基础理论，并创立了贝叶斯统计理论，对于统计决策函数、统计推断、统计的估算等做出了贡献.1763年发表了这方面的论著，对于现代概率论和数理统计都有很重要的作用。贝叶斯的另一著作《机会的学说概论》发表于1758年。贝叶斯所采用的许多术语被沿用至今。贝叶斯（ReverendThomasBayes1702-1761）学派奠基性的作品是贝叶斯的论文“关于几率性问题求解的评论”。或许是他自己感觉到他的学说还有不完善的地方，这一论文在他生前并没有发表，而是在他死后，由他的朋友发表的。著名的数学家拉普拉（Laplace,P.S）用贝叶斯的方法导出了重要的“相继律”，贝叶斯的方法和理论逐渐被人理解和重视起来。但由于当时贝叶斯方法在理论和实际应用中还存在很多不完善的地方，因而在十九世纪并未被普遍接受。二十世纪初，意大利的菲纳特（B.deFinetti）及其英国的杰弗（Jeffreys,H.）都对贝叶斯学派的理论作出重要的贡献。第二次世界大战后，瓦尔德（Wald,A.）提出了统计的决策理论，在这一理论中，贝叶斯解占有重要的地位；信息论的发展也对贝叶斯学派做出了新的贡献。1958年英国最悠久的统计杂志Biometrika全文重新刊登了贝叶斯的论文，20世纪50年代，以罗宾斯（Robbins,H.）为代表，提出了经验贝叶斯方法和经典方法相结合，引起统计界的广泛注意，这一方法很快就显示出它的优点，成为很活跃的一个方向。在这里值得一提的是，八十年代以后，人工智能的发展，尤其是机器学习、数据挖掘的兴起，为贝叶斯理论的发展和应用提供了更为广阔的空间。1.2贝叶斯基本原理贝叶斯辨识方法的基本思想是把所要估计的参数看作随机变量，然后设法通过观测与该参数有关联的其他变量，以此来推断这个参数。核心思想为：()()()()PBAPAPABPB,(,)()(,)()PBACPACPABCPBC由已知的先验概率()PBA求出()PAB设μ是描述某一动态系统的模型，θ是模型μ的参数：(),(),(1),(1),,(1),(1),(0),(0)kzkukzkukzuzuZ1(1),(1),(2),(2),,(1),(1),(0),(0)kzkukzkukzuzuZ1(),(),kkzkukZZ如果系统的输出变量()zk在参数θ及其历史记录1kZ下的概率是已知的，并假定u(k)是一个已知的确定性序列，可将其忽略，那么利用Bayes公式可得：1111(()|,)(|)(|)(|(),)(()|)kkkkkpzkpppzkpzkZZZZZ(1)当k=1时，有：00100((1)|,)(|)(|)(|(1),)((1)|)pzzpzppzzpzzZ(2)原则上说，根据(1)式可求得θ的后验概率密度函数，但实际上这是有困难的，只有在参数θ和数据之间的关系是线性的，噪声又是高斯分析的情况下，才有可能得到(1)式的解析解。求得参数θ的后验概率密度后，就可以利用它进一步求的参数θ的估计值。常用的方法有两种，一种是极大后验参数估计方法，另一种是条件期望参数估计方法，这两种方法统称为贝叶斯方法。二、贝叶斯参数估计方法2.1极大后验参数估计方法极大后验参数估计方法就是把后验概率密度函数(|)kpZ，达到极大值作为估计准则。问题的提出：(|)log(|)kkMaxporMaxpZZ需满足：ˆlog(|)0MAPkpZ(3)将(2)式代入(3)可得：11ˆˆlog(()|,)log(|)0MAPMAPkkpzkpZZ(4)若让式（4）第一项取0，则对应的估计就是极大似然估计。可见，极大后验估计比极大似然估计多考虑了参数θ的先验概率知识。通常情况下，若参数θ的先验概率已知，则极大后验估计优于极大似然估计，即精度更好。但是通常由于没有参数θ的先验概率知识，加之参数θ的后验概率密度函数计算难度较大，因此极大似然估计仍比极大后验估计使用得更普遍。2.2条件期望参数估计方法条件期望参数估计方法就是以参数θ的条件数学期望作为参数估计值。ˆ|(|)kkEpdZZ(5)可以看出上式的物理意义是用随机变量的均值作为它的估计值。三、Bayes应用于系统辨识3.1最小二乘模型的Bayes参数辨识残差表示为：()()Tkkzkθ其中：11()[(1),,(),(1),()][,,,,,]TTnmkzkzknukukmaabb利用两种方法都需先求θ的后验概率密度函数，利用贝叶斯公式，得：111((),)()()(())kkkkPzkZPZPZPzkZ(6)假定的后验分布是高斯分布，在样本Zi情况下，均值为ˆi，协方差阵为Pi。其中0≤i≤k-1。由上述分析可得式(7)、(8)、(9)。/2100000(2)1ˆˆ(|)exp2NTpPPZ(7)/2111111(2)1ˆˆ(|)exp2NTkkkkkpPPZ(8)()()Tkkzkθ(9)假定的先验分布是高斯分布，在、Zk-1已知情况下，z(k)的均值为Tk，方差为2，则有：212211(()|,)exp()22TkkpzkzkZ(10)将(7)、(8)、(10)代入(6)式，可得：21111211ˆˆ(|)exp()22TTkkkkkpNormzkPZ(11)其中Norm是与θ无关的数，如下式：/2211(2)(()|)2detNkkNormpzkPZ(12)将(11)式两边取对数后分别对θ求导，得：1112log(|)1ˆ()0TkkkkkpzkPZ(13)由(13)式计算得：111112211ˆ()TkkkkkkPPzk(14)令11121TkkkkPP，则有：1121TkkkkkPPIP(15)1121ˆˆ()TkkkkkPzk(16)令21kkkKP，则得出：11ˆˆˆ()TkkkkkKzk(17)将式(15)右乘Pk-1，可得：11TkkkkkPPKP(18)将Pk代入Kk表达式，可得：121kkkTkkkPKP(19)配方法求参数的过程：222211111111111121112log(|)()()()1ˆˆˆˆTTkkkTTTTTTkkkkkkkkkkpConstzkzkzkPPPPZ2111112111111221111111221ˆˆ()11ˆˆ()()11ˆˆ()()TTkkkkTTkkkkkkTkkkkkkkkkkkConstzkPPzkPzkPConstPzkPPPPzkPP3.2贝叶斯估计收敛性利用21kkkKP和式(17)有：1111112ˆˆˆ()(){[()]}ˆˆ()[()]ˆ(1)[()]ˆ(1)[()]TTTkkkkkkkTTTkkkkkkTTkkkkTTkkkkkzkzkKzkzkKzkKzkPzk因为()Tkkzk，且：222ˆˆ()[(())][][()]ˆˆˆˆ[()()][()()]TTTkkkkkkkTTTTTkkkkkkkkkkkkkkEEDzkEEEEP由上两式可知，经过辨识过程的数据处理后，残差按照期望的方向迅速减小。四、MATLAB仿真4.1.仿真模型模型的选取如图1所示。图中v(k)是服从N(0,1)分布的不相关随机噪声。ddV(k)y(k)u(k)e(k)z(k)++N(z-1)G(z-1)图1.系统模型111()()()BzGzAz111()()()DzNzCz1121112112()11.50.7()()1.00.5()10.2AzzzCzBzzzDzzz上述模型状态方程的推导如下：11()()()()()ukGzvkNzzk可得：1111()()()()()()()BzDzukvkzkAzCz因为A(z-1)=C(z-1)，于是：111()()()()()()ukBzvkDzzkAz将Az-1)、B(z-1)、C(z-1)的表达式代入上式：121212()(1.00.5)()(10.2)()(11.50.7)ukzzvkzzzkzz1212121.0()0.5()()()0.2()()1.5()0.7()ukzukzvkvkzvkzzkzkzzkz1.0(1)0.5(2)()(1)0.2(2)()1.5(1)0.7(2)ukukvkvkvkzkzkzk4.2程序流程图如图2所示，模型结构选取如下方式121213()(1)(2)(1)(2)()(1)((2)zkazkakbukbukvkdvkdvk工作间清零，给M序列的长度L赋值用4位移位寄存器产生幅值为1输入信号产生随机噪声信号并计算其方差画出输入信号杆线图及随机噪声图形产生输出采样信号给被识别的ϴ和P赋初值按照公式计算K(k)按照公式计算K(k)按照公式计算K(k)计算系统的实际输出响应及模型响应计算被识别参数的收敛情况满足要求画出被识别参数、收敛情况、输出采样、系统实际输出、模型输出等数值停机NY程序流程图图2贝叶斯辨识结果如表1：参数1a2a1b2b1d2d3d真值-1.50.71.00.51.0-1.00.2估计值-1.50.71.00.51.0-1.00.2表1五、仿真结果及分析1.仿真结果如图3-6：2.结果分析：仿真结果表明，Bayes最小二乘辨识，递推到第9步时，参数辨识的结果达到稳定状态，即a1=-1.5000,a2=0.7000,b1=1.0000,b2=0.5000,d1=1.0000,d2=-1.0000,d3=0.2000。此时，辨识参数的相对变化量E0.000000005。可见，此例的Bayes最小二乘递推算法的辨识结果与增广最小二乘递推算法的辨识结果是一致的。程序中，预先把参数的估计值矩阵、参数的收敛状况矩阵及系统和模型的响应矩阵的初始值均设定为零，所