Programação Evolutiva em R

Continuando os posts relacionados à programação de algoritmos evolutivos em R, vou mostrar como implementar a programação evolutiva em R. Se você não conhece R tire uns minutos para ler o artigo “O que você precisa saber sobre linguagem R para implementar algoritmos evolutivos”. Para outras implementações visite o “Semper Evolutionis”.

Então mãos a obra!

A programação evolutiva (PE) foi proposta por Fogel em 1966 [1] e desde então tem sido aplicada com sucesso principalmente na resolução de problema de otimização numérica. Seu algoritmo é relativamente simples como podemos ver no pseudocódigo a seguir:

O algoritmo da PE é bastante semelhante a outros algoritmos evolutivos. E como evolutivo, deve-se criar uma população inicial aleatória no domínio especificado pelo problema. Porém, aqui há uma diferença se comparado com, por exemplo, algoritmos genéticos (AGs). Na PE um indivíduo é formado pela sua posição no espaço de busca x mais um vetor de desvios padrão n. O vetor de desvios padrão será a base para modificar a posição x no espaço de busca. O código em R para inicializar a população é:

init.population <- function(pop.size,dimension,lb,ub){
   pop <- matrix(runif(pop.size*dimension),nrow=pop.size)
   pop <- lb + pop * (ub-lb)
   std.dev <- matrix(runif(pop.size*dimension),nrow=pop.size)
   return(list(pop = pop,std.dev = std.dev))
}

Ao entrar no laço cada indivíduo deve gerar um filho de acordo com a equação abaixo, na qual a nova posição no espaço sera formada pela posição atual mais o desvio padrão multiplicado por um número aleatório de Cauchy ou pode-se usar também um número aleatório normal com média zero e desvio padrão igual a 1, sendo que o desvio padrão também será mutado. Vale destacar que Nj(0,1) refere-se a um número aleatório para cada gene enquanto que N(0,1) refere-se a um único número aleatório por indivíduo. Nas constantes tau e tau’, d indica a quantidade de genes ou a dimensão do problema.

O gráfico abaixo mostra o que seria um número aleatório normal com essas características, ou seja, de acordo com essa função haverá uma maior probabilidade de se gerar números aleatórios mais ou menos entre -2 e 2.

Após avaliar os filhos (calcular a aptidão dos novos indivíduos) deve-se fazer a união entre os filhos e os genitores com o objetivo de realizar um torneio, ou seja, cada indivíduo sera comparado com outros q indivíduos selecionados aleatoriamente dentro da população unida. Para cada indivíduo computa-se quantas vitórias cada um teve. A nova geração será formada pelos indivíduos que obtiveram o maior número de vitórias. Segue o código em R para o restante do algoritmo:

EP <- function(func,pop.size, dimension,lb,ub,num.it,q){
   tau  <- sqrt(2 * sqrt(pop.size))^-1
   taul <- sqrt(2 * pop.size)^-1
   pop <- init.population(pop.size, dimension,lb,ub) 
   x <- pop$pop
   std.dev <- pop$std.dev
   fit <- apply(x,1,func)
   for(i in 1:num.it){
      #Cria filhos usando Cauchy mutation
      n1 <- matrix(rcauchy(pop.size*dimension),nrow=pop.size) 
      xl <- x + std.dev * n1
      #Verifica os limites de cada gene dos filhos 
      idx <- which(xl < lb, arr.ind=TRUE)
      xl[idx] <- lb[idx[,2]]
      idx <- which(xl > ub, arr.ind=TRUE)()
      xl[idx] <- ub[idx[,2]]
      #Cria novos desvios    
      n1 <- rnorm(pop.size, mean = 0, sd = 1)
      n.norm <- rnorm(1, mean = 0, sd = 1)
      std.devl <- std.dev * exp(taul*n.norm + tau*n1)
      fitl <- apply(xl,1,func)
      #Faz a união entre pais e filhos
      tmp.pop <- rbind(xl,x)
      tmp.fit <- c(fitl,fit)
      tmp.std <- rbind(std.devl,std.dev)
      
      #Computa as vitórias de cada individuo
      win <- rep(NA, pop.size*2)
      for(j in 1:(2*pop.size)){ 
         idx <- sample(1:(2*pop.size),q)
         win[j] <- length(which(tmp.fit[j] < tmp.fit[idx]))
      }
      #Une as duas populações e ordena pelo numero de vitorias
      #Em seguida trunca-se pelo tamanho da população
      merge <- cbind(win,tmp.fit,tmp.pop,tmp.std)
      merge <- merge[order(merge[,1],decreasing = TRUE),]
      x <- merge[1:pop.size,3:(dimension+2)]
      fit <- merge[1:pop.size,2]
      std.dev <- merge[1:pop.size,(dimension+3):(dimension*2 + 2)]
   }
   return(list(pop=x,fit=fit,eta=std.dev))
}

O próximo código é o script de teste onde se pretende minimizar a função de Schwefel cujo mínimo é -12569,49.

dimension <- 30
max.it <- 9000
pop.size <- 100
q <- 10
func <- Schwefel
low <- -500
high <- 500
lb <- rep(low,dimension)
ub <- rep(high,dimension)
r <- EP(func,pop.size,dimension,lb,ub,max.it,q)
print(min(r$fit))

O código para função de Schwefel é:

Schwefel <- function(x){
   y <- sum(-x*sin(sqrt(abs(x))))
   return(y)
}

Espero que tenham gostado de mais este post e nos encontramos em breve com outro algoritmos evolutivo ou de enxame.

[1] L. J. Fogel, A. J. Owens, and M. J. Walsh, Artificial Intelligence Through Simulated Evolution. New York: Wiley, 1966.