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

前言

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

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

---

matlab代码

function PAES_ZDT1_Complete()% PAES算法求解ZDT1问题 - 完整实现% Pareto Archived Evolution Strategy for ZDT1 Problemclc; 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_candidatecurrent_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) == 0fprintf('迭代 %d/%d: 档案大小 = %d\n', iter, params.max_iterations, size(archive, 1));endendelapsed_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);% 辅助函数gif n > 1g = 1 + 9 * sum(x(2:end)) / (n - 1);elseg = 1;end% 第二个目标函数h = 1 - sqrt(f1 / g);f2 = g * h;objectives = [f1, f2];endfunction mutated_solution = mutate(solution, mutation_rate)% 高斯变异操作mutated_solution = solution + mutation_rate * randn(size(solution));% 边界处理 - 确保变量在[0,1]范围内mutated_solution = max(0, min(1, mutated_solution));endfunction result = dominates(obj1, obj2)% 判断obj1是否支配obj2 (最小化问题)% 支配条件：obj1在所有目标上不劣于obj2，且至少在一个目标上严格优于obj2all_better_or_equal = all(obj1 <= obj2);at_least_one_better = any(obj1 < obj2);result = all_better_or_equal && at_least_one_better;endfunction grid_coords = calculate_grid_coordinates(objectives, archive_obj, grid_divisions)% 计算解在自适应网格中的坐标if size(archive_obj, 1) == 1grid_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));endfunction crowding = calculate_crowding(grid_coords, archive_obj, grid_divisions)% 计算指定网格的拥挤度（该网格中解的数量）if size(archive_obj, 1) <= 1crowding = 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;endendendfunction [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;  % 被档案中的解支配，拒绝endend% 3.2 检查候选解是否支配档案中的解dominated_indices = [];for i = 1:size(archive_obj, 1)if dominates(candidate_obj, archive_obj(i, :))dominated_indices = [dominated_indices, i];endend% 如果候选解支配档案中的某些解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_crowdingaccept = true;new_current = candidate;new_archive = [new_archive; candidate];new_archive_obj = [new_archive_obj; candidate_obj];endendendfunction [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) > 1remove_idx = max_crowding_indices(randi(length(max_crowding_indices)));elseremove_idx = max_crowding_indices(1);end% 移除选中的解archive(remove_idx, :) = [];archive_obj(remove_idx, :) = [];endnew_archive = archive;new_archive_obj = archive_obj;endfunction 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);endavg_convergence = mean(min_distances);fprintf('与真实前沿的平均距离: %.6f\n', avg_convergence);fprintf('\n算法运行完成！\n');endend

运行结果

代码简析

相关引用：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;  % 统计同网格解数量endend
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, :) = [];  % 移除对应目标值endnew_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_candidatecurrent_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问题并逼近真实帕累托前沿。