粒子群优化算法应用练习

最新推荐文章于 2024-11-09 21:51:15 发布

q浅夜

最新推荐文章于 2024-11-09 21:51:15 发布

阅读量1.3k

点赞数 25

文章标签：算法

本文链接：https://blog.csdn.net/qianyeguiji/article/details/135482532

版权

粒子群优化算法应用练习

负荷经济分配

参考文献：杨俊杰, 周建中, 吴玮, 等.改进粒子群优化算法在负荷经济分配中的应用[J].电网技术, 2005, 29(12): 1-4.

负荷经济分配（Economic Dispatch，ED）问题是机组组合问题的一个子优化问题。ED 问题是指当给定某运行时段机组的开停机计划后，在满足各种约束的条件下在运行机组间分配负荷，使电力系统的总运行费用最小或发电公司的利润最大。

优化目标

电力系统的总运行费用 F 最小
$\min F = \sum\limits^M_{i=1}F(P_i)$
其中，F(P_i) 为机组 i 的费用特性函数。

可以用分段二次曲线表示该类机组的混合费用特性曲线，每段曲线表示燃烧不同的燃料：
$F(P_i)=\begin{cases} a_{i,1}P_i^2+b_{i,1}P_i+c_{i,1}\ \ \ \ (P_{i,\min}\leq P_i \leq P_{i,1}) \\ a_{i,2}P_i^2+b_{i,2}P_i+c_{i,2}\ \ \ \ (P_{i,1} < P_i <P_{i,2}) \\ ...\\a_{i,k}P_i^2+b_{i,k}P_i+c_{i,k}\ \ \ \ (P_{i,k-1} < P_i \leq P_{i, \max}) \end{cases}$
式中，a_i,j 、b_i,j 、c_i,j（j = 1, 2, …, k）为机组 i 在 j 类燃料下的费用系数；P_i,j （j = 1, 2, …, k-1）为机组 i 在出力区 j 的出力上限值。

约束条件

（1）爬坡约束条件

在任意两个连续的运行时间段内改变机组 i 的出力必须满足下列约束条件：
$\max \lbrace P_{i,\min},P_{i,0}-R_{Di}\rbrace \leq P_i \leq \min \lbrace P_{i, \max},P_{i,0}+R_{Ui}\rbrace$
式中，P_i 为机组 i 当前时段的出力；P_i,0 为机组 i 在前一时段的出力；P_i,min 、P_i,max 分别为机组 i 的最小和最大出力；R_Di 、R_Ui 分别为机组 i 爬坡限制的下限值和上限值。

（2）机组出力限制区约束

确定机组的出力范围一般需综合考虑机组的技术允许最低出力、避开振动区及出力限制线等因素，机组的可行出力区可分为多段：
$\begin{cases} P_{i,\min} \leq P_i \leq P_{i,1}' \\ P_{i,j-1}' \leq P_i \leq P_{i,j}' \ \ \ \ \ \ \ \ (j=2,3,...,n_i) \\ P_{i,n_i}' \leq P_i \leq P_{i,\max} \end{cases}$
式中，P‘_i,j（j = 2, 3, …, n_i）为机组 i 在第 j 段禁止出力区的上限；j、n_i 为机组 i 禁止出力区的符号和数量。

（3）系统负荷平衡约束
$\sum\limits^M_{i=1}P_i=P_D+P_L$
式中，M为运行机组台数；P_D为系统负荷需求；P_L为损失电量。

代码练习

选择了文中较简单的案例练习：只考虑系统负荷平衡约束，且只考虑一种燃料、单个出力区。简化为较为简单的分配问题。

机组参数如下：
在这里插入图片描述

共10台机组，令粒子维数为10，每一维代表该机组的负荷值。

由于每个机组的负荷上限不同，在初始化和边界条件处理时对每一维依次进行计算。

对于负荷平衡约束条件，则设计一定的惩罚，在适应度函数中体现。

主函数：

clear all;
close all;
clc;
%% 初始化
N = 60; % 粒子个数
D = 10; % 粒子维数
T = 150; % 最大迭代次数
c1 = 2; % 学习因子1
c2 = 2; % 学习因子2
Wmax = 0.8; % 惯性权重最大值
Wmin = 0.4; % 惯性权重最小值
Pmax = [72 70 64 61 72 71 73 73 143 143]; % 每个机组的允许负荷最大值
Pmin = 0;
Vmax = 10; % 速度最大值
Vmin = -10; % 速度最小值
%% 初始化种群个体
x = ones(N, D);
for i = 1 : N
    for j = 1 : D
        x(i, j) = x(i, j)*rand*(Pmax(j) - Pmin);
    end
