请作为数学建模工程师，解释我给出的代码，把下面的代码转化成数学建模形式，并进行学术化分析，要求保持段落间的逻辑关系。这一部分工作完成后，请作为数学建模编程人员，整理代码中的驱动条件、限制条件，仔细整理出一份详细的仿真计算策略，形成数学建模问题求解的规划策略研究部分。
%% 无人艇最优巡航路径规划 - MPA优化版
clear all; close all; clc;
% 设置随机种子以确保可重复性
rng(42);
%% 系统参数定义
% 网格参数
grid_size = 500; % 网格大小(m)
map_rows = 25; % 行数
map_cols = 48; % 列数
target_pos = [1, 24]; % 目标点位置(行,列)
% 防御区域参数
warning_dist_min = 10; % 预警区最小距离(网格单元) 5km->10
warning_dist_max = 20; % 预警区最大距离(网格单元) 10km->20
intercept_dist_min = 6; % 拦截区最小距离(网格单元) 3km->6
intercept_dist_max = 10;% 拦截区最大距离(网格单元) 5km->10
% 无人艇参数
uav1_speed = 20; % 速度(节)
uav1_range = 400; % 续航(海里)
uav1_detect_radius = 4; % 探测半径(网格单元) 2km->4
uav1_strike_radius = 8; % 打击半径(网格单元) 4km->8
uav1_detect_type = 'semicircle'; % 半圆形探测区域
uav2_speed = 30; % 速度(节)
uav2_range = 300; % 续航(海里)
uav2_detect_radius = 2; % 探测半径(网格单元) 1km->2
uav2_detect_type = 'circle'; % 圆形探测区域
% 进攻船只参数
ship_speed_kn = 10; % 进攻船速(节)，假设值
ship_speed_ms = ship_speed_kn * 0.5144; % 转换为米/秒
time_step_sec = 300; % 时间步长(秒)，5分钟
step_distance = ship_speed_ms * time_step_sec / grid_size; % 每步移动网格数
% 威胁模型参数
k_threat = 10; % 威胁系数
% 路径规划参数
simulation_steps = 100; % 路径长度(时间步数量)
%% MPA优化参数
SearchAgents_no = 10; % 种群数量
Max_iter = 20; % 最大迭代次数
dim = 6; % 优化变量维度 [uav1_row, uav1_col, uav2_row, uav2_col, weight1, weight2]
lb = [1, 1, 1, 1, 0.1, 0.1]; % 下界
ub = [map_rows, map_cols, map_rows, map_cols, 10, 10]; % 上界
%% 加载数据
disp('加载实际数据...');
% 加载小船路径和权重
if exist('100paths_with_weights.mat', 'file')
load('100paths_with_weights.mat', 'all_paths', 'path_weights');
fprintf('成功加载小船路径数据: %d条路径\n', length(all_paths));
% 预处理：确保路径点坐标在合理范围内
for i = 1:length(all_paths)
path = all_paths{i};
path(path(:,1) < 1, 1) = 1;
path(path(:,1) > map_rows, 1) = map_rows;
path(path(:,2) < 1, 2) = 1;
path(path(:,2) > map_cols, 2) = map_cols;
all_paths{i} = path;
end
else
error('小船路径数据文件100paths_with_weights.mat未找到');
end
% 加载无人机路径
if exist('UAV_path.mat', 'file')
load('UAV_path.mat', 'optimal_path1', 'optimal_path2');
fprintf('成功加载无人机巡逻路径\n');
else
error('无人机路径数据文件UAV_path.mat未找到');
end
% 加载水深数据
filename = '附件1 发现概率与海域水深分布表.xlsx';
depth_sheet = '水深分布表';
if exist(filename, 'file')
depth_data = readmatrix(filename, 'Sheet', depth_sheet, 'Range', 'B2:AX26');
valid_grids = depth_data > 2; % 水深大于2米的区域
fprintf('成功加载水深数据\n');
else
error('水深数据文件未找到');
end
%% 初始化和预计算
disp('计算距离地图和威胁热力图...');
% 计算每个网格到目标点的欧氏距离
[X, Y] = meshgrid(1:map_cols, 1:map_rows);
dist_map = sqrt((X - target_pos(2)).^2 + (Y - target_pos(1)).^2) * grid_size; % 实际距离(m)
% 构建无人机覆盖区域
drone_coverage = false(map_rows, map_cols);
% 无人机1覆盖
if exist('optimal_path1', 'var')
for i = 1:size(optimal_path1, 1)
pos = optimal_path1(i, :);
r = round(pos(1)); c = round(pos(2));
if r >= 1 && r <= map_rows && c >= 1 && c <= map_cols
% 无人机探测半径为2km(4个网格)
for dr = max(1, r-4):min(map_rows, r+4)
for dc = max(1, c-4):min(map_cols, c+4)
if norm([dc-c, dr-r]) <= 4
drone_coverage(dr, dc) = true;
end
end
end
end
end
end
% 无人机2覆盖
if exist('optimal_path2', 'var')
for i = 1:size(optimal_path2, 1)
pos = optimal_path2(i, :);
r = round(pos(1)); c = round(pos(2));
if r >= 1 && r <= map_rows && c >= 1 && c <= map_cols
% 无人机探测半径为1.5km(3个网格)
for dr = max(1, r-3):min(map_rows, r+3)
for dc = max(1, c-3):min(map_cols, c+3)
if norm([dc-c, dr-r]) <= 3
drone_coverage(dr, dc) = true;
end
end
end
end
end
end
%% 计算初始威胁地图（静态）
threat_map_static = zeros(map_rows, map_cols);
num_paths = length(all_paths);
path_lengths = zeros(num_paths, 1); % 各路径长度
for i = 1:num_paths
path = all_paths{i};
path_lengths(i) = size(path, 1);
for j = 1:path_lengths(i)
r = path(j,1); c = path(j,2);
if r >= 1 && r <= map_rows && c >= 1 && c <= map_cols
d = dist_map(r, c);
if d > 0
threat_val = k_threat / (d^2) * path_weights(i);
threat_map_static(r, c) = threat_map_static(r, c) + threat_val;
end
end
end
end
% 无人机覆盖区域威胁降低90%
threat_map_static(drone_coverage) = threat_map_static(drone_coverage) * 0.1;
%% 为MPA目标函数创建参数结构体
simParams.dist_map = dist_map;
simParams.valid_grids = valid_grids;
simParams.all_paths = all_paths;
simParams.path_weights = path_weights;
simParams.num_paths = num_paths;
simParams.path_lengths = path_lengths;
simParams.simulation_steps = simulation_steps;
simParams.step_distance = step_distance;
simParams.threat_map_static = threat_map_static;
simParams.uav1_detect_radius = uav1_detect_radius;
simParams.uav2_detect_radius = uav2_detect_radius;
simParams.uav1_detect_type = uav1_detect_type;
simParams.uav2_detect_type = uav2_detect_type;
simParams.map_rows = map_rows;
simParams.map_cols = map_cols;
simParams.k_threat = k_threat;
simParams.warning_min_m = 5000;
simParams.warning_max_m = 10000;
simParams.intercept_min_m = 3000;
simParams.intercept_max_m = 5000;
%% 定义目标函数（用于MPA优化）
function final_threat = mpa_objective_function(x, params)
% 解析优化变量
uav1_start = [x(1), x(2)];
uav2_start = [x(3), x(4)];
weight1 = x(5);
weight2 = x(6);
% 检查位置有效性
if ~is_position_valid(uav1_start, params.dist_map, params.valid_grids, ...
params.warning_min_m, params.warning_max_m) || ...
~is_position_valid(uav2_start, params.dist_map, params.valid_grids, ...
params.intercept_min_m, params.intercept_max_m)
final_threat = 1e10; % 无效位置给予高惩罚
return;
end
% 初始化船只状态
ship_positions = cell(params.num_paths, 1);
threat_ratios = ones(params.num_paths, 1);
for i = 1:params.num_paths
path = params.all_paths{i};
ship_positions{i} = [path(1,1), path(1,2), 1];
end
% 初始化无人艇路径
uav1_path = zeros(params.simulation_steps, 2);
uav2_path = zeros(params.simulation_steps, 2);
uav1_path(1, :) = uav1_start;
uav2_path(1, :) = uav2_start;
% 初始方向
uav1_direction = [0, 1];
uav2_direction = [0, 1];
% 方向向量（右, 下, 左, 上）
directions = [0, 1; 1, 0; 0, -1; -1, 0];
% 生成路径的每一步
threat_history = zeros(params.simulation_steps, 1);
for step = 1:params.simulation_steps
% 更新进攻船只位置
active_threat = 0;
for i = 1:params.num_paths
pos_info = ship_positions{i};
current_index = pos_info(3);
if current_index < params.path_lengths(i)
move_progress = min(params.step_distance, params.path_lengths(i) - current_index);
new_index = floor(current_index + move_progress);
path = params.all_paths{i};
new_pos = path(new_index, :);
ship_positions{i} = [new_pos(1), new_pos(2), new_index];
end
r = ship_positions{i}(1); c = ship_positions{i}(2);
d = params.dist_map(r, c);
if d > 0
current_threat = params.k_threat / (d^2) * params.path_weights(i) * threat_ratios(i);
active_threat = active_threat + current_threat;
end
end
threat_history(step) = active_threat;
% 无人艇移动
if step > 1
% UAV1: 在预警区巡航
[uav1_next, uav1_dir] = move_agent_smart(uav1_path(step-1, :), uav1_direction, ...
params.threat_map_static, directions, params.map_rows, params.map_cols, params.dist_map, ...
params.valid_grids, params.uav1_detect_radius, params.uav1_detect_type, ...
params.warning_min_m, params.warning_max_m, weight1);
uav1_path(step, :) = uav1_next;
uav1_direction = uav1_dir;
% UAV2: 在拦截区巡航
[uav2_next, uav2_dir] = move_agent_smart(uav2_path(step-1, :), uav2_direction, ...
params.threat_map_static, directions, params.map_rows, params.map_cols, params.dist_map, ...
params.valid_grids, params.uav2_detect_radius, params.uav2_detect_type, ...
params.intercept_min_m, params.intercept_max_m, weight2);
uav2_path(step, :) = uav2_next;
uav2_direction = uav2_dir;
end
% 更新覆盖区域
[uav1_detect, ~] = compute_coverage(uav1_path(step, :), uav1_direction, ...
params.uav1_detect_radius, params.uav1_detect_type, params.map_rows, params.map_cols);
[uav2_detect, ~] = compute_coverage(uav2_path(step, :), uav2_direction, ...
params.uav2_detect_radius, params.uav2_detect_type, params.map_rows, params.map_cols);
% 察打一体无人艇执行打击
if any(uav1_detect(:) > 0)
for i = 1:params.num_paths
pos = [ship_positions{i}(1), ship_positions{i}(2)];
if uav1_detect(pos(1), pos(2)) && rand() < 0.8 % 80%打击成功率
threat_ratios(i) = threat_ratios(i) * 0.1; % 威胁降低90%
end
end
end
% 高速巡逻艇执行拦截
if any(uav2_detect(:) > 0)
for i = 1:params.num_paths
pos = [ship_positions{i}(1), ship_positions{i}(2)];
if uav2_detect(pos(1), pos(2)) && rand() < 0.7 % 70%拦截成功率
threat_ratios(i) = threat_ratios(i) * 0.2; % 威胁降低80%
end
end
end
end
final_threat = threat_history(end);
end
%% 使用MPA优化无人艇位置和参数
disp('使用海洋捕食者算法(MPA)优化无人艇配置...');
fobj = @(x) mpa_objective_function(x, simParams);
[Top_predator_fit, Top_predator_pos, Convergence_curve] = MPA(SearchAgents_no, Max_iter, lb, ub, dim, fobj);
% 解析最优解
uav1_start = [Top_predator_pos(1), Top_predator_pos(2)];
uav2_start = [Top_predator_pos(3), Top_predator_pos(4)];
weight1 = Top_predator_pos(5);
weight2 = Top_predator_pos(6);
fprintf('MPA优化完成! 最优威胁值: %.4f\n', Top_predator_fit);
fprintf('UAV1初始位置: [%d, %d], 权重: %.2f\n', uav1_start(1), uav1_start(2), weight1);
fprintf('UAV2初始位置: [%d, %d], 权重: %.2f\n', uav2_start(1), uav2_start(2), weight2);
%% 使用最优配置运行完整仿真
disp('使用最优配置运行完整仿真...');
% 初始化船只状态
ship_positions = cell(num_paths, 1);
threat_ratios = ones(num_paths, 1);
for i = 1:num_paths
path = all_paths{i};
ship_positions{i} = [path(1,1), path(1,2), 1];
end
% 初始化无人艇路径
uav1_path = zeros(simulation_steps, 2);
uav2_path = zeros(simulation_steps, 2);
uav1_path(1, :) = uav1_start;
uav2_path(1, :) = uav2_start;
% 存储覆盖区域
uav1_coverage = false(map_rows, map_cols);
uav2_coverage = false(map_rows, map_cols);
% 存储打击记录
strike_positions = [];
intercept_positions = [];
% 初始方向
uav1_direction = [0, 1];
uav2_direction = [0, 1];
% 方向向量（右, 下, 左, 上）
directions = [0, 1; 1, 0; 0, -1; -1, 0];
% 总威胁跟踪
threat_history = zeros(simulation_steps, 1);
% 生成路径的每一步
for step = 1:simulation_steps
% 更新进攻船只位置
active_threat = 0;
for i = 1:num_paths
pos_info = ship_positions{i};
current_index = pos_info(3);
if current_index < path_lengths(i)
move_progress = min(step_distance, path_lengths(i) - current_index);
new_index = floor(current_index + move_progress);
path = all_paths{i};
new_pos = path(new_index, :);
ship_positions{i} = [new_pos(1), new_pos(2), new_index];
end
r = ship_positions{i}(1); c = ship_positions{i}(2);
d = dist_map(r, c);
if d > 0
current_threat = k_threat / (d^2) * path_weights(i) * threat_ratios(i);
active_threat = active_threat + current_threat;
end
end
threat_history(step) = active_threat;
% 无人艇移动
if step > 1
% UAV1: 在预警区巡航
[uav1_next, uav1_dir] = move_agent_smart(uav1_path(step-1, :), uav1_direction, ...
threat_map_static, directions, map_rows, map_cols, dist_map, valid_grids, ...
uav1_detect_radius, uav1_detect_type, 5000, 10000, weight1);
uav1_path(step, :) = uav1_next;
uav1_direction = uav1_dir;
% UAV2: 在拦截区巡航
[uav2_next, uav2_dir] = move_agent_smart(uav2_path(step-1, :), uav2_direction, ...
threat_map_static, directions, map_rows, map_cols, dist_map, valid_grids, ...
uav2_detect_radius, uav2_detect_type, 3000, 5000, weight2);
uav2_path(step, :) = uav2_next;
uav2_direction = uav2_dir;
end
% 更新覆盖区域
[uav1_detect, uav1_cells] = compute_coverage(uav1_path(step, :), uav1_direction, ...
uav1_detect_radius, uav1_detect_type, map_rows, map_cols);
uav1_coverage = uav1_coverage | uav1_cells;
[uav2_detect, uav2_cells] = compute_coverage(uav2_path(step, :), uav2_direction, ...
uav2_detect_radius, uav2_detect_type, map_rows, map_cols);
uav2_coverage = uav2_coverage | uav2_cells;
% 察打一体无人艇执行打击
if any(uav1_detect(:) > 0)
for i = 1:num_paths
pos = [ship_positions{i}(1), ship_positions{i}(2)];
if uav1_detect(pos(1), pos(2)) && rand() < 0.8 % 80%打击成功率
threat_ratios(i) = threat_ratios(i) * 0.1; % 威胁降低90%
strike_positions = [strike_positions; pos, step];
end
end
end
% 高速巡逻艇执行拦截
if any(uav2_detect(:) > 0)
for i = 1:num_paths
pos = [ship_positions{i}(1), ship_positions{i}(2)];
if uav2_detect(pos(1), pos(2)) && rand() < 0.7 % 70%拦截成功率
threat_ratios(i) = threat_ratios(i) * 0.2; % 威胁降低80%
intercept_positions = [intercept_positions; pos, step];
end
end
end
end
%% 创建网格地图展示
disp('创建网格地图展示...');
% 归一化威胁图用于可视化
threat_normalized = (threat_map_static - min(threat_map_static(:))) / ...
(max(threat_map_static(:)) - min(threat_map_static(:)));
% 创建地图
create_enhanced_grid_map(uav1_path, uav2_path, threat_normalized, ...
warning_dist_min, warning_dist_max, intercept_dist_min, intercept_dist_max, ...
target_pos, map_rows, map_cols, ...
uav1_detect_radius, uav2_detect_radius, uav1_coverage, uav2_coverage, ...
strike_positions, intercept_positions, threat_history, time_step_sec, ...
Top_predator_fit, Convergence_curve);
%% 威胁统计与输出
final_threat = threat_history(end);
uav1_cover_percent = sum(uav1_coverage(:)) / (map_rows * map_cols) * 100;
uav2_cover_percent = sum(uav2_coverage(:)) / (map_rows * map_cols) * 100;
disp(' ');
fprintf('=== 威胁统计 ===\n');
fprintf('最终剩余威胁: %.2f\n', final_threat);
fprintf('UAV1覆盖率: %.1f%%\n', uav1_cover_percent);
fprintf('UAV2覆盖率: %.1f%%\n', uav2_cover_percent);
fprintf('联合覆盖率: %.1f%%\n', sum(uav1_coverage(:) | uav2_coverage(:)) / (map_rows*map_cols)*100);
%% ============= 辅助函数 ============= %
function valid = is_position_valid(pos, dist_map, valid_grids, min_dist, max_dist)
r = round(pos(1));
c = round(pos(2));
[rows, cols] = size(dist_map);
% 检查位置是否在网格范围内
if r < 1 || r > rows || c < 1 || c > cols
valid = false;
return;
end
% 检查水深和距离限制
d = dist_map(r, c);
valid = (d >= min_dist) && (d <= max_dist) && valid_grids(r, c);
end
function [detect_grid, coverage] = compute_coverage(pos, direction, radius, detect_type, rows, cols)
% 计算覆盖区域，返回逻辑矩阵
detect_grid = false(rows, cols);
[X, Y] = meshgrid(1:cols, 1:rows);
dist = sqrt((X - pos(2)).^2 + (Y - pos(1)).^2);
if strcmp(detect_type, 'semicircle')
% 计算方向向量角度
if all(direction == [0,0])
angle = 0; % 默认向右
else
angle = atan2(direction(1), direction(2));
end
% 计算各点相对于中心的角度
angles = atan2(Y - pos(1), X - pos(2)) - angle;
angles = mod(angles + pi, 2*pi) - pi; % 标准化到[-π,π]
% 创建半圆覆盖区域 (前向180度)
detect_grid = (dist <= radius) & (abs(angles) <= pi/2);
elseif strcmp(detect_type, 'circle')
% 圆形覆盖
detect_grid = dist <= radius;
end
% 边界处理
detect_grid = detect_grid & (X >= 1) & (X <= cols) & (Y >= 1) & (Y <= rows);
coverage = detect_grid;
end
function [new_pos, new_dir] = move_agent_smart(current_pos, current_dir, threat_map, ...
directions, rows, cols, dist_map, valid_grids, detect_radius, detect_type, ...
min_dist_m, max_dist_m, weight)
% 智能移动代理：选择最大威胁覆盖的方向
best_score = -inf;
new_pos = current_pos; % 默认不移动
new_dir = current_dir; % 默认方向不变
% 尝试所有可能方向 + 当前位置
candidate_positions = [current_pos; current_pos + directions];
[num_candidates, ~] = size(candidate_positions);
for i = 1:num_candidates
pos = candidate_positions(i, :);
r = round(pos(1));
c = round(pos(2));
% 检查位置有效性
if r < 1 || r > rows || c < 1 || c > cols
continue;
end
% 检查水深和任务区域
dist_val = dist_map(r, c);
if ~valid_grids(r, c) || dist_val < min_dist_m || dist_val > max_dist_m
continue;
end
% 确定方向 (如果移动)
if i == 1 % 当前位置
dir = current_dir;
else
dir = directions(i-1, :);
end
% 计算在该位置的探测覆盖率
[~, coverage] = compute_coverage([r, c], dir, detect_radius, detect_type, rows, cols);
coverage_score = sum(threat_map(coverage), 'all');
% 距离保持分数 (保持在任务区域中心附近)
ideal_dist = (min_dist_m + max_dist_m)/2;
dist_score = 1/(1 + abs(dist_val - ideal_dist)/1000); % 归一化
% 总分数 = 覆盖威胁 + 距离保持
score = coverage_score + weight * dist_score;
% 更新最佳位置
if score > best_score
best_score = score;
new_pos = [r, c];
new_dir = dir;
end
end
end
function create_enhanced_grid_map(uav1_path, uav2_path, threat_normalized, ...
warn_min, warn_max, int_min, int_max, ...
target, rows, cols, uav1_radius, uav2_radius, ...
uav1_coverage, uav2_coverage, strike_pos, intercept_pos, ...
threat_history, time_step, best_fit, conv_curve)
% 增强的网格地图展示
fig = figure('Position', [100, 100, 1200, 900], 'Name', '无人艇巡航路径与威胁覆盖');
% 主路径图
subplot(2,3,[1,2,4,5]);
hold on;
axis equal;
grid on;
box on;
xlim([0.5, cols+0.5]);
ylim([0.5, rows+0.5]);
set(gca, 'YDir', 'reverse', 'FontSize', 10);
xticks(1:cols);
yticks(1:rows);
xlabel('列号', 'FontSize', 12);
ylabel('行号', 'FontSize', 12);
title(sprintf('无人艇巡航路径 (最终威胁: %.2f)', threat_history(end)), 'FontSize', 14);
% 绘制威胁热力图背景
imagesc(threat_normalized, 'AlphaData', threat_normalized > 0);
colormap(jet);
colorbar('southoutside');
caxis([0, 1]);
% 标记防御区域
[X, Y] = meshgrid(1:cols, 1:rows);
dist_to_target = sqrt((X - target(2)).^2 + (Y - target(1)).^2);
warn_zone = (dist_to_target >= warn_min) & (dist_to_target <= warn_max);
int_zone = (dist_to_target >= int_min) & (dist_to_target <= int_max);
[warn_y, warn_x] = find(warn_zone);
[int_y, int_x] = find(int_zone);
plot(warn_x, warn_y, 'y.', 'MarkerSize', 4, 'DisplayName', '预警区');
plot(int_x, int_y, 'g.', 'MarkerSize', 4, 'DisplayName', '拦截区');
% 标记目标点
plot(target(2), target(1), 'rp', 'MarkerSize', 25, 'MarkerFaceColor', 'r');
text(target(2), target(1), ' 目标', 'Color', 'r', 'FontSize', 14, ...
'VerticalAlignment', 'middle', 'FontWeight', 'bold');
% 绘制无人艇路径
plot(uav1_path(:, 2), uav1_path(:, 1), 'b-', 'LineWidth', 1.5);
plot(uav1_path(:, 2), uav1_path(:, 1), 'b.', 'MarkerSize', 8);
plot(uav1_path(1, 2), uav1_path(1, 1), 'bo', 'MarkerSize', 12, 'MarkerFaceColor', 'b');
plot(uav1_path(end, 2), uav1_path(end, 1), 'bd', 'MarkerSize', 12, 'MarkerFaceColor', 'b');
plot(uav2_path(:, 2), uav2_path(:, 1), 'm-', 'LineWidth', 1.5);
plot(uav2_path(:, 2), uav2_path(:, 1), 'm.', 'MarkerSize', 8);
plot(uav2_path(1, 2), uav2_path(1, 1), 'mo', 'MarkerSize', 12, 'MarkerFaceColor', 'm');
plot(uav2_path(end, 2), uav2_path(end, 1), 'md', 'MarkerSize', 12, 'MarkerFaceColor', 'm');
% 绘制打击和拦截位置
if ~isempty(strike_pos)
plot(strike_pos(:,2), strike_pos(:,1), 'rx', 'MarkerSize', 10, 'LineWidth', 2, 'DisplayName', '打击点');
end
if ~isempty(intercept_pos)
plot(intercept_pos(:,2), intercept_pos(:,1), 'kx', 'MarkerSize', 10, 'LineWidth', 1, 'DisplayName', '拦截点');
end
% 添加图例
legend_items = {'目标点', '预警区', '拦截区', ...
'察打艇路径', '巡逻艇路径', '打击点', '拦截点'};
valid_items = isgraphics(findobj(gca, 'Type', 'line'));
legend(legend_items(1:min(length(valid_items), length(legend_items))), 'Location', 'southoutside', 'NumColumns', 4, 'FontSize', 9);
% 创建子图：威胁覆盖图
subplot(2,3,3);
hold on;
imagesc(threat_normalized);
colormap(jet);
colorbar;
% 覆盖区域边界
[b_y1, b_x1] = find(bwperim(uav1_coverage));
[b_y2, b_x2] = find(bwperim(uav2_coverage));
plot(b_x1, b_y1, 'b.', 'MarkerSize', 4, 'DisplayName', 'UAV1覆盖');
plot(b_x2, b_y2, 'm.', 'MarkerSize', 4, 'DisplayName', 'UAV2覆盖');
plot(target(2), target(1), 'rp', 'MarkerSize', 20, 'MarkerFaceColor', 'r');
title('最终威胁分布与覆盖', 'FontSize', 12);
axis equal;
xlim([0.5, cols+0.5]);
ylim([0.5, rows+0.5]);
set(gca, 'YDir', 'reverse', 'FontSize', 9);
xlabel('列号');
ylabel('行号');
legend('Location', 'southoutside', 'Orientation', 'horizontal');
% 创建子图：威胁减少过程
subplot(2,3,6);
time_vec = (0:length(threat_history)-1) * time_step / 60; % 转换为分钟
plot(time_vec, threat_history, 'b-', 'LineWidth', 2);
grid on;
xlabel('时间 (分钟)', 'FontSize', 10);
ylabel('实时威胁指数', 'FontSize', 10);
title('威胁随时间减少过程', 'FontSize', 12);
% 添加网格线
min_time = 0;
max_time = time_vec(end);
time_ticks = linspace(min_time, max_time, 6);
set(gca, 'XTick', time_ticks);
% 添加MPA收敛曲线
yyaxis right
plot(1:length(conv_curve), conv_curve, 'r-', 'LineWidth', 1.5);
ylabel('MPA收敛曲线', 'FontSize', 10);
title(sprintf('MPA优化结果: 最优威胁=%.2f', best_fit), 'FontSize', 12);
legend('威胁变化', 'MPA收敛', 'Location', 'southeast');
% 保存图像
saveas(fig, 'unmanned_boat_optimal_paths_MPA.png');
disp('巡航路径图已保存为 unmanned_boat_optimal_paths_MPA.png');
end
%% ============= MPA算法实现 ============= %
function [Top_predator_fit,Top_predator_pos,Convergence_curve]=MPA(SearchAgents_no,Max_iter,lb,ub,dim,fobj)
Top_predator_pos=zeros(1,dim);
Top_predator_fit=inf;
stepsize=zeros(SearchAgents_no,dim);
fitness=inf(SearchAgents_no,1);
% 初始化种群
Prey=initialization(SearchAgents_no,dim,ub,lb);
Xmin=repmat(lb,SearchAgents_no,1);
Xmax=repmat(ub,SearchAgents_no,1);
Iter=0;
FADs=0.2;
P=0.5;
while Iter<Max_iter
%------------------- Detecting top predator -----------------
for i=1:size(Prey,1)
Flag4ub=Prey(i,:)>ub;
Flag4lb=Prey(i,:)<lb;
Prey(i,:)=(Prey(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
fitness(i,1)=fobj(Prey(i,:));
if fitness(i,1)<Top_predator_fit
Top_predator_fit=fitness(i,1);
Top_predator_pos=Prey(i,:);
end
end
%------------------- Marine Memory saving -------------------
if Iter==0
fit_old=fitness; Prey_old=Prey;
end
Inx=(fit_old<fitness);
Indx=repmat(Inx,1,dim);
Prey=Indx.*Prey_old+~Indx.*Prey;
fitness=Inx.*fit_old+~Inx.*fitness;
fit_old=fitness; Prey_old=Prey;
%------------------------------------------------------------
Elite=repmat(Top_predator_pos,SearchAgents_no,1); %(Eq. 10)
CF=(1-Iter/Max_iter)^(2*Iter/Max_iter);
RL=0.05*levy(SearchAgents_no,dim,1.5); %Levy random number vector
RB=randn(SearchAgents_no,dim); %Brownian random number vector
for i=1:size(Prey,1)
for j=1:size(Prey,2)
R=rand();
%------------------ Phase 1 (Eq.12) -------------------
if Iter<Max_iter/3
stepsize(i,j)=RB(i,j)*(Elite(i,j)-RB(i,j)*Prey(i,j));
Prey(i,j)=Prey(i,j)+P*R*stepsize(i,j);
%--------------- Phase 2 (Eqs. 13 & 14)----------------
elseif Iter>Max_iter/3 && Iter<2*Max_iter/3
if i>size(Prey,1)/2
stepsize(i,j)=RB(i,j)*(RB(i,j)*Elite(i,j)-Prey(i,j));
Prey(i,j)=Elite(i,j)+P*CF*stepsize(i,j);
else
stepsize(i,j)=RL(i,j)*(Elite(i,j)-RL(i,j)*Prey(i,j));
Prey(i,j)=Prey(i,j)+P*R*stepsize(i,j);
end
%----------------- Phase 3 (Eq. 15)-------------------
else
stepsize(i,j)=RL(i,j)*(RL(i,j)*Elite(i,j)-Prey(i,j));
Prey(i,j)=Elite(i,j)+P*CF*stepsize(i,j);
end
end
end
%------------------ Detecting top predator ------------------
for i=1:size(Prey,1)
Flag4ub=Prey(i,:)>ub;
Flag4lb=Prey(i,:)<lb;
Prey(i,:)=(Prey(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
fitness(i,1)=fobj(Prey(i,:));
if fitness(i,1)<Top_predator_fit
Top_predator_fit=fitness(i,1);
Top_predator_pos=Prey(i,:);
end
end
%---------------------- Marine Memory saving ----------------
if Iter==0
fit_old=fitness; Prey_old=Prey;
end
Inx=(fit_old<fitness);
Indx=repmat(Inx,1,dim);
Prey=Indx.*Prey_old+~Indx.*Prey;
fitness=Inx.*fit_old+~Inx.*fitness;
fit_old=fitness; Prey_old=Prey;
%---------- Eddy formation and FADs effect (Eq 16) -----------
if rand()<FADs
U=rand(SearchAgents_no,dim)<FADs;
Prey=Prey+CF*((Xmin+rand(SearchAgents_no,dim).*(Xmax-Xmin)).*U);
else
r=rand(); Rs=size(Prey,1);
stepsize=(FADs*(1-r)+r)*(Prey(randperm(Rs),:)-Prey(randperm(Rs),:));
Prey=Prey+stepsize;
end
Iter=Iter+1;
Convergence_curve(Iter)=Top_predator_fit;
end
end
function Positions=initialization(SearchAgents_no,dim,ub,lb)
Boundary_no= size(ub,2); % numnber of boundaries
% If the boundaries of all variables are equal and user enter a signle
% number for both ub and lb
if Boundary_no==1
Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;
end
% If each variable has a different lb and ub
if Boundary_no>1
for i=1:dim
ub_i=ub(i);
lb_i=lb(i);
Positions(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
end
end
end
function [z] = levy(n,m,beta)
num = gamma(1+beta)*sin(pi*beta/2); % used for Numerator
den = gamma((1+beta)/2)*beta*2^((beta-1)/2); % used for Denominator
sigma_u = (num/den)^(1/beta);% Standard deviation
u = random('Normal',0,sigma_u,n,m);
v = random('Normal',0,1,n,m);
z =u./(abs(v).^(1/beta));
end
博客内容仅提供了一个链接，未包含其他关键信息。链接为http://wangyalei.iteye.com/blog/744383 。
扫一扫
2233
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
