数学建模学习（111）：改进遗传算法（引入模拟退火、轮盘赌和网格搜索）求解JSP问题-EW帮帮网

文章目录

一、车间调度问题
- 1.1目前处理方法
- 1.2简单案例
二、基于改进遗传算法求解车间调度
三、本文代码创新点
四、附上完整代码

一、车间调度问题

车间作业调度问题（JSP）是生产调度领域的经典问题之一，广泛应用于制造业、物流等领域。在一个典型的车间中，多个工件需要在多台机器上按特定的顺序进行加工。每个工件的加工步骤及其顺序是预先确定的，不同工件的加工顺序可能不同。 车间作业调度问题的目标是确定工件在各机器上的加工顺序，以最小化所有工件的完成总时间（即 makespan）。

1.1目前处理方法

车间作业调度问题因其高度的复杂性和计算上的挑战性，被归类为NP难问题。为了解决这一问题，学术界和工业界已经提出了多种方法，涵盖了精确算法、启发式算法以及元启发式算法。
精确算法，如分支定界法和动态规划，理论上能够找到全局最优解，但由于其高计算复杂度，通常不适用于大规模问题。
启发式算法，例如优先规则，如最短加工时间优先（SPT）和最早截止时间优先（EDD），通过简单的规则快速生成可行解。尽管这些方法计算速度快，但解的质量往往无法得到保证。
元启发式算法，包括遗传算法（GA）、模拟退火（SA）和粒子群优化（PSO），在处理复杂优化问题时表现出色。遗传算法通过模拟自然选择和遗传机制来不断优化解的种群；模拟退火算法通过模拟物质冷却过程中的能量变化来避免局部最优；粒子群优化算法则通过模拟鸟群觅食行为在多维空间中搜索全局最优解。
近年来，混合算法也逐渐受到关注，例如结合遗传算法和禁忌算法的混合算法，能够利用不同算法的优势，提高解的质量和求解效率。

1.2简单案例

我们假设有三个机器（M1, M2, M3）和三个工件（J1, J2, J3），每个工件已有不同的加工步骤和时间。这是一个简单的示例：
假设加工步骤和时间如下：
J1: M1(3) → M2(2) → M3(2)
J2: M2(4) → M3(1) → M1(3)
J3: M3(2) → M1(4) → M2(3)
以时间为x轴，每个机器为y轴，利用Python求解，每种工件代表1种颜色，可绘制甘特图如下所示：
请添加图片描述
绘图相应代码如下：

import matplotlib.pyplot as plt
import random
plt.rcParams['font.sans-serif'] = ['SimHei']
# 创建一个表示任务的列表，每个任务用一个字典表示
tasks = [
    {"Task": "J1", "Machine": "M1", "Start": 0, "Finish": 3},
    {"Task": "J1", "Machine": "M2", "Start": 3, "Finish": 5},
    {"Task": "J1", "Machine": "M3", "Start": 5, "Finish": 7},
    {"Task": "J2", "Machine": "M2", "Start": 0, "Finish": 4},
    {"Task": "J2", "Machine": "M3", "Start": 4, "Finish": 5},
    {"Task": "J2", "Machine": "M1", "Start": 5, "Finish": 8},
    {"Task": "J3", "Machine": "M3", "Start": 0, "Finish": 2},
    {"Task": "J3", "Machine": "M1", "Start": 2, "Finish": 6},
    {"Task": "J3", "Machine": "M2", "Start": 6, "Finish": 9},
]
# 创建一个新的图形
fig, gnt = plt.subplots()

# 设置Y轴
gnt.set_ylim(0, 40)
gnt.set_xlim(0, 12)
# 设置Y轴标签
gnt.set_yticks([10, 20, 30])
gnt.set_yticklabels(['M1', 'M2', 'M3'])
# 设置X轴标签
gnt.set_xlabel('时间')
gnt.set_ylabel('机器')
# 生成不同的颜色，每个工件一个颜色
task_colors = {}
for task in tasks:
    if task["Task"] not in task_colors:
        task_colors[task["Task"]] = (random.random(), random.random(), random.random())
