Improved Beam Search for Optimizing No-Wait Flowshops With Release Times

Production management of perishable goods is highly complex and requires well-informed decisions in corresponding stages. In such production environments, scheduling problems with time constraints are of high relevance to ensure the timely flow of the work-in-process material and goods. This study introduces the no-wait flowshop scheduling problem with release times (NWFSP-RT) to help advance decision support systems in the food production industry. For this purpose, an original mixed-integer linear programming (MILP) formulation is proposed for minimizing makespan. A BS algorithm, and an improved variant, the local search-based Beam Search (BSLS) algorithms are developed to solve the NWFSP-RT problem. Extensive numerical analysis is conducted to analyze the performance of the algorithms in solving this highly intractable extension of the scheduling problems. We showed that BSLS effectively avoids early convergence and local optimality while dismissing non-promising search directions within a partial enumeration solution approach. The statistical analysis confirmed that the improved BS algorithm performs better in terms of solution quality. Applications of the developed heuristic are worthwhile research topics to pursue in solving other complex optimization problems.


I. INTRODUCTION
No-wait scheduling problems (NWSP) have received recent recognition due to the wide industrial applications in steel manufacturing [1], metal, plastic, chemical, and food industries [2] as well as the modern manufacturing environments [3]. Given the importance of scheduling in the design of decision support systems [4], case-specific formulations, and effective solution algorithms are of high significance to help modernize the production systems. Recent studies developed state-of-the-art solution approaches [5]- [9] to facilitate the industry-scale application of the NWSPs, which are rather complex. Many other studies extended the NWSP, for example through including setup times [10], machine maintenance [11], no-idle-time situation [12], and distributed The associate editor coordinating the review of this manuscript and approving it for publication was Chi-Tsun Cheng . production environment [13], among the other examples (see [14]), to address case-specific industrial needs.
In reality, there are production situations where a particular process can begin after a certain time, the release time [15]. The extant scheduling predominantly considered a release time equal to zero to simplify the problem. Reference [16] was the first to include the release time constraints in a single-machine production environment. Since then, many studies focused on developing solution algorithms to solve the single-machine scheduling with release times; Genetic Algorithm (GA; [17]), learning-based Branch-andbound of [18], [19], Simulated Annealing of [20], Ant Colony Optimization of [21], and the hybrid GA developed by [22] are the seminal examples. Some studies extended the basic single-machine scheduling with release times, for example through developing a two-agent single-machine problem with release time [23], considering expiration time in addition to the release times [24], and including sequence-dependent setup times [25].
Two-machine and flowshop scheduling problems with time restrictions, however, are relatively less developed. Reference [26] is one of the first examples that developed solution algorithms to solve the two-machine scheduling with time constraints. Addressing the same problem, [27] proposed a branch and bound algorithm that yields better solutions compared to the earlier study. Reference [28] further improved the solution quality to the two-machine scheduling problem with time constraints in a reasonable computational time. Most recently, [29] analyzed the complexities involved in the scheduling problems with release time constraints in a service application area. From the limited studies of flowshop scheduling problems with time constraints, reference [30] was the first to address permutation flowshop scheduling with release times. Proposing eight solution algorithms, they showed that the local search-based GA outperforms the other algorithms, including the constructive heuristic they put forward. Reference [31] solved the flowshop scheduling problem with release dates and sequence-dependent setup times in a multi-objective solution scheme. Reference [32] proposed heuristics for solving flowshops with delivery times, which is a variety of job release dates. These studies did not consider the no-wait production situation that is necessary for some industrial applications, as well as the advanced manufacturing environments featured by just-in-time production.
Taking the canned-food and refreshment production as an example, the cooking operation has to be followed by canning in a no-wait flowshop situation such that the freshness of the products can be guaranteed. Besides, due to the operational characteristics of the consumer goods industry [33], [34], more particularly its time-sensitive nature, the release time constraints are necessary to optimize the supply chains more effectively. The no-wait flowshop scheduling problem with release time (NWFSP-RT) is, therefore, of high relevance in the optimization of production operations in canned-food and other similar industries. Although the NWFSP-RT problem has practical applications in many industries, it has not been studied thus far. Recognizing the mentioned research gap, this study sought to propose effective and efficient algorithms to solve this scheduling extension. The main contributions and novelty of this paper are as follows.
• This study is the first to address the NWFSP-RT problem in the literature, which helps expand the boundaries of research in the field of scheduling theory.
• To facilitate the industry applications of the NWFSP-RT problem, an original mixed-integer linear programming (MILP) formulation to this problem is developed.
• Given the inherent complexities involved in the NWFSP-RT problem, a Beam Search (BS) and a local searchbased improvement of BS (BS LS ) are developed to solve the proposed problem effectively and efficiently.
• The basic BS algorithm has the problem of premature convergence and getting trapped in local optimality. To handle this issue, the proposed BS LS algorithm integrates a local search-based mechanism to the partial-enumeration solution approach to help avoid early convergence and local optimality. Overall, this work contributes to the advancement of decision support systems in the food production industry meeting the practical requirements for solving this computationally intractable scheduling extension.
The rest of this manuscript begins with a review of NWFSP's literature in Section 2. The proposed mathematical formulation is presented in Section 3. Section 4 elaborates on the developed solution algorithms. Extensive numerical analysis is provided in Section 5 to examine the performance of the proposed algorithms. Section 6 concludes this research by reflecting on the major takeaways and directions for future research.

