Python优化算法16——鲸鱼优化算法(WOA)

科研里面优化算法都用的多，尤其是各种动物园里面的智能仿生优化算法，但是目前都是MATLAB的代码多，python几乎没有什么包，这次把优化算法系列的代码都从底层手写开始。

需要看以前的优化算法文章可以参考：Python优化算法_阡之尘埃的博客-CSDN博客

鲸鱼优化算法（Whale Optimization Algorithm, WOA）是一种基于自然启发的优化算法，由S. Mirjalili及其同事在2016年提出。该算法模拟了座头鲸的捕食行为，特别是座头鲸通过“气泡网捕食”策略来捕捉猎物的过程。

基本概念

WOA的核心思想是通过模拟座头鲸的捕食过程来进行搜索和优化。鲸鱼在捕猎时会围绕猎物游动并产生气泡网，迫使猎物聚集。这一行为被用来设计搜索策略，使算法能够有效地找到全局最优解。

算法流程

初始化：

随机生成一组初始解，这些解被称为鲸鱼个体。

适应度评估：

计算每个鲸鱼个体的适应度值，根据优化问题的目标函数来评估解的质量。

气泡网攻击：

模拟鲸鱼围绕猎物进行的螺旋形游动，将这一过程用于局部搜索。
在这个阶段，鲸鱼个体通过螺旋上升的方式逐步逼近猎物（当前最优解）。

包围猎物：

模拟鲸鱼围绕猎物收缩包围圈，通过缩小搜索空间实现全局搜索。
在这个阶段，鲸鱼个体按照某种概率在当前解附近进行局部开发。

随机搜索：

在某些情况下，鲸鱼个体会随机搜索以增强探索能力，避免陷入局部最优。

更新最优解：

根据适应度信息更新全局最佳解，以指导鲸鱼个体的下一步搜索。

迭代：

重复上述过程，直到满足停止条件，如达到最大迭代次数或找到满意的解。

优势与应用

鲸鱼优化算法具有以下优势：

平衡探索与开发：通过气泡网攻击、包围猎物和随机搜索三种行为，WOA在全局搜索和局部开发之间实现了良好的平衡。
简单高效：算法结构简单，易于实现，适用于各种连续和离散优化问题。

由于这些优势，WOA在工程设计优化、机器学习参数优化、图像处理等多个领域得到了广泛应用。像其他群体智能算法一样，WOA的性能可能依赖于参数设置和具体问题特征，因此在实际应用中需要根据问题性质进行调整和优化。

原理不多介绍了，直接看代码就好。

代码实现

导入包

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 
import seaborn as sns
import warnings
import copyplt.rcParams ['font.sans-serif'] ='SimHei'               #显示中文
plt.rcParams ['axes.unicode_minus']=False               #显示负号
warnings.filterwarnings('ignore')
plt.rcParams['font.family'] = 'DejaVu Sans'

只给代码不给使用案例就都是钓鱼的。我这里给出代码，也要给使用案例，先采用一些简单的优化算法常用的测试函数。由于都优化算法需要测试函数，我们先都定义好常见的23个函数：

