三种智能优化算法解决十个无约束单目标优化问题

摘要

本文主要讲述的是解决十个无约束单目标优化问题的过程和结果,即求这些函数在给定区间内的最小值。
首先说明了用遗传算法解决的思路和结果,但是其中有两个函数无法得到理想结果,于是使用粒子群算法重新计算,但是无法解决的问题反而变成了四个,于是又用差分演化算法尝试解决,结果发现遗传算法无法解决的问题差分演化也没办法解决,但是另外八个问题的结果精度上比遗传算法高了很多。最后经过反复尝试各种修改参数和代码调试,在把迭代次数从5000上调到50000之后,发现无法解决的两个函数才表现出一定的稳定收敛。

关键词

无约束单目标优化、连续优化、最小化、遗传算法、粒子群算法、差分演化算法

一 目标问题概述

Wenyin Gong. 2019. Test functions for single-objective optimization. 1, 1, Article 1 (January 2019), 2 pages.

1. Rastrigin Function

2. Ackley Function

3. Sphere Function

4. Rosenbrock Function

5. Beale Function

6. Schaffer Function

7. Styblinski-Tang Function

8. Bukin Function

9. Himmelblaus Function

10. Cross in tray Function

程序设计思路

程序实现的思路主要是以下几个部分:

  1. 输入部分(决定解决的问题)
  2. 确定变量维度以及每个维度的上下界(根据选择的问题)
  3. 初始化参数
  4. 算法迭代
  5. 输出结果

为提高可读性,本文的程序运行结果数据以及代码统一放在文末,并会有对比分析。如果有兴趣可以前往参考。

遗传算法

用遗传算法(Genetic Algorithm)解决,程序算法迭代部分如下:

  1. 选择。采用君主方案选择,即将结果排序后,选出最好的那一个个体,然后奇数位个体全部替换成这个最好的,偶数位个体全部与最好的个体进行交叉。
  2. 变异。随机变异。
  3. 新的种群和上一代种群合并,并排序,留下前一半的个体作为新的种群进行迭代。
  4. 记录排序后种群的第一个个体的函数值

迭代运行5000次,使用遗传算法能够确认八个函数的迭代结果会比较稳定,误差也在可接受范围内,但是Rosenbrock函数和Bukin函数无法稳定到可以推测出的最小值,考虑到可能是这种算法不适合解决这两个问题,所以我们尝试其他的优化算法。

粒子群算法

由于Rosenbrock函数和Bukin函数使用遗传算法没法收敛,我们把初始化参数和算法迭代换成粒子群算法,其迭代的步骤如下:

  1. 遍历每个个体
  2. 对每个个体,如果其目标函数值要比历史最佳要好,更新这个个体历史最佳函数值,以及将这个个体的参数更新
  3. 再对每个个体,根据其个体历史最佳函数值,决定是否更新全局最佳函数值以及更改全局最佳个体参数
  4. 更新速度和位置
  5. 边界处理
  6. 记录这一代最佳函数值

然后对十个函数进行了重新测试,迭代次数为5000次,测试结果有点出人意料,Rosenbrock函数多次出现了直接等于0.0的结果(即直接出了最小值,而不是接近),但是不稳定,多次测试还是会出现其他结果。Bukin函数则还是不收敛。而之前没有问题的Rastrigin函数和Styblinski-Tang函数却也无法收敛了。于是我们再尝试差分演化算法。

差分演化算法

我们将初始化参数和算法迭代改成差分演化算法,其迭代步骤如下:

  1. 变异算子随着迭代进行深入,每次都相应改变
  2. 随机变异
  3. 随机交叉
  4. 边界处理
  5. 选择
  6. 记录最好的函数值

差分演化算法迭代次5000次,遗传算法中没办法得到结果的Rosenbrock Function和Bukin Function还是没办法得到结果,其他的问题能很好的解决,并且和粒子群算法一样,精度极高,这个在文末数据分析里会对比讨论。

