nlsqcon.m

基于matlab的约束非线性规划算法库

function [Xsol, Residual, Aset, Abounds, term, Isave, Iter] = ...
    NLSQcon(ffun,cfun,X0,Teq,Tineq,LB,UB,Maxiter,epsrel,epsabs,epscon,params)
%Call: [Xsol, Residual, Aset, Abounds, term, Isave, Iter] = ...
%           NLSQcon(ffun,cfun,X0,Teq,Tineq,LB,UB,Maxiter,epsrel,epsabs,epscon,{m,t,y})
%If no parameters "params" is used the call is done like
%           NLSQcon(ffun,cfun,X0,Teq,Tineq,LB,UB,Maxiter,epsrel,epsabs,epscon,{:})
%with {:} to indicate an empty string,
%where 'ffun' and 'cfun' are the names of the user defined functions
%that defines the residual vector f(x) and the constraint function
%c(x) respectively or their corresponding Jacobian matrix at the point x.
%f(x) should be m*1 and corresponding Jacobian is m*n
%c(x) should be (Teq+Tineq)*1 and corresponding Jacobian (Teq+Tineq)*n
%The problem to solve is (Note the different usage of indeces i and j!)
%          min 0.5*norm2(f(x)) over x that is n*1
%             subject to
%          ci(x) =  0 i=1,2,....Teq
%          cj(x) >= 0 j=Teq+1,....Teq+Tineq
%          LB(j) <= x(j) <= UB(j) j=1,2,......,n
%X0 the starting point should be feasible w.r.t. the bounds. If not,
%the starting point is adjusted.
%Maxit is the maximum number of allowed iterations until termination.
%To be a well-posed problem Teq should be <= n and m >= n.
%The user-written functions 'ffun' and 'cfun' should be written as
%Example
%----------------for example-----------
%function [f, J] = fmukai(x, ctrl,m,t,y)
%n=max(size(x));
%for i=1:m
%    f(i)=modell(x)-y(i);
%end
%f=f(:);
%The Jacobian is set to the empty matix to indicate that it should
%be computed using forward differences
%J=[];
%----------------for example------------
%function [h, A] = cmukai(x, ctrl,m,t,y)
%h(1)=-x(1)^2-x(2)^2-x(3)^2+5;
%h(2)=-(x(4)-3)^2-x(5)^2+1;
%h=h(:);
%n=max(size(x));
%if ctrl > 0
%   A=zeros(2,n);
%   A(1,1)=-2*x(1);
%   A(1,2)=-2*x(2);
%   A(1,3)=-2*x(3);
%   A(2,4)=-2*(x(4)-3);
%   A(2,5)=-2*x(5);
%else
%   A=[];
%end
%----------------end example-----------

kappa = 0.1*ones(4,1); %Initial sequence of weights in the merit function
Wkm1 = kappa(1);
alphakm1 = 1;
betakm1 = 1;
k=0;
phase=1; % 1 = Gauss-Newton
Aset=[]; % holds indeces of general constraints in current working set
%Abounds indicates if some bounds are attained at the current point
%Check input parameters
[m, n, X, LB, UB, f, c, Aset, Abounds, exit] = CheckIn(ffun,cfun,X0,Teq,Tineq,LB,UB,Maxiter,params);
term = exit;
Isave = [];
while ( term == 0) & (k < Maxiter)
       [Aset, Abounds, T, GNlambda, p, Jc, b, betak, lenJp, cTc, JpTf, Jp] = Direction(phase,m,n,X,f,c,k,Teq,...
          Tineq,LB,UB,Aset,Abounds,ffun,cfun,params);
       if isempty(b)
           bTb = 0;
       else
           bTb = dot(b,b);
       end
       fTf = dot(f,f);
       term = TermCrit(p,betak,f,epsrel,epsabs,cTc,epscon);
       if term == 0
          if phase == 2 %Newton
              alpha = 1;
              [f dummy] = feval(ffun,X+p,0,params{:});
              [c dummy] = feval(cfun,X+p,0,params{:});
              Wk = -1; %signal inside PrintInfo
          else
              [kappa, Wk] = Weights(kappa,Wkm1,cTc,JpTf,lenJp);
              [alphaArm, fArm, cArm] = Steplength(Aset,Teq,Tineq,f,c,b,Abounds,X,LB,UB,p,Jc,JpTf,Wk,ffun,cfun,...
              alphakm1,params);
          %
              [alphaINM, fINM, cINM] = SteplengthINM(Aset,Teq,Tineq,f,c,b,Abounds,X,LB,UB,p,Jc,JpTf,Jp,Wk,ffun,cfun,...
              alphakm1,params);
          %alphaArm
          %fArmSqr=dot(fArm,fArm)
          %cArmSqr=dot(cArm,cArm)
          %alphaINM
          %fINMsqr=dot(fINM,fINM)
          %cINMsqr=dot(cINM,cINM)
          %pause
              alpha = alphaArm;
              f = fArm;
              c = cArm;
          end % of if phase == 2
          Isave = PrintInfo(Isave,k,Aset,Abounds,X,fTf,bTb,Teq,Tineq,m,n,p,betak,betakm1,Wk,alpha);
          X = X + alpha*p;
          [phase, term, X, Aset, Abounds, alphakm1, Wkm1, betakm1] = ...
             Analysis(k,X,n,p,f,c,Teq,Tineq,LB,UB,alpha,alphakm1,Wk,Wkm1,betak,betakm1,phase,term);
          k = k+1;
      end % of if ~term
end % of while
Xsol = X;
Residual = f;
Iter = k;
if isempty(Aset)
    bTb=0;
else
    b=c(Aset);
    bTb=dot(b,b);
end
Isave = PrintInfo(Isave,k,Aset,Abounds,X,f'*f,bTb,Teq,Tineq,m,n,p,betak,betakm1,Inf,Inf);
disp('No. of iterations: current and Maxiter')
[Iter Maxiter]
disp('       fTf          bTb        Weight       betak       speed         pTp            alpha')
format short e
disp(Isave(:,2:8))