R 语言手写遗传算法（待完善）

R
启发式优化算法
mcga
Author

Rui

Published

September 8, 2022

我太懒了，文字部分回头补上

遗传算法思想

步骤

实例

目标函数

\[ f(x,y) = \frac{\sin(\sqrt{x^2+y^2})}{\sqrt{x^2+y^2}}+\exp{\{\frac{\cos(2\pi x)+\cos(2\pi y)}{2}}\}-2.71289 \]

Code 
a <- b <- seq(-2, 2, length=101)

f <- function(a, b) {
  term1 <- sin(sqrt(a^2 + b^2)) / sqrt(a^2 + b^2)
  term2 <- exp(0.5*(cos(2*pi*a) + cos(2*pi*b)))
  term3 <- -2.71289
  
  term1 + term2 + term3
}

z <- outer(a, b, f)

library(plot3D)
persp3D(a, b, z, contour = TRUE, zlim= c(-4, 2), image = FALSE)

随机化初始值

假设该函数中每个变量的定义域都为 [-10, 10] 若选取 100 个染色体，那么每个染色体具有 2 个基因且与上述 2 个变量一一对应。思路是创建一个 \(100\times 2\) 的随机矩阵，每行代表一个染色体，每一列代表一种基因。使用 runif 函数产生 [-10, 10] 的随机数向量，再放入 \(100\times 2\) 的矩阵中。

Code 
n <- 100 # 染色体个数
D <- 2 # 变量个数（维度），也是基因数
bond <- 10 # 基因界限，限制每个基因的取值范围[-10, 10]

set.seed(1)
data <- matrix(runif(n*D, -bond, bond), nrow = n, ncol = D, byrow = TRUE)
head(data)
##           [,1]      [,2]
## [1,] -4.689827 -2.557522
## [2,]  1.457067  8.164156
## [3,] -5.966361  7.967794
## [4,]  8.893505  3.215956
## [5,]  2.582281 -8.764275
## [6,] -5.880509 -6.468865

适应度函数计算

稍微改动一下目标函数：

\[ f(x,y) = \frac{\sin(\sqrt{x^2+y^2})+\varepsilon}{\sqrt{x^2+y^2}+\varepsilon}+\exp{\{\frac{\cos(2\pi x)+\cos(2\pi y)}{2}}\}-2.71289 \]

Note

\(\varepsilon\) 是一个较小的正数，这么做是为了防止分母为零，同时保证当 x 和 y 趋于 0 时，极限为 1。

适应度函数在目标函数的基础上稍作改动，目的是保证适应度函数为正。对于不同的优化问题（最大化、最小化）、目标函数的具体形式，适应度函数可以做不同改变，也可以取目标函数的倒数、绝对值、自然对数等。

\[ F(c) = \frac{\sin(\sqrt{x^2+y^2})+\varepsilon}{\sqrt{x^2+y^2}+\varepsilon}+\exp{\{\frac{\cos(2\pi x)+\cos(2\pi y)}{2}}\}+10 \]

Code 
################################### 目标函数 ###################################

fn1 <- function (x) {
  # 将这个复杂的函数拆解成
  term1 <- (sin(sqrt(sum(x^2))) + 0.001) / (sqrt(sum(x^2)) + 0.001) 
  term2 <- exp(0.5*sum(cos((2*pi*x))))
  
  # 将各部分合并起来得到最终的目标函数
  term1 + term2 
}
# 这是向量化的目标函数（x是向量）


################################## 适应度函数 ##################################

fitness_function <- function (x) {
  outcome = fn1(x) + 10 # 适应度函数应该为正
  
  return(outcome)
}

选择算子

使用轮盘赌方法，先计算每个染色体的适应度，再按照适应度计算概率，被选中的概率与适应度函数成正比。即适应度越强的染色体越容易被选中，而适应度越差的染色体越不容易被选中，以此来模拟”适者生存”。

\[ P_i=\frac{F_{i}(c)}{Total}=\frac{F_{i}(c)}{\sum_{i=1}^{n}F_{i}(c)} \]

Code 
################################### 选择算子 ###################################

