PAES算法求解 ZDT1 双目标优化问题

前言

提醒：
文章内容为方便作者自己后日复习与查阅而进行的书写与发布，其中引用内容都会使用链接表明出处（如有侵权问题，请及时联系）。
其中内容多为一次书写，缺少检查与订正，如有问题或其他拓展及意见建议，欢迎评论区讨论交流。

内容由AI辅助生成，仅经笔者审核整理，请甄别食用。

matlab代码

function PAES_ZDT1_Complete()
    % PAES算法求解ZDT1问题 - 完整实现
    % Pareto Archived Evolution Strategy for ZDT1 Problem
    
    clc; clear; close all;
    
    fprintf('=== PAES算法求解ZDT1问题 ===\n');
    fprintf('开始优化...\n\n');
    
    % ==================== 算法参数设置 ====================
    params.max_iterations = 3*10000;    % 最大迭代次数
    params.archive_size = 100;        % 档案最大容量
    params.grid_divisions = 10;       % 网格划分数
    params.mutation_rate = 0.1;       % 变异率
    params.dimension = 30;            % 决策变量维数
    params.num_objectives = 2;        % 目标函数数量
    
    % ==================== 初始化 ====================
    % 随机生成初始解
    current_solution = rand(1, params.dimension);  % ZDT1变量范围[0,1]
    archive = current_solution;
    
    % 计算初始解的目标函数值
    current_obj = ZDT1_objective(current_solution);
    archive_obj = current_obj;
    
    % ==================== 主优化循环 ====================
    tic;
    for iter = 1:params.max_iterations
        % 1. 变异产生候选解
        candidate = mutate(current_solution, params.mutation_rate);
        candidate_obj = ZDT1_objective(candidate);
        
        % 2. 评估候选解
        [accept_candidate, new_current, archive, archive_obj] = ...
            evaluate_candidate(current_solution, current_obj, ...
                             candidate, candidate_obj, ...
                             archive, archive_obj, params);
        
        % 3. 更新当前解
        if accept_candidate
            current_solution = new_current;
            current_obj = ZDT1_objective(current_solution);
        end
        
        % 4. 维护档案大小
        if size(archive, 1) > params.archive_size
            [archive, archive_obj] = maintain_archive(archive, archive_obj, params);
        end
        
        % 5. 显示进度
        if mod(iter, 2000) == 0
            fprintf('迭代 %d/%d: 档案大小 = %d\n', iter, params.max_iterations, size(archive, 1));
        end
    end
    
    elapsed_time = toc;
    
    % ==================== 结果输出和可视化 ====================
    fprintf('\n=== 优化完成 ===\n');
    fprintf('总耗时: %.2f秒\n', elapsed_time);
    fprintf('最终档案大小: %d\n', size(archive, 1));
    fprintf('f1范围: [%.4f, %.4f]\n', min(archive_obj(:, 1)), max(archive_obj(:, 1)));
    fprintf('f2范围: [%.4f, %.4f]\n', min(archive_obj(:, 2)), max(archive_obj(:, 2)));
    
    % 绘制结果
    plot_results(archive_obj);
    
    % ==================== 嵌套函数定义 ====================
    
    function objectives = ZDT1_objective(x)
        % ZDT1测试函数
        % f1(x) = x1
        % f2(x) = g(x) * h(f1, g)
        % g(x) = 1 + 9/(n-1) * sum(xi, i=2 to n)
        % h(f1, g) = 1 - sqrt(f1/g)
        
        n = length(x);
        
        % 第一个目标函数
        f1 = x(1);
        
        % 辅助函数g
        if n > 1
            g = 1 + 9 * sum(x(2:end)) / (n - 1);
        else
            g = 1;
        end
        
        % 第二个目标函数
        h = 1 - sqrt(f1 / g);
        f2 = g * h;
        
        objectives = [f1, f2];
    end
    
    function mutated_solution = mutate(solution, mutation_rate)
        % 高斯变异操作
        mutated_solution = solution + mutation_rate * randn(size(solution));
        
        % 边界处理 - 确保变量在[0,1]范围内
        mutated_solution = max(0, min(1, mutated_solution));
    end
    
    function result = dominates(obj1, obj2)
        % 判断obj1是否支配obj2 (最小化问题)
        % 支配条件：obj1在所有目标上不劣于obj2，且至少在一个目标上严格优于obj2
        
        all_better_or_equal = all(obj1 <= obj2);
        at_least_one_better = any(obj1 < obj2);
        
        result = all_better_or_equal && at_least_one_better;
    end
    
    function grid_coords = calculate_grid_coordinates(objectives, archive_obj, grid_divisions)
        % 计算解在自适应网格中的坐标
        
        if size(archive_obj, 1) == 1
            grid_coords = zeros(1, length(objectives));
            return;
        end
        
        % 计算目标空间的边界
        min_obj = min(archive_obj, [], 1);
        max_obj = max(archive_obj, [], 1);
        
        % 避免除零错误
        range_obj = max_obj - min_obj;
        range_obj(range_obj == 0) = 1;
        
        % 计算归一化坐标
        normalized = (objectives - min_obj) ./ range_obj;
        
        % 计算网格坐标
        grid_coords = floor(normalized * grid_divisions);
        
        % 确保坐标在有效范围内
        grid_coords = max(0, min(grid_divisions - 1, grid_coords));
    end
    
    function crowding = calculate_crowding(grid_coords, archive_obj, grid_divisions)
        % 计算指定网格的拥挤度（该网格中解的数量）
        
        if size(archive_obj, 1) <= 1
            crowding = 1;
            return;
        end
        
        % 计算档案中所有解的网格坐标
        crowding = 0;
        for i = 1:size(archive_obj, 1)
            current_grid = calculate_grid_coordinates(archive_obj(i, :), archive_obj, grid_divisions);
            if isequal(current_grid, grid_coords)
                crowding = crowding + 1;
            end
        end
    end
    
    function [accept, new_current, new_archive, new_archive_obj] = ...
        evaluate_candidate(current, current_obj, candidate, candidate_obj, archive, archive_obj, params)
        % 根据PAES规则评估候选解是否被接受
        
        accept = false;
        new_current = current;
        new_archive = archive;
        new_archive_obj = archive_obj;
        
        % 情况1: 候选解支配当前解
        if dominates(candidate_obj, current_obj)
            accept = true;
            new_current = candidate;
            % 将候选解添加到档案
            new_archive = [new_archive; candidate];
            new_archive_obj = [new_archive_obj; candidate_obj];
            return;
        end
        
        % 情况2: 当前解支配候选解
        if dominates(current_obj, candidate_obj)
            return;  % 拒绝候选解
        end
        
        % 情况3: 两解互不支配，需要进一步判断
        
        % 3.1 检查候选解是否被档案中的解支配
        for i = 1:size(archive_obj, 1)
            if dominates(archive_obj(i, :), candidate_obj)
                return;  % 被档案中的解支配，拒绝
            end
        end
        
        % 3.2 检查候选解是否支配档案中的解
        dominated_indices = [];
        for i = 1:size(archive_obj, 1)
            if dominates(candidate_obj, archive_obj(i, :))
                dominated_indices = [dominated_indices, i];
            end
        end
        
        % 如果候选解支配档案中的某些解
        if ~isempty(dominated_indices)
            accept = true;
            new_current = candidate;
            % 从档案中移除被支配的解
            keep_indices = setdiff(1:size(archive, 1), dominated_indices);
            new_archive = new_archive(keep_indices, :);
            new_archive_obj = new_archive_obj(keep_indices, :);
            % 添加候选解到档案
            new_archive = [new_archive; candidate];
            new_archive_obj = [new_archive_obj; candidate_obj];
            return;
        end
        
        % 3.3 候选解与档案中所有解互不支配
        if size(archive, 1) < params.archive_size
            % 档案未满，直接接受
            accept = true;
            new_current = candidate;
            new_archive = [new_archive; candidate];
            new_archive_obj = [new_archive_obj; candidate_obj];
        else
            % 档案已满，根据拥挤度判断
            candidate_grid = calculate_grid_coordinates(candidate_obj, [archive_obj; candidate_obj], params.grid_divisions);
            current_grid = calculate_grid_coordinates(current_obj, [archive_obj; candidate_obj], params.grid_divisions);
            
            candidate_crowding = calculate_crowding(candidate_grid, [archive_obj; candidate_obj], params.grid_divisions);
            current_crowding = calculate_crowding(current_grid, [archive_obj; candidate_obj], params.grid_divisions);
            
            % 如果候选解所在网格的拥挤度更小，则接受
            if candidate_crowding < current_crowding
                accept = true;
                new_current = candidate;
                new_archive = [new_archive; candidate];
                new_archive_obj = [new_archive_obj; candidate_obj];
            end
        end
    end
    
    function [new_archive, new_archive_obj] = maintain_archive(archive, archive_obj, params)
        % 当档案超过容量限制时，移除拥挤度最高的解
        
        while size(archive, 1) > params.archive_size
            % 计算所有解的拥挤度
            crowding_values = zeros(size(archive, 1), 1);
            for i = 1:size(archive, 1)
                grid_coords = calculate_grid_coordinates(archive_obj(i, :), archive_obj, params.grid_divisions);
                crowding_values(i) = calculate_crowding(grid_coords, archive_obj, params.grid_divisions);
            end
            
            % 找到拥挤度最高的解的索引
            [~, max_crowding_indices] = max(crowding_values);
            
            % 如果有多个解具有相同的最高拥挤度，随机选择一个
            if length(max_crowding_indices) > 1
                remove_idx = max_crowding_indices(randi(length(max_crowding_indices)));
            else
                remove_idx = max_crowding_indices(1);
            end
            
            % 移除选中的解
            archive(remove_idx, :) = [];
            archive_obj(remove_idx, :) = [];
        end
        
        new_archive = archive;
        new_archive_obj = archive_obj;
    end
    
    function plot_results(archive_obj)
        % 绘制优化结果对比图
        
        figure('Position', [100, 100, 800, 600]);
        
        % 绘制PAES找到的解
        scatter(archive_obj(:, 1), archive_obj(:, 2), 60, 'ro', 'filled', 'MarkerEdgeColor', 'k');
        hold on;
        
        % 绘制ZDT1的真实Pareto前沿
        f1_true = linspace(0, 1, 1000);
        f2_true = 1 - sqrt(f1_true);
        plot(f1_true, f2_true, 'b-', 'LineWidth', 2.5);
        
        % 图形美化
        xlabel('f_1', 'FontSize', 14, 'FontWeight', 'bold');
        ylabel('f_2', 'FontSize', 14, 'FontWeight', 'bold');
        title('PAES算法求解ZDT1问题结果对比', 'FontSize', 16, 'FontWeight', 'bold');
        legend({'PAES找到的解', '真实Pareto前沿'}, 'FontSize', 12, 'Location', 'northeast');
        grid on;
        grid minor;
        
        % 设置坐标轴
        xlim([0, 1]);
        ylim([0, 1]);
        
        % 添加文本信息
        text(0.05, 0.95, sprintf('解的数量: %d', size(archive_obj, 1)), ...
             'FontSize', 12, 'BackgroundColor', 'white', 'EdgeColor', 'black');
        
        % 计算并显示一些性能指标
        fprintf('\n=== 性能分析 ===\n');
        
        % 计算分布均匀性（相邻解之间的平均距离）
        if size(archive_obj, 1) > 1
            [~, sort_idx] = sort(archive_obj(:, 1));
            sorted_obj = archive_obj(sort_idx, :);
            distances = sqrt(sum(diff(sorted_obj).^2, 2));
            avg_distance = mean(distances);
            std_distance = std(distances);
            fprintf('解的平均间距: %.4f\n', avg_distance);
            fprintf('间距标准差: %.4f\n', std_distance);
        end
        
        % 计算收敛性（与真实前沿的平均距离）
        min_distances = zeros(size(archive_obj, 1), 1);
        for i = 1:size(archive_obj, 1)
            f1_current = archive_obj(i, 1);
            f2_true_at_f1 = 1 - sqrt(f1_current);
            min_distances(i) = abs(archive_obj(i, 2) - f2_true_at_f1);
        end
        avg_convergence = mean(min_distances);
        fprintf('与真实前沿的平均距离: %.6f\n', avg_convergence);
        
        fprintf('\n算法运行完成！\n');
    end

