智能算法|Python实现麻雀优化算法(SSA)

背景介绍

麻雀搜索算法[1](Sparrow Search Algorithm, SSA)于2020年提出，主要通过模仿麻雀的觅食行为和反捕食行为实现位置寻优，以找到部分NP问题的局部最优值。

在该算法的预设中，麻雀种群内部被分为发现者和跟随者两种角色，同时模仿真实的捕食情景，增加了麻雀的危险预警机制。

问题定义

下面以一个2维平面搜索问题为例，对SSA进行介绍。

假设我们需要解决的问题是计算给定范围内 $x_1\in [lb,ub],x_2 \in [lb,ub]$ ，两个数字 $x_1$ 和 $x_2$ 的平方和最小值。其中lb为搜索空间的下界，ub为搜索空间的上界。

将每只麻雀视为2维平面上的一个点，该点的横纵坐标 $x_1$ 、 $x_2$ 即为该麻雀的位置。当前位置好坏通过计算适应度 $fit = x^2+y^2$ 的大小来评价。

则上述的问题可以抽象成如下的数学公式：

$\mathop{min}\limits_{x,y} fit = x^2+y^2$

麻雀搜索算法实现

针对上述的问题，寻优的目标具有最小 $fit$ 的麻雀位置 $(x_1,x_2)$ 。则具体的搜索过程按照如下步骤进行。

1. 麻雀种群初始化

假设该种群共有麻雀数量为 $pop=50$ ，则该种群可用如下矩阵表示：

$X=\left[ \begin{matrix} {{x}_{1,1}} & {{x}_{1,2}} & \cdots & {{x}_{1,d}} \\ {{x}_{2,1}} & {{x}_{2,2}} & \cdots & {{x}_{2,d}} \\ \vdots & \vdots & \ddots & \vdots \\ {{x}_{n,1}} & {{x}_{n,2}} & \cdots & {{x}_{n,d}} \\ \end{matrix} \right]$

其中，d=2，n=50。

该过程对应的 Python 代码如下：

#载入所需的包
import numpy as np
import random
#初始化麻雀种群
def initial(pop, dim, ub, lb):
    X = np.zeros([pop, dim])
    for i in range(pop):
        X[i,:] = np.random.uniform(low=lb[0], high=ub[0], size=(1, dim))
    return X, lb, ub

则初始种群所对应的适应度 $F_X$ 为：

$F_X=\left[ \begin{matrix} f({{x}_{1,1}} & {{x}_{1,2}} & \cdots & {{x}_{1,d}}) \\ f({{x}_{2,1}} & {{x}_{2,2}} & \cdots & {{x}_{2,d}}) \\ \vdots & \vdots & \ddots & \vdots \\ f({{x}_{n,1}} & {{x}_{n,2}} & \cdots & {{x}_{n,d}}) \\ \end{matrix} \right]$

'''定义适应度函数'''
def fun(X):
    O = 0
    for i in X:
        O += i ** 2
    return O

'''计算适应度函数'''
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

将初始化得到的麻雀种群按照适应度的大小进行排序，则得到的具有最优适应度值的麻雀为 X[0,:]。

'''适应度排序'''
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

接下来，对初始生成的麻雀种群按照发现者更新、追随者更新、危险预警的公式进行位置更新。

2. 麻雀种群更新

这群麻雀中有N只麻雀,每代选取种群中位置最好的PN只麻雀作为发现者,剩余的N-PN只麻雀作为跟随者。

2.1 发现者更新

发现者的更新公式如下：

