An improved ant colony optimization for solving the flexible job shop scheduling problem with multiple time constraints

1.   Introduction

The flexible job shop scheduling problem was originally developed by Brucker and Schlie ^[1], and since then the problem has received a lot of attention. In traditional shop scheduling, the machines for all operations of a job are unique, and only the sequence of operations on each machine needs to be determined. And flexible job shop scheduling eliminates the limitation of resource uniqueness ^[2], and each operation of a job can be processed on multiple machines. Scheduling requires not only sequencing the operation but also making decisions about the machines that will process the operation.

FJSP started with the main research objective of optimizing the maximum completion time. Later on, as the research continued, more objectives came into focus. There are optimization objectives regarding machine load that was studied in earlier years ^[3,4]. Different jobs in the flexible job shop are processed on the same machine requiring job setup that consumes a lot of time. Many scholars have considered the loading of workers and the change of tools, as well as the setup time that exists for machining on the same machine ^[5,6,7,8]. In addition, in flexible job shop scheduling, the same job often needs to be processed on different machines. This will greatly reduce the feasibility of the scheduling solution if the actual execution does not deliver the job or material to the corresponding station on time. Some scholars have added constraints on transport routes and transportation time considering the transport processes that exist in the middle of job processing in different machines ^[9,10,11,12]. In recent years, some scholars have also considered both transportation time and setup time ^[7,13]. The demand for punctuality in manufacturing orders is very high, and companies need to find a suitable production completion point according to the delivery date. A heterogeneous DAFJSP with the total cost and the earliness/tardiness as the objectives is studied ^[14]. The above literature has studied FJSP under different objectives, however, these papers are relatively scarce for considering multiple time constraints of setup time, transportation time, and delivery time simultaneously. ^[2] studied the flexible job shop scheduling problem with multiple time constraints. The paper treats transportation time and setup time as independent time factors, with the objectives of minimizing maximum completion time, minimizing total setup time, and minimizing total transportation time. It does not take into account the load capacity of individual machine and the delivery time of orders, which are important references for actual production. In ^[13], transportation and setup time are considered, while minimizing the maximum completion time, total delay, and total energy consumption. However, it ignores the limitations of the workload of each machine in a realistic shop.

As the demand for customized manufacturing continues to grow, multi-variety and low volume have become a very important aspect of modern manufacturing companies ^[15]. The problem studied in this paper is abstracted from a customized furniture production shop for panels. In this customized production shop, the diversity of products produced increases the setup time and transportation time. In this case, the scheduling shop that considers only the processing time schedule will have a large error with the real production completion time. Another requirement that cannot be ignored in custom processing orders is the delivery time. Different orders have different delivery time. Producing the product too early takes up storage space, while producing it too late violates the order requirements, so the right production completion time is important for customized production. Therefore, based on this customized production shop this paper considers processing time, setup time, transportation time, and order delivery time as independent time factors for FJSP.

The study of the multi-objective flexible job shop scheduling problem (MOFJSP) is of great engineering significance, but the difficulty of this type of problem lies in dealing with the existence of multi-objective conflicts. The individual solutions found in this problem are not directly comparable and there are non-dominated solutions. For any two solutions S₁, and S₂. if all solutions of S₁ are better than S₂, it means S₁ dominates S₂. If there exists a solution within S₁ that is not dominated by other solutions, it means that S₁ is a non-dominated solution. The set of these non-dominated solutions is the Pareto Front. Zhang et al. ^[4] use a weighted summation method to construct the fitness function in solving MOFJSP, but the weights of each objective are difficult to determine. Non-dominated sorting is one of the more widely used sorting methods in recent years ^[16]. Therefore, this paper addresses these problems by proposing a combination of the entropy weight method and non-dominated sorting to solve the multi-objective problem.

