遗传算法(GA)求解带约束的流水车间(flow shop)调度问题——python

算法分析与设计课程的期末大作业，问题如下：

Problem 2: job scheduling problem.

The production line has 8 machines, and 28 jobs are ready to be processed. The process time of each job to the machine is given as follows, the number in the table is the process time of each job on each machine:

	machine1	machine2	machine3	machine4	machine5	machine6	machine7	machine8
Job1	11	16	15	4	15	22	17	22
Job2	9	9	15	3	21	4	16	2
Job3	15	15	19	15	22	15	17	8
Job4	22	33	16	19	21	21	27	9
Job5	14	28	14	12	12	17	8	21
Job6	7	18	19	11	10	23	13	19
Job7	18	19	20	5	29	16	9	10
Job8	8	32	32	16	1	12	8	5
Job9	23	22	14	18	22	11	11	5
Job10	9	17	27	45	12	11	8	6
Job11	21	15	11	10	10	13	17	29
Job12	8	8	12	12	12	16	16	16
Job13	7	7	27	7	7	7	19	9
Job14	11	11	11	12	25	33	11	11
Job15	19	5	5	5	7	8	21	21
Job16	21	12	11	3	3	3	5	5
Job17	7	7	4	4	12	12	16	16
Job18	8	1	1	8	8	1	1	1
Job19	11	22	22	33	33	11	11	22
Job20	11	35	33	18	6	6	6	2
Job21	2	2	33	1	1	51	22	2
Job22	15	15	15	5	5	5	5	5
Job23	19	2	19	2	2	2	19	2
Job24	18	18	38	2	2	22	38	3
Job25	10	9	9	29	1	1	1	22
Job26	50	1	1	5	2	1	44	1
Job27	35	7	7	12	12	35	1	5
Job28	22	1	1	29	11	29	1	1

We suppose:

Each machine can process one job at one time.

The process of each job at each machine cannot be spitted.

Each job has to be processed along machines based on sequence 1,2,3,, ,8  (for example, any job cannot be processed on machine 2 before machine 1, it cannot be processed on machine 3 before machines 1 or 2).

For Jobs 3,5,7,9,11,14 you can choose only 3 job to process, and the others can be ignored.

You need give a process sequence of jobs, such that can minimize the final complete time (at last machine) of last job.

This process sequence cannot be changed along all the machines. For example, once the sequence is 1-2-3-…-27-28, then, for each machine, the job 1 is processed first, then job 2,…, and job28 is the last process job on each machine.

程序设计思路：

一、编码

28个job需要按照确定的时间依次在8台machine上完成8道工序，job间的加工顺序一旦确定在每台machine上都保持一致，同时需要满足job 3,5,7,9,11,14中只有3个被加工，这是一个带有约束的flow shop问题。在进行染色体编码时可以采用如下方式：

每条染色体由31个{1,2,…,28}间的整数排列而成，其中前28个数为head，可以理解为不考虑约束情况下包含28个job的flow shop的加工顺序，后3个数为tail，tail为{3,5,7,9,11,14}中的任意3个不重复的数。则把每条染色体tail部分包含的数从head部分去除将会对应一个原带有约束的flow shop问题的解，如：染色体[8,3,4,…,9,7,6,…,11,28,10,…,3,7,11]表示[8,4,…,9,6,…,28,10…]的加工顺序。

二、初始化种群（initial）

根据种群规模pop_size，将{1,2,…,28}随机打乱pop_size次以获得head，从{3,5,7,9,11,14}中任取3个数，重复操作pop_size次以获得tail，将head和tail拼接，获得初始种群每个个体的染色体。

三、交叉（cross over）

父代种群中的两条染色体（parent1，parent2），以p_cr的概率进行交叉产生子代（child1, child2）。为了保证子代的可行性，可以设计如下交叉方式：

对于head部分，采取单点交叉，即令两条父代染色体的head部分从某一任意索引后的基因片段互换，然后分别遍历两条父代染色体head部分的基因，将没有被包含进子代染色体中的基因依次加入，直到子代染色体的head长度也为28，形成两个子代染色体的head，如下图（以1到7的序列示意）：

对于tail部分，任取1~3间的数i，若parent1的tail中没有parent2的tail部分第i位对应的基因(gene2)，则将parent1 tail部分的第i个基因(gene1)替换位gene2，反之保持不变；同理，若parent2的tail部分没有gene1，则进行替换，反之不替换，如：

四、变异（mutate）

对每个交叉产生的子代染色体child，以p_mu的概率发生变异，变异规则为：对于head部分，任意生成两个1到28间的索引值，使索引值区间内的基因片段逆序排列，对于tail部分，使其内的任意一位基因替换为{3,5,7,9,11,14}中不被包含在原来的tail内的数，如（head以1到7的序列为例）：

