复现了 《混合差分进化与调度算法》-王凌 中 求解置换流水车间调度问题(PFSP)的混合差分进化算法。

算法基本流程如下：

 局部搜索的基本流程如下：

基于LOV规则进行编码： 

 局部搜索后基于LOV规则调整编码：

算法的c++源代码如下：

DE.h:

#pragma once
#include <vector>
using namespace std;

class DE
{
private:
    int DE_Pop_Size;  //种群规模，Pop_Size ≥ 4
    int DE_Max_Iter;  //最大迭代次数
    int DE_N;  //解的维度(工件数)
    int DE_M;  //加工机器数
    float DE_F;  //变异操作中的缩放比例因子
    float DE_CR;  //交叉操作中的接收交叉概率
    float DE_Lower;  //取值下界
    float DE_Upper;  //取值上界
    vector<vector<float>> DE_Population;  //DE种群
    vector<vector<float>> DE_Temp_Population;  //临时种群，保存上一代信息
    vector<vector<int>> DE_Solution;  //解码后的DE种群(工件排序)
    vector<double> DE_Makespan_Set;  //种群中各个体的完工时间
    vector<float> DE_Bestit;  //种群中的最优个体
    vector<int> DE_Best_Solution;  //最优个体对应的排序
    int DE_Bestit_Pos;  //最优个体在种群中的位置
    double DE_Min_Makespan;  //最优个体的完工时间
    vector<vector<double>> DE_T;  //加工时间矩阵(各工件在各机器上的加工时间)  

public:
    DE(const vector<vector<double>>& process_time, int population_size, int max_iteration, float F, float CR);  //构造函数
    void optimize();  //调用DE算法进行优化
    void init();  //初始化
    void mutation_crossover_selection();  //变异、交叉、选择
    void local_search();  //局部搜索
    vector<vector<double>> calculate_makespan(const vector<int>& solution);  //计算完工时间矩阵，完工时间 makespan = schedule[N-1][M-1]
    vector<int> LOV_decode(const vector<float>& DE_Individual);  //根据LOV规则进行解码
    void LOV_encode(const vector<int>& solution, vector<float>& DE_Individual); //根据LOV规则调整编码
    vector<int> insert_operation(vector<int> solution, int insertPos, int selectPos);  //插入操作
    vector<int> get_best_solution();  //输出最优工件排序
    double get_min_makespan();  //输出最小完工时间
};

DE.cpp:

#include "DE.h"
#include <vector>
#include <numeric>
#include <algorithm>
#include <unordered_map>
#include <iostream>

using namespace std;

//构造函数
DE::DE(const vector<vector<double>>& process_time, int population_size, int max_iteration, float F, float CR) :
    DE_T(process_time), DE_Pop_Size(population_size), DE_Max_Iter(max_iteration), 
    DE_N(process_time.size()), DE_M(process_time[0].size()), 
    DE_F(F), DE_CR(CR), DE_Lower((float)0.000), DE_Upper((float)4.000), 
    DE_Bestit_Pos(0), DE_Min_Makespan(DBL_MAX)
{
    DE_Population.resize(DE_Pop_Size, vector<float>(DE_N));
    DE_Temp_Population.resize(DE_Pop_Size, vector<float>(DE_N));
    DE_Solution.resize(DE_Pop_Size, vector<int>(DE_N));
    DE_Makespan_Set.resize(population_size,DBL_MAX);
    DE_Bestit.reserve(DE_N);
    DE_Best_Solution.reserve(DE_N);
}

//调用DE算法进行优化
void DE::optimize() {
    //1.初始化
    init();

    //2.迭代
    for (int i = 0; i < DE_Max_Iter; i++)
    {
        //变异、交叉、选择
        mutation_crossover_selection();
        //局部搜索
        local_search();
    }
}

//初始化
void DE::init() {    
    //随机生成DE初始种群个体
    float random;
    for (int i = 0; i < DE_Pop_Size; i++)
    {
        for (int j = 0; j < DE_N; j++)
        {
            random = rand() % 1000 / (float)1000;
            DE_Population[i][j] = DE_Lower + random * (DE_Upper - DE_Lower);
        }
    }
    
    //初始种群个体解码，得到工件排序
    for (int i = 0; i < DE_Pop_Size; i++)
    {
        DE_Solution[i] = LOV_decode(DE_Population[i]);
    }

    //计算初始种群适应度
    for (int i = 0; i < DE_Pop_Size; i++)
    {
        DE_Makespan_Set[i] = calculate_makespan(DE_Solution[i])[DE_N - 1][DE_M - 1];
    }

    //寻找初始种群中的最优个体并记录
    vector<double>::iterator itMin = min_element(DE_Makespan_Set.begin(), DE_Makespan_Set.end());
    DE_Bestit_Pos = distance(DE_Makespan_Set.begin(), itMin);
    DE_Bestit = DE_Population[DE_Bestit_Pos];
    DE_Best_Solution = DE_Solution[DE_Bestit_Pos];
    DE_Min_Makespan = DE_Makespan_Set[DE_Bestit_Pos];
}

