all.m

来自「用matlab开发的系统辩识原代码。利用最小二乘法实现辩识.」· M 代码 · 共 92 行

z=zeros(10,100);zs=zeros(10,100);zm=zeros(10,100);zmd=zeros(10,100); %输出采样矩阵、不考虑噪声时系统输出矩阵、不考虑噪声时模型输出矩阵、模型输出矩阵的大小
z(3)=0;z(2)=0; z(1)=0; zs(3)=0;zs(2)=0; zs(1)=0;zm(3)=0; zm(2)=0; zm(1)=0;zmd(3)=0; zmd(2)=0; zmd(1)=0; %输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出的初值
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(10,10); %直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.00000000005; %相对误差E=0.000000005
c=[c0,zeros(10,14)]; %被辨识参数矩阵的初始值及大小
e=zeros(10,15); %相对误差的初始值及大小
for k=4:100; %开始求K 
    z(k)=1.5*z(k-1)-0.7*z(k-2)+0.1*z(k-3)+u(k-1)+0.5*u(k-2)-0.1*u(k-3)+v(k)-v(k-1)+0.2*v(k-2)+0.1*v(k-3); %系统在M序列输入下的输出采样信号    
    h1=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3),v(k),v(k-1),v(k-2),v(k-3)]'; %为求K(k)作准备
    x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K
    d1=z(k)-h1'*c0; c1=c0+k1*d1; %辨识参数c 
    zs(k)=1.5*z(k-1)-0.7*z(k-2)+0.1*z(k-3)+u(k-1)+0.5*u(k-2)-0.1*u(k-3); %系统在M序列的输入下的输出响应
    zm(k)=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3)]*[c1(1);c1(2);c1(3);c1(4);c1(5);c1(6)]; %模型在M序列的输入下的输出响应
    zmd(k)=h1'*c1; %模型在M序列的输入下的输出响应
e1=c1-c0; e2=e1./c0; %求参数误差的相对变化
    e(:,k)=e2;
    c0=c1; %给下一次用
    c(:,k)=c1; %把递推出的辨识参数c 的列向量加入辨识参数矩阵 
    p1=p0-k1*h1'*p0; %find p(k)
    p0=p1; %给下次用
    if e2<=E break; %若收敛情况满足要求，终止计算
    end %判断结束
end %循环结束
c, e, %显示被辨识参数及参数收敛情况
z, zs, zm %显示输出采样值、系统实际输出值、模型输出值
%分离变量
a1=c(1,:); a2=c(2,:);             b1=c(3,:);                     b2=c(4,:); %分离出a1、 a2、 b1、 b2
d1=c(5,:);                      d2=c(6,:);                        d3=c(7,:); %分离出d1、 d2、 d3
ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); %分离出a1、 a2、 b1、 b2的收敛情况
ed1=e(5,:); ed2=e(6,:); ed3=e(7,:); %分离出d1 、d2 、d3的收敛情况
figure(2); %画第二个图形
i=1:100;plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+') %画出各个被辨识参数
title('Parameter Identification with Recursive Least Squares Method') %标题
figure(3); i=1:100; %画出第三个图形
plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+') %画出各个参数收敛情况
title('Identification Precision') %标题
figure(4); subplot(4,1,1); %画出第四个图形，第一个子图
i=1:100;plot(i,zs(i),'r') %画出被辨识系统在没有噪声情况下的实际输出响应
subplot(4,1,2); i=1:100;plot(i,z(i),'g') %第二个子图，画出被辨识系统的采样输出响应
subplot(4,1,3); i=1:100;plot(i,zm(i),'b') %第三个子图，画出模型含有噪声的输出响应
subplot(4,1,4); i=1:100;plot(i,zs(i),'b') %第四个子图，画出模型去除噪声后的输出响应





