基于NSGA-II的多目标柔性车间调度算法

一、算法实现

%% 基于NSGA-II的多目标柔性车间调度算法
% 问题：柔性作业车间调度问题（FJSP）
% 目标：最小化最大完工时间（makespan）和总机器负载（total machine load）
% 算法：NSGA-II（非支配排序遗传算法II）

clear; clc; close all;
warning('off', 'all');

%% 1. 问题定义与参数设置
disp('=== 多目标柔性车间调度系统 ===');
disp('初始化参数...');

% 调度问题参数
problem = struct();
problem.n_jobs = 6;           % 工件数量
problem.n_machines = 5;       % 机器数量
problem.n_operations = [3, 4, 3, 2, 4, 3];  % 每个工件的工序数
problem.max_ops = max(problem.n_operations); % 最大工序数

% 生成随机加工时间矩阵（工序×机器）
% 如果机器不可用，加工时间为inf
problem.processing_times = cell(problem.n_jobs, 1);
for i = 1:problem.n_jobs
    problem.processing_times{i} = zeros(problem.n_operations(i), problem.n_machines);
    for j = 1:problem.n_operations(i)
        % 每道工序可在部分机器上加工
        available_machines = randperm(problem.n_machines);
        n_available = randi([ceil(problem.n_machines*0.6), problem.n_machines]);
        available_machines = available_machines(1:n_available);
        
        for k = 1:problem.n_machines
            if ismember(k, available_machines)
                problem.processing_times{i}(j,k) = randi([5, 30]); % 加工时间5-30
            else
                problem.processing_times{i}(j,k) = inf; % 不可用
            end
        end
    end
end

% NSGA-II算法参数
params = struct();
params.pop_size = 100;        % 种群大小
params.max_gen = 200;         % 最大迭代次数
params.pc = 0.9;              % 交叉概率
params.pm = 0.1;              % 变异概率
params.n_obj = 2;             % 目标函数数量
params.eta_c = 20;            % 交叉分布指数
params.eta_m = 20;            % 变异分布指数
params.k = 2;                 % 锦标赛选择参数

% 显示问题信息
fprintf('柔性车间调度问题:\n');
fprintf('  工件数量: %d\n', problem.n_jobs);
fprintf('  机器数量: %d\n', problem.n_machines);
fprintf('  总工序数: %d\n', sum(problem.n_operations));
fprintf('  种群大小: %d\n', params.pop_size);
fprintf('  最大代数: %d\n', params.max_gen);

%% 2. 初始化种群
disp('初始化种群...');

% 编码方案：两层编码
% 第一层：工序顺序编码（OS）
% 第二层：机器选择编码（MS）

population = initialize_population(problem, params);

% 计算初始种群的适应度
for i = 1:params.pop_size
    population(i).fitness = evaluate_fitness(population(i), problem);
    population(i).rank = 0;
    population(i).crowding_distance = 0;
end

% 显示初始种群信息
fprintf('初始种群已生成，大小为: %d\n', params.pop_size);

%% 3. NSGA-II主循环
disp('开始NSGA-II优化...');

% 存储每代的最佳解
best_solutions = cell(params.max_gen, 1);
pareto_fronts = cell(params.max_gen, 1);
convergence = zeros(params.max_gen, 2); % 存储两个目标的收敛情况

for gen = 1:params.max_gen
    fprintf('第 %d 代/%d...\n', gen, params.max_gen);
    
    % 3.1 非支配排序
    population = non_dominated_sort(population);
    
    % 3.2 计算拥挤度距离
    population = crowding_distance_assignment(population);
    
    % 3.3 选择操作（锦标赛选择）
    parents = tournament_selection(population, params);
    
    % 3.4 交叉操作
    offspring = crossover(parents, problem, params);
    
    % 3.5 变异操作
    offspring = mutation(offspring, problem, params);
    
    % 3.6 计算子代适应度
    for i = 1:length(offspring)
        offspring(i).fitness = evaluate_fitness(offspring(i), problem);
        offspring(i).rank = 0;
        offspring(i).crowding_distance = 0;
    end
    
    % 3.7 合并父代和子代
    combined_pop = [population, offspring];
    
    % 3.8 环境选择（精英保留）
    population = environmental_selection(combined_pop, params);
    
    % 3.9 记录当前代的最佳解
    [best_solutions{gen}, pareto_fronts{gen}] = get_best_solutions(population);
    
    % 3.10 记录收敛情况
    convergence(gen, :) = get_convergence_metrics(population);
    
    % 3.11 显示进度
    if mod(gen, 20) == 0
        fprintf('  第 %d 代: makespan=%.2f, total_load=%.2f\n', ...
            gen, convergence(gen,1), convergence(gen,2));
    end