end
v = rand(N, D)*(Vmax - Vmin) + Vmin;
%% 初始化个体最优位置和最优值
p = x;
pbest = ones(N, 1);
for i = 1 : N
    pbest(i) = func(x(i, :), D);
end
%% 初始化全局位置和最优值
g = ones(1, D);
gbest = eps;
for i = 1 : N
    if (pbest(i) > gbest)
        g = p(i, :);
        gbest = pbest(i);
    end
end
gb = ones(1, T);
%% 迭代
for i = 1 : T
    for j = 1 : N
        %% 更新个体最优位置和最优值
        if (func(x(j, :), D) > pbest(j))
            p(j, :) = x(j, :);
            pbest(j) = func(x(j, :), D);
        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, :));
        x(j, :) = x(j, :) + v(j, :);
        %% 边界条件处理
        for ii = 1 : D
            if (v(j, ii) > Vmax) | (v(j, ii) < Vmin)
                v(j, ii) = rand*(Vmax - Vmin) + Vmin;
            end
            if (x(j, ii) > Pmax(ii)) | (x(j, ii) < Pmin)
                x(j, ii) = rand*(Pmax(ii) - Pmin);
            end
        end
    end
    %% 记录历代全局最优
    gb(i) = gbest;
end

g;
gbest;
figure
plot(gb);
xlabel('迭代次数')
ylabel('适应度值')
title('适应度进化曲线')

适应度函数：

%% 适应度函数
function value = func(x, D)
C = [0.0036 1.6358 22.734;
    0.0042 1.5519 23.150;
    0.0025 1.6728 21.680;
    0.0089 1.3690 17.797;
    0.0065 1.4409 21.150;
    0.0033 1.6290 22.490;
    0.0022 1.5615 24.670;
    0.0022 1.5615 24.670;
    0.0023 1.5499 42.384;
    0.0026 1.3807 46.881]; % 费用系数
value = 0;
for i = 1 : D
    val = C(i, 1)*x(i)^2 + C(i, 2)*x(i) + C(i, 3);
    value = value + val;
end
fmax = 600; % 总负荷需求
f = 0; 
for i = 1: D
    f = f + x(i);
end
if (f > fmax)
    value = value*0.9; % 设置惩罚
end

适应度进化曲线：
在这里插入图片描述

分配结果：

在这里插入图片描述

与文中结果基本相近。

易腐商品最优订货批量与定价

参考文献：田志友, 蒋录全, 吴瑞明. 易腐商品最优订货批量与定价及其粒子群优化解[J]. 系统工程理论与实践, 2005(03): 46-51.

从广义上讲，超过正常销售期后市场价值有明显降低的商品均可归属易腐商品。销售者需要在销售期初，综合考虑市场需求的波动、顾客的消费偏好，以及此类商品的销售期长短，制定合理的采购批量和售价，以确保实现其利润最大化目的。

优化目标

（推导过程省略）

当期初订货量为 s 、商品成本为 c、正常销售期内商品定价为 ω 、销售期过后商品按处理价 φ 进行低价清仓时，期望收入：
$R_s(w)=s(w-c)-(w-\phi)\sum\limits^{s-1}_{m=0}(s-m)P_w(m)$
其中，P_w(m)为正常销售期内，定价为 ω 时，顾客商品需求为 m 的概率分布
$P_w(m)=C_{m+\alpha-1}^m×[\frac {\beta t[1-F(w)]} {1+\beta t[1-F(w)]}]^m ×[\frac 1 {1+\beta t [1-F(w)]}]^\alpha \ ,\ n=0,1,...,\infty;m=0,1,...,\infty$
令其中参数 α = 3，β = 2

F(w)为正态分布函数，令其参数 μ = 10，σ = 1

约束条件

订货量 s 取值范围为[1，20]，定价 ω 取值范围[6, 12]。

代码练习

设商品单位成本 c = 6，超过正常销售期后的处理价 φ = 5。

粒子维数为2维，分别为订货量和商品定价。

由于订货量是整数，故在初始化、更新速度和位置、边界条件处理时，使用round()函数进行取整。

主函数：