II. LITERATURE REVIEW
Reference [35] introduced Flowshop Scheduling Problems (FSP). Inspired by this seminal work of the scheduling theory, many of the FSP-related articles focused on extending the mathematical formulation for various industry applications. As a prime example, the no-wait constraint is of particular relevance for the applications in steel manufacturing, chemical, food, and modern industries [1]- [3]. This section reviews the no-wait flowshop scheduling extensions published in the Web of Science database, considering machine features, i.e. two-and m-machine variants, job-related constraints, as well as the studied performance indicator as suggested by [36].
The NWFSPs was first time introduced by [37]. Working on this seminal flowshop extension, the early works were mostly focused on developing general heuristics [38]- [46], branch-and-bound methods [47], [48], and other heuristic approaches [49] to solve the two-machine NWFSPs. Reference [50] showed that the no-wait flowshops with a makespan objective function is equivalent to the traveling salesman problem and that the problem with fixed processing times is solvable in polynomial time. References [51], [52] developed polynomial-time algorithms to solve the twomachine NWFSPs considering availability constraints. References [53] and [54] included uncertain setup times and due date constraints to two-machine NWFSPs, respectively. Reference [55] introduced convex resource-dependent processing times, and, more recently, [56], [57] integrated the convex resource allocation and common flow allowance, respectively, into the two-machine NWFSPs considering the learning effect.
NWFSP with more than two machines is more realistic, hence, has received wider recognition in the academic literature. Job-related constraints for case-specific needs are developed to further facilitate the real-world applications of m-machine NWFSPs, among which, the setup times are studied the most [10], [58]- [61]. Sequence-dependent [62]- [65] and group/family setup times [66]- [68] are the other examples of this type. Besides, [69] simultaneously considered precedence constraints and setup times. The mixed-no-wait situation was examined by [70].
References [71] and [72] considered a makespan constraint while minimizing total completion time. Reference [73] minimized makespan in NWFSPs considering a threshold for a mean completion time constraint. NWFSP for two-batch processing machines was studied by [74]. NWFSP was also adapted to address lot-streaming problems [75]. Reference [76] considered NWFSPs in a distributed production environment. No-wait situation under hybrid flowshops with parallel machines was investigated by [77]- [79].
The other studies solved m-machine NWFSPs considering various performance indicators. The makespan criterion, i.e. the maximum completion time of all jobs, which is concerned with better utilization of the resources, is considered the most for solving NWFSPs [80]- [90]. Total completion time, i.e. the total time required to finish the jobs was considered by [2], [91], [92] to solve NWFSPs. References [93]- [95] solved NWFSPs considering total flow time, which is a measure to address the work-in-process inventory. Other studies addressed service-oriented performance indicators, like total tardiness [96], [97], the number of tardy jobs [98], and maximum lateness [99], [100]. Sum of the makespan measure, and total weighted tardiness [101], as well as makespan and total completion time [101], are among the multi-objective NWFSPs solved considering a priori preference articulation. Alternatively, some other studies developed multi-objective optimization algorithms under a posterior preference articulation to solve the NWFSP addressing conflicting objectives, like makespan and total tardiness [102], [103]; makespan and maximum lateness [104].
Overall, m-machine NWFSP considering release times and the makespan criterion has not been studied thus far. Besides, no studies explored the BS algorithm for solving the NWSPs. The present study sought to address the mentioned gaps, developing a MILP formulation and an improved BS algorithm to help facilitate the case-specific applications of NWFSP-RT.

