[演算法] 粒子群演算法概念 Particle Swarm Optimization, PSO

在看過螞蟻演算法、基因演算法後，本篇應該是筆者介紹演算法中的最後一篇，有鑑於並非長期耕耘在這個領域，故大多數的介紹，讀者可以當成聽聽故事，或許內容部分不見得完全正確，不過筆者就是盡一份心力，試著能讓欲嘗試這領域的新手快速理解，建立一個簡單的概念即可。

同前面的演算法，本篇也是請人工智慧快速整理出一份代碼，而我們直接針對代碼部分做進一步的了解。

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

# 定義函數 f(x, y) = x^2 + y^2
def f(x, y):
    return x**2 + y**2

# PSO 之重要參數
num_particles = 30                  # 粒子數量
num_dimensions = 2                  # 搜索空間的維度
max_iterations = 100                # 最大疊代次數
inertia = 0.5                       # 慣性權重
c1 = 1.5                            # 個體最優影響因子
c2 = 1.5                            # 全局最優影響因子

# 規範搜索空間範圍
x_min, x_max = -10, 10
y_min, y_max = -10, 10

# 初始化粒子群
particles = []
for _ in range(num_particles):
    position = np.array([random.uniform(x_min, x_max), random.uniform(y_min, y_max)])
    velocity = np.array([random.uniform(-1, 1), random.uniform(-1, 1)])
    best_position = np.copy(position)
    best_fitness = f(position[0], position[1])
    particles.append({
        'position': position,
        'velocity': velocity,
        'best_position': best_position,
        'best_fitness': best_fitness
    })

# 初始化全局最優
global_best_position = particles[0]['best_position']
global_best_fitness = particles[0]['best_fitness']

# 主循環
for iteration in range(max_iterations):
    for particle in particles:
        # 更新速度
        r1 = np.random.random(num_dimensions)  # 隨機數
        r2 = np.random.random(num_dimensions)
        inertia_component = inertia * particle['velocity']
        cognitive_component = c1 * r1 * (particle['best_position'] - particle['position'])
        social_component = c2 * r2 * (global_best_position - particle['position'])
        particle['velocity'] = inertia_component + cognitive_component + social_component
        
        # 更新位置
        particle['position'] += particle['velocity']
        
        # 保持在搜索空間範圍內
        particle['position'][0] = np.clip(particle['position'][0], x_min, x_max)
        particle['position'][1] = np.clip(particle['position'][1], y_min, y_max)
        
        # 計算適應度
        current_fitness = f(particle['position'][0], particle['position'][1])
        
        # 更新個體最優
        if current_fitness < particle['best_fitness']:
            particle['best_fitness'] = current_fitness
            particle['best_position'] = np.copy(particle['position'])
        
        # 更新全局最優
        if current_fitness < global_best_fitness:
            global_best_fitness = current_fitness
            global_best_position = np.copy(particle['best_position'])

print(f"全局最優位置：{global_best_position}, 全局最優適應度：{global_best_fitness}")

# 繪製粒子位置
x = [p['position'][0] for p in particles]
y = [p['position'][1] for p in particles]
plt.scatter(x, y, color='blue', label='Particles')
plt.scatter(global_best_position[0], global_best_position[1], color='red', label='Global Best')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Particle Swarm Optimization')
plt.legend()
plt.show()

在粒子群演算法中，有分成所謂的個體最優和全體最優概念。

起初，如果是第一次接觸到這個概念或許會覺得有點特別，為什麼要分個體最優和全體最優，兩者又有何差異。

我們回顧一下螞蟻演算法和基因演算法，在螞蟻演算法裡面，我們給定初始的螞蟻數量是 10，在基因演算法中，儘管我們假定一個人染色體上的基因數量是 5，但重點在於我們在每一代每一代的人數也都設定為 10 人，當作觀察目標。而在粒子群演算法之中，給定最初始的觀察粒子數量是 30。

回顧螞蟻演算法疊代過程中，我們不會在意螞蟻現在的位置是在哪裡，因為每一隻螞蟻的重點是去完成整趟旅程，並且紀錄下這趟旅程的總距離，最終檢視是否為最佳解，也就是所謂的全局最佳。

