前言

提醒：
文章內容為方便作者自己后日復習與查閱而進行的書寫與發布，其中引用內容都會使用鏈接表明出處（如有侵權問題，請及時聯系）。
其中內容多為一次書寫，缺少檢查與訂正，如有問題或其他拓展及意見建議，歡迎評論區討論交流。

內容由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?\frac{{\sum }_{i=2}^{30}{x}_i}{29}\text{（輔助函數，平衡多變量影響）}\\ f_2(x)& =g(x)?\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}?\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(?i,A_i\le B_i\right)\wedge \left(?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}=?\text{normalized}?\text{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;  % 統計同網格解數量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問題并逼近真實帕累托前沿。