%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序各符号定义列表&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%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('实际数据','估计数据');
