粒子群算法求解多元函数最值问题

粒子群算法求解多元函数最值问题

一、简介

多元函数极值&最值问题通常使用导数/偏导数进行推导，这里尝试使用启发式算法进行求解近似最优解。选用粒子群算法进行求解，粒子群算法模仿鸟群觅食行为，核心思想是通过向距离食物最近的鸟集聚，不断更新速度和位置以达到最优解，即表现不好的个体通过向表现好的个体学习使得自身往好的方向转变，这里存在一个前提：所有鸟知道距离食物的远近，距离食物最近包含两部分：当前最近和历史最近。标准粒子群算法适合求解函数极值问题，在TSP、背包问题上多用混合型粒子群算法。算法详细介绍可参考【粒子群算法研究】。

二、粒子群算法求解

多元函数：y = ($x_1$+$x_2$)/($x_1^2$+$x_2^2$+1)

粒子群编程求解

import pylab as mpl 
import random import
matplotlib.pyplot as plt 
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 添加这条可以让图形显示中文

def fitness(x):
    """
    计算适应值：y = (x1+x2)/(x1**2+x2**2+1)
    """
    x1 = x[0]
    x2 = x[1]
    y = (x1+x2)/(x1**2+x2**2+1)
    return y

def update(x, p, g, w, c1, c2, v_low, v_high, vi, d, x_low, x_up):
    '''
    更新速度和位置
    '''
    # 更新速度
    vi = w * vi + c1 * random.uniform(0, 1) * (p-x) + c2 * random.uniform(0, 1) * (g-x)
    for j in range(d):
        if vi[j] >v_high:
            vi[j] = v_high
        elif vi[j] < v_low:
            vi[j] = v_low

    # 更新位置
    x = x +vi
    for j in range(d):
        if x[j] > x_up[j]:
            x[j] = x_up[j]
        elif x[j] < x_low[j]:
            x[j] = x_low[j]
    return x


if __name__ == '__main__':
    # 参数
    birdNum = 20  # 粒子数量
    w = 0.2  # 惯性因子
    c1 = 0.4  # 自我认知因子
    c2 = 0.4  # 社会认知因子
    # pBest, pLine = 0, []  # 当前最优值、当前最优解，（自我认知部分）
    # gBest, gLine = 0, []  # 全局最优值、全局最优解，（社会认知部分）
    v_low = -2 #速度约束
    v_high =2
    iterMax = 100  # 迭代次数
    iterI = 1  # 当前迭代次数
    bestfit = -100000  # 记录每代最优值
    #y = (x1+x2)/(x1**2+x2**2+1)
    d = 2 #函数变量个数
    x_low = [-10, -10]#变量范围
    x_up = [10, 10]

    #随机生成初始解
    xs = np.array([np.array([random.uniform(x_low[0],x_up[0]),random.uniform(x_low[1],x_up[1])])
	for i in range(birdNum)])
    v = np.array([random.uniform(v_low,v_high) for i in range(birdNum)])
    fits = np.array([fitness(x) for x in xs])

    fit = max(fits)
    x = xs[np.argmax(fits)]

    gBest = pBest = fit # 全局最优值、当前最优值
    gSolve = pSolve = x  # 全局最优解、当前最优解
    bestfit = []
    bestfit.append(gBest)

    while iterI <= iterMax:  # 迭代开始
        for i in range(birdNum):
            xs[i] = update(xs[i], pSolve, gSolve, w, c1, c2, v_low, v_high, v[i], d, x_low, x_up)
            fits[i] = fitness(xs[i])

        pBest, pSolve = max(fits), xs[np.argmax(fits)]

        if max(fits) >= gBest:
            gBest, gSolve = max(fits), xs[np.argmax(fits)]

        bestfit.append(gBest)
        print(iterI, gBest)  # 打印当前代数和最佳适应度值
        print(gSolve, pSolve)
        iterI += 1  # 迭代计数加一

    #迭代图
    plt.plot(range(1,len(bestfit)+1), bestfit)
    plt.xlabel("迭代次数")
    plt.ylabel("最优值")
    plt.show()

结果
$y$=0.7068,$x_1=0.8701,x_2=0.8287$