'''F1函数'''
def F1(X):Results=np.sum(X**2)return Results'''F2函数'''
def F2(X):Results=np.sum(np.abs(X))+np.prod(np.abs(X))return Results'''F3函数'''
def F3(X):dim=X.shape[0]Results=0for i in range(dim):Results=Results+np.sum(X[0:i+1])**2return Results'''F4函数'''
def F4(X):Results=np.max(np.abs(X))return Results'''F5函数'''
def F5(X):dim=X.shape[0]Results=np.sum(100*(X[1:dim]-(X[0:dim-1]**2))**2+(X[0:dim-1]-1)**2)return Results'''F6函数'''
def F6(X):Results=np.sum(np.abs(X+0.5)**2)return Results'''F7函数'''
def F7(X):dim = X.shape[0]Temp = np.arange(1,dim+1,1)Results=np.sum(Temp*(X**4))+np.random.random()return Results'''F8函数'''
def F8(X):Results=np.sum(-X*np.sin(np.sqrt(np.abs(X))))return Results'''F9函数'''
def F9(X):dim=X.shape[0]Results=np.sum(X**2-10*np.cos(2*np.pi*X))+10*dimreturn Results'''F10函数'''
def F10(X):dim=X.shape[0]Results=-20*np.exp(-0.2*np.sqrt(np.sum(X**2)/dim))-np.exp(np.sum(np.cos(2*np.pi*X))/dim)+20+np.exp(1)return Results'''F11函数'''
def F11(X):dim=X.shape[0]Temp=np.arange(1,dim+1,+1)Results=np.sum(X**2)/4000-np.prod(np.cos(X/np.sqrt(Temp)))+1return Results'''F12函数'''
def Ufun(x,a,k,m):Results=k*((x-a)**m)*(x>a)+k*((-x-a)**m)*(x<-a)return Resultsdef F12(X):dim=X.shape[0]Results=(np.pi/dim)*(10*((np.sin(np.pi*(1+(X[0]+1)/4)))**2)+\np.sum((((X[0:dim-1]+1)/4)**2)*(1+10*((np.sin(np.pi*(1+X[1:dim]+1)/4)))**2)+((X[dim-1]+1)/4)**2))+\np.sum(Ufun(X,10,100,4))return Results'''F13函数'''
def Ufun(x,a,k,m):Results=k*((x-a)**m)*(x>a)+k*((-x-a)**m)*(x<-a)return Resultsdef F13(X):dim=X.shape[0]Results=0.1*((np.sin(3*np.pi*X[0]))**2+np.sum((X[0:dim-1]-1)**2*(1+(np.sin(3*np.pi*X[1:dim]))**2))+\((X[dim-1]-1)**2)*(1+(np.sin(2*np.pi*X[dim-1]))**2))+np.sum(Ufun(X,5,100,4))return Results'''F14函数'''
def F14(X):aS=np.array([[-32,-16,0,16,32,-32,-16,0,16,32,-32,-16,0,16,32,-32,-16,0,16,32,-32,-16,0,16,32],\[-32,-32,-32,-32,-32,-16,-16,-16,-16,-16,0,0,0,0,0,16,16,16,16,16,32,32,32,32,32]])bS=np.zeros(25)for i in range(25):bS[i]=np.sum((X-aS[:,i])**6)Temp=np.arange(1,26,1)Results=(1/500+np.sum(1/(Temp+bS)))**(-1)return Results'''F15函数'''
def F15(X):aK=np.array([0.1957,0.1947,0.1735,0.16,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246])bK=np.array([0.25,0.5,1,2,4,6,8,10,12,14,16])bK=1/bKResults=np.sum((aK-((X[0]*(bK**2+X[1]*bK))/(bK**2+X[2]*bK+X[3])))**2)return Results'''F16函数'''
def F16(X):Results=4*(X[0]**2)-2.1*(X[0]**4)+(X[0]**6)/3+X[0]*X[1]-4*(X[1]**2)+4*(X[1]**4)return Results'''F17函数'''
def F17(X):Results=(X[1]-(X[0]**2)*5.1/(4*(np.pi**2))+(5/np.pi)*X[0]-6)**2+10*(1-1/(8*np.pi))*np.cos(X[0])+10return Results'''F18函数'''
def F18(X):Results=(1+(X[0]+X[1]+1)**2*(19-14*X[0]+3*(X[0]**2)-14*X[1]+6*X[0]*X[1]+3*X[1]**2))*\(30+(2*X[0]-3*X[1])**2*(18-32*X[0]+12*(X[0]**2)+48*X[1]-36*X[0]*X[1]+27*(X[1]**2)))return Results'''F19函数'''
def F19(X):aH=np.array([[3,10,30],[0.1,10,35],[3,10,30],[0.1,10,35]])cH=np.array([1,1.2,3,3.2])pH=np.array([[0.3689,0.117,0.2673],[0.4699,0.4387,0.747],[0.1091,0.8732,0.5547],[0.03815,0.5743,0.8828]])Results=0for i in range(4):Results=Results-cH[i]*np.exp(-(np.sum(aH[i,:]*((X-pH[i,:]))**2)))return Results'''F20函数'''
def F20(X):aH=np.array([[10,3,17,3.5,1.7,8],[0.05,10,17,0.1,8,14],[3,3.5,1.7,10,17,8],[17,8,0.05,10,0.1,14]])cH=np.array([1,1.2,3,3.2])pH=np.array([[0.1312,0.1696,0.5569,0.0124,0.8283,0.5886],[0.2329,0.4135,0.8307,0.3736,0.1004,0.9991],\[0.2348,0.1415,0.3522,0.2883,0.3047,0.6650],[0.4047,0.8828,0.8732,0.5743,0.1091,0.0381]])Results=0for i in range(4):Results=Results-cH[i]*np.exp(-(np.sum(aH[i,:]*((X-pH[i,:]))**2)))return Results'''F21函数'''
def F21(X):aSH=np.array([[4,4,4,4],[1,1,1,1],[8,8,8,8],[6,6,6,6],[3,7,3,7],\[2,9,2,9],[5,5,3,3],[8,1,8,1],[6,2,6,2],[7,3.6,7,3.6]])cSH=np.array([0.1,0.2,0.2,0.4,0.4,0.6,0.3,0.7,0.5,0.5])Results=0for i in range(5):Results=Results-(np.dot((X-aSH[i,:]),(X-aSH[i,:]).T)+cSH[i])**(-1)return Results'''F22函数'''
def F22(X):aSH=np.array([[4,4,4,4],[1,1,1,1],[8,8,8,8],[6,6,6,6],[3,7,3,7],\[2,9,2,9],[5,5,3,3],[8,1,8,1],[6,2,6,2],[7,3.6,7,3.6]])cSH=np.array([0.1,0.2,0.2,0.4,0.4,0.6,0.3,0.7,0.5,0.5])Results=0for i in range(7):Results=Results-(np.dot((X-aSH[i,:]),(X-aSH[i,:]).T)+cSH[i])**(-1)return Results'''F23函数'''
def F23(X):aSH=np.array([[4,4,4,4],[1,1,1,1],[8,8,8,8],[6,6,6,6],[3,7,3,7],\[2,9,2,9],[5,5,3,3],[8,1,8,1],[6,2,6,2],[7,3.6,7,3.6]])cSH=np.array([0.1,0.2,0.2,0.4,0.4,0.6,0.3,0.7,0.5,0.5])Results=0for i in range(10):Results=Results-(np.dot((X-aSH[i,:]),(X-aSH[i,:]).T)+cSH[i])**(-1)return Results
把他们的参数设置都用字典装起来
Funobject = {'F1': F1,'F2': F2,'F3': F3,'F4': F4,'F5': F5,'F6': F6,'F7': F7,'F8': F8,'F9': F9,'F10': F10,'F11': F11,'F12': F12,'F13': F13,'F14': F14,'F15': F15,'F16': F16,'F17': F17,'F18': F18,'F19': F19,'F20': F20,'F21': F21,'F22': F22,'F23': F23}
Funobject.keys()#维度，搜索区间下界，搜索区间上界，最优值
Fundim={'F1': [30,-100,100],'F2': [30,-10,10],'F3': [30,-100,100],'F4': [30,-10,10],'F5': [30,-30,30],'F6': [30,-100,100],'F7': [30,-1.28,1.28],'F8': [30,-500,500],'F9':[30,-5.12,5.12],'F10': [30,-32,32],'F11': [30,-600,600],'F12': [30,-50,50],'F13': [30,-50,50],'F14': [2,-65,65],'F15':[4,-5,5],'F16': [2,-5,5],'F17':[2,-5,5],'F18': [2,-2,2],'F19': [3,0,1],'F20': [6,0,1],'F21':[4,0,10],'F22': [4,0,10],'F23': [4,0,10]}