五、适应度计算（fitness）利用动态规划的思想可以求解每个解对应的加工顺序的完工时间，令seq表示某染色体代表的加工顺序，dij,i=1,2…8, j=1,2,…,28表示第 j 个job在第 i 台machine上的加工时间，Tij,  i=1,2,…,8, j=1,2,…,25表示第i个machine完成seq中第j个job的时间，则有：

T8,25为最大完工时间，将最大完工时间的倒数作为染色体的适应度。

六、种群选择（selection）

采用精英保留法与轮盘赌结合的方式筛选种群，从每一代交叉变异后的种群中选出一定数量的适应度最高的个体直接进入下一代，接着采用轮盘赌的方式从未被选中的个体中选择，直到下一代种群规模达到预设值，适应度越高的个体被选中的几率越大。

七、进化（evolution）

进化过程如下：

八、求解

设定种群规模为80，交叉概率为0.8，变异概率0.01，保留精英解数目为20，进化500代，记录第一代中完工时间的最小值为546，最后一代中完工时间的最小值为473，缩短了13.4%，最优加工顺序为{21, 15, 6, 12, 17, 1, 19, 24, 5, 18, 2, 28, 23, 14, 13, 25, 3, 4, 26, 8, 27, 10, 16, 20, 22}

两种加工顺序的甘特图如下：

可以发现，优化后的甘特图排程明显更加紧密，空置时间较少。

python代码如下，水平有限大家多多交流哇：

from datetime import time
import pandas as pd
import numpy as np
import random
import heapq
import matplotlib.pyplot as plt

class GA_TSP(object):

    pop_size = 50    # 种群规模
    p_cr = 0.8       # 交换概率
    p_mu = 0.01      # 变异概率
    constraints=[3,5,7,9,11,14]    # 约束条件
    job_num=28                     # job数量
    machine_num=8                  # 机器数量
    elite_number=20                # 每一代中精英解的数量
    minimize_complete_time=[]      # 最小完工时间
    best_solve=[]                  # 最佳作业排序
    work_time = np.array(pd.read_csv('工作时间.csv', header=None, engine='python'))

    def __init__(self,popsize,p_cr,p_mu,iteration,elite_number):
        self.pop_size=popsize
        self.p_cr=p_cr
        self.p_mu=p_mu
        self.iteration=iteration
        self.elite_number=elite_number

    def initial(self):
        pop=[]

        for i in range(self.pop_size):
            gene1=np.arange(1,self.job_num+1)
            np.random.shuffle(gene1)
            gene1=gene1.tolist()
            gene2=random.sample(self.constraints,3)
            gene1.extend(gene2)
            pop.append(gene1)

        return np.array(pop)        #初始化种群

    def cross_over(self,parent1,parent2):
        child1=[]
        child2=[]

        if np.random.rand()>self.p_cr:
            return np.zeros(self.job_num+3),np.zeros(self.job_num+3)
        else:
            index1=np.random.randint(0,self.job_num)
            index2=np.random.randint(self.job_num,self.job_num+3)
            for i in range(index1,self.job_num):
                child1.append(parent2[i])
            for i in range(self.job_num):
                if parent1[i] not in child1:
                    child1.append(parent1[i])
            for i in range(self.job_num,self.job_num+3):
                child1.append(parent1[i])
            if parent2[index2] not in child1[-3:]:
                child1[index2]=parent2[index2]

            for i in range(index1,self.job_num):
                child2.append(parent1[i])
            for i in range(self.job_num):
                if parent2[i] not in child2:
                    child2.append(parent2[i])
            for i in range(self.job_num,self.job_num+3):
                child2.append(parent2[i])
            if parent1[index2] not in child2[-3:]:
                child2[index2]=parent1[index2]

            return np.array(child1),np.array(child2)    #交叉

    def mutate(self,gene):

        if np.random.rand()>self.p_mu:
            return gene
        else:
            index1=np.random.randint(0,self.job_num)
            index2=np.random.randint(index1,self.job_num)
            index3=np.random.randint(self.job_num,self.job_num+3)
            gene_slice=gene[index1:index2].copy()
            gene_slice_reverse=gene_slice[::-1]

            for i in range(len(gene_slice)):
                gene[index1+i]=gene_slice_reverse[i]
            constraints_rem=self.constraints.copy()

            for i in gene[-3:]:
                constraints_rem.remove(i)
            gene_insert=random.choice(constraints_rem)
            gene[index3]=gene_insert

            return gene                                #变异

    def fitness(self,pop):
        c_time=[]
        fitness_all=[]

        for i in range(pop.shape[0]):
            pop_cal=pop[i].copy().tolist()
            pop_cal_head=pop_cal[:self.job_num].copy()
            pop_cal_tail=pop_cal[-3:].copy()

            for j in pop_cal_tail:
                pop_cal_head.remove(j)
            ctime_dp=np.zeros((self.machine_num,self.job_num-3))

            for ii in range(self.job_num-3):
                for iii in range(ii+1):
                    ctime_dp[0][ii]=ctime_dp[0][ii]+self.work_time[pop_cal_head[iii]-1][0]
            for ii in range(1,self.machine_num):
                for iii in range(ii+1):
                    ctime_dp[ii][0]=ctime_dp[ii][0]+self.work_time[pop_cal_head[0]-1][iii]
            for ii in range(1,self.job_num-3):
                for iii in range(1,self.machine_num):
                    ctime_dp[iii][ii]=max(ctime_dp[iii-1][ii],ctime_dp[iii][ii-1])+\
                                      self.work_time[pop_cal_head[ii]-1][iii]

            c_time.append(ctime_dp[self.machine_num-1][self.job_num-4])
            fitness_all.append(1/ctime_dp[self.machine_num-1][self.job_num-4])
        index=c_time.index(min(c_time))
        best_solve_iter=pop[index].copy()

        return min(c_time),fitness_all,best_solve_iter      #动态规划求解最大完工时间，计算适应度

    def selection(self,pop,fitness):
        pop_after_select=[]
        pop_before_select=pop.copy().tolist()
        index1=list(map(fitness.index,heapq.nlargest(self.elite_number,fitness)))

        for i in index1:
            pop_after_select.append(pop_before_select[i])
        pop_before_select=[pop_before_select[i] for i in range(len(pop_before_select)) \
                           if i not in index1]
        fitness=[fitness[i] for i in range(len(fitness)) \
                           if i not in index1]
        probability=np.array(fitness)/np.array(fitness).sum()
        index2=np.random.choice(np.arange(len(fitness)),size=self.pop_size-self.elite_number,\
                                replace=False,p=probability)

        for i in index2:
            pop_after_select.append(pop_before_select[i])

        return np.array(pop_after_select)                  #筛选种群：精英保留法+轮盘赌

    def evolution(self):
        population=self.initial()

        for i in range(self.iteration):
            for ii in range(self.pop_size):
                iii=np.random.randint(0,self.pop_size-1)
                if ii != iii:
                    child1,child2=self.cross_over(population[ii],population[iii])
                    if child1[0] != 0:
                        mut_child1=self.mutate(child1)
                        mut_child2=self.mutate(child2)
                        population=np.vstack((population,mut_child1))
                        population=np.vstack((population,mut_child2))

            best_time_iter,fitness_iter,best_solve_iter=self.fitness(population)
            self.minimize_complete_time.append(best_time_iter)
            self.best_solve.append(best_solve_iter)
            population=self.selection(population,fitness_iter)

        return self.minimize_complete_time,self.best_solve   #进化

