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

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

2、域，与各个领域的专业知识相给 合，形成了一个个新的交叉学科分支。关于系统辨识的含义，早在 1962 年 Zacleh 曾作如下定义：“根据系统的输入和 输出，在指定的一类系统中确定一个相被 辨识系统等价的系统” 。根据这个定义， 在系统辨识中必须确定三方面的问题；第 一，必须指定一类系统即根据先验信息 确定系统模型的类型。第二，必须规定一 类插入信号。例如正弦信号、阶跃信号、 脉冲信号、白噪声、伪随机信号等。而且 这些信号从时域考虑，必须能持续地激励 系统的所有状态；从频域考虑，输入信号 的频带能覆盖系统的频带宽度。第三，必 须规定 “系统等价 ”的含义及其度量准则。 2 线性系统的辨识问题描

3、述 考虑如下线性系统： z k a1z k 1 L L an z k na b1u k 1 L L bn u k nb e k(1)其中，u(k)为系统激励信号，y(k)为系 统输出，e(k)为模型噪声。其系统模型如 图 1 所示：图1 SIS o的系统模型结构图 其中G(z-1)是系统函数模型， 有色噪声系统模型，e(k)为白噪声 系统函数为N(z-1)的系统后的输出N(z-1)为 v(k)经过通常1式中:11 a1z1 . 1bizib2za：2D z 1C z 1(2)L anaz na% z nb(3)1C z 1C|D z 1 d1z 1c?zd2z 2LCncL bzzndnc则系

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

5、状态，二是信号频带能覆盖系 统的频带宽度。除此之外还要求信号有可 重复性，不能是不可重复的随机噪声，因 此我们通常选择 M 序列或逆 M 序列作为 输入。准则函数 因为本文主要探讨最小二乘类辨识 方法，在此选取准则函数Jek2zkhT k2k1k1（8） 使准则函数 J min 的 估计值记做 ）LS， 称作参数 的最小二乘估计值。在式（7）中，令k=1,2,3,L,可构成线性方程组：zL k H LT k eL k（9）z1e1z2e2zLM,eLLM式中zLeLz 0 Lz 1 nau 0 Lu 1 nbz 1 Lz 2 nau 1 Lu 2 nbHLMMMMz L 1 Lz L nau

6、L 1 Lu L nb10)准则函数相应变为：L2LT2TekzkhT kzL HLzL H Lk1k1（11） 极小化 J ,求得参数 的估计值，将使 模型更好的预报系统的输出。辨识算法 常用的最小二乘类辨识方法有下列 三种：最小二乘法，递推最小二乘法和广义最小二乘法。2.4.1 最小二乘法设使得J min，则有）一 zLHlT zLHl0LS1-1-1-1-（12）展幵上式，并根据以下两个向量微分 公式：TTax a xA为对称阵xTAx 2xTAx（13）得正则方程： H【Hl）ls H【Zl（14）当h【Hl为正则阵时，有I hTHl 1hTzl（15）2且有|）ls 2H【Hl 0，

7、所以满足式（15） 的s唯一使得Jmin，这种通过极小化式（11）计算）ls的方法称作最小二乘法。而 且可以证明，当噪声 e（k）是均值为0的高 斯白噪声时，可实现无偏估计。2.4.2递推最小二乘算法为了减少计算量，减少数据在计算机中占用的内存，并实时辨识出系统动态特 性，我们常利用最小二乘法的递推形式。 下面我们来推导递推最小二乘算法的原 理。首先，算法写为将式(11)的最小二乘一次完成)WLST 1 THLTHL H LTzLL h i hTi1P L HLTzLLh i z ii116)定义P 1 kHkTHki1Hkh i hT ik1 hihi117)式中HkhThTMhThThTM

8、hTk1(18)HLh(i)是一个列向量，也就是置，P(k)是一个方阵，他的维式中， 的第 i 行的倒 数取决于未知参数的个数，假设未知参数19)的个数是n,则P(k)的维数是nX n.o 由式17可得P(k)的递推关系为：k1h i hTi11P 1 k 1i h k hTh k hT kzkzkTz 1 ,z 2 ,L ,z k 1T1 , z 2 ,L ,z k20)则)kHkT 1HHkThii121)此可1)kk1h i z ii122)由式 19 和 22 可得) k HkT Hk HkT zkPkkhii1PkPk)k1)P 1 k 1 k 1 h k1 T )P 1 k h k