end

%% 4. 结果分析
disp('分析优化结果...');

% 4.1 获取最终Pareto前沿
final_pareto = get_final_pareto_front(population);

% 4.2 选择折中解（基于TOPSIS方法）
compromise_solution = select_compromise_solution(final_pareto);

% 4.3 解码最佳调度方案
best_schedule = decode_schedule(compromise_solution, problem);

% 4.4 计算性能指标
metrics = calculate_performance_metrics(population, final_pareto, convergence);

%% 5. 可视化结果
disp('生成可视化结果...');

% 5.1 收敛曲线
plot_convergence(convergence, params);

% 5.2 Pareto前沿
plot_pareto_front(final_pareto, compromise_solution);

% 5.3 甘特图
plot_gantt_chart(best_schedule, problem);

% 5.4 调度性能分析
plot_schedule_analysis(best_schedule, problem);

% 5.5 算法性能对比
plot_algorithm_comparison(metrics);

%% 6. 输出结果
disp('输出调度方案...');
print_schedule_results(best_schedule, problem, metrics);

fprintf('\n=== NSGA-II优化完成 ===\n');
fprintf('最佳调度方案:\n');
fprintf('  最大完工时间: %.2f\n', best_schedule.makespan);
fprintf('  总机器负载: %.2f\n', best_schedule.total_load);
fprintf('  机器利用率: %.2f%%\n', best_schedule.machine_utilization * 100);

二、种群初始化模块

function population = initialize_population(problem, params)
    % 初始化种群
    
    population = struct();
    
    for i = 1:params.pop_size
        % 1. 工序顺序编码（OS）
        % 基于工件的编码，每个工件号出现次数等于其工序数
        OS = [];
        for j = 1:problem.n_jobs
            OS = [OS, j * ones(1, problem.n_operations(j))];
        end
        OS = OS(randperm(length(OS))); % 随机排列
        
        % 2. 机器选择编码（MS）
        % 为每道工序选择可用的机器
        MS = zeros(1, length(OS));
        
        % 创建工件-工序计数器
        job_op_counter = ones(1, problem.n_jobs);
        
        for pos = 1:length(OS)
            job_id = OS(pos);
            op_id = job_op_counter(job_id);
            
            % 获取该工序可用的机器
            available_machines = find(problem.processing_times{job_id}(op_id, :) < inf);
            
            if isempty(available_machines)
                error('工序 %d-%d 没有可用机器!', job_id, op_id);
            end
            
            % 随机选择一台可用机器
            MS(pos) = available_machines(randi(length(available_machines)));
            
            % 更新工序计数器
            job_op_counter(job_id) = job_op_counter(job_id) + 1;
        end
        
        % 存储个体
        population(i).OS = OS;
        population(i).MS = MS;
        population(i).fitness = zeros(1, params.n_obj);
        population(i).rank = 0;
        population(i).crowding_distance = 0;
    end
end

三、适应度评估模块

function fitness = evaluate_fitness(individual, problem)
    % 评估个体的适应度（目标函数值）
    % 目标1：最小化最大完工时间（makespan）
    % 目标2：最小化总机器负载（total machine load）
    
    % 解码调度方案
    schedule = decode_schedule(individual, problem);
    
    % 计算目标函数值
    fitness = zeros(1, 2);
    
    % 目标1：最大完工时间
    fitness(1) = schedule.makespan;
    
    % 目标2：总机器负载
    fitness(2) = schedule.total_load;
end

