%{
程序功能：
1、变度量法算法（DFP）求解无约束问题
2、调用文件夹下Newton的子函数：nfx,ndfx,ndfx2,vectorLength
3、z3=A(:,:,i)\b;%计算当前d的值
矩阵计算可能存在奇异值
4、请根据不同的目标函数，设置精度、迭代次数、初始迭代值。
5、迭代初始点的选取很重要
Name：李承霖
Num:35020181152203
Date:2018.11.07
%}
%初始化程序
function SearchDFP()
clear
clc
syms x1 x2
format rat
eps=0.001; %精度要求
k=1;
n=10;%设置最大迭代次数、
x=zeros(n,2);
% x(1,:)=[1 1];
%设置初始迭代点
x(1,:)=input('Please input initial point matrix:');
d=zeros(n,2);
dfx1=zeros(n,2);%一阶导数值fx
A=zeros(2,2,n);%二阶导数矩阵A
% beta=zeros(n,1);
alpha=zeros(n,1);%搜索步长
v=ndfx(x1,x2);%一阶导数
z1=ndfx2(x1,x2);%二阶导数
H=zeros(2,2,n);%置初始矩阵
H(:,:,1)=eye(2,2);%初始化为单位矩阵
s=zeros(n,2);
y=zeros(n,2);
for i=1:n
f=subs(v,x1,x(i,1));
dfx1(i,:)=subs(f,x2,x(i,2));%此处计算一阶导数的值
%此处精度设置要慎重
% m=sqrt( dfx1(i,1)^2+dfx1(i,2)^2); %取矢量长度
m=vectorLength(dfx1(i,:));
if( m<eps)
break;%程序终止
end
z2=subs(z1,x1,x(i,1));
% A=subs(z2,x2,x(i,2));
A(:,:,i)=subs(z2,x2,x(i,2));
%DFP 核心算法
%------------------------------------
g=dfx1(i,:);%
G=A(:,:,i);
if(i==1)
d(i,:)=-H(:,:,i)*dfx1(i,:)';
else
%第二轮循环
k=i-1;
s(k,:)=x(i,:)-x(i-1,:);
y(k,:)=dfx1(i,:)-dfx1(i-1,:);
t1=y(k,:)*H(:,:,k)*y(k,:)';
t2=s(k,:)*y(k,:)';
m1=H(:,:,k)*y(k,:)'*y(k,:)*H(:,:,k);
m2=s(k,:)'*s(k,:);
H(:,:,i)=H(:,:,i-1)-m1/t1+m2/t2;
m3=-H(:,:,i)*dfx1(i,:)'; %小心奇异矩阵
d(i,:)=m3';
end
alpha(i)=-d(i,:)*g'/(d(i,:)*G*d(i,:)');
x(i+1,:)=x(i,:)+alpha(i)*d(i,:) ; %第一轮循环
%BFGS 核心算法
%------------------------------------
%第一轮循环
%{
g=dfx1(i,:);%
G=A(:,:,i);
if(i==1)
d(i,:)=-inv(B(:,:,i))*dfx1(i,:)';
else
%第二轮循环
k=i-1;
s(k,:)=x(i,:)-x(i-1,:);
y(k,:)=dfx1(i,:)-dfx1(i-1,:);
t1=s(k,:)*B(:,:,k)*s(k,:)';
t2=y(k,:)*s(k,:)';
m1=B(:,:,k)*s(k,:)'*s(k,:)*B(:,:,k);
m2=y(k,:)'*y(k,:);
B(:,:,i)=B(:,:,i-1)-m1/t1+m2/t2;
m3=-inv(B(:,:,i))*dfx1(i,:)'; %小心奇异矩阵
d(i,:)=m3';
end
alpha(i)=-d(i,:)*g'/(d(i,:)*G*d(i,:)');
x(i+1,:)=x(i,:)+alpha(i)*d(i,:) ; %第一轮循环
%}
%------------------------------------
%FR 核心算法
%------------------------------------
%{
%下面计算一维搜索步长
g=dfx1(i,:);
G=A(:,:,i);
if(i==1)
beta(i)=0;
d(i,:)=-g;
else
beta(i)=abs( vectorLength(dfx1(i,:))^2/vectorLength(dfx1(i-1,:))^2 );
d(i,:)=-dfx1(i,:)+beta(i)*d(i-1,:);
end
d1=d(i,:);
alpha(i)=-d1*g'/(d1*G*d1');
x(i+1,:)=x(i,:)+alpha(i)*d1;
%}
%------------------------------------
%{
%下面是newton搜索无约束的核心算法
%计算二阶导数方程的解
b=-dfx1(i,:)';
z1=ndfx2(x1,x2);
z2=subs(z1,x1,x(i,1));
% A=subs(z2,x2,x(i,2));
A(:,:,i)=subs(z2,x2,x(i,2));
z3=A(:,:,i)\b;%计算当前d的值
d(i,:)=z3';
x(i+1,:)=x(i,:)+d(i,:); %第二轮循环
%此处精度设置要慎重
if(abs(x(i+1,1)-x(i,1))<eps && abs(x(i+1,2)-x(i,2))<eps)
break;%程序终止
end
%}
end
disp('The result is:')
% d,x,fx,A
d=d(1:i,:)
x=x(1:i,:)
dfx1=dfx1(1:i,:)
A=A(:,:,1:i)
alpha=alpha(1:i,:)
H=H(:,:,1:i)
% beta=beta(1:i,:)
fprintf('程序尝试迭代次数：%d\n',i)
fprintf('最优解：%f ,%f\n',x(i,:))
fprintf('最优目标函数值：y=%.3f\n',nfx(x(i,1),x(i,2)))
end
function y= nfx( x1,x2 )
%UNTITLED8 此处显示有关此函数的摘要
% 此处显示详细说明
%目标函数
y=2*x1^2-2*x1*x2+x2^2+2*x1-2*x2;
%
% y=3*x1*x1/2+x2*x2/2-x1*x2-2*x1;
% y=(x1^2)/3+ exp(x2);%目标函数,要避免导函数分母除0的情况出现。
end
function y=ndfx(x1,x2)
y(1)=diff(nfx(x1,x2),x1);
y(2)=diff(nfx(x1,x2),x2); %对目标函数取一阶偏微分
end
function y=ndfx2(x1,x2)
%z=diff(nfx(x1,x2),x1)+diff(nfx(x1,x2),x2); %对目标函数取二阶偏微分
syms y
% y=zeros(2,2);
% sym(y)
% z=diff((nfx(x1,x2)),x1);
% y2=diff((nfx(x1,x2)),x2);
z=ndfx(x1,x2);
y(1,1)=diff(z(1),x1);%二阶
y(1,2)=diff(z(1),x2);%二阶
y(2,1)=diff(z(2),x1);%二阶
y(2,2)=diff(z(2),x2);%二阶
% y=[diff(nfx(x1,x2),x1,2),diff(y1,x2);diff(y2,x1),diff(nfx(x1,x2),x2,2)];
% y=[diff(z(1),x1),diff(z(1),x2);diff(z(2),x1), diff(z(2),x2)];%构造二阶导数矩阵
y=[y(1,1),y(1,2);y(2,1),y(2,2)];
end
%{
程序功能：
1、计算一个矢量的长度
%}
function y=vectorLength(x)
n=length(x);
sum=0;
for i=1:n
sum=sum+x(i)*x(i);
end
y=sqrt(sum);
end