clear all;
close all;
clc;
%% 初始化
N = 150; % 粒子个数
D = 2; % 粒子维数
T = 100; % 最大迭代次数
c1 = 1.5; % 学习因子1
c2 = 1.5; % 学习因子2
Wmax = 0.8; % 惯性权重最大值
Wmin = 0.4; % 惯性权重最小值
X1max = 20; % 最大订货量
X1min = 1; % 最小订货量
X2max = 12; % 最大定价
X2min = 6; % 最小定价
V1max = 0.3;
V1min = 0;
V2max = 0.1;
V2min = 0;
%% 初始化种群个体
x = ones(N, D);
v = ones(N, D);
for i = 1 : N
    x(i, 1) = round(rand*(X1max - X1min) + X1min); % 取整
    x(i, 2) = rand*(X2max - X2min) + X2min;
    v(i, 1) = round(rand*(V1max - V1min) + V1min); % 取整
    v(i, 2) = rand*(V2max - V2min) + V2min;
end
%% 初始化个体最优位置和最优值
p = x;
pbest = ones(N, 1);
for i = 1 : N
    pbest(i) = func(x(i, :));
end
%% 初始化全局最优位置和最优值
g = ones(1, D);
gbest = eps;
for i = 1 : N
    if pbest(i) > gbest
        g = p(i, :);
        gbest = pbest(i);
    end
end
gb = ones(1, T);
%% 迭代
for i = 1 : T
    for j = 1 : N
        %% 更新个体最优位置和最优值
        if (func(x(j, :)) > pbest(j))
            p(j, :) = x(j, :);
            pbest(j) = func(x(j, :));
        end
        %% 更新全局最优位置和最优值
        if (pbest(j) > gbest)
            g = p(j, :);
            gbest = pbest(j);
        end
        %% 计算动态惯性权重值
        w = Wmax - (Wmax - Wmin)*i/T;
        %% 更新位置和速度值
        v(j, 1) = round(w*v(j, 1)+c1*rand*(p(j, 1)-x(j, 1))+c2*rand*(g(1)-x(j, 1)));
        v(j, 2) = w*(v(j, 2))+c1*rand*(p(j, 2)-x(j, 2))+c2*rand*(g(2)-x(j, 2));
        x(j, :) = x(j, :) + v(j, :);
        %% 边界条件处理
        if (v(j, 1) > V1max) | (v(j, 1) < V1min)
            v(j, 1) = round(rand*(V1max - V1min) + V1min); 
        end
        if (v(j, 2) > V2max) | (v(j, 2) < V2min)
            v(j, 2) = rand*(V2max - V2min) + V2min;
        end
        if (x(j, 1) > X1max) | (x(j, 1) < X1min)
            x(j, 1) = round(rand*(X1max - X1min) + X1min);
        end
        if (x(j, 2) > X2max) | (x(j, 2) < X2min)
            x(j, 2) = rand*(X2max - X2min) + X2min;
        end
    end
    %% 记录历代全局最优值
    gb(i) = gbest;
end

g;
figure
plot(gb);
xlabel('迭代次数')
ylabel('适应度值')
title('适应度进化曲线')

适应度函数：

%% 适应度函数
function value = func(x)
c = 6; % 单位成本
phi = 5; % 超过正常销售期后的处理价
value = 0;

    function res = p(m, w) % 并愿意以当前定价购买一单位商品的潜在人次的分布概率
    alpha = 3;
    beta = 2; % gamma分布参数
    mu = 10;
    sigma = 1; % 正态分布参数
    C = nchoosek(m + alpha - 1, m); % 组合数
    F = normcdf(w, mu, sigma); % 顾客感知价值的累积分布函数值（正态分布）
    res = C* (beta*(1-F)/(1+beta*(1-F)))^m * (1/(1+beta*(1-F)))^alpha;
    end

for m = 0 : (x(1)-1)
    value = value + (x(1)-m)*p(m,x(2));
end
value = x(1)*(x(2)-c)-(x(2)-phi)*value;
end

适应度进化曲线：

在这里插入图片描述

优化结果：

订货量为7，定价9.1714。与文中结果一致。
在这里插入图片描述

武器-目标分配

参考文献：高尚, 杨静宇.武器-目标分配问题的粒子群优化算法[J].系统工程与电子技术, 2005, 27(7): 1250-1252.

有 n 个目标 T₁，T₂，…，T_n，迎击武器分布于 m 个武器平台 W₁，W₂，…，W_m，第 i（ i = 1，2，…，m）个武器平台最多可使用 r_i 个武器，对目标 T_j 最多可使用 S_j 个武器，武器平台 W_i 迎击目标 T_j 的概率为 p_ij（ i = 1，2，…，m；j = 1，2，…，n），武器最佳分配以分配迎击武器迎击全部目标的失败概率最小为目标。

