Algoritmos de otimização populacional: Método Nelder-Mead (NM)

Conteúdo:

1. Introdução
2. Descrição do algoritmo
3. Resultados dos testes

1. Introdução

O método Nelder-Mead foi desenvolvido em 1965 por John Nelder e Roger Mead. Eles procuravam um método de otimização que pudesse funcionar com funções que não possuíam derivadas ou que não tinham fórmulas analíticas para suas derivadas. Eles também queriam desenvolver um método que fosse simples de implementar e eficiente para uso em computadores daquela época. As pesquisas os levaram à ideia de usar um simplex, um poliedro no espaço de parâmetros da função.

A história do desenvolvimento do método começou com o trabalho de John Nelder e seus colegas no laboratório de ciências da computação em Oxford. Eles enfrentaram o problema de otimização de funções que não possuíam derivadas analíticas ou eram muito complexas para calcular. Métodos de otimização tradicionais, como os métodos de gradiente, não eram aplicáveis nesses casos. Em vez disso, Nelder e Mead propuseram um novo método, que se baseava na busca iterativa da solução ótima no espaço de parâmetros da função.

O método Nelder-Mead foi chamado de "método simplex" e foi publicado no artigo "A Simplex Method for Function Minimization" no jornal "The Computer Journal" em 1965. Este método foi aceito pela comunidade científica e passou a ser amplamente utilizado em diversas áreas que requerem otimização de funções.

O simplex é um conjunto de pontos que formam um poliedro, onde cada ponto é um conjunto de valores de parâmetros da função a ser otimizada. A ideia é alterar e mover o simplex no espaço de parâmetros para encontrar o valor ótimo da função.

O método Nelder-Mead (método simplex de Nelder-Mead) pertence à classe de algoritmos de otimização irrestrita. É um algoritmo determinístico que não requer o uso de derivadas da função e pode operar com funções que possuem vários mínimos locais.

2. Descrição do algoritmo

O método Nelder-Mead não é considerado um algoritmo populacional porque utiliza apenas um agente de otimização, representando o simplex, que por sua vez consiste em um conjunto de pontos no espaço de busca, o que por si só pode ser considerado uma população. No entanto, usaremos vários agentes, cada um com seu próprio simplex, portanto, essa implementação pode ser considerada como um algoritmo populacional.

Assim, para simplificar o entendimento das operações com o simplex, suponhamos que é necessário otimizar uma função de duas variáveis, ou seja, a dimensão do espaço z = 2, então o simplex consistirá de z + 1 = 3 vértices. A alteração do simplex é uma operação sobre o pior ponto desses três. Vamos denominar esses pontos como Best, Good, Worst, sendo Good o segundo ponto após Worst (no caso de a dimensão ser maior que dois, Good será sempre o segundo após Worst da lista ordenada de vértices).

A seguir, com esses três vértices, é necessário realizar cálculos e ajustar o Worst em relação aos demais vértices. Movemos o Worst com base no centro geométrico dos vértices, excluindo-o, ou seja, calculamos o centroide e movemos o Worst em relação a esse centroide.

Em cada iteração do método de Nelder-Mead, executamos os seguintes passos:

1. Ordenação dos vértices do simplex pelos valores da função objetivo, ou seja, 

f(P0) ⩾ f(P1) ⩾ ... ⩾ f(P(n-1))

onde:

• P0, P1 ... P(n-1) são pontos do simplex
• n é o número de vértices do simplex
  

2. Cálculo do centroide: calcular o valor médio correspondente às coordenadas de todos os vértices do simplex, exceto o pior vértice. Suponha que o vetor Xo seja o centro de massa de todos os vértices, exceto X(n-1), então: 

Xo = (X0 + X1 + ... + X(n-2)) / (n-1)

onde:

• X0, X1 ... X(n-2) são os vetores de coordenadas correspondentes de todos os vértices, exceto o pior

3. Reflexão (Reflection): operação de reflexão do pior vértice em relação ao centroide. O vetor Xr do novo vértice refletido é calculado pela seguinte fórmula:

Xr = Xo + α (Xo - Xw)

onde:

• Xw é o vetor do pior vértice, correspondendo ao vértice P(n-1)
  
• α é o coeficiente de reflexão, normalmente igual a 1
  

4. Avaliação do novo vértice: calcula-se o valor da função objetivo neste ponto. Se o novo vértice tiver uma função objetivo melhor do que o pior vértice do simplex, ela pode substituí-lo, ou seja, a reflexão substitui o original. A operação de reflexão permite explorar uma nova área do espaço de parâmetros, refletindo-se do centroide do simplex.

5. Expansão (Expansion): usada para explorar ainda mais o espaço de parâmetros quando a operação de reflexão produziu um vértice melhor do que o melhor, sob a suposição de que um movimento adicional na direção da reflexão permitirá melhorar ainda mais sua posição. A operação de expansão permite explorar uma área maior do espaço de parâmetros do que a operação de reflexão. Ela aumenta a distância entre o centroide e o ponto refletido, o que pode ajudar a encontrar novos mínimos locais da função objetivo. O vetor de expansão Xe pode ser calculado da seguinte maneira:

Xe = Xo + γ(Xr - Xo)

onde:

• γ é o coeficiente de expansão, geralmente maior que 1, por padrão 2

