论文复现|An Iterated Greedy Heuristic for Mixed No-Wait Flowshop Problems|流水车间调度，迭代贪婪IG算法

这篇文章是IEEE TRANSACTIONS ON CYBERNETICS（TCYB）顶刊上的文章，主要是使用改进的IG求解MIXED零等待问题，也就是部分机器需要零等待，部分不需要这么个问题，伪代码有一小部分没看懂，就没有实现，其他的大体框架都实现完成，可以直接运行。有问题欢迎留言或者联系某站  lvqaq的数学建模

import numpy as np
import random
from matplotlib import pyplot as plt

def decode(p, maxLen, pi, Mw, Mv, d1, d2):
    """
    SMC Algorithm
    p 为加工时间矩阵，其中需要包含虚拟工件
    maxLen 为工件数
    pi 为调度序列
    Mw 为FSP约束的机器号
    Mv 为NWFSP约束的机器组的列表
    d1，d2 为文中的d1 d2
    """
    pi = [0] + pi
    maxLen += 1
    n = len(pi)
    m = p.shape[1]  # 机器数
    # print(f'机器数{m}')
    maxV = len(Mv)      # NW机器组的块数
    C = np.zeros((maxLen, m+1))
    S = np.zeros((maxLen, m+1))
    for k in range(1, n): # 多一个虚拟工件
        v = 0       # 这里伪码中写的是1，v表示第几块NW机器组
        j = 1       # 这里伪码写的是1，j表示机器的下角标，按理解机器号应该从1开始，下角标从0开始
        while j <= m:
            # print(j, Mw, m)
            if j in Mw:
                S[k,j] = max(C[k,j-1], C[k-1,j])
                C[k,j] = S[k,j]+p[pi[k],j-1]
                j += 1
                # print(f'j in Mw:{j}')
                continue
            if j == Mv[v][0]:
                # print(f'j in Mv:{j}')
                if C[k,j-1] - C[k-1,j] <= d1[v][pi[k-1]][pi[k]]:
                    a = 0
                else:
                    a = C[k, j-1] - C[k-1, j] - d1[v][pi[k-1],pi[k]]
                # print(C)
                C[k, j] = C[k-1, j] + d1[v][pi[k-1],pi[k]] + p[pi[k],j-1] + a
                lv = len(Mv[v])
                # print(d1[v][pi[k-1],pi[k]], p[pi[k],j-1], a)
                # print(C)
                # print(S)
                # print('\n\n')
                # print(k, j)
                # print(C[k-1,j+lv-1])
                C[k, j+lv-1] = C[k-1,j+lv-1] + d2[v][pi[k-1],pi[k]] + a
                j += lv
                if v+1 < maxV:
                    v += 1
                # print(f'j in Mv:{j}')
    return C.max()


def NEH(k, p, Mv, Mw, d1, d2):
    num_jobs, num_machines = p.shape
    num_jobs -= 1       # 考虑NEH时，不考虑虚拟工件
    # print(num_jobs)
    job_indices = np.argsort(-np.sum(p, axis=1))
    print(job_indices)
    sequence = [job_indices[0]]

    for i in range(1, num_jobs):
        best_sequence = None
        best_makespan = float('inf')

        for j in range(len(sequence) + 1):
            new_sequence = sequence[:j] + [job_indices[i]] + sequence[j:]
            # print(new_sequence, j, len(new_sequence))
            makespan = decode(p, num_jobs, new_sequence, Mw, Mv, d1, d2)

            if makespan < best_makespan:
                best_makespan = makespan
                best_sequence = new_sequence

        sequence = best_sequence

    return sequence, best_makespan


def cal_d2(p, Mv):
    n = p.shape[0]      # 工件数 包括虚拟工件
    Mv = [[element - 1 for element in sublist] for sublist in Mv]
    d2 = []
    # arr = np.zeros((n, n))
    for subMv in Mv:
        # print(subMv)
        arr = np.zeros((n, n))
        for i in range(n):
            for j in range(1, n):
                if i == j:
                    continue
                maxTemp = []
                for h in subMv:
                    sum = 0

                    for t in range(h, subMv[-1]+1):
                        sum += p[j, t] - p[i, t]
                    else:
                        # print(sum)
                        sum += p[i, h]
                    maxTemp.append(sum)
                    # print(f'i:{i} j:{j} h:{h} sum:{sum}')
                arr[i, j] = max(maxTemp)
        d2.append(arr)
    return d2


