f1 <- function(x) { return( x[1]^2 + x[2]^2 ) } f2 <- function(x) { return( (x[1]-1)^2 + (x[2]-1)^2 ) } bif <- function(x) { return(c(f1(x), f2(x))) } rnd_sol <- function(f, mu) { pop <- array(dim=c(mu, 2)) fpop <- array(dim=c(mu, 2)) for(i in 1:mu) { x <- runif(2, min = -5, max = 5) pop[i,] <- x fpop[i,] <- f(x) } return(list(pop=pop, fitness_pop=fpop)) } moead_d <- function(f, mu, sigma, n_eval_max) { # directions lambda <- array(dim=c(mu, 2)) for(i in 0:(mu-1)) { angle = i*pi/((mu-1) * 2) lambda[i+1,] <- c(cos(angle), sin(angle)) } # neighbor directions neighbor_dir <- array(dim=c(mu, 3)) neighbor_dir[1,] <- c(1, 2, 3) for(i in 2:(mu-1)) { neighbor_dir[i,] <- c(i-1, i, i+1) } neighbor_dir[mu,] <- c(mu - 2, mu - 1, mu) # initial population res <- rnd_sol(f, mu) pop <- res$pop fitness_pop <- res$fitness_pop n_eval <- mu while (n_eval < n_eval_max) { # pick at random one direction i <- sample(1:mu, size = 1) # select at random one solution from the neighboring direction j <- sample(neighbor_dir[i,], size = 1) x <- pop[j,] # random mutation y <- x + rnorm(n = 2, mean = 0, sd = sigma) # evaluation fy <- f(y) n_eval <- n_eval + 1 # local replacement for all neighboring directions for(k in neighbor_dir[j,]) { gk <- lambda[k,1] * fitness_pop[k,1] + lambda[k,2] * fitness_pop[k,2] gy <- lambda[k,1] * fy[1] + lambda[k,2] * fy[2] if (gy < gk) { # better, so replace pop[k, ] <- y fitness_pop[k,] <- fy } } } return(fitness_pop) } main <- function() { res0 <- rnd_sol(bif, 20) res <- moead_d(bif, mu = 20, sigma = 0.1, n_eval_max = 1000) res1 <- moead_d(bif, mu = 20, sigma = 0.1, n_eval_max = 100000) plot(res[,1], res[,2], type ="p", ylim=c(0, 1), xlim=c(0, 1), pch = "+") points(res1[,1], res1[,2], pch = "+", col = "green") }