# 绘制任务
for task in tasks:
    machine_map = {"M1": 10, "M2": 20, "M3": 30}
    start = task["Start"]
    finish = task["Finish"]
    machine = machine_map[task["Machine"]]
    color = task_colors[task["Task"]]
    gnt.broken_barh([(start, finish - start)], (machine - 5, 9), facecolors=color)
# 显示甘特图
plt.show()

二、基于改进遗传算法求解车间调度

2.1车间调度背景介绍

在一个车间内，有10台机器，每台机器负责一道加工步骤。这些机器需要完成10个工件的加工，每个工件都有10个加工步骤。不同工件的加工步骤顺序各不相同且顺序不可更改。每个加工步骤由对应的机器在一定时间内完成。你需要确定各个工件在不同机器上的加工顺序，以最小化所有工件完成加工所需的总时间。具体的不同工件加工顺序与加工时长如下表所示:
在这里插入图片描述
对于每个工件，其加工顺序和每台机器上的加工耗时一定，拿第一行举例。工件1的加工顺序为：M1-M2-M4-M3-M5-M6-M8-M4-M7-M10。其花费时间固定，分别为28-33-11-48-18-19-86-64-65-90。我们要做的就是对机器进行调度，什么时候加工第几个工件的第几个工序。

2.2遗传算法介绍

遗传算法（Genetic Algorithm，GA）是一类模拟自然选择和遗传机制的进化算法，主要用于优化和搜索问题。它最早由约翰·霍兰德（John Holland）在20世纪70年代提出。遗传算法通过模拟自然界生物进化的过程，包括选择、交叉（杂交）和变异等操作，逐步优化问题的解。

2.2.1基本流程

初始化种群：随机生成一定数量的初始解，称为个体，这些个体构成初始种群。
适应度评估：计算每个个体的适应度值，适应度值越高表示该个体越适合问题的求解。
选择：根据个体的适应度值，选择一些个体作为父代。常用的选择方法有轮盘赌选择、锦标赛选择等。
交叉：通过交叉操作（即杂交）生成新的个体，交叉操作将两个父代个体的部分基因交换，从而生成子代个体。
变异：对新生成的个体进行变异操作，即随机改变个体的一部分基因，以增加种群的多样性。
生成新种群：用子代个体替换部分或全部父代个体，形成新一代种群。
重复迭代：重复上述步骤，直到满足终止条件，如达到最大迭代次数或种群中的最佳适应度值达到预期目标。

2.2.2遗传算法的基本操作和公式

选择操作：选择操作通过适应度值选择个体。轮盘赌选择法的概率计算公式：
$P_i = \frac{f_i}{\sum_{j=1}^{N} f_j}$
其中，(P_i) 是第 (i) 个个体被选中的概率，(f_i) 是第 (i) 个个体的适应度值，(N) 是种群的大小。
交叉操作：交叉操作将两个父代个体的基因部分交换，生成两个新的个体。单点交叉的公式：
$\begin{aligned} &\text{子代1} = (\text{父代1}_1, \text{父代1}_2, \ldots, \text{父代1}_c, \text{父代2}_{c+1}, \ldots, \text{父代2}_n) \\ &\text{子代2} = (\text{父代2}_1, \text{父代2}_2, \ldots, \text{父代2}_c, \text{父代1}_{c+1}, \ldots, \text{父代1}_n) \end{aligned}$
其中，(c) 是交叉点的位置，(\text{父代1}) 和 (\text{父代2}) 是两个父代个体，(n) 是个体的基因长度。
变异操作：变异操作随机改变个体的一部分基因。对于一个基因序列 (x = (x_1, x_2, \ldots, x_n))，其变异公式可以表示为：
$x_i' = \begin{cases} \text{随机值} & \text{如果 } \text{随机概率} < \text{变异概率} \\ x_i & \text{否则} \end{cases}$
其中，(x_i’) 是变异后的基因值，变异概率通常设为一个小值，如0.01或0.1。

2.2.3遗传算法的优势

