BPSO离散二进制粒子群算法解决背包问题

痛定思痛。 2022-11-15 12:26 120阅读 0赞

## 算法原理 ##

粒子群算法原理 见链接  
[粒子群算法求解无约束优化问题 源码实现][Link 1]  
[粒子群算法求解带约束优化问题 源码实现][Link 2]

![20200702203245749.png_pic_center][] ![202103291336007.jpg_pic_center][]

离散粒子群  
      基本的粒子群算法是在连续域中搜索函数极值的有力工具，在此之后，提出了离散二进制粒子群，将离散问题空间映射到连续粒子运动空间，并适当修改粒子群算法来求解，在计算上仍保留经典粒子群算法速度-位置更新运算法则，粒子在状态空间的取值和变化只限于0和1两个值，而速度的每一维Vij代表位置每一位xij的取值为1的可能性，，因此在连续粒子群中的vij更新公式仍保持不变，但是Pbest和Gbest只在【0,1】的取值，其位置更新等式为：

![20210401141415703.png_pic_center][]

式中：r是0-1之间的随机数。

爱学习，爱玉酱。下面开始讲解算例。  
![20210520153651109.png_pic_center][]

## 背包问题 ##

背包问题。有N件物品和一个容量为V的背包。第i件物品的体积是c(i)，价值是w(i)。求解将哪些物品放入背包可使物品的体积总和不超过背包的容量，且价值总和最大。假设物品数量为10，背包的容量为300。每件物品的体积为\[95，75，23， 73，50， 22，6，57，89， 98\]， 价值为\[89， 59，19， 43,100，72，44，16，7, 64\]。

## python求解 ##

