# A Mutualism Quantum Genetic Algorithm to Optimize the Flow Shop Scheduling with Pickup and Delivery Considerations

Authors: Jinwei Gu; Manzhan Gu; Xingsheng Gu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A mutualism quantum genetic algorithm (MQGA) is proposed for an integrated supply chain scheduling with the materials pickup, flow shop scheduling, and the finished products delivery. The objective is to minimize the makespan, that is, the arrival time of the last finished product to the customer. In MQGA, a new symbiosis strategy named mutualism is proposed to adjust the size of each population dynamically by regarding the mutual influence relation of the two subpopulations. A hybridQ-bit coding method and a local speeding-up method are designed to increase the diversity of genes, and a checking routine is carried out to ensure the feasibility of each solution; that is, the total physical space of each delivery batch could not exceed the capacity of the vehicle. Compared with the modified genetic algorithm (MGA) and the quantum-inspired genetic algorithm (QGA), the effectiveness and efficiency of the MQGA are validated by numerical experiments.


---

## Body

## 1. Introduction

The coordination of logistics activities in a supply chain has received a lot of attention recently. From the manufacturer’s point of view, the important problem in the supply chain scheduling is the coordination of the three stages including material supply, production scheduling, and product delivery. Traditionally, research on scheduling generally focuses on the models with various machine setting, job characteristics, and performance measures [1], and the transportation arrangements of the materials and the finished products are normally ignored in these models.In the past decades, many researchers studied the machine scheduling problems with transportation under consideration. Some of the results consider the transportation of the jobs between machines in the flow shop model, which incorporates transport times when the jobs are transferred from one machine to another. Maggu et al. [2–4] firstly considered the scheduling problem with transportation between machines. In addition, Langston [5] designed algorithms for the problem of planning and coordinating movement within a deterministic flow shop system and analyzed the worst-case performances of the algorithms. Haouari and Ladhari [6] presented an effective branch and bound algorithm for the model proposed by Langston [5]. Hurink and Knust [7] considered the flow shop problems with transportation and a single robot with the objective of minimizing the makespan. They derived complexity results for the problems. Hurink and Knust [8] developed a tabu search algorithm for the problem with time windows. Hurink and Knust [9] extended some of the results in [8] to the job-shop scheduling models. For the model with a single robot, Lee and Strusevich [10] studied the two machine flow shop and open shop problems. They presented a best possible approximation algorithm for each of the two problems with some constraints. Naderi et al. [11] proposed simulated annealing algorithms for flow shop problems with the objective of minimizing the total weighted tardiness and makespan. For more new results on the scheduling problems with transportation between machines, please see [12–15], among others.Another type of scheduling model with transportation focuses on the delivery of finished jobs to customers. Lee and Chen [16] considered some scheduling models that incorporate the delivery decisions of the finished jobs. They researched the computational complexity of some problems and proposed polynomial or pseudopolynomial algorithms for them. Chang and Lee [17] extended one of the models in [16] to the environment in which each job has an individual amount of space. Li et al. [18] studied a single-machine scheduling problem with routing decisions of a vehicle that serves customers at different addresses. In addition, Chen and Vairaktarakis [19] researched eight problems with production and distribution under consideration and either designed a polynomial time algorithm or proved the NP-hardness for each of the problems. Geismar et al. [20] extended the Chen and Vairaktarakis model to the problem with a short shelf life and no inventory in the process. In addition, Soukhal et al. [21] investigated the scheduling models with constraints on both transportation and buffer capacities to minimize the makespan. They proved new complexity results for special cases of the problems considered. We refer to, for example, [22–25], for more recent results on the model with delivery operations for finished products.This paper considers the scheduling model that integrates the pickup of materials, flow shop scheduling, and the delivery of finished jobs. In this model, the material warehouse, the factory, and the customer are located at three different places. There are two vehicles (namely, conveyor and truck) each with a limited capacity. One vehicle (conveyor) travels between the factory and the warehouse for material transportation, and the other vehicle (truck) travels between the customer and the factory for finished products delivery. This model applies to many situations in supply chain business activities. For example, a shoes manufacturer purchases materials from India and arranges the production in China. Finally, the finished products are delivered to the USA.To the best of our knowledge, research on this model includes Hall and Potts [26], Li and Ou [27], and Wang and Cheng [28]. Hall and Potts [26] provided algorithms for models in which jobs are produced on machines and formed in batches for delivery. Li and Ou [27] studied the problem with a vehicle traveling between a warehouse and the factory. They proved the problem is NP-hard in strong sense and developed polynomial time and effective heuristic algorithms, respectively, for special case and the general setting. Wang and Cheng [28] studied the complexity of the model in which the warehouse, the factory, and the customer locations are different. They proved the problem is strongly NP-hard and proposed a heuristic with a tight bound of 2.The model considered in this paper is different from that in [26]. Each job in our model has its own transportation time between the warehouse and the factory and between the factory and the customer. Our model also differs from Li and Ou’s model in which the locations of the warehouse and the customer are the same. Furthermore, this paper considers flow shop during the production scheduling, while the model in [28] considers processing jobs in a single machine.Quantum genetic algorithms (QGA) are heuristic search techniques inspired from the principles of survival of the fittest in natural genetic evolution and quantum theory. They are known to be efficient in a large search space, without explicitly requiring additional information (such as convexity or derivative information) about the objective. In addition, theQ-bit encoding method has the superposition and probability expression characteristics, and each individual can express more states to increase the diversity of populations. For these reasons, in the last few years, they have been applied to many combinatorial problems, including scheduling and vehicle routing applications that are partially related to our problem.In order to solve the integrated supply chain scheduling model with the materials pickup, flow shop scheduling, and the finished products delivery, a technique called mutualism quantum genetic algorithm (MQGA) is developed, which differs from traditional QGA method in two aspects.The first aspect is the encoding method. With the encoding method in [29], a chromosome including two segments, the job sequence and the job-to-batch assignment, is designed, and this encoding method could equally and uniquely represent all possible solutions. Furthermore, in order to improve the encoding performance, a local speeding-up method is introduced to perform the deep search operation. After that, a checking routine is carried out to ensure the feasibility of each solution, that is, ensure that the total physical space of each delivery batch could not exceed the capacity of the vehicle.The second aspect is as follows: based on the relationship between two mutually symbiotic species, a population growth model (namely, mutualism) is introduced to improve the performance of the algorithm. In this model, a symbiosis strategy is developed to dynamically adjust the sizes of two subpopulations. This strategy helps to increase the diversity of genes and avoid premature convergence.The remainder of this paper is organized as follows. In Section2, the problem including the assumptions and notations is introduced. In Section 3, the two species growth model and the corresponding evolutionary strategy are proposed. The mechanism of MQGA is introduced in Section 4, and the experimental results are discussed in Section 5. Section 6 contains some concluding remarks.


## 2. Problem Description

In this paper, it is assumed that the material warehouse, the factory, and the customer are located at different addresses, and there is a set of jobs,J
=
J
1
,
J
2
,
…
,
J
n, processed in a flow shop scheduling, which consists of m machines M
=
M
1
,
M
2
,
…
,
M
m. All unprocessed jobs are initially located at the supplier’s warehouse and need to be transported to the factory for processing by a vehicle (conveyor). The conveyor is located at the warehouse at the beginning and is available for transportation of jobs from warehouse to the factory and comes back. Each pickup journey takes a constant amount of time, and in a trip the conveyor loads a limited number of jobs due to the limitation of its space capacity. In addition, the finished jobs need to be delivered to a customer by another vehicle (truck), which is initially located at the factory. Because of the constraint of its space capacity, the truck transports limited finished jobs to the customer and returns to the factory. Each trip also takes a constant amount of time. The objective is to minimize the makespan, that is, the arrival time of the last finished job to the customer. Figure 1 gives a description of the integrated supply chain problem including the pickup of materials, flow shop scheduling, and the delivery of finished products.Figure 1
Schematic of the cooperative problem of flow shop production with pickup and delivery transportations.


### 2.1. Assumptions

Before the introduction of the notations, a number of assumptions are given as follows:(i)
all the facilities including conveyor, machines, and truck are available from time zero;
(ii)
there is no idle time between any consecutive two pickup journeys;
(iii)
the time of loading and unloading jobs is included in the pickup time and delivery time;
(iv)
the time of transporting jobs between machines is negligible;
(v)
the storage or buffer capacities between successive machines are unlimited;
(vi)
there is no priority among jobs;
(vii)
machine failure is not considered.


### 2.2. Notations

The following notations are used in the problem model:c
i
: physical space of jobJ
i;
S
1
S
2: physical space capacity of the conveyor (truck);
T
11
T
12: travel time from the warehouse to the factory (from the factory to the warehouse);
T
21
T
22: travel time from the factory to the customer (from the customer to the factory);
l
1
i: the arrival time of J
i to the factory;
l
2
i: the arrival time of J
i to the customer;
p
i
k: processing time of J
i on machine M
k;
s
k
o: the earliest start time of the oth position job on machine M
k;
u
1
j: departure time of the jth batch from the warehouse to the workshop;
u
2
j: departure time of the jth batch from the workshop to the customer;
q
i
k
o: equal to 1 if job J
i is scheduled on the oth position on machine M
k; 0, otherwise;
L
j: the latest completion time of the jobs in the jth delivery batch;
U
i
j: equal to 1 if job J
i is in the jth pickup batch; 0, otherwise;
V
i
j: equal to 1 if job J
i is in the jth delivery batch; 0, otherwise;
M: very large positive constant.


### 2.3. The Optimization Model

Tang and Gong [30] proposed a mixed integer programming (MIP) model for a single-machine batch scheduling with pickup and delivery transportations. Based on it, this paper studies a different integrated supply chain problem, which consists of the pickup of materials, flow shop scheduling, and the delivery of finished products. For this problem, a MIP model is constructed below. This model is to decide the following policy: (
1
) how and when the materials are assigned to an appropriate delivery batch, (
2
) what the scheduling ofn jobs on a series of m machines is, (
3
) how to dispatch the finished product as soon as possible to satisfy the customer’s demand:(1)
min
⁡
max
⁡
1
≤
i
≤
n
l
2
i
(2)
Subject
to
∑
o
=
1
n
q
i
k
o
=
1
,
i
=
1,2
,
…
,
n
,
k
=
1,2
,
…
,
m
,
(3)
∑
i
=
1
n
q
i
k
o
=
1
,
o
=
1,2
,
…
,
n
,
k
=
1,2
,
…
,
m
,
(4)
s
k
o
+
∑
i
=
1
n
p
i
k
q
i
k
o
≤
s
k
,
o
+
1
,
k
=
1,2
,
…
,
m
,
o
=
1,2
,
…
,
n
-
1
(5)
s
k
o
+
∑
i
=
1
n
p
i
k
q
i
k
o
≤
s
k
+
1
,
o
k
=
1,2
,
…
,
m
-
1
,
o
=
1,2
,
…
,
n
,
(6)
s
1
o
≥
l
1
i
-
M
1
-
q
i
1
o
,
i
=
1,2
,
…
,
n
,
o
=
1,2
,
…
,
n
,
(7)
L
j
≥
s
m
o
+
p
i
m
-
M
2
-
V
i
j
-
q
i
m
o
,
i
=
1,2
,
…
,
n
,
o
=
1,2
,
…
,
n
,
(8)
∑
i
=
1
n
c
i
U
i
j
≤
S
1
,
j
=
1,2
,
…
,
n
,
(9)
∑
i
=
1
n
c
i
V
i
j
≤
S
2
,
j
=
1,2
,
…
,
n
,
(10)
∑
j
=
1
n
U
i
j
=
1
,
i
=
1,2
,
…
,
n
,
(11)
∑
j
=
1
n
V
i
j
=
1
,
i
=
1,2
,
…
,
n
,
(12)
∑
i
=
1
n
U
i
(
j
+
1
)
≤
M
∑
i
=
1
n
U
i
j
,
j
=
1,2
,
…
,
n
-
1
,
(13)
∑
i
=
1
n
V
i
(
j
+
1
)
≤
M
∑
i
=
1
n
V
i
j
,
j
=
1,2
,
…
,
n
-
1
,
(14)
u
21
=
L
1
,
(15)
u
2
j
≥
L
j
,
j
=
2,3
,
…
,
n
,
(16)
u
1
j
≥
u
1
(
j
-
1
)
+
T
11
+
T
12
,
j
=
2,3
,
…
,
n
,
(17)
u
2
j
≥
u
2
(
j
-
1
)
+
T
21
+
T
22
,
j
=
2,3
,
…
,
n
,
(18)
l
1
i
≥
u
1
j
+
T
11
-
M
1
-
U
i
j
,
i
=
1,2
,
…
,
n
,
j
=
1,2
,
…
,
n
,
(19)
l
2
i
≥
u
2
j
+
T
21
-
M
1
-
V
i
j
,
i
=
1,2
,
…
,
n
,
j
=
1,2
,
…
,
n
,
(20)
l
2
i
≥
0
,
i
=
1,2
,
…
,
n
,
(21)
L
j
≥
0
,
j
=
1,2
,
…
,
n
,
(22)
u
1
j
≥
0
,
u
2
j
≥
0
,
j
=
1,2
,
…
,
n
,
(23)
s
k
o
≥
0
,
k
=
1,2
,
…
,
m
,
o
=
1,2
,
…
,
n
,
(24)
U
i
j
,
V
i
j
∈
{
0,1
}
,
i
=
1,2
,
…
,
n
,
j
=
1,2
,
…
,
n
,
(25)
q
i
k
o
∈
{
0,1
}
,
i
=
1,2
,
…
,
n
,
k
=
1,2
,
…
,
m
,
o
=
1,2
,
…
,
n
.Objective function (1) represents the arrival time of the last finished product to the customer. Constraints (2) ensure that job J
i should be placed in one and only one position on each machine. Constraints (3) guarantee that there is one and only one job processed on the oth position on each machine. Constraints (4) define that the earliest start time of the job on the (
o
+
1
)th position should not be earlier than the completion time of the oth position job on the same machine. Constraints (5) ensure that, for the same job, the earliest start time on machine M
k
+
1 should not be earlier than its completion time on machine M
k. Constraints (6) define that the start processing time of a job should not be earlier than its release time. Constraints (7) indicate that the ready time for delivering the oth job should not be earlier than its completion time on the last machine. Constraints (8) and (9) define that the number of jobs in each batch should satisfy the space capacity of the conveyor or the truck. Constraints (10) and (11) determine each job to be assigned to one and only one batch. Constraints (12) and (13) indicate that if no job is assigned to batch j, then no job will be assigned to batch (
j
+
1
). Constraints (14) and (15) define the departure time u
1
j. Constraints (16) and (17) define that the departure time of the jth batch should not be earlier than the returning time of the (
j
-
1
)th batch. Constrains (18) and (19) specify the nonnegativity of l
1
i and l
2
i. Equations (20)–(25) define the variables.