全局搜索能力强：遗传算法通过选择、交叉和变异等操作，可以在较大范围内进行搜索，能够较好地避免陷入局部最优。
适用范围广：遗传算法适用于各种优化问题，包括离散问题和连续问题。
易于并行化：由于种群中的个体可以并行计算适应度值，遗传算法易于实现并行计算，提升计算效率。

2.2.4遗传算法的不足

计算开销大：遗传算法需要多次迭代，且每次迭代需要计算多个个体的适应度值，计算开销较大。
参数选择复杂：遗传算法的性能对参数（如种群大小、交叉率、变异率等）依赖较大，参数选择不当可能影响算法效果。
收敛速度慢：在某些情况下，遗传算法的收敛速度较慢，可能需要较多的迭代次数才能找到较优的解。

2.3讲解本文思路及代码

本文改进遗传算法的代码更侧重于结合多种优化算法解决问题。利用网格搜索找到最优参数，并且在搜索过程中结合模拟退火算法进行局部优化，探索了更大范围的解空间，寻找最优解的可能性更大，利用上述显著优势尝试求解JSP问题。
Step1：初始化。创建一个初始种群，每个个体代表一个可能的作业顺序。通过随机生成的方式确保种群的多样性。

def createPop(Jm, popSize):
    pop = []
    for i in range(popSize):
        pop.append(createInd(Jm))
    return pop

Step2：选择。使用轮盘赌选择方法，根据个体的适应度（即完工时间）概率性地选择个体。适应度越高（即完工时间越短），被选择的概率越大。

def roulette_wheel_selection(fitness):
    total_fitness = sum(fitness)
    pick = random.uniform(0, total_fitness)
    current = 0
    for i, fit in enumerate(fitness):
        current += fit
        if current > pick:
            return i

Step3：交叉。应用自定义的交叉操作，将父代个体组合生成新的后代。通过选择两个父代个体，交换它们的一部分基因片段，以继承父母双方的特征。

def cross(A, B):
    n = len(A)
    r1 = np.random.randint(n)
    r2 = np.random.randint(n)
    rl, rr = min(r1, r2), max(r1, r2)
    if rl == rr:
        return A, B
    bt = copy.deepcopy(B)
    afinal = copy.deepcopy(A)
    for i in range(rl, rr + 1):
        bt.remove(A[i])
    k = 0
    for i in range(n):
        if i < rl or i > rr:
            afinal[i] = bt[k]
            k += 1
    at = copy.deepcopy(A)
    bfinal = copy.deepcopy(B)
    for i in range(rl, rr + 1):
        at.remove(B[i])
    k = 0
    for i in range(n):
        if i < rl or i > rr:
            bfinal[i] = at[k]
            k += 1
    return afinal, bfinal

Step4：变异。通过随机交换作业操作引入变异，以保持基因多样性。变异操作有助于防止算法陷入局部最优。

def mutate(individual, mutation_rate):
    if random.random() < mutation_rate:
        index1, index2 = random.sample(range(len(individual)), 2)
        individual[index1], individual[index2] = individual[index2], individual[index1]
    return individual

Step5：解码。使用详细的调度算法解码每个个体，以评估其完工时间，确保所有作业约束得到满足。解码过程包括为每个作业在特定机器上安排开始和结束时间。

def decode(J, P, s):
    n, m = J.shape
    T = [[[0]] for _ in range(m)]
    C = np.zeros((n, m))
    k = np.zeros(n, dtype=int)
    for job in s:
        if job < 0 or job >= n:
            print(f"Invalid job index: {job}")
            continue
        if k[job] < 0 or k[job] >= m:
            print(f"Invalid operation index for job {job}: {k[job]}")
            continue
        machine = J[job, k[job]] - 1
        process_time = P[job, k[job]]
        last_job_finish = C[job, k[job] - 1] if k[job] > 0 else 0
        start_time = max(last_job_finish, T[machine][-1][-1])
        insert_index = len(T[machine])
        for i in range(1, len(T[machine])):
            gap_start = max(last_job_finish, T[machine][i - 1][-1])
            gap_end = T[machine][i][0]
            if gap_end - gap_start >= process_time:
                start_time = gap_start
                insert_index = i
                break
        end_time = start_time + process_time
        C[job, k[job]] = end_time
        T[machine].insert(insert_index, [start_time, job, k[job], end_time])
        k[job] += 1
    M = [[] for _ in range(m)]
    for machine in range(m):
        for j in T[machine][1:]:
            M[machine].append(j[1])
    return T, M, C

