# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[xmin,fmin,m,ntype]=nelmed(...
```function [xmin,fmin,m,ntype]=nelmed(...
F,x0,dx,epsx,epsf,M,ifpr,varargin)
% [xmin,fmin,m,ntype]=nelmed(...
%             F,x0,dx,epsx,epsf,M,ifpr,varargin)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function performs multidimensional
% unconstrained function minimization using the
% direct search procedure developed by
% J. A. Nelder and R. Mead. The method is
% described in various books such as:
% 'Nonlinear Optimization', by M. Avriel
%
% F        - objective function of the form
%            F(x,p1,p2,...) where x is vector
%            in n space and p1,p2,... are any
%            auxiliary parameters needed to
%            define F
% x0       - starting vector to initiate
%            the search
% dx       - initial polyhedron side length
% epsx     - convergence tolerance on x
% epsf     - convergence tolerance on
%            function values
% M        - function evaluation limit to
%            terminate search
% ifpr     - when this parameter equals one,
%            different stages in the search
%            are printed
% varargin - variable length list of parameters
%            which can be passed to function F
% xmin     - coordinates of the smallest
%            function value
% fmin     - smallest function value found
% m        - total number of function
% ntype    - a vector containing
%            [ninit,nrefl,nexpn,ncontr,nshrnk]
%            which tells the number of reflect-
%            ions, expansions, contractions,and
%            shrinkages performed
%
% User m functions called: objective function
%                  named in the argument list

if isempty(ifpr), ifpr=0; end
if isempty(M), M=500; end;
if isempty(epsf), epsf=1e-5; end
if isempty(epsx), epsx=1e-5; end

% Initialize the simplex array
x0=x0(:); n=length(x0); N=n+1; f=zeros(1,N);
x=repmat(x0,1,N)+[zeros(n,1),dx*eye(n,n)];
for k=1:N
f(k)=feval(F,x(:,k),varargin{:});
end

ninit=N; nrefl=0; nexpn=0; ncontr=0;
nshrnk=0; m=N;

Erx=realmax; Erf=realmax;
alpha=1.0; % Reflection coefficient
beta= 0.5; % Contraction coefficient
gamma=2.0; % Expansion coefficient

% Top of the minimization loop

while Erx>epsx | Erf>epsf

[f,k]=sort(f); x=x(:,k);

% Exit if maximum allowable number of
%function values is exceeded
if m>M, xmin=x(:,1); fmin=f(1); return; end

% Generate the reflected point and
% function value
c=sum(x(:,1:n),2)/n; xr=c+alpha*(c-x(:,N));
fr=feval(F,xr,varargin{:}); m=m+1;
nrefl=nrefl+1;
if ifpr==1, fprintf(' :RFL \n'); end

if fr<f(1)
% Expand and take best from expansion
% or reflection
xe=c+gamma*(xr-c);
fe=feval(F,xe,varargin{:});
m=m+1; nexpn=nexpn+1;
if ifpr==1, fprintf(' :EXP \n'); end

if fr<fe
% The reflected point was best
f(N)=fr; x(:,N)=xr;
else
% The expanded point was best
f(N)=fe; x(:,N)=xe;
end

elseif fr<=f(n)  % In the middle zone
f(N)=fr; x(:,N)=xr;

else
% Reflected point exceeds second the
% highest value so either use contraction
% or shrinkage
if fr<f(N)
xx=xr; ff=fr;
else
xx=x(:,N); ff=f(N);
end

xc=c+beta*(xx-c);
fc=feval(F,xc,varargin{:});
m=m+1; ncontr=ncontr+1;

if fc<=ff
% Accept the contracted value
x(:,N)=xc; f(N)=fc;
if ifpr==1, fprintf(' :CNT \n'); end

else
% Shrink the simplex toward
% the best point
x=(x+repmat(x(:,1),1,N))/2;
for j=2:N
f(j)=feval(F,x(:,j),varargin{:});
end
m=m+n; nshrnk=nshrnk+n;
if ifpr==1, fprintf(' :SHR \n'); end
end
end

% Evaluate parameters to check convergence
favg=sum(f)/N; Erf=sqrt(sum((f-favg).^2)/n);
xcent=sum(x,2)/N; xdif=x-repmat(xcent,1,N);
Erx=max(sqrt(sum(xdif.^2)));

end % Bottom of the optimization loop

xmin=x(:,1); fmin=f(1);
ntype=[ninit,nrefl,nexpn,ncontr,nshrnk]; ```