最终解决

到了这个时候,我大概能猜到一点问题的所在了,通过差分算法,即使这两个函数的结果还是有部分波动,但是它们确实在一个贴近结果(0.0)的附近波动,所以方法应该是没错了,但是精度还是有问题。

而我之前由于在个人电脑上运行,考虑到出结果的时间,同时也通过观察迭代过程图发现很早就能收敛到较低的范围,而越往后提升基本不大了,没想过用一个特别大的迭代数,因为要等的时间确实太久了。

在我把迭代次数从5000次改成50000次之后,精度提高了一个数量级,这说明了对于优化算法迭代次数将很大程度决定最终的效果。但是通过分析数据,Rosenbrock函数和Bukin函数相对于其他函数,更难获得最小值,这和函数本身的性质以及算法都有关系,由于篇幅原因就不进行深入了。

最后,很感谢《智能优化技术》这门课的老师,龚文引教授对我的启发以及帮助,学习完这门课收获很多!同时由于我本人水平有限,文中难免出现疏漏之处,还望谅解!

数据分析

关于遗传算法、粒子群算法、差分演化算法在迭代次数为5000次时解决十个问题分别得到的结果对比:

对比可以发现,

  1. 遗传算法和差分演化算法在解决Rosenbrock函数和Bukin函数的时候效果都不佳,但是差分演化算法还是表现的要比遗传算法要好。
  2. 粒子群算法无法解决Rastrigin函数和Bukin函数。但是在解决Rosenbrock函数、Styblinski-Tang函数以及Himmelblau’s函数时出现了一个奇怪的现象,它会有时候收敛到最小值(另外两个算法很难解决的Rosenbrock函数在这里直接给出了精确值),但是有时候会收敛到另一个固定的值附近。
  3. 在正常解决的函数里,粒子群算法和差分演化算法要比遗传算法精确很多很多,粒子群算法和差分演化算法在解决的函数里表现的精度很接近。

完整代码

首先是模块funclib.py,这个模块用来实现十个优化问题,同时用只一个函数调用,提供输出菜单中显示的序号就行。

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
62
# funclib.py

from math import *

def rastrigin_func(X):
A = 10.0
sum = 0.0
for xi in X:
sum += (xi ** 2 - A * cos(2 * pi * xi))
return sum + A * len(X)

def ackley_func(X):
return -20 * exp(-0.2 * ((0.5 * (X[0] ** 2 + X[1] ** 2)) ** 0.5)) - exp(0.5 * (cos(2 * pi * X[0]) + cos(2 * pi * X[1]))) + e + 20

def sphere_func(X):
return sum([xi ** 2 for xi in X])

def rosenbrock_func(X):
sum = 0.0
for i in range(0, len(X) - 1):
sum += (100 * (X[i + 1] - X[i] ** 2) ** 2 + (1 - X[i]) ** 2)
return sum

def beale_func(X):
return (1.5 - X[0] + X[0] * X[1]) ** 2 + (2.25 - X[0] + X[0] * (X[1] ** 2)) ** 2 + (2.625 - X[0] + X[0] * (X[1] ** 3)) ** 2

def schaffer_func(X):
return 0.5 + (cos(sin(abs(X[0] ** 2 - X[1] ** 2))) ** 2 - 0.5) / (1 + 0.001 * (X[0] ** 2 + X[1] ** 2)) ** 2

def styblinski_tang_func(X):
return sum([(xi ** 4 - 16 * (xi ** 2) + 5 * xi) / 2 for xi in X])

def bukin_func(X):
return 100 * (abs(X[1] - 0.01 * (X[0] ** 2)) ** 0.5) + 0.01 * abs(X[0] + 10)

def himmelblaus_func(X):
return (X[0] ** 2 + X[1] - 11) ** 2 + (X[0] + X[1] ** 2 - 7) ** 2

