function [x,fval,exitflag,output] = fminsearch2(funfcn,x,options,varargin)
%FMINSEARCH Multidimensional unconstrained nonlinear minimization (Nelder-Mead).
% X = FMINSEARCH(FUN,X0) starts at X0 and attempts to find a local minimizer
% X of the function FUN. FUN is a function handle. FUN accepts input X and
% returns a scalar function value F evaluated at X. X0 can be a scalar, vector
% or matrix.
%
% X = FMINSEARCH(FUN,X0,OPTIONS) minimizes with the default optimization
% parameters replaced by values in the structure OPTIONS, created
% with the OPTIMSET function. See OPTIMSET for details. FMINSEARCH uses
% these options: Display, TolX, TolFun, MaxFunEvals, MaxIter, FunValCheck,
% PlotFcns, and OutputFcn.
%
% X = FMINSEARCH(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a
% structure with the function FUN in PROBLEM.objective, the start point
% in PROBLEM.x0, the options structure in PROBLEM.options, and solver
% name 'fminsearch' in PROBLEM.solver. The PROBLEM structure must have
% all the fields.
%
% [X,FVAL]= FMINSEARCH(...) returns the value of the objective function,
% described in FUN, at X.
%
% [X,FVAL,EXITFLAG] = FMINSEARCH(...) returns an EXITFLAG that describes
% the exit condition of FMINSEARCH. Possible values of EXITFLAG and the
% corresponding exit conditions are
%
% 1 Maximum coordinate difference between current best point and other
% points in simplex is less than or equal to TolX, and corresponding
% difference in function values is less than or equal to TolFun.
% 0 Maximum number of function evaluations or iterations reached.
% -1 Algorithm terminated by the output function.
%
% [X,FVAL,EXITFLAG,OUTPUT] = FMINSEARCH(...) returns a structure
% OUTPUT with the number of iterations taken in OUTPUT.iterations, the
% number of function evaluations in OUTPUT.funcCount, the algorithm name
% in OUTPUT.algorithm, and the exit message in OUTPUT.message.
%
% Examples
% FUN can be specified using @:
% X = fminsearch(@sin,3)
% finds a minimum of the SIN function near 3.
% In this case, SIN is a function that returns a scalar function value
% SIN evaluated at X.
%
% FUN can also be an anonymous function:
% X = fminsearch(@(x) norm(x),[1;2;3])
% returns a point near the minimizer [0;0;0].
%
% If FUN is parameterized, you can use anonymous functions to capture the
% problem-dependent parameters. Suppose you want to optimize the objective
% given in the function myfun, which is parameterized by its second argument c.
% Here myfun is an M-file function such as
%
% function f = myfun(x,c)
% f = x(1)^2 + c*x(2)^2;
%
% To optimize for a specific value of c, first assign the value to c. Then
% create a one-argument anonymous function that captures that value of c
% and calls myfun with two arguments. Finally, pass this anonymous function
% to FMINSEARCH:
%
% c = 1.5; % define parameter first
% x = fminsearch(@(x) myfun(x,c),[0.3;1])
%
% FMINSEARCH uses the Nelder-Mead simplex (direct search) method.
%
% See also OPTIMSET, FMINBND, FUNCTION_HANDLE.
% Reference: Jeffrey C. Lagarias, James A. Reeds, Margaret H. Wright,
% Paul E. Wright, "Convergence Properties of the Nelder-Mead Simplex
% Method in Low Dimensions", SIAM Journal of Optimization, 9(1):
% p.112-147, 1998.
% Copyright 1984-2006 The MathWorks, Inc.
% $Revision: 1.21.4.14 $ $Date: 2006/12/15 19:27:41 $
% Modified initial run choice to improve performance. Based on
% adaptive-one-factor-at-a-time experimental design, to aim the simplex
% in the best direction.
% crfoster@mit.edu
defaultopt = struct('Display','notify','MaxIter','200*numberOfVariables',...
'MaxFunEvals','200*numberOfVariables','TolX',1e-4,'TolFun',1e-4, ...
'FunValCheck','off','OutputFcn',[],'PlotFcns',[]);
% If just 'defaults' passed in, return the default options in X
if nargin==1 && nargout <= 1 && isequal(funfcn,'defaults')
x = defaultopt;
return
end
if nargin<3, options = []; end
% Detect problem structure input
if nargin == 1
if isa(funfcn,'struct')
[funfcn,x,options] = separateOptimStruct(funfcn);
else % Single input and non-structure
error('MATLAB:fminsearch:InputArg','The input to FMINSEARCH should be either a structure with valid fields or consist of at least two arguments.');
end
end
if nargin == 0
error('MATLAB:fminsearch:NotEnoughInputs',...
'FMINSEARCH requires at least two input arguments');
end
% Check for non-double inputs
if ~isa(x,'double')
error('MATLAB:fminsearch:NonDoubleInput', ...
'FMINSEARCH only accepts inputs of data type double.')
end
n = numel(x);
numberOfVariables = n;
printtype = optimget(options,'Display',defaultopt,'fast');
tolx = optimget(options,'TolX',defaultopt,'fast');
tolf = optimget(options,'TolFun',defaultopt,'fast');
maxfun = optimget(options,'MaxFunEvals',defaultopt,'fast');
maxiter = optimget(options,'MaxIter',defaultopt,'fast');
funValCheck = strcmp(optimget(options,'FunValCheck',defaultopt,'fast'),'on');
% In case the defaults were gathered from calling: optimset('fminsearch'):
if ischar(maxfun)
if isequal(lower(maxfun),'200*numberofvariables')
maxfun = 200*numberOfVariables;
else
error('MATLAB:fminsearch:OptMaxFunEvalsNotInteger',...
'Option ''MaxFunEvals'' must be an integer value if not the default.')
end
end
if ischar(maxiter)
if isequal(lower(maxiter),'200*numberofvariables')
maxiter = 200*numberOfVariables;
else
error('MATLAB:fminsearch:OptMaxIterNotInteger',...
'Option ''MaxIter'' must be an integer value if not the default.')
end
end
switch printtype
case 'notify'
prnt = 1;
case {'none','off'}
prnt = 0;
case 'iter'
prnt = 3;
case 'final'
prnt = 2;
case 'simplex'
prnt = 4;
otherwise
prnt = 1;
end
% Handle the output
outputfcn = optimget(options,'OutputFcn',defaultopt,'fast');
if isempty(outputfcn)
haveoutputfcn = false;
else
haveoutputfcn = true;
xOutputfcn = x; % Last x passed to outputfcn; has the input x's shape
% Parse OutputFcn which is needed to support cell array syntax for OutputFcn.
outputfcn = createCellArrayOfFunctions(outputfcn,'OutputFcn');
end
% Handle the plot
plotfcns = optimget(options,'PlotFcns',defaultopt,'fast');
if isempty(plotfcns)
haveplotfcn = false;
else
haveplotfcn = true;
xOutputfcn = x; % Last x passed to plotfcns; has the input x's shape
% Parse PlotFcns which is needed to support cell array syntax for PlotFcns.
plotfcns = createCellArrayOfFunctions(plotfcns,'PlotFcns');
end
header = ' Iteration Func-count min f(x) Procedure';
% Convert to function handle as needed.
funfcn = fcnchk(funfcn,length(varargin));
% Add a wrapper function to check for Inf/NaN/complex values
if funValCheck
% Add a wrapper function, CHECKFUN, to check for NaN/complex values without
% having to change the calls that look like this:
% f = funfcn(x,varargin{:});
% x is the first argument to CHECKFUN, then the user's function,
% then the elements of varargin. To accomplish this we need to add the
% user's function to the beginning of varargin, and change funfcn to be
% CHECKFUN.
varargin = {funfcn, varargin{:}};
funfcn = @checkfun;
end
n = numel(x);
% Initialize parameters
rho = 1; chi = 2; psi = 0.5; sigma = 0.5;
onesn = ones(1,n);
two2np1 = 2:n+1;
one2n = 1:n;
% Set up a simplex near the initial guess.
xin = x(:); % Force xin to be a column vector
v = zeros(n,n+1); fv = zeros(1,n+1);
v(:,1) = xin; % Place input guess in the simplex! (credit L.Pfeffer at Stanford)
x(:) = xin; % Change x to the form expected by funfcn
fv(:,1) = funfcn(x,varargin{:});
func_evals = 1;
itercount = 0;
how = '';
% Initial simplex setup continues later
% Initialize the output and plot functions.
if haveoutputfcn || haveplotfcn
[xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,v(:,1),xOutputfcn,'init',itercount, ...
func_evals, how, fv(:,1),varargin{:});
if stop
[x,fval,exitflag,output] = cleanUpInterrupt(xOutputfcn,optimValues);
if prnt > 0
disp(output.message)
end
return;
end
end
% Print out initial f(x) as 0th iteration
if prnt == 3
disp(' ')
disp(header)
disp(sprintf(' %5.0f %5.0f %12.6g %s', itercount, func_evals, fv(1), how));
elseif prnt == 4
clc
formatsave = get(0,{'format','formatspacing'});
format compact
format short e
disp(' ')
disp(how)
v
fv
func_evals
end
% OutputFcn and PlotFcns call
if haveoutputfcn || haveplotfcn
[xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,v(:,1),xOutputfcn,'iter',itercount, ...
func_evals, how, fv(:,1),varargin{:});
if stop % Stop per user request.
[x,fval,exitflag,output] = cleanUpInterrupt(xOutputfcn,optimValues);
if prnt > 0
disp(output.message)
end
return;
end
end
% Continue setting up the initial simplex.
% The inital improment suggested L.Pfeffer L.Pfeffer at Stanford
% includede setting 5 percent deltas for non-zero terms and 0.00025
% for zero elements of x
% This has been modified by crfoster@mit.edu to include the number of
% variables in setting up the deltas. This is based on an
% adaptive-one-factor-at-a-time experiment to set up the initial simplex.
% The following deltas try to maintain a similar distance between the
% data points and the centroid of the design. For n=2 the values are
% identical to L.Pfeffer, but they increase to .2 as the number of
% variables increase
usual_delta = .2*.5*(n-n*.5+.5^n-1)/(n*(.5-1)^2);
zero_term_delta = .001*usual_delta/.2;
% This performs an adaptive-one-factor-at-a-time experiment based on the
% work of Dan Frey at MIT. The goal is to 'aim' the simplex in the right
% direction with these n+1 values. Thus it begins at the given starting
% point and walks towards the higher values.
% This change was run on the 35 routines provided by
% More, Garbow, Hilltrom's "Testing Unconstrained Optimization", ACM
% Transactions on Mathematical Software, Vol 7, No.1 March 1981 p17-41
% It solved one problem than the native could not, and converged on two
% others where the native could not. The change is inexpensive in terms of
% cost of runs and performs remarkably better than the original.
y = xin;
if y(1) ~= 0
y(1) = (1 + usual_delta)*y(1);
else
y(1) = zero_term_delta;
end
v(:,2) = y;
x(:) = y; f = funfcn(x,varargin{:});
fv(1,2) = f;
oldloc = 1;
for j = 2:n
if fv(1,oldloc)<=fv(1,j)
y = v(:,oldloc);
else
y = v(:,j);
oldloc=j;
end
if y(j) ~= 0
y(j) = (1 + usual_delta)*y(j);
else
y(j) = zero_term_delta;
end
v(:,j+1) = y;
x(:) = y; f = funfcn(x,varargin{:});
fv(1,j+1) = f;
end
% sort so v(1,:) has the lowest function value
[fv,j] = sort(fv);
v = v(:,j);
how = 'initial simplex';
itercount = itercount + 1;
func_evals = n+1;
if prnt == 3
disp(sprintf(' %5.0f %5.0f %12.6g %s', itercount, func_evals, fv(1), how))
elseif prnt == 4
disp(' ')
disp(how)
v
fv
func_evals
end
% OutputFcn and PlotFcns call
if haveoutputfcn || haveplotfcn
[xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,v(:,1),xOutputfcn,'iter',itercount, ...
func_evals, how, fv(:,1),varargin{:});
if stop % Stop per user request.
[x,fval,exitflag,output] = cleanUpInterrupt(xOutputfcn,optimValues);
if prnt > 0
disp(output.message)
end
return;
end
end
exitflag = 1;
% Main algorithm: iterate until
% (a) the maximum coordinate difference between the current best point and the
% other points in the simplex is less than or equal to TolX. Specifically,
% until max(||v2-v1||,||v2-v1||,...,||v(n+1)-v1||) <= TolX,
% where ||.|| is the infinity-norm, and v1 holds the
% vertex with the current lowest value; AND
% (b) the corresponding difference in function values is less than or equal
% to TolFun. (Cannot use OR instead of AND.)
% The iteration stops if the maximum number of iterations or function evaluations
% are exceeded
while func_evals < maxfun && itercount < maxiter
if max(abs(fv(1)-fv(two2np1))) <= tolf && ...
max(max(abs(v(:,two2np1)-v(:,onesn)))) <= tolx
break
end
% Compute the reflection point
% xbar = average of the n (NOT n+1) best points
xbar = sum(v(:,one2n), 2)/n;
xr = (1 + rho)*xbar - rho*v(:,end);
x(:) = xr; fxr = funfcn(x,varargin{:});
func_evals = func_evals+1;
if fxr < fv(:,1)
% Calculate the expansion point
xe = (1 + rho*chi)*xbar - rho*chi*v(:,end);
x(:) = xe; fxe = funfcn(x,varargin{:});
func_evals = func_evals+1;
if fxe < fxr
v(:,end) = xe;
fv(:,end) = fxe;
how = 'expand';
else
v(:,end) = xr;
fv(:,end) = fxr;
how = 'reflect';
end
else % fv(:,1) <= fxr
if fxr < fv(:,n)
v(:,end) = xr;
fv(:,end) = fxr;
how = 'reflect';
else % fxr >= fv(:,n)
% Perform contraction
if fxr < fv(:,end)
% Perform an outside contraction
xc = (1 + psi*rho)*xbar - psi*rho*v(:,end);
x(:) = xc; fxc = funfcn(x,varargin{:});
func_evals = func_evals+1;
if fxc <= fxr
v(:,end) = xc;
fv(:,end) = fxc;
how = 'contract outside';
else
% perform a shrink
how = 'shrink';
end
else
% Perform an inside contraction
xcc = (1-psi)*xbar + psi*v(:,end);
x(:) = xcc; fxcc = funfcn(x,varargin{:});
func_evals = func_evals+1;
if fxcc < fv(:,end)
v(:,end) = xcc;
fv(:,end) = fxcc;
how = 'contract inside';
else
% perform a shrink
how = 'shrink';
end
end
if strcmp(how,'shrink')
for j=two2np1
v(:,j)=v(:,1)+sigma*(v(:,j) - v(:,1));
x(:) = v(:,j); fv(:,j) = funfcn(x,varargin{:});
end
func_evals = func_evals + n;
end
end
end
[fv,j] = sort(fv);
v = v(:,j);
itercount = itercount + 1;
if prnt == 3
disp(sprintf(' %5.0f %5.0f %12.6g %s', itercount, func_evals, fv(1), how))
elseif prnt == 4
disp(' ')
disp(how)
v
fv
func_evals
end
% OutputFcn and PlotFcns call
if haveoutputfcn || haveplotfcn
[xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,v(:,1),xOutputfcn,'iter',itercount, ...
func_evals, how, fv(:,1),varargin{:});
if stop % Stop per user request.
[x,fval,exitflag,output] = cleanUpInterrupt(xOutputfcn,optimValues);
if prnt > 0
disp(output.message)
end
return;
end
end
end % while
x(:) = v(:,1);
fval = fv(:,1);
if prnt == 4,
% reset format
set(0,{'format','formatspacing'},formatsave);
end
output.iterations = itercount;
output.funcCount = func_evals;
output.algorithm = 'Nelder-Mead simplex direct search';
% OutputFcn and PlotFcns call
if haveoutputfcn || haveplotfcn
callOutputAndPlotFcns(outputfcn,plotfcns,x,xOutputfcn,'done',itercount, func_evals, how, fval, varargin{:});
end
if func_evals >= maxfun
msg = sprintf(['Exiting: Maximum number of function evaluations has been exceeded\n' ...
' - increase MaxFunEvals option.\n' ...
' Current function value: %f \n'], fval);
if prnt > 0
disp(' ')
disp(msg)
end
exitflag = 0;
elseif itercount >= maxiter
msg = sprintf(['Exiting: Maximum number of iterations has been exceeded\n' ...
' - increase MaxIter option.\n' ...
' Current function value: %f \n'], fval);
if prnt > 0
disp(' ')
disp(msg)
end
exitflag = 0;
else
msg = ...
sprintf(['Optimization terminated:\n', ...
' the current x satisfies the termination criteria using OPTIONS.TolX of %e \n' ...
' and F(X) satisfies the convergence criteria using OPTIONS.TolFun of %e \n'], ...
tolx, tolf);
if prnt > 1
disp(' ')
disp(msg)
end
exitflag = 1;
end
output.message = msg;
%--------------------------------------------------------------------------
function [xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,x,xOutputfcn,state,iter,...
numf,how,f,varargin)
% CALLOUTPUTANDPLOTFCNS assigns values to the struct OptimValues and then calls the
% outputfcn/plotfcns.
%
% state - can have the values 'init','iter', or 'done'.
% For the 'done' state we do not check the value of 'stop' because the
% optimization is already done.
optimValues.iteration = iter;
optimValues.funccount = numf;
optimValues.fval = f;
optimValues.procedure = how;
xOutputfcn(:) = x; % Set x to have user expected size
stop = false;
% Call output functions
if ~isempty(outputfcn)
switch state
case {'iter','init'}
stop = callAllOptimOutputFcns(outputfcn,xOutputfcn,optimValues,state,varargin{:}) || stop;
case 'done'
callAllOptimOutputFcns(outputfcn,xOutputfcn,optimValues,state,varargin{:});
otherwise
error('MATLAB:fminsearch:InvalidState', ...
'Unknown state in CALLOUTPUTANDPLOTFCNS.')
end
end
% Call plot functions
if ~isempty(plotfcns)
switch state
case {'iter','init'}
stop = callAllOptimPlotFcns(plotfcns,xOutputfcn,optimValues,state,varargin{:}) || stop;
case 'done'
callAllOptimPlotFcns(plotfcns,xOutputfcn,optimValues,state,varargin{:});
otherwise
error('MATLAB:fminsearch:InvalidState', ...
'Unknown state in CALLOUTPUTANDPLOTFCNS.')
end
end
%--------------------------------------------------------------------------
function [x,FVAL,EXITFLAG,OUTPUT] = cleanUpInterrupt(xOutputfcn,optimValues)
% CLEANUPINTERRUPT updates or sets all the output arguments of FMINBND when the optimization
% is interrupted.
x = xOutputfcn;
FVAL = optimValues.fval;
EXITFLAG = -1;
OUTPUT.iterations = optimValues.iteration;
OUTPUT.funcCount = optimValues.funccount;
OUTPUT.algorithm = 'golden section search, parabolic interpolation';
OUTPUT.message = 'Optimization terminated prematurely by user.';
%--------------------------------------------------------------------------
function f = checkfun(x,userfcn,varargin)
% CHECKFUN checks for complex or NaN results from userfcn.
f = userfcn(x,varargin{:});
% Note: we do not check for Inf as FMINSEARCH handles it naturally.
if isnan(f)
error('MATLAB:fminsearch:checkfun:NaNFval', ...
'User function ''%s'' returned NaN when evaluated;\n FMINSEARCH cannot continue.', ...
localChar(userfcn));
elseif ~isreal(f)
error('MATLAB:fminsearch:checkfun:ComplexFval', ...
'User function ''%s'' returned a complex value when evaluated;\n FMINSEARCH cannot continue.', ...
localChar(userfcn));
end
%--------------------------------------------------------------------------
function strfcn = localChar(fcn)
% Convert the fcn to a string for printing
if ischar(fcn)
strfcn = fcn;
elseif isa(fcn,'inline')
strfcn = char(fcn);
elseif isa(fcn,'function_handle')
strfcn = func2str(fcn);
else
try
strfcn = char(fcn);
catch
strfcn = '(name not printable)';
end
end
