系统辨识第五章作业分析.docx
《系统辨识第五章作业分析.docx》由会员分享,可在线阅读,更多相关《系统辨识第五章作业分析.docx(26页珍藏版)》请在冰点文库上搜索。
![系统辨识第五章作业分析.docx](https://file1.bingdoc.com/fileroot1/2023-7/11/92f90e39-259d-4081-9552-bf6212c141ec/92f90e39-259d-4081-9552-bf6212c141ec1.gif)
系统辨识第五章作业分析
摘要
系统辨识是描述各种各样系统运动规律的一种方法论,是研究系统的一种有效工具。
利用这个工具可以对我们要研究的系统进行定量描述。
随着现代控制理论的迅速发展,系统辨识得到迅速而蓬勃发展,并已经成功运用与多种工程应用领域。
但针对有色噪声干扰系统,传统的辨识方法不能得到良好的参数估计,而工程上大多系统都为有色噪声干扰系统。
有色噪声干扰系统的一类系统为广义输出误差模型,本文对白噪声和有色噪声两种噪声干扰下运用最小二乘法,递推最小二乘法分别比较噪声的不同,以及进一步得出最小二乘法的适用性。
进一步研究广义最小二乘法对有色噪声系统辨识的改进。
关键词:
系统辨识;白噪声;有色噪声;最小二乘;递推最小二乘;广义最小二乘
ABSTRACT
thelawofmotionofthesystemidentificationistodescribethevarioussystemandaneffectivetoolforthesystem.Italsocanbeusedtostudyquantitativedescriptionofthesystem.Withtherapiddevelopmentofmoderncontroltheory,thedevelopmentofsystemidentificationhasbeenrapidandvigorous,andhasbeensuccessfullyappliedtoavarietyofengineeringapplications.However,thetraditionalidentificationmethodscan'tgetgoodparameterestimationinengineeringsystemsforcolorednoise.Colorednoisejammingsystemisforgeneralizedoutputerrormodel.undertwokindsofnoiseofwhitenoiseandcolornoise,byusingtheleastsquaremethodandtherecursiveleastsquaremethod,respectivelytocomparethedifferentnoise,andtofurtherconcludedthattheapplicabilityoftheleast-squaremethod.Tofurtherstudiedtheimprovementofgeneralizedleastsquaresidentificationofcolorednoisesystem.
KEYWORDS:
systemidentification,whitenoise,colorednoise.,leastsquares;,recursiveleastsquares,thegeneralizedleastsquares
目录
摘要1
第一章系统辨识2
1.1辨识的定义2
1.2辨识的内容和步骤3
1.3辨识目的3
1.4先验知识3
1.5数据预处理4
1.6模型结构辨识4
1.7模型参数辨识4
1.8模型检验4
第二章题目详解5
2.1.1最小二乘法5
2.1.2递推法8
2.2.1最小二乘法12
2.2.2递推法13
2.3广义最小二乘法15
附录15
第一章系统辨识
1.1辨识的定义
Zadeh对辨识的定义(1962年)辨识就是在输入和输出数据的基础上,从一组给定的模型类中,确定一个与所测系统等价的模型。
Liung对辨识的的定义(1978年)系统辩识有三个要素——数据、模型类和(等价)准则。
系统辩识是按照一个准则,在模型类中选择一个与数据拟合得最好的模型。
Liung认为,实际系统的复杂性很难找到一个适用的模型与之等价。
因此,系统辩识的任务只是要求从输入输出数据出发,找到一个与实际系统相逼近的模型。
该定义体现了逼近的观点。
系统辨识是在输入和输出观测的基础上,从指定的一类模型中确定一个与系统等价的模型。
1.2辨识的内容和步骤
简单地说,
辨识就是一种从观测到的含有噪声的输入输出数据中提取数学模型的方法。
根据现场的情况,辨识可以离线进行,也可以在线进行;辨识的内容主要包括四个方面:
实验设计,模型结构辨识,模型参数辨识,模型检验。
1.3辨识目的
明确模型应用的最终目的很重要,因为它将决定模型的类型,精度要求和采用哪种辨识方法等问题。
1.4先验知识
对于给定的过程进行辨识之前,要通过一些手段对过程有一些了解粗略的掌握一些过程和数据,这些先验知识对实验设计将其指导作用。
1.5数据预处理
输入输出数据通常含有直流或交流成分,任何方法都无法消除对系统精度的影响。
此外,高频成分也是不利的。
因此需对输入数据进行零均值化和剔除高频成分预处理。
1.6模型结构辨识
模型结构辨识包括模型验前结构和模型结构参数两部分。
1.7模型参数辨识
当模型结构确定之后,就需要进行模型参数辨识。
方法很多,最基本,应用最广的是最小二乘法。
但是最小二乘法也有一些重大缺陷,比如过程辨识受到有色噪声污染时,他几乎比能适应。
1.8模型检验
模型检验是辨识过程不可或缺的步骤之一。
但是,它没有一般方法可循。
它和模型结构问题密切相关。
如果模型结构不合适,模型检验是不能通过的。
第二章题目详解
2.1.1最小二乘法
一个单输入-单输出线性定常系统的差分方程
式中
为输入信号,
为理论上的输出值。
通过观测得到,在观测过程中往往附加有随机干扰。
观测值用
表示:
ξ(k)=v(k)+Σ
则
;
如果u(k)也有测量误差,则在
ξ(k)中应包含这一测量误差。
则
原题设为二阶系统,且
的真值已知,输入
已知,求
,并且由输入输出序列得出参数估值,其程序已在matlab中实现。
实验结果如下图:
参数
a1
a2
b1
b2
真实值
1.642
0.715
0.39
0.35
估计值
1.6008
0.6483
0.4143
0.3309
实验结果分析:
由于所用的输出观测值有白噪声成分,所以辨识结果有误差。
2.1.2递推法
图1.2最小二乘递推算法辨识的Malab6.0程序流程图
Matlab仿真结果:
c=
Columns1through12
1.642001.6420000000000
0.715000.7150000000000
0.390000.3689000000000
0.350000.2298000000000
Columns13through24
000000000000
000000000000
000000000000
000000000000
Columns25through30
000000
000000
000000
000000
e=
Columns1through12
000000000000
000000000000
00-0.0540000000000
00-0.3433000000000
Columns13through24
000000000000
000000000000
000000000000
000000000000
Columns25through30
000000
000000
000000
000000
结果分析:
仿真结果表明,大约递推到第五步时,参数辨识的结果基本达到稳定状态。
此时,参数的相对变化量E
0.000000005。
从整个辨识过程来看,精度的要求直接影响辨识的速度。
虽然最终的精度可以达到很小,但开始阶段的相对误差会很大。
2.2.1最小二乘法
题设将白噪声改为有色噪声,
;有色噪声可以看做是白噪声经过成型滤波器后得到的。
Matlab仿真结果:
参数
a1
a2
b1
b2
真实值
1.642
0.715
0.39
0.35
估计值
1.2875
0.3702
0.3462
0.2773
结果分析:
由于所用的输出观测值有有色噪声成分,所以辨识结果有误差,且用最小二乘法计算得到的参数估计误差较大。
2.2.2递推法
结果:
结果分析:
第一问和第二问的区别只是噪声不同,但是从仿真结果图像来看,都是从第五步开始达到稳定状态的,但有色噪声参数辨识显然误差更大更明显。
2.3广义最小二乘法
广义最小二乘法辨识的计算步骤如下:
(1)应用已得到的输入和输出数据u(k)和y(k)
(2)计算残差e(k),并用其代替
(3)计算和
和
(4)应用得到的
和
,
,再用最小二乘法重新估计
,得
的第2次估值
。
然后按步骤⑵计算残差e(k)。
重新估计
,得到估值
。
再按步骤⑶计算
和
,按步骤⑷求
的第3次估计
。
重复上循环步骤,直到
的估值
收敛为止。
附录
程序
%1.1二阶系统的最小二乘一次完成算法辨识程序,在光盘中的文件名:
FLch3LSeg1.m
u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.474,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177];%系统辨识的输入信号为一个周期的M序列
y=zeros(1,30);%定义输出观测值的长度
w=0.1*randn(1,30)
fork=3:
30
y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+w(k);%用理想输出值作为观测值
end
subplot(3,1,1)%画三行一列图形窗口中的第一个图形
stem(u)%画输入信号u的径线图形
subplot(3,1,2)%画三行一列图形窗口中的第二个图形
i=1:
1:
30;%横坐标范围是1到16,步长为1
plot(i,y)%图形的横坐标是采样时刻i,纵坐标是输出观测值z,图形格式为连续曲线
subplot(3,1,3)%画三行一列图形窗口中的第三个图形
stem(y),gridon%画出输出观测值z的径线图形,并显示坐标网格
u,y%显示输入信号和输出观测信号
%L=30%数据长度
HL=[-y
(2)-y
(1)u
(2)u
(1);-y(3)-y
(2)u(3)u
(2);-y(4)-y(3)u(4)u(3);-y(5)-y(4)u(5)u(4);-y(6)-y(5)u(6)u(5);-y(7)-y(6)u(7)u(6);-y(8)-y(7)u(8)u(7);-y(9)-y(8)u(9)u(8);-y(10)-y(9)u(10)u(9);-y(11)-y(10)u(11)u(10);-y(12)-y(11)u(12)u(11);-y(13)-y(12)u(13)u(12);-y(14)-y(13)u(14)u(13);-y(15)-y(14)u(15)u(14);-y(16)-y(15)u(16)u(15);-y(17)-y(16)u(17)u(16);-y(18)-y(17)u(18)u(17);-y(19)-y(18)u(19)u(18);-y(20)-y(19)u(20)u(19);-y(21)-y(20)u(21)u(20);-y(22)-y(21)u(22)u(21);-y(23)-y(22)u(23)u(22);-y(24)-y(23)u(24)u(23);-y(25)-y(24)u(25)u(24);-y(26)-y(25)u(26)u(25);-y(27)-y(26)u(27)u(26);-y(28)-y(27)u(28)u(27);-y(29)-y(28)u(29)u(28)]%给样本矩阵HL赋值
Y=[y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16);y(17);y(18);y(19);y(20);y(21);y(22);y(23);y(24);y(25);y(26);y(27);y(28);y(29);y(30)]%给样本矩阵zL赋值
%CalculatingParameters
c1=HL'*HL;c2=inv(c1);c3=HL'*Y;c=c2*c3%计算并显示
a1=c
(1),a2=c
(2),b1=c(3),b2=c(4)%从中分离出并显示a1、a2、b1、b2
%1.2RLS递推最小二乘辨识
c0=[1.6420.7150.390.35]';%被辨识参数的初始值
p0=10^6*eye(4,4);%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=0.000000005;%取相对误差E=0.000000005
c=[c0,zeros(4,29)];%被辨识参数矩阵的初始值及大小
e=zeros(4,30);%相对误差的初始值及大小
u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.474,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177];%系统辨识的输入信号为一个周期的M序列
y=zeros(1,30);%定义输出观测值的长度
w=0.5*randn(1,30);
fork=3:
30
y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+w(k);%用理想输出值作为观测值
end
fork=3:
30;%开始求K
h1=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';
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;%给下次用
ife2<=E
break;%如果参数收敛情况满足要求,终止计算
end%小循环结束
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:
30;%横坐标从1到30
plot(i,a1,'r',i,a2,':
',i,b1,'g',i,b2,':
')%画出a1,a2,b1,b2的各次辨识结果
title('ParameterIdentificationwithRecursiveLeastSquaresMethod')%图形标题
figure(3);%第三个图形
i=1:
30;%横坐标从1到30
plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:
')%画出a1,a2,b1,b2的各次辨识结果的收敛情况
title('IdentificationPrecision')%图形标题
%2.1二阶系统的最小二乘一次完成算法辨识程序,在光盘中的文件名:
FLch3LSeg1.m
u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.474,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177];%系统辨识的输入信号为一个周期的M序列
y=zeros(1,30);%定义输出观测值的长度
w=0.1*randn(1,30);
fork=3:
30
W(k)=w(k)+1.642*w(k-1)+0.715*w(k-2)
y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+W(k);%用理想输出值作为观测值
end
subplot(3,1,1)%画三行一列图形窗口中的第一个图形
stem(u)%画输入信号u的径线图形
subplot(3,1,2)%画三行一列图形窗口中的第二个图形
i=1:
1:
30;%横坐标范围是1到16,步长为1
plot(i,y)%图形的横坐标是采样时刻i,纵坐标是输出观测值z,图形格式为连续曲线
subplot(3,1,3)%画三行一列图形窗口中的第三个图形
stem(y),gridon%画出输出观测值z的径线图形,并显示坐标网格
u,y%显示输入信号和输出观测信号
%L=30%数据长度
HL=[-y
(2)-y
(1)u
(2)u
(1);-y(3)-y
(2)u(3)u
(2);-y(4)-y(3)u(4)u(3);-y(5)-y(4)u(5)u(4);-y(6)-y(5)u(6)u(5);-y(7)-y(6)u(7)u(6);-y(8)-y(7)u(8)u(7);-y(9)-y(8)u(9)u(8);-y(10)-y(9)u(10)u(9);-y(11)-y(10)u(11)u(10);-y(12)-y(11)u(12)u(11);-y(13)-y(12)u(13)u(12);-y(14)-y(13)u(14)u(13);-y(15)-y(14)u(15)u(14);-y(16)-y(15)u(16)u(15);-y(17)-y(16)u(17)u(16);-y(18)-y(17)u(18)u(17);-y(19)-y(18)u(19)u(18);-y(20)-y(19)u(20)u(19);-y(21)-y(20)u(21)u(20);-y(22)-y(21)u(22)u(21);-y(23)-y(22)u(23)u(22);-y(24)-y(23)u(24)u(23);-y(25)-y(24)u(25)u(24);-y(26)-y(25)u(26)u(25);-y(27)-y(26)u(27)u(26);-y(28)-y(27)u(28)u(27);-y(29)-y(28)u(29)u(28)]%给样本矩阵HL赋值
Y=[y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16);y(17);y(18);y(19);y(20);y(21);y(22);y(23);y(24);y(25);y(26);y(27);y(28);y(29);y(30)]%给样本矩阵zL赋值
%CalculatingParameters
c1=HL'*HL;c2=inv(c1);c3=HL'*Y;c=c2*c3%计算并显示
a1=c
(1),a2=c
(2),b1=c(3),b2=c(4)%从中分离出并显示a1、a2、b1、b2
%2.2RLS递推最小二乘辨识
c0=[1.6420.7150.390.35]';%被辨识参数的初始值
p0=10^6*eye(4,4);%直接给出初始状态P0,即一个充分大的实数单位矩阵
E=0.000000005;%取相对误差E=0.000000005
c=[c0,zeros(4,29)];%被辨识参数矩阵的初始值及大小
e=zeros(4,30);%相对误差的初始值及大小
u=[1.147,0.201,-0.787,-1.589,-1.052,0.866,1.152,1.573,0.626,0.433,-0.958,0.81,-0.044,0.947,-1.474,-0.719,-0.086,-1.099,1.45,1.151,0.485,1.633,0.043,1.326,1.706,-0.340,0.89,1.444,1.177];%系统辨识的输入信号为一个周期的M序列
y=zeros(1,30);%定义输出观测值的长度
w=0.5*randn(1,30);
fork=3:
30
W(k)=w(k)+1.642*w(k-1)+0.715*w(k-2)
y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+W(k);%用理想输出值作为观测值
end
fork=3:
30;%开始求K
h1=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';
x=h1'*p0*h1+1;
x1=inv(x);%开始求K(k)
k1=p0*h1*x1;%求出K的值
d1=y(k)-h1'*c0;
c1=c0+k1*d