def corssintray_func(X):
return -0.0001 * ((abs(sin(X[0]) * sin(X[1]) * exp(abs(100 - (X[0] ** 2 + X[1] ** 2) ** 0.5 / pi))) + 1) ** 0.1)

def obj_func(func_i, X):
if func_i == '0':
return rastrigin_func(X)
elif func_i == '1':
return ackley_func(X)
elif func_i == '2':
return sphere_func(X)
elif func_i == '3':
return rosenbrock_func(X)
elif func_i == '4':
return beale_func(X)
elif func_i == '5':
return schaffer_func(X)
elif func_i == '6':
return styblinski_tang_func(X)
elif func_i == '7':
return bukin_func(X)
elif func_i == '8':
return himmelblaus_func(X)
else:
return corssintray_func(X)

每个算法文件中需要导入的其他模版。

1
2
3
4
5
import numpy as np
# 通过obj_func来调用函数,提供对应函数序号和自变量数组作为参数
from funclib import obj_func
# 绘图库
import matplotlib.pyplot as plt

程序的输入部分,由用户选择要测试哪个函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 注意从0开始
message = "Choose the index of function to calculate:\n \
0. Rastrigin function\n \
1. Ackley function\n \
2. Sphere function\n \
3. Rosenbrock function\n \
4. Beale function\n \
5. Schaffer function\n \
6. Styblinski_tang function\n \
7. Bukin function\n \
8. Himmelblaus function\n \
9. Cross_in_tray function\n"

func_to_cal = input(message)
while func_to_cal > '9' and func_to_cal < '0':
func_to_cal = input(message)

根据前面的用户选的函数,设定自变量的维度和上下限。

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
#解的维度
pop_dim = 2
if func_to_cal in "0236":
pop_dim = 10

#解空间的上下界
#如果是2维则有x,y之分
#x_top, x_bottom, y_top, y_bottom
#如果是n维则只有一个上下界
#xi_top, xi_bottom

if func_to_cal in "0236":
if func_to_cal == '0':
xi_top = 5.12
xi_bottom = -5.12
elif func_to_cal == '2':
xi_top = 100
xi_bottom = -100
elif func_to_cal == '3':
xi_top = 30
xi_bottom = -30
else:
xi_top = 5
xi_bottom = -5
else:
if func_to_cal == '1':
x_top = y_top = 5
x_bottom = y_bottom = -5
elif func_to_cal == '4':
x_top = y_top = 4.5
x_bottom = y_bottom = -4.5
elif func_to_cal == '5':
x_top = y_top = 100
x_bottom = y_bottom = -100
elif func_to_cal == '7':
x_top = -5
x_bottom = -15
y_top = 3
y_bottom = -3
elif func_to_cal == '8':
x_top = y_top = 5
x_bottom = y_bottom = -5
else:
x_top = y_top = 10
x_bottom = y_bottom = -10

遗传算法的参数初始化。

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
# 种群规模
pop_num = 100
# 随机初始化种群
pop = np.random.rand(pop_num, pop_dim)
# 选择、交叉、变异后的新种群
new_pop = np.random.rand(pop_num, pop_dim)

# 确定维度
if pop_dim == 2 :
pop[:,0] *= (x_top - x_bottom)
pop[:,0] += np.repeat(np.array([x_bottom]), pop_num, axis=0)
pop[:,1] *= (y_top - y_bottom)
pop[:,1] += np.repeat(np.array([y_bottom]), pop_num, axis=0)
else:
for i in range(0, pop_dim):
pop[:,i] *= (xi_top - xi_bottom)
pop[:,i] += np.repeat(np.array([xi_bottom]), pop_num, axis=0)

# 交叉概率
porb_cro = 0.8
# 变异概率
prob_mut = 0.1
# 迭代次数
iterations = 5000
# 历代最优
trace = np.zeros(iterations, dtype=float)

# 计算初代的结果
results = np.zeros(pop_num, dtype=float)
for i in range(0, pop_num):
results[i] = obj_func(func_to_cal, pop[i])

