以Ras函數（Rastrigin's Function）為目標函數，求其在x1,x2∈[-5,5]上的最小值。這個函數對模擬退火、進化計算等算法具有很強的欺騙性，因為它有非常多的局部最小值點和局部最大值點，很容易使算法陷入局部最優，而不能得到全局最優解。如下圖所示，該函數只在(0,0)處存在全局最小值0。

Python代碼實現

import numpy as np
import matplotlib.pyplot as plt




# 目標函數定義
def ras(x):
    y = 20 + x[0] ** 2 + x[1] ** 2 - 10 * (np.cos(2 * np.pi * x[0]) + np.cos(2 * np.pi * x[1]))
    return y




# 參數初始化
w = 1.0
c1 = 1.49445
c2 = 1.49445


maxgen = 200  # 進化次數
sizepop = 20  # 種群規模


# 粒子速度和位置的范圍
Vmax = 1
Vmin = -1
popmax = 5
popmin = -5


# 產生初始粒子和速度
pop = 5 * np.random.uniform(-1, 1, (2, sizepop))
v = np.random.uniform(-1, 1, (2, sizepop))


fitness = ras(pop)  # 計算適應度
i = np.argmin(fitness)  # 找最好的個體
gbest = pop  # 記錄個體最優位置
zbest = pop[:, i]  # 記錄群體最優位置
fitnessgbest = fitness  # 個體最佳適應度值
fitnesszbest = fitness[i]  # 全局最佳適應度值


# 迭代尋優
t = 0
record = np.zeros(maxgen)
while t < maxgen:


    # 速度更新
    v = w * v + c1 * np.random.random() * (gbest - pop) + c2 * np.random.random() * (zbest.reshape(2, 1) - pop)
    v[v > Vmax] = Vmax  # 限制速度
    v[v < Vmin] = Vmin


    # 位置更新
    pop = pop + 0.5 * v
    pop[pop > popmax] = popmax  # 限制位置
    pop[pop < popmin] = popmin


    '''
    # 自適應變異
    p = np.random.random()             # 隨機生成一個0~1內的數
    if p > 0.8:                          # 如果這個數落在變異概率區間內，則進行變異處理
        k = np.random.randint(0,2)     # 在[0,2)之間隨機選一個整數
        pop[:,k] = np.random.random()  # 在選定的位置進行變異 
    '''


    # 計算適應度值
    fitness = ras(pop)


    # 個體最優位置更新
    index = fitness < fitnessgbest
    fitnessgbest[index] = fitness[index]
    gbest[:, index] = pop[:, index]


    # 群體最優更新
    j = np.argmin(fitness)
    if fitness[j] < fitnesszbest:
        zbest = pop[:, j]
        fitnesszbest = fitness[j]


    record[t] = fitnesszbest  # 記錄群體最優位置的變化


    t = t + 1


# 結果分析
print(zbest)


plt.plot(record, 'b-')
plt.xlabel('generation')
plt.ylabel('fitness')
plt.title('fitness curve')
plt.show()

結果為

[0.99699579 0.00148844]

可以知道求解的點非最小值，算法陷入了局部最小值。

刪除自適應變異部分的注釋，運行后結果如下，可以看出收斂到全局最優解。

[0.00022989 0.00014612]

Matlab有個自帶的粒子群優化函數particleswarm也可以使用。本例的代碼如下：

y = @(x) 20 + x(1).^2 + x(2).^2 - 10*(cos(2*pi*x(1))+cos(2*pi*x(2)));
rng default
options = optimoptions('particleswarm','SwarmSize',200,'HybridFcn',@fmincon,'MaxIterations',200, 'Display','iter');
lb = [-5 -5];     % 這是變量的下限
ub = [5 5];       % 這是變量的上限
[x,fval,exitflag,output] = particleswarm(y,length(lb),lb,ub,options);

結果如下