select_function <- function (data) {
  
  # 轮盘赌选择染色体（总数不变）
  # 计算概率
  f <- apply(data, MARGIN = 1, FUN = fitness_function) # 得到所有染色体的适应度（向量）
  P <- f/sum(f) # 概率
  
  # 依据概率对染色体抽样
  temp <- sample(1:n, size = n, replace = TRUE, prob = P)
  afterselect <- data[temp, ]
  
  return (afterselect)
}
Note

一些细节的实现因人而异，比如有人提出在每次选择过程中应该强制将适应度函数最高的染色体 100% 保留下来，而将其他染色体进行轮盘赌选择。这里我选择最原始的轮盘赌方法。

交叉算子

在生物学中，交叉是指生殖的一个术语。两条染色体被随机选择并通过数学运算进行匹配。在本例中使用单点交叉。

单点交叉意味着两个亲本的基因被一个交叉线交换。

选择用于交叉的染色体数量是由交叉率（Pc）控制的，介于 0-1 之间。例如确定 Pc = 0.25，这意味着大约 25% 的染色体将成为亲本进行交叉。

为了确定交叉线的位置，需要生成一个 1 到 D-1 之间的随机数，D-1 为交叉的点位数。

Code 
################################## 交叉算子 ####################################

cross_function <- function (afterselect, Pc) {
  
  # 随机选择亲本
  cross_rand <- sample(c(FALSE, TRUE), size = n, replace = TRUE, prob = c(1-Pc, Pc))
  cross_data <- afterselect[cross_rand, ]
  uncross_data <- afterselect[!cross_rand, ]
  
  # 用于保存子代染色体的空矩阵
  r <- nrow(cross_data)
  c <- ncol(cross_data)
  offspring <- matrix(vector(length = r*c), nrow = r, ncol = c) 
  
  # 交叉
  for (i in nrow(cross_data)) {
    
    r1 <- sample(1:n, size = 1)
    r2 <- sample(1:n, size = 1)
      
      while (r1 == r2) {
        r1 <- sample(1:n, size = 1)
        r2 <- sample(1:n, size = 1)
      } # 确保两个随机数不一样
    
    # 随机交叉点位
    cross_position <- sample(x = 1:(D-1), size = 1) 
    
    # 染色体交叉
    offspring[i, 1:cross_position] <- data[r1, 1:cross_position]
    offspring[i, (cross_position+1):D] <- data[r2, (cross_position+1):D]
    
  }
  
  # 交叉后的子代与未被选做亲本的染色体组成新的染色体矩阵
  aftercross <- rbind(offspring, uncross_data)
}

注：

  • 大部分教程只是阐述了交叉的大致过程，并没有说明细节。比如：

    • 被选做亲本的每一条染色体，是与其他被选做亲本的染色体任意配对并交叉的吗？如果是这样，那么一共有 \(C_k^2=\frac{k(k-1)}{2}\) 种组合（k 为被选做亲本的染色体数）。

    • 其次，每种染色体组合交叉后产生 1 个子代还是 2 个子代？如果产生 2 个子代，那么交叉后一共产生 \(k(k-1)\) 个子代。这意味着什么呢？要知道 k 还是随着染色体总个数 n 的增加而增加的，如此一来 n 随着迭代次数呈指数型增长，大约 20 次迭代后开始难以计算了！

  • 所以应该要保证迭代过程中染色体总数保持不变（或者相对稳定）。在参考 多起点遗传算法 后，认为被选做亲本的染色体 cross_data 从所有染色体 data 中随机抽选配对、交叉，并且只产生 1 个子代（父亲的前半段与母亲的后半段染色体拼接，而不是反过来），如此一共得到 k 个子代，再加上未被选做亲本的染色体（n-k 个），最终可以保证迭代过程中染色体总数不变。

变异算子

又称”基因突变”，本例使用随机突变。

Code 
##################################### 变异算子 #################################