优化目标

若分配了武器平台 W_i 迎击目标 T_i，则 x_ij = 1，否则 x_ij = 0

分配迎击武器迎击全部目标的失败概率最小：
$\min E = \sum\limits^n_{j=1}\prod^m_{i=1}(1-p_{ij}x_{ij})$

约束条件

（1）第 i（ i = 1，2，…，m）个武器平台最多可使用 r_i 个武器
$\sum\limits^n_{j=1}x_{ij}\leq r_i\ ,\ \ \ \ i=1,2,...,m$
（2）对目标 T_j 最多可使用 S_j 个武器
$\sum\limits^m_{i=1}x_{ij}\leq s_j\ , \ \ \ \ j=1,2,...,n$

代码练习

使用离散粒子群算法。对于约束条件（2），以罚函数形式在适应度函数中体现；对于约束条件（1），在初始化和边界条件处理分别进行约束。

使用数据：

在这里插入图片描述

主函数：

clear all;
close all;
clc;
%% 初始化
N = 500; % 粒子个数
D = 6; % 粒子维数(打击目标数)
m = 4; % 武器平台数
T = 100; % 最大迭代次数
c1 = 1.5; % 学习因子1
c2 = 1.5; % 学习因子2
Wmax = 0.8; % 惯性权重最大值
Wmin = 0.4; % 惯性权重最小值
Vmax = 2; % 速度最大值
Vmin = -2; % 速度最小值
r = [2 1 2 1]; % 每个武器平台最多使用几个武器
%% 初始化种群个体
x = zeros(m, D, N);
for i = 1 : N
    for j = 1 : m
        x(j, :, i) = func1(D, r(j));
    end
end
v = rand(m, D, N)*(Vmax - Vmin) + Vmin;
%% 初始化个体最优位置和最优值
p = x;
pbest = ones(N, 1);
for i = 1 : N
    pbest(i) = func2(x(:, :, i), m, D);
end
%% 初始化全局最优位置和最优值
g = ones(m, D, 1);
gbest = inf;
for i = 1 : N
    if (pbest(i) < gbest)
        g = p(:, :, i);
        gbest = pbest(i);
    end
end
gb = ones(1, T);
%% 迭代
for i = 1 : T
    for j = 1 : N
        %% 更新个体最优位置和最优值
        if (func2(x(:, :, j), m, D) < pbest(j))
            p(:, :, j) = x(:, :, j);
            pbest(j) = func2(x(:, :, j), m, D);
        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 i1 = 1 : m
            for i2 = 1 : D
                if (v(i1, i2, j) > Vmax) | (v(i1, i2, j) < Vmin)
                    v(i1, i2, j) = rand*(Vmax-Vmin) + Vmin;
                end
            end
        end
        %% 更新位置值
        vx(:, :, j) = 1./(1+exp(-v(:, :, j)));
        for j1 = 1 : m
            n = 0; % 计1的个数
            while (n == 0) || (n > r(j1))
                n = 0; 
                for j2 = 1 : D
                    if vx(j1, j2, j) > rand
                        x(j1, j2, j) = 1;
                    else
                        x(j1, j2, j) = 0;
                    end
                    n = n + x(j1, j2, j);
                end
            end
        end
    end
    %% 记录历代全局最优值
    gb(i) = gbest;
end

g;
gbest;
figure
plot(gb);
xlabel('迭代次数')
ylabel('适应度值')
title('适应度进化曲线')

粒子群初始化函数func1：

%% 初始化函数
%% 约束：武器平台最多可使用r个武器
function res = func1(D, r)
res = zeros(1, D);
f = 0; 
i = 1;
while (f == 0) || (f > r) 
    if i > D
        i = 1;
        f = 0;
    end
    res(i) = randi([0, 1]);
    if res(i) == 1
        f = f + 1;
    end
    i = i + 1;
end

适应度函数func2：

%% 适应度函数
function res = func2(x, m, D)
p = [0.5 0.7 0.7 0.8 0.4 0.8;
    0.4 0.8 0.5 0.7 0.6 0.9;
    0.8 0.7 0.7 0.7 0.6 0.5;
    0.5 0.7 0.6 0.8 0.7 0.9]; % 迎击概率表
res = 0;
E = 0;
for i = 1 : D
    e = 1;
    for j = 1 : m
        e = e *(1 - p(j, i)*x(j, i));
    end
    E = E + e;
