Implementado Algoritmos Genéticos em R

A primeira pergunta que pode passar pela sua cabeça é “por que em R?”. Bom, vários são os motivos. Vou tentar enumerar os principais deles que consigo me lembrar. Primeiro, é uma linguagem que contém inúmeras facilidades para se trabalhar com vetores e matrizes, estruturas de dados estas que são a base para a implementação de algoritmos evolutivos e de enxame. Segundo, é uma linguagem que possui mais de 10.000 pacotes em seu repositório (CRAN) em 2017, isso mesmo, mais de dez mil pacotes, que variam de aplicações em ecologia, passando por equações diferenciais, aprendizagem de máquina, dentre muitas outros. Terceiro, é uma linguagem pronta para estatística, ou seja, já possui todas as funções necessárias para realizar comparações estatísticas entre algoritmos. E quarta, mas não menos importante, é completamente free.

Antes de começar a mostrar a implementação, é importante que o leitor esteja seguro que conhece quatro importantes conceitos de R que serão amplamente usados neste artigo: vetores, matrizes, listas e funções de grupo. Acessar vetores e matrizes em R é muito semelhante ao que se faz em linguagens como C e Java, porém com as mesmas facilidades apresentadas pelo Matlab, tais como, operações diretas com vetores e matrizes, e indexação lógica. Listas são como objetos que armazenam diferentes tipos de dados. Por exemplo, uma lista pode ser formada por três estruturas de dados diferentes, tais como, um vetor, uma matriz e uma variável inteira, não necessariamente nessa ordem. Muito embora seja possível trabalhar com R de forma orientada a objetos, neste artigo iremos trabalhar de forma procedimental. Finalmente, funções de grupos são funções que permitem ao programador executar operações em vetores e matrizes sem a necessidade de utilizar laços tais como for e while. Na verdade, o uso de funções de grupo são preferíveis ao uso de laços, pois elas normalmente são implementadas de forma otimizada. Além disso, estou assumindo também que o leitor já possua um certo conhecimento em algoritmos.

Então vamos começar!

Neste momento você deve estar se perguntando, mas o que é um Algoritmo Genético (AG)? Formalmente é uma meta-heurística baseada no processo de seleção natural na qual o indivíduo que melhor se adapta ao ambiente é o que sobrevive.

What?

Calma. Vamos passo a passo.

Imagine que em um deserto existem duas espécies de animais: pequenos répteis, que vou chamar de lagartixas, e uma ave de rapina, vamos supor que seja uma águia. Se as lagartixas forem, por exemplo, de cor avermelhada, e a velocidade de reprodução delas não for muito grande, em pouco tempo as águias irão extinguí-las. Agora imagine que por uma mutação genética uma lagartixa nasceu bege, ou seja, mais próximo da cor da areia do deserto. Se essa lagartixa não se mexer enquanto esta sobre a areia, a águia não conseguirá vê-la. Na verdade eu não sei se a águia irá ou não vê-la, mas vamos supor que não. Nesse caso, a nova lagartixa terá mais chances de sobreviver, e como consequência de gerar novos descendentes. Todos nós sabemos que ao gerar descendentes, estes terão parte do material genético dos geradores. Isso significa que há uma probabilidade de que as lagartixas descendentes da lagartixa bege sejam da mesma cor. Assim, essa mutação irá equilibrar as coisas entre as lagartixas e as águias, até que uma águia sofra uma mutação que permita visualizar sua presa.

Então, tomando como base o exemplo das lagartixas e águias, podemos começar a entender conceitos como seleção, cruzamento e mutação, que são as principais operações dos AGs. O processo de seleção natural é aquele em que os mais aptos tem maior probabilidade de sobreviver e transmitir seus genes. Essa adaptação se deu, no exemplo, pela mutação quando a lagartixa mudou de cor. Já, quando as lagartixas geraram novos descendentes através do compartilhamento genético houve o que se conhece como cruzamento.

Neste ponto, já temos quase toda a teoria necessária para começarmos a entender como funciona um AG. Mas calma, ainda precisamos de mais um pouco de teoria. Na verdade, precisamos definir como simular esse processo em um algoritmo que possa ser executado em um computador.