function schedule = decode_schedule(individual, problem)
    % 解码调度方案
    
    OS = individual.OS;
    MS = individual.MS;
    
    n_machines = problem.n_machines;
    n_jobs = problem.n_jobs;
    
    % 初始化机器时间表
    machine_timetable = cell(n_machines, 1);
    for m = 1:n_machines
        machine_timetable{m} = []; % 每个机器存储[开始时间, 结束时间, 工件, 工序]
    end
    
    % 初始化工件进度
    job_progress = zeros(1, n_jobs); % 每个工件已完成的工序数
    job_completion = zeros(1, n_jobs); % 每个工件的完成时间
    
    % 按照工序顺序解码
    for pos = 1:length(OS)
        job_id = OS(pos);
        op_id = job_progress(job_id) + 1;
        machine_id = MS(pos);
        
        % 获取加工时间
        processing_time = problem.processing_times{job_id}(op_id, machine_id);
        
        if isinf(processing_time)
            error('无效的机器选择!');
        end
        
        % 计算工序的开始时间
        % 约束1：同一工件的工序顺序约束
        job_ready_time = job_completion(job_id);
        
        % 约束2：机器的可用时间
        machine_schedule = machine_timetable{machine_id};
        
        if isempty(machine_schedule)
            % 机器空闲
            start_time = job_ready_time;
            end_time = start_time + processing_time;
        else
            % 找到机器上的合适时间窗口
            start_time = max(job_ready_time, machine_schedule(end, 2));
            
            % 检查是否有更早的空闲时间窗口
            for i = 1:size(machine_schedule, 1)-1
                gap_start = machine_schedule(i, 2);
                gap_end = machine_schedule(i+1, 1);
                gap_duration = gap_end - gap_start;
                
                if gap_duration >= processing_time && gap_start >= job_ready_time
                    if gap_start < start_time
                        start_time = gap_start;
                    end
                end
            end
            
            end_time = start_time + processing_time;
        end
        
        % 更新机器时间表
        machine_timetable{machine_id} = [machine_timetable{machine_id}; 
                                         start_time, end_time, job_id, op_id];
        machine_timetable{machine_id} = sortrows(machine_timetable{machine_id}, 1);
        
        % 更新工件进度
        job_progress(job_id) = op_id;
        job_completion(job_id) = end_time;
    end
    
    % 计算调度指标
    schedule = struct();
    schedule.OS = OS;
    schedule.MS = MS;
    schedule.machine_timetable = machine_timetable;
    schedule.job_completion = job_completion;
    
    % 最大完工时间
    schedule.makespan = max(job_completion);
    
    % 总机器负载
    total_load = 0;
    for m = 1:n_machines
        if ~isempty(machine_timetable{m})
            machine_load = sum(machine_timetable{m}(:,2) - machine_timetable{m}(:,1));
            total_load = total_load + machine_load;
        end
    end
    schedule.total_load = total_load;
    
    % 机器利用率
    schedule.machine_utilization = total_load / (schedule.makespan * n_machines);
    
    % 平均流程时间
    schedule.mean_flow_time = mean(job_completion);
    
    % 最大机器负载
    machine_loads = zeros(1, n_machines);
    for m = 1:n_machines
        if ~isempty(machine_timetable{m})
            machine_loads(m) = sum(machine_timetable{m}(:,2) - machine_timetable{m}(:,1));
        end
    end
    schedule.max_machine_load = max(machine_loads);
    
    % 机器负载均衡指标
    schedule.load_balance = std(machine_loads) / mean(machine_loads);
end

四、NSGA-II核心操作模块

function population = non_dominated_sort(population)
    % 快速非支配排序
    
    n_pop = length(population);
    
    % 初始化支配关系
    for i = 1:n_pop
        population(i).domination_set = [];
        population(i).dominated_count = 0;
        population(i).rank = inf;
    end
    
    % 计算支配关系
    for i = 1:n_pop
        for j = i+1:n_pop
            % 比较个体i和j
            if dominates(population(i).fitness, population(j).fitness)
                population(i).domination_set = [population(i).domination_set, j];
                population(j).dominated_count = population(j).dominated_count + 1;
            elseif dominates(population(j).fitness, population(i).fitness)
                population(j).domination_set = [population(j).domination_set, i];
                population(i).dominated_count = population(i).dominated_count + 1;
            end
        end
    end
    
    % 第一前沿
    front = [];
    for i = 1:n_pop
        if population(i).dominated_count == 0
            population(i).rank = 1;
            front = [front, i];
        end
    end
    
    % 后续前沿
    current_rank = 1;
    while ~isempty(front)
        next_front = [];
        
        for i = front
            for j = population(i).domination_set
                population(j).dominated_count = population(j).dominated_count - 1;
                
                if population(j).dominated_count == 0
                    population(j).rank = current_rank + 1;
                    next_front = [next_front, j];
                end
            end
        end
        
        front = next_front;
        current_rank = current_rank + 1;
    end
end

function result = dominates(f1, f2)
    % 判断f1是否支配f2（最小化问题）
    % f1支配f2当且仅当：对于所有目标，f1 <= f2，且至少一个目标f1 < f2
    
    result = all(f1 <= f2) && any(f1 < f2);
end

