Development of a Multi-Objective Optimization Model Using NSGA-II for Flow Shop Production Machine Scheduling to Support Manufacturing Industry

Production machine scheduling in the industrial world, whether in manufacturing, production, or agro-industry, plays a crucial role in decision-making [1]. Companies strive to achieve the most effective and efficient scheduling to increase productivity while minimizing total cost and time [2]. One model that can be applied to made-to-order manufacturing environments is the flow shop scheduling model [3]. Production machine scheduling problems typically involve the arrangement and management of jobs to be processed on a series of machines. One of the difficulties in arranging and managing jobs across available machines is the difficulty of finding appropriate techniques to create an optimal production machine scheduling model that meets all established scheduling criteria.

The scheduling model developed in this study currently assumes ideal conditions where all machines are continuously available and each worker has homogeneous skills. To improve the practicality and robustness of the solution, sensitivity analysis and robustness testing are recommended. Sensitivity analysis can be performed by varying critical parameters, such as machine processing time, breakdown probability, or worker productivity, and then evaluating the changes in objectives, such as makespan or total tardiness.

In addition, robustness testing can be performed through Monte Carlo simulations incorporating random distributions on critical variables, or through robust optimization methods, where the optimization objective is expanded to minimize the worst-case impact of disruptions. In this way, the solution is not only nominally optimal but also robust to uncertainty, thus providing a clear trade-off between efficiency (makespan) and schedule stability. Integrating robustness metrics into multi-objective models, for example as additional objectives in NSGA-II or AOF, will yield Pareto solutions that reflect the balance between performance and practicality under real-world production conditions.

In flow shop scheduling, there are a number of jobs, each of which has the same sequence of tasks on the machines. A schedule can be modeled as a flow shop scheduling problem if the sequence of tasks is aligned [4]. Research on production machine scheduling has been conducted using the Johnson and Campbell algorithm [5]. Flow shop scheduling has evolved from single-objective (optimization with a single function) to multi-objective (optimization with multiple objective functions) [6]. In multi-objective cases, a set of optimal solutions will be produced, known as Pareto-optimal solutions and multi-objective evolutionary algorithm framework [7, 8].

Several studies have developed models related to flow shop production machine scheduling, using both single-objective and multi-objective optimization. Research on production machine scheduling with a single objective function that optimizes makespan value has been conducted by Ying et al. [9], Xu et al. [10] and Karacan et al. [11]. Research on production machine scheduling conducted by Kurt [12] relates to flow shop scheduling with 2 objective functions, namely completion time and tardiness, but has a fairly long computation time. Nejjarou et al. [13] and Ezugwu [14] studied production machine scheduling with 2 objective functions, namely total weighted completion time and makespan [15]. We studied two objective functions, namely makespan and total flow time [16, 17]. Research on production machine scheduling with two objective functions: makespan and total tardiness has been conducted by Wang et al. [18] and Mousavi et al. [19]. In general, all of the above studies demonstrated good computational performance, capable of mathematically formulating objective functions for flow shop scheduling. However, most of the methods used were only capable of solving dependent problems, meaning they could only be used on specific problems depending on the type of problem (heuristic). Research on optimization and mathematical patterns has been carried out by research teams [20, 21]. Research on multi-objective was also conducted by Hutagalung and Azlan [22] using the AHP method.

In the case of multi-objective optimization, the Non-Dominated Sorting Genetic Algorithm for multi-objective Optimization: NSGA-II algorithm is a group of Metaheuristic Algorithms that have been tested for reliability compared to other multi-objective optimizations. NSGA-II is a development method of Genetic Algorithm (GA) and NSGA. Compared with GA and NSGA, NSGA-II is distinguished by the use of crowding distance operators to produce better Pareto optimal solutions. Research using NSGA-II has been conducted by Labidi et al. [23] and Ransikarbum et al. [24], who studied multi-objective for vehicle routing problems, Nielsen et al. [25] and Araújo et al. [26] conducted a study of multi-objective optimization cases in portfolio management cases and Preuß et al. [27] conducted a study that was able to create a new method called omni optimizer which was adopted from NSGA-II for both single and multi-objective optimization cases.

NSGA-II (Non-dominated Sorting Genetic Algorithm II) is a multi-objective evolutionary algorithm used to solve optimization problems with more than one objective function (multi-objective optimization problems). NSGA-II was introduced by Deb and Tiwari [28] as an improvement on the previous NSGA algorithm.

