NSGA-II 求解 JSP+AGV调度 三目标联合优化问题

import numpy as np
import random
import matplotlib.pyplot as plt
import math
from itertools import chain

agv_num = 2

####工件的加工信息表格（Fisher and Thompson 6x6 instance mt06）####
# job_info=np.array([[2,1,0,3,1,6,3,7,5,3,4,6],
#                    [1,8,2,5,4,10,5,10,0,10,3,4],
#                    [2,5,3,4,5,8,0,9,1,1,4,7],
#                    [1,5,0,5,2,5,3,3,4,8,5,9],
#                    [2,9,1,3,4,5,5,4,0,3,3,1],
#                    [1,3,3,3,5,9,0,10,4,4,2,1]])

####工件的加工信息表格（Fisher and Thompson 5x6 instance mt06）####
# job_info=np.array([[2,1,0,3,1,6,3,7,5,3,4,6],
#                    [1,8,2,5,4,10,5,10,0,10,3,4],
#                    [2,5,3,4,5,8,0,9,1,1,4,7],
#                    [1,5,0,5,2,5,3,3,4,8,5,9],
#                    [2,9,1,3,4,5,5,4,0,3,3,1]])

####工件的加工信息表格（Fisher and Thompson 7x6 instance mt06）####
job_info=np.array([[2,1,0,3,1,6,3,7,5,3,4,6],
                   [1,8,2,5,4,10,5,10,0,10,3,4],
                   [2,5,3,4,5,8,0,9,1,1,4,7],
                   [1,5,0,5,2,5,3,3,4,8,5,9],
                   [2,9,1,3,4,5,5,4,0,3,3,1],
                   [2,9,1,3,4,5,5,4,0,3,3,1],
                   [1,3,3,3,5,9,0,10,4,4,2,1]
                   ])

job_num = job_info.shape[0] #工件数量  6
operation_num = int(job_info.shape[1]/2) #单个工件的工序数量  6

####工件的延迟信息表格 ####
# delay_info = np.array ([35,70,105,140,175,210])

####工件的延迟信息表格 ####
delay_info = np.array ([35,70,105,140,175,180,190])


####工件的利润信息表格 ####
# money_info = np.array([20,30,40,50,60,70])

####工件的利润信息表格 ####
money_info = np.array([20,30,40,50,60,70,80])

####AGV两点位运动时间表格####
agv_move_info=np.array([[0,4,6,8,14,12,10],
                        [10,0,3,5,11,9,7],
                        [12,15,0,3,9,7,9],
                        [14,17,15,0,7,9,11],
                        [8,11,9,7,0,3,5],
                        [6,9,7,9,15,0,3],
                        [4,7,9,11,17,15,0]])

mach_num = agv_move_info.shape[0]-1 # machine的数量  6


####调度方案生成####

def init(N_size):
    p_list = []
    pop = []

    for i in range(0, job_num):
        for q in range(operation_num):
            p_list.append(i)

    for j in range(0, N_size):
        pop.append(random.sample(p_list, job_num * operation_num))

    pop = np.array(pop)

    return pop

####agv运输时间计算####
def time_compute(start,end):
    return agv_move_info[start+1][end+1]


####完工时间计算####

# 最后呈现两张图
# 图1：machine的甘特图
# 图2：agv的甘特图 unloading / waiting / loading