III. MATHEMATICAL FORMULATION
Following the well-known three-field notation suggested by [105], the NWFSP-RT problem can be denoted as F m nwt, r j C max , with F m determining the flowshop production environment; C max specifying makespan as the objective value, and nwt, r j pointing out the no-wait and release time settings. The F m nwt, r j C max problem is formulated assuming m available machines for processing n jobs. A deterministic release time, r j , is considered to impose a constraint on the start of job's processing. That is, the job j cannot be processed before the defined ready time. A processing time, P j,i , is also associated with jobs j ∈ {1, 2, . . . , n} when they are processed on machine i ∈ {1, 2, . . . , m}. Assuming zero-buffer capacity between the machines, a job can be released by the machine i only when the machine i + 1 is available to seize it. In this situation, the process i + 1 begins immediately after completion of the process i. The remainder of assumptions is as the following.
• Jobs cannot be interrupted or re-assigned once a machine starts processing them.
• Each machine can only process one job at a time, for the jobs being processed by one machine at a time.
• Jobs are assumed to follow the same processing order.
• Assuming that machines do not breakdown, no maintenance operations are considered in between the production operations. The following indices, parameters, and decision variables are defined to provide the MILP formulation. The problem is to find a permutation of n jobs, to be processed in m machines such that the makespan value, shown in Equation 1, is minimized.
The objective function is subject to the following constraints: n k=1 x j,k = 1; j ∈ {1, . . . , n} n j=1 Constraint (2) ensures that the job in the position k cannot be processed on the first machine (i = 1) until the release time is reached. Constraint (3) is defined to guarantee that the completion time associated with the job with a specific position in the job sequence is greater than the summation of the completion time of the jobs in the previous positions and the processing time of the job situated in the current position. Equation (4) helps calculate the job completion time. Given that there is no buffer zone for the work-inprocess items, the completion time of a job on the machine i is equal to its completion time on the machine i − 1 plus the respective processing time. Constraint (5) is a guarantee that the completion time of a job on a specific machine is greater than the respective processing time. The maximum completion time, makespan, is calculated using constraint (6). According to Equations (7)(8), only one job can be situated at one specific position in the job sequence. Finally, the integrality constraints (9-10) specifies the type of decision variables.

