GLS_ARARX.zip

%《系统辨识理论与方法》递推广义最小二乘算法仿真实例 P157 %(除最后生成的表格格式不对外，其他调试通过) clear; format shortG M='The RGLS algorithm for the CARMA model' FF=1; %The Forgetting factor FF=\Lambda=1 sigma=6; % sigma=1,6 PlotLength=3000; length1=PlotLength + 100; na = 2; nb = 2; nc=1; n = na+nb+nc; a=[1 -1.6 0.8]; b=[0 0.412 0.309]; c=[1 0.8]; d=[1]; ac=[1, a(2)+c(2) a(3)+a(2)*c(2) a(3)*c(2)]; %AZ(z)=A(z)C(z)，用于计算噪信比 par0=[a(2:na+1),b(2:nb+1),c(2:nc+1)]'; p0=1e6; P=eye(n)*p0; r=1; for i=na+nb+1:n %噪声协方差矩阵部分的初值直接取单位阵 P(i,i)=1; end par1=ones(n,1)/p0; % The parameter estimation vector %——Compute the noise-to-signal ratio sy=f_integral(a,b); sv=f_integral(ac,d); delta_ns=sqrt(sv/sy)*100*sigma; [sy,sv,delta_ns] %——Generate the input-output data rand('state',1); randn('state',1); %rng(1,'v5uniform'); rng(1,'v5normal'); u=(rand(length1,1)-0.5)*sqrt(12); v = randn(length1,1)*sigma; t0=20; y=ones(t0,1)/p0; e=y; e1=e; %设置输出、中间变量、估计中间变量初值 for t=n:length1 %得到输出真值 e(t)=-par0(na+nb+1:n)'*e(t-1:-1:t-nc) + v(t); y(t)=par0(1:na+nb)'*[-y(t-1:-1:t-na); u(t-1:-1:t-nb)]+e(t); end % ——Compute RGLS estimate jj=0;j1=0; for t=t0:length1 jj=jj+1; varphi=[-y(t-1:-1:t-na); u(t-1:-1:t-nb); -e1(t-1:-1:t-nc)]; L=P*varphi/(FF+varphi'*P*varphi); P=(P-L*varphi'*P)/FF; par1=par1+L*(y(t)-varphi'*par1); e1(t)=y(t)-par1(1:na+nb)'*varphi(1:na+nb); %刷新中间变量估计值 delta=norm(par1-par0)/norm(par0); ls(jj,:)=[jj, par1', delta]; if (jj==100)||(jj==200)||(jj==500)||mod(jj,1000)==0 j1 = j1+1; ls_100(j1,:)=[jj, par1', delta*100]; end if jj==PlotLength break end end ls_100(j1+1,:)=[0, par0', 0]; fprintf('\n%s\n','$t$ & $a_1$ & $a_2$ & $b_1$ & $b_2$ & $\delta\ (\%)\ \ $ \\\hline'); fprintf('%5d &%10.5f &%10.5f &%10.5f &%10.5f &%10.5f &%10.5f\\\\\n',ls_100); figure(1); k=(22:5:PlotLength-1)'; plot(ls(k,1), ls(k,n+2)); if sigma==1 data1=[ls(:,1), ls(:,n+2)]; save data1 data1 else %sigma==6 load data1 z0=[data1, ls(:,n+2)]; figure(2); plot(k,z0(k,2),'k',k,z0(k,3),'b') axis([0, 3000, 0, 1.02]); text(1000,0.31,'{\it\sigma}^2 = 1.00^2') text(250,0.43,'{\it\sigma}^2 = 6.00^2') line([755,1010],[0.027,0.275]) xlabel('\it t');ylabel('{\it\delta}'); end