In multi-objective algorithm-based research (e.g., NSGA-II, MOEA/D, AOF), computational complexity and real-time performance are very important aspects, especially if the final target is implementation in industry. Without this analysis, it is difficult to assess the feasibility of adoption because the algorithm's performance in the real world is greatly influenced by time and resource constraints. In conducting Real-Time Performance experiments, this can be done by measuring real-world performance, recording: (1) Average iteration time (ms): how long it takes one generation to complete, (2) Total runtime (s): total time until convergence, and (3) Scalability: how the runtime increases with larger populations/generations.

Advantages of NSGA-II [29, 30]: (1) Efficient time complexity (O(MN²), M = number of objectives, N = population size), (2) Eliteness mechanism, (3) Maintains solution diversity with crowding distance.

In this study, the comparison was only conducted between NSGA-II and AOF. Benchmarks against other widely used multi-objective algorithms, such as MOEA/D or SPEA2, were not conducted, so the superiority of NSGA-II in this context is limited to the selected test cases.

To demonstrate the universal superiority of NSGA-II, additional benchmarks against other multi-objective algorithms, such as MOEA/D or SPEA2, are needed, which could be a potential direction for future research.

Based on the above problems, it is necessary to conduct research to analyze the multi-objective NSGA-II algorithm in flow shop production machine scheduling to optimize 2 objective functions, namely makespan and total tardiness, so as to provide a set of alternative solutions for decision makers.

3.1 Mathematical model

The formulation of the makespan minimization problem with a mathematical model is carried out by Rego et al. [32] and Chen et al. [33] the following:

$\mathrm{M}(\mathrm{S})=\mathrm{C}_{\mathrm{m}, \mathrm{n}}$            (1)

$C_{m, n}=\max \left\{C_{m-1, n}, C_{m, n-1}\right\}+t_{m, n}$        (2)

where,

M(S) = makespan (time to process all jobs to completion)

C_m,n = completion time of job m on machine n

t_m,n = processing time for job m on machine n

m = total number of jobs in the schedule

n = total number of machines in the schedule

In addition to formulating makespan, Huang et al. [34] formulated the problem of minimizing total tardiness with the following mathematical model:

$\mathrm{T}=\sum_{\mathrm{i}=0}^{\mathrm{n}} \max \left[0, \mathrm{C}_{\mathrm{i}}-\mathrm{d}_{\mathrm{i}}\right]$              (3)

where,

T = total tardiness (the total time a job can be completed earlier than its due date)

C_i = the complete time for all machines to complete a process on job i

Therefore, Eqs. (1)-(3) can be modeled as a multi-objective problem with the following mathematical model:

$\begin{aligned} & \operatorname{Min} \cdot \mathrm{Z}_{\mathrm{M}}=\max \left\{\mathrm{C}_{\mathrm{m}-1, \mathrm{n}}, \mathrm{C}_{\mathrm{m}, \mathrm{n}-1}\right\}+\mathrm{t}_{\mathrm{m}, \mathrm{n}} \\ & \operatorname{Min} \cdot \mathrm{Z}_{\mathrm{T}}=\sum_{i=0}^n \max \left[0, C_i-d_i\right]\end{aligned}$             (4)

where,

m = total number of jobs in the schedule

n = total number of machines in the schedule

t_m,n = processing time for the mth job on the nth machine

$C_i$  = complete time for all machines to complete a process on job i (i = 1, ..., m)

3.2 Implementation of NSGA-II stages

The NSGA-II multi-objective optimization used to find a solution to the mathematical model of Equation 4 is divided into several stages as follows:

(1) Population Initialization

The population is initialized randomly, but must still comply with the constraints of the mathematical model.

(2) Non-Dominated Sort

The initialized population is then sorted based on non-dominance using the Fast Non-dominated Sorting algorithm defined by Deb as follows:

(a) For each individual p in population P, perform the following steps:

(i) Initialize Sp =$\varnothing$. This set will contain all individuals dominated by p.

(ii) Initialize n_p = 0. This indicates the number of individuals that dominate individual p.

(iii) For each individual q in P

·If p dominates q, add q to the set Sp. $\left(S_p=S_p \cup q\right)$

·If q dominates p, increment the dominance counter n_p. (n_p=n_p+1)

(iv) If n_p = 0, this means no individual dominates p, so p is in the first front; Assign p a rank of one (p_rank = 1). Add p to the first front $\left(F_1=F_1 \cup\{p\}\right)$

(b) Initialize the front counter with one (i = 1)

(c) As long as the i-th front is not empty (Ƒ_i ≠$\varnothing$), do the following:

(i) Q = $\varnothing$