IV. PROPOSED SOLUTION ALGORITHM
Inspired by the cluster search approach of [106] applied in the artificial intelligence context, reference [107] introduced the BS method in its current form by including the searchfiltering mechanism. Since then, BS has been successfully applied to solve several scheduling extensions. Adding a recovery mechanism to the original BS, [108] solved the single-machine scheduling problem with release times, outperforming the then-best-performing solution algorithms. Applying the same solution approach, [109] minimized the maximum completion time of the unrelated parallel machines scheduling problem. Reference [110] showed that both the basic and recovery BS algorithms effectively minimize the total earliness-tardiness in single-machine scheduling with no idle-time in large-scale problems. BS was integrated with the Ant Colony Optimization [111] and the Genetic Algorithm developed by [22] to solve the open shop scheduling and single-machine scheduling problems with release dates, respectively. Reference [25] successfully tested the BS algorithms to solve the single-machine scheduling problem with release dates and sequence-dependent setup times. Most recently, [112] tackled the permutation flowshops using the BS algorithm, minimizing total flow time. Overall, the BS algorithms have demonstrated a strong track record to solve complex scheduling problems, particularly those with time constraints.
The BS method is similar to the renowned branch and bound algorithm, except that it avoids full enumeration of all solution nodes to shorten the computational time. Besides, BS does not include the back-tracking mechanism, where the algorithm cannot search among the solutions in the previous levels once the nodes are branched. Filtering rules, including the maximum number of allowed nodes per branch, α, and per layer, β, are used to keep the best-performing partial solutions and avoid searching less promising solution areas.
BS structure consists of nodes, π L, i , branched into L layers. For a given node, π L, i , L demonstrates the number of fixed jobs, the associated layer, and i identifies the various job assignment possibilities. π 0 demonstrates the nodes with no fixed job positions at level 0, and π L−1, i identifies the root node for π L, i . At any layer, the set θ L comprises the jobs that are yet to be assigned to the sequence. Of the n − L + 1 possible partial solutions for branching a given route node, only the α best alternatives with the smallest idle-time, to be calculated using Equation 11, will be considered for further computations. The partial solutions under node i will be saved in set B L,i in level L and B L = B L,1 ∪ . . . ∪ B L,i represent the collection of all B L,i in the layer. Of the B L alternatives, β nodes with the best objective function value of the partial solutions, to be calculated using Equation 12 will be considered as activate nodes for the rest of the computational steps.
I π L, i and F(π L, i ) demonstrate the idle-time and makespan of the partial solution π L, p , respectively, based on which nodes at the branch and layer filtration steps get sorted.
The basic BS algorithm may result in a local optimal solution. Reducing the likelihood of being trapped in local optima, local search approaches help improve the BS algorithm's performance for tackling optimization problems with complex solution space. To avoid early convergence and being trapped in local optimality, the developed BS LS considers two local search procedures. In the first approach, two positions from the job sequence are randomly selected for the exchange, swap procedure, whether it has resulted in any improvements or not. In the insertion procedure, a job is randomly extracted from the sequence and inserted into all the possible positions until a better fitness value is obtained. In this procedure, the number of executions has a direct impact on the quality of the solution and the computational time. Considering n × m × k function evaluations in this study, different values of k is considered to conduct the experiments on the problems with different workloads. The pseudocode of the computational procedure is provided in Figure 1.
The seven-step procedure is described below.
Step 1. Set the beam parameters α, and β(and the local search threshold, if necessary) considering the problem size and the computational needs. Set L = 0, and π 0 equal to an empty vector. Save all the available jobs in the set θ 0 .
Step 2. Branch the root node π L,i into L + n − 1 nodes, assigning one of the available jobs from the set θ L to fix the L th position at the sequence. The branching possibilities may be less in the presence of pre-defined ordering conditions. Set L ← L + 1. Step 3. Calculate the idle-time, I π L, i , of every θ L job's positioning in the partial sequence; sort them on an ascending order, and select α from the top of each list; save the resulted partial solution with a newly-fixed job position in the set B L+1,i . Continue this procedure until all the root nodes in the layer L are branched.
Step 4. Apply random local search and replace B L,i with B L,i , if the objective function value associated with B L,i is worse than the resulted neighborhood partial solution. This step will be skipped for the basic BS algorithm.
Step 5. Given B L = B L,1 ∪ . . . ∪ B L,i as the collective set of new nodes, if the number of B L members is less-than-or-equal-to β, branch all of them; otherwise, calculate the respective F(π L, i ), and select the best β alternatives as active nodes for branching. Discard the remainder of the partial solutions.
Step 6. Return to Step 2, and continue until θ L is empty (L < n − 1). At this stage, each node represents a complete solution.
Step 7. Calculate the objective function value of all θ n−1,i nodes; the solution associated with the best node is the desired (near-) optimum solution.
An illustrative example, including n = 6 jobs and m = 3 machines, is provided in Table 1 to clarify the computational procedure. In this example, α = 2 and β = 8 are considered as the BS algorithm parameters to go through the solution process. Figure 2 elaborates on the computational steps on a level-by-level basis. Given n = 6 jobs and the fact that there are no pre-ordering conditions, π 0 is branched into six nodes. Of the possible choices in level L = 1, two with the shortest idle time, I π 1,6 = 24 and I π 1,4 = 27, π 1,6 and π 1,4 are selected as active nodes. Considering π 1,6 and π 1,4 as root nodes for the L = 2 nodes, a total of 10 (5 for each root node) nodes have resulted. Keeping two of the options with the lowest idle time from each of the two branches, the beam set in level 2 is B 2 = π 2,2 , π 2,4 , π 2,6 , π 2,10 . Given that the number of active nodes is still less than the predefined threshold, β = 8, all the B 2 nodes will be branched. Continuing the same procedure to level 4, a total of 24 nodes, B 4 = π 4,1 , . . . , π 2,24 , are resulted where the non-active nodes are first truncated based on the I π . Then, β = 8 with the smallest makespan value, π 4,1 , π 4,2 , π 4,5 , π 4,7 , π 4,9 , π 4,13 , π 4,19 and π 4,23 are selected as root nodes for the next stage, where the fifth job position will be fixed. At the final level of the BS algorithm, L = 6 in this example, the comparison amongst the complete solutions are presented in Table 2.

