[未答复] 采用BFGS算法进行优化无法求解

[复制链接]
秋天的尾影 发表于 2021-10-11 21:27:54
10 财富积分

本帖最后由 秋天的尾影 于 2021-10-17 11:42 编辑

采用BFGS函数进行优化,所用代码之前论坛里发过,应该问题不大,但目前遇到的问题是,我采用下述两个算例进行验证时,算例一能够返回正确结果,但算例二直接报错,经检查发现应该是gk=feval(gfun,x0)  中的feval出问题了,其输出结果与Bk矩阵不匹配导致计算出错。此外,我采用龚纯老师书中的代码也试了一下,两个算例一直是busy,试了一下其他的算例,有些能够解出来,但大部分都是busy。
算例一:
syms t s;
f=100*(s^2-t)^2+(s-1)^2;
df=hessian(f);
x0=[0;0];
[x,val,k]=BFGS('f','df',x0);


算例二:
syms t s r;
f=100*(s^2-t)^2+(s-1)^2+(r-1)^2;
df=hessian(f);
x0=[0;0;0];
[x,val,k]=BFGS('f','df',x0)


BFGS函数
function [x,val,k]=BFGS(fun,gfun,x0)
maxk=500;
rho=0.55;sigma=0.4;epsilon=1e-5;
k=0; n=length(x0);
Bk=eye(n) %Bk=feval('Hess',x0);
while(k<maxk)
    gk=feval(gfun,x0)  
    if(norm(gk)<epsilon),break;end
    dk=-Bk\gk;
    m=0;mk=0;
    while(m<20)
        newf=feval(fun,x0+rho^m*dk);
        oldf=feval(fun,x0);
        if(newf<oldf+sigma*rho^m*gk'*dk)
            mk=m;break;
        end
        m=m+1;
    end
    x=x0+rho^mk*dk;
    sk=x-x0;yk=feval(gfun,x)-gk;
    if(yk'*sk>0)
        Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);
    end
    k=k+1;x0=x;
end
val=feval(fun,x0);


龚纯老师书中的代码
算例一:
syms t s;
var=[t s];
f=100*(s^2-t)^2+(s-1)^2;
x0=[0,0];
[x,minf]=minBFGS(f,x0,var)



function [x,minf] = minBFGS(f,x0,var,eps)
format long;
if nargin == 3
    eps = 1.0e-6;
end
x0 = transpose(x0);
n = length(var);
syms l;
H = eye(n,n);
gradf = jacobian(f,var);
v0  = Funval(gradf,var,x0);
p = -H*transpose(v0);
k = 0;

while 1
    v  = Funval(gradf,var,x0);
    tol = norm(v);
    if tol<=eps
        x = x0;
        break;
    end   
    y = x0 + l*p;
    yf = Funval(f,var,y);
    [a,b] = minJT(yf,0,0.1);
    xm = minHJ(yf,a,b);
    x1 = x0 + xm*p;
    vk = Funval(gradf,var,x1);
    tol = norm(vk);
    if tol<=eps
        x = x1;
        break;
    end
    if k+1==n
        x0 = x1;
        continue;
    else
        dx = x1 - x0;
        dgf = vk - v;
        dgf = transpose(dgf);
        dxT = transpose(dx);
        dgfT = transpose(dgf);
        mdx = dx*dxT;
        mdgf = dgf*dgfT;
        H = H + (1+dgfT*(H*dgf)/(dxT*dgf))*mdx/(dxT*dgf)- ...
            (dx*dgfT*H + H*dgf*dxT)/(dxT*dgf);
        p = -H*transpose(vk);
        k = k+1;
        x0 = x1;
    end
end
minf = Funval(f,var,x);
format short;


您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表