## 2.1. Assumptions

Before the introduction of the notations, a number of assumptions are given as follows:(i)
all the facilities including conveyor, machines, and truck are available from time zero;
(ii)
there is no idle time between any consecutive two pickup journeys;
(iii)
the time of loading and unloading jobs is included in the pickup time and delivery time;
(iv)
the time of transporting jobs between machines is negligible;
(v)
the storage or buffer capacities between successive machines are unlimited;
(vi)
there is no priority among jobs;
(vii)
machine failure is not considered.


## 2.2. Notations

The following notations are used in the problem model:c
i
: physical space of jobJ
i;
S
1
S
2: physical space capacity of the conveyor (truck);
T
11
T
12: travel time from the warehouse to the factory (from the factory to the warehouse);
T
21
T
22: travel time from the factory to the customer (from the customer to the factory);
l
1
i: the arrival time of J
i to the factory;
l
2
i: the arrival time of J
i to the customer;
p
i
k: processing time of J
i on machine M
k;
s
k
o: the earliest start time of the oth position job on machine M
k;
u
1
j: departure time of the jth batch from the warehouse to the workshop;
u
2
j: departure time of the jth batch from the workshop to the customer;
q
i
k
o: equal to 1 if job J
i is scheduled on the oth position on machine M
k; 0, otherwise;
L
j: the latest completion time of the jobs in the jth delivery batch;
U
i
j: equal to 1 if job J
i is in the jth pickup batch; 0, otherwise;
V
i
j: equal to 1 if job J
i is in the jth delivery batch; 0, otherwise;
M: very large positive constant.


## 2.3. The Optimization Model

Tang and Gong [30] proposed a mixed integer programming (MIP) model for a single-machine batch scheduling with pickup and delivery transportations. Based on it, this paper studies a different integrated supply chain problem, which consists of the pickup of materials, flow shop scheduling, and the delivery of finished products. For this problem, a MIP model is constructed below. This model is to decide the following policy: (
1
) how and when the materials are assigned to an appropriate delivery batch, (
2
) what the scheduling ofn jobs on a series of m machines is, (
3
) how to dispatch the finished product as soon as possible to satisfy the customer’s demand:(1)
min
⁡
max
⁡
1
≤
i
≤
n
l
2
i
(2)
Subject
to
∑
o
=
1
n
q
i
k
o
=
1
,
i
=
1,2
,
…
,
n
,
k
=
1,2
,
…
,
m
,
(3)
∑
i
=
1
n
q
i
k
o
=
1
,
o
=
1,2
,
…
,
n
,
k
=
1,2
,
…
,
m
,
(4)
s
k
o
+
∑
i
=
1
n
p
i
k
q
i
k
o
≤
s
k
,
o
+
1
,
k
=
1,2
,
…
,
m
,
o
=
1,2
,
…
,
n
-
1
(5)
s
k
o
+
∑
i
=
1
n
p
i
k
q
i
k
o
≤
s
k
+
1
,
o
k
=
1,2
,
…
,
m
-
1
,
o
=
1,2
,
…
,
n
,
(6)
s
1
o
≥
l
1
i
-
M
1
-
q
i
1
o
,
i
=
1,2
,
…
,
n
,
o
=
1,2
,
…
,
n
,
(7)
L
j
≥
s
m
o
+
p
i
m
-
M
2
-
V
i
j
-
q
i
m
o
,
i
=
1,2
,
…
,
n
,
o
=
1,2
,
…
,
n
,
(8)
∑
i
=
1
n
c
i
U
i
j
≤
S
1
,
j
=
1,2
,
…
,
n
,
(9)
∑
i
=
1
n
c
i
V
i
j
≤
S
2
,
j
=
1,2
,
…
,
n
,
(10)
∑
j
=
1
n
U
i
j
=
1
,
i
=
1,2
,
…
,
n
,
(11)
∑
j
=
1
n
V
i
j
=
1
,
i
=
1,2
,
…
,
n
,
(12)
∑
i
=
1
n
U
i
(
j
+
1
)
≤
M
∑
i
=
1
n
U
i
j
,
j
=
1,2
,
…
,
n
-
1
,
(13)
∑
i
=
1
n
V
i
(
j
+
1
)
≤
M
∑
i
=
1
n
V
i
j
,
j
=
1,2
,
…
,
n
-
1
,
(14)
u
21
=
L
1
,
(15)
u
2
j
≥
L
j
,
j
=
2,3
,
…
,
n
,
(16)
u
1
j
≥
u
1
(
j
-
1
)
+
T
11
+
T
12
,
j
=
2,3
,
…
,
n
,
(17)
u
2
j
≥
u
2
(
j
-
1
)
+
T
21
+
T
22
,
j
=
2,3
,
…
,
n
,
(18)
l
1
i
≥
u
1
j
+
T
11
-
M
1
-
U
i
j
,
i
=
1,2
,
…
,
n
,
j
=
1,2
,
…
,
n
,
(19)
l
2
i
≥
u
2
j
+
T
21
-
M
1
-
V
i
j
,
i
=
1,2
,
…
,
n
,
j
=
1,2
,
…
,
n
,
(20)
l
2
i
≥
0
,
i
=
1,2
,
…
,
n
,
(21)
L
j
≥
0
,
j
=
1,2
,
…
,
n
,
(22)
u
1
j
≥
0
,
u
2
j
≥
0
,
j
=
1,2
,
…
,
n
,
(23)
s
k
o
≥
0
,
k
=
1,2
,
…
,
m
,
o
=
1,2
,
…
,
n
,
(24)
U
i
j
,
V
i
j
∈
{
0,1
}
,
i
=
1,2
,
…
,
n
,
j
=
1,2
,
…
,
n
,
(25)
q
i
k
o
∈
{
0,1
}
,
i
=
1,2
,
…
,
n
,
k
=
1,2
,
…
,
m
,
o
=
1,2
,
…
,
n
.Objective function (1) represents the arrival time of the last finished product to the customer. Constraints (2) ensure that job J
i should be placed in one and only one position on each machine. Constraints (3) guarantee that there is one and only one job processed on the oth position on each machine. Constraints (4) define that the earliest start time of the job on the (
o
+
1
)th position should not be earlier than the completion time of the oth position job on the same machine. Constraints (5) ensure that, for the same job, the earliest start time on machine M
k
+
1 should not be earlier than its completion time on machine M
k. Constraints (6) define that the start processing time of a job should not be earlier than its release time. Constraints (7) indicate that the ready time for delivering the oth job should not be earlier than its completion time on the last machine. Constraints (8) and (9) define that the number of jobs in each batch should satisfy the space capacity of the conveyor or the truck. Constraints (10) and (11) determine each job to be assigned to one and only one batch. Constraints (12) and (13) indicate that if no job is assigned to batch j, then no job will be assigned to batch (
j
+
1
). Constraints (14) and (15) define the departure time u
1
j. Constraints (16) and (17) define that the departure time of the jth batch should not be earlier than the returning time of the (
j
-
1
)th batch. Constrains (18) and (19) specify the nonnegativity of l
1
i and l
2
i. Equations (20)–(25) define the variables.


## 3. The Mutualism Strategy for Population Growth

As the long evolutionary process of nature, the relationship between living beings is complicated. The phenomenon that two species live together is generally referred to as the symbiosis. This paper studies a kind of symbiosis called mutualism, which refers to two species living together and depending on each other over a long period of time. The nutrition of one species is the food source of the other one.Mutualism brings two advantages to our algorithm. Firstly, there is no need to design the fitness function, and the fitness value of an individual is obtained by the cooperation between two populations, which could reduce the dependence on the domain knowledge. Secondly, different from the traditional framework in which the evolution of populations depends on the fitness values of individuals, mutualism incorporates the cooperative behavior between an individual and its surroundings, so as to postpone the premature convergence and improve the convergence speed.


### 3.1. The Mutualism Population Growth Model

A differential equation of mutualism population growth model is introduced in this subsection. Assume two species A and B live in the same environment and they do not take each other as food, but the existence of one species can promote the population growth of the other one. For example, algae and fungi are two species living in lichen. Algae provide nutrient to fungi through photosynthesis, and fungi offer algae water and inorganic substance. If algae and fungi are separated in lichen, both of them will die.Letx
(
t
) and y
(
t
), respectively, be the numbers of individuals of species A and B at the tth generation. Define the maximal individual number of species A (B) by K
1
K
2, and the cooperation degree is a
1
a
2. In this way, the differential equation of the cooperation system can be written as follows:(26)
d
x
(
t
)
d
t
=
r
1
1
-
x
t
K
2
+
a
1
y
t
K
1
x
t
,
(27)
d
y
(
t
)
d
t
=
r
2
1
+
a
2
x
t
K
2
-
y
t
K
1
y
t
,where r
i
i
=
1,2 denotes the accrual rate, that is, birth rate subtracting death rate. In time interval t
,
t
+
Δ
t, define x
t
+
Δ
t
-
x
t
=
r
1
·
Δ
t
·
x
t as the modified number of species A. So Δ
x
/
Δ
t
=
r
1
x
(
t
). Once the population number reaches the maximal value, species A will stop growing, and the accrual rate r
1 will become zero. 1
-
x
(
t
)
/
K
2 denotes the block function on the population size of species A. Suppose the total amount of food is 1; then x
(
t
)
/
K
2 is the amount of food consumed by species A. When species A and B live together in the same environment, they cooperate with each other and bring direct advantages for population growth, and α
1
y
(
t
)
/
K
1 denotes the influence degree of species B on the population size of species A. We have similar explanation for (27).


### 3.2. The Mutualism Strategy

In the standard genetic algorithm, the fitness function will be given as an input, and the fitness is easy to reach the peak value. However, in the real environment the adaptability of a species is dynamic and is affected by its surroundings. Therefore, in MQGA, the fitness value of an individual is replaced by its cooperation degree with individuals in other species, and the details on the computation of the cooperation degree are given below.


