1. Introduction
Improving production efficiency is one of the final aims of the manufacturing industry. Production scheduling can meet that need, so it is used in manufacturing systems extensively. Considering the differences in factory conditions and job processing requirements, the scheduling problem has been divided into several kinds, such as flow shop [
1,
2], hybrid flow shop [
3,
4], job shop [
5,
6], open shop [
7,
8], and so on. The flow shop problem (FSP) claims that each job should be processed in the same order through every stage, and the hybrid flow shop refers to parallel machines at one or more stages. A job shop environment permits jobs to pass some steps or repeat them at a given machine route, and the processing order of the job in an open shop can be arbitrary. It can be concluded that the flow shop problem is the basic one, then the others can be seen as its variants. In other words, the FSP reflects some common characteristics of production scheduling problems.
In recent years, flow shop scheduling problems with constraints have become research hotspots owing to the actual technology demand. Typically, no-wait attribute, which stipulates all the jobs should be processed from the first stage to the last without any interruption, fits lots of production industries. For instance, in steel processing, maintaining the temperature is a necessary condition, so it calls for no-wait demand. Once you wait, you need more time and energy to heat to the desired temperature, or accept degradation. As early as 1971, Callahan [
9] considered that the steel industry required a dependent processing and the no-wait constraint. Tang and Song [
10] modeled the mill roll annealing process as a no-wait hybrid flow shop scheduling problem. Höhn and Jacobs [
11] also considered continuous casting in the multistage production process a variant of no-wait flowshop scheduling. Yuan and Jing [
12] explained and affirmed the importance of the no-wait flowshop problem to steel manufacturing. Similarly, Reklaits [
13] pointed out that the no-wait constraint is also suitable for the chemical industry due to its instability. Hall [
14] emphasized that a similar situation also occurs in plastic molding and silver-ware industries, as well as the food processing industry [
15]. In addition, there are also many other practical problems that can be solved by the assistance of being modeled as no-wait scheduling problems, such as controlling health care management costs [
16], developing a metro train real-time control system [
17], and so on.
In this paper, the NWFSP with makespan criterion is considered because of its universality. It minimizes the maximum completion time of the jobs (
), maximizing the efficiency of machine usage [
18]. The problem can be described as
when the symbolic notation proposed by Graham [
19] is applied. The first field represents the production environment, which is a flow shop with
stages. The second field shows the constraints and “
“ is for “no-wait constraint”. The performance criteria occupy the third field.
The research on NWFSP began around the 1970s. Garey and Johnson [
20] proved that the NWFSP is a strongly NP-hard problem when the number of machines is no less than three. At the same time, the problem was documented as NP-hard for more than two machines by Papadimitriou and Kanellakis [
21]. Due to the complexity of NWFSP and the wide range of its applications, many heuristic algorithms have emerged or been tried by Bonney and Gundry [
22], King and Spachis [
23], Gangadharan and Rajendran [
24], Laha [
25], and so on, which aimed to make the solution as close as possible to the optimal within an acceptable time. In addition, exact methods [
26] like branch-and-bound [
27] and the column generation method [
28] are widely used to obtain optimal solutions for small-scale problems. Metaheuristic algorithms, whose search mechanisms are inspired by phenomena in various fields of the real world are also active for solving NWFSP. Hall [
14] and Allahverdi [
29] jointly detailed the research progress of NWFSP before 2016. A typical example is the discrete particle swarm optimization algorithm [
30]. In recent years, various improved or emerging metaheuristic algorithms have been used to solve the NWFSP. For example, a hybrid ant colony optimization algorithm (HACO [
31]), a discrete water wave optimization algorithm (DWWO [
32]), a factorial based particle swarm optimization with a population adaptation mechanism (FPAPSO [
33]), a single water wave optimization (SWWO [
34]), a hybrid biogeography-based optimization with variable neighborhood search (HBV [
35]), a quantum-inspired cuckoo co-search (QCCS [
36]) algorithm, and an improved discrete migrating birds optimization algorithm (IDMBO [
37]) have also been used with good results. In 2021, Lai used the Discrete Wolf Pack Algorithm [
38] to solve the NWFSP problem. Each algorithm has its own advantages, which makes it perform better than others in certain applications. In order to give full play to the advantages of each algorithm and enhance the generality of the algorithm, hyper-heuristic (HH) [
39] and matheuristic [
40] algorithms have appeared. The emergence of a new generation of large-scale mathematical programming optimizers such as CPLEX [
41] and Gurobi [
42] has greatly improved the performance of matheuristic algorithms that combine heuristic algorithms and exact algorithms. In recent years, matheuristics [
43] have been widely used in scheduling problems, including NWFSP. Since NWFSP can be transformed into the ATSP [
44], it is effective to incorporate more mathematical computations into its search mechanism.
Existing NWFSP researches mostly focus on designing algorithms from the perspective of instances rather than problems. The heuristic rules or the component mechanisms and parameter selections of the metaheuristics are determined by the quality of the test results, so as to balance the local search and the global search. In this way, the algorithm is dependent on the instances, and the performance of the algorithm will be different for different instances. This paper aims to explore the heuristic rules from the perspective of the problem, so that the rules show consistent effects on the instances, thereby enhancing the stability of the algorithm. Existing neighborhood structures for local search focus on insertion, swap, destruction-construction, etc., and do not divide the already searched area and the unsearched area, so the same sequence cannot avoid being tested repeatedly. Likewise, good structures are not necessarily preserved after multiple iterations. In this paper, by expanding the search domain gradually, the search process is kept in the unsearched potential areas, and the repeated calculation is reduced. The experimental results in
Section 4.1 show that such a search method can also effectively preserve good structures. Exact algorithms are limited by the complexity and scale of the problem, so they are less applicable than heuristics and metaheuristics. Using heuristic rules to guide the search direction of the exact algorithm can improve the solving efficiency and increase the application scope of the exact algorithm. The experiments in
Section 4.2 demonstrate that by fully exploiting the advantages of both methods, better results can be obtained than the current mainstream metaheuristic algorithms for solving NWFSP problems.
The rest of this article is organized as follows.
Section 2 describes the problem of NWFSP.
Section 3 proposes two heuristic algorithms for improving the solutions of NWFSP.
Section 4 discusses the performance of the proposed matheuristic algorithm on two benchmarks.
Section 5 summarizes the contributions of this article and provides some suggestions for future research directions.
2. No-Wait Flow Shop Scheduling Problem
Compared to FSP, NWFSP permits no interruption when a job is processed from the first stage to the last one. The details are expressed as follows.
Typically, the NWFSP includes stages, preparing for jobs to follow the same given route, where the and are determined. It also restricts only one machine at each stage, so it contains machines overall. Furthermore, the processing time of job on machine is also known and may be named as . It is natural to consider every as the element of row , column in processing time matrix . One job should be processed on only one machine at the same time and it is merely allowed to be transferred to the next stage until the current work is finished; similarly, each machine cannot deal with more than one job at any time. Here, one only considers the situation of ignoring the setup time and transfer time.
2.1. Problem Description
As we known, makespan, the complete time of the last job in a sequence, is one of the indicators that directly measures economic benefits. Let denote the minimum makespan and its expression is shown in Equation (1). The meanings of related symbols are as follows.
The job set:
The machine set:
any sequence
the set of all possible sequence
the job in the sequence
end time of job on machine i
is the delay of the job when it is arranged after job ,
the tail value of the job .
the sum of the processing times of all the jobs on the first machine
Owing to the no-wait constraint, the delay value between two adjacent jobs is constant and can be calculated by matrix
. As shown in
Figure 1,
has been divided into three parts and can be calculated from Equation (4). Obviously,
is a constant, so the change of
is only affected by two parts. Then, the simplified objective function was formulated in Equation (5).
2.2. Equivalent Model
Considering jobs as cities and the delay value as the distance between two cities, the NWFSP can be transformed into the ATSP. After adding virtual city/job “0”and “
n + 1”, the model can be described as follows. For convenience, “0” and “
n + 1” are counted as the same city, and “0” represents the starting city/job, and “
n + 1” represents the ending city/job. The city/job set is
.
3. The Proposed Heuristic Algorithms
The NWFSP can be transformed into the ATSP, so the path-exchange-strategy used in the Lin-Kernighan–Helsgaun (LKH) [
45] algorithm is also suitable for solving NWFSP. Based on this idea, we can divide the neighborhood and construct heuristic algorithms.
3.1. Basic Neighborhood Structures
Considering the search efficiency fully, the search solution domain is initially limited to three levels, that is, Neighborhood
,
,
. The characteristics of the three basic neighborhoods are given in
Table 1. The
optimum means that the solution cannot be better by exchanging
areas, which are arbitrarily divided in the sequence.
The
neighborhood is shown in
Figure 2 as an example, where
and
are the two area lines and the three regions are
,
,
.
3.2. Selected Potential Sub-Neighborhoods
In order to improve the search efficiency and advance the search process in more potential sub-domains, three rules are proposed and three sub-domains are selected, namely
,
, and
. The details are shown in
Table 2. Two heuristic algorithms will be carried out along these three neighborhoods.
Considering that the current solution sequence is almost derived from the last region-exchange operation, the combinations containing the adjacent area number pair (that is “01”, “12”, “23”, “34”) as a part do not temporarily need to test again.
From the experimental results, along the evolutionary direction, the job in last position is easier to have coincidence with the optimal solution because of being affected by values. Then is the first position. When the evolution process has advanced to a certain level, the first and last job can be fixed.
The lower the neighborhood level, the higher the search priority. The algorithm flow chart is shown in
Figure 3. “US” is the abbreviation of “Update successfully” and “UF” is the abbreviation of “Update failed”. For ease of representation, the structure in the dashed box in
Figure 3 is referred to as “PHM”.
3.3. Two Heuristic Algorithms
Based on the mechanism in
Figure 3, it is easy to propose two heuristic algorithms according to whether the
value characteristic is highlighted. The two heuristic algorithms HG1 and HG2 differ only in the first step, that is, the input feasible solutions are different. The details are shown in
Figure 4. The advantage of HG1 is randomness. HG2 always takes the job with the smallest
value as the final job to start the search, which is more in line with the characteristics of NWFSP and can speed up the algorithm convergence.
In order to improve the utilization of historical data and improve the efficiency of neighborhood search. According to the different characteristics of the neighborhoods of , and , corresponding search algorithms are provided one by one. They are vmcmp2, vmcmp3, and vmcmp4.
3.3.1. Local Evolutionary Algorithm in Z2
Searching the only manner {102} can improve evolution speed. Without the change of
, the evaluation speed is further improved. Let the size of
characterize the degree of improvement of the new sequence; the larger the value, the better the improvement.
In
Figure 2, it can be seen that for the same
, when
and
move,
,
,
,
are changed, but
,
remains unchanged, which can be regarded as constants. Because the first job is stored in
, once the adjustment is made according to {102}, the
value of the new sequence must change, and the value of
will also change. As can be seen from Equation (15), in order to make full use of historical data and reduce repeated calculating,
is regarded as two components. The first part and the second part are respectively represented by functions
and
, as shown in Equations (16) and (17).
The two-dimensional array is constructed, and the calculation result of Equation (17) is stored in the row and the column of . Such an operation can ensure that after the sequence is updated, the intermediate result only needs to recalculate part of the data. For the convenience of description, several definitions are given below.
Definition 1 (Point). Adjacent job pair makes up a point. If is a point,
is the predecessor job of
, and
is the successor job of
.
Definition 2 (Single-point Attribute (SPA)). The function value calculated by a point only or with information of a stable job. For example, function is the single-point attribute of the point
, where
is the stable job.
Definition 3 (Double-point Attribute (DPA)). The function value is calculated by two points and reflects the relationship of four related jobs. For example, function is the double-point attribute of point
and point
.
When traversing , function is associated with the attributes of job which belongs to the second point . But in essence, is still a job in . If the index number of the first job in the sequence is 0, then the index range of is (denote the index of in the sequence as , then ). Function is seen as the single-point attributes of the jobs whose index ranges from 1 to . Once is given, we can directly calculate the function . It can be seen from the characteristics of the double-point attribute that if the two points do not change, then their DPAs do not change. That is, as long as the con1 is met, the saved values in need not be updated, and the calculation speed can be increased.
The con1 is described as follows.
First, the jobs and remain in the original pair order, that is, in the sequence , and is still on the left side of .
Second, to find the value of the job and job , that is, in , x cannot be equal to or in the last operation; cannot be equal to or .
The local optimal algorithm (Algorithm 1) for
is described as follows.
Algorithm 1. vmcmp2(con1)// Search {102} only |
Step 1. compute all the SPAs in CS, that is, all the |
Step 2. compute all the DPAs in CS, that is, all the |
Step 3. //Save the best positions to which provide the |
Step 4. if |
L1 = and L2 = |
CS=CS ({102})//Perform operation {102} on CS to update CS |
update (); //Recalculate all the f values |
update (, ); // Recalculate g values that do not meet con1 |
jump to Step 3 |
end if |
else// CS is the best of Z2 now |
Stop searching in Z2 |
end else |
3.3.2. Local Evolutionary Algorithm in Z3
Since the first job and the last job are constant, the point attribute can be used with the historical data to optimize the algorithm.
Figure 5 shows the evolution form of “0213”, from which Equation (20) can be obtained.
Where
represents the sum of the delay values of all the jobs from the job
to the job
in the current sequence. Equation (20) is the sum of the expressions in the three brackets, and the three sub-expressions, which can be seen as
DPAs are in a rotation-symmetric relationship. Parentheses 1 and 2 can share the function
. Parenthesis 3 is denoted as function
.
and
are considered as functions of the job numbers and the function values are saved. See Equations (21) and (22) for details.
Compute and save each
value to the
row,
column of the two-dimensional array
. Assume point1
, point2
, point3
, then Equation (20) can be replaced by Equation (23) to reduce double counting.
In contrast to the search in , there is only one combination to try, so compute all and choose the biggest improvement. If all is computed but only updated at most once, the utilization rate of historical data is too low. Then Multiple Iterations in One step (MIO) is proposed to improve this situation.
Traverse the
Z3 neighborhood of
CS, the specific method is to calculate by Formula 23, and save all the values greater than 0 in
. In
Figure 6,
is a multiset that sorts the results according to its
values, from largest to smallest, and itRS7 is the iterator of
.
.
Suppose there is a NWFSP with
, and the
is shown in
Figure 6. As shown in
Figure 7, its current solution is
.
Take out the element pointed to by the iterator
, that is
.Then the evolution process can be seen in
Figure 8. The indices of points (1, 2), (3, 4), (10, 11) in the sequence are respectively recorded as
. In
Figure 6,
means
.
Because of the constant last job, the
is constant too. Then the essence of every evolution manner is the change of boundary values. Combining
Figure 6 and
Figure 7, we can draw the conclusions of
Table 3.
It is easy to see from
Table 3 that no matter how many iterations are completed, as long as the changed path of boundary values still exist, evolution in this direction will succeed.
It can be seen from
Figure 9 that after two iterations, the corresponding point of
cannot be found in
at this time, so it cannot be evolved in this way. Next, derive the conditions
con2 whose paths still exist after iterations, as shown in
Table 4.
Traverse
and determine whether each element meets the
con2 in turn. If one is satisfied, it is applicable and the iteration can be completed directly. Therefore,
is calculated only once, and there is a chance to iterate multiple times. The algorithm (Algorithm 2) vmcmp3(
con2) suitable for
is described as follows.
Algorithm 2. vmcmp3(con2)
|
Step 1. Input , record Loc (j, ) and save them to vector k1, Let its jth element, k1(j) to save the location of job j in |
Step 2. flag = 0// Set the initial value of the flag |
Step 3. Compute and save , , then compute and save them to multiset DY0 when . |
Step 4. if (DY0.size()! = 0) // if DY0 is not empty |
for itRS7 = DY0.begin() to DY0.end ()-1)// Traverse DY0 by itRS7 |
if(con2) //If con2 is met |
CS= update (itRS7, {2013}, CS);// According the locations |
//provided by DY0, give CS a {2013} operation |
flag=1// Updated successfully |
update (Loc (j, CS)); // Save the new correspondence of location //and job number |
end if(con2) |
end for |
end if |
else |
flag=2 |
end else |
if flag==1 // When update successfully |
jump to Step2 |
end if |
else // When update failed |
Stop searching in Z3 |
end else |
3.3.3. Local Evolutionary Algorithm in Z4
Experiments show that the combination {03214} has a higher hit rate. Its main feature is that the first and last jobs are not changed. The
Figure 10 shows the evolution form.
Still, start with the expression of
, specific as shown in Equations (24)–(26).
Similarly, using Equation (27) to uniformly calculate the two parts separated by parentheses in Equation (26) can simplify the calculation. Let
be
here.
Compute all the results by Equation (27) and save them in multiset .
Use two iterators,
and
to traverse
. Then select all the elements that meet the condition of
con3, and store them in multiset
. See
Figure 11 and
Figure 12 for details.
collects all successful evolution paths. Its iterator is
itRS9 in
Figure 12 (
). Suppose the variables saved by a single element in
are
and
, then the conditions for
con3 to be established are shown in
Table 5.
Similarly, use MIO Strategy to construct an algorithm (Algorithm 3) to speed up the search process.
Algorithm 3. vmcmp4(con3, con4) |
Step 1. Input , record Loc (j, ) and save them to vector k1, Let its jth element, k1(j) to save the location of job j in |
Step 2. flag = 0// Set the initial value of the flag |
Step 3. for = 0 to n-4 // The union of the ranges L1 and L2 |
for = 3 to n − 1// The union of the ranges L3 and L4 |
cmpRS2(,)//compute the DPAs of (, , , ) and save them in RS2 |
end for // End the inner loop |
end for// End outer loop |
Step 4. it1begin = RS2.begin() |
it2begin = it1 + 1 // Assign initial values to iterators it1 and it2 |
for it1 = it1begin to it1end |
for it2 = it2begin to it2end // Double loop for traverse RS2 |
if (con3(it1, it2)) //If the elements pointed to by it1 and it2 satisfy con3 |
JSF. insert (it1, it2)//Save the evolution path provided |
//by the elements which are pointed to by it1 and it2 in JSF |
end if |
end for// End the inner loop |
end for// End outer loop |
Step 5. if (JSF. size ()0)//If JSF is not empty, it means that can be |
//updated in its neighborhood Z4 |
for itJSF9 = JSF. begin () to JSF. end ()−1 //Traverse multiset JSF |
//which is the neighborhood Z4 of |
if (con4(itJSF9))//If the element pointed to by itJSF9 satisfies the con4 |
CS = update (JSF9, {03214}, CS)//Update CS as specified by JSF9 |
update (Loc (j, CS)); // Save the new correspondence of |
//location and job number |
flag = 1; // Assign 1 to the flag bit |
end if |
end for |
end if |
else |
flag = 2 |
end else |
if flag == 1// If there is an update in this round |
jump to Step2 |
end if |
else// This round of searching this in Z4 has not been updated |
Stop searching in Z4 |
end else |
Let the condition that satisfies MIO in Z4 be con4. If con4 in the algorithm is to be established, two conditions must be met, as shown in
Table 6.
3.4. Experimental Results of the Two Heuristic Algorithms
Use HG1 and HG2 to solve Tailard (TA) benchmark. Then record the calculation results in
Table 7. The heuristic algorithms were coded in C++. All the experiments executed in Windows 10 on a desktop PC with an 11th Gen Intel(R) Core (TM) i7-1165G7 (2.80GHz) processor and 16.0 GB RAM. The calculation formula of
RPI (Relative Percentage Increase) and
are shown in Equations (28) and (29), where
refers to the optimal
of instance a1.
refers to the result of algorithm
i1 to solve the instance
a1.
In
Table 7, HG2 has a smaller
value than HG1 in all scale instances; its fluctuation range is from 0.8% to 2.58%, the maximum
value of 2.58% appears on the small-scale instances (scale is 20 × 10), while its
value on the largest-scale instances (500 × 20) is only 1.28%. On the whole, it performs better for large-scale instances. The average
is 1.30%, and the fluctuation is not large, indicating that the algorithm has strong local search ability, and can quickly converge to a local optimal solution regardless of the scale, with stable performance.
It can be seen from
Table 7 that out of 120 instances, HG1 only exceeds HG2 in 6 instances (the part in bold in the table). Therefore, for the goal of approximating the optimal value, the effect of the HG2 algorithm constructed by introducing the
value is much better than that of HG1, and even HG2 directly obtains the optimum on some small-scale instances (the part in italics in the table). Taking “TA 05” as an example, its Gantt chart is shown in
Figure 13. If the result of HG1 is imported into Gurobi, another set of optimal solution sequences can be obtained. This Gantt chart is shown in
Figure 14. Therefore, the use of different heuristic algorithms is effective for the breadth search. The same job-blocks in the HG1 sequence and the HG2 sequence can also be considered as high-quality structures by comparing them with the optimal solution sequence. Therefore, the heuristic algorithms can easily construct multiple high-quality solutions. Arrange the jobs according to their
values from small to large, take them as the last job in turn, and then perform the PHM operation.
Similarly, we used HG1 and HG2 to solve Reeves (REC) benchmark. The results are recorded in
Table 8.
On this benchmark, the values obtained by HG2 are also better than that obtained by HG1. Except for REC11, HG1 does not get a better than HG2, and HG2 even gets the optimal solution on REC03. Obviously, the effect of HG2 is better than that of HG1, so HG2 was selected to participate in the calculation in subsequent studies.