example3_8.m

这是一个MATLAB6。5 优化计算研究例题的所有源码

function [x,fct,ier,nsub]=minq(gam,c,G,xu,xo,prt,xx);
%%% 初始化 %%%
if prt>0, printlevel=prt, end;
convex=0;
n=size(G,1);
% 最大迭代次数
maxit=3*n;       	
	% 设置初始点
if nargin<7,
      xx=zeros(n,1);
end;
% 初始点位于界限内
xx=max(xu,min(xx,xo));
% 规范化低阶子问题
hpeps=100*eps;
G=G+spdiags(hpeps*diag(G),0,n,n);

K=logical(zeros(n,1));	
if issparse(G), L=speye(n); else L=eye(n); end;
dd=ones(n,1);	
free=logical(zeros(n,1)); 
nfree=0;
nfree_old=-1;

fct=inf; 		% 函数最优值
nsub=0;			% 子空间的数目
unfix=1;		% 允许变量自由变化
nitref=0;		
improvement=1;		

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 主循环：进行坐标和子空间搜索
while 1,
   if prt>1, disp('enter main loop'); end;
   if norm(xx,inf)==inf, error('infinite xx in minq.m'); end;
   g=G*xx+c;
   fctnew=gam+0.5*xx'*(c+g);
   if ~improvement,
    % good termination 
    if prt, 
      disp('terminate: no improvement in coordinate search'); 
    end;
    ier=0; break; 
   elseif nitref>nitrefmax,
    % good termination 
    if prt, disp('terminate: nitref>nitrefmax'); end;
    ier=0; break; 
    elseif nitref>0 & nfree_old==nfree & fctnew >= fct,
    % good termination 
    if prt, 
      disp('terminate: nitref>0 & nfree_old==nfree & fctnew>=fct'); 
    end;
    ier=0; break; 
    elseif nitref==0,
    x=xx;
    fct=min(fct,fctnew);
    if prt>1, fct, end;
    if prt>2, X=x', fct, end;
    else % more accurate g and hence f if nitref>0
    x=xx;
    fct=fctnew;
    if prt>1, fct, end;
    if prt>2, X=x', fct, end;    
    end;
    if nitref==0 & nsub==maxit, 
     if prt,
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'); 
      disp('!!!!!           minq          !!!!!'); 
      disp('!!!!! incomplete minimization !!!!!'); 
      disp('!!!!!   too many iterations   !!!!!'); 
      disp('!!!!!     increase maxit      !!!!!'); 
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
    else
      disp('iteration limit exceeded');
    end;
    ier=99;
    break;
    end;

  % 坐标搜索
  count=0; 	
  k=0;     	
  while 1,
    while count<=n,
      count=count+1;
      if k==n, k=0; end;
      k=k+1;
      if free(k) | unfix, break; end;
    end;
    if count>n, 
      break; 
    end;
    q=G(:,k);
    alpu=xu(k)-x(k); alpo=xo(k)-x(k); % bounds on step

    % 寻找步长大小
    [alp,lba,uba,ier]=getalp(alpu,alpo,g(k),q(k));
    if ier,
      x=zeros(n,1); 
      if lba, x(k)=-1; else x(k)=1; end;
      if prt, 
        gTp=g(k),pTGp=q(k),quot=pTGp/(norm(p,1)^2*norm(G(:),inf))
        disp('minq: function unbounded below in coordinate direction'); 
        disp('      unbounded direction returned'); 
        disp('      possibly caused by roundoff'); 
      end;
      if prt>1, 
        disp('f(alp*x)=gam+gam1*alp+gam2*alp^2/2, where'); 
        gam1=c'*x
        gam2=x'*(G*x)
        ddd=diag(G);
        min_diag_G=min(ddd)
        max_diag_G=max(ddd)
      end;
      return;
    end;
    xnew=x(k)+alp;
    if prt & nitref>0,
      xnew,alp
    end;
    
    if lba | xnew<=xu(k),
      % 下界约束
      if prt>2, disp([num2str(k), ' at lower bound']); end;
      if alpu~=0,
        x(k)=xu(k);
        g=g+alpu*q;
        count=0;
      end;
      free(k)=0;
    elseif uba | xnew>=xo(k),
      % 上界约束
      if prt>2, disp([num2str(k), ' at upper bound']); end;
      if alpo~=0,
        x(k)=xo(k);
        g=g+alpo*q;
        count=0;
      end;
      free(k)=0;
    else
      % no bound active
      if prt>2, disp([num2str(k), ' free']); end;
      if alp~=0,
        if prt>1 & ~free(k), 
          unfixstep=[x(k),alp], 
        end;
        x(k)=xnew;
        g=g+alp*q;
        free(k)=1;
      end;
    end;

  end;
  % end of coordinate search

  nfree=sum(free);
  if (unfix & nfree_old==nfree),
    g=G*x+c;
    nitref=nitref+1;
    if prt>0,
      disp('optimum found; iterative refinement tried');
    end;
  else
    nitref=0;
  end;
  nfree_old=nfree;
  gain_cs=fct-gam-0.5*x'*(c+g);
  improvement=(gain_cs>0 | ~unfix);

  if prt, 
    nfree=pr01('csrch ',free);
  end; 
  if prt, gain_cs, end;
  if prt>2, X=x', end;

  % 子空间搜索
  xx=x; 
  if ~improvement | nitref>nitrefmax,
    % optimal point found - nothing done
  elseif nitref>nitrefmax,
  elseif nfree==0,
    if prt>0,
      disp('no free variables - no subspace step taken')
    end;
    unfix=1;
  else
    minqsub; 
    if ier, return; end;
  end;

  if prt>0, 
    nfree=pr01('ssrch ',free);
    disp(' ');
    if unfix & sum(nfree)<n,
      disp('bounds may be freed in next csearch'); 
    end;
  end; 
end;
% end of main loop
if prt>0, 
  fct
  disp;
end;