#### 3.2.1. Expert Evaluation Method

Given two individualsi and j, respectively, from populations A and B, a method is proposed to evaluate the cooperation degree of individual i. Let o
b
i and o
b
j be the objective values of i and j. Denote, respectively, by o
b
min
⁡, o
b
max
⁡, and ob
¯ the minimum value, the maximum value, and the average objective value of all individuals in A and B. By cooperating with individual j, define Scor
e
i
j as the cooperative degree of individual i.Case 1.
Considero
b
i
∈
o
b
min
⁡
,
ob
¯, o
b
j
∈
o
b
min
⁡
,
ob
¯. In this case, each of the two individuals has a high survival rate and the “matching” ability to cooperate. They can offer nutrition to each other, and the collaborative effect is the best. So, the cooperative degree of individual i is set to be “4.”Case 2.
Considero
b
i
∈
o
b
min
⁡
,
ob
¯, o
b
j
∈
ob
¯
,
o
b
max
⁡. In this case, compared to individual j, individual i has a smaller objective value, which indicates that individual i is more capable of growing and developing. If the two individuals work together, individual i cannot obtain enough nutrition from individual j, and the cooperation is a “mismatch.” Therefore, the cooperative degree of individual i is set to be “2.”Case 3.
Considero
b
i
∈
ob
¯
,
o
b
max
⁡, o
b
j
∈
o
b
min
⁡
,
ob
¯. Known from Case 2, in this case individual j is more capable of growing and developing, and individual j will absorb more nutrition from individual i. That will eventually impede the development and growth of individual i. Therefore, the cooperative degree of individual i is defined to be “0.”Case 4.
Considero
b
i
∈
ob
¯
,
o
b
max
⁡, o
b
j
∈
ob
¯
,
o
b
max
⁡. In this case, the two individuals have the similar low levels of growing, and a long period of time is needed for the development of the individuals. But they have the “matching” abilities and can cooperate with each other. Therefore, the cooperative degree of individual i is set to be “1.” The four cases are summarized in Table 1.Table 1
The cooperation degree of individuali (with individual j) Score
i
j.
Condition
Score
o
b
i
∈
[
ob
min
⁡
,
ob
¯
], ob
j
∈
[
ob
min
⁡
,
ob
¯
]
4
ob
i
∈
[
ob
min
⁡
,
ob
¯
], ob
j
∈
[
ob
¯
,
ob
max
⁡
]
2
ob
i
∈
[
ob
¯
,
ob
max
⁡
], ob
j
∈
[
ob
min
⁡
,
ob
¯
]
0
ob
i
∈
[
ob
¯
,
ob
max
⁡
], ob
j
∈
[
ob
¯
,
ob
max
⁡
]
1


#### 3.2.2. Cooperative Degree of Individual

Given two populations A and B, denote byc
Fitnes
s
i the cooperative degree of individual i in A. Firstly choose some individuals from B to form a delegation X, and X can be chosen by ways including random selection, greedy selection, or tournament. Then let individual i cooperate with each individual in X with the expert evaluation method. The fitness of i is calculated as follows:(28)
c
Fitnes
s
i
=
∑
j
∈
X
Scor
e
i
j
,where Scor
e
i
j denotes the cooperative degree of individual i by cooperating with individual j.


#### 3.2.3. Cooperative Degree of a Population

The cooperative degree of a population stands for the cooperative level of the population, and it can be calculated as follows:(29)
a
1
=
∑
i
∈
Pop
1
c
Fitnes
s
i
∑
i
∈
Pop
2
c
Fitnes
s
i
,
a
2
=
1
a
1
,where Pop1 (Pop2) is the delegation of population A (B). If a
1
>
1, A is more cooperative and obtains more advantage or nutrition from B. Conversely, if a
2
>
1, B is more cooperative than A. The two parameters a
1 and a
2 show the mutual influence degree on population sizes.


#### 3.2.4. Adjustment of Population Sizes

In our mutualism population growth model, assume(30)
r
1
=
K
1
-
x
(
t
)
s
1
,
r
2
=
K
2
-
y
(
t
)
s
2
,where s
1
,
s
2 are constants. Based on (26) and (27), at each generation t, the population sizes of A and B will be modified by Δ
x and Δ
y, respectively, which are given as follows: (31)
Δ
x
=
K
1
-
x
s
1
1
-
x
K
2
+
a
1
y
K
1
x
,
Δ
y
=
K
2
-
y
s
2
1
+
a
2
x
K
2
-
y
K
1
y
.The newly added individuals are generated randomly, and the deleted individuals are the ones with the smallest cooperative degree values.If the cooperative degree of a population is high, its population size will increase. Thus the growth of a population size mostly depends on its ability to coordinate. If the ability is strong, the size will become bigger and bigger till the maximum value. Otherwise, the size will become small and eventually extinct. In order to survive, each population will try its best to gain good genes and cooperate with others. In this way, the quality of solution is improved and the performance of algorithm is enhanced.


### 3.3. Framework of the Mutualism Strategy for Population Growth

A description of the mutualism strategy for population growth is given in Figure2.Figure 2
Description of the mutualism strategy for population growth.


## 3.1. The Mutualism Population Growth Model

A differential equation of mutualism population growth model is introduced in this subsection. Assume two species A and B live in the same environment and they do not take each other as food, but the existence of one species can promote the population growth of the other one. For example, algae and fungi are two species living in lichen. Algae provide nutrient to fungi through photosynthesis, and fungi offer algae water and inorganic substance. If algae and fungi are separated in lichen, both of them will die.Letx
(
t
) and y
(
t
), respectively, be the numbers of individuals of species A and B at the tth generation. Define the maximal individual number of species A (B) by K
1
K
2, and the cooperation degree is a
1
a
2. In this way, the differential equation of the cooperation system can be written as follows:(26)
d
x
(
t
)
d
t
=
r
1
1
-
x
t
K
2
+
a
1
y
t
K
1
x
t
,
(27)
d
y
(
t
)
d
t
=
r
2
1
+
a
2
x
t
K
2
-
y
t
K
1
y
t
,where r
i
i
=
1,2 denotes the accrual rate, that is, birth rate subtracting death rate. In time interval t
,
t
+
Δ
t, define x
t
+
Δ
t
-
x
t
=
r
1
·
Δ
t
·
x
t as the modified number of species A. So Δ
x
/
Δ
t
=
r
1
x
(
t
). Once the population number reaches the maximal value, species A will stop growing, and the accrual rate r
1 will become zero. 1
-
x
(
t
)
/
K
2 denotes the block function on the population size of species A. Suppose the total amount of food is 1; then x
(
t
)
/
K
2 is the amount of food consumed by species A. When species A and B live together in the same environment, they cooperate with each other and bring direct advantages for population growth, and α
1
y
(
t
)
/
K
1 denotes the influence degree of species B on the population size of species A. We have similar explanation for (27).


## 3.2. The Mutualism Strategy

In the standard genetic algorithm, the fitness function will be given as an input, and the fitness is easy to reach the peak value. However, in the real environment the adaptability of a species is dynamic and is affected by its surroundings. Therefore, in MQGA, the fitness value of an individual is replaced by its cooperation degree with individuals in other species, and the details on the computation of the cooperation degree are given below.


### 3.2.1. Expert Evaluation Method

Given two individualsi and j, respectively, from populations A and B, a method is proposed to evaluate the cooperation degree of individual i. Let o
b
i and o
b
j be the objective values of i and j. Denote, respectively, by o
b
min
⁡, o
b
max
⁡, and ob
¯ the minimum value, the maximum value, and the average objective value of all individuals in A and B. By cooperating with individual j, define Scor
e
i
j as the cooperative degree of individual i.Case 1.
Considero
b
i
∈
o
b
min
⁡
,
ob
¯, o
b
j
∈
o
b
min
⁡
,
ob
¯. In this case, each of the two individuals has a high survival rate and the “matching” ability to cooperate. They can offer nutrition to each other, and the collaborative effect is the best. So, the cooperative degree of individual i is set to be “4.”Case 2.
Considero
b
i
∈
o
b
min
⁡
,
ob
¯, o
b
j
∈
ob
¯
,
o
b
max
⁡. In this case, compared to individual j, individual i has a smaller objective value, which indicates that individual i is more capable of growing and developing. If the two individuals work together, individual i cannot obtain enough nutrition from individual j, and the cooperation is a “mismatch.” Therefore, the cooperative degree of individual i is set to be “2.”Case 3.
Considero
b
i
∈
ob
¯
,
o
b
max
⁡, o
b
j
∈
o
b
min
⁡
,
ob
¯. Known from Case 2, in this case individual j is more capable of growing and developing, and individual j will absorb more nutrition from individual i. That will eventually impede the development and growth of individual i. Therefore, the cooperative degree of individual i is defined to be “0.”Case 4.
Considero
b
i
∈
ob
¯
,
o
b
max
⁡, o
b
j
∈
ob
¯
,
o
b
max
⁡. In this case, the two individuals have the similar low levels of growing, and a long period of time is needed for the development of the individuals. But they have the “matching” abilities and can cooperate with each other. Therefore, the cooperative degree of individual i is set to be “1.” The four cases are summarized in Table 1.Table 1
The cooperation degree of individuali (with individual j) Score
i
j.
Condition
Score
o
b
i
∈
[
ob
min
⁡
,
ob
¯
], ob
j
∈
[
ob
min
⁡
,
ob
¯
]
4
ob
i
∈
[
ob
min
⁡
,
ob
¯
], ob
j
∈
[
ob
¯
,
ob
max
⁡
]
2
ob
i
∈
[
ob
¯
,
ob
max
⁡
], ob
j
∈
[
ob
min
⁡
,
ob
¯
]
0
ob
i
∈
[
ob
¯
,
ob
max
⁡
], ob
j
∈
[
ob
¯
,
ob
max
⁡
]
1


### 3.2.2. Cooperative Degree of Individual

Given two populations A and B, denote byc
Fitnes
s
i the cooperative degree of individual i in A. Firstly choose some individuals from B to form a delegation X, and X can be chosen by ways including random selection, greedy selection, or tournament. Then let individual i cooperate with each individual in X with the expert evaluation method. The fitness of i is calculated as follows:(28)
c
Fitnes
s
i
=
∑
j
∈
X
Scor
e
i
j
,where Scor
e
i
j denotes the cooperative degree of individual i by cooperating with individual j.


### 3.2.3. Cooperative Degree of a Population

The cooperative degree of a population stands for the cooperative level of the population, and it can be calculated as follows:(29)
a
1
=
∑
i
∈
Pop
1
c
Fitnes
s
i
∑
i
∈
Pop
2
c
Fitnes
s
i
,
a
2
=
1
a
1
,where Pop1 (Pop2) is the delegation of population A (B). If a
1
>
1, A is more cooperative and obtains more advantage or nutrition from B. Conversely, if a
2
>
1, B is more cooperative than A. The two parameters a
1 and a
2 show the mutual influence degree on population sizes.


### 3.2.4. Adjustment of Population Sizes

In our mutualism population growth model, assume(30)
r
1
=
K
1
-
x
(
t
)
s
1
,
r
2
=
K
2
-
y
(
t
)
s
2
,where s
1
,
s
2 are constants. Based on (26) and (27), at each generation t, the population sizes of A and B will be modified by Δ
x and Δ
y, respectively, which are given as follows: (31)
Δ
x
=
K
1
-
x
s
1
1
-
x
K
2
+
a
1
y
K
1
x
,
Δ
y
=
K
2
-
y
s
2
1
+
a
2
x
K
2
-
y
K
1
y
.The newly added individuals are generated randomly, and the deleted individuals are the ones with the smallest cooperative degree values.If the cooperative degree of a population is high, its population size will increase. Thus the growth of a population size mostly depends on its ability to coordinate. If the ability is strong, the size will become bigger and bigger till the maximum value. Otherwise, the size will become small and eventually extinct. In order to survive, each population will try its best to gain good genes and cooperate with others. In this way, the quality of solution is improved and the performance of algorithm is enhanced.


## 3.2.1. Expert Evaluation Method