import numpy as np
    import matplotlib.pyplot as plt
    from tqdm import tqdm#进度条设置
    import matplotlib.pyplot as plt
    import matplotlib; matplotlib.use('TkAgg')
    from pylab import *
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
    
    class BackPack(object):
        def __init__(self):
            self.Weight = [95,75,23,73,50,22,6,57,89,98]  # 物品体积
            self.Value = [89, 59, 19, 43, 100, 72, 44, 16, 7, 64]  # 物品价值
            self.N = 100  # 群体粒子个数
            self.D = len(self.Weight)  # 粒子维数
            self.T = 100  # 最大迭代次数
            self.c1 = 1.5  # 学习因子1
            self.c2 = 1.5  # 学习因子2
            self.w = 1  # 惯性因子，一般取1
            self.V_max = 10  # 速度最大值
            self.V_min = -10  # 速度最小值
            self.Weight_max = 300  # 背包容量
            self.afa = 10  # 惩罚系数
            self.G=100 #迭代次数
    
        #初始化种群
        def init_x(self):
            """ :return: 随机生成的种群(二维list) """
            X = np.random.choice([0, 1], size=(self.N, self.D))
            return X
        #初始化速度
        def init_v(self):
            """ :return: 速度 """
            V = np.random.random(size=(self.N, self.D))  ##10维
            return V
    
    
        # 适应度值
        def fitness_func(self,X):
            """ :param X: 种群 :return: """
            """计算粒子的的适应度值，也就是目标函数值，X 的维度是 size * 10 """
            size = np.shape(X)[0] #种群个数
            result = np.zeros((size, 1))
            for i in range(size):  # 遍历每一个粒子
                fit = sum(X[i] * self.Value)  # 物品价值
                TotalSize = sum(X[i] * self.Weight)  # 物品重量
                if TotalSize <= self.Weight_max:  # 如果小于限制重量
                    fit = fit
                    result[i] = fit
                else:
                    #fit = fit - self.afa * (TotalSize - self.Weight_max) # 对fit进行惩罚
                    fit=0
                    result[i] = fit
            return result  # 我们要求result越大
    
        # 速度更新公式
        def velocity_update(self,V, X, pbest, gbest):
            """ 根据速度更新公式更新每个粒子的速度 种群size=100 :param V: 粒子当前的速度矩阵，size*10 的矩阵 :param X: 粒子当前的位置矩阵，size*10 的矩阵 :param pbest: 每个粒子历史最优位置，size*10 的矩阵 :param gbest: 种群历史最优位置，1*10 的矩阵 """
            r1 = np.random.random((self.N, 1)) #(粒子个数,1)
            r2 = np.random.random((self.N, 1))
            V = self.w * V + self.c1 * r1 * (pbest - X) + self.c2 * r2 * (gbest - X)  # 直接对照公式写就好了
            # 防止越界处理
            V[V < self.V_min] = np.random.random() * (self.V_max - self.V_min) + self.V_min
            V[V > self.V_max] = np.random.random() * (self.V_max - self.V_min) +self.V_min
            return V
    
        # 位置更新公式
        def position_update(self,X, V):
            """ 根据公式更新粒子的位置 :param X: 粒子当前的位置矩阵，维度是 size*10 :param V: 粒子当前的速度举着，维度是 size*10 return 更新后的种群 """
    
            for i in range(self.N):  # 遍历每一个粒子
                # 修改速度为sigmoid形式
    
                V[i, :] = 1. / (1 + np.exp(-np.array(V[i, :])))
                for j in range(self.D):  # 遍历粒子中的每一个元素
                    rand = np.random.random()  # 生成 0-1之间的随机数
                    if V[i, j] > rand:
                        X[i, j] = 1
                    else:
                        X[i, j] = 0
                #对当前个体进行限制，放在超重
                while np.sum(X[i]*self.Weight)>self.Weight_max:#如果当前粒子超重
                    X[i]=np.random.choice([0, 1], size=(1, self.D))
            return X
    
        def update_pbest(self,X, fitness, pbest, pbestfitness, m):
            """ 更新个体最优 :param X: 当前种群 :param fitness: 当前每个粒子的适应度 :param pbest: 更新前的个体最优解 :param pbestfitness: 更新前的个体最优适应度 :param m: 粒子数量 :return: 更新后的个体最优解、个体最优适应度 """
    
            for i in range(m):
                if fitness[i] > pbestfitness[i]:
                    pbest[i] = X[i]
                    pbestfitness[i] = fitness[i]
            return pbest, pbestfitness
    
        def update_gbest(self,pbest, pbestfitness, gbest, gbestfitness, m):
            """ 更新全局最优解 :param pbest: 粒子群 :param pbestfitness: 个体适应度(个体最优) :param gbest: 全局最优解 :param gbestfitness: 全局最优解 :param m: 粒子数量 :return: gbest全局最优，g对应的全局最优解 """
            for i in range(m):
                if pbestfitness[i] > gbestfitness:
                    gbest = pbest[i]
                    gbestfitness = pbestfitness[i]
            return gbest, gbestfitness
    
    
        def main(self):
            fitneess_value_list = []  # 记录每次迭代过程中的种群适应度值变化
            x = self.init_x()  # 初始化x
            v = self.init_v()  # 初始化v
            # 计算种群各个粒子的初始适应度值
            p_fitness = self.fitness_func(x)
            # 计算种群的初始最优适应度值
            g_fitness = p_fitness.max()
            # 讲添加到记录中
            fitneess_value_list.append(g_fitness)
            # 初始的个体最优位置和种群最优位置
            pbest = x
            gbest = x[p_fitness.argmax()]  #
    
            # 接下来就是不断迭代了
            for i in tqdm(range(self.G)):
                pbest=pbest.copy()#必须加，不然会被篡改值，造成结果错
                p_fitness= p_fitness.copy()#必须加，不然会被篡改值，造成结果错
                gbest=gbest.copy()#必须加，不然会被篡改值，造成结果错
                g_fitness=g_fitness.copy()#必须加，不然会被篡改值，造成结果错
                v = self.velocity_update(v, x, pbest, gbest)  # 更新速度
                x = self.position_update(x, v)  # 更新位置
                p_fitness2 = self.fitness_func(x)  # 计算子代各个粒子的适应度
    
                # 更新每个粒子的历史最优位置
                pbest, p_fitness=self.update_pbest(x, p_fitness2, pbest, p_fitness, self.N)
    
                #更新群体的最优位置
                gbest, g_fitness=self.update_gbest(pbest, p_fitness, gbest, g_fitness, self.N)
    
                # 记录最优迭代结果
                fitneess_value_list.append(g_fitness)
    
            print("最优适应度是：%.5f" % fitneess_value_list[-1])
            print("最优解是", gbest)
            print('最优解对应的体积',np.sum(gbest*self.Weight))
            print('最优解对应的价值',np.sum(gbest*self.Value))
    
            plt.plot(fitneess_value_list,label='迭代曲线')
            plt.xlabel('迭代次数')
            plt.ylabel('适应度')
            plt.legend()
            plt.show()
    
    if __name__=='__main__':
        pso = BackPack()
        pso.main()

![20210520153448841.png_pic_center][] ![20210520153546115.png_pic_center][]

## MATLAB求解 ##

