matlab基于最小二乘全局化算法遺傳算法的參數(shù)識別

1、最小二乘法：%遞推公式，更新p0=p2;for n=2:N-1%遞推最小二乘法     K0=p0*X(n,:)'*inv(1+X(n,:)*p0*X(n,:)');%計算     Theta_abs=Theta_abs+K0*(Y(n)-X(n,:)*Theta_abs);%計算估計值Theta；     p3=p0-K0*X(n,:)*p0;%計算P     p0=p3;    

2、60;  %誤差平方和最小    Y1=X(n,:)*Theta_abs;%遞推值    J=(Y(n,:)-Y1)*(Y(n,:)-Y1)'    if (J<err)%設(shè)定平方誤差最小，跳出循環(huán)        break;    end;end對于引進后移算子假定在初始條件0時z變換得到ARX模型有：；為均值為0的噪聲項上式可以改寫為上式改寫為最小二乘格式（3）對于（3）式的次觀測構(gòu)成一個線性方程組即.取極小化準則函數(shù)

3、，極小化，求得參數(shù)的估計值，表示為了確定使準則最小的條件，將該式對各參數(shù)求導，并令其結(jié)果等于零：，只要矩陣是滿秩的，則是正定的，使準則為極小的條件得到滿足，最小二乘估計的遞推算法（RLS）最小二乘法，不僅占用大量內(nèi)存，而且不適合于在線辨識，為了解決這個問題，把它轉(zhuǎn)化為遞推算法：若令,則加權(quán)遞推最小二乘（RWLS）：e(k)為有色噪聲，v(k)為白噪聲。取， 當噪聲為有色噪聲時，采用增廣最小二乘法：其思路采用CARMA模型。在實際應用中噪聲v(k)有兩種形式：1 matlab最小二乘法擬合a,Jm=lsqcurvefit(fun,a0,x,y);fun不支持句柄函數(shù)a0為最優(yōu)化的初始值，fun為

4、數(shù)據(jù)原型函數(shù)。x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)； lb x ubx,norm,res,ef,out,lam,jac=lsqcurvefit(F,x0,t,y,v1,v2,opt,P1,P2,.)其中輸出變量的含義為：1） x ： 最優(yōu)解2） norm ： 誤差的平方和3）res: 誤差向量4） ef ： 程序結(jié)束時的狀態(tài)指示：· >0：收斂· 0：函數(shù)調(diào)用次數(shù)或迭代次數(shù)達到最大值（該值在options中指定）· <0：不收斂5) out: 包含以下數(shù)據(jù)的一個結(jié)構(gòu)變量· funcCount 函數(shù)

5、調(diào)用次數(shù)· iterations 實際迭代次數(shù)· cgiterations 實際PCG迭代次數(shù)（大規(guī)模計算用）· algorithm 實際使用的算法· stepsize 最后迭代步長（中等規(guī)模計算用）· firstorderopt 一階最優(yōu)條件滿足的情況（大規(guī)模計算用）6) lam：上下界所對應的Lagrange乘子7) jac：結(jié)果（x點）處的雅可比矩陣輸入?yún)?shù)其中輸入變量的含義為：· x0為初始解（缺省時程序自動取x0=0）· t,y: 擬合數(shù)據(jù)· v1,v2: 參數(shù)待求x的上下界· options：