It is well known that the traditional FJSP has proven to be an NP-hard problem. The FJSP, which considers multiple time constraints, is more complex and therefore requires an efficient algorithm to solve the problem. At present, many scholars in this field have used different optimization algorithms to solve FJSP. For example, genetic algorithms (GA) ^[17,18,19], ant colony optimization (ACO) ^[20], particle swarm optimization (PSO) ^[21], simulated annealing (SA) ^[22], imperialist competitive algorithm (ICA) ^[23], African buffalo optimization (ABO) ^[24,25], whale optimization algorithm (WOA) ^[26,27]. Ant colony optimization is a swarm intelligence algorithm. It is mainly used to solve traveling salesman problem (TSP) problem, assignment scheduling problem, and job-shop scheduling problem, and achieve better experimental results. Its search process uses distributed computing, with multiple individuals performing parallel computation at the same time, greatly improving the computational power and operational efficiency of the algorithm. In this paper, we choose ant colony optimization to solve the problem by combining the problem characteristics and algorithm features. Ant colony optimization has the characteristic of positive feedback. If the algorithm starts with a suboptimal solution, the positive feedback will make the suboptimal solution dominate quickly and make the algorithm fall into the local optimum, and it is difficult to jump out of the local optimum. Therefore, the use of ant colony optimization is more demanding for the initial solution.

Based on the characteristics of the ant colony optimization, this paper designs a distributed coding method and three ways to generate the initial solution to improve the quality of the initial solution. The search is performed using an improved ant colony optimization for operation sequencing. The scheduling set is then divided into two populations by combining the entropy weight method and the non-dominated sorting. Finally, mutation and closeness operations are proposed for the machine assignment of the population to improve the quality and diversity of the population, respectively.

In summary, there are relatively few studies that consider multiple time constraints on setup time, transportation time, and delivery time simultaneously. Therefore, this paper proposes an improved multi-objective ant colony optimization to solve the flexible job shop scheduling problem with multiple time constraints. The main contributions of this paper are as follows.

1) A flexible job shop that simultaneously considers multiple time constraints such as transportation time, setup time and delivery time, and takes into account machine load capacity, is studied from the perspective of actual production needs. Better guidance is provided for actual production.

2) To solve this problem effectively, an improved ant colony algorithm is proposed. The distributed coding approach designed in this paper expands the search range of the solution. Then the proposed three initialization strategies can generate high-quality initial solutions. In addition, the improved ant colony algorithm enhances the search accuracy of the algorithm. Finally, the population diversity is further improved by non-dominated selection combined with mutation and closeness operations.

3) The feasibility of the proposed algorithm in solving flexible job shop scheduling under multiple time constraints is verified through extensive experiments.

The rest of this paper is as follows. Section 2 describes the flexible job shop scheduling problem with multiple time constraints in detail and provide an instance of a small-scale for understanding. Section 3 details the overall process of the algorithm in this paper. In Section 4, we compare the algorithm proposed in this paper with other algorithms on 28 benchmark instances. Section 5 presents our conclusions and plans for future research.

2.   Problem formulation

MaOFJSP is described as follows. There are n operations and m machines in a shop. Each job consists of h_i operations that must be processed in a specified sequence. Each operation can only be processed by one machine, which is any of the optional sets of machines. The processing time for each operation depends on which machine is processing the operation. When an operation of a job completes, it is moved to the next operation. There may be transportation time and setup time between the completion time of one operation and the start time of another consecutive operation. If two adjacent operations of the same job are processed on different machines, there is a transfer process, which generates additional transportation time. Setup time exists on the machine when two consecutive operations on the same machine are different types of jobs. Note that each machine needs to be set up for a time before the first job can be processed. Transportation time and setup time are zero when two adjacent operations of the same job are processed on the same machine. The problem is to assign operations to machines and find the sequence of operations on all machines to create a feasible plan that minimizes makespan, total workload (TW), workload of critical machine (WCM), and penalties of earliness/tardiness (PET). The following assumptions exist in this paper for MaOFJSP.

1) Each machine can only process one operation at a time, and once the machine starts processing this operation it cannot be stopped.

2) At the moment 0 all machines are available for work and at the moment 0 all jobs are released at the same time.

3) The first job transportation time of the operation is 0, but its setup time is not 0.

4) Each operation requires only one machine for processing.

5) Each job has a fixed sequence of operations to be processed, and after each operation is completed the next operation will be processed.

