最优化原理与方法 BFGS算法及PRP算法

实验内容：

编制无约束最优化程序（BFGS算法，PRP算法）求解下列无约束问题 （说明：一维搜索可调用实验一编写的程序）
$$minf(x)=2x_1^2-2x_1{x}_2+x_2^2+2x_1-2x_2$$

---

实验步骤：

1. 理解和分析算法
2. 编写程序
3. 调试程序
4. 求解无约束问题

---

BGFS.m文件

clear all
clc
x=sym('x',[2,1]);
syms y;

y=2*x(1)^2-2*x(1)*x(2)+x(2)^2+2*x(1)-2*x(2);
m_h=zeros(2,1);   %初始x后
B_h=eye(2);          %初始B后
k=1;                     %计数器

while (fanshu(df(m_h,y))>0.001)
           d=B_h\(-df(m_h,y))';     %求解这个线性方程组
           alpha=p2c3(m_h,d);  %一维搜索求得alpha
           m_f=m_h;
           m_h=m_f+alpha*d;
           s=m_h-m_f;
           y_t=(-df(m_h,y))'-(-df(m_f,y))';
           B_f=B_h;
           B_h=B_f-(B_f*s*s'*B_f)/(s'*B_f*s)+(y_t*y_t')/(y_t'*s);
           k=k+1;
end
solver = m_h

其中fanshu、df和p2c3下面将以.m文件实现

fanshu.m

function y=fanshu(x)
y=sqrt(x(1)^2+x(2)^2);

df.m

function y=df(m,f)
x=sym('x',[2,1]);
syms y1 y2
y1=diff(f,x(1));
y2=diff(f,x(2));
t1=subs(y1,x,m);

t2=subs(y2,x,m);

y(1)=double(t1);
y(2)=double(t2);

p2c3.m

%%两点三次插值

%初始步长α记为a=1;步长缩减因子ρ记为rho
%将点μ,μ1,μ2记为m,m1,m2，将函数值φ,φ1,φ2记为f,f1,f2
% 导数值φ',φ1',φ2'记为ff,f11,f22
% 
% clear all;
% clc
function  m=p2c3(x,d)   
a=1;rho=0.1;m1=0;epsilon=0.001;
f1=fun(m1,x,d); f11=fun_t(m1,x,d);

if f11>0
    a=-abs(a);
else
    a=abs(a);
end

m2=m1+a; f2=fun(m2,x,d); f22=fun_t(m2,x,d);

while (f11*f22)>0
    a=2*a;
    m1=m2;
    f1=f2;
    f11=f22;
    m2=m1+a; f2=fun(m2,x,d); f22=fun_t(m2,x,d);
end

z=3*(f2-f1)/(m2-m1)-f11-f22;
w=sign(m2-m1)*sqrt(z^2-f11*f22);
m=m1+(m2-m1)*(1-(f22+w+z)/(f22-f11+2*w));
f=fun(m,x,d); ff=fun_t(m,x,d);

while abs(ff)>epsilon||abs(ff)==epsilon
        a=rho*a;m1=m; f1=f; f11=ff;
        if f11>0
            a=-abs(a);
        else
            a=abs(a);
        end
        
        m2=m1+a; f2=fun(m2,x,d); f22=fun_t(m2,x,d);
        
        while f11*f22>0
            a=2*a;
            m1=m2;
            f1=f2;
            f11=f22;
            m2=m1+a; f2=fun(m2,x,d); f22=fun_t(m2,x,d);
        end
        
        z=3*(f2-f1)/(m2-m1)-f11-f22;
        w=sign(m2-m1)*sqrt(z^2-f11*f22);
        m=m1+(m2-m1)*(1-(f22+w+z)/(f22-f11+2*w));
        f=fun(m,x,d); ff=fun_t(m,x,d);
end
double(m);

当运行BFGS.m文件时，即可得出无约束问题的解

solver =

    0.0000
    1.0000

---

PRP.m文件

clear all
clc
x=sym('x',[2,1]);
syms y;

y=2*x(1)^2-2*x(1)*x(2)+x(2)^2+2*x(1)-2*x(2);
m_f=zeros(2,1);     %初始x前
d=(-df(m_f,y))';      %初始搜索方向
k=1;                       %计数器

if (fanshu(df(m_f,y))<0.001)
    solver=m_f;         
else
     alpha=gold(m_f,d);   %一维搜索求alpha
     m_h=m_f+alpha*d;       %  x后
     belta=0;              %初始值     
     while (fanshu(df(m_h,y))>0.001)
           belta=(df(m_h,y)*(df(m_h,y)-df(m_f,y))')/(df(m_f,y)*df(m_f,y)');
           d=(-df(m_h,y))'+belta*d;
           alpha=gold(m_h,d);  %一维搜索求得alpha
           m_f=m_h;
           m_h=m_f+alpha*d;
           k=k+1;
     end
    solver=m_h
end

其中df与fanshu在上面已经实现，接下来实现gold.m文件

gold.m

%%0.618法求精确一维搜索
%在这我们将α记为a,将αl,αr记为al,ar
%在这我们将φl,φr记为fl,fr
%μ记为m
% clear all
% clc 
function  m=gold(x,d)   
t=(sqrt(5)-1)/2;
[b1,b2]=s_field(x,d);
epsilon=0.0001;
al=b1+(1-t)*(b2-b1);
ar=b1+t*(b2-b1);
fl=fun(al,x,d); fr=fun(ar,x,d);

while  abs(b2-b1)>epsilon
    if fl<fr
        b2=ar;ar=al;fr=fl;
        al=b1+(1-t)*(b2-b1);
        fl=fun(al,x,d);
    else
        b1=al;al=ar;fl=fr;
        ar=b1+t*(b2-b1);
        fr=fun(ar,x,d);
    end
end

if al<ar
    m=al;
else
    m=ar;
end

当运行PRP.m文件时，即可得出无约束问题的解

solver =

    0.0001
    0.9999