end
%% 对每个目标最多使用1个武器
afa = 0;
for i = 1 : D
    n = 0;
    for j = 1 : m
        n = n + x(j, i);
    end
    if n > 1
        afa = afa + 1;
    end
end
res = E + 10*afa;

适应度进化曲线：

在这里插入图片描述

优化结果：

在这里插入图片描述

武器1迎击目标4和目标6，武器2迎击目标2，武器3迎击目标1和3，武器4迎击目标5。

迎击全部目标的失败概率为1.4。

与文章结果一致。

除了上面的结果外，还出现了其他优化结果：

在这里插入图片描述

武器1迎击目标3和目标4，武器2迎击目标6，武器3迎击目标1和2，武器4迎击目标5。

在这里插入图片描述

武器1迎击目标2和目标4，武器2迎击目标6，武器3迎击目标1和3，武器4迎击目标5。

迎击全部目标的失败概率均同样为1.4。

工作车间调度

参考文献：谷峰, 陈华平, 卢冰原. 粒子群算法在柔性工作车间调度中的应用[J]. 系统工程, 2005, 23(9): 20-23.

车间配有 m 台机器，要加工 n 种工件，每种工件都由 n_j 个工序组成，用 O_ij 代表工件 J_j 上的第 i 个工序。

调度过程满足以下假设：

（1）同一时刻同一台机器只能加工一个零件；

（2）每个工件在某一时刻只能在一台机器上加工，不能中途中断每一个操作；

（3）同一工件的工序之间有先后约束，不同工件的工序之间没有先后约束；

（4）不同工件具有相同的优先级。

优化目标

假设在可行调度中确定每个工序的开始时间为 s_ij，每道工序的加工时间为 t_ij，使总完工时间 C_max 最小：
$\min C_{\max}=\min \lbrace (s_{ij}+t_{ij})\rbrace$

约束条件

（1）同一工件工序的前后约束

用 t_jk 和 c_jk 分别表示工件 j 在机器 k 上的加工时间和完工时间，如果工件在机器 h 上的加工工序先于机器 k ，则有
$c_{jk}-t_{jk} \geq c_{jh}$
（2）工件的非堵塞约束（保证工序的完整性，不能中断操作）

若工件 i 先于工件 j 到达机器 k，则有
$c_{jk}-c_{ik} \geq t_{jk}$

代码练习

采用如下表的数据：

在这里插入图片描述

初始化时，构建 2*8 矩阵为每一个粒子的结构，第一行 X_process 表示 8 道工序，第二行 X_machine 表示各工序所对应作业的机器，比如
$X_{process}: \ \ \ 1\ 3\ 1\ 2\ 2\ 1\ 3 \ 2 \\X_{machine}: \ \ 1\ 3\ 2\ 1\ 4\ 2\ 4\ 2$
进行粒子位置更新操作时，先进行常规的速度更新计算和位置更新计算，然后对计算后的 X_process 按数值从小到大排序，根据该顺序调换原来 X_process 中的数；对于 X_machine，则采取向上取整操作，并进行边界条件处理。

例如，经过某一次迭代得到结果
$X_{process}: \ \ \ \ \ \ 1\ \ \ \ 3\ \ \ \ 1\ \ \ \ 2\ \ \ \ 2\ \ \ \ 1\ \ \ \ 3\ \ \ \ 2\ \ \ \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 3.2\ 3.9\ 6.4\ 5.8\ 4.3\ 1.5\ 3.6\ 2.4\\X_{machine}: \ \ 1\ \ \ \ 3\ \ \ \ 2\ \ \ \ 1\ \ \ \ 4\ \ \ \ 2\ \ \ \ 4\ \ \ \ 2$
将计算后的 X_process 按数值从小到大排序，1.5（对应初始值1）< 2.4(2) < 3.2(1) < 3.6(3) < 3.9(3) < 4.3(2) < 5.8(2) < 6.4(1)，则对应的迭代后的 X_process 为 [1，2，1，3，3，2，2，1]。

这种编码方式与迭代方式可以解决工件工序的顺序约束。在计算适应度时，按照已有的调度安排，计算总完工时间，也不会产生工序堵塞现象。

但是，这样只能对某一特定的调度情况（几种待加工工件、分别有几道工序、有几台机器可以使用）进行计算，一旦发生了变化，就要重新编码。

主函数：