9、 hT k k 1T)1 P k h k z k hT k k 1hkzk引进增益矩阵K(k) ，定义 K k P k h k23)(24)(23 )步写为)k)k 1 K k z k hT k)k 1(25)接下来可以进一步把式(20)写为(30)k 1 h k hT k 1(26)式T 111A CCTA 1 A 1C Icta 1c 1cta 1将式(26)演变成P k P k 1Pk1hkhTkPk1hTkPk1hk 1 1P k 1 h k hT kI 十P k 1hT k P k 1 h k 1(27)将上式代入式(24)，整理后可得K k Pk 1hk hTkPk 1hk 1 1

10、(28)综合式25、27和28可得最小二乘递推参数估计算法RLS)T)kk 1K k z khT k k 1K k P k 1hkhTkPk1hk 1 1P kI Kk hT k P k 12.4.3广义最小二乘法设SIS0系统采用如下模型:B z1 u k(30)假定模型阶次na, nb和nc已知，用广 义最小二乘法可以得到无偏一致估计。令zf k C z 1 z k1uf k C z u k(31)a1,a2,L ,ana,b1,b2,L ,bnbfT hf kZf k 1 ,L , Zf k n& , Uf k 1 ,L , Uf k n(32)将模型化为最小二乘格式：zf k h； k

11、 v k(33)由于v(k)是白噪声，所以用最小二乘 可以获得参数e的无偏估计，由于噪声模 型C(z1)未知，还需要用迭代的方法来求得 c(z1)。令1e kv kC z 1(34)e k c1, c2 ,L ,cnc T he k e k 1 ,L , e k nc (35) 这样就把噪声模型也转变为最小二 乘格式： e k heT k e v k( 36)由于上式中的噪声已为白噪声，所以 用最小二乘也可获得参数9 e的无偏估计， 但是数据向量中依然含有不可测的噪声 量 e k 1 ,L , e k nc ，可用相应的估计值来代 替，置 he k e) k 1 ,L , e) k nc ，其

12、中 k 0 时，按照e) k z k hT k )(37) 计算，式中T h k z k 1 ,L , z k na ,u k 1 ,L ,u k nb38)综上所述，广义最小二乘法可归纳为) ) T ) k k 1 K f k zf k hTf k k 1K f k P k 1 hf k hTf k Pf k 1 hf k 1Pf kPf k 1ekI K f k hTf k ) e k 1 Ke kKekPek 1hekheTkPek1 hek 1Pe kI Ke k heT k Pe k 111(39)3仿真研究已知系统模型：x(k)(k-1)+(k-2)= u(k-1)+(k-2),

