Implementando PSO em R

Seguindo a linha de implementações de algoritmos evolutivos e de enxame em R, vamos implementar um PSO (Particle Swarm Optimization), em português conhecido como Otimização por Nuvem de Partículas. Para não criar confusão com nome em português ou inglês, irei sempre me referir ao algoritmo como PSO.

Antes de começarmos, se você não sabe nada de R então sugiro a leitura do artigo “O que você precisa saber sobre linguagem R para implementar algoritmos evolutivos”. Caso queira entender o funcionamento dos Algoritmos Genéticos antes de começar o PSO, então leia o artigo “Implementando Algoritmos Genéticos em R”.

Então vamos começar…

O PSO não é exatamente um algoritmo evolutivo, mas sim um algoritmo de enxame. Essa linha é comumente chamada de inteligência de enxame ou inteligência coletiva, abrangendo algoritmos tais como PSO, Colônia de Formigas (Ant Colony), Colonia de Abelhas Artificias (Artificial Bee Colony), Alimentação Bacteriana (Bacterial Foraging), dentre outros. A principal diferença desses algoritmos com relação aos evolutivos é que os de enxame simulam o comportamento de sociedades de indivíduos, normalmente no processo de busca por alimentação.

Proposto por J. Kenndy e R. Eberhart em 1995 em seu artigo “Particle Swarm Optimization” [1], o PSO tornou-se um dos algoritmos mais importantes no mundo da otimização contínua, possuindo também variações para otimização multiobjetivo, na qual dois ou mais objetivos conflitantes devem ser otimizados. Seu objetivo é simular o comportamento de pássaros em busca de alimentos e embora inicialmente criado para simulação comportamental, rapidamente foi aplicado na otimização numérica.

Bom, vamos ao algoritmo.

Ao inicializarmos o enxame, três importantes matrizes devem ser criadas: o enxame propriamente dito (S), que representa a posição das partículas no espaço de busca, o histórico da melhor posição por onde a partícula já passou (P) e a velocidade (V), sendo que fitness é o vetor que contem a avaliação de Se P.fit é o vetor que contém o fitness da melhor posição já visita por cada partícula. Com essas estruturas definidas cria-se ainda o vetor g que armazena a melhor posição até o momento e g.fit que é seu fitness.

Após o processo de inicialização começa-se as iterações do algoritmo, que equivale às gerações dos algoritmos genéticos. Enquanto o critério de parada não é alcançado, atualiza-se a matriz de velocidade usando w (peso de inércia), c1 e c2 (constantes aceleradoras), e, r1 e r2 (matrizes de variáveis aleatórias entre 0 e 1). Observe que a velocidade utiliza os dois melhores fitness encontrados para direcioná-la para em seguida atualizar as posições. Após a atualização da posição de cada partícula verifica-se se o melhor fitness (g.fit) foi melhorado e se as novas posições são melhores do que o histórico das partículas.

Já acabou? Sim, o algoritmo é bastante simples. Então vamos por a mão na massa em R.

A função de inicialização é mostrada a seguir, na qual são inicializadas as matrizes S, P e V. Um detalhe que não foi tratado no algoritmo são os parâmetros v.min e v.max que tem como objetivo controlar a velocidade das partículas, pois elas podem crescer indefinidamente, o que gera a chamada explosão de velocidade, que faz com que o algoritmo perca a capacidade de exploitation, ou seja, de busca local.

O resto da implementação é mostrada na figura abaixo, na qual faz-se uma alocação de memória para as matrizes e se inicializam as estruturas necessárias. Já nas iterações do algoritmo, observe como a atualização tanto das velocidades quanto da posição das partículas é feito em instruções matriciais. Em seguida verifica-se se nenhuma dimensão da partícula ultrapassou os limites inferior e superior, respectivamente, fazendo a mesma coisa na matriz de velocidades. Aqui é importante introduzir um novo parâmetro da função which() que é o arr.ind = TRUE. Esse parâmetro devolve a posição da matriz em que a condição é verdadeira. Se fosse usada a instrução which(S < lb), sem o parâmetro arr.ind, seriam mostrados os valores que ferem essa condição e não os índices. Os índices também formam uma matriz, e por lb e ub serem vetores, é mandatório usar a segunda posição dos índices que contém a coluna que não obedece a condição (idx[,2]). Na verdade, a operação S[idx] <- lb[idx[,2]] obriga que a quantidade de elementos seja a mesma em ambos os lados. A operação é feita se a quantidade não for igual, mas irá gerar um warning que não é facilmente detectado durante a execução, e consequentemente, irá gerar resultados errados. Já quando os limites da velocidade são verificados pode-se utilizar a indexação lógica simples, pois v.min e v.max são escalares.

O script de teste é mostrado a seguir. Nele são definidas os parâmetros de execução do PSO, sendo w um número menor que 0.9, e, c1 e c2 normalmente igual a 2.

Ainda é possível implementar um decaimento linear para a variável de inércia w. Para isso é necessário trocar w por dois novos parâmetros w.min e w.max, normalmente, 0.9 e 0.4, respectivamente, e computar w através da equação:

O PSO é um algoritmo bastante simples de se implementar, especialmente em R. Espero que tenham gostado e comecem a praticar esse conhecimento o mais rápido possível. Para saber mais…

[1] Kennedy, J.; Eberhart, R. (1995). “Particle Swarm Optimization”. Proceedings of IEEE International Conference on Neural Networks. IV. pp. 1942–1948.

[2] Shi, Y.; Eberhart, R.C. (1998). “A modified particle swarm optimizer”. Proceedings of IEEE International Conference on Evolutionary Computation. pp. 69–73.

[3] Kennedy, J. (1997). “The particle swarm: social adaptation of knowledge”. Proceedings of IEEE International Conference on Evolutionary Computation. pp. 303–308.