# 由于是求最小值,所以升序排列结果,并记录原来的索引位置
sorted_indices = np.argsort(results)
results.sort()

# 根据目标函数值对种群进行排序,按照函数值升序排列
sorted_pop = np.zeros((pop_num, pop_dim), dtype=float)
for i in range(0, pop_num):
sorted_pop[i] = pop[sorted_indices[i]].copy()

pop = sorted_pop.copy()

遗传算法的迭代过程。

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
62
63
#开始迭代遗传
for iter in range(0, iterations):
#君主方案进行选择
emper = pop[0].copy()
#每次交叉点的个数
noPoint = round(pop_dim * porb_cro)
#交叉基因的位置
poPoint = np.random.randint(0, pop_dim, size=(int(pop_num / 2), noPoint))
#交叉操作
new_pop = pop.copy()
for i in range(0, int(pop_num / 2)):
new_pop[2 * i] = emper.copy()
for j in range(0, noPoint):
new_pop[2 * i][poPoint[i][j]] = new_pop[2 * i + 1][poPoint[i][j]]
new_pop[2 * i + 1][poPoint[i][j]] = emper[poPoint[i][j]]

#变异操作
for i in range(0, pop_num):
for j in range(0, pop_dim):
choice = np.random.rand()
if choice < prob_mut:
if pop_dim == 2:
if j == 0:
new_pop[i][j] = np.random.rand() * (x_top - x_bottom) + x_bottom
else:
new_pop[i][j] = np.random.rand() * (y_top - y_bottom) + y_bottom
else:
new_pop[i][j] = np.random.rand() * (xi_top - xi_bottom) + xi_bottom

#子种群按照函数值升序排列
new_results = np.zeros(pop_num, dtype=float)
for i in range(0, pop_num):
new_results[i] = obj_func(func_to_cal, new_pop[i])

new_sorted_indices = np.argsort(new_results)
new_results.sort()

#获得新的按照函数值升序的种群
new_sorted_pop = np.zeros((pop_num, pop_dim), dtype=float)
for i in range(0, pop_num):
new_sorted_pop[i] = new_pop[new_sorted_indices[i]].copy()
new_pop = new_sorted_pop.copy()

#合并两代种群并排序
#函数值合并
total_results = np.hstack((results, new_results))
#种群合并
total_pop = np.vstack((pop, new_pop))
#获得合并后的函数值升序索引列表
total_sorted_indives = np.argsort(total_results)
total_results.sort()
#将合并后的种群按升序排列
sorted_total_pop = np.zeros((pop_num * 2, pop_dim), dtype=float)
for i in range(0, 2 * pop_num):
sorted_total_pop[i] = total_pop[total_sorted_indives[i]].copy()
total_pop = sorted_total_pop.copy()

# 获得前一半函数值和种群
results = total_results[0:pop_num].copy()
pop = total_pop[0:pop_num].copy()

# 记录本代最佳的函数值
trace[iter] = results[0]

粒子群算法的参数初始化。

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
# 种群数量
pop_num = 100
# 迭代次数
iterations = 50000
# 两个速度更新参数
c1 = 1.5
c2 = 1.5
# 速度变化权重
w = 0.8
# 根据维度设置速度上限和下限
if pop_dim == 2:
v1_max = x_top / 2.0
v1_min = x_bottom / 2.0
v2_max = y_top / 2.0
v2_min = y_bottom / 2.0
else:
v_max = xi_top / 2.0
v_min = xi_bottom / 2.0
# 随机初始化种群和其对应速度
pop = np.random.rand(pop_num, pop_dim)
pop_v = np.random.rand(pop_num, pop_dim)

if pop_dim == 2 :
pop[:,0] *= (x_top - x_bottom)
pop[:,0] += np.repeat(np.array([x_bottom]), pop_num, axis=0)
pop[:,1] *= (y_top - y_bottom)
pop[:,1] += np.repeat(np.array([y_bottom]), pop_num, axis=0)