clear all;
close all;
clc;
%% 初始化
N = 40; % 粒子个数
D = 8; % 粒子维数
T = 200; % 最大迭代个数
c1 = 1.5; % 学习因子1
c2 = 1.5; % 学习因子2
Wmax = 0.8; % 惯性权重最大值
Wmin = 0.4; % 惯性权重最小值
Xmax = 3; % 位置最大值
Xmin = 1; % 位置最小值
Vmax = 1; % 速度最大值
Vmin = -1; % 速度最小值
m = 4; % 机器数
%% 初始化种群个体
x = ones(2, D, N);
for i = 1 : N
    x1 = [1 1 1 2 2 2 3 3];
    randIndex_x1 = randperm(D);
    x1 = x1(randIndex_x1); % 随机排序
    x(1, :, i) = x1;
    x2 = randi([1, m], 1, D);
    x(2, :, i) = x2;
end
v = rand(2, D, N)*(Vmax - Vmin) + Vmin;
%% 初始化最优个体和最优值
p = x;
pbest = ones(N, 1);
for i = 1 : N
    pbest(i) = func(x(:, :, i), D, m);
end
%% 初始化全局最优位置和最优值
g = ones(2, D, 1);
gbest = inf;
for i = 1 : N
    if (pbest(i) < gbest)
        g = p(:, :, i);
        gbest = pbest(i);
    end
end
gb = ones(1, T);
%% 迭代
for i = 1 : T
    for j = 1 : N
        %% 更新个体最优位置和最优值
        if (func(x(:, :, j), D, m) < pbest(j))
            p(:, :, j) = x(:, :, j);
            pbest(j) = func(x(:, :, j), D, m);
        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));
        x1 = x(1, :, j) + v(1, :, j);
        x1_s = sort(x1,'ascend'); % 对x1进行升序排序
        x_tmp = x;
        for k = 1 : D
            col = find(x1 == x1_s(k));
            x(1, k, j) = x_tmp(1, col, j);
        end
        x(2, :, j) = ceil(x(2, :, j) + v(2, :, j)); % 向上取整
        %% 边界条件处理
        for ii = 1 : D
            if (v(1, ii, j) > Vmax) | (v(1, ii, j) < Vmin)
                v(1, ii, j) = rand*(Vmax - Vmin) + Vmin;
            end
            if (v(2, ii, j) > Vmax) | (v(2, ii, j) < Vmin)
                v(2, ii, j) = rand*(Vmax - Vmin) + Vmin;
            end
            if (x(2, ii, j) > m) | (x(2, ii, j) < 1)
                x(2, ii, j) = randi([1, m]);
            end
        end
    end
    %% 记录历代全局最优值
    gb(i) = gbest;
end

g;
gbest;
figure
plot(gb);
xlabel('迭代次数')
ylabel('适应度值')
title('适应度进化曲线')

适应度函数：

%% 适应度函数
function tmax = func(x, D, m)
t(:, :, 1) = [1 4 6 2;
              4 1 1 3;
              3 2 2 3]; % 第1个工件的各工序在各机器上的加工时间
t(:, :, 2) = [1 4 3 3;
              4 8 7 1;
              7 2 3 2]; % 第2个工件的各工序在各机器上的加工时间
t(:, :, 3) = [3 2 2 7;
              9 3 6 1;
              0 0 0 0]; % 第3个工件的各工序在各机器上的加工时间
time = zeros(1, m); % 4台机器的工作时间
n = ones(1, 3); % 3个工件在下面循环中为第几道工序
for ii = 1 : D
    f = x(1, ii); % 加工的是哪个工件
    mm = x(2, ii); % 在哪个机器上加工
    nn = n(f); % 该工件的第几道工序
    time(mm) = time(mm) + t(nn, mm, f);
    n(f) = n(f) + 1;
end
tmax = max(time);

适应度进化曲线：
在这里插入图片描述

优化结果：
在这里插入图片描述

机器1：O₂₁ → O₁₁

机器2：O₃₁ → O₂₃

机器3：O₁₂

机器4：O₁₁ → O₂₂ → O₃₂

总完工时间：4min

还有其他几种优化结果，总完工时间均为4min。

鸢尾花分类规则提取

参考文献：高亮, 高海兵, 周驰, 等. 基于粒子群优化算法的模式分类规则获取[J]. 华中科技大学学报(自然科学版), 2004, 32(11): 24-26.

规则提取：

令 x_il, x_iu 分别表示规则中的分类变量 x_i 的下界和上界。设存在 k 个分类变量，即 x = [x₁, x₂, …, x_k]，规则的一般形式可表示如下：

