带精英策略的非支配排序遗传算法（NSGA-II）C++实现

带精英策略的非支配排序遗传算法（NSGA-II）C++实现

一、算法原理与实现思路

1.1 NSGA-II核心思想

NSGA-II（Non-dominated Sorting Genetic Algorithm II）是一种经典的多目标优化算法，通过非支配排序和拥挤度计算实现种群进化，结合精英保留策略提高收敛性。

核心步骤：

  1. 非支配排序：将种群划分为不同前沿等级
  2. 拥挤度计算：保持种群多样性
  3. 二元锦标赛选择：基于前沿等级和拥挤度
  4. 模拟二进制交叉（SBX）和多项式变异
  5. 精英保留：合并父代与子代种群

1.2 C++实现优势

二、完整C++实现代码

2.1 头文件与数据结构

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <random>
#include <fstream>
#include <chrono>
#include <omp.h>

using namespace std;

// 个体结构体
struct Individual {
    vector<double> variables;  // 决策变量
    vector<double> objectives; // 目标函数值
    vector<int> dominated;     // 被支配的个体索引
    int dominate_count;        // 支配其他个体的数量
    int rank;                 // 非支配前沿等级
    double crowding_distance;  // 拥挤度距离
    
    Individual(int n_var, int n_obj) 
        : variables(n_var), objectives(n_obj), 
          dominated(vector<int>()), dominate_count(0), 
          rank(0), crowding_distance(0.0) {}
};

// 算法参数
struct Parameters {
    int pop_size = 100;        // 种群大小
    int max_generations = 500;  // 最大迭代次数
    double crossover_prob = 0.9;// 交叉概率
    double mutation_prob = 0.1; // 变异概率
    double eta_c = 15.0;        // 交叉分布指数
    double eta_m = 20.0;        // 变异分布指数
    int n_var = 10;             // 变量个数
    int n_obj = 2;              // 目标函数个数
    vector<double> var_min;     // 变量下界
    vector<double> var_max;     // 变量上界
};

2.2 核心算法实现

class NSGA2 {
public:
    NSGA2(Parameters params) : params(params) {
        // 初始化变量边界
        if (params.var_min.empty()) {
            params.var_min = vector<double>(params.n_var, 0.0);
        }
        if (params.var_max.empty()) {
            params.var_max = vector<double>(params.n_var, 1.0);
        }
        
        // 初始化随机数生成器
        rng = mt19937(chrono::steady_clock::now().time_since_epoch().count());
    }
    
    void run() {
        // 初始化种群
        initializePopulation();
        
        // 主循环
        for (int gen = 0; gen < params.max_generations; gen++) {
            // 评估种群
            evaluatePopulation();
            
            // 非支配排序
            vector<vector<int>> fronts = nonDominatedSort(population);
            
            // 计算拥挤度
            calculateCrowdingDistance(fronts);
            
            // 选择父代
            vector<Individual> parents = selection(population, fronts);
            
            // 交叉和变异
            vector<Individual> offspring = generateOffspring(parents);
            
            // 合并种群
            vector<Individual> combined = population;
            combined.insert(combined.end(), offspring.begin(), offspring.end());
            
            // 环境选择
            population = environmentalSelection(combined, params.pop_size);
            
            // 输出进度
            if (gen % 50 == 0) {
                cout << "Generation " << gen 
                     << ", First Front Size: " << fronts[0].size() << endl;
            }
        }
        
        // 最终结果
        evaluatePopulation();
        vector<vector<int>> final_fronts = nonDominatedSort(population);
        pareto_front = extractParetoFront(final_fronts[0]);
    }
    
    // 获取Pareto前沿
    vector<Individual> getParetoFront() const {
        return pareto_front;
    }

private:
    Parameters params;
    vector<Individual> population;
    vector<Individual> pareto_front;
    mt19937 rng;
    
