Getting an infinite loop in Nelder-Mead code, help please

[Godric]

I've been making a go at writing out algorithms myself in MATLAB rather than using pre-existing code. I've just attempted to write the Nelder-Mead optimization method and I have hit a stumbling block where the code is now looping infinitely. It seems to have something to do with sorting the vertices again after reflecting/expanding/contracting/shrinking. As whenever I cancel the operation this is where it has gotten to. Perhaps I made a mistake in my sorting function? It's probably an obvious mistake. A huge thank you to anyone who can help me!

Here's my code:

function [solution] = NMS()
clc
%Initialize Simplex
N = 15; %number of vertices
dim = 30; %dimensions
LB = -100; %lower bounds
UB = 100; %upper bounds
errorgoal = 0.01;
goal = 0;
imax = 200; %maximum iterations
objective_function = @(x) sum(x.^2); %this is just a test function

vertices=zeros(N,dim);
vert_eval=zeros(N,1);

%intializing vertices
for i=1:N
    for j=1:dim
     vertices(i,j)=LB+(UB-LB).*rand(1);
    end
    vert_eval(i)=objective_function(vertices(i,:));
end

%Simplex manipulation constants
RE=1; %reflection coefficient
EX=2; %expansion coefficient
CO=0.5; %outside contraction coefficient
CI=-0.5; %inside contraction coefficient

%Sort the simplex
SortVertices = SortEm(vertices,vert_eval,objective_function,N);

for j=1:N
vert_eval(j)=objective_function(SortVertices(j,:));
endwhile(abs(vert_eval(1)-goal) > errorgoal)
%set fcount
fcount = N;

%Find centroid
centroid = CalcCentroid(SortVertices,dim,N);
cent_eval = objective_function(centroid);

%Calculate Reflection
SortVertices(N,:);
ref = (1+RE).*centroid - RE.*SortVertices(N,:);
ref_eval = objective_function(ref);
fcount = fcount + 1;

if vert_eval(1) <= ref_eval < vert_eval(N) %Reflect
     if fcount==imax
         break;
     end
     SortVertices(N,:)=ref;
  
elseif ref_eval < vert_eval(1) %Expand
      if fcount==imax
         break;
      end
     exp = (1+EX).*centroid - EX.*SortVertices(N,:);
     ex_eval = objective_function(exp);
     f_count = f_count+1;
  
     if ex_eval < ref_eval
         SortVertices(N,:)=exp;
      
     else
         SortVertices(N,:)=ref;
      
     end
elseif vert_eval(N-1) <= ref_eval < vert_eval(N) %Outside contrac
     if fcount==imax
         break;
     end
     con = (1+CO).*centroid - CO.*SortVertices(N,:);
     co_eval = objective_function(con);
     f_count = f_count + 1;
  
     if co_eval < ref_eval
         SortVertices(N,:)=con;
      
     else %shrink
         if fcount >= imax - N-1
             break;
          
         else
             for j=2:N
                 SortVertices(j,:) = SortVertices(1)-(SortVertices(j) - SortVertices(1))./2;
                 vert_eval(j)=objective_function(SortVertices(j,:));
             end
         end
     end
  
elseif ref_eval > vert_eval(N) %Inside contrac
     if fcount==imax
         break;
     end
     coi = (1+CI).*centroid - CI.*SortVertices(N,:);
     co_eval = objective_function(coi);
     f_count = f_count + 1;
   
     if co_eval < ref_eval
         SortVertices(N,:)=coi;
      
     else %shrink
         if fcount >= imax - N-1
             break
          
         else
             for j=2:N
                 SortVertices(j,:) = SortVertices(1)-(SortVertices(j) - SortVertices(1))./2;
                 vert_eval(j)=objective_function(SortVertices(j,:));
             end
         end
     end
  
end
  vertices = SortVertices;
  SortVertices = SortEm(vertices,vert_eval,objective_function,N);
  for j=1:N
   vert_eval(j)=objective_function(SortVertices(j,:));
  end
end
solution = vert_eval
end

%Vertex sorting
function [SortVertices] = SortEm(vertices,vert_eval,objective_function,N)
for m=1:N
        for i=1:N-1
            if(vert_eval(i)>vert_eval(i+1))
                temp=vertices(i,:);
                vertices(i,:)= vertices(i+1,:);
                vertices(i+1,:)= temp;
            else
                continue;
            end
        end  
      for k=1:N
       vert_eval(k)=objective_function(vertices(k,:));
      end
end
   SortVertices=vertices;
end

%Calculate the centroid
function [centroid] = CalcCentroid(SortVertices,dim,N)
SumVert=zeros(1,dim);
i=1:dim;
for k=1:N-1
  SumVert=SumVert+SortVertices(k,i);
end
centroid=SumVert./(N-1);
end

 

[berkeman]

Maybe add some debug output statements to show you what it is doing?

 

[Godric]

So it appears the worst value just switches between two values every loop and never gets better. I am not certain what could be causing this effect.

The infinite looping was a dumb mistake where the f_count was getting reset, but the actual algorithm is still not working because of the constant back and forth.

Edit: I think I solved my issues.