《MATLAB智能优化算法》学习笔记一

一、0-1背包问题概述

有$n$件物品和一个载重量为$C$的背包。第$i$件物品的质量是$w_i$，其价值是$v_i$。求解在不超过背包载重量的情况下，将哪些物品装入背包可使价值总和最大。

因此，0-1背包问题的数学模型如下：

式中，$x_i$为0-1决策变量，表示物品$i$是否被装包，如果是，则$x_i=1$，否则，$x_i=0$。

接下来以一个实例讲解上述0-1背包问题。假设现有1个背包和5个物品，背包的最大载重量为$6kg$，这5个物品的质量和价值如表所示：

枚举所有可行的装包方案：

从上表可以看出方案8的装包物品总价值是最大的。

二、算法设计

这是$n$为5的情况下，可以把所有情况列举出来，但是当$n$稍微大一点时，装包的方案就很多，一一列举出来很费力。此时枚举法已经不适合用来解决此问题了，因此我将会使用智能优化算法中的遗传算法（GA）来求解该问题。首先要给位明确的一点是，使用GA求解0-1背包问题并不能保证找到最有的装包方案。

（1）编码

编码采用01编码的形式。如选择装入物品2、3、5对应的编码就是“01101”，这个染色体（个体）就为 [0 1 1 0 1]。

（2）解码

与编码恰好相反。染色体（个体）[1 0 1 0 0] 对应的是装入背包的物品是1号和3号。

（3）约束处理

约束处理的步骤如下：

