function mckinnon_nedler % Derivative-free Nedler-Mead algorithm lambda_1 = (1+sqrt(33))/8; lambda_2 = (1-sqrt(33))/8; % initial simplex X = [[0;0], [lambda_1^1; lambda_2^1], [lambda_1^2; lambda_2^2]]; N = size(X,2)-1; % evaluate the function F = zeros(N+1,1); for i=1:N+1, F(i) = f(X(:,i)); end % sort the points according to the function values [X,F]=sort_points(X,F); vizualize0; HL=[]; Nmax = 3000; for k=1:Nmax, [X,F]=nedler_step(X,F); fbar = mean(F); if(fbar < 1.0E-12), break; end printf('It: %3d, f: %15.6e\n', k, fbar); HL=vizualize1(X,F,HL); drawnow; pause(0.1); end function [X,F]=sort_points(X,F) [F,I]=sort(F); X=X(:,I); function [X,F]=nedler_step(X,F) N=length(F)-1; xbar = mean(X(:,1:N),2); xref = @(t) xbar + t*(X(:,N+1)-xbar); x_m1 = xref(-1); f_m1 = f(x_m1); if (F(1) <= f_m1) && (f_m1 < F(N)) X(:,N+1) = x_m1; F(N+1) = f_m1; [X,F] = sort_points(X,F); return; elseif f_m1 < F(1) x_m2 = xref(-2); f_m2 = f(x_m2); if f_m2 < f_m1 X(:,N+1) = x_m2; F(N+1) = f_m2; [X,F] = sort_points(X,F); return; else X(:,N+1) = x_m1; F(N+1) = f_m1; [X,F] = sort_points(X,F); return; end elseif f_m1 >= F(N) if (F(N) <= f_m1) && (f_m1 < F(N+1)) x_mh = xref(-0.5); f_mh = f(x_mh); if f_mh <= f_m1 X(:,N+1) = x_mh; F(N+1) = f_mh; [X,F] = sort_points(X,F); return; end else x_ph = xref(0.5); f_ph = f(x_ph); if f_ph < F(N+1) X(:,N+1) = x_ph; F(N+1) = f_ph; [X,F] = sort_points(X,F); return; end end for i=2:N+1, X(:,i) = (X(:,1)+X(:,i))/2; F(i) = f(X(:,i)); end [X,F] = sort_points(X,F); return; end function vizualize0() clf; hold on; x1=linspace(-1.5,1.5,100); y1=linspace(-1,2,100); [x2,y2]=meshgrid(x1,y1); f2=zeros(size(x2)); for i=1:length(x2(:)), f2(i) = f([x2(i),y2(i)]); end plot(0,-1/2,'Color','r','Marker','*','MarkerSize',10); fmax = max(f2(:)); fmin = 0; p = 2; z=linspace(fmin^(1/p),fmax^(1/p),50); contour(x2,y2,f2,z.^p); hold off; HA=gca; function [HL]=vizualize1(X,F,HL) delete(HL); hold on; N = length(F)-1; HL=line(X(1,[1:N+1,1]),X(2,[1:N+1,1]),'Color','k','Marker','x'); hold off; function r=f(x) %tau = 1.0; theta = 15.0; phi = 10.0; %tau = 2.0; theta = 6.0; phi = 60.0; tau = 3.0; theta = 6.0; phi = 400.0; if x(1)<0, r = theta*phi*abs(x(1))^tau + x(2) + x(2)^2; else r = theta*x(1)^tau + x(2) + x(2)^2; end