参考：《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]之间的随机数
生成随机数： 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]]
若要生成制定范围的随机数，表达式： r = l b + ( u b − l b ) ∗ r a n d o m ( ) r = lb+(ub-lb)*random() r=lb+(ub−lb)∗random()
# 若要生成制定范围的随机数，表达式：$r = lb+(ub-lb)*random()$
RandValue2 = 0+(4-0)*np.random.random([1,5])
print("生成随机数：", RandValue2) # 生成1*5个[0,4]之间的随机数
生成随机数： [[0.73625711 1.71444157 0.90705669 0.87469253 2.35687254]]
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,:])
(100, 5)
[ 1.0648139 4.97723729 -1.21941472 1.24288407 -0.57282862]
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)
fitness: 5
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)
X: [[ 1 -1 1 -1]
[ 1 -1 1 -1]]
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()
最优适应度值： [2.28139543e-171]
最优解[x1,x2]： [-2.22236523e-86 4.22788920e-86]
扫一扫
1026
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