Segundo Norvig, um AG é um algoritmo de busca local, no qual a geração de um novo estado não depende seu histórico, mas somente do seu estado atual. Conceitualmente isso esta correto. No entanto, no mundo da computação evolutiva, os AGs são conhecidos como algoritmos de busca global, pois seu objetivo é fazer uma busca geral no que conhecemos como espaço de busca. Mas em que isso muda os conceitos envolvidos em AGs? Boa pergunta! No mundo da computação evolutiva também existe o conceito de busca local, que é a busca por soluções em uma região próxima ao estado atual. É desses conceitos que derivam os conceitos conhecidos em inglês como Exploration e Exploitation. O primeiro refere-se a capacidade de busca global, enquanto que o segundo refere-se a capacidade busca local. Nesse contexto, é interessante que o algoritmo tenha uma boa capacidade de Exploration em grande parte de sua execução; e uma boa capacidade de Exploitation na parte final de sua execução, já que imagina-se que no final o algoritmo já se encontre próximo da solução ótima.

Mas que diabos esse AG esta buscando?

AGs buscam soluções para um problema. Imagine que você quer encontrar o menor caminho saindo de uma cidade qualquer, passando por várias cidades, como por exemplo, São Luis, Bacabal, Santa Inês, Imperatriz, Açailandia, Carolina, Santa Luzia e Arari. Não seria inteligente, por exemplo, fazer a rota inicial São Luis — Açailandia — Santa Inês — Imperatriz — Arari, pois se gastaria muito combustível e tempo. Por outro lado, a sequencia inicial São Luis —Arari — Santa Inês — Santa Luzia — Açailandia seria mais adequada. Nesse exemplo, o espaço de busca seria a combinação de todas as cidades. Matematicamente, seria a combinação das 8 cidades 8-a-8, onde não pode haver repetições. Esse problema é conhecido como o problema do caixeiro viajante e nesse exemplo em particular gera um espaço de busca de 40320 possibilidades de combinações.

Algoritmicamente, sabendo o aspecto que se deseja avaliar, com por exemplo a distância percorrida, é possível encontrar a solução para o problema das cidades usando força bruta, ou seja, testando todas as possibilidades. Devemos lembrar que o exemplo tem apenas 8 cidades. Quando o problema passa a ter mais cidade, 21 por exemplo, o espaço de busca aumenta para 51.090.942.171.709.440.000 possibilidades. Literalmente cinquenta e um quintilhões de possibilidades. Agora imagine que você tem um processador de 3 Ghz, o que dá aproximadamente 3.000.000.000 de operações por segundo (aqui não estou considerando pipeline, paralelismo e outras características que podem deixar o processamento mais rápido), para terminar a avaliação de todas as possibilidades são necessários aproximadamente 4.730.643 horas ou 197 mil dias ou 540 anos, isto é, inviável.

Agora sabendo que podemos precisar de um AG, é necessário definir o que seria um indivíduo. Seria uma lagartixa? Isso, mas precisamos de algo mais preciso ou mais formal. No caso das cidades, um indivíduo seria uma sequencia de cidades a serem visitadas. Se as cidades forem numeradas na sequencia 1-São Luis, 2-Bacabal, 3-Santa Inês, 4-Imperatriz, 5-Açailandia, 6-Carolina, 7-Santa Luzia e 8-Arari, o vetor a seguir é um individuo válido, isto é, uma possível solução.

Ah entendi. Então se uma solução é um indivíduo, cada cidade é um gene? Exatamente! Desse conceito derivamos os AGs binários, cujos genes são formados apenas por 0s e/ou 1s, e os AGs em código real, no qual os genes podem ser números inteiros ou reais. E a partir da representação podemos introduzir os operadores de cruzamento e de mutação.

O operador de cruzamento mais comum é o cruzamento simples. Nesse cruzamento, dois genitores e um ponto de corte são escolhidos. Em seguida são gerados dois filhos trocando-se os genes de cada pai. Considerando então dois pais com código real, poderíamos ter algo como apresentado na figura a seguir (o mesmo processo serve para indivíduos com código binário). Não estamos usando o exemplo das cidades porque esse exemplo em particular exige a utilização de operadores especiais que serão vistos em outro post.

Com o cruzamento feito, agora podemos mostrar o que é a mutação. No caso de um AG em código real, a mutação mais simples, também chamada de uniforme, é feita gerando-se um número aleatório dentro do domínio de cada gene, isto é, um número é gerado considerando os limites inferior e superior do gene. No caso do AG ser binário, basta trocar 0 por 1 e vice-versa, como mostrado a seguir.

Agora já sabemos como funcionam os principais operadores, portanto já podemos apresentar o algoritmo e começar a colocar a mão na massa. Mas calma, primeiro precisamos definir qual o problema que vamos resolver. Suponhamos que o problema a ser resolvido é minimizar a função:

Essa função é bem conhecida no mundo da computação evolutiva, chamada de Rosenbrock, sendo bastante difícil de solucionar em grandes dimensões. Em duas dimensões, o objetivo é encontrar o ponto mínimo da figura a seguir:

Considerando a função Rosenbrock em duas dimensões, resolver esse problema é relativamente fácil. Mas você não acabou de dizer que é difícil? Sim, em grandes dimensões. O gráfico mostrado representa apenas duas dimensões, ou variáveis. Além disso, como pode ser observado, essa função é genérica, na qual d representa uma dimensão qualquer, normalmente 30, tornando bastante difícil o processo de encontrar sua solução ótima. O domínio de cada gene no nosso código deve estar no intervalo [-5,5], em testes com o código você pode definir outros limites. O mínimo dessa função é um vetor de 1s, ou seja, [1,…,1]. Como estamos considerando uma dimensão igual a 30, o mínimo será um vetor de 30 posições.

Rosenbrock <- function(x){
  d <- length(x)
  xi <- x[1:(d-1)]
  xnext <- x[2:d]
  sum <- sum(100*(xnext-xi^2)^2 + (xi-1)^2)
  return(sum)
}

Com esses dados nas mãos já podemos começar a nossa implementação. Eu costumo usar o RStudio, que é uma IDE open source que tem o apoio da comunidade de R, porém essa IDE não é obrigatória, pois existem outras como, por exemplo, a Microsoft R (antes conhecida como Revolution R) que é free para estudantes.

Vamos ao código? Let’s do this!

Vamos começar apresentando o algoritmo canônico do AG. Para colocar logo a mão na massa irei explicando o algoritmo juntamente com o código em R

O código para testar o AG é apresentado a seguir, no qual carregam-se todas as funções que serão usadas pelo código através da instrução source. Esses códigos são os mesmos apresentados no restante deste artigo e possuem o mesmo nome das funções que implementam. Em seguida definem-se a dimensão do problema (dim), a quantidade de iterações (it) que é o critério de parada do AG, os limites inferior e superior (lb e ub), a função que será otimizada (func), a quantidade de execuções (execs), um vetor que receberá uma lista de resultados e um vetor que conterá a melhor solução encontrada em cada execução. Para o AG ainda são passadas: a probabilidade de cruzamento (pc), a probabilidade de mutação (pm) e o tipo de seleção que será usado (sel = 1 indica roulette wheel; sel = 2 indica tournament)

A função GA é muito similar ao algoritmo canônico, mas com dois parâmetros a mais, que são o tamanho do torneio (para seleção via torneio) e se o elitismo será ou não feito. Como o parâmetros t.size não estão sendo passado pelo código de teste, adota-se o padrão definido no código (t.size = 4) do AG propriamente dito (função GA no código).

No AG canônico quando uma nova população é criada ela substitui a população corrente. Em alguns casos, a nova população não contém nenhum individuo melhor que a população corrente, o que faz com que esse individuo seja perdido. O que por um lado pode diminuir um pouco a pressão seletiva, mas por outro lado, pode fazer a população retroceder em termos de evolução. Quando o elitismo é usado, garante-se que pelo menos o melhor individuo passe de uma geração (iteração) para outra.

Quando a função GA é chamada, o primeiro passo do algoritmo é executado, que é inicializar a população de modo aleatório, mas dentro do domínio de cada gene. Isso é feito pela função init.population(). Observe que apesar de precisar criar vários indivíduos e cada um com 30 genes, nenhuma instrução de laço é utilizada. Note também a função apply que aplica a função Rosenbrock para cada linha da matriz criada. Esse processo é chamado de cálculo de aptidão ou fitness, e quando feita em toda população determina-se a função de Avalia População() do AG canônico.

O próximo passo é entrar no laço onde as gerações são executadas, sendo a primeira função a de seleção. O processo de seleção visa escolher dentro da população quais individuos devem ser cruzados. Isso pode ser feito de diversas formas, sendo que na maioria deles utiliza-se como base a aptidão dos indivíduos. A função selection() apenas escolhe se a seleção será feita via roleta (roulette wheel) ou competição (tournament).

Vamos começar estudando o algoritmo de seleção por roleta, no qual para cada indivíduo é dada um probabilidade de ser escolhido. Essa probabilidade é obtida através da divisão de cada aptidão pela soma de todas as aptidões da população. Em seguida é computada a probabilidade cumulativa de cada indivíduo de modo a poder selecionar os indivíduos por sorteio. É como se para cada indivíduo fosse dado um pedaço de uma roleta de acordo com sua aptidão. Assim, quanto maior a aptidão, maior o pedaço da roleta. O processo é mostrado no algoritmo a seguir.

