递推最小二乘估计及模型阶次辨识

实验二递推最小二乘估计(RLS)及模型阶次辨识(F-Test）1实验方案设计1.1生成输入数据和噪声用M序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。1.2过程仿真辨识模型的形式取)()()()()(11kekuzBkzzA，为方便起见，取nnnba即nnnnzbabzbzBzaaazazA...1)(...1)(22112211用M序列作为辨识的输入信号。1.3递推遗忘因子法数据长度L取534，初值1000010000100001)0(001.0)0(P1.4计算损失函数、噪声标准差损失函数)()1()()]1(ˆ)()([)1()(2khkPkhkkhkzkJkJ噪声标准差dim)(ˆLLJ1.6F-Test定阶法计算模型阶次统计量t)22,2(~222)1()1()()1,(nLFnLnJnJnJnnt其中，)(J为相应阶次下的损失函数值，L为所用的数据长度，n为模型的估计阶次。若atnnt)1,(，拒绝00:nnH，若atnnt)1,(，接受00:nnH，其中t为风险水平下的阀值。这时模型的阶次估计值可取1n。1.6计算噪信比和性能指标噪信比22ye参数估计平方相对偏差iiiniiiˆ~,~1221参数估计平方根偏差iiiniiniiˆ~,)()~(21221222编程说明M序列中，M序列循环周期取15124pN，时钟节拍t=1Sec，幅度1a，特征多项式为1)(56sssF。白噪声循环周期为32768215。)(sG采样时间0T设为1Sec，5.0,1,7.0,5.12121bbaa。3源程序清单3.1正态分布白噪声生成函数functionv=noise(N)%生成正态分布N(0,sigma)%生成N个[01]均匀分布随机数A=179;x0=11;M=2^15;fork=1:Nx2=A*x0;x1=mod(x2,M);v1=x1/(M+1);v(:,k)=v1;x0=x1;endaipi=v;sigma=1;%标准差fork=1:length(aipi)ksai=0;fori=1:12temp=mod(i+k,length(aipi))+1;ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-6);endend3.2M序列生成函数function[NprM]=createM(n,a)%生成长度为n的M序列,周期为Np，周期数为rx=[1111];%初始化初态fori=1:ny=x;x(2:4)=y(1:3);x(1)=xor(y(1),y(4));U(i)=2*y(4)-a;endM=U*a;lenx=length(x);Np=2^lenx-1;r=n/Np;end3.3加权最小二乘递推算法函数function[Aes,Bes,Error]=RLS(na,nb,Z,U,f)%Aes、Bes为参数估计值，na、nb为模型阶次，Z、U为输出输入数据，f为加权因子N=na+nb;n_max=length(Z);X=0.001.*ones(N,1);%初始估计值P=10^5.*eye(N);%初始Pe=0.0001;stop=1;%误差要求，循环停止信号n=N;Error=zeros(n_max,1);while(stop==1&&n=n_max)H=[];%新的数据向量fori=1:naH=[H;-Z(n-i)];endforj=1:nbH=[H;U(n-j)];endK=P*H*inv(H'*P*H+f);%计算增益矩阵X_past=X;X=X+K*(Z(n)-H'*X);%计算新的估计值P=P-K*K'*(H'*P*H+f);%计算下次递推用到的Ptemp=abs((X-X_past)./X_past);%相对误差stop=sum(temp)=e;%判断精度Error(n)=Z(n)-H'*X;n=n+1;endAes=X(1:na)';Bes=X(na+1:N)';3.4方差函数unctioncc=fangcha(bb)%UNTITLED6Summaryofthisfunctiongoeshere%Detailedexplanationgoesherea=mean(bb);QP=0;fori=1:1:534s(i)=(bb(i)-a)^2;QP=s(i)+QP;endcc=QP/534;end3.5主函数unction[]=leastsquares()L=534;%M序列的周期,四级移位寄存器生成M序列,作为输入信号u(k)ex=60;%在图像中展示的数据个数a=1;aa1=-1.5;aa2=0.7;bb1=1;bb2=0.5;%提前规定的a,b,c,d[Npru]=createM(L,a);%生成M序列figure(1);%画第1个图形：u(k)stem(u(1:ex)),grid;%以径的形式显示出部分输入信号并给图形加上网格xlabel('k')%标注横轴变量ylabel('输入信号')%标注纵轴变量title(['四级移位寄存器生成M序列输入信号(前',int2str(ex),'位)'])%图形标题axis([160-1.51.5])z(2)=0;z(1)=0;%取z的前两个初始值为零y=z;v=noise(L);%生成白噪声lamat=0.1;fork=3:L;%循环变量从3到Ly(k)=-aa1*z(k-1)-aa2*z(k-2)+bb1*u(k-1)+bb2*u(k-2);z(k)=y(k)+lamat*v(k);%给出辨识输出采样信号endov=fangcha(v);%计算噪声方差oy=fangcha(y);%计算信号方差yita=sqrt(oy\ov);%计算噪信比%用最小二乘递推算法辨识参数：a,b,c,de0=[0.0010.0010.0010.001]';%被辨识参数的初始值采用直接取方式，取一个充分小的实向量p0=10^7*eye(4,4);%初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵E=1e-10;%相对误差E参考值取0.000000005e=[e0,zeros(4,L-1)];%被辨识参数矩阵的初始值及大小eee=zeros(4,L);%相对误差的初始值及大小n=0;%用于统计递推次数fork=3:L;%开始递推运算hk=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';%求h(k)K=p0*hk*inv(hk'*p0*hk+1);%求K(k)e1=e0+K*[z(k)-hk'*e0];%求θ(k)ee=(K*(z(k)-hk'*e0));%求相对误差eee(:,k)=ee;%把当前相对变化的列向量加入误差矩阵的最后一列e0=e1;%新获得的参数作为下一次递推的旧参数e(:,k)=e1;%把当前所辨识参数的c1列向量加入辨识参数矩阵的最后一列pk=p0-K*K'*[hk'*p0*hk+1];%求p(k)值p0=pk;%把当前的p(k)值给下次用n=n+1;%完成一次递推，统计值加1ifee==Ebreak;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束a1=e(1,:);a2=e(2,:);b1=e(3,:);b2=e(4,:);ae1=eee(1,:);ae2=eee(2,:);be1=eee(3,:);be2=eee(4,:);%分离参数figure(2);%画第2个图形i=1:ex;%横坐标从1到Lplot(i,a1(i),'r',i,a2(i),'m',i,b1(i),'c',i,b2(i),'g')%画出a，b，c，d的各次辨识结果gridon;xlabel('k')ylabel('辨识参数')%标注纵轴变量title('最小二乘各次递推参数估计值')%图形标题str1=',理论值为：';legend(['a1',str1,num2str(aa1)],['a2',str1,num2str(aa2)],['b1',str1,num2str(bb1)],['b2',str1,num2str(bb2)]);figure(3);%画第3个图形i=1:ex;%横坐标从1到Lplot(i,ae1(i),'r',i,ae2(i),'g',i,be1(i),'b',i,be2(i),'r:')%画出a,b,c,d的各次辨识结果的收敛情况gridon;xlabel('k')%标注横轴变量ylabel('参数误差')%标注纵轴变量legend('a1','a2','b1','b2');title('参数的误差收敛情况')%图形标题d=1;Us=u(11:L);%舍弃前10个数据Zs=v(11:L);arfa=0.01;%置信度if(arfa==0.05)%由自由度确定百分点T_arfa=5.3;elseif(arfa==0.05)T_arfa=4.6;elseif(arfa==0.025)T_arfa=3.69;elseif(arfa==0.05)T_arfa=2.99;elseif(arfa==0.1)T_arfa=2.3;elseT_arfa=0.1;endnorder_max=10;%最大阶数J=zeros(1,norder_max);%各阶残差方差norder=d;%起始阶数stop=0;%终止信号N1=0;PEs=zeros(norder_max,norder_max);%存储RLS估计结果while(stop==0&&norder=norder_max)[Aes,Bes,Error]=RLS(norder,norder,Zs,Us,0.9);%RLS参数估计theta=[Aes,Bes]';thetalength=length(theta);PEs(1:thetalength,norder)=theta;J(norder)=Error'*Error;%残差方差ifnorder1t=((J(norder-1)-J(norder))/J(norder))*((L-2*norder-2)/2);ift=T_arfaN1=norder-1;N2=N1-d;endendnorder=norder+1;enddisp(['TheestimatedorderbyF_testis',num2str(N1),',Theηis',num2str(yita)])end4曲线打印图1驱动序列：M序列图2递推函数估计值图像图3参数误差收敛图像5结果分析根据输出TheestimatedorderbyF_testis2,Theηis0.20874可以得出该模型估计阶数为2，噪信比为0.20874。6实验体会通过这次试验，又验证了另外一种M序列的产生方法的算法，并且取得了预期试验结果。通过实验，对递推最小二乘法有了进一步了解，对课本的理论的知识有了更深一层的认识。电子13级梁建勋（201305051）4、她们宁可做一时的女王，不愿一世的平庸。5、男人插足叫牛逼，女人插足叫小三。6、你要成佛成仙，我跟你去，你要下十八层地狱，我也跟你去。你要投胎，我不答应！7、忘川之畔，与君长相憩，烂泥之中，与君发相缠。寸心无可表，唯有魂一缕。燃起灵犀一炉，枯骨生出曼陀罗。8、如果还有机会的话，我一定会让你回到我的身边，我不想让你和别人结婚。9、每个人心里都有脆弱的一面，如果放大这种脆弱的话，没人想活。10、我做了一个很伟大的决定，看你这么可怜，又没有朋友，我们做朋友吧！11、我不该只是等待，我应该去寻找。12、哪怕再花上七十年，七百年，我想我肯定会找到他！13、恶鬼：你敢打我！夏冬青：你都要吃我了，我还不能打你啊！14、人活着就会