系统辨识大作业

系统辨识作业

学院：信息与控制工程学院专业班级：

姓名：

学号：

2011年4月15日

一、（1） 题目要求

分析：

对于初速度为零，从初始位置下落，下落距离等于重力加速度乘以下落时间t 的平方的二分之一，因此，此处使用最小二乘法的基本算法，即可求得g 的估计值，同时经画图可在坐标上描出对应的点，通过其曲线对结果进行相应的验证。 作业程序： clear all close all clc %%%% 第一题（1）程序 %%%% t=[0 1 2 3 4 5 6]'; L=[0 8.94 20.05 50.65 72.19 129.85 171.56]'; T=0.5*t.^2;

A=polyfit(T,L,1); z=polyval(A,T); figure plot(T,z); figure plot(T,L,'b+') g=A(1)%辨识参数 hold on plot(T,z,'r'); xlabel('T/s^2'); ylabel('L/m'); title('最小二乘法求重力加速度g'); hold off 运行结果：

a.最小二乘法拟合曲线图

b ．参数辨识： g =9.5793 即求得重力加速度值为9.5793。

（2）

题目要求

设两国军备竞赛模型为

x(k)=ax(k-1)+by(k-1)+f

y(k)=cx(k-1)+dy(k-1)+g

分析

同上题，仍然利用最小二乘法，只不过此题不是单个数，而是矩阵。 作业程序

●伊朗和伊拉克

clear all

close all

clc

%%%% 第一题（2）1程序 %%%%

y=[909 1123 2210 2247 2204 2303 2179 2675 ];

x=[2891 3982 8801 11230 12178 9867 9165 5080 ];

for k=2:8

h(k,:)=[x(k-1) y(k-1) 1];

z(k,:)=[x(k) y(k)];

end

est=inv(h'*h)*h'*z; % 算出a、b、f、c、d、g估计值，为三行两列的矩阵

a=est(1,1),b=est(2,1),f=est(3,1) ,c=est(1,2),d=est(2,2),g=est(3,2 )%分别取到a、b、f、c、d、g plot([1972:1979],x,'r');hold on%画出实际与各自的估计曲线

hg=h*est;

abf=(hg(:,[1]))';

plot([1972:1979],abf,'b');hold on plot([1972:1979],y,'r');hold on

hg=h*est;

cdg=(hg(:,[2]))';

plot([1972:1979],cdg,'g'); xlabel('年代');

ylabel('军事费用/百万美元');%写出坐标表头

legend('伊朗实际值','伊朗估计值','伊拉克实际值','伊拉克估计值');

title('伊朗与伊拉克军事费用对比'); hold off

北约和华约

clear all

close all

clc

x=[216478 211146 212267 210525 205717 212009 215988 218561 255411 233957 ];

y=[112893 115020 117169 119612 121461 123561 125498 127185 129000 131595];

for k=2:10

h(k,:)=[x(k-1) y(k-1) 1];

z(k,:)=[x(k) y(k)];

end

est=inv(h'*h)*h'*z;%用最小二乘法求出估计值

a=est(1,1),b=est(2,1),f=est(3,1)

,c=est(1,2),d=est(2,2),g=est(3,2

)

plot([1972:1981],x,'r');hold on

hg=h*est;

abf=(hg(:,[1]))';%取guji矩阵的第一列即a，b ,f的值

plot([1972:1981],abf,'b');hold on plot([1972:1981],y,'g');hold on

cdg=(hg(:,[2]))';

plot([1972:1981],cdg,'b');hold on xlabel('年代'); ylabel('军事费用/百万美元');%写出坐标表头

legend('北约实际值','北约估计值','华约克实际值','华约克估计值');

title('北约与华约事费用对比'); hold off

最小二乘法辨识结果

结果分析

由以上对比曲线可以看出，经过最小二乘法估计得到的数据与实际数据之间虽然存在区别，但是基本符合要求, 故可取该组辨识参数数据。

二、 题目要求

考虑理想数学模型为

选择如下的辨识模型进行增广递推最小二乘参数辨识。

() 1.5(1)0.7(2)(1)0.5(2) 1.2()(1)0.2(2)z k z k z k u k u k v k v k v k +-+-=-+-+--+-式中，)(k V 是服从正态分布的白噪声)1,0(N 。输入信号采用4阶M 序列，其幅值为1。

1212123()(1)(2)(1)(2)()(1)(2)

z k a z k a z k bu k b u k c v k c v k c v k +-+-=-+-++-+-

给出各参数的辨识曲线和辨识误差曲线。 分析

应该使用增广递推最小二乘算法。

作业程序

clear all

close all

clc

%%%% 第二题程序 %%%%

N=15;

%4位移位寄存器产生的M序列的周期

y1=1;y2=1;y3=1;y4=0;

for i=1:N;

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

%白噪声的产生

A=19;x0=12;M=1024;

for k=1:N

x=A*x0;

x1=mod(x,M);

v(k)=x1/512;

x0=x1;

end

%辨识主程序

z=zeros(7,N);zs=zeros(7,N);zm=ze ros(7,N);zmd=zeros(7,N);

z(1)=0;z(2)=0;

zs(1)=0; zs(2)=0;

zm(1)=0; zm(2)=0;

zmd(1)=0;zmd(2)=0;

theta=zeros(7,1);

p0=10^6*eye(7,7);

the=[theta,zeros(7,14)];

e=zeros(7,15);

for k=3:N; z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1 )+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k -2);

