最小二乘辨识方法的优劣比较

1、最小二乘辨识方法的优劣比较摘要：本文系统的探讨了三种最小二乘类辨识方法的原理和性能，并对各种方法在各种不同的环境下进行了MATLAB仿真，仿真结果证明：最小二乘法不适合实时处理，在同等情况下，递推最小二乘的辨识速度较快，但在有色噪声干扰下效果不理想，广义最小二乘法的辨识效果最好，且不受噪声是否有色的影响，但是费时最多。关键词：最小二乘辨识速度MATLAB仿真1引言系统辨识是一门介于现代控制理论和系统理论的边缘学科它将现代控制论的平滑、滤波、预测和参数估计理论，以及系统论的系统分析方法和建模思想应用于自然科学、社会科学和工程实践中的各个领域，与各个领域的专业知识相给合，形成了一个个新的交叉学科分

2、支。关于系统辨识的含义，早在1962年Zacleh曾作如下定义：“根据系统的输入和输出，在指定的一类系统中确定一个相被辨识系统等价的系统”。根据这个定义，在系统辨识中必须确定三方面的问题；第一，必须指定一类系统即根据先验信息确定系统模型的类型。第二，必须规定一类插入信号。例如正弦信号、阶跃信号、脉冲信号、白噪声、伪随机信号等。而且这些信号从时域考虑，必须能持续地激励系统的所有状态；从频域考虑，输入信号的频带能覆盖系统的频带宽度。第三，必须规定“系统等价”的含义及其度量准则。2线性系统的辨识2.1 问题描述考虑如下线性系统：z(k)+az(k-1)+LL+az(k-n)=bu(k-1)+LL+b

3、u(k-n)+e(k)(1)1na1nbab其中，u(k)为系统激励信号，y(k)为系统输出，e(k丿为模型噪声。其系统模型如图1所示:图1SISO的系统模型结构图其中G(z-i)是系统函数模型，N(z-i丿为有色噪声系统模型，e(k)为白噪声v(k)经过系统函数为N(z-1)的系统后的输出。通常(、BQ-1/(、DQ-1/2)Gz-i=(),Nz-i=()Az-1Cz-1A(z-i)=BC-J=1+az-i+az-2+L+az-nana+bz-nbnb12bz-i+bz-2+L123)CC-J=1+cz-i+cz-2+L4)5)()12Dz-i丿=dz-i+dz-2+L+bz-ndI12()

4、B(z-1)()D(z-1)()则系统可表示为：z卜右v(k)设样本和参数集为：h(k)=-z(k-1),-z(k-2),-z(k-n),u(k-1),u(k-2),u(k-n)T(6)0=a,a,a,b,b,bt12n12nh(k)为可观测的量，差分方程可写为最小二乘形式z(k)=hr(k)0+e(k)(7)如何系统噪声e(k)存在的情况下从该方程中正确的解出0，即是系统辨识的任务。为了求出0，我们面临三大问题：一是输入信号的选择，二是判决准则的选取，三是辨识算法的选择，下面一一探讨。2.2 选择输入为了准确辨识系统参数，我们对输入信号有两大要求，一是信号要能持续的激励系统所有状态，二是信号

5、频带能覆盖系统的频带宽度。除此之外还要求信号有可重复性，不能是不可重复的随机噪声，因此我们通常选择M序列或逆M序列作为输入。2.3 准则函数因为本文主要探讨最小二乘类辨识方法，在此选取准则函数J(0)=艺e(k)2=艺z(k)-hT(k)02(8)k=1k=1使准则函数J(0)=min的0估计值记做0，称作参数0的最小二乘估计值。LS在式(7)中，令k=1,2,3,L可构成线性方程组：z(k)=Ht(k)0+e(k)(9)LLLz(1)1e(1)1z(2)e(2)z=,e=LMLM式中_z(L)_e（L）_-z(0)L-z(1-n)u(0)Lu(1-n)1-z(1)Lz(2-n)u(1)Lu(

6、2-nb)H=abLMMMM_-z(L-1)L-z(L-n)u(L-1)Lau(L-n)b准则函数相应变为：J)=e(k)2=z(k)-hT(k)&2=(zk=1k=110)11)-H0)T(z-H6)LLL极小化J（6）,求得参数6的估计值，将使模型更好的预报系统的输出。2.4 辨识算法常用的最小二乘类辨识方法有下列三种：最小二乘法，递推最小二乘法和广义最小二乘法。2.4.1最小二乘法设6使得J（6）=min，贝V有aj(6)a66lsa6=6(z-H6)t(z-H6)=0LLLL展开上式，并根据以下两个向量微分公式：a()aTx丿=aTaxCxtAx)=2xtAdxA为对称阵得正贝

