完整word版版系统辨识最小二乘法大作业.docx
- 文档编号:11594818
- 上传时间:2023-03-19
- 格式:DOCX
- 页数:18
- 大小:276.36KB
完整word版版系统辨识最小二乘法大作业.docx
《完整word版版系统辨识最小二乘法大作业.docx》由会员分享,可在线阅读,更多相关《完整word版版系统辨识最小二乘法大作业.docx(18页珍藏版)》请在冰豆网上搜索。
完整word版版系统辨识最小二乘法大作业
西北工业大学
系统辩识大作业
题目:
最小二乘法系统辨识
一、问题重述:
用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数
离散化有
z^4-3.935z^3+5.806z^2-3.807z+0.9362
----------------------------------------------=
z^4-3.808z^3+5.434z^2-3.445z+0.8187
噪声的成形滤波器
离散化有
4.004e-010z^3+4.232e-009z^2+4.066e-009z+3.551e-010
-----------------------------------------------------------------------------=
z^4-3.808z^3+5.434z^2-3.445z+0.8187
采样时间0.01s
要求:
1.用Matlab写出程序代码;
2.画出实际模型和辨识得到模型的误差曲线;
3.画出递推算法迭代时各辨识参数的变化曲线;
最小二乘法:
在系统辨识领域中,最小二乘法是一种得到广泛应用的估计方法,可用于动态,静态,线性,非线性系统。
在使用最小二乘法进行参数估计时,为了实现实时控制,必须优化成参数递推算法,即最小二乘递推算法。
这种辨识方法主要用于在线辨识。
MATLAB是一套高性能数字计算和可视化软件,它集成概念设计,算法开发,建模仿真,实时实现于一体,构成了一个使用方便、界面友好的用户环境,其强大的扩展功能为各领域的应用提供了基础。
对于一个简单的系统,可以通过分析其过程的运动规律,应用一些已知的定理和原理,建立数学模型,即所谓的“白箱建模”。
但对于比较复杂的生产过程,该建模方法有很大的局限性。
由于过程的输入输出信号一般总是可以测量的,而且过程的动态特性必然表现在这些输入输出数据中,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。
这种建模方法就称为系统辨识。
把辨识建模称作“黑箱建模”。
系统辨识又分为参数辨识和阶次辨识,在本文中只讨论参数辨识问题
最小二乘递推算法所用的模型:
Z(k)=B(
)u(k)+v(k)
最小二乘递推算法为:
是服从N
分布的不相关随机噪声。
,
,
考虑如图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)+
(3.69)
仿真对象选择如下的模型结构
(3.68)
其中,
是服从正态分布的白噪声N
。
输入信号采用4位移位寄存器产生的M序列,幅度为1。
按式(3.69)构造h(k);加权阵取单位阵
;利用式(3.61)计算K(k)、
和P(k),计算各次参数辨识的相对误差,精度满足要求式(3.67)后停机。
最小二乘递推算法辨识的Malab6.0程序流程图:
最小二乘递推算法辨识程序
clear%清理工作间变量
L=35;%M序列的周期
y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值
fori=1:
L;%开始循环,长度为L
x1=xor(y3,y4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出
y(i)=y4;%取出第四个移位寄存器的幅值为"0"和"1"的输出信号,即M序列
ify(i)>0.5,u(i)=-1;%如果M序列的值为"1",辨识的输入信号取“-1”
elseu(i)=1;%如果M序列的值为"0",辨识的输入信号取“1”
end%小循环结束
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备
end%大循环结束,产生输入信号u
figure
(1);%第一个图形
stem(u),gridon%显示出输入信号径线图并给图形加上网格
z
(2)=0;z
(1)=0;z(3)=0;z(4)=0;%设z的前四个初始值为零
fork=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.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4);%输出采样信号
end
%RLS递推最小二乘辨识
c0=[0.0010.0010.0010.0010.0010.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(9,9);%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=0.000000005;%取相对误差E=0.000000005
c=[c0,zeros(9,24)];%被辨识参数矩阵的初始值及大小
e=zeros(9,25);%相对误差的初始值及大小
fork=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;%给下次用
ife2<=Ebreak;%如果参数收敛情况满足要求,终止计算
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.935b3=5.806b4=-3.807b5=0.9362。
此时,参数的相对变化量E
0.000000005。
从整个辨识过程来看,精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小,但开始阶段的相对误差会很大,从相对误差图形中可看出,参数的最大相对误差会达到三位数
增广最小二乘法算法所用的模型:
Z(k)=B(
)u(k)+D(
)e(k)
增广最小二乘法算法为:
模型结构选用如下形式:
增广最小二乘算法流程图如下页图所示
增广最小二乘辨识程序
clear
L=55;%4位移位寄存器产生的M序列的周期
y1=1;y2=1;y3=1;y4=0;%4个移位寄存器的输出初始值
fori=1:
L;
x1=xor(y3,y4);%第一个移位寄存器的输入信号
x2=y1;%第二个移位寄存器的输入信号
x3=y2;%第三个移位寄存器的输入信号
x4=y3;%第四个移位寄存器的输入信号
y(i)=y4;%第四个移位寄存器的输出信号,M序列,幅值"0"和"1",
ify(i)>0.5,u(i)=-1;%M序列的值为"1"时,辨识的输入信号取“-1”
elseu(i)=1;%M序列的值为"0"时,辨识的输入信号取“1”
end
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号准备
end
figure
(1);%画第一个图形
subplot(2,1,1);%画第一个图形的第一个子图
stem(u),gridon%画出M序列输入信号
v=randn(1,60);%产生一组60个正态分布的随机噪声
subplot(2,1,2);%画第一个图形的第二个子图
plot(v),gridon;%画出随机噪声信号
u,v%显示输入信号和噪声信号
z=zeros(13,55);%输出采样矩阵
z
(2)=0;z
(1)=0;z(3)=0;z(4)=0;
%输出采样的初值
c0=[0.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye(13,13);%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=0.00000000005;%相对误差E=0.000000005
c=[c0,zeros(13,54)];%被辨识参数矩阵的初始值及大小
e=zeros(13,55);%相对误差的初始值及大小
fork=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.806*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];%findp(k)
p0=p1;%给下次用
ife2<=Ebreak;%若收敛情况满足要求,终止计算
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、d3d4
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、d3d4的收敛情况
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('误差曲线')%标题
参数变化曲线
相对误差曲线
仿真结果表明,递推到第十步时,参数辨识的结果基本达到稳定状态,即a1=-3.808,a2=5.434,a3=-3.445,a4=-0.8187,b1=1,b2=-3.935b3=5.806b4=-3.807b5=0.9362,d1=0.0000,d2=0.0000,d3=0.0000,d4=0.0000。
此时,辨识参数的相对变化量E
0.000000005。
与最小二乘递推算法相比,增广最小二乘递推算法虽然考虑了噪声模型,同样具有速度快、辨识结果精确的特点。
广义最小二乘法所用的模型是:
Z(k)=B(
)u(k)+
e(k)
其中我们选定模型参数:
a1=1.5,a2=-0.7,b1=1,b2=0.5,c1=0,c2=0
广义最小二乘递推算法的计算步骤:
1.给定初始条件
2.利用式
,计算
和
;
3.利用式
,构造
;
4.利用式
递推计算
;
5.利用
和
计算
;
6.根据
来构造
;
7.利用
广义最小二乘法程序代码
clear%清理工作间变量
L=55;%M序列的周期
y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值
fori=1:
L;%开始循环,长度为L
x1=xor(y3,y4);%第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”
x2=y1;%第二个移位寄存器的输入是第一个移位寄存器的输出
x3=y2;%第三个移位寄存器的输入是第二个移位寄存器的输出
x4=y3;%第四个移位寄存器的输入是第三个移位寄存器的输出
y(i)=y4;%取出第四个移位寄存器的幅值为"0"和"1"的输出信号,即M序列
ify(i)>0.5,u(i)=-1;%如果M序列的值为"1",辨识的输入信号取“-1”
elseu(i)=1;%如果M序列的值为"0",辨识的输入信号取“1”
end%小循环结束
y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备
end%大循环结束,产生输入信号u
z
(2)=0;z
(1)=0;%设z的前四个初始值为零
fork=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.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
pf0=10^6*eye(4,4);%直接给出初始状态P0,即一个充分大的实数单位矩阵
ce0=[0.0010.001]';
pe0=eye(2,2);
c=[c0,zeros(4,44)];%被辨识参数矩阵的初始值及大小
ce=[ce0,zeros(2,44)];
e=zeros(4,45);%相对误差的初始值及大小
ee=zeros(2,45);
fork=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。
这种对象利用最小二乘法已可获得无偏一致估计,但广义最小二乘法同时又能给出噪声模型的参数估计。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 完整 word 系统 辨识 最小二乘法 作业