If x₁ ∈ [x_1l, x_1u] and x₂ ∈ [x_2l, x_2u] and … and x_k ∈ [x_kl, x_ku], then x ∈ c_i.

优化目标

确定每个粒子的规则正确率指标 α，β。其中 α 为满足规则条件部分的当前类的数目与当前类的总数目的比例；β 为满足规则条件部分的非当前类的数目与当前类总数目的比例。

检验规则的正确率指标：若 α >= α₀ 且 β <= β₀，则该规则满足正确率指标。

对于满足正确率指标的粒子，以规则中的分类变量的范围 σ 为目标函数。

总结如下：
$\ α \geq \alpha_0 \ AND\ \beta \leq \beta_0 \ \ \ THEN \\ \sigma = \sum\limits^k_{i=1}(x_{iu}-x_{il})/k \\ ELSE\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ f=2-(\alpha-\beta) \\ END \ \ IF\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$

约束条件

观察鸢尾花数据的特点，初始化时初始为 [0, 5] 的随机数。

代码练习

将每个粒子构造为 3×2×4 的数组，3表示3种鸢尾花，2表示分类变量的上界和下界，4表示鸢尾花的4种属性。

在进行边界处理时，如果粒子的位置超出了一定范围（代码中设定为 [0, 10]），或下界超过了上界，都对该粒子的位置重新初始化。

Iris数据集共150条，代码中选取了三种鸢尾花的各前30条，共90条数据，参与适应度函数的计算。

主函数：

clear all;
close all;
clc;
%% 初始化
N = 100; % 粒子个数
D = 4; % 粒子维数
T = 200; % 最大迭代次数
c1 = 1.5; % 学习因子1
c2 = 1.5; % 学习因子2
Wmax = 0.9; % 惯性权重最大值
Wmin = 0.4; % 惯性权重最小值
Vmax = 2; % 速度最大值
Vmin = -2; % 速度最小值
%% 初始化种群个体
x = ones(3, 2, D, N);
x(3, 1, :, :) = x(3, 1, :, :)*rand*5;
x(3, 2, :, :) = x(3, 2, :, :)*rand*5 + 5;
v = rand(3, 2, D, N)*(Vmax - Vmin) + Vmin;
%% 初始化个体最优位置和最优值
p = x;
pbest = ones(N, 1);
for i = 1 : N
    pbest(i) = func(x(:, :, :, i));
end
%% 初始化全局最优位置和最优值
g = ones(3, 2, D, 1);
gbest = inf;
for i = 1 : N
    if(pbest(i) < gbest)
        g = p(:, :, :, i);
        gbest = pbest(i);
    end
end
gb = ones(1, T);
%% 迭代
for i = 1 : T
    for j = 1 : N
        %% 更新个体最优位置和最优值
        if func(x(:, :, :, j) < pbest(j))
            p(:, :, :, j) = x(:, :, :, j);
            pbest(j) = func(x(:, :, :, j));
        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));
        x(:,:,:,j) = x(:,:,:,j) +v(:,:,:,j);
        %% 边界条件处理
        for ii = 1 : D
            for m = 1 : 3
                for k = 1 : 2
                    if (v(m, k, ii, j) > Vmax) | (v(m, 1, ii, j)) < Vmin
                        v(m, k, ii, j) = rand*(Vmax - Vmin) + Vmin;
                    end
                    if (x(m, k, ii, j) > 10) | (x(m, k, ii, j) < 0)
                        x(m, k, ii, j) = rand*10;
                    end
                end
                if (x(m, 1, ii, j) >= x(m, 2, ii, j))
                    x(m, 1, ii, j) = rand*5;
                    x(m, 2, ii, j) = 5 + rand*5;
                end
            end
        end
    end
    %% 记录历代全局最优值
    gb(i) = gbest;
end
g;
figure
plot(gb);
xlabel('迭代次数')
ylabel('适应度值')
title('适应度进化曲线')

适应度函数：

