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

最小二乘辨识方法的优劣比较摘 要：本文系统的探讨了三种最小二乘类辨识方法的原理和性能，并对各种方法在各种不同的环境下进行了MATLAB仿真，仿真结果证明：最小二乘法不适合实时处理，在同等情况下，递推最小二乘的辨识速度较快，但在有色噪声干扰下效果不理想，广义最小二乘法的辨识效果最好，且不受噪声是否有色的影响，但是费时最多。关键词：最小二乘 辨识速度 MATLAB仿真 1 引言系统辨识是一门介于现代控制理论和系统理论的边缘学科它将现代控制论的平滑、滤波、预测和参数估计理论，以及系统论的系统分析方法和建模思想应用于自然科学、社会科学和工程实践中的各个领域，与各个领域的专业知识相给合，形成了一个个新的交叉学科分支。关于系统辨识的含义，早在1962年Zacleh曾作如下定义：“根据系统的输入和输出，在指定的一类系统中确定一个相被辨识系统等价的系统”。根据这个定义，在系统辨识中必须确定三方面的问题；第一，必须指定一类系统即根据先验信息确定系统模型的类型。第二，必须规定一类插入信号。例如正弦信号、阶跃信号、脉冲信号、白噪声、伪随机信号等。而且这些信号从时域考虑，必须能持续地激励系统的所有状态；从频域考虑，输入信号的频带能覆盖系统的频带宽度。第三，必须规定“系统等价”的含义及其度量准则。2 线性系统的辨识2.1 问题描述考虑如下线性系统： （1）其中，u(k)为系统激励信号，y(k)为系统输出，e(k)为模型噪声。其系统模型如图1所示：图1 SISO的系统模型结构图其中G(z-1)是系统函数模型，N(z-1)为有色噪声系统模型，e(k)为白噪声v(k)经过系统函数为N(z-1)的系统后的输出。通常 （2）式中： （3） （4）则系统可表示为： （5）设样本和参数集为： （6）h(k)为可观测的量， 差分方程可写为最小二乘形式 （7）如何系统噪声e(k)存在的情况下从该方程中正确的解出，即是系统辨识的任务。为了求出，我们面临三大问题：一是输入信号的选择，二是判决准则的选取，三是辨识算法的选择，下面一一探讨。2.2 选择输入 为了准确辨识系统参数，我们对输入信号有两大要求，一是信号要能持续的激励系统所有状态，二是信号频带能覆盖系统的频带宽度。除此之外还要求信号有可重复性，不能是不可重复的随机噪声，因此我们通常选择M序列或逆M序列作为输入。2.3 准则函数因为本文主要探讨最小二乘类辨识方法，在此选取准则函数 （8）使准则函数的估计值记做，称作参数的最小二乘估计值。在式(7)中，令k=1,2,3,L，可构成线性方程组： （9）式中 （10）准则函数相应变为： （11）极小化,求得参数的估计值，将使模型更好的预报系统的输出。2.4 辨识算法常用的最小二乘类辨识方法有下列三种：最小二乘法，递推最小二乘法和广义最小二乘法。2.4.1 最小二乘法设使得，则有 （12）展开上式，并根据以下两个向量微分公式： （13）得正则方程： （14）当为正则阵时，有 （15）且有，所以满足式（15）的唯一使得，这种通过极小化式（11）计算的方法称作最小二乘法。而且可以证明，当噪声e(k)是均值为0的高斯白噪声时，可实现无偏估计。2.4.2 递推最小二乘算法为了减少计算量，减少数据在计算机中占用的内存，并实时辨识出系统动态特性，我们常利用最小二乘法的递推形式。下面我们来推导递推最小二乘算法的原理。首先，将式(11)的最小二乘一次完成算法写为 （16）定义 （17）式中 （18）式中，h(i)是一个列向量，也就是HL的第i行的倒置，P(k)是一个方阵，他的维数取决于未知参数的个数，假设未知参数的个数是n，则P(k)的维数是nn.。由式17可得P(k)的递推关系为： （19）设 （20）则 （21）由此可得： （22）由式19和22可得 （23）引进增益矩阵K(k)，定义 （24）式（23）可以进一步写为 （25）接下来可以进一步把式(20)写为 （26）利用矩阵反演公式将式（26）演变成 （27）将上式代入式(24)，整理后可得 (28)综合式25、27和28可得最小二乘递推参数估计算法RLS （29）2.4.3 广义最小二乘法设SISO系统采用如下模型： （30）假定模型阶次na，nb和nc已知，用广义最小二乘法可以得到无偏一致估计。令 （31） （32）将模型化为最小二乘格式： （33）由于v(k)是白噪声，所以用最小二乘可以获得参数的无偏估计，由于噪声模型C(z-1)未知，还需要用迭代的方法来求得C(z-1)。令 （34）置 （35）这样就把噪声模型也转变为最小二乘格式： （36）由于上式中的噪声已为白噪声，所以用最小二乘也可获得参数e的无偏估计，但是数据向量中依然含有不可测的噪声量，可用相应的估计值来代替，置，其中k 0时，按照 （37）计算，式中 （38）综上所述，广义最小二乘法可归纳为 （39）3 仿真研究已知系统模型:x(k)-1.5x(k-1)+0.7x(k-2)= u(k-1)+0.5u(k-2), y(k)=x(k)+(k),(k)=(k)， u、x、y、分别为模型输入、模型输出、测量输出、干扰噪声。输入u为逆m序列：信号幅值a=1、寄存器位数为n=5(信号长度N=2n-1=31)，重复周期数q=40。为噪信比调整因子，噪信比定义为：NSR=,x、分别为模型输出x和噪声的均方差（标准差），有两种模型：(1)为白噪声，(2)为有色噪声，噪声模型为：(k)=e(k)+0.5e(k-1)+0.9(k-1)-0.95(k-2),e(k)为白噪声。定义辨识误差值：，其中：为独立的实验次数,为模型真值，为估计值。3.1 产生输入和输出数据选择自相关特性好的逆m序列作为输入。利用MATLAB产生寄存器位数n=5，每周期长为31，重复周期数q=40的逆m序列，并将其作为输入得到系统输出。绘出一个周期的输入输出图形分别如图2和图3所示。 图2 寄存器位数为5的逆m序列 图3 逆m序列经过系统的输出3.2 产生系统噪声为了后面能较好的区分每种辨识方法的性能，我们分别在输出中叠加白噪声和有色噪声。取NSR=20%，用同一噪声源产生两种噪声模型，分别绘制白噪声、用相同噪声模型产生的有色噪声和不同噪声影响下的系统输出的曲线，如图4、图5图6和图7所示。 图4 方差为0.2的高斯白噪声 图5 白噪声影响下的系统输出 图6 白噪声经过噪声模型后产生的有色噪声 图7 有色噪声影响下的系统输出3.3 最小二乘辨识模型辨识为较好的研究最小二乘辨识模型的性能，作者分别在不同的噪声模型下，用不同的噪信比影响系统输出，利用输入输出数据对系统进行辨识。分别采用白噪声模型和有色噪声模型，取NSR=0%、5%、10%、15%、20%、25%、30%、35%、40%、45%、50%，每种工况下取独立试验次数N=50（每次独立产生噪声），数据序列取前1024点，用最小二乘法辨识模型，分别画出NSR曲线（图8和图9）。图中的纵坐标（辨识误差）是50次辨识误差的均值。 图8 白噪声对各系数的辨识精度影响 图9 有色噪声对各系数的辨识精度影响由图可见： 在白噪声影响下，各系数的辨识误差都很小，欲辨识参数为a1=-1.5，a2=0.7，b1=2，b2=0.5，即使在噪信比为50%的情况下，四个参数的辨识误差都在10-3数量级，相对误差非常小，均可视为无偏估计，与理论相符。 在有色噪声影响下，各系数的辨识误差相对白噪声影响偏大，当噪信比达到50%时，其中，a1、a2和b1的辨识误差都在0.020.04之间，相对误差在10%左右，b2的辨识误差甚至达到0.16以上，相对误差达到30%以上。综上所述：在只有白噪声影响下，最小二乘辨识法可以达到无偏估计，但是在有色噪声影响下辨识结果的相对误差较大。最小二乘法只适合用于只有白噪声影响下的系统辨识，对于有色噪声影响下的系统，我们应该寻求更好的辨识方法。3.4 递推最小二乘辨识模型辨识分别采用白噪声模型和有色噪声模型，取NSR=10%、40%，用递推最小二乘法辨识模型参数，对比画出各参数辨识结果随递推次数变化的曲线。为了对比研究，我们在同一组u、x序列下，用同一白噪声源产生给定噪信比的白噪声和有色噪声干扰。欲辨识参数为a1=-1.5，a2=0.7，b1=2，b2=0.5，设定在两种辨识情况下，前后两次辨识误差小于0.000000005时，结束仿真，当设定NSR=0.1时，本次仿真循环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所示 图14白噪声影响下参数辨识结果 图15有色噪声影响下参数辨识结果 图16 白噪声影响下的参数辨识误差 图17有色噪声影响下的参数辨识误差（NSR=40%时）从仿真结果图中我们可以看到，当噪信比较小时（如NSR=0.1）时，在白噪声影响下，递推最小二乘能正确辨识系统参数，且辨识曲线较平稳。而在有色噪声影响下，辨识结果有一点误差，但辨识曲线尚较平稳。当噪信比较大（如NSR=0.4）时，不管是在白噪声还是有色噪声影响下，辨识曲线的波动都较大，且辨识误差都比噪信比小时的辨识误差有所增加。综上所述：递推最小二乘法只适合与噪信比较小时的白噪声影响下的系统辨识，对于有色噪声影响下的系统辨识或者噪信比较大时的白噪声影响下的系统辨识，我们应该寻找更好的辨识方法。3.5 广义最小二乘辨识模型3.5.1 广义最小二乘与递推最小二乘的性能比较采用有色噪声模型，取NSR=10%、30%，分别用递推最小二乘和广义最小二乘递推法辨识系统参数，为了对比研究，我们在同一组u、y序列下进行辨识试验，并画出各参数辨识结果随次数变化的曲线。图18和图19是取NSR=10%时分别采用递推最小二乘法和广义最小二乘法的辨识结果，图20和图21分别是其相对辨识误差。 图18 递推最小二乘法的辨识结果 图19广义最小二乘法的辨识结果 图20 递推最小二乘法的辨识精度 图21 广义最小二乘法的辨识精度保持其他参数不变，取NSR=30%，再次运行程序可得图22图25。 图22 递推最小二乘法的辨识结果 图23广义最小二乘法的辨识结果 图24 递推最小二乘法的辨识精度 图25 广义最小二乘法的辨识精度 由图可见：两种方法在较小噪信比的有色噪声影响下均能较好的辨识出系统参数，把相对误差曲线图锁定在（-0.02，0.02）的区间上进行观察可得：广义最小二乘法能更快更平坦的接近辨识结果，辨识效果更好。当有色噪声强度增加到30%时，两种方法的辨识结果均出现了不同程度的失真，尤其是递推最小二乘法，结果失真严重且辨识曲线波动较大，广义最小二乘法的辨识结果虽有失真，但仍比较接近真实结果，且辨识曲线平坦。因此我们可以得出结论：由于引入了对噪声的估计，在噪信比较大的有色噪声影响下，广义最小二乘法仍然能够得到较好的辨识结果。3.5.2 广义最小二乘的辨识性能与噪声模型的最高阶次在上面的广义最小二乘法的辨识过程中，我们取的噪声模型的最高阶次为4，如果改变噪声模型的最高阶次结果会有何改变呢？为此，作者编写程序得到了在噪信比为30%的有色噪声影响下，广义最小二乘法进行辨识时，噪声模型阶次依次从2开始以2为步长递增到8时的辨识结果曲线，如图26图29所示。 图26 噪声模型阶次为2时的辨识结果 图27噪声模型阶次为3时的辨识结果 图28 噪声模型阶次为4时的辨识结果 图29 噪声模型阶次为6时的辨识结果由图可见，在噪声模型最高阶次过小（为2）的情况下，辨识结果精度较低且在辨识初期波动较大，随着模型最高阶次的增加，辨识结果精度增加，且辨识过程曲线平坦，但增加到一定程度后辨识曲线无明显区别。但是噪声模型的阶次越高，辨识过程的运算量就越大，因此我们在实际仿真时要综合考虑，选取噪声模型阶次能较好的辨识出系统参数即可，不必一味追求过高的模型阶次。4 结论从本文的介绍和仿真结果可以看出，系统辨识常用的三种最小二乘类辨识方法中，递推最小二乘法的算法简单，能减少计算量，减少数据在计算机中占用的内存，并实时辨识出系统动态特性。但只能在噪信比较小的白噪声影响下，才能精确的辨识出系统参数。广义最小二乘法算法较为复杂，运算量较大，但是能在系统噪声为有色噪声时，仍能较好的辨识出系统参数。实际仿真中，如果我们合理的选取噪声模型的最高阶次，可以达到较好的辨识效果，且运算量也较小。附源程序：clearclose all;L=500;x1=1;x2=1;x3=1;x4=0; %移位寄存器初值S=1;%方波初值un=zeros(L,1);for i=1:L+4 IM=xor(S,x4); %进行异或运算，产生逆M序列 if IM=0 u(i)=-1; %fixed else u(i)=1; end un(i,1)=u(i); S=not(S);%产生方波 M(i)=xor(x1,x4); %进行异或运算，产生M序列 x4=x3;x3=x2;x2=x1;x1=M(i); %寄存器移位endz(2)=0;z(1)=0; %取z的前二个初始值为零for k=3:L+4; 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.001 0.001 0.001 0.001 ; %a1 a2 b1 b2给出被辨识参数的初始值p0=106*eye(4,4);E=0.000000005;%相对误差E=0.000000005c=c0,zeros(4,499);%被辨识参数矩阵的初始值及大小e=zeros(4,500);%相对误差的初始值及大小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;%求被辨识参数c e1=c1-c0;%求参数当前值与上一次的值的差值 e2=e1./c0 ;%求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1;%新获得的参数作为下一次递推的旧参数 c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1;%求出 p(k)的值 p0=p1;%给下次用 if abs(e2)=E break;%若参数收敛满足要求，终止计算 end%小循环结束end%大循环结束% c %显示被辨识参数% e %显示辨识结果的收敛情况% e2;a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);i=1:500;plot(i,a1,r,i,a2,-,i,b1,g:,i,b2,k-.,LineWidth,2) %画出a1，a2，a3,b1，b2的各次辨识结果title(各参数最小二乘法辨识结果)%图形标题legend(a1,a2,b1,b2,0);grid on;figure;for i=1:L H1(i,1)=z(i); H1(i,2)=u(i); end estimate1=inv(H1*H1)*H1*z(2:501); D1