用Python实现人工鱼群算法

引言

前几天做了一个“华为杯”数学建模,用到了一个新算法,以此记录一下,也顺带解读一下用Python的实现方法。

PS:顺带一说,这一次的数模题目真的是难做。而且Bug贼多。


鱼群算法

简单介绍

关于鱼群算法的介绍巴拉巴拉的,网上多的不能再多了。我也就不赘述了。

主要说说,鱼群算法属于启发式算法的一种。我觉得自己对于启发式算法有一种,偏爱(其实就是普通算法数学上搞不懂啦XD)。同样属于启发式算法的还有:遗传算法、粒子群算法、禁忌算法、蚁群算法巴拉巴拉,哎呀,也是多了去了。

我感觉启发式算法的好,就好在不用从数学上去解决一个问题,而是从模型的角度,把问题抽象成各种能被算法识别的模型,然后剩下的就交给电脑和调参了。

鱼群算法就是四个四种行为:觅食行为、聚群行为、追尾行为和随机行为。关于整个算法的介绍,这里有一篇文章,我觉得还可以,介绍的很详细。详见这里

Python实现

今天的重点也不是讲解算法的原理,主要是想讲一讲算法的实现。一入Python深似海,所以还是选择用Python实现吧。其实用Matlab也可以,但是我还是钟爱Python。本篇使用的代码源自这个博客。后面还有根据今年“华为杯”数学建模B题目小改的算法。

人工智障鱼个体

既然是写程序,就要有写程序的思路。现在有一群鱼,每个鱼都有相同的行为和一些参数,所以首当其中,肯定是面向对象编程,一条鱼就是一个对象。所以先来定义一个鱼的类。

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
import numpy as np
import copy
import ObjFunction

class AFSIndividual:
"""人工智障鱼个体"""
def __init__(self, vardim, bound):
'''
vardim: 状态的维度
bound: 变量的上下限
'''
self.vardim = vardim
self.bound = bound

def generate(self):
'''
随机初始化本鱼
'''
len = self.vardim
rnd = np.random.random(size=len)
self.chrom = np.zeros(len)
self.velocity = np.random.random(size=len)
for i in xrange(0, len):
self.chrom[i] = self.bound[0, i] + \
(self.bound[1, i] - self.bound[0, i]) * rnd[i]
self.bestPosition = np.zeros(len)
self.bestFitness = 0.

def calculateFitness(self):
'''
计算本鱼的目标函数
'''
self.fitness = ObjFunction.GrieFunc(
self.vardim, self.chrom, self.bound)

其中,init函数接收了两个变量,第一个vardim是你的状态数有多少个,我感觉只能是一个数字,然后生成一个vardim维的行向量。

generate函数就是把这条鱼实例的n个状态随机初始化,还要保证在bound边界之内。for循环之中是启发式中很经典的初始化赋值过程。

calculateFitness函数就显得简单了很多,就是根据目标函数,计算这条鱼当前状态的适应值。

在这里我觉得,每条鱼都应该有四种行为,毕竟每个鱼都是依靠这四种行为更新自己的状态的,但是那个博主把四种行为写在了鱼群类里面。

人工智障鱼群

创建类

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
import numpy as np
from AFSIndividual import AFSIndividual
import random
import copy
import matplotlib.pyplot as plt


class ArtificialFishSwarm:

"""class for ArtificialFishSwarm"""

def __init__(self, sizepop, vardim, bound, MAXGEN, params):
'''
sizepop: 鱼群中有多少条人工智障鱼
vardim: 鱼有多少个状态
bound: 每个状态的上下限,2行n列,第一行是下限,第二行是上限
MAXGEN: 最大迭代次数
params: 需要的一些参数[visual, step, delta, trynum]依次为视觉范围,最大步长,拥挤度因子,尝试次数
'''
self.sizepop = sizepop
self.vardim = vardim
self.bound = bound
self.MAXGEN = MAXGEN
self.params = params
self.population = []
self.fitness = np.zeros((self.sizepop, 1))
self.trace = np.zeros((self.MAXGEN, 2))
self.lennorm = 6000

在这个里面他还有两个量,用来保存所有鱼的适应度fitness,还有一个不知道是啥的trace,可能和他的目标函数有关。还有一个lennorm也不知道是啥,我暂时认定为对视野的一个补偿值。

四种行为

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
class ArtificialFishSwarm:
# ···这里省去上面的那些了
def evaluation(self, x):
'''
evaluation the fitness of the individual
'''
x.calculateFitness()

def forage(self, x):
'''
artificial fish foraging behavior
'''
newInd = copy.deepcopy(x)
found = False
for i in xrange(0, self.params[3]):
indi = self.randSearch(x, self.params[0])
if indi.fitness > x.fitness:
newInd.chrom = x.chrom + np.random.random(self.vardim) * self.params[1] * self.lennorm * (
indi.chrom - x.chrom) / np.linalg.norm(indi.chrom - x.chrom)
newInd = indi
found = True
break
if not (found):
newInd = self.randSearch(x, self.params[1])
return newInd