Given two individualsi and j, respectively, from populations A and B, a method is proposed to evaluate the cooperation degree of individual i. Let o
b
i and o
b
j be the objective values of i and j. Denote, respectively, by o
b
min
⁡, o
b
max
⁡, and ob
¯ the minimum value, the maximum value, and the average objective value of all individuals in A and B. By cooperating with individual j, define Scor
e
i
j as the cooperative degree of individual i.Case 1.
Considero
b
i
∈
o
b
min
⁡
,
ob
¯, o
b
j
∈
o
b
min
⁡
,
ob
¯. In this case, each of the two individuals has a high survival rate and the “matching” ability to cooperate. They can offer nutrition to each other, and the collaborative effect is the best. So, the cooperative degree of individual i is set to be “4.”Case 2.
Considero
b
i
∈
o
b
min
⁡
,
ob
¯, o
b
j
∈
ob
¯
,
o
b
max
⁡. In this case, compared to individual j, individual i has a smaller objective value, which indicates that individual i is more capable of growing and developing. If the two individuals work together, individual i cannot obtain enough nutrition from individual j, and the cooperation is a “mismatch.” Therefore, the cooperative degree of individual i is set to be “2.”Case 3.
Considero
b
i
∈
ob
¯
,
o
b
max
⁡, o
b
j
∈
o
b
min
⁡
,
ob
¯. Known from Case 2, in this case individual j is more capable of growing and developing, and individual j will absorb more nutrition from individual i. That will eventually impede the development and growth of individual i. Therefore, the cooperative degree of individual i is defined to be “0.”Case 4.
Considero
b
i
∈
ob
¯
,
o
b
max
⁡, o
b
j
∈
ob
¯
,
o
b
max
⁡. In this case, the two individuals have the similar low levels of growing, and a long period of time is needed for the development of the individuals. But they have the “matching” abilities and can cooperate with each other. Therefore, the cooperative degree of individual i is set to be “1.” The four cases are summarized in Table 1.Table 1
The cooperation degree of individuali (with individual j) Score
i
j.
Condition
Score
o
b
i
∈
[
ob
min
⁡
,
ob
¯
], ob
j
∈
[
ob
min
⁡
,
ob
¯
]
4
ob
i
∈
[
ob
min
⁡
,
ob
¯
], ob
j
∈
[
ob
¯
,
ob
max
⁡
]
2
ob
i
∈
[
ob
¯
,
ob
max
⁡
], ob
j
∈
[
ob
min
⁡
,
ob
¯
]
0
ob
i
∈
[
ob
¯
,
ob
max
⁡
], ob
j
∈
[
ob
¯
,
ob
max
⁡
]
1


## 3.2.2. Cooperative Degree of Individual

Given two populations A and B, denote byc
Fitnes
s
i the cooperative degree of individual i in A. Firstly choose some individuals from B to form a delegation X, and X can be chosen by ways including random selection, greedy selection, or tournament. Then let individual i cooperate with each individual in X with the expert evaluation method. The fitness of i is calculated as follows:(28)
c
Fitnes
s
i
=
∑
j
∈
X
Scor
e
i
j
,where Scor
e
i
j denotes the cooperative degree of individual i by cooperating with individual j.


## 3.2.3. Cooperative Degree of a Population

The cooperative degree of a population stands for the cooperative level of the population, and it can be calculated as follows:(29)
a
1
=
∑
i
∈
Pop
1
c
Fitnes
s
i
∑
i
∈
Pop
2
c
Fitnes
s
i
,
a
2
=
1
a
1
,where Pop1 (Pop2) is the delegation of population A (B). If a
1
>
1, A is more cooperative and obtains more advantage or nutrition from B. Conversely, if a
2
>
1, B is more cooperative than A. The two parameters a
1 and a
2 show the mutual influence degree on population sizes.


## 3.2.4. Adjustment of Population Sizes

In our mutualism population growth model, assume(30)
r
1
=
K
1
-
x
(
t
)
s
1
,
r
2
=
K
2
-
y
(
t
)
s
2
,where s
1
,
s
2 are constants. Based on (26) and (27), at each generation t, the population sizes of A and B will be modified by Δ
x and Δ
y, respectively, which are given as follows: (31)
Δ
x
=
K
1
-
x
s
1
1
-
x
K
2
+
a
1
y
K
1
x
,
Δ
y
=
K
2
-
y
s
2
1
+
a
2
x
K
2
-
y
K
1
y
.The newly added individuals are generated randomly, and the deleted individuals are the ones with the smallest cooperative degree values.If the cooperative degree of a population is high, its population size will increase. Thus the growth of a population size mostly depends on its ability to coordinate. If the ability is strong, the size will become bigger and bigger till the maximum value. Otherwise, the size will become small and eventually extinct. In order to survive, each population will try its best to gain good genes and cooperate with others. In this way, the quality of solution is improved and the performance of algorithm is enhanced.


## 3.3. Framework of the Mutualism Strategy for Population Growth

A description of the mutualism strategy for population growth is given in Figure2.Figure 2
Description of the mutualism strategy for population growth.


## 4. The Mutualism Quantum Genetic Algorithm

The genetic algorithm (GA) proposed by Holland has achieved a great success for solving complex combinational problem. To improve GA, Han and Kim [31] recently introduced a quantum-inspired genetic algorithm, which attracted many researchers’ attention. However, the experiment results are still not satisfactory for large-scale problems. In this section, a method named mutualism quantum genetic algorithm (MQGA) is introduced to solve the integrated supply chain problem including the pickup of materials, the flow shop scheduling, and the delivery of finished jobs. The symbiosis evolutionary strategy is embedded in the quantum genetic algorithm (QGA).


### 4.1. HybridQ-Bit Coding

The encoding in [29] is adopted in this paper, as it has obtained the success on the coordinated scheduling problem with batching machine scheduling and two-stage transportation. Their encoding divided a chromosome into two segments, the job sequence and the job-to-batch assignment, which can represent all possible solutions. Besides, the Q-bit representation is combined in our algorithm for its advantage of linear superposition of solutions and increase of diversity of genes. The details of the encoding are illustrated as follows.Firstly, based on the concept and principles of quantum computing, the smallest unit of information stored in a two-state quantum computer is called a Q-bit. A Q-bit may be in the “1” state, “0” state, or any superposition of the two. A Q-bit can be represented as ψ
=
α
0
+
β
1, where α and β are complex numbers that specify the probability amplitudes of the corresponding states. α
2 gives the probability that the Q-bit will be found in the “0” state and β
2 gives the probability that the Q-bit will be found in the “1” state. Normalization of the state to unity guarantees α
2
+
β
2
=
1. Therefore, in MQGA, a population P
Q
t
=
P
1
t
,
…
,
P
N
t in Q-bit representation is initially randomly generated, and an individual is defined as a string of m
Q-bits α
1
β
1
α
2
β
2
…
…
α
m
β
m, where α
i
2
+
β
i
2
=
1, i
=
1,2
,
…
,
m.Secondly, the individual is converted into binary representation in the following way. For each Q-bit, let η
i be a random number generated in [
0,1
]. If α
i
2
>
η
i, the Q-bit is converted to be 1; otherwise it is 0.Thirdly, the Q-bit representation needs to be changed into two segments, the job sequence and the job-to-batch assignment, which include the comprehensive information of production and delivery. Thus we give a converting rule to form three parts, flow shop part (the code size is lo
g
2
n
+
1
×
n), material pickup part (the code size is 1
×
n), and product delivery part (the code size is 1
×
n), in which we set the last Q-bit of part 2 and part 3 as “1” to ensure the last job belongs to the last batch.For example, given a 3-job, 3-machine problem, let 15Q-bits represent a code, and every three Q-bits form one group. So a code contains five groups, which are then separated into 3 parts. The first three groups are defined as Part 1, the fourth group is Part 2, and the fifth group is Part 3.For Part 1, suppose the binary representation is[
0
1
1
∣
1
0
1
∣
1
1
0
], which is already converted from Q-bit representation; then get [
3
5
6
]. This job permutation is regarded as a random key representation. If values of two elements in the representation are different, let the smaller one denote the job with smaller index; otherwise, let the first one denote the job with smaller index. So, the above random key representation is corresponding to job permutation [
1
2
3
].For Part 2 and Part 3, suppose the binary representations are[
1
0
1
] and [
0
0
1
], which represent the job-to-batch assignment of the conveyor and the truck, respectively. The final code is the combination of Part 1, Part 2, and Part 3. An example is given in Figure 3.Figure 3
The encoding and decoding mode ofQ-bit.The final code for another problem with 8 jobs and 2 machines is presented in (32) which is the individual structure of the MQGA, where the first line represents a permutation of 8 jobs and the second (third) line denotes the job-to-batch assignment in the conveyor (truck): (32)
5
8
2
1
3
6
4
7
1
0
1
0
0
0
1
1
0
0
1
0
1
0
0
1
.The rule is that the jobs in the positions after each “1” to the next “1” belong to a new batch. Therefore, in the pickup transportation, job 5 is in the first batch, job 8 and job 2 belong to the second batch, job 1, job 3, job 6, and job 4 are in the third batch, and job 7 is in the last batch. The third line is the delivery batch of truck. Similarly, in the delivery transportation, job 5, job 8, and job 2 are in the same batch, job 1 and job 3 belong to the next batch, and job 6, job 4, and job 7 are in the last batch. Figure4 illustrates the job sequence and the job-to-batch assignment.Figure 4
Job sequence and the job-to-batch assignment.Assume the processing times of the 8 jobs on the two machines are given in Table2. The pickup time and returning time of the conveyor are 4 time units and 2 time units, respectively. The delivery time and returning time of truck are 8 time units and 4 time units, respectively. Figure 5 gives the visual assignment of jobs to be transported and produced.Table 2
The processing times of the 8 jobs on 2 machines.
Processing time
J
1
J
2
J
3
J
4
J
5
J
6
J
7
J
8
M
1
3
1
2
1
2
3
3
4
M
2
4
1
2
3
4
3
1
2Figure 5
The decoding process and the Gantt chart.


### 4.2. Further Processing Mechanisms Solutions

#### 4.2.1. The Speeding-Up Method of Inserting Neighborhood

When applying local search techniques to search better solutions, the neighborhood structure is very important, since it will directly affect the results. For the flow shop problem, the results in [32–34] show that the inserting operation is better than exchanging when searching in the local neighborhood. Therefore, this paper employs the inserting operation during the process of search in the neighborhood. The inserting neighborhood of a permutation is a set of solutions, that is, the permutations generated by deleting a job from the original position and inserting the job into a new position.Denote byπ
(
i
) the job in the ith position in the permutation π, and define ω
π
(
i
) as the subpermutation after removing job π
(
i
). Assuming job π
(
i
) is inserted in the hth position of ω
i, let ω
π
(
i
)
(
h
) be the new permutation and σ the position where job π
(
i
) is inserted to achieve the minimal makespan; that is, (33)
ω
π
i
σ
=
min
⁡
1
≤
h
≤
n
ω
π
i
h
.Take the two-machine flow shop problem as an example. LetC
max
⁡
π
;
M
1
,
M
2 be the makespan for a given permutation π, and let T
M
1
π
j and T
M
2
π
j, respectively, be the earliest possible completion times of job π
(
j
) on machines M
1 and M
2 in theforward channel. In addition, define T
M
1
′
π
j, T
M
2
′
π
j, respectively, denoting the earliest possible completion times of job π
(
j
) on machines M
1 and M
2 in thebackward channel. The iterative formulas are given as follows:(34)
T
M
1
π
j
=
T
M
1
π
j
-
1
+
p
1
π
(
j
)
,
j
=
1,2
,
…
,
n
,
T
M
2
π
j
=
max
⁡
T
M
1
π
j
,
T
M
2
π
j
-
1
+
p
2
π
(
j
)
,
j
=
1,2
,
…
,
n
,
T
M
1
π
0
=
T
M
2
π
0
=
0
,
T
M
2
′
π
j
=
T
M
2
′
π
j
+
1
+
p
2
π
(
j
)
,
j
=
n
-
1
,
n
-
2
,
…
,
1
,
T
M
1
′
π
j
=
max
⁡
T
M
2
′
π
j
,
T
M
1
′
π
j
+
1
+
p
1
π
(
j
)
,
j
=
n
-
1
,
n
-
2
,
…
,
1
,
T
M
1
′
π
n
=
T
M
2
′
π
n
=
0
,
C
max
⁡
π
=
T
M
2
π
n
=
T
M
1
′
π
1
,
C
max
⁡
ω
π
(
i
)
h
;
M
1
,
M
2
=
max
⁡
T
M
1
ω
h
-
1
+
p
1
π
i
+
T
M
1
′
ω
h
,
T
M
2
ω
h
-
1
+
p
2
π
i
+
T
M
2
′
ω
h
,
T
M
1
ω
h
-
1
+
p
1
π
i
+
p
2
π
i
+
T
M
2
′
ω
h
.For the flow shop problem with more than two machines, the details of obtaining the local optimal solution are given in Procedure1.Procedure 1:The speeding-up method of inserting neighborhood.
Begin
i
←
1
While (i
≤
n) do
Begin
Obtain the subpermutationω
π
(
i
) by removing job π
(
i
) in the original permutation π.
For   k
=
1,2
,
…
,
m
-
1
Calculate T
M
k
(
ω
π
(
i
)
(
j
)
), T
M
k
+
1
(
ω
π
(
i
)
(
j
)
), T
M
k
+
1
′
(
ω
π
(
i
)
(
j
)
) and T
M
k
′
(
ω
π
(
i
)
(
j
)
).
End
For   h
=
1,2
,
…
,
n and h
≠
i
For   k
=
1,2
,
…
,
m
-
1
Calculate C
max
⁡
(
ω
π
(
i
)
(
h
)
;
M
k
,
M
k
+
1
).
End
CalculateC
max
⁡
(
ω
π
(
i
)
(
h
)
)
=
C
max
⁡
(
ω
π
(
i
)
(
h
)
;
M
m
-
1
,
M
m
).
End
Denote byσ the position where job π
(
i
) is inserted to get the minimal objective value.
i
=
i
+
1;
End
End
Output the local optimal individual.