h=[-z(k-1),-z(k-2),u(k-1),u(k-2) ,v(k),v(k-1),v(k-2)]';

x=inv(h'*p0*h+1);

q=p0*h*x;

d1=z(k)-h'*theta;

theta1=theta+q*d1;

zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k -1)+0.5*u(k-2);

zm(k)=[-z(k-1),-z(k-2),u(k-1),u( k-2)]*[theta1(1);theta1(2);theta 1(3);theta1(4)];

zmd(k)=h'*theta1;

e(:,k)=theta1-theta;

theta=theta1;

the(:,k)=theta1;

p1=p0-q*q'*[h'*p0*h+1];

p0=p1;

end

figure(1);

i=1:N;

plot(i,the(1,:),'r',i,the(2,:),' r:',i,the(3,:),'b',i,the(4,:),'b :',i,the(5,:),'g',i,the(6,:),'g: ',i,the(7,:),'g+')

title('辨识曲线')

figure(2);

i=1:N;

plot(i,e(1,:),'r',i,e(2,:),'r:', i,e(3,:),'b',i,e(4,:),'b:',i,e(5 ,:),'g',i,e(6,:),'g:',i,e(7,:),' r+')

title('辨识误差曲线')

运行程序，得到以下曲线

得到参数辨识结果

三、

题目要求

完善试验2，给出由相关分析法给出的脉冲响应函数辨识曲线，在试验2的基础上，给出拟合二阶及三阶系统的最小二乘辨识参数曲线，及由脉冲响应函数反求传函()

G z。

分析

利用带遗忘因子的递推算法。

作业程序

clear all

close all

clc

%%%% 第三题程序 %%%%

n=2;Ts=2;np=15; %n阶次，Ts采样时间

U=UY(61:300,1); %输入初值U

Y=UY(61:300,2); %输入初值Y

f=zeros(120,4);

f(:,1)=-1*Y(61:180);

f(:,2)=-1*Y(60:179);

f(:,3)=U(61:180);

f(:,4)=U(60:179);

yy=Y(62:181);

q=inv(f'*f)*f'*yy;

a=q(1:2);

b=q(3:4); G=tf(b',[1 a'],2) %传递函数GLS(z)

yc=zeros(n+np,1); %脉冲响应

uc=zeros(n+np,1);

uc(n+1)=1/Ts;

for i=1:np-1

yc(n+i)=-1*[yc(i+n-1:-1:i)]'*a+[ uc(i+n-1:-1:i)]'*b;

end

yimp=yc(n+1:n+np);

figure(1);

stem(0:Ts:(np-1)*Ts,2*yimp,'r') hold on

plot(0:Ts:(np-1)*Ts,2*yimp,'r^-' )

impulse(G)

程序运行结果：

（n=2）

Transfer function:

0.3933 z + 0.2639

-------------------------------

z^2 - 0.9544 z + 0.2004

Sampling time: 2

得到辨识曲线如图所示：

（n=3）

Transfer function:

0.3954 z^2 + 0.323 z + 0.05386

