Viewing the Jupyter Notebooks from nbviewer is encouraged because GitHub is still not fully integrated with the Jupyter Notebook: http://nbviewer.jupyter.org/github/suyunu/Flow-Shop-Scheduling/blob/master/ga-fss.ipynb
In this project, we tried to solve Flow Shop Scheduling Problem (FSSP) with Genetic Algorithm (GA). Before I start doing anything on the problem, I made a literature survey and found these 2 papers:
These 2 papers have done lots of good optimization tests for the parameters and obtained good results. So, I took pieces’ form both papers:
There are n machines and m jobs. Each job contains exactly n operations. The i-th operation of the job must be executed on the i-th machine. No machine can perform more than one operation simultaneously. For each operation of each job, execution time is specified. Operations within one job must be performed in the specified order. The first operation gets executed on the first machine, then (as the first operation is finished) the second operation on the second machine, and so until the n-th operation. Jobs can be executed in any order, however. Problem definition implies that this job order is exactly the same for each machine. The problem is to determine the optimal such arrangement, i.e. the one with the shortest possible total job execution makespan.
I used a simple permutation representation.
The string "ABCDEF" represents a job sequence where "Job A" is processed first, then "Job B" is processed, and so on. If the string "ABCADE" is generated by genetic operators such as crossover and mutation, this string is not a feasible solution of the flow shop scheduling problem because the job "A" appears twice in the string and the job "F" does not appear. Therefore the string used in the flow shop scheduling problem should be the permutation of given jobs.
In this part I will explain the genetic operators such as crossover and mutation as well as the selection mechanism for the flow shop scheduling problem.
I randomly generated $N_{pop}$ number of solutions for the initial population. Suprisingly after lots of test, $N_{pop} = 3$ gave the best results.
The usual method (for maximization) is to measure relative fitness as the ratio of the value of a given chromosome to the population mean. However, in minimization problems, we need to modify this so that low-valued chromosomes are the "good" ones. One approach is to define fitness as $(v_{max} - v)$, but of course we are unlikely to know what $v_{max}$ is. We could use the largest value found so far as a surrogate for $v_{max}$, but in practice this method was not very good at distinguishing between chromosomes. In essence, if the population consists of many chromosomes whose fitness values are relatively close, the resolution between "good" and "bad" ones is not of a high order. This approach was therefore abandoned, in favour of a simple ranking mechanism. Sorting the chromosomes is a simple task, and updating the list is also straightforward. The selection of parents was then made in accordance with the probability distribution
$ p(k) = \frac{2k}{M(M+1)} $
where $k$ is the $k^{th}$ chromosome in ascending order of fitness (i.e. descending order of makespan). This implies that the median value has a chance of $\frac{1}{M}$ of being selected, while the $M^{th}$ (the fittest) has a chance of $\frac{2}{M+1}$, roughly twice that of the median.
We choose $N_{pop}$ pairs of strings from the current population according to this selection probability. So at the end we will have $N_{pop}$ number of children reproduced from those selected pairs.
According to Murata et al. Two-point crossover gives overall better results than other crossover techniques. Also having a crossover probability $P_c = 1$ gives the best results. So we applied crossover to every parent.
In Two-point crossover, two points are randomly selected for dividing one parent. The jobs outside the selected two points are always inherited from one parent to the child, and the other jobs are placed in the order of their appearance in the other parent.
Again according to Murata et al. Shift change gives overall better results than other mutation techniques. Also having a mutation probability $P_m = 1$ gives the best results. So we applied mutation to every child.
In shift change mutation, a job at one position is removed and put at another position. Then all other jobs shifted accordingly. The two positions are randomly selected.
We applied all the operations to $N_{pop}$ pairs of strings (parents). So we ended up with $N_{pop}$ number of children. As a last step we randomly remove one string from the current population and add the best string in the previous population to the current one. Then we continue our processes with this newly generated population.
Total number of evaluations/generations used as a stopping condition. The total number of generations was specifed as to be inversely proportional to the population size $N_{pop}$. For example, when the total number of evaluations was $10,000$, the total number of generations was specified as $\frac{10,000}{N_{pop}}$
I choose the number of generation to stop the program as $10000$
import numpy as np
import math
import time
import random
import itertools
import queue
import pandas as pd
from IPython.display import display, Markdown
Reading input from document and initializing variables. You can change dataset variable to 1, 2 or 3 to use different input sets.
# Dataset number. 1, 2 or 3
dataset = "2"
if dataset == "1":
optimalObjective = 4534
elif dataset == "2":
optimalObjective = 920
else:
optimalObjective = 1302
filename = "data" + dataset + ".txt"
f = open(filename, 'r')
l = f.readline().split()
# number of jobs
n = int(l[0])
# number of machines
m = int(l[1])
# ith job's processing time at jth machine
cost = []
for i in range(n):
temp = []
for j in range(m):
temp.append(0)
cost.append(temp)
for i in range(n):
line = f.readline().split()
for j in range(int(len(line)/2)):
cost[i][j] = int(line[2*j+1])
f.close()
We described all this functions above, except how do we calculate the objective function, which is actually an essential part of this code.
We want to find out the time when the last machine finishes processing of the last job in the list, which is makespan.
We first defined a variable called time to follow the times of the processed jobs and to calculate when a job will be finished when it enters a machine. Then we put those calculations into a priority queue, which sorts the job-machine pairs according to the completion time. In each iteration we take the next item in the priority queue, which has 3 variables:
Then we check two processes:
def initialization(Npop):
pop = []
for i in range(Npop):
p = list(np.random.permutation(n))
while p in pop:
p = list(np.random.permutation(n))
pop.append(p)
return pop
def calculateObj(sol):
qTime = queue.PriorityQueue()
qMachines = []
for i in range(m):
qMachines.append(queue.Queue())
for i in range(n):
qMachines[0].put(sol[i])
busyMachines = []
for i in range(m):
busyMachines.append(False)
time = 0
job = qMachines[0].get()
qTime.put((time+cost[job][0], 0, job))
busyMachines[0] = True
while True:
time, mach, job = qTime.get()
if job == sol[n-1] and mach == m-1:
break
busyMachines[mach] = False
if not qMachines[mach].empty():
j = qMachines[mach].get()
qTime.put((time+cost[j][mach], mach, j))
busyMachines[mach] = True
if mach < m-1:
if busyMachines[mach+1] == False:
qTime.put((time+cost[job][mach+1], mach+1, job))
busyMachines[mach+1] = True
else:
qMachines[mach+1].put(job)
return time
def selection(pop):
popObj = []
for i in range(len(pop)):
popObj.append([calculateObj(pop[i]), i])
popObj.sort()
distr = []
distrInd = []
for i in range(len(pop)):
distrInd.append(popObj[i][1])
prob = (2*(i+1)) / (len(pop) * (len(pop)+1))
distr.append(prob)
parents = []
for i in range(len(pop)):
parents.append(list(np.random.choice(distrInd, 2, p=distr)))
return parents
def crossover(parents):
pos = list(np.random.permutation(np.arange(n-1)+1)[:2])
if pos[0] > pos[1]:
t = pos[0]
pos[0] = pos[1]
pos[1] = t
child = list(parents[0])
for i in range(pos[0], pos[1]):
child[i] = -1
p = -1
for i in range(pos[0], pos[1]):
while True:
p = p + 1
if parents[1][p] not in child:
child[i] = parents[1][p]
break
return child
def mutation(sol):
pos = list(np.random.permutation(np.arange(n))[:2])
if pos[0] > pos[1]:
t = pos[0]
pos[0] = pos[1]
pos[1] = t
remJob = sol[pos[1]]
for i in range(pos[1], pos[0], -1):
sol[i] = sol[i-1]
sol[pos[0]] = remJob
return sol
def elitistUpdate(oldPop, newPop):
bestSolInd = 0
bestSol = calculateObj(oldPop[0])
for i in range(1, len(oldPop)):
tempObj = calculateObj(oldPop[i])
if tempObj < bestSol:
bestSol = tempObj
bestSolInd = i
rndInd = random.randint(0,len(newPop)-1)
newPop[rndInd] = oldPop[bestSolInd]
return newPop
# Returns best solution's index number, best solution's objective value and average objective value of the given population.
def findBestSolution(pop):
bestObj = calculateObj(pop[0])
avgObj = bestObj
bestInd = 0
for i in range(1, len(pop)):
tObj = calculateObj(pop[i])
avgObj = avgObj + tObj
if tObj < bestObj:
bestObj = tObj
bestInd = i
return bestInd, bestObj, avgObj/len(pop)
I've tried 3 different population sizes for 3 datasets. You can find the results below. Obviously and suprisingly population size = 3 gives the best results.
popSize1 = [3,3,3,3,3,5,5,5,5,5,10,10,10,10,10]
objVal1 = [4534,4534,4534,4534,4534,4534,4534,4534,4534,4534,4534,4534,4534,4534,4534]
avgObj1 = [4929.33,4840.33,5002.00,5154.33,4836.00,5066.80,5356.20,5663.20,4850.20,5050.60,5530.40,5839.50,5671.60,5725.70,5459.10]
gap1 = [0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
cpuTime1 = [24.80,25.36,26.22,26.35,25.44,43.37,43.84,43.72,44.11,43.64,86.21,86.90,82.90,80.21,68.50]
popSize2 = [3,3,3,3,3,5,5,5,5,5,10,10,10,10,10]
objVal2 = [920,929,920,920,937,931,937,937,926,920,937,931,929,938,937]
avgObj2 = [1001.00,1002.33,1055.00,964.67,1028.67,1052.40,1185.20,1074.00,1076.00,1063.60,1177.00,1166.00,1139.80,1182.50,1181.40]
gap2 = [0.00,0.98,0.00,0.00,1.85,1.20,1.85,1.85,0.65,0.00,1.85,1.20,0.98,1.96,1.85]
cpuTime2 = [41.74,42.15, 45.23,40.75,39.75,64.36,64.04,64.76,67.21,66.68,123.27,122.97,122.32, 122.29,122.42]
popSize3 = [3,3,3,3,3,5,5,5,5,5,10,10,10,10,10]
objVal3 = [1302,1302,1302,1302,1302,1302,1302,1302,1302,1323,1302,1302,1302,1302,1323]
avgObj3 = [1363.00, 1390.00,1341.00,1302.00, 1346.33, 1521.20,1403.40,1398.40,1497.00,1453.00, 1610.50, 1619.90,1654.90,1600.30,1678.80]
gap3 = [0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,1.61,0.00,0.00,0.00,0.00,1.61]
cpuTime3 = [57.17, 57.09,56.90,57.42,56.91,94.06,93.96, 93.88,93.97,93.81,182.90,189.32,199.76,200.02,193.59]
dSet = dSet1 + dSet2 + dSet3
popSize = popSize1 + popSize2 + popSize3
objVal = objVal1 + objVal2 + objVal3
avgObj = avgObj1 + avgObj2 + avgObj3
gap = gap1 + gap2 + gap3
cpuTime = cpuTime1 + cpuTime2 + cpuTime3
dfDict1 = {'Pop Size':popSize1, 'Obj Val':objVal1, 'Pop Avg Obj':avgObj1, '%Gap':gap1, 'CPU Time (s)':cpuTime1 }
ds1 = pd.DataFrame(dfDict1)
ds1 = ds1[['Pop Size', 'Obj Val', 'Pop Avg Obj', '%Gap', 'CPU Time (s)']]
dfDict2 = {'Pop Size':popSize2, 'Obj Val':objVal2, 'Pop Avg Obj':avgObj2, '%Gap':gap2, 'CPU Time (s)':cpuTime2 }
ds2 = pd.DataFrame(dfDict2)
ds2 = ds2[['Pop Size', 'Obj Val', 'Pop Avg Obj', '%Gap', 'CPU Time (s)']]
dfDict3 = {'Pop Size':popSize3, 'Obj Val':objVal3, 'Pop Avg Obj':avgObj3, '%Gap':gap3, 'CPU Time (s)':cpuTime3 }
ds3 = pd.DataFrame(dfDict3)
ds3 = ds3[['Pop Size', 'Obj Val', 'Pop Avg Obj', '%Gap', 'CPU Time (s)']]
dfDict = {'Pop Size':popSize, 'Obj Val':objVal, 'Pop Avg Obj':avgObj, '%Gap':gap, 'CPU Time (s)':cpuTime }
ds = pd.DataFrame(dfDict)
ds = ds[['Pop Size', 'Obj Val', 'Pop Avg Obj', '%Gap', 'CPU Time (s)']]
display(Markdown("_**Dataset1:**_"))
display(ds1)
print()
display(Markdown("_**Dataset2:**_"))
display(ds2)
print()
display(Markdown("_**Dataset3:**_"))
display(ds3)
Dataset1:
Pop Size | Obj Val | Pop Avg Obj | %Gap | CPU Time (s) | |
---|---|---|---|---|---|
0 | 3 | 4534 | 4929.33 | 0.0 | 24.80 |
1 | 3 | 4534 | 4840.33 | 0.0 | 25.36 |
2 | 3 | 4534 | 5002.00 | 0.0 | 26.22 |
3 | 3 | 4534 | 5154.33 | 0.0 | 26.35 |
4 | 3 | 4534 | 4836.00 | 0.0 | 25.44 |
5 | 5 | 4534 | 5066.80 | 0.0 | 43.37 |
6 | 5 | 4534 | 5356.20 | 0.0 | 43.84 |
7 | 5 | 4534 | 5663.20 | 0.0 | 43.72 |
8 | 5 | 4534 | 4850.20 | 0.0 | 44.11 |
9 | 5 | 4534 | 5050.60 | 0.0 | 43.64 |
10 | 10 | 4534 | 5530.40 | 0.0 | 86.21 |
11 | 10 | 4534 | 5839.50 | 0.0 | 86.90 |
12 | 10 | 4534 | 5671.60 | 0.0 | 82.90 |
13 | 10 | 4534 | 5725.70 | 0.0 | 80.21 |
14 | 10 | 4534 | 5459.10 | 0.0 | 68.50 |
Dataset2:
Pop Size | Obj Val | Pop Avg Obj | %Gap | CPU Time (s) | |
---|---|---|---|---|---|
0 | 3 | 920 | 1001.00 | 0.00 | 41.74 |
1 | 3 | 929 | 1002.33 | 0.98 | 42.15 |
2 | 3 | 920 | 1055.00 | 0.00 | 45.23 |
3 | 3 | 920 | 964.67 | 0.00 | 40.75 |
4 | 3 | 937 | 1028.67 | 1.85 | 39.75 |
5 | 5 | 931 | 1052.40 | 1.20 | 64.36 |
6 | 5 | 937 | 1185.20 | 1.85 | 64.04 |
7 | 5 | 937 | 1074.00 | 1.85 | 64.76 |
8 | 5 | 926 | 1076.00 | 0.65 | 67.21 |
9 | 5 | 920 | 1063.60 | 0.00 | 66.68 |
10 | 10 | 937 | 1177.00 | 1.85 | 123.27 |
11 | 10 | 931 | 1166.00 | 1.20 | 122.97 |
12 | 10 | 929 | 1139.80 | 0.98 | 122.32 |
13 | 10 | 938 | 1182.50 | 1.96 | 122.29 |
14 | 10 | 937 | 1181.40 | 1.85 | 122.42 |
Dataset3:
Pop Size | Obj Val | Pop Avg Obj | %Gap | CPU Time (s) | |
---|---|---|---|---|---|
0 | 3 | 1302 | 1363.00 | 0.00 | 57.17 |
1 | 3 | 1302 | 1390.00 | 0.00 | 57.09 |
2 | 3 | 1302 | 1341.00 | 0.00 | 56.90 |
3 | 3 | 1302 | 1302.00 | 0.00 | 57.42 |
4 | 3 | 1302 | 1346.33 | 0.00 | 56.91 |
5 | 5 | 1302 | 1521.20 | 0.00 | 94.06 |
6 | 5 | 1302 | 1403.40 | 0.00 | 93.96 |
7 | 5 | 1302 | 1398.40 | 0.00 | 93.88 |
8 | 5 | 1302 | 1497.00 | 0.00 | 93.97 |
9 | 5 | 1323 | 1453.00 | 1.61 | 93.81 |
10 | 10 | 1302 | 1610.50 | 0.00 | 182.90 |
11 | 10 | 1302 | 1619.90 | 0.00 | 189.32 |
12 | 10 | 1302 | 1654.90 | 0.00 | 199.76 |
13 | 10 | 1302 | 1600.30 | 0.00 | 200.02 |
14 | 10 | 1323 | 1678.80 | 1.61 | 193.59 |
# Number of population
Npop = 3
# Probability of crossover
Pc = 1.0
# Probability of mutation
Pm = 1.0
# Stopping number for generation
stopGeneration = 10000
# Start Timer
t1 = time.clock()
# Creating the initial population
population = initialization(Npop)
# Run the algorithm for 'stopGeneration' times generation
for i in range(stopGeneration):
# Selecting parents
parents = selection(population)
childs = []
# Apply crossover
for p in parents:
r = random.random()
if r < Pc:
childs.append(crossover([population[p[0]], population[p[1]]]))
else:
if r < 0.5:
childs.append(population[p[0]])
else:
childs.append(population[p[1]])
# Apply mutation
for c in childs:
r = random.random()
if r < Pm:
c = mutation(c)
# Update the population
population = elitistUpdate(population, childs)
#print(population)
#print(findBestSolution(population))
# Stop Timer
t2 = time.clock()
# Results Time
bestSol, bestObj, avgObj = findBestSolution(population)
print("Population:")
print(population)
print()
print("Solution:")
print(population[bestSol])
print()
print("Objective Value:")
print(bestObj)
print()
print("Average Objective Value of Population:")
print("%.2f" %avgObj)
print()
print("%Gap:")
G = 100 * (bestObj-optimalObjective) / optimalObjective
print("%.2f" %G)
print()
print("CPU Time (s)")
timePassed = (t2-t1)
print("%.2f" %timePassed)
Population: [[11, 2, 4, 8, 5, 9, 6, 12, 10, 1, 7, 0, 3], [11, 2, 4, 6, 12, 9, 10, 8, 1, 7, 5, 0, 3], [0, 8, 3, 9, 11, 2, 4, 6, 12, 10, 1, 7, 5]] Solution: [11, 2, 4, 6, 12, 9, 10, 8, 1, 7, 5, 0, 3] Objective Value: 920 Average Objective Value of Population: 1003.67 %Gap: 0.00 CPU Time (s) 43.08