权重改进的粒子群算法

在粒子群算法中,惯性权重w是最重要的参数。

自适应权重法

自适应权重法的基本策略是:对于具有高适应度的粒子,减小惯性权重以增加局部搜索能力;对于适应度较低的粒子,增加惯性权重以增加全局搜索能力。

依据早熟收敛程度和适应度值进行调整

设定粒子 p_i 的适应度值为 f_i,最优粒子适应度为 f_m,粒子群平均适应度为 f_{avg},优于平均适应度的粒子适应度平均值为f_{avg}^′,定义 Δ=|f_m−f_{avg}^′ |

如果 f_i≥f_{avg},那么

w=w−(w−w_{min})⋅|(f_i−f_{avg}^′)/(f_m−f_{avg}^′ )|

如果 f_{avg}≤f_i,则不改变权重。

如果 f_i ,则

w=1.5−1/(1+k_1⋅exp⁡(−k_2⋅Δ) )

其中k_1 控制w的上限,k_2 用于控制调节能力。

根据全局最优点的距离进行调整

w=w_{min} −\frac{(w_{max} −w_{min} )×(f−f_{min} )}{f_{avg}−f_{min} }, f≤f_{avg}

w=w_{max} ,f>f_{avg}

随机权重法

w=w^′+σ⋅N(0,1)

w^′=w_{min} +(w_{max} −w_{min} )×rand(0,1)

线性递减权重法

w=w_{max} −(t×(w_{max} −w_{min} ))/t_{max}

其中t为迭代次数。

分布估计算法解决旅行商问题(TSP)

在用分布估计算法解决旅行商问题时,结构与传统的分布估计算法相似,只不过是把概率向量换成了“概率矩阵”而已:

  1. 通过概率矩阵生成解。
  2. 评估解。
  3. 选择优势群体。
  4. 更新概率矩阵。
  5. 重复以上4步直至迭代结束。

这里说的“概率矩阵”记录了上一代优势群体中,“城市对”出现的次数(或与城市对出现次数成正比的一个数)。“城市对”是指路线中相邻的两个城市,这两个城市不分先后。假设概率矩阵第i行第j列的元素为 p_{ij},它代表在优势群体中城市i和城市j相邻的次数(即i−j和j−i发生的次数)。

经过测试(于2018-11-09),这种算法非常容易陷入局部最优,效果不是很好。

下面是程序运行结果:

每代最短距离变化图

每代平均距离变化图

最短路线图

主程序

主程序如下:

function main()
pos = [randperm(20); randperm(20)];
% pos = load('berlin52.txt'); % 7542
% pos = pos(:, 2:3);
% pos = pos';

popsize = 300;                  % 种群规模
maxgen = 100;                   % 最大迭代次数
citysize = size(pos, 2);        % 城市数量
p = ones(citysize, citysize);   % 概率矩阵
sn = ceil(popsize * 0.5);       % 优势群体规模
dm = makedm(pos);               % 距离矩阵

bestlists = zeros(maxgen, citysize);    % 记录每代最优解
bestfits = zeros(1, maxgen);            % 记录每代最优解适应度值
avgfits = zeros(1, maxgen);             % 记录每代平均适应度值

for gen = 1:maxgen
    pop = makepop(popsize, citysize, p);    % 制作种群

    fitness = callength(pop, dm);
    [bfit, bfiti] = min(fitness);
    bestlists(gen, :) = pop(bfiti, :);
    bestfits(1, gen) = bfit;
    avgfits(1, gen) = sum(fitness) / popsize;

    spop = selection(pop, sn, fitness);     % 选择优势群体
    p = makep(spop);                        % 更新概率矩阵
end
fitness = callength(bestlists, dm);
[bfit, bfiti] = min(fitness);
fprintf("最短距离: %f\n", bfit);
end

制作种群函数

对于每个个体,第一个城市是随机生成的,之后每个城市都根据上一个城市进行轮盘赌法选择。有时候上一代优势群体中可能不会出现所有可能的城市对,这意味着该城市不会出现在候选城市中,但是路径必须经过每一个城市,应该处理这种情况,这种情况也可以通过给概率矩阵设置最小值以避免。

