最小二乘法.docx
- 文档编号:30684621
- 上传时间:2023-08-19
- 格式:DOCX
- 页数:19
- 大小:183.27KB
最小二乘法.docx
《最小二乘法.docx》由会员分享,可在线阅读,更多相关《最小二乘法.docx(19页珍藏版)》请在冰豆网上搜索。
最小二乘法
最小二乘法
一、最小二乘法概述
最小二乘法是1795年高斯在预测星体运行轨道最先提出的,它奠定了最小二乘估计理论的基础.到了20世纪60年代瑞典学者Austron把这个方法用于动态系统的辨识中,在这种辨识方法中,首先给出模型类型,在该类型下确定系统模型的最优参数。
我们可以将所研究的对象按照对其了解的程度分成白箱、灰箱和黑箱。
于其内部结构、机制只了解一部分,对于其内部运行规律并不十分清楚,这样的研究对象通常称之为“灰箱”;如果我们对于研究对象的内部结构、内部机制及运行规律均一无所知的话,则把这样的研究对象称之为“黑箱”。
研究灰箱和黑箱时,将研究的对象看作是一个系统,通过建立该系统的模型,对模型参数进行辨识来确定该系统的运行规律。
对于动态系统辨识的方法有很多,但其中应用最广泛,辨识效果良好的就是最小二乘辨识方法,研究最小二乘法在系统辨识中的应用具有现实的、广泛的意义。
应用最小二乘法对系统模型参数进行辨识的方法有离线辨识和在线辨识两种离线辨识是在采集到系统模型所需全部输入输出数据后,用最小二乘法对数据进行集中处理,从而获得模型参数的估计值;而在线辨识是一种在系统运行过程中进行的递推辨识方法,所应用的数据是实时采集的系统输入输出数据,应用递推算法对参数估计值进行不断修正,以取得更为准确的参数估计值。
假设一个SISO系统如下图所示:
图1SISO系统结构图
其离散传递函数为:
(1)
输入输出的关系为:
(2)
进一步,我们可以得到:
(3)
其中,扰动量
为均值为0,不相关的白噪声。
将式(3)写成差分方程的形式:
(4)
令
则式(4)可以写为:
(5)
将上述式子扩展到N个输入、输出观测值{
},k=1,2,…,N+n。
将其代入到式(5)中,写成矩阵的形式为:
(6)
其中,
取泛函
为
最小二乘法原理既是使
最小,对其求极值得:
由此可得系统的最小二乘法估计值为:
二、选取模型处理数据后并用最小二乘法进行辨识
选取K=-8,T=50,n=4的系统传递函数;
MATLAB仿真程序如下:
1、首先对系统进行离散化,将离散化后的数据保存以备之后用最小二乘法辨识:
离散程序:
clc;
closeall;
clearall;
DT=1;
ST=1000;
LP=ST/DT;
K=-8;
T=50;
n=4;
x(1:
4)=0;
fori=1:
LP
u0=rand();
x
(1)=x
(1)+(K*u0-x
(1))/T*DT;
x(2:
4)=x(2:
4)+(x(1:
3)-x(2:
4))/T*DT;
y(i)=x(4);
u(i)=u0;
t(i)=i*DT;
end
plot(t,y,'k-')
holdon;
legend('系统输出')
save('lisanhua.mat');
MATLAB输出波形如下:
2、对离散化后保存的数据进行最小二乘法辨识,辨识程序如下:
clc;
clearall;
load('lisanhua.mat')
[Nm]=size(y);
ifN N=m; end n=1; X=zeros(N-n,2*n); Y=zeros(N-n,1); forj=1: N-n fori=1: n X(j,i)=-y(j+n-i); X(j,i+n)=u(j+n-i); end Y(j,1)=y(j); end xita=inv(X'*X)*X'*Y; fortt=1: N-n y_cap(n+tt)=xita'*X(tt,: )'; end forj=1: N-n fori=1: n x_cap(i)=-y_cap(j+n-i); x_cap(n+i)=u(j+n-i); end y_cat(n+j)=xita'*x_cap'; end plot(t,y,'r-') holdon; plot(t,y_cat,'b--') holdon; legend('系统输出','辨识输出') 上述的程序是在n=1时的辨识输出,此时系统输出和辨识输出的波形如下: 当n=2时,此时系统输出和辨识输出的波形如下: 当n=3时,此时系统输出和辨识输出的波形如下: 当n=4和5时,此时系统输出和辨识输出的波形如下: 由上面的曲线可以看出,当n=1、2、3的时候,辨识出来的曲线能够和原系统的输出曲线有较好的拟合程度,当n继续增加时,辨识曲线渐渐偏离原系统的输出曲线。 在使用最小二乘法进行系统辨识时,系统扰动不能为阶跃扰动,应采用非周期的方波或者随机信号。 在使用最小二乘法进行辨识时,其核心在于怎样构建x矩阵和y向量。 由仿真结果可知此次辨识结果较为理想,由最小二乘法拟合的曲线能较好的反应实际系统的模型。 在随机情况下,利用最小二乘法时,并不要求观测数据提供其概率统计方法的信息,而其估计结果,却有相当好的统计特性。 但它具有两方面的缺陷: 一是当模型噪声是有色噪声时,最小二乘估计不是无偏、一致估计;二是随着数据的增长,将出现所谓的“数据饱和”现象。 水箱建模 一、加载水箱数据,并对数据进行数据预处理,将数据进行平滑归零处理,加载水箱数据后输出波形为: clc; clearall; closeall; load('y_tt.mat') plot(tt,u);title('水箱数据输入'); holdon; plot(tt,y);title('水箱数据输出'); holdon; dt=tt (2)-tt (1) dt的输出值为2,说明计算步长为2,再根据水箱的数据输出波形 可知,当波形在800和2000附近输出值较为稳定,故可以取水箱的输入输出值在800到2000之间的数据进行数据归零、平滑等预处理。 数据预处理的MATLAB程序如下: clc; clearall; closeall; load('y_tt.mat') plot(tt,u);title('水箱数据输入'); holdon; plot(tt,y);title('水箱数据输出'); holdon; dt=tt (2)-tt (1) u=u(400: 1000); y=y(400: 1000); t=tt(400: 1000)-400*dt;%时间归零 long=length(y); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%起始点归零%%%%%%%%%%%%%%%%%%%%%%%%%% y1=0; u1=0; fori=1: 5 y1=y1+y(i); u1=u1+u(i); end y1=y1/5;%起始点均值 u1=u1/5;%起始点均值 fori=1: long y11(i)=y(i)-y1; u11(i)=u(i)-u1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%起始点归零%%%%%%%%%%%%%%%%%%%%%%%%%% fori=3: long-2 y11(i)=(y11(i-2)+y11(i-1)+y11(i)+y11(i+1)+y11(i+2))/5; end plot(t,u11,'k');title('预处理后输入'); holdon; plot(t,y11,'k');title('预处理后输出'); holdon; 预处理后的输出波形如下图: 可知,当波形在100和800附近输出值较为稳定,故可以取水箱的输入输出值在100到800之间的数据再进行数据归零处理。 MATLAB程序如下: u=u11(53: 400); y=y11(53: 400); t=tt(53: 400)-53*dt; plot(t,u,'k');title('处理数据输入'); holdon; plot(t,y,'k');title('处理数据输出'); holdon; 波形如下图: 波形仍然有点波动,在程序里再次进行光滑处理: u=u11(53: 400); y=y11(53: 400); t=tt(53: 400)-53*dt; long1=length(y); fornum=1: 50 y1=0; u1=0; N=5; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Æðʼµã¹éÁã%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fori=1: N y1=y1+y(i); u1=u1+u(i); end y1=y1/N; u1=u1/N; fori=1: long1 y12(i)=y(i)-y1; u12(i)=u(i)-u1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Æðʼµã¹éÁã%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fori=3: long1-2 y12(i)=(y12(i-2)+y12(i-1)+y12(i)+y12(i+1)+y12(i+2))/N; y(i)=y12(i); end end plot(t,u,'k');title('处理数据输入'); holdon; plot(t,y12,'k');title('处理数据输出'); holdon; save('shuixiangData.mat') 最终波形为: 二、将处理后的数据保存下来以用于系统辨识,辨识程序如下: clc; clearall; closeall; load('shuixiangData.mat'); [lp,m]=size(y); ifm>lp lp=m; end dt=t (2)-t (1); [ts,Mp,fai,tr,tp,ys,text]=value(y,dt); k=ys Point040=ys*0.4; Point080=ys*0.8; err(1: lp)=0; fori=1: lp err(i)=abs(y(i)-Point040); end [e40,pos1]=min(err); T40=pos1*dt; fori=1: lp err(i)=abs(y(i)-Point080); end [e80,pos2]=min(err); T80=(pos2+1)*dt; order=[1,2,3,4,5,6,7,8,9,10]; err_n(1: 10)=0; forj=1: 10 n1=(1.075*T40)/(T80-T40)+0.5; err_n(j)=abs(sqrt(order(j))-n1); end [en,pos_n]=min(err_n); n=pos_n T=(T40+T80)/(2.16*n) a=exp(-dt/T); b=1-a; x(1: n)=0; y1=[0]; fori=1: lp-1 x (1)=a*x (1)+b*ys*1; ifn>1 x(2: n)=a*x(2: n)+b*x(1: n-1); end ifn==1 y1=[y1x (1)]; else y1=[y1x(n)]; end end plot(t,y,'k-'); holdon; plot(t,y1,'r--'); 系统输出和辨识输出如下图: 参数辨识输出为: k= 0.4800 n= 1 T= 170.3704 三、用穷举法对系统数据进行参数寻优,MATLAB程序如下: DTA_L=0.1; DTA_H=1; D_DTA=0.05; Ti_L=140; Ti_H=240; D_Ti=1; Qmin=10^40; DTAbest=1; Tibest=200; fori=DTA_L: D_DTA: DTA_H forj=Ti_L: D_Ti: Ti_H [Q,DT,LP,t13,u13,y13]=O_obj(i,j,y12); if(Q Qmin=Q; DTAbest=i; Tibest=j; end end end [Q,DT,LP,t13,u13,y13]=O_obj(DTAbest,Tibest,y12); DTAbest Tibest plot(t13,y13) holdon; plot(t13,u13) holdon; 最终的输出曲线如下图: 参数输出为: DTAbest= 0.5000 Tibest= 188
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最小二乘法