系统辨识实验答案相关分析法

实验一 辩识离散线性系统脉冲响应特性的相关分析法

实验步骤及问题

1.首先采用批量算法，步骤如下： （1）首先产生M序列。

通过Matlab软件编程产生M序列，程序如下。 x1=1;x2=1;x3=1;x4=1; m=15;

for i=1:m

y4=x4;y3=x3;y2=x2;y1=x1; x4=y3;x3=y2;x2=y1;x1=xor(y3,y4); if y4==0

u(i)=-1; else u(i)=y4; end end m=u %grapher i1=i; k=1:1:i1;

subplot(3,1,1)

plot(k,u,k,u,'rx') xlabel('k')

ylabel('m序列')

title('移位寄存器产生的m序列') 产生的M序列如下图所示：

（2）系统输出

y(k) + a1y(k-1) + a2y(k-2) = b1u(k-1-d) + b2u(k-2-d)

式中参数值为a1=-0.9,a2=0.5,b1=1.1,b2=0.5,d=0时产生的图形如下:

（3）采用建议的系统参数a和b，观察冲激响应曲线的真值，并估计系统的整定时间Ts；改变系统参数a和b，查看其对结果的影响。 利用批量算法求脉冲响应。

脉冲响应估计值为：1.1928 2.2057 2.6252 1.9906 1.2154 0.9030 0.6983 0.6772 0.8593 1.2405 1.4121 1.5508 1.5512 1.4251 1.3490 图形如下图所示：

当改变a, b参数值时，使a1??1,a2?0.5,b1?1,b2?0.5。此时的脉冲响应估计值为：1.2900 2.2734 2.7686 2.3477 1.6035 1.1689 0.8643 0.7402 0.8232 1.2314 1.4990 1.7188 1.7344 1.6143 1.5176 脉冲响应图形如下图所示：

当改变参数值后，从图形中可以看出脉冲响应的图形形状没有大的改变，但是脉冲响应估计

值发生变化。

（4）采用建议的系统参数a、b和噪声幅度，选取不同的n,k0值，得到不同的M序列周期值，观察结果是否收敛于冲激响应曲线的真值，并分析原因，总结Np的取值规律，判断它与Ts的取值有何关系？

当取n=5, k0=3时产生的图形如下：

从此图可以看出，当n和k0的值改变时，由于m序列是周期性的，系统输出与脉冲响应估计值也具有一定的周期性，结果最后收敛于冲激响应曲线的真值。Np的取值与n有关满足

Np =2n-1 。为了使所得冲激响应能够在Tp=NpΔ之内结束，防止Ru(τ)曲线上出现重叠现象，因此要求Tp=NpΔ>Ts 。

（5）观察M序列的递推算法辩识结果，如果总计算点数不正好等于Np的整数倍，冲激响

应曲线的辩识结果会怎样？

利用递推算法求脉冲响应 图形如下所示：

观察运用递推算法所得到的脉冲响应估计图形，所得到的脉冲响应估计值为：0 0 0.4000 1.2480 2.4027 3.3792 3.4101 2.2806 0.4633 0.9540 -1.5634 -1.3414 -0.8367 -0.3349 0.2585。可以看出与运用批量算法所得到的脉冲响应估计值不同，但是两种算法的辨识结果基本上是一致的，只是递推算法可以用于在线辨识，就这一点来说，递推算法优于批量算法。

如果总计算点数不正好等于Np的整数倍，由于不是在一个完整的周期取值，所得的图形形状以及得到的脉冲响应估计值数据，都会影响系统辨识结果。

（6）观察辩识得到的冲激响应曲线的光滑度，分析导致曲线不光滑的原因？

由于M序列的取值不是连续的 输出也不是连续的，因此曲线不是光滑的。

（7） 观察L序列的递推算法和批量算法辩识结果，同M序列的递推算法和批量算法辩识结果相比较，对两种序列的优缺点作出简要分析。

M序列含有直流成份，这将造成对辨识对象的“净扰动”，这通常是不希望的，而L序列可以很好的克服这一缺点，是一种比M序列更理想的伪随机序列。

1.批量算法程序 x1=1;x2=1;x3=1;x4=1; m=31; for i=1:m

y4=x4;y3=x3;y2=x2;y1=x1; x4=y3;x3=y2;x2=y1;x1=xor(y3,y4); if y4==0 u(i)=-1; else u(i)=y4; end end m=u %grapher i1=i; k=1:1:i1;

subplot(3,1,1) plot(k,u,k,u,'rx')

xlabel('k')

ylabel('m序列')

title('移位寄存器产生的m序列')

a1=-0.9;a2=0.5;b1=1.1;b2=0.5;Y(1)=0;Y(2)=0; for K=3:1:32

Y(K)=0.9*Y(K-1)-0.5*Y(K-2)+1.1*u(K-1)+0.5*u(K-2) end

subplot(3,1,2) k=1:1:i1+1; plot(k,Y,k,Y,'rx') xlabel('k') ylabel('输出') title('系统输出') H=1/32;

p=ones(31);q=eye(31); R=p+q for j=1:31