end

运行结果

代码简析

相关引用：PAES (Pareto Archived Evolution Strategy)优化方法简介

以下结合数学公式和代码逻辑，对PAES求解ZDT1问题的完整实现进行拆解讲解，帮助理解算法如何通过“变异-评估-存档维护”流程逼近帕累托前沿。

一、核心数学公式与算法逻辑

PAES（Pareto Archived Evolution Strategy）的核心是通过高斯变异探索解空间，用存档（archive）维护非支配解，并基于“支配关系”和“网格拥挤度”筛选解，最终逼近ZDT1的真实帕累托前沿$f_2=1-\sqrt{f_1}$。

二、代码模块与公式对应解析

1. ZDT1目标函数（数学定义与实现）

ZDT1 双目标优化问题简介

ZDT1的两个目标函数：
$$\begin{array}{rl}{f}_1(x)& =x_1\text{（第一个目标，直接取第一个决策变量）}\\ g(x)& =1+9\cdot \frac{{\sum }_{i=2}^{30}{x}_i}{29}\text{（辅助函数，平衡多变量影响）}\\ f_2(x)& =g(x)\cdot \left(1-\sqrt{\frac{f_1(x)}{g(x)}}\right)\text{（第二个目标，与}{f}_1\text{和}g\text{相关）}\end{array}$$

