526互联

粒子群算法(Particle Swarm Optimization, PSO)

发布时间 2023-09-17 22:34:13作者: 自倚修行

Particle Swarm Optimization

算法原理参考: https://zhuanlan.zhihu.com/p/404198434

Question

使用PSO算法计算函数$ f(x) = x_1^2 + 3 x_2^2 - x_1 + 2 x_2 - 5 $ 在 \(x \in [-100, 100]^2\) 上的最小值

Code

import numpy as np
from typing import List

Parameters

dim = 2      # varible dimension
N = 200       # number of partials 
max_v = 20   # the upper bound of velocity
max_iter = 50 # max iter
boundary = np.array([[-100, 100], [-100, 100]]) # [[dimK Lower Bound, dimK Upper Bound]]

c1 = 0.9     # inertia weight
c2 = 2       # neighborhood accelerate coefficient
c3 = 2       # global best accelerate coefficient

param = {
    'c1' : c1,
    'c2' : c2,
    'c3' : c3,
    'boundary': boundary
}


## compute the f
def f(t: List[float]) -> float:
    
    # x = t[0]
    # y = t[1]
    # z =((1*np.cos((1+1)*x+1))+(2*np.cos((2+1)*x+2))+(3*np.cos((3+1)*x+3))+ \
    #     (4*np.cos((4+1)*x+4))+(5*np.cos((5+1)*x+5)))*((1*np.cos((1+1)*y+1))+ \
    #     (2*np.cos((2+1)*y+2))+(3*np.cos((3+1)*y+3))+(4*np.cos((4+1)*y+4))+(5*np.cos((5+1)*y+5)))

    z = t[0]**2 + 3 * t[1]**2 - t[0] + 2 * t[1] - 5
    return z

PSO algorithm

## limit position
def limit_p(P, boundary):
    N = P.shape[0]
    # for each dimension
    def fun(a):
        if a < boundary[i][0]:
            return boundary[i][0]
        elif a > boundary[i][1]:
            return boundary[i][1]
        else:
            return a
    

    for i in range(boundary.shape[0]):
        P[:, i] = np.array(list(map(fun, P[:, i])))

    return P


## limit velocity
def limit_v(V, boundary):
    N = V.shape[0]
    # for each dimension
    for i in range(boundary.shape[0]):
        # the max velocity may not proceed 10% boundary size
        lim = 0.1 * (boundary[i][1] - boundary[i][0])
        V[:, i] = np.array(list(map(min, [lim] * N, V[:, i])))
    return V


## update the inertia weight
def update_param(param: dict, max_iter: int) -> dict:
    param['c1'] = 0.9 - 0.4 / max_iter * param['c1']
    return param


## update position and velocity
def update_p_v(P, V, pBest, gBest, param):
    V = param['c1'] * V + \
        param['c2'] * (pBest - P) * np.random.random() + \
        param['c3'] * (gBest - P) * np.random.random()
    V = limit_v(V, param['boundary'])
    P = limit_p(P + V, param['boundary'])
    return P, V


## compute nBest gBest
def compute_pB_gB(P, f, pBest, gBest):
    y = np.array(list(map(f, P)))
    for i in range(len(y)):
        if f(pBest[i]) > y[i]:
            pBest[i] = P[i]
    if min(y) < f(gBest):
        gBest = P[np.argmin(y)]
    return pBest, gBest

P = np.random.random(size=[N, dim]) * 200 - 100     # initial positions
V = np.random.random(size=[N, dim]) * 20 - 10       # initial velocties
pBest = P                                           # initial personal best
gBest = P[0]                                        # initial the fake global best

## iteration
# gBest_his = []
for i in range(max_iter):
    pBest, gBest = compute_pB_gB(P, f, pBest, gBest)
    P, V = update_p_v(P, V, pBest, gBest, param)
    param = update_param(param, max_iter)

    print("Iter {}, Best value = {:.3f}".format(i, f(gBest)))
    # gBest_his.append(f(gBest))

print("The min_f reached at", gBest)
# from matplotlib import pyplot as plt
# plt.plot(gBest_his)
# plt.show()

Iter 0, Best value = 137.194
Iter 1, Best value = -5.293
Iter 2, Best value = -5.293
Iter 3, Best value = -5.293
Iter 4, Best value = -5.293
Iter 5, Best value = -5.293
Iter 6, Best value = -5.293
Iter 7, Best value = -5.293
Iter 8, Best value = -5.293
Iter 9, Best value = -5.293
Iter 10, Best value = -5.293
Iter 11, Best value = -5.293
Iter 12, Best value = -5.293
Iter 13, Best value = -5.293
Iter 14, Best value = -5.336
Iter 15, Best value = -5.336
Iter 16, Best value = -5.336
Iter 17, Best value = -5.336
Iter 18, Best value = -5.336
Iter 19, Best value = -5.439
Iter 20, Best value = -5.439
Iter 21, Best value = -5.439
Iter 22, Best value = -5.500
Iter 23, Best value = -5.500
Iter 24, Best value = -5.500
Iter 25, Best value = -5.500
Iter 26, Best value = -5.581
Iter 27, Best value = -5.581
Iter 28, Best value = -5.581
Iter 29, Best value = -5.581
Iter 30, Best value = -5.581
Iter 31, Best value = -5.581
Iter 32, Best value = -5.581
Iter 33, Best value = -5.581
Iter 34, Best value = -5.581
Iter 35, Best value = -5.581
Iter 36, Best value = -5.581
Iter 37, Best value = -5.581
Iter 38, Best value = -5.581
Iter 39, Best value = -5.582
Iter 40, Best value = -5.582
Iter 41, Best value = -5.582
Iter 42, Best value = -5.582
Iter 43, Best value = -5.582
Iter 44, Best value = -5.582
Iter 45, Best value = -5.582
Iter 46, Best value = -5.582
Iter 47, Best value = -5.582
Iter 48, Best value = -5.583
Iter 49, Best value = -5.583
The min_f reached at [ 0.4941146  -0.33489041]