7、方程：（HTH）6LLLS）当HTH为正贝阵时，有6LL=HTzLL=(HTH)-1HTzLSLLLL12)13)14)15)LSa2j（6）且有6ls=2HtH>0，所以满足式（15）的6唯一使得J（6）=min，这LLLS种通过极小化式（11）计算6的方法称作最小二乘法。而且可以证明，当噪声e（k）是均值LS为0的高斯白噪声时，可实现无偏估计。2.4.2递推最小二乘算法为了减少计算量，减少数据在计算机中占用的存，并实时辨识出系统动态特性，我们常利用最小二乘法的递推形式。下面我们来推导递推最小二乘算法的原理。首先，将式(11)的最小二乘一次完成算法写为9)=(HtHLHtz=P(L)H

8、tzWLSLLLLLL=Yh(i)hT(i)|T|Fh(i)ZC)i=1i=1定义P-1(k)=HTH=Ykh(i)hT(i)kki=1P-1(k-1)=HTH=Yh(i)hT(i)k-1k-1i=1式中H=k-hT(1)hT(2)MhT(k)H=k-1-hT(1)hT(2)MhT(k-1)式中，h(i丿是一个列向量，也就是址的第i行的倒置，P(k丿是一个方阵,于未知参数的个数，假设未知参数的个数是n则P(k丿的维数是nXn.o由式17可得P(k)的递推关系为：P-1(k)=匸h(i)hT(i)+h(k)hT(k)i=1=P-1(k-1)+h(k)hT(k)设z=z(l),z(2),L,z(k

9、-1)TJk-1匚z=z(l),z(2),L,z(k)Tk则)()9(k-1)=(HtHLHtzk-1k-1k-1k-1=P(k-1)匸h(i)z(i)i=1由此可得：P-1(k-1)9(k-1)=Yh(i)z(i)16)17)(18)他的维数取决19)20)21)22)i=124)由式19和22可得0(k)=(HrHLHrz=P(k)fh(i)z(i)kkkkI=P(k)P-1(k-1)0(k-1)+h(k)z(k)=P(k)(P-1(k)-h(k)hr(k)0(k-1)+h(k)z(k)=0(k-1)+P(k)h(k)Iz(k)-hr(k)0(k-1)1引进增益矩阵K(k)，定义K(k)=

10、P(k)h(k)式(23)可以进一步写为0(k)=0(k-1)+K(k)Iz(k)-hr(k)0(k-1)1接下来可以进一步把式(20)写为P(k)=I24)25)P-1(k-1)+h(k)hT(kj利用矩阵反演公式(A+CCT)_1=A-1A-1CVi+CTA-1CCTA-1将式(26)演变成P(k)=P(k-1)-P(k-1)h(k)hT(k)P(k-1)hr(k)P(k-1)h(k)+1_P(k-1)h(k)hr(k)-hr(k)P(k-1)h(k)+1-1-1P(k-1)将上式代入式(24)，整理后可得K(k)=P(k-1)h(k)hr(k)P(k-1)h(k)+1_综合式25、)27

11、和28可得最小二乘递推参数)古计算法RLS0(k)=d(k-1)+K(k)z(k)-hr(k)0(k-1)K(k)=P(k-1)h(k)P(k)=-1hT(k)P(k-1)h(k)+1I-K(k)hT(k)1P(k-1)-126)27)(28)29)2.4.3广义最小二乘法设SISO系统采用如下模型：A(z-1)z(k)=B(z-1)u(k)+1)v(k)30)假定模型阶次n,nb和n已知，用广义最小二乘法可以得到无偏一致估计。令abcz(k)=C(z-1)z(k)u(k)=CQ-Ju(k)32)33)0=a,a,L,a,b,b,L,bt12n12nabh(k)=-z(k-1),L,-z(k-

12、n),u(k-1),L,u(k-n)Tlf1-ffaffb将模型化为最小二乘格式：z(k)=hT(k)9+v(k)由于v(k丿是白噪声，所以用最小二乘可以获得参数。的无偏估计，由于噪声模型C(z-1)未知，还需要用迭代的方法来求得C(z-i)。令34)e(k)=(Z)v(k)置35)0(k)=c,c,L,ctVe12nch(k)=-e(k-1),L,-e(k-n)Tc36)这样就把噪声模型也转变为最小二乘格式：e(k)=hT(k)0+v(k)ee由于上式中的噪声已为白噪声，所以用最小二乘也可获得参数0e的无偏估计，但是数据向量中依然含有不可测的噪声量-e(k-1),L,-e(k-作)，可用相应