Fundim字典里面装的是对应这个函数的，维度，搜索区间下界，搜索区间上界。这样写好方便我们去遍历测试所有的函数。

鲸鱼优化算法

终于到了算法的主代码阶段了：

import numpy as np
import random
import math
import copydef 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 Xdef 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 Xdef CaculateFitness(X,fun):'''计算种群的所有个体的适应度值'''pop = X.shape[0]fitness = np.zeros([pop, 1])for i in range(pop):fitness[i] = fun(X[i, :])return fitnessdef SortFitness(Fit):'''适应度值排序''''''输入为适应度值输出为排序后的适应度值，和索引'''fitness = np.sort(Fit, axis=0)index = np.argsort(Fit, axis=0)return fitness,indexdef SortPosition(X,index):'''根据适应度值对位置进行排序'''Xnew = np.zeros(X.shape)for i in range(X.shape[0]):Xnew[i,:] = X[index[i],:]return Xnewdef WOA(pop, dim, lb, ub, MaxIter, fun):'''鲸鱼优化算法''''''输入：pop:为种群数量dim:每个个体的维度ub:为个体上边界信息，维度为[1,dim]lb:为个体下边界信息，维度为[1,dim]fun:为适应度函数接口MaxIter:为最大迭代次数输出：GbestScore:最优解对应的适应度值GbestPositon:最优解Curve:迭代曲线'''X = initialization(pop,ub,lb,dim) #初始化种群fitness = CaculateFitness(X, fun)  # 计算适应度值fitness, sortIndex = SortFitness(fitness)  # 对适应度值排序X = SortPosition(X, sortIndex)  # 种群排序GbestScore = copy.copy(fitness[0])#记录最优适应度值GbestPositon = np.zeros([1,dim])GbestPositon[0,:] = copy.copy(X[0, :])#记录最优解Curve = np.zeros([MaxIter, 1])for t in range(MaxIter):print('第'+str(t)+'次迭代')Leader = X[0, :]  # 领头鲸鱼a = 2 - t * (2 / MaxIter)  # 线性下降权重2 - 0for i in range(pop):r1 = random.random()r2 = random.random()A = 2 * a * r1 - aC = 2 * r2b = 1l = 2 * random.random() - 1 #[-1,1]之间的随机数for j in range(dim):p = random.random()if p < 0.5:if np.abs(A) >= 1:#寻找猎物rand_leader_index = min(int(np.floor(pop * random.random() + 1)), pop - 1)#随机选择一个个体X_rand = X[rand_leader_index, :]D_X_rand = np.abs(C * X_rand[j] - X[i, j])X[i, j] = X_rand[j] - A * D_X_randelif np.abs(A) < 1:#包围猎物D_Leader = np.abs(C * Leader[j] - X[i, j])X[i, j] = Leader[j] - A * D_Leaderelif p >= 0.5:#气泡网攻击distance2Leader = np.abs(Leader[j] - X[i, j])X[i, j] = distance2Leader * np.exp(b * l) * np.cos(l * 2 * math.pi) + Leader[j]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[0,:] = copy.copy(X[0, :])Curve[t] = GbestScorereturn GbestScore, GbestPositon, Curve

