2014_系统辨识理论及应用_萧德云_2

【实例简介】
萧德云的系统辨识理论及应用，本文讲的很全，不过感觉本文偏重辨识算法，输入信号的设计、最优实验设计没讲，不得不说有点缺陷，其次，感觉文笔一般，没有太让人想看的欲望，丁锋出了更全的系统辨识，不过暂时没搞到全的，好像也只出了第一、三册。 由于文件大，分1,2两部分，这是第2部分。
第11章增广UD分解辨识算法 285 6,(k)和1(k)对应的损失函数J,(k),r=0,1,2,…,n和V(k),r=0,1,2,…, n-1。当选定可能的最高阶次n时,算法可同时获得1至n阶的模型辨识结果,包括 1至n阶模型的损失函数和参数估计向量。模型式(11,4.2)的辨识结果是我们所关 注的,模型式(11.4.5)的辨识结果没有什么实用价值,可以不予关注。 例11.3考虑如下SlSO线性差分方程模型 z(k)=1.5z(k-1)—0.7z(k—2)+1.0u(k-1)+0.5l(k一2)+Av(k) (11.4.22) 其中,z(k)、u(k)是模型的输入和输出变量;υ(k)是均值为零、方差为1的白噪声,λ 为噪声标准差。当模型输入选用特征多项式为F(s)=ss⊕1、幅度为1的M序 列,噪声标准差取λ=0.3(噪信比约为20.5%)时,可以仿真获得模型的输入和输出 数据(数据长度L=300)。假设模型可能的最高阶次n=4,根据式(11.4.4)和 式(11.4.9)定义,利用仿真数据构造如下的增广数据向量和信息压缩阵 -x(k-4),(k-4),一x(k-3),k-3),一x(k一2),(一2),一z(-1),a(k-1) C4(300) (iof(i) 0.01840.0118-0.03900.02530.0270-0.0002-0.00250.00120.0012 0.01180.0144-0.02740.01900.02120.0002-0.00340.00190.0018 0.0390 0.0274 0.1018-0.0428-0.0974 0.0256 0.0320-0,0038-0.0036 0.02530.0190-0.04280.05020.0114001880.01570,00140,0014 0.02700.02120,09740.01140.1434-0.04270.09960.02990.0291 0.00020.00020.02560.01880.04270.05010.01020.01810,0178 0.0025-0.00340.03200.0157-0.09960.01020.1219-0.0574-0.0561 0.00120.0019-0.00380.00140.02990.0181-0.05740.04180.0376 0.00120.00180.00360.00140.02910.0178-0.05610.03760.0368 (11.4.23) 利用AUDI辨识算法式(11.4.10)和UDU分解 MATLAB程序式(11.3.4),可 以获得如下的参数辨识矩阵U4(300)和损失函数矩阵D4(300) 0000-0.0056-0.90200.00410.68610.02020.01900,00390.0314 01.00001.00720.00550.51530.00180.01720.00980.0489 01.0000-0.0082-1.48980.04340.72710.03270.0974 01.00001.02270.02120.48740.01440.0369 U4(300)= 01.00000.03091.51760.05350.7889 b00 1.0232-0.01840.4831 0000 0.0314 1.5224 01.00001.0217 01.0000 D2(300)=diag[0.00020.00340.00120.00330.03650.00330.03650.00330.0368 矩阵U(300)和D4(300)第9、7、5和3列元素分别为4、3、2和1阶模型的参数估 计向量和损失函数 286 系统辨识理论及应用 0(300)=[a4(300),b(300),a4(300),b1(300,a2(300),b2(300)0a1(300),b1(300) 4阶模型 =[0.0314,0.0489,0.0974,0.0369,0.7889,0.4831,-1.5224,1.0217] J4(300)=27.14 (11.4.24) 62(300)=[a3(300),b2(300),a2(300),2(300)a1(300),b(300) 3阶模型 [-0.0190,-0.0172,0.7271,0.4874,-1.5176,1.0232] J3(300)=27.41 (11.4.25) 02(300)=[a2(300),b2(300)a2(300,(3001 2阶模型 =[0.6861,0.5153,-1.4898,1.0227]T (11.4.26) J2(300)=27.43 61(300)=[a1(300),b1(300)=[一0.9020,1.0072 1阶模型 (11.4.27) J1(300)=825.20 损失函数随阶次变化的曲线如图11.2所示,显然模型阶次大于2后,损失函数 变化不显著,所以模型阶次应该取2。 6000 (564886) 5000 4000 3000 2000 1000 825.20) (2743 2741) (2714) 阶次 图11.2损失函数随模型阶次变化的曲线 最终获得的辨识模型为 z(k)=1.4898z(k-1)-0.6861z(k-2)+1.0227(k-1) 0.5153a(k-2)+0.3008v(k) 该辨识模型与仿真模型式(11.4.22)很接近,说明AUDI算法的辨识效果还是不 错的。 第11章增广UD分解辨识算法 87 仿真例表明,1~4阶模型参数估计值都包含在矩阵U4(300)中,模型过参数的估计 值一般很小(理论上应趋于0),如a4(300)=0.0314,b1(300)=0.0489,a3(300)= 0.0974,b3(300)=0.0369。依据矩阵D4(300)中的J(30),=1,2,3,4,利用第 10章讨论的F-Test定阶法,可以判定模型的阶次,本例确定为2阶,与仿真模型的 阶次是一致的。 另外,如果模型式(11.4.1)输入和输出的参数个数不相等,即a,i=1,2,…,n, b,=1,2,…,n,n≠n,则增广数据向量gn(k),n=max(na,n)中的一些数据将不 以输入输出对成对出现。比如,n=2,n=4,这时n=max(na,n)=4,增广数据向 量定义为q(k)=[(k-4),(k-3),一z(k-2),t(k-2),一z(k-1),(k-1) z(k)],矩阵U4(k)和D4(k)的元素组成也会发生相应变化。 11.5增广UD分解最小二乘辨识算法 考虑如下的ARX模型(最小二乘模型) A(2)2(k)=B(2)u(k)+u(k A(z1)=1+ a2z1+…+ B(z1)=b1z+b2x1+…十b 其中,n(k)、x(k)分别是模型的输入和输出变量;v(k)是均值为零的不相关随机噪 声。若模型阶次取ηnn,利用AUD算法式(11.4.10),完全可以解决该模型 的辨识问题。不过,只能像例11.3那样处理过程是批处理的,不方便实时在线应用。 下面讨论基于增广UD分解的最小二乘递推辨识算法,记作AUDI-RLS(AUDl recursive least squares) 11.51 AUDERLS算法 根据式(11.4.9)信息压缩阵的定义,ARX模型的信息压缩阵Cn(k)可写成如下 的递推形式 Cn(k)=[S,(k-1)+g(k)g(k)] (11.5.2) 利用矩阵反演公式(见附录E.3),式(11.5.2)演算成 Cn(k)=Cn(k-1)一 k-1)e(k) (k)Cm(k-1 1+g(k)C(k-1)g(k) (11.5.3) 根据定理11.1,对Cn(k)和C(k-1)分别进行UDUr分解,式(11.5.3)写成 C(k)=U(k)D(k)UT(k) U(k-1D(k-lUT(k-1) U(k-1)D(k-1)U(k-1)gn()g(k)U(k-1)Dn(一1)Un(k-1) 1+g(k)U(k-1)D,(k-1)U(k-1)P() 288 系统辨识理论及应用 =Un(k-1)D(k-1)U(k-1)Un(k-1)g,(k)B()U(k一1) 1+ fm(k)8n(k) U,(k-1)D,(k-1)g(k)g(R) (k) U(k-1) 其中 gn(k)=[g1(k),g2(k),…,gN(k)]=Dn(k一1)∫n(k) f(k)=[f1(k),f2(k),…,fN(k)=U(k-1)gn(k)(11.5.5) R()=1+厂(k)gn(k) 式中,N=2n+1。再次根据定理11.1,可令 D,(k-1)8,(k)gi(k) BN(k)=U(k)D,(kUT(k) (11.5.6) 则式(11.5.4)写成 C(k)=U(k)D, (k)UT(k) U,(k-1U(k)D (k)U(kUT(k-1) 显然有 D(k)=d(k) (11.5.8) U Ck)=U(k-lU(k) 如果能由Un(k-1)和D(k-1)求得Un(k)和Dn(k),则根据式(11.5.8)即可实 现Un(k)和D(k)的递推计算,由此也就构成了AUDI-RIS算法。 引入N维基向量 e;=[0,…,1,…,0] (11.5.9) 第t元素为1,其他元素为0 并记 (k)u13(k) (k) (k) U(k) N-DN(k) la12(h)u13(k) UI(k) 1 u2(k) 1(k) (11.5.10) (k) k) [a1(k) (k) (k D,(k)= diagld,(k) d2(k) dN-1(k) dN(k)I D(k)= diagld, (k) d2(k) n-i(k) dN(k)] 第11章增广UD分解辨识算法 289 .·:···.:·::··:··::;·“;;;“+,++++,+,a 那么式(11.5.6)可写成 饭,(b)正(k)d1(k)=>d1(k-1)e日一MN (11.5.11) ga (k)g(k) BN(k) 把上式写成更一般的形式 ∑n1(k)aT(h)d,(k)一∑d(k-1)ee-MN-,(1.5.12) 其中,=0,1,…,N-1,且 M gN=,(k)g一(k) BN-(k) (11.5.13) BN-i(k) gn(k),BN-,(k)=1+2f(k)g, (k) =1 式中,当j=0时,gN(k)=g(k)。显然,式(11.5,11)是式(11.5.12)取j=0时的特 例。当j=1时,式(11.5.12)写成 ∑a(k)n(k)d(k)=∑a(k-1)ee-M1(1.5.14) 式中 MN-=gN-1()8N-(k) BN-1(k) (11.5.15) 比较式(11.5.11)和式(11.5.14),可得 MN-1=uN(k)uk(k)dN(k)-dN(k-1)eNeN+En rinck) BN(k) 11.5.16) 根据式(11.5.13)中gN-(k)的定义,矩阵MN-1第N行和第N列的元素均为 0,由此可以导出 (k)=dN(k-1) gN BN(k (k) gN(R)g (k) 1,2,…,N-1 (11,5.17) dN(kBN(k) N() 又根据式(11.5.5)和式(11.5.13),式(11.5.17)进一步演变成 dOck) (k-1)8-(k BAck) u,n (k) N(k)g: (k) (11.5.18) BN-I(k i=1,2,“,N-1 b NN 同理,矩阵MN-,+1第N-j行和第N一j列的元素均为0,与推导式(11.5.17) 和式(11.5.18)类似,U(k)和Dn(k)第N一j列元素可表示为 290 系统辨识理论及应用 a、,(k)=k1(k-1)B=1(k BN-(k) ui(N-n)(k)= fN- (k)8,(k R(6),i=1,2,…,N-1 (11.5.19) (N-)(N-j k)=1 0,1,2,…,N-1 其中,fN-,(k)和g(k)由式(11.5.5)确定,BN-,(k)由式(11.5.13)定义 式(11.5.19)是Un(k)和Dn(k)的递推计算公式,注意到式(11.5.8),即可实现 Un(k)和Dn(k)的递推计算。进一步利用式(11.5.8),并对式(11.5,19)进行下脚标 置换,也就是N-j→j,且j=0,1,2, =1,2,…,N,由此便可构成如下的 AUDI-RLS辨识算法。 AUDI-RLS辨识算法 P(k)=[-z(k-n),(k-n) -1),u(k-1),-z(k) f(k)=U(k-1)gn(k),gn(k)=D(k-1)fn(k) (k)=1+∑f()g(k) (k) f, (k)g (k 月-1(k) u1(k)=1 (11.5.20) n(k)=(k-1)+∑v1(k-1)(k),1(k)=1 d(k)=d(k-1)8-1(k) 月(k i=1,2,…,N;j=1,2,…,N;N=2n+1 上述AUD-RLS辨识算法是一种依时间k的递推计算结构,选定可能的最高模 型阶次n,利用k时刻的增广数据向量φ,(k)和k-1时刻的U(k-1)和D(k-1), 便可递推计算Un(k)和D(),获得1至n阶模型的辨识结果,包括损失函数和模型 参数估计值。 11.52AUD卜RLS算法分析 通过对AUD-RLS辨识算法式(11.5.20)的分析,可以得到以下几个重要的关 系,它们对进一步认识 AUDI-RLS辨识算法有益。 (1)根据式(11.4.2)、式(11.4.4)和式(11.5.5),fn(k)奇数行元素可以表示为 f2r+1(k)=-[z(k)一h(k)日,(k-1)] 一z,(k),r=0,1,2,…,n (11.5.21) 式中,,(k为式(11.4.2)第r阶模型新息。该式表明∫n(k)奇数行元素等于模型新 第11章增广UD分解辨识算法 291 息的负数 (2)根据式(11.5,5),gn(k)奇数行元素可以表示为 g2+1(k)=-z,(k)/J,(k-1),r=0,1,2,…,n (11.5.22) 该式表明gn(k)奇数行元素与式(11.4.2)第r阶模型新息和损失函数有关。 (3)根据式(1.5.13),B(k),=1,2,…,N可以表示为 B(k)=1+∑f1(k)E(k) 1+ 6(k) (k) (k) u2(k) J。(k-1)+v。(k-1)+1(k-1)+v1(k一1 z7-1/2(k) JG;=1)/2(k-1) j为奇数 (11.5.23) l222(k) V,2/(k-1),j为偶数 式中,z,(k)为式(11.4.2)第r阶模型新息,,(k-1)为式(11.4.5)第r阶模型新息。 该式表明A(k)与低于(-2)/2or(-1)/2阶模型新息和损失函数有关。 (4)矩阵D(k)奇数对角元素是式(11.4.2)模型损失函数的倒数,有如下递推计 算形式 d2r+1(k)= d2H(k-1)B2-+1(k) ,r=0,1,2,…,n 11.5.24) B2r (k) 式中,dx+1(k)为式(11.4.2)第r阶模型k时刻损失函数的倒数。根据式(11.5.5), 上式可进一步写成 J,(k-1)B2,+(k) B2n() =1,(k-D)/1、f(kR)+1(k) +∑f(k)g1(k) J,(k-1) (k) 1+∑f(k)g1(k =J,(k-1)+z,(k)E,(k) (11.5,25) 其中 z,(k) e,(k)= ,r=0,1,2,…,n (11.5.26) ∑f(k)g(k) 式中,ε,(k)为式(11.4.2)第r阶模型残差。式(11.5.25)与第5章讨论的最小二乘 损失函数递推计算 292 系统辨识理论及应用 想曹mtt世世+世+量+++++是++,+. J(k)=J(k-1)+z(k)E(k) (11.5.27) e(k)= z(k) 1+h(k)P(k-1)h(k) 具有类似的结构,不过式(11.5.27)的残差E(k)是用k时刻以前的数据信息更新的, 而式(11.5.25)的残差ε,(k)要用到低于r阶的模型信息进行修正。 5)矩阵U(k)奇数列元素是模型式(11.4.2)的参数估计值,第r阶第i个模型 参数估计值具有如下的递推计算形式 L(2r+1) (k)=ua+1,(k-1)+∑un(k-1)unx(k),r=1,2,…,n (11.5.28) 根据式(11.5.5)和n(k)的定义,上式可进一步写成 k)=tg+n(k-1)+K2+1yz,(k) ∑u(k-1)g(D) (11.5,29) Icrt B,(k) 1,2,…,n 式中,z,(k)为式(11.4.2)第r阶模型新息,K2,+为第r阶第i个模型参数估计值 增益。式(11.5.29)与第5章讨论的最小二乘参数递推估计 0(k)=日(k-1)+K(k)z(k) (11.5.30) 也具有类似的结构,不过式(11.530)中的增益矩阵K(k)要用k时刻以前的数据信 息更新,而式(11.5.29)增益K+1的更新与低于r阶的模型信息有关。 11.53 AUDFRFF算法 当 AUDIRLS辨识算法应用于时变系统时,为了跟踪模型参数的变化,需要引 人遗忘因子k。这时式(11.5.2)写成 Cn(k)= Lu(k),(k-1)+on()T(k)J-I (11.5.31) 利用矩阵反演公式(见附录E.3),式(11.5.31)演算成 C,(k) (k)C.(k一1)_Cn(k-1)甲p,(k)P(k)Cn(k-1) u(k)+ pi(k)Cn(k-1)n(k) (11.5,32) 与推导 AUDI-RIS辨识算法式(11.5.20)一样,只要更换两个式子(见式(11.5.33) 和式(11.5.34)),便可构成基于增广UD分解的遗忘因子递推辨识算法,记作 AUDI-RFF(AUDl- recursive forgetting factor),如式(11.5.35)所示,相应的 MATLAB程序也只要在 AUDI-RLS算法程序的基础上稍作修改即可,即 B(k)=1+∑f(k)g,(k)→月(k)=(k)+∑f(k)g(k)(1.1.33 【实例截图】
【核心代码】