13、的估计值来代替,置h(k)=，其中k0时，e(k)=0k0时，按照37)计算，式中38)h(k)=-z(k-1),L，一z(k-n),u(k-1),L,u(k-n)_|Tab综上所述，广义最小二乘法可归纳为0(k)=0(k-1)+K(k)z(k)-hT(k)0(k-1)fLffK(k)=P(k-1)h(k)hT(k)P(k-1)h(k)+1_ffff(k)=LI-K(k)hT(k)1P(k-1)0(k)=0(k-1)+K(k)Le)(k)-hT(k)0(k-1)1eeeeeK(k)=P(k-1)h(k)LhT(k)P(k-1)h(k)+1eeeeeeP(k)=LeeeeeeI-K(k)hT(k

14、)1P(k-1)eee-1-139)3仿真研究已知系统模型:x(k)-l5x(k-l)+07x(k-2)=u(k-l)+O5u(k-2),y(k)=x(k)+v(k),v(k)=ay(k),u、x、y、V分别为模型输入、模型输出、测量输出、干扰噪声。输入u为逆m序列：信号幅值a=1、寄存器位数为n=5(信号长度N=2n-l=31),重复周期数q=40。a为噪信比调整因gao子，噪信比定义为：NSR=iX100%=X100%,o、g分别为模型输出x和噪声VooxVxx的均方差(标准差)，Y有两种模型：(1)Y为白噪声，(2)Y为有色噪声，噪声模型为：Y(k)=e(k)+0.5e(k-1)+0.9

15、Y(k-1)-0.95Y(k-2),e(k)为白噪声。定义辨识误差值：5=1卜£牛甲x100%，其中：N为独立的实验次数,9°为模型真值，f为估计值。N町0i=103.1 产生输入和输出数据选择自相关特性好的逆m序列作为输入。利用MATLAB产生寄存器位数n=5，每周期长为31，重复周期数q=40的逆m序列，并将其作为输入得到系统输出。绘出一个周期的图3逆m序列经过系统的输出毬悅静挪产主曲帆甲團输入输出图形分别如图2和图3所示。图2寄存器位数为5的逆m序列3.2 产生系统噪声为了后面能较好的区分每种辨识方法的性能，我们分别在输出中叠加白噪声和有色噪声。取NSR=20%，用同

16、一噪声源产生两种噪声模型，分别绘制白噪声、用相同噪声模型产图4方差为0.2的高斯白噪声图5白噪声影响下的系统输出图7有色噪声影响下的系统输出图6白噪声经过噪声模型后产生的有色噪声3.3 最小二乘辨识模型辨识的均值。图8白噪声对各系数的辨识精度影响图9有色噪声对各系数的辨识精度影响为较好的研究最小二乘辨识模型的性能，作者分别在不同的噪声模型下，用不同的噪信比影响系统输出，利用输入输出数据对系统进行辨识。V分别采用白噪声模型和有色噪声模型，取NSR=O%、5%、10%、15%、20%、25%、30%、35%、40%、45%、50%，每种工况下取独立试验次数N=50（每次独立产生噪声），数据序列取前

17、1024点，用最小二乘法辨识模型，分别画出NSR5曲线（图8和图9）。图中的纵坐标（辨识误差）是50次辨识误差由图可见： 在白噪声影响下，各系数的辨识误差都很小，欲辨识参数为al=-1.5,a2=0.7,b1=2,b2=0.5，即使在噪信比为50%的情况下，四个参数的辨识误差都在10-3数量级，相对误差非常小，均可视为无偏估计，与理论相符。 在有色噪声影响下，各系数的辨识误差相对白噪声影响偏大，当噪信比达到50%时，其中，al、a2和bl的辨识误差都在0.020.04之间，相对误差在10%左右，b2的辨识误差甚至达到0.16以上，相对误差达到30%以上。综上所述：在只有白噪声影响下，最小二乘辨