function population = crowding_distance_assignment(population)
    % 计算拥挤度距离
    
    n_pop = length(population);
    n_obj = length(population(1).fitness);
    
    % 按等级分组
    max_rank = max([population.rank]);
    
    for rank = 1:max_rank
        % 获取当前前沿的个体索引
        front_indices = find([population.rank] == rank);
        n_front = length(front_indices);
        
        if n_front == 0
            continue;
        end
        
        % 初始化拥挤度距离
        for idx = front_indices
            population(idx).crowding_distance = 0;
        end
        
        % 对每个目标计算拥挤度
        for obj = 1:n_obj
            % 按当前目标值排序
            [~, sorted_idx] = sort([population(front_indices).fitness]);
            sorted_indices = front_indices(sorted_idx);
            
            % 边界个体的拥挤度设为无穷大
            population(sorted_indices(1)).crowding_distance = inf;
            population(sorted_indices(end)).crowding_distance = inf;
            
            % 计算中间个体的拥挤度
            f_min = population(sorted_indices(1)).fitness(obj);
            f_max = population(sorted_indices(end)).fitness(obj);
            
            if abs(f_max - f_min) < 1e-10
                continue; % 避免除零
            end
            
            for i = 2:n_front-1
                idx = sorted_indices(i);
                idx_prev = sorted_indices(i-1);
                idx_next = sorted_indices(i+1);
                
                population(idx).crowding_distance = population(idx).crowding_distance + ...
                    (population(idx_next).fitness(obj) - population(idx_prev).fitness(obj)) / (f_max - f_min);
            end
        end
    end
end

function parents = tournament_selection(population, params)
    % 锦标赛选择
    
    n_pop = length(population);
    parents = [];
    
    while length(parents) < n_pop
        % 随机选择k个个体
        candidates = randperm(n_pop, params.k);
        
        % 选择最优个体
        best_idx = candidates(1);
        for i = 2:params.k
            idx = candidates(i);
            
            % 首先比较等级
            if population(idx).rank < population(best_idx).rank
                best_idx = idx;
            elseif population(idx).rank == population(best_idx).rank
                % 等级相同则比较拥挤度距离
                if population(idx).crowding_distance > population(best_idx).crowding_distance
                    best_idx = idx;
                end
            end
        end
        
        parents = [parents, population(best_idx)];
    end
end

五、遗传操作模块

function offspring = crossover(parents, problem, params)
    % 交叉操作
    
    n_parents = length(parents);
    offspring = [];
    
    for i = 1:2:n_parents-1
        if rand() < params.pc
            % 执行交叉
            parent1 = parents(i);
            parent2 = parents(i+1);
            
            % 1. 工序顺序编码（OS）交叉 - 基于工件的顺序交叉
            child1_OS = crossover_OS(parent1.OS, parent2.OS, problem);
            child2_OS = crossover_OS(parent2.OS, parent1.OS, problem);
            
            % 2. 机器选择编码（MS）交叉 - 两点交叉
            [child1_MS, child2_MS] = crossover_MS(parent1.MS, parent2.MS);
            
            % 创建子代个体
            child1 = struct();
            child1.OS = child1_OS;
            child1.MS = child1_MS;
            child1.fitness = zeros(1, params.n_obj);
            child1.rank = 0;
            child1.crowding_distance = 0;
            
            child2 = struct();
            child2.OS = child2_OS;
            child2.MS = child2_MS;
            child2.fitness = zeros(1, params.n_obj);
            child2.rank = 0;
            child2.crowding_distance = 0;
            
            offspring = [offspring, child1, child2];
        else
            % 不交叉，直接复制
            offspring = [offspring, parents(i), parents(i+1)];
        end
    end
end

function child_OS = crossover_OS(parent1_OS, parent2_OS, problem)
    % 基于工件的顺序交叉
    
    n_genes = length(parent1_OS);
    
    % 随机选择两个交叉点
    points = sort(randperm(n_genes, 2));
    start_point = points(1);
    end_point = points(2);
    
    % 初始化子代
    child_OS = zeros(1, n_genes);
    
    % 复制父代1的片段
    child_OS(start_point:end_point) = parent1_OS(start_point:end_point);
    
    % 从父代2填充剩余位置
    pos = 1;
    for i = 1:n_genes
        if pos == start_point
            pos = end_point + 1;
        end
        
        if pos > n_genes
            break;
        end
        
        gene = parent2_OS(i);
        
        % 检查该工件在子代中出现的次数
        job_count_in_child = sum(child_OS == gene);
        job_count_required = problem.n_operations(gene);
        
        if job_count_in_child < job_count_required
            child_OS(pos) = gene;
            pos = pos + 1;
        end
    end
end