代码实现（嵌套函数 ZDT1_objective）：

function objectives = ZDT1_objective(x)
    f1 = x(1);  % 直接取第一个变量
    g = 1 + 9 * sum(x(2:end))/(length(x)-1);  % 计算g(x)
    f2 = g * (1 - sqrt(f1/g));  % 计算第二个目标
    objectives = [f1, f2];  % 输出目标向量
end

2. 高斯变异（探索新解的数学逻辑）

PAES通过高斯变异生成候选解，公式：
$$x_i^{\prime }=x_i+\text{mutation\_rate}\cdot \mathcal{N}(0,1)$$
其中$\mathcal{N}(0,1)$是标准正态分布随机数，mutation_rate 控制变异幅度（代码中设为 0.1）。

代码实现（嵌套函数 mutate）：

function mutated_solution = mutate(solution, mutation_rate)
    % 高斯变异：给原解加正态分布噪声
    mutated_solution = solution + mutation_rate * randn(size(solution));  
    mutated_solution = max(0, min(1, mutated_solution));  % 边界截断（保证在[0,1]）
end

3. 支配关系判断（多目标优化的核心逻辑）

Pareto 最优解（Pareto Optimal Solution）简介

若解$A$在所有目标上不差于解$B$，且至少一个目标严格优于$B$，则称$A$支配$B$，公式：
$$\text{dominates}(A,B)=\left(\forall i,A_i\le B_i\right)\wedge \left(\exists j,A_j<B_j\right)$$