//变异、交叉、选择
void DE::mutation_crossover_selection() {
    DE_Temp_Population.assign(DE_Population.begin(), DE_Population.end());  //存放上一代种群
    vector<int> Temp(DE_N);  //临时解(工件排序)

    for (int i = 0; i < DE_Pop_Size; i++)
    {
        //变异和交叉
        int r1 = rand() % DE_Pop_Size, r2 = rand() % DE_Pop_Size;
        while (r1 == r2 || r1 == i || r2 == i) {
            r1 = rand() % DE_Pop_Size;
            r2 = rand() % DE_Pop_Size;
        }
        int L = 0, k = rand() % DE_N;
        do
        {
            DE_Population[i][k] = DE_Temp_Population[i][k] +
                DE_F * (DE_Bestit[k] - DE_Population[i][k]) +
                DE_F * (DE_Temp_Population[r1][k] - DE_Population[r2][k]);
            if (DE_Population[i][k] < DE_Lower) {
                DE_Population[i][k] = 2 * DE_Lower - DE_Population[i][k];
            }
            else if (DE_Population[i][k] > DE_Upper) {
                DE_Population[i][k] = 2 * DE_Upper - DE_Population[i][k];
            }
            k = (k + 1) % DE_N;
            L++;
        } while ((rand() % 1000 / (float)1000) < DE_CR && L < DE_N);

        //选择
        Temp = LOV_decode(DE_Population[i]);  //当前个体对应的工件排序
        double Temp_Makespan = calculate_makespan(Temp)[DE_N - 1][DE_M - 1];  //当前个体的完工时间
        if (Temp_Makespan < DE_Makespan_Set[i]) {
            DE_Solution[i] = Temp;
            DE_Makespan_Set[i] = Temp_Makespan;
        }
        else {
            DE_Population[i] = DE_Temp_Population[i];
        }
        if (Temp_Makespan < DE_Min_Makespan) {
            //更新最优个体
            DE_Bestit = DE_Population[i];  
            DE_Best_Solution = DE_Solution[i];
            DE_Min_Makespan = Temp_Makespan;
            DE_Bestit_Pos = i;
        }
    }
}

//对种群中的最优个体进行局部搜索
void DE::local_search() {
    int u = rand() % DE_N, v = rand() % DE_N;
    while (u == v) {
        v = rand() % DE_N;
    }
    vector<int> Temp = insert_operation(DE_Best_Solution, u, v);
    vector<int> Temp1(DE_N);
    double Makespan_Temp = calculate_makespan(Temp)[DE_N - 1][DE_M - 1];

    for (int i = 1; i < DE_N * (DE_N - 1); i++)
    {
        u = rand() % DE_N, v = rand() % DE_N;
        while (u == v) {
            v = rand() % DE_N;
        }
        Temp1 = insert_operation(DE_Best_Solution, u, v);        
        double Makespan_Temp1 = calculate_makespan(Temp1)[DE_N - 1][DE_M - 1];
        if (Makespan_Temp1 < Makespan_Temp) {
            Temp = Temp1;
            Makespan_Temp = Makespan_Temp1;
        }
    }

    //如果局部搜索产生的新解比原来好，则接受
    if (Makespan_Temp < DE_Min_Makespan) {        
        DE_Best_Solution = Temp;
        DE_Min_Makespan = Makespan_Temp;
        LOV_encode(DE_Best_Solution, DE_Bestit);
        //在种群中更新局部搜索后的更优个体
        DE_Population[DE_Bestit_Pos] = DE_Bestit;
        DE_Solution[DE_Bestit_Pos] = DE_Best_Solution;
        DE_Makespan_Set[DE_Bestit_Pos] = DE_Min_Makespan;
    }
}