mutate_function <- function(aftercross, mutation_rate) {
  
  r <- sample(
    x = c(0, 1), 
    size = dim(aftercross)[1] * dim(aftercross)[2], 
    replace = TRUE,
    prob = c(1-mutation_rate, mutation_rate)
  )
  
  mutation_position <- matrix(r, ncol = D, nrow = n, byrow = TRUE) # 基因突变点位
  
  mutation_num <- sum(r) # 基因突变个数
  mutation_rand <- runif(n = mutation_num, min = -bond, max = bond) # 基因突变随机数
  r[r==1] <- mutation_rand
  mutation <- matrix(r, ncol = D, nrow = n, byrow = TRUE) 
  
  aftercross[mutation_position == 1] <- 0
  
  aftermutation <- aftercross + mutation # 发生基因突变
  
  return (aftermutation)
}

停止条件

  1. 人为设置迭代次数

  2. 人为设置运行时间

  3. 若两次迭代的适应度函数增量的绝对值小于事先给定的阈值则停止

本文选择第 1 和第 3 种方法，并这一部分将嵌入主体代码中，参见下一节。

完整实现

以上各代码块为了方便展示，每个函数都只分别考虑算法的某一部分。而下面的实现方式则是按照算法流程将各个函数嵌套起来，形式上略有不同但是输出结果不变。

准备函数

Code 
################################### 目标函数 ###################################
fn1 <- function (x) {
  # 将这个复杂的函数拆解成3部分
  term1 <- (sin(sqrt(sum(x^2))) + 0.001) / (sqrt(sum(x^2)) + 0.001) 
  # 分母加上0.001是为了防止分母为0，分子同样加上0.001是为了保证x趋于零向量时，极限依然为1
  term2 <- exp(0.5*sum(cos((2*pi*x))))
  # 将3个部分合并起来得到最终的函数
  term1 + term2 
}
# 这是向量化的目标函数


################################## 适应度函数 ##################################
fitness_function <- function (x) {
  outcome = fn1(x) + 10 # 适应度函数应该为正
  
  return(outcome)
} 


################################### 选择算子 ###################################
select_function <- function (data) {
  
  # 轮盘赌选择染色体（总数不变）
  # 计算概率
  f <- apply(data, MARGIN = 1, FUN = fitness_function) # 得到所有染色体的适应度（向量）
  P <- f/sum(f) # 概率
  
  # 依据概率对染色体抽样
  temp <- sample(1:n, size = n, replace = TRUE, prob = P)
  afterselect <- data[temp, ]
  
  return (afterselect)
}  


################################### 交叉算子####################################
cross_function <- function (data, Pc) {
  
  # 得到经选择后的染色体
  afterselect <- select_function(data)
  
  # 随机选择亲本
  cross_rand <- sample(c(FALSE, TRUE), size = n, replace = TRUE, prob = c(1-Pc, Pc))
  cross_data <- afterselect[cross_rand, ]
  uncross_data <- afterselect[!cross_rand, ]
  
  r <- nrow(cross_data)
  c <- ncol(cross_data)
  offspring <- matrix(vector(length = r*c), nrow = r, ncol = c)
  
  for (i in nrow(cross_data)) {
    
    r1 <- sample(1:n, size = 1)
    r2 <- sample(1:n, size = 1)
      
      while (r1 == r2) {
        r1 <- sample(1:n, size = 1)
        r2 <- sample(1:n, size = 1)
      } # 确保两个随机数不一样
    
    # 随机交叉点位
    cross_position <- sample(x = 1:(D-1), size = 1) 
    
    # 染色体交叉
    offspring[i, 1:cross_position] <- data[r1, 1:cross_position]
    offspring[i, (cross_position+1):D] <- data[r2, (cross_position+1):D]
    
  }
  
  aftercross <- rbind(offspring, uncross_data)
  
  return (aftercross)
}


################################## 变异算子 ####################################
mutate_function <- function(data, mutation_rate) {
  
  # 得到经交叉后的染色体
  aftercross <- cross_function(data, Pc)
  
  r <- sample(
    x = c(0, 1), 
    size = dim(aftercross)[1] * dim(aftercross)[2], 
    replace = TRUE,
    prob = c(1-mutation_rate, mutation_rate)
  )
  
  mutation_position <- matrix(r, ncol = D, nrow = n, byrow = TRUE) # 基因突变点位
  mutation_num <- sum(r) # 基因突变个数
  mutation_rand <- runif(n = mutation_num, min = -bond, max = bond) # 基因突变随机数
  r[r==1] <- mutation_rand
  mutation <- matrix(r, ncol = D, nrow = n, byrow = TRUE) 
  
  aftercross[mutation_position == 1] <- 0
  
  aftermutation <- aftercross + mutation # 发生基因突变
  
  return (aftermutation)
}