6) When two adjacent operations of a job are processed on different machines, the transportation time is determined by the two machines. The transportation time between the machines is known.

7) The setup time before processing is different for different operations on different machines. The setup time before processing on different machines is also known for each operation.

For convenience of description, we define the notations in Table 1.

The four optimization objectives of this paper are as follows:

1) Minimization of makespan (f₁):

2) Minimization of TW (f₂):

3) Minimization of WCM (f₃):

4) Minimization of PET (f₄):

To facilitate understanding, we now provide an instance of a small-scale flexible job shop scheduling problem. Table 2 shows an instance of a flexible job shop scheduling problem. The instance consists of 3 jobs and 4 machines. Table 2 indicates the processing time and setup time required to handle different operations for different jobs on each machine. Table 3 indicates the transportation time between the different machines.

Figure 1(a) shows the scheduling Gantt chart without setup time and transportation time. In Figure 1(a) 101 is the first operation of the job. Operation 101 corresponds to 4 on the y-axis and 1 on the x-axis, indicating that the time required for operation 101 to be processed on machine 4 is 1-time unit. Figure 1(b) shows the scheduling Gantt chart considering setup time and transportation time. The brown box in Figure 1(b) represents the setup time and the light blue box represents the transportation time. The light blue 301 box indicates that the transportation time required for job 3 to be transported from machine 4 to machine 2 is the 1-time unit. As can be seen in Figure 1(b), the scheduling Gantt with setup time and transportation time has a maximum completion time of 5 more than the scheduling Gantt without setup time and transportation time. In the actual production process, the extra 5 make sense. The feasibility of the resulting plan is low if setup time and transportation time are not taken into account.

3.   The proposed IACO for MaOFJSP

In this paper, we propose an improved ant colony optimization to solve MaOFJSP. The framework of IACO is shown in Figure 2. The machine assignment is first determined and the solution space is decomposed into multiple mutually independent regions. Then an improved ant colony optimization is used to search for the optimal scheduling solution for each optimization objective in this region. In addition, a method combining entropy power method and non-dominated ranking method is presented for finding non-dominated solutions. Finally, mutation and closeness operations are proposed for the machine assignment part to improve the search accuracy. The algorithm is described in detail in Sections 3.1–3.6.

3.1.   Chromosome coding

L Coding is the representation of a feasible solution in simple chromosomal form, converting the feasible solution into an array form that can be processed by a computer. An effective coding approach can better express the relationship between individuals and feasible solutions. FJSP contains two subproblems: machine assignment and operation sequencing. The two-stage integer coding method is one of the most common coding methods for FJSP. A sequence of two ends of length $L = \sum O_{i j}$ is used to represent machine assignment and operation sequencing, respectively. However, this coding is not sufficient for the search of individual solution spaces and can lead to an inefficient start of the search. In this paper, we propose a distributed coding approach, which is different from the traditional two-stage integer coding of the search framework. This coding approach first determines the machine assignment and decomposes the solution space into multiple mutually independent regions. This is then combined with an ant colony optimization to sort the operations in order to adequately search the solution space allocated to this machine. The final set of non-dominated solutions is obtained by continuously adjusting the machine assignment and operation sequencing.

1) Machine assignment: For machine selection, the length L of the chromosome is the total number of operations. One machine is assigned to each operation. The gene strings of the machine part are represented as integers. Each integer is generated by the three methods of generating initial solutions proposed in this paper, and no illegal solutions are generated. As shown in Figure 3, the gene string is 1-3-2-3-2-2-1. The 1 in the first light green box means that operation O₁₁ selects the first machine to be processed from among the three machines that can be selected. Operation O₁₃ selects the second machine for processing, and so on.

2) Operation sequencing: When sequencing the operations, the length L of the gene string is the total number of operations. The number of integer occurrences represents the total number of operations for this artifact. The order of occurrence of an integer is the order of operation. As shown in Figure 3, the gene strings are 1-1-2-2-1-3-3 and their processing order is O₁₁-O₁₂-O₂₁-O₂₂-O₁₃-O₃₁-O₃₂.

3.2.   Chromosome decoding

