【数值优化】启发搜索之群体智能（4）-- 菌群

细菌觅食算法（Bacterial Foraging Optimization，BFO）在2002年被K.M.Passino提出，它是模仿Eeoli大肠杆菌在人体肠道内吞噬食物的行为而提出一种群智能优化算法。

菌群算法过程

细菌觅食算法主要包括三层循环，外层是驱散操作，中间层是繁殖操作，内层是趋化操作。算法的核心是内层的趋化性操作，它对应着细菌在寻找食物过程中所采取的方向选择策略（包括翻转和前进），对算法的收敛性有着极其重要的影响。
第一个循环：驱散
该操作是为了提高算法的全局搜索能力，因为当一个问题的解空间存在多个极值点时，细菌的群聚性就会使得算法容易陷入局部极值。这个过程是用新的个体来代替原有的个体，不同于复制操作，驱散操作是按照一定的概率P而发生的，当某个细菌符合迁移的条件时，该细菌将被随机分配到解空间中去。
第二个循环：繁殖
该操作主要模拟了细菌个体优胜劣汰的繁殖过程，该操作具体为按照适应度值对所有细菌进行排序，将较差的一半细菌清除，用较好的一半细菌代替，保证细菌总量不变。
第三个循环：驱化
该操作主要模拟了细菌的运动过程，包括向前移动和转向移动两个过程，每当细菌完成一次翻转后，检查适应度值是否改变，若适应度得到改善，细菌将沿同一方向继续移动若干步，如此循环直至适应度不再改善，或达到设定的移动步数临界值，此过程定义为前进

菌群算法的代码实现

import numpy as np
import random
import copy
import matplotlib.pyplot as plt
import math

#========================================================
#  辅助工具
#========================================================

class BFOIndividual:
    '''
    模拟细菌个体
    '''
    def __init__(self,  vardim, bound):
        '''
        个体初始化
        vardim: 变量维数
        bound: 变量取值范围
        '''
        self.vardim = vardim
        self.bound = bound
        self.fitness = 0.
        self.trials = 0.

    def generate(self):
        '''
        生成随机性
        '''
        len = self.vardim
        rnd = np.random.random(size=len)
        self.chrom = np.zeros(len)
        for i in range(0, len):
            self.chrom[i] = self.bound[0, i] + (self.bound[1, i] - self.bound[0, i]) * rnd[i]

    def GrieFunc(self, vardim, x, bound):
        '''
        差分算法
        '''
        s1 = 0.
        s2 = 1.
        for i in range(1, vardim + 1):
            s1 = s1 + x[i - 1] ** 2
            s2 = s2 * math.cos(x[i - 1] / math.sqrt(i))
        y = (1. / 4000.) * s1 - s2 + 1
        y = 1. / (1. + y)
        return y
    
    def calculateFitness(self):
        '''
        计算适应度
        '''
        self.fitness = self.GrieFunc(self.vardim, self.chrom, self.bound)

#========================================================
#  核心算法
#========================================================