代码实现（嵌套函数 dominates）：

function result = dominates(obj1, obj2)
    all_better_or_equal = all(obj1 <= obj2);  % 所有目标不差
    at_least_one_better = any(obj1 < obj2);    % 至少一个目标更优
    result = all_better_or_equal && at_least_one_better;  % 同时满足则支配
end

4. 网格拥挤度（维护解分布性的工具）

为避免解过度集中，PAES用网格划分量化拥挤度：

1. 计算目标空间边界：$\min \_obj=\min (\text{archive\_obj},[],1)$，$\max \_obj=\max (\text{archive\_obj},[],1)$
2. 归一化目标值：$\text{normalized}=\frac{\text{objectives}-\min \_obj}{\max \_obj-\min \_obj}$（避免除零，若范围为0则设为1）
3. 计算网格坐标：$\text{grid\_coords}=\lfloor \text{normalized}\cdot \text{grid\_divisions}\rfloor$
4. 拥挤度定义：同一网格内的解数量，数量越多则拥挤度越高。

代码实现（嵌套函数 calculate_crowding）：

function crowding = calculate_crowding(grid_coords, archive_obj, grid_divisions)
    crowding = 0;
    for i = 1:size(archive_obj, 1)
        current_grid = calculate_grid_coordinates(archive_obj(i, :), archive_obj, grid_divisions);
        if isequal(current_grid, grid_coords)
            crowding = crowding + 1;  % 统计同网格解数量
        end
    end
