2003版系统辨识最小二乘法大作业
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
西北工业大学系统辩识大作业
题目:最小二乘法系统辨识
一、 问题重述:
用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数
离散化有
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 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。对
4324326.51411.5320120232320
Y s s s s G U s s s s ++++==
++++432
120120232320
E N W s s s s ==
++++
于一个简单的系统
,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。但对于比较复杂的生产过程 ,该建模方法有很大的局限性。由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。系统辨识又分为参数辨识和阶次辨识 ,在本文中只讨论参数辨识问题
最小二乘递推算法所用的模型: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)
+ e (k ) 图1 最小二乘递推算法辨识实例结构图
y (k )
u (k )
z (k )
v (k )
)(1-z N
)(1-z G
其中,)(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 列向量加入辨识参数矩阵的最后一列