黏菌算法的Python实现_黏菌优化算法python

参考：《Python智能优化算法：从原理到代码实现与应用》

1 黏菌算法的基本原理

核心思想：迭代、替换最优解、精英策略

2 准备工作

2.1 种群初始化

import numpy as np
RandValue = np.random.random()
print("生成随机数：", RandValue) # 生成一个[0,1]之间的随机数

# 生成多个随机数
RandValue1 = np.random.random([3,4])
print("生成随机数：", RandValue1) # 生成3*4个[0,1]之间的随机数
1
2
3
4
5
6
7

生成随机数： 0.77163342604771
生成随机数： [[0.23492356 0.30544602 0.76447861 0.25265406]
[0.44425183 0.61387201 0.97107584 0.86323364]
[0.9268611  0.92894365 0.26571302 0.79201158]]
1
2
3
4

若要生成制定范围的随机数，表达式：$r=lb+(ub-lb)\ast random()$

# 若要生成制定范围的随机数，表达式：$r = lb+(ub-lb)*random()$
RandValue2 = 0+(4-0)*np.random.random([1,5])
print("生成随机数：", RandValue2) # 生成1*5个[0,4]之间的随机数
1
2
3

生成随机数： [[0.73625711 1.71444157 0.90705669 0.87469253 2.35687254]]
1

2.2 编写初始化函数

def initialization(pop,ub,lb,dim):
    ''' 初始化函数
        pop:种群数量,dim:每个个体的维度,ub:每个个体的变量上界，维度为[dim],
        lb:每个个体的变量下界，维度为[dim],X:输出的种群，维度为[pop,dim]
    '''
    X = np.zeros([pop,dim]) # 声明空间
    X = lb+(ub-lb)*np.random.random([pop,dim])  # 方法1
    
    # for i in range(pop):
    #     for j in range(dim):
    #         X[i,j] = lb[j] + (ub[j]-lb[j])*np.random.random() # 方法2

    print(X.shape)
    return X

pop = 100
dim = 5
ub = np.array([5,5,5,5,5])
lb = np.array([-5,-5,-5,-5,-5])
X = initialization(pop,ub,lb,dim)
print(X[1,:])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

(100, 5)
[ 1.0648139   4.97723729 -1.21941472  1.24288407 -0.57282862]
1
2

2.2 适应度函数（优化问题的目标函数）

def fun(x):
    '''适应度函数'''
    '''X为输入的一个个体,维度为[1,dim],fitness 为输出的适应度值'''
    fitness = np.sum(x**2)
    return fitness
x = np.array([1,2])
fitness = fun(x)
print("fitness:",fitness)
1
2
3
4
5
6
7
8

fitness: 5
1

2.3 边界检查和约束函数

