以下代码在运行时出错:无法解析名称 'flights.distance',出错位置在第502行、427行和108行
%% 碳交易机制下基于多目标的机队优化配置 - 最终修正版
clear variables;
close all;
clc;
tic; % 开始计时
fprintf('程序开始运行...\n\n');
%% 1. 数据初始化
% 机型数据
aircraft = struct();
aircraft.k = [1, 2, 3]; % 机型编号
aircraft.num = [36, 29, 10]; % 可用数量
aircraft.weight = [21500, 42400, 25000]; % 空重(kg)
aircraft.seats = [76, 140, 90]; % 座位数
aircraft.range = [2800, 6200, 3700]; % 最大航程(km) - 修正此处
aircraft.vk = [860, 956, 956]; % 巡航速度(km/h)
aircraft.ratio = [0.071, 0.073, 0.155]; % 燃油效率参数
% 航线数据
flights = struct();
flights.id = 1:28; % 航班编号 - 修正此处
% 航班运营数据(完整客座率数据)
flights.distance = [533, 1241, 1633, 494, 862, 1772, 521, 1709, 663, 1417, 813, 795, 1444, 1357, 2272, 835, 563, 1159, 1392, 211, 1123, 2697, 867, 1137, 1048, 1286, 388, 1363];
flights.price = [730, 1180, 1415, 520, 650, 1710, 695, 998, 460, 625, 795, 795, 1045, 827, 1357, 645, 790, 635, 735, 560, 1360, 1195, 905, 685, 855, 740, 1175, 1180];
flights.plf = [0.835, 0.846, 0.835, 0.836, 0.858, 0.835, 0.871, 0.84, 0.9, 0.93, 0.835, 0.58, 0.85, 0.835, 0.81, 0.835, 0.835, 0.835, 0.835, 0.835, 0.88, 0.87, 0.86, 0.81, 0.9,0.835, 0.93, 0.835];
% 碳交易参数
mu = 0.7; % 燃油排放系数(kg CO₂/kg燃油)
free_quota_ratio = 0.8; % 免费碳配额比例
carbon_price = 67.9; % 碳交易价格(元/吨CO₂)
spill_cost = 500; % 旅客溢出成本(元/人)
min_seats = 100; % 最低座位要求
% 遗传算法参数
pop_size = 200; % 种群大小
max_gen = 500; % 最大迭代次数
pm_max = 0.5; % 最大变异概率
pm_min = 0.05; % 最小变异概率
pm_mid = 0.15; % 中间变异概率
alpha = 0.5; % 自适应参数
beta = 0.5; % 自适应参数
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
seating_capacity = aircraft.seats(k) * flights.plf(i);
valid_aircraft_mask(i, k) = (flights.distance(i) <= aircraft.range(k)) && ...
(seating_capacity >= min_seats);
end
end
fprintf('完成!\n');
% 计算初始碳排放
fprintf('计算初始碳排放...');
initial_emission = 0;
[~, max_seats_idx] = max(aircraft.seats);
for i = 1:n_flights
disi = flights.distance(i);
vk_max = aircraft.vk(max_seats_idx);
ratio_max = aircraft.ratio(max_seats_idx);
% 重量计算(使用实际客座率)
passengers_weight = 100 * (aircraft.seats(max_seats_idx) * flights.plf(i));
Wk = aircraft.weight(max_seats_idx) + 50 * aircraft.seats(max_seats_idx) + passengers_weight; % 修正此处
emission_i = mu * (1 - exp(-disi * ratio_max / (10 * vk_max))) * Wk;
initial_emission = initial_emission + emission_i;
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);
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, spill_cost);
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;
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, spill_cost);
% 方案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, spill_cost);
% 方案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), flights.plf(i), carbon_price, spill_cost, aircraft.seats(k)*flights.plf(i)); % 修正此处
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, spill_cost);
% 输出结果对比
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);
%% 5. 可视化分析
% 适应度进化曲线
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;
% 多样性曲线
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;
% 收益-排放对比图
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;
% 机型使用频率对比图
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'}, 'FontSize', 10);
% 优化机型分配热力图
figure('Name', '优化机型分配热力图', 'NumberTitle', 'off', 'Position', [300, 300, 900, 500]);
assignment_matrix = decode_assignment(best_individual, n_flights, n_aircraft);
imagesc(assignment_matrix);
colormap(jet(3)); % 3种颜色对应3种机型
% 定制化colorbar
c = colorbar('Ticks', [0.5, 1.5, 2.5], ...
'TickLabels', {'机型1', '机型2', '机型3'});
c.Label.String = '分配机型';
c.Label.FontSize = 12;
% 添加客座率标注
for i = 1:n_flights
for k = 1:n_aircraft
if assignment_matrix(i, k) == 1
% 显示客座率百分比
text(k, i, sprintf('%.0f%%', flights.plf(i)*100), ...
'HorizontalAlignment', 'center', ...
'Color', [1, 1, 1], ...
'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);
fprintf('\n所有可视化图表生成完成!\n');
%% 6. 函数定义
% ===== 初始化种群函数 =====
function population = initialize_population(pop_size, n_flights, n_aircraft, valid_aircraft_mask, flights, aircraft, ~)
population = zeros(pop_size, n_flights * n_aircraft);
avg_demand = mean(aircraft.seats) * mean(flights.plf); % 平均需求参考值
for p = 1:pop_size
for i = 1:n_flights
valid_aircraft = find(valid_aircraft_mask(i, :));
if ~isempty(valid_aircraft)
if rand() < 0.8
seating_capacities = aircraft.seats(valid_aircraft) .* flights.plf(i);
diffs = abs(seating_capacities - avg_demand);
[~, min_idx] = min(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, plf, carbon_price, spill_cost, demand)
fuel_consumption_rate = 0.05;
actual_passengers = min(seats * plf, demand);
spilled_passengers = max(0, demand - actual_passengers);
total_weight = weight + 50 * seats + 100 * actual_passengers;
fuel_consumption = fuel_consumption_rate * distance * total_weight / 1000;
fuel_cost = fuel_consumption * 6.274;
carbon_emission = fuel_consumption * 3.15 / 1000;
carbon_cost = carbon_emission * carbon_price;
other_costs = 0.2 * distance * seats;
spill_penalty = spilled_passengers * spill_cost;
cost = fuel_cost + carbon_cost + other_costs + spill_penalty;
end
% ===== 计算适应度函数 =====
function [fitness, revenue_arr, emission_arr] = calculate_fitness(population, n_flights, n_aircraft, ...
flights, aircraft, mu, free_quota, carbon_price, valid_aircraft_mask, spill_cost)
pop_size = size(population, 1);
fitness = zeros(pop_size, 1);
revenue_arr = zeros(pop_size, 1);
emission_arr = zeros(pop_size, 1);
all_net_revenue = zeros(pop_size, 1);
all_emission = zeros(pop_size, 1);
avg_demand = mean(aircraft.seats) * mean(flights.plf);
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 && valid_aircraft_mask(i, k)
pi = flights.price(i);
Capk = aircraft.seats(k);
plf = flights.plf(i);
disi = flights.distance(i);
actual_passengers = Capk * plf;
revenue_i = pi * actual_passengers;
total_revenue = total_revenue + revenue_i;
Wk = aircraft.weight(k) + 50 * Capk + 100 * actual_passengers;
vk_current = aircraft.vk(k);
ratio_current = aircraft.ratio(k);
emission_i = mu * (1 - exp(-disi * ratio_current / (10 * vk_current))) * Wk;
total_emission = total_emission + emission_i;
cost = calculate_operation_cost(disi, aircraft.weight(k), Capk, plf, carbon_price, spill_cost, avg_demand);
total_cost = total_cost + cost;
end
end
end
total_emission_ton = total_emission / 1000;
if total_emission_ton > free_quota
carbon_cost = (total_emission_ton - free_quota) * carbon_price;
else
carbon_cost = 0;
end
total_cost = total_cost + carbon_cost;
all_net_revenue(p) = total_revenue - total_cost;
all_emission(p) = total_emission_ton;
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
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
% ===== 自适应变异概率函数 =====
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, ~)
avg_demand = mean(aircraft.seats) * mean(flights.plf);
for i = 1:size(new_population, 1)
if rand() < pm(i)
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)
if rand() < 0.8
seating_capacities = aircraft.seats(valid_aircraft) .* flights.plf(flight_idx);
diffs = abs(seating_capacities - avg_demand);
[~, min_idx] = min(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;
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