5.1. No caso de uma operação de expansão: calcula-se a função de adaptação para o vértice Xe e, se for melhor do que o vértice Xr, então ela pode substituir a pior vértice Xw do simplex. Após a execução da expansão, voltamos ao passo 1.

É importante observar que o coeficiente de expansão γ deve ser escolhido com cuidado. Se for escolhido muito alto, pode levar à expansão do simplex em uma grande área do espaço de parâmetros, o que pode retardar a convergência ou resultar na perda de informações sobre extremos locais. Se for escolhido muito baixo, pode não explorar adequadamente novas áreas do espaço de parâmetros.

6. Contração (Contraction): usada para calcular o vértice Xc, quando a operação de reflexão produziu um vértice pior ou igual ao bom ponto Xb (segundo após o pior), ou seja, é realizada na esperança de encontrar uma posição melhor dentro do simplex. Ela é realizada da seguinte forma:

Xc = Xo + β(Xw - Xo)

onde:

• β - coeficiente de contração, geralmente entre 0 e 1, padrão 0.5

6.1.  No caso de uma operação de contração: calcula-se o valor da função objetivo para o vértice Xc. Se o novo vértice tiver uma função objetivo melhor do que o pior vértice do simplex, então ela pode substituir o pior. Após a execução da expansão, voltamos ao passo 1.

A operação de contração permite deslocar o simplex mais próximo ao pior ponto, o que pode ajudar na exploração da vizinhança desse ponto. Como na operação de expansão, é importante escolher o coeficiente de contração β com cautela. Se for escolhido muito alto, isso pode levar a uma contração excessiva do simplex, resultando na perda de informações sobre mínimos locais. Se for escolhido muito baixo, a contração pode ser insuficiente para explorar a vizinhança do pior ponto.

7. Encolhimento (Shrinkage): a última operação no método de Nelder-Mead é realizada quando nenhuma das operações anteriores (reflexão, expansão, contração) leva a uma melhoria no valor da função objetivo. A operação de encolhimento permite reduzir o tamanho do simplex para estreitar a área de busca e focar ao redor do melhor ponto.

Como vemos, os autores têm quatro operações com o simplex. Ao trabalhar com simplexes, é necessário ter em mente - para iniciar a otimização, é necessário calcular a função de adaptação para todos os pontos do simplex. O número de pontos do simplex é igual a "z+1", onde z é a dimensão do espaço de busca. Por exemplo, se temos uma dimensão de busca de 1000, então devemos calcular a função de adaptação 1001 vezes para começar as operações com o simplex. Enquanto que com um algoritmo populacional comum com uma população de 50 agentes, poderíamos realizar 20 épocas, para este algoritmo será realizada apenas uma época e, ao mesmo tempo, será gasta a maior parte do número limitado de iterações.

A operação "Shrinkage" é a mudança de todos os pontos para o melhor, após o qual será necessário calcular a função de adaptação 1000 vezes, ou seja, isso significa que a operação "Shrinkage" é muito custosa. Segundo a ideia dos autores, essa operação deve ser chamada quando o algoritmo está preso, quando o movimento do pior ponto do simplex não traz melhorias na solução. No entanto, meus experimentos mostraram que essa operação é muito custosa em termos de recursos computacionais e, além disso, é inútil e autodestrutiva para o algoritmo de otimização, pois a convergência de todos os pontos do simplex-agente a um único extremo local significaria ficar preso e parar a exploração do espaço de busca. Por isso, decidi desistir da operação "Shrinkage" para a implementação prática do algoritmo.

Assim, usaremos as primeiras 3 operações com o simplex, que podem ser visualizadas no diagrama 1:

Figura 1. As três operações principais do método de Nelder-Mead.

Como o método de busca simplex pode enfrentar uma escolha infeliz de vértices iniciais, sendo determinístico, o resultado de sua otimização pode ser insatisfatório. Experimentos com este algoritmo apenas confirmaram as preocupações: ele funciona aceitavelmente apenas em um espaço limitado de coordenadas.

Por causa de problemas com a convergência, os vértices do simplex simplesmente se deslocam com um valor fixo. Imagine que você está em pernas de pau tentando sentar-se em um banco, mas ele está muito baixo para isso. A mesma situação ocorre com o simplex em mudança. Para que o algoritmo funcione perfeitamente, é necessário ter muita sorte para que os vértices iniciais do simplex estejam no lugar certo. Isso é semelhante à situação com as pernas de pau: para sentar-se com segurança, é necessária uma banqueta alta. O algoritmo converge relativamente bem em funções monodimensionais suaves, como uma parábola. No entanto, nossas tarefas práticas são muito mais complexas e geralmente cheias de "armadilhas" locais, especialmente em tarefas de negociação, que são por natureza discretas.

Então surge a questão: o que fazer quando o simplex fica "preso" para sempre em um extremo local? A lógica do algoritmo de Nelder-Mead não prevê maneiras de escapar dessa "armadilha". Após cada operação, seja contração ou expansão, o algoritmo retorna à reflexão, tentando superar o pior ponto. Visualmente, isso parece um "piscar". Para resolver esse problema, tentei dar ao simplex a oportunidade de deixar a "armadilha" local e continuar a busca em novas áreas do espaço. Para isso, usei voos de Levy, que em alguns casos ajudaram a "reviver" a população, como descrito em artigos anteriores. No entanto, vale ressaltar que os voos de Levy nem sempre são úteis e sua aplicação depende da lógica do algoritmo de otimização. Em nosso caso, foi um experimento, e os resultados de melhoria não foram garantidos.

