Python优化算法22——自适应变异麻雀搜索优化算法(AMSSA)

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

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

之前写过黏菌优化算法的文章，现在有很多新的黏菌优化算法，都是进行了一些改进。本次带来的是自适应变异麻雀搜索算法，当然也会和普通的麻雀优化进行对比。

这是文章的摘要：

针对麻雀搜索算法前期易陷入局部极值点、后期寻优精度不高等问题，提出一种

自适应变异麻雀搜索算法（AMSSA）。先通过猫映射混沌序列初始化种群，增强初始种群的随机

性、遍历性，提高算法的全局搜索能力；再引入柯西变异和Tent 混沌扰动，拓展局部搜索能力，使

陷入局部极值点的个体跳出限制继续搜索；最后，提出探索者-跟随者数量自适应调整策略，利用各

阶段探索者和跟随者数量的改变增强算法前期的全局搜索能力和后期的局部深度挖掘能力，提高算

法的寻优精度。选取16 个基准函数和Wilcoxon 检验进行验证，实验结果表明：所提算法与其他算

法相比，寻优精度、收敛速度和稳定性都取得较大提升。

文章和理论看看就行，主要是代码。本文都是python实现。文章和代码文件都在文末有获取方式。

代码实现

导入包

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 
import seaborn as sns
import warnings
import copy

plt.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=0
    for i in range(dim):
        Results=Results+np.sum(X[0:i+1])**2
    return 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*dim
    return 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)))+1
    return Results
 
'''F12函数'''
def Ufun(x,a,k,m):
    Results=k*((x-a)**m)*(x>a)+k*((-x-a)**m)*(x<-a)
    return Results
 
def 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 Results
 
def 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/bK
    Results=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])+10
    return 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=0
    for 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=0
    for 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=0
    for 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=0
    for 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=0
    for 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 copy


''' 改进点：cat混沌种群初始化函数 '''
def initial(pop, dim, ub, lb):
    X = np.zeros([pop, dim])
    x=np.random.random()
    y=np.random.random()
    N=1
    a=1
    b=1
    for i in range(pop):
        for j in range(dim):
            x1=np.mod(x+b*y,N)
            y1=np.mod(a*x+a*b*y,N)
            X[i, j] = x1*(ub[j] - lb[j]) + lb[j]
            x=x1;
            y=y1;
    
    return X,lb,ub
            
'''边界检查函数'''
def BorderCheck(X,ub,lb,pop,dim):
    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 PDUpdate(X,PDNumber,ST,Max_iter,dim,t):
    X_new  = copy.copy(X)
    R2 = random.random()
    for j in range(PDNumber):
        if R2<ST:
            #改进点：改进探索者位置更新公式
            X_new[j,:] = X[j,:]*(2/np.exp(4*t/(np.random.random()*Max_iter)**2))
        else:
            X_new[j,:] = X[j,:] + np.random.randn()*np.ones([1,dim])
    return X_new
        
'''麻雀加入者更新'''            
def JDUpdate(X,PDNumber,pop,dim):
    X_new = copy.copy(X)
    for j in range(PDNumber+1,pop):
         if j>(pop - PDNumber)/2 + PDNumber:
             X_new[j,:]= np.random.randn()*np.exp((X[-1,:] - X[j,:])/j**2)
         else:
             #产生-1，1的随机数
             A = np.ones([dim,1])
             for a in range(dim):
                 if(random.random()>0.5):
                     A[a]=-1       
         AA = np.dot(A,np.linalg.inv(np.dot(A.T,A)))
         X_new[j,:]= X[1,:] + np.abs(X[j,:] - X[1,:])*AA.T
           
    return X_new                    
            
'''危险更新'''   
def SDUpdate(X,pop,SDNumber,fitness,BestF):
    X_new = copy.copy(X)
    Temp = range(pop)
    RandIndex = random.sample(Temp, pop)
    SDchooseIndex = RandIndex[0:SDNumber]
    for j in range(SDNumber):
        if fitness[SDchooseIndex[j]]>BestF:
            X_new[SDchooseIndex[j],:] = X[0,:] + np.random.randn()*np.abs(X[SDchooseIndex[j],:] - X[1,:])
        elif fitness[SDchooseIndex[j]] == BestF:
            K = 2*random.random() - 1
            X_new[SDchooseIndex[j],:] = X[SDchooseIndex[j],:] + K*(np.abs( X[SDchooseIndex[j],:] - X[-1,:])/(fitness[SDchooseIndex[j]] - fitness[-1] + 10E-8))
    return X_new
              
    

