LM拟合优化算法的MATLAB实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序各符号定义列表&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& %x:初始输入序列，为0到6.2以0.2为步长的序列 %y:对应输出序列，[102.225,99.815,-21.585,-35.099,2.523,-38.865,-39.020,89.147,125.249,-63.405 % -183.606,-11.287,197.627,98.355,-131.977,-129.887,52.596,101.193,5.412,-20.805 % ,6.549,-40.176,-71.425,57.366,153.032,5.301,-183.830,-84.612,159.602,155.021 % -73.318,-146.955] %Y_objection:拟合函数，Y_objection=aCoefficient*cos(bCoefficient*x)+bCoefficient*sin(aCoefficient*x) %aCoefficient:a系数 %bCoefficient:b系数 %fRemainder:残差，fRemainder=y-Y_objection=y-aCoefficient*cos(bCoefficient*x)-bCoefficient*sin(aCoefficient*x) %fRemainderNew:新残差，fRemainderNew=y-Y_objection=y-aCoefficient*cos(bCoefficient*(x+hlm))-bCoefficient*sin(aCoefficient*(x+hlm)) %Jacobi:Jacobi矩阵，Jacobi=[-cos(bx)-bx*cos(ax),-sin(ax)+ax*sin(bx)] %gMatrix:过程矩阵，gMatrix=Jacobi'*fRemainder %AMatrix:过程矩阵，AMatrix=Jacobi'*Jacobi %uParameter:过程参数,uParameter=t*max{diag(AMatrix)},t为用户自定义值 %hlm:过程参数,hlm=-inv(uParameter*eye(2)+AMatrix)*gMatrix %Q_parameter:过程参数Q_parameter=(Fx-Fhlm)/L0-Lhlm %Fx=0.5norm((fRemainder))^2,Fhlm=0.5norm((fRemainderNew))^2 %L0-Lhlm=0.5hlm'*(u_parameter*hlm-g_matrix) %初始化各数值 x=0:0.2:6.2; y=[102.225,99.815,-21.585,-35.099,2.523,-38.865,-39.020,89.147,125.249,-63.405, ..., -183.606,-11.287,197.627,98.355,-131.977,-129.887,52.596,101.193,5.412,-20.805, ..., 6.549,-40.176,-71.425,57.366,153.032,5.301,-183.830,-84.612,159.602,155.021, ..., -73.318,-146.955]; x=x'; y=y'; criterion1=10^(-10); criterion2=10^(-10); t=1; kkParameter=0; kmax=20; for aParameter=80:120 for bParameter=80:120 %每次赋给a系数b系数一个初值，在初值范围查找满足条件的最小值 %初始化 kParameter=0; vParameter=2; abParameter=[aParameter;bParameter]; aCoefficient=abParameter(1,1);%a系数初值 bCoefficient=abParameter(2,1);%b系数初值 for i=1:32 %残差初始化 fRemainder(i,1)= y(i,1)-aCoefficient*cos(bCoefficient*x(i,1))-bCoefficient*sin(aCoefficient*x(i,1)); %Jacobi矩阵初始化 Jacobi(i,:)=[-cos(bCoefficient*x(i,1))-bCoefficient*x(i,1)*cos(aCoefficient*x(i,1)), ..., -sin(aCoefficient*x(i,1))+aCoefficient*x(i,1)*sin(bCoefficient*x(i,1))]; end gMatrix=Jacobi'*fRemainder;%g矩阵初始化 AMatrix=Jacobi'*Jacobi;%A矩阵初始化 uParameter=t*max(diag(AMatrix));%u参数初始化 stopCritetion=(norm(gMatrix)<=criterion1); %while循环中找寻满足当前条件的最小值 while((stopCritetion==0) && (kParameter<=kmax)) kParameter=kParameter+1; hlm=-inv(uParameter*eye(2)+AMatrix)*gMatrix; if(norm(hlm)<=(criterion2*(norm(abParameter)+criterion2))) stopCritetion=1; break; else abParameterNew=abParameter+hlm; aCoefficienNew=abParameterNew(1,1); bCoefficienNew=abParameterNew(2,1); for i=1:32 fRemainderNew(i,1)=y(i,1)-aCoefficienNew*cos(bCoefficienNew*x(i,1))-bCoefficienNew*sin(aCoefficienNew*x(i,1)); end QParameter=((norm(fRemainder))^2-(norm(fRemainderNew))^2)/(hlm'*(uParameter*hlm-gMatrix)); if (QParameter>0) abParameter=abParameterNew; aCoefficient=abParameter(1,1); bCoefficient=abParameter(2,1); for i=1:32 fRemainder(i,1)= y(i,1)-aCoefficient*cos(bCoefficient*x(i,1))-bCoefficient*sin(aCoefficient*x(i,1)); Jacobi(i,:)=[-cos(bCoefficient*x(i,1))-bCoefficient*x(i,1)*cos(aCoefficient*x(i,1)), ..., -sin(aCoefficient*x(i,1))+aCoefficient*x(i,1)*sin(bCoefficient*x(i,1))]; end AMatrix=Jacobi'*Jacobi; gMatrix=Jacobi'*fRemainder; stopCritetion=(norm(gMatrix)<criterion1); uParameter=uParameter*max(1/3,1-(2*QParameter-1)^3); vParameter=2; else uParameter=uParameter*vParameter; vParameter=2*vParameter; end end end %记录每次找寻到的最小值 kkParameter=kkParameter+1; Optimal(kkParameter,:)=[aCoefficient ,bCoefficient , norm(fRemainder)]; end end %minValue存储最小残差范数，abCoefficient存储ab系数的数组 [minValue,abCoefficient]=min(Optimal(:,3)); aCoefficient=Optimal(abCoefficient,1) bCoefficient=Optimal(abCoefficient,2) minValue %计算估计值 for i=1:32 Y_objection(i,1) = bCoefficient*cos(bCoefficient*x(i,1))+bCoefficient*sin(aCoefficient*x(i,1)); end %画图 plot(x,y,'r');hold on; plot(x,Y_objection,'b');hold on; legend('实际数据','估计数据');