Vamos para a parte mais interessante: escrever o código.

Após cada iteração, é necessário saber qual operação foi realizada com o pior vértice do simplex. Para isso, é conveniente usar um enumerador. Além das três operações principais com o vértice do simplex, adicionei mais uma, "nenhuma", para o caso de surgir uma ideia de como complementar o algoritmo. O enumerador chamaremos de E_SimplexOperation, cujos campos são:

• none: não prevê operações
• reflection: operação "reflexão"
• expansion: operação "expansão"
• contraction: operação "contração"
  

//——————————————————————————————————————————————————————————————————————————————
enum E_SimplexOperation
{
  none        = 0,
  reflection  = 1,
  expansion   = 2,
  contraction = 3
};
//——————————————————————————————————————————————————————————————————————————————

Para descrever um vértice do simplex, precisaremos da estrutura S_Point, que contém os seguintes campos:

• Init(int coords): inicialização do ponto, recebe o argumento "coords"
• c []: array de coordenadas do ponto em um espaço multidimensional, o tamanho do array é determinado no método "Init"
• f: valor da função de adaptação do vértice
  

//——————————————————————————————————————————————————————————————————————————————
struct S_Point
{
  void Init (int coords)
  {
    ArrayResize (c, coords);
    f = -DBL_MAX;
  }
  double c [];
  double f;
};
//——————————————————————————————————————————————————————————————————————————————

Para o simplex, usaremos a estrutura S_Simplex, que representa um simplex simples para cada agente e contém os seguintes campos:

• Init(int coords): inicialização do simplex, recebe o argumento "coords" - quantidade de coordenadas
• S_Point p []: array de vértices do simplex, o tamanho do array é determinado no método "Init"
• c []: array de coordenadas do centroide do simplex, o tamanho do array é determinado no método "Init".
• S_Point Xr: ponto de reflexão
• S_Point Xe: ponto de expansão
• S_Point Xc: ponto de contração
• indG: índice de um bom ponto (Good) no simplex
• indW: índice do pior ponto no simplex
• E_SimplexOperation: tipo de operação do simplex

//——————————————————————————————————————————————————————————————————————————————
struct S_Simplex
{
  void Init (int coords)
  {
    ArrayResize (p, coords + 1);

    for (int i = 0; i <= coords; i++)
    {
      p [i].Init (coords);
    }

    ArrayResize (c, coords);

    Xr.Init (coords);
    Xe.Init (coords);
    Xc.Init (coords);

    operation = none;
  }

  S_Point p [];
  double c  []; //centroid

  S_Point Xr;  //reflection point
  S_Point Xe;  //expansion point
  S_Point Xc;  //expansion point

  int indG;    //index of good point
  int indW;    //index of worst point

  E_SimplexOperation operation; //type of simplex operation
};
//——————————————————————————————————————————————————————————————————————————————

Para o agente de otimização, definiremos a estrutura S_Agent, que contém os seguintes campos:

• Init(int coords): método de inicialização do agente, recebe o argumento "coords", que indica a quantidade de coordenadas. O método altera o tamanho do array "c" para "coords" usando a função ArrayResize. Então, o campo "f" é definido para o valor -DBL_MAX, que é o valor inicial para a função de avaliação do agente. Em seguida, o método "Init" é chamado para o campo "s", que representa o simplex do agente.
• c []: array de coordenadas do agente em um espaço multidimensional para interação com o programa externo. O tamanho do array é definido no método "Init".
• f: valor da função de adaptação do agente.
• s: simplex do agente, representado pela estrutura S_Simplex. O simplex é inicializado chamando o método "Init" para o campo "s".

//——————————————————————————————————————————————————————————————————————————————
struct S_Agent
{
  void Init (int coords)
  {
    ArrayResize (c, coords);
    f = -DBL_MAX;

    s.Init (coords);
  }

  double c  []; //coordinates
  double f;     //fitness
  S_Simplex s;  //agent simplex
};
//——————————————————————————————————————————————————————————————————————————————

Definiremos a classe do algoritmo do método de busca simplex C_AO_NMm, que contém os seguintes campos:

• cB []: array das melhores coordenadas.
• fB: valor da função de adaptação das melhores coordenadas.
• a []: array de agentes, representados pela estrutura S_Agent.
• rangeMax []: array dos valores máximos do intervalo de busca para cada coordenada.
• rangeMin []: array dos valores mínimos do intervalo de busca para cada coordenada.
• rangeStep []: array dos passos de busca para cada coordenada.
• Init(const int coordsP, const int popSizeP, const double reflectionCoeffP, const double expansionCoeffP, const double contractionCoeffP): método de inicialização do algoritmo, recebe argumentos que definem a quantidade de coordenadas, o tamanho da população, e os coeficientes para as operações de reflexão, expansão e contração, e realiza a inicialização de todos os campos da classe.
• Moving(): método que realiza o movimento dos agentes no algoritmo.
• Revision(): método que realiza a revisão dos agentes no algoritmo.
• Sorting, CalcCentroid, Reflection, Expansion, Contraction, Flying: realizam operações com os simplexes.
• SeInDiSp, RNDfromCI, Scale: são executados para gerar e transformar valores em intervalos aceitáveis com o passo definido.
  

