%% 船舶航线优化主程序
clear; clc; close all;
%% 参数设置
% 船舶参数(Suezmax型)
ship.L_pp = 228.9; % 垂线间长(m)
ship.B = 32.26; % 船宽(m)
ship.T = 14.471; % 吃水(m)
ship.Cb = 0.803; % 方形系数
ship.D = 8.5; % 螺旋桨直径(m)
ship.sfoc = 104; % 燃油消耗率(g/kWh)
ship.V = 12; % 设计航速(kn)
ship.V_max = 16; % 最大航速(kn)
ship.V_min = 10; % 最小航速(kn)
ship.P_max = 32000; % 最大主机功率(kW)
ship.rho_sea = 1026; % 海水密度(kg/m^3)
ship.rho_air = 1.225; % 空气密度(kg/m^3)
ship.A_T = 1918; % 水面上的最大横向面积(m^2)
% 算法参数
population_size = 50; % 种群大小
max_generations = 100; % 最大迭代次数
pc = 0.8; % 交叉概率
pm = 0.1; % 变异概率
T0 = 100; % 初始温度
alpha = 0.95; % 退火系数
start_time = datetime(2025,2,12,18,0,0);
% 航线参数
start_port = [35.2391, 119.7921]; % 日照港
end_port = [-28.5202, 114.3272]; % 杰拉尔顿港
max_segments = 30; % 最大航段数
lat_range = [-30, 55];
lon_range = [110, 160];
%% 气象数据加载与预处理
% 读取NetCDF文件
op_ncfile = 'C:\Users\dell\Downloads\41c65edbd5a653750e1d93b2a34fe6a8.nc';
wv_ncfile = 'C:\Users\dell\Downloads\abbf6737917f8802584d2b89ac9cf3ca.nc';
% 获取维度信息
op_lon = ncread(op_ncfile, 'longitude'); % 风场文件经度向量
op_lat = ncread(op_ncfile, 'latitude'); % 风场文件纬度向量
time = ncread(op_ncfile, 'valid_time'); % 时间(小时数,从1900-01-01起)
wv_lon = ncread(wv_ncfile, 'longitude'); % 浪场文件经度向量
wv_lat = ncread(wv_ncfile, 'latitude'); % 浪场文件纬度向量
% 正确时间转换(基于NC文件元数据)
time_ref = datetime(1970,1,1, 'TimeZone', 'UTC'); % 明确指定UTC时区
date_vec = time_ref + seconds(time); % 保持UTC时区
% 提取2025年1-2月数据
start_date = datetime(2025,2,12, 'TimeZone', 'UTC');
end_date = datetime(2025,2,27,23,59,59, 'TimeZone', 'UTC');
time_mask = (date_vec >= start_date) & (date_vec <= end_date);
time_idx = find(time_mask);
% 读取气象变量(优化维度处理)
% 注意:NetCDF维度顺序通常为(longitude, latitude, time)
u10 = permute(ncread(op_ncfile, 'u10', [1,1,1], [Inf, Inf, Inf]), [2,1,3]); % 调整为(lat,lon,time)
v10 = permute(ncread(op_ncfile, 'v10', [1,1,1], [Inf, Inf, Inf]), [2,1,3]);
swh = permute(ncread(wv_ncfile, 'swh', [1,1,1], [Inf, Inf, Inf]), [2,1,3]);
cdww = permute(ncread(wv_ncfile, 'cdww', [1,1,1], [Inf, Inf, Inf]), [2,1,3]);
% 创建纬度-经度网格
[LonGrid, LatGrid] = meshgrid(op_lon, op_lat);
% 构建气象数据结构体
weather_data = struct(...
'op_lon', op_lon, ... % 经度向量 [nLon×1]
'op_lat', op_lat, ... % 纬度向量 [nLat×1]
'wv_lon', wv_lon, ... % 经度向量 [nLon×1]
'wv_lat', wv_lat, ... % 纬度向量 [nLat×1]
'LonGrid', LonGrid, ... % 经度网格 [nLat×nLon]
'LatGrid', LatGrid, ... % 纬度网格 [nLat×nLon]
'time', date_vec(time_idx), ... % 时间向量 [nTime×1]
'u10', u10(:,:,time_idx), ... % 风速U分量 [nLat×nLon×nTime]
'v10', v10(:,:,time_idx), ... % 风速V分量 [nLat×nLon×nTime]
'swh', swh(:,:,time_idx), ... % 有效波高 [nLat×nLon×nTime]
'cdww', cdww(:,:,time_idx)); % 波浪阻力系数 [nLat×nLon×nTime]
%% 初始化种群
population = initialize_population(population_size, start_port, end_port, max_segments, lat_range, lon_range);
%% 混合遗传算法主循环
best_fitness = inf;
best_route = [];
best_fuel = Inf;
best_time = Inf;
fitness_history = zeros(max_generations, 1);
temperature = T0;
for gen = 1:max_generations
% 计算适应度
fitness = zeros(population_size, 1);
fuel_consumptions = zeros(population_size, 1);
total_times = zeros(population_size, 1);
for i = 1:population_size
[fitness(i), fuel_consumptions(i), total_times(i)] = calculate_fitness(population{i}, ship, weather_data, start_time);
end
% 记录最佳个体
[min_fitness, idx] = min(fitness);
current_fuel = fuel_consumptions(idx);
current_time = total_times(idx);
if min_fitness < best_fitness
best_fitness = min_fitness;
best_route = population{idx};
best_fuel = current_fuel;
best_time = current_time;
end
fitness_history(gen) = best_fitness;
% 选择操作
selected = tournament_selection(population, fitness, population_size/2);
% 交叉操作
offspring = hybrid_crossover(selected, pc, start_port, end_port);
% 变异操作
offspring = mutation(offspring, pm, lat_range, lon_range);
% 合并种群
combined_population = [population; offspring];
offspring_fitness = zeros(length(offspring), 1);
for i = 1:length(offspring)
offspring_fitness(i) = calculate_fitness(offspring{i}, ship, weather_data, start_time);
end
combined_fitness = [fitness; offspring_fitness];
% 模拟退火选择
population = simulated_annealing_selection(combined_population, combined_fitness, population_size, temperature);
% 更新温度
temperature = alpha * temperature;
% 显示进度
fprintf('Generation %d: Fitness=%.2f, Fuel=%.2fkg, Time=%.2fh (Current) | Best: Fitness=%.2f, Fuel=%.2fkg, Time=%.2fh\n',...
gen, min_fitness, current_fuel, current_time, best_fitness, best_fuel, best_time);
% 更新温度
temperature = alpha * temperature;
end
%% 结果可视化
plot_route(best_route, start_port, end_port, weather_data);
figure;
plot(1:max_generations, fitness_history);
xlabel('Generation');
ylabel('Best Fitness');
title('Convergence Curve');
%% 初始化种群函数
function population = initialize_population(pop_size, start_port, end_port, max_segments, lat_range, lon_range)
population = cell(pop_size, 1); % 初始化为列向量
for i = 1:pop_size
num_segments = randi([2, max_segments]);
route = zeros(num_segments+1, 2);
route(1,:) = start_port;
route(end,:) = end_port;
% 生成中间点
for j = 2:num_segments
while true
new_lat = lat_range(1) + (lat_range(2)-lat_range(1))*rand();
new_lon = lon_range(1) + (lon_range(2)-lon_range(1))*rand();
if ~isequal([new_lat, new_lon], start_port) && ~isequal([new_lat, new_lon], end_port)
route(j,:) = [new_lat, new_lon];
break;
end
end
end
population{i} = route;
% 按经度排序中间点,确保东向航行
route(2:end-1,:) = sortrows(route(2:end-1,:), 2);
end
end
%% 适应度计算函数
function [fitness, total_fuel, total_time] = calculate_fitness(route, ship, weather_data, start_time)
% 检查路径有效性
if isempty(route) || size(route,1) < 2 || any(isnan(route(:)))
fitness = Inf;
return;
end
total_distance = 0;
total_time = 0;
total_fuel = 0;
% 计算每个航段的参数
for i = 1:size(route,1)-1
% 计算航段距离 (大圆距离)
[dist, ~] = lldistkm(route(i,:), route(i+1,:));
if isnan(dist) % 如果距离计算失败
fitness = Inf;
return;
end
total_distance = total_distance + dist;
% 获取航段中点位置和时间的气象数据
mid_point = [(route(i,1)+route(i+1,1))/2, (route(i,2)+route(i+1,2))/2];
segment_time = start_time + seconds(total_time*3600);
segment_time.TimeZone = 'UTC'; % 强制设置时区
% 双线性插值获取气象数据
[u10, v10, swh, cdww] = get_weather_data(mid_point, segment_time, weather_data);
% 计算航向角 (0-360度,正北为0)
heading = calculate_heading(route(i,:), route(i+1,:));
% 计算相对风向和波向
wind_dir = atan2d(v10, u10); % 风向(度)
wind_speed = sqrt(u10^2 + v10^2);
rel_wind_angle = wrapTo360(wind_dir - heading);
% 计算船舶实际速度 (考虑风浪影响)
[V_actual, fuel_consumption] = calculate_ship_performance(ship, wind_speed, rel_wind_angle, swh, cdww);
% 计算航段时间和燃油消耗
segment_time_hours = dist / (V_actual * 1.852); % 海里转公里
total_time = total_time + segment_time_hours;
total_fuel = total_fuel + fuel_consumption * segment_time_hours;
% 添加航速有效性检查
if V_actual <= 0 || isnan(V_actual)
fitness = Inf;
return;
end
% 添加时间计算保护
segment_time_hours = max(0, dist / (V_actual * 1.852));
end
% 适应度函数 (加权和)
fuel_weight = 0.6; % 燃油消耗权重
fitness = fuel_weight * total_fuel + (1-fuel_weight) * total_time;
fprintf('Segment %d: V=%.2fkn, Fuel=%.2fkg/h, Time=%.2fh\n',...
i, V_actual, fuel_consumption, segment_time);
end
%% 气象数据获取函数
function [u10, v10, swh, cdww] = get_weather_data(point, time, weather_data)
% 双线性插值获取气象数据
% 风场数据 (0.25°分辨率)
[u10, v10] = bilinear_interp(point(1), point(2), time, ...
weather_data.op_lat, weather_data.op_lon, weather_data.time, ...
weather_data.u10, weather_data.v10);
% 浪场数据 (0.5°分辨率)
[swh, cdww] = bilinear_interp(point(1), point(2), time, ...
weather_data.wv_lat, weather_data.wv_lon, weather_data.time, ...
weather_data.swh, weather_data.cdww);
end
%% 双线性插值法函数
function [varargout] = bilinear_interp(lat, lon, time, lat_grid, lon_grid, time_grid, varargin)
% 通用双线性插值函数
varargout = cell(1, nargin-6);
% 找到最近的4个空间点和2个时间点
[lat_idx, lon_idx, time_idx, weights] = find_nearest_points(lat, lon, time, lat_grid, lon_grid, time_grid);
for i = 1:length(varargin)
var = varargin{i};
% 空间插值
spatial_interp = (1-weights(1))*(1-weights(2)) * var(lat_idx(1), lon_idx(1), time_idx(1)) + ...
(1-weights(1))*weights(2) * var(lat_idx(1), lon_idx(2), time_idx(1)) + ...
weights(1)*(1-weights(2)) * var(lat_idx(2), lon_idx(1), time_idx(1)) + ...
weights(1)*weights(2) * var(lat_idx(2), lon_idx(2), time_idx(1));
% 时间插值
if length(time_idx) > 1
spatial_interp2 = (1-weights(1))*(1-weights(2)) * var(lat_idx(1), lon_idx(1), time_idx(2)) + ...
(1-weights(1))*weights(2) * var(lat_idx(1), lon_idx(2), time_idx(2)) + ...
weights(1)*(1-weights(2)) * var(lat_idx(2), lon_idx(1), time_idx(2)) + ...
weights(1)*weights(2) * var(lat_idx(2), lon_idx(2), time_idx(2));
varargout{i} = (1-weights(3)) * spatial_interp + weights(3) * spatial_interp2;
else
varargout{i} = spatial_interp;
end
end
end
%% 最近的网格点寻找函数
function [lat_idx, lon_idx, time_idx, weights] = find_nearest_points(lat, lon, time, lat_grid, lon_grid, time_grid)
% 找到最近的网格点和插值权重
% 纬度
[~, lat_idx(1)] = min(abs(lat_grid - lat));
if lat_grid(lat_idx(1)) < lat
lat_idx(2) = min(lat_idx(1)+1, length(lat_grid));
else
lat_idx(2) = max(lat_idx(1)-1, 1);
end
weights(1) = (lat - lat_grid(lat_idx(1))) / (lat_grid(lat_idx(2)) - lat_grid(lat_idx(1)));
% 经度
[~, lon_idx(1)] = min(abs(lon_grid - lon));
if lon_grid(lon_idx(1)) < lon
lon_idx(2) = min(lon_idx(1)+1, length(lon_grid));
else
lon_idx(2) = max(lon_idx(1)-1, 1);
end
weights(2) = (lon - lon_grid(lon_idx(1))) / (lon_grid(lon_idx(2)) - lon_grid(lon_idx(1)));
% 时间
time_diff = abs(time_grid - time);
[~, time_idx(1)] = min(time_diff);
if time_grid(time_idx(1)) < time
time_idx(2) = min(time_idx(1)+1, length(time_grid));
else
time_idx(2) = max(time_idx(1)-1, 1);
end
if length(time_idx) > 1
weights(3) = seconds(time - time_grid(time_idx(1))) / seconds(time_grid(time_idx(2)) - time_grid(time_idx(1)));
end
end
%% 角度计算函数
function heading = calculate_heading(point1, point2)
% 计算两点之间的航向角 (0-360度,正北为0)
lat1 = deg2rad(point1(1));
lon1 = deg2rad(point1(2));
lat2 = deg2rad(point2(1));
lon2 = deg2rad(point2(2));
dLon = lon2 - lon1;
y = sin(dLon) * cos(lat2);
x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dLon);
heading = wrapTo360(rad2deg(atan2(y, x)));
end
%% 船舶油耗及速度计算函数
function [V_actual, fuel_consumption] = calculate_ship_performance(ship, wind_speed, rel_wind_angle, swh, cdww)
% 计算考虑风浪影响的船舶实际速度和燃油消耗
% 计算风阻力
R_wind = 0.5 * ship.rho_air * wind_speed^2 * ship.A_T * abs(cosd(rel_wind_angle));
% 计算波浪阻力 (简化模型)
R_wave = 0.5 * ship.rho_sea * (swh/2)^2 * ship.B * ship.T * (1 + 4 * (ship.T/ship.L_pp)^2);
% 计算静水阻力
R_water = 0.5 * ship.rho_sea * ship.L_pp * ship.B * cdww;
V_knots = ship.V; % 初始使用设计航速
V_ms = V_knots * 0.5144; % 节转m/s
% 总静水阻力
R_total = R_wind + R_wave + R_water;
% 考虑风浪影响的速度损失 (经验公式)
speed_loss = 0.05 * (swh/4) + 0.02 * wind_speed/10;
V_actual = max(ship.V_min, ship.V * (1 - speed_loss));
% 计算燃油消耗 (简化模型)
V_actual_ms = V_actual * 0.5144; % 节转m/s
power_required = R_total * V_actual_ms / 0.6; % 推进效率0.6
fuel_consumption = power_required * ship.sfoc / 1e6; % kg/h
end
%% 选择算子
function selected = tournament_selection(population, fitness, num_selected)
% 锦标赛选择
selected = cell(num_selected, 1);
for i = 1:num_selected
% 随机选择2个个体进行竞争
candidates = randperm(length(population), 2);
if fitness(candidates(1)) < fitness(candidates(2))
selected{i} = population{candidates(1)};
else
selected{i} = population{candidates(2)};
end
end
end
%% 交叉算子
function offspring = hybrid_crossover(parents, pc, start_port, end_port)
offspring = cell(length(parents), 1); % 确保子代数量与父代相同
i = 1;
while i <= length(parents)
if i+1 > length(parents)
% 处理奇数情况:最后一个父代直接复制
offspring{i} = parents{i};
break;
end
if rand() < pc
parent1 = parents{i};
parent2 = parents{i+1};
% 检查父代路径有效性
if size(parent1,1) < 3 || size(parent2,1) < 3
offspring{i} = parent1;
offspring{i+1} = parent2;
i = i + 2;
continue;
end
min_len = min(size(parent1,1), size(parent2,1));
if min_len <= 2
offspring{i} = parent1;
offspring{i+1} = parent2;
i = i + 2;
continue;
end
% 随机交叉点
crossover_point = randi([2, min_len-1]);
% 生成子代
child1 = [parent1(1:crossover_point,:); parent2(crossover_point+1:end,:)];
child2 = [parent2(1:crossover_point,:); parent1(crossover_point+1:end,:)];
% 确保子代有效性
child1(1,:) = start_port;
child1(end,:) = end_port;
child2(1,:) = start_port;
child2(end,:) = end_port;
offspring{i} = child1;
offspring{i+1} = child2;
else
% 不交叉时直接复制父代
offspring{i} = parents{i};
offspring{i+1} = parents{i+1};
end
i = i + 2; % 移动到下一对父代
end
end
%% 变异算子
function mutated_pop = mutation(population, pm, lat_range, lon_range)
mutated_pop = cell(size(population));
for i = 1:length(population)
if rand() < pm
route = population{i};
if size(route,1) <= 2
mutated_pop{i} = route;
continue;
end
idx = randi([2, size(route,1)-1]); % 随机选择中间点
new_lat = lat_range(1) + (lat_range(2)-lat_range(1))*rand();
new_lon = lon_range(1) + (lon_range(2)-lon_range(1))*rand();
route(idx,:) = [new_lat, new_lon];
mutated_pop{i} = route;
else
mutated_pop{i} = population{i};
end
end
end
%% 模拟退火
function new_population = simulated_annealing_selection(population, fitness, pop_size, temperature)
% 模拟退火选择
[~, sorted_idx] = sort(fitness);
new_population = cell(pop_size, 1);
% 保留最优个体
new_population{1} = population{sorted_idx(1)};
% 对其他个体进行退火选择
for i = 2:pop_size
% 随机选择两个个体
candidates = randperm(length(population), 2);
delta = fitness(candidates(2)) - fitness(candidates(1));
% 接受概率
if delta < 0
new_population{i} = population{candidates(2)};
else
p = exp(-delta / temperature);
if rand() < p
new_population{i} = population{candidates(2)};
else
new_population{i} = population{candidates(1)};
end
end
end
end
%% 绘图函数
function plot_route(route, start_port, end_port, weather_data)
% 绘制优化后的航线
figure;
% 绘制地图背景
worldmap([-40 60], [100 170]);
load coastlines;
plotm(coastlat, coastlon);
% 绘制航线
plotm(route(:,1), route(:,2), 'r-o', 'LineWidth', 2);
plotm(start_port(1), start_port(2), 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g');
plotm(end_port(1), end_port(2), 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b');
% 添加标题和图例
title('Optimized Ship Route from Dalian to Gladstone');
legend('Route', 'Dalian', 'Gladstone', 'Location', 'best');
% 添加气象数据示例 (可选)
hold on;
[LON, LAT] = meshgrid(weather_data.op_lon, weather_data.op_lat);
quiverm(LAT(1:5:end,1:5:end), LON(1:5:end,1:5:end), ...
weather_data.u10(1:5:end,1:5:end,1), weather_data.v10(1:5:end,1:5:end,1), 'k');
hold off;
end
%% 距离及角度计算函数
function [km, arc] = lldistkm(latlon1, latlon2)
% 计算两点之间的大圆距离
lat1 = latlon1(1);
lon1 = latlon1(2);
lat2 = latlon2(1);
lon2 = latlon2(2);
R = 6371; % 地球半径(km)
lat1 = deg2rad(lat1);
lon1 = deg2rad(lon1);
lat2 = deg2rad(lat2);
lon2 = deg2rad(lon2);
delta_lat = lat2 - lat1;
delta_lon = lon2 - lon1;
a = sin(delta_lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta_lon/2)^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
km = R * c;
arc = rad2deg(c);
end
那你把这个代码修改一下发给我
最新发布