%% 适应度函数
function value = func(x)
C = [5.1 3.5 1.4 0.2; % Iris-setosa
4.9 3.0 1.4 0.2;
4.7 3.2 1.3 0.2;
4.6 3.1 1.5 0.2;
5.0 3.6 1.4 0.2;
5.4 3.9 1.7 0.4;
4.6 3.4 1.4 0.3;
5.0 3.4 1.5 0.2;
4.4 2.9 1.4 0.2;
4.9 3.1 1.5 0.1;
5.4 3.7 1.5 0.2;
4.8 3.4 1.6 0.2;
4.8 3.0 1.4 0.1;
4.3 3.0 1.1 0.1;
5.8 4.0 1.2 0.2;
5.7 4.4 1.5 0.4;
5.4 3.9 1.3 0.4;
5.1 3.5 1.4 0.3;
5.7 3.8 1.7 0.3;
5.1 3.8 1.5 0.3;
5.4 3.4 1.7 0.2;
5.1 3.7 1.5 0.4;
4.6 3.6 1.0 0.2;
5.1 3.3 1.7 0.5;
4.8 3.4 1.9 0.2;
5.0 3.0 1.6 0.2;
5.0 3.4 1.6 0.4;
5.2 3.5 1.5 0.2;
5.2 3.4 1.4 0.2;
4.7 3.2 1.6 0.2;
7.0 3.2 4.7 1.4; % Iris-versicolor
6.4 3.2 4.5 1.5;
6.9 3.1 4.9 1.5;
5.5 2.3 4.0 1.3;
6.5 2.8 4.6 1.5;
5.7 2.8 4.5 1.3;
6.3 3.3 4.7 1.6;
4.9 2.4 3.3 1.0;
6.6 2.9 4.6 1.3;
5.2 2.7 3.9 1.4;
5.0 2.0 3.5 1.0;
5.9 3.0 4.2 1.5;
6.0 2.2 4.0 1.0;
6.1 2.9 4.7 1.4;
5.6 2.9 3.6 1.3;
6.7 3.1 4.4 1.4;
5.6 3.0 4.5 1.5;
5.8 2.7 4.1 1.0;
6.2 2.2 4.5 1.5;
5.6 2.5 3.9 1.1;
5.9 3.2 4.8 1.8;
6.1 2.8 4.0 1.3;
6.3 2.5 4.9 1.5;
6.1 2.8 4.7 1.2;
6.4 2.9 4.3 1.3;
6.6 3.0 4.4 1.4;
6.8 2.8 4.8 1.4;
6.7 3.0 5.0 1.7;
6.0 2.9 4.5 1.5;
5.7 2.6 3.5 1.0;
6.3 3.3 6.0 2.5; % Iris-virginica
5.8 2.7 5.1 1.9;
7.1 3.0 5.9 2.1;
6.3 2.9 5.6 1.8;
6.5 3.0 5.8 2.2;
7.6 3.0 6.6 2.1;
4.9 2.5 4.5 1.7;
7.3 2.9 6.3 1.8;
6.7 2.5 5.8 1.8;
7.2 3.6 6.1 2.5;
6.5 3.2 5.1 2.0;
6.4 2.7 5.3 1.9;
6.8 3.0 5.5 2.1;
5.7 2.5 5.0 2.0;
5.8 2.8 5.1 2.4;
6.4 3.2 5.3 2.3;
6.5 3.0 5.5 1.8;
7.7 3.8 6.7 2.2;
7.7 2.6 6.9 2.3;
6.0 2.2 5.0 1.5;
6.9 3.2 5.7 2.3;
5.6 2.8 4.9 2.0;
7.7 2.8 6.7 2.0;
6.3 2.7 4.9 1.8;
6.7 3.3 5.7 2.1;
7.2 3.2 6.0 1.8;
6.2 2.8 4.8 1.8;
6.1 3.0 4.9 1.8;
6.4 2.8 5.6 2.1;
7.2 3.0 5.8 1.6]; % 鸢尾花数据集
value = 0;
a0 = 0.99;
b0 = 0.01; % 正确率指标阈值
n = 0; % 满足正确率指标的样本个数
for i = 1 : 90
    flag = 0; 
    m = 0; % 表示哪一种鸢尾花
    if i <= 30
        m = 1;
    elseif (i > 30) & (i <= 60)
        m = 2;
    else 
        m = 3;
    end
    for j = 1 : 4 % 对该样本的四种属性依次进行比较
        if (abs(x(m, 1, j)-C(i,j)) <= 1.5) & (abs(x(m, 2, j)-C(i,j)) <= 1.5)
            flag = 1;
        else
            flag = 0;
        end
    end
    if flag == 1 % 说明该样本满足该粒子表示的分类准则
        n = n + 1;
    end
end
a = n/90;
b = 1 - n/90; % 正确率指标
if (a >= a0) & (b <= b0)
    for m = 1 : 3
        for i = 1 : 4
            value = value + x(m, 2, i) - x(m, 1, i);
        end
    end
    value = value/12;
else
    value = 6 - (a - b);
end