function pop = makepop(popsize, citysize, p)
% 制作种群
% popsize       input  种群规模
% citysize      input  城市数量
% p             input  概率矩阵
% pop           output 种群
pop = zeros(popsize, citysize);

for i = 1:popsize
    for j = 1:citysize
        if j == 1   % 随机产生第一个城市
            pop(i, j) = randi(citysize);
        else        % 使用轮盘赌法选择下一个城市
            cityp = p(pop(i, j-1), :);      % 找出备选城市
            cityp(pop(i, 1:j-1)) = 0;       % 已经访问过的概率为0
            sumcityp = sum(cityp);          % 开始轮盘赌法
            cityp = cityp ./ sumcityp;
            cityp = cumsum(cityp);
            index = find(cityp >= rand);

            % 上一代优势群体中没有城市对记录,则找一个城市填充
            if isempty(index)
                tempcity = 1:1:citysize;
                tempcity(pop(i, 1:j-1)) = 0;
                tempindex = find(tempcity > 0);
                index = tempindex(1);
            end
            pop(i, j) = index(1);
        end
    end
end
end

更新概率矩阵函数

根据上面对概率矩阵的定义更新概率矩阵:

function p = makep(spop)
% 更新概率矩阵
% spop          input  优势群体
% np            output 概率矩阵
[popsize, citysize] = size(spop);
count = zeros(citysize, citysize);   % 记录出现次数
spop = [spop spop(:, 1)];
for i = 1:popsize
    for j = 1:citysize
        a = spop(i, j);
        b = spop(i, j+1);
        count(a, b) = count(a, b) + 1;
        count(b, a) = count(b, a) + 1;
    end
end
sumcount = sum(count, 2);
p = count ./ sumcount;
end

一些其他函数

选择函数:

function spop = selection(pop, sn, fitness)
% 选择
% pop       input  种群
% sn        input  优势群体规模
% fitness   input  种群适应度值
% spop      output 优势群体
[~, index] = sort(fitness);
index = index(1: sn);
spop = pop(index, :);
end

计算路径长度函数:

function len = callength(pop, dm)
% 求路径距离
% pop           input  种群
% dm            input  距离矩阵
% len           output 距离
[NP, D] = size(pop);
pop = [pop pop(:,1)];
sum = zeros(NP, 1);
for i = 1:NP
    for j = 1:D
        sum(i,1) = sum(i,1) + dm(pop(i, j), pop(i, j+1));
    end
end
len = sum;
end

制作距离矩阵函数:

function dm = makedm(pos)
% 制作距离矩阵
% pos       input  城市坐标
% dm        output 距离矩阵
[~, len] = size(pos);
deltax = repmat(pos(1,:)', 1, len) - repmat(pos(1,:), len, 1);
deltay = repmat(pos(2,:)', 1, len) - repmat(pos(2,:), len, 1);
dm = round((deltax .^ 2 + deltay .^ 2) .^ 0.5);
end

粒子群算法求函数最小值

求下面函数的最小值:

f(x) = \sum_{i=1}^{30} (x_i^2 + x_i – 6 )

程序运行结果如下:

函数最小值: -182.160634

每代最优个体适应值变化图

主函数

主函数首先初始化种群,对于第1代种群,个体极值和全局极值都在本代种群中;之后进行迭代,每次迭代根据公式更新速度和位置,并更新个体极值和全局极值,重复此过程直至迭代结束。

function main()
popsize = 50;       % 种群规模
birdsize = 30;      % 粒子数量
w = 0.5;            % 惯性权重
c1 = 1.0;           % 认知因子
c2 = 2.0;           % 社会因子
maxgen = 100;       % 最大迭代次数

% 初始化
x = randn(popsize, birdsize);
v = randn(popsize, birdsize);

% 初始化pid,pgd
fitness = calfitness(x);
pid = x;
pidfit = fitness;
[bfit, bfiti] = min(fitness);
pgd = x(bfiti, :);
pgdfit = bfit;

% 记录每代最优值
bestpidfit = zeros(popsize, 1);

for gen = 1:maxgen
    % 更新速度和位置
    v = w .* v + c1 .* rand .* (pid - x) + ... 
        c2 .* rand .* (repmat(pgd, popsize, 1) - x);
    x = x + v;

    % 更新pid,pgd
    fitness = calfitness(x);
    index = find(fitness < pidfit);
    pid(index, :) = x(index, :);
    pidfit(index, 1) = fitness(index, 1);
    [bfit, bfiti] = min(fitness);
    bestpidfit(gen, 1) = bfit;
    if bfit < pgdfit
        pgd = x(bfiti, :);
        pgdfit = bfit;
    end
end
fprintf("函数最小值: %f\n", pgdfit);
figure(1);
plot(1:maxgen, bestpidfit);
title("每代最优适应度值变化曲线");
end

适应值函数

function fitness = calfitness(x)
% 计算适应度值
% f = sum(x^2+x-6)
% x         input  种群
% fitness   output 适应度值
x = x .^ 2 + x - 6;
fitness = sum(x, 2);
end

粒子群算法简单介绍

原理

粒子群算法(也称粒子群优化算法(particle swarm optimization, PSO)),模拟鸟群随机搜索食物的行为。粒子群算法中,每个优化问题的潜在解都是搜索空间中的一只鸟,叫做“粒子”。所有的粒子都有一个由被优化的函数决定的适应值(fitness value),每个粒子还有一个速度决定它们“飞行”的方向和距离。

粒子群算法初始化为一群随机的粒子(随机解),然后根据迭代找到最优解。每一次迭代中,粒子通过跟踪两个极值来更新自己:第1个是粒子本身所找到的最优解,这个称为个体极值;第2个是整个种群目前找到的最优解,这个称为全局极值。也可以不用整个种群,而是用其中的一部分作为粒子的邻居,称为局部极值。

假设在一个D维搜索空间中,有N个粒子组成一个群落,其中第i个粒子表示为一个D维的向量:

X_i = (x_{i1}, x_{i2}, \ldots, x_{iD}), i=1,2,\ldots ,N

第i个粒子的速度表示为:

V_i = (v_{i1}, v_{i2}, \ldots, v_{iD}), i=1,2, \dots, N

还要保存每个个体的已经找到的最优解p_{best},和一个整个群落找到的最优解g_{best}

第i个粒子根据下面的公式更新自己的速度和位置:

v_{id} = w \times v_{id} + c_1 r_1 (p_{id}-x_{id}) + c_2 r_2 (p_{gd}-x_{id}) ,(1)

x_{id} = x_{id} + v_{id}

其中,p_{id} 是个体已知最优解,p_{gd} 是种群已知最优解,w 为惯性权重,c_1, c_2 为学习因子(或加速常数 acceleration constant),r_1,r_2 是[0,1]范围内的随机数。

式(1)由三部分组成:

  1. 惯性或动量部分,反应粒子的运动习惯。
  2. 认知部分,粒子有向自身历史最佳位置逼近的优势。
  3. 社会部分,粒子有向群体或领域历史最佳位置逼近的趋势。

特点

  • 粒子群算法是一种高效的并行搜索算法。
  • 粒子群算法保留了基于种群的全局搜索策略,操作模型比较简单,避免了复杂的遗传操作。
  • 保留了每个粒子的个体历史极值。
  • 算法在后期收敛速度缓慢。
  • 粒子群算法对种群大小不十分敏感。

遗传算法求解混合流水车间调度问题(HFSP)三:算法实现二

初始化种群的函数:

function pop = initpop(popsize, piecesize)
% 初始化种群
% popsize           input  种群规模
% piecesize         input  工件数量
% pop               output 种群
pop = zeros(popsize, piecesize);
for i = 1:popsize
    pop(i, :) = randperm(piecesize);
end
end

计算适应度值的函数:

function fitness = calfitness(objvalue)
% 计算适应度值
% objvalue      input  目标函数值
% fitness       output 适应度值
fitness = 1 ./ objvalue;
end

制作加工流程矩阵的函数,该矩阵行代表机器,按照工序进行排序(假设工序1有两个机器,则第1行为工序1机器1,第2行为工序1机器2,第3行为工序2机器1……),列代表时间:

function gantt = makegantt(ptr, per, equsize)
% 制作加工流程矩阵
% ptr       input  工件时间记录
% per       input  工件设备记录
% equsize   input  工序设备数量
% gantt     output 加工流程矩阵
finaltime = max(max(ptr));
[piecesize, prosize] = size(per);
cumsumequ = cumsum(equsize);
gantt = zeros(sum(equsize), finaltime);
for pro = 1:prosize
    for i = 1:piecesize
        if pro == 1
            equ = per(i, pro);
        else
            equ = cumsumequ(pro - 1) + per(i, pro);
        end
        starttime = ptr(i, pro*2-1);
        endtime = ptr(i, pro*2);
        gantt(equ, starttime:endtime) = i; 
    end
end
end

生成两个随机数的函数:

function [rl, ru] = makerlru(maxnum)
% 制作两个随机整数, rl < ru
% maxnum        input  最大数
% rl            output 小随机数
% ru            output 大随机数
r1 = randi(maxnum);
r2  =randi(maxnum);
while r2 == r1
    r2 = randi(maxnum);
end
rl = min([r1, r2]);
ru = max([r1, r2]);
end

遗传算法求解混合流水车间调度问题(HFSP)二:算法实现一

遗传算法的设计

  • 编码:对工件进行优先级编码,编码越小,优先级越高。
  • 解码:按照工件优先级进行生产,求出整体完工时间。
  • 目标函数值:整体完工时间。
  • 适应度值:目标函数越小,适应度值越大。
  • 选择:按照轮盘赌方法进行选择。
  • 交叉:按照交叉概率,选择两个父代个体(父1和父2),交叉生成两个子代个体(子1和子2)。生成一个随机的交换空间,交换空间内,子1基因来自于父2,子2基因来自于父1;交换空间外,子1基因来自于父1,子2基因来自于父2。注意任意两个零件优先级不能相等。
  • 变异:按照变异概率,随机交换两个基因。
  • 终止条件:迭代固定次数。

计算结果

在本例中,有6个工件,有3道工序,每道工序有2台机器,下面是执行结果:

最优序列:
3 4 6 2 1 5

加工流程图

上图中,第1、2行是第1工序的2台设备,第3、4行是第2工序的2台设备,第5、6行是第3工序的两台设备,纵轴代表时间。按照最优序列3 4 6 2 1 5赋予每个零件优先级,一共用时25.

主函数

首先定义问题的参数:

piecetime = [2 4 6; ...             % 设备加工时间
    4 9 2; 4 2 8; 9 5 6; 5 2 7; 9 4 3];
equsize = [2 2 2];                  % 每个工序设备数目
piecesize = size(piecetime, 1);     % 工件数目
prosize = size(piecetime, 2);       % 工序数目

在本例中,一共有6个设备,3道工序,设备必须按照1-2-3的工序进行加工,每个工序有2台机器。一共有6个工件。piecetime中行代表工件,列代表工序,内容是加工时间,比如第1行第1列,表示工件1在工序1加工时间为2.

定义遗传算法的参数:

popsize = 20;       % 种群规模
cr = 0.6;           % 交叉概率
mr = 0.05;          % 变异概率
maxgen = 100;       % 迭代次数

进行迭代:

pop = initpop(popsize, piecesize);
for gen = 1:maxgen
    [objvalue, ptr, per] = calobjvalue(pop, piecetime, equsize);    
    fitness = calfitness(objvalue);     % 计算适应度值
    pop = selection(pop, fitness);      % 选择
    pop = crossover(pop, cr);           % 交叉
    pop = mutation(pop, mr);            % 变异
end

主函数全部代码如下:

function main()
piecetime = [2 4 6; ...             % 设备加工时间
    4 9 2; 4 2 8; 9 5 6; 5 2 7; 9 4 3];
equsize = [2 2 2];                  % 每个工序设备数目
piecesize = size(piecetime, 1);     % 工件数目
prosize = size(piecetime, 2);       % 工序数目

popsize = 20;       % 种群规模
cr = 0.6;           % 交叉概率
mr = 0.05;          % 变异概率
maxgen = 100;       % 迭代次数

bestobjvalue = zeros(1, maxgen);
bestpop = zeros(maxgen, piecesize);
avgobjvalue = zeros(1, maxgen);
bestptr = cell(1, maxgen);
bestper = cell(1, maxgen);

pop = initpop(popsize, piecesize);
for gen = 1:maxgen
    [objvalue, ptr, per] = calobjvalue(pop, piecetime, equsize);

    [bobjvalue, bindex] = min(objvalue);
    bestobjvalue(1, gen) = bobjvalue;
    bestpop(gen, :) = pop(bindex, :);
    avgobjvalue(1, gen) = sum(objvalue) / popsize;
    bestptr{1, gen} = ptr{1, bindex};
    bestper{1, gen} = per{1, bindex};

    fitness = calfitness(objvalue);     % 计算适应度值
    pop = selection(pop, fitness);      % 选择
    pop = crossover(pop, cr);           % 交叉
    pop = mutation(pop, mr);            % 变异
end

[~, finalindex] = min(bestobjvalue);
finalptr = bestptr{1, finalindex};
finalper = bestper{1, finalindex};

fprintf("最优序列:\n");
disp(bestpop(finalindex, :));

gantt = makegantt(finalptr, finalper, equsize);

figure(1);
imagesc(gantt);
colorbar;
title("加工流程图");


figure(2);
plot(1:maxgen, bestobjvalue);
title("最优时间变化图");
xlabel("代数"); ylabel("最优时间");

figure(3);
plot(1:maxgen, avgobjvalue);
title("平均时间变化图");
xlabel("代数"); ylabel("平均时间");

end

计算目标函数值函数

在第一道工序中,所有的工件同时等待被加工,则按照优先级进行加工;在第二道和之后的工序中,由于上一道工序中工件完工时间不同,上一道工序先加工完的工件先进行本工序加工。

function [objvalue, ptr, per] = calobjvalue(pop, piecetime, equsize)
% 计算目标函数值
% pop           input  种群
% piecetime     input  工件加工时间
% equsize       input  每个工序设备数量
% objvalue      output 目标函数值(加工时间)
% ptr           output 工件加工时间记录,cell
% per           output 工件加工设备记录,cell
[popsize, piecesize] = size(pop);
prosize = size(equsize, 2);
objvalue = zeros(popsize, 1);
ptr = cell(1, popsize);
per = cell(1, popsize);
for i = 1:popsize
    pieceweight = pop(i, :);

    % 设备状态序列
    % [工序1设备1 工序1设备2 工序2设备1 工序2设备2 ……]
    % 记录当前设备使用结束时间,默认为0表示未开始
    equstu = zeros(1, sum(equsize));

    % 对设备状态序列的工序分隔符
    % 大于等于当前设备最小值的索引是当前设备所处的工序
    % [2 3 5] 工序1有2台设备 工序2有1台设备 工序3有2台设备
    prosep = cumsum(equsize);

    % 工件时间记录,记录每个工件每个工序的开始时间和结束时间
    % 行表示工件,相邻两列表示开始加工时间和停止加工时间
    % [1 2 2 3; 4 5 6 7]
    % 表示工件1第1工序加工时间为1-2,第2工序加工时间为2-3
    % 工件2第1工序加工时间为4-5,第2工序加工时间为6-7
    piecetimerecord = zeros(piecesize, prosize*2);

    % 工件设备记录,记录每个工件在工序中的加工设备
    % 行数表示工件,列表示该零件在每个工序加工设备
    % [1 2; 2 1]
    % 表示工件1在第1工序加工设备为1,第2工序加工设备为2
    % 工件2在第1工序加工设备为2,第2工序加工设备为1
    pieceequrecord = zeros(piecesize, prosize);

    % 对每一道工序
    %   如果是第1道工序,对工件按优先级排序
    %     其余工序上上一道工序完工时间对工件排序
    %   对排序后的每一件工件
    %     对该工序中可用机器按使用结束时间排序
    %     使用使用结束时间最小的机器
    %     加工开始时间为max{设备使用结束时间,零件上一工序完工时间}
    %     加工结束时间=加工开始时间+加工时间
    %     更新各个状态和记录矩阵
    for pro = 1:prosize
        if(pro == 1)
            [~, piecelist] = sort(pieceweight);
        else
            tempendtime = piecetimerecord(:, (pro-1)*2);
            tempendtime = tempendtime';
            [~, piecelist] = sort(tempendtime);
        end
        for pieceindex = 1:length(piecelist)
            piece = piecelist(pieceindex);
            equtimelist = equstu(prosep(pro)-equsize(pro)+1:prosep(pro));
            [equtime, equlist] = sort(equtimelist);
            equ = equlist(1);
            if pro == 1
                piecestarttime = 0;
            else
                piecestarttime = piecetimerecord(piece, pro*2-2);
            end
            starttime = max(equtime(1), piecestarttime) + 1;
            endtime = starttime + piecetime(piece, pro) - 1;
            equstuindex = prosep(pro)-equsize(pro)+equ;
            equstu(equstuindex) = endtime;
            piecetimerecord(piece, pro*2-1) = starttime;
            piecetimerecord(piece, pro*2) = endtime;
            pieceequrecord(piece, pro) = equ;
        end
    end
    objvalue(i, 1) = max(max(piecetimerecord));
    ptr{1, i} = piecetimerecord;
    per{1, i} = pieceequrecord;
end
end

选择函数

使用轮盘赌方法进行选择:

function spop = selection(pop, fitness)
% 轮盘赌选择
% pop       input  种群
% fitness   input  适应度值
% spop      output 选择后生成的种群
[popsize, piecesize] = size(pop);
spop = zeros(popsize, piecesize);
sumfit = sum(fitness);
fitness = fitness ./ sumfit;
fitness = cumsum(fitness);
r = rand(1, popsize);
r = sort(r);
j = 1;
for i = 1:popsize
    while fitness(j) < r(i)
        j = j + 1;
    end
    spop(i, :) = pop(j, :);
end

% 由于上面轮盘赌方法特殊性,一个个体在相邻位置多次重复,故随机排序
rr = randperm(popsize);
spop(:, :) = spop(rr, :);
end

交叉函数

按照交叉概率,选择两个父代个体(父1和父2),交叉生成两个子代个体(子1和子2)。生成一个随机的交换空间,交换空间内,子1基因来自于父2,子2基因来自于父1;交换空间外,子1基因来自于父1,子2基因来自于父2。注意任意两个零件优先级不能相等。

function cpop = crossover(pop, cr)
% 交叉
% pop       input  种群
% cr        input  交叉概率
% cpop      output 交叉后种群
[popsize, piecesize] = size(pop);
cpop = pop;

if mod(popsize,2) ~= 0
    nn = popsize - 1;
else
    nn = popsize;
end

% 父代father mother, 子代son daughter
% 在rl:ru中,son继承mother,daughter继承father
% 其余位置son继承father,daughter继承mother
for i = 1:2:nn
    if rand > cr
        continue;
    end
    [rl, ru] = makerlru(piecesize);

    father = pop(i, :);
    mother = pop(i+1, :);
    if father == mother
        continue;
    end
    son = zeros(1, piecesize);
    daughter = zeros(1, piecesize);
    son(rl:ru) = mother(rl:ru);
    daughter(rl:ru) = father(rl:ru);

    j = 1;
    for k = 1:piecesize
        if k >= rl && k <= ru
            continue;
        end
        while ~isempty(find(son == father(j), 1))
            j = j + 1;
        end
        son(k) = father(j);
    end

    j = 1;
    for k = 1:piecesize
        if k >= rl && k <= ru
            continue;
        end
        while ~isempty(find(daughter == mother(j), 1))
            j = j + 1;
        end
        daughter(k) = mother(j);
    end

    cpop(i, :) = son;
    cpop(i+1, :) = daughter;
end
end

变异函数

随机交换染色体中两个位置的基因:

function mpop = mutation(pop, mr)
% 变异,交换两个随即位置的基因
% pop       input  种群
% mr        input  变异概率
% mpop      output 变异后种群
[popsize, piecesize] = size(pop);
mpop = pop;
for i = 1:popsize
    if rand > mr
        continue;
    end
    r1 = randi(piecesize);
    r2 = randi(piecesize);
    temp  = mpop(i, r1);
    mpop(i, r1) = mpop(i, r2);
    mpop(i, r2) = temp;
end
end

遗传算法求解混合流水车间调度问题(HFSP)一:问题介绍

混合流水车间调度问题(Hybrid Flow-shop Scheduling Problem, HFSP)是车间调度中的一类经典问题。混合流水车间调度问题,在一道工序有一台或多台机器,工件的加工需要满足一定的工艺顺序。

假设和约束

  1. 一个工件在一道工序上被任意一个机器加工。
  2. 一个机器在某一时刻只能空闲或加工一个工件。
  3. 工件必须按照加工工序顺序进行加工。
  4. 同一道工序中机器都相同。
  5. 工件加工过程不允许中断。
  6. 如果多个工件同时需要被加工,则按优先级顺序进行加工。

需要解决的问题

确定零件的加工优先级,以使整体加工时间最短。

解决思路

在第一道工序中,所有的工件同时等待被加工,则按照优先级进行加工;在第二道和之后的工序中,由于上一道工序中工件完工时间不同,上一道工序先加工完的工件先进行本工序加工。求出整体完工时间作为目标函数值,运用遗传算法求解以使目标函数值最小。

分布估计算法求解0-1背包问题二

一些其他函数

重量计算函数:

function wgtsum = weightsum(pop, weights)
% 计算种群的重量
% pop           input  种群
% weights       input  重量向量
% wgtsum        output 种群重量
popsize = size(pop, 1);
wgtsum = zeros(popsize, 1);
for i = 1:popsize
    wgtsum(i, 1) = weightsumv(pop(i, :), weights);
end
end
function wgtsum = weightsumv(stuffs, weights)
% 计算一个个体的重量
% stuffs        input  物品序列
% weights       input  重量向量
% wgtsum        output 个体重量
wgtsum = sum(weights(stuffs ~= 0));
end

收益计算函数:

function pftsum = profitssum(pop, profits)
% 计算种群收益
% pop       input  种群
% profits   input  收益向量
% pftsum    output 种群收益
popsize = size(pop, 1);
pftsum = zeros(popsize, 1);
for i = 1:popsize
    pftsum(i, 1) = sum(profits(pop(i, :) ~= 0));
end
end

分布估计算法求解0-1背包问题一

0-1背包问题是:有一个固定容量的背包,和固定种类的物品,每种物品只有一件。每件物品有各自的价值和重量,求解哪些物品放入背包可以使价值总和最大,且不超过背包容量。

本例中用分布估计算法求解0-1背包问题结果如下:

每代最优选择收益图

每代平均收益图

本例中的算例在下面下载:

0-1背包问题算例下载 people.sc.fsu.edu

0-1背包问题算例下载 本站

可以看到,分布估计算法可能在很靠前的迭代中就能得到很好的解,但是由于该算法不会保留上一代最优解,因此该解很可能丢失。从平均收益图来看,种群的整体收益平均值还是增加的。

主程序

主程序主要有以下四步:

for gen = 1:maxgen
    pop = makepop(popsize, stuffsize, p);               % 制作种群
    pop = capacitylimit(pop, capacity, weights, p);     % 限制重量
    spop = selection(pop, sn, profits);     % 选择优势个体
    p = makep(spop, p, alpha);              % 更新概率向量
end

根据概率向量随机采样得到种群,并限制种群中个体的重量(不能超过背包容量),之后选择优势个体,并根据优势个体更新概率向量。进行下次迭代。

主程序如下:

function main()
capacity = load("bag\P08\p08_c.txt");       % 背包容量
bks = load("bag\P08\p08_s.txt");            % 最优解
bks = bks';
weights = load("bag\P08\p08_w.txt");        % 重量
weights = weights';
profits = load("bag\P08\p08_p.txt");        % 收益
profits = profits';

popsize = 100;                  % 群体规模
maxgen = 50;                    % 迭代次数
stuffsize = length(weights);    % 物品数量
p = ones(1, stuffsize) .* 0.5;  % 概率向量
alpha = 1;                      % 学习速率
sn = 0.7;                       % 优势个体数量
sn = ceil(popsize * sn);

bestselection = zeros(maxgen, stuffsize);   % 记录每代最优选择
avgweights = zeros(1, maxgen);              % 记录每代平均收益

for gen = 1:maxgen
    pop = makepop(popsize, stuffsize, p);               % 制作种群
    pop = capacitylimit(pop, capacity, weights, p);     % 限制重量

    wgtsum = weightsum(pop, weights);
    [~, index] = max(wgtsum);
    bestselection(gen, :) = pop(index, :);
    avgweights(1, gen) = sum(wgtsum) / popsize;

    spop = selection(pop, sn, profits);     % 选择优势个体
    p = makep(spop, p, alpha);              % 更新概率向量
end
wgtsum = weightsum(bestselection, weights);
[~, index] = max(wgtsum);
figure(1);
plot(1:1:maxgen, wgtsum');
title("每代最优选择收益图");
figure(2);
plot(1:1:maxgen, avgweights);
title("每代平均收益图");
end

制作种群函数

制作种群函数根据概率向量p进行随机采样,直至达到种群规模。概率向量p中的一项代表在该位置上取1的概率:

function pop = makepop(popsize, stuffsize, p)
% 初始化种群,但没有限制重量
% popsize       input  种群规模
% stuffsize     input  物品数目
% p             input  概率向量
% pop           output 构造的种群
pop = zeros(popsize, stuffsize);
for i =1:popsize
    pop(i, :) = makepopv(stuffsize, p);
end
end

function pop = makepopv(stuffsize, p)
% 初始化个体,没有限制重量
% stuffsize     input  物品数目
% p             input  概率向量
% pop           output 构造的个体
tpop = rand(1, stuffsize);
pop = zeros(1, stuffsize);
for j = 1:stuffsize
    if tpop(1, j) <= p(1, j)
        pop(1, j) = 1;
    end
end
end

限制重量函数

如果种群中某一个个体重量超过背包容量,则重新生成该个体:

function npop = capacitylimit(pop, capacity, weights, p)
% 限制重量
% pop           input  种群
% capacity      input  背包容量
% weights       input  重量
% p             input  概率向量
% npop          output 新种群
npop = pop;
[popsize, stuffsize] = size(pop);
for i = 1:popsize
    wgtsum = weightsumv(npop(i, :), weights);
    while wgtsum > capacity
        npop(i, :) = makepopv(stuffsize, p);
        wgtsum = weightsumv(npop(i, :), weights);
    end
end
end

选择优势个体函数

从种群中选择指定数量的优势个体:

function spop = selection(pop, sn, profits)
% 选择
% pop           input  种群
% sn            input  选择数量
% profits       input  收益向量
% spop          output 选择的种群
pftsum = profitssum(pop, profits);
pftsum = pftsum';
[~, index] = sort(pftsum, 'descend');
index = index(1: sn);
spop = pop(index, :);
end

更新概率向量函数

该函数更新概率向量:

function np = makep(pop, p, alpha)
% 更新概率向量
% pop           input  种群
% p             input  概率向量
% alpha         input  学习速率
% np            output 更新后的概率向量
popsize = size(pop, 1);
np = (1 - alpha) .* p + alpha .* sum(pop) ./ popsize;
end

几种分布估计算法介绍

PBIL算法

用以解决二进制编码的、变量无关的优化问题。

在PBIL算法中,表示解空间的概率模型是一个概率向量:

p(x)= (p(x_1), p(x_2), p(x_3), \ldots ,p(x_n))

其中p(x_i)表示在位置i上取值为1的概率。

PBIL算法过程如下:

  1. 在每一代中,通过概率向量p(x)随机产生M个个体。
  2. 计算M个个体的适应值。
  3. 选择最优的N个个体来更新p(x),N<=M。
  4. 迭代,直至结束。

更新概率向量的方法如下:

p_{l+1}(x) = (1-\alpha)p_l(x)+\alpha \frac{1}{N} \sum_{k=1}^{N}x_l^k

其中l表示代数,α为学习速率。

UMDA算法

与PBIL算法的区别在于概率向量的更新:

p_{l+1}(x) = \frac{1}{N} \sum_{k=1}^{N}x_l^k