end

5. 存档维护（保留优质解的逻辑）

外部存档（External Archive）机制

• 添加解：若候选解不被支配，且存档未满则直接加入；若存档已满，比较“候选解所在网格”与“当前解所在网格”的拥挤度，拥挤度高的网格移除解。
• 移除解：当存档超过容量（archive_size），持续移除“拥挤度最高网格”中的解，直到容量达标。

代码实现（嵌套函数 maintain_archive）：

function [new_archive, new_archive_obj] = maintain_archive(archive, archive_obj, params)
    while size(archive, 1) > params.archive_size
        % 计算所有解的拥挤度
        crowding_values = zeros(size(archive, 1), 1);
        for i = 1:size(archive, 1)
            grid_coords = calculate_grid_coordinates(archive_obj(i, :), archive_obj, params.grid_divisions);
            crowding_values(i) = calculate_crowding(grid_coords, archive_obj, params.grid_divisions);
        end
        
        % 移除拥挤度最高的解（随机选一个拥挤度最高的）
        [~, max_crowding_indices] = max(crowding_values);
        remove_idx = max_crowding_indices(randi(length(max_crowding_indices)));
        
        archive(remove_idx, :) = [];  % 移除解
        archive_obj(remove_idx, :) = [];  % 移除对应目标值
    end
    new_archive = archive;
    new_archive_obj = archive_obj;
end

6. 主循环（算法流程的串联）

PAES通过“变异→评估→更新→维护存档”的循环迭代优化：

for iter = 1:params.max_iterations
    % 1. 变异产生候选解
    candidate = mutate(current_solution, params.mutation_rate);
    candidate_obj = ZDT1_objective(candidate);
    
    % 2. 评估候选解（支配关系+存档规则）
    [accept_candidate, new_current, archive, archive_obj] = evaluate_candidate(...);
    
    % 3. 更新当前解
    if accept_candidate
        current_solution = new_current;
        current_obj = ZDT1_objective(current_solution);
    end
    
    % 4. 维护存档大小
    if size(archive, 1) > params.archive_size
        [archive, archive_obj] = maintain_archive(archive, archive_obj, params);
    end
end

三、算法流程总结

PAES通过以下步骤求解ZDT1：

1. 初始化：随机生成初始解，计算目标值并初始化存档。
2. 变异：用高斯变异生成候选解，探索新的解空间。
3. 评估：基于支配关系判断是否接受候选解，若接受则更新当前解和存档。
4. 维护存档：通过网格拥挤度移除冗余解，保证存档解分布均匀。
5. 循环迭代：重复“变异-评估-维护”，直到达到最大迭代次数。

四、关键可视化与指标

1. 真实帕累托前沿：f2_true = 1 - sqrt(f1_true)，代码中用蓝色线绘制。
2. 算法找到的解：存档中的非支配解，用红色点绘制，分布越接近蓝色线且均匀，说明算法性能越好。
3. 性能指标：
   • 平均间距：排序后相邻解的平均欧氏距离，反映分布均匀性。
   • 与真实前沿的平均距离：解到$f_2=1-\sqrt{f_1}$的垂直距离均值，反映收敛性。

该代码完整实现了PAES的核心逻辑，通过数学公式（如支配关系、高斯变异）与工程实现（如网格拥挤度、存档维护）的结合，有效求解ZDT1问题并逼近真实帕累托前沿。