该代码在运行时发生以下错误,第101行:aircraft.vk 不是标量或单元素数组,请检查数据!
%% 碳交易机制下基于多目标的机队优化配置 - MATLAB遗传算法实现(优化版)
clear variables;
close all;
clc;
tic; % 开始计时
fprintf('程序开始运行...\n\n');
%% 1. 数据初始化
% 机型数据(4种机型)
aircraft = struct();
aircraft.k = [1, 2, 3, 4]; % 机型编号
aircraft.num = [6, 6, 6, 3]; % 可用数量
aircraft.weight = [42400, 129000, 80127, 78250]; % 空重(kg)
aircraft.seats = [156, 324, 226, 132]; % 座位数
aircraft.range = [4310, 12000, 7300, 6000]; % 航程(km)
aircraft.vk = [800,860,925,900] ; %巡航速度(km/h)
% 航线数据(43个航班)
flights = struct();
flights.id = 1:43; % 航班编号
% 旅客需求(示例数据)
flights.demand = [132,134,108,156,200,180,150,140,130,120,...
110,100,90,80,170,160,150,140,130,120,...
110,100,90,80,170,160,150,140,130,120,...
110,100,90,80,170,160,150,140,130,120,110,100,90];
% 飞行距离(km,示例数据)
flights.distance = [2349,3425,1234,2567,3456,2345,...
3456,1234,2345,3456,4567,5678,...
2345,3456,4567,1234,2345,3456,...
4567,5678,1234,2345,3456,4567,...
5678,1234,2345,3456,4567,5678,...
1234,2345,3456,4567,5678,1234,...
2345,3456,4567,5678,1234,2345,2756];
% 平均票价(元,示例数据)
flights.price = [1260,2020,1500,1800,2200,1900,...
1700,1600,1500,1400,1300,1200,...
1100,1000,2100,2000,1900,1800,...
1700,1600,1500,1400,1300,1200,...
1100,1000,2100,2000,1900,1800,...
1700,1600,1500,1400,1300,1200,...
1100,1000,2100,2000,1900,1800,1300];
% 碳交易参数
mu = 0.7; % 燃油排放系数(kg CO₂/kg燃油)
free_quota_ratio = 0.8; % 免费碳配额比例
carbon_price = 50; % 碳交易价格(元/吨CO₂)
% 遗传算法参数(优化后)
pop_size = 200; % 增大种群大小
max_gen = 500; % 增加迭代次数
pm_max = 0.5; % 提高最大变异概率
pm_min = 0.05; % 提高最小变异概率
pm_mid = 0.15; % 中间变异概率
alpha = 0.5; % 降低alpha值
beta = 0.5; % 提高beta值
elite_ratio = 0.1; % 精英保留比例
fprintf('数据初始化完成!\n');
%% 2. 预处理数据
n_flights = length(flights.id);
n_aircraft = length(aircraft.k);
% 预计算满足条件的机型矩阵(性能优化关键)
fprintf('预计算机型可用性矩阵...');
valid_aircraft_mask = false(n_flights, n_aircraft);
for i = 1:n_flights
for k = 1:n_aircraft
% 检查机型是否满足航程和座位要求
valid_aircraft_mask(i, k) = (flights.distance(i) <= aircraft.range(k)) && ...
(aircraft.seats(k) >= flights.demand(i));
end
end
fprintf('完成!\n');
% 计算初始碳排放总量(使用最大机型执飞所有航班)
fprintf('计算初始碳排放...');
initial_emission = 0;
max_seats_idx = find(aircraft.seats == max(aircraft.seats), 1);
% 使用方案二:更严谨地处理 ratio 和 vk 的标量转换
if isfield(aircraft, 'ratio')
ratio_tmp = aircraft.ratio;
if numel(ratio_tmp) == 1
ratio = ratio_tmp;
else
error('aircraft.ratio 不是标量或单元素数组,请检查数据!');
end
else
ratio = 0.8; % 合理默认值
end
if isfield(aircraft, 'vk')
vk_tmp = aircraft.vk;
if numel(vk_tmp) == 1
vk = vk_tmp;
else
error('aircraft.vk 不是标量或单元素数组,请检查数据!');
end
else
vk = 800; % 合理默认值,单位 km/h,可根据实际情况调整
end
for i = 1:n_flights
disi = flights.distance(i);
% 这里使用新公式计算飞机重量:W_i = EW_i + 50*SN_i + 100*Q ,Q就是flights.demand(i)
Wk = aircraft.weight(max_seats_idx) + 50 * aircraft.seats(max_seats_idx) + 100 * flights.demand(i);
initial_emission = initial_emission + mu * (1 - exp(-disi * ratio / (10 * vk))) * Wk;
end
initial_emission = initial_emission / 1000; % 转换为吨
free_quota = initial_emission * free_quota_ratio;
fprintf('完成! 初始排放: %.2f吨, 免费配额: %.2f吨\n', initial_emission, free_quota);
%% 3. 遗传算法主程序(优化后)
fprintf('\n开始遗传算法优化...\n种群大小: %d, 最大迭代次数: %d\n', pop_size, max_gen);
% 初始化种群
population = initialize_population(pop_size, n_flights, n_aircraft, valid_aircraft_mask, flights, aircraft);
% 记录每一代的最佳适应度和平均适应度
best_fitness = zeros(max_gen, 1);
avg_fitness = zeros(max_gen, 1);
% 记录种群多样性
diversity_history = zeros(max_gen, 1);
% 记录最佳个体历史
best_individual_history = zeros(max_gen, n_flights * n_aircraft);
% 预分配内存(性能优化)
revenue_arr = zeros(pop_size, 1);
emission_arr = zeros(pop_size, 1);
for gen = 1:max_gen
% 计算适应度(使用新的适应度函数)
[fitness, revenue_arr, emission_arr] = calculate_fitness(population, n_flights, n_aircraft, ...
flights, aircraft, mu, free_quota, carbon_price, valid_aircraft_mask);
% 记录种群多样性
diversity_history(gen) = calculate_diversity(population);
% 记录最佳个体
[max_fit, max_idx] = max(fitness);
best_fitness(gen) = max_fit;
best_individual = population(max_idx, :);
best_individual_history(gen, :) = best_individual;
avg_fitness(gen) = mean(fitness);
% 选择操作(使用锦标赛选择)
new_population = selection(population, fitness, pop_size, elite_ratio);
% 交叉操作
new_population = crossover(new_population, pop_size, n_flights, n_aircraft, valid_aircraft_mask);
% 计算自适应变异概率
pm = adaptive_mutation_probability(fitness, pm_max, pm_min, pm_mid, alpha, beta);
% 变异操作(增强变异强度)
new_population = mutation(new_population, pm, n_flights, n_aircraft, valid_aircraft_mask, flights, aircraft);
% 更新种群
population = new_population;
% 每10代显示一次迭代进度
if mod(gen, 10) == 0 || gen == 1
fprintf('迭代 %d/%d: 最佳适应度 = %.4f, 平均适应度 = %.4f, 多样性 = %.4f\n', ...
gen, max_gen, max_fit, avg_fitness(gen), diversity_history(gen));
end
end
% 获取最终最佳个体
[~, best_gen] = max(best_fitness);
best_individual = best_individual_history(best_gen, :);
fprintf('\n遗传算法优化完成! 耗时 %.2f 秒\n', toc);
%% 4. 结果分析
fprintf('\n===== 结果分析 =====\n');
% 计算最终方案
[~, final_revenue, final_emission] = calculate_fitness(best_individual, n_flights, n_aircraft, ...
flights, aircraft, mu, free_quota, carbon_price, valid_aircraft_mask);
% 方案1: 最大机型执飞所有航班
base_individual = zeros(1, n_flights * n_aircraft);
for i = 1:n_flights
[~, max_seats_idx] = max(aircraft.seats);
pos = (i-1)*n_aircraft + max_seats_idx;
base_individual(pos) = 1;
end
[~, base_revenue, base_emission] = calculate_fitness(base_individual, n_flights, n_aircraft, ...
flights, aircraft, mu, free_quota, carbon_price, valid_aircraft_mask);
% 方案2: 成本最小化方案
cost_individual = zeros(1, n_flights * n_aircraft);
for i = 1:n_flights
min_cost = inf;
best_aircraft = 0;
% 获取当前航班可用的机型(使用预计算的矩阵)
valid_k = find(valid_aircraft_mask(i, :));
for k_idx = 1:length(valid_k)
k = valid_k(k_idx);
cost = calculate_operation_cost(flights.distance(i), aircraft.weight(k), ...
aircraft.seats(k), carbon_price);
if cost < min_cost
min_cost = cost;
best_aircraft = k;
end
end
if best_aircraft > 0
pos = (i-1)*n_aircraft + best_aircraft;
cost_individual(pos) = 1;
end
end
[~, cost_revenue, cost_emission] = calculate_fitness(cost_individual, n_flights, n_aircraft, ...
flights, aircraft, mu, free_quota, carbon_price, valid_aircraft_mask);
% 输出结果对比
fprintf('\n===== 方案对比 =====\n');
fprintf('方案1 (最大机型): 收益 = %.2f元, 排放 = %.2f吨\n', base_revenue, base_emission);
fprintf('方案2 (成本最小): 收益 = %.2f元, 排放 = %.2f吨\n', cost_revenue, cost_emission);
fprintf('方案3 (遗传算法): 收益 = %.2f元, 排放 = %.2f吨\n', final_revenue, final_emission);
% 计算优化效果
revenue_decrease1 = (base_revenue - final_revenue) / base_revenue * 100;
emission_decrease1 = (base_emission - final_emission) / base_emission * 100;
revenue_improve2 = (final_revenue - cost_revenue) / cost_revenue * 100;
emission_improve2 = (cost_emission - final_emission) / cost_emission * 100;
fprintf('\n===== 优化效果 =====\n');
fprintf('方案3 vs 方案1: 收益变化 %.2f%%, 排放减少 %.2f%%\n', revenue_decrease1, emission_decrease1);
fprintf('方案3 vs 方案2: 收益提升 %.2f%%, 排放减少 %.2f%%\n', revenue_improve2, emission_improve2);
% 绘制适应度进化曲线
figure('Name', '适应度进化曲线', 'NumberTitle', 'off');
subplot(2,1,1);
plot(1:max_gen, best_fitness, 'b-', 'LineWidth', 1.5);
hold on;
plot(1:max_gen, avg_fitness, 'r--', 'LineWidth', 1.5);
title('适应度进化曲线', 'FontSize', 14);
xlabel('迭代次数', 'FontSize', 12);
ylabel('适应度', 'FontSize', 12);
legend('最佳适应度', '平均适应度', 'Location', 'southeast');
grid on;
set(gca, 'FontSize', 10);
% 绘制多样性曲线
subplot(2,1,2);
plot(1:max_gen, diversity_history, 'g-', 'LineWidth', 1.5);
title('种群多样性变化', 'FontSize', 14);
xlabel('迭代次数', 'FontSize', 12);
ylabel('多样性指数', 'FontSize', 12);
grid on;
set(gca, 'FontSize', 10);
% 绘制收益-排放对比图
figure('Name', '收益-排放对比', 'NumberTitle', 'off');
scatter([base_emission, cost_emission, final_emission], ...
[base_revenue, cost_revenue, final_revenue], 120, ...
[1 0 0; 0 0.8 0; 0 0 1], 'filled');
text(base_emission, base_revenue, ' 方案1 (最大机型)', 'FontSize', 10, 'VerticalAlignment', 'bottom');
text(cost_emission, cost_revenue, ' 方案2 (成本最小)', 'FontSize', 10, 'VerticalAlignment', 'top');
text(final_emission, final_revenue, ' 方案3 (遗传算法)', 'FontSize', 10, 'VerticalAlignment', 'bottom');
title('收益-排放对比', 'FontSize', 14);
xlabel('碳排放量 (吨)', 'FontSize', 12);
ylabel('运营收益 (元)', 'FontSize', 12);
grid on;
set(gca, 'FontSize', 10);
%% 5. 新增可视化部分(优化后的热力图等)
% ===== 1. 机型使用频率对比图 =====
figure('Name', '机型使用频率对比', 'NumberTitle', 'off', 'Position', [100, 100, 800, 400]);
% 解码三种方案的机型分配
base_assigned = decode_assignment(base_individual, n_flights, n_aircraft);
cost_assigned = decode_assignment(cost_individual, n_flights, n_aircraft);
final_assigned = decode_assignment(best_individual, n_flights, n_aircraft);
% 统计各机型使用次数
base_counts = sum(base_assigned);
cost_counts = sum(cost_assigned);
final_counts = sum(final_assigned);
% 绘制对比条形图
bar_data = [base_counts; cost_counts; final_counts]';
bar_handle = bar(bar_data, 'grouped');
title('机型使用频率对比', 'FontSize', 14);
xlabel('机型', 'FontSize', 12);
ylabel('使用次数', 'FontSize', 12);
legend({'方案1 (最大机型)', '方案2 (成本最小)', '方案3 (遗传算法)'}, 'Location', 'best', 'FontSize', 10);
grid on;
% 添加数值标签
for i = 1:length(bar_handle)
x = bar_handle(i).XEndPoints;
y = bar_handle(i).YEndPoints;
text(x, y, string(y), 'HorizontalAlignment','center',...
'VerticalAlignment','bottom', 'FontSize', 9, 'FontWeight','bold');
end
set(gca, 'XTickLabel', {'机型1','机型2','机型3','机型4'}, 'FontSize', 10);
% ===== 2. 优化效果雷达图 =====
figure('Name', '多目标优化效果对比', 'NumberTitle', 'off', 'Position', [200, 200, 700, 600]);
% 准备数据(归一化)
metrics = [
base_revenue, base_emission;
cost_revenue, cost_emission;
final_revenue, final_emission
];
% 归一化处理(越大越好)
norm_metrics = metrics ./ max(metrics);
% 雷达图参数
angles = linspace(0, 2*pi, 3);
angles = angles(1:end-1); % 只需要两个角度(收益和排放)
rlims = [0 1.2]; % 设置雷达图半径范围
% 创建极坐标图
polaraxes;
hold on;
% 绘制雷达图网格
max_val = max(max(norm_metrics));
grid_vals = linspace(0, max_val, 5);
for i = 1:length(grid_vals)
polarplot(angles, grid_vals(i)*ones(1, length(angles)), 'k:', 'LineWidth', 0.5);
end
% 绘制各方案雷达图
colors = ['r', 'g', 'b'];
line_styles = {'-', '--', '-.'};
for i = 1:3
data = [norm_metrics(i, :), norm_metrics(i, 1)]; % 闭合数据
polarplot([angles, angles(1)], data, ...
'Color', colors(i), 'LineStyle', line_styles{i}, ...
'LineWidth', 2, 'Marker', 'o', 'MarkerSize', 6, ...
'MarkerFaceColor', colors(i));
end
% 设置角度标签
thetaticks(rad2deg(angles));
thetaticklabels({'经济收益', '环境效益'});
rticks(grid_vals);
rticklabels(arrayfun(@(x) sprintf('%.1f', x), grid_vals, 'UniformOutput', false));
% 添加图例和标题
legend({'方案1 (最大机型)', '方案2 (成本最小)', '方案3 (遗传算法)'}, ...
'Location', 'southoutside', 'Orientation', 'horizontal', 'FontSize', 10);
title('多目标优化效果雷达图', 'FontSize', 14);
set(gca, 'FontSize', 10);
hold off;
% ===== 3. 优化后的热力图(核心部分) =====
figure('Name', '优化机型分配热力图', 'NumberTitle', 'off', 'Position', [300, 300, 900, 500]);
% 解码遗传算法优化后的分配方案
assignment_matrix = decode_assignment(best_individual, n_flights, n_aircraft);
% 绘制热力图
imagesc(assignment_matrix);
colormap(jet(4)); % 使用4种颜色对应4种机型
% 定制化colorbar
c = colorbar('Ticks', [0.375, 1.125, 1.875, 2.625], ...
'TickLabels', {'机型1', '机型2', '机型3', '机型4'});
c.Label.String = '分配机型';
c.Label.FontSize = 12;
% 添加注释(需求满足率)
for i = 1:n_flights
for k = 1:n_aircraft
if assignment_matrix(i, k) == 1
% 计算需求满足率
satisfy_rate = min(1, flights.demand(i)/aircraft.seats(k));
% 根据满足率选择文本颜色
if satisfy_rate < 0.7
text_color = [1, 0.8, 0.8]; % 浅红色
elseif satisfy_rate > 0.95
text_color = [0.8, 1, 0.8]; % 浅绿色
else
text_color = [1, 1, 1]; % 白色
end
% 添加文本标注
text(k, i, sprintf('%.0f%%', satisfy_rate*100), ...
'HorizontalAlignment', 'center', ...
'Color', text_color, ...
'FontSize', 9, ...
'FontWeight', 'bold');
end
end
end
% 添加标题和标签
title_str = {
'优化后航班-机型分配热力图'
sprintf('总收益: %.0f元 | 碳排放: %.1f吨 | 优化效果: +%.1f%%收益, -%.1f%%排放', ...
final_revenue, final_emission, revenue_improve2, emission_improve2)
};
title(title_str, 'FontSize', 14);
xlabel('机型', 'FontSize', 12);
ylabel('航班编号', 'FontSize', 12);
% 设置坐标轴
set(gca, 'FontSize', 10, 'YTick', 1:n_flights, 'YTickLabel', flights.id);
% 添加航班分组参考线
for i = 10:10:n_flights
line(xlim, [i+0.5 i+0.5], 'Color', 'k', 'LineWidth', 1.5, 'LineStyle', '--');
end
% 高亮显示需求满足率低的航班
for i = 1:n_flights
k = find(assignment_matrix(i, :));
if ~isempty(k)
satisfy_rate = min(1, flights.demand(i)/aircraft.seats(k));
if satisfy_rate < 0.7
rectangle('Position', [k-0.5, i-0.5, 1, 1], ...
'EdgeColor', 'r', 'LineWidth', 1.5, 'LineStyle', '-');
end
end
end
fprintf('\n所有可视化图表生成完成!\n');
%% 5. 函数定义(优化后)
% ===== 初始化种群函数 =====
function population = initialize_population(pop_size, n_flights, n_aircraft, valid_aircraft_mask, flights, aircraft)
population = zeros(pop_size, n_flights * n_aircraft);
for p = 1:pop_size
for i = 1:n_flights
% 获取当前航班可用的机型
valid_aircraft = find(valid_aircraft_mask(i, :));
if ~isempty(valid_aircraft)
% 80%概率选择座位数最接近需求的机型,20%概率随机选择
if rand() < 0.8
seat_diffs = abs(aircraft.seats(valid_aircraft) - flights.demand(i));
[~, min_idx] = min(seat_diffs);
best_aircraft = valid_aircraft(min_idx);
else
best_aircraft = valid_aircraft(randi(length(valid_aircraft)));
end
pos = (i-1)*n_aircraft + best_aircraft;
population(p, pos) = 1;
else
% 如无满足需求机型,选择满足航程的随机机型
valid_range = find(flights.distance(i) <= aircraft.range);
if ~isempty(valid_range)
best_aircraft = valid_range(randi(length(valid_range)));
pos = (i-1)*n_aircraft + best_aircraft;
population(p, pos) = 1;
end
end
end
end
fprintf('种群初始化完成! 种群大小: %d\n', pop_size);
end
% ===== 计算运营成本函数 =====
function cost = calculate_operation_cost(distance, weight, seats, carbon_price)
% 燃油消耗模型 (kg/km)
fuel_consumption_rate = 0.05; % 燃油消耗率 (kg/km/ton)
total_weight = weight + seats * 75; % 总重量 (kg)
% 总燃油消耗 (kg)
fuel_consumption = fuel_consumption_rate * distance * total_weight / 1000;
fuel_cost = fuel_consumption * 5; % 燃油价格5元/kg
% 碳排放成本(基于燃油消耗)
carbon_emission = fuel_consumption * 3.15 / 1000; % 吨CO2(3.15吨CO2/吨燃油)
carbon_cost = carbon_emission * carbon_price; % 碳交易成本
% 其他运营成本(机组、维护等)
other_costs = 0.2 * distance * seats; % 其他运营成本
cost = fuel_cost + carbon_cost + other_costs;
end
% ===== 计算适应度函数(优化后) =====
function [fitness, revenue_arr, emission_arr] = calculate_fitness(population, n_flights, n_aircraft, ...
flights, aircraft, mu, free_quota, carbon_price, valid_aircraft_mask)
pop_size = size(population, 1);
fitness = zeros(pop_size, 1);
revenue_arr = zeros(pop_size, 1);
emission_arr = zeros(pop_size, 1);
% 预计算常量
ratio = 0.8;
vk = 800;
% 预先计算整个种群的净收益和排放
all_net_revenue = zeros(pop_size, 1);
all_emission = zeros(pop_size, 1);
for p = 1:pop_size
total_revenue = 0;
total_emission = 0;
total_cost = 0;
for i = 1:n_flights
for k = 1:n_aircraft
pos = (i-1)*n_aircraft + k;
if population(p, pos) == 1
% 检查分配的机型是否有效
if ~valid_aircraft_mask(i, k)
continue; % 跳过无效分配
end
pi = flights.price(i);
Di = flights.demand(i);
Capk = aircraft.seats(k);
% 收入计算
revenue_i = pi * min(Di, Capk);
total_revenue = total_revenue + revenue_i;
% 碳排放计算
disi = flights.distance(i);
Wk = aircraft.weight(k) + Capk * 75;
emission_i = mu * (1 - exp(-disi * ratio / (10 * vk))) * Wk;
total_emission = total_emission + emission_i;
% 运营成本
cost = calculate_operation_cost(disi, aircraft.weight(k), Capk, carbon_price);
total_cost = total_cost + cost;
end
end
end
% 碳交易成本
if total_emission > free_quota
carbon_cost = (total_emission - free_quota) * carbon_price / 1000; % 转换为元
else
carbon_cost = 0;
end
total_cost = total_cost + carbon_cost;
% 净利润计算
net_revenue = total_revenue - total_cost;
% 保存当前个体的结果
all_net_revenue(p) = net_revenue;
all_emission(p) = total_emission / 1000; % 转换为吨
end
% 计算整个种群的最大最小值
max_revenue = max(all_net_revenue);
min_revenue = min(all_net_revenue);
max_emission = max(all_emission);
min_emission = min(all_emission);
% 归一化处理(避免除零)
revenue_range = max_revenue - min_revenue + eps;
emission_range = max_emission - min_emission + eps;
% 权重系数
lambda1 = 0.7; % 经济目标权重
lambda2 = 0.3; % 环境目标权重
for p = 1:pop_size
% 归一化目标值
norm_revenue = (all_net_revenue(p) - min_revenue) / revenue_range;
norm_emission = (max_emission - all_emission(p)) / emission_range; % 排放越小越好
% 加权适应度函数
fitness(p) = lambda1 * norm_revenue + lambda2 * norm_emission;
revenue_arr(p) = all_net_revenue(p);
emission_arr(p) = all_emission(p);
end
end
% ===== 选择操作(锦标赛选择) =====
function new_population = selection(population, fitness, pop_size, elite_ratio)
% 精英选择
elite_num = floor(pop_size * elite_ratio);
[~, elite_idx] = sort(fitness, 'descend');
elite_idx = elite_idx(1:elite_num);
new_population = population(elite_idx, :);
% 锦标赛选择
for i = (elite_num+1):pop_size
% 随机选择3个个体
candidates = randperm(size(population,1), 3);
[~, best_idx] = max(fitness(candidates));
new_population(i, :) = population(candidates(best_idx), :);
end
end
% ===== 交叉操作(单点交叉) =====
function new_population = crossover(new_population, pop_size, n_flights, n_aircraft, valid_aircraft_mask)
cross_prob = 0.8;
for i = 1:2:(pop_size-1)
if rand() < cross_prob
% 随机选择航班交叉点(在航班边界上)
cross_point = randi(n_flights-1) * n_aircraft;
% 执行交叉
temp = new_population(i, cross_point+1:end);
new_population(i, cross_point+1:end) = new_population(i+1, cross_point+1:end);
new_population(i+1, cross_point+1:end) = temp;
% 修复交叉可能引起的约束问题
new_population(i, :) = repair_constraints(new_population(i, :), n_flights, n_aircraft, valid_aircraft_mask);
new_population(i+1, :) = repair_constraints(new_population(i+1, :), n_flights, n_aircraft, valid_aircraft_mask);
end
end
end
% ===== 自适应变异概率函数(移除delta参数) =====
function pm = adaptive_mutation_probability(fitness, pm_max, pm_min, pm_mid, alpha, beta)
f_min = min(fitness);
f_max = max(fitness);
% 归一化适应度
f_norm = (fitness - f_min) / (f_max - f_min + eps);
% 计算自适应概率(基于指数函数)
pm = pm_mid + (pm_max - pm_mid) * exp(-alpha * f_norm) - ...
(pm_mid - pm_min) * exp(-beta * (1 - f_norm));
% 限制在合理范围内
pm = max(pm_min, min(pm_max, pm));
end
% ===== 变异操作(增强变异强度) =====
function new_population = mutation(new_population, pm, n_flights, n_aircraft, valid_aircraft_mask, flights, aircraft)
for i = 1:size(new_population, 1)
if rand() < pm(i)
% 变异多个航班(最多30%的航班)
num_mutations = randi([1, ceil(n_flights*0.3)]);
flights_to_mutate = randperm(n_flights, num_mutations);
for j = 1:num_mutations
flight_idx = flights_to_mutate(j);
start_pos = (flight_idx-1)*n_aircraft + 1;
end_pos = flight_idx*n_aircraft;
% 重置当前航班分配
new_population(i, start_pos:end_pos) = 0;
% 使用预计算的机型掩码选择新机型
valid_aircraft = find(valid_aircraft_mask(flight_idx, :));
if ~isempty(valid_aircraft)
% 80%概率选择最优机型,20%概率随机选择
if rand() < 0.8
% 选择座位数最接近需求的机型
seat_diffs = abs(aircraft.seats(valid_aircraft) - flights.demand(flight_idx));
[~, min_idx] = min(seat_diffs);
selected = valid_aircraft(min_idx);
else
selected = valid_aircraft(randi(length(valid_aircraft)));
end
new_population(i, start_pos + selected - 1) = 1;
else
% 无有效机型时选择满足航程的随机机型
valid_range = find(flights.distance(flight_idx) <= aircraft.range);
if ~isempty(valid_range)
selected = valid_range(randi(length(valid_range)));
new_population(i, start_pos + selected - 1) = 1;
end
end
end
end
end
end
% ===== 约束修复函数 =====
function chromosome = repair_constraints(chromosome, n_flights, n_aircraft, valid_aircraft_mask)
for i = 1:n_flights
start_idx = (i-1)*n_aircraft + 1;
end_idx = i*n_aircraft;
gene = chromosome(start_idx:end_idx);
% 检查是否已分配且有效的机型
assigned = find(gene);
if ~isempty(assigned)
k = assigned(1);
if ~valid_aircraft_mask(i, k)
gene(k) = 0; % 移除无效分配
end
end
% 确保每个航班有且仅有一个有效分配
if sum(gene) ~= 1
valid_aircraft = find(valid_aircraft_mask(i, :));
gene = zeros(1, n_aircraft);
if ~isempty(valid_aircraft)
% 随机选择有效机型
selected = valid_aircraft(randi(length(valid_aircraft)));
gene(selected) = 1;
else
% 无有效机型时选择第一个满足航程的机型
% (实际应用中应避免此情况)
valid_range = find(flights.distance(i) <= aircraft.range);
if ~isempty(valid_range)
selected = valid_range(1);
gene(selected) = 1;
end
end
chromosome(start_idx:end_idx) = gene;
end
end
end
% ===== 解码分配方案 =====
function assignment_matrix = decode_assignment(individual, n_flights, n_aircraft)
assignment_matrix = zeros(n_flights, n_aircraft);
for i = 1:n_flights
start_idx = (i-1)*n_aircraft + 1;
end_idx = i*n_aircraft;
gene = individual(start_idx:end_idx);
k = find(gene);
if ~isempty(k)
assignment_matrix(i, k) = 1;
end
end
end
% ===== 计算种群多样性(新增) =====
function div = calculate_diversity(pop)
% 计算种群汉明距离
n = size(pop, 1);
total_dist = 0;
% 如果种群太小,直接返回0
if n < 2
div = 0;
return;
end
for i = 1:n-1
for j = i+1:n
% 计算两个个体之间的差异
dist = sum(pop(i,:) ~= pop(j,:));
total_dist = total_dist + dist;
end
end
% 计算平均距离
pair_count = n*(n-1)/2;
div = total_dist / pair_count;
end