网站导航
> 技术文档 > PAES算法求解 ZDT1 双目标优化问题

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

网友: 技术文档 2025-09-03 11:37:28

前言

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

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

文章目录

  • 前言
  • matlab代码
  • 代码简析
      • 一、核心数学公式与算法逻辑
      • 二、代码模块与公式对应解析
        • 1. ZDT1目标函数（数学定义与实现）
        • 2. 高斯变异（探索新解的数学逻辑）
        • 3. 支配关系判断（多目标优化的核心逻辑）
        • 4. 网格拥挤度（维护解分布性的工具）
        • 5. 存档维护（保留优质解的逻辑）
        • 6. 主循环（算法流程的串联）
      • 三、算法流程总结
      • 四、关键可视化与指标

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\'); endend

运行结果
在这里插入图片描述
在这里插入图片描述

代码简析

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

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

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

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

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

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

ZDT1 双目标优化问题简介

ZDT1的两个目标函数：
f 1 ( x ) = x 1 （第一个目标，直接取第一个决策变量） g ( x ) = 1 + 9 ⋅ ∑ i = 2 30 x i 29 （辅助函数，平衡多变量影响） f 2 ( x ) = g ( x ) ⋅ ( 1 − f 1( x ) g ( x ) ) （第二个目标，与 f 1 和 g 相关） \\begin{align*} f_1(x) &= x_1 \\quad \\text{（第一个目标，直接取第一个决策变量）} \\\\ g(x) &= 1 + 9 \\cdot \\frac{\\sum_{i=2}^{30} x_i}{29} \\quad \\text{（辅助函数，平衡多变量影响）} \\\\ f_2(x) &= g(x) \\cdot \\left(1 - \\sqrt{\\frac{f_1(x)}{g(x)}}\\right) \\quad \\text{（第二个目标，与} f_1 \\text{和} g \\text{相关）} \\end{align*} f1​(x)g(x)f2​(x)​=x1​（第一个目标，直接取第一个决策变量）=1+9⋅29∑i=230​xi​​（辅助函数，平衡多变量影响）=g(x)⋅(1−g(x)f1​(x)​ ​)（第二个目标，与f1​和g相关）​

代码实现（嵌套函数 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 ′ = x i + mutation_rate ⋅ N ( 0 , 1 ) x_i\' = x_i + \\text{mutation\\_rate} \\cdot \\mathcal{N}(0, 1) xi′​=xi​+mutation_rate⋅N(0,1)
其中 N(0,1) \\mathcal{N}(0, 1) 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 A A在所有目标上不差于解 B B B，且至少一个目标严格优于 B B B，则称 A A A支配 B B B，公式：
dominates ( A , B ) = ( ∀ i , A i ≤ B i ) ∧ ( ∃ j , A j < B j ) \\text{dominates}(A, B) = \\left( \\forall i, A_i \\leq B_i \\right) \\land \\left( \\exists j, A_j < B_j \\right) dominates(A,B)=(∀i,Ai​≤Bi​)∧(∃j,Aj​<Bj​)

代码实现（嵌套函数 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 ⁡ _ o b j = min ⁡ ( archive_obj , [ ] , 1 ) \\min\\_obj = \\min(\\text{archive\\_obj}, [], 1) min_obj=min(archive_obj,[],1)， max ⁡ _ o b j = max ⁡ ( archive_obj , [ ] , 1 ) \\max\\_obj = \\max(\\text{archive\\_obj}, [], 1) max_obj=max(archive_obj,[],1)
  2. 归一化目标值： normalized = objectives − min ⁡ _ o b j max ⁡ _ o b j − min ⁡ _ o b j \\text{normalized} = \\frac{\\text{objectives} - \\min\\_obj}{\\max\\_obj - \\min\\_obj} normalized=max_obj−min_objobjectives−min_obj​（避免除零，若范围为0则设为1）
  3. 计算网格坐标： grid_coords = ⌊ normalized ⋅ grid_divisions ⌋ \\text{grid\\_coords} = \\lfloor \\text{normalized} \\cdot \\text{grid\\_divisions} \\rfloor grid_coords=⌊normalized⋅grid_divisions⌋
  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 endend
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); endend

三、算法流程总结

PAES通过以下步骤求解ZDT1：

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

四、关键可视化与指标

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

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

网络标签：

相关问题