function [child1_MS, child2_MS] = crossover_MS(parent1_MS, parent2_MS)
    % 机器选择编码的两点交叉
    
    n_genes = length(parent1_MS);
    
    % 随机选择两个交叉点
    points = sort(randperm(n_genes, 2));
    start_point = points(1);
    end_point = points(2);
    
    % 执行交叉
    child1_MS = parent1_MS;
    child2_MS = parent2_MS;
    
    child1_MS(start_point:end_point) = parent2_MS(start_point:end_point);
    child2_MS(start_point:end_point) = parent1_MS(start_point:end_point);
end

function offspring = mutation(offspring, problem, params)
    % 变异操作
    
    for i = 1:length(offspring)
        if rand() < params.pm
            % 执行变异
            individual = offspring(i);
            
            % 随机选择变异类型
            mutation_type = randi([1, 3]);
            
            switch mutation_type
                case 1
                    % 工序顺序交换变异
                    individual.OS = swap_mutation_OS(individual.OS);
                case 2
                    % 机器选择变异
                    individual.MS = machine_mutation(individual.MS, individual.OS, problem);
                case 3
                    % 同时变异
                    individual.OS = swap_mutation_OS(individual.OS);
                    individual.MS = machine_mutation(individual.MS, individual.OS, problem);
            end
            
            offspring(i) = individual;
        end
    end
end

function mutated_OS = swap_mutation_OS(OS)
    % 工序顺序交换变异
    
    n_genes = length(OS);
    
    % 随机选择两个不同的位置
    pos1 = randi(n_genes);
    pos2 = randi(n_genes);
    
    while pos1 == pos2
        pos2 = randi(n_genes);
    end
    
    % 交换基因
    mutated_OS = OS;
    mutated_OS([pos1, pos2]) = mutated_OS([pos2, pos1]);
end

function mutated_MS = machine_mutation(MS, OS, problem)
    % 机器选择变异
    
    n_genes = length(MS);
    
    % 随机选择一个位置进行变异
    pos = randi(n_genes);
    
    % 获取该位置对应的工件和工序
    job_id = OS(pos);
    
    % 计算这是该工件的第几道工序
    job_ops = find(OS == job_id);
    op_idx = find(job_ops == pos);
    
    % 获取该工序可用的机器
    available_machines = find(problem.processing_times{job_id}(op_idx, :) < inf);
    
    if length(available_machines) > 1
        % 随机选择另一台可用机器
        current_machine = MS(pos);
        available_machines = setdiff(available_machines, current_machine);
        
        if ~isempty(available_machines)
            new_machine = available_machines(randi(length(available_machines)));
            mutated_MS = MS;
            mutated_MS(pos) = new_machine;
        else
            mutated_MS = MS;
        end
    else
        mutated_MS = MS;
    end
end

参考代码 基于NSGA-2的求解多目标柔性车间调度算法 www.youwenfan.com/contentcnt/160664.html

六、环境选择模块

function new_population = environmental_selection(combined_pop, params)
    % 环境选择（精英保留策略）
    
    n_combined = length(combined_pop);
    n_pop = params.pop_size;
    
    % 非支配排序
    combined_pop = non_dominated_sort(combined_pop);
    
    % 计算拥挤度距离
    combined_pop = crowding_distance_assignment(combined_pop);
    
    % 按等级和拥挤度排序
    [~, sorted_idx] = sortrows([[combined_pop.rank]', -[combined_pop.crowding_distance]']);
    
    % 选择前pop_size个个体
    new_population = combined_pop(sorted_idx(1:n_pop));
end

function [best_solutions, pareto_front] = get_best_solutions(population)
    % 获取当前代的最佳解（Pareto前沿）
    
    % 找到第一前沿的个体
    front1_indices = find([population.rank] == 1);
    pareto_front = population(front1_indices);
    
    % 选择拥挤度最大的个体作为最佳解
    if ~isempty(pareto_front)
        [~, best_idx] = max([pareto_front.crowding_distance]);
        best_solutions = pareto_front(best_idx);
    else
        best_solutions = [];
        pareto_front = [];
    end
end

function convergence = get_convergence_metrics(population)
    % 获取收敛指标
    
    % 找到第一前沿
    front1_indices = find([population.rank] == 1);
    
    if isempty(front1_indices)
        convergence = [inf, inf];
        return;
    end
    
    front1 = population(front1_indices);
    fitness_values = [front1.fitness];
    
    % 计算两个目标的最小值
    convergence = [min(fitness_values(1,:)), min(fitness_values(2,:))];
end

七、结果分析与可视化模块