//——————————————————————————————————————————————————————————————————————————————
class C_AO_NMm
{
  //----------------------------------------------------------------------------
  public: double cB [];  //best coordinates
  public: double fB;     //FF of the best coordinates
  public: S_Agent a [];  //agent

  public: double rangeMax  []; //maximum search range
  public: double rangeMin  []; //manimum search range
  public: double rangeStep []; //step search

  public: void Init (const int    coordsP,            //coordinates number
                     const int    popSizeP,           //population size
                     const double reflectionCoeffP,   //Reflection coefficient
                     const double expansionCoeffP,    //Expansion coefficient
                     const double contractionCoeffP); //Contraction coefficient

  public: void Moving   ();
  public: void Revision ();

  //----------------------------------------------------------------------------
  private: int    coords;        //coordinates number
  private: int    popSize;       //population size
  private: int    simplexPoints; //simplex points


  private: double reflectionCoeff;  //Reflection coefficient
  private: double expansionCoeff;   //Expansion coefficient
  private: double contractionCoeff; //Contraction coefficient

  private: bool   revision;

  private: S_Point pTemp [];
  private: int     ind   [];
  private: double  val   [];

  private: void   Sorting      (S_Point &p []);
  private: void   CalcCentroid (S_Simplex &s, int indW);
  private: void   Reflection   (S_Agent &agent, int indW);
  private: void   Expansion    (S_Agent &agent);
  private: void   Contraction  (S_Agent &agent, int indW);
  private: void   Flying       (S_Agent &agent, int indW);

  private: double SeInDiSp     (double In, double InMin, double InMax, double Step);
  private: double RNDfromCI    (double min, double max);
  private: double Scale        (double In, double InMIN, double InMAX, double OutMIN, double OutMAX,  bool revers);
};
//——————————————————————————————————————————————————————————————————————————————

Para a inicialização do algoritmo, usaremos o método "Init" da classe C_AO_NMm, que realiza a inicialização de todos os campos da classe e cria os arrays necessários para o funcionamento do algoritmo de otimização.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Init (const int    coordsP,            //coordinates number
                     const int    popSizeP,           //population size
                     const double reflectionCoeffP,   //reflection coefficient
                     const double expansionCoeffP,    //expansion coefficient
                     const double contractionCoeffP)  //contraction coefficient
{
  MathSrand ((int)GetMicrosecondCount ()); // reset of the generator
  fB       = -DBL_MAX;
  revision = false;

  coords        = coordsP;
  popSize       = popSizeP;
  simplexPoints = coords + 1;

  reflectionCoeff  = reflectionCoeffP;
  expansionCoeff   = expansionCoeffP;
  contractionCoeff = contractionCoeffP;

  ArrayResize (pTemp, simplexPoints);
  ArrayResize (ind,   simplexPoints);
  ArrayResize (val,   simplexPoints);

  ArrayResize (a, popSize);
  for (int i = 0; i < popSize; i++)
  {
    a [i].Init (coords);
  }

  ArrayResize (rangeMax,  coords);
  ArrayResize (rangeMin,  coords);
  ArrayResize (rangeStep, coords);
  ArrayResize (cB,        coords);
}
//——————————————————————————————————————————————————————————————————————————————