def cal_d1(d2, p, Mv):
    n = p.shape[0]  # 工件数 包括虚拟工件
    Mv = [[element - 1 for element in sublist] for sublist in Mv]
    d1 = []
    # arr = np.zeros((n, n))
    for index, subMv in enumerate(Mv):
        # print(subMv)
        arr = np.zeros((n, n))
        for i in range(n):
            for j in range(1, n):
                if i == j:
                    continue
                arr[i, j] = d2[index][i, j]
                sum1 = 0
                for h in range(subMv[1], subMv[-1] + 1):
                    sum1 += p[i, h]
                sum2 = 0
                for h in range(subMv[0], subMv[-1] + 1):
                    sum2 += p[j, h]
                arr[i, j] = arr[i, j] + sum1 - sum2
        d1.append(arr)
    return d1


def Insert(L, position, x):
    """
    把x插入到position的位置
    其中position需要满足  0<= position <= len(L)
    """
    L = L.copy()
    if not isinstance(x, list):
        x = [x]
    if position == 0:
        return x + L
    if position == len(L):
        return L + x
    return L[:position] + x + L[position:]


def Swap(L, position1, position2):
    L = L.copy()
    L[position1], L[position2] = L[position2], L[position1]
    return L


def MDR(pi, r, p, Mv, Mw, d1, d2):
    pi = pi.copy()
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    B = []
    r = min(r, num_jobs-1)
    for _ in range(r):
        index = random.randint(0, len(pi)-1)
        # print(f'index {index}')
        B.append(pi[index])
        pi.remove(B[-1])
    # print(f'B:{B}')
    # print(pi)
    for i in B:
        best_pi = None
        best_makespan = float('inf')
        for index in range(len(pi)+1):
            new_pi = Insert(pi, index, i)       # 把B[i]插入pi中
            makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
            if makespan < best_makespan:
                best_makespan = makespan
                best_pi = new_pi
        else:
            pi = best_pi
    # print(pi)
    return pi, best_makespan


def CINS(pi, p, Mv, Mw, d1, d2):
    """
    CINS型领域结构
    返回领域结果，同时返回是否更改
    """
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    makespan_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
    pi_s = pi.copy()
    pi_t = pi.copy()
    for job in pi_s:
        pi_t.remove(job)
        best_pi = None
        best_makespan = float('inf')
        for index in range(num_jobs+1):
            new_pi = Insert(pi_t, index, job)
            makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
            if makespan < best_makespan:
                best_makespan = makespan
                best_pi = new_pi
        else:
            pi_t = pi.copy()
    else:
        if best_makespan < makespan_star:
            pi_star = best_pi
            flag = True
        else:
            flag = False
    return pi_star, flag


def GINS(pi, p, Mv, Mw, d1, d2):
    """
    GINS 型领域结构
    """
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    makespan_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
    pi_s = pi.copy()
    pi_t = pi.copy()
    random.shuffle(pi_s)
    for job in pi_s:
        pi_t.remove(job)
        best_pi = None
        best_makespan = float('inf')
        for index in range(num_jobs+1):
            new_pi = Insert(pi_t, index, job)
            makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
            if makespan < best_makespan:
                best_makespan = makespan
                best_pi = new_pi
        else:
            pi_t = best_pi
    else:
        if decode(p, num_jobs, pi_t, Mw, Mv, d1, d2) < makespan_star:
            pi_star = pi_t
            flag = True
        else:
            flag = False
    return pi_star, flag


def CPINS(pi, p, Mv, Mw, d1, d2):
    """
    CPINS型领域结构
    返回领域结果，同时返回是否更改
    """
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    makespan_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
    pi_s = pi.copy()
    for job_index in range(len(pi_s)-1):
        pi_t = pi.copy()
        adjacent_job = pi_s[job_index:job_index+2].copy()
        pi_t.remove(adjacent_job[0])
        pi_t.remove(adjacent_job[1])
        best_pi = None
        best_makespan = float('inf')
        for index in range(num_jobs-1):
            new_pi = Insert(pi_t, index, adjacent_job)
            makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
            if makespan < best_makespan:
                best_makespan = makespan
                best_pi = new_pi

    else:
        if best_makespan < makespan_star:
            pi_star = best_pi
            flag = True
        else:
            flag = False
    return pi_star, flag