def single_scheduling_time(scheduling_plan,show,count):  # 传入某一个个体的scheduling plan

    t = 0

    agv_state = list()
    makespan_list = job_num * [0]
    makespan_list = job_num * [0]
    delay_list = job_num * [0]
    profit_list = job_num * [0]

    # AGV初始状态标定
    for i in range(agv_num):
        agv_state.append(["waiting", -1])  # 第一个参数：区分AGV两种状态"loading"、"unloading"、"waiting"
        # 第二个参数：区分AGV的目的地，-1表示“L/U”，0表示machine 1...

    machine_timestamp = list()  # 6xN maxtrix
    for j in range(mach_num):
        machine_timestamp.append([])
        machine_timestamp[j].append(0)  # 第一个时刻为0

    agv_timestamp = list()  # 2xN matrix
    for j in range(agv_num):
        agv_timestamp.append([])
        agv_timestamp[j].append(0)  # 第一个时刻为0

    tag = np.zeros((job_num), dtype=int)  # 记录某一个job的第几次步骤  #[0,0,0,0,0,0] 编号0代表第一个步骤，编号5代表最后一个步骤

    for i in range(scheduling_plan.shape[0]):
        mach_p = job_info[scheduling_plan[i]][tag[scheduling_plan[i]] * 2]  # 该工序由mach_p号进行加工

        #         print("mach_p",mach_p)

        time_p = job_info[scheduling_plan[i]][tag[scheduling_plan[i]] * 2 + 1]  # 加工该工序所用的时间

        #         print("time_p",time_p)

        order_executiving(mach_p, time_p, scheduling_plan[i], tag, agv_state, agv_timestamp, machine_timestamp,
                          makespan_list, delay_list)

        tag[scheduling_plan[i]] += 1

    # 计算优化函数1：makespan 记录最大的完工时间

    #     print("makespan_list",makespan_list)

    makespan = max(makespan_list)

    # 计算优化函数2：mean_tardiness 记录平均延迟时间

    #     print("delay_list",delay_list)

    mean_tardiness = round(np.mean(np.array(delay_list)), 1)

    # 计算优化函数3：profit 记录利润

    profit_list = profit_cal(delay_list)
    profit = round(np.sum(np.array(profit_list)), 1)

    # 添加时间刻度,用于绘制甘特图

    max_machine_list = list()
    for j in range(len(machine_timestamp)):
        max_machine_list.append(max(machine_timestamp[j]))
    max_machine_t = max(max_machine_list)

    max_agv_list = list()
    for j in range(len(agv_timestamp)):
        max_agv_list.append(max(agv_timestamp[j]))
    max_agv_t = max(max_agv_list)

    if show == 1:
        plot_fig_agv(max_agv_t,agv_timestamp,count)
        plot_fig_machine(max_machine_t,machine_timestamp,count)

    return makespan, mean_tardiness, profit


#### 指数递减的利润计算函数 ####
def profit_cal(delay_list):
    profit = []

    num = 0

    for i in delay_list:
        p = money_info[num] * math.exp(-abs(i / 1000))

        profit.append(p)

        num += 1

    return profit


#### 交叉、变异算子 ####
crossover_mutation_rate = 0.3

#### 交叉算子 POX交叉####

def POX_improve(Father_1, Father_2):

    children_1 = np.zeros(job_num * operation_num, dtype=int)
    children_2 = np.zeros(job_num * operation_num, dtype=int)
    # 建立工件集合
    Job_set = [i_ for i_ in range(job_num)]  # 有job_num个工件
    J_set_1 = np.random.choice(Job_set, 2, replace=False)
    J_set_2 = np.random.choice(Job_set, 2, replace=False)
    index_J2_of_F1_set = []
    index_J1_of_F2_set = []
    for i in range(len(J_set_1)):
        index_J1_of_F1 = np.where(Father_1 == J_set_1[i])[0]
        index_J1_of_F2 = np.where(Father_2 == J_set_1[i])[0]
        index_J2_of_F1 = np.where(Father_1 == J_set_2[i])[0]
        index_J2_of_F2 = np.where(Father_2 == J_set_2[i])[0]
        index_J2_of_F1_set.append(index_J2_of_F1)
        index_J1_of_F2_set.append(index_J1_of_F2)
        for i_1 in range(np.size(index_J1_of_F1)):
            a_1 = index_J1_of_F1[i_1]
            a_2 = index_J2_of_F2[i_1]
            children_1[a_1] = Father_1[a_1]
            children_2[a_2] = Father_2[a_2]
            # 第二部将F1的转接到C2
    # 将二维数组转化为一维数组
    index_delet_1 = list(chain.from_iterable(index_J2_of_F1_set))
    index_delet_2 = list(chain.from_iterable(index_J1_of_F2_set))
    Father_1_deleted = np.delete(Father_1, index_delet_1)
    Father_2_deleted = np.delete(Father_2, index_delet_2)
    index_0_insert_1 = np.where(children_2 == 0)[0]
    index_0_insert_2 = np.where(children_1 == 0)[0]
    for j in range(np.size(Father_1_deleted)):
        a_3 = index_0_insert_1[j]
        a_4 = index_0_insert_2[j]
        children_2[a_3] = Father_1_deleted[j]
        children_1[a_4] = Father_2_deleted[j]
    return children_1, children_2

#### 变异算子 ###