    // 初始化种群
    void initializePopulation() {
        population.clear();
        uniform_real_distribution<double> dist(0.0, 1.0);
        
        for (int i = 0; i < params.pop_size; i++) {
            Individual ind(params.n_var, params.n_obj);
            for (int j = 0; j < params.n_var; j++) {
                double r = dist(rng);
                ind.variables[j] = params.var_min[j] + 
                                 r * (params.var_max[j] - params.var_min[j]);
            }
            population.push_back(ind);
        }
    }
    
    // 评估种群（示例函数，需替换为实际问题）
    void evaluatePopulation() {
        #pragma omp parallel for
        for (int i = 0; i < population.size(); i++) {
            // 示例目标函数：ZDT1
            double f1 = population[i].variables[0];
            double sum = 0.0;
            for (int j = 1; j < params.n_var; j++) {
                sum += population[i].variables[j];
            }
            double g = 1.0 + 9.0 * sum / (params.n_var - 1);
            double f2 = g * (1.0 - sqrt(f1 / g));
            
            population[i].objectives[0] = f1;
            population[i].objectives[1] = f2;
        }
    }
    
    // 非支配排序
    vector<vector<int>> nonDominatedSort(vector<Individual>& pop) {
        int n = pop.size();
        vector<vector<int>> fronts;
        vector<int> dom_count(n, 0);
        vector<vector<int>> dominated(n);
        
        // 计算支配关系
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                if (i == j) continue;
                
                if (dominates(pop[i], pop[j])) {
                    dominated[i].push_back(j);
                } else if (dominates(pop[j], pop[i])) {
                    dom_count[i]++;
                }
            }
            