for i=1:32-j M(i,i+j-1)=u(j) end end

for t=1:31 Z(t,1)=Y(1,t)

从图4和图5所得图形与a1,a2,b1,b2的真值相比较可知，当遗忘因子较小时对参数的跟踪效果较好。

（3）在有噪声，系统参数稳定时（即模型不时变），采用不同的辩识阶次n和时滞延时d，用RLS法求取系统参数和拟合度J，观察系统的辨识结果，然后列表、作图，估计系统的阶次n’和时延r’，判断它们与真值n=2，d=1是否符合？ 由于当n变得比真正的阶次n0大时J几乎停止减小。因此，确定阶次的方法就

?,从下图可以看出当是取最后一次明显地陡降后的n值作为合适的阶次估计值nn>2时，J几乎停止减小，所以可以估计系统阶次n’=2。

（4）参数估计中为什么要采用对输入、输出实测值U’和Y’去除均值后得到的动态分量u’、y’而不用U’和Y’本身？

因为对输入、输出实测值U’和Y’去除均值后得到的动态分量u’、y’可以增加估计参数的精度。

(5) 程序中在进行参数估计之前安排若干步生成u’、y’的计算？如果省去会对结果有什么影响？

如果省去这些计算的话，会影响参数估计的精度。

（6) IV法可否不经过RLS法而直接启动？在IV法中若取不同的滤波器系数ρ

和q，估计的精度是否有差别？为什么？

IV估计的递推算法在形式上与RLS算法相似，因其对初值设定敏感，一般宜先以RLS算法启动，待约50步后再切换到IV法。当选用不同的滤波器参数时估计得精度有差别。

(7) 在有噪声和正确设定阶次及延时(n=2，d=1)的情况下，选取相同的遗忘因

子α=0.995，对三种估计方法的精度进行比较，分析三种估计方法产生误差的原因各有哪些？

RLS法：当噪声比较小时，辨识精度较高；当噪声比较大时，收敛点可能

不唯一，参数估计值往往是有偏的；但是运用此方法时，数据要充分多，否则辨识精度明显下降；另外噪声模型阶次不宜过高，也会影响精度，所以误差的原因主要是数据量的大小以及噪声模型阶次。

IV法：辨识效果较好，能适应较大的范围的噪声特性。对初值 P(0) 敏感，

可用LS法起步，否则没有可靠的收敛性；对数据 u(k) 在 z(k) 的直流分量敏感，对z(k)的直流分量不敏感。所以误差的主要原因是初值P(0)以及数据计算过程产生的误差。

ELS法：ELS是极大似然法的一种近似形式，当数据长度较小时，辨识精

度可能优于极大似然法，但数据长度较大时，精度低于极大似然法。 MATLAB程序: clc

Np=15;r=4;

x1=1;x2=1;x3=1;x4=1; m_length=r*Np; a=1;

for i=1:1:m_length y4=x4;y3=x3;y2=x2;y1=x1; x4=y3;x3=y2;x2=y1; x1=xor(y3,y4); if y4==0 M(i)=-a; else M(i)=a; end end figure;

i=1:1:m_length; plot(i,M);

noise=zeros(1,m_length); for i=1:1:m_length

temp=noise+0.5*rands(1,m_length); noise=temp; end

noise=noise/12;

n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5; s_u0=0.2;s_y0=0.2; u0=ones(1,m_length)*s_u0; u=M+u0+noise; figure;

i=1:1:m_length; plot(i,u);

y(1)=0;y(2)=0;y(3)=0; y(1)=s_y0+y(1)+noise(1); y(2)=s_y0+y(2)+noise(2); y(3)=s_y0+y(3)+noise(3); for k=4:m_length

y(k)=b1*u(k-1-d)+b2*u(k-2-d)-a1*y(k-1)-a2*y(k-2); y(k)=s_y0+y(k)+noise(k); end figure;

i=1:1:m_length; plot(i,y); c=0.5;

p=(c^2)*eye(2*n); sita(:,3)=[0,0,0,0]'; alf=0.995; sum_l=0;sum_y=0; for k=1:1:m_length sum_l=sum_l+u(k); sum_y=sum_y+y(k); end

t_l0=sum_l/m_length;t_y0=sum_y/m_length; for k=1:1;m_length t_l(k)=u(k)-t_l0; t_y(k)=y(k)-t_y0; end

for k=1:1:m_length t_l(k)=u(k)-t_l0; t_y(k)=y(k)-t_y0; end figure;

i=1:1:m_length; plot(i,t_y); for k=4:1:m_length

fai=[-t_y(k-1),-t_y(k-2),t_l(k-1-d),t_l(k-2-d)]'; w=p*fai;

dan=1/(alf+fai'*w);

sita(:,k)=sita(:,k-1)+dan*w*(t_y(k)-fai'*sita(:,k-1)); p=(p-dan*w*w')/alf; end figure;

i=3:1:m_length;

plot(i,sita(1,3:m_length),i,sita(2,3:m_length),i,sita(3,3:m_length),i,sita(4,3:m_length))

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