Step6：模拟退火。为了增强搜索能力，本文结合模拟退火方法，对个体解进行微调，允许偶尔接受较差解，以跳出局部最优。模拟退火通过**逐渐降低“温度”**来减少接受较差解的概率，从而达到全局优化的效果。

def simulated_annealing(individual, J, P, T0, alpha, max_iter):
    current_solution = copy.deepcopy(individual)
    Tmax, _, current_fitness = decode(J, P, current_solution)
    best_solution = copy.deepcopy(current_solution)
    best_fitness = current_fitness.max()
    T = T0
    for i in range(max_iter):
        new_solution = mutate(copy.deepcopy(current_solution), 1.0)
        Tmax, _, new_fitness = decode(J, P, new_solution)
        delta_fitness = new_fitness.max() - current_fitness.max()
        if delta_fitness < 0 or random.random() < math.exp(-delta_fitness / T):
            current_solution = new_solution
            current_fitness = new_fitness
            if current_fitness.max() < best_fitness:
                best_solution = copy.deepcopy(current_solution)
                best_fitness = current_fitness.max()
        T *= alpha
    return best_solution, best_fitness

Step7：算法执行过程。通过初始化种群并不断执行选择、交叉、变异和解码操作，结合模拟退火优化，遗传算法逐步优化种群中的个体。每一代后，种群会逐步接近最优解，最终输出最优的调度方案。以下是主函数部分，展示了算法的执行过程：

# 主函数
if __name__ == "__main__":
    J, P = load_data('05.txt')
    pop_size = 50
    max_iter = 5000
    w = 0.5
    c1 = 1.5
    c2 = 1.5
    best_position, best_fitness = pso(J, P, pop_size, max_iter, w, c1, c2)
    print(f'Best Solution: {best_position}')
    print(f'Best Fitness: {best_fitness}')
    T, M, C = decode(J, P, best_position)
    drawGantt(T)
    plt.show()

2.4算法执行结果：

经过长时间运行，本文得到最少用时为 893s, 由网格搜索得到的具体参数如下所示：
在这里插入图片描述

三、本文代码创新点

创新1：引入混合优化算法。本文结合了遗传算法、模拟退火和网格搜索三种优化方法，提升了搜索效率和解的质量。遗传算法利用选择、交叉和变异操作进行全局搜索，模拟退火通过局部搜索进一步优化解，而网格搜索则用于寻找最佳的参数组合。
创新2：采用网格搜索。param_grid定义了一系列参数范围，通过遍历所有参数组合，自动化寻找最佳参数设置，确保算法在参数空间内找到最优解。
创新3：引入模拟退火策略。在遗传算法的基础上，引入模拟退火策略，有效地跳出局部最优解，进一步提升解的质量。通过控制退火过程的温度，平衡了探索与开发。
创新4：自适应交叉和变异操作。交叉和变异操作的实现更加灵活，能够根据个体特性进行调整。变异操作通过概率控制增加了种群的多样性，帮助避免过早收敛。
创新5：甘特图绘制。通过rawGantt函数，使用不同颜色表示不同工件的操作，并在甘特图上显示详细的操作信息，提供了直观的调度结果可视化。
创新6：基于轮盘赌选择。采用轮盘赌选择机制，增加了适应度高的个体被选中的概率，确保了优良基因的传递，提高了算法的收敛速度和效果。

四、附上完整代码

import copy
import numpy as np
import matplotlib.pyplot as plt
import random
import math
#------------------------------------------------------------------------------------------------------------------