pop_v[:,0] *= (v1_max - v1_min)
pop_v[:,0] += np.repeat(np.array([v1_min]), pop_num, axis=0)
pop_v[:,1] *= (v2_max - v2_min)
pop_v[:,1] += np.repeat(np.array([v2_min]), pop_num, axis=0)
else:
for i in range(0, pop_dim):
pop[:,i] *= (xi_top - xi_bottom)
pop[:,i] += np.repeat(np.array([xi_bottom]), pop_num, axis=0)

pop_v[:,i] *= (v_max - v_min)
pop_v[:,i] += np.repeat(np.array([v_min]), pop_num, axis=0)

# 初始化个体最优位置和最优值
pop_best = pop.copy()
pop_best_obj = np.zeros(pop_num)
for i in range(0, pop_num):
pop_best_obj[i] = obj_func(func_to_cal, pop[i])

# 初始化全局最优位置和最优值
glo_indi_best = pop[0].copy()
glo_indi_best_obj = np.Infinity

# 根据初代随机值第一次更新全局最佳位置和函数值
for i in range(0, pop_num):
if pop_best_obj[i] < glo_indi_best_obj:
glo_indi_best = pop_best[i].copy()
glo_indi_best_obj = pop_best_obj[i]
# 用来记录历代最佳
trace = np.zeros(iterations)

粒子群算法的迭代过程。

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
# 开始迭代
for iter in range(0, iterations):
for indi in range(0, pop_num):
# 更新个体最优位置和最优值
if obj_func(func_to_cal, pop[indi]) < obj_func(func_to_cal, pop_best[indi]):
pop_best[indi] = pop[indi].copy()
pop_best_obj[indi] = obj_func(func_to_cal, pop[indi])

# 更新全局最优位置和最优值
if pop_best_obj[indi] < glo_indi_best_obj:
glo_indi_best = pop_best[indi].copy()
glo_indi_best_obj = pop_best_obj[indi]

# 更新位置和速度
pop_v[indi] = w * pop_v[indi]
pop_v[indi] += c1 * np.random.rand() * (pop_best[indi] - pop[indi])
pop_v[indi] += c2 * np.random.rand() * (glo_indi_best - pop[indi])

pop[indi] = pop[indi] + pop_v[indi]

# 边界处理
if pop_dim == 10:
for i in range(0, pop_dim):
if pop_v[indi][i] > v_max or pop_v[indi][i] < v_min:
pop_v[indi][i] = np.random.rand() * (v_max - v_min) + v_min
if pop[indi][i] > xi_top or pop[indi][i] < xi_bottom:
pop[indi][i] = np.random.rand() * (xi_top - xi_bottom) + xi_bottom
else:
if pop_v[indi][0] > v1_max or pop_v[indi][0] < v1_min:
pop_v[indi][0] = np.random.rand() * (v1_max - v1_min) + v1_min
if pop[indi][0] > x_top or pop[indi][0] < x_bottom:
pop[indi][0] = np.random.rand() * (x_top - x_bottom) + x_bottom

if pop_v[indi][1] > v2_max or pop_v[indi][1] < v2_min:
pop_v[indi][1] = np.random.rand() * (v2_max - v2_min) + v2_min
if pop[indi][1] > y_top or pop[indi][1] < y_bottom:
pop[indi][1] = np.random.rand() * (y_top - y_bottom) + y_bottom

# 记录本代最佳函数值
trace[iter] = glo_indi_best_obj

差分演化算法的参数初始化。

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
# 种群数量
pop_num = 100
# 迭代次数
iterations = 50000
# 初始变异算子
mut_coe_0 = 0.4
# 交叉算子
cro_coe = 0.2

# 原始种群
raw_pop = np.random.rand(pop_num, pop_dim)
# 变异种群
mut_pop = np.zeros((pop_num, pop_dim))
# 选择种群
cho_pop = np.zeros((pop_num, pop_dim))

