%% TSP问题：多种算法对比分析% 功能：对比暴力求解、最邻近算法、蚁群算法和intlinprog求解器clear; clc; close all;%% 参数设置n_cities = 6; % 城市数量（建议4-6，确保所有算法都能运行）%% 生成城市坐标rng(42); % 设置随机种子以保证结果可重现if n_cities == 6cities_x = [0, 1, 2, 1, 0, 3];cities_y = [0, 0, 0, 1, 1, 2];elseif n_cities == 5cities_x = [0, 2, 4, 2, 0,];cities_y = [0, 0, 2, 2, 2];else % 随机分布 cities_x = rand(n_cities, 1) * 10; cities_y = rand(n_cities, 1) * 10;end%% 计算距离矩阵distance_matrix = zeros(n_cities, n_cities);for i = 1:n_cities for j = 1:n_cities if i ~= j distance_matrix(i,j) = sqrt((cities_x(i)-cities_x(j))^2 + ... (cities_y(i)-cities_y(j))^2); else distance_matrix(i,j) = inf; % 避免自循环 end endendfprintf('=== TSP问题设置 ===\n');fprintf('城市数量: %d\n', n_cities);fprintf('城市坐标:\n');for i = 1:n_cities fprintf(' 城市%d: (%.2f, %.2f)\n', i, cities_x(i), cities_y(i));end%% 方法1：暴力求解（精确算法）fprintf('\n=== 暴力求解 ===\n');if n_cities <= 8 % 暴力求解限制在8个城市以内 [best_route_brute, best_distance_brute, time_brute] = solve_brute_force(distance_matrix);else fprintf('城市数过多，跳过暴力求解\n'); best_route_brute = []; best_distance_brute = NaN; time_brute = NaN;end%% 方法2：最邻近算法（贪心算法）fprintf('\n=== 最邻近算法 ===\n');[best_route_nn, best_distance_nn, time_nn] = solve_nearest_neighbor(distance_matrix);%% 方法3：蚁群算法（启发式算法）fprintf('\n=== 蚁群算法 ===\n');[best_route_aco, best_distance_aco, time_aco] = solve_aco(distance_matrix);%% 方法4：intlinprog求解器fprintf('\n=== intlinprog求解器 ===\n');if license('test', 'Optimization_Toolbox') && n_cities <= 6 [best_route_exact, best_distance_exact, time_exact, exact_success] = solve_tsp_practical(distance_matrix); if ~exact_success fprintf('intlinprog求解失败\n'); best_route_exact = []; best_distance_exact = NaN; time_exact = NaN; endelse if ~license('test', 'Optimization_Toolbox') fprintf('未安装MATLAB优化工具箱，跳过intlinprog求解\n'); else fprintf('城市数过多，跳过intlinprog求解\n'); end best_route_exact = []; best_distance_exact = NaN; time_exact = NaN;end%% 结果可视化figure('Position', [100, 100, 1400, 600]);% 子图1：城市分布和路径对比subplot(1, 3, 1);plot(cities_x, cities_y, 'ro', 'MarkerSize', 12, 'LineWidth', 2);hold on;% 标注城市编号for i = 1:n_cities text(cities_x(i)+0.1, cities_y(i)+0.1, num2str(i), 'FontSize', 12, 'FontWeight', 'bold');end% 绘制各种算法结果legend_items = {'城市', '城市编号'};% 最邻近算法（蓝色）if ~isempty(best_route_nn) route_x_nn = cities_x([best_route_nn, best_route_nn(1)]); route_y_nn = cities_y([best_route_nn, best_route_nn(1)]); plot(route_x_nn, route_y_nn, 'b-', 'LineWidth', 1.5); legend_items{end+1} = sprintf('最邻近(%.2f)', best_distance_nn);end% 蚁群算法（红色）if ~isempty(best_route_aco) route_x_aco = cities_x([best_route_aco, best_route_aco(1)]); route_y_aco = cities_y([best_route_aco, best_route_aco(1)]); plot(route_x_aco, route_y_aco, 'r-.', 'LineWidth', 2); legend_items{end+1} = sprintf('蚁群算法(%.2f)', best_distance_aco);end% intlinprog（绿色）if ~isempty(best_route_exact) route_x_exact = cities_x([best_route_exact, best_route_exact(1)]); route_y_exact = cities_y([best_route_exact, best_route_exact(1)]); plot(route_x_exact, route_y_exact, 'g--', 'LineWidth', 1.5); legend_items{end+1} = sprintf('intlinprog(%.2f)', best_distance_exact);end% 暴力求解（品红）if ~isempty(best_route_brute) route_x_brute = cities_x([best_route_brute, best_route_brute(1)]); route_y_brute = cities_y([best_route_brute, best_route_brute(1)]); plot(route_x_brute, route_y_brute, 'm:', 'LineWidth', 2); legend_items{end+1} = sprintf('暴力求解(%.2f)', best_distance_brute);endxlabel('X坐标');ylabel('Y坐标');title(sprintf('TSP问题多种算法求解结果对比 (%d个城市)', n_cities));legend(legend_items, 'Location', 'best');grid on;% 子图2：算法性能对比（距离）subplot(1, 3, 2);algorithms = {'最邻近', '蚁群算法', 'intlinprog', '暴力求解'};distances = [best_distance_nn, best_distance_aco, best_distance_exact, best_distance_brute];times = [time_nn, time_aco, time_exact, time_brute];% 只显示有效的算法valid_distances = distances(~isnan(distances));if ~isempty(valid_distances) best_known = min(valid_distances); valid_algorithms = {}; valid_distances_plot = []; valid_times_plot = []; if ~isnan(best_distance_nn) valid_algorithms{end+1} = '最邻近'; valid_distances_plot(end+1) = best_distance_nn; valid_times_plot(end+1) = time_nn; end if ~isnan(best_distance_aco) valid_algorithms{end+1} = '蚁群算法'; valid_distances_plot(end+1) = best_distance_aco; valid_times_plot(end+1) = time_aco; end if ~isnan(best_distance_exact) valid_algorithms{end+1} = 'intlinprog'; valid_distances_plot(end+1) = best_distance_exact; valid_times_plot(end+1) = time_exact; end if ~isnan(best_distance_brute) valid_algorithms{end+1} = '暴力求解'; valid_distances_plot(end+1) = best_distance_brute; valid_times_plot(end+1) = time_brute; end bar(1:length(valid_distances_plot), valid_distances_plot, 'FaceColor', [0.2, 0.6, 0.8]); set(gca, 'XTick', 1:length(valid_distances_plot), 'XTickLabel', valid_algorithms); ylabel('最短距离'); title('算法结果质量对比'); grid on; % 在柱状图上标注数值 for i = 1:length(valid_distances_plot) text(i, valid_distances_plot(i) + max(valid_distances_plot)*0.02, ... sprintf('%.2f', valid_distances_plot(i)), ... 'HorizontalAlignment', 'center', 'FontWeight', 'bold'); endend% 子图3：时间对比subplot(1, 3, 3);if exist('valid_times_plot', 'var') && ~isempty(valid_times_plot) bar(1:length(valid_times_plot), valid_times_plot, 'FaceColor', [0.8, 0.4, 0.2]); set(gca, 'XTick', 1:length(valid_times_plot), 'XTickLabel', valid_algorithms); ylabel('计算时间 (秒)'); title('算法计算时间对比'); grid on; % 在柱状图上标注数值 for i = 1:length(valid_times_plot) text(i, valid_times_plot(i) + max(valid_times_plot)*0.02, ... sprintf('%.3fs', valid_times_plot(i)), ... 'HorizontalAlignment', 'center', 'FontWeight', 'bold'); endend%% 性能对比表fprintf('\n=== 性能对比结果 ===\n');fprintf('城市数量: %d\n', n_cities);fprintf('--------------------------------------------------------\n');fprintf('%-15s %-15s %-15s %-15s\n', '算法', '最优距离', '计算时间(s)', '相对误差(%)');fprintf('--------------------------------------------------------\n');% 找到最佳解作为基准valid_distances = distances(~isnan(distances));if ~isempty(valid_distances) best_known = min(valid_distances); if ~isnan(best_distance_nn) error_nn = abs(best_distance_nn - best_known) / best_known * 100; fprintf('%-15s %-15.2f %-15.4f %-15.2f\n', '最邻近算法', best_distance_nn, time_nn, error_nn); end if ~isnan(best_distance_aco) error_aco = abs(best_distance_aco - best_known) / best_known * 100; fprintf('%-15s %-15.2f %-15.4f %-15.2f\n', '蚁群算法', best_distance_aco, time_aco, error_aco); end if ~isnan(best_distance_exact) error_exact = abs(best_distance_exact - best_known) / best_known * 100; fprintf('%-15s %-15.2f %-15.4f %-15.2f\n', 'intlinprog', best_distance_exact, time_exact, error_exact); end if ~isnan(best_distance_brute) error_brute = abs(best_distance_brute - best_known) / best_known * 100; fprintf('%-15s %-15.2f %-15.4f %-15.2f\n', '暴力求解', best_distance_brute, time_brute, error_brute); endendfprintf('--------------------------------------------------------\n');%% 暴力求解算法function [best_route, best_distance, time_taken] = solve_brute_force(distance_matrix) tic; n_cities = size(distance_matrix, 1); % 生成所有可能的排列（固定起点为城市1） if n_cities > 1 perms_list = perms(2:n_cities); else perms_list = []; end best_distance = inf; best_route = []; % 检查所有可能的路径 if isempty(perms_list) best_route = [1, 2]; best_distance = distance_matrix(1, 2) + distance_matrix(2, 1); else for i = 1:size(perms_list, 1) route = [1, perms_list(i, :)]; distance = calculate_total_distance(route, distance_matrix); if distance < best_distance best_distance = distance; best_route = route; end end end time_taken = toc; fprintf('最优距离: %.2f，耗时: %.4f秒\n', best_distance, time_taken); fprintf('最优路径: '); fprintf('%d ', best_route); fprintf('\n');end%% 最邻近算法function [best_route, best_distance, time_taken] = solve_nearest_neighbor(distance_matrix) tic; n_cities = size(distance_matrix, 1); best_route = nearest_neighbor_tsp(distance_matrix); best_distance = calculate_total_distance(best_route, distance_matrix); time_taken = toc; fprintf('最优距离: %.2f，耗时: %.6f秒\n', best_distance, time_taken); fprintf('最优路径: '); fprintf('%d ', best_route); fprintf('\n');endfunction route = nearest_neighbor_tsp(distance_matrix) n_cities = size(distance_matrix, 1); route = zeros(1, n_cities); % 从第一个城市开始 current_city = 1; route(1) = current_city; visited = false(1, n_cities); visited(current_city) = true; % 依次选择最近的未访问城市 for i = 2:n_cities min_distance = inf; next_city = 0; for j = 1:n_cities if ~visited(j) && j ~= current_city && isfinite(distance_matrix(current_city, j)) if distance_matrix(current_city, j) < min_distance min_distance = distance_matrix(current_city, j); next_city = j; end end end if next_city > 0 && next_city <= n_cities route(i) = next_city; visited(next_city) = true; current_city = next_city; else % 如果找不到有效城市，选择第一个未访问的城市 for j = 1:n_cities if ~visited(j) route(i) = j; visited(j) = true; current_city = j; break; end end end endend%% 蚁群算法function [best_route, best_distance, time_taken] = solve_aco(distance_matrix) tic; n_cities = size(distance_matrix, 1); max_iter = 30; n_ants = 15; alpha = 1; beta = 5; rho = 0.1; Q = 100; [best_route, best_distance, ~] = ant_colony_optimization(... distance_matrix, n_ants, max_iter, alpha, beta, rho, Q); time_taken = toc; fprintf('最优距离: %.2f，耗时: %.4f秒\n', best_distance, time_taken); fprintf('最优路径: '); fprintf('%d ', best_route); fprintf('\n');endfunction [best_route, best_distance, distance_history] = ant_colony_optimization(... distance_matrix, n_ants, max_iter, alpha, beta, rho, Q) n_cities = size(distance_matrix, 1); % 初始化信息素矩阵 pheromone = ones(n_cities, n_cities); % 初始化最佳解 best_distance = inf; best_route = []; distance_history = zeros(max_iter, 1); % 主循环 for iter = 1:max_iter % 每只蚂蚁构建解 routes = cell(n_ants, 1); distances = zeros(n_ants, 1); for ant = 1:n_ants routes{ant} = construct_solution(distance_matrix, pheromone, ... alpha, beta, n_cities); distances(ant) = calculate_total_distance(routes{ant}, distance_matrix); end % 更新全局最佳解 [min_dist, min_idx] = min(distances); if iter>1 if min_dist==best_distance rho=max(0.7*(1-(iter/max_iter)),0.1); else rho=0.1; end end if min_dist < best_distance && min_dist > 0 best_distance = min_dist; best_route = routes{min_idx}; end % 更新信息素 update_pheromone(pheromone, routes, distances, rho, Q, n_cities); distance_history(iter) = best_distance; end % 确保返回有效解 if isempty(best_route) best_route = 1:n_cities; best_distance = calculate_total_distance(best_route, distance_matrix); endendfunction route = construct_solution(distance_matrix, pheromone, alpha, beta, n_cities) % 随机选择起始城市 current_city = randi(n_cities); route = current_city; unvisited = setdiff(1:n_cities, current_city); % 逐个城市构建路径 while ~isempty(unvisited) % 计算转移概率 probabilities = zeros(length(unvisited), 1); for i = 1:length(unvisited) next_city = unvisited(i); if current_city >= 1 && current_city <= n_cities && ... next_city >= 1 && next_city <= n_cities tau = pheromone(current_city, next_city)^alpha; eta = (1/distance_matrix(current_city, next_city))^beta; probabilities(i) = tau * eta; end end % 归一化概率 if sum(probabilities) > 0 probabilities = probabilities / sum(probabilities); else probabilities = ones(size(probabilities)) / length(probabilities); end % 轮盘赌选择下一个城市 r = rand(); cumulative_prob = 0; selected_idx = 1; for i = 1:length(probabilities) cumulative_prob = cumulative_prob + probabilities(i); if r <= cumulative_prob selected_idx = i; break; end end % 移动到选中的城市 current_city = unvisited(selected_idx); route = [route, current_city]; unvisited(selected_idx) = []; end % 验证路径 if length(route) ~= n_cities || any(route < 1) || any(route > n_cities) route = 1:n_cities; % 返回默认路径 endendfunction update_pheromone(pheromone, routes, distances, rho, Q, n_cities) % 信息素挥发 pheromone = (1 - rho) * pheromone; % 根据蚂蚁路径添加信息素 n_ants = length(routes); for ant = 1:n_ants if distances(ant) > 0 delta_pheromone = Q / distances(ant); route = routes{ant}; if ~isempty(route) && all(route >= 1) && all(route <= n_cities) n = length(route); for i = 1:n from_city = route(i); to_city = route(mod(i, n) + 1); if from_city >= 1 && from_city <= n_cities && ... to_city >= 1 && to_city <= n_cities pheromone(from_city, to_city) = pheromone(from_city, to_city) + delta_pheromone; pheromone(to_city, from_city) = pheromone(to_city, from_city) + delta_pheromone; end end end end end % 确保对角线为0 pheromone(logical(eye(n_cities))) = 0;end%% intlinprog求解器（实用版本）function [best_route, best_distance, time_taken, success] = solve_tsp_practical(distance_matrix) success = false; best_route = []; best_distance = inf; time_taken = 0; tic; try n = size(distance_matrix, 1); % 多次独立求解 for attempt = 1:10 [f, Aeq, beq, lb, ub, intcon, var_idx] = build_basic_tsp_model(distance_matrix); options = optimoptions('intlinprog', ... 'Display', 'off', ... 'MaxTime', 10, ... 'RelativeGapTolerance',0.1); [x, fval, exitflag] = intlinprog(f, intcon, [], [], Aeq, beq, lb, ub, options); if exitflag > 0 % 验证解的有效性 route = extract_tsp_route(x, var_idx, n); if ~isempty(route) && length(route) == n distance = calculate_total_distance(route, distance_matrix); if distance < best_distance best_distance = distance; best_route = route; success = true; end end end end time_taken = toc; if success fprintf('最优距离: %.2f，耗时: %.4f秒\n', best_distance, time_taken); fprintf('最优路径: '); fprintf('%d ', best_route); fprintf('\n'); else fprintf('求解失败，总耗时: %.4f秒\n', time_taken); end catch ME fprintf('求解器错误: %s\n', ME.message); time_taken = toc; success = false; endendfunction [f, Aeq, beq, lb, ub, intcon, var_idx] = build_basic_tsp_model(distance_matrix) n = size(distance_matrix, 1); % 创建变量索引 var_count = n * (n - 1); f = zeros(var_count, 1); var_idx = zeros(n, n); var_counter = 0; for i = 1:n for j = 1:n if i ~= j var_counter = var_counter + 1; var_idx(i, j) = var_counter; f(var_counter) = distance_matrix(i, j); end end end % 约束条件 Aeq = []; beq = []; % 出度约束 for i = 1:n aeq = zeros(1, var_count); for j = 1:n if i ~= j && var_idx(i, j) > 0 aeq(var_idx(i, j)) = 1; end end Aeq = [Aeq; aeq]; beq = [beq; 1]; end % 入度约束 for j = 1:n aeq = zeros(1, var_count); for i = 1:n if i ~= j && var_idx(i, j) > 0 aeq(var_idx(i, j)) = 1; end end Aeq = [Aeq; aeq]; beq = [beq; 1]; end % 变量边界和整数约束 lb = zeros(var_count, 1); ub = ones(var_count, 1); intcon = 1:var_count;endfunction route = extract_tsp_route(x, var_idx, n) route = []; % 构建邻接矩阵 adj_matrix = zeros(n, n); threshold = 0.5; for i = 1:n for j = 1:n if i ~= j && var_idx(i, j) > 0 && var_idx(i, j) <= length(x) if x(var_idx(i, j)) >= threshold adj_matrix(i, j) = 1; end end end end % 提取路径 current = 1; route = [current]; visited = false(1, n); visited(current) = true; for step = 1:(n-1) next_city = find(adj_matrix(current, :) == 1); if ~isempty(next_city) && ~visited(next_city) current = next_city; route = [route, current]; visited(current) = true; else route = []; return; end endend%% 计算路径总距离function total_distance = calculate_total_distance(route, distance_matrix) if isempty(route) total_distance = inf; return; end n = length(route); if n == 0 total_distance = inf; return; end total_distance = 0; for i = 1:n from_city = route(i); to_city = route(mod(i, n) + 1); % 回到起点 if from_city >= 1 && from_city <= size(distance_matrix, 1) && ... to_city >= 1 && to_city <= size(distance_matrix, 2) dist = distance_matrix(from_city, to_city); if isfinite(dist) total_distance = total_distance + dist; else total_distance = inf; return; end else total_distance = inf; return; end endend
10个月宝宝每天需要喝多少奶粉？