Decoding is the process of converting chromosomal sequences into a number of scheduleable schedules. An effective decoding approach can lead to better results. For both sequences of machine assignment and operation sequencing, we first decode the machine part and then decode the operation sequence using full insertion decoding. For the machine assignment part, it is first converted into a machine selection matrix $J_m(i, j)$, a machine processing time matrix $T(i, j)$, and a machine setup time $H M(i, j)$. We use the four operations O₁₁, O₁₂, O₂₁, O₂₂ of Tables 2 and 3 to decode the following equations.

$J_m(i, j)$ represents the set of optional machines of $O_{i j}$. $T(i, j)$ denotes the processing time of $O_{i j}$ by different optional machines. $H M(i, j)$ represents the setup time for processing $O_{i j}$ by different machines. These three time matrices are interconnected. After decoding these three time matrices, the machine selects the chromosome part. The machine for the corresponding operation is selected from the $J_m(i, j)$ matrix based on the individual integer values of the chromosome sequence from left to right.

The operation part uses full insertion decoding to reduce the maximum completion time as much as possible. The principle of this method is to insert operations in scheduled machine idle time slots. The insertion operation will advance the start time of the operation, thus allowing the total completion time to be reduced. The insertion operation is performed as follows: When h operations are processed on a certain machine M, (h-1) idle time slots will be generated. When $O_{i j}$ is not the first operation of the job and $O_{i j}$ has already been processed by machine M. Find the processing machine corresponding to $O_{i j}$ and find whether the free processing segment on that machine satisfies the early insertion of $O_{i j}$. The specific steps are as follows.

1) If the previous operation of this machine idle time period is the same job as $O_{i j}$, the setup time and transportation time will be 0.

2) Calculate CT (CJ_Oi(j-1) + transport time) with ET (completion time for the operation before the gap T₁ + setup time). Select Max (CT, ET) and $O_{i j}$ processing time on this machine and add to get Q. If Q is less than T₂ (the start time of the operation after the gap) then the insertion is satisfied, as in type (a) in Figure 4, and $O_{i j}$ is inserted into the gap. If Q is greater than T₂ then the insertion is not satisfied, as in type (b) in Figure 4, and the search for the next gap continues. If all gaps are not satisfied, $O_{i j}$ is arranged in the original order. For each completed step, the start time and completion time of the operation will be calculated and recorded. Then the value of each target is calculated according to Eqs (1)–(4).

3.3.   Initial population

The initialization of the population is a very important step in the algorithm. A good initial solution has an important impact on the speed and quality of the algorithm. The distributed encoding approach proposed in this paper requires more quality and diversity in the initial solution assignment by machines. Therefore, in this paper, three initialization methods are proposed to improve the quality and diversity of the initial solutions in response to this requirement.

1) Random selection method: We randomly generate a machine selection sequence. In the machine sequence, we randomly select a machine from among the machines that can be selected by the operation.

2) Least time selection method: Each operation may have more than one machine available for processing. From these, we select the machine that takes the least amount of time to process the operation.

3) Earliest start time selection method: The machine with the smallest start time is searched for in the operation of the optional machine.

Because of the need to balance the quality and diversity of the initial solutions, we search for these three different initialization methods using different probabilities. A random selection is made in 50% of a population, 30% choose the machine with the least processing operation time and 20% choose the machine with the least operation can start time.

3.4.   Ant colony optimization based local search

Ant colony optimization is a probabilistic algorithm used to find the optimal path in the graph. According to the characteristics of FJSP, an improved ACO is proposed in this paper as the local search and calculation of the objective value.

In this part of the algorithm, the neighborhoods corresponding to the individuals in the machine assignment population are searched. For β objectives, β ant colonies are generated to search for a local-optimal solution on each objective. The specific steps are as follows:

1) Initialize the pheromone concentration table, and the concentration is all PC_init. The pheromone concentration table is shown in Figure 5, the number of rows is equal to the number of jobs and the number of columns is equal to the number of operations. The value in each grid indicates the pheromone of the selected job when machining the corresponding operation. For example, the value in grid (i, j) is to select job i as the jth processing operation.