V. NUMERICAL ANALYSIS
This study solves the proposed MILP formulation via the LINGO 12.0 software to obtain the optimal solution for small-scale test instances. The proposed BS and BS LS algorithms were coded on Visual C++ using a personal computer with the following specs: Intel(R) Core (TM) i5-3337U (1.80GHz) processor, 8GB of RAM, and a Windows 10 operating system. This section continues by elaborating on the test instances and the performance measure. The heuristic algorithm parameters calibration is explained next. Finally, the numerical experiments are provided, where the algorithms' performance over the small, medium, and large test instances are compared and validated using statistical analysis.

A. TEST INSTANCES
Four different groups of test instances are considered to conduct numerical tests. To account for the release times, Equation (13) suggested by [113], is used to generate random  where γ = 0.05, and p j, 1 is the total processing time of all jobs assigned to the first machine.
The second group of the test instances, the REC benchmark, was initially proposed by [115], which is available in the OR-library. This dataset consists of three instances from the seven configurations made of n = {20, 30, 50, 70} jobs and m = {5, 10, 15, 20} machines, each of which including three distinct examples. The second dataset, therefore, is made of 3 × 7 =21 instances.
Proposed by [116], the TA benchmark is featured by n = {20, 50, 100, 200, 500} jobs and m = {5, 10, 20} machines, resulting in a total of 10 × 12 = 120 test instances, where the processing time is a randomly generated integers from the [1,99] interval. Besides, the extended benchmark (ETA) introduced by [81] is considered to analyze the performance of the algorithms in solving very-large-scale test instances. The extension is made of n = {1000, 1500, 2000} jobs, and m = 20 machines, each of which including ten distinct instances, which result in a total of 10 × 3 = 30 instances. Overall, the extensive numerical analysis in the present study comprises a total of 480 + 21 + 120 + 30 = 651 test instances.

B. PARAMETERS CALIBRATION
Four levels are considered for each of the algorithm parameters to conduct the calibration experiments. Given the search width α ∈ {2, 3, 4, 5} and cluster width β ∈ {200, 400, 600, 800} as the BS parameters, the best parameter combination is selected considering the smallest relative percentage deviation (RPD), calculated by Equation (14). This performance measure is also used for analyzing the final test results, where C Best is the best makespan value using all the solution algorithms.
Random  considered for the calibration test, respectively. Given α × β = 4 × 4 = 16 possible combinations of BS parameters, the performance values for each of the categories mentioned above are recorded in Table 3. Overall, for fixed α, larger β results in better solution quality, and larger α, β results in longer computational time; these are shown in Figure 3. It is also observed that the solution time increases sharply for the β values larger than 400. Overall, α = 2 and β = 600 are decided for the final experimental tests over all the test instances. To limit the computational time of the local search procedure, k = 0.1 and 0.01 are considered in test instances with n ≤ 400 and n ≥ 500 workloads, respectively.