(ii) For each individual p in front Ƒ_i

For each individual q in Sp

n_q = n_q - 1, decrement the dominance count for individual q

(iii) If n_q = 0, then no individual in the next front dominates q, so q_rank = i + 1. Add q to Q $(Q=Q \cup q)$.

(iv) Increment the front counter by one (i=i+1)

(v) Q becomes the next front (Ƒ_i = Q)

(3) Crowding Distance

Crowding distance is used as a comparison between two individuals in the same front, so that the resulting solution can represent the overall Pareto-optimal solution.

The method for calculating crowding distance is defined by Zheng et al. [35] and Zheng et al. [36] as follows:

(a) For each individual in front Ƒ_i, where n is the number of individuals

(b) Initialize the distance to zero for all individuals (Ƒ_i(d_j)=0), where j represents the jth individual in front Ƒ_i

(c) For each objective function m

‪(i) Sort the individuals in front Ƒ_i i based on the objective value m

‪(ii) Assign infinite distances to the first and last individuals (I(d₁)=$\infty$ and I(d_n)=$\infty$)

‪(iii) For k ranging from 2 to (n-1)

$I\left(d_{\mathrm{k}}\right)=I\left(d_{\mathrm{k}}\right)+\frac{I(k+1) \cdot m-I(k-1) \cdot m}{\sum_m^{max }-\sum_m^{min } 1}$

I(k).m is the objective value m for individual k in I

(4) Selection

The roulette wheel selection method used is Roulette Wheel Selection. In the roulette wheel method, the parent chromosomes to be crossed over are selected based on their fitness values.

The fitness value of each chromosome is divided by the total fitness value of all chromosomes in the population. Each chromosome is considered a piece of the roulette wheel, with the size of the piece proportional to its fitness value. A target value between zero and one is randomly assigned. The roulette wheel is then spun N times, where N is the number of individuals or chromosomes in the population. In each spin, the chromosome with a fitness value below the target value is selected to become the parent for the next generation.

(5) Genetic Operators

The genetic operator used is Precedence Preservative Crossover (PPX) for crossover, and the mutation operators are removing and insert mutations. PPX can be explained as follows:

(a) A new string is randomly constructed from the alleles of the parent strings.

(b) A random number 1 or 2 is used to select the parent.

(c) f 1 is derived from the leftmost allele of the first parent, if 2 is derived from the leftmost allele of the second parent. The selected allele is then removed from both parents.

(d) The process continues until the characters in both parents are exhausted.

For example, two parents ABCDEF and CABFDE, and the random number 121122, will produce the new string ACBDFE.

In remove and insert mutations, one locus is randomly selected and the character at that position is deleted. A new locus is selected, and the previously deleted character is inserted.

(6) Recombination

The crossover and mutation population is then combined with the parent population, which is then selected using Non-Dominated Sort and Crowding Distance to obtain the next generation population.

3.3 Research results

The data used in this study is taken from Eva Instances: http://www.upv.es/gio/rruiz. The data consists of jobs ranging from 50 to 350 and machines ranging from 10 to 50.

Several samples will be taken from each job and machine. The data will then be tabulated using the format shown in Table 1.

Table 1. Production machine scheduling data table

3.3.1 Optimizing production machine scheduling using AOF

The AOF method is a method for generating solutions that minimize makespan and total tardiness. In principle, this method applies single-objective optimization, such as the Genetic Algorithm, to multi-objective problems.

The two objective functions in Eq. (4) are combined into a new objective function using the addition operator (+), and each function is assigned a weight between 0 and 1, with the sum of the two weights necessarily equal to 1. The new objective function is the result of combining the two objective functions.

$C_\mathrm{m, n}=\max \left\{\mathrm{{C}_\mathrm{m-1, n}}, \mathrm{{C}_\mathrm{m, n-1}}\right\}+\mathrm{{t}_{m, n}}$             (5)

Eq. (4) is: $\operatorname{Min} Z=\alpha^*\left(\max \left\{\mathrm{C}_{\mathrm{i}-\mathrm{1}, \mathrm{j}}, \mathrm{C}_{\mathrm{i}, \mathrm{j}-1}\right\}+\mathrm{t}_{\mathrm{i}, \mathrm{j}}\right)+(1-\alpha) * \sum_{i=0}^n \max \left[0, C_i-d_i\right]$ where $\alpha$ is a weight with a value between 0 and 1 .

To optimize Eq. (4), a Genetic Algorithm is used, with the determination of the weight $\alpha$ carried out through research with variations of $\alpha$ = 0.1, 0.5 and 0.9, the data for which is tabulated in Tables 2 to 4.