Como se pode observar, cada individuo recebe um pedaço maior da roleta se tiver uma aptidão maior. Isso faz com que o método da roleta seja adequado para problemas de maximização. Por outro lado, um artifício matemático permite que um problema de maximização seja transformado em um problema de minimização, pois Max f(x) = Min -f(x). Nesse caso, o código em R para seleção via roleta fica como mostrado logo abaixo.

Observe no código da roleta como a probabilidade é computada a partir do valor negativo do fitness. Já na seleção por torneio, o processo de seleção é bem mais simples como mostra a função tournament(). Na função, quatro individuos são escolhidos, pois t.size = 4, para competirem entre si. Aquele que dentre os quatro escolhidos tiver a menor aptidão (melhor) é selecionado para o cruzamento.

Em seguida toda a população intermediária escolhida é passada para a função de cruzamento (simple.crossover()) propriamente dita. No entanto, nem todos serão cruzados. Isso vai ser feito de acordo com a probabilidade de cruzamento (Pc), que normalmente é um valor entre 0.6 e 0.99, indicando a porcentagem média de indivíduos que serão cruzados por iteração. Especificamente nesta implementação estou usando Pc = 0.6, que implica que em média 60% dos indivíduos selecionados serão cruzados. Se o cruzamento será feito em ordem de posicionamento no vetor, aleatoriamente, ou de uma outra forma é uma decisão do programador. Eu decidi implementar a escolha aleatória dos genitores como pode ser visto na função simple.crossover().

Após realizar o cruzamento será necessário realizar a mutação, que é baseada na probabilidade de mutação (Pm), sendo normalmente um valor entre 0.01 e 0.05, ou seja, entre 1% e 5% de probabilidade de ser mutado. Aqui estou usando Pm = 0.05. A função uniform.mutation() realiza a mutação em todos os genes sem usar uma única instrução de laço. Após sua execução, utiliza-se a função apply() para avaliar a população resultante.

Após a mutação chegamos no ponto onde o elitismo deve ou não ser utilizado. Como eu passei o parâmetro na função GA (elitism = TRUE), o mesmo será utilizado. No trecho de código apenas verifico se o melhor pai (geração corrente) é melhor que o melhor filho (nova geração). Em caso afirmativo, eu substituo o pior filho (que apresenta a pior aptidão) pelo melhor pai. Caso contrário não é necessário realizar a substituição. Assim, o processo é repetido até atingir o critério de parada.

Vamos executar o código em 4 versões: Roleta com elitismo, roleta sem elitismo, torneio com elitismo e torneio sem elitismo e ver os resultados para a função Rosenbrock em 5 execuções de cada, e usando os parâmetros previamente definidos.

Nitidamente, o elitismo fez diferença nos dois métodos de seleção, isto é, em ambos casos os resultados com elitismo foram melhores. Já em uma comparação superficial, o AG usando o método do torneio com elitismo apresentou os melhores resultados. Para uma comparação adequada métodos estatísticos deveriam ser utilizados, mas isso é assunto para outro post.

Bem chegamos ao fim. Espero que tenham gostado e que esse código seja útil para seus estudos e aprendizagem. A seguir indico alguns livros que podem ajudar a aumentar seus conhecimentos em AGs.

Este foi meu primeiro contato com Algoritmos Genéticos. Os primeiros capítulos são bem detalhados e mostram como um AG binário funciona passo a passo. Em seguida Michalewicz discute os motivos pelos quais os AGs funcionam e discute variações no algoritmo. A parte II é direcionada para a otimização numérica, na qual ele apresenta AGs em código real e como lidar com problemas com restrições. Possui um anexo com a implementação de um AG com código real em C.

O livro de Ricardo Linden discute o funcionamento de AGs de maneira bastante aprofundada, descrevendo diferentes métodos de seleção, cruzamento e mutação. Há um capítulo sobre a teoria de schemas e como eles estão presentes na evolução dos AGs binários. Além disso, ele discute os problemas que os operadores genéticos podem acarretar em determinados tipos de problemas. O livro todo apresenta implementações em Java que podem ser usadas para entender melhor os conceitos por ele apresentados.

Os códigos funcionando podem ser encontrados no endereço: https://gitlab.com/omar.carmona/real-coded-genetic-algorithm

Se usar os códigos por favor citar o seguinte artigo: “ Cortes, O. A. C. and Silva, J. C. (2019), Unconstrained numerical optimization using real-coded genetic algorithms: a study case using benchmark functions in R from Scratch. Revista Brasileira De Computação Aplicada, 11(3), 1–11”.