//计算完工时间矩阵，完工时间 makespan = schedule[N-1][M-1]
vector<vector<double>> DE::calculate_makespan(const vector<int>& solution) {
    vector<vector<double>> schedule(DE_N, vector<double>(DE_M, 0));  //调度矩阵，存放各机器上的调度结果
    //计算第1个工件：sche(J1,1) = T(J1,1); sche(J1,k) = sche(J1,k-1) + T(J1,k)
    schedule[0][0] = DE_T.at(solution[0] - 1).at(0);   //T[solution[0]-1][0];
    for (int j = 1; j < DE_M; j++)
    {
        schedule[0][j] = schedule[0][j - 1] + DE_T.at(solution[0] - 1).at(j);  //T[solution[0]-1][j]
    }
    //计算第2至n个工件：sche(Ji,1) = sche(Ji-1,1) + T(Ji,1); sche(Ji,k) = max{sche(Ji-1,k), sche(Ji,k-1)} + T(J1,k)
    for (int i = 1; i < DE_N; i++)
    {
        schedule[i][0] = schedule[i - 1][0] + DE_T.at(solution[i] - 1).at(0);  //T[solution[i]-1][0]
        for (int j = 1; j < DE_M; j++)
        {
            schedule[i][j] = max(schedule[i - 1][j], schedule[i][j - 1]) + DE_T.at(solution[i] - 1).at(j);  //T[solution[i]-1][j]
        }
    }
    return schedule;
}

//根据LOV规则进行解码（旧，当个体存在相同分量时，由于哈希冲突，所得排列会出错）
//例：{1.36, 3.85, *2.55, 0.63, *2.55, 0.82} -> {2, *1, 5, *1, 6, 4}
//vector<int> DE::LOV_decode(const vector<float>& DE_Individual) {
//    int n = DE_Individual.size();
//    vector<int> solution(n);
//    vector<float> temp(DE_Individual);  //创建个体副本
//    sort(temp.rbegin(), temp.rend());  //对副本进行降序排列
//    unordered_map<float, int> record;  //记录temp[i]对应的下标
//
//    //求解中间序列φ
//    for (int i = 0; i < n; i++)
//    {
//        record[temp[i]] = i + 1;
//    }
//    for (int i = 0; i < n; i++)
//    {
//        solution[i] = record[DE_Individual[i]];
//    }
//
//    //求解最终加工次数π，π(φ(i)) = i
//    vector<int> temp_solution = solution;
//    for (int i = 0; i < n; i++)
//    {
//        solution[temp_solution[i] - 1] = i + 1;;
//    }
//    return solution;
//}

//根据LOV规则进行解码（新）
//例：{1.36, 3.85, 2.55, 0.63, 2.68, 0.82} -> {4, 1, 3, 6, 2, 5} -> {2, 5, 3, 1, 6, 4}
vector<int> DE::LOV_decode(const vector<float>& DE_Individual) {
    int n = DE_Individual.size();
    vector<pair<float, int>> temp(n);
    for (int i = 0; i < n; i++) {
        temp[i].first = DE_Individual[i];  //first保存个体信息
        temp[i].second = i + 1;  //second保存对应位置
    }
    sort(temp.begin(), temp.end(), [&](auto a, auto b) {
        return a.first > b.first;  //根据first降序排列
        });
    vector<int> solution(n);
    for (int i = 0; i < n; i++) {
        solution[i] = temp[i].second;  //排序后的second即所求工件排列
    }
    return solution;
}

//根据LOV规则调整编码
//例：{2, 5, 1, 3, 6, 4} -> {3, 1, 4, 6, 2, 5} -> {2.55, 3.85, 1.36, 0.63, 2.68, 0.82}
void DE::LOV_encode(const vector<int>& solution, vector<float>& DE_Individual) {
    int n = DE_Individual.size();
    vector<int> temp(n);  //创建中间序列
    unordered_map<int, float> record;
    //求解中间序列
    for (int i = 0; i < n; i++)
    {
        temp[solution[i] - 1] = i + 1;
    }
    //重新调整个体编码
    sort(DE_Individual.rbegin(), DE_Individual.rend());
    for (int i = 0; i < n; i++)
    {
        record[i + 1] = DE_Individual[i];
    }
    for (int i = 0; i < n; i++)
    {
        DE_Individual[i] = record[temp[i]];
    }
}