class BFO:
    '''
    菌群觅食算法
    '''
    def __init__(self, sizepop, vardim, bound, params):
        '''
        sizepop: 种群规模
        vardim: 变量维数
        bound: 变量取值范围
        param: 其他参数(列表)
        '''
        self.sizepop = sizepop
        self.vardim = vardim
        self.bound = bound
        self.bestPopulation = []
        self.accuFitness = np.zeros(self.sizepop)
        self.population = []  # 种群
        self.fitness = np.zeros(self.sizepop)  # 适应性评分
        self.params = params
        self.history = np.zeros((self.params[0] * self.params[1] * self.params[2], 2))  # 演化历史

    def initialize(self):
        '''
        初始化菌群
        '''
        for i in range(0, self.sizepop):
            ind = BFOIndividual(self.vardim, self.bound)
            ind.generate()
            self.population.append(ind)

    def evaluation(self):
        '''
        迭代
        '''
        for i in range(0, self.sizepop):
            self.population[i].calculateFitness()
            self.fitness[i] = self.population[i].fitness

    def sortPopulation(self):
        '''
        种群排序
        '''
        sortedIdx = np.argsort(self.accuFitness)
        newpop = []
        newFitness = np.zeros(self.sizepop)
        for i in range(0, self.sizepop):
            ind = self.population[sortedIdx[i]]
            newpop.append(ind)
            self.fitness[i] = ind.fitness
        self.population = newpop

    def chemotaxls(self):
        '''
        趋化行为
        '''
        for i in range(0, self.sizepop):
            tmpInd = copy.deepcopy(self.population[i])
            self.fitness[i] += self.communication(tmpInd)
            Jlast = self.fitness[i]
            rnd = np.random.uniform(low=-1, high=1.0, size=self.vardim)
            phi = rnd / np.linalg.norm(rnd)
            tmpInd.chrom += self.params[4] * phi
            for k in range(0, self.vardim):
                if tmpInd.chrom[k] < self.bound[0, k]:
                    tmpInd.chrom[k] = self.bound[0, k]
                if tmpInd.chrom[k] > self.bound[1, k]:
                    tmpInd.chrom[k] = self.bound[1, k]
            tmpInd.calculateFitness()
            m = 0
            while m < self.params[3]:
                if tmpInd.fitness < Jlast:
                    Jlast = tmpInd.fitness
                    self.population[i] = copy.deepcopy(tmpInd)
                    # print m, Jlast
                    tmpInd.fitness += self.communication(tmpInd)
                    tmpInd.chrom += self.params[4] * phi
                    for k in range(0, self.vardim):
                        if tmpInd.chrom[k] < self.bound[0, k]:
                            tmpInd.chrom[k] = self.bound[0, k]
                        if tmpInd.chrom[k] > self.bound[1, k]:
                            tmpInd.chrom[k] = self.bound[1, k]
                    tmpInd.calculateFitness()
                    m += 1
                else:
                    m = self.params[3]
            self.fitness[i] = Jlast
            self.accuFitness[i] += Jlast

    def communication(self, ind):
        '''
        细菌间交流行为
        '''
        Jcc = 0.0
        term1 = 0.0
        term2 = 0.0
        for j in range(0, self.sizepop):
            term = 0.0
            for k in range(0, self.vardim):
                term += (ind.chrom[k] -
                         self.population[j].chrom[k]) ** 2
            term1 -= self.params[6] * np.exp(-1 * self.params[7] * term)
            term2 += self.params[8] * np.exp(-1 * self.params[9] * term)
        Jcc = term1 + term2

        return Jcc

    def reproduction(self):
        '''
        繁殖行为
        '''
        self.sortPopulation()
        newpop = []
        for i in range(0, self.sizepop // 2):
            newpop.append(self.population[i])
        for i in range(self.sizepop // 2, self.sizepop):
            self.population[i] = newpop[i - self.sizepop // 2]

    def eliminationAndDispersal(self):
        '''
        分散与淘汰
        '''
        for i in range(0, self.sizepop):
            rnd = random.random()
            if rnd < self.params[5]:
                self.population[i].generate()

    def printResult(self):
        '''
        绘制菌群迭代过程
        '''
        x = np.arange(0, self.t)
        y1 = self.history[:, 0]
        y2 = self.history[:, 1]
        plt.plot(x, y1, 'r', label='optimal value')
        plt.plot(x, y2, 'g', label='average value')
        plt.xlabel("Iteration")
        plt.ylabel("function value")
        plt.title("Baterial clony foraging algorithm for function optimization")
        plt.legend()
        plt.show()

    def main(self):
        '''
        主流程:蜂群搜索过程
        '''
        # 迭代次数
        self.t = 0
        # 种群初始化
        self.initialize()
        # 评估适应性得分
        self.evaluation()
        # 评估迭代结果
        bestIndex = np.argmin(self.fitness)
        self.best = copy.deepcopy(self.population[bestIndex])
        # 三层循环
        for i in range(0, self.params[0]):
            for j in range(0, self.params[1]):
                for k in range(0, self.params[2]):
                    self.t += 1
                    self.chemotaxls()
                    self.evaluation()
                    best = np.min(self.fitness)
                    bestIndex = np.argmin(self.fitness)
                    if best < self.best.fitness:
                        self.best = copy.deepcopy(self.population[bestIndex])
                    self.avefitness = np.mean(self.fitness)
                    self.history[self.t - 1, 0] = self.best.fitness
                    self.history[self.t - 1, 1] = self.avefitness
                    print("Generation %d: optimal function value is: %f; average function value is %f" % (self.t, self.history[self.t - 1, 0], self.history[self.t - 1, 1]))
                self.reproduction()
            self.eliminationAndDispersal()

        # 打印最优结果
        print("Optimal function value is: %f; " % self.history[self.t - 1, 0])
        # 打印最优个体参数信息
        print("Optimal solution is:")
        print(self.best.chrom)
        # 绘制历史
        self.printResult()

if __name__ == "__main__":
    bound = np.tile([[-600], [600]], 25)
    bfo = BFO(60, 25, bound, [2, 2, 50, 4, 50, 0.25, 0.1, 0.2, 0.1, 10])
    bfo.main()