#### 4.2.2. A Modified Mechanism of the Solution

Because the total physical space of all jobs in a batch may exceed the capacity of the vehicle, our encoding could not ensure a feasible solution, and a modified mechanism is needed to update the solution, which is an important step, not only in encoding process, but also in the steps of crossover and mutation, as an unfeasible solution may affect the search process in local neighborhood.Let TVol be the capacity of the vehicle andVol
B
j the total space of all jobs in batch B
j. For simplicity, denote by J
j
,
1
,
J
j
,
2
,
…
J
j
,
b
j the jobs in batch B
j, where b
j is the number of jobs in B
j. The modified mechanism is described in Procedure 2.Procedure 2:The modified mechanism of a solution.
i
=
1
For batch B
j
IfVol
(
B
j
)
>
TVol,
Find the jobJ
j
,
b in B
j such that
Vol
(
J
j
,
1
+
J
j
,
2
+
⋯
+
J
j
,
b
)
≤
TVol and Vol
(
J
j
,
1
+
J
j
,
2
+
⋯
+
J
j
,
s
+
J
j
,
b
+
1
)
>
TVol.
LetB
j
=
{
J
j
,
1
,
…
,
J
j
,
b
} and B
j
+
1
=
B
j
+
1
∪
{
J
j
,
b
+
1
,
…
,
J
j
,
b
j
}. Update the index of jobs in B
j
+
1.
i
=
i
+
1.
End


### 4.3. Other Operations

The operations introduced below are performed at each generation of the quantum genetic algorithm (QGA). The fitness of each individual can be calculated by the mutualism strategy in Section3.2. In order to ensure the individual with good genes can be chosen in the offspring, the roulette selection is applied here, which could maintain the diversity of genes. In the crossover operation, two decoded individuals are selected as parents to operate single-point crossover. See Figure 6 for the illustration of this operation, in which the corresponding genes of Q-bit parents are swapped to produce two new offspring.Figure 6
Crossover operation of MQGA.In the mutation operation, a NOT Gate is used as the mutation operator. Firstly, select individuals with mutation probabilityP
m; then randomly generate a mutation position. Quantum rotation gate is used to maintain diversity of genes, which could change the probability amplitude of quantum states. See more details in our previous work [35]. Besides, a catastrophe operator is used to avoid premature convergence. If the solution does not change in two consecutive generations, it could be regarded to be trapped in local optimal solution, and the local best solution will be reserved while others will be replaced by solutions generated randomly.


### 4.4. The Main Procedure of MQGA

In this subsection, the details of the MQGA are introduced in Algorithm1.Algorithm 1:MQGA (mutualism quantum genetic algorithm).
Begin
t
←
0
Step  1. Set parameters: crossover probability P
c, mutation probability P
m, iterative
generation GN, two subpopulations sizesPS
1, PS
2.
Step  2. Initialize subpopulations in hybrid Q-bit representation.
While (t
≠
GN) do
Begin
Step 3. If (subpopulation is not empty)
Apply mutualism strategy, rescale the size of subpopulations according to population growth model.
End.
Step  4. Perform selection, crossover, mutation operations in each subpopulation;
Step  5. If (catastrophe condition is satisfied)
Perform catastrophe operation.
Else
Apply quantum rotation operation.
End
Step  6. Calculate objective value and fitness value of each individual in the current generation,
and update the best solution if possible.
t
←
t
+
1
End
End
Output the global best result.


## 4.1. HybridQ-Bit Coding

The encoding in [29] is adopted in this paper, as it has obtained the success on the coordinated scheduling problem with batching machine scheduling and two-stage transportation. Their encoding divided a chromosome into two segments, the job sequence and the job-to-batch assignment, which can represent all possible solutions. Besides, the Q-bit representation is combined in our algorithm for its advantage of linear superposition of solutions and increase of diversity of genes. The details of the encoding are illustrated as follows.Firstly, based on the concept and principles of quantum computing, the smallest unit of information stored in a two-state quantum computer is called a Q-bit. A Q-bit may be in the “1” state, “0” state, or any superposition of the two. A Q-bit can be represented as ψ
=
α
0
+
β
1, where α and β are complex numbers that specify the probability amplitudes of the corresponding states. α
2 gives the probability that the Q-bit will be found in the “0” state and β
2 gives the probability that the Q-bit will be found in the “1” state. Normalization of the state to unity guarantees α
2
+
β
2
=
1. Therefore, in MQGA, a population P
Q
t
=
P
1
t
,
…
,
P
N
t in Q-bit representation is initially randomly generated, and an individual is defined as a string of m
Q-bits α
1
β
1
α
2
β
2
…
…
α
m
β
m, where α
i
2
+
β
i
2
=
1, i
=
1,2
,
…
,
m.Secondly, the individual is converted into binary representation in the following way. For each Q-bit, let η
i be a random number generated in [
0,1
]. If α
i
2
>
η
i, the Q-bit is converted to be 1; otherwise it is 0.Thirdly, the Q-bit representation needs to be changed into two segments, the job sequence and the job-to-batch assignment, which include the comprehensive information of production and delivery. Thus we give a converting rule to form three parts, flow shop part (the code size is lo
g
2
n
+
1
×
n), material pickup part (the code size is 1
×
n), and product delivery part (the code size is 1
×
n), in which we set the last Q-bit of part 2 and part 3 as “1” to ensure the last job belongs to the last batch.For example, given a 3-job, 3-machine problem, let 15Q-bits represent a code, and every three Q-bits form one group. So a code contains five groups, which are then separated into 3 parts. The first three groups are defined as Part 1, the fourth group is Part 2, and the fifth group is Part 3.For Part 1, suppose the binary representation is[
0
1
1
∣
1
0
1
∣
1
1
0
], which is already converted from Q-bit representation; then get [
3
5
6
]. This job permutation is regarded as a random key representation. If values of two elements in the representation are different, let the smaller one denote the job with smaller index; otherwise, let the first one denote the job with smaller index. So, the above random key representation is corresponding to job permutation [
1
2
3
].For Part 2 and Part 3, suppose the binary representations are[
1
0
1
] and [
0
0
1
], which represent the job-to-batch assignment of the conveyor and the truck, respectively. The final code is the combination of Part 1, Part 2, and Part 3. An example is given in Figure 3.Figure 3
The encoding and decoding mode ofQ-bit.The final code for another problem with 8 jobs and 2 machines is presented in (32) which is the individual structure of the MQGA, where the first line represents a permutation of 8 jobs and the second (third) line denotes the job-to-batch assignment in the conveyor (truck): (32)
5
8
2
1
3
6
4
7
1
0
1
0
0
0
1
1
0
0
1
0
1
0
0
1
.The rule is that the jobs in the positions after each “1” to the next “1” belong to a new batch. Therefore, in the pickup transportation, job 5 is in the first batch, job 8 and job 2 belong to the second batch, job 1, job 3, job 6, and job 4 are in the third batch, and job 7 is in the last batch. The third line is the delivery batch of truck. Similarly, in the delivery transportation, job 5, job 8, and job 2 are in the same batch, job 1 and job 3 belong to the next batch, and job 6, job 4, and job 7 are in the last batch. Figure4 illustrates the job sequence and the job-to-batch assignment.Figure 4
Job sequence and the job-to-batch assignment.Assume the processing times of the 8 jobs on the two machines are given in Table2. The pickup time and returning time of the conveyor are 4 time units and 2 time units, respectively. The delivery time and returning time of truck are 8 time units and 4 time units, respectively. Figure 5 gives the visual assignment of jobs to be transported and produced.Table 2
The processing times of the 8 jobs on 2 machines.
Processing time
J
1
J
2
J
3
J
4
J
5
J
6
J
7
J
8
M
1
3
1
2
1
2
3
3
4
M
2
4
1
2
3
4
3
1
2Figure 5
The decoding process and the Gantt chart.


## 4.2. Further Processing Mechanisms Solutions

### 4.2.1. The Speeding-Up Method of Inserting Neighborhood

When applying local search techniques to search better solutions, the neighborhood structure is very important, since it will directly affect the results. For the flow shop problem, the results in [32–34] show that the inserting operation is better than exchanging when searching in the local neighborhood. Therefore, this paper employs the inserting operation during the process of search in the neighborhood. The inserting neighborhood of a permutation is a set of solutions, that is, the permutations generated by deleting a job from the original position and inserting the job into a new position.Denote byπ
(
i
) the job in the ith position in the permutation π, and define ω
π
(
i
) as the subpermutation after removing job π
(
i
). Assuming job π
(
i
) is inserted in the hth position of ω
i, let ω
π
(
i
)
(
h
) be the new permutation and σ the position where job π
(
i
) is inserted to achieve the minimal makespan; that is, (33)
ω
π
i
σ
=
min
⁡
1
≤
h
≤
n
ω
π
i
h
.Take the two-machine flow shop problem as an example. LetC
max
⁡
π
;
M
1
,
M
2 be the makespan for a given permutation π, and let T
M
1
π
j and T
M
2
π
j, respectively, be the earliest possible completion times of job π
(
j
) on machines M
1 and M
2 in theforward channel. In addition, define T
M
1
′
π
j, T
M
2
′
π
j, respectively, denoting the earliest possible completion times of job π
(
j
) on machines M
1 and M
2 in thebackward channel. The iterative formulas are given as follows:(34)
T
M
1
π
j
=
T
M
1
π
j
-
1
+
p
1
π
(
j
)
,
j
=
1,2
,
…
,
n
,
T
M
2
π
j
=
max
⁡
T
M
1
π
j
,
T
M
2
π
j
-
1
+
p
2
π
(
j
)
,
j
=
1,2
,
…
,
n
,
T
M
1
π
0
=
T
M
2
π
0
=
0
,
T
M
2
′
π
j
=
T
M
2
′
π
j
+
1
+
p
2
π
(
j
)
,
j
=
n
-
1
,
n
-
2
,
…
,
1
,
T
M
1
′
π
j
=
max
⁡
T
M
2
′
π
j
,
T
M
1
′
π
j
+
1
+
p
1
π
(
j
)
,
j
=
n
-
1
,
n
-
2
,
…
,
1
,
T
M
1
′
π
n
=
T
M
2
′
π
n
=
0
,
C
max
⁡
π
=
T
M
2
π
n
=
T
M
1
′
π
1
,
C
max
⁡
ω
π
(
i
)
h
;
M
1
,
M
2
=
max
⁡
T
M
1
ω
h
-
1
+
p
1
π
i
+
T
M
1
′
ω
h
,
T
M
2
ω
h
-
1
+
p
2
π
i
+
T
M
2
′
ω
h
,
T
M
1
ω
h
-
1
+
p
1
π
i
+
p
2
π
i
+
T
M
2
′
ω
h
.For the flow shop problem with more than two machines, the details of obtaining the local optimal solution are given in Procedure1.Procedure 1:The speeding-up method of inserting neighborhood.
Begin
i
←
1
While (i
≤
n) do
Begin
Obtain the subpermutationω
π
(
i
) by removing job π
(
i
) in the original permutation π.
For   k
=
1,2
,
…
,
m
-
1
Calculate T
M
k
(
ω
π
(
i
)
(
j
)
), T
M
k
+
1
(
ω
π
(
i
)
(
j
)
), T
M
k
+
1
′
(
ω
π
(
i
)
(
j
)
) and T
M
k
′
(
ω
π
(
i
)
(
j
)
).
End
For   h
=
1,2
,
…
,
n and h
≠
i
For   k
=
1,2
,
…
,
m
-
1
Calculate C
max
⁡
(
ω
π
(
i
)
(
h
)
;
M
k
,
M
k
+
1
).
End
CalculateC
max
⁡
(
ω
π
(
i
)
(
h
)
)
=
C
max
⁡
(
ω
π
(
i
)
(
h
)
;
M
m
-
1
,
M
m
).
End
Denote byσ the position where job π
(
i
) is inserted to get the minimal objective value.
i
=
i
+
1;
End
End
Output the local optimal individual.