主体部分

Code 

ga_alg <- function (n, D, bond, Pc, mutation_rate, maxiter, theta) {
  
  # 随机初始化染色体矩阵
  set.seed(1)
  data <- matrix(runif(n*D, -bond, bond), nrow = n, ncol = D, byrow = TRUE)

  # 监视迭代过程中的适应度
  fit <- c() # 用来保存迭代中的适应度函数,用以监视算法迭代过程
  fitness_value <- apply(data, MARGIN = 1, FUN = fitness_function)
  fit[1] <- mean(fitness_value) # 保存初始值的适应度

  # 监视迭代过程中的适应度增量
  delta_fit <- c() # 用来保存迭代中的适应度增量，用以监视算法迭代过程

  # 监视迭代过程中的染色体总数
  chrom_num <- c()
  chrom_num[1] <- nrow(data)

  # 监视迭代过程中的染色体种类
  chrom_variety <- c()
  chrom_variety[1] <- nrow(unique(data))

  # 开始迭代
  for (i in 1:maxiter) {
  
  
    # 计算迭代次数
    iteration <- i
  
    # 遗传算法
    data <- mutate_function(data, mutation_rate)
  
    # 计算适应度
    fitness_value <- apply(data, MARGIN = 1, FUN = fitness_function)
    fit[i+1] <- mean(fitness_value)
  
    # 计算适应度增量
    delta_fit[i] <- abs(fit[i+1] - fit[i])
  
    # 计算经过选择、交叉、变异后染色体总数
    chrom_num[i+1] <- nrow(data)
  
    # 计算经过选择、交叉、变异后染色体种类
    chrom_variety[i+1] <- nrow(unique(data))
  
    # 算法停止条件
    if (!is.null(delta_fit[i])) {
      if (delta_fit[i] < theta) {break}
    }  
  }
  
  # 计算适应度函数最大值
  temp <- apply(data, MARGIN = 1, FUN = fitness_function)
  optim <- max(temp)
  
  # 计算最优解
  solution <- data[which.max(temp), ]
  
  # 将结果汇总
  out <- list(
    "data" = data, 
    "optim" = optim,
    "solution" = solution,
    "iteration" = iteration, 
    "fitness" = fit,
    "deltafitness" = delta_fit, 
    "chrom_num" = chrom_num, 
    "chrom_variety" = chrom_variety
  )
  
  return (out)
}

在计算适应度函数部分 fit[*] <- mean(fitness_value) ，也就是说在每轮迭代中选择所有染色体适应度的均值（mean()）来代表本次迭代的适应度水平。如果使用 min()，max() 或者 median() 效果会如何呢？

Note

主体部分也抽象为一个函数

准备参数并运行

Code 
n <- 100 # 染色体个数
D <- 2 # 基因个数
bond <- 10 # 基因界限，限制每个基因的取值范围[-10, 10]
Pc <- 0.2 # 亲本选择率
mutation_rate <- 0.01 # 变异率
maxiter <- 500 # 最大迭代次数
theta <- 1e-7 # 适应度函数增量阈值，本例中各变量都为连续型实数
Note

为了方便展示，本文以 2 维函数为例。实际上变量维数可以拓展到多维。

Code 
ga <- ga_alg(n, D, bond, Pc, mutation_rate, maxiter, theta)

最大值：

Code 
ga$optim
## [1] 13.71828

最优解：

Code 
ga$solution
## [1] 0 0

可视化

Code 
################################ 可视化迭代过程 ################################
par(mfrow = c(2, 1))
plot(
  x = 0:ga$iteration, y = ga$fitness, 
  xlab = "iteration", ylab = "fitness",
  type = "l", main = "适应度"
)

plot(
  x = 1:ga$iteration, y = ga$deltafitness, 
  xlab = "iteration", ylab = "delta fitness",
  type = "l", main = "适应度增量"
)