def mut_job(DNA):
    DNA = DNA[:job_num * operation_num]
    DNA_1 = DNA[:job_num * operation_num].copy()
    mut_two_of_position = np.random.choice([i_ for i_ in range(job_num * operation_num)], 2, replace=False)
    a = mut_two_of_position[0]
    b = mut_two_of_position[1]
    DNA[b] = DNA_1[a]
    DNA[a] = DNA_1[b]
    return DNA

####AGV调度规则####


def order_executiving(mach_p, time_p, order, tag, agv_state, agv_timestamp, machine_timestamp, makespan_list,
                      delay_list):  # 参数1：目的地，参数2：job的第几个工序

    agv_distance = []

    dp = job_info[order][(tag[order] - 1) * 2]  # 此时工件所处的位置 demand point

    if (tag[order] == 0):  # 若是job的第1个工序

        for i in range(agv_num):  # 计算哪个AGV可以最快到达L/U
            agv_distance.append(time_compute(agv_state[i][1], -1) + agv_timestamp[i][-1])
        exe_agv = np.where(agv_distance == np.min(agv_distance))[0][0]

        #         print("exe_agv",exe_agv)

        agv_move_unloading(exe_agv, agv_state[exe_agv][1], -1, agv_timestamp,
                           agv_state)  # 将AGV从current position 移动到 demand point(L/U)

        makespan_list[order] = agv_timestamp[exe_agv][-1]  # 记录job开始加工的时间

    else:  # 若不是job的第1个工序

        for i in range(agv_num):  # 计算哪个AGV可以最快到达demand point
            agv_distance.append(time_compute(agv_state[i][1], dp) + agv_timestamp[i][-1])
        exe_agv = np.where(agv_distance == np.min(agv_distance))[0][0]

        #         print("exe_agv",exe_agv)

        agv_move_unloading(exe_agv, agv_state[exe_agv][1], dp, agv_timestamp,
                           agv_state)  # 将AGV从current position 移动到 demand point

    agv_state[exe_agv][1] = mach_p  # 更改current point

    if (tag[order] == 0):  # 若是job的第1个工序,无需check

        agv_timestamp[exe_agv].append(agv_timestamp[exe_agv][-1])  # 将waiting_time=0添加至agv_timestamp中

        agv_move_loading(exe_agv, -1, mach_p, agv_timestamp, agv_state)  # 将AGV从 demand point(L/U) 移动到 request point

    else:

        if (machine_timestamp[dp][-1] <= agv_timestamp[exe_agv][-1]):  # check whether the job is ready--yes

            agv_timestamp[exe_agv].append(agv_timestamp[exe_agv][-1])  # 将waiting_time=0添加至agv_timestamp中

            agv_move_loading(exe_agv, dp, mach_p, agv_timestamp, agv_state)  # 将AGV从 demand point 移动到 request point

        else:  # check whether the job is ready--no

            agv_move_waiting(exe_agv, dp, agv_timestamp, machine_timestamp, agv_state)

            agv_move_loading(exe_agv, dp, mach_p, agv_timestamp, agv_state)  # 将AGV从 demand point 移动到 request point

    # load the job on the machine

    job_loading(exe_agv, mach_p, time_p, agv_timestamp, machine_timestamp)

    # 计算优化函数1/2：makespan_list 记录各个job完工的时间，delay_list 记录各个job的延迟时间

    if (tag[order] == operation_num - 1):  # 如果这是job的最后一个工序，记录下完工所需时间

        makespan_list[order] = machine_timestamp[mach_p][-1] - makespan_list[order]

        delay_list[order] = abs(machine_timestamp[mach_p][-1] - delay_info[order])


def agv_move_unloading(exe_agv, start, end, agv_timestamp, agv_state):
    agv_state[exe_agv][0] = "unloading"  # agv状态为 unloading 运行中

    agv_timestamp[exe_agv].append(time_compute(start, end) + agv_timestamp[exe_agv][-1])  # 将时间刻度添加至agv_timestamp中


def agv_move_loading(exe_agv, start, end, agv_timestamp, agv_state):
    agv_state[exe_agv][0] = "loading"  # agv状态为 loading 运行中

    agv_timestamp[exe_agv].append(time_compute(start, end) + agv_timestamp[exe_agv][-1])  # 将时间刻度添加至agv_timestamp中


def agv_move_waiting(exe_agv, dp, agv_timestamp, machine_timestamp, agv_state):
    agv_state[exe_agv][0] = "waiting"  # agv状态为 waiting 中

    agv_timestamp[exe_agv].append(
        machine_timestamp[dp][-1] - agv_timestamp[exe_agv][-1] + agv_timestamp[exe_agv][-1])  # 将时间刻度添加至agv_timestamp中