            if (dom_count[i] == 0) {
                pop[i].rank = 0;
                if (fronts.empty()) fronts.push_back(vector<int>());
                fronts[0].push_back(i);
            }
        }
        
        // 构建前沿
        int front_index = 0;
        while (!fronts[front_index].empty()) {
            vector<int> next_front;
            for (int i : fronts[front_index]) {
                for (int j : dominated[i]) {
                    dom_count[j]--;
                    if (dom_count[j] == 0) {
                        pop[j].rank = front_index + 1;
                        next_front.push_back(j);
                    }
                }
            }
            front_index++;
            if (!next_front.empty()) {
                fronts.push_back(next_front);
            }
        }
        
        return fronts;
    }
    
    // 判断个体i是否支配个体j
    bool dominates(const Individual& a, const Individual& b) {
        bool better = false;
        for (int i = 0; i < params.n_obj; i++) {
            if (a.objectives[i] > b.objectives[i]) {
                return false;
            }
            if (a.objectives[i] < b.objectives[i]) {
                better = true;
            }
        }
        return better;
    }
    
    // 计算拥挤度
    void calculateCrowdingDistance(const vector<vector<int>>& fronts) {
        for (const auto& front : fronts) {
            int size = front.size();
            if (size == 0) continue;
            
            // 初始化拥挤度
            for (int idx : front) {
                population[idx].crowding_distance = 0.0;
            }
            
            // 对每个目标函数计算
            for (int obj_idx = 0; obj_idx < params.n_obj; obj_idx++) {
                // 按目标函数值排序
                vector<int> sorted = front;
                sort(sorted.begin(), sorted.end(), 
                    int a, int b { 
                        return population[a].objectives[obj_idx] < 
                               population[b].objectives[obj_idx]; 
                    });
                
                // 设置边界点的拥挤度为无穷大
                population[sorted[0]].crowding_distance = numeric_limits<double>::max();
                population[sorted.back()].crowding_distance = numeric_limits<double>::max();
                
                // 计算中间点的拥挤度
                double min_obj = population[sorted[0]].objectives[obj_idx];
                double max_obj = population[sorted.back()].objectives[obj_idx];
                double range = max_obj - min_obj;
                if (range < 1e-10) continue;
                
                for (int i = 1; i < size - 1; i++) {
                    double dist = population[sorted[i+1]].objectives[obj_idx] - 
                                 population[sorted[i-1]].objectives[obj_idx];
                    population[sorted[i]].crowding_distance += dist / range;
                }
            }
        }
    }
    
    // 选择操作（二元锦标赛）
    vector<Individual> selection(const vector<Individual>& pop, 
                                const vector<vector<int>>& fronts) {
        vector<Individual> parents;
        uniform_int_distribution<int> dist(0, pop.size() - 1);
        
        while (parents.size() < pop.size()) {
            int idx1 = dist(rng);
            int idx2 = dist(rng);
            const Individual& a = pop[idx1];
            const Individual& b = pop[idx2];
            
            // 选择规则：优先选择前沿等级低的，相同则选择拥挤度大的
            if (a.rank < b.rank || 
                (a.rank == b.rank && a.crowding_distance > b.crowding_distance)) {
                parents.push_back(a);
            } else {
                parents.push_back(b);
            }
        }
        
        return parents;
    }
    
    // 生成子代
    vector<Individual> generateOffspring(const vector<Individual>& parents) {
        vector<Individual> offspring;
        uniform_real_distribution<double> prob_dist(0.0, 1.0);
        
        for (int i = 0; i < params.pop_size; i += 2) {
            int idx1 = i % parents.size();
            int idx2 = (i + 1) % parents.size();
            const Individual& p1 = parents[idx1];
            const Individual& p2 = parents[idx2];
            
            // 交叉
            Individual c1(params.n_var, params.n_obj);
            Individual c2(params.n_var, params.n_obj);
            
            if (prob_dist(rng) < params.crossover_prob) {
                crossover(p1, p2, c1, c2);
            } else {
                c1 = p1;
                c2 = p2;
            }
            
            // 变异
            mutate(c1);
            mutate(c2);
            
            offspring.push_back(c1);
            if (offspring.size() < params.pop_size) {
                offspring.push_back(c2);
            }
        }
        
        return offspring;
    }
    
    // SBX交叉
    void crossover(const Individual& p1, const Individual& p2, 
                  Individual& c1, Individual& c2) {
        uniform_real_distribution<double> dist(0.0, 1.0);
        
        for (int i = 0; i < params.n_var; i++) {
            double u = dist(rng);
            double beta;
            if (u <= 0.5) {
                beta = pow(2.0 * u, 1.0 / (params.eta_c + 1.0));
            } else {
                beta = pow(1.0 / (2.0 * (1.0 - u)), 1.0 / (params.eta_c + 1.0));
            }
            
            double x1 = p1.variables[i];
            double x2 = p2.variables[i];
            
            c1.variables[i] = 0.5 * ((1 + beta) * x1 + (1 - beta) * x2);
            c2.variables[i] = 0.5 * ((1 - beta) * x1 + (1 + beta) * x2);
            
            // 边界处理
            c1.variables[i] = max(params.var_min[i], 
                                 min(params.var_max[i], c1.variables[i]));
            c2.variables[i] = max(params.var_min[i], 
                                 min(params.var_max[i], c2.variables[i]));
        }
    }
    
    // 多项式变异
    void mutate(Individual& ind) {
        uniform_real_distribution<double> dist(0.0, 1.0);
        
        for (int i = 0; i < params.n_var; i++) {
            if (dist(rng) < params.mutation_prob) {
                double u = dist(rng);
                double delta;
                if (u < 0.5) {
                    delta = pow(2.0 * u, 1.0 / (params.eta_m + 1.0)) - 1.0;
                } else {
                    delta = 1.0 - pow(2.0 * (1.0 - u), 1.0 / (params.eta_m + 1.0));
                }
                
                ind.variables[i] += delta * (params.var_max[i] - params.var_min[i]);
                
                // 边界处理
                ind.variables[i] = max(params.var_min[i], 
                                      min(params.var_max[i], ind.variables[i]));
            }
        }
    }
    
    // 环境选择
    vector<Individual> environmentalSelection(vector<Individual> combined, int pop_size) {
        vector<vector<int>> fronts = nonDominatedSort(combined);
        vector<Individual> new_pop;
        
        for (const auto& front : fronts) {
            if (new_pop.size() + front.size() <= pop_size) {
                // 整个前沿加入
                for (int idx : front) {
                    new_pop.push_back(combined[idx]);
                }
            } else {
                // 部分加入，按拥挤度排序
                vector<int> sorted = front;
                sort(sorted.begin(), sorted.end(), 
                    int a, int b {
                        return combined[a].crowding_distance > 
                               combined[b].crowding_distance;
                    });
                
                int remaining = pop_size - new_pop.size();
                for (int i = 0; i < remaining; i++) {
                    new_pop.push_back(combined[sorted[i]]);
                }
                break;
            }
        }
        
        return new_pop;
    }
    
    // 提取Pareto前沿
    vector<Individual> extractParetoFront(const vector<int>& indices) {
        vector<Individual> front;
        for (int idx : indices) {
            front.push_back(population[idx]);
        }
        return front;
    }
};