z=zeros(10,10);zs=zeros(10,10);zm=zeros(10,10);zmd=zeros(10,10); %输出采样矩阵、不考虑噪声时系统输出矩阵、不考虑噪声时模型输出矩阵、模型输出矩阵的大小
z(3)=0;z(2)=0; z(1)=0; zs(3)=0;zs(2)=0; zs(1)=0;zm(3)=0; zm(2)=0; zm(1)=0;zmd(3)=0; zmd(2)=0; zmd(1)=0; %输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出的初值
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(10,10); %直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.00000000005; %相对误差E=0.000000005
c=[c0,zeros(10,14)]; %被辨识参数矩阵的初始值及大小
e=zeros(10,15); %相对误差的初始值及大小
for k=4:100; %开始求K 
    z(k)=1.5*z(k-1)-0.7*z(k-2)+0.1*z(k-3)+u(k-1)+0.5*u(k-2)-0.1*u(k-3)+v(k)-v(k-1)+0.2*v(k-2)+0.1*v(k-3); %系统在M序列输入下的输出采样信号    
    h1=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3),v(k),v(k-1),v(k-2),v(k-3)]'; %为求K(k)作准备
    x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K
    d1=z(k)-h1'*c0; c1=c0+k1*d1; %辨识参数c 
    zs(k)=1.5*z(k-1)-0.7*z(k-2)+0.1*z(k-3)+u(k-1)+0.5*u(k-2)-0.1*u(k-3); %系统在M序列的输入下的输出响应
    zm(k)=[-z(k-1),-z(k-2),-z(k-3),u(k-1),u(k-2),u(k-3)]*[c1(1);c1(2);c1(3);c1(4);c1(5);c1(6)]; %模型在M序列的输入下的输出响应
    zmd(k)=h1'*c1; %模型在M序列的输入下的输出响应
e1=c1-c0; e2=e1./c0; %求参数误差的相对变化
    e(:,k)=e2;
    c0=c1; %给下一次用
    c(:,k)=c1; %把递推出的辨识参数c 的列向量加入辨识参数矩阵 
    p1=p0-k1*h1'*p0; %find p(k)
    p0=p1; %给下次用
    if e2<=E break; %若收敛情况满足要求，终止计算
    end %判断结束
end %循环结束
c, e, %显示被辨识参数及参数收敛情况
z, zs, zm %显示输出采样值、系统实际输出值、模型输出值
%分离变量
a1=c(1,:); a2=c(2,:);     a3=c(3,:);      b1=c(4,:);b2=c(5,:);b3=c(6,:);   %分离出a1、 a2、 b1、 b2
d1=c(7,:); d2=c(8,:);d3=c(9,:);d4=c(10,:);      %分离出d1、 d2、 d3
ea1=e(1,:); ea2=e(2,:); ea3=e(3,:); eb1=e(4,:);eb2=e(5,:);eb3=e(6,:); %分离出a1、 a2、 b1、 b2的收敛情况
ed1=e(7,:); ed2=e(8,:); ed3=e(9,:);ed4=e(10,:); %分离出d1 、d2 、d3的收敛情况
figure(2); %画第二个图形
i=1:100;
plot(i,a1,'r',i,a2,'r:',i,a3,'r',i,b1,'b',i,b2,'b:',i,b3,'b',i,d1,'g',i,d2,'g:',i,d3,'g+',i,d4,'g:') %画出各个被辨识参数
title('Parameter Identification with Recursive Least Squares Method') %标题
figure(3); i=1:100; %画出第三个图形
plot(i,ea1,'r',i,ea2,'r:',i,ea3,'r',i,eb1,'b',i,eb2,'b:',i,eb3,'b',i,ed1,'g',i,ed2,'g:',i,ed3,'g+',i,ed4,'g') %画出各个参数收敛情况
title('Identification Precision') %标题
figure(4); subplot(4,1,1); %画出第四个图形，第一个子图
i=1:100;plot(i,zs(i),'r') %画出被辨识系统在没有噪声情况下的实际输出响应
subplot(4,1,2); i=1:100;plot(i,z(i),'g') %第二个子图，画出被辨识系统的采样输出响应
subplot(4,1,3); i=1:100;plot(i,zm(i),'b') %第三个子图，画出模型含有噪声的输出响应
subplot(4,1,4); i=1:100;plot(i,zmd(i),'b') %第四个子图，画出模型去除噪声后的输出响应