def job_loading(exe_agv, mach_p, time_p, agv_timestamp, machine_timestamp):
    if (agv_timestamp[exe_agv][-1] <= machine_timestamp[mach_p][-1]):  # AGV到达时，机器正在加工

        machine_timestamp[mach_p].append(machine_timestamp[mach_p][-1])  # 加载一个值为0的waiting time 方便后续甘特图
        machine_timestamp[mach_p].append(time_p + machine_timestamp[mach_p][-1])  # 将时间刻度添加至machine_timestamp中

    else:

        machine_timestamp[mach_p].append(agv_timestamp[exe_agv][-1])  # 1. 加载AGV到达的时刻
        machine_timestamp[mach_p].append(time_p + machine_timestamp[mach_p][-1])  # 2. 加载机器处理job结束的时刻


#### 三目标函数快速非支配算法 ####

def fast_non_dominated_sort(values1, values2, values3):
    S = [[] for i in range(0, len(values1))]  # S记录被P支配的集合
    n = [0 for i in range(0, len(values1))]  # n记录p被支配的个数
    front = [[]]  # 返回前沿

    #     rank = [0 for i in range(0, len(values1))]  # 生成values1(种群大小)个0的一维列表,rank指非支配序

    for p in range(0, len(values1)):  # p为种群中的每个个体,遍历种群
        S[p] = []  # 初始化当前个体p支配的集合
        n[p] = 0  # 初始化当前个体p被支配的个数
        for q in range(0, len(values1)):  # q为种群中的每个个体,遍历种群(两层遍历达到每两两都进行比较)
            if values1[p] == values1[q] and values2[p] == values2[q] and values3[p] == values3[q]:
                continue
            else:
                if values1[p] <= values1[q] and values2[p] <= values2[q] and values3[p] <= values3[q]:  # 证明p支配q
                    if q not in S[p]:
                        S[p].append(q)  # 添加被p支配的集合S[p]
                elif values1[q] <= values1[p] and values2[q] <= values2[p] and values3[q] <= values3[p]:  # 证明p被q支配
                    n[p] = n[p] + 1  # 使其p被支配个数n[p]加1
        if n[p] == 0:  # 如果p不被其他群体支配
            front[0].append(p)

    i = 0
    while (front[i] != []):  # 通过判断前一列序列是否为空(保证不为空) i = 0(即第二步)
        Q = []  # 初始化另一个集Q
        for p in front[i]:  # 遍历当前列序列里面的每个个体
            for q in S[p]:  # 考察它所支配的个体集s[p]
                n[q] = n[q] - 1  # 集合S[p]中的每个个体q的(被支配个数)n[q]减1
                if n[q] == 0:
                    Q.append(q)
        i = i + 1  # 继续下一序列集
        front.append(Q)  # 将其添加到前沿序列

    del front[len(front) - 1]  # 最后一个序列什么都没有添加元素,即(front[i]==0)循环结束,故删除最后一个空列表
    return front  # 返回快速非支配排序后的结果,rank为其所在二维列表中所在一维列表的下标


#### 三目标函数的拥挤度排序 （矩阵变换）####
##[0,1,2,3] front
##[0,0,0,0] crowding distance
##[ , , , ] values 1
##[ , , , ] values 2
##[ , , , ] values 3

