SS.rar_ scatter search_scatter search_ss

function [x,y,X,Y] = simplex(func,x,opts) % % SIMPLEX - multidimensional unconstrained non-linear optimsiation % % X = SIMPLEX(FUNC,X) finds a local minumum of a function, via a function % handle FUNC, starting from an initial point X. The local minimum is % located via the Nelder-Mead simplex algorithm [1], which does not require % any gradient information. % % [X,Y] = SIMPLEX(FUNC,X) also returns the value of the function, Y, at % the local minimum, X. % % X = SIMPLEX(FUNC,X,OPTS) allows the optimisation parameters to be % specified via a structure, OPTS, with members % % opts.Chi - Parameter governing expansion steps % opts.Delta - Parameter governing size of initial simplex. % opts.Gamma - Parameter governing contraction steps. % opts.Rho - Parameter governing reflection steps. % opts.Sigma - Parameter governing shrinkage steps. % opts.MaxIter - Maximum number of optimisation steps. % opts.MaxFunEvals - Maximum number of function evaluations. % opts.TolFun - Stopping criterion based on the relative change in % value of the function in each step. % opts.TolX - Stopping criterion based on the change in the % minimiser in each step. % % OPTS = SIMPLEX() returns a structure containing the default optimisation % parameters, with the following values: % % opts.Chi = 2 % opts.Delta = 0.01 % opts.Gamma = 0.5 % opts.Rho = 1 % opts.Sigma = 0.5 % opts.MaxIter = 200 % opts.MaxFunEvals = 1000 % opts.TolFun = 1e-6 % opts.TolX = 1e-6 % % X = SIMPLEX(FUNC,X,OPTS, P1, P2, ...) allows addinal parameters to be % passed to the function to be minimised. % % [X,Y,XX,YY] = SIMPLEX(FUNC, X) also returns in XX all of the values of % X evaluated during the optimisation process and in YY the corresponding % values of the function. % % References: % % [1] J. A. Nelder and R. Mead, "A simplex method for function % minimization", Computer Journal, 7:308-313, 1965. % % File : simplex.m % % Date : Friday 5th January 2007 % % Author : Dr Gavin C. Cawley % % Description : Simple implementation of the Nelder-Mead simplex optimisation % algorithm [1]. Similar to the fminsearch routine from the % MATLAB optimisation toolbox. % % References : [1] J. A. Nelder and R. Mead, "A simplex method for function % minimization", Computer Journal, 7:308-313, 1965. % % History : 05/01/2007 - v1.00 % % Copyright : (c) Dr Gavin C. Cawley, January 2007. % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA % % use default optimisation parameters if none given % if nargin < 3 % % func='fcolville' % opts.Chi = 2; % opts.Delta = 0.01; % opts.Gamma = 0.5; % opts.Rho = 1; % opts.Sigma = 0.5; % opts.MaxIter = 500; % opts.MaxFunEvals = 1000; % opts.TolFun = 1e-6; % opts.TolX = 1e-6; % % end % return structure containing default optimisation parameters if nargin == 0 x = opts; return end % get initial parameters x = x(:); n = length(x); x = repmat(x', n+1, 1); y = zeros(n+1, 1); % form initial simplex for i=1:n x(i,i) = x(i,i) + opts.Delta; y(i) = feval(func,x(i,:)); end y(n+1) = feval(func,x(n+1,:)); X = x; Y = y; count = n+1; format = ' % 4d % 4d % 12f %s\n'; fprintf(1, '\n Iteration Func-count min f(x) Procedure\n\n'); fprintf(1, format, 1, count, min(y), 'initial'); % iterative improvement for i=2:opts.MaxIter % order [y,idx] = sort(y); x = x(idx,:); % reflect centroid = mean(x(1:end-1,:)); x_r = centroid + opts.Rho*(centroid - x(end,:)); y_r = feval(func,x_r); count = count + 1; X = [X ; x_r]; Y = [Y ; y_r]; if y_r >= y(1) & y_r < y(end-1) % accept reflection point x(end,:) = x_r; y(end) = y_r; fprintf(1, format, i, count, min(y), 'reflect'); else if y_r < y(1) % expand x_e = centroid + opts.Chi*(x_r - centroid); y_e = feval(func,x_e); count = count + 1; X = [X ; x_e]; Y = [Y ; y_e]; if y_e < y_r % accept expansion point x(end,:) = x_e; y(end) = y_e; fprintf(1, format, i, count, min(y), 'expand'); else % accept reflection point x(end,:) = x_r; y(end) = y_r; fprintf(1, format, i, count, min(y), 'reflect'); end else % contract shrink = 0; if y(end-1) <= y_r & y_r < y(end) % contract outside x_c = centroid + opts.Gamma*(x_r - centroid); y_c = feval(func,x_c); count = count + 1; X = [X ; x_c]; Y = [Y ; y_c]; if y_c <= y_r % accept contraction point x(end,:) = x_c; y(end) = y_c; fprintf(1, format, i, count, min(y), 'contract outside'); else shrink = 1; end else % contract inside x_c = centroid + opts.Gamma*(centroid - x(end,:)); y_c = feval(func,x_c); count = count + 1; X = [X ; x_c]; Y = [Y ; y_c]; if y_c <= y(end) % accept contraction point x(end,:) = x_c; y(end) = y_c; fprintf(1, format, i, count, min(y), 'contract inside'); else shrink = 1; end end if shrink % shrink for j=2:n+1 x(j,:) = x(1,:) + opts.Sigma*(x(j,:) - x(1,:)); y(j) = feval(func,x(j,:)); count = count + 1; X = [X ; x(j,:)]; Y = [Y ; y(j)]; end fprintf(1, format, i, count, min(y), 'shrink'); end end end % evaluate stopping criterion if max(abs(min(x) - max(x))) < opts.TolX fprintf(1, 'optimisation terminated sucessfully (TolX criterion)\n'); break; end if abs(max(y) - min(y))/max(abs(y)) < opts.TolFun fprintf(1, 'optimisation terminated sucessfully (TolFun criterion)\n'); break; end end if i == opts.MaxIter fprintf(1, 'Warning : maximim number of iterations exceeded\n'); end % update model structure [y, idx] = min(y); x = x(idx,:); % bye bye...