#-------------------------------------------------------------------------------------
def createInd(J):
    #创建个体
    n = J.shape[0]
    s = []
    Jm = J.copy()
    while not np.all(Jm == 0):
        I = np.random.randint(0, n)
        M = Jm[I, 0]
        if M != 0:
            s.append(I)
            b = np.roll(Jm[I, :], -1)
            b[-1] = 0
            Jm[I, :] = b
    return s

def createPop(Jm, popSize):
    #创建种群
    pop = []
    for i in range(popSize):
        pop.append(createInd(Jm))
    return pop

def decode(J, P, s):
    #解码
    n, m = J.shape
    T = [[[0]] for _ in range(m)]
    C = np.zeros((n, m))
    k = np.zeros(n, dtype=int)
    for job in s:
        machine = J[job, k[job]] - 1
        process_time = P[job, k[job]]
        last_job_finish = C[job, k[job] - 1] if k[job] > 0 else 0
        start_time = max(last_job_finish, T[machine][-1][-1])
        insert_index = len(T[machine])
        for i in range(1, len(T[machine])):
            gap_start = max(last_job_finish, T[machine][i - 1][-1])
            gap_end = T[machine][i][0]
            if gap_end - gap_start >= process_time:
                start_time = gap_start
                insert_index = i
                break
        end_time = start_time + process_time
        C[job, k[job]] = end_time
        T[machine].insert(insert_index, [start_time, job, k[job], end_time])
        k[job] += 1
    M = [[] for _ in range(m)]
    for machine in range(m):
        for j in T[machine][1:]:
            M[machine].append(j[1])
    return T, M, C

def drawGantt(timelist):
    #绘制甘特图
    T = timelist.copy()
    plt.rcParams['font.sans-serif'] = ['SimHei']
    fig, ax = plt.subplots(figsize=(10, 6))
    color_map = {}
    for machine in T:
        for task_data in machine[1:]:
            job_idx, operation_idx = task_data[1], task_data[2]
            if job_idx not in color_map:
                color_map[job_idx] = (random.random(), random.random(), random.random())
    for machine_idx, machine_schedule in enumerate(T):
        for task_data in machine_schedule[1:]:
            start_time, job_idx, operation_idx, end_time = task_data
            color = color_map[job_idx]
            ax.barh(machine_idx, end_time - start_time, left=start_time, height=0.4, color=color)
            label = f'{job_idx}-{operation_idx}'
            ax.text((start_time + end_time) / 2, machine_idx, label, ha='center', va='center', color='white', fontsize=10)
    ax.set_yticks(range(len(T)))
    ax.set_yticklabels([f'machine{i + 1}' for i in range(len(T))])
    plt.xlabel("时间")
    plt.title("JSP甘特图")
    l = []
    for job_idx, color in color_map.items():
        l.append(plt.Rectangle((0, 0), 1, 1, color=color, label=f'{job_idx}'))
    plt.legend(handles=l, title='工件')

def cross(A, B):
    #交叉操作
    n = len(A)
    r1 = np.random.randint(n)
    r2 = np.random.randint(n)
    rl, rr = min(r1, r2), max(r1, r2)
    if rl == rr:
        return A, B
    bt = copy.deepcopy(B)
    afinal = copy.deepcopy(A)
    for i in range(rl, rr + 1):
        bt.remove(A[i])
    k = 0
    for i in range(n):
        if i < rl or i > rr:
            afinal[i] = bt[k]
            k += 1
    at = copy.deepcopy(A)
    bfinal = copy.deepcopy(B)
    for i in range(rl, rr + 1):
        at.remove(B[i])
    k = 0
    for i in range(n):
        if i < rl or i > rr:
            bfinal[i] = at[k]
            k += 1
    return afinal, bfinal

