Implementando a Evolução Diferencial em R

A evolução diferencial foi proposta por Storn e Price em 1997 [1]. Seu funcionamento é muito similar a dos AGs, porém com os operadores invertidos. Enquanto nos AGs executa-se o cruzamento e depois a mutação, na evolução diferencial (ED) o processo é inverso, ou seja, primeiro realiza-se a mutação para em seguida executar o cruzamento. O algoritmo é apresentado a seguir:

Assim como nas outras meta-heurísticas que apresentei, na ED uma população inicial aleatória de tamanho n e dimensão d é criada com o código a seguir, que basicamente é o mesmo dos algoritmos genéticos (AGs).

init.population <- function(func,lb,ub,pop.size,dimension){
   pop <- matrix(runif(pop.size*dimension),nrow=pop.size)
   fitness <- rep(NA,pop.size)
   pop <- lb + pop*(ub-lb)
   fitness <- apply(pop,1,func)
   return(list(pop = pop, fit = fitness))
}

Em seguida entra-se no laço até atingir o critério de parada. Dentro desse laço, três indivíduos são selecionados aleatoriamente. Em seguida a diferença entre dois dos indivíduos (idx1 e idx2) é multiplicado pelo fator F sendo adicionado a um terceiro indivíduo (idx3). A ordem não precisa ser necessariamente essa, ou seja, poderia ser a diferença entre 3 e 2 sendo adiciona a 1, por exemplo, como será feito no código. Essa operação é a criação do vetor de diferenças, que também é chamada de mutação. Além disso, essa operação dá nome à estratégia utilizada no processo de otimização. No caso apresentado, a estratégia é chamada DE/Rand/1, pois a operação de seleção é aleatória e gera um vetor v para cada seleção de 3 indivíduos.

Outra estratégia que pode ser utilizada é a DE/Best/1, na qual o elemento de indice 3 sempre será o melhor indivíduo da população, enquanto idx1 e idx2, continuam aleatórios. Aqui incentivo o leitor a buscar por outras estratégias que podem ser usadas na ED.

O código a seguir ilustra em R como as duas estratégias mencionadas podem ser implementadas. Note que na estratégia DE/Best/1, o melhor individuo é replicado em um vetor (best), facilitando assim a operação vetorial sem a necessidade de uma instrução for.

if(strategy == 1){ #DE/Rand/1
   idx <- matrix(sample(1:pop.size,3*pop.size,replace = TRUE),nrow =   
          pop.size)
}
else { #DE/Best/1
   idx <- matrix(sample(1:pop.size,2*pop.size, replace = TRUE), 
          nrow = pop.size)
   best <- rep(which.min(fit),pop.size)
   idx <- cbind(best,idx)
}
v <- pop[idx[,1],] + F * (pop[idx[,2],] - pop[idx[,3],])

Em seguida para cada dimensão do vetor um valor aleatório é sorteado, se esse valor for menos do que o crossover rate (taxa de cruzamento) o gene do novo indivíduo virá do vetor de diferenças, caso contrário virá do elemento da população que estiver sendo mutado. Particularmente em R, costumo criar a matriz inteira de CRs, isto é, para todos os genes da população, e em seguida formar a nova população usando indexação lógica, como mostrado no código abaixo:

r <- matrix(runif(pop.size*dimension), nrow = pop.size)
idx <- r < pc
new.pop[idx] <- v[idx]
new.pop[!idx] <- pop[!idx]
new.fit <- apply(new.pop,1,func)

E para finalizar é necessário identificar quem dos novos indivíduos são melhores do que a população atual, o que também é feito através da indexação lógica:

idx <- new.fit < fit
fit[idx] <- new.fit[idx]
pop[idx,] <- new.pop[idx,]

E pronto, esta finalizado o algoritmo. O código a seguir é da função completa de ED. Adiciono a ela somente um pedaço de código no qual verifico se os vetores de diferenças estão dentro dos limites inferior e superior de cada gene.

DE <- function(func,lb,ub,pop.size,dimension,iterations,F, pc = 
               0.6,strategy = 1){
  tmp.pop <- init.population(func,lb,ub,pop.size,dimension)
  pop <- tmp.pop$pop
  fit <- tmp.pop$fit
  new.pop <- matrix(rep(NA, pop.size * dimension), nrow = pop.size)
  new.fit <- c()
  for(i in 1:iterations){
    if(strategy == 1){
       idx <- matrix(sample(1:pop.size, 3 * pop.size, replace = 
              TRUE), nrow = pop.size)
    }
    else {
      idx <- matrix(sample(1:pop.size, 2 * pop.size, replace = 
             TRUE), nrow = pop.size)
      best <- rep(which.min(fit),pop.size)
      idx <- cbind(best,idx)
    }
    
    #Verifies lower and upper bounds
    v <- pop[idx[,1],] + F * (pop[idx[,2],] - pop[idx[,3],])
    idx <- which(v < lb,arr.ind=TRUE)
    v[idx] <- lb[idx[,2]]
    idx <- which(v > ub,arr.ind=TRUE)
    v[idx] <- ub[idx[,2]]
    
    #Crossover
    r <- matrix(runif(pop.size*dimension), nrow = pop.size)
    idx <- r < pc
    new.pop[idx] <- v[idx]
    new.pop[!idx] <- pop[!idx]
    
    #Update population 
    new.fit <- apply(new.pop,1,func)
    idx <- new.fit < fit
    fit[idx] <- new.fit[idx]
    pop[idx,] <- new.pop[idx,]
  }
  return(list(pop = pop, fit = fit))
}

O código de teste esta logo abaixo o qual pode selecionar tanto a função de Schwefel quanto a de Rosenbrock, bastando comentar a linha correta e estabelecer os novos limites inferiores e superiores, respectivamente:

source('DE.R')
source('init.population.R')
#source('./../Numerical-Benchmarks/Rosenbrock.R')
source('./../Numerical-Benchmarks/Schwefel.R')
dim <- 30
it <- 2000
pop.size <- 50
lb <- rep(-500,dim)
ub <- rep(500,dim)
func <- Schwefel
execs <- 20
F <- 0.3
pc <- 0.6
#DE/Rand/1 = 1
#DE/Best/1 = 2
strategy <- 2
result <-  vector("list",execs)
best <- rep(NA,execs)
for(i in 1:execs){
   cat("exec = ",i,"\n")
   result[[i]] <- DE(func,lb,ub,pop.size,dim,it,F, pc=0.6,strategy)
   best[i] <- min(result[[i]]$fit)
}
print(min(best))

A função de Schwefel foi apresentada no artigo “Programação Evolutiva em R” enquanto que a Rosenbrock foi mostrada no artigo “Implementando AGs em R”.

[1] STORN, R. and PRICE, K. Differential evolution — A simple and efficient heuristic for global optimization over continuous spaces. J. Global Optimization, v. 11, pp. 341–359, 1997.