Table 2. Production machine scheduling data table

Table 3. AOF research results instance I_150_30

Table 4. AOF research results instance I_250_50

AOF weights can be dynamically adjusted. This approach leverages the distribution of solutions on the Pareto front, allowing the weights to reflect the trade-off balance between objectives that naturally emerges from the evolutionary algorithm. Furthermore, weights can be determined based on the decision maker's explicit preferences, for example through the integration of an interactive decision support system (DSS). This mechanism ensures that the selected solution is not only computationally optimal but also more closely aligned with practical needs and industry strategy. The dynamic weight approach also offers significant advantages. By adjusting weights based on the distribution of the Pareto front, the method is able to represent the trade-off between objectives in a more natural and adaptive manner.

To analyze the results of the AOF study above, normality and homogeneity tests were first conducted to determine the appropriate test for comparing the average results. The normality tests used were the Kolmogorov-Smirnov and Shapiro-Wilk tests, while the homogeneity test used Levene's test.

Research results with a normal distribution and homogeneous variances were compared using the ANOVA test, while those with a normal distribution and non-homogeneous variances were compared using the Welch test. For non-normally distributed results, the Kruskal-Wallis test was used. The statistical tests above were calculated using SPSS software, with the analysis results tabulated in the table below.

The results of the normality test in Table 5 with the Kolmogorov-Smirnov and Shapiro-Wilk tests show that the results of the AOF research for ZM, ZT and Generation at $\alpha$ = 0.1, 0.5 and 0.9 are normally distributed with a significance level of 0.05.

Table 5. Results of AOF normality test

The results of the homogeneity test in Table 5 using Levene's test indicate that the AOF research results for ZM, ZT, and Generation at α = 0.1, 0.5, and 0.9 have homogeneous variances with a significance level of 0.05 (the Sig. column in Table 6 is greater than 0.05).

Table 6. AOF research results instance I_350_50

Based on the results of the two tests above, the AOF experimental results are normally distributed with homogeneous variances. Therefore, the ANOVA test was used to compare the average values of ZM, ZT, and Generation at α = 0.1, 0.5, and 0.9 with the test results tabulated in the Table 7.

Table 7. Results of AOF homogeneity test

The ANOVA test results show that the average values of ZM, ZT, and Generations do not differ significantly for α=0.1, 0.5, and 0.9 at a significance level of 0.05 (the Sig. column in Table 7 is greater than 0.05). This means that changes in α do not significantly affect changes in the average values of the first objective function (ZM), the second objective function (ZT), or the number of generations. The Genetic Algorithm for AOF optimization applies the same PPX and Remove and Insert operators as the Genetic Operators for the NSGA-II stage with parameters tabulated in Table 8.

Table 8. Results of ANOVA AOF test

To show the population changes related to the increase in generations, Figure 3 displays the Genetic Algorithm population for production machine scheduling at the time of population initialization, the 50th, 125th, 250th and 294th generations. It turns out that the value of the 1st objective function (ZM) and the value of the 2nd objective function (ZT) are always better than the previous generation.

31.png

32.png

33.png

34.png