防止变量超出可行变量范围
v a l = { u b , v a l > u b l b , v a l < l b v a l , 其它 val=\left\{

$$\begin{array}{l} ub, val>ub \\ lb, val<lb \\ val, 其它 \end{array}$$ \right. val=⎩ ⎨ ⎧​ub,val>ublb,val<lbval,其它​

def BorderCheck(X,ub,lb,pop,dim):
    '''边界检查函数'''
    '''pop:种群数量,dim:每个个体的维度,ub:每个个体的变量上界，维度为[dim],
        lb:每个个体的变量下界，维度为[dim],X:输入的种群，维度为[pop,dim]'''
    # # 方法1
    # lb = lb*np.ones([pop,dim])
    # ub = ub*np.ones([pop,dim])
    # X = np.min([np.max([X,lb],axis=0),ub],axis=0) #axis用来指定在哪个维度上比较

    # 方法2
    for i in range(pop):
        for j in range(dim):
            if X[i,j] > ub[j]:
                X[i,j] = ub[j]
            if X[i,j] <lb[j]:
                X[i,j] = lb[j]
    return X

x= np.array([(1,-2,3,-4),
            (1,-2,3,-4)])
ub = np.array([1,1,1,1])
lb = np.array([-1,-1,-1,-1])
dim = 4
pop = 2
X = BorderCheck(x, ub,lb,pop,dim)
print("X:",X)
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

X: [[ 1 -1  1 -1]
 [ 1 -1  1 -1]]
1
2

3 黏菌算法Python实现

import numpy as np
import copy 
from matplotlib import pyplot as plt
'''适应度函数'''
def fun(X):
        O=X[0]**2 + X[1]**2
        return O

'''黏菌优化算法求解x1^2 + x2^2的最小值'''

'''主函数 '''
#设置参数
pop = 50 #种群数量
MaxIter = 100 #最大迭代次数
dim = 2 #维度
lb = -10*np.ones(dim) #下边界
ub = 10*np.ones(dim)#上边界
#适应度函数选择
fobj = fun

def initialization(pop,ub,lb,dim):
    ''' 黏菌种群初始化函数'''
    '''
    pop:为种群数量
    dim:每个个体的维度
    ub:每个维度的变量上边界，维度为[dim,1]
    lb:为每个维度的变量下边界，维度为[dim,1]
    X:为输出的种群，维度[pop,dim]
    '''
    X = np.zeros([pop,dim]) #声明空间
    for i in range(pop):
        for j in range(dim):
            X[i,j]=(ub[j]-lb[j])*np.random.random()+lb[j] #生成[lb,ub]之间的随机数
    
    return X
     
def BorderCheck(X,ub,lb,pop,dim):
    '''边界检查函数'''
    '''
    dim:为每个个体数据的维度大小
    X:为输入数据，维度为[pop,dim]
    ub:为个体数据上边界，维度为[dim,1]
    lb:为个体数据下边界，维度为[dim,1]
    pop:为种群数量
    '''
    for i in range(pop):
        for j in range(dim):
            if X[i,j]>ub[j]:
                X[i,j] = ub[j]
            elif X[i,j]<lb[j]:
                X[i,j] = lb[j]
    return X


def CaculateFitness(X,fun):
    '''计算种群的所有个体的适应度值'''
    pop = X.shape[0]
    fitness = np.zeros([pop, 1])
    for i in range(pop):
        fitness[i] = fun(X[i, :])
    return fitness


def SortFitness(Fit):
    '''适应度排序'''
    '''
    输入为适应度值
    输出为排序后的适应度值，和索引
    '''
    fitness = np.sort(Fit, axis=0)
    index = np.argsort(Fit, axis=0)
    return fitness,index

def SortPosition(X,index):
    '''根据适应度对位置进行排序'''
    Xnew = np.zeros(X.shape)
    for i in range(X.shape[0]):
        Xnew[i,:] = X[index[i],:]
    return Xnew

def SMA(pop,dim,lb,ub,MaxIter,fun):
    '''黏菌优化算法'''
    '''
    输入：
    pop:为种群数量
    dim:每个个体的维度
    ub:为个体上边界信息，维度为[1,dim]
    lb:为个体下边界信息，维度为[1,dim]
    fun:为适应度函数接口
    MaxIter:为最大迭代次数
    输出：
    GbestScore:最优解对应的适应度值
    GbestPositon:最优解
    Curve:迭代曲线
    '''
    z = 0.03 #位置更新参数
    X = initialization(pop,ub,lb,dim) #初始化种群
    fitness = CaculateFitness(X,fun) #计算适应度值
    fitness,sortIndex = SortFitness(fitness) #对适应度值排序
    X = SortPosition(X,sortIndex) #种群排序
    GbestScore = copy.copy(fitness[0])
    GbestPositon = copy.copy(X[0,:])
    Curve = np.zeros([MaxIter,1])
    W = np.zeros([pop,dim]) #权重W矩阵
    for t in range(MaxIter):
        worstFitness = fitness[-1]
        bestFitness = fitness[0]
        S=bestFitness-worstFitness+ 10E-8 #当前最优适应度于最差适应度的差值，10E-8为极小值，避免分母为0；
        for i in range(pop):
            if i<pop/2: #适应度排前一半的W计算
                W[i,:]= 1+np.random.random([1,dim])*np.log10((bestFitness-fitness[i])/(S)+1)
            else:#适应度排后一半的W计算
                W[i,:]= 1-np.random.random([1,dim])*np.log10((bestFitness-fitness[i])/(S)+1)
        #惯性因子a,b
        tt = -(t/MaxIter)+1
        if tt!=-1 and tt!=1:
            a = np.math.atanh(tt)
        else:
            a = 1
        b = 1-t/MaxIter
        #位置更新
        for i in range(pop):
            if np.random.random()<z:
                X[i,:] = (ub.T-lb.T)*np.random.random([1,dim])+lb.T #公式（1.4）第一个式子
            else:
                p = np.tanh(abs(fitness[i]-GbestScore))
                vb = 2*a*np.random.random([1,dim])-a
                vc = 2*b*np.random.random([1,dim])-b
                for j in range(dim):
                    r = np.random.random()
                    A = np.random.randint(pop)
                    B = np.random.randint(pop)
                    if r<p:
                        X[i,j] = GbestPositon[j] + vb[0,j]*(W[i,j]*X[A,j]-X[B,j]) #公式（1.4）第二个式子
                    else:
                        X[i,j] = vc[0,j]*X[i,j]         #公式(1.4)第三个式子
        
        X = BorderCheck(X,ub,lb,pop,dim) #边界检测       
        fitness = CaculateFitness(X,fun) #计算适应度值
        fitness,sortIndex = SortFitness(fitness) #对适应度值排序
        X = SortPosition(X,sortIndex) #种群排序
        if(fitness[0]<=GbestScore): #更新全局最优
            GbestScore = copy.copy(fitness[0])
            GbestPositon = copy.copy(X[0,:])
        Curve[t] = GbestScore
    
    return GbestScore,GbestPositon,Curve

GbestScore,GbestPositon,Curve = SMA(pop,dim,lb,ub,MaxIter,fobj)
   
print('最优适应度值：',GbestScore)
print('最优解[x1,x2]：',GbestPositon)

#绘制适应度曲线
plt.figure(1)
plt.plot(Curve,'r-',linewidth=2)
plt.xlabel('Iteration',fontsize='medium')
plt.ylabel("Fitness",fontsize='medium')
plt.grid()
plt.title('SMA',fontsize='large')
plt.show()
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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

最优适应度值： [2.28139543e-171]
最优解[x1,x2]： [-2.22236523e-86  4.22788920e-86]
1
2