'''改进麻雀搜索算法'''
def AMSSA(pop,dim,lb,ub,Max_iter,fun):
    ST = 0.6 #预警值
    PD = 0.7 #发现者的比列，剩下的是加入者
    SD = 0.2 #意识到有危险麻雀的比重
    PDNumber = int(pop*PD) #发现者数量
    SDNumber = int(pop*SD) #意识到有危险麻雀数量
    X,lb,ub = initial(pop, dim, ub, lb) #初始化种群
    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([Max_iter,1])
    for i in range(Max_iter):
        #改进点：比例系数改进
        ST = 0.6 #预警值
        b = 0.7 #比例系数
        k = 0.1 #扰动因子
        r = b*(np.tan(-np.pi*i/(4*Max_iter) + np.pi/4))-k*np.random.random();
        #防止r为负数或者大于1的数
        if r<0:
            r=0.1
        if r>1:
            r=1
        PDNumber=int(np.round(pop*r))
        SDNumber=int(np.round(pop*(1-r)))
        BestF = fitness[0]
        
        X = PDUpdate(X,PDNumber,ST,Max_iter,dim,i)#发现者更新
        
        X = JDUpdate(X,PDNumber,pop,dim) #加入者更新
        
        X = SDUpdate(X,pop,SDNumber,fitness,BestF) #危险更新
        
        X = BorderCheck(X,ub,lb,pop,dim) #边界检测
        
        fitness = CaculateFitness(X,fun) #计算适应度值
        #改进点：tent扰动和柯西扰动
        Avgf=np.mean(fitness)
        x0=np.random.random()
        for j in range(pop):
            if fitness[j]<Avgf:#柯西变异
                Temp=np.zeros([1,dim])
                Temp[0,:]=copy.copy(X[j,:])
                r=np.tan((np.random.random()-0.5)*np.pi)
                Temp[0,:]=Temp[0,:]+r*Temp[0,:]
                for jj in range(dim):
                    if Temp[0,jj]>ub[jj]:
                        Temp[0,jj]=ub[jj]
                    if Temp[0,jj]<lb[jj]:
                        Temp[0,jj]=lb[jj]
                fitTemp=fun(Temp[0]).mean()
                #print(fitTemp,fitness.shape)

                if fitTemp<fitness[j]:
                    fitness[j]=fitTemp
                    X[j,:]=copy.copy(Temp)
            else:#Tent 扰动
                if x0<0.5:
                    tentV=2*x0+np.random.random()/pop
                else:
                    tentV=2*(1-x0)+np.random.random()/pop
                Temp=np.zeros([1,dim])
                Temp[0,:]=copy.copy(X[j,:])
                Temp[0,:]=Temp[0,:]+tentV*Temp[0,:]
                for jj in range(dim):
                    if Temp[0,jj]>ub[jj]:
                        Temp[0,jj]=ub[jj]
                    if Temp[0,jj]<lb[jj]:
                        Temp[0,jj]=lb[jj]
                fitTemp=fun(Temp[0]).mean()
                #print(fitTemp,fitness.shape)
                if fitTemp<fitness[j]:
                    fitness[j]=fitTemp
                    X[j,:]=copy.copy(Temp)
                x0=tentV
                
        fitness,sortIndex = SortFitness(fitness) #对适应度值排序
        X = SortPosition(X,sortIndex) #种群排序
        if(fitness[0]<=GbestScore): #更新全局最优
            GbestScore = copy.copy(fitness[0])
            GbestPositon[0,:] = copy.copy(X[0,:])
        Curve[i] = GbestScore
    
    return GbestScore,GbestPositon,Curve

普通麻雀优化(SSA)：