To strengthen the validity of the results, clarification of the parameter tuning method is needed. One commonly used approach is grid search, where parameter combinations are systematically tested to identify the best configuration based on performance indicators (e.g., hypervolume or spread). By trying several parameter combinations (e.g., PopSize = {50, 100, 200}, Crossover = {0.8, 0.9, 1.0}, Mutation = {0.05, 0.1, 0.2}. These parameter combinations are presented in Table 9.

Table 9. Grid search parameter NSGA-II

Another, more adaptive alternative is automatic parameter tuning, such as adaptive parameter control or self-adaptive evolutionary strategies, which dynamically adjust parameter values during the evolution process. Implementing these tuning methods makes experimental results more credible, replicable, and relevant for real-world applications in industrial environments.

3.3.2 Optimizing production machine scheduling using NSGA-II

To efficiently obtain the best set of Pareto-optimal solutions, experiments are needed to determine the smallest number of generations required to obtain them.

Applying Pareto solutions in the context of production scheduling presents its own challenges, as optimization results generated by multi-objective algorithms are often only displayed in tables or static graphs.

Therefore, integrating Pareto solutions into a Decision Support System (DSS) or interactive tool is a strategic step to ensure that the resulting solutions are not only theoretical but also operational. Future research will attempt to utilize an Integrated DSS.

Through a DSS, Pareto solutions, which represent trade-offs between objectives (e.g., minimizing makespan, total delay, and increasing machine utilization), can be visualized in interactive graphs such as scatter plots or parallel coordinate plots. With this display, decision-makers can understand the distribution patterns of solutions and select the alternative that best aligns with company priorities.

The quality of the Pareto-optimal solution sets from different generations is compared based on the average estimate of the hypervolume indicator.

The hypervolume indicator is estimated using a Monte Carlo approach [37-39], which normalizes all objective function values between 0 and 1, then generates a random set of objective function values and tests each random objective function value to see whether it is dominated by one of the Pareto-optimal solutions.

The research for hypervolume estimation was conducted three times for each instance in the same generation, the results of which are tabulated in Table 10.

Table 10. Genetic algorithm parameters for AOF

Table 11. Results of NSGA-II hypervolume estimation research

The normality test using SPSS on the data in Table 11 is tabulated in Table 12. The normality test assesses the normality of data distribution. This test is the most commonly used for parametric statistical analysis. Normally distributed data is a prerequisite for parametric tests. For data that does not have a normal distribution, non-parametric tests are used for analysis.

To show population changes related to increasing generations, Figure 4 displays the NSGA-II population at the time of population initialization, the 100th, 250th, 500th, 750th and 1000th generations.

Figure 4 shows that the 20 solutions in the population are increasingly evenly distributed, approaching the Pareto-optimal solution as the number of generations increases. As generations increase, the values of the first and second objective functions decrease simultaneously. It can be seen that of the 20 resulting solutions, not all are Pareto-optimal. In this example, only 12 Pareto-optimal solutions are on front 1 (symbolized by a blue '*'), while the other 8 solutions are on front 2 (symbolized by a red '+').

41.png

42.png

43.png

44.png

45.png

46.png

The population obtained in the 1000th generation of NSGA-II is tabulated in Table 12.

Table 12. Normality test results of NSGA-II hypervolume estimation

Table 13 is a table of NSGA II parameters, which consists of: Population size = 20. Advantages, fast computation per generation. Generations = 1000. Many generations good for allowing convergence. Pc = 0.9 common and reasonable value. Emphasizes the exploitation of gene combinations. Pm = 0.1 somewhat high for some representations (mutation is typically 0.001–0.1 depending on the representation).

Table 13. NSGA-II parameters

Table 14. NSGA-II population at the 1000th generation

Table 14 is the result of the NSGA-II population at the 1000th generation.

3.4 Comparison of AOF and NSGA-II

Solutions to multi-objective problems can be compared using dominance.

The i-th solution is declared better than the j-th solution if it dominates the j-th solution. The i-th solution dominates the j-th solution if all the objective function values of the i-th solution are no worse than the j-th solution, and at least one objective function value of the i-th solution is better than the objective function of the j-th solution.

Table 15 shows that a dominance comparison is performed between the AOF Genetic Algorithm solution for α = 0.1 and the NSGA-II solution. The symbol 'M' is used to indicate dominance, 'D' =100 to indicate dominance, and '-' to indicate non-dominance of the solutions.

Table 15. Comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_150_30

Figure 5 is a 2-dimensional image of the shape of comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_150_30.

5.png

Figure 5. Comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_150_30

Figure 6 is a 2-dimensional image of the shape of comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_250_50.

6.png

Figure 6. Comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_250_50

Figure 7 is a 2-dimensional image of the shape of comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_350_30.

7.png

Figure 7. Comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_350_30

The comparison of the dominance of AOF solutions with NSGA-II solutions in Tables 16 and 17 shows that there are always AOF solutions dominated by NSGA-II solutions, while no AOF solutions dominate NSGA-II solutions. This indicates that the Pareto-optimal NSGA-II solution successfully provides a better solution than the AOF solution.

Table 16. Comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_250_50

Table 17. Comparison of the dominance of the AOF solution with the NSGA-II solution for instance I_350_30

Based on the comparison above, it is clear that the population of Pareto-optimal NSGA-II solutions is superior, successfully dominating 100% of the AOF system solutions. Similarly, the comparison of NSGA-II solutions with AOF solutions shows that the population of Pareto-optimal NSGA-II solutions is superior, successfully dominating 100% of the AOF solutions.

The population of Pareto-optimal solutions can be used to provide managers with several alternative solutions. Pareto-optimal solutions are optimal multi-objective optimization solutions with trade-offs across all objectives, so all Pareto-optimal solutions are equally good. Therefore, managers still need to make informed decisions about which solution to implement in the company.