13、y(k)=x(k)+ v (k), v (k)= a= (k), u、x、y、v分别为模型输入、模型输出、 测量输出、干扰噪声。输入u为逆m序列： 信号幅值a=1、寄存器位数为n=5(信号长 度N=2n-仁31),重复周期数q=40。a为噪 信比调整因子，噪信比定义为： NSRd 100% 100%, (Tx、分别为模型输出xxx和噪声v的均方差(标准差)，Y有两种 模型：(1) Y为白噪声，(2) 丫为有色噪声， 噪声模型为：丫 (k)=e(k)+(k-1)+ 丫 (k-1) y (k-2),e(k)为白噪 声。定义辨识误差值：3= F杲100%，其N i 1中：N为独立的实验次数,0为模型

14、真值，- 为估计值。产生输入和输出数据选择自相关特性好的逆 m序列作为 输入。利用MATLAB产生寄存器位数n=5, 每周期长为31，重复周期数q=40的逆m序列，并将其作为输入得到系统输出。绘出一个周期的输入输出图形分别如图 2和 图3所示。图2 寄存器位数为 5的逆m序列 图3逆m序列经过系统的输出3.2 产生系统噪声为了后面能较好的区分每种辨识方 法的性能，我们分别在输出中叠加白噪声 和有色噪声。取NSR=20%用同一噪声源 产生两种噪声模型，分别绘制白噪声、用 相同噪声模型产生的有色噪声和不同噪 声影响下的系统输出的曲线，如图4、图5图6和图7所示。叮电1皐出1*(wn 图4方差为的高

15、斯白噪声 图5白噪声影响下的系统输出医唏直VTF肆.IED 2CD XDXX mi TlXhSDD 1XD图6白噪声经过噪声模型后产生的有色噪声图7有色噪声影响下的系统输出3.3 最小二乘辨识模型辨识为较好的研究最小二乘辨识模型的 性能，作者分别在不同的噪声模型下，用 不同的噪信比影响系统输出，利用输入输 出数据对系统进行辨识。V分别采用白噪 声模型和有色噪声模型，取NSR=O% 5%、10%、15%、20%、25%、30%、35%、40%、 45%、50%，每种工况下取独立试验次数N=50 （每次独立产生噪声），数据序列取 前1024点，用最小二乘法辨识模型，分 别画出NSRY曲线（图8和图

16、9）。图中的纵坐标（辨识误差） EIDE DI DH HZ 口二 UJ 3B- 4 D.C- EL SiMlff 拠 此 . *iri( *是 50次辨识误差图8白噪声对各系数的辨识精度影响图9有色噪声对各系数的辨识精度影响 由图可见： 在白噪声影响下，各系数的辨识误差 都很小，欲辨识参数为 a1 = , a2=, b仁2, b2=，即使在噪信比为 50%的 情况下，四个参数的辨识误差都在 10-3数量级，相对误差非常小，均可 视为无偏估计，与理论相符。 在有色噪声影响下，各系数的辨识误 差相对白噪声影响偏大，当噪信比达 到50%时，其中，al、a2和bl的辨 识误差都在之间，相对误差在 10

17、% 左右，b2的辨识误差甚至达到以上， 相对误差达到 30%以上。综上所述：在只有白噪声影响下，最 小二乘辨识法可以达到无偏估计，但是在 有色噪声影响下辨识结果的相对误差较 大。最小二乘法只适合用于只有白噪声影 响下的系统辨识，对于有色噪声影响下的 系统，我们应该寻求更好的辨识方法。3.4 递推最小二乘辨识模型辨识V分别采用白噪声模型和有色噪声模 型，取NSR=10% 40%,用递推最小二乘 法辨识模型参数，对比画出各参数辨识结 果随递推次数变化的曲线。为了对比研 究，我们在同一组 u、x 序列下，用同一 白噪声源丫产生给定噪信比的白噪声和有 色噪声干扰。欲辨识参数为 a1=，a2=，b1=2

18、，b2=， 设定在两种辨识情况下，前后两次辨识误 差小于时，结束仿真，当设定NSR=寸，本 次仿真循环 35 次时结束仿真，仿真结果 见图 10图 13 所示。图10白噪声影响下参数辨识结果图11有色噪声影响下参数辨识结果hBH int理谓柚nJ讨sID13JO33XI削畔勺素4-lfl图12白噪声影响下的参数辨识误差图13有色噪声影响下的参数辨识误差NSR=10%寸）修改参数NSR=,其他条件不变，欲辨 识参数为a仁，a2=, b仁2, b2=,再次运 行仿真，循环35次以后结束仿真，仿真 结果见图14图仃所示图14白噪声影响下参数辨识结果图15有色噪声影响下参数辨识结果yMil.乘赠虑底伸

19、豎此iIBHZ!2KS-“刻&* i *潭*卢更克图16白噪声影响下的参数辨识误差图17有色噪声影响下的参数辨识误差NSR=40% 寸）从仿真结果图中我们可以看到，当噪 信比较小时（如NSR弓时，在白噪声影响 下，递推最小二乘能正确辨识系统参数， 且辨识曲线较平稳。而在有色噪声影响 下，辨识结果有一点误差，但辨识曲线尚较平稳。当噪信比较大（如 NSR=时，不 管是在白噪声还是有色噪声影响下，辨识 曲线的波动都较大，且辨识误差都比噪信 比小时的辨识误差有所增加。综上所述：递推最小二乘法只适合与 噪信比较小时的白噪声影响下的系统辨 识，对于有色噪声影响下的系统辨识或者 噪信比较大时的白噪声影响下的

20、系统辨 识，我们应该寻找更好的辨识方法。3.5 广义最小二乘辨识模型广义最小二乘与递推最小二乘的性 能比较v采用有色噪声模型，取 NSR=10% 30%，分别用递推最小二乘和广义最小二 乘递推法辨识系统参数，为了对比研究， 我们在同一组 u、y 序列下进行辨识试验， 并画出各参数辨识结果随丫次数变化的曲线。图18和图19是取NSR=10%寸分别 采用递推最小二乘法和广义最小二乘法 的辨识结果，图 20 和图 21 分别是其相对 辨识误差。M3M 1IT也購亍谢耳卜二建#|1出幻肇if 帆图18递推最小二乘法的辨识结果图19广义最小二乘法的辨识结果iir ft 电廈* Wh 二-r图20递推最小

21、二乘法的辨识精度图21广义最小二乘法的辨识精度保持其他参数不变，取 NSR=30%再次运行程序可得图22图25M3N FT h色第邛事卜二港.1出细肇島wHrf k国卒第 rr： i二鼻话谨坤护 * jm rr 母m n AJTI .Ibl rfc丄 i Ml划即聞旅图22递推最小二乘法的辨识结果图23广义最小二乘法的辨识结果图24递推最小二乘法的辨识精度图25广义最小二乘法的辨识精度由图可见：两种方法在较小噪信比的 有色噪声影响下均能较好的辨识出系统 参数，把相对误差曲线图锁定在（， ）的 区间上进行观察可得：广义最小二乘法能 更快更平坦的接近辨识结果，辨识效果更 好。当有色噪声强度增加到

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

23、开始以 2 为步长递增到 8 时的辨识结果曲线，如图26图29所示和邑唏rr-： * i i etfi w is w ji/i：Jr #. * .H i l rmi 占 QQ JT* tJTI .Ibl rfcu i M i Ml4tlK图26 噪声模型阶次为2时的辨识结果图27噪声模型阶次为3时的辨识结果董it 0色唏議看# rr.xi怡二鼻诵冷潅見统只2報 *.4 JI ri母 EQJT* tJTI ilhl rfc 丄 i M MlQ-C * ,H i I rmi 当 QQ JT* tJTI ulhl rh 丄 i IB ” Ml土社片图28噪声模型阶次为4时的辨识结果 图29噪声模型阶

24、次为6时的辨识结果 由图可见，在噪声模型最高阶次过小 （为 2）的情况下，辨识结果精度较低且 在辨识初期波动较大，随着模型最高阶次 的增加，辨识结果精度增加，且辨识过程 曲线平坦，但增加到一定程度后辨识曲线 无明显区别。但是噪声模型的阶次越高， 辨识过程的运算量就越大，因此我们在实 际仿真时要综合考虑，选取噪声模型阶次 能较好的辨识出系统参数即可，不必一味 追求过高的模型阶次。4 结论 从本文的介绍和仿真结果可以看出， 系统辨识常用的三种最小二乘类辨识方 法中，递推最小二乘法的算法简单，能减 少计算量，减少数据在计算机中占用的内 存，并实时辨识出系统动态特性。但只能 在噪信比较小的白噪声影响下

25、，才能精确 的辨识出系统参数。广义最小二乘法算法 较为复杂，运算量较大，但是能在系统噪 声为有色噪声时，仍能较好的辨识出系统 参数。实际仿真中，如果我们合理的选取 噪声模型的最高阶次，可以达到较好的辨 识效果，且运算量也较小。%移%进附源程序：clearclose all;L=500;x1=1;x2=1;x3=1;x4=0; 位寄存器初值S=1;%方波初值un=zeros(L,1);for i=1:L+4IM=xor(S,x4); 行异或运算，产生逆 M 序列 if IM=0u(i)=-1; %fixed elseu(i)=1; end un(i,1)=u(i);S=not(S);%产生方波M

26、(i)=xor(x1,x4);%进行异或运算，产生 M 序列 x4=x3;x3=x2;x2=x1;x1=M(i);%寄存器移位 end z(2)=0;z(1)=0;%取 z 的前二个初始值为零 for k=3:L+4; z(k)=*z(k-1)*z(k-2)+*u(k-1)+*u(k-2)+nor mrnd(0,1,1,1);%给出理想的辨识输出采样信号 zn(k,1)=z(k);endc0= ; %a1 a2 b1 b2 给出被辨识参数 的初始值p0=10W*eye(4,4);E=;%相对误差E=c=c0,zeros(4,499);%被辨识参数矩阵的初 始值及大小e=zeros(4,500)

27、;%相对误差的初始值及大小for k=3:500; %开始求 K h1=-z(k-1),-z(k-2),u(k-1),u(k-2);x=h1*p0*h1+1; x1=inv(x); % 开始求 K(k) k1=p0*h1*x1;% 求出 K 的值 d1=z(k)-h1*c0; c1=c0+k1*d1;% 求 被 辨识参数 ce1=c1-cO;%求参数当前值与上一次的 值的差值e2=e1./c0 ;%求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向 量加入误差矩阵的最后一列cO=c1;%新获得的参数作为下一次递 推的旧参数c(:,k)=c1;%把辨识参数 c列向量加入 辨识参数矩阵的最后一列p1=p0-k1*k1*h1*p0*h1+1;% 求 出p(k)的值pO=p1;%给下次用if abs(e2)v=E break;%若参数收敛满足要求，终止计算end%小循环结束en d%大循环结