MS_student_m.rar_MS算法_MS辨识_skin5nm_最小二乘法

%--------------------------------------------------------------------------------- % LS_CAR.m Book, Chapter 3 for A(z)y(t)=B(z)u(t) +v(t) , for given A(z) and B(z) %Example 3.2.1, August 15 Sunday, 2004 % sigma=0.1, 1.0, LS_CAR_sigma0110.* k0=17 % st=40, FF=1, sigma=0.10, 1.00 fprintf('\n----------------------------------------------------------------------'); %---------------------------------------------------------------------* % Filename: LS_CAR.m for the ARX models * % A(z)y(t)=B(z)u(t)+v(t) * % The RLS algorithm * % The forgetting factor FF=\lambda * % The noise variance sigma^2, sigma=0.10 and 1.00 * % Feng Ding * % Ryerson University, Toronto, Canada * % January 29, 2009, Thursday 0:30 am * %---------------------------------------------------------------------* clear; format short g M='The RLS algorithm for the ARX model' %---------------------------------------------------------------------- for st=40 %1, 7, 18, 20,22,40 %Set initial state for random variable save st st clear;clf;format short g load st Method = 'CAR-RLS'; FF=0.98; % Forgetting factor sigma=1; % sigma=0.1, 1.0 % variance of noise figure(1); length1=3100; na = 2; nb = 2; dd=1; % dd: delay =1 a=[-1.6, 0.8]; b=[0.412, -0.309]; d=[1]; %a=[1, 0.412, 0.309]; b=[0, 0.6804,0.6303]; par0=[a,b]'; roots([1,a]); k0=16; %k0=16 n=length(par0); PP = eye(n)*1e6; RR = 1; h = ones(n,1)*1e-6; par1 = ones(n,1)*1e-6; % Compute the noise-to-signal ratio-------------------------- a1=[1,a]; b1=[0,b]; sy=f_integral(a1,b1); sv=f_integral(a1,d); delta_ns = sqrt(sv/sy)*100*sigma; ns = [sv,sy,delta_ns] fprintf('$\\sigma_v^2=%6.2f^2$, $\\delta_{\\ns}=%6.2f%s \n',sigma,delta_ns,'\%$'); %----Generate the input-output data----------------------------------------- %rand('state',1); randn('state',1); %u=(rand(length1,1)- 0.5)*sqrt(12); v=randn(length1,1); rand('state',40); % randn('state',1); eta=f_rn(length1); u=eta(:,1); v=eta(:,2); Gz=tf(b1,a1,1); Gn=tf(d,a1,1); y=lsim(Gz,u) + lsim(Gn,v * sigma); %for k = n:length1 % y(k)=-par0(1:na)'*y(k-1:-1:k-na) + par0(na+1:n)'*u(k-1:-1:k-nb)+v(k)*sigma; %end %----RLS---------------------------------------------- jj = 0; j1 = 0; j2 = 0; for k = 20:length1 jj=jj+1; h = [-y(k-1:-1:k-na); u(k-1:-1:k-nb)]; %delay = 1 yy=y(k); [par1,delta,PP,RR] = f_ls(par0,par1,h,yy,FF,'SG',PP,RR); % [par1,delta,PP,RR] = f_ls(par0,par1,h,yy,0.7,'SG',PP,RR); ls(jj,:)=[jj, par1',delta]; if (jj==100)|(jj==200)|(jj==300)|mod(jj,500)==0 j1 = j1+1; ls_100(j1,:)=[jj, par1', delta*100]; end ls_100(j1+1,:)=[ 0, par0', 0]; end fprintf('\nst = %d, Method = %s, ($\\sigma^2=%5.2f^2$, $\\delta_{\\ns}=%6.2f\\%s$)', st, Method, sigma,delta_ns,'%'); fprintf('\n %s \n','$k$ & $a_1$ & $a_2$ & $b_1$ & $b_2$ & $\delta\ (\%)$ \\'); fprintf('%5d &%10.5f &%10.5f &%10.5f &%10.5f &%10.5f \\\\\n',ls_100'); % st=1:40, fprintf('%d %9.5f %9.2f \\\\\n',st,ls_100(9,n+2),sigma); %figure(1) jk=(17:10:2990)'; plot(ls(jk,1), ls(jk,n+2)); end fprintf('=====================================================================\n') if FF==1 z1=[ls(:,1), ls(:,n+2)]; save z1 z1 elseif FF==0.99 load z1 z1=[z1, ls(:,n+2)]; save z1 z1 else %sigma==1.0 load z1 z0=[z1, ls(:,n+2)]; figure(2) k=(k0:10:2990)'; %%plot(z0(jk,1),z0(jk,2),'k',z0(jk,1),z0(jk,3),'b',z0(jk,1),z0(jk,4),'m') plot(z0(k,1),z0(k,2),'k',z0(k,1),z0(k,3),'b',z0(k,1),z0(k,4),'k') axis([0, 3000, 0, 0.33]) text(600,0.058,'{\it\sigma^2} = 1.00^2') text(600,0.13,'{\it\sigma^2} = 0.10^2') line([247,620],[0.024,0.119]) % legend('{\it\sigma} = 0.10','{\it\sigma} = 1.00') end xlabel('\it t'); ylabel('{\it \delta}');