-------------------------------------------------- z^3 - 0.8075 z^2 + 0.08785 z + 0.004068

得到的传递函数为： Transfer function: 0.3882 z + 0.2645 ----------------------- z^2 - 0.9641 z + 0.2137 三阶系统的脉冲响应函数

四、 题目要求

查资料，比较最小二乘法及其改进算法的性能，优缺点。

解答

最小二乘法：最小二乘的思想就是寻找一个θ的估计值θ?，使得各次测量的

),1(m i Z i =与由估计θ?确定的量测估计θ??i i H Z =之差的平方和最小，由于此方

法兼顾了所有方程的近似程度，使整体误差达到最小，因而对抑制误差是有利的。但是最小二乘一般方法的估计精度不够高，这是由于对各个测量数据同等对待，而各次测量数据一般不会在相同的条件下获得，造成测量数据的置信度不变较大，当新数据源源而来时，将出现以下问题：1.数据增加，要求计算机的存储空间增加2、每增加一组数据，即作一次求逆，导致计算量增加，难以用于在线辨识。只有当模型的噪声项是独立的随机变量时，普通最小二乘法才能得到真实参数的无偏估计值，否则所得到的估计值是有偏的。

因此，最小二乘法有很多改进算法，虽然没有一个是完美的，但是能够适应不同的情况，条件，对应选择不同的算法，其各自的性能及优缺点如下：

广义最小二乘法的优点是计算精度高，估计的效果比较好，是无偏估计，但缺点是计算量大，其收敛是比较缓慢的，为了得到准确的参数估值，往往需要进行多次迭代计算，另外，对于循环的循环性没有给出证明，并非总是收敛于最优估值上。由于一般情况下，系统信噪比比较低，准则函数为非单值函数（即存在多个局部极小值），如果初值给的不合理，用GLS 方法得到的将是局部极小值，若想得到总体最优解，初值应接近该最优值。

递推算法实现了实时控制，减少了计算量和存储量，但未解决最小二乘法的有偏估计问题。矩形窗（限定记忆）RLS 方法需要保留一定的数据存储量，此存储量大小取决于矩形窗宽度，因而在应用范围上有一定程度的限制。

辅助变量法的计算较简单，估计是无偏估计，但计算精度较低； 增广矩阵法能保证精度和收敛，算法简单，可同时得到参数和噪声模型的估计，工程应用效果很好但计算量也较大。

多级最小二乘法把复杂的辨识问题分成三个阶段来处理，分三步获得系统模型参数和噪声模型参数的估计计算量小，速度快，而且每个阶段只用到简单的最小二乘法，省去了广义最小二乘法的迭代过程，简化了计算，并且可以得到参数的一致性无偏估计，计算精度比辅助变量法高。但是这种方法也存在着与广义最小二乘法相类似的收敛问题，而且精度难提高，工程应用不便。

加权最小二乘法可对不同置信度的测量值采用加权的办法分别对待，置信度高的，权重取得大些；置信度低的，权重取的小些。但加权最小二乘法仅能用于事先能估计方程误差 对参数估计的影响。

快速多阶段最小二乘法大大减少计算时间。

五、 题目要求

时变系统的参数估计问题。

已知系统 y(k)+a(k)y(k-1)=b(k)u(k-1)+v(k) 其中a(k)、 b(k)具有以下数值：

a(k)＝0.8、 b(k)＝0.5 当0300k ≤< a(k)＝0.6、 b(k)＝0.3 当300k ≥

v(k)为零均值白噪声，输入()u k 采用四级逆M 序列，系统用数字计算机模拟，利用渐消记忆法估计时变系统参数：

()[(),()k a k b k θ=

分别取0.9、0.95、0.99，试说明ρ取值不同对参数估计的影响。

作业程序

clear all close all clc

%%%% 第五题程序 %%%%

a=[1 0.8]'; b=[0.5]'; d=1; %对象参数

na=length(a)-1;

nb=length(b)-1; %na 、nb 为A 、B 阶次 L=2000; %仿真长度

uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)

yk=zeros(na,1); %输出初值 x1=1; x2=1; x3=1;

x4=0;%移位寄存器初值 S=1;%方波初值 for k=1:L

IM=xor(S,x4);%进行异或运算，产生逆M 序列

if IM==0 u(k)=-1;

else

u(k)=1; end

S=not(S);%产生方波

M(k)=xor(x3,x4);%进行异或运算，产生M 序列

x4=x3;x3=x2;x2=x1;x1=M(k);%寄存器移位 end

xi=randn(L,1);

thetae_1=zeros(na+nb+1,1); %thet ae 初值

P=10^6*eye(na+nb+1);

lambda=0.99; %遗忘因子范围[0.9 1] for k=1:L

if k==301

a=[1 0.6]';b=[0.3]'; %对象参数突变 end

theta(:,k)=[a(2:na+1);b]; %对象参数真值

phi=[-yk;uk(d:d+nb)];

y(k)=phi'*theta(:,k)+xi(k); %采集输出数据

%遗忘因子递推最小二乘法

K=P*phi/(lambda+phi'*P*phi) thetae(:,k)=thetae_1+K*(y(k)-phi '*thetae_1);

P=(eye(na+nb+1)-K*phi')*P/lambda ;

%更新数据

thetae_1=thetae(:,k); for i=d+nb:-1:2 uk(i)=uk(i-1); end

uk(1)=u(k); for i=na:-1:2 yk(i)=yk(i-1); end

yk(1)=y(k); end

subplot(1,2,1)

plot([1:L],thetae(1:na,:)); hold on;

plot([1:L],theta(1:na,:),'k:'); xlabel('k'); ylabel('参数估计a'); legend('a'); axis([0 L -2 2]); subplot(1,2,2)

plot([1:L],thetae(na+1:na+nb+1,:));

hold on;

plot([1:L],theta(na+1:na+nb+1,:) ,'k:'); xlabel('k'); ylabel('参数估计b'); legend('b'); axis([0 L -0.5 2]);

取不同ρ,得到不同的参数估计1)遗忘因子ρ=0.9

2)遗忘因子ρ=0.95