2) Determine the operation sequencing of each ant according to the roulette method. When determining the jth processing operation, the job set is first to be processed, then the fitness value $P_{i j}^k$ of each job is calculated by Eq (8) and composed a turntable. Each job is distributed on the turntable according to the $P_{i j}^k$. Each time the turntable rotates, the point indicated by the arrow is the selected job.

where i is the last determined job, and j is the job in the job set allowed^k for jth processing, $\tau_{i j}(k)$ is the jth pheromone concentration from job i to j, $\eta_{i s}(k)$ is the reciprocal of the distance $d_{i j}^k$ from job i to j, as shown in Eq (9):

where S_J is the allowable start time of the job, S_M is the allowable start time of the machine. T₂ is the transportation time of the job and T₃ is the setup time of the job. In fact, $d_{i j}^k$ reflects the size of the machine idle gap that will be generated by job j at time k. If the gap is 0, it means that scheduling machining job j at this time can improve local productivity. Since 0 cannot be the denominator, $d_{i j}^k$ is taken to be 0.1 at minimum.

Since each job has multiple operations, each job can be selected multiple times by ants, and when the maximum number of selections is reached (the number of operations for the job), the job is removed from the allowed^k set. This ensures that when the allowed^k set is empty, the processing order of all operations is determined, resulting in a complete scheduling solution. While determining the processing sequence of an operation, it is also necessary to calculate information such as the start time and end time of processing the operation, and with this information, the scheduling plan can guide the actual production. For each machine and each job, set its initial allowable start time S_Mk and S_Ji to 0. Whenever an operation is determined, the start time ts and the end time te are calculated for that operation, then the allowable start time S_M of the corresponding machine and the allowable start time S_J of the workpiece are updated. The equations for ts and te are shown below.

3) When the ant selects all the jobs, the objective value of the corresponding scheduling scheme of the ant is calculated according to the Eqs (1)–(4). The time information of each operation has been calculated at the same time in the selection process, and there is no need to repeat the encoding and decoding.

4) Record the scheduling scheme corresponding to the optimal individual in the ant colony.

5) Update pheromone. In this paper, the pheromone update of ant colony optimization with elite strategy is used to update the pheromone secreted during the traversal of ants. Each completed cycle gives an additional pheromone increment to the path searched by the elite ant. This approach accelerates the increase of pheromones on the optimal path, which leads the subsequent ants to find the global optimal path and the global optimal solution quickly. The pheromone update with elite strategy is as follows.

$\Delta \tau_{i j}(k)$ denotes the sum of the pheromone increments left by route (i, j) during the traversal of the ant. $\Delta \tau_{i j}(k)$ denotes the sum of the pheromone increments left by the elite ant pathway (i, j). $\rho$ is the pheromone evaporation rate. $\omega^{\varphi}$ is the number of elite ants.

Based on the proposed paper when giving additional pheromone increments to elite ants, it may lead to the pheromone stacking on a certain path and falling into a locally optimal solution. We control the amount of pheromones on each path by setting a maximum value and a minimum value for the pheromone values on the ant search path. Limit the amount of pheromones to [$\tau_{\min }, \tau_{\max }$]. When the pheromone value of (i, j) is greater than $\tau_{\max }$, $\tau_{i j}(k) = \tau_{\max }$. When the pheromone value of (i, j) is less than $\tau_{\max }, \tau_{i j}(k) = \tau_{\min }$. This results in better control of the amount of pheromones on each path.

6) Judge whether the ant colony optimization is terminated. If the algorithm is terminated, the current optimal solution will be output. If not, return to step 2) continue to execute.

3.5.   Non-dominated sorting selection

Due to the existence of four objectives in this paper, the obtained objective values cannot be compared with each other. And each objective is not able to obtain an optimal solution by constraining each other. Thus the scheduling set obtained by the IACO algorithm above is a non-dominated sorting selection ^[16]. First, the four target values are assigned weights by the entropy weighting method, and the comprehensive score of this scheduling scheme is calculated. After the combined score of each solution is obtained, the solutions are then ranked non-dominated and the solution set is divided into different dominated layers. The solutions of different dominating layers are selected in turn up to the F_n layer. When the number of selected solutions is greater than or equal to N_pop, the population consisting of the first N_pop solutions is treated as F₁, and the population consisting of the remaining solutions is F_2n. The machine assignment information of the solutions in F₁ and F_2n is extracted to form population F₁ and population F₂. This process is shown in Figure 6(a).

