东北大学系统建模作业 1(2)

系统建模方法作业

姓名：冯天宇

学号：1001258

班级：双控3班

1.考虑如下系统

y () 1.5(1)0.7(2)(3)0.5(4)()k y k y k u k u k k ξ--+-=-+-+

式中，()k ξ为白噪声。

取初值6(0)10P I =，ˆ(0)0θ=。分别选择M 序列和方差为1的正态分布白噪声作为输入信号()u k ，采用递推最小二乘算法进行参数估计，迭代L=400步停止计算。

要求

i) 给出基本迭代公式；

ii) 画出程序流程框图；

iii) 画出输入输出数据曲线、参数估计曲线、误差曲线；

提示：

产生长度为L 方差为1的正态分布白噪声，相应的MATLAB 命令为 randn （L,1）。

i)迭代公式:

1ˆˆˆ()(1)()[()()(1)]1()(1)()()(1)()()()[()()](1)T T T k k k z k k k k k k k k k k k I k k k θθθ-⎧=-+--⎪⎪⎡⎤=--+⎨⎢⎥Λ⎣⎦⎪

⎪=--⎩

K h K P h h P h P K h P

在计算过程中，增益矩阵K(k)可能产生误差，经过算法的迭代，造成误差传

递和累积，最后将影响辨识算法的精度。为此，可采用下面的计算式，以截断误

差的传递，保证辨识精度。

1()(1)()()()(1)()()T T k k k k k k k k ⎡⎤=---+⎢⎥Λ⎣

⎦P P K K h P h

ii)程序流程框图：

N

y(k)

工作间清零

产生输出采样信号

给被辨识参数θ和P赋初值

计算P(k)

计算被辨识参数的相对变化量

参数收敛满足要求？

停机

计算K(k)

计算θ(k)

第四个移位寄存器的输出取反，并将幅值变为1得到辨识系统的输入信号样本值给M序列的长度L和移位寄存器的输入赋初始值

画出被辨识参数θ的各次递推估计值图形

分离参数

画出被辨识参数θ的相对误差的图形

画出辨识的输入信号曲线图形

图最小二乘递推算法辨识的Malab程序流程框图

Y

iii）程序、输入输出数据曲线、参数估计曲线、误差曲线、结果分析：

（1）当输入信号u（k）为M序列时，程序：

clear %清理工作间变量

L=400; % M序列的周期

m1=1;m2=1;m3=1;m4=0; %四个移位寄存器的输出初始值

for i=1:L; %开始循环，长度为L

x1=xor(m3,m4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“异或”

x2=m1; %第二个移位寄存器的输入是第一个移位寄存器的输出x3=m2; %第三个移位寄存器的输入是第二个移位寄存器的输出x4=m3; %第四个移位寄存器的输入是第三个移位寄存器的输出m(i)=m4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列

if m(i)>0.5,u(i)=-1; %如果M序列的值为"1", 辨识的输入信号取“-1”

else u(i)=1; %如果M序列的值为"0", 辨识的输入信号取“1”

end %小循环结束

m1=x1;m2=x2;m3=x3;m4=x4; %为下一次的输入信号做准备

end %大循环结束，产生输入信号u

y(4)=0;y(3)=0;y(2)=0;y(1)=0; %设y的前四个初始值为零

N

for k=5:400 %循环变量从5到400

y(k)=1.5*y(k-1)-0.7*y(k-2)+u(k-3)+0.5*u(k-4); %输出采样信号

end

figure(1), subplot(2,1,1), stem(u),

title('input')

subplot(2,1,2), stem(y), hold on

k=1:400;

plot(k,y)

title('output')

%RLS递推最小二乘辨识

c0=[0 0 0 0]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵c=[c0,zeros(4,399)]; %被辨识参数矩阵的初始值及大小

e=zeros(4,400); %相对误差的初始值及大小

for k=5:400; %开始求K

h1=[-y(k-1),-y(k-2),u(k-3),u(k-4)]';

x=h1'*p0*h1+1;

x1=inv(x); %开始求K(k)

k1=p0*h1*x1; %求出K的值

d1=y(k)-h1'*c0; c1=c0+k1*d1; %求被辨识参数c

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

e2=e1./c0; %求参数的相对变化

e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列

c0=c1; %新获得的参数作为下一次递推的旧参数

c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*[h1'*p0*h1+1]; %求出p(k)的值

p0=p1; %给下次用

end %循环结束

c,e %显示被辨识参数及其误差(收敛)情况

%分离参数

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:400; %横坐标从1到400

plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果

title('Parameter Identification with Recursive Least Squares Method') %图形标题

figure(3); %第三个图形

i=1:400; %横坐标从1到400

plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题