Code 


############################### 可视化染色体 ###################################
par(mfrow = c(2, 1))
plot(
  x = c(0:ga$iteration), y = ga$chrom_num, 
  xlab = "iteration", ylab = "num of chroms",
  type = "l", main = "染色体总数"
)

plot(
  x = c(0:ga$iteration), y = ga$chrom_variety, 
  xlab = "iteration", ylab = "variety of chroms",
  type = "l", main = "染色体种类"
)

Code 


############################### 可视化解的范围 #################################
par(mfrow = c(2, 1))

plot(density(ga$data[ , 1]), xlab = "", ylab = "", main = "")
rug(jitter(ga$data[ , 1]))

plot(density(ga$data[ , 2]), xlab = "", ylab = "", main = "")
rug(jitter(ga$data[ , 2]))

Note

当然，也可以将可视化部分也抽象为函数方便重复调用

最后一幅图可以看出两个变量的取值都集中分布在 0 附近。

使用 mcga 实现遗传算法

mcga 包是一个遗传算法快速的工具包，主要解决实值优化的问题，用于求解多维函数的最小值。它使用的变量值表示基因序列，而不是字节码，因此不需要编解码的处理。mcga 实现了遗传算法的交配和突变的操作，并且可以进行大范围和高精度的搜索空间的计算，算法的主要缺点是使用了 256 位的一元字母表。mcga 包中用于实现遗传算法的函数是 mcga()，具体信息可以使用 ?mcga 查看。

安装

安装 mcga 包之前需要先安装 GA 包。GA 包也是一个强大的遗传算法包，它能够编码、遗传算法以及变种、可视化、统计描述等一系列完整功能，之后一节再来介绍。

install.packages("GA")
install.packages("mcga")

参数说明

输入：

  • popsize：个体数量，即染色体数目

  • chsize：基因数量，限参数的数量

  • crossprob：交叉概率，默认为 1.0

  • mutateprob：突变概率，默认为 0.01

  • elitism：精英数量，直接复制到下一代的染色体数目，默认为 1

  • minval：随机生成初始种群的下边界值

  • maxval：随机生成初始种群的上边界值

  • maxiter：繁殖次数，即循环次数，默认为 10

  • evalFunc：适应度函数，用于给个体进行评价

输出：

  • population 经过迭代后排序过的种群结果

  • costs 结果迭代后每一个个体（染色体）对应的适应度函数值

Code 
library(mcga)

fn2 <- function (x) {
  # 将这个复杂的函数拆解成3部分
  term1 <- (sin(sqrt(sum(x^2))) + 0.001) / (sqrt(sum(x^2)) + 0.001) 
  # 分母加上0.001是为了防止分母为0，分子同样加上0.001是为了保证x趋于零向量时，极限依然为1
  term2 <- exp(0.5*sum(cos((2*pi*x))))
  
  
  # 将3个部分合并起来得到最终的函数
  - term1 - term2 + 2.71289 
}
# 因为mcga求解最小值，所以在原来的目标函数前加负号


m <- mcga(
  popsize = 100, 
  chsize = 2, 
  minval = -10,              
  maxval = 10, 
  maxiter = 500, 
  crossprob = 0.2, 
  mutateprob = 0.01, 
  evalFunc = fn2
)
Code 
# 最优个体
cat("Best Chromosome :\n")
## Best Chromosome :
print(m$population[1, ])
## [1]  1.957547e-19 -1.357197e-14
Code 
# 最优适应度函数值
cat("Cost: ", m$costs[1], "\n")
## Cost:  -1.005392

参考

[1] https://mp.weixin.qq.com/s/RQblfcxXsorXChE732Pn5w

[2] https://mp.weixin.qq.com/s/Bmymma5f-znKjpd86OmtYg

[3] https://www.bilibili.com/video/BV1cL411t7LF/?spm_id_from=333.788.recommend_more_video.29&vd_source=f9f3226ec290b0512a03337b4cb90df3

[4] http://blog.fens.me/algorithm-ga-r/

https://www.bilibili.com/video/BV1WK4y1U7af?vd_source=f9f3226ec290b0512a03337b4cb90df3