function final_pareto = get_final_pareto_front(population)
    % 获取最终的Pareto前沿
    
    % 非支配排序
    population = non_dominated_sort(population);
    
    % 获取第一前沿
    front1_indices = find([population.rank] == 1);
    final_pareto = population(front1_indices);
    
    % 按第一个目标排序
    [~, sorted_idx] = sort([final_pareto.fitness]);
    final_pareto = final_pareto(sorted_idx);
end

function compromise_solution = select_compromise_solution(pareto_front)
    % 使用TOPSIS方法选择折中解
    
    if isempty(pareto_front)
        compromise_solution = [];
        return;
    end
    
    n_solutions = length(pareto_front);
    n_obj = length(pareto_front(1).fitness);
    
    % 构建决策矩阵
    decision_matrix = zeros(n_solutions, n_obj);
    for i = 1:n_solutions
        decision_matrix(i, :) = pareto_front(i).fitness;
    end
    
    % 标准化决策矩阵
    norm_matrix = decision_matrix ./ sqrt(sum(decision_matrix.^2, 1));
    
    % 理想解和负理想解（最小化问题）
    ideal_solution = min(norm_matrix, [], 1);
    negative_ideal_solution = max(norm_matrix, [], 1);
    
    % 计算距离
    d_plus = sqrt(sum((norm_matrix - ideal_solution).^2, 2));
    d_minus = sqrt(sum((norm_matrix - negative_ideal_solution).^2, 2));
    
    % 计算相对接近度
    closeness = d_minus ./ (d_plus + d_minus);
    
    % 选择最接近理想解的方案
    [~, best_idx] = max(closeness);
    compromise_solution = pareto_front(best_idx);
end

function plot_convergence(convergence, params)
    % 绘制收敛曲线
    
    figure('Position', [100, 100, 1200, 500]);
    
    subplot(1,2,1);
    plot(1:params.max_gen, convergence(:,1), 'b-', 'LineWidth', 2);
    xlabel('迭代次数');
    ylabel('最大完工时间');
    title('最大完工时间收敛曲线');
    grid on;
    
    subplot(1,2,2);
    plot(1:params.max_gen, convergence(:,2), 'r-', 'LineWidth', 2);
    xlabel('迭代次数');
    ylabel('总机器负载');
    title('总机器负载收敛曲线');
    grid on;
    
    sgtitle('NSGA-II收敛性能');
end

function plot_pareto_front(pareto_front, compromise_solution)
    % 绘制Pareto前沿
    
    if isempty(pareto_front)
        return;
    end
    
    figure('Position', [200, 200, 800, 600]);
    
    % 提取目标值
    f1 = [pareto_front.fitness];
    makespan = f1(1:2:end);
    total_load = f1(2:2:end);
    
    % 绘制Pareto前沿
    scatter(makespan, total_load, 50, 'b', 'filled');
    hold on;
    
    % 绘制折中解
    if ~isempty(compromise_solution)
        scatter(compromise_solution.fitness(1), compromise_solution.fitness(2), ...
                100, 'r', 'filled', '^', 'LineWidth', 2);
        text(compromise_solution.fitness(1), compromise_solution.fitness(2), ...
             ' 折中解', 'FontSize', 12, 'FontWeight', 'bold');
    end
    
    % 连接Pareto点
    [sorted_makespan, sort_idx] = sort(makespan);
    sorted_load = total_load(sort_idx);
    plot(sorted_makespan, sorted_load, 'b--', 'LineWidth', 1);
    
    xlabel('最大完工时间 (makespan)');
    ylabel('总机器负载 (total load)');
    title('Pareto最优前沿');
    grid on;
    legend('Pareto解', '折中解', 'Location', 'best');
    
    % 添加统计信息
    text(0.05, 0.95, sprintf('Pareto解数量: %d', length(pareto_front)), ...
         'Units', 'normalized', 'FontSize', 10, 'BackgroundColor', 'white');
end