def CPSNS(pi, p, Mv, Mw, d1, d2):
    """
   CPSNS 型领域结构
   返回领域结果，同时返回是否更改
   """
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    makespan_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
    pi_t = pi.copy()
    best_pi = None
    best_makespan = decode(p, num_jobs, pi_t, Mw, Mv, d1, d2)
    for i in range(num_jobs-1):
        new_pi = Swap(pi_t, i, i+1)
        makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
        if makespan < best_makespan:
            best_makespan = makespan
            best_pi = new_pi
    if best_makespan < makespan_star:
        pi_star = best_pi
        flag = True
    else:
        flag = False
    return pi_star, flag


def GPINS(pi, p, Mv, Mw, d1, d2):
    """
    GPINS型领域结构
    返回领域结果，同时返回是否更改
    """
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    makespan_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
    pi_s = pi.copy()
    pi_t = pi.copy()
    for job_index in range(len(pi_s)-1):
        adjacent_job = pi_s[job_index:job_index+2].copy()
        pi_t.remove(adjacent_job[0])
        pi_t.remove(adjacent_job[1])
        best_pi = None
        best_makespan = float('inf')
        for index in range(num_jobs-1):
            new_pi = Insert(pi_t, index, adjacent_job)
            makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
            if makespan < best_makespan:
                best_makespan = makespan
                best_pi = new_pi
        pi_t = new_pi

    else:
        if best_makespan < makespan_star:
            pi_star = best_pi
            flag = True
        else:
            flag = False
    return pi_star, flag


def CSNS(pi, p, Mv, Mw, d1, d2):
    """
    CSNS 型领域结构
    返回领域结果，同时返回是否更改
    """
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    makespan_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
    pi_t = pi.copy()
    best_pi = None
    best_makespan = decode(p, num_jobs, pi_t, Mw, Mv, d1, d2)
    for i in range(num_jobs):
        for j in range(i+1, num_jobs):
            new_pi = Swap(pi_t, i, j)
            makespan = decode(p, num_jobs, new_pi, Mw, Mv, d1, d2)
            if makespan < best_makespan:
                best_makespan = makespan
                best_pi = new_pi
    if best_makespan < makespan_star:
        pi_star = best_pi
        flag = True
    else:
        flag = False
    return pi_star, flag


def select_neighborhood(k, pi, p, Mv, Mw, d1, d2):
    new_pi = pi.copy()
    if k == 1:
        new_pi, flag0 = CINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag1 = GINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag2 = CSNS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag3 = CINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag4 = CSNS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag5 = CSNS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag6 = CSNS(new_pi, p, Mv, Mw, d1, d2)
        flag = any([flag0, flag1, flag2, flag3, flag4, flag5, flag6])
    elif k == 2:
        new_pi, flag0 = CINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag1 = CSNS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag2 = CINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag3 = GINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag4 = GINS(new_pi, p, Mv, Mw, d1, d2)
        flag = any([flag0, flag1, flag2, flag3, flag4])
    elif k == 3:
        new_pi, flag0 = CPINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag1 = GPINS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag2 = CPSNS(new_pi, p, Mv, Mw, d1, d2)
        flag = any([flag0, flag1, flag2])
    elif k == 4:
        new_pi, flag0 = CPSNS(new_pi, p, Mv, Mw, d1, d2)
        new_pi, flag1 = GPINS(new_pi, p, Mv, Mw, d1, d2)
        flag = any([flag0, flag1])
    else:
        raise ValueError
    return new_pi.copy(), flag     # 若flag为True表明更新了最优解


def VND(pi, p, Mv, Mw, d1, d2):
    pi = pi.copy()
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    pi_star = pi.copy()
    i = 1
    while i <= 4:
        Cmax_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
        pi_t, flag = select_neighborhood(i, pi_star, p, Mv, Mw, d1, d2)
        Cmax_t = decode(p, num_jobs, pi_t, Mw, Mv, d1, d2)
        if Cmax_t < Cmax_star:
            pi_star = pi_t.copy()
            i = 1
        else:
            i += 1
    return pi_star