%%%%%%%%%%%%%%%%%离散粒子群算法解决0-1背包问题%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear all;              %清除所有变量
    close all;              %清图
    clc;                    %清屏
    N=100;                  %群体粒子个数
    D=10;                   %粒子维数
    c1=1.5;                 %学习因子1
    c2=1.5;                 %学习因子2
    Wmax=0.8;               %惯性权重最大值
    Wmin=0.4;               %惯性权重最小值
    Vmax=10;                %速度最大值
    Vmin=-10;               %速度最小值
    V = 300;                             %背包容量
    C = [95,75,23,73,50,22,6,57,89,98];  %物品体积
    W = [89,59,19,43,100,72,44,16,7,64]; %物品价值
    afa = 2;                             %惩罚函数系数
    %%%%%%%%%%%%%%%%初始化种群个体（限定位置和速度）%%%%%%%%%%%%%%%%
    x=randint(N,D);        %随机获得二进制编码的初始种群
    v=rand(N,D) * (Vmax-Vmin)+Vmin;
    %%%%%%%%%%%%%%%%%%初始化个体最优位置和最优值%%%%%%%%%%%%%%%%%%%
    p=x;
    pbest=ones(N,1);
    for i=1:N
        pbest(i)= func(x(i,:),C,W,V,afa);
    end
    %%%%%%%%%%%%%%%%%%%初始化全局最优位置和最优值%%%%%%%%%%%%%%%%%%
    g=ones(1,D);
    gbest=-10000;
    for i=1:N
        if(pbest(i)>gbest)
            g=p(i,:);
            gbest=pbest(i);
        end
    end
    gb=ones(1,T);
    %%%%%%%%%%%按照公式依次迭代直到满足精度或者迭代次数%%%%%%%%%%%%%
    while 1>0 %无限迭代
        for j=1:N
            %%%%%%%%%%%%%%更新个体最优位置和最优值%%%%%%%%%%%%%%%%%
            if (func4(x(j,:),C,W,V,afa)>pbest(j))
                p(j,:)=x(j,:);
                pbest(j)=func(x(j,:),C,W,V,afa);
            end
            %%%%%%%%%%%%%%%%更新全局最优位置和最优值%%%%%%%%%%%%%%%
            if(pbest(j)>gbest)
                g=p(j,:);
                gbest=pbest(j);
            end
            %%%%%%%%%%%%%%%%计算动态惯性权重值%%%%%%%%%%%%%%%%%%%%
            w=Wmax-(Wmax-Wmin)*i/T;
            %%%%%%%%%%%%%%%%%跟新位置和速度值%%%%%%%%%%%%%%%%%%%%%
            v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
                +c2*rand*(g-x(j,:));
            %%%%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%
            for ii=1:D
                if (v(j,ii)>Vmax)  |  (v(j,ii)< Vmin)
                    v(j,ii)=rand * (Vmax-Vmin)+Vmin;
                end
            end
            vx(j,:)=1./(1+exp(-v(j,:)));
            for jj=1:D
                if vx(j,jj)>rand
                    x(j,jj)=1;
                else
                    x(j,jj)=0;
                end
            end
        end
        %%%%%%%%%%%%%%%%%%%%记录历代全局最优值%%%%%%%%%%%%%%%%%%%%%
        gb(i)=gbest;
        if g.*C<=300 %如果满足条件 退出
            break
    
    
    end
    g;                       %最优个体
    figure
    plot(gb)
    xlabel('迭代次数');
    ylabel('适应度值');
    title('适应度进化曲线')
    
    
    %%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%
    function result = func(f,C,W,V,afa)
    fit = sum(f.*W);
    TotalSize = sum(f.*C);
    if TotalSize <= V
        fit = fit;
    else
        fit = fit - afa * (TotalSize - V);
    end
    result = fit;

![在这里插入图片描述][20210401192449829.jpg_pic_center]  
作者：电气 余登武  
![20210401192555481.jpg_pic_center][]

[Link 1]: https://blog.csdn.net/KOBEYU652453/article/details/115298847
[Link 2]: https://blog.csdn.net/KOBEYU652453/article/details/115300479
[20200702203245749.png_pic_center]: https://img-blog.csdnimg.cn/20200702203245749.png#pic_center
[202103291336007.jpg_pic_center]: /images/20221022/d749fe48d7df4c02a900546f2a4337db.png
[20210401141415703.png_pic_center]: /images/20221022/b6d288c22c4646229209edda1f9da3bc.png
[20210520153651109.png_pic_center]: /images/20221022/1dab5f8cc97a4f22bded05a4ca66cbf4.png
[20210520153448841.png_pic_center]: https://img-blog.csdnimg.cn/20210520153448841.png#pic_center
[20210520153546115.png_pic_center]: /images/20221022/cc90eca86fff4d80b10f567ee93263aa.png
[20210401192449829.jpg_pic_center]: /images/20221022/9737b4ac1eed4f248232497b6960e8ed.png
[20210401192555481.jpg_pic_center]: /images/20221022/64120b9c1d614695874d58ebdb8f99ea.png