# 原始种群赋值
if pop_dim == 2 :
raw_pop[:,0] *= (x_top - x_bottom)
raw_pop[:,0] += np.repeat(np.array([x_bottom]), pop_num, axis=0)
raw_pop[:,1] *= (y_top - y_bottom)
raw_pop[:,1] += np.repeat(np.array([y_bottom]), pop_num, axis=0)
else:
for i in range(0, pop_dim):
raw_pop[:,i] *= (xi_top - xi_bottom)
raw_pop[:,i] += np.repeat(np.array([xi_bottom]), pop_num, axis=0)

# 群体的每个个体的初始目标函数值
objs = np.zeros(pop_num)
for i in range(0, pop_num):
objs[i] = obj_func(func_to_cal, raw_pop[i])

# 每一代的最佳函数值(包含初始的一代)
trace = np.zeros((iterations + 1))
# 初始一代最佳函数值
trace[0] = objs.min()

差分演化算法的迭代过程。

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
62
63
# 迭代
for iter in range(0, iterations):
# 变异
# 自适应变异算子
lamda = exp(1.0 - float(iterations) / (float(iterations) - iter))
mut_coe = mut_coe_0 * (2 ** lamda)

for indi in range(0, pop_num):
r1 = np.random.randint(0, pop_num)
while r1 == indi:
# 确保不是当前个体
r1 = np.random.randint(0, pop_num)
r2 = np.random.randint(0, pop_num)
while r2 == indi or r2 == r1:
# 确保不是当前个体且不与r1相等
r2 = np.random.randint(0, pop_num)
r3 = np.random.randint(0, pop_num)
while r3 == indi or r3 == r1 or r3 == r2:
# 确保不是当前个体且不与r1,r2相等
r3 = np.random.randint(0, pop_num)
mut_pop[indi] = raw_pop[r1] + mut_coe * (raw_pop[r2] - raw_pop[r3])

# 交叉

for n in range(0, pop_num):
# 交叉点
# r = np.random.randint(0, pop_dim)
for i in range(0, pop_dim):
cro_r = np.random.rand()
if cro_r < cro_coe:
cho_pop[n][i] = mut_pop[n][i].copy()
else:
cho_pop[n][i] = raw_pop[n][i].copy()

# 边界条件处理
for i in range(0, pop_num):
if pop_dim == 2:
if cho_pop[i][0] > x_top or cho_pop[i][0] < x_bottom:
cho_pop[i][0] = np.random.rand() * (x_top - x_bottom) + x_bottom
if cho_pop[i][1] > y_top or cho_pop[i][1] < y_bottom:
cho_pop[i][1] = np.random.rand() * (y_top - y_bottom) + y_bottom
else:
for j in range(0, pop_dim):
if cho_pop[i][j] > xi_top or cho_pop[i][j] < xi_bottom:
cho_pop[i][j] = np.random.rand() * (xi_top - xi_bottom) + xi_bottom

# 选择
# 选择种群的函数值
cur_objs = np.random.rand(pop_num)
for i in range(0, pop_num):
cur_objs[i] = obj_func(func_to_cal, cho_pop[i])

# 如果选择种群的函数值好一些就赋值给原始种群
for i in range(0, pop_num):
if cur_objs[i] < objs[i]:
raw_pop[i] = cho_pop[i].copy()

# 再重新计算原始种群的函数值
for i in range(0, pop_num):
objs[i] = obj_func(func_to_cal, raw_pop[i])

# 记录这一代最好的函数值
trace[iter + 1] = objs.min()

结果的输出。

1
2
3
4
5
print(trace[-1])
# 这里是绘图过程
x = range(0, iterations + 1)
plt.scatter(x, trace)
plt.show()

关于引用

如果有需要转载文章或者使用源代码,只需要在文首加上:
原作者:Yaron Ho
作者Blog:https://hoyyy.me

十分感谢。