riv01.m

本人用matlab编写的rels(最小二乘法）

%RIV(2阶)
fidin1=fopen('test1.txt');                               % 打开test1.txt文件              
fidout1=fopen('data1.txt','w');                          % 创建data1.txt文件 
while ~feof(fidin1)                                     % 判断是否为文件末尾                
   tline=fgetl(fidin1);                                 % 从文件读行    
   fprintf(fidout1,'%s\n\n',tline);                     % 把此行数据写入文件data1.txt 
end 
fclose(fidout1); 
u=importdata('data1.txt');      % u
fidin2=fopen('test2.txt');                               % 打开test2.txt文件              
fidout2=fopen('data2.txt','w');                          % 创建data2.txt文件 
while ~feof(fidin2)                                     % 判断是否为文件末尾                
   tline=fgetl(fidin2);                                 % 从文件读行    
   fprintf(fidout2,'%s\n\n',tline);                     % 把此行数据写入文件data.txt 
end 
fclose(fidout2); 
y=importdata('data2.txt');      % y
%RIV,RLS启动
C0=[0.0001 0.0001 0.0001 0.0001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
P0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵
r=0.99;
C=[C0,zeros(4,299)];%被辨识参数矩阵的初始值及大小
e=zeros(4,300);%相对误差的初始值及大小
for k=3:100; %开始求K 
x=[-y(k-1),-y(k-2),u(k-1),u(k-2)]'; m=r+x'*P0*x; n=inv(m); K=P0*x*n;%求K的值
C1=C0+K*(y(k)-x'*C0);%求被辨识参数C
P1=(P0-K*x'*P0)/r;%求出 P(k)的值
P0=P1;
e1=C1-C0;%求参数当前值与上一次的值的差值
e2=e1./C0;%求参数的相对变化
e(:,k)=e2; 
C0=C1;%新获得的参数作为下一次递推的旧参数
C(:,k)=C1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 
end%RLS启动结束
C1
z=zeros(300,1);
for k=101:300; %开始求K 
x=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';z1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';
m=r+x'*P0*z1; n=inv(m); K=P0*z1*n;%求K的值
C1=C0+K*(y(k)-x'*C0);%求被辨识参数C
P1=(P0-K*x'*P0)/r;%求出 P(k)的值
P0=P1;
z1(k)=z1'*C1;%求z(k)
e1=C1-C0;%求参数当前值与上一次的值的差值
e2=e1./C0;%求参数的相对变化
e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 
C0=C1;%新获得的参数作为下一次递推的旧参数
C(:,k)=C1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列 
end%循环结束
%分离参数
C1
a1=C(1,:); a2=C(2,:) ; b1=C(3,:) ; b2=C(4,:) ; ea1=e(1,:) ; ea2=e(2,:) ; eb1=e(3,:) ; eb2=e(4,:) ; 
subplot(2,1,1)
i=1:300;%横坐标从1到300
plot(i,a1,'r',i,a2,'r:',i,b1,'g',i,b2,'g:') %画出a1，a2，b1，b2的各次辨识结果
title('被辨识参数每次递推估计值的图形')
subplot(2,1,2)
i=1:300; %横坐标从1到300
plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况
title('被辨识参数的相对误差变化的图形')