function [ xmin, ynewlo, icount, numres, ifault ] = nelmin ( fn, n, start, ... reqmin, step, konvge, kcount ) %*****************************************************************************80 % %% NELMIN minimizes a function using the Nelder-Mead algorithm. % % Discussion: % % This routine seeks the minimum value of a user-specified function. % % Simplex function minimisation procedure due to Nelder+Mead(1965), % as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with % subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976, % 25, 97) and Hill(1978, 27, 380-2) % % The function to be minimized must have the form: % % function value = fn ( x ) % % where X is a vector, and VALUE is the (scalar) function value. % The name of this function must be passed as the argument FN. % % This routine does not include a termination test using the % fitting of a quadratic surface. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 27 October 2011 % % Author: % % Original FORTRAN77 version by R ONeill. % MATLAB version by John Burkardt. % % Reference: % % John Nelder, Roger Mead, % A simplex method for function minimization, % Computer Journal, % Volume 7, 1965, pages 308-313. % % R ONeill, % Algorithm AS 47: % Function Minimization Using a Simplex Procedure, % Applied Statistics, % Volume 20, Number 3, 1971, pages 338-345. % % Parameters: % % Input, real value = FN ( X ), the name of the MATLAB function which % evaluates the function to be minimized, preceded by an "@" sign. % % Input, integer N, the number of variables. % % Input, real START(N). On input, a starting point % for the iteration. On output, this data may have been overwritten. % % Input, real REQMIN, the terminating limit for the variance % of function values. % % Input, real STEP(N), determines the size and shape of the % initial simplex. The relative magnitudes of its elements should reflect % the units of the variables. % % Input, integer KONVGE, the convergence check is carried out % every KONVGE iterations. % % Input, integer KCOUNT, the maximum number of function % evaluations. % % Output, real XMIN(N), the coordinates of the point which % is estimated to minimize the function. % % Output, real YNEWLO, the minimum value of the function. % % Output, integer ICOUNT, the number of function evaluations. % % Output, integer NUMRES, the number of restarts. % % Output, integer IFAULT, error indicator. % 0, no errors detected. % 1, REQMIN, N, or KONVGE has an illegal value. % 2, iteration terminated because KCOUNT was exceeded without convergence. % xmin = []; ynewlo = []; icount = 0; numres = 0; ccoeff = 0.5; ecoeff = 2.0; eps = 0.001; rcoeff = 1.0; % % Check the input parameters. % if ( reqmin <= 0.0 ) ifault = 1; return end if ( n < 1 ) ifault = 1; return end if ( konvge < 1 ) ifault = 1; return end jcount = konvge; dn = n; nn = n + 1; dnn = nn; del = 1.0; rq = reqmin * dn; % % Initial or restarted loop. % while ( 1 ) for i = 1 : n p(i,nn) = start(i); end y(nn) = fn ( start ); icount = icount + 1; for j = 1 : n x = start(j); start(j) = start(j) + step(j) * del; for i = 1 : n p(i,j) = start(i); end y(j) = fn ( start ); icount = icount + 1; start(j) = x; end % % The simplex construction is complete. % % Find highest and lowest Y values. YNEWLO = Y(IHI) indicates % the vertex of the simplex to be replaced. % ylo = y(1); ilo = 1; for i = 2 : nn if ( y(i) < ylo ) ylo = y(i); ilo = i; end end % % Inner loop. % while ( 1 ) if ( kcount <= icount ) break end ynewlo = y(1); ihi = 1; for i = 2 : nn if ( ynewlo < y(i) ) ynewlo = y(i); ihi = i; end end % % Calculate PBAR, the centroid of the simplex vertices % excepting the vertex with Y value YNEWLO. % for i = 1 : n z = 0.0; for j = 1 : nn z = z + p(i,j); end z = z - p(i,ihi); pbar(i) = z / dn; end % % Reflection through the centroid. % for i = 1 : n pstar(i) = pbar(i) + rcoeff * ( pbar(i) - p(i,ihi) ); end ystar = fn ( pstar ); icount = icount + 1; % % Successful reflection, so extension. % if ( ystar < ylo ) for i = 1 : n p2star(i) = pbar(i) + ecoeff * ( pstar(i) - pbar(i) ); end y2star = fn ( p2star ); icount = icount + 1; % % Check extension. % if ( ystar < y2star ) for i = 1 : n p(i,ihi) = pstar(i); end y(ihi) = ystar; % % Retain extension or contraction. % else for i = 1 : n p(i,ihi) = p2star(i); end y(ihi) = y2star; end % % No extension. % else l = 0; for i = 1 : nn if ( ystar < y(i) ) l = l + 1; end end if ( 1 < l ) for i = 1 : n p(i,ihi) = pstar(i); end y(ihi) = ystar; % % Contraction on the Y(IHI) side of the centroid. % elseif ( l == 0 ) for i = 1 : n p2star(i) = pbar(i) + ccoeff * ( p(i,ihi) - pbar(i) ); end y2star = fn ( p2star ); icount = icount + 1; % % Contract the whole simplex. % if ( y(ihi) < y2star ) for j = 1 : nn for i = 1 : n p(i,j) = ( p(i,j) + p(i,ilo) ) * 0.5; xmin(i) = p(i,j); end y(j) = fn ( xmin ); icount = icount + 1; end ylo = y(1); ilo = 1; for i = 2 : nn if ( y(i) < ylo ) ylo = y(i); ilo = i; end end continue % % Retain contraction. % else for i = 1 : n p(i,ihi) = p2star(i); end y(ihi) = y2star; end % % Contraction on the reflection side of the centroid. % elseif ( l == 1 ) for i = 1 : n p2star(i) = pbar(i) + ccoeff * ( pstar(i) - pbar(i) ); end y2star = fn ( p2star ); icount = icount + 1; % % Retain reflection? % if ( y2star <= ystar ) for i = 1 : n p(i,ihi) = p2star(i); end y(ihi) = y2star; else for i = 1 : n p(i,ihi) = pstar(i); end y(ihi) = ystar; end end end % % Check if YLO improved. % if ( y(ihi) < ylo ) ylo = y(ihi); ilo = ihi; end jcount = jcount - 1; if ( 0 < jcount ) continue end % % Check to see if minimum reached. % if ( icount <= kcount ) jcount = konvge; z = 0.0; for i = 1 : nn z = z + y(i); end x = z / dnn; z = 0.0; for i = 1 : nn z = z + ( y(i) - x )^2; end if ( z <= rq ) break end end end % % Factorial tests to check that YNEWLO is a local minimum. % for i = 1 : n xmin(i) = p(i,ilo); end ynewlo = y(ilo); if ( kcount < icount ) ifault = 2; break end ifault = 0; for i = 1 : n del = step(i) * eps; xmin(i) = xmin(i) + del; z = fn ( xmin ); icount = icount + 1; if ( z < ynewlo ) ifault = 2; break end xmin(i) = xmin(i) - del - del; z = fn ( xmin ); icount = icount + 1; if ( z < ynewlo ) ifault = 2; break end xmin(i) = xmin(i) + del; end if ( ifault == 0 ) break end % % Restart the procedure. % for i = 1 : n start(i) = xmin(i); end del = eps; numres = numres + 1; end return end