(完整word版)2003版系统辨识最小二乘法大作业

西北工业大学系统辩识大作业
题目：最小二乘法系统辨识
一、 问题重述：
用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数
离散化有
z^4 - 3.935 z^3 + 5.806 z^2 - 3.807 z + 0.9362
---------------------------------------------- =
z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187
噪声的成形滤波器
离散化有
4.004e-010 z^3 + 4.232e-009 z^2 + 4.066e-009 z + 3.551e-010
-----------------------------------------------------------------------------
=
z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187
采样时间0.01s
要求：1.用Matlab 写出程序代码；
2.画出实际模型和辨识得到模型的误差曲线；
3.画出递推算法迭代时各辨识参数的变化曲线；
最小二乘法:
在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态 ,静态 , 线性 ,非线性系统。

在使用最小二乘法进行参数估计时 ,为了实现实时控制 ,必须优化成参数递推算法 ,即最小二乘递推算法。

这种辨识方法主要用于在线辨识。

MATLAB 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。

对
于一个简单的系统 ,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。

但对于比较复杂的生产过程 ,该建模方法有很大的局限性。

由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。

这种建模方法就称为系统辨识。

把辨识建模称作“黑箱建模”。

系统辨识又分为参数辨识和阶次辨识 ,在本文中只讨论参数辨识问题
最小二乘递推算法所用的模型：Z(k)=B(
)u(k)+v(k)
最小二乘递推算法为
:
)(k v 是服从N )1,0(分布的不相关随机噪声。

)(1-z G )()(11--=z A z B ，)(1-z N )
()(11
--=z C z D ，
考虑如图 1示仿真对象 ,系统的差分方程为
z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+)(k v (3.69) 仿真对象选择如下的模型结构
)
()4(5)3()2(3)1()()4(4)3(3)2()1()(42121k v k u b k u b k u b k u b k u b k z a k z a k z a k z a k z +-+-+-+-+=-+-+-+-+
(3.68)
其中，)(k v 是服从正态分布的白噪声N )1,0(。

输入信号采用4位移位寄存器产生的M 序列，
幅度为1。

按式(3.69)构造h (k )；加权阵取单位阵I Λ L ；利用式(3.61)计算K (k )、)(ˆk θ
和P (k )，计算各次参数辨识的相对误差，精度满足要求式(3.67)后停机。

最小二乘递推算法辨识的Malab6.0程序流程图：
最小二乘递推算法辨识程序
clear %清理工作间变量
L=35; % M序列的周期
y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值
for i=1:L;%开始循环，长度为L
x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出
y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列
if y(i)>0.5,u(i)=-1; %如果M序列的值为"1", 辨识的输入信号取“-1”
else u(i)=1; %如果M序列的值为"0", 辨识的输入信号取“1”
end %小循环结束
y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备
end %大循环结束，产生输入信号u
figure(1); %第一个图形
stem(u),grid on %显示出输入信号径线图并给图形加上网格
z(2)=0;z(1)=0; z(3)=0;z(4)=0;%设z的前四个初始值为零
for k=5:25; %循环变量从5到25
z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.8 06*u(k-2)-3.807*u(k-3)+0.9362*u(k-4); %输出采样信号
end
%RLS递推最小二乘辨识
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(9,9); %直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.000000005; %取相对误差E=0.000000005
c=[c0,zeros(9,24)]; %被辨识参数矩阵的初始值及大小
e=zeros(9,25); %相对误差的初始值及大小
for k=5:25; %开始求K
h1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4)]';
x=h1'*p0*h1+1; x1=inv(x); %开始求K(k)
k1=p0*h1*x1;%求出K的值
d1=z(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; %给下次用
if e2<=E break; %如果参数收敛情况满足要求，终止计算
end %小循环结束
end %大循环结束
c,e%显示被辨识参数及其误差(收敛)情况
%分离参数
a1=c(1,:); a2=c(2,:);a3=c(3,:); a4=c(4,:);b1=c(5,:); b2=c(6,:); b3=c(7,:); b4=c(8,:); b5=c(9,:);
ea1=e(1,:); ea2=e(2,:);ea3=e(3,:); ea4=e(4,:);eb1=e(5,:); eb2=e(6,:); eb3=e(7,:); eb4=e(8,:); eb5=e(9,:);
figure(2); %第二个图形
i=1:25; %横坐标从1到25
plot(i,a1,'r',i,a2,':',i,a3,'r',i,a4,'k',i,b1,'y',i,b2,':',i,b3,'m',i,b4,':',i, b5,'g') %画出a1，a2，a3,a4,b1，b2,b3,b4,b5的各次辨识结果
title('参数变化曲线') %图形标题
figure(3); %第三个图形
i=1:25; %横坐标从1到25
plot(i,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'c:',i,eb3,'r',i, eb4,':',i,eb5,'g') %画出a1，a2，a3,a4,b1，b2,b3,b4,b5的各次辨识结果的收敛情况title('误差曲线') %图形标题
参数变化图
相对误差图
仿真结果表明，大约递推到第十五步时，参数辨识的结果基本达到稳定状态，即a1=-3.808, a2=5.434, a3=-3.445, a4= -0.8187, b1=1, b2=-3.935 b3=5.806 b4=-3.807 b5=0.9362。

此时，参数的相对变化量E 0.000000005。

从整个辨识过程来看，精度的要求直接影响辨识的速度。

虽然最终的精度可以达到很小，但开始阶段的相对误差会很大，从相对误差图形中可看出，参数的最大相对误差会达到三位数
增
广
最
小
二
乘
法
算
法
所
用
的
模
型
：
Z(k)=B()u(k)+D()e(k)
增广最小二乘法算法为：
模型结构选用如下形式：
)4)-v(k *4d 3)-v(k *3d 2)-v(k *2d 1)-v(k *1d )4(5)3()2(3)1()()4(4)3(3)2()1()(42121++++-+-+-+-+=-+-+-+-+k u b k u b k u b k u b k u b k z a k z a k z a k z a k z
增广最小二乘算法流程图如下页图所示
增广最小二乘辨识程序
clear
L=55; %4位移位寄存器产生的M序列的周期
y1=1;y2=1;y3=1;y4=0; %4个移位寄存器的输出初始值
for i=1:L;
x1=xor(y3,y4); %第一个移位寄存器的输入信号
x2=y1; %第二个移位寄存器的输入信号
x3=y2; %第三个移位寄存器的输入信号
x4=y3; %第四个移位寄存器的输入信号
y(i)=y4; %第四个移位寄存器的输出信号，M序列, 幅值"0"和"1"，
if y(i)>0.5,u(i)=-1; %M序列的值为"1"时,辨识的输入信号取“-1”
else u(i)=1; %M序列的值为"0"时,辨识的输入信号取“1”
end
y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号准备
end
figure(1); %画第一个图形
subplot(2,1,1); %画第一个图形的第一个子图
stem(u),grid on %画出M序列输入信号
v=randn(1,60); %产生一组60个正态分布的随机噪声
subplot(2,1,2); %画第一个图形的第二个子图
plot(v),grid on; %画出随机噪声信号
u ,v %显示输入信号和噪声信号
z=zeros(13,55); %输出采样矩阵
z(2)=0; z(1)=0;z(3)=0; z(4)=0;
%输出采样的初值
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(13,13); %直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.00000000005; %相对误差E=0.000000005
c=[c0,zeros(13,54)]; %被辨识参数矩阵的初始值及大小
e=zeros(13,55); %相对误差的初始值及大小
for k=5:55; %开始求K
z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.8 06*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+4.004e-010*v(k-1)+4.232e-009*v(k-2)+4.066e -009*v(k-3)+3.551e-010*v(k-4); %系统在M序列输入下的输出采样信号
h1=[-z(k-1),-z(k-2),-z(k-3),-z(k-4),u(k),u(k-1),u(k-2),u(k-3),u(k-4),v(k-1),v(k -2),v(k-3),v(k-4)]'; %为求K(k)作准备
x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %K
d1=z(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]; %find p(k)
p0=p1; %给下次用
if e2<=E break; %若收敛情况满足要求，终止计算
end %判断结束
end %循环结束
c, e, %显示被辨识参数及参数收敛情况
%分离变量
a1=c(1,:); a2=c(2,:);a3=c(3,:); a4=c(4,:);b1=c(5,:); b2=c(6,:); b3=c(7,:); b4=c(8,:); b5=c(9,:); %分离出a1，a2，a3,a4,b1，b2,b3,b4,b5
d1=c(10,:); d2=c(11,:); d3=c(12,:); d4=c(13,:);%分离出d1、 d2、 d3 d4
ea1=e(1,:); ea2=e(2,:);ea3=e(3,:); ea4=e(4,:);eb1=e(5,:); eb2=e(6,:); eb3=e(7,:); eb4=e(8,:); eb5=e(9,:); %分离出a1，a2，a3,a4,b1，b2,b3,b4,b5的收敛情况
ed1=e(10,:); ed2=e(11,:); ed3=e(12,:);ed4=e(13,:); %分离出d1 、d2 、d3 d4的收敛情况
figure(2); %画第二个图形
i=1:55;plot(i,a1,'r',i,a2,':',i,a3,'r',i,a4,'k',i,b1,'y',i,b2,':',i,b3,'m',i,b4 ,':',i,b5,'g',i,d1,'k',i,d2,'g:',i,d3,'m',i,d4,'r') %画出各个被辨识参数
title('参数变化曲线') %标题
figure(3); i=1:55; %画出第三个图形
plot(i,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'c:',i,eb3,'r',i, eb4,':',i,eb5,'g',i,ed1,'y',i,ed2,'g:',i,ed3,'k',i,ed4,'m') %画出各个参数收敛情况
title('误差曲线') %标题
参数变化曲线
相对误差曲线
仿真结果表明，递推到第十步时，参数辨识的结果基本达到稳定状态，即a 1=-3.808, a 2=5.434,
a 3=-3.445, a 4= -0.8187,
b 1=1, b 2=-3.935 b 3=5.806 b 4=-3.807 b 5=0.9362，d 1=0.0000, d 2=0.0000, d 3=0.0000，d4=0.0000。

此时，辨识参数的相对变化量E ≤0.000000005。

与最
小二乘递推算法相比，增广最小二乘递推算法虽然考虑了噪声模型，同样具有速度快、辨识结果精确的特点。

广义最小二乘法所用的模型是：
Z(k)=B(
)u(k)+e(k)
⎩⎨
⎧=-+-++-+-=-+-+)
()2()1()()
()2()1()2()1()(2
1
2
1
2
1
k k e c k e c k e k e k u b k u b k z a k z a k z υ 其中我们选定模型参数：a1=1.5,a2=-0.7,b1=1,b2=0.5,c1=0,c2=0
广义最小二乘递推算法的计算步骤：
1.
给定初始条件 ⎪⎪⎩
⎪⎪
⎨⎧====I )0(P )0(ˆ)
(I )0(P )()0(ˆe
e 2
f 0θθ为充分大的数充分小的实向量a a ε
2.
利用式)
()(z )()()(z )(1
1
k z C k u k z C k z f
f
--==，计算)(k z f 和)(k u f
；
3.
利用式⎩⎨
⎧------==τ
τ
)]
(,),1(),(,),1([)(],,,,,[f
f
f
f
f
1
1
b
a
n n n k u k u n k z k z k b b a a b
a
h θ，构造)(f k h ；
4.
利用式[]
⎪⎪⎩
⎪
⎪⎨⎧--=+--=--+-=-)1()]()([)(1)()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆ1
k k k k k k k k k k k k k z k k k f
f f f f
f f f f f f
f
f
P h K P h P h h P K h K τ
τ
τ
θθθI 递推计算)(ˆk θ
； 5. 利用)(ˆ)()()(ˆk k k z k e
θ
τ
h -=和 τ
)](,),1(),(,),1([)(b
a
n k u k u n k z k z k ------= h 计算)(ˆk e
； 6. 根据τ
)](ˆ,),1(ˆ[)(c
e n k e
k e
k ----= h 来构造)(k e h ；
7.
利用[]⎪⎩
⎪
⎨⎧--=+--=--+-=-)1()]()([)(1)()1()()()1()()]1(ˆ)()(ˆ)[()1(ˆ)(ˆ1
k k k k k k k k k k k k k e k k k e
e e e e
e e e e e e
e
e
e
e
P h K P h P h h P K h K τ
τ
τ
θθθI
广义最小二乘法程序代码
clear %清理工作间变量
L=55; % M序列的周期
y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值
for i=1:L;%开始循环，长度为L
x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出
y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列
if y(i)>0.5,u(i)=-1; %如果M序列的值为"1", 辨识的输入信号取“-1”
else u(i)=1; %如果M序列的值为"0", 辨识的输入信号取“1”
end %小循环结束
y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备
end %大循环结束，产生输入信号u
z(2)=0;z(1)=0; %设z的前四个初始值为零
for k=3:45; %循环变量从5到45
z(k)=1.5*z(k-1)-0.7*z(k-2)+1*u(k-1)+0.5*u(k-2); %输出采样信号
end
%RGLS广义最小二乘辨识
c0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量pf0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵
ce0=[0.001 0.001]';
pe0=eye(2,2);
c=[c0,zeros(4,44)]; %被辨识参数矩阵的初始值及大小
ce=[ce0,zeros(2,44)];
e=zeros(4,45); %相对误差的初始值及大小
ee=zeros(2,45);
for k=3:45; %开始求K
zf(k)=z(k)+0*z(k-1)+0*z(k-2);
uf(k)=u(k)+0*u(k-1)+0*u(k-2);
hf1=[-zf(k-1),-zf(k-2),uf(k-1),uf(k-2)]';
x=hf1'*pf0*hf1+1; x1=inv(x); %开始求K(k)
k1=pf0*hf1*x1;%求出K的值
d1=zf(k)-hf1'*c0; c1=c0+k1*d1; %求被辨识参数c
e1=c1-c0; %求参数当前值与上一次的值的差值
e2=e1./c0; %求参数的相对变化
e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1; %新获得的参数作为下一次递推的旧参数
c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列
pf1=pf0-k1*hf1'*pf0; %求出 p(k)的值
pf0=pf1; %给下次用
h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';
ee(k)=z(k)-h1'*c1;
he1=[-ee(k-1),-ee(k-2)]';
x=he1'*pe0*he1+1; x1=inv(x);
k1=pe0*he1*x1;
d1=ee(k)-he1'*ce0;
ce1=ce0+k1*d1;
pe1=pe0-k1*he1'*pe0;
ce0=ce1;
ce(:,k)=ce1;
pe0=pe1;
end %大循环结束
%显示被辨识参数及其误差(收敛)情况
%分离参数
a1=c(1,:); a2=c(2,:);b1=c(3,:);b2=c(4,:); c1=ce(1,:);c2=ce(2,:);
ea1=e(1,:); ea2=e(2,:);eb1=e(3,:); eb2=e(4,:);
figure(2); %第二个图形
i=1:45; %横坐标从1到25
plot(i,a1,'r',i,a2,':',i,b1,'b',i,b2,':',i,c1,'b',i,c2,':') %画出a1，a2,b1，b2,c1,c2的各次辨识结果
title('参数变化曲线') %图形标题
figure(3); %第三个图形
i=1:45; %横坐标从1到25
plot(i,ea1,'r',i,ea2,'k:',i,eb1,'g',i,eb2,'c:') %画出a1，a2，b1，b2,的各次辨识结果的收敛情况
title('误差曲线') %图形标题
结果分析：广义最小二乘法的基本思想是基于对数据先进行一次滤波预处理，然后利用普通的最小二乘法对滤波后的数据进行辨识。

辨识结果a1=1.5016,a2=-0.7010,b1=1.0120,b2=0.4986,c1=0.0040,c2=0.0027与最小二乘法递推算法相比较，可以发现它们的结果基本上是一致的。

这是因为本例的仿真对象相当于
=1。

这种对象利用最小二乘法已可获得无偏一致估计，但广义最小二乘法同时又能
给出噪声模型的参数估计。