def MIG(r0, p, Mv, Mw, d1, d2):
    delta_r = 5
    num_jobs, num_machines = p.shape
    num_jobs -= 1
    r = r0
    x = 0
    pi, NEH_makespan = NEH(None, p, Mv, Mw, d1, d2)
    print(f'NEH makespan is {NEH_makespan}')

    # pi = list(range(1, 16))
    # random.shuffle(pi)
    # random_makespan = decode(p, 15, pi, Mw, Mv, d1, d2)
    # print(f'random makespan is {random_makespan}')

    pi, _ = CINS(pi, p, Mv, Mw, d1, d2)
    pi_star = pi.copy()
    g, gen = 0, 50
    history = []
    while g < gen:
        print(g, r, x)
        g += 1
        pi_, _ = MDR(pi_star, r, p, Mv, Mw, d1, d2)
        pi__ = VND(pi_, p, Mv, Mw, d1, d2)
        Cmax_pi__ = decode(p, num_jobs, pi__, Mw, Mv, d1, d2)
        Cmax_pi = decode(p, num_jobs, pi, Mw, Mv, d1, d2)
        Cmax_pi_star = decode(p, num_jobs, pi_star, Mw, Mv, d1, d2)
        if  Cmax_pi__ < Cmax_pi:
            pi = pi__.copy()
            x = 0
            r = r0
            if Cmax_pi__ < Cmax_pi_star:
                pi_star = pi
        else:
            alpha = np.random.rand()
            if alpha < np.exp(-((Cmax_pi__ - Cmax_pi)/Cmax_pi)*100):
                pi = pi__.copy()
            x += 1
            if x > 10:
                r = min(r+1, r0+delta_r)
        history.append(Cmax_pi_star)
    plt.plot(history)
    print(Cmax_pi_star)
    plt.show()
    return Cmax_pi_star





# p = np.array([[3,5,6],[6,3,2],[1,3,3],[4,2,4]])
# # print(p.shape)
# p = np.array([[0,0,0,0,0], [3,6,1,4,1],[5,3,3,2,1],[6,2,3,4,1]])
# print(p)
# pi = [1,2,3]
# Mw = [1,5]
# Mv = [[2,3,4]]
# # m = 4
# # n = 4
# d2 = cal_d2(p, Mv)
# d1 = cal_d1(d2, p, Mv)
# C = decode(p,3,[2,1,3],Mw,Mv,d1,d2)
# # print(decode(p,len(pi),[1,3,2],Mw,Mv,d1,d2))
# print(C)
# print(MDR(pi, 1, p, Mv, Mw, d1, d2))

# print(VND([1,2,3], p, Mv, Mw, d1, d2))

data = np.array([
    [63, 540, 420, 216, 963, 144, 143, 30, 196, 38],
    [28, 80, 360, 112, 428, 52, 55, 50, 784, 31],
    [84, 420, 480, 133, 294, 112, 205, 55, 785, 130],
    [60, 20, 0, 3, 60, 12, 10, 23, 83, 11],
    [60, 90, 180, 78, 105, 27, 15, 59, 57, 29],
    [240, 80, 420, 40, 96, 88, 26, 24, 97, 84],
    [60, 80, 180, 100, 132, 68, 44, 27, 136, 136],
    [30, 90, 120, 87, 57, 66, 200, 34, 123, 31],
    [105, 420, 480, 49, 98, 119, 110, 46, 80, 42],
    [32, 120, 300, 48, 244, 32, 96, 14, 84, 37],
    [120, 120, 60, 14, 98, 52, 13, 19, 130, 21],
    [200, 300, 480, 60, 220, 90, 144, 27, 427, 323],
    [90, 270, 660, 180, 954, 90, 363, 54, 850, 38],
    [105, 420, 480, 91, 434, 84, 128, 23, 62, 95],
    [60, 120, 540, 48, 600, 120, 315, 24, 424, 20]
])
zero_row = np.zeros((1, data.shape[1]), dtype=data.dtype)  # 创建一行全0的行
p = np.vstack((zero_row, data))  # 将全0行与原数组垂直堆叠
# print(new_data)
Mv = [[4,5,6,7,8]]
Mw = [1,2,3,9,10]



d2 = cal_d2(p, Mv)
d1 = cal_d1(d2, p, Mv)
print('ready')

# print(decode(p, 15,
#              [2,4,11,15,3,13,14,7,8,12,1,6,9,10,5],
#              Mw, Mv, d1, d2))
#
print(MIG(1, p, Mv, Mw, d1, d2))