其实优化算法差不多都是这个流程，边界函数，适应度函数排序，然后寻优过程等等。

OPT_algorithms = {'WOA':WOA}
OPT_algorithms.keys()

简单使用

我们选择F8来测试，先看看F8函数三维的情况：

'''F8绘图函数'''
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3Ddef F8(X):Results=np.sum(-X*np.sin(np.sqrt(np.abs(X))))return Resultsdef F8Plot():fig = plt.figure(1) #定义figureax = Axes3D(fig) #将figure变为3dx1=np.arange(-500,500,10) #定义x1，范围为[-500,500],间隔为10x2=np.arange(-500,500,10) #定义x2，范围为[-500,500],间隔为10X1,X2=np.meshgrid(x1,x2) #生成网格nSize = x1.shape[0]Z=np.zeros([nSize,nSize])for i in range(nSize):for j in range(nSize):X=[X1[i,j],X2[i,j]] #构造F8输入X=np.array(X) #将格式由list转换为arrayZ[i,j]=F8(X)  #计算F8的值#绘制3D曲面# rstride:行之间的跨度  cstride:列之间的跨度# rstride:行之间的跨度  cstride:列之间的跨度# cmap参数可以控制三维曲面的颜色组合ax.plot_surface(X1, X2, Z, rstride = 1, cstride = 1, cmap = plt.get_cmap('rainbow'))ax.contour(X1, X2, Z, zdir='z', offset=-1000)#绘制等高线ax.set_xlabel('X1')#x轴说明ax.set_ylabel('X2')#y轴说明ax.set_zlabel('Z')#z轴说明ax.set_title('F8_space')plt.show()F8Plot()

然后我们使用优化算法来寻优，自定义好所有的参数：

#设置参数
pop = 30 #种群数量
MaxIter = 200#最大迭代次数
dim = 30 #维度
lb = -500*np.ones([dim, 1]) #下边界
ub = 500*np.ones([dim, 1])#上边界
#选择适应度函数
fobj = F8
#原始算法
GbestScore,GbestPositon,Curve = WOA(pop,dim,lb,ub,MaxIter,fobj) 
#改进算法print('------原始算法结果--------------')
print('最优适应度值：',GbestScore)
print('最优解：',GbestPositon)

其实f8测试函数的最小值是-418.9829*30。所以可以看到这些结果还是有些差距的，，但是我那别的函数测试了一下，效果还是挺好的，WOA效果是很不错的。可能是再F8这个测试函数还是哪个表现不出来。

自己使用解决实际问题的时候只需要替换fobj这个目标函数的参数就可以了。

这个函数就如同上面所有的自定义的测试函数一样，你只需要定义输入的x，经过1系列实际问题的计算逻辑，返回的适应度值就可以。

绘制适应度曲线

#绘制适应度曲线
plt.figure(figsize=(6,2.7),dpi=128)
plt.plot(Curve,'b-',linewidth=2)
plt.xlabel('Iteration',fontsize='medium')
plt.ylabel("Fitness",fontsize='medium')
plt.grid()
plt.title('WOA',fontsize='large')
plt.legend(['WOA'], loc='upper right')
plt.show()

我这里是对数轴，但是也收敛了，在100左右基本就到0了，WOA还是很厉害的。