function plot_gantt_chart(schedule, problem)
    % 绘制甘特图
    
    figure('Position', [300, 300, 1000, 600]);
    
    colors = lines(problem.n_jobs); % 为每个工件分配不同颜色
    
    hold on;
    
    % 绘制每个机器的调度
    for m = 1:problem.n_machines
        machine_schedule = schedule.machine_timetable{m};
        
        if ~isempty(machine_schedule)
            for i = 1:size(machine_schedule, 1)
                start_time = machine_schedule(i, 1);
                end_time = machine_schedule(i, 2);
                job_id = machine_schedule(i, 3);
                op_id = machine_schedule(i, 4);
                
                % 绘制矩形
                rectangle('Position', [start_time, m-0.4, end_time-start_time, 0.8], ...
                          'FaceColor', colors(job_id,:), 'EdgeColor', 'k', 'LineWidth', 1);
                
                % 添加文本标签
                text(start_time + (end_time-start_time)/2, m, ...
                     sprintf('J%d-O%d', job_id, op_id), ...
                     'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', ...
                     'FontSize', 8, 'FontWeight', 'bold', 'Color', 'white');
            end
        end
    end
    
    % 设置图形属性
    ylabel('机器');
    xlabel('时间');
    title('柔性车间调度甘特图');
    ylim([0.5, problem.n_machines+0.5]);
    xlim([0, schedule.makespan * 1.1]);
    grid on;
    
    % 添加图例
    legend_handles = [];
    legend_labels = {};
    for j = 1:min(problem.n_jobs, 10) % 最多显示10个工件的图例
        h = fill(NaN, NaN, colors(j,:));
        legend_handles = [legend_handles, h];
        legend_labels = [legend_labels, sprintf('工件 %d', j)];
    end
    legend(legend_handles, legend_labels, 'Location', 'eastoutside');
    
    % 添加统计信息
    text(0.02, 0.98, sprintf('最大完工时间: %.2f', schedule.makespan), ...
         'Units', 'normalized', 'FontSize', 10, 'BackgroundColor', 'white');
    text(0.02, 0.94, sprintf('总机器负载: %.2f', schedule.total_load), ...
         'Units', 'normalized', 'FontSize', 10, 'BackgroundColor', 'white');
    text(0.02, 0.90, sprintf('机器利用率: %.1f%%', schedule.machine_utilization*100), ...
         'Units', 'normalized', 'FontSize', 10, 'BackgroundColor', 'white');
    
    hold off;
end

八、性能评估模块

