基于距离浓度的人工免疫算法（DCAIA）求函数最值（python实现）

基于距离浓度的人工免疫算法（DCAIA）求函数最值（python实现）

免疫算法是人工免疫系统形式中应用最广法的一个部分，是兴起于2012年左右的机器学习算法，大多数的免疫算法被用于实现全局优化问题，包括组合优化、函数优化、神经网络优化等。

一、免疫算法理论

免疫算法是模仿生物免疫机制，结合基因的进化机理，人工构造出的一种新型智能优化算法。采用群体搜索策略，通过迭代计算，最终以较大的概率得到问题的最优解。相比于其他算法，免疫算法利用自身产生多样性和维持机制的特点，保证了群体多样性，克服了‘早熟’问题，可以得到全局最优解。具有自适应性，随机性并行性，全局收敛性，种群多样性等特点。

二、基于距离浓度算法的免疫机理

免疫系统是高度进化的生物系统，它的功能是针对外来抗原的输入产生识别并消灭抗原的抗体，本文提出的算法（DCAIA），借鉴了如下生物免疫特性和机理：

（1）抗体对抗原的识别是基于抗原的表位与抗体的受体的受体之间的匹配，抗体与抗原的亲和度决定了两者的匹配程度。

（2）免疫系统中同一种抗体的浓度是受到限制的，当某种抗体数量增加，达到一定浓度时，系统要对其进行抑制，同样对浓度低的抗体要进行繁殖以保证其多样性。

（3）抗体在繁殖过程中会产生高频变异，即克隆变异。以便可以识别种类繁多的抗原。

（4）抗体在识别抗原过程中，可以产生记忆，使得该抗体更快地识别抗原。在识别过程中，通过克隆变异和学习记忆，提高抗体对抗原的匹配度，亲和度成熟。

三、免疫算法对应关系

抗原 优化问题

抗体 优化问题的可行解

亲和度 可行解的质量

细胞活化 免疫选择

细胞分化 个体克隆

亲和度成熟 变异

动态维持平衡 种群刷新

四、免疫算法算子

编码方式主要为实数编码和离散编码。本文章利用基于距离浓度的人工免疫算法求解函数最值，采用实数编码。

（1）亲和度评价算子

表征抗体与抗原的结合度。通常是一个函数 $$aff(x)$$，函数的输入为一个抗体（可行解）。

（2）抗体浓度评价算子

抗体浓度通常定义为：
$$den(a{b}_i)=\frac{1}{N}\sum _{j=1}^{N}S(a{b}_i,a{b}_j)$$
式中N为种群规模， $$S(a{b}_i,a{b}_j)$$表示抗体间相似度：

