Due to their complex combinatorial nature, designing metaheuristics for multi-objective scheduling problems is a big challenge. This task generally involves a proper selection of a solution representation scheme that suits both the structure of the problem and the search mechanism of a metaheuristic. However, working on the solution representation using the plain diversification and intensification strategies of a metaheuristic are usually not sufficient for achieving non-dominated solutions that are close enough to the optimal Pareto front. Therefore, hybridization with a suitable local search heuristic is usually employed. Nonetheless, the design of such local search heuristics requires a proper selection of neighborhood search functions that promote moves towards the optimal Pareto front, which is another challenge.
Scatter search is a population-based metaheuristic that was introduced by Glover [
39]. It utilizes problem-specific neighborhood search and solution-recombination functions, which provide intensification and diversification mechanisms, to iteratively update a set of reference solutions, denoted as
. Solution quality and diversity are essential characteristics for solutions to be included in
. In multi-objective optimization problems, the set
can partially include the non-dominated solutions that are found and updated throughout iterations. The first application of the scatter search to multi-objective optimization problems was for the proctor assignment problem by Martí et al. [
40]. For non-linear optimization problems, early competitive results for multi-objective scatter search implementations include Nebro et al. [
41] and Beausoleil [
42].
3.1. Outline of the Developed Metaheuristic
As outlined in
Figure 1, the proposed MOSS approach consists of two phases, namely initialization and improvement. In the initialization phase, the set of reference solutions,
, is constructed as follows. A random solution is generated using the construction heuristic described in
Section 3.2. This is followed by the application of a local search heuristic based on tabu search (TS), which is outlined in
Section 3.5. The new, improved solution is then checked for similarity with other solutions that have been added to the set
in previous iterations. The degree of dissimilarity between two solutions,
and
, is measured using a distance function, denoted as
, which is defined in
Section 3.3. The distance between the new, improved solution and every solution included in
up to the current iteration is evaluated and checked against a predefined distance threshold value,
. The new solution will be added to
only if
is exceeded for all solutions currently in that set. If a number of improved solutions failed to be added to
in
times, the
threshold value is reduced to half of its value. This is done to ensure that the initialization procedure generates the required number of solutions in the reference set (
). This construction approach tries to maintain both the quality and diversity of the initial solutions included in the reference set
.
In the improvement phase, a total of
iterations are conducted. In each iteration, a set of trial solutions, denoted as
, is constructed. This is done by randomly selecting two solutions from the reference set
and then generating a new solution by combining some of their features, as described in
Section 3.4. The combined solution is then improved via TS and used to update the non-dominated set
. Within the steps of TS either during the construction or the improvement phase, an incumbent solution is compared with all found non-dominated solutions that are included in a set denoted as
. If a non-dominated solution is found, it will be added to that set.
The newly generated solutions in set
are then used to update the reference set
. The size of set
is limited to twice the size of the set
, which in turn should not exceed the specified
. The main idea behind this updating mechanism is to maintain a set of diverse and high-quality solutions in the set
, which are used in generating new trial solutions. For the set of non-dominated solutions, no size restriction is imposed. If the generated solutions in set
do not change the solutions in both sets,
and
, a diversification procedure is applied, as explained in
Section 3.6. The following subsections explain the algorithmic details of the proposed MOSS metaheuristic’s main components.
3.2. Solution Representation and Construction
The operations-based (OB) representation scheme is commonly used in the applications of metaheuristic approaches to scheduling problems [
4,
26,
31,
45]. It appears in the form of a permutation that represents the sequence by which operations will be scheduled according to an interpretation algorithm, which generally chooses the earliest start times for the operations based on the given sequence. Despite its successful implementations and reported competitive results in the literature, the OB representation is many-to-one, meaning that more than one solution representation for a specific instance can result in the same generated schedule. As indicated by Abdelmaguid [
4], such redundancy decreases the efficiency of the search operators of a metaheuristic. Alternatively, in this paper, we propose a representation based on vectors of permutations that correspond to the workstations visiting orders for jobs and the processing sequences on machines. It is referred to as a vectors of permutations-based (VPB) representation. The proposed representation is based on the model used in [
29] for the MOSP.
To illustrate the proposed VPB chromosome structure, consider a simple instance of the studied problem for which the parameters are given in
Table 1,
Table 2 and
Table 3. In this instance, there are five workstations, indexed according to
. A machine is labeled as
, where
i is the assigned index of the machine in workstation
w. The set of all machines is denoted as
.
Table 1 presents the structure of the system of the sample instance in terms of the machines that constitute each workstation and the ready times of the machines. There are six jobs to be scheduled, indexed according to
. The assigned priority values for the jobs (
), their release times (
), and the sets of required workstations (
) are listed in
Table 2. Lastly, the processing times, denoted as
, for every job
and machine
, are provided in
Table 3.
The proposed VPB chromosome structure, denoted as
, is constituted of two strings—
and
—and a vector of integer values, denoted as
, representing locations in the string
. The two strings are composed of sub-strings of fixed sizes, and each sub-string is composed of elements that refer to the operations to be scheduled.
Figure 2 represents the chromosome structure of a solution to the considered sample instance of the studied problem. Each operation is represented by a letter that refers to its job label, as given in
Table 2, followed by a number that refers to its processing workstation, as provided in
Table 1. As shown in
Figure 2, the two strings that constitute the chromosome structure include
(the upper string), which corresponds to the orders by which jobs visit the workstations, and
(the lower string), which corresponds to the operations’ processing sequences in the workstations. Meanwhile, the vector
defines locations for separators that divide the workstations’ sub-strings into divisions representing the processing sequences on the machines.
For
, each sub-string corresponds to a job, and its size equals the number of operations of that job. To facilitate the visual identification of sub-strings, alternating grayscale colors are used in
Figure 2. For instance, for the job labeled with the letter E, there are four elements in its sub-string. The order of the elements (operations) given in the job labeled with the letter C, for instance, indicates that job C will visit workstations in the order of 2, followed by 5, followed by 3. Accordingly, the upper string shown in
Figure 2 can be directly interpreted into the vector of permutations expressed as
=
=
=
,
,
,
,
,
, as in [
29]. Having sub-strings with fixed sizes in the proposed chromosome structure entails that moves or exchanges of the elements as part of the search procedure of a metaheuristic will be conducted only within a sub-string; meanwhile, inter sub-string moves or exchanges are not allowed.
In
, fixed-size sub-strings are defined to represent the processing sequences of operations on the machines that belong to workstations. Each sub-string corresponds to a workstation. Since a workstation may contain more than one machine, these sub-strings are divided by what is referred to as machine separators, which are defined in the vector
. If workstation
contains the set
of machines, there will be
machine separators in its sub-string. The processing sequences of the operations on the machines in set
are defined using the divisions of the sub-string of workstation
w. For instance, in
Figure 2, the first division of the sub-string of workstation 4 contains the operation
, which will be processed on the first machine,
. Meanwhile, the second division contains operation
, followed by operation
, which will both be processed on machine
in that order. The processing sequences provided in
Figure 2 for the sample DMOSP instance can be easily interpreted with the vector of vectors of permutations expressed as
=
=
, where
=
=
,
=
=
,
=
=
,
=
=
, and
=
=
, as in [
29]. Similar to the upper string, moves or exchanges as part of the search procedures of a metaheuristic must be conducted only on elements within a sub-string, and inter sub-string moves are not allowed. In addition, moves on the machine separators in the form of changing positions also need to be considered. Such moves allow for changing the sets of assigned operations for the machines in a workstation.
The decoding procedure of the proposed solution-representation scheme is based on assigning the earliest start time for each operation based on the given permutations by both strings. The pseudo-code in Algorithm 1 outlines its main steps.
Algorithm 1 Chromosome decoder for the proposed solution-representation scheme. |
- 1:
function VPB Chromosome Decoder - 2:
Input: DMOSP instance and a VPB chromosome () - 3:
Output: Schedule - 4:
Construct = based on - 5:
Construct = based on and , where - 6:
= - 7:
= - 8:
= first operation in - 9:
= empty schedule - 10:
do - 11:
= - 12:
for every , - 13:
= first operation in - 14:
if and = for then - 15:
= - 16:
Set start time of , = - 17:
Set completion time of , = + - 18:
Add to - 19:
= = - 20:
Remove from - 21:
= first operation in or if is empty - 22:
Remove from - 23:
- 24:
break for loop - 25:
end if - 26:
end for - 27:
while () - 28:
if is empty return - 29:
else return INFEASIBLE SCHEDULE - 30:
end function
|
Algorithm 1 starts by converting the given VPB chromosome into a solution pair, , as demonstrated earlier. Then, three sets of variables are initialized. For every machine, m, the variable refers to the time at which machine m is ready to process its next operation. For every job, j, the variable refers to the time at which the job can start its next operation, while the variable refers to the job’s next operation to be scheduled, based on its given order of visiting the workstations. Then, the main loop of the algorithm constructs the schedule by iteratively scanning the next operation to be scheduled based on the processing order for each machine. If the next operation to be scheduled for job is found to be the same as the next operation to be scheduled for machine m, that operation is added to the schedule after its start and completion times are determined, and the values of its corresponding machine’s and job’s and are updated. The loop continues until all operations are scheduled. The condition of infeasibility is identified, and no schedule is returned if the search loop is terminated with remaining unscheduled operations.
Figure 3 shows the Gantt chart of the schedule generated for the solution representation shown in
Figure 2 when Algorithm 1 is applied. Note that a feasible schedule generated via Algorithm 1 is semi-active, meaning that the start time of an operation cannot be set earlier without changing the processing order within its assigned machine and/or changing its job’s visiting order to its required workstations [
46]. The set of semi-active schedules includes optimal solutions for many objective functions, including the ones considered in this paper. In addition, the following proposition states an important property for the relationship between the input chromosomes of Algorithm 1 and the generated schedules.
Proposition 1. For any instance of the DMOSP, a schedule generated via Algorithm 1 has one and only one input VPB chromosome.
Proof. The proof follows from the observation that a given chromosome
is interpreted to one and only one solution pair,
since, as demonstrated in
Figure 2, the divisions of
directly correspond to the order by which jobs visit workstations, which results in one and only one
, and the divisions of
, as defined according to the set of positions
, directly represent the processing sequences of operations on machines, which results in one and only one
. Therefore, any change in
,
or
will result in a solution pair,
≠
. Let schedules
and
be generated via Algorithm 1 from solution pairs
and
, respectively, where
≠
. Since Algorithm 1 is a semi-active schedule generator that schedules operations in their orders provided in a solution pair using the earliest start time for each operation, as stated in step 16, and by definition, any operation in a semi-active schedule cannot start earlier without changing the sequence to process the operations [
46],
. Therefore, a schedule generated via Algorithm 1 cannot have more than one input VPB chromosome. □
The property stated in Proposition 1 represents an advantage of the proposed VPB chromosome structure compared to the OB one in terms of eliminating redundancy from the chromosome domain to the solution space. Theoretically, this can help enhance the search efficiency by avoiding search moves that change the chromosome but do not change the solution. However, some guidance is needed for the search moves of the VPB representation to avoid generating chromosomes that result in infeasible solutions.
To construct solutions to be used within the proposed MOSS, the construction heuristic developed in [
4] is used as is. That heuristic is based on applying different rules for selecting machines and selecting jobs to be scheduled at each iteration. Such scheduling rules tend to either improve the makespan or improve the MWFT. In addition, a random selection rule is used to allow some variability in the generated initial population. The encoding of a generated solution to its corresponding VPB chromosome is straightforward since the processing orders of operations of each job and on each machine and the jobs’ visiting order to their required workstations are clearly defined in every schedule.
3.3. Measuring Dissimilarity between Two Solutions
An essential characteristic of the scatter search metaheuristic is maintaining diversity among the solutions in the reference set
along with their quality throughout both the initialization and the improvement phases. Accordingly, before a solution to set
is added, its dissimilarity with existing solutions in this set is checked. This is done here by measuring the degree of dissimilarity, which is defined for any two solutions
and
using a distance function denoted as
, which is provided in Abdelmaguid [
29]. A solution can be added only when its evaluated distances with all existing solutions in set
exceed a predefined limit.
The distance function
is based on comparing the vectors of permutations of the two solutions,
and
, and counting the minimum number of moves that are needed to convert one solution to another. To illustrate how the distance function
is evaluated, consider the two sample solutions shown in
Figure 4 for the Bi-DMOSP instance described in
Section 3.2. Determining the minimum number of moves required to convert the two strings of the second solution’s chromosome involves the determination of the common sub-strings in both chromosomes and then choosing the moves that will result in the minimal number of common sub-strings in the resulting strings. Details are provided in [
29] for a computationally efficient algorithm that can be used for that task.
Figure 5 illustrates the minimal moves that are required to convert the chromosome of the second solution into the first one’s. The first move (labeled with the circled number 1 in
Figure 5) corresponds to exchanging the positions of operations D1 and D3 in the jobs’ string. Meanwhile, the move labeled with the number 3 corresponds to moving operation F1 after the position of the machine separator while shifting operation E1 backward in the workstations’ string. Such a move changes the sequence of processing for the operations of the machines of the workstation, and it changes the machine assignment for operation F1. Other moves are conducted similarly. As a result of the illustrated minimal moves, the distance between both solutions equals 5 in this case.
This distance function’s advantage is that it does not depend on the objective functions’ values, enabling the search procedures to explore diverse areas in the search space. To explain, consider the two sample solutions presented in
Figure 4. The objective functions’ values for the first solution,
, are
and
, and for the second solution,
, they are
and
. Accordingly, the differences between the objective functions’ values for both solutions are not a small quantity compared to the value of
. Therefore, if the distance based on the differences between the objective functions’ values is utilized, both solutions would have been added to the reference set. Nevertheless, in this case, one of the solutions would easily reach the same structure as the other one in a few neighborhood search moves. Therefore, using the value of
will ensure true diversity as related to the neighborhood search moves that are utilized in the local search procedure based on tabu search.
3.4. Solution Recombination
Solution recombination is an essential operator in the search procedures of modern metaheuristics. It combines some features from parent solutions to generate a new child solution. This operator helps explore more regions in the search space between or even beyond the parent solutions. This operator appears in evolutionary metaheuristics under the name of a crossover operator [
47]. In swarm-intelligence metaheuristics, the mechanism of solution recombination is conducted via position-update functions [
48]. Path relinking is another form of solution-recombination operators [
39].
In the current implementation of the scatter-search metaheuristic, the solution-recombination operator works as follows. First, two different parent solutions are selected randomly from the reference set, where all solutions have equal probabilities of being selected. One selected parent is labeled the leader (), and the other is labeled the follower (). These labels are used to identify the order by which genes will be inherited through the child solution ().
Figure 6 demonstrates the gene-inheriting mechanism using solutions for the sample instance presented in
Section 3.2. In this mechanism, genes are inherited for each sub-string independently for both strings constituting each chromosome. For each sub-string, a number,
, is randomly generated from a uniform distribution. If
is greater than a predefined recombination threshold value, denoted as
, the genes from the leader parent’s sub-string will be copied as-is to the child solution.
When
, genes in the child solution’s sub-string are inherited from both parents’ sub-strings as follows. First, an integer number,
, is randomly generated from a discrete uniform distribution, where
ℓ is the length (number of genes) of the sub-string.
Figure 6 represents the randomly generated positions using a small arrow underneath the sub-strings. Then, genes from the leader solution’s sub-string are copied to the child solution’s sub-string, starting from the first gene up to the gene in location
. Afterward, starting from location
, the remaining genes are inherited from the follower parent using the same sequence they appear in that parent. Positions of the machine separators are inherited as-is from the parent having the larger number of inherited genes in the corresponding workstation’s sub-string, where tie-breaking is done randomly.
Since the resultant child solution can be infeasible, the recombination operator is applied a maximum of five times on the same selected parent solutions until a feasible child solution is generated. If no feasible solution is obtained after these five iterations, either the leader or the follower solution is copied to the child solution. In this case, both parents have equal copying probabilities.
3.5. A Novel Bi-Objective Tabu Search
Local improvements to solutions generated during the initialization and the improvement phases are conducted using a proposed novel bi-objective tabu search (Bi-TS) approach. This Bi-TS utilizes the two neighborhood search functions developed by Abdelmaguid [
29] to improve the makespan, denoted as
and
. To improve the
, two similar neighborhood search functions, denoted as
and
, are introduced here.
The neighborhood search function
implements moves of operations within
permutations. It performs moves on critical operations that, by definition, cannot be moved from their current positions in the permutations without affecting the current makespan [
46]. Accordingly, a critical path is defined as a sequence of consecutive critical operations from start to finish. First,
selects an operation belonging to the critical path, and then it investigates possible moves within the processing sequence of the job to which the selected operation belongs. The feasibility of a move and the quality of the resultant solution are derived by employing simple mathematical relationships. These relationships reduce the required computational time by avoiding the detailed evaluations of the operations’ exact start and finish times. Similarly, the neighborhood search function
implements moves of critical operations within and between
permutations for the machines belonging to the same workstation.
The new functions
and
follow the same structure of operation removal-reinsertion mechanism, as in
and
, respectively. This similarity allows for a unified tabu structure and tabu list for all neighborhood search moves. In
and
, the feasibility conditions proven by Abdelmaguid [
29] for moves conducted in
and
are employed. Unlike
and
, they are applied to critical and non-critical operations. The
objective values as a result of applying
and
are calculated by modifying the solution’s network representation, which was provided in Abdelmaguid [
29]. Since this is a time-consuming process, non-critical operations that belong to the first half of a job or machine sequence are excluded from the search. This exclusion is decided after it is noticed that the value of
generally does not improve by moving operations scheduled early in a given solution.
The search structure of the developed Bi-TS approach is demonstrated via the pseudo-code in Algorithm 2. It starts in steps 4 to 7 by initializing the search and iteration variables.
Algorithm 2 Bi-objective Tabu search for the Bi-DMOSP |
- 1:
function Bi-TS() - 2:
Input: MOSP solution: - 3:
Output: Improved solution: - 4:
Start with empty tabu list - 5:
= , = , = - 6:
1 - 7:
0 - 8:
while do - 9:
- 10:
If then BestMoves = ordered list of at most non-tabu neighborhood moves from all possible applications of and on resulting in—and ordered according to—the lowest estimates - 11:
Else BestMoves = ordered list of at most non-tabu neighborhood moves from all possible applications of and on resulting in—and ordered according to—the lowest estimates - 12:
If then and - 13:
If then NSMove = randomly selected move from BestMoves - 14:
and - 15:
Else NSMove = the first (best) move from BestMoves - 16:
= resultant solution of applying NSMove to - 17:
NSMov = inverse of NSMove that is needed to transform back to - 18:
If NSMov is not tabu, then update the tabu list by adding NSMov - 19:
Else update the tabu list by adding NSMove - 20:
- 21:
- 22:
If ( and ) or ( and ) then = , = , = , Check and update set based on - 23:
Else If ( and ) or ( and ) Check and update set based on If is added to set then = , = , = , End If - 24:
End If - 25:
If then = + 1 - 26:
If then = + 1 - 27:
= , , - 28:
end while - 29:
return - 30:
end function
|
Within the main loop, the algorithm keeps a record of the number of times each objective value is improved as a result of applying the neighborhood search functions. This is done by updating two variables, and , for recording improvements in the makespan and the MWFT objective values, respectively. These two variables are used to probabilistically determine which set of neighborhood search functions will be applied in each iteration: either or . This is done in step 10 by calling the function, which is a Boolean function that returns if a randomly generated number within the range from 0 to 1 is found to be less than the passed value p and returns otherwise. Accordingly, the ordered list is used to store the best non-tabu neighborhood moves based on the current solution . The maximum number of neighborhood moves that can be stored in the list is denoted , which is one of the control parameters of the developed Bi-TS.
From the list, either the first (best) neighborhood move or a random one is selected. Random move selection is permitted to allow for a more diversified search whenever the repeated selection of the best moves does not improve the dominance of the incumbent solution for iterations. If a random selection of moves is decided, it is repeated for at least consecutive iterations of the Bi-TS until an improvement occurs. Accordingly, this diversification mechanism is controlled via the two parameters and .
The tabu list is updated in step 18 by adding the inverse of the selected move, which is the move that will bring the solution back to its original structure before applying the selected move. If the inverse move is already included in the tabu list, the best move itself is added instead, as shown in step 19. The tabu list has a fixed size, denoted as , which cannot be exceeded. If the tabu list has moves included, the oldest move is removed.
The dominance of the incumbent solution is checked starting at step 20, when the function is applied to the difference between the objective values of the best found and the incumbent solutions. The function returns −1, 0, or +1, depending on the sign of the passed argument. If the incumbent solution is found to dominate the best found solution, a replacement is made, as shown in step 22. In such a case, the global non-dominated set of solutions, , is checked for an update. Otherwise, in step 23, if the incumbent and best found solutions do not dominate each other, set is checked for an update via an investigation of the dominance of the incumbent solution with all the solutions therein. If the incumbent solution is added to set , it will replace the current best found solution.
3.6. Updating and Diversifying Solutions in the Reference Set
During the improvement phase of the developed MOSS, the set is used to update the reference set such that does not exceed the specified fixed value of the parameter . This is done by first removing all solutions from set and adding them to set . Then, an ordered list, , is generated based on solutions in in which solutions are ordered based on their dominance relations. That is, if dominates , will be ordered first, and vice versa. If the two solutions do not dominate each other, their order will be arbitrary.
In the given solutions’ order in
, solutions are iteratively added to set
until its size reaches
. To maintain the diversity of the solutions in set
, before a solution is added, its distance from every already added solution is evaluated as described in
Section 3.3. This distance must exceed the predefined distance threshold value (
). Otherwise, the solution will not be added. After all solutions in
are examined, if the size of the resultant set
is found to be less than
, additional solutions will be generated and added to it using the same procedure described in the construction process, as in
Section 3.2.
The developed reference set-updating procedure aims to maintain both solution quality and diversity throughout iterations. Furthermore, to avoid long-term entrapment in local optimal solutions, a diversification strategy is annexed to the reference set-updating process. In this strategy, the percentage changes in the summations of the values of the makespan and the MWFT for all solutions in set
are evaluated before and after every improvement cycle. If these percentage changes are found to be less than 2%, and the set of non-dominated solutions (set
) has not been changed by adding new solutions during the improvement cycle, the solutions in the reference set
will be diversified. The diversification is achieved by replacing one-quarter of the solutions with newly constructed ones using the construction procedure described in
Section 3.2. New randomly constructed solutions are improved using tabu search, as usual. In addition, at most one-quarter of the solutions in the reference set
are replaced with solutions selected arbitrarily from the non-dominated set
.