### 4.2.2. A Modified Mechanism of the Solution

Because the total physical space of all jobs in a batch may exceed the capacity of the vehicle, our encoding could not ensure a feasible solution, and a modified mechanism is needed to update the solution, which is an important step, not only in encoding process, but also in the steps of crossover and mutation, as an unfeasible solution may affect the search process in local neighborhood.Let TVol be the capacity of the vehicle andVol
B
j the total space of all jobs in batch B
j. For simplicity, denote by J
j
,
1
,
J
j
,
2
,
…
J
j
,
b
j the jobs in batch B
j, where b
j is the number of jobs in B
j. The modified mechanism is described in Procedure 2.Procedure 2:The modified mechanism of a solution.
i
=
1
For batch B
j
IfVol
(
B
j
)
>
TVol,
Find the jobJ
j
,
b in B
j such that
Vol
(
J
j
,
1
+
J
j
,
2
+
⋯
+
J
j
,
b
)
≤
TVol and Vol
(
J
j
,
1
+
J
j
,
2
+
⋯
+
J
j
,
s
+
J
j
,
b
+
1
)
>
TVol.
LetB
j
=
{
J
j
,
1
,
…
,
J
j
,
b
} and B
j
+
1
=
B
j
+
1
∪
{
J
j
,
b
+
1
,
…
,
J
j
,
b
j
}. Update the index of jobs in B
j
+
1.
i
=
i
+
1.
End


## 4.2.1. The Speeding-Up Method of Inserting Neighborhood

When applying local search techniques to search better solutions, the neighborhood structure is very important, since it will directly affect the results. For the flow shop problem, the results in [32–34] show that the inserting operation is better than exchanging when searching in the local neighborhood. Therefore, this paper employs the inserting operation during the process of search in the neighborhood. The inserting neighborhood of a permutation is a set of solutions, that is, the permutations generated by deleting a job from the original position and inserting the job into a new position.Denote byπ
(
i
) the job in the ith position in the permutation π, and define ω
π
(
i
) as the subpermutation after removing job π
(
i
). Assuming job π
(
i
) is inserted in the hth position of ω
i, let ω
π
(
i
)
(
h
) be the new permutation and σ the position where job π
(
i
) is inserted to achieve the minimal makespan; that is, (33)
ω
π
i
σ
=
min
⁡
1
≤
h
≤
n
ω
π
i
h
.Take the two-machine flow shop problem as an example. LetC
max
⁡
π
;
M
1
,
M
2 be the makespan for a given permutation π, and let T
M
1
π
j and T
M
2
π
j, respectively, be the earliest possible completion times of job π
(
j
) on machines M
1 and M
2 in theforward channel. In addition, define T
M
1
′
π
j, T
M
2
′
π
j, respectively, denoting the earliest possible completion times of job π
(
j
) on machines M
1 and M
2 in thebackward channel. The iterative formulas are given as follows:(34)
T
M
1
π
j
=
T
M
1
π
j
-
1
+
p
1
π
(
j
)
,
j
=
1,2
,
…
,
n
,
T
M
2
π
j
=
max
⁡
T
M
1
π
j
,
T
M
2
π
j
-
1
+
p
2
π
(
j
)
,
j
=
1,2
,
…
,
n
,
T
M
1
π
0
=
T
M
2
π
0
=
0
,
T
M
2
′
π
j
=
T
M
2
′
π
j
+
1
+
p
2
π
(
j
)
,
j
=
n
-
1
,
n
-
2
,
…
,
1
,
T
M
1
′
π
j
=
max
⁡
T
M
2
′
π
j
,
T
M
1
′
π
j
+
1
+
p
1
π
(
j
)
,
j
=
n
-
1
,
n
-
2
,
…
,
1
,
T
M
1
′
π
n
=
T
M
2
′
π
n
=
0
,
C
max
⁡
π
=
T
M
2
π
n
=
T
M
1
′
π
1
,
C
max
⁡
ω
π
(
i
)
h
;
M
1
,
M
2
=
max
⁡
T
M
1
ω
h
-
1
+
p
1
π
i
+
T
M
1
′
ω
h
,
T
M
2
ω
h
-
1
+
p
2
π
i
+
T
M
2
′
ω
h
,
T
M
1
ω
h
-
1
+
p
1
π
i
+
p
2
π
i
+
T
M
2
′
ω
h
.For the flow shop problem with more than two machines, the details of obtaining the local optimal solution are given in Procedure1.Procedure 1:The speeding-up method of inserting neighborhood.
Begin
i
←
1
While (i
≤
n) do
Begin
Obtain the subpermutationω
π
(
i
) by removing job π
(
i
) in the original permutation π.
For   k
=
1,2
,
…
,
m
-
1
Calculate T
M
k
(
ω
π
(
i
)
(
j
)
), T
M
k
+
1
(
ω
π
(
i
)
(
j
)
), T
M
k
+
1
′
(
ω
π
(
i
)
(
j
)
) and T
M
k
′
(
ω
π
(
i
)
(
j
)
).
End
For   h
=
1,2
,
…
,
n and h
≠
i
For   k
=
1,2
,
…
,
m
-
1
Calculate C
max
⁡
(
ω
π
(
i
)
(
h
)
;
M
k
,
M
k
+
1
).
End
CalculateC
max
⁡
(
ω
π
(
i
)
(
h
)
)
=
C
max
⁡
(
ω
π
(
i
)
(
h
)
;
M
m
-
1
,
M
m
).
End
Denote byσ the position where job π
(
i
) is inserted to get the minimal objective value.
i
=
i
+
1;
End
End
Output the local optimal individual.


## 4.2.2. A Modified Mechanism of the Solution

Because the total physical space of all jobs in a batch may exceed the capacity of the vehicle, our encoding could not ensure a feasible solution, and a modified mechanism is needed to update the solution, which is an important step, not only in encoding process, but also in the steps of crossover and mutation, as an unfeasible solution may affect the search process in local neighborhood.Let TVol be the capacity of the vehicle andVol
B
j the total space of all jobs in batch B
j. For simplicity, denote by J
j
,
1
,
J
j
,
2
,
…
J
j
,
b
j the jobs in batch B
j, where b
j is the number of jobs in B
j. The modified mechanism is described in Procedure 2.Procedure 2:The modified mechanism of a solution.
i
=
1
For batch B
j
IfVol
(
B
j
)
>
TVol,
Find the jobJ
j
,
b in B
j such that
Vol
(
J
j
,
1
+
J
j
,
2
+
⋯
+
J
j
,
b
)
≤
TVol and Vol
(
J
j
,
1
+
J
j
,
2
+
⋯
+
J
j
,
s
+
J
j
,
b
+
1
)
>
TVol.
LetB
j
=
{
J
j
,
1
,
…
,
J
j
,
b
} and B
j
+
1
=
B
j
+
1
∪
{
J
j
,
b
+
1
,
…
,
J
j
,
b
j
}. Update the index of jobs in B
j
+
1.
i
=
i
+
1.
End


## 4.3. Other Operations

The operations introduced below are performed at each generation of the quantum genetic algorithm (QGA). The fitness of each individual can be calculated by the mutualism strategy in Section3.2. In order to ensure the individual with good genes can be chosen in the offspring, the roulette selection is applied here, which could maintain the diversity of genes. In the crossover operation, two decoded individuals are selected as parents to operate single-point crossover. See Figure 6 for the illustration of this operation, in which the corresponding genes of Q-bit parents are swapped to produce two new offspring.Figure 6
Crossover operation of MQGA.In the mutation operation, a NOT Gate is used as the mutation operator. Firstly, select individuals with mutation probabilityP
m; then randomly generate a mutation position. Quantum rotation gate is used to maintain diversity of genes, which could change the probability amplitude of quantum states. See more details in our previous work [35]. Besides, a catastrophe operator is used to avoid premature convergence. If the solution does not change in two consecutive generations, it could be regarded to be trapped in local optimal solution, and the local best solution will be reserved while others will be replaced by solutions generated randomly.


## 4.4. The Main Procedure of MQGA

In this subsection, the details of the MQGA are introduced in Algorithm1.Algorithm 1:MQGA (mutualism quantum genetic algorithm).
Begin
t
←
0
Step  1. Set parameters: crossover probability P
c, mutation probability P
m, iterative
generation GN, two subpopulations sizesPS
1, PS
2.
Step  2. Initialize subpopulations in hybrid Q-bit representation.
While (t
≠
GN) do
Begin
Step 3. If (subpopulation is not empty)
Apply mutualism strategy, rescale the size of subpopulations according to population growth model.
End.
Step  4. Perform selection, crossover, mutation operations in each subpopulation;
Step  5. If (catastrophe condition is satisfied)
Perform catastrophe operation.
Else
Apply quantum rotation operation.
End
Step  6. Calculate objective value and fitness value of each individual in the current generation,
and update the best solution if possible.
t
←
t
+
1
End
End
Output the global best result.


## 5. Experiment

In this paper, all algorithms are programmed with MATLAB language, and all the computations are conducted on a Pentium PC 1.66 GHZ with 512 MB memory. In order to evaluate the performance of MQGA, an extensive set of instances with different characteristics are generated based on the flow shop benchmark problems, including FT benchmark and ABZ benchmark. In order to test the relative large-scale problem, three test instances are designed, where there are 50 jobs and 5 machines, 50 jobs and 10 machines, and 100 jobs and 5 machines, seen in the Appendix. Table3 shows the transporting times of conveyor and truck in different problems. Table 4 shows the physical spaces of the 100 jobs. These test instances are named N(FT06), N(FT10), N(FT20), N(ABZ3), N(ABZ6), N(ABZ7), N(ABZ8), N(ABZ9), N(TA1), N(TA2), and N(TA3).Table 3
Parameters of the conveyor and truck.
Problem
Parameter
C1
Tr1
R1
C2
Tr2
R2
N(FT06)
6
10
3
10
20
10
N(FT10)
6
40
20
10
100
30
N(FT20)
6
20
10
10
60
30
N(ABZ3)
6
100
100
10
300
100
N(ABZ6)
6
100
100
10
300
100
N(ABZ7)
6
30
30
10
100
100
N(ABZ8)
6
30
30
10
100
100
N(ABZ9)
6
30
30
10
100
100
N(TA1)
6
30
30
10
100
100
N(TA2)
6
100
30
10
300
200
N(TA3)
6
100
30
10
300
200
Note:C1 and C2 are the capacities of conveyor and truck, respectively. Tr1 and Tr2 are the pickup and delivery times. R1 and R2 are the returning times of conveyor and truck.Table 4
Physical space of jobs.
Jobs 1–20
2 6 1 0.3 3 1.3 2 3.3 1.2 4 2.8 3.2 0.7 0.1 0.9 2 6 1 0.3 3
Jobs 21–40
1.3 2 3.3 1.2 4 2.8 3.2 0.7 0.1 0.9 2 6 1 0.3 3 1.3 2 3.3 1.2 4
Jobs 41–60
2.8 3.2 0.7 0.1 0.9 3 4 3 6 1 2 6 1 0.3 3 1.3 2 3.3 1.2 4
Jobs 61–80
2.8 3.2 0.7 0.1 0.9 2 6 1 0.3 3 1.3 2 3.3 1.2 4 2.8 3.2 0.7 0.1 0.9
Jobs 81–100
2 6 1 0.3 3 1.3 2 3.3 1.2 4 2.8 3.2 0.7 0.1 0.9 3 4 3 6 1