def randSearch(self, x, searLen):
'''
artificial fish random search behavior
'''
ind = copy.deepcopy(x)
ind.chrom += np.random.uniform(-1, 1,
self.vardim) * searLen * self.lennorm
for j in xrange(0, self.vardim):
if ind.chrom[j] < self.bound[0, j]:
ind.chrom[j] = self.bound[0, j]
if ind.chrom[j] > self.bound[1, j]:
ind.chrom[j] = self.bound[1, j]
self.evaluation(ind)
return ind

def huddle(self, x):
'''
artificial fish huddling behavior
'''
newInd = copy.deepcopy(x)
dist = self.distance(x)
index = []
for i in xrange(1, self.sizepop):
if dist[i] > 0 and dist[i] < self.params[0] * self.lennorm:
index.append(i)
nf = len(index)
if nf > 0:
xc = np.zeros(self.vardim)
for i in xrange(0, nf):
xc += self.population[index[i]].chrom
xc = xc / nf
cind = AFSIndividual(self.vardim, self.bound)
cind.chrom = xc
cind.calculateFitness()
if (cind.fitness / nf) > (self.params[2] * x.fitness):
xnext = x.chrom + np.random.random(
self.vardim) * self.params[1] * self.lennorm * (xc - x.chrom) / np.linalg.norm(xc - x.chrom)
for j in xrange(0, self.vardim):
if xnext[j] < self.bound[0, j]:
xnext[j] = self.bound[0, j]
if xnext[j] > self.bound[1, j]:
xnext[j] = self.bound[1, j]
newInd.chrom = xnext
self.evaluation(newInd)
# print "hudding"
return newInd
else:
return self.forage(x)
else:
return self.forage(x)

def follow(self, x):
'''
artificial fish following behivior
'''
newInd = copy.deepcopy(x)
dist = self.distance(x)
index = []
for i in xrange(1, self.sizepop):
if dist[i] > 0 and dist[i] < self.params[0] * self.lennorm:
index.append(i)
nf = len(index)
if nf > 0:
best = -999999999.
bestIndex = 0
for i in xrange(0, nf):
if self.population[index[i]].fitness > best:
best = self.population[index[i]].fitness
bestIndex = index[i]
if (self.population[bestIndex].fitness / nf) > (self.params[2] * x.fitness):
xnext = x.chrom + np.random.random(
self.vardim) * self.params[1] * self.lennorm * (self.population[bestIndex].chrom - x.chrom) / np.linalg.norm(self.population[bestIndex].chrom - x.chrom)
for j in xrange(0, self.vardim):
if xnext[j] < self.bound[0, j]:
xnext[j] = self.bound[0, j]
if xnext[j] > self.bound[1, j]:
xnext[j] = self.bound[1, j]
newInd.chrom = xnext
self.evaluation(newInd)
# print "follow"
return newInd
else:
return self.forage(x)
else:
return self.forage(x)
def distance(self, x):
'''
return the distance array to a individual
'''
dist = np.zeros(self.sizepop)
for i in xrange(0, self.sizepop):
dist[i] = np.linalg.norm(x.chrom - self.population[i].chrom) / 6000
return dist

evaluation函数就是用来评价每个鱼,计算每个鱼的适应度。

forage就是觅食行为,其中的可以对应之前那个讲解原理的博客看,理解起来会容易些。randSearch就是随机行为。huddle就是群聚行为。follow就是追尾行为。distance就是计算每条鱼和自己的距离用的,在有些行为中,需要在视野范围内的鱼的状态,就用这个函数来计算视野范围内的鱼。

迭代过程和结果输出

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
class ArtificialFishSwarm:
# ···这里省去上面的那些了
def solve(self):
'''
evolution process for afs algorithm
'''
self.t = 0
self.initialize()
# evaluation the population
for i in xrange(0, self.sizepop):
self.evaluation(self.population[i])
self.fitness[i] = self.population[i].fitness
best = np.max(self.fitness)
bestIndex = np.argmax(self.fitness)
self.best = copy.deepcopy(self.population[bestIndex])
self.avefitness = np.mean(self.fitness)
self.trace[self.t, 0] = (1 - self.best.fitness) / self.best.fitness
self.trace[self.t, 1] = (1 - self.avefitness) / self.avefitness
print("Generation %d: optimal function value is: %f; average function value is %f" % (
self.t, self.trace[self.t, 0], self.trace[self.t, 1]))
while self.t < self.MAXGEN - 1:
self.t += 1
# newpop = []
for i in xrange(0, self.sizepop):
xi1 = self.huddle(self.population[i])
xi2 = self.follow(self.population[i])
if xi1.fitness > xi2.fitness:
self.population[i] = xi1
self.fitness[i] = xi1.fitness
else:
self.population[i] = xi2
self.fitness[i] = xi2.fitness
best = np.max(self.fitness)
bestIndex = np.argmax(self.fitness)
if best > self.best.fitness:
self.best = copy.deepcopy(self.population[bestIndex])
self.avefitness = np.mean(self.fitness)
self.trace[self.t, 0] = (1 - self.best.fitness) / self.best.fitness
self.trace[self.t, 1] = (1 - self.avefitness) / self.avefitness
print("Generation %d: optimal function value is: %f; average function value is %f" % (
self.t, self.trace[self.t, 0], self.trace[self.t, 1]))