ga=GA_TSP(80,0.8,0.01,500,20)
solution_time,solution_schedule=ga.evolution()

worst_schedule=list(solution_schedule[0].copy())
worst_schedule_head=worst_schedule[:28]
worst_schedule_tail=worst_schedule[-3:]
for i in worst_schedule_tail:
    worst_schedule_head.remove(i)                           #第一代的排班方案

best_schedule=list(solution_schedule[-1].copy())
best_schedule_head=best_schedule[:28]
best_schedule_tail=best_schedule[-3:]
for i in best_schedule_tail:
    best_schedule_head.remove(i)                            #最后一代的排班方案


def gantt(n,m,seq,wt):

    t = np.zeros(n)
    for i in range(m):
        for j in range(n):
            if j == 0:
                t[j]=t[j]+wt[seq[j]-1][i]
                plt.barh(y=i+1,left=t[j]-wt[seq[j]-1][i],height=1,width=wt[seq[j]-1][i])
            else:
                t[j]=max(t[j],t[j-1])+wt[seq[j]-1][i]
                plt.barh(y=i+1,left=t[j]-wt[seq[j]-1][i],height=1,width=wt[seq[j]-1][i])
            if (wt[seq[j]-1][i] != 0):
                pass

    plt.plot([t[-1], t[-1]], [-0.5, m + 1], color=(0, 0, 0))
    x_tickss = [i for i in range(0, int(t[-1]), int((t[-1] / 4) // 10) * 10)]

    if (t[-1] - x_tickss[-1] < (int((t[-1] / 4) // 10) * 10) / 4):
        x_tickss.pop()
    x_tickss.append(t[-1])
    plt.title('n=%d,m=%d,need_time=%.3f' % (n, m, t[-1]))
    plt.xticks(x_tickss)
    plt.show()

    return (t)                                             #绘制甘特图

gantt(ga.job_num-3,ga.machine_num,worst_schedule_head,ga.work_time)
gantt(ga.job_num-3,ga.machine_num,best_schedule_head,ga.work_time)