18、识法可以达到无偏估计，但是在有色噪声影响下辨识结果的相对误差较大。最小二乘法只适合用于只有白噪声影响下的系统辨识，对于有色噪声影响下的系统，我们应该寻求更好的辨识方法。3.4递推最小二乘辨识模型辨识v分别采用白噪声模型和有色噪声模型，取NSR=10%。、40%，用递推最小二乘法辨识模型参数，对比画出各参数辨识结果随递推次数变化的曲线。为了对比研究，我们在同一组u、x序列下，用同一白噪声源卩产生给定噪信比的白噪声和有色噪声干扰。欲辨识参数为a1=-1.5,a2=0.7,b1=2,b2=0.5,设定在两种辨识情况下，前后两次辨识误差小于0.000000005时，结束仿真，当设定NSR=0.1时，本

19、次仿真循环35次时结束仿真，仿真结果见图10图13所示。图10白噪声影响下参数辨识结果图11有色噪声影响下参数辨识结果图12白噪声影响下的参数辨识误差图13有色噪声影响下的参数辨识误差(NSR=10%时)修改参数NSR=0.4,其他条件不变，欲辨识参数为a1=-1.5,a2=0.7,b1=2,b2=0.5,再次运行仿真，循环35次以后结束仿真，仿真结果见图14图17所示Hcwjir菲弐酸祎t吾序戡曲讣砧舉凤谨越威嚴旳芟R.iET£-IJR也图14白噪声影响下参数辨识结果图15有色噪声影响下参数辨识结果图16白噪声影响下的参数辨识误差图17有色噪声影响下的参数辨识误差（NSR=40%时

20、）从仿真结果图中我们可以看到，当噪信比较小时（如NSR=0.1）时，在白噪声影响下，递推最小二乘能正确辨识系统参数，且辨识曲线较平稳。而在有色噪声影响下，辨识结果有一点误差，但辨识曲线尚较平稳。当噪信比较大（如NSR=0.4）时，不管是在白噪声还是有色噪声影响下，辨识曲线的波动都较大，且辨识误差都比噪信比小时的辨识误差有所增加。综上所述：递推最小二乘法只适合与噪信比较小时的白噪声影响下的系统辨识，对于有色噪声影响下的系统辨识或者噪信比较大时的白噪声影响下的系统辨识，我们应该寻找更好的辨识方法。3.5 广义最小二乘辨识模型3.5.1广义最小二乘与递推最小二乘的性能比较v采用有色噪声模型，取NSR

21、=10%、30%，分别用递推最小二乘和广义最小二乘递推法辨识系统参数，为了对比研究，我们在同一组u、y序列下进行辨识试验，并画出各参数辨识结果随Y次数变化的曲线。图18和图19是取NSR=10%时分别采用递推最小二乘法和广义最小二乘法的辨识结果，图20和图21分别是其相对辨识误差。IrtW色垮戶昨鳞下超梅小二曲注金1003DD400HDDBOOTW!BDC900ICEDHEM-litl,包#1惠审T"迟HHteS炳豪曲Hiif?it車looyjDmormormsddsnoiocd图18递推最小二乘法的辨识结果图19广义最小二乘法的辨识结果图20递推最小二乘法的辨识精度图21广义最小二

22、乘法的辨识精度保持其他参数不变，取NSR=30%，再次运行程序可得图22图25。讯匸氏wrt斉色曙齊卡科下*通枯龜沪二出扫it法评.宰鼻Hsp-otiffiF厂反轨丰二总通出乐旺輕曲車ido3co<sn-®ounosac*rwaccxioiocd疇声妍愛“昂厲阶衣图22递推最小二乘法的辨识结果图23广义最小二乘法的辨识结果图24递推最小二乘法的辨识精度图25广义最小二乘法的辨识精度由图可见：两种方法在较小噪信比的有色噪声影响下均能较好的辨识出系统参数，把相对误差曲线图锁定在（-0.02，0.02）的区间上进行观察可得：广义最小二乘法能更快更平坦的接近辨识结果，辨识效果更好。当有