print("Optimal function value is: %f; " % self.trace[self.t, 0])
print "Optimal solution is:"
print self.best.chrom
self.printResult()

def printResult(self):
'''
plot the result of afs algorithm
'''
x = np.arange(0, self.MAXGEN)
y1 = self.trace[:, 0]
y2 = self.trace[:, 1]
plt.plot(x, y1, 'r', label='optimal value')
plt.plot(x, y2, 'g', label='average value')
plt.xlabel("Iteration")
plt.ylabel("function value")
plt.title("Artificial Fish Swarm algorithm for function optimization")
plt.legend()
plt.show()

solve就是整个迭代过程,其中就是一次次按照四种行为更新每条鱼的状态,然后求每条鱼的适应度,把最好的记录下来,然后继续下一次迭代,直到最后满足条件,就停止。

printResult这个函数把最后的他自己计算的一个值trace的整个迭代过程给画出来了(话说这个值到底干啥用的我还是没搞明白)。

主函数

1
2
3
4
if __name__ == "__main__":
bound = np.tile([[-600], [600]], 25)
afs = AFS(60, 25, bound, 500, [0.001, 0.0001, 0.618, 40])
afs.solve()

目标函数

这个目标函数是他遗传算法中用的目标函数,在这里也可以用用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import math


def GrieFunc(vardim, x, bound):
"""
Griewangk function
"""
s1 = 0.
s2 = 1.
for i in range(1, vardim + 1):
s1 = s1 + x[i - 1] ** 2
s2 = s2 * math.cos(x[i - 1] / math.sqrt(i))
y = (1. / 4000.) * s1 - s2 + 1
y = 1. / (1. + y)
return y

其实就是一个y=f(x)函数,把x给进去计算yy就代表了适应度,f(x)就是目标函数,也叫适应度函数。


如何取状态变量

这个问题,是用启发式算法的关键。和大多数启发式算法一样,都需要将实际问题转化成一个个状态变量,这样就可以方便交给电脑去迭代,寻找一组最优的状态变量。

对于上面的那个数学公式问题,可以就把x作为状态变量,有多少个x,就有多少个状态变量。

对于今年“华为杯”的B题第二问,我知道你们不会去查的,所以我还是简述一下问题。

要在十二座城市之间搭建电网,每两个城市之间的电网价值是两个城市的人口和能传输的总容量的一个函数。一共能搭建十六根电线,每两个城市之间只能有一根,问能达到的最大电网价值是多少。

在解这样的问题中,要设计一种状态变量,变量数要尽可能的少,还可以根据这些变量计算出各种需要的东西。

这个题目中,我们需要计算每根线的价值,网络中有没有没有被连接到的城市和信息孤岛。

我们的解法是,为每一根线设计了两个变量,分别表示这根线两端的城市。这样,我们用了三十二个变量把十六根线表示出来了。而用着三十二个变量就可以求解处总价值和网络是否全联通。

目标函数就是总价值+是否全联通*惩罚系数。关于这个可以多看看启发式算法,一般叫做代价函数,包括你需要的值,和不满足约束产生的惩罚量。


离散化问题

用这个算法的时候,我就发现了离散化的问题。这个算法里面的状态是连续的,不是离散的。在计算适应度函数的时候就会出问题。

还是刚刚那个题,题目中的状态变量,我们设定为每根线两端代表的城市,然后用1~12表示十二个城市,所以变量的上下限是112,但是计算目标函数的时候,只能是整数。

在这里有几种做法,我用的一种是,在迭代过程中让他保持连续,但是在计算代价函数的时候,将每个状态按照一定规则取整。这个一定的规则可以是四舍五入向下取整向上取整等,很多方法都行。但是一旦决定,就要重新确定一下状态的取值范围,要保证状态能取得每一个值在取整的过程中都有相同的概率。

从上面可以看出,如果是向下取整,12就只有在状态变量刚好为12的时候取得到,所以应该将上限改为13,这样每一个值的概率都一样。


后记

不知道是因为接触过一些启发式算法的原因还是怎么,看鱼群算法的时候觉得还比较好理解。以前也接触过粒子群算法、遗传算法和Hopfield神经网络。粒子群和鱼群有点像。

不过这个代码的原博客,那个博主写了好多得启发式算法,而且都是Python实现的,有机会我再拿几个来,给大家讲解讲解。