基因演算法中，我們可以看到從一開始給定 10 位父代染色體，並且透過一代一代的遺傳，試圖想從基因的觀點找尋到最佳的染色體特徵，同樣的我們想要找到的也僅只於全局的最佳，並沒有所謂的個體最佳的概念。

那麼，究竟何謂個體最佳呢？

粒子群演算法與前述兩個方法不同的之處在於，我們給定初始的 30 個粒子，這 30 個粒子個別的歷史移動軌跡都需要紀錄保留下來。也因為每個移動軌跡都被記錄，可以很容易想像：以第 1 個粒子來說，他可能總計移動 100 次，那麼在這個 100 次裡面，每次的移動不僅僅自己和過往的自己比較，同時也可以視其他粒子的當前過往的歷史軌跡，給予一個最適合的移動方向。而自己和過往的自己比較的這件事情就叫做個體最佳。

也就是說，與螞蟻演算法的差異是，這 30 個粒子有點像是一個團隊去搜索地震的倖存者，我們在意的不是 30 個粒子過往的移動路線是否相似。我們的目的就是讓每個人都分散式去搜尋各種可能，意思是，非常有可能倖存者不是僅僅一位，這對於每個人在搜尋時，他自己過往的移動路線上一樣的重要，不能僅僅朝向已經搜索到倖存者的那位隊員前進。因為有可能在你身邊，也同樣存在著倖存者。因此，自己過往的資訊也就非常關鍵，需要被納入考慮。

與基因演算法不同的是，基因演算法適合組合的問題，也就是說，我們可以發現，粒子群演算法適合處理的是連續型的問題，讀者可以自行體會兩者之間的差異。

簡單來說，可以把演算法歸納看成連續與不連續型的用法根本差異。而與螞蟻演算法, 基因演算法不同，例子群演算法特別適合處理這種連續性的問題。

r1 = np.random.random(num_dimensions)  # 隨機數
r2 = np.random.random(num_dimensions)
inertia_component = inertia * particle['velocity']
cognitive_component = c1 * r1 * (particle['best_position'] - particle['position'])
social_component = c2 * r2 * (global_best_position - particle['position'])
particle['velocity'] = inertia_component + cognitive_component + social_component

以上這幾行可以說是本演算法的核心，我們想問的問題是，當前個粒子，究竟應該如何移動？移動到哪？

本次代碼我們探討的是一個二維度的問題，所以你會看到我們以下的向量維度多數為二：

>>> r1 = np.random.random(num_dimensions)
>>> r2 = np.random.random(num_dimensions)
>>> inertia_component = inertia * particle[‘velocity’]
>>> cognitive_component = c1 * r1 * (particle[‘best_position’] — particle[‘position’])
>>> social_component = c2 * r2 * (global_best_position — particle[‘position’])
>>> particle[‘velocity’] = inertia_component + cognitive_component + social_component
>>> r1
array([0.5813269 , 0.45378116])
>>> r2
array([0.73456962, 0.2201673 ])
>>> inertia_component
array([-5.40754058e-13, -4.80791677e-12])
>>> cognitive_component
array([0., 0.])
>>> social_component
array([1.88226192e-11, 5.33599185e-13])
>>> particle[‘velocity’]
array([ 1.82818652e-11, -4.27431759e-12])
>>> inertia # 慣性權重 0.5
>>> c1 # 個體最優影響因子1.5
>>> c2 # 全局最優影響因子1.5

在以上的代碼中，由 cognitive_component = array([0., 0.])，我們可以看到當前粒子的個體最優位置即在目前的位置。然而，並不是意味著就不搜尋了，意思是，並非只有這個因子完全決定粒子的前進。除此之外，我們還有粒子本身的慣性 inertia_component，指的是由上個時間點位置到當前位置的速度，我們也採信。除此之外，還剩下一個重要因子就是全局最優的移動 social_component，也就是說，我們同樣也會傾向相信要往目前 30 個粒子移動過程中，大家探索到的那個最佳解位置方向移動。