def crowding_distance(values1, values2, values3, front):  # 注意:此处front是一维列表形式,每个层次序号

    l = len(front)

    distance = [0 for i in range(0, l)]  # 初始化该层个体距离集合,且初值为0

    values1_p = []
    values2_p = []
    values3_p = []

    for i in range(l):
        values1_p.append(values1[i])
        values2_p.append(values2[i])
        values3_p.append(values3[i])

    matrix = np.array([front, distance, values1_p, values2_p, values3_p])

    ## 按照函数1进行拥挤度赋值(选择排序)

    # 排序

    for m in range(l - 1):
        tag = np.argmin(matrix[2, m:]) + m

        matrix[:, [m, tag]] = matrix[:, [tag, m]]

    #  赋值

    matrix[1, 0] = matrix[1, l - 1] = 444444444  # 首尾赋值无穷

    for k in range(1, l - 1):
        matrix[1, k] = (matrix[2, k + 1] - matrix[2, k - 1]) / (max(values1) - min(values1)) + matrix[1, k]

    ## 按照函数2进行拥挤度赋值(选择排序)

    # 排序

    for m in range(l - 1):
        tag = np.argmin(matrix[3, m:]) + m

        matrix[:, [m, tag]] = matrix[:, [tag, m]]

    #  赋值

    matrix[1, 0] = matrix[1, l - 1] = 444444444  # 首尾赋值无穷

    for k in range(1, l - 1):
        matrix[1, k] = (matrix[3, k + 1] - matrix[3, k - 1]) / (max(values2) - min(values2)) + matrix[1, k]

    ## 按照函数3进行拥挤度赋值(选择排序)

    # 排序

    for m in range(l - 1):
        tag = np.argmin(matrix[4, m:]) + m

        matrix[:, [m, tag]] = matrix[:, [tag, m]]

    #  赋值

    matrix[1, 0] = matrix[1, l - 1] = 444444444  # 首尾赋值无穷

    for k in range(1, l - 1):
        matrix[1, k] = (matrix[4, k + 1] - matrix[4, k - 1]) / (max(values3) - min(values3)) + matrix[1, k]

    # 输出按照拥挤度由大到小排序的sorted_final

    for m in range(l - 1):
        tag = np.argmax(matrix[1, m:]) + m

        matrix[:, [m, tag]] = matrix[:, [tag, m]]

    return list(map(int, matrix[0, :]))  # 返回改层拥挤距离集 ，如[68，154，58，110]
####用甘特图显示AGV状态####

def plot_fig_agv(max_t, ts_list,k):
    M = len(ts_list)
    fig, ax = plt.subplots()

    ax.set_ylim(5, 15+10*M)
    ax.set_xlim(0, max_t + 30)

    label_list = []
    y_loc = []
    for m in range(M):
        label_list.append('AGV-'+str(m))
        y_loc.append(15+m*10)

    plt.yticks(y_loc, label_list)

    for m in range(M):
        barh_list = []
        color_list = []
        if len(ts_list[m])!=0:
            for i in range(len(ts_list[m])):
                if i ==0:
                    barh_list.append((0, ts_list[m][i]))
                else:
                    barh_list.append((ts_list[m][i-1], ts_list[m][i]-ts_list[m][i-1]))
                if i%3==0:
                    color_list.append('green') # loading
                elif i%3==1:
                    color_list.append('red') # unloading
                elif i%3==2:
                    color_list.append('white') # waiting
        else:
            barh_list = [(0, 0)]
            color_list = 'red'
            # 设置条形数据
        ax.broken_barh(barh_list, (10 + 10*m, 6), facecolor=color_list)
    path = 'C:/Users/15820/Desktop/' + str(k) + '-agv.jpg'
    plt.savefig(path, dpi=750, bbox_inches='tight')
    plt.show()
# plot_fig_agv(max_agv_t, agv_timestamp)


# In[10]:


####用甘特图显示机器状态####
from matplotlib import pyplot as plt

def plot_fig_machine(max_t, ts_list,k):
    M = len(ts_list)
    fig, ax = plt.subplots()

    ax.set_ylim(5, 15+10*M)
    ax.set_xlim(0, max_t + 30)

    label_list = []
    y_loc = []
    for m in range(M):
        label_list.append('machine-'+str(m))
        y_loc.append(15+m*10)

    plt.yticks(y_loc, label_list)

    for m in range(M):
        barh_list = []
        color_list = []
        if len(ts_list[m])!=0:
            for i in range(len(ts_list[m])):
                if i ==0:
                    barh_list.append((0, ts_list[m][i]))
                else:
                    barh_list.append((ts_list[m][i-1], ts_list[m][i]-ts_list[m][i-1]))
                if i%2==0:
                    color_list.append('black') # working
                else :
                    color_list.append('white') # waiting
        else:
            barh_list = [(0, 0)]
            color_list = 'red'
            # 设置条形数据
        ax.broken_barh(barh_list, (10 + 10*m, 6), facecolor=color_list)
    path = 'C:/Users/15820/Desktop/' + str(k) + '-machine.jpg'
    plt.savefig(path, dpi=750, bbox_inches='tight')
    plt.show()

###### 正式优化 ######
import time

start = time.perf_counter()

N_size = 80  # 种群数量 80
max_gen = 100  # 迭代数量 500
mutation_rate = 0.1 # 变异概率
count = 0