2.3 主函数与结果输出

// 输出Pareto前沿到文件
void writeParetoFront(const vector<Individual>& front, const string& filename) {
    ofstream outfile(filename);
    if (!outfile.is_open()) {
        cerr << "Error opening file: " << filename << endl;
        return;
    }
    
    // 写入表头
    outfile << "Variable1";
    for (int i = 1; i < front[0].variables.size(); i++) {
        outfile << ",Variable" << (i+1);
    }
    for (int i = 0; i < front[0].objectives.size(); i++) {
        outfile << ",Objective" << (i+1);
    }
    outfile << "\n";
    
    // 写入数据
    for (const auto& ind : front) {
        for (int i = 0; i < ind.variables.size(); i++) {
            if (i > 0) outfile << ",";
            outfile << ind.variables[i];
        }
        for (int i = 0; i < ind.objectives.size(); i++) {
            outfile << "," << ind.objectives[i];
        }
        outfile << "\n";
    }
    
    outfile.close();
}

int main() {
    // 设置参数
    Parameters params;
    params.pop_size = 100;
    params.max_generations = 500;
    params.n_var = 10;  // ZDT1问题有10个变量
    params.n_obj = 2;
    
    // 创建并运行算法
    NSGA2 nsga2(params);
    auto start = chrono::high_resolution_clock::now();
    nsga2.run();
    auto end = chrono::high_resolution_clock::now();
    
    // 计算运行时间
    auto duration = chrono::duration_cast<chrono::milliseconds>(end - start);
    cout << "NSGA-II completed in " << duration.count() << " ms" << endl;
    
    // 获取并输出结果
    vector<Individual> pareto = nsga2.getParetoFront();
    cout << "Pareto front size: " << pareto.size() << endl;
    
    // 输出到文件
    writeParetoFront(pareto, "pareto_front.csv");
    
    return 0;
}

三、关键算法优化

3.1 并行计算加速

// 并行评估种群
void evaluatePopulation() {
    #pragma omp parallel for schedule(dynamic)
    for (int i = 0; i < population.size(); i++) {
        // 目标函数计算
        double f1 = population[i].variables[0];
        double sum = 0.0;
        for (int j = 1; j < params.n_var; j++) {
            sum += population[i].variables[j];
        }
        double g = 1.0 + 9.0 * sum / (params.n_var - 1);
        double f2 = g * (1.0 - sqrt(f1 / g));
        
        population[i].objectives[0] = f1;
        population[i].objectives[1] = f2;
    }
}

3.2 内存优化

// 使用指针减少数据复制
void environmentalSelection(vector<Individual>& combined, int pop_size) {
    // ... 使用指针操作而非复制整个个体
}

3.3 性能监控

// 在run()方法中添加
auto start_gen = chrono::high_resolution_clock::now();
// ... 主循环代码
auto end_gen = chrono::high_resolution_clock::now();
auto gen_duration = chrono::duration_cast<chrono::milliseconds>(end_gen - start_gen);
cout << "Generation " << gen << " took " << gen_duration.count() << " ms" << endl;

四、应用案例：ZDT1测试函数

4.1 ZDT1问题描述

4.2 结果可视化（Python脚本）

import matplotlib.pyplot as plt
import pandas as pd

# 读取结果
df = pd.read_csv('pareto_front.csv')

