24小时热门版块排行榜

返回列表

bingbing5380

铁虫 (正式写手)

应助: 1 (幼儿园)
金币: 72.7
散金: 450
红花: 3
帖子: 331
在线: 144.1小时
虫号: 640275
注册: 2008-10-29
专业: 宏观经济管理与战略

[求助] 《精通matlab最优化计算》（龚纯著）中的BFGS算法为何一直busy？（代码见贴正文)

《精通matlab最优化计算》（龚纯著）中的BFGS算法为何一直busy？（代码见贴正文)
但如果去掉 if k+1==n 这个条件，可以出结果，但是与书上例子的结果有出入。
请教大牛，指点迷津一下？非常感谢。

附--拟牛顿算法BFGS代码：
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;

function [minx,maxx] = minJT(f,x0,h0,eps)
format long;
if nargin == 3
eps = 1.0e-6;
end

x1 = x0;
k = 0;
h = h0;
while 1
x4 = x1 + h;
k = k+1;
f4 = subs(f, findsym(f),x4);
f1 = subs(f, findsym(f),x1);
if f4 < f1
      x2 = x1;
      x1 = x4;
      f2 = f1;
      f1 = f4;
      h = 2*h;
else
      if k==1
         h = -h;
         x2 = x4;
         f2 = f4;
      else
         x3 = x2;
         x2 = x1;
         x1 = x4;
         break;
      end
end
end

minx = min(x1,x3);
maxx = x1+x3 - minx;
format short;

function [x,minf] = minHJ(f,a,b,eps)
format long;
if nargin == 3
eps = 1.0e-6;
end

l = a + 0.382*(b-a);
u = a + 0.618*(b-a);
k=1;
tol = b-a;

while tol>eps && k<100000
fl = subs(f , findsym(f), l);
fu = subs(f , findsym(f), u);
if fl > fu
      a = l;
      l = u;
      u = a + 0.618*(b - a);
else
      b = u;
      u = l;
      l = a + 0.382*(b-a);
end
k = k+1;
tol = abs(b - a);
end
if k == 100000
disp('找不到最小值！');
x = NaN;
minf = NaN;
return;
end
x = (a+b)/2;
minf = subs(f, findsym(f),x);
format short;

回复此楼