def load_data(path):
    #载入操作
    with open(path, 'r') as file:
        lines = file.readlines()
    job_num, machines_num = map(int, lines[0].split())
    J = np.zeros((job_num, len(lines[1].split()) // 2), dtype=int)
    P = np.zeros((job_num, len(lines[1].split()) // 2), dtype=int)
    for i in range(1, len(lines)):
        data = list(map(int, lines[i].split()))
        for j in range(len(data)):
            if j % 2 == 0:
                J[i - 1][j // 2] = data[j]
            else:
                P[i - 1][j // 2] = data[j]
    return J, P

def mutate(individual, mutation_rate):
    #变异操作
    if random.random() < mutation_rate:
        index1, index2 = random.sample(range(len(individual)), 2)
        individual[index1], individual[index2] = individual[index2], individual[index1]
    return individual

def roulette_wheel_selection(fitness):
    #轮盘赌
    total_fitness = sum(fitness)
    pick = random.uniform(0, total_fitness)
    current = 0
    for i, fit in enumerate(fitness):
        current += fit
        if current > pick:
            return i

def simulated_annealing(individual, J, P, T0, alpha, max_iter):
    #模拟退火
    current_solution = copy.deepcopy(individual)
    Tmax, _, current_fitness = decode(J, P, current_solution)
    best_solution = copy.deepcopy(current_solution)
    best_fitness = current_fitness.max()
    T = T0
    for i in range(max_iter):
        new_solution = mutate(copy.deepcopy(current_solution), 1.0)
        Tmax, _, new_fitness = decode(J, P, new_solution)
        delta_fitness = new_fitness.max() - current_fitness.max()
        if delta_fitness < 0 or random.random() < math.exp(-delta_fitness / T):
            current_solution = new_solution
            current_fitness = new_fitness
            if current_fitness.max() < best_fitness:
                best_solution = copy.deepcopy(current_solution)
                best_fitness = current_fitness.max()
        T *= alpha
    return best_solution, best_fitness

param_grid = {
    #网格搜索
    'pop_size': [300, 400, 500],
    'mutation_rate': [ 0.5, 0.7, 0.9],
    'T0': [ 1500, 2000],
    'alpha': [0.95, 0.97, 0.99],
    'max_iter': [50, 100, 150, 200]
}

# 加载数据
J, P = load_data('05.txt')
n, m = J.shape

# 初始化最优参数和最优结果
best_params = {}
best_Cmax = float('inf')

# 网格搜索循环
for pop_size in param_grid['pop_size']:
    for mutation_rate in param_grid['mutation_rate']:
        for T0 in param_grid['T0']:
            for alpha in param_grid['alpha']:
                for max_iter in param_grid['max_iter']:
                    print(f'Testing parameters: pop_size={pop_size}, mutation_rate={mutation_rate}, T0={T0}, alpha={alpha}, max_iter={max_iter}')
                    pop = createPop(J, pop_size)
                    Tmax, _, C = decode(J, P, pop[0])
                    fitness = [C.max()]
                    Cmax = C.max()
                    bestInd = copy.deepcopy(pop[0])
                    for i in range(1, pop_size):
                        T_, _, C = decode(J, P, pop[i])
                        if C.max() < Cmax:
                            Tmax = T_
                            Cmax = C.max()
                            bestInd = copy.deepcopy(pop[i])
                        fitness.append(C.max())
                    g = 0
                    gen = n * m
                    while g < gen:
                        g += 1
                        newInd = []
                        newFitness = []
                        l = 0
                        while l < pop_size / 2:
                            l += 1
                            tm = roulette_wheel_selection(fitness)
                            l1, l2 = cross(pop[tm], bestInd)
                            T1, _, C1 = decode(J, P, l1)
                            newInd.append(l1)
                            newFitness.append(C1.max())
                            if C1.max() < Cmax:
                                Tmax = T1
                                Cmax = C1.max()
                                bestInd = copy.deepcopy(l1)
                            T2, _, C2 = decode(J, P, l2)
                            newInd.append(l2)
                            newFitness.append(C2.max())
                            if C2.max() < Cmax:
                                Tmax = T2
                                Cmax = C2.max()
                                bestInd = copy.deepcopy(l2)
                        newpop = pop + newInd
                        newFit = fitness + newFitness
                        newId = np.array(newFit).argsort()[:pop_size]
                        pop = copy.deepcopy([newpop[i] for i in newId])
                        fitness = [newFit[i] for i in newId]
                        for i in range(pop_size):
                            pop[i] = mutate(pop[i], mutation_rate)
                            Ind = copy.deepcopy(pop[i])
                            Tt, _, Ct = decode(J, P, Ind)
                            fitness[i] = Ct.max()
                            if Ct.max() < Cmax:
                                Tmax = Tt
                                Cmax = Ct.max()
                                bestInd = copy.deepcopy(Ind)
                        # 模拟退火
                        for i in range(pop_size):
                            pop[i], fitness[i] = simulated_annealing(pop[i], J, P, T0, alpha, max_iter)
                            if fitness[i] < Cmax:
                                Tmax, _, _ = decode(J, P, pop[i])
                                Cmax = fitness[i]
                                bestInd = copy.deepcopy(pop[i])
                        print(f'第{g}代, Cmax={Cmax}')
                    if Cmax < best_Cmax:
                        best_Cmax = Cmax
                        best_params = {
                            'pop_size': pop_size,
                            'mutation_rate': mutation_rate,
                            'T0': T0,
                            'alpha': alpha,
                            'max_iter': max_iter
                        }

print(f'Best parameters: {best_params}')
print(f'Best Cmax: {best_Cmax}')

# 使用最优参数重新运行算法
pop_size = best_params['pop_size']
mutation_rate = best_params['mutation_rate']
T0 = best_params['T0']
alpha = best_params['alpha']
max_iter = best_params['max_iter']
pop = createPop(J, pop_size)
Tmax, _, C = decode(J, P, pop[0])
fitness = [C.max()]
Cmax = C.max()
bestInd = copy.deepcopy(pop[0])
for i in range(1, pop_size):
    T_, _, C = decode(J, P, pop[i])
    if C.max() < Cmax:
        Tmax = T_
        Cmax = C.max()
        bestInd = copy.deepcopy(pop[i])
    fitness.append(C.max())
g = 0
gen = n * m
chistory = []
while g < gen:
    g += 1
    newInd = []
    newFitness = []
    l = 0
    while l < pop_size / 2:
        l += 1
        tm = roulette_wheel_selection(fitness)
        l1, l2 = cross(pop[tm], bestInd)
        T1, _, C1 = decode(J, P, l1)
        newInd.append(l1)
        newFitness.append(C1.max())
        if C1.max() < Cmax:
            Tmax = T1
            Cmax = C1.max()
            bestInd = copy.deepcopy(l1)
        T2, _, C2 = decode(J, P, l2)
        newInd.append(l2)
        newFitness.append(C2.max())
        if C2.max() < Cmax:
            Tmax = T2
            Cmax = C2.max()
            bestInd = copy.deepcopy(l2)
    newpop = pop + newInd
    newFit = fitness + newFitness
    newId = np.array(newFit).argsort()[:pop_size]
    pop = copy.deepcopy([newpop[i] for i in newId])
    fitness = [newFit[i] for i in newId]
    for i in range(pop_size):
        pop[i] = mutate(pop[i], mutation_rate)
        Ind = copy.deepcopy(pop[i])
        Tt, _, Ct = decode(J, P, Ind)
        fitness[i] = Ct.max()
        if Ct.max() < Cmax:
            Tmax = Tt
            Cmax = Ct.max()
            bestInd = copy.deepcopy(Ind)
    # 模拟退火
    for i in range(pop_size):
        pop[i], fitness[i] = simulated_annealing(pop[i], J, P, T0, alpha, max_iter)
        if fitness[i] < Cmax:
            Tmax, _, _ = decode(J, P, pop[i])
            Cmax = fitness[i]
            bestInd = copy.deepcopy(pop[i])
    print(f'第{g}代, Cmax={Cmax}')
    chistory.append(Cmax)
plt.plot(chistory)
plt.show()
drawGantt(Tmax)
plt.show()

数学建模学习（111）：改进遗传算法（引入模拟退火、轮盘赌和网格搜索）求解JSP问题