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

实验二 递推最小二乘估计(RLS)

及模型阶次辨识(F-Test）

1 实验方案设计

1.1 生成输入数据和噪声

用M序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。 生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。 1.2 过程仿真

辨识模型的形式取A(z1)z(k)B(z1)u(k)e(k)，为方便起见，取nanbn 即

A(z)1a1z1a2a2...anznB(z)1b1zb2a...bnz用M序列作为辨识的输入信号。 1.3 递推遗忘因子法

(0)0.0011001数据长度L取534，初值P(0)000012n

0000 10011.4 计算损失函数、噪声标准差

ˆ(k1)]2[z(k)h(k)损失函数J(k)J(k1)

h(k)P(k1)h(k)ˆ噪声标准差J(L)

Ldim

1.6 F-Test 定阶法计算模型阶次

统计量t

J(n)J(n1)L2n2~F(2,L2n2)

J(n1)2t(n,n1)其中，J()为相应阶次下的损失函数值，L为所用的数据长度，n为模型

的估计阶次。若t(n,n1)ta，拒绝H0:nn0，若t(n,n1)ta，接受H0:nn0，其中t为风险水平下的阀值。这时模型的阶次估计值可取n1。 1.6 计算噪信比和性能指标

e2噪信比 2y~2i~ˆ参数估计平方相对偏差1,iii i1i2n参数估计平方根偏差2()ii12ni12n~22()i~ˆ ,iii2 编程说明

M序列中，M序列循环周期取Np24115，时钟节拍t=1Sec，幅度a1，

G(s)采样时间T0特征多项式为F(s)s6s51。白噪声循环周期为21532768。

设为1Sec，a11.5, a20.7, b11, b20.5。

3 源程序清单

3.1 正态分布白噪声生成函数

function v=noise(N) %生成正态分布N(0,sigma) %生成N个[0 1]均匀分布随机数 A=179; x0=11; M=2^15; for k=1:N x2=A*x0;

x1=mod(x2,M); v1=x1/(M+1); v(:,k)=v1; x0=x1; end

aipi=v;

sigma=1; %标准差

for k=1:length(aipi) ksai=0; for i=1:12

temp=mod(i+k,length(aipi))+1; ksai=ksai+aipi(temp); end

v(k)=sigma*(ksai-6); end end

3.2 M序列生成函数

function [Np r M]=createM(n,a) %生成长度为n的M序列,周期为Np，周期数为r x=[1 1 1 1]; %初始化初态 for i=1:n y=x;

x(2:4)=y(1:3);

x(1)=xor(y(1),y(4)); U(i)=2*y(4)-a; end M=U*a;

lenx=length(x); Np=2^lenx-1; r=n/Np; end

3.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); %初始P

e=0.0001; stop=1; %误差要求，循环停止信号 n=N;

Error=zeros(n_max,1); while(stop==1&&n<=n_max)

H=[]; %新的数据向量 for i=1:na

H=[H;-Z(n-i)]; end

for j=1:nb

H=[H;U(n-j)]; end

K=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); %计算下次递推用到的P temp=abs((X-X_past)./X_past); %相对误差 stop=sum(temp)>=e; %判断精度 Error(n)=Z(n)-H'*X; n=n+1; end

Aes=X(1:na)'; Bes=X(na+1:N)';

3.4 主函数

clear %清理工作间变量

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 [Np r u]=createM(L,a); %生成M序列 figure(1); %画第1个图形：u(k)

stem(u(1:ex)),grid; %以径的形式显示出部分输入信号并给图形加上网格 xlabel('k') %标注横轴变量

ylabel('输入信号') %标注纵轴变量

title(['四级移位寄存器生成M序列输入信号(前',int2str(ex),'位)']) %图形标题 z(2)=0;z(1)=0; %取z的前两个初始值为零

y=z;

v=noise(L);% 生成白噪声 lamat=0.1;

for k=3:L; %循环变量从3到L

y(k)=-aa1*z(k-1)-aa2*z(k-2)+bb1*u(k-1)+bb2*u(k-2); z(k)=y(k)+lamat*v(k); %给出辨识输出采样信号 end

ov=fangcha(v); %计算噪声方差 oy=fangcha(y); %计算信号方差 yita=sqrt(oy\\ov); %计算噪信比

%用最小二乘递推算法辨识参数：a,b,c,d

e0=[0.001 0.001 0.001 0.001]'; %被辨识参数的初始值采用直接取方式，取一个充分小的实向量

p0=10^7*eye(4,4); %初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵 E=1e-10; %相对误差E参考值取0.000000005

e=[e0,zeros(4,L-1)]; %被辨识参数矩阵的初始值及大小 eee=zeros(4,L); %相对误差的初始值及大小 n=0; %用于统计递推次数 for k=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; %完成一次递推，统计值加1

if ee==E break; %若参数收敛满足要求，终止计算 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到L

plot(i,a1(i),'r',i,a2(i),'m',i,b1(i),'c',i,b2(i),'g') %画出a，b，c，d的各次辨识结果 grid on;

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到L

plot(i,ae1(i),'r',i,ae2(i),'g',i,be1(i),'b',i,be2(i),'r:') %画出a,b,c,d的各次辨识结果的收敛情况 grid on;

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; else T_arfa=0.1; end

norder_max=10; %最大阶数

J=zeros(1,norder_max); %各阶残差方差 norder=d; %起始阶数 stop=0; %终止信号 N=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; %残差方差 if norder>1

t=((J(norder-1)-J(norder))/J(norder))*((L-2*norder)/2);

if t<=T_arfa N1=norder-1; N2=N1-d; end end

norder=norder+1; end

disp(['The estimated order by F_test is ',num2str(N1),',The η is ',num2str(yita)]);

4 曲线打印

图1 驱动序列：M序列

图2 递推函数估计值图像

图3 参数误差收敛图像

5 结果分析

根据输出The estimated order by F_test is 2,The η is 0.20874 可以得出该模型估计阶数为2，噪信比为0.20874。

6 实验体会

通过本次试验，我不仅学到了噪声生成和系统辨识的方法，也复习了Matlab编程的方法。在这次试验中，无论是专业课知识还是动手能力上都得到了更高层次提升。