# 绘制Pareto前沿
plt.figure(figsize=(10, 6))
plt.scatter(df['Objective1'], df['Objective2'], c='blue', alpha=0.6)
plt.xlabel('f1')
plt.ylabel('f2')
plt.title('NSGA-II Results on ZDT1 Problem')
plt.grid(True)

# 绘制理论Pareto前沿
x = [0.01*i for i in range(101)]
y = [1 - (xi**0.5) for xi in x]
plt.plot(x, y, 'r--', label='Theoretical Pareto Front')
plt.legend()
plt.savefig('pareto_front.png')
plt.show()

五、性能评估与参数调优

5.1 性能指标

// 计算超体积指标
double calculateHypervolume(const vector<Individual>& front, 
                           const vector<double>& ref_point) {
    // 按第一个目标排序
    vector<Individual> sorted = front;
    sort(sorted.begin(), sorted.end(), 
        const Individual& a, const Individual& b {
            return a.objectives[0] < b.objectives[0];
        });
    
    double hv = 0.0;
    double prev_y = ref_point[1];
    
    for (const auto& ind : sorted) {
        double width = ref_point[0] - ind.objectives[0];
        double height = prev_y - ind.objectives[1];
        hv += width * height;
        prev_y = ind.objectives[1];
    }
    
    return hv;
}

5.2 参数调优建议

参数 推荐值 影响 调优策略
种群大小 50-200 多样性vs计算量 问题复杂度↑ → 种群大小↑
交叉概率 0.8-0.95 探索能力 保持较高值
变异概率 1/n_var 开发能力 随维度增加而减小
分布指数η_c 10-20 交叉分布 值越大，子代越接近父代
分布指数η_m 20-100 变异幅度 值越大，变异越温和

参考代码 cpp 带精英策略的非支配排序遗传算法（NSGA-II） www.youwenfan.com/contentcst/160549.html

六、扩展功能

6.1 约束处理

// 在评估函数中添加约束违反度
void evaluateWithConstraints(Individual& ind) {
    // 计算目标函数
    // ...
    
    // 计算约束违反度
    double constraint_violation = 0.0;
    for (auto& c : constraints) {
        double val = c.evaluate(ind.variables);
        if (val > 0) {
            constraint_violation += val;
        }
    }
    
    // 惩罚目标函数
    for (int i = 0; i < ind.objectives.size(); i++) {
        ind.objectives[i] += 1e6 * constraint_violation;
    }
}

6.2 动态目标函数

// 运行时更换目标函数
void setObjectiveFunction(function<void(Individual&)> func) {
    objective_func = func;
}

// 在评估函数中使用
void evaluatePopulation() {
    for (auto& ind : population) {
        objective_func(ind);
    }
}

6.3 混合优化策略

// 在进化过程中加入局部搜索
void localSearch(Individual& ind) {
    // 使用Nelder-Mead或其他局部搜索方法
    // 改进当前解
}

// 在生成子代后应用
void generateOffspring(...) {
    // ... 交叉变异后
    if (shouldApplyLocalSearch()) {
        localSearch(offspring[i]);
    }
}

七、编译与运行

7.1 编译命令

# 使用g++编译
g++ -std=c++11 -O3 -fopenmp nsga2.cpp -o nsga2

# 使用CMake
mkdir build
cd build
cmake ..
make

7.2 运行示例

./nsga2

# 输出示例
Generation 0, First Front Size: 18
Generation 50, First Front Size: 32
Generation 100, First Front Size: 35
...
NSGA-II completed in 1250 ms
Pareto front size: 42

八、总结

本实现提供了NSGA-II算法的完整C++解决方案，具有以下特点：

  1. 高效实现：使用STL容器和现代C++特性
  2. 并行计算：OpenMP加速种群评估
  3. 模块化设计：易于扩展新的选择、交叉和变异算子
  4. 工业级应用：适合嵌入式系统和实时优化
  5. 全面功能：支持约束处理、动态目标函数和混合优化

应用场景：

通过本实现，开发者可以快速将NSGA-II算法应用于各种多目标优化问题，并根据具体需求进行定制和扩展。

专注于matlab/simulink,电子电路,编程

关于我们湘ICP备19023172号-4