$x_{i,j}^{t+1}=\left\{ \begin{align} & \begin{matrix} x_{i,j}^{t}\cdot \exp \left( \frac{-i}{\alpha \cdot ite{{r}_{\max }}} \right) & if\text{ }{{R}_{2}}<ST \\ \end{matrix} \\ & \begin{matrix} x_{i,j}^{t}+Q & if\text{ }{{R}_{2}}<ST \\ \end{matrix} \\ \end{align} \right.$

其中， $iter_{max}$ 是预先设定的最大迭代次数， $\alpha$ 为(0，1]中的均匀随机数，Q为一个标准正态分布随机数。

'''麻雀发现者勘探更新'''
def PDUpdate(X, PDNumber, ST, Max_iter, dim):
    X_new = copy.copy(X)
    R2 = random.random()
    for p in range(PDNumber):
        for j in range(dim)
            if R2 < ST:
                X_new[p, j] = X[p, j] * np.exp(-p / (random.random() * Max_iter))
            else:
               X_new[p, j] = X[p, j] + np.random.randn()
    return X_new

2.2 追随者更新

追随者的更新公式如下：
$x_{i,j}^{t+1}=\left\{ \begin{align} & \begin{matrix} Q\cdot \exp \left( \frac{x_{worst}^{t}-x_{i,j}^{t}}{{{i}^{2}}} \right) & if\text{ }i>n/2 \\ \end{matrix} \\ & \begin{matrix} x_{best}^{t}-\frac{1}{d}\sum\limits_{j=1}^{d}{\left| x_{i,j}^{t}-x_{best}^{t} \right|\cdot rand\left( \left\{ -1,1 \right\} \right)} & otherwise \\ \end{matrix} \\ \end{align} \right.$
其中， $x_{best}^{t}$ 为当前种群中具有最优适应度的麻雀， $x_{worst}^{t}$ 为当前种群中具有最差适应度的麻雀。

'''麻雀加入者勘探更新'''
def JDUpdate(X, PDNumber, pop, dim):
    X_new = copy.copy(X)
    # 产生-1，1的随机数
    A = np.ones([dim, 1])
    for a in range(dim):
        if (random.random() > 0.5):
            A[a] = -1
    for i in range(PDNumber + 1, pop):
        for j in range(dim):
            if i > (pop - PDNumber) / 2 + PDNumber:
                X_new[i, j] = np.random.randn() * np.exp((X[-1, j] - X[i, j]) / i ** 2)
            else:
                AA = np.mean(np.abs(X[i, :] - X[0, :])*A)
                X_new[i, j] = X[0, j] - AA
    return X_new

2.3 危险预警

危险预警的更新公式如下：
$x_{i,j}^{t+1}=\left\{ \begin{align} & \begin{matrix} x_{i,j}^{best}+\beta \cdot \left( x_{i,j}^{t}-x_{i,j}^{best} \right) & {{f}_{i}}\ne {{f}_{g}} \\ \end{matrix} \\ & \begin{matrix} x_{i,j}^{t}+K\cdot \left( \frac{x_{i,j}^{t}-x_{i,j}^{worst}}{\left| {{f}_{i}}-{{f}_{worst}} \right|+\varepsilon } \right) & {{f}_{i}}={{f}_{g}} \\ \end{matrix} \\ \end{align} \right.$

其中 $\beta$ 为符合正态分布的随机数， $K$ 为[-1,1]之间的随机数， $\varepsilon$ 为一个较小的数字，防止分母为0。 ${f}_{worst}$ 为具有当前种群的最差适应度值。 ${f}_{g}$ 为全局最优适应度值。

'''危险更新'''
def SDUpdate(X, pop, SDNumber, fitness, BestF):
    X_new = copy.copy(X)
    dim = X.shape[1]
    Temp = range(pop)
    RandIndex = random.sample(Temp, pop)
    SDchooseIndex = RandIndex[0:SDNumber]
    for i in range(SDNumber):
        for j in range(dim):
            if fitness[SDchooseIndex[i]] > BestF:
                X_new[SDchooseIndex[i], j] = X[0, j] + np.random.randn() * np.abs(X[SDchooseIndex[i], j] - X[0, j])
            elif fitness[SDchooseIndex[i]] == BestF:
                K = 2 * random.random() - 1
                X_new[SDchooseIndex[i], j] = X[SDchooseIndex[i], j] + K * (
                        np.abs(X[SDchooseIndex[i], j] - X[-1, j]) / (fitness[SDchooseIndex[i]] - fitness[-1] + 10E-8))
    return X_new

3. 完整的麻雀优化算法

3.1 麻雀算法流程图

麻雀算法流程图

3.1 完整的麻雀算法

import copy
import random
import numpy as np

''' Tent种群初始化函数 '''
def initial(pop, dim, ub, lb):
    X = np.zeros([pop, dim])
    for i in range(pop):
        for j in range(dim):
            X[i, j] = np.random.rand() * (ub[j] - lb[j]) + lb[j]
    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] = np.random.rand() * (ub[j] - lb[j]) + lb[j]
            elif X[i,j]<lb[j]:
                X[i, j] = np.random.rand() * (ub[j] - lb[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):
    X_new = copy.copy(X)
    R2 = random.random()
    for p in range(PDNumber):
        for j in range(dim):
            if R2 < ST:
                X_new[p, j] = X[p, j] * np.exp(-p / (random.random() * Max_iter))
            else:
               X_new[p, j] = X[p, j] + np.random.randn()
    return X_new
        
'''麻雀加入者更新'''            
def JDUpdate(X, PDNumber, pop, dim):
    X_new = copy.copy(X)
    # 产生-1，1的随机数
    A = np.ones([dim, 1])
    for a in range(dim):
        if (random.random() > 0.5):
            A[a] = -1
    for i in range(PDNumber + 1, pop):
        for j in range(dim):
            if i > (pop - PDNumber) / 2 + PDNumber:
                X_new[i, j] = np.random.randn() * np.exp((X[-1, j] - X[i, j]) / i ** 2)
            else:
                AA = np.mean(np.abs(X[i, :] - X[0, :])*A)
                X_new[i, j] = X[0, j] - AA
    return X_new
            
'''危险更新'''   
def SDUpdate(X, pop, SDNumber, fitness, BestF):
    X_new = copy.copy(X)
    dim = X.shape[1]
    Temp = range(pop)
    RandIndex = random.sample(Temp, pop)
    SDchooseIndex = RandIndex[0:SDNumber]
    for i in range(SDNumber):
        for j in range(dim):
            if fitness[SDchooseIndex[i]] > BestF:
                X_new[SDchooseIndex[i], j] = X[0, j] + np.random.randn() * np.abs(X[SDchooseIndex[i], j] - X[0, j])
            elif fitness[SDchooseIndex[i]] == BestF:
                K = 2 * random.random() - 1
                X_new[SDchooseIndex[i], j] = X[SDchooseIndex[i], j] + K * (
                        np.abs(X[SDchooseIndex[i], j] - X[-1, j]) / (fitness[SDchooseIndex[i]] - fitness[-1] + 10E-8))
    return X_new
              
    

'''麻雀搜索算法'''
def Tent_SSA(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):
        BestF = fitness[0]
        
        X = PDUpdate(X,PDNumber,ST,Max_iter,dim)#发现者更新

        X = JDUpdate(X,PDNumber,pop,dim) #加入者更新

        X = SDUpdate(X,pop,SDNumber,fitness,BestF) #危险更新

        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[i] = GbestScore
    return GbestScore,GbestPositon,Curve

3.2 结果展示

结果	取值
最优值	1.03409414e-08
最优位置	(0.00048558,0.00023213)

迭代过程图

Tips

文内公式部分参考：stronghorse
文内代码部分参考：智能算法研学社（Jack旭）

最后编辑于：2023.01.28 10:34:07

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 200,841评论 5赞 472
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 84,415评论 2赞 377
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 147,904评论 0赞 333
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,051评论 1赞 272
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,055评论 5赞 363
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,255评论 1赞 278
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,729评论 3赞 393
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,377评论 0赞 255
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,517评论 1赞 294
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,420评论 2赞 317
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,467评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,144评论 3赞 317
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,735评论 3赞 303
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,812评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,029评论 1赞 256
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,528评论 2赞 346
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,126评论 2赞 341