time   calls  line
                  1 function [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = fmincon(FUN,X,A,B,Aeq,Beq,LB,UB,NONLCON,options,varargin)
                  2 %FMINCON finds a constrained minimum of a function of several variables.
                  3 %   FMINCON attempts to solve problems of the form:
                  4 %    min F(X)  subject to:  A*X  <= B, Aeq*X  = Beq (linear constraints)
                  5 %     X                     C(X) <= 0, Ceq(X) = 0   (nonlinear constraints)
                  6 %                              LB <= X <= UB        (bounds)
                  7 %    
                  8 %   FMINCON implements four different algorithms: interior point, SQP, active 
                  9 %   set, and trust region reflective. Choose one via the option Algorithm: 
                 10 %   for instance, to choose SQP, set OPTIONS = optimset('Algorithm','sqp'), 
                 11 %   and then pass OPTIONS to FMINCON. 
                 12 %                                                           
                 13 %   X = FMINCON(FUN,X0,A,B) starts at X0 and finds a minimum X to the 
                 14 %   function FUN, subject to the linear inequalities A*X <= B. FUN accepts 
                 15 %   input X and returns a scalar function value F evaluated at X. X0 may be
                 16 %   a scalar, vector, or matrix. 
                 17 %
                 18 %   X = FMINCON(FUN,X0,A,B,Aeq,Beq) minimizes FUN subject to the linear 
                 19 %   equalities Aeq*X = Beq as well as A*X <= B. (Set A=[] and B=[] if no 
                 20 %   inequalities exist.)
                 21 %
                 22 %   X = FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB) defines a set of lower and upper
                 23 %   bounds on the design variables, X, so that a solution is found in 
                 24 %   the range LB <= X <= UB. Use empty matrices for LB and UB
                 25 %   if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; 
                 26 %   set UB(i) = Inf if X(i) is unbounded above.
                 27 %
                 28 %   X = FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) subjects the minimization
                 29 %   to the constraints defined in NONLCON. The function NONLCON accepts X 
                 30 %   and returns the vectors C and Ceq, representing the nonlinear 
                 31 %   inequalities and equalities respectively. FMINCON minimizes FUN such 
                 32 %   that C(X) <= 0 and Ceq(X) = 0. (Set LB = [] and/or UB = [] if no bounds
                 33 %   exist.)
                 34 %
                 35 %   X = FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS) minimizes with 
                 36 %   the default optimization parameters replaced by values in the structure
                 37 %   OPTIONS, an argument created with the OPTIMSET function. See OPTIMSET
                 38 %   for details. For a list of options accepted by FMINCON refer to the
                 39 %   documentation.
                 40 %  
                 41 %   X = FMINCON(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a
                 42 %   structure with the function FUN in PROBLEM.objective, the start point
                 43 %   in PROBLEM.x0, the linear inequality constraints in PROBLEM.Aineq
                 44 %   and PROBLEM.bineq, the linear equality constraints in PROBLEM.Aeq and
                 45 %   PROBLEM.beq, the lower bounds in PROBLEM.lb, the upper bounds in 
                 46 %   PROBLEM.ub, the nonlinear constraint function in PROBLEM.nonlcon, the
                 47 %   options structure in PROBLEM.options, and solver name 'fmincon' in
                 48 %   PROBLEM.solver. Use this syntax to solve at the command line a problem 
                 49 %   exported from OPTIMTOOL. The structure PROBLEM must have all the fields.
                 50 %
                 51 %   [X,FVAL] = FMINCON(FUN,X0,...) returns the value of the objective 
                 52 %   function FUN at the solution X.
                 53 %
                 54 %   [X,FVAL,EXITFLAG] = FMINCON(FUN,X0,...) returns an EXITFLAG that 
                 55 %   describes the exit condition of FMINCON. Possible values of EXITFLAG 
                 56 %   and the corresponding exit conditions are listed below. See the
                 57 %   documentation for a complete description.
                 58 %   
                 59 %   All algorithms:
                 60 %     1  First order optimality conditions satisfied.
                 61 %     0  Too many function evaluations or iterations.
                 62 %    -1  Stopped by output/plot function.
                 63 %    -2  No feasible point found.
                 64 %   Trust-region-reflective, interior-point, and sqp:
                 65 %     2  Change in X too small.
                 66 %   Trust-region-reflective:
                 67 %     3  Change in objective function too small.
                 68 %   Active-set only:
                 69 %     4  Computed search direction too small.
                 70 %     5  Predicted change in objective function too small.
                 71 %   Interior-point and sqp:
                 72 %    -3  Problem seems unbounded.
                 73 %
                 74 %   [X,FVAL,EXITFLAG,OUTPUT] = FMINCON(FUN,X0,...) returns a structure 
                 75 %   OUTPUT with information such as total number of iterations, and final 
                 76 %   objective function value. See the documentation for a complete list.
                 77 %
                 78 %   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = FMINCON(FUN,X0,...) returns the 
                 79 %   Lagrange multipliers at the solution X: LAMBDA.lower for LB, 
                 80 %   LAMBDA.upper for UB, LAMBDA.ineqlin is for the linear inequalities, 
                 81 %   LAMBDA.eqlin is for the linear equalities, LAMBDA.ineqnonlin is for the
                 82 %   nonlinear inequalities, and LAMBDA.eqnonlin is for the nonlinear 
                 83 %   equalities.
                 84 %
                 85 %   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD] = FMINCON(FUN,X0,...) returns the 
                 86 %   value of the gradient of FUN at the solution X.
                 87 %
                 88 %   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = FMINCON(FUN,X0,...) 
                 89 %   returns the value of the exact or approximate Hessian of the Lagrangian
                 90 %   at X. 
                 91 %
                 92 %   Examples
                 93 %     FUN can be specified using @:
                 94 %        X = fmincon(@humps,...)
                 95 %     In this case, F = humps(X) returns the scalar function value F of 
                 96 %     the HUMPS function evaluated at X.
                 97 %
                 98 %     FUN can also be an anonymous function:
                 99 %        X = fmincon(@(x) 3*sin(x(1))+exp(x(2)),[1;1],[],[],[],[],[0 0])
                100 %     returns X = [0;0].
                101 %
                102 %   If FUN or NONLCON are parameterized, you can use anonymous functions to
                103 %   capture the problem-dependent parameters. Suppose you want to minimize 
                104 %   the objective given in the function myfun, subject to the nonlinear 
                105 %   constraint mycon, where these two functions are parameterized by their 
                106 %   second argument a1 and a2, respectively. Here myfun and mycon are 
                107 %   MATLAB file functions such as
                108 %
                109 %        function f = myfun(x,a1)      
                110 %        f = x(1)^2 + a1*x(2)^2;       
                111 %                                      
                112 %        function [c,ceq] = mycon(x,a2)
                113 %        c = a2/x(1) - x(2);
                114 %        ceq = [];
                115 %
                116 %   To optimize for specific values of a1 and a2, first assign the values 
                117 %   to these two parameters. Then create two one-argument anonymous 
                118 %   functions that capture the values of a1 and a2, and call myfun and 
                119 %   mycon with two arguments. Finally, pass these anonymous functions to 
                120 %   FMINCON:
                121 %
                122 %        a1 = 2; a2 = 1.5; % define parameters first
                123 %        options = optimset('Algorithm','interior-point'); % run interior-point algorithm
                124 %        x = fmincon(@(x) myfun(x,a1),[1;2],[],[],[],[],[],[],@(x) mycon(x,a2),options)
                125 %
                126 %   See also OPTIMSET, OPTIMTOOL, FMINUNC, FMINBND, FMINSEARCH, @, FUNCTION_HANDLE.
                127 
                128 %   Copyright 1990-2010 The MathWorks, Inc.
                129 %   $Revision: 1.1.6.15 $  $Date: 2010/10/08 17:17:22 $
                130 
  0.15     860  131 defaultopt = struct( ... 
                132     'Algorithm','trust-region-reflective', ...
                133     'AlwaysHonorConstraints','bounds', ...
                134     'DerivativeCheck','off', ...
                135     'Diagnostics','off', ...
                136     'DiffMaxChange',Inf, ...
                137     'DiffMinChange',0, ...
                138     'Display','final', ...
                139     'FinDiffType','forward', ...    
                140     'FunValCheck','off', ...
                141     'GradConstr','off', ...
                142     'GradObj','off', ...
                143     'HessFcn',[], ...
                144     'Hessian',[], ...    
                145     'HessMult',[], ...
                146     'HessPattern','sparse(ones(numberOfVariables))', ...
                147     'InitBarrierParam',0.1, ...
                148     'InitTrustRegionRadius','sqrt(numberOfVariables)', ...
                149     'LargeScale','on', ...
                150     'MaxFunEvals',[], ...
                151     'MaxIter',[], ...
                152     'MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
                153     'MaxProjCGIter','2*(numberOfVariables-numberOfEqualities)', ...    
                154     'MaxSQPIter','10*max(numberOfVariables,numberOfInequalities+numberOfBounds)', ...
                155     'NoStopIfFlatInfeas','off', ...
                156     'ObjectiveLimit',-1e20, ...
                157     'OutputFcn',[], ...
                158     'PhaseOneTotalScaling','off', ...
                159     'PlotFcns',[], ...
                160     'PrecondBandWidth',0, ...
                161     'RelLineSrchBnd',[], ...
                162     'RelLineSrchBndDuration',1, ...
                163     'ScaleProblem','obj-and-constr', ...
                164     'SubproblemAlgorithm','ldl-factorization', ...
                165     'TolCon',1e-6, ...
                166     'TolConSQP',1e-6, ...    
                167     'TolFun',1e-6, ...
                168     'TolGradCon',1e-6, ...
                169     'TolPCG',0.1, ...
                170     'TolProjCG',1e-2, ...
                171     'TolProjCGAbs',1e-10, ...
                172     'TolX',[], ...
                173     'TypicalX','ones(numberOfVariables,1)', ...
                174     'UseParallel','never' ...
                175     );
                176 
                177 % If just 'defaults' passed in, return the default options in X
           860  178 if nargin==1 && nargout <= 1 && isequal(FUN,'defaults') 
                179    X = defaultopt;
                180    return
                181 end
                182 
  0.01     860  183 if nargin < 10 
                184     options = [];
                185     if nargin < 9
                186         NONLCON = [];
                187         if nargin < 8
                188             UB = [];
                189             if nargin < 7
                190                 LB = [];
                191                 if nargin < 6
                192                     Beq = [];
                193                     if nargin < 5
                194                         Aeq = [];
                195                     end
                196                 end
                197             end
                198         end
                199     end
                200 end
                201 
           860  202 problemInput = false; 
           860  203 if nargin == 1 
                204     if isa(FUN,'struct')
                205         problemInput = true;
                206         [FUN,X,A,B,Aeq,Beq,LB,UB,NONLCON,options] = separateOptimStruct(FUN);
                207     else % Single input and non-structure.
                208         error('optim:fmincon:InputArg','The input to FMINCON should be either a structure with valid fields or consist of at least four arguments.' );
                209     end
                210 end
                211 
           860  212 if nargin < 4 && ~problemInput 
                213   error('optim:fmincon:AtLeastFourInputs','FMINCON requires at least four input arguments.')
                214 end
                215 
  0.01     860  216 if isempty(NONLCON) && isempty(A) && isempty(Aeq) && isempty(UB) && isempty(LB) 
                217    error('optim:fmincon:ConstrainedProblemsOnly', ...
                218          'FMINCON is for constrained problems. Use FMINUNC for unconstrained problems.')
                219 end
                220 
                221 % Check for non-double inputs
                222 % SUPERIORFLOAT errors when superior input is neither single nor double;
                223 % We use try-catch to override SUPERIORFLOAT's error message when input
                224 % data type is integer.
  0.01     860  225 try 
  0.01     860  226     dataType = superiorfloat(X,A,B,Aeq,Beq,LB,UB); 
                227 catch ME
                228     if strcmp(ME.identifier,'MATLAB:datatypes:superiorfloat')
                229        dataType = 'notDouble';
                230     end
                231 end
                232 
           860  233 if ~strcmp(dataType,'double') 
                234     error('optim:fmincon:NonDoubleInput', ...
                235         'FMINCON only accepts inputs of data type double.')
                236 end
                237 
           860  238 if nargout > 4 
                239    computeLambda = true;
           860  240 else  
           860  241    computeLambda = false; 
           860  242 end 
                243 
           860  244 activeSet = 'medium-scale: SQP, Quasi-Newton, line-search'; 
           860  245 sqp = 'sequential quadratic programming'; 
           860  246 trustRegionReflective = 'trust-region-reflective'; 
           860  247 interiorPoint = 'interior-point'; 
                248 
           860  249 XOUT = X(:); 
           860  250 numberOfVariables=length(XOUT); 
                251 % Check for empty X
           860  252 if numberOfVariables == 0 
                253    error('optim:fmincon:EmptyX','You must provide a non-empty starting point.');
                254 end
                255 
  0.07     860  256 display = optimget(options,'Display',defaultopt,'fast'); 
  0.04     860  257 flags.detailedExitMsg = ~isempty(strfind(display,'detailed')); 
           860  258 switch display 
  0.01     860  259     case {'off','none'} 
  0.01     860  260         verbosity = 0; 
                261     case {'notify','notify-detailed'}
                262         verbosity = 1;
                263     case {'final','final-detailed'}
                264         verbosity = 2;
                265     case {'iter','iter-detailed'}
                266         verbosity = 3;
                267     case 'testing'
                268         verbosity = 4;
                269     otherwise
                270         verbosity = 2;
                271 end
                272 
                273 % Set linear constraint right hand sides to column vectors
                274 % (in particular, if empty, they will be made the correct
                275 % size, 0-by-1)
           860  276 B = B(:); 
           860  277 Beq = Beq(:); 
                278 
                279 % Check for consistency of linear constraints, before evaluating
                280 % (potentially expensive) user functions 
                281 
                282 % Set empty linear constraint matrices to the correct size, 0-by-n
           860  283 if isempty(Aeq) 
  0.01     860  284     Aeq = reshape(Aeq,0,numberOfVariables); 
           860  285 end 
           860  286 if isempty(A) 
           860  287     A = reshape(A,0,numberOfVariables);    
           860  288 end 
                289 
           860  290 [lin_eq,Aeqcol] = size(Aeq); 
  0.01     860  291 [lin_ineq,Acol] = size(A); 
                292 % These sizes checks assume that empty matrices have already been made the correct size
           860  293 if Aeqcol ~= numberOfVariables 
                294    error('optim:fmincon:WrongNumberOfColumnsInAeq','Aeq must have %i column(s).',numberOfVariables)
                295 end
           860  296 if lin_eq ~= length(Beq) 
                297     error('optim:fmincon:AeqAndBeqInconsistent', ...
                298         'Row dimension of Aeq is inconsistent with length of beq.')
                299 end
           860  300 if Acol ~= numberOfVariables 
                301    error('optim:fmincon:WrongNumberOfColumnsInA','A must have %i column(s).',numberOfVariables)
                302 end
           860  303 if lin_ineq ~= length(B) 
                304     error('optim:fmincon:AeqAndBinInconsistent', ...
                305         'Row dimension of A is inconsistent with length of b.')
                306 end
                307 % End of linear constraint consistency check
                308 
  0.09     860  309 LargeScaleFlag = strcmpi(optimget(options,'LargeScale',defaultopt,'fast'),'on');  
  0.07     860  310 Algorithm = optimget(options,'Algorithm',defaultopt,'fast');  
                311 % Option needed for processing initial guess
  0.02     860  312 AlwaysHonorConstraints = optimget(options,'AlwaysHonorConstraints',defaultopt,'fast');  
                313 
                314 % Read in and error check option TypicalX
  0.10     860  315 [typicalx,ME] = getNumericOrStringFieldValue('TypicalX','ones(numberOfVariables,1)', ... 
                316     ones(numberOfVariables,1),'a numeric value',options,defaultopt);
  0.02     860  317 if ~isempty(ME) 
                318     throw(ME)
                319 end
  0.05     860  320 checkoptionsize('TypicalX', size(typicalx), numberOfVariables); 
           860  321 chckdOpts.TypicalX = typicalx; 
                322 
                323 % Determine algorithm user chose via options. (We need this now
                324 % to set OUTPUT.algorithm in case of early termination due to 
                325 % inconsistent bounds.) 
                326 % This algorithm choice may be modified later when we check the 
                327 % problem type and supplied derivatives
           860  328 algChoiceOptsConflict = false; 
           860  329 if strcmpi(Algorithm,'active-set') 
                330     OUTPUT.algorithm = activeSet;
           860  331 elseif strcmpi(Algorithm,'sqp') 
                332     OUTPUT.algorithm = sqp;
  0.01     860  333 elseif strcmpi(Algorithm,'interior-point') 
           860  334     OUTPUT.algorithm = interiorPoint; 
                335 elseif strcmpi(Algorithm,'trust-region-reflective')
                336     if LargeScaleFlag
                337         OUTPUT.algorithm = trustRegionReflective;
                338     else
                339         % Conflicting options Algorithm='trust-region-reflective' and
                340         % LargeScale='off'. Choose active-set algorithm.
                341         algChoiceOptsConflict = true; % warn later, not in early termination
                342         OUTPUT.algorithm = activeSet;
                343     end
                344 else
                345     error('optim:fmincon:InvalidAlgorithm',...
                346         ['Invalid choice of option Algorithm for FMINCON. Choose ''sqp'',', ...
                347         ' ''interior-point'', ''trust-region-reflective'', or ''active-set''.']);
                348 end    
                349 
  0.15     860  350 [XOUT,l,u,msg] = checkbounds(XOUT,LB,UB,numberOfVariables); 
  0.02     860  351 if ~isempty(msg) 
                352    EXITFLAG = -2;
                353    [FVAL,LAMBDA,GRAD,HESSIAN] = deal([]);
                354    
                355    OUTPUT.iterations = 0;
                356    OUTPUT.funcCount = 0;
                357    OUTPUT.stepsize = [];
                358    if strcmpi(OUTPUT.algorithm,activeSet) || strcmpi(OUTPUT.algorithm,sqp)
                359        OUTPUT.lssteplength = [];
                360    else % trust-region-reflective, interior-point
                361        OUTPUT.cgiterations = [];
                362    end
                363    if strcmpi(OUTPUT.algorithm,interiorPoint) || strcmpi(OUTPUT.algorithm,activeSet) || ...
                364       strcmpi(OUTPUT.algorithm,sqp)
                365        OUTPUT.constrviolation = [];
                366    end
                367    OUTPUT.firstorderopt = [];
                368    OUTPUT.message = msg;
                369    
                370    X(:) = XOUT;
                371    if verbosity > 0
                372       disp(msg)
                373    end
                374    return
                375 end
           860  376 lFinite = l(~isinf(l)); 
           860  377 uFinite = u(~isinf(u)); 
                378 
                379 % Create structure of flags and initial values, initialize merit function
                380 % type and the original shape of X.
           860  381 flags.meritFunction = 0; 
  0.01     860  382 initVals.xOrigShape = X; 
                383 
  0.07     860  384 diagnostics = isequal(optimget(options,'Diagnostics',defaultopt,'fast'),'on'); 
  0.07     860  385 funValCheck = strcmpi(optimget(options,'FunValCheck',defaultopt,'fast'),'on'); 
  0.07     860  386 flags.grad = strcmpi(optimget(options,'GradObj',defaultopt,'fast'),'on'); 
                387 
                388 % Notice that defaultopt.Hessian = [], so the variable "hessian" can be empty
  0.05     860  389 hessian = optimget(options,'Hessian',defaultopt,'fast');  
                390 % If calling trust-region-reflective with an unavailable Hessian option value, 
                391 % issue informative error message
           860  392 if strcmpi(OUTPUT.algorithm,trustRegionReflective) && ... 
                393         ~( isempty(hessian) || strcmpi(hessian,'on') || strcmpi(hessian,'user-supplied') || ...
                394            strcmpi(hessian,'off') || strcmpi(hessian,'fin-diff-grads')  )
                395     error('optim:fmincon:BadTRReflectHessianValue', ...
                396         ['Value of Hessian option unavailable in trust-region-reflective algorithm. Possible\n' ...
                397          ' values are ''user-supplied'' and ''fin-diff-grads''.'])
                398 end
                399 
           860  400 if ~iscell(hessian) && ( strcmpi(hessian,'user-supplied') || strcmpi(hessian,'on') ) 
                401     flags.hess = true;
           860  402 else 
  0.01     860  403     flags.hess = false; 
           860  404 end 
                405 
           860  406 if isempty(NONLCON) 
           860  407    flags.constr = false; 
                408 else
                409    flags.constr = true;
                410 end
                411 
                412 % Process objective function
           860  413 if ~isempty(FUN)  % will detect empty string, empty matrix, empty cell array 
                414    % constrflag in optimfcnchk set to false because we're checking the objective, not constraint
  0.35     860  415    funfcn = optimfcnchk(FUN,'fmincon',length(varargin),funValCheck,flags.grad,flags.hess,false,Algorithm); 
                416 else
                417    error('optim:fmincon:InvalidFUN', ...
                418          ['FUN must be a function handle;\n', ...
                419           ' or, FUN may be a cell array that contains function handles.']);
                420 end
                421 
                422 % Process constraint function
           860  423 if flags.constr % NONLCON is non-empty 
                424    flags.gradconst = strcmpi(optimget(options,'GradConstr',defaultopt,'fast'),'on');
                425    % hessflag in optimfcnchk set to false because hessian is never returned by nonlinear constraint 
                426    % function
                427    %
                428    % constrflag in optimfcnchk set to true because we're checking the constraints
                429    confcn = optimfcnchk(NONLCON,'fmincon',length(varargin),funValCheck,flags.gradconst,false,true);
           860  430 else 
           860  431    flags.gradconst = false;  
           860  432    confcn = {'','','','',''}; 
           860  433 end 
                434 
           860  435 [rowAeq,colAeq] = size(Aeq); 
                436 
                437 % Look at problem type and supplied derivatives, and determine if need to
                438 % switch algorithm
  0.01     860  439 if strcmpi(OUTPUT.algorithm,activeSet) || strcmpi(OUTPUT.algorithm,sqp) 
                440     % Check if Algorithm was originally set to 'trust-region-reflective'
                441     % and then changed to 'active-set' because LargeScale was 'off'.
                442     % Do not warn if Algorithm was set to 'sqp' (defensive code).
                443     if algChoiceOptsConflict && strcmpi(OUTPUT.algorithm,activeSet)
                444         % Active-set algorithm chosen as a result of conflicting options
                445         warning('optim:fmincon:NLPAlgLargeScaleConflict', ...
                446             ['Options LargeScale = ''off'' and Algorithm = ''trust-region-reflective'' conflict.\n' ...
                447             'Ignoring Algorithm and running active-set algorithm. To run trust-region-reflective, set\n' ...
                448             'LargeScale = ''on''. To run active-set without this warning, use Algorithm = ''active-set''.'])
                449     end
                450     if issparse(Aeq) || issparse(A)
                451         warning('optim:fmincon:ConvertingToFull', ...
                452             'Cannot use sparse matrices with %s algorithm: converting to full.',Algorithm)
                453     end
                454     if flags.hess % conflicting options
                455         flags.hess = false;
                456         warning('optim:fmincon:HessianIgnored', ...
                457             ['An analytic Hessian cannot be used when OPTIONS.Algorithm = ''%s''.\n' ...
                458             'OPTIONS.Hessian will be ignored (user-supplied Hessian will not be used).'], Algorithm);
                459         if isequal(funfcn{1},'fungradhess')
                460             funfcn{1}='fungrad';
                461         elseif  isequal(funfcn{1},'fun_then_grad_then_hess')
                462             funfcn{1}='fun_then_grad';
                463         end
                464     end
  0.01     860  465 elseif strcmpi(OUTPUT.algorithm,trustRegionReflective) 
                466     if isempty(NONLCON) && isempty(A) && isempty(Aeq) && flags.grad
                467         % if only l and u then call sfminbx
                468     elseif isempty(NONLCON) && isempty(A) && isempty(lFinite) && isempty(uFinite) && flags.grad ...
                469             && colAeq > rowAeq
                470         % if only Aeq beq and Aeq has more columns than rows, then call sfminle
                471     else
                472         warning('optim:fmincon:SwitchingToMediumScale', ...
                473             ['Trust-region-reflective algorithm does not solve this type of problem, ' ...
                474             'using active-set algorithm. You could also try the interior-point or ' ...
                475             'sqp algorithms: set the Algorithm option to ''interior-point'' ', ...
                476             'or ''sqp'' and rerun. For more help, see %s in the documentation.'], ...
                477             addLink('Choosing the Algorithm','choose_algorithm'))
                478         if isequal(funfcn{1},'fungradhess')
                479             funfcn{1}='fungrad';
                480             warning('optim:fmincon:HessianIgnored', ...
                481                 ['Active-set algorithm is a Quasi-Newton method and does not use\n' ...
                482                 'analytic Hessian. Hessian flag in options will be ignored.'])
                483         elseif  isequal(funfcn{1},'fun_then_grad_then_hess')
                484             funfcn{1}='fun_then_grad';
                485             warning('optim:fmincon:HessianIgnored', ...
                486                 ['Active-set algorithm is a Quasi-Newton method and does not use\n' ...
                487                 'analytic Hessian. Hessian flag in options will be ignored.'])
                488         end
                489         flags.hess = false;
                490         OUTPUT.algorithm = activeSet; % switch to active-set
                491     end
                492 end
                493 
           860  494 lenvlb = length(l); 
           860  495 lenvub = length(u); 
                496 
           860  497 shiftedX0 = false; % boolean that indicates if initial point was shifted 
  0.01     860  498 if strcmpi(OUTPUT.algorithm,activeSet) 
                499    %
                500    % Ensure starting point lies within bounds
                501    %
                502    i=1:lenvlb;
                503    lindex = XOUT(i)<l(i);
                504    if any(lindex)
                505       XOUT(lindex)=l(lindex)+1e-4; 
                506       shiftedX0 = true;
                507    end
                508    i=1:lenvub;
                509    uindex = XOUT(i)>u(i);
                510    if any(uindex)
                511       XOUT(uindex)=u(uindex);
                512       shiftedX0 = true;
                513    end
                514    X(:) = XOUT;
           860  515 elseif strcmpi(OUTPUT.algorithm,trustRegionReflective) 
                516    %
                517    % If components of initial x not within bounds, set those components  
                518    % of initial point to a "box-centered" point
                519    %
                520    arg = (u >= 1e10); arg2 = (l <= -1e10);
                521    u(arg) = inf;
                522    l(arg2) = -inf;
                523    xinitOutOfBounds_idx = XOUT < l | XOUT > u;
                524    if any(xinitOutOfBounds_idx)
                525        shiftedX0 = true;
                526        XOUT = startx(u,l,XOUT,xinitOutOfBounds_idx);
                527        X(:) = XOUT;
                528    end
           860  529 elseif strcmpi(OUTPUT.algorithm,interiorPoint) 
                530     % Variables: fixed, finite lower bounds, finite upper bounds
  0.24     860  531     xIndices = classifyBoundsOnVars(l,u,numberOfVariables,true); 
                532 
                533     % If honor bounds mode, then check that initial point strictly satisfies the
                534     % simple inequality bounds on the variables and exactly satisfies fixed variable
                535     % bounds.
  0.02     860  536     if strcmpi(AlwaysHonorConstraints,'bounds') || strcmpi(AlwaysHonorConstraints,'bounds-ineqs') 
           860  537         violatedFixedBnds_idx = XOUT(xIndices.fixed) ~= l(xIndices.fixed); 
           860  538         violatedLowerBnds_idx = XOUT(xIndices.finiteLb) <= l(xIndices.finiteLb); 
           860  539         violatedUpperBnds_idx = XOUT(xIndices.finiteUb) >= u(xIndices.finiteUb); 
  0.01     860  540         if any(violatedLowerBnds_idx) || any(violatedUpperBnds_idx) || any(violatedFixedBnds_idx) 
  0.04     860  541             XOUT = shiftInitPtToInterior(numberOfVariables,XOUT,l,u,Inf); 
           860  542             X(:) = XOUT; 
           860  543             shiftedX0 = true; 
  0.01     860  544         end 
  0.01     860  545     end 
                546 else % SQP
                547     % Classify variables: finite lower bounds, finite upper bounds
                548     xIndices = classifyBoundsOnVars(l,u,numberOfVariables,false);
                549     
                550     % SQP always honors the bounds. Check that initial point
                551     % strictly satisfies the bounds on the variables.
                552     violatedLowerBnds_idx = XOUT(xIndices.finiteLb) < l(xIndices.finiteLb);
                553     violatedUpperBnds_idx = XOUT(xIndices.finiteUb) > u(xIndices.finiteUb);
                554     if any(violatedLowerBnds_idx) || any(violatedUpperBnds_idx)
                555         finiteLbIdx = find(xIndices.finiteLb);
                556         finiteUbIdx = find(xIndices.finiteUb);
                557         XOUT(finiteLbIdx(violatedLowerBnds_idx)) = l(finiteLbIdx(violatedLowerBnds_idx));
                558         XOUT(finiteUbIdx(violatedUpperBnds_idx)) = u(finiteUbIdx(violatedUpperBnds_idx));
                559         X(:) = XOUT;
                560         shiftedX0 = true;
                561     end
                562 end
                563 
                564 % Display that x0 was shifted in order to honor bounds
           860  565 if shiftedX0 
           860  566     if verbosity >= 3 
                567         if strcmpi(OUTPUT.algorithm,interiorPoint) 
                568             fprintf(['Your initial point x0 is not strictly between bounds lb and ub; FMINCON\n' ... 
                569                 ' shifted x0 to strictly satisfy the bounds.\n'])
                570             fprintf('\n');
                571         else
                572             fprintf(['Your initial point x0 is not between bounds lb and ub; FMINCON\n' ...
                573                 ' shifted x0 to satisfy the bounds.\n'])
                574             % Trust-region-reflective and active-set display a blank line right before 
                575             % iterative display; the newer algorithms, such as sqp, don't - we add one here.
                576             if strcmpi(OUTPUT.algorithm,sqp)  
                577                 fprintf('\n');
                578             end
                579         end
                580     end
           860  581 end 
                582     
                583 % Evaluate function
  0.01     860  584 initVals.g = zeros(numberOfVariables,1); 
           860  585 HESSIAN = []; 
                586 
           860  587 switch funfcn{1} 
           860  588 case 'fun' 
           860  589    try 
  0.10     860  590       initVals.f = feval(funfcn{3},X,varargin{:}); 
                591    catch userFcn_ME
                592         optim_ME = MException('optim:fmincon:ObjectiveError', ...
                593             'Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.');
                594         userFcn_ME = addCause(userFcn_ME,optim_ME);
                595         rethrow(userFcn_ME)
                596    end
                597 case 'fungrad'
                598    try
                599       [initVals.f,initVals.g(:)] = feval(funfcn{3},X,varargin{:});
                600    catch userFcn_ME
                601         optim_ME = MException('optim:fmincon:ObjectiveError', ...
                602             'Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.');
                603         userFcn_ME = addCause(userFcn_ME,optim_ME);
                604         rethrow(userFcn_ME)
                605    end
                606 case 'fungradhess'
                607    try
                608       [initVals.f,initVals.g(:),HESSIAN] = feval(funfcn{3},X,varargin{:});
                609    catch userFcn_ME
                610         optim_ME = MException('optim:fmincon:ObjectiveError', ...
                611             'Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.');
                612         userFcn_ME = addCause(userFcn_ME,optim_ME);
                613         rethrow(userFcn_ME)
                614    end
                615 case 'fun_then_grad'
                616    try
                617       initVals.f = feval(funfcn{3},X,varargin{:});
                618    catch userFcn_ME
                619         optim_ME = MException('optim:fmincon:ObjectiveError', ...
                620             'Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.');
                621         userFcn_ME = addCause(userFcn_ME,optim_ME);
                622         rethrow(userFcn_ME)
                623    end
                624    try
                625       initVals.g(:) = feval(funfcn{4},X,varargin{:});
                626    catch userFcn_ME
                627         optim_ME = MException('optim:fmincon:GradientError', ...
                628             'Failure in initial user-supplied objective gradient function evaluation. FMINCON cannot continue.');
                629         userFcn_ME = addCause(userFcn_ME,optim_ME);
                630         rethrow(userFcn_ME)
                631    end
                632 case 'fun_then_grad_then_hess'
                633    try
                634       initVals.f = feval(funfcn{3},X,varargin{:});
                635    catch userFcn_ME
                636         optim_ME = MException('optim:fmincon:ObjectiveError', ...
                637             'Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.');
                638         userFcn_ME = addCause(userFcn_ME,optim_ME);
                639         rethrow(userFcn_ME)
                640    end
                641    try
                642       initVals.g(:) = feval(funfcn{4},X,varargin{:});
                643    catch userFcn_ME
                644         optim_ME = MException('optim:fmincon:GradientError', ...
                645             'Failure in initial user-supplied objective gradient function evaluation. FMINCON cannot continue.');
                646         userFcn_ME = addCause(userFcn_ME,optim_ME);
                647         rethrow(userFcn_ME)
                648    end
                649    try
                650       HESSIAN = feval(funfcn{5},X,varargin{:});
                651    catch userFcn_ME
                652         optim_ME = MException('optim:fmincon:HessianError', ...
                653             'Failure in initial user-supplied objective Hessian function evaluation. FMINCON cannot continue.');            
                654         userFcn_ME = addCause(userFcn_ME,optim_ME);
                655         rethrow(userFcn_ME)
                656    end
                657 otherwise
                658    error('optim:fmincon:UndefinedCallType','Undefined calltype in FMINCON.');
                659 end
                660 
                661 % Check that the objective value is a scalar
           860  662 if numel(initVals.f) ~= 1 
                663    error('optim:fmincon:NonScalarObj','User supplied objective function must return a scalar value.')
                664 end
                665 
                666 % Evaluate constraints
  0.01     860  667 switch confcn{1} 
           860  668 case 'fun' 
                669     try
                670         [ctmp,ceqtmp] = feval(confcn{3},X,varargin{:});
                671     catch userFcn_ME
                672         if strcmpi('MATLAB:maxlhs',userFcn_ME.identifier)
                673                 error('optim:fmincon:InvalidHandleNonlcon', ...
                674                     ['The constraint function must return two outputs; ' ... 
                675                      'the nonlinear inequality constraints and ' ...
                676                      'the nonlinear equality constraints.'])
                677         else
                678             optim_ME = MException('optim:fmincon:NonlconError', ...
                679                 'Failure in initial user-supplied nonlinear constraint function evaluation. FMINCON cannot continue.');
                680             userFcn_ME = addCause(userFcn_ME,optim_ME);
                681             rethrow(userFcn_ME)
                682         end
                683     end
                684     initVals.ncineq = ctmp(:);
                685     initVals.nceq = ceqtmp(:);
                686     initVals.gnc = zeros(numberOfVariables,length(initVals.ncineq));
                687     initVals.gnceq = zeros(numberOfVariables,length(initVals.nceq));
           860  688 case 'fungrad' 
                689    try
                690       [ctmp,ceqtmp,initVals.gnc,initVals.gnceq] = feval(confcn{3},X,varargin{:});
                691    catch userFcn_ME
                692        optim_ME = MException('optim:fmincon:NonlconError', ...
                693            'Failure in initial user-supplied nonlinear constraint function evaluation. FMINCON cannot continue.');           
                694        userFcn_ME = addCause(userFcn_ME,optim_ME);
                695        rethrow(userFcn_ME)
                696    end
                697    initVals.ncineq = ctmp(:);
                698    initVals.nceq = ceqtmp(:);
  0.01     860  699 case 'fun_then_grad' 
                700     try
                701         [ctmp,ceqtmp] = feval(confcn{3},X,varargin{:});
                702     catch userFcn_ME
                703         optim_ME = MException('optim:fmincon:NonlconError', ...
                704             'Failure in initial user-supplied nonlinear constraint function evaluation. FMINCON cannot continue.');
                705         userFcn_ME = addCause(userFcn_ME,optim_ME);
                706         rethrow(userFcn_ME)
                707     end
                708     initVals.ncineq = ctmp(:);
                709     initVals.nceq = ceqtmp(:);
                710     try
                711         [initVals.gnc,initVals.gnceq] = feval(confcn{4},X,varargin{:});
                712     catch userFcn_ME
                713         optim_ME = MException('optim:fmincon:NonlconFunOrGradError', ...
                714             'Failure in initial user-supplied nonlinear constraint gradinet function evaluation. FMINCON cannot continue.');
                715         userFcn_ME = addCause(userFcn_ME,optim_ME);
                716         rethrow(userFcn_ME)
                717     end
  0.01     860  718 case '' 
                719    % No nonlinear constraints. Reshaping of empty quantities is done later
                720    % in this file, where both cases, (i) no nonlinear constraints and (ii)
                721    % nonlinear constraints that have one type missing (equalities or
                722    % inequalities), are handled in one place
           860  723    initVals.ncineq = []; 
           860  724    initVals.nceq = []; 
           860  725    initVals.gnc = []; 
  0.01     860  726    initVals.gnceq = []; 
                727 otherwise
                728    error('optim:fmincon:UndefinedCalltype','Undefined calltype in FMINCON.');
                729 end
                730 
           860  731 non_eq = length(initVals.nceq); 
           860  732 non_ineq = length(initVals.ncineq); 
                733 
                734 % Make sure empty constraint and their derivatives have correct sizes (not 0-by-0):
           860  735 if isempty(initVals.ncineq) 
           860  736     initVals.ncineq = reshape(initVals.ncineq,0,1); 
           860  737 end 
           860  738 if isempty(initVals.nceq) 
           860  739     initVals.nceq = reshape(initVals.nceq,0,1); 
           860  740 end 
           860  741 if isempty(initVals.gnc) 
           860  742     initVals.gnc = reshape(initVals.gnc,numberOfVariables,0); 
           860  743 end 
           860  744 if isempty(initVals.gnceq) 
           860  745     initVals.gnceq = reshape(initVals.gnceq,numberOfVariables,0); 
  0.01     860  746 end 
           860  747 [cgrow,cgcol] = size(initVals.gnc); 
           860  748 [ceqgrow,ceqgcol] = size(initVals.gnceq); 
                749 
           860  750 if cgrow ~= numberOfVariables || cgcol ~= non_ineq 
                751    error('optim:fmincon:WrongSizeGradNonlinIneq', ...
                752          'Gradient of nonlinear inequality constraints must have size %i-by-%i.', ...
                753          numberOfVariables,non_ineq)
                754 end
  0.01     860  755 if ceqgrow ~= numberOfVariables || ceqgcol ~= non_eq 
                756    error('optim:fmincon:WrongSizeGradNonlinEq', ...
                757          'Gradient of nonlinear equality constraints must have size %i-by-%i.', ...
                758          numberOfVariables,non_eq)
                759 end
                760 
           860  761 if diagnostics > 0 
                762    % Do diagnostics on information so far
                763    diagnose('fmincon',OUTPUT,flags.grad,flags.hess,flags.constr,flags.gradconst,...
                764       ~LargeScaleFlag,options,defaultopt,XOUT,non_eq,...
                765       non_ineq,lin_eq,lin_ineq,l,u,funfcn,confcn,initVals.f,initVals.g,HESSIAN, ...
                766       initVals.ncineq,initVals.nceq,initVals.gnc,initVals.gnceq);
                767 end
                768 
                769 % call algorithm
  0.02     860  770 if strcmpi(OUTPUT.algorithm,activeSet) % active-set 
                771     defaultopt.MaxIter = 400; defaultopt.MaxFunEvals = '100*numberofvariables'; defaultopt.TolX = 1e-6;
                772     defaultopt.Hessian = 'off';
                773     problemInfo = []; % No problem related data
                774     [X,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN]=...
                775         nlconst(funfcn,X,l,u,full(A),B,full(Aeq),Beq,confcn,options,defaultopt, ...
                776         chckdOpts,verbosity,flags,initVals,problemInfo,varargin{:});
           860  777 elseif strcmpi(OUTPUT.algorithm,trustRegionReflective) % trust-region-reflective 
                778    if (isequal(funfcn{1}, 'fun_then_grad_then_hess') || isequal(funfcn{1}, 'fungradhess'))
                779       Hstr = [];
                780    elseif (isequal(funfcn{1}, 'fun_then_grad') || isequal(funfcn{1}, 'fungrad'))
                781       n = length(XOUT); 
                782       Hstr = optimget(options,'HessPattern',defaultopt,'fast');
                783       if ischar(Hstr) 
                784          if isequal(lower(Hstr),'sparse(ones(numberofvariables))')
                785             Hstr = sparse(ones(n));
                786          else
                787             error('optim:fmincon:InvalidHessPattern', ...
                788                   'Option ''HessPattern'' must be a matrix if not the default.')
                789          end
                790       end
                791       checkoptionsize('HessPattern', size(Hstr), n);
                792    end
                793    
                794    defaultopt.MaxIter = 400; defaultopt.MaxFunEvals = '100*numberofvariables'; defaultopt.TolX = 1e-6;
                795    defaultopt.Hessian = 'off';
                796    % Trust-region-reflective algorithm does not compute constraint
                797    % violation as it progresses. If the user requests the output structure,
                798    % we need to calculate the constraint violation at the returned
                799    % solution.
                800    if nargout > 3
                801        computeConstrViolForOutput = true;
                802    else
                803        computeConstrViolForOutput = false;
                804    end
                805    if isempty(Aeq)
                806       [X,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = ...
                807          sfminbx(funfcn,X,l,u,verbosity,options,defaultopt,computeLambda,initVals.f,initVals.g, ...
                808          HESSIAN,Hstr,flags.detailedExitMsg,computeConstrViolForOutput,varargin{:});
                809    else
                810       [X,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = ...
                811          sfminle(funfcn,X,sparse(Aeq),Beq,verbosity,options,defaultopt,computeLambda,initVals.f, ...
                812          initVals.g,HESSIAN,Hstr,flags.detailedExitMsg,computeConstrViolForOutput,varargin{:});
                813    end
           860  814 elseif strcmpi(OUTPUT.algorithm,interiorPoint) 
           860  815     defaultopt.MaxIter = 1000; defaultopt.MaxFunEvals = 3000; defaultopt.TolX = 1e-10; 
           860  816     defaultopt.Hessian = 'bfgs'; 
           860  817     mEq = lin_eq + non_eq + nnz(xIndices.fixed); % number of equalities 
                818     % Interior-point-specific options. Default values for lbfgs memory is 10, and 
                819     % ldl pivot threshold is 0.01
  2.91     860  820     [options_ip,optionFeedback] = getIpOptions(options,numberOfVariables,mEq,flags.constr,defaultopt,10,0.01);  
                821 
 206.33     860  822     [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options_ip.HessFcn, ... 
                823         initVals.f,initVals.g,initVals.ncineq,initVals.nceq,initVals.gnc,initVals.gnceq,HESSIAN, ...
                824         xIndices,options_ip,optionFeedback,varargin{:});
                825 else % sqp
                826     defaultopt.MaxIter = 400; defaultopt.MaxFunEvals = '100*numberofvariables'; 
                827     defaultopt.TolX = 1e-6; defaultopt.Hessian = 'bfgs';
                828     % Validate options used by sqp
                829     [options,optionFeedback] = getSQPOptions(options,defaultopt,numberOfVariables);
                830     optionFeedback.detailedExitMsg = flags.detailedExitMsg;
                831     % Call algorithm
                832     [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = sqpLineSearch(funfcn,X,full(A),full(B),full(Aeq),full(Beq), ...
                833         full(l),full(u),confcn,initVals.f,full(initVals.g),full(initVals.ncineq),full(initVals.nceq), ...
                834         full(initVals.gnc),full(initVals.gnceq),xIndices,options,verbosity,optionFeedback,varargin{:});
                835 end