$$\begin{gathered} S(a{b}_i,a{b}_j)=\{\begin{array}{l}1,dis(a{b}_i,a{b}_j)\le \delta \\ 0,dis(a{b}_i,a{b}_j)>\delta \end{array} \\ 式中dis()函数表示欧氏距离，\delta 为浓度阈值 \end{gathered}$$

（3）激励度计算算子

抗体鼓励度是对抗体质量的最终评价结果，通常亲和度大、浓度低的抗体会得到较大的鼓励度。

本文使用以下公式：
$$sim(a{b}_i)=\alpha \times aff(a{b}_i)-\beta \times den(a{b}_i)$$

（4）免疫选择算子

免疫选择算子根据抗体的激励度确定哪些抗体进入克隆选择操作。

（5）克隆算子

克隆算子将选中的抗体个体进行复制。

（6）变异算子

变异算子对克隆算子结果进行变异操作，以产生亲和度突变，实现局部搜索。变异算子在免疫算法中产生有潜力的新抗体，实现区域搜索的重要算子，对算法影响很大。本文采用以下公式：
$$\begin{gathered} T_m(a{b}_{i,j,m})=\{\begin{array}{l}a{b}_{i,j,m}+(rand-0.5)\ast \epsilon ,rand<p_m\\ a{b}_{i,j,m},其他\end{array} \\ \epsilon 表示定义的邻域范围，p_m为变异概率 \end{gathered}$$

（7）克隆抑制算子

克隆抑制算子用于对经过变异后的克隆体再进行选择，抑制亲和度低的抗体。

五、python实例

例（1）

计算函数$$f(x)=\sum _{i=1}^n{x}_i^2(-20\le x_i\le 20)$$的最小值，其中个体x的维数n=10，这是一个简单的平方和函数，只有一个极小点x=（0,0,…,0），理论上最小值分f(0,0,…,0)=0。

#亲和度：与目标函数相关
def affinity(x):
    return np.sum(x**2)

#激励度：根据亲和度和浓度计算，影响变异
def motivation(x,num):
    ND = np.zeros((1,num)) #临时数组
    nd = np.zeros((1,num)) #浓度
    for i in range(num):
        for j in range(num):
            tmp = np.sqrt(np.sum((x[:,i]-x[:,j])**2))
            if(tmp > detas):
                ND[0,i] = 0
            else:
                ND[0,i] = 1                
        nd[0,i] = np.sum(ND,axis=1)/num
    ## 激励度
    MSLL = np.zeros((1,num))
    for i in range(num):
        tmp = affinity(x[:,i])
        MSLL[0,i] = tmp
    mtv = alphe * MSLL - beta * nd
    return mtv

#变异函数
def variation(x , gen):
    temp = num//2
    first_x = np.zeros((D,temp))
    for i in range(temp):
        tx=x[:,i].reshape((10,1))
        tx_rpt=np.tile(tx,10)
        randt=np.random.rand(1)
        #以pm为概率变异
        if randt > pm:
            for g in range(D):
                for h in range(D):
                    tx_rpt[g,h] = tx_rpt[g,h] + (np.random.rand(1) - 0.5) * (defalut / (gen+1))
        for j in range(D):
            for k in range(D):
                if (tx_rpt[j,k] > top) or (tx_rpt[j,k] < bottom):
                    tx_rpt[j,k] = np.random.rand(1) * (top - bottom) + bottom
        tx_rpt[:,0] = tx[:,0]
        mvt = motivation(tx_rpt,10)
        index = np.argsort(mvt).flatten()
        temp = tx_rpt[:,index]
        first_x[:,i] = temp[:,0]
    return first_x

import numpy as np
import matplotlib.pyplot as plt

D = 10 #维度
num = 100 #个体数量
top = 20 #取值上限
bottom = -20 #取值下限
alphe = 3 #激励度系数
beta = 1 #激励度系数
detas = 0.3 #浓度阈值
pm = 0.6   #变异几率
trace = [] #记录过程
defalut = top 
x = np.random.rand(D,num) * (top - bottom) + bottom  #10*100初始种群

gen = 1
iteration = 300 #最大迭代代数
while gen<iteration:
    
    mtv = motivation(x,num) 
    #按激励度排序
    index = np.argsort(mtv).flatten()
    sort_x = x[:,index]
    #变异
    sort_x[:,:50] = variation(sort_x[:,:50],gen)
    sort_x[:,50:] = np.random.rand(D,50) * (top - bottom) + bottom
    gen = gen+1
    print("\r%d/%d" % (gen,iteration),end="")
        
    x = sort_x
    #计算亲和度
    MSLL = np.zeros((1,num))
    for i in range(num):
        MSLL[0,i] = affinity(x[:,i])
    #按亲和度排序
    index = np.argsort(MSLL).flatten()
    x = x[:,index]
    #记录优化过程，此时第一个为做接近最值的结果
    trace.append(affinity(x[:,0]))
    
#作图
xx = range(gen-1)
plt.plot(xx , trace)
print(trace[-1])
plt.show()

例（2）

求函数$$f(x,y)=5sin(xy)+x^2+y^2$$的最小值，其中x的范围是[-4,4]，y的取值范围是[-4,4]，这是一个有多个局部极值的函数。把上例的亲和度函数更改，参数略作修改即可，详细注释参考上面。

#亲和度：与目标函数相关
def affinity(x):
    y = x[1]
    x = x[0]
    res = 5 * np.sin( x*y ) + x*x + y*y
    return res

#激励度：根据亲和度和浓度计算，影响变异
def motivation(x,num):
    ND = np.zeros((1,num)) #临时
    nd = np.zeros((1,num)) #浓度
    for i in range(num):
        for j in range(num):
            tmp = np.sqrt(np.sum((x[:,i]-x[:,j])**2))
            if(tmp > detas):
                ND[0,i] = 0
            else:
                ND[0,i] = 1                
        nd[0,i] = np.sum(ND,axis=1)/num
    ## 激励度
    MSLL = np.zeros((1,num))
    for i in range(num):
        tmp = affinity(x[:,i])
        MSLL[0,i] = tmp
    mtv = alphe * MSLL - beta * nd
    return mtv
#变异
def variation(x , gen):
    temp = num//2
    first_x = np.zeros((D,temp))
    for i in range(temp):
        tx=x[:,i].reshape((D,1))
        tx_rpt=np.tile(tx,10)
        randt=np.random.rand(1)
        if randt > pm:
            for g in range(D):
                for h in range(D):
                    tx_rpt[g,h] = tx_rpt[g,h] + (np.random.rand(1) - 0.5) * (defalut / (gen+1))
        for j in range(D):
            for k in range(D):
                if (tx_rpt[j,k] > top) or (tx_rpt[j,k] < bottom):
                    tx_rpt[j,k] = np.random.rand(1) * (top - bottom) + bottom
        tx_rpt[:,0] = tx[:,0]
        mvt = motivation(tx_rpt,10)
        index = np.argsort(mvt).flatten()
        temp = tx_rpt[:,index]
        first_x[:,i] = temp[:,0]
    return first_x

import numpy as np
import matplotlib.pyplot as plt
D = 2 #维度
num = 100 #个体数量
top = 4
bottom = -4
alphe = 3
beta = 1
detas = 0.3 #浓度阈值
pm = 0.6   #变异几率
trace = [] 
defalut = top
x = np.random.rand(D,num) * (top - bottom) + bottom  #10*100

gen = 1
trace = []
iteration = 80
while gen<iteration:
    
    mtv = motivation(x,num) 
    #激励度排序
    index = np.argsort(mtv).flatten()
    sort_x = x[:,index]
    sort_x[:,:50] = variation(sort_x[:,:50],gen)

    sort_x[:,50:] = np.random.rand(D,50) * (top - bottom) + bottom
    gen = gen+1

    print("\r%d/%d" % (gen,iteration),end="")
    x = sort_x    
    MSLL = np.zeros((1,num))
    for i in range(num):
        MSLL[0,i] = affinity(x[:,i])
    index = np.argsort(MSLL).flatten()
    x = x[:,index]
    trace.append(affinity(x[:,0]))
    
xx = range(gen-1)
plt.plot(xx , trace)
print("\n最值：")
print(trace[-1])
plt.show()

参考：
[1].刘韬.《人工免疫系统及其数据挖掘应用研究》.中国矿业大学出版社.
[2].https://zhuanlan.zhihu.com/p/112727537