//插入操作 (需满足 insertPos != selectPos 且两者均 ≥ 0，均 ＜ solution.size)
//例： insert_operation({1, 2, 3, 4}, 0, 2) ---> {3, 1, 2, 4}
vector<int> DE::insert_operation(vector<int> solution, int insertPos, int selectPos) {
    if (insertPos > selectPos) {
        for (int i = selectPos; i < insertPos; i++) {
            swap(solution[i], solution[i + 1]);
        }
    }
    else {
        for (int i = selectPos; i > insertPos; i--) {
            swap(solution[i], solution[i - 1]);
        }
    }
    return solution;
}

//输出最优工件排序
vector<int> DE::get_best_solution() {
    return DE_Best_Solution;
}

//输出最小完工时间
double DE::get_min_makespan() {
    return DE_Min_Makespan;
}

main.cpp:

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <numeric>
#include <chrono>
#include <unordered_map>
#include "DE.h"

using namespace std;

extern void printVect(const vector<int>& v);
extern void printVect(const vector<vector<double>>& v);

int main() {
    srand((unsigned)time(NULL));
    chrono::steady_clock::time_point t1 = chrono::steady_clock::now();

    //**测试代码**//
    //初始化测试案例(Taillard-PFSP, ta001)，设定DE算法参数
    vector<vector<double>> matrixT(20, vector<double>(5, 0));
    matrixT[0] = { 54, 79, 16, 66, 58 };
    matrixT[1] = { 83, 3, 89, 58, 56 };
    matrixT[2] = { 15, 11, 49, 31, 20 };
    matrixT[3] = { 71, 99, 15, 68, 85 };
    matrixT[4] = { 77, 56, 89, 78, 53 };
    matrixT[5] = { 36, 70, 45, 91, 35 };
    matrixT[6] = { 53, 99, 60, 13, 53 };
    matrixT[7] = { 38, 60, 23, 59, 41 };
    matrixT[8] = { 27, 5, 57, 49, 69 };
    matrixT[9] = { 87, 56, 64, 85, 13 };
    matrixT[10] = { 76, 3, 7, 85, 86 };
    matrixT[11] = { 91, 61, 1, 9, 72 };
    matrixT[12] = { 14, 73, 63, 39, 8 };
    matrixT[13] = { 29, 75, 41, 41, 49 };
    matrixT[14] = { 12 ,47 ,63 ,56 ,47 };
    matrixT[15] = { 77 ,14 ,47 ,40 ,87 };
    matrixT[16] = { 32 ,21 ,26 ,54 ,58 };
    matrixT[17] = { 87 ,86 ,75 ,77 ,18 };
    matrixT[18] = { 68 , 5 ,77 ,51 ,68 };
    matrixT[19] = { 94 ,77 ,40 ,31 ,28 };   
    int population = 50, max_iteration = 100;
    float F = 0.5, CR = 0.1;
    
    //调用DE算法进行优化求解
    DE* de = new DE(matrixT, population, max_iteration, F, CR);
    de->optimize();
    
    //打印输入输出结果
    cout << "**** 输入 ****" << endl;
    cout << "工件数：" << matrixT.size() << "  机器数：" << matrixT[0].size() << endl << endl << "加工时间：" << endl;
    printVect(matrixT);
    cout << "算法参数：" << endl << "种群规模：" << population << "  迭代次数：" << max_iteration;
    cout << "  F：" << F << "  CR：" << CR << endl << endl;
    cout << "**** 输出 ****" << endl << "最优工件排序：";
    printVect(de->get_best_solution());
    cout << "对应完工时间：" << de->get_min_makespan() << endl;

    //计算程序运行时间
    chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
    //秒
    double dr_s = std::chrono::duration<double>(t2 - t1).count();
    //毫秒级
    double dr_ms = std::chrono::duration<double, std::milli>(t2 - t1).count();
    //微妙级
    //double dr_us = std::chrono::duration<double, std::micro>(t2 - t1).count();
    //纳秒级
    //double dr_ns = std::chrono::duration<double, std::nano>(t2 - t1).count();
    cout << endl << "程序运行时间：" << dr_ms << " ms" << endl;

    system("pause");
    return 0;
}

//打印数组
void printVect(const vector<int>& v) {
    for (size_t i = 0; i < v.size(); i++)
    {
        cout << v[i] << " ";
    }
    cout << endl;
}

//打印数组
void printVect(const vector<vector<double>>& v) {
    for (size_t i = 0; i < v.size(); i++)
    {
        for (size_t j = 0; j < v[0].size(); j++)
        {
            cout << v[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}

测试结果：

参考文献：

[1]    王凌等. 混合差分进化与调度算法[M]. 清华大学出版社, 2012.