(1）将已经装进包中的物品按照性价比（性价比＝价值／质量）由低到高进行排序。

(2）按照第（1）步排序结果，取走排在第1位的物品后，检验此时是否满足背包的载重量约束。如果满足约束，则将染色体中基因位上的数字1改为数字0．此时染色体初步修复完毕；如果可以不满足约束，则先将染色体中基因位上的数字1改为数字0．然后继续取走排在第2位的物品，并再次检验此时的染色体是否满足背句的载重量约束。循环往复，一直到满足背包的载重量约束为止．染色体初步修复完毕。

(3）在第（2）步已经得到满足背包载重量约束的染色体，但此时背包可能还有剩余空间。因此，将此时未装包的物品按照性价比从高到低排序，然后按照该顺序依次将物品装进包中。在装包过程中，将不满足约束的物品不装包，将满足约束的物品装包，并将染色体中基因位上的数字0改为数字1。一直遍历到最后一个未装包的物品为止，染色体最终修复完毕。

（4）适应度函数

只有对一条染色体进行合适的评价,才能保证GA在搜索过程中方向不”跑偏”。

适应度函数为物品价值之和，即

适应度值越大的染色体越优，反之越劣质。

（5）种群初始化

假设种群的数目为$NIND$，那么需要初始化$NIND$个个体（染色体）。因此只需要理解如何初始化一个个体就可以以此类推了。

假设有$n$个物品，那么一个染色体的长度为$n$，即这个染色体有$n$个数字组成（每个数字或是1或是0）。

初始化一个个体步骤如下：

①随机生成$n$个数字（每个数字或是1或是0），将此时的个体命名为Individual。

②检验Individual是否满足背包的载重约束，如果满足约束，个体Individual初始化为完毕；如果不满足约束，则对个体Individual进行约束处理，约束处理结束后，个体Individual初始化完毕。

按照一个个体初始化的方法，对$NIND$个个体全部进行初始化，初始化结束后，即完成对父代种群Chrom的初始化。

（6）选择操作

初始化种群后，种群中每个个体的适应度值都会存在差异。按道理我们应该选择适应度值大的个体来进行后续操作，但其实这样做很可能会导致种群在后续的进化中停滞不前，陷入局部最优。

所以我们不仅要多选择适应度值大的个体，还要兼顾适应度值小的个体，防止陷入局部最优。

在这里我们采用轮盘赌策略。

由于在自然界中动物也不是百分百孕育出后代，所以我们在选择个体时也不必要选择$NIND$个，可能选择出$Nsel=n\times GGAP$个（$GGAP$称为代沟，是一个大于0小于等于1的随机数）。因此，我们只要转动$Nsel$次转盘就能选出子代种群SelCh。

在轮盘赌转盘上，每个个体对应一个被选中的概率。假设第$i$个个体的适应度值为$Fitness_i$，那么其被选中的概率为：

假如有4个个体，每个个体被选中的概率分别为：25%，40%，21%，14%。轮盘赌转盘如图所示：

（7）交叉操作

选择出子代种群后，我们要种群向着适应度值增大的方向进化，所以需要改变个体上的基因，因此我们的第一个操作就是交叉操作。我们看下面这个例子：

个体1 ： 1 0 0 1 1 0 1 1

个体2 ： 0 1 0 0 0 1 0 1

随机选择两个交叉位置，$a_1=3,a_2=5$,那么两个位置中间那段染色体就会进行交换。

交换后：

个体1 ： 1 0 0 0 0 0 1 1

个体2 ： 0 1 0 1 1 1 0 1

由于自然界中交叉是概率性事件且交叉概率较大，所以我们在编写代码时也要给交叉操作一个概率并且取值较大。

当种群中所有个体都需要交叉操作时，我们就会将种群中$Nsel$个个体分成$Nsel/2$组（如果$Nsel$为奇数，则在分组时不考虑最后一个个体）

（8）变异操作

这就是模拟自然界中的基因变异，我们看下面这个例子：

变异前个体 ： 1 0 1 0 1 0 1 1

随机产生两个变异位置a和b，如a=2，b=5，则变异操作就是将两位置中间的基因进行逆序排列。

变异后个体 ： 1 1 0 1 0 0 1 1

由于自然界中变异是概率性事件且变异概率较小，所以我们在编写代码时也要给变异操作一个概率并且取值较小。

（9）重组操作

经过上述选择操作、交叉操作和变异操作得到$Nsel$个个体，由于父代种群数目为$NIND$，所以还需要$NIND-Nsel$个个体才可以对父代种群Chrom进行更新。这$NIND-Nsel$个个体就从父代种群Chrom中找出适应度值排在前$NIND-Nsel$位的$NIND-Nsel$个个体，然后添加到子代种群SelCh中。至此新的父代种群Chrom更新完毕，作为下一次选择操作的种群。

三、求解流程

四、MATLAB程序实现

（1）判断函数：

1
2
3
4
5
6
7
8
9
10
11
%% 判断一个个体是否满足背包的载重量约束，1表示满足，0表示不满足
%输入Individual：          个体
%输入w：                   各个物品的质量
%输入cap：                 背包的载重量
%输出flag：                表示一个个体是否满足背包的载重量约束，1表示满足，0表示不满足
function flag=judge_individual(Individual,w,cap)
pack_item= Individual==1;                   %判断第i个位置上的物品是否装包，1表示装包，0表示未装包
w_pack=w(pack_item);                        %找出装进背包中物品的质量
total_w=sum(w_pack);                        %计算装包物品的总质量
flag= total_w<=cap;                         %如果装包物品的总质量小于等于背包的载重量约束，则为1，否则为0
end

（2）约束处理函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
%% 对违反约束的个体进行修复
%输入Individual：          个体
%输入w：                   各个物品的质量
%输入p：                   各个物品的价值
%输入cap：                 背包的载重量
%输出Individual：          修复后的个体
function Individual=repair_individual(Individual,w,p,cap)
%% 判断一个个体是否满足背包的载重量约束，1表示满足，0表示不满足
flag=judge_individual(Individual,w,cap);
%% 只有不满足约束的个体才进行修复
if flag==0
    %% 初步修复
    pack_item=find(Individual==1);                  %找出装进背包中物品的序号
    num_pack=numel(pack_item);                      %装进背包中物品的总数目
    w_pack=w(pack_item);                            %找出装进背包中物品的质量
    total_w=sum(w_pack);                            %计算装包物品的总质量
    p_pack=p(pack_item);                            %找出装进背包中物品的价值
    ratio_pack=p_pack./w_pack;                      %计算装进背包中物品的性价比=价值/质量
    [~,rps_index]=sort(ratio_pack);                 %将已经装进包中的物品按照性价比（性价比=价值?质量）由低到高进行排序
    %% 按照rps_index顺序，依次将物品从背包中移除
    for i=1:num_pack
        remove_item=pack_item(rps_index(i));        %被移除的物品的序号
        %如果移除该物品后满足背包的载重量约束，则将该物品对应的基因位改为0，然后终止循环
        if (total_w-w_pack(rps_index(i)))<=cap
            total_w=total_w-w_pack(rps_index(i));   %装包中物品总质量减少
            Individual(remove_item)=0;              %将该物品对应的基因位改为0
            break;
        else
            %如果移除该物品后依然不满足背包的载重量约束，则也要将该物品对应的基因位改为0，然后继续移除其它物品
            total_w=total_w-w_pack(rps_index(i));   %装包中物品总质量减少
            Individual(remove_item)=0;              %将该物品对应的基因位改为0
        end
    end
    %% 进一步修复
    unpack_item=find(Individual==0);                %找出此时未装进背包中物品的序号
    num_unpack=numel(unpack_item);                  %此时未装进背包中物品的总数目
    w_unpack=w(unpack_item);                        %找出此时未装进背包中物品的质量
    p_unpack=p(unpack_item);                        %找出此时未装进背包中物品的价值
    ratio_unpack=p_unpack./w_unpack;                %计算此时未装进背包中物品的性价比=价值/质量
    [~,rups_index]=sort(ratio_unpack,'descend');    %将此时未装进包中的物品按照性价比（性价比=价值?质量）由高到低进行排序
    %% 按照rups_index顺序，依次将物品装包
    for j=1:num_unpack
        pack_wait=unpack_item(rups_index(i));       %待装包物品编号
        %如果装包该物品后满足背包的载重量约束，则将该物品对应的基因位改为1，然后继续装包其它物品
        if (total_w+w_unpack(rups_index(i)))<=cap
            total_w=total_w+w_unpack(rups_index(i));%装包中物品总质量增加
            Individual(pack_wait)=1;                %将该物品对应的基因位改为1
        else
            %如果装包该物品后不满足背包的载重量约束，则终止循环
            break;
        end
    end
end
end

（3）编码函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
%% 编码，生成满足约束的个体
%输入n：                   物品数目
%输入w：                   各个物品的质量
%输入p：                   各个物品的价值
%输入cap：                 背包的载重量
%输出Individual：          满足背包载重量约束的个体
function Individual=encode(n,w,p,cap)
Individual=round(rand(1,n));                                %随机生成n个数字（每个数字是0或1）
flag=judge_individual(Individual,w,cap);                    %判断Individual是否满足背包的载重量约束，1表示满足，0表示不满足
%% 如果flag为0，则需要修复个体Individual。否则，不需要修复
if flag==0
    Individual=repair_individual(Individual,w,p,cap);       %修复个体Individual
end
end

（4）种群初始化函数：

1
2
3
4
5
6
7
8
9
10
11
12
%% 初始化种群
%输入NIND：                种群大小
%输入n：                   物品数目
%输入w：                   各个物品的质量
%输入p：                   各个物品的价值
%输入cap：                 背包的载重量
%输出Chrom：               初始种群
function Chrom=InitPop(NIND,N,w,p,cap)
Chrom=zeros(NIND,N);                    %用于存储种群
for i=1:NIND
    Chrom(i,:)=encode(N,w,p,cap);       %编码，生成满足约束的个体
end

（5）目标函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%% 计算单个染色体的装包物品总价值和总重量
%输入n：                      物品数目
%输入Individual：             个体
%输入p：                      各个物品价值
%输入w：                      各个物品质量
%输出sumP：                   该个体的装包物品总价值
%输出sumW：                   该个体的装包物品总重量
function [sumP,sumW]=Individual_P_W(n,Individual,p,w)
sumP=0;
sumW=0;
for i=1:n
    %如果为1，则表示物品被装包
    if Individual(i)==1
        sumP=sumP+p(i);
        sumW=sumW+w(i);  
    end
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
%% 计算种群中每个染色体的物品总价值
%输入Chrom：               种群
%输入p：                   各个物品的价值
%输入w：                   各个物品的质量
%输出Obj：                 种群中每个个体的物品总价值
function Obj=Obj_Fun(Chrom,p,w)
NIND=size(Chrom,1);        %种群大小
n=size(Chrom,2);           %物品数目
Obj=zeros(NIND,1);
for i=1:NIND
    Obj(i,1)=Individual_P_W(n,Chrom(i,:),p,w);
end
end

（6）轮盘赌选择操作函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
%% 选择操作
%输入Chrom：               种群
%输入FitnV：               适应度值
%输入GGAP：                代沟
%输出SelCh：               被选择的个体
function SelCh=Select(Chrom,FitnV,GGAP)
NIND=size(Chrom,1);                 %种群数目
Nsel=NIND*GGAP;
total_FitnV=sum(FitnV);             %所有个体的适应度之和
select_p=FitnV./total_FitnV;        %计算每个个体被选中的概率
select_index=zeros(Nsel,1);         %储存被选中的个体序号
%对select_p进行累加操作，c(i)=sum(select_p(1:i))
%如果select_p=[0.1,0.2,0.3,0.4]，则c=[0.1,0.3,0.6,1]
c=cumsum(select_p);
%% 循环NIND次，选出NIND个个体
for i=1:Nsel
    r=rand;                         %0~1之间的随机数
    index=find(r<=c,1,'first');     %每次被选择出的个体序号
    select_index(i,1)=index;
end
SelCh=Chrom(select_index,:);        %被选中的个体

（7）交叉操作函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
%% 交叉操作
%输入SelCh：               被选择的个体
%输入Pc：                  交叉概率
%输出SelCh：               交叉后的个体
function SelCh=Crossover(SelCh,Pc)
[NSel,n]=size(SelCh);                               %n为染色体长度
for i=1:2:NSel-mod(NSel,2)
    if Pc>=rand %交叉概率Pc
        cross_pos=unidrnd(n);                       %随机生成一个1~N之间的交叉位置
        
        cross_Selch1=SelCh(i,:);                    %第i个进行交叉操作的个体
        cross_Selch2=SelCh(i+1,:);                  %第i+1个进行交叉操作的个体
        
        cross_part1=cross_Selch1(1:cross_pos);      %第i个进行交叉操作个体的交叉片段
        cross_part2=cross_Selch2(1:cross_pos);      %第i+1个进行交叉操作个体的交叉片段
        
        cross_Selch1(1:cross_pos)=cross_part2;      %用第i+1个个体的交叉片段替换掉第i个个体交叉片段
        cross_Selch2(1:cross_pos)=cross_part1;      %用第i个个体的交叉片段替换掉第i+1个个体交叉片段
        
        SelCh(i,:)=cross_Selch1;                    %更新第i个个体
        SelCh(i+1,:)=cross_Selch2;                  %更新第i+1个个体
    end
end

（8）变异操作函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
%% 变异操作
%输入SelCh：               被选择的个体
%输入Pm：                  变异概率
%输出SelCh：               变异后的个体
function SelCh=Mutate(SelCh,Pm)                     
[NSel,n]=size(SelCh);                               %n为染色体长度
for i=1:NSel
    if Pm>=rand
        R=randperm(n);                              %随机生成1~n的随机排列
        pos1=R(1);                                  %第1个变异位置
        pos2=R(2);                                  %第2个变异位置
        
        left=min([pos1,pos2]);                      %更小的那个值作为变异起点
        right=max([pos1,pos2]);                     %更大的那个值作为变异终点
        
        mutate_Selch=SelCh(i,:);                    %第i个进行变异操作的个体
        mutate_part=mutate_Selch(right:-1:left);    %进行变异操作后的变异片段
        mutate_Selch(left:right)=mutate_part;       %将mutate_Selch上的第left至right位上的片段进行替换
        
        SelCh(i,:)=mutate_Selch;                    %更新第i个进行变异操作的个体
    end
end

（9）重组操作函数：

1
2
3
4
5
6
7
8
9
10
%% 重插入子代的新种群
%输入Chrom：               父代种群
%输入SelCh：               子代种群
%输入Obj：                 父代适应度
%输出Chrom：               重组后得到的新种群
function Chrom=Reins(Chrom,SelCh,Obj)
NIND=size(Chrom,1);
NSel=size(SelCh,1);
[~,index]=sort(Obj,'descend');
Chrom=[Chrom(index(1:NIND-NSel),:);SelCh];

（10）主函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
tic
clear
clc
%% 创建数据
%各个物品的质量，单位kg
w=[80,82,85,70,72,70,82,75,78,45,49,76,45,35,94,49,76,79,84,74,76,63,...
    35,26,52,12,56,78,16,52, 16,42,18,46,39,80,41,41,16,35,70,72,70,66,50,55,25, 50,55,40];
%各个物品的价值，单位元
p=[200,208,198,192,180,180,168,176,182,168,187,138,184,154,168,175,198,...
    184,158,148,174,135, 126,156,123,145,164,145,134,164,134,174,102,149,134,...
    156,172,164,101,154,192,180,180,165,162,160,158,155, 130,125];    
cap=1000;                                       %每个背包的载重量为1000kg
n=numel(p);                                     %物品个数
%% 参数设置
NIND=500;                                       %种群大小
MAXGEN=500;                                     %迭代次数
Pc=0.9;                                         %交叉概率
Pm=0.08;                                        %变异概率
GGAP=0.9;                                       %代沟
%% 初始化种群
Chrom=InitPop(NIND,n,w,p,cap);
%% 优化
gen=1;
bestIndividual=Chrom(1,:);                      %初始将初始种群中一个个体赋值给全局最优个体
bestObj=Individual_P_W(n,bestIndividual,p,w);   %计算初始bestIndividual的物品总价值
BestObj=zeros(MAXGEN,1);                        %记录每次迭代过程中的最优适应度值
while gen<=MAXGEN
    %% 计算适应度
    Obj=Obj_Fun(Chrom,p,w);                     %计算每个染色体的物品总价值
    FitnV=Obj;                                  %适应度值=目标函数值=物品总价值
    %% 选择
    SelCh=Select(Chrom,FitnV,GGAP);
    %% 交叉操作
    SelCh=Crossover(SelCh,Pc);
    %% 变异
    SelCh=Mutate(SelCh,Pm);
    %% 重插入子代的新种群
    Chrom=Reins(Chrom,SelCh,Obj);
    %% 将种群中不满足载重量约束的个体进行约束处理
    Chrom=adjustChrom(Chrom,w,p,cap);
    %% 记录每次迭代过程中最优目标函数值
    [cur_bestObj,cur_bestIndex]=max(Obj);       %在当前迭代中最优目标函数值以及对应个体的编号
    cur_bestIndividual=Chrom(cur_bestIndex,:);  %当前迭代中最优个体
    %如果当前迭代中最优目标函数值大于等于全局最优目标函数值，则进行更新
    if cur_bestObj>=bestObj
        bestObj=cur_bestObj;
        bestIndividual=cur_bestIndividual;
    end
    BestObj(gen,1)=bestObj;                     %记录每次迭代过程中最优目标函数值
    %% 打印每次迭代过程中的全局最优解
    disp(['第',num2str(gen),'次迭代的全局最优解为：',num2str(bestObj)]);
    %% 更新迭代次数
    gen=gen+1 ;
end
%% 画出迭代过程图
figure;
plot(BestObj,'LineWidth',1);
xlabel('迭代次数');
ylabel('目标函数值（物品总价值）');
%% 最终装进包中的物品序号
pack_item=find(bestIndividual==1);
%% 计算最优装包方案的物品总价值和总重量
[bestP,bestW]=Individual_P_W(n,bestIndividual,p,w);
toc

五、实验结果展示

1. 《MATLAB智能优化算法》（曹旺著） ↩