其实看到这里差不多就可以去把这个优化算法的函数拿去使用了，演示结束了，但是由于我们这里还需要对它的性能做一些测试，我们会把它在所有的测试函数上都跑一遍，这个时间可能是有点久的。

所有函数都测试一下

准备存储评价结果的数据框

functions = list(Funobject.keys())
algorithms = list(OPT_algorithms.keys())
columns = ['Mean', 'Std', 'Best', 'Worth']
index = pd.MultiIndex.from_product([functions, algorithms], names=['function_name', 'Algorithm_name'])
df_eval = pd.DataFrame(index=index, columns=columns)
df_eval.head()

索引和列名称都建好了，现在就是一个个跑，把值放进去就行了。

准备存储迭代图的数据框

df_Curve=pd.DataFrame(columns=index)
df_Curve

自定义训练函数

#定义训练函数
def train_fun(fobj_name=None,opt_algo_name=None, pop=30,MaxIter=200,Iter=30,show_fit=False):fundim=Fundim[fobj_name]  ; fobj=Funobject[fobj_name]dim=fundim[0]lb = fundim[1]*np.ones([dim, 1]) ; ub = fundim[2]*np.ones([dim, 1])opt_algo=OPT_algorithms[opt_algo_name]GbestScore_one=np.zeros([Iter])GbestPositon_one=np.zeros([Iter,dim])Curve_one=np.zeros([Iter,MaxIter])for i in range(Iter):GbestScore_one[i],GbestPositon_one[i,:],Curve_oneT =opt_algo(pop,dim,lb,ub,MaxIter,fobj)Curve_one[i,:]=Curve_oneT.Toneal_Mean=np.mean(GbestScore_one) #计算平均适应度值oneal_Std=np.std(GbestScore_one)#计算标准差oneal_Best=np.min(GbestScore_one)#计算最优值oneal_Worst=np.max(GbestScore_one)#计算最差值oneal_MeanCurve=Curve_one.mean(axis=0) #求平均适应度曲线#储存结果df_eval.loc[(fobj_name, opt_algo_name), :] = [oneal_Mean,oneal_Std, oneal_Best,oneal_Worst]df_Curve.loc[:,(fobj_name,opt_algo_name)]=oneal_MeanCurve#df_Curve[df_Curve.columns[(fobj_name,opt_algo_name)]] = oneal_MeanCurveif show_fit:print(f'{fobj_name}函数的{opt_algo_name}算法的平均适应度值是{oneal_Mean}，标准差{oneal_Std}，最优值{oneal_Best}，最差值{oneal_Worst}')

训练测试

#设置参数
pop = 30#种群数量
MaxIter = 100 #代次数
Iter = 30 #运行次数

计算，遍历所有的测试函数

#所有函数，所有算法全部一次性计算
for fobj_name in list(Funobject.keys()):for opt_algo_name in OPT_algorithms.keys():try:train_fun(fobj_name=fobj_name,opt_algo_name=opt_algo_name, pop=pop,MaxIter=MaxIter,Iter=Iter)print(f'{fobj_name}的{opt_algo_name}算法完成')except Exception as e: # 使用 except 来捕获错误print(f'{fobj_name}的{opt_algo_name}算法报错了：{e}') # 打印错误信息

查看计算出来的评价指标

df_eval

由于这里大部分的测试函数最优值都是零，我们可以看到。woa在很多函数上基本是可以找到接近最优值的，效果虽然不如SSA，SMA这些，但是还是不错的。和蝴蝶优化，海鸥优化算是一个等级的，中等效果。

画出迭代图

colors = ['darkorange', 'limegreen', 'lightpink', 'deeppink', 'red', 'cornflowerblue', 'grey']
markers = ['^', 'D', 'o', '*', 'X', 'p', 's']def plot_log_line(df_plot, fobj_name, step=10, save=False):plt.figure(figsize=(6, 3), dpi=128)for column, color, marker in zip(df_plot.columns, colors, markers):plt.semilogy(df_plot.index[::step], df_plot[column][::step].to_numpy(), color=color, marker=marker, label=column, markersize=4, alpha=0.7)plt.xlabel('Iterations')plt.ylabel('f')plt.legend(loc='best', fontsize=8)if save:plt.savefig(f'./图片/{fobj_name}不同迭代图.png', bbox_inches='tight')plt.show()# 使用示例
# plot_log_line(your_dataframe, 'example_plot')
for fobj_name in df_Curve.columns.get_level_values(0).unique():df1=df_Curve[fobj_name]print(f'{fobj_name}的不同算法效果对比：')plot_log_line(df1,fobj_name,5,False)   #保存图片-True