6、包含算法控制參數(shù)的結(jié)構(gòu)LineSearchType 線搜索方法(cubicpoly,quadcubic(缺省值)opt=optimset(oldopts,newopts)可以設(shè)定的參數(shù)比較多，對lsqnonlin和lsqcurvefit，常用的有以下一些參數(shù)：Diagnostics 是否顯示診斷信息（ 'on' 或'off）Display 顯示信息的級別（'off' ， 'iter' ， 'final，'notify）LargeScale 是否采用大規(guī)模算法（ 'on' 或'off）缺省值為onMa

7、xIter 最大迭代次數(shù)TolFun 函數(shù)計算的誤差限TolX 決策變量的誤差限Jacobian 目標函數(shù)是否采用分析Jacobi矩陣（'on' ，'off）MaxFunEvals 目標函數(shù)最大調(diào)用次數(shù)LevenbergMarquardt 搜索方向選用LM法(on), GN法(off,缺省值)LineSearchType 線搜索方法(cubicpoly,quadcubic(缺省值)LargeScale: on | off LevenbergMarquardt: on | off 例子1用matlab實現(xiàn)對1 首先編寫m文件function f =lsq(x,xdata)

8、f=x(1)*sin(xdata)+0.5*x(2)*sin(2*xdata)./(1-x(2)2*sin(xdata).2).0.52 利用lsqcurvefit函數(shù)調(diào)用m文件m=4;k=0.4o=0:0.01*pi:2*pi;xdata=o;ydata=m*sin(o)+0.5*k*sin(2*o)./(1-(k2*sin(o).2).0.5;x0 = 0; 0;x,resnorm = lsqcurvefit(lsq,x0,xdata,ydata)結(jié)果得到：x(1)= 4.0000;x(2)=0.4000;resnorm = 6.3377e-0162 nlinfitm=4;k=0.4;o=

9、0:0.005*pi:2*pi;xdata=o;ydata=m*sin(o)+0.5*k*sin(2*o)./(1-(k2*sin(o).2).0.5;x0 = 0; 0;beta = nlinfit(xdata,ydata,lsq,x0)例子1.1 用fminunc函數(shù)；k1=13;k2=1.3;k3=9.1;xdata=0:pi/100:pi;ydata= k1*exp(k2*xdata)+k3*sin(xdata);F =(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata);x,fr = fminunc(F,0 0 0)x,fr =

10、 fminsearch(F,1 1 1)例子1.2 用遺傳算法的參數(shù)識別k1=13;k2=1.3;k3=9.1;xdata=0:pi/100:pi;ydata= k1*exp(k2*xdata)+k3*sin(xdata);F =(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata);x,fval = ga(F,3,10;1;8, 20;10;15)%ee=norm(E);%使用差平方和最小原則；或者使用sum(abs(E);%ee=norm(E)/sqrt(n);%使用rms準則例子1.3 利用multistart方法k1=13;k2=1

11、.3;k3=9.1;xdata=0:pi/100:pi;ydata= k1*exp(k2*xdata)+k3*sin(xdata);F =(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata);ms=MultiStart;opts=optimset('Algorithm','interior-point', 'LargeScale','off');problem=createOptimProblem('fmincon','x0',10,1,8,&

12、#39;objective',F,'lb',1,0,1,'ub',20,10,15,'options',opts);xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,200)例子1.4 利用globalsearchk1=13;k2=1.3;k3=9.1;xdata=0:pi/100:pi;ydata= k1*exp(k2*xdata)+k3*sin(xdata);F =(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata);gs = G

13、lobalSearch('Display','iter');opts=optimset('Algorithm','interior-point');problem=createOptimProblem('fmincon','x0',10,1,8,'objective',F,'lb',1,0,1,'ub',20,10,15,'options',opts);xming,fming,flagg,outptg,manyminsg = run(g

14、s,problem)例子1.5 利用multistart和lsqcurvefitk1=13;k2=1.3;k3=9.1;xdata=0:pi/100:pi;ydata= k1*exp(k2*xdata)+k3*sin(xdata);ms=MultiStart;opts=optimset('Algorithm', 'trust-region-reflective');problem=createOptimProblem('lsqcurvefit','x0',10,1,8,'xdata',xdata, 'yda

15、ta',ydata,'objective',myfun,'lb',1,0,1,'ub',20,10,15,'options',opts);xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,100)function y=myfun(x,xdata)y=x(1).*exp(x(2).*xdata)+x(3).*sin(xdata);例子1.6利用multistart和lsqnonlin(8s)k1=13;k2=1.3;k3=9.1;xdata=0:pi/100:pi;ydata

16、= k1*exp(k2*xdata)+k3*sin(xdata);F =(x)x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata;ms=MultiStart;opts=optimset('Algorithm', 'trust-region-reflective');problem=createOptimProblem('lsqnonlin','x0',10,1,8,'xdata',xdata, 'ydata',ydata,'objective'

17、,F,'lb',1,0,1,'ub',20,10,15,'options',opts);xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,100)1.7 利用matlabpool parallel加速tic;k1=13;k2=1.3;k3=9.1;xdata=0:pi/100:pi;ydata= k1*exp(k2*xdata)+k3*sin(xdata);F =(x)x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata;ms=MultiStart('

18、Display','iter','UseParallel','always');opts=optimset('Algorithm', 'trust-region-reflective');matlabpool open 2problem=createOptimProblem('lsqnonlin','x0',10,1,8,'xdata',xdata, 'ydata',ydata,'objective',F,'lb'

19、;,1,0,1,'ub',20,10,15,'options',opts);xminm,fminm,flagm,outptm,manyminsm=run(ms,problem,100)time=toc;%47.3541s函數(shù)模型 function y=myfun(x,xdata)y=x(1)*exp(x(2).*xdata)+x(3)*sin(xdata);例子2（多輸入變量）function F=myfun(k,xdata)F = k(1)*xdata(:,1)+k(2)*xdata(:,2)+k(3)*xdata(:,3);xdata = 3.6 7.7 9.3; 4.1 8.6 2.8; 1.3 7.9 10.0'ydata = 16.5 150.6 263.1'k0 = 0