The basic concept of Pareto theory is domination. Domination determines the relationship between two solutions. A condition for a solution X₁ to completely dominate another solution X₂ is that all objectives of X₁ f_m(X₁) (m = 1, 2, …, N) are all better than X₂. If both solutions directly have a better objective than the other one, the two do not dominate each other. The process of delineating the different non-dominated layers is shown in Figure 6(b). The three solutions in Rank1, i, i-1, and i+1 are not dominated by each other.

The basic idea of the entropy weighting method is to determine the objective weights based on the magnitude of the variability of the indicators. According to the definition of information entropy, for a certain index, the entropy value can be used to judge the dispersion degree of a certain index. The smaller its information entropy value is, the greater the dispersion degree of the indicator, and the greater the influence of the indicator on the comprehensive evaluation. The specific steps of the entropy method are as follows.

Normalization of the target values using Eq (16).

where x_ij is the original value of the jth indicator data for the ith object, y_ij is the jth indicator value for the ith object, i = 1, 2, ..., n, n = 4 × N_pop.

Equation (17) is used to calculate the share of the ith object in this indicator under the jth indicator Y_ij:

Calculate the entropy value e_j and the information utility value d_j for the jth indicator with the formulas shown in Eqs (18) and (19).

The weights W_j of the values of j indicators can be obtained as shown in Eq (20).

Calculate the composite score F, as shown in Eq (21).

3.6.   Mutation and closeness

Mutation operator: we design the mutation operator for the individuals in F₁ to increase the search scope. As shown in Figure 7, let the mutation position is 1, 3, 5 and 6, and the values are randomly selected and transformed into any number of optional machines at this position. After mutation, the value in position 1 is transformed into 3 and value in position 3 is transformed into 2, etc., to generate a new individual in the machine assignment population.

Closeness operator: The closeness operation is performed on the individuals in F_2n to approach the individuals in F₁. First, calculating the Hamming distance between the individuals in F_2n and F₁. The Hamming distance is the number of different corresponding positions in the two individuals. As shown in Figure 8, the information of position 1, 2 and 4 of individual 1 and individual 2 is different, so the Hamming distance between them is equal to 3. Then, each individual in F_2n randomly approaches to one nearest individual in F₁, approaches the machine part in the way of multi-point cross mutation, randomly generates multiple positions of cross, replaces the individual information in F₁ to the corresponding individual in F_2n to be approached, and the information in other positions of individuals in F₁ remains unchanged.

The next generation with the information of machine assignment is composed by merging the two population after mutation and closeness.

4.   Experimental results

4.1.   Test instances

We use Matlab 2016a to perform simulation experiments for this algorithm, running on an Intel(R) Core(TM) i5-7300HQ CPU @ 2.50 GHz and 2.50 GHz. In this paper, we use the Brandimarte dataset ^[28] containing 10 instances (MK01-MK10) and the Dauzere-Peres and Paulli dataset ^[29] containing 18 instances (01a–18a) for our comparative experiments. Where the data in the MK sets are small-medium scale problems, while the data in (01a–18a) are large-scale problems. The transportation and setup times in the instances are randomly generated. The equations for generating the start delivery date and end delivery date are

where $\theta_i$ is the random number. The range of values in Mk01-Mk04 is [1, 1.2], the range of values in Mk06 is [0.8, 1], the value range in MK10 is [1.6, 1.8], and the value range in MK10 is [1.6, 1.8]. γ is the delivery interval, Mk1-Mk10 is 10, 01a-18a is 100. All detailed data for this paper can be found in the following websites: https://pan.baidu.com/s/16_zjfTGb6A2JMEsLVRELrA code: ysf1.

4.2.   Performance metrics

