系统辨识实验报告

序号16

系统辨识理论及Matlab仿真

上 机 实 验 报 告

姓 名： 班 级： 学 号： 指导老师：

2016 年 1 月

一、上机题目

设单输入—单输出系统的差分方程为

z（k）+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k) v(k)=c1v(k)+c2v(k-1)+c3v(k-2)

取真值a1=1.6、a2=0.7、b1=10、b2=0.4、c1=0.9、c2=1.2和c3=0.2,输入信号采用4阶M序列，幅值为1.当v(k)的均值为0，方差分别为0.1和0.5的高斯噪声时，分别用一般最小二乘法、递推最小二乘法和增广递推最小二乘法估计参数θ。并通过对三种方法的辨识结果的分析和比较，说明上述三种参数辨识方法的优缺点。

解：

一般最小二乘法：

Matlab仿真程序：

clear all close all clc

randn('seed',100);

v=0.1*randn(1,60);%产生一组高斯分布的随机噪声 %%%%M序列产生程序%%%%%%% L=60;

y1=1;y2=1;y3=1;y4=0;%移位寄存器初值 for i=1:L; x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4;

if y(i)>0.5,u(i)=-1; else u(i)=1; end

y1=x1;y2=x2;y3=x3;y4=x4; end figure(1); stem(u),grid on

title('输入信号M序列')

%%%%最小二乘辨识程序%%%%%%%%%% z=zeros(1,16);%定义输出观测值的长度 for k=3:16

z(k)=-1.6*z(k-1)-0.7*z(k-2)+u(k-1)+0.4*u(k-2)+1*v(k)+1.6*v(k-1)+0.7*v(k-2);%观测值 end figure(2) plot([1:16],z) title('输出观测值') figure(3) stem(z),grid on

title('输出观测值z的经线图形')

H=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);

-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)];%给出样本系数矩阵

Z=[z(3); z(4); z(5); z(6); z(7); z(8); z(9); z(10); z(11); z(12); z(13); z(14); z(15); z(16);]%给出样本观测矩阵

c=inv(H'*H)*H'*Z;%计算参数

al=c(1),a2=c(2),b1=c(3),b2=c(4)%分离参数

程序运行结果： Z =

-0.7165 -0.2873 -0.3262 -0.5818 1.7445 -1.0598 1.8351 -2.6803 3.7050 -2.8528 0.8821 -1.2964 1.7166 -2.6295 al =1.3684 a2 =0.4863 b1 =0.9853 b2 =0.1820

递推最小二乘法：

Matlab仿真程序：

%%%%递推最小二乘算法程序%%%%%%% z(2)=0;z(1)=0;

%观测值由理想输出值加噪声 for k=3:60;%循环变量从3到15

z(k)=-1.6*z(k-1)-0.7*z(k-2)+u(k-1)+0.4*u(k-2)+1*v(k)+1.6*v(k-1)+0.7*v(k-2); end

%%%%RLS递推最小二乘辨识%%%%%%%% c0=[0.001 0.001 0.001 0.001 ]'; p0=10^3*eye(4,4); E=0.000000005;%相对误差

c=[c0,zeros(4,59)];%被辨识参数矩阵初始值 e=zeros(4,60);%相对误差的初始值 lamt=1; for k=3:60;

h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; k1=p0*h1*inv(h1'*p0*h1+1*lamt);%求出k的值 new=z(k)-h1'*c0;

c1=c0+k1*new;%被辨识参数c p1=1/lamt*(eye(4)-k1*h1')*p0;

e1=(c1-c0)./c0;%求参数当前值与上一次的值的差值

e(:,k)=e1;%把当前相对变化的列向量加入误差矩阵的最后一列 c(:,k)=c1;%把辨识参数c列向量加入便是参数矩阵的最后一列 c0=c1;%新的参数作为下一次递推的旧参数 p0=p1;

if norm(e1)<=E%若参数收敛满足要求，终值计算 break; end end %分离参数

a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:); ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:); figure(2); i=1:60;

plot(i,a1,'k',i,a2,'b',i,b1,'r',i,b2,'g') legend('a1','a2','b1','b2'); title('递推最小参数辨识') figure(3); i=1:60;

plot(i,ea1,'k',i,ea2,'b',i,eb1,'r',i,eb2,'g') legend('a1','a2','b1','b2');

title('辨识精度')

程序运行结果：

增广递推最小二乘法：

Matlab仿真程序：

%%%%%%增广递推最小二乘辨识%%%%%%%%%% z(2)=0;z(1)=0;

theat0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初值 p0=10^4*eye(7,7);%初始状态P0

theat=[theat0,zeros(7,59)];%被辨识参数矩阵的初始值 for k=3:60;

z(k)=-1.6*z(k-1)-0.7*z(k-2)+u(k-1)+0.4*u(k-2)+1*v(k)+1.6*v(k-1)+0.7*v(k-2); h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';

x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; d1=z(k)-h1'*theat0;

theat1=theat0+k1*d1;%辨识参数c theat0=theat1;%给下一次用

theat(:,k)=theat1;%把辨识参数c列向量加入辨识参数矩阵 p1=p0-k1*k1'*[h1'*p0*h1+1]; p0=p1;%给下次用 end

%%%%分离变量%%%%%%%%%%%

a1=theat(1,:);a2=theat(2,:);b1=theat(3,:);b2=theat(4,:);c1=theat(5,:);c2=theat(6,:);c3=theat(7,:); i=1:60; figure(2); plot(i,z) figure(3)

plot(i,a1,'r',i,a2,'b',i,b1,'k',i,b2,'y',i,c1,'g',i,c2,'c',i,c3,'m')%画出各个被辨识参数 title('增广递推最小二乘辨识算法')

程序运行结果：

三种最小二乘辨识方法的优缺点： 一般最小二乘算法：

优点：适用于理论分析，一次完成计算。

缺点：①数据量越多，系统参数估计的精度就越高，为了获得满意的辨识结果，矩阵阶数常常取的很大，从而造成了计算量和需要的存储量也很大，所以多用于离线辨识。

②每增加一次观测量，都必须重新计算。

递推最小二乘算法：

优点：当矩阵维数增加时，能相比于一般最小二乘辨识较少的计算量，为实时地辨识出动态系统的特性。

增广递推最小二乘算法：

优点：当噪声不是白噪声时，一般最小二乘法参数估计不是无偏估计，增广最小二乘法通过先对数据进行一次滤波预处理，进而进行最小二乘法辨识，最终获得最小二乘法的无偏、一致估计。

教 师 评 语 成 绩

年 月 日

本文来源：https://www.bwwdw.com/article/9myt.html