3)遗忘因子 =0.99

分析

由上述辨识曲线可以看出，遗忘因子在0.9到1之间变化时，其对应的参数估计曲线随着遗忘因子的增大而趋于稳定，波动减小，但相比而言，0.9处

最快达到稳定值。

六、

题目要求

广义最小二乘法和辅助变量法的区别

解答

这道题与第四道题中部分问题类似，重复部分不再赘述。

★广义最小二乘法：

其基本思想在于对数据先进行一次白化滤波处理，然后利用基本的最小二乘法对滤波后的数据进行辨识。如滤波模型选的合适，对数据进行了较好的白化处理，则直接基本的最小二乘法就能得到无偏一致估计。所用的滤波模型实际上是一种动态模型，经过几次迭代调整后，便可对数据进行较好的白化处理。其收敛速度比较慢，需要经过多次迭代计算，才能得到较准确的参数估计值。一般情况下，经过多次迭代后，估计值便会收敛到稳态值。但在某些情况下（如噪声比较低时）存在局部极小值，估计值不一定收敛到准则函数的全局极小值上。为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。其过于复杂。

★辅助变量法（工具变量法）：

广义LS算法过于复杂，因此寻求简便一些的一致无偏估计算法很有实用价值。辅助变量法最小二乘法便是这样的算法。

某一个变量与模型中随机解释变量高度相关，但却不与随机误差项相关，那么就可以用此变量与模型中相应回归系数的一个一致估计量，这个变量就称为工具变量，这种估计方法就叫工具变量法。

辅助变量法的辨识精度叫最小二乘法高，尤其适当造声势有色噪声，且噪声结构不好确定时，辅助变量法更显示出它的优势，因为这种方法无需知道噪声结构能获得系统的无偏一致估计。

缺点：工具变量法的关键是选择一个有效的工具变量，由于工具变量选择中的困难，工具变量法本身存在两方面不足：

一是由于工具变量不是惟一的，因而工具变量估计量有一定的任意性；

二是由于误差项实际上是不可观测的，因而要寻找严格意义上与误差项无关而与所替代的随机解释变量高度相关的变量事实上是困难的。

参考书目：

《系统辨识理论及应用》国防工业出版社

《自适应控制理论及应用》哈尔滨工业大学出版社

《系统辨识与自适应控制Matlab编程》北京航空航天大学出版社