23、色噪声强度增加到30%时，两种方法的辨识结果均出现了不同程度的失真，尤其是递推最小二乘法，结果失真严重且辨识曲线波动较大，广义最小二乘法的辨识结果虽有失真，但仍比较接近真实结果，且辨识曲线平坦。因此我们可以得出结论：由于引入了对噪声的估计，在噪信比较大的有色噪声影响下，广义最小二乘法仍然能够得到较好的辨识结果。3.5.2广义最小二乘的辨识性能与噪声模型的最高阶次在上面的广义最小二乘法的辨识过程中，我们取的噪声模型的最高阶次为4，如果改变噪声模型的最高阶次结果会有何改变呢？为此，作者编写程序得到了在噪信比为30%的有色噪声影响下，广义最小二乘法进行辨识时，噪声模型阶次依次从2开始以2为步长递增到

24、8时的辨识结果曲线，如图26图29所示。图27噪声模型阶次为3时的辨识结果图26噪声模型阶次为2时的辨识结果图28噪声模型阶次为4时的辨识结果HEM-stt可削nita-t-"迟HHteS炳豪曲Hiif?it車图29噪声模型阶次为6时的辨识结果由图可见，在噪声模型最高阶次过小（为2）的情况下，辨识结果精度较低且在辨识初期波动较大，随着模型最高阶次的增加，辨识结果精度增加，且辨识过程曲线平坦，但增加到一定程度后辨识曲线无明显区别。但是噪声模型的阶次越高，辨识过程的运算量就越大，因此我们在实际仿真时要综合考虑，选取噪声模型阶次能较好的辨识出系统参数即可，不必一味追求过高的模型阶次。4结论

25、从本文的介绍和仿真结果可以看出，系统辨识常用的三种最小二乘类辨识方法中，递推最小二乘法的算法简单，能减少计算量，减少数据在计算机中占用的存，并实时辨识出系统动态特性。但只能在噪信比较小的白噪声影响下，才能精确的辨识出系统参数。广义最小二乘法算法较为复杂，运算量较大，但是能在系统噪声为有色噪声时，仍能较好的辨识出系统参数。实际仿真中，如果我们合理的选取噪声模型的最高阶次，可以达到较好的辨识效果且运算量也较小。附源程序：clearcloseall;L=500;%移位寄存器初值x1=1;x2=1;x3=1;x4=0;s=l;%方波初值un=zeros(L,1);fori=1:L+4fork=3:L+

26、4;IM=xor(s,x4);ifIM=0%进行异或运算，产生逆M序列u(i)=-1;%fixedelseu(i)=1;endun(i,1)=u(i);S=not(S);%产生方波M(i)=xor(x1,x4);x4=x3;x3=x2;x2=x1;x1=M(i);%进行异或运算，产生M序列%寄存器移位endz(2)=0;z(1)=0;%取z的前二个初始值为零z(k)=1.5*z(k-1)-0.7*z(k-2)+1.0*u(k-1)+0.5*u(k-2)+normrnd(0,1,1,1);%给出理想的辨识输出采样信号zn(k,1)=z(k);endc0=0.0010.0010.0010.001&

27、#39;%ala2blb2给出被辨识参数的初始值p0=10"6*eye(4,4);E=0.000000005;%相对误差E=0.000000005c=c0,zeros(4,499);%被辨识参数矩阵的初始值及大小e=zeros(4,500);%相对误差的初始值及大小fork=3:500;%开始求Khl二z(k1),z(k2),u(k1),u(k2)；x二h1'*p0*h1+1;xl二inv(x);%开始求K(k)kl二p0*h1*x1;%求出K的值d1=z(k)h1'*c0;c1=c0+k1*d1;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2

28、=e1./c0;%求参数的相对变化e(:,k)=e2;%把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列p1=p0k1*k1'*h1'*p0*h1+1;%求出p(k)的值p0二p1;%给下次用ifabs(e2)<=Ebreak;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束%c%显示被辨识参数%e%显示辨识结果的收敛情况%e2;a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);i=1:500;plot(i,al,'r',i,a2,'',i,bl,'g:',i,b2,'k-.','LineWidth',2)%画出al,a2,a3,bl,b2的各次辨识结果title('