import numpy as np
import copy
import random
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 SSA(pop,dim,lb,ub,Max_iter,fun):
    '''麻雀搜索算法'''
    '''
    输入：
    pop:为种群数量
    dim:每个个体的维度
    ub:为个体上边界信息，维度为[1,dim]
    lb:为个体下边界信息，维度为[1,dim]
    fun:为适应度函数接口
    MaxIter:为最大迭代次数
    输出：
    GbestScore:最优解对应的适应度值
    GbestPositon:最优解
    Curve:迭代曲线
    '''
    
    ST = 0.8 #预警值
    PD = 0.2 #发现者的比列，剩下的是加入者
    SD = 0.1 #意识到有危险麻雀的比重
    PDNumber = int(pop*PD) #发现者数量
    SDNumber = int(pop*SD) #意识到有危险麻雀数量
    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([Max_iter,1])
    for t in range(Max_iter):
        print("第"+str(t)+"次迭代")
        BestF = copy.copy(fitness[0])
        Xworst = copy.copy(X[-1,:])
        Xbest = copy.copy(X[0,:])
        '''发现者位置更新'''
        R2 = np.random.random()
        for i in range(PDNumber):
            if R2<ST:
                X[i,:] = X[i,:]*np.exp(-i/(np.random.random()*Max_iter))
            else:
                X[i,:] = X[i,:] + np.random.randn()*np.ones([1,dim])
        
        X = BorderCheck(X,ub,lb,pop,dim) #边界检测   
        fitness = CaculateFitness(X,fun) #计算适应度值
        
        bestII=np.argmin(fitness)
        Xbest = copy.copy(X[bestII,:])
        '''加入者位置更新'''
        for i in range(PDNumber+1,pop):
             if i>(pop - PDNumber)/2 + PDNumber:
                 X[i,:]= np.random.randn()*np.exp((Xworst - X[i,:])/i**2)
             else:
                 #产生-1，1的随机数
                 A = np.ones([dim,1])
                 for a in range(dim):
                     if(np.random.random()>0.5):
                         A[a]=-1       
                 AA = np.dot(A,np.linalg.inv(np.dot(A.T,A)))
                 X[i,:]= X[0,:] + np.abs(X[i,:] - GbestPositon)*AA.T
        
        X = BorderCheck(X,ub,lb,pop,dim) #边界检测   
        fitness = CaculateFitness(X,fun) #计算适应度值
        '''意识到危险的麻雀更新'''
        Temp = range(pop)
        RandIndex = random.sample(Temp, pop) 
        SDchooseIndex = RandIndex[0:SDNumber]#随机选取对应比列的麻雀作为意识到危险的麻雀
        for i in range(SDNumber):
            if fitness[SDchooseIndex[i]]>BestF:
                X[SDchooseIndex[i],:] = Xbest + np.random.randn()*np.abs(X[SDchooseIndex[i],:] - Xbest)
            elif fitness[SDchooseIndex[i]] == BestF:
                K = 2*np.random.random() - 1
                X[SDchooseIndex[i],:] = X[SDchooseIndex[i],:] + K*(np.abs( X[SDchooseIndex[i],:] - X[-1,:])/(fitness[SDchooseIndex[i]] - fitness[-1] + 10E-8))
        
        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] = GbestScore
    
    return GbestScore,GbestPositon,Curve

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

两个算法都放入字典，后面批量化测试。

OPT_algorithms = {'AMSSA':AMSSA,'SSA':SSA}
OPT_algorithms.keys()

简单使用

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

'''F14绘图函数'''
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def F14Plot():
    fig = plt.figure(1) #定义figure
    ax = Axes3D(fig) #将figure变为3d
    x1=np.arange(-65,65,2) #定义x1，范围为[-65,65],间隔为2
    x2=np.arange(-65,65,2) #定义x2，范围为[-65,65],间隔为2
    X1,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]] #构造F14输入
            X=np.array(X) #将格式由list转换为array
            Z[i,j]=F14(X)  #计算F14的值
    #绘制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=0)#绘制等高线
    ax.set_xlabel('X1')#x轴说明
    ax.set_ylabel('X2')#y轴说明
    ax.set_zlabel('Z')#z轴说明
    ax.set_title('F14_space')
    plt.show()

F14Plot()

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

#设置参数
pop = 30 #种群数量
MaxIter = 200#最大迭代次数
dim = 2 #维度
lb = -100*np.ones([dim, 1]) #下边界
ub = 100*np.ones([dim, 1])#上边界
#选择适应度函数
fobj = F14
#原始算法
GbestScore,GbestPositon,Curve = AMSSA(pop,dim,lb,ub,MaxIter,fobj) 
#改进算法

print('------原始算法结果--------------')
print('最优适应度值：',GbestScore)
print('最优解：',GbestPositon)

其实f14测试函数的最小值是0。所以可以看到这个算法几乎达到了最优，效果看起来还行。

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

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

绘制适应度曲线

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

我这里是对数轴，但是也收敛了，30轮左右最小值。

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

所有函数都测试一下

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

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.T
    
    oneal_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_MeanCurve
    if 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

由于这里大部分的测试函数最优值都是零，我们可以看到。两个算法在寻优的最优值情况都是差不多的，其中AMSSA会稍微好一点点,基本能找到比SSA更小的情况，当然也不一定，例如F6和F5,效果还不一定有原始版的好。

画出迭代图

评价一个优化算法肯定不能只从最优的情况来看，还要看他训练轮数寻优收敛的时间轮数，来综合评价。

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