C. TEST RESULTS
The experimental test begins with comparing the performance of the algorithms with the optimum solution obtained VOLUME 8, 2020 by an exact approach. For this purpose, the MILP formulation is solved using LINGO, considering the first two groups of VRF instances, shown in Tables 4-5. Evidently, the optimal solution can be obtained for the problems with n = 10 workloads, where the average solution time is approximately 201.3 seconds, which is significantly longer than that of the proposed heuristics. Considering ten independent runs for each of the instances, it is observed that the BS  algorithm obtains the optimum solutions in 7 out of the 40 cases, while the BS LS algorithm nailed it in a total of 28 attempts.
Given the extensive computational time required to solve the instances with n = 20 jobs using LINGO, the solution time is limited to 3600 seconds, and the obtained results are compared with those of the BS and BS LS algorithms. Offering a significantly shorter computational time when compared to the MILP approach, the BS LS algorithm outperforms both the MILP model and BS algorithms in terms of solution quality in all of the instances. It is observed that when m = 5 machines are considered, the average RPD of the solutions obtained by the BS algorithm is about the same, while BS performs better than the exact approach when a larger number of machines are included in the instance.
The remainder of the VRF benchmark is solved using the BS and BS LS algorithms. The initial evaluation shows that BS LS yields the best solution in all of the experiments. According to the performance overview for the small-scale (S-VRF) and large-scale (L-VRF) instances in Tables 4-5, BS LS outperforms BS in terms of solution quality. However, BS LS requires longer computational time, where, expectedly, the difference in computational time becomes significant VOLUME 8, 2020  when the problem size increases (Figures 4-5). The paired t-test, one-tail test, and the respective p-value shown in Table 6 confirm the superiority of BS LS considering the hard benchmarks developed by [114]. The obtained C max , as well as the average RPD values, are recorded in Tables 15∼21 in the Appendix.     Numerical analysis continues using the REC benchmark for the analysis. The overall performance measures of the BS and BS LS algorithms are provided in Table 7. On this basis, Figure 6 visualizes the effect of the problem size on the performance of the algorithms in terms of the solution quality and the computational time. The statistical analysis in Table 8 confirms that BS LS performs significantly better compared to the BS algorithm. The recorded average RPD and CPU times are provided in Table 9 for reference in the future studies, where BS LS yields the best solution in all of the instances.
As a final step to our experiments, the renowned TA benchmark, and its extended version developed by [81] are used to enrich the numerical analysis. The average RPD and CPU time over the test instances are provided in Table 10 and visualized in Figures 7-8. Taking the statistical results in Table 11 as a proof, the local searchbased approach, BS LS performs better in both the original and extended TA benchmarks. BS LS yields the best solutions in all of the instances recorded in Tables 12-13. An important observation is that the difference in computational time becomes significant for the instances with n = 500 and larger.      Table 14 provides the overall analysis of the average RPD and computational time of the BS and BS LS algorithms for 651 test cases, over four standard datasets. It can be concluded that the BS LS is comparatively a better solution algorithm for solving the F m nwt, r j C max problem concerning solution quality. In terms of the computational time, however, BS LS is more time-extensive when compared to the basic BS. A thorough cost-benefit trade-off would help the practitioners to find the extent to which the local search procedure should be involved in the BS LS applications in solving complex industry-scale problems.

VI. CONCLUSION
This study contributed to the understudied topic of flowshop scheduling with no-wait and time constraints, which have wide applications in the chemical and food industries. Given the intractable nature of the F m nwt, r j C max problem, two VOLUME 8, 2020  effective heuristics are put forward to solve the proposed scheduling extension.
Ignoring the non-promising search directions within a partial enumeration solution approach, the BS heuristic builds on the good partial solutions to search for (near-) optimum solutions. Employing two kinds of filtration rules, it is shown that BS produces comparably good solutions in a slight fraction of time the full-enumeration approach needs to solve  the problem. To avoid early convergence and local optimality and further improve the solution quality, a local search approach was integrated into the BS. It is proven that the proposed improvement to BS makes a more effective heuristic, obtaining better solutions in all of the test instances. Although BS LS is computationally extensive, it produces more dependable outcomes for very-large industry-scale problems, while the computational time can be controlled by limiting the local-search computations. Overall, BS-based solution approaches appeared to be viable options for tackling optimization problems with complex search space.
This research can be extended in several directions. No-wait and time constraints can be simultaneously addressed in other production environments.
Service-oriented performance indicators such as the total number of tardy jobs can be studied to analyze the effect of time constraints in prioritizing unprecedented demands. Including other case-specific production situations like mixed-blocking is another worthwhile research direction to contribute to the development of NWFSP-RT. Considering the freshness of this scheduling extension, we feel that metaheuristic algorithms can be helpful to improve the bestfound solutions to the NWFSP-RT. Finally, applications of the improved BS are promising research directions to pursue in solving other complex optimization problems.

APPENDIX
See Tables 15-21.   TABLE 19. Algorithms performance and the best-found solutions for the VRF301-375 instances. PEI-YU LIN received the master's degree in industrial engineering and management from the National Taipei University of Technology, Taiwan. She is currently an Industrial Engineer at Good Young Company Ltd. Her research interests include the operations research and production scheduling.