### 5.1. Parameter Settings

Two parameters including crossover rateP
c and mutation rate P
m are discussed in this section. Assume the two rates are set with the following levels:(35)
P
m
:
0.2,0.4,0.6,0.8,1
;
P
c
:
0.2,0.4,0.6,0.8,1
.In this way, there are a total of 25 combinations. For each combination, MGA [13] is used to test the 11 testing instances in Table 3. After the experiment, the data is transformed to be the relative percentage deviation (RPD), which is calculated as follows:(36)
RPD
=
W
T
A
-
W
T
ref
W
T
ref
×
100
%
,where W
T
A is the corresponding objective value of the solution achieved by proposed algorithm and W
T
ref is the minimum objective value under multiple simulation of the same test problem. In addition, for a combination of P
c and P
m values, the analysis of variance (AOV) model is employed to eliminate the different RPD values caused by the block (each test problem is regarded as a block). The main experiment results of the AOV are shown in Table 5.Table 5
AOV results for the experiment on tuning the parameters.
Source
Sum of squares
Df
Mean square
F value
A:P
c
0.0083
4
0.0021
3.3334
B:P
m
0.0072
4
0.0018
4.3732
Interactions AB
0.0438
16
0.0029
7.3679According to the theory of statistical analysis, the greater theF value is, the more significantly the corresponding factors affect. Table 5 shows the interaction between the two parameters has a significant effect on the performance of the algorithm. Based on the experiment, we find the best value of P
c is 0.8 and the best value of P
m is 0.2, since they could find the most times of optimal objective value.


### 5.2. Simulation Work

#### 5.2.1. The Detailed Simulation Result of N(TA1) Problem

In order to validate the performance of MQGA, MGA [13] and QGA [36] are compared. Firstly problem N(TA1) is chosen as the representative instance to analyze the performance of algorithms. Three combinations of iterative generations and population sizes are 500 × 10 (i.e., the maximum number of iterations is 500, population size is 10, and subpopulation size is 5), 1000 × 10, and 2000 × 20. For each combination, the results are presented in Table 6, and more results are given in Table 7, where “BV” and “WV,” respectively, denote the best and worst objective values in the 10 run times. “AV” is the average objective value.Table 6
Detailed comparison results of MQGA, QGA, and MGA on N(TA1) problem.
RUN
500 × 10
1000 × 10
2000 × 20
MQGA
MGA
QGA
MQGA
MGA
QGA
MQGA
MGA
QGA
1
3166
3329
3199
3139
3326
3194
3152
3308
2
3182
3336
3204
3166
3328
3193
3147
3327
3
3165
3382
3229
3154
3351
3209
3154
3345
4
3172
3277
3218
3162
3266
3210
3147
3210
5
3164
3307
3220
3151
3323
3208
3148
3285
6
3179
3330
3215
3165
3319
3170
3155
3324
7
3164
3267
3222
3144
3253
3222
3128
3253
8
3145
3371
3221
3138
3218
3191
3167
3313
9
3154
3297
3226
3163
3293
3221
3140
3305
10
3163
3298
3206
3140
3285
3145
3135
3276Table 7
Comparison results of MQGA, QGA, and MGA on N(TA1) problem.
Size
MQGA
MGA
QGA
BV
WV
AV
BV
WV
AV
BV
WV
AV
500 × 10
3155
3192
3165.4
3267
3382
3325.7
3199
3229
3216
1000 × 10
3138
3166
3150.2
3218
3351
3296.2
3145
3222
3196.3
2000 × 20
3128
3167
3147.3
3210
3345
3294.6
3131
3226
3186.3Known from the results in Table6, for each calculation, MQGA always obtains a better solution than that obtained by MGA and QGA. When the iteration is 500, the range of objective values of MQGA is (3155, 3192) and (3267, 3382) and (3199, 3229) for MGA and QGA. When the iteration is 2000, the ranges of objective values are, respectively, (3128, 3167), (3210, 3345), and (3131, 3226). Therefore, the search ability of MQGA increases faster than that of MGA and QGA as the iteration increases, and MGA and QGA are more likely to be trapped in local optimal solutions.The convergence curves of MQGA, MGA, and QGA obtained for each combination are depicted in Figure7, which indicates the convergence speed of MQGA is faster than that of MGA and QGA, especially for the combination with larger iteration.Figure 7
The convergence curves of MQGA, QGA, and MGA.The solutions distributions of MQGA, QGA, and MGA are given in Figure8, which shows that there is no overlap in the solution space of each algorithm. Moreover, the change in the solutions achieved by MQGA is smaller than that of MGA and QGA, which implies that the quality of solution obtained by MQGA is better than that of MGA and QGA.Figure 8
Distributions of solutions of MQGA, QGA, and MGA.


#### 5.2.2. The Comparison Results of MQGA with MGA and QGA

The maximum number of iterations is 1000, and the size of population is 100 (in MQGA, there are two subpopulations, and the size of each subpopulation is 50). To be fair, each algorithm runs 10 times for each instance. The experiment results are shown in Table8, in which each item is an average value. The “Improved%” column represents the percentage difference between the average objective values obtained by the current algorithm with MGA and is calculated by the following formula:(37)
Improved
%
=
(
Ob
¯
-
MGA
¯
)
MGA
¯
×
100
%
,where Ob
¯ represents the objective value achieved by the current algorithm and MGA
¯ represents the objective value achieved by MGA. Boldface and italic indicate the best results for each problem.Table 8
Computational results of MQGA and QGA, MGA.
Problem
Approach
Ob
¯
Improved%
AG
¯
AT
¯
N(FT06)(6 * 6)
MGA
103.3
—
73.3
13
QGA
103.3
0
33.6
30
MQGA
103
0.48%
9.8
113
N(FT10)(10 * 10)
MGA
1282.2
—
300.8
29
QGA
1280.7
0.12%
97.3
98
MQGA
1273.6
0.31%
76.2
163
N(FT20)(20 * 5)
MGA
1330
—
26.7
23
QGA
1328.6
0.09%
18.3
36
MQGA
1323.1
0.32%
36.7
133
N(ABZ3)(10 * 10)
MGA
2088.6
—
430.8
29
QGA
2079.4
0.44%
146.6
101
MQGA
2071.3
0.83%
83.4
164
N(ABZ6)(10 * 10)
MGA
1883.2
—
363.6
29
QGA
1874.4
0.47%
111.3
101
MQGA
1861.4
1.16%
123.3
164
N(ABZ7)(15 * 20)
MGA
1276.2
—
140.3
68
QGA
1233.8
3.32%
104.3
172
MQGA
1228.2
3.76%
86.8
310
N(ABZ8)(15 * 20)
MGA
1290.3
—
114.4
68
QGA
1277.2
1.03%
17.7
172
MQGA
1268.4
1.71%
19.8
310
N(ABZ9)(15 * 20)
MGA
1273.3
—
143.9
68
QGA
1238.3
1.18%
29.7
171
MQGA
1241.3
2.3%
11.3
310
N(TA1)(50 * 5)
MGA
3286.2
—
71.7
84
QGA
3173.3
3.37%
17
426
MQGA
3134.2
4.02%
21.1
778
N(TA2)(50 * 10)
MGA
3346.6
—
91.2
130
QGA
4917
11.33%
9.1
480
MQGA
4803
13.4%
19.3
814
N(TA3)(100 * 5)
MGA
10706.3
—
63.3
171
QGA
9289.33
13.23%
16.2
632
MQGA
9126.4
14.76%
17.8
963
Note:AG
¯ represents the converged generation and AT
¯ represents the average computation time.In the performance of solution, from Table8, MQGA is superior to MGA and QGA for all the test instances, and QGA always outperforms MGA. The improvement becomes more obvious as the scale of the problem increases. For example, for medium-scale problem N(FT20), the optimal ability of MQGA improves 0.32%. For the large-scale problem N(TA3), the MQGA improves 14.76%. The reason is that when the number of jobs increases the transportation of vehicle will become the bottleneck for MGA and QGA. Therefore, MGA and QGA are not applicable to the relative large-scale problem instances. In addition, in the performance of convergence, MQGA can obtain a better solution than MGA and QGA with less iterative generations. For the relative large-scale problems like N(TA1), N(TA2), and N(TA3), MGA and QGA converge faster to the local optimal solution. However, MQGA could jump out of the local optimum with more iterative generations. Finally, in the performance of time, MQGA needs more computation time than MGA and QGA but obtains an improved solution.


## 5.1. Parameter Settings

Two parameters including crossover rateP
c and mutation rate P
m are discussed in this section. Assume the two rates are set with the following levels:(35)
P
m
:
0.2,0.4,0.6,0.8,1
;
P
c
:
0.2,0.4,0.6,0.8,1
.In this way, there are a total of 25 combinations. For each combination, MGA [13] is used to test the 11 testing instances in Table 3. After the experiment, the data is transformed to be the relative percentage deviation (RPD), which is calculated as follows:(36)
RPD
=
W
T
A
-
W
T
ref
W
T
ref
×
100
%
,where W
T
A is the corresponding objective value of the solution achieved by proposed algorithm and W
T
ref is the minimum objective value under multiple simulation of the same test problem. In addition, for a combination of P
c and P
m values, the analysis of variance (AOV) model is employed to eliminate the different RPD values caused by the block (each test problem is regarded as a block). The main experiment results of the AOV are shown in Table 5.Table 5
AOV results for the experiment on tuning the parameters.
Source
Sum of squares
Df
Mean square
F value
A:P
c
0.0083
4
0.0021
3.3334
B:P
m
0.0072
4
0.0018
4.3732
Interactions AB
0.0438
16
0.0029
7.3679According to the theory of statistical analysis, the greater theF value is, the more significantly the corresponding factors affect. Table 5 shows the interaction between the two parameters has a significant effect on the performance of the algorithm. Based on the experiment, we find the best value of P
c is 0.8 and the best value of P
m is 0.2, since they could find the most times of optimal objective value.


## 5.2. Simulation Work

### 5.2.1. The Detailed Simulation Result of N(TA1) Problem

In order to validate the performance of MQGA, MGA [13] and QGA [36] are compared. Firstly problem N(TA1) is chosen as the representative instance to analyze the performance of algorithms. Three combinations of iterative generations and population sizes are 500 × 10 (i.e., the maximum number of iterations is 500, population size is 10, and subpopulation size is 5), 1000 × 10, and 2000 × 20. For each combination, the results are presented in Table 6, and more results are given in Table 7, where “BV” and “WV,” respectively, denote the best and worst objective values in the 10 run times. “AV” is the average objective value.Table 6
Detailed comparison results of MQGA, QGA, and MGA on N(TA1) problem.
RUN
500 × 10
1000 × 10
2000 × 20
MQGA
MGA
QGA
MQGA
MGA
QGA
MQGA
MGA
QGA
1
3166
3329
3199
3139
3326
3194
3152
3308
2
3182
3336
3204
3166
3328
3193
3147
3327
3
3165
3382
3229
3154
3351
3209
3154
3345
4
3172
3277
3218
3162
3266
3210
3147
3210
5
3164
3307
3220
3151
3323
3208
3148
3285
6
3179
3330
3215
3165
3319
3170
3155
3324
7
3164
3267
3222
3144
3253
3222
3128
3253
8
3145
3371
3221
3138
3218
3191
3167
3313
9
3154
3297
3226
3163
3293
3221
3140
3305
10
3163
3298
3206
3140
3285
3145
3135
3276Table 7
Comparison results of MQGA, QGA, and MGA on N(TA1) problem.
Size
MQGA
MGA
QGA
BV
WV
AV
BV
WV
AV
BV
WV
AV
500 × 10
3155
3192
3165.4
3267
3382
3325.7
3199
3229
3216
1000 × 10
3138
3166
3150.2
3218
3351
3296.2
3145
3222
3196.3
2000 × 20
3128
3167
3147.3
3210
3345
3294.6
3131
3226
3186.3Known from the results in Table6, for each calculation, MQGA always obtains a better solution than that obtained by MGA and QGA. When the iteration is 500, the range of objective values of MQGA is (3155, 3192) and (3267, 3382) and (3199, 3229) for MGA and QGA. When the iteration is 2000, the ranges of objective values are, respectively, (3128, 3167), (3210, 3345), and (3131, 3226). Therefore, the search ability of MQGA increases faster than that of MGA and QGA as the iteration increases, and MGA and QGA are more likely to be trapped in local optimal solutions.The convergence curves of MQGA, MGA, and QGA obtained for each combination are depicted in Figure7, which indicates the convergence speed of MQGA is faster than that of MGA and QGA, especially for the combination with larger iteration.Figure 7
The convergence curves of MQGA, QGA, and MGA.The solutions distributions of MQGA, QGA, and MGA are given in Figure8, which shows that there is no overlap in the solution space of each algorithm. Moreover, the change in the solutions achieved by MQGA is smaller than that of MGA and QGA, which implies that the quality of solution obtained by MQGA is better than that of MGA and QGA.Figure 8
Distributions of solutions of MQGA, QGA, and MGA.


