最小二乘法

最小二乘法 假设一个SISO系统如下图所示: 图1 SISO系统结构图 其离散传递函数为: ?1?2?n?1bz?bz???bzB(z)12n G ( z ) ? ? (1) ?1A(z)1?a1z?1?a2z?2???anz?n?1 输入输出的关系为: u(k)?G(z?1)?e(k)?y(k) (2) 进一步，我们可以得到: 3) y(k)?A(z?1)?u(k)?B(z?1)?e(k) ( 其中，扰动量e(k)为均值为0，不相关的白噪声。 将式(3)写成差分方程的形式: y(k)??a1y(k?1)?a2y(k?2)???any(k?n)?b1u(k?1)?b2u(k?2)??bnu(k?n)?e(k) (4) 令 ? (k)?[?y(k?1)?y(k?2)??y(k?n)u(k?1)u(k?2)?u(k?n)]T ??[a1a2?anb1b2?bn] 则式(4)可以写为: y(k)??T(k)??e(k) (5) 将上述式子扩展到N个输入、输出观测值{u(k),y(k)}，k=1,2,?，N+n。 将其代入到式(5)中，写成矩阵的形式为: Y????e (6) 其中， Y?[y(n?1)y(n?2)?y(n?N)]T??y(1)u(n) ??y(2)u(n?1) ??? ??y(N)u(n?N)?u(1)??u(2)???????u(N)?e?[e(n?1)e(n?2)?e(n?N)]T?y(n?1)? ?y(n)??y(n?1)?y(n) ?????? ???y(n?N)?y(n?N?1) 取泛函J(?)为 N J(?)??(Y???)2??e2(n?i)?eT?e?(Y???)T(Y???) ?J??[(Y???)T(Y???)]?0????i?1 ??(?T?)T?TY 选取K=-8，T=50，n=4的系统传递函数; G?s?? ?8?50s?1?4 MATLAB仿真程序如下: 离散程序: clc; close all; clear all; DT=1; ST=1000; LP=ST/DT; K=-8; T=50; n=4; x(1:4)=0; for i=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-') hold on; legend('系统输出') save ('lisanhua.mat'); MATLAB输出波形如下: clear all; load('lisanhua.mat') [N m]=size(y); if N<m N=m; end n=1; X=zeros(N-n,2*n); Y=zeros(N-n,1); for j=1:N-n for i=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; for tt=1:N-n y_cap(n+tt)=xita'*X(tt,:)'; end for j=1:N-n for i=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-') hold on; plot(t,y_cat,'b--') hold on; legend('系统输出','辨识输出') 上述的程序是在n=1时的辨识输出，此时系统输出和辨识输出的波形如下: 当n=2时，此时系统输出和辨识输出的波形如下: 当n=3时，此时系统输出和辨识输出的波形如下: 当n=4和5时，此时系统输出和辨识输出的波形如下: 由上面的曲线可以看出，当n=1、2、3的时候，辨识出来的曲线能够和原系统的输出曲线有较好的拟合程度，当n继续增加时，辨识曲线渐渐偏离原系统的输出曲线。 水箱建模 一、 加载水箱数据，并对数据进行数据预处理，将数据进行平滑归零处理，加载水箱数据 后输出波形为: clc; clear all; close all; load('y_tt.mat') plot(tt,u);title('水箱数据输入'); hold on; plot(tt,y);title('水箱数据输出'); hold on; dt=tt(2)-tt(1) dt的输出值为2，说明计算步长为2，再根据水箱的数据输出波形 可知，当波形在800和2000附近输出值较为稳定，故可以取水箱的输入输出值在800到2000之间的数据进行数据归零、平滑等预处理。数据预处理的MATLAB程序如下: clc; clear all; close all; load('y_tt.mat') plot(tt,u);title('水箱数据输入'); hold on; plot(tt,y);title('水箱数据输出'); hold on; 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; for i=1:5 y1=y1+y(i); u1=u1+u(i); end y1=y1/5; % 起始点均值 u1=u1/5; % 起始点均值 for i=1:long y11(i)=y(i)-y1; u11(i)=u(i)-u1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%起始点归零%%%%%%%%%%%%%%%%%%%%%%%%%% for i=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('预处理后输入'); hold on; plot(t,y11,'k');title('预处理后输出'); hold on; 预处理后的输出波形如下图: 可知，当波形在100和800附近输出值较为稳定，故可以取水箱的输入输出值在100到800之间的数据再进行数据归零处理。MATLAB程序如下: u=u11(53:400); y=y11(53:400); t=tt(53:400)-53*dt; plot(t,u,'k');title('处理数据输入'); hold on; plot(t,y,'k');title('处理数据输出'); hold on; 波形如下图: 波形仍然有点波动，在程序里再次进行光滑处理: u=u11(53:400); y=y11(53:400); t=tt(53:400)-53*dt; long1=length(y); for num=1:50 y1=0; u1=0; N=5; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%?eê?μ?1éá?%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:N y1=y1+y(i); u1=u1+u(i); end y1=y1/N; u1=u1/N; for i=1:long1 y12(i)=y(i)-y1; u12(i)=u(i)-u1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%?eê?μ?1éá?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=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('处理数据输入'); hold on; plot(t,y12,'k');title('处理数据输出'); hold on; save ('shuixiangData.mat') 最终波形为: 二、将处理后的数据保存下来以用于系统辨识,辨识程序如下: clc; clear all; close all; load ('shuixiangData.mat'); [lp,m]=size(y); if m>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; for i=1:lp err(i)=abs(y(i)-Point040); end [e40,pos1]=min(err); T40=pos1*dt; for i=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; for j=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]; for i=1:lp-1 x(1)=a*x(1)+b*ys*1; if n>1 x(2:n)=a*x(2:n)+b*x(1:n-1); end if n==1 y1=[y1 x(1)]; else y1=[y1 x(n)]; end end plot(t,y,'k-'); hold on; 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; for i=DTA_L:D_DTA:DTA_H for j=Ti_L:D_Ti:Ti_H [Q,DT,LP,t13,u13,y13]=O_obj(i,j,y12); if(Q<Qmin) 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) hold on; plot(t13,u13) hold on; 最终的输出曲线如下图: 参数输出为: DTAbest = 0.5000 Tibest = 188