An Enhanced Two-Level Metaheuristic Algorithm with Adaptive Hybrid Neighborhood Structures for the Job-Shop Scheduling Problem

For solving the job-shop scheduling problem (JSP), this paper proposes a novel two-level metaheuristic algorithm, where its upper-level algorithm controls the input parameters of its lower-level algorithm. The lower-level algorithm is a local search algorithm searching for an optimal JSP solution within a hybrid neighborhood structure. To generate each neighbor solution, the lower-level algorithm randomly uses one of two neighbor operators by a given probability. The upper-level algorithm is a population-based search algorithm developed for controlling the five input parameters of the lower-level algorithm, i.e., a perturbation operator, a scheduling direction, an ordered pair of two neighbor operators, a probability of selecting a neighbor operator, and a start solution-representing permutation. Many operators are proposed in this paper as options for the perturbation and neighbor operators. Under the control of the upper-level algorithm, the lower-level algorithm can be evolved in its input-parameter values and neighborhood structure. Moreover, with the perturbation operator and the start solution-representing permutation controlled, the two-level metaheuristic algorithm performs like a multistart iterated local search algorithm. The experiment’s results indicated that the two-level metaheuristic algorithm outperformed its previous variant and the two other high-performing algorithms in terms of solution quality.


Introduction
Production scheduling is an important tool for controlling and optimizing workloads in an industrial production system. It is a decision-making process which involves assigning jobs to machines on a timetable. e job-shop scheduling problem (JSP) is one of well-known production scheduling problems. Such a problem is also defined as a much complex optimization problem both in theoretical and practical aspects. e objective of JSP is commonly to find a feasible schedule which completes all jobs by the shortest makespan (note that makespan stands for the length of the schedule). Approximation algorithms, such as [1][2][3][4][5][6], have been developed for JSP since the problem cannot be optimally solved in a reasonable (polynomial) amount of time [7,8]. e two-level metaheuristic algorithm of [9] is one of the high-performing approximation algorithms for JSP. is algorithm is, as its name implies, a combination between its upper-level algorithm (named UPLA) and its lower-level algorithm (named LOLA). Its mechanism is that UPLA controls the input parameters of LOLA, and LOLA then searches for an optimal schedule. Due to successful results of [9], this paper aims at developing an enhanced two-level metaheuristic algorithm for JSP. To do so, the two-level metaheuristic algorithm proposed in this paper has been changed from its original variant [9] in both levels.
e lower-level algorithm proposed in this paper, named LOSAP, is a local search algorithm exploring in a probabilistic-based hybrid neighborhood structure. To generate each neighbor solution, LOSAP randomly uses one from the two predetermined neighbor operators (i.e., the first and second operators) by a preassigned probability of selecting the first operator. A high value of this probability leads the hybrid neighborhood structure to be more likely similar to the first operator's neighborhood structure, and vice versa. Previous successful applications of randomly using one from two different operators have been found in [10,11]. Major differences between LOSAP and LOLA [9] are briefly presented as follows. e LOLA's search ability is mainly based on its special solution space, which is a solution space of parameterized-active schedules. e LOSAP's search ability is, however, mainly based on its hybrid neighborhood structure since it searches in an ordinary solution space of semiactive schedules. LOSAP also has many proposed operators as the options for its perturbation and neighbor operators; most of them were not used in LOLA. In addition, LOSAP uses a different criterion from LOLA on accepting a new best-found solution.
Although LOLA and LOSAP have many differences from each other as above-mentioned, they still share a common weakness. e weakness is that no single combination of inputparameter values performs best for all instances. In other words, a combination of input-parameter values performing best for an instance may not perform as best for another instance. For each instance, each algorithm has a specific best combination of input-parameter values; however, it cannot be foreknown without doing an experiment on the being-considered instance.
is weakness is, in fact, a common weakness of most metaheuristic algorithms. To overcome such a weakness, past researches, e.g., [4,9,[12][13][14][15], developed an upper-level metaheuristic algorithm to control input parameters of a problem-solving metaheuristic algorithm in a lower level. e upper-level algorithm proposed in this paper, named MUPLA, is a modification of UPLA (i.e., the upperlevel algorithm of [9]). MUPLA is a population-based metaheuristic algorithm designed to be a parameter controller for LOSAP. us, its population consists of a number of combinations of the LOSAP's input-parameter values which are evolved (updated) over iterations. For short, let a combination of LOSAP's input-parameter values be called a parameter-value combination. Each parameter-value combination contains specific values of the perturbation operator, the scheduling direction, the ordered pair of two neighbor operators, the probability of selecting a neighbor operator, and the start solution-representing permutation. A major change of MUPLA from UPLA [9] is that each parameter-value combination has a different start solutionrepresenting permutation from those of the others. As a consequence, the MUPLA combined with LOSAP acts as a multistart iterated local search algorithm.
is is a large upgrade because the UPLA combined with LOLA [9] is just an iterated local search algorithm. e remainder of this paper is divided into four main sections. Section 2 provides an overview of the relevant publications of the research topic. Section 3 describes the proposed two-level metaheuristic algorithm in both levels. us, the lower-level algorithm (LOSAP) and the upperlevel algorithm (MUPLA) are described in Sections 3.1 and 3.2, respectively. en, Section 4 presents the results and discussions on the evaluation of the two-level metaheuristic algorithm's performance. Section 5 finally summarizes the findings and recommendations.

Literature Review
e job-shop scheduling problem (JSP) comes with n given jobs J 1 , J 2 ,. . ., J n and m given machines M 1 , M 2 ,. . ., M m . Each job J i is composed of a sequence of m given operations O i1 , O i2 ,. . ., O im as a chain of precedence constraints. To process each job J i , O ij (where j = 1, 2, . . ., m − 1) is defined as an immediate-preceding operation of O ik (where k = j + 1); thus, O ij must be finished before O ik can start. In addition, each operation must be processed on a preassigned machine with a predetermined processing time. Each machine cannot process more than one operation at a time, and it cannot be stopped or paused during processing an operation. All n given jobs arrive at the time 0, and all m given machines are also available at the time 0. A schedule is feasible if it completely allocates all n jobs under all the given constraints. An optimal schedule is a feasible schedule which minimizes the makespan, i.e., the total amount of time required to complete all jobs. Excellent reviews about JSP are available in [7,8,16,17].
In JSP, feasible schedules can be alternatively constructed in forward or backward (reverse) directions. A forward schedule is a schedule constructed in the forward direction, while a backward (reverse) schedule is a schedule constructed in the backward direction. In other words, a forward schedule is a schedule in which all jobs J i (where i = 1, 2, . . ., n) are constructed forward from O i1 to O im , while a backward schedule is a schedule in which all jobs J i are constructed backward from O im to O i1 . Although the forward scheduling is commonly used for the makespan criterion, the backward scheduling as an alternative has been applied in many researches, e.g., [18][19][20]. Besides the schedule's classification based on the scheduling directions, the feasible schedules can be classified based on their allowable delay times, e.g., semiactive schedules and active schedules [17,21]. A feasible schedule is defined as a semiactive schedule if no operation can be started earlier without altering an operation sequence on any machine. A semiactive schedule is then defined as an active schedule if no operation can be started earlier without delaying any other operation or violating any precedence constraint.
Many approximation algorithms have been developed based on metaheuristic algorithms for solving JSP. Iterated local search, a well-known type of metaheuristic algorithms, has also been applied for JSP [22,23]. In general, an iterated local search algorithm is a single-solution-based local search technique which can search for a global optimal solution. During an exploration, it uses a neighbor operator repeatedly to find a local optimum and then uses a perturbation operator to escape the just-found local optimum (note that a perturbation operator stands for an operator that generates a new start solution by largely modifying a found local optimal solution [24,25]). It has also been found that some researches such as [26][27][28] enhanced their iterated local search algorithms by adding multistart properties.
In the iterated local search and related algorithms, there are three operators commonly used as a neighbor operator and a perturbation operator. ese three operators are the common swap, insert, and inverse operators [29]. Some iterated local search algorithms, such as [30,31], use the 2 Complexity common swap operator or the common insert operator multiple times as their perturbation operators. e definitions of the three common operators are given as follows (let u and v are random integers from 1 to the number of all members in the permutation, and v ≠ u): (i) e common swap operator is to swap between the two members in the u th and v th positions of a permutation. (ii) e common insert operator is to remove a member from the u th position of a permutation and then insert it back at the v th position. (iii) e common inverse operator is to inverse the sequence of all members from the u th to v th positions of a permutation.
e two-level metaheuristic algorithm of [9] can be classified as an adaptive iterated local search algorithm for JSP. Its upper-level algorithm (named UPLA) controls the input parameters of its lower-level algorithm, and its lowerlevel algorithm (named LOLA) then searches for an optimal schedule. us, the two-level metaheuristic algorithm can adapt itself for every single JSP instance. e development of the two-level metaheuristic algorithm of [9] followed in the successes of the previous researches of [4,[12][13][14][15] in using a metaheuristic algorithm to control parameters of another metaheuristic algorithm.
In [9], LOLA is a local search algorithm exploring in a solution space of parameterized-active schedules. Its input parameters (i.e., an acceptable idle-time limit, a scheduling direction, a perturbation operator, and a neighbor operator) are controlled by UPLA. UPLA is a population-based metaheuristic algorithm searching in a real-number search space. In this view's point, UPLA is similar to the other population-based algorithms, such as particle swarm optimization [32], differential evolution [33], fish swarm [34], and cuckoo search [35]. However, the evolving procedure of the UPLA's population is different from those of the others mentioned. e UPLA's population consists of the combinations of input-parameter values of LOLA. For a parameter-value combination, each parameter's value is iteratively changed by a sum of two changeable opposite-direction vectors. e first vector's direction is toward the memorized best-found value, whereas the second vector's direction is away from. e magnitudes of these two vectors are generated randomly between zeros to their given maximum values.
e first vector's maximum magnitude (0.05) is usually larger than the second vector's maximum magnitude (0.01). However, if the parameter's value equals the memorized best-found value, the maximum magnitudes of both vectors then equal the same value (0.01).

Method
Section 3 describes the procedure of the proposed two-level metaheuristic algorithm in both levels. In this section, the lower-level algorithm (LOSAP) is described in Section 3.1, and the upper-level algorithm (MUPLA) is described in Section 3.2.

Proposed Lower-Level Algorithm.
e lower-level algorithm proposed in this paper, named LOSAP, is a local search algorithm searching for an optimal solution in a probabilisticbased hybrid neighborhood structure. Its framework is similar to those of the other local search algorithms [36][37][38][39]. However, its neighborhood structure is generated based on the two predetermined operators, i.e., the first and second neighbor operators. By a given probability, LOSAP randomly uses one from the two predetermined operators in order to generate a neighbor solution-representing permutation. is means that based on the given probability, LOSAP can switch between the two given operators anytime during its exploration. In this paper, many operators are proposed as the options for being the LOSAP's neighbor operators (note that the successes of randomly using one from two neighbor operators have been found in different algorithms, e.g., [10,11]).
LOSAP generates a hybrid neighborhood structure between the first operator's neighborhood structure and the second operator's neighborhood structure. e hybridization is controlled by the probability of selecting the first neighbor operator as a LOSAP's input parameter (the probability of selecting the second neighbor operator, as a complement, is the unity minus the probability of selecting the first operator). Note that the higher the probability of selecting the first operator, the more likely the hybrid neighborhood structure is like the first operator's neighborhood structure. It equally means that the lower the probability of selecting the first operator, the more likely the hybrid neighborhood structure is like the second operator's neighborhood structure. At boundaries, the mentioned probability's values of 1.00 and 0.00 make the hybrid neighborhood structure be identical to the first operator's neighborhood structure and the second operator's neighborhood structure, respectively. e probability of selecting the first neighbor operator is not the only LOSAP's input parameter. LOSAP has total five input parameters consisting of the perturbation operator, the scheduling direction, the ordered pair of the first and second neighbor operators, the probability of selecting the first neighbor, and the start operation-based permutation. LOSAP provides many options for setting each input parameter in order that LOSAP with proper parameter values can perform well for every single instance. Below, Section 3.1.1 describes Algorithm 1 as the decoding procedure used by LOSAP. In detail, Algorithm 1 is the procedure of transforming each solutionrepresenting permutation generated by LOSAP into a semiactive schedule. Section 3.1.2 then describes Algorithm 2 as the procedure of LOSAP.
3.1.1. Decoding Procedure. LOSAP searches in a solution space of semiactive schedules, and it uses an operation-based permutation [40][41][42][43] to represent a semiactive schedule. An operation-based permutation has been used to represent a semiactive schedule in many researches, such as [20,43,44]. For an n-job/m-machine instance, an operation-based Step 1. Receive an operation-based permutation and a scheduling direction (forward or backward) from LOSAP (Algorithm 2).
Step 2. If the scheduling direction received in Step 1 is forward, then the precedence relations of the operations of each job are unchanged. However, if it is backward, then the precedence relations of the operations must be reversed by using Steps 2.1 and 2.2. Step Step 2.2. Assign the precedence relations of all operations O ij (where j � 1, 2, . . ., m) of each job J i in ascending order of j values. us, O i1 must be finished before O i2 can start, O i2 must be finished before O i3 can start, and so on.
Step 3. If the scheduling direction received in Step 1 is forward, then the operation-based permutation received in Step 1 is unchanged. However, if it is backward, then the order of all members in the operation-based permutation must be reversed. For example, the permutation (3, 2, 3, 1, 1, 2) with a backward scheduling direction must be changed into (2, 1, 1, 3, 2, 3).
Step 4. Transform the permutation taken from Step 3 by changing the number i in its j th occurrence (from left to right) into the operation O ij (i � 1, 2, . . ., n and j � 1, 2, . . ., m). For example, the permutation (3, 2, 3, 1, 1, 2) must be transformed into (O 31 Step 5. Transform the permutation taken from Step 4 into a semiactive schedule by using Steps 5.1 to 5.3.
Step 5.1. Let Φ t be the partial schedule of the t scheduled operations. us, Φ 0 is empty. Now, let t ⟵ 1.
Step 5.2. Let O L ⟵ the leftmost as-yet-unscheduled operation in the permutation. en, create Φ t by scheduling O L into Φ t-1 at its earliest possible start time on its preassigned machine (the earliest possible start time of O L is the maximum between the finished time of its immediate-preceding operation in its job and the finished time of the current last-scheduled operation on its machine).
Step 5.3. If t < mn, then let t ⟵ t + 1 and repeat from Step 5.2. Otherwise, go to Step 6.
Step 6. If the scheduling direction received in Step 1 is forward, then return Φ mn (which is a completed forward semiactive schedule) as the final result. However, if it is backward, then modify Φ mn to satisfy the original precedence relations by using Steps 6.1 and 6.2. Step Step 6.2. Turn the schedule modified from Step 6.1 back to front in order that the last-finished operation in the schedule becomes the first-started operation, and so on. After that, let the schedule be started at the time 0. en, return the schedule modified in this step (which is a completed backward semiactive schedule) as the final result. ALGORITHM 1: e procedure of decoding an operation-based permutation into a semiactive schedule.
Step 2. Generate an initial P 0 by using PTBT on P.
Step 3. Execute Algorithm 1 by inputting SD and P 0 in order to receive S 0 .
Step 4. Find a local optimal schedule by using Steps 4.1 to 4.4.
Step 4.2. Randomly generate u from U[0, 1). If u ≤ PROB, then generate P 1 by using NO1 on P 0 ; otherwise, generate P 1 by using NO2 on P 0 .
Step 4.3. Execute Algorithm 1 by inputting SD and P 1 in order to receive S 1 .
Step 4.4. Update P 0 , S 0 , and t L by using Steps 4.4.1 to 4.4.3.
Step 5. Return P 0 and S 0 as the final (best-found) operation-based permutation and the final (best-found) schedule, respectively. ALGORITHM 2: e procedure of LOSAP. 4 Complexity maximum between the finished time of its immediate-preceding operation in its job and the finished time of the current last-scheduled operation on its machine. Algorithm 1 is the solution-decoding procedure used by LOSAP. As the options, it can transform an operation-based permutation into a forward semiactive schedule (i.e., a semiactive schedule constructed by the forward direction) or a backward semiactive schedule (i.e., a semiactive schedule constructed by the backward direction). us, Algorithm 1 requires assigning values of the two input parameters, i.e., an operation-based permutation and a scheduling direction. Remind that m and n represent the number of all machines and the number of all jobs, respectively, in the being-considered JSP instance.
us, mn (i.e., m multiplied by n) represents the number of all operations.
As mentioned above, Algorithm 1 is the solutiondecoding method used by LOSAP, and Algorithm 2 in Section 3.1.2 is the LOSAP's procedure. us, it can be said that Algorithm 1 is a component of Algorithm 2.
3.1.2. LOSAP Procedure. LOSAP, as shown in Algorithm 2, has the five input parameters whose values need to be assigned, i.e., PTBT, SD, PNO ≡ (NO1, NO2), PROB, and P. e definitions of these LOSAP's input parameters are given below: (i) PTBT and P stand for the perturbation operator and the start operation-based permutation, respectively. LOSAP uses PTBT on P in order to generate its initial best-found operation-based permutation. (ii) SD stands for the scheduling direction (forward or backward) of all solutions generated by LOSAP. If SD is selected to be forward, LOSAP transforms each generated operation-based permutation into a forward schedule. Otherwise, LOSAP transforms each permutation into a backward schedule. (iii) PNO ≡ (NO1, NO2) is the ordered pair of the first neighbor operator (called NO1) and the second neighbor operator (called NO2). (iv) PROB is the probability of selecting the first neighbor operator (NO1). Consequently, the probability of selecting the second neighbor operator (NO2) is the unity minus PROB.
Besides the definitions of the input parameters mentioned above, the other abbreviations used in LOSAP (i.e., Algorithm 2) are defined below: (i) P 0 stands for the current best-found operationbased permutation. As mentioned above, the initial P 0 is generated by using PTBT on P. (ii) S 0 , which is decoded from P 0 , stands for the current best-found schedule. In addition, Makespan (S 0 ) stands for the makespan of S 0 . (iii) P 1 , which is generated by using NO1 or NO2 on P 0 , stands for the current neighbor operation-based permutation. (iv) S 1 , which is decoded from P 1 , stands for the current neighbor schedule. In addition, Makespan (S 1 ) stands for the makespan of S 1 .
(v) m and n are the number of machines and the number of jobs, respectively, in the being-considered JSP instance. us, mn is the number of operations. e LOSAP's procedure given in Algorithm 2 is efficient but simple. In brief, LOSAP starts its procedure by using PTBT on P to generate an initial P 0 ; after that, LOSAP transforms P 0 into S 0 . en, LOSAP starts its repeated loop by using PROB to randomly select either NO1 or NO2. LOSAP uses the selected operator (i.e., either NO1 or NO2) on P 0 to generate P 1 , and LOSAP then transforms P 1 into S 1 . If the criterion of accepting a new best-found permutation is satisfied, then LOSAP updates P 0 ⟵ P 1 and S 0 ⟵ S 1 . Finally, LOSAP determines whether to continue the repeated loop's next round or stop its procedure. Note that S 0 and S 1 are always generated in the forward direction if SD is selected as forward at the beginning; otherwise, they are always generated in the backward direction.
In LOSAP, there are many optional operators for PTBT and PNO. ese optional operators are modified from the common swap, insert, and inverse operators [29] (the definitions of the common operators are given in Section 2). In detail, the five LOSAP's options for PTBT consist of the nmedium swap operator, the n-large swap operator, the nmedium inverse operator, the n-large insert operator, and the n-medium insert operator. LOSAP also provides the four options for PNO ≡ (NO1, NO2), consisting of (1-small inverse, 1-medium insert), (1-large swap, 1-large insert), (1medium swap, 1-medium insert), and (1-small swap, 1-small insert). ese optional operators are defined as follows. e number in front of the hyphen sign (-) indicates the number of repeated uses of the operator mentioned in back of the hyphen sign. For example, the 1-small swap operator is to use the small swap operator once on a permutation, and the n-medium inverse operator is to use the medium inverse operator n times on a permutation.
In addition to the above paragraph, the words small, medium, and large in the names of the optional operators are used to restrict the value of v in its distance from u as explained below (remind that u is a random integer within [1, mn]): (i) For all operators with small, v is a random integer within [u − 4, u + 4] (note that the small swap in [45] means the swap of two adjacent members, and it thus differs from the small swap in LOSAP). (ii) For all operators with medium, v is a random integer within [u − (mn/5), u + (mn/5)]. (iii) For all operators with large, v is a random integer within [1, mn]. It means that the operators with large are identical to the common operators.
After generating v successfully, its value must be verified whether it can be used or not. If the generated value of v is outside of [1, mn] or equal to u, then it must be randomly regenerated within the given same range.
is procedure must be repeated until receiving the value of v within [1, mn] and unequal to u.

Complexity
Remind that an addition of more optional operators into LOSAP may not always give a benefit for the two-level metaheuristic algorithm. Moreover, it sometimes makes the two-level metaheuristic algorithm harder to find a better solution. As shown in Algorithm 2, LOSAP does not include the n-small swap, n-small insert, and n-large inverse operators as the options for PTBT. e reason of excluding the n-small swap and n-small insert operators is that they make a too-small change into P; by using them, the two-level metaheuristic algorithm can hardly escape a local optimum. In contrast, a change from the n-large inverse operator is almost as large as a change from a re-initialization. For PNO, the 1-medium inverse and 1-large inverse operators are not used as the options because, as neighbor operators, they make a too-large change into P 0 .
LOLA [9] and LOSAP both are lower-level algorithms of their own two-level metaheuristic algorithms; however, there are many differences between them. e LOLA's search ability is mainly based on its solution space of parameterizedactive schedules. By contrast, LOSAP uses an ordinary solution space of semiactive schedules; thus, its search ability is mainly based on its probabilistic-based hybrid neighborhood structure. Most optional perturbation and neighbor operators of LOSAP are different from those of LOLA; the n-large insert operator is the only optional perturbation operator found in both LOLA and LOSAP. Another difference between LOLA and LOSAP is in their criteria of accepting a new best-found solution. LOLA accepts only a better neighbor solution, while LOSAP accepts a nonworsening neighbor solution. LOSAP uses this acceptance criterion to escape from a shoulder (i.e., a flat area of search space adjacent to a downhill edge [46]). In addition, LOSAP will not reset t L to 0 for any sideways move in order to avoid an endless loop when finding a flat local minimum.

Proposed Upper-Level
Algorithm. MUPLA is the upperlevel algorithm of the proposed two-level algorithm. e purpose of MUPLA is to evolve the LOSAP's input-parameter values so that LOSAP can return its best performance on every single JSP instance. MUPLA is a population-based search algorithm specifically developed for being a parameter tuner. e MUPLA's population contains N combinations of the LOSAP's input-parameter values, i.e., C 1 (t), C 2 (t),. . ., C N (t). In short, let a parameter-value combination stand for a combination of the LOSAP's input-parameter values. In the population, each parameter-value combination is adjusted over iterations by the MUPLA's specific evolutionary process.
Let C i (t) ≡ (c 1i (t), c 2i (t), c 3i (t), c 4i (t), p i (t)) represent the i th parameter-value combination (where i = 1, 2, . . ., N) in the MUPLA's population at the t th iteration. ese c 1i (t), c 2i (t), c 3i (t), and c 4i (t) are real numbers representing the values of the perturbation operator (PTBT), the scheduling direction (SD), the ordered pair of the first and second neighbor operators (PNO), and the probability of selecting the first neighbor operator (PROB) of LOSAP, respectively. In addition, p i (t) represents the start operation-based permutation of LOSAP. Note that p i (t) is an operationbased permutation, not a real number like others. Table 1 presents the transformation from c 1i (t), c 2i (t), c 3i (t), c 4i (t), and p i (t) into the values of PTBT, SD, PNO, PROB, and P of LOSAP, respectively. e abbreviations used in MUPLA (i.e., Algorithm 3) are defined below: (i) Score (C i (t)) stands for the performance score of C i (t). Note that the lower the performance score, the better the performance. Between any two parameter-value combinations, the combination with a lower performance score is the better one. (ii) P fi (t) stands for the final (best-found) operationbased permutation of the LOSAP using the inputparameter values decoded from C i (t). (iii) S fi (t) stands for the final (best-found) schedule of the LOSAP using the input-parameter values decoded from C i (t). In addition, Makespan (S fi (t)) stands for the makespan of S fi (t). (iv) C best ≡ (c 1best , c 2best , c 3best , c 4best , p best ) stands for the best parameter-value combination ever-found by the population. In addition, Score (C best ) stands for the performance score of C best (note that C best in definition is similar to the global best position of PSO [32]). (v) S best stands for the best schedule ever-found by the population.
e procedure of MUPLA is presented in detail in Algorithm 3. In addition, it is also presented in a form of flow chart in Figure 1. In brief, MUPLA starts its procedure by assigning t ⟵ 1 and Score (C best ) ⟵ + ∞. en, MUPLA randomly generates C i (t), where i = 1, 2, . . ., N. After that, MUPLA processes its repeated loop as follows. For each C i (t), MUPLA decodes it into the LOSAP input-parameter values and then runs the LOSAP using these input-parameter values to receive P fi (t) and S fi (t); then, MUPLA assigns Score (C i (t)) ⟵ Makespan (S fi (t)). If MUPLA finds any C i (t) whose Score (C i (t)) is less than or equal to Score (C best ), then it updates C best ⟵ C i (t), Score (C best ) ⟵ Score (C i (t)), and S best ⟵ S fi (t). After that, MUPLA generates C i (t + 1), where i = 1, 2, . . ., N, by using its specific evolutionary process (as shown in Step 3 of Algorithm 3), and it then assigns t ⟵ t + 1. Finally, MUPLA determines whether to continue its repeated loop's next round or stop its procedure.
A main difference of MUPLA from the other population-based algorithms such as [32][33][34][35] is in its specific evolutionary process (i.e., the procedure of adjusting its population) given by the equation in Step 3.3. In the equation in Step 3.3, each c ji (t + 1) is updated from c ji (t) by the sum of two opposite-direction vectors. e first vector is toward the best-found value, whereas the second vector is away from. e equation in Step 3.3 used in MUPLA is slightly modified from that in UPLA [9]. e modification is to reduce the first vector's maximum magnitude by a half (from 0.05 to 0.025) in cases that c ji (t) differs from c jbest . e purpose of this modification is to delay the population from getting stuck in a local optimum. In a preliminary experiment, the proposed two-level metaheuristic algorithm performed better on average after reducing the first vector's maximum magnitude as mentioned. 6 Complexity Besides the modification in the equation in Step 3.3, there are larger other changes of MUPLA from UPLA [9]. One change is that, unlike UPLA, each C i (t) of the MUPLA's population includes a start operation-based permutation p i (t). By this change mentioned, MUPLA can add a multistart property into LOSAP. Consequently, the combination of MUPLA and LOSAP becomes a multistart iterated local search algorithm.
is is a large upgrade because the combination of UPLA and LOLA in [9] is just an iterated local search algorithm. Another change is in its criterion of accepting a new best-found parameter-value combination. UPLA accepts only a better parameter-value combination, while MUPLA accepts a nonworsening parameter-value combination.

Results and Discussions
e performance of the proposed two-level metaheuristic algorithm was evaluated via the experiment conducted in this research. Section 4 then presents results of the Table 1: Transformation from C i (t) into the LOSAP's input-parameter values.
Step 2. Evaluate Score (C i (t)), and update C best and S best by using Steps 2.1 to 2.5.
Step 2.2. Transform C i (t) into the values of PTBT, SD, PNO, PROB, and P of LOSAP by the relationships shown in Table 1.
Step 2.3. Execute Algorithm 2 (LOSAP) by inputting the PTBT, SD, PNO, PROB, and P taken from Step 2.2 in order to receive P fi (t) and S fi (t). en, let Score (C i (t)) ⟵ Makespan (S fi (t)).
Step 2.4. If Score (C i (t)) ≤ Score (C best ), then let C best ⟵ C i (t), Score (C best ) ⟵ Score (C i (t)), and S best ⟵ S fi (t).
Step 2.5. If i < N, then let i ⟵ i + 1 and repeat from Step 2.2; otherwise, go to Step 3.
Step 3.3. If t mod 25 � 0, then randomly generate c ji (t + 1) from U[0, 1) (where j � 1, 2, 3, 4); otherwise, generate c ji (t + 1) by the following equation. Let u 1 and u 2 be randomly generated from U[0, 1). c ji (t + 1) � experiment. In this section, let MUPLA stand for the whole two-level metaheuristic algorithm (i.e., the MUPLA combined with LOSAP). e reason is that MUPLA uses LOSAP as its component when it solves JSP. For a comparison purpose, the MUPLA's results were compared to those of TS, GA, and UPLA taken from [5,6], and [9], respectively. Let TS stand for the taboo search algorithm developed by Nowicki and Smutnicki [5], and let GA stand for the genetic algorithm developed by Piroozfard et al. [6]. In addition, let UPLA in this section stand for the UPLA combined with LOLA. e performance of MUPLA was evaluated on 53 benchmark instances. ese 53 instances consisted of the ft01, ft10, and ft20 instances from [47], the la01 to la40 instances from [48], and the orb01 to orb10 instances from [49]. e 53 mentioned instances were chosen because they covered all instances used by [5,6,9] in their experiments. In detail, the TS's performance was evaluated in [5] on the 43 first-mentioned instances, i.e., ft06, ft10, ft20, and all 40 la instances. e GA's performance was evaluated in [6] on the 42 instances, i.e., ft06, ft20, and all 40 la instances. In addition, the UPLA's performance was evaluated in [9] on all the 53 instances.
ere were two main objectives of the experiment. e first objective was to evaluate the best performance of which MUPLA was capable on its solution quality without a computational time's limitation. To achieve this objective, MUPLA was run on an extremely long computational time (i.e., 5,000 iterations in this paper) for each trial unless it could find the known optimal solution. e reason was to ensure that MUPLA already reached the best solution as possible as it could. e second objective of the experiment was to evaluate the MUPLA's search performance over iterations. To do so, the solution's convergence rates of MUPLA were plotted. A finding from the solution's convergence rates was to identify a proper maximum iteration for MUPLA. It was very important because a proper maximum iteration could balance the solution quality and the consumed computational time. Remind that the higher the number of iterations used, the higher the computational time consumed. e parameter settings of MUPLA used in the experiment are shown below: (i) e population of MUPLA consisted of three parameter-value combinations (i.e., N = 3 in Algorithm 3). (ii) e stopping criterion of MUPLA was that either the 5,000 th iteration (i.e., t = 5,000 in Algorithm 3) was reached or the known optimal solution value was found (note that solution and solution value in this paper represent schedule and makespan, respectively). (iii) MUPLA was coded in C# and executed on an Intel ® CoreTM i5 CPU processor M580 @ 2.67 GHz with RAM of 6 GB (2.3 GB usable). (iv) For each instance, MUPLA was executed five trials with different random-seed numbers.
e results of the experiment using the above settings are presented in Table 2. In the table, the MUPLA's results on each instance consisted of the best-found solution value over five trials, the average of the first trial's best-found solution value to the fifth trial's best-found solution value, the average number of used iterations, and the average computational time consumed. In the purpose of comparison, Table 2 presents the best-found solution values of TS, GA, and UPLA taken from [5,6], and [9], respectively. Remind that in the table, the TS's results are shown on only the ft06, ft10, ft20, and all 40 la instances, and the GA's results are shown on only the ft06, ft20, and all 40 la instances. e abbreviations used in Table 2 are defined as follows: (i) Instance and Opt represent the name of each instance and its known optimal solution value given by literature (e.g., [2]), respectively. (ii) In each instance, Best represents the best-found solution value of each algorithm. Best of MUPLA was taken from the five trials in this experiment.
Execute LOSAP to receive P fi (t) and S fi (t).
Transform C i (t) into LOSAP's input-parameter values.
Update C best and S best .

No
Let i ← 1.

No
Yes i = N? Stop?
Based on the results in Table 2, Section 4.1 evaluates the best search performance of which MUPLA was capable without a computational time's limitation. Section 4.2 then evaluates the MUPLA's search performance over iterations.

Evaluation on Search Performance without Computational Time's Limitation.
As just mentioned, Section 4.1 aims at revealing the best search performance of which MUPLA was capable on its solution quality without a computational time's restriction. To achieve the aim, MUPLA was executed on an extremely long computational time (i.e., 5,000 iterations) unless MUPLA could find the known optimal solution. e purpose of this setting was to ensure that MUPLA already reached its best solution as possible as it could. In the performance evaluation, the % BDs of MUPLA were compared to the %BDs of TS [5], GA [6], and UPLA [9] shown in Table 2. Remind that the %BDs of TS, GA, and UPLA were calculated from the best-found solution values given by their original articles. e use of these %BDs of TS, GA, and UPLA can be defined as a limitation of this experiment because TS, GA, and UPLA might find improving solutions if they were run on more computational times than those used in their articles. On the other hand, they might not find any improving solution although much more computational times were provided; it usually happens to any metaheuristic algorithm when its search gets stuck in a local optimum.
Discussions on the %BDs shown in Table 2 are given hereafter along with using one-sided paired t-tests to compare the mean %BDs. Let the significance level (α) be equal to 0.1. As shown in Table 2, MUPLA performed better than TS [5] in terms of %BD (remind that the performances were compared on only the 43 instances). Of the total 43 instances, the average %BD of MUPLA was 0.03%, while the average %BD of TS was 0.06%. MUPLA returned better %BDs than TS on five instances (i.e., la21, la24, la27, la37, and la40), while TS returned a better %BD than MUPLA on only one instance (i.e., la29). On the 37 remaining instances, they both returned equal %BDs. In addition, there were 41 instances where MUPLA returned %BDs of 0.00%, while there were only 37 instances where TS returned % BDs of 0.00%. It means that from all 43 instances, MUPLA found optimal solutions on 41 instances, while TS found optimal solutions on only 37 instances. With α = 0.1, an one-sided paired t-test indicated that the mean %BD of MUPLA was significantly better than the mean %BD of TS (p value = 0.066).
MUPLA also outperformed GA [6] in terms of %BD (remind that the performances were compared on only the 42 instances). Over the total 42 instances, the average %BD of MUPLA was 0.03%, while the average %BD of GA was 0.35%. MUPLA returned better %BDs than GA on 13 instances, and they both returned equal %BDs on the 29 remaining instances. is means that GA could not return a better %BD than MUPLA on any instance. Moreover, MUPLA returned %BDs of 0.00% on 40 instances, while GA returned %BDs of 0.00% on only 29 instances. is means that from the 42 instances, MUPLA found the optimal solutions on 40 instances, while GA found the optimal solutions on only 29 instances. In detail, MUPLA failed to find optimal solutions on la29 and la40, while GA failed on la29, la40, and the 11 other instances. Although they both failed on la29 and la40, MUPLA performed much better on these two instances. On la29 and la40, respectively, the MUPLA's %BDs were 0.95% and 0.16%, while the GA's % BDs were 2.34% and 0.90%. An one-sided paired t-test with α = 0.1 indicated that the mean %BD of MUPLA was significantly better than the mean %BD of GA (p value = 0.001).
Based on %BDs in Table 2, MUPLA outperformed UPLA [9]. From the total 53 instances, the average %BD of MUPLA was 0.02%, while the average %BD of UPLA was 0.24%. MUPLA returned better %BDs than UPLA on 13 instances, and they both returned the same %BDs on 40 instances. It means that there were no instances where UPLA returned better %BDs than MUPLA. In addition, MUPLA returned the %BDs of 0.00% on 51 instances, while UPLA returned the % BDs of 0.00% on only 40 instances. It means that MUPLA could find optimal solutions on 51 instances, while UPLA could find optimal solutions on only 40 instances. In detail, MUPLA failed to find optimal solutions on only la29 and la40, while UPLA failed on la29, la40, and the 11 other instances. Based on the results from Table 2, an one-sided paired t-test with α = 0.1 indicated that the mean %BD of MUPLA was significantly better than the mean %BD of UPLA (p value = 0.002).
A summary of the above discussions on %BDs is given below: In discussion on %ADs, Table 2 shows that over all 53 instances, MUPLA returned %ADs of 0.00% on 45 instances. It means that MUPLA found the optimal solution by every single trial for the 45 instances. On the eight remaining instances, %ADs of MUPLA were also very low. Notice that the highest %AD, belonging to la29, was only 1.02%. is can be interpreted that MUPLA could perform very well in overall, not only in its best trial. When compared with % BDs, %ADs were equal or almost equal to %BDs for all instances. is emphasizes that MUPLA kept its high performance with good consistency from one trial to another.

Evaluation on Search Performance over Iterations.
In this section, the solution's convergence rates of MUPLA were plotted in order to evaluate the MUPLA's search performance over iterations (remind that a solution's convergence rate is usually presented by a scatter plot showing the relationship between the solution quality and the number of iterations used). As shown in Figures 2 and 3 respectively, % AD-over-iteration plots (i.e., the plots of %ADs over iterations) and %BD-over-iteration plots (i.e., the plots of %BDs over iterations) were drawn for the eight instances, i.e., la21, la24, la27, la29, la38, la40, orb2, and orb6. ese eight instances were selected because their %ADs in Table 2 were greater than 0% (it means that for each mentioned instance, at least one trial of MUPLA could not find the known optimal solution). e total number of iterations in each figure was only 1,500 iterations (not 5,000 iterations) because %ADs and %BDs, in general, were changed hardly after the 1,500 th iteration.
To simplify the analysis, the pieces of information in Figures 2 and 3 were combined altogether and then put into Figure 4. At each iteration, the %ADs and %BDs on eight instances from Figures 2 and 3 were averaged into the avg % AD and avg %BD, respectively, in Figure 4. ey resulted in the avg-%AD-over-iteration plot and the avg-%BD-over-iteration plot in Figure 4. For a comparison purpose, Figure 4 also presents the plots of TS's final avg %BD, GA's final avg % BD, and UPLA's final avg %BD. ese plots represent the average %BDs of the final solutions on eight instances of TS, GA, and UPLA (i.e., 0.31%, 1.38%, and 1.41%, respectively). ey were used to identify the lowest numbers of iterations of which MUPLA returned better solutions than the final solutions of TS, GA, and UPLA.
In Figure 4, the avg-%AD-over-iteration and avg-% BD-over-iteration plots could be divided into three periods based on their reduction rates. e first period roughly started from the first iteration to the 160 th iteration, where the avg %AD and avg %BD were reduced sharply. e second period roughly started from the 160 th iteration to the 700 th iteration, where the avg %AD and avg %BD were reduced relatively quickly but slower than the first period's rate. e third period roughly started from the 700 th iteration onwards, where the avg %AD and avg %BD were reduced very slowly. Since the higher number of iterations results in the higher computational time, the proper maximum iteration can help MUPLA to balance the solution quality and the computational time. Based on the two given plots, the MUPLA's maximum Complexity iteration should be set within the 160 th to 700 th iterations. At extremes, the 160 th iteration should be used when highly concerned on the computational time, while the 700 th iteration should be used when highly concerned on the solution quality. For any other point within the 160 th to 700 th iterations, the higher possibility to find a better  For more analysis on Figure 4, the avg-%AD-over-iteration and avg-%BD-over-iteration plots were compared to the final avg %BDs of TS, GA, and UPLA. For each algorithm, the final avg %BD is the average of the %BDs of the final solutions on eight instances taken from Table 2. In Figure 4, the avg-%BD-over-iteration plot of MUPLA was lower than the final avg %BDs of UPLA, GA, and TS at the eighth, ninth, and 484 th iterations, respectively. It means that the MUPLA's best trial (among the total five trials) at its eighth, ninth, and 484 th iterations provided better solutions than the final solutions of UPLA, GA, and TS, respectively. Moreover, the avg-%AD-over-iteration plot was lower than the final avg %BDs of UPLA and GA at the 35 th and 39 th iterations, respectively. It can be interpreted roughly that MUPLA could usually return the better solutions than the final solutions of UPLA and GA within the first 40 iterations. e findings mentioned in this paragraph fulfilled those of the previous paragraph. One finding was that MUPLA could return acceptable-quality solutions within its first 10 iterations and then high-quality solutions within its first 500 iterations. Moreover, MUPLA was still able to find improving solutions after the 500 th iteration.
In addition, Figure 4 also reveals an enhancement percentage of the LOSAP's search performance received from MUPLA. Remind that at the MUPLA's first iteration (as the start), the LOSAP's input-parameter values were generated full-randomly. It means that at the MUPLA's first iteration, LOSAP performed its own search performance without any support from the upper-level algorithm yet. In Figure 4, the avg %AD at the first iteration was 4.80%. After passing through the MUPLA's evolutionary process to the 100 th iteration, the avg %AD became 1.00%. At the 500 th iteration, the avg %AD was then reduced to 0.56%. e improvement of the avg %AD has been a strong evidence of the enhancement in the LOSAP's search performance received from MUPLA. is can be concluded in the same way when explaining by the avg-%BD-over-iteration plot.

Conclusion
e proposed two-level metaheuristic algorithm is composed of the upper-level algorithm and the lower-level algorithm named MUPLA and LOSAP, respectively. In the upper level, MUPLA is a population-based search algorithm developed for controlling the LOSAP's input parameters. In the lower level, LOSAP is a local search algorithm with a probabilistic-based hybrid neighborhood structure. e work relation between MUPLA and LOSAP is given in brief as follows: MUPLA starts a repeated loop by generating the LOSAP's input-parameter values; then, LOSAP uses these input-parameter values to return its best-found JSP solution. e best-found JSP solution becomes a feedback sent back to MUPLA for improving the LOSAP's input-parameter values. MUPLA then starts the next round of the repeated loop. e MUPLA combined with LOSAP performs like an adaptive multistart iterated local search algorithm. e experiment's results indicated that the proposed two-level metaheuristic algorithm (i.e., the MUPLA combined with LOSAP) outperformed its original variant (i.e., the UPLA combined with LOLA) and the two other high-performing algorithms in terms of solution quality. Based on the convergence rates, the maximum iteration of the two-level metaheuristic algorithm was recommended to be set within the 160 th to the 700 th iterations. A future study should be conducted to enhance the performance of the two-level metaheuristic algorithm for JSP and related other problems. Another interesting future study is to modify MUPLA (uncombined with LOSAP) for being a lower-level algorithm. Consequently, the two-level metaheuristic algorithm using MUPLA in both levels will be possibly developed.

Data Availability
e data used to support the findings of this study are available from the author upon request.

Conflicts of Interest
e author declares that there are no conflicts of interest regarding the publication of this paper.