### 5.2.2. The Comparison Results of MQGA with MGA and QGA

The maximum number of iterations is 1000, and the size of population is 100 (in MQGA, there are two subpopulations, and the size of each subpopulation is 50). To be fair, each algorithm runs 10 times for each instance. The experiment results are shown in Table8, in which each item is an average value. The “Improved%” column represents the percentage difference between the average objective values obtained by the current algorithm with MGA and is calculated by the following formula:(37)
Improved
%
=
(
Ob
¯
-
MGA
¯
)
MGA
¯
×
100
%
,where Ob
¯ represents the objective value achieved by the current algorithm and MGA
¯ represents the objective value achieved by MGA. Boldface and italic indicate the best results for each problem.Table 8
Computational results of MQGA and QGA, MGA.
Problem
Approach
Ob
¯
Improved%
AG
¯
AT
¯
N(FT06)(6 * 6)
MGA
103.3
—
73.3
13
QGA
103.3
0
33.6
30
MQGA
103
0.48%
9.8
113
N(FT10)(10 * 10)
MGA
1282.2
—
300.8
29
QGA
1280.7
0.12%
97.3
98
MQGA
1273.6
0.31%
76.2
163
N(FT20)(20 * 5)
MGA
1330
—
26.7
23
QGA
1328.6
0.09%
18.3
36
MQGA
1323.1
0.32%
36.7
133
N(ABZ3)(10 * 10)
MGA
2088.6
—
430.8
29
QGA
2079.4
0.44%
146.6
101
MQGA
2071.3
0.83%
83.4
164
N(ABZ6)(10 * 10)
MGA
1883.2
—
363.6
29
QGA
1874.4
0.47%
111.3
101
MQGA
1861.4
1.16%
123.3
164
N(ABZ7)(15 * 20)
MGA
1276.2
—
140.3
68
QGA
1233.8
3.32%
104.3
172
MQGA
1228.2
3.76%
86.8
310
N(ABZ8)(15 * 20)
MGA
1290.3
—
114.4
68
QGA
1277.2
1.03%
17.7
172
MQGA
1268.4
1.71%
19.8
310
N(ABZ9)(15 * 20)
MGA
1273.3
—
143.9
68
QGA
1238.3
1.18%
29.7
171
MQGA
1241.3
2.3%
11.3
310
N(TA1)(50 * 5)
MGA
3286.2
—
71.7
84
QGA
3173.3
3.37%
17
426
MQGA
3134.2
4.02%
21.1
778
N(TA2)(50 * 10)
MGA
3346.6
—
91.2
130
QGA
4917
11.33%
9.1
480
MQGA
4803
13.4%
19.3
814
N(TA3)(100 * 5)
MGA
10706.3
—
63.3
171
QGA
9289.33
13.23%
16.2
632
MQGA
9126.4
14.76%
17.8
963
Note:AG
¯ represents the converged generation and AT
¯ represents the average computation time.In the performance of solution, from Table8, MQGA is superior to MGA and QGA for all the test instances, and QGA always outperforms MGA. The improvement becomes more obvious as the scale of the problem increases. For example, for medium-scale problem N(FT20), the optimal ability of MQGA improves 0.32%. For the large-scale problem N(TA3), the MQGA improves 14.76%. The reason is that when the number of jobs increases the transportation of vehicle will become the bottleneck for MGA and QGA. Therefore, MGA and QGA are not applicable to the relative large-scale problem instances. In addition, in the performance of convergence, MQGA can obtain a better solution than MGA and QGA with less iterative generations. For the relative large-scale problems like N(TA1), N(TA2), and N(TA3), MGA and QGA converge faster to the local optimal solution. However, MQGA could jump out of the local optimum with more iterative generations. Finally, in the performance of time, MQGA needs more computation time than MGA and QGA but obtains an improved solution.


## 5.2.1. The Detailed Simulation Result of N(TA1) Problem

In order to validate the performance of MQGA, MGA [13] and QGA [36] are compared. Firstly problem N(TA1) is chosen as the representative instance to analyze the performance of algorithms. Three combinations of iterative generations and population sizes are 500 × 10 (i.e., the maximum number of iterations is 500, population size is 10, and subpopulation size is 5), 1000 × 10, and 2000 × 20. For each combination, the results are presented in Table 6, and more results are given in Table 7, where “BV” and “WV,” respectively, denote the best and worst objective values in the 10 run times. “AV” is the average objective value.Table 6
Detailed comparison results of MQGA, QGA, and MGA on N(TA1) problem.
RUN
500 × 10
1000 × 10
2000 × 20
MQGA
MGA
QGA
MQGA
MGA
QGA
MQGA
MGA
QGA
1
3166
3329
3199
3139
3326
3194
3152
3308
2
3182
3336
3204
3166
3328
3193
3147
3327
3
3165
3382
3229
3154
3351
3209
3154
3345
4
3172
3277
3218
3162
3266
3210
3147
3210
5
3164
3307
3220
3151
3323
3208
3148
3285
6
3179
3330
3215
3165
3319
3170
3155
3324
7
3164
3267
3222
3144
3253
3222
3128
3253
8
3145
3371
3221
3138
3218
3191
3167
3313
9
3154
3297
3226
3163
3293
3221
3140
3305
10
3163
3298
3206
3140
3285
3145
3135
3276Table 7
Comparison results of MQGA, QGA, and MGA on N(TA1) problem.
Size
MQGA
MGA
QGA
BV
WV
AV
BV
WV
AV
BV
WV
AV
500 × 10
3155
3192
3165.4
3267
3382
3325.7
3199
3229
3216
1000 × 10
3138
3166
3150.2
3218
3351
3296.2
3145
3222
3196.3
2000 × 20
3128
3167
3147.3
3210
3345
3294.6
3131
3226
3186.3Known from the results in Table6, for each calculation, MQGA always obtains a better solution than that obtained by MGA and QGA. When the iteration is 500, the range of objective values of MQGA is (3155, 3192) and (3267, 3382) and (3199, 3229) for MGA and QGA. When the iteration is 2000, the ranges of objective values are, respectively, (3128, 3167), (3210, 3345), and (3131, 3226). Therefore, the search ability of MQGA increases faster than that of MGA and QGA as the iteration increases, and MGA and QGA are more likely to be trapped in local optimal solutions.The convergence curves of MQGA, MGA, and QGA obtained for each combination are depicted in Figure7, which indicates the convergence speed of MQGA is faster than that of MGA and QGA, especially for the combination with larger iteration.Figure 7
The convergence curves of MQGA, QGA, and MGA.The solutions distributions of MQGA, QGA, and MGA are given in Figure8, which shows that there is no overlap in the solution space of each algorithm. Moreover, the change in the solutions achieved by MQGA is smaller than that of MGA and QGA, which implies that the quality of solution obtained by MQGA is better than that of MGA and QGA.Figure 8
Distributions of solutions of MQGA, QGA, and MGA.


## 5.2.2. The Comparison Results of MQGA with MGA and QGA

The maximum number of iterations is 1000, and the size of population is 100 (in MQGA, there are two subpopulations, and the size of each subpopulation is 50). To be fair, each algorithm runs 10 times for each instance. The experiment results are shown in Table8, in which each item is an average value. The “Improved%” column represents the percentage difference between the average objective values obtained by the current algorithm with MGA and is calculated by the following formula:(37)
Improved
%
=
(
Ob
¯
-
MGA
¯
)
MGA
¯
×
100
%
,where Ob
¯ represents the objective value achieved by the current algorithm and MGA
¯ represents the objective value achieved by MGA. Boldface and italic indicate the best results for each problem.Table 8
Computational results of MQGA and QGA, MGA.
Problem
Approach
Ob
¯
Improved%
AG
¯
AT
¯
N(FT06)(6 * 6)
MGA
103.3
—
73.3
13
QGA
103.3
0
33.6
30
MQGA
103
0.48%
9.8
113
N(FT10)(10 * 10)
MGA
1282.2
—
300.8
29
QGA
1280.7
0.12%
97.3
98
MQGA
1273.6
0.31%
76.2
163
N(FT20)(20 * 5)
MGA
1330
—
26.7
23
QGA
1328.6
0.09%
18.3
36
MQGA
1323.1
0.32%
36.7
133
N(ABZ3)(10 * 10)
MGA
2088.6
—
430.8
29
QGA
2079.4
0.44%
146.6
101
MQGA
2071.3
0.83%
83.4
164
N(ABZ6)(10 * 10)
MGA
1883.2
—
363.6
29
QGA
1874.4
0.47%
111.3
101
MQGA
1861.4
1.16%
123.3
164
N(ABZ7)(15 * 20)
MGA
1276.2
—
140.3
68
QGA
1233.8
3.32%
104.3
172
MQGA
1228.2
3.76%
86.8
310
N(ABZ8)(15 * 20)
MGA
1290.3
—
114.4
68
QGA
1277.2
1.03%
17.7
172
MQGA
1268.4
1.71%
19.8
310
N(ABZ9)(15 * 20)
MGA
1273.3
—
143.9
68
QGA
1238.3
1.18%
29.7
171
MQGA
1241.3
2.3%
11.3
310
N(TA1)(50 * 5)
MGA
3286.2
—
71.7
84
QGA
3173.3
3.37%
17
426
MQGA
3134.2
4.02%
21.1
778
N(TA2)(50 * 10)
MGA
3346.6
—
91.2
130
QGA
4917
11.33%
9.1
480
MQGA
4803
13.4%
19.3
814
N(TA3)(100 * 5)
MGA
10706.3
—
63.3
171
QGA
9289.33
13.23%
16.2
632
MQGA
9126.4
14.76%
17.8
963
Note:AG
¯ represents the converged generation and AT
¯ represents the average computation time.In the performance of solution, from Table8, MQGA is superior to MGA and QGA for all the test instances, and QGA always outperforms MGA. The improvement becomes more obvious as the scale of the problem increases. For example, for medium-scale problem N(FT20), the optimal ability of MQGA improves 0.32%. For the large-scale problem N(TA3), the MQGA improves 14.76%. The reason is that when the number of jobs increases the transportation of vehicle will become the bottleneck for MGA and QGA. Therefore, MGA and QGA are not applicable to the relative large-scale problem instances. In addition, in the performance of convergence, MQGA can obtain a better solution than MGA and QGA with less iterative generations. For the relative large-scale problems like N(TA1), N(TA2), and N(TA3), MGA and QGA converge faster to the local optimal solution. However, MQGA could jump out of the local optimum with more iterative generations. Finally, in the performance of time, MQGA needs more computation time than MGA and QGA but obtains an improved solution.


## 6. Conclusions and Future Research

This paper studies an integrated scheduling with the materials pickup, flow shop scheduling, and the finished products delivery. The objective is to find a coordinated schedule to minimize the arrival time of the last completed product to the customer. In order to solve the problem, a biologically inspired quantum genetic algorithm is proposed with a new mutualism strategy. The experiment results demonstrate that MQGA can find a satisfactory solution with an acceptable amount of computation time.Future research could address problems with multiple customers or multiple transport vehicles or different shop environments, including flexible scheduling and job-shop. Problems with other performance measures, including minimum mean tardiness, and multimeasures should also be studied.



---

*Source: 387082-2015-03-19.xml*