pop = init(N_size)  # 构建初始种群 array [80x36]

# print(pop)

gen_n = 0

while (gen_n < max_gen):

    show = 0
    #### 交叉，变异 ####
    pop2 = pop[:]

    while (len(pop2) != 2 * N_size):  # 使产生的种群大小为2N
        m = random.randint(0, N_size - 1)
        n = random.randint(0, N_size - 1)
        child_1,child_2 = POX_improve(pop[m],pop[n]) # 交叉算子
        if random.random() < mutation_rate: # 变异算子
            child_1 = mut_job(child_1)
        if random.random() < mutation_rate:
            child_2 = mut_job(child_2)

        pop2 = np.r_[pop2, [child_1]]
        pop2 = np.r_[pop2, [child_2]]
    #### 精英选择 ####
    function1_values = [0 for i in range(2 * N_size)]  # 记录种群function1_values
    function2_values = [0 for i in range(2 * N_size)]  # 记录种群function2_values
    function3_values = [0 for i in range(2 * N_size)]  # 记录种群function2_values

    for i in range(pop2.shape[0]):
        function1_values[i], function2_values[i], function3_values[i] = single_scheduling_time(pop2[i],show,count)

    non_dominated_sorted_solution = fast_non_dominated_sort(function1_values, function2_values, [-1*i for i in function3_values])

    crowding_distance_solution = []  # 记录按照拥挤度改层之后的索引

    for i in range(0, len(non_dominated_sorted_solution)):
        crowding_distance_solution.append(
            crowding_distance(function1_values, function2_values, function3_values, non_dominated_sorted_solution[i]))

    new_solution = []  # 初始化新解集

    for i in range(0, len(crowding_distance_solution)):  # 遍历拥挤距离集层

        l = len(new_solution)

        if len(crowding_distance_solution[i]) + l <= N_size:

            for k in crowding_distance_solution[i]:
                new_solution.append(k)

        else:

            for k in crowding_distance_solution[i]:

                new_solution.append(k)

                if (len(new_solution) == N_size):
                    break

            break

    for i in range(N_size):
        pop[i] = pop2[new_solution[i]]

    #### 如果这是最后一轮迭代 则输出结果 ####

    if gen_n == max_gen - 1:

        function1_values = [0 for i in range(N_size)]
        function2_values = [0 for i in range(N_size)]
        function3_values = [0 for i in range(N_size)]

        for i in range(N_size):
            function1_values[i], function2_values[i],function3_values[i] = single_scheduling_time(pop[i],show,count)

        non_dominated_sorted_solution = fast_non_dominated_sort(function1_values, function2_values,[-1*i for i in function3_values])

        # 输出结果
        result = []

        show = 1

        for valuez in non_dominated_sorted_solution[0]:  # 遍历最优解集合

            result.append(pop[valuez])

        result_final = list(set([tuple(t) for t in result]))

        i = 1

        output1=[]
        output2=[]
        output3=[]

        for pre in result_final:

            pre = list(pre)

            count = count + 1

            fun1, fun2, fun3 = single_scheduling_time(np.array(pre),show,count)
            if((fun1 in output1) or (fun2 in output2) or (fun3 in output3)):
                continue
            output1.append(fun1)
            output2.append(fun2)
            output3.append(fun3)

            print("第", i, "个最优的解:(", fun1, ",", fun2, ",", fun3,")")

            i = i + 1

            if(i>5):

                break

    gen_n += 1

# for i in range(pop.shape[0]):
#     print(pop[i])
print(output1)
print(output2)
print(output3)
from mpl_toolkits.mplot3d import Axes3D

function1 = [i for i in function1_values]
function2 = [j for j in function2_values]
function3 = [k for k in function3_values]
# print("function2",function2)


# 绘制散点图
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(function1, function2, function3)

# 添加坐标轴(顺序是Z, Y, X)
ax.set_zlabel('profit', fontdict={'size': 15, 'color': 'red'})
ax.set_ylabel('Tardiness', fontdict={'size': 15, 'color': 'red'})
ax.set_xlabel('Makespan', fontdict={'size': 15, 'color': 'red'})
plt.show()

# plt.xlabel('Makespan', fontsize=10)
# plt.ylabel('Tardiness', fontsize=10)
# plt.clabel('profit', fontsize=10)
# plt.scatter(function1, function2,function3)
# plt.show()