系统辨识报告

一、
题目要求：
（1）一钢球以零初速由电视天线上自由下落，实验测得下落时间和距
度g。

分析解答：
分析可得，我们已知重力加速度乘以下落时间t的平方除以二等于其降落距离，因此，此处使用最小二乘法的基本算法，即可求得g的估计值，同时经画图可在坐标上描出对应的点，通过其曲线对结果进行相应的验证。

（1）Matlab程序如下：
t=[1 2 3 4 5 6]';
T=0.5*t.^2;
l=[8.94 20.05 50.65 72.19 129.85 171.56]';
H=[ones(6,1),T];
Z=inv(H'*H)*H'*l;
l0=Z(1),g=Z(2)
L=l0+g*T;
plot(T,L,'r',T,l,'b+');hold on
xlabel('t/s');ylabel('l/m');
title('下落时间和距离关系');
hold off
结果：a =
9.7729
（2）做出曲线如下图所示：
图一：小球下落时间与距离关系曲线图
（2）
题目要求：
设两国军备竞赛模型为
x(k)=ax(k-1)+by(k-1)+f
y(k)=cx(k-1)+dy(k-1)+g
式中x(k)和y(k)为两国军事费用（单位：百万美元），已知数据如下：
试用最小二乘法确定模型参数a,b,c,d,f和g.
分析解答：
仍然可直接使用最小二乘法进行估计计算.
1.伊朗和伊拉克：
Matlab程序如下：
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
guji=inv(h'*h)*h'*z; %
算出a、b、f、c、d、g估计值，为三行两列的矩阵
a=guji(1,1),b=guji(2,1),f=g uji(3,1),c=guji(1,2),d=guji
(2,2),g=guji(3,2)%分别取到a、b、f、c、d、g
plot([1:8],x,'r');hold on%画出实际与各自的估计曲线
hg=h*guji;
abf=(hg(:,[1]))';
plot([1:8],abf,'b');hold on
plot([1:8],y,'r');hold on hg=h*guji;
cdg=(hg(:,[2]))';
plot([1:8],cdg,'g');hold on xlabel('k（1-8代表1972到1980）'); ylabel('军事费用/百万美元');%写出坐标表头
legend('伊朗实际值','伊朗
估计值','伊拉克实际值','伊拉克
估计值');
title('伊朗与伊拉克军事费用对比'); hold off
a =
0.5481
b =
-0.2034
f =
4.4475e+003
c =
-0.0306
d =
0.7443
g =
987.1565
：
图二：伊朗与伊拉克军事费用比较图2、北约和华约：
Matlab程序如下：
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
guji=inv(h'*h)*h'*z;%用最小二乘法求出估计值
a=guji(1,1),b=guji(2,1),f=guji(3,1),c=guj i(1,2),d=guji(2,2),g=guji(3,2)
plot([1:10],x,'r');hold on
hg=h*guji;
abf=(hg(:,[1]))';%取guji矩阵的第一列即a，b ,f的值
plot([1:10],abf,'b');hold on
plot([1:10],y,'g');hold on
cdg=(hg(:,[2]))';
plot([1:10],cdg,'b');hold on xlabel('k（1-8代表1972到1980）'); ylabel('军事费用/百万美元');%写出坐标表头
legend('北约实际值','北约估计值','华约克实际值','华约克估计值');
title('北约与华约事费用对比
'); hold off
a =
0.1540
b =
1.7438 f =
-2.5468e+004 c = 0.0181 d =
0.9675 g =
2.0756e+003
图三：华约与北约军事费用对比图
据之间虽然存在区别，但是 二：
题目要求：
考虑理想数学模型为
式中，)(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 +-+-=-+-++-+-() 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 +-+-=-+-+--+-
分析解答：
分析：应该使用增广递推最小二乘算法
n=60;
%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
figure(1);subplot(2,1,1); stem(u),grid on
title('M序列') % M序列%白噪声的产生
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=zeros(7,n);zm d=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);theta1(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('辨识误差')
运行程序，得到如下曲线：
产生的M序列信号
递推增广最小二乘辨识方法
辨识误差
三：题目要求：
完善试验2，给出由相关分析法给出的脉冲响应函数辨识曲线，在试验2的基础上，给出拟合二阶及三阶系统的最小二乘辨识参数曲线，及由脉冲响应函数反求传函()
G z。

1、由实验二已经求得，相关分析法得到的脉冲响应函数辨识曲线
相关分析法所得脉冲响应函数辨识曲线
2、二阶脉冲响应辨识：
Matlab程序：
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)
得到的传递函数为：
Transfer function:
0.3882 z + 0.2645
-----------------------
z^2 - 0.9641 z + 0.2137
Sampling time: 2
得到辨识曲线如图所示：
改为三阶系统的传递函数
运行结果：
0.3888 z^2 + 0.2581 z + 0.002418
G(z)=-------------------------------------
z^3 - 0.9819 z^2 + 0.2452 z - 0.01913
三阶系统的脉冲响应函数
四：比较最小二乘法及其改进算法的性能、优缺点。

答：从最小二乘的定义中得知，最小二乘的思想就是寻找一个θ的估计值θˆ，使
得各次测量的),1(m i Z i =与由估计θˆ确定的量测估计θˆˆi i H Z =之差的平方
和最小，由于此方法兼顾了所有方程的近似程度，使整体误差达到最小，因
而对抑制误差是有利的。

但是。

最小二乘一般方法的估计精度不够高，这是
由于对各个测量数据同等对待，而各次测量数据一般不会在相同的条件下获
得，造成测量数据的置信度不变较大。

因此，最小二乘法有很多改进算法，虽然没有一个是完美的，但是能够适应不同的情况，条件，对应选择不同的
算法，其各自的性能及优缺点如下：
1) 基本的最小二乘法（LS ，RLS ）
对低噪声水平最有效，参数估计可很快收敛到真值，所需计算量相对较少。

对于有色噪声的情况只能得到有偏估计。

几乎不需要特殊的先验假设，在递
推算法中只要求给定 和 ，RLS 具有可靠的收敛性。

2）广义最小二乘法（GLS）
一般可给出好的结果，但计算量较大，信噪比较小时收敛不到真值。

可能出现局部最小值的情况。

3）多级最小二乘法（MSLS）
分三步获得系统模型参数和噪声模型参数的估计，计算量小，速度快，但精度难提高，工程应用不便。

4）增广最小二乘法（RELS）
算法简单，可同时得到参数和噪声模型的估计，工程应用效果很好。

5）加权最小二乘法可对不同置信度的测量值采用加权的办法分别对待，置信度高的，权重取得大些；置信度低的，权重取的小些。

但加权最小二乘法仅能用于事先能估计方程误差对参数估计的影响。

6）递推最小二乘法的可以对计算量大、存储大的数据进行批处理，在线辨识，可以减少数据存储量，避免矩阵求逆，从而减少计算量。

矩形窗（限定记忆）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 当0300
k
≤<
a(k)＝0.6、b(k)＝0.3 当300
k≥
v(k)为零均值白噪声，输入()
u k采用四级逆M序列，系统用数字计算机模拟，利用渐消记忆法估计时变系统参数：
()[(),()
k a k b k
θ=
分别取0.9、0.95、0.99，试说明取值不同对参数估计的影响。

分析解答如下：
分析题意可知，此题应该使用带遗忘因子的递推最小二乘法
L=100;
%设置仿真长度
x1=1,x2=1,x3=1,x4=0;
%四位移位寄存器产生M序列s=1;
for k=1:L
im=xor(s,x4); if(im==0)
u(k)=1; else
u(k)=1; end
s=not(s);
m(k)=xor(x3,x4);
x4=x3;x3=x2;x2=x1;x1=m(k); end %M 序列产生结束 stairs(u);grid; v=randn(L); num=v; k1=k; y(1)=0; y=zeros(1,L); for k=2:L %输出响应 if k<300 y(k)=-0.8*y(k-1)+0.5*u(k-1)+v(k-1); else y(k)=-0.6*y(k-1)+0.3*u(k-1)+v(k-1); end
end
L=1000;
Z=zeros(3,1);
P=10^6*eye(3);
%设置初值
Q=zeros(3,n);
r=0.98;
for k=1:L
K=P*Q/(r+Q'*P*Q);
Q=[-y((k-1):-1:(k-na))
u((k-1):-1:(k-nb))];
C(:,k)=Z+K*(y(k)-Q*Z);
P=(eye(3)-K*Q)*P/r;
Z=C(:,k);
end
k=1:L
plot(k,C);
hold on
1）遗忘因子 =0.9
2）遗忘因子ρ=0.95
2）遗忘因子ρ=0.99
分析：
由上述辨识曲线可以看出，当遗忘因子增大时，其估计所得参数曲线变化并不明显，但是能够看出估计所得参数的变化则越平滑；λ越小，则跟踪越好，但参数会出现一定的跳跃。

但遗忘因子越大，辨识精度及速度都会降低。

但
是由图1、图2、图3可以得到相应的结论。

六、广义最小二乘法和辅助变量法的区别
★辅助变量法（工具变量法）：
某一个变量与模型中随机解释变量高度相关，但却不与随机误差项相关，那么就可以用此变量与模型中相应回归系数的一个一致估计量，这个变量就称为工具变量，这种估计方法就叫工具变量法。

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

缺点：工具变量法的关键是选择一个有效的工具变量，由于工具变量选择中的困难，工具变量法本身存在两方面不足：
一是由于工具变量不是惟一的，因而工具变量估计量有一定的任意性；
二是由于误差项实际上是不可观测的，因而要寻找严格意义上与误差项无关而与所替代的随机解释变量高度相关的变量事实上是困难的。

★广义最小二乘法：
其基本思想在于对数据先进行一次白化滤波处理，然后利用基本的最小二乘法对滤波后的数据进行辨识。

如滤波模型选的合适，对数据进行了较好的白化处理，则直接基本的最小二乘法就能得到无偏一致估计。

所用的滤波模型实际上是一种动态模型，经过几次迭代调整后，便可对数据进行较好的白化处理。

其收敛速度比较慢，需要经过多次迭代计算，才能得到较准确的参数估计值。

一般情况下，经过多次迭代后，估计值便会收敛到稳态值。

但在某些情况下（如噪声比较低时）存在局部极小值，估计值不一定收敛到准则函数的全局极小值上。

为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。

如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。