模拟退火算法

求解函数

以求解以下这么一个函数为例子，实现代码为R语言

f(x)=x∗sin(10∗π∗x)+2x∈[−1,2] f(x)=x∗sin(10∗π∗x)+2x∈[−1,2]

其函数图像为：

求解流程与概念

初始解的产生

又称做状态产生函数，通常由两部分组成：
1）第一次产生候选解的分布函数。
2）第一次产生候选解不满足条件的情况下，再次产生解的分布函数。
下面采用的是标准差为2的正态分布、标准差为3的正态分布。其中要注意的是状态产生函数(领域函数)的出发点应该是尽可能保证产生的候选解遍布全部的解空间。

劣解接受概率（Metropolis准则）

又称做状态接受函数，这里是模拟退火法这个名字的来源，模仿固体退火原理，随着温度（迭代次数上升）的下降，能量逐渐稳定。即劣解的接受概率p p逐渐下降，其公式为：

p={1,e −E(x new )−E(x old )T  , E(x new )<E(x old )E(x new )≥E(x old )  p={1,E(xnew)<E(xold)e−E(xnew)−E(xold)T,E(xnew)≥E(xold)

从公式中可以看出这是求最小值时的分布，因为当新的解小于旧解时百分百接受。又可以看出当 E(x new )≥E(x old ) E(xnew)≥E(xold)时，p恒在【0，1】之间。其概率变化情况：

关于初始温度，初温越大，获得高质量解的机率越大，但花费较多的计算时间，原理可以从上面的图中或者结合下面温度随时间变化进行理解、解释，温度越高对于劣解的接受能力越大（搜索范围越大），下降速度越慢（求解速度下降）。
关于初始解如何设定一般有以下方法：
（1）按照某一概率分布抽样一组K个解，以K个解的目标值的方差为初温（下面使用的时正态分布）；。
（2）随机产生一组K个状态，确定两两状态间的最大目标差值，根据差值，利用一定的函数确定初温。
（3）利用经验公式。

温度更新

温度更新函数：

T 0 =T 0 ln(1+t)  T0=T0ln(1+t)

温度可以看出不管初始温度如何，都有一个先加热后散热的过程。即先提高搜索范围，再慢慢的减少搜索范围。至于为什么会在差不多时刻大家都把温度降低，我是这么进行理解的，他们只是在不断的接近，而不会重叠，同一时刻初始温度越大，其当前时刻温度也越大，都接近在同一刻收敛是因为他们总是在以一个倍数在进行减少。

代码

算法流程图：

代码：

# 定义目标函数
sloveFun <- function(x){
  (x*sin(10*pi * x) + 2)
}

limitX <- c(-1, 2)
#1. 初始化以及参数调整
stmp <- runif(100, limitX[1], limitX[2]) # 用于求初始温度设置
t0 <- var(sloveFun(stmp)) # 初始温度
s0 <- runif(1, limitX[1], limitX[2])  # 初始解
iters   <- 1000           # 迭代次数
ccnt    <- 200            # 终止条件，连续ccnt个新解都没有接受时终止算法
hisBest <- limitX[2]      # 默认取上界
ccntVc  <- NULL

#2.迭代
for (i in 1:iters) {
  #2.1 在s0附近产生新解，但又能包含定义内所有值
  s1 <- rnorm(1, mean = s0, sd = 2)
  while (s1 < limitX[1] || s1 > limitX[2]){ 
    s1 <- rnorm(1, mean = s0, sd=3)
  }
  #2.2 计算能量增量
  delta_t <- sloveFun(s0) - sloveFun(s1)

  #2.3 根据增量判断是否接受解
  if (delta_t > 0) { # 接受更好解
    s0 <- s1
    ccntVc <- c(ccntVc, 1)
  }else{ # 较差的解以一定概率接受
    p  <- exp(delta_t / t0)
    p1 <- sample(x = c(0, 1), size = 1, prob = c(1-p, p)) # 以一定概率生成 0,1：0 不接受劣解,1 接受劣解 
    if( p1 == 1 ) {
      s0 <- s1
      ccntVc <- c(ccntVc, 1)
    } else {
      ccntVc <- c(ccntVc, 0)
    }
  }

  #2.4 更新最优值
  hisBest <- ifelse(sloveFun(s1) < sloveFun(hisBest), s1, hisBest)
  hisBest <- ifelse(sloveFun(s0) < sloveFun(hisBest), s0, hisBest)

  #2.5 更新温度
  t0 <- t0/log(1 + i) # 温度随时间的推移而下降

  #2.6 检查终止条件
  ccntVcLen <- length(ccntVc)
  if (ccntVcLen > ccnt && # 模拟次数大于ccnt
      sum(ccntVc[(ccntVcLen - ccnt +1):ccntVcLen]) == 0) { # 连续ccnt次没有接受新解，ccntVc末尾200个全为0
    cat("连续", ccnt, "次没有接受新解，算法终止！","\n")
    cat("迭代次数为:", i, "\n")
    cat("最优解为:", hisBest, "\n")
    break()
  }
}
# 连续 200 次没有接受新解，算法终止！ 
# 迭代次数为: 257 
# 最优解为: 1.848313 

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• 52
• 53
• 54
• 55
• 56
• 57
• 58
• 59
• 60
• 61

在看了游皓麟一书的某些代码块后，越发感觉利用R语言学习算法的效率之高，代码简短，感觉就像可以运行的流程图。是一个学习的快速工具。

同问题粒子群算法求解：
http://blog.csdn.net/qq_27755195/article/details/62216762
遗传算法求解：
http://blog.csdn.net/qq_27755195/article/details/56597467