Normalmente, usamos o método "Moving" para mover os agentes no espaço de busca, mas para o algoritmo NM, manteremos apenas a função de posicionamento aleatório inicial das vértices dos simplexes dos agentes no espaço de busca.
Para criar um vértice aleatório dos simplexes, usaremos RNDfromCI e a ajustaremos ao intervalo necessário com o passo definido pelo método SeInDiSp.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Moving ()
{
  if (!revision)
  {
    int cnt = 0;

    for (int i = 0; i < popSize; i++)
    {
      for (int p = 0; p < simplexPoints; p++)
      {
        cnt++;
        for (int c = 0; c < coords; c++)
        {
          a [i].s.p [p].c [c] = RNDfromCI (rangeMin [c], rangeMax [c]);
          a [i].s.p [p].c [c] = SeInDiSp  (a [i].s.p [p].c [c], rangeMin [c], rangeMax [c], rangeStep [c]);
        }
      }
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

A lógica principal do movimento das vértices dos simplexes dos agentes será aplicada no método "Revision".

Se a revisão ainda não foi realizada "if (!revision)", é necessário ordenar as vértices no simplex, construir o centroide e realizar a operação de reflexão para as piores vértices dos simplexes dos agentes. Esta é a primeira operação com as vértices e sempre é "reflexão". Aqui também podemos atualizar a solução global, pois neste ponto temos conhecimento sobre o valor de adaptação para cada vértice. A definição do valor da variável "revision" para "true" e a saída do método.

A partir da próxima iteração, conheceremos a função de adaptação não para as vértices dos simplexes, mas para os agentes de otimização específicos. E de acordo com as operações realizadas sobre os simplexes, designar soluções conhecidas dos agentes aos respectivos vértices do simplex. Precisamos atualizar a solução global com os melhores valores de adaptação dos agentes.

A seguir, verificamos qual operação foi realizada nas vértices na iteração anterior. Consideramos que os vértices dos simplexes são vetores.

Se a operação realizada foi "reflexão":

• atribuir ao vetor Xr o valor de adaptação do agente
• comparar a adaptação do vetor Xr com Xw, e se o valor for melhor, atualizar o vetor Xw
• comparar a adaptação do vetor Xr com Xb, e se o valor for melhor, realizar "expansão"
• comparar a adaptação do vetor Xr com Xg (a segunda após a pior vértice), e se o valor for pior, realizar "contração"
• se houver melhoria no vértice, realizar "reflexão", com ordenação prévia e construção do centroide, caso contrário, executar "Flying"
  

Se a operação realizada foi "expansão":

• atribuir ao vetor Xe o valor de adaptação do agente
• comparar a adaptação do vetor Xe com Xw, e se o valor for melhor, atualizar o vetor Xw
  
• realizar "reflexão", com ordenação prévia e construção do centroide
  

Se a operação realizada foi "contração":

• atribuir ao vetor Xс o valor de adaptação do agente
  
• comparar a adaptação do vetor Xс com Xw, e se o valor for melhor, atualizar o vetor Xw, caso contrário, executar "Flying"
• se houve melhoria no vértice, realizar "reflexão", com ordenação prévia e construção do centroide
  

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Revision ()
{
  //----------------------------------------------------------------------------
  if (!revision)
  {
    //sort agent simplex points by FF value and assign BGW
    for (int i = 0; i < popSize; i++)
    {
      Sorting (a [i].s.p);
    }

    //calculate the simplex centroid
    for (int i = 0; i < popSize; i++)
    {
      a [i].s.indG = simplexPoints - 2;
      a [i].s.indW = simplexPoints - 1;

      CalcCentroid (a [i].s, a [i].s.indW);
    }

    //assign the next type of operation - reflection
    for (int i = 0; i < popSize; i++)
    {
      Reflection (a [i], a [i].s.indW);
      a [i].s.operation = reflection;
    }

    //save the best point of the agents’ simplexes as a global solution 
    for (int i = 0; i < popSize; i++)
    {
      if (a [i].s.p [0].f > fB)
      {
        fB = a [i].s.p [0].f;
        ArrayCopy (cB, a [i].s.p [0].c, 0, 0, WHOLE_ARRAY);
      }
    }

    revision = true;
    return;
  }

  //----------------------------------------------------------------------------
  if (revision)
  {
    int pos = -1;

    for (int i = 0; i < popSize; i++)
    {
      if (a [i].f > fB) pos = i;
    }

    if (pos != -1)
    {
      fB = a [pos].f;
      ArrayCopy (cB, a [pos].c, 0, 0, WHOLE_ARRAY);
    }
  }

  //----------------------------------------------------------------------------
  for (int i = 0; i < popSize; i++)
  {
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    //if there was reflection ++++++++++++++++++++++++++++++++++++++++++++
    if (a [i].s.operation == reflection)
    {
      a [i].s.Xr.f = a [i].f;
      bool needUpd = false;

      //------------------------------------------------------------------------
      if (a [i].s.Xr.f > a [i].s.p [a [i].s.indW].f)  //Xr > Xw                 |---w--Xr--g----------b---|
      {
        a [i].s.p [a [i].s.indW] = a [i].s.Xr;        //replace Xw with Xr
        needUpd = true;
      }
      
      //------------------------------------------------------------------------
      if (a [i].s.Xr.f > a [i].s.p [0].f)             //if Xr > Xb            |---w----g----------b---Xr|
      {
        Expansion (a [i]);                            //perform expansion
        continue;
      }

      //------------------------------------------------------------------------
      if (a [i].s.Xr.f <= a [i].s.p [a [i].s.indG].f) //if Xr < Xg            |---w---Xr--g----------b--|
      {
        Contraction (a [i], a [i].s.indG);            //perform contraction
        continue;
      }
      
      if (needUpd)
      {
        Sorting (a [i].s.p);
        a [i].s.indG = simplexPoints - 2;
        a [i].s.indW = simplexPoints - 1;
        CalcCentroid (a [i].s, a [i].s.indW);
        Reflection   (a [i],   a [i].s.indW);
      }
      else Flying (a [i], a [i].s.indW);

      continue;
    }

    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    //if there was expansion +++++++++++++++++++++++++++++++++++++++++++++
    if (a [i].s.operation == expansion)
    {
      a [i].s.Xe.f = a [i].f;

      if (a [i].s.Xe.f > a [i].s.Xr.f) a [i].s.p [a [i].s.indW] = a [i].s.Xe;
      else                             a [i].s.p [a [i].s.indW] = a [i].s.Xr;

      Sorting (a [i].s.p);
      a [i].s.indG = simplexPoints - 2;
      a [i].s.indW = simplexPoints - 1;
      CalcCentroid (a [i].s, a [i].s.indW);
      Reflection   (a [i],   a [i].s.indW);

      continue;
    }

    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    //if there was contraction +++++++++++++++++++++++++++++++++++++++++++
    if (a [i].s.operation == contraction)
    {
      a [i].s.Xc.f = a [i].f;

      if (a [i].s.Xc.f > a [i].s.p [a [i].s.indW].f)
      {
        a [i].s.p [a [i].s.indW] = a [i].s.Xc;

        Sorting (a [i].s.p);
        a [i].s.indG = simplexPoints - 2;
        a [i].s.indW = simplexPoints - 1;
        CalcCentroid (a [i].s, a [i].s.indW);
        Reflection   (a [i],   a [i].s.indW);
      }
      else Flying (a [i], a [i].s.indW);

      continue;
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

Para calcular o centroide, escreveremos o método CalcCentroid, que calcula o valor médio das coordenadas em um dado simplex "s" até o índice indW especificado.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::CalcCentroid (S_Simplex &s, int indW)
{
  double summ = 0.0;

  for (int c = 0; c < coords; c++)
  {
    summ = 0.0;

    for (int p = 0; p < simplexPoints; p++)
    {
      if (p != indW) summ += s.p [p].c [c];
    }

    s.c [c] = summ / coords;
  }
}
//——————————————————————————————————————————————————————————————————————————————

A operação "reflexão" será executada pelo método "Reflection" para o vértice de índice "indW" do simplex do agente "agent".

Dentro do ciclo, o cálculo do valor refletido Xr é feito pela fórmula Xr = Xo + reflectionCoeff * (Xo - Xw). Aqui, reflectionCoeff é o coeficiente de reflexão, que define o grau de desvio do vértice refletido em relação ao vértice inicial.
Em seguida, a função SeInDiSp é aplicada ao valor de Xr para garantir que esteja dentro dos limites de faixa aceitáveis rangeMin[c] e rangeMax[c], com o passo rangeStep[c]. O resultado é salvo em agent.s.Xr.c[c].
Configurar a operação "reflexão" para o agente agent.s.operation, indicando que a operação "reflexão" foi realizada para esse agente.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Reflection (S_Agent &agent, int indW)
{
  double Xo;
  double Xr;
  double Xw;

  for (int c = 0; c < coords; c++)
  {
    Xo = agent.s.c [c];
    Xw = agent.s.p [indW].c [c];

    //Xr = Xo + RNDfromCI (0.0, reflectionCoeff) * (Xo - Xw);
    Xr = Xo + reflectionCoeff * (Xo - Xw);

    agent.s.Xr.c [c] = SeInDiSp  (Xr, rangeMin [c], rangeMax [c], rangeStep [c]);
    agent.c      [c] = agent.s.Xr.c [c];
  }

  agent.s.operation = reflection;
}
//——————————————————————————————————————————————————————————————————————————————

A operação "expansão" será executada pelo método "Expansion" para o agente "agent".

Dentro do ciclo, o cálculo do valor expandido Xe é feito pela fórmula Xe = Xo + expansionCoeff * (Xr - Xo). Aqui, expansionCoeff é o coeficiente de expansão, que define o grau de aumento do desvio do vértice expandido em relação ao vértice inicial.
Em seguida, a função SeInDiSp é aplicada ao valor de Xe para garantir que esteja dentro dos limites de faixa aceitáveis rangeMin[c] e rangeMax[c], com o passo rangeStep[c]. O resultado é salvo em agent.s.Xe.c[c].
Configurar a operação "expansão" para o agente agent.s.operation, indicando que a operação "expansão" foi realizada para esse agente.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Expansion (S_Agent &agent)
{
  double Xo;
  double Xr;
  double Xe;

  for (int c = 0; c < coords; c++)
  {
    Xo = agent.s.c    [c];
    Xr = agent.s.Xr.c [c];

    //Xe = Xo + RNDfromCI (0.0, expansionCoeff) * (Xr - Xo);
    Xe = Xo + expansionCoeff * (Xr - Xo);

    agent.s.Xe.c [c] = SeInDiSp  (Xe, rangeMin [c], rangeMax [c], rangeStep [c]);
    agent.c      [c] = agent.s.Xe.c [c];
  }

  agent.s.operation = expansion;
}
//——————————————————————————————————————————————————————————————————————————————

A operação "contração" será realizada pelo método "Contraction" para o vértice de índice "indW" do simplex do agente "agent".

Dentro do ciclo, o cálculo do valor comprimido Xc é feito pela fórmula Xc = Xo + contractionCoeff * (Xw - Xo). Aqui, contractionCoeff é o coeficiente de contração, que define o grau de redução do desvio do vértice comprimido em relação ao vértice inicial.
Em seguida, a função SeInDiSp é aplicada ao valor de Xc para garantir que esteja dentro dos limites de faixa aceitáveis rangeMin[c] e rangeMax[c], com o passo rangeStep[c]. O resultado é salvo em agent.s.Xc.c[c].
Configurar a operação "contraction" para o agente agent.s.operation, indicando que a operação "contração" foi realizada para esse agente.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Contraction (S_Agent &agent, int indW)
{
  double Xo;
  double Xw;
  double Xc;

  for (int c = 0; c < coords; c++)
  {
    Xo = agent.s.c    [c];
    Xw = agent.s.p [indW].c [c];

    //Xc = Xo + RNDfromCI (0.0, contractionCoeff) * (Xw - Xo);
    Xc = Xo + contractionCoeff * (Xw - Xo);

    agent.s.Xc.c [c] = SeInDiSp  (Xc, rangeMin [c], rangeMax [c], rangeStep [c]);
    agent.c      [c] = agent.s.Xc.c [c];
  }

  agent.s.operation = contraction;
}
//——————————————————————————————————————————————————————————————————————————————

Para os voos de Lévy, aplicamos o método "Flying", que realiza a operação "voo" para o agente "agent" com o objetivo de empurrar o vértice do simplex para novos locais do espaço de busca com uma distribuição aleatória que concentra a probabilidade mais próxima da melhor solução global e diminui a probabilidade até zero nas fronteiras do espaço explorado.

Após a aplicação do método "Flying", configurar a operação "reflection" para o agente agent.s.operation, para que nas iterações subsequentes a lógica do algoritmo seja executada como se uma operação de reflexão tivesse sido realizada.
A função SeInDiSp é aplicada ao valor de Xr para garantir que ele esteja dentro dos limites de faixa aceitáveis rangeMin[c] e rangeMax[c], com o passo rangeStep[c]. O resultado é salvo em agent.s.Xr.c[c].

Assim, o método "Flying" realiza a operação "voo" para o agente, alterando suas coordenadas com base nas coordenadas atuais da melhor solução global e valores gerados aleatoriamente. Isso permite que o agente explore novas áreas do espaço de soluções e potencialmente encontre soluções mais ótimas.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_NMm::Flying (S_Agent &agent, int indW)
{
  for (int c = 0; c < coords; c++)
  {
    double r1 = RNDfromCI (-1.0, 1.0);
    r1 = r1 > 0.5 ? 1.0 : -1.0;
    double r2 = RNDfromCI (1.0, 20.0);
    r2 = pow (r2, -2.0);

    if (r1 > 0.0) Xr = cB [c] + Scale (r2, 0.0, 1.0, 0.0, rangeMax [c] - cB [c], false);
    else          Xr = cB [c] - Scale (r2, 0.0, 1.0, 0.0, cB [c] - rangeMin [c], false);

    //Xr = RNDfromCI (rangeMin [c], rangeMax [c]);

    agent.s.Xr.c [c] = SeInDiSp  (Xr, rangeMin [c], rangeMax [c], rangeStep [c]);
    agent.c      [c] = agent.s.Xr.c [c];
  }

  agent.s.operation = reflection;
}
//——————————————————————————————————————————————————————————————————————————————

3. Resultados dos testes

Impressão dos resultados do algoritmo IWD na bancada de testes:

C_AO_NMm:5;1.0;2.0;0.5
=============================
5 Rastrigin's; Func runs 10000 result: 77.96849465506082
Score: 0.96607
25 Rastrigin's; Func runs 10000 result: 61.68192953373487
Score: 0.76427
500 Rastrigin's; Func runs 10000 result: 50.62713783555849
Score: 0.62730
=============================
5 Forest's; Func runs 10000 result: 0.9552378700292385
Score: 0.54033
25 Forest's; Func runs 10000 result: 0.45877475613538293
Score: 0.25951
500 Forest's; Func runs 10000 result: 0.09902427974784639
Score: 0.05601
=============================
5 Megacity's; Func runs 10000 result: 5.6
Score: 0.46667
25 Megacity's; Func runs 10000 result: 2.304
Score: 0.19200
500 Megacity's; Func runs 10000 result: 0.40280000000000005
Score: 0.03357
=============================
All score: 3.90573

Os resultados do trabalho do algoritmo mostram resultados bastante dignos na função suave de Rastrigin.

Na visualização, podemos notar sinais de degeneração da população e concentração de agentes em um único lugar. A situação é um pouco resgatada pelo uso de voos de Lévy, o que é evidente por pontos individuais correspondentes às coordenadas. Os gráficos de convergência não são muito estáveis; eu apresentei a visualização apenas dos dois primeiros testes para cada função, porque a execução do terceiro teste com 1000 variáveis é tão longa que não há possibilidade de gravação em GIF.

  NMm na função de teste Rastrigin.

  NMm na função de teste Forest.

  NMm na função de testeMegacity.

O método de Nelder-Mead, apesar de suas numerosas desvantagens, ocupou o honroso nono lugar.

#	AO	Description	Rastrigin			Rastrigin final	Forest			Forest final	Megacity (discrete)			Megacity final	Final result
			10 p (5 F)	50 p (25 F)	1000 p (500 F)		10 p (5 F)	50 p (25 F)	1000 p (500 F)		10 p (5 F)	50 p (25 F)	1000 p (500 F)
1	SDSm	stochastic diffusion search M	0.99809	1.00000	0.69149	2.68958	0.99265	1.00000	1.00000	2.99265	1.00000	1.00000	1.00000	3.00000	100.000
2	SSG	saplings sowing and growing	1.00000	0.92761	0.51630	2.44391	0.72120	0.65201	0.83760	2.21081	0.54782	0.61841	0.99522	2.16146	77.683
3	DE	differential evolution	0.98295	0.89236	0.51375	2.38906	1.00000	0.84602	0.65510	2.50112	0.90000	0.52237	0.12031	1.54268	73.099
4	HS	harmony search	0.99676	0.88385	0.44686	2.32747	0.99148	0.68242	0.37529	2.04919	0.71739	0.71842	0.41338	1.84919	70.623
5	IWO	invasive weed optimization	0.95828	0.62227	0.27647	1.85703	0.70170	0.31972	0.26613	1.28755	0.57391	0.30527	0.33130	1.21048	48.250
6	ACOm	ant colony optimization M	0.34611	0.16683	0.15808	0.67103	0.86147	0.68980	0.64798	2.19925	0.71739	0.63947	0.05579	1.41265	47.387
7	MEC	mind evolutionary computation	0.99270	0.47345	0.21148	1.67763	0.60244	0.28046	0.21324	1.09615	0.66957	0.30000	0.26045	1.23002	44.049
8	SDOm	spiral dynamics optimization M	0.69840	0.52988	0.33168	1.55996	0.59818	0.38766	0.37600	1.36183	0.35653	0.15262	0.25842	0.76758	40.289
9	NMm	Nelder-Mead method M	0.88433	0.48306	0.45945	1.82685	0.46155	0.24379	0.21903	0.92437	0.46088	0.25658	0.16810	0.88555	39.660
10	COAm	cuckoo optimization algorithm M	0.92400	0.43407	0.24120	1.59927	0.57881	0.23477	0.13842	0.95200	0.52174	0.24079	0.17001	0.93254	37.830
11	FAm	firefly algorithm M	0.59825	0.31520	0.15893	1.07239	0.50637	0.29178	0.41704	1.21519	0.24783	0.20526	0.35090	0.80398	33.139
12	ABC	artificial bee colony	0.78170	0.30347	0.19313	1.27829	0.53379	0.14799	0.11177	0.79355	0.40435	0.19474	0.13859	0.73768	29.766
13	BA	bat algorithm	0.40526	0.59145	0.78330	1.78002	0.20664	0.12056	0.21769	0.54488	0.21305	0.07632	0.17288	0.46225	29.499
14	CSS	charged system search	0.56605	0.68683	1.00000	2.25289	0.13961	0.01853	0.13638	0.29452	0.07392	0.00000	0.03465	0.10856	27.930
15	GSA	gravitational search algorithm	0.70167	0.41944	0.00000	1.12111	0.31390	0.25120	0.27812	0.84322	0.42609	0.25525	0.00000	0.68134	27.807
16	BFO	bacterial foraging optimization	0.67203	0.28721	0.10957	1.06881	0.39364	0.18364	0.17298	0.75026	0.37392	0.24211	0.18841	0.80444	27.542
17	EM	electroMagnetism-like algorithm	0.12235	0.42928	0.92752	1.47915	0.00000	0.02413	0.29215	0.31628	0.00000	0.00527	0.10872	0.11399	19.002
18	SFL	shuffled frog-leaping	0.40072	0.22021	0.24624	0.86717	0.19981	0.02861	0.02221	0.25063	0.19565	0.04474	0.06607	0.30646	13.200
19	MA	monkey algorithm	0.33192	0.31029	0.13582	0.77804	0.09927	0.05443	0.07482	0.22852	0.15652	0.03553	0.10669	0.29874	11.777
20	FSS	fish school search	0.46812	0.23502	0.10483	0.80798	0.12730	0.03458	0.05458	0.21647	0.12175	0.03947	0.08244	0.24366	11.332
21	IWDm	intelligent water drops M	0.26459	0.13013	0.07500	0.46972	0.28358	0.05445	0.05112	0.38916	0.22609	0.05659	0.05054	0.33322	10.423
22	PSO	particle swarm optimisation	0.20449	0.07607	0.06641	0.34696	0.18734	0.07233	0.18207	0.44175	0.16956	0.04737	0.01947	0.23641	8.426
23	RND	random	0.16826	0.09038	0.07438	0.33302	0.13381	0.03318	0.03949	0.20648	0.12175	0.03290	0.04898	0.20363	5.054
24	GWO	grey wolf optimizer	0.00000	0.00000	0.02093	0.02093	0.06514	0.00000	0.00000	0.06514	0.23478	0.05789	0.02545	0.31812	1.000

Considerações finais

O algoritmo de Nelder-Mead é pouco aplicável em tarefas modernas e complexas de otimização, pois exige o cálculo da função de adaptação para cada vértice do simplex na fase inicial, o que anula sua capacidade de aplicação em espaços de busca multidimensionais. Daí que o algoritmo não possa ser recomendado para uso, apesar de seus resultados relativamente bons e de sua posição na tabela de classificação.

O algoritmo SDS foi removido da tabela, deixando apenas sua versão modificada.

Figura 2. Gradiente de cor dos algoritmos de acordo com os testes correspondentes.

Figura 3. Histograma dos resultados dos testes dos algoritmos (em uma escala de 0 a 100, quanto maior, melhor,

com um script no arquivo para calcular a tabela de classificação com base na metodologia deste artigo).

Vantagens e desvantagens da versão modificada do algoritmo de Nelder-Mead (NMm):

Vantagens:

1. Pequeno número de parâmetros externos.
2. Resultados razoáveis em funções suaves com poucas variáveis.
   

Desvantagens:

1. Implementação complexa.
2. Dificuldade de uso com aplicativos de terceiros devido à necessidade de cálculo prévio da adaptação para todos os vértices do simplex.
3. Baixa velocidade de operação.
4. Muito pobre em escalabilidade.
5. Alta variação nos resultados.
6. Tendência para ficar preso em extremos locais.

Note-se que as Vantagens e Desvantagens se referem especificamente à versão modificada da busca simplex, enquanto a versão comum do NM nem vale a pena ser incluída na tabela devido aos resultados muito fracos.

O artigo inclui um arquivo com versões atualizadas dos códigos dos algoritmos descritos em artigos anteriores. O autor do artigo não assume responsabilidade pela precisão absoluta na descrição dos algoritmos canônicos, muitos dos quais foram modificados para melhorar as capacidades de busca. As conclusões e opiniões apresentadas nos artigos se baseiam nos resultados dos experimentos realizados.