function metrics = calculate_performance_metrics(population, pareto_front, convergence)
    % 计算算法性能指标
    
    metrics = struct();
    
    % 1. 超体积指标（Hypervolume）
    if ~isempty(pareto_front)
        % 提取目标值
        f_values = [pareto_front.fitness];
        f1 = f_values(1:2:end);
        f2 = f_values(2:2:end);
        
        % 参考点（比最差点稍差）
        ref_point = [max(f1)*1.1, max(f2)*1.1];
        
        % 计算超体积
        metrics.hypervolume = calculate_hypervolume([f1; f2]', ref_point);
    else
        metrics.hypervolume = 0;
    end
    
    % 2. 间距指标（Spacing）
    if length(pareto_front) > 1
        metrics.spacing = calculate_spacing(pareto_front);
    else
        metrics.spacing = 0;
    end
    
    % 3. 世代距离（Generational Distance）
    if ~isempty(pareto_front)
        % 使用第一代作为参考
        metrics.generational_distance = calculate_generational_distance(pareto_front);
    else
        metrics.generational_distance = 0;
    end
    
    % 4. 收敛性指标
    metrics.convergence_gen = find_convergence_generation(convergence);
    
    % 5. Pareto解的数量
    metrics.n_pareto_solutions = length(pareto_front);
    
    % 6. 目标函数范围
    if ~isempty(pareto_front)
        f_values = [pareto_front.fitness];
        metrics.makespan_range = [min(f_values(1:2:end)), max(f_values(1:2:end))];
        metrics.load_range = [min(f_values(2:2:end)), max(f_values(2:2:end))];
    else
        metrics.makespan_range = [0, 0];
        metrics.load_range = [0, 0];
    end
end

function hv = calculate_hypervolume(points, ref_point)
    % 计算超体积指标
    
    if isempty(points)
        hv = 0;
        return;
    end
    
    % 按第一个目标排序
    [~, idx] = sort(points(:,1));
    points = points(idx, :);
    
    % 计算超体积
    hv = 0;
    n_points = size(points, 1);
    
    for i = 1:n_points
        if i == 1
            width = ref_point(1) - points(i,1);
        else
            width = points(i-1,1) - points(i,1);
        end
        height = ref_point(2) - points(i,2);
        hv = hv + width * height;
    end
end

function spacing = calculate_spacing(pareto_front)
    % 计算间距指标
    
    n_solutions = length(pareto_front);
    
    if n_solutions < 2
        spacing = 0;
        return;
    end
    
    % 计算每个解到其他解的最小距离
    d = zeros(n_solutions, 1);
    
    for i = 1:n_solutions
        min_dist = inf;
        for j = 1:n_solutions
            if i ~= j
                dist = norm(pareto_front(i).fitness - pareto_front(j).fitness);
                if dist < min_dist
                    min_dist = dist;
                end
            end
        end
        d(i) = min_dist;
    end
    
    % 计算平均距离
    d_bar = mean(d);
    
    % 计算间距
    spacing = sqrt(sum((d - d_bar).^2) / (n_solutions - 1));
end

function gd = calculate_generational_distance(pareto_front)
    % 计算世代距离
    
    % 这里简化处理，使用理想点作为参考
    ideal_point = [min([pareto_front.fitness]), min([pareto_front.fitness])];
    
    n_solutions = length(pareto_front);
    distances = zeros(n_solutions, 1);
    
    for i = 1:n_solutions
        distances(i) = norm(pareto_front(i).fitness - ideal_point);
    end
    
    gd = mean(distances);
end

function conv_gen = find_convergence_generation(convergence)
    % 找到收敛的代数
    
    n_gen = size(convergence, 1);
    threshold = 0.001; % 收敛阈值
    
    for gen = 2:n_gen
        improvement1 = abs(convergence(gen,1) - convergence(gen-1,1)) / convergence(gen-1,1);
        improvement2 = abs(convergence(gen,2) - convergence(gen-1,2)) / convergence(gen-1,2);
        
        if improvement1 < threshold && improvement2 < threshold
            conv_gen = gen;
            return;
        end
    end
    
    conv_gen = n_gen; % 未收敛
end

九、调度结果输出模块

function print_schedule_results(schedule, problem, metrics)
    % 输出调度结果
    
    fprintf('\n=== 调度方案详情 ===\n\n');
    
    % 1. 基本调度信息
    fprintf('调度性能指标:\n');
    fprintf('  最大完工时间: %.2f\n', schedule.makespan);
    fprintf('  总机器负载: %.2f\n', schedule.total_load);
    fprintf('  机器利用率: %.2f%%\n', schedule.machine_utilization * 100);
    fprintf('  平均流程时间: %.2f\n', schedule.mean_flow_time);
    fprintf('  最大机器负载: %.2f\n', schedule.max_machine_load);
    fprintf('  负载均衡指标: %.4f\n', schedule.load_balance);
    
    fprintf('\n算法性能指标:\n');
    fprintf('  Pareto解数量: %d\n', metrics.n_pareto_solutions);
    fprintf('  超体积指标: %.4f\n', metrics.hypervolume);
    fprintf('  间距指标: %.4f\n', metrics.spacing);
    fprintf('  世代距离: %.4f\n', metrics.generational_distance);
    fprintf('  收敛代数: %d\n', metrics.convergence_gen);
    
    % 2. 机器调度详情
    fprintf('\n机器调度详情:\n');
    for m = 1:problem.n_machines
        machine_schedule = schedule.machine_timetable{m};
        
        if isempty(machine_schedule)
            fprintf('  机器 %d: 空闲\n', m);
        else
            fprintf('  机器 %d: ', m);
            total_time = 0;
            for i = 1:size(machine_schedule, 1)
                start_time = machine_schedule(i, 1);
                end_time = machine_schedule(i, 2);
                job_id = machine_schedule(i, 3);
                op_id = machine_schedule(i, 4);
                
                fprintf('J%d-O%d[%.1f-%.1f] ', job_id, op_id, start_time, end_time);
                total_time = total_time + (end_time - start_time);
            end
            fprintf('(利用率: %.1f%%)\n', total_time/schedule.makespan*100);
        end
    end
    
    % 3. 工件完成时间
    fprintf('\n工件完成时间:\n');
    for j = 1:problem.n_jobs
        fprintf('  工件 %d: %.2f\n', j, schedule.job_completion(j));
    end
    
    % 4. 编码信息
    fprintf('\n最优个体编码:\n');
    fprintf('  工序顺序: ');
    fprintf('%d ', schedule.OS);
    fprintf('\n  机器选择: ');
    fprintf('%d ', schedule.MS);
    fprintf('\n');
    
    % 5. 建议与优化
    fprintf('\n优化建议:\n');
    if schedule.machine_utilization < 0.7
        fprintf('  ⚠️ 机器利用率较低，考虑减少机器数量或增加工件\n');
    end
    
    if schedule.load_balance > 0.3
        fprintf('  ⚠️ 机器负载不均衡，考虑调整调度策略\n');
    end
    
    if schedule.makespan > 100
        fprintf('  ⚠️ 最大完工时间较长，考虑优化关键路径\n');
    end
    
    fprintf('\n调度完成!\n');
end

十、扩展功能：多算法对比

function compare_algorithms(problem, params)
    % 对比不同算法的性能
    
    algorithms = {'NSGA-II', 'MOEA/D', 'SPEA2'};
    n_algorithms = length(algorithms);
    
    results = cell(n_al