Neldermead

function [ x_opt, has ] = neldermead( x, function_handle)
% Define algorithm constants

%x – vector n+1 by 2
rho = 1; % rho > 0 %koef reflection
xi = 2; % xi > max(rho, 1) %koef expansion
gam = 0.5; % 0 < gam < 1
sig = 0.5; % 0 < sig < 1

tolerance = 1e-8;
max_feval = 250;

% Initialization
[temp,n_dim] = size(x);

index = 1:n_dim+1;

[f] = evaluate(x,function_handle);
%[f] = feval(function_handle,x);
n_feval = n_dim+1;

[f,index] = sort(f); % sorts values from smallest to largest…
x = x(index,:);

% Begin Nelder Mead iterations
deltaf=1000;
it=0;
while deltaf >10e-9
% compute midpoint of simplex opposite worst point
x_bar = sum(x(1:n_dim,:)) / n_dim; %mean
% compute reflection point
x_r = (1+rho)*x_bar – rho*x(n_dim+1,:);

f_r = feval(function_handle,x_r);
n_feval = n_feval+1;

if (f(1) <= f_r & f_r <= f(n_dim)) % accept the point
x(n_dim+1,:) = x_r;
f(n_dim+1) = f_r;

elseif (f_r < f(1)) % test for possible expansion
x_e = (1+rho*xi)*x_bar – rho*xi*x(n_dim+1,:);
f_e = feval(function_handle,x_e);
n_feval = n_feval+1;

if (f_e<f_r) % accept the expanded point
x(n_dim+1,:) = x_e;
f(n_dim+1 ) = f_e;
else
x(n_dim+1,:) = x_r;
f(n_dim+1 ) = f_r;
end
elseif (f(n_dim) <= f_r & f_r < f(n_dim+1)) % outside contraction
x_c = (1+rho*gam)*x_bar – rho*gam*x(n_dim+1,:);
f_c = feval(function_handle,x_c); n_feval = n_feval+1;

if (f_c <= f_r) % accept the contracted point
x(n_dim+1,:) = x_c;
f(n_dim+1 ) = f_c;

else
[x,f] = shrink(x,function_handle,sig); n_feval = n_feval+n_dim;

end

else % f_r must be >= f(n_dim+1), try an inside contraction
x_c = (1-gam)*x_bar + gam*x(n_dim+1,:);
f_c = feval(function_handle,x_c); n_feval = n_feval+1;

if (f_c < f(n_dim+1)) % accept the contracted point
x(n_dim+1,:) = x_c;
f(n_dim+1 ) = f_c;
else
[x,f] = shrink(x,function_handle,sig);
n_feval = n_feval+n_dim;
end

end

% resort the points. note that we are not implementing the usual
% Nelder-Mead tie-breaking rules (when f(1) = f(2) or f(n_dim) =
% f(n_dim+1)…
[f,index] = sort(f); % sorts values from smallest to largest…
x = x(index,:)

% test for convergence
deltaf = abs(f(1)-f(end));
end

it=it+1;
x_opt = x(1,:);
end % function nelder_mead

function [f] = evaluate(x,function_handle);

[temp,n_dim] = size(x);

f = zeros(1,n_dim+1);

for i=1:n_dim+1
f(i) = feval(function_handle,x(i,:));
end
end

function [x,f] = shrink(x,function_handle,sig)
% In the worst case, we need to shrink the simplex along each edge towards
% the current “best” point. This is quite expensive, requiring n_dim new
% function evaluations.

[temp,n_dim] = size(x);

x1 = x(1,:);
f(1) = feval(function_handle,x1);

for i=2:n_dim+1
x(i,:) = sig*x(i,:) + (1-sig)*x(1,:);
f(i) = feval(function_handle,x(i,:));
end
end

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s