To evaluate the IACO proposed in this paper, we use the set coverage (SC) ^[30], Hypervolume (HV) ^[30], the metrics of the number of non-dominated solutions (NONS), and the convergence of the algorithm are used to evaluate the algorithm proposed in this paper.

The SC is described as follows:

A and B are the two solution sets that we are to compare. SC(A, B) represents the proportion of solutions in B that are dominated by solutions in A. When SC(A, B) = 1 means that the solutions in B are all dominated by A. When SC(A, B) = 0 means that none of the solutions in B is dominated by the solutions of A. It is worth noting that SC(A, B) is not necessarily equal to 1- SC(B, A). If SC(A, B) > SC(B, A) then it means that the number of solutions in A dominates the number of solutions in B, and the solution set A is better than the solution set B.

HV is used to represent the volume of the hypercube enclosed by the individuals in the solution set and the reference points in the target space. The convergence and distribution of the solution set S are evaluated by calculating the HV value, which is defined in Eq (25).

S is the approximate solution set of the Pareto optimal frontier. P is the reference point corresponding to the Pareto frontier. v(s, P) denotes the hypervolume of the space formed between the solution s and the reference point P in the non-prevailing solution set S. Larger HV represents better diversity of the solution set S.

4.3.   Parameter settings

The algorithm proposed in this paper has three key parameters: evaporation rate (E), mutation probability (M), and pheromone factor weight (P). In order to choose the best combination of parameters for the algorithm, the currently widely used design of experiment (DOE) method of Taguchi ^[31] is used on a medium-sized problem (MK01). The parameter level table as in Table 4 has three parameters, each with three levels, and a total of nine different combinations of the three parameters can be seen in the orthogonal table of parameter combinations in Table 5. The HV (Hypervolume) index is used to evaluate different combinations of parameters, and the accuracy of HV is highly dependent on the choice of reference points. The reference set for this experiment is the optimal set filtered by the algorithm under MK01 instance using the first three coefficient combinations after non-dominated sorting. Figure 9 shows the parameter level trends chart. Based on the results of DOE, the optimal values of the parameters are set to E = 0.2, M = 0.4, and P = 2.The other parameter values of the algorithm were set as in Table 6.

4.4.   Comparison with other classical algorithms

In order to evaluate the effectiveness of the proposed algorithm in this paper more comprehensively, we compare IACO with NSGA-II ^[32] and IGA ^[15]. The reasons for choosing these two algorithms as comparison algorithms are: (ⅰ) These algorithms are multi-objective optimization algorithms, which are easy to reproduce. (ⅱ) Both of these algorithms are very classical and have been proven to be effective. Considering the randomness of the algorithm. For each instance, each algorithm was run 10 times independently with a stopping condition of 200 iterations. We take the average of 10 times results for comparison and analysis.

The sc-based results obtained for IACO, IGA, and NSGA-II are given in Table 7. As can be seen from Table 7, IACO obtains a larger number of non-dominated solutions in the instance than the other algorithms. The bold black font in the table represents the better data between the two. In almost all instances, SC(HA, A1) > SC(A1, HA), SC(HA, A2) > SC(A2, HA). For example, in Mk01, Mk03, Mk04, MK05, Mk06, Mk07, Mk09, Mk10, 01a, 02a, 03a, 04a, 05a, 06a, 07a, 08a, 12a, 14a, 15a, 18a. According to the box plot in Figure 10, these relationships are even more evident. The red line in the graph is the mean value line and the small blue squares are the outliers.

Table 8 shows the non-based results from IACO, IGA, and NSGA-II. In Table 8, IACO is much better at finding non-dominated solutions than the other two algorithms, except for the two examples MK03 and MK08. According to the box plot in Figure 11, these relationships are even more evident. The red line in the graph represents the mean value line. To ensure more transparent comparison results, Pareto front plots were drawn to demonstrate the performance of IACO, NSGA-II and IGA. To confirm our experiments, Figure 12 shows the Pareto front plots for the three algorithms based on the MK02, MK07, 03a instances. These three instances were chosen because they are three examples of different scales. We can clearly see the ability of the algorithm to find non-dominated solutions from the comparison of these three instances of different scales. Each colored dot in the Figure 12 represents a non-dominated solution, and from the figure we can clearly see that IACO has many more dots than the other two algorithms. Thus it shows that the algorithm proposed in this paper is stronger than the other two algorithms in finding non-dominated solutions on this instance.

To find out the specific cause of SC(HA, A1) > SC(A1, HA), SC(HA, A2) > SC(A2, HA), We further compare the convergence of IACO and IGA, and NSGA-Ⅱ. As shown in Table 9, the value in the table is the optimal value for each objective in the first of 10 experiments, where f₁ represents the makespan, f₂ represents the total workload, f₃ represents the workload of the critical machine, and f₄ represents the penalties of earliness/tardiness. As can be seen from Table 9, the f₁ of IACO is better than IGA, and NSGA-Ⅱ on 20 instances; The f₂ of IACO is better than IGA, and NSGA-Ⅱ on 24 instances; The f₃ of IACO is better than IGA, and NSGA-Ⅱ on 24 instances; The f₄ of IACO is better than IGA, and NSGA-Ⅱ on 25 instances. All four target values of IACO in these instances outperform the other two algorithms. For example, on MK05, MK09, MK10, 01a, 02a, 03a, 06a, 08a, 09a, 12a, 14a, 15a, 18a. The convergence curves of the three algorithms for the four objectives on instance MK07 and 03a are given in Figures 13 and 14. We can see from the Figure 13 that IACO outperforms the other two algorithms in terms of convergence on the other three objectives, except for f₁ where it performs slightly worse than IGA. From the Figure 14, we can see that the IACO proposed in this paper is better than the other two algorithms in terms of convergence of all four objectives in 03a instance. To demonstrate our experimental process, we give the Gantt charts obtained by applying the algorithm proposed in this paper at three different scale instances. The light blue boxes in the Gantt chart represent transportation time and the light brown boxes represent setup time. The numbers above the different boxes represent different job numbers, and the different operation numbers are indicated in the order in which they appear. The results show that the IACO structure is more suitable for minimizing the makespan, workload of the critical machine, and penalties of earliness/tardiness.

In summary, these results show that the proposed IACO has better performance than IGA and NSGA-Ⅱ in solving the MaOFJSP.

5.   Conclusions and future work

In this paper, an improved ant colony optimization is proposed to solve the multi-time-constrained multi-objective flexible job shop scheduling problem (MaOFJSP) with transportation time, setup time and delivery time. The four objectives the makespan, total workload (TW), workload of the critical machine (WCM), and penalties of earliness/tardiness (PET) were also optimized. Combining the problem characteristics, we design a distributed coding approach. The different machine neighborhoods are searched by iteratively updating the machine distribution part in the early stage of the algorithm. The local search is later performed by the improved ant colony optimization in this paper for the different neighborhoods assigned by the machine. The resulting scheduling set is sorted and filtered using entropy weight method and non-dominated sorting. In addition, we again proposed mutation and closeness operations for the machine assignment in order to improve the diversity of the population. Finally, the algorithm is evaluated and the effectiveness of the proposed algorithm is verified by experimenting with 28 benchmark instances of different sizes.

The study in this paper shows that scheduling with multiple time constraints will take more time to complete production. And this extra time is not negligible. Therefore, we hope that the research in this paper can be a good guide for the managers of shop production. In the future, our research will focus on the application of multi-objective FJSP in dynamic production scheduling environments. In addition, we can also use some emerging machine learning techniques to further improve the performance of the algorithm, such as generating high-quality initial populations by a two-stage method with the help of recurrent neural networks (RNN). Adjusting the algorithm parameters by dynamic control of reinforcement learning (RL) allows the algorithm to be self-adjusting.

Acknowledgments

This paper presents work funded by the National Natural Science Foundation of China (Nos. U1904167, 51905494), Innovative Research Team (in Science and Technology) at the University of Henan Province (No. 21IRTSTHN018), Henan Province Philosophy and Social Sciences Planning Project (No. 2019BZX017), and Zhengzhou Science and Technology Collaborative Innovation Project (21ZZXTCX19).

Conflict of interest

The authors declare no conflict of interest.