Multi-objective optimization of simultaneous buffer and service rate allocation in manufacturing systems based on a data-driven hybrid approach

The challenge presented by simultaneous buffer and service rate allocation in manufacturing systems represents a difficult non-deterministic polynomial problem. Previous studies solved this problem by iteratively utilizing a generative method and an evaluative method. However, it typically takes a long computation time for the evaluative method to achieve high evaluation accuracy, while the satisfactory solution quality realized by the generative method requires a certain number of iterations. In this study, a data-driven hybrid approach is developed by integrating a tabu search–non-dominated sorting genetic algorithm II with a whale optimization algorithm–gradient boosting regression tree to maximize the throughput and minimize the average buffer level of a manufacturing system subject to a total buffer capacity and total service rate. The former algorithm effectively searches for candidate simultaneous allocation solutions by integrating global and local search strategies. The prediction models built by the latter algorithm efficiently evaluate the candidate solutions. Numerical examples demonstrate the efficacy of the proposed approach. The proposed approach improves the solution efficiency of simultaneous allocation, contributing to dynamic production resource reconfiguration of manufacturing systems.


Introduction
Uncertain factors inherent in manufacturing systems, such as varying service rates and blocking, may cause a loss in their system performance.One option to decrease this loss is to allocate additional buffers and utilize machines with high service rates (Jiao et al., 2018;Zeid et al., 2021).However, too many buffers will increase system redundancy, leading to longer sojourn time and higher cost, whereas an extremely high service rate will increase machine idle time (Smith, 2018).Therefore, it is necessary to carefully configure the buffers and service rates of machines applied in manufacturing systems.Numerous studies have separately investigated the buffer allocation problem (Gao et al., 2020;Shi & Gershwin, 2016;Weiss et al., 2019) or service rate allocation problem (Shaaban & McNamara, 2009;D. Song et al., 1998) in manufacturing systems.Numerical results have shown that a minor improvement in the optimization of buffer or service rate allocation may contribute significantly to cost savings, increased profits, and improved working efficiency (Ng et al., 2017).However, these studies did not explore superior system performance compared to separate optimization owing to the positive relationship between buffers and service rates (Nahas et al., 2014;Su et al., 2017).The simultaneous allocation problem is more difficult than the sole allocation problem, particularly in the case of multi-objective optimization (Xi et al., 2021).To date, limited research has been conducted on the simultaneous buffer and service rate allocation problem (SBSRAP) for multi-objective optimization.(Cruz, 2009) maximized system throughput while simultaneously reducing the total buffer size and total service rate using a genetic algorithm.(Ng et al., 2017) proposed a hybrid algorithm to maximize the system throughput and minimize the average waiting time.(Smith, 2018) optimized the SBSRAP to achieve threshold throughput using a hybrid approach employing sequential quadratic programming and the generalized expansion method.Although previous studies have proposed algorithms to achieve sufficient solution quality, the computation time required to do so has not been further discussed; as a result, the tradeoff between solution quality and computation efficiency has not been addressed.Furthermore, previous scholars (Gao et al., 2019(Gao et al., , 2021;;Oljira et al., 2020;Z. Song & Moon, 2019) typically utilized simulation and analytical methods as evaluative methods to calculate system performance and determine better candidate solutions for SBSRAP.However, simulation methods require long setting and simulation times, whereas analytical methods still require considerable time when executing a large number of iterations to search with acceptable accuracy (Shi & Gershwin, 2016;Weiss et al., 2019).Data-driven methods were an option to improve evaluation efficiency (Gao & Liu, 2023;Tsadiras et al., 2013).However, according to our best knowledge, there have been no previous studies related to performance prediction of manufacturing systems considering simultaneous buffer and service rate allocation (SBSRA).
This article therefore proposes a data-driven hybrid approach in which the tabu search-non-dominated sorting genetic algorithm II (TS-NSGA-II) is integrated with the whale optimization algorithm-gradient boosting regression tree (WOA-GBRT) to maximize the throughput and minimize the average buffer level (ABL) of a manufacturing system subject to a total buffer capacity and total service rate.The TS-NSGA-II is applied to obtain SBSRA solutions of sufficient quality.Furthermore, performance prediction models are developed based on the WOA-GBRT to decrease the evaluation time of the throughput and ABL.
The contributions of this study are as follows: • Performance prediction models constructed by the WOA-GBRT efficiently predict the throughput and ABL, providing a new possibility in the field of rapid performance evaluation of manufacturing systems.• The proposed TS-NSGA-II balances global and local search abilities, improving the solution quality for the SBSARP in manufacturing systems.
The remainder of this article is structured as follows.Section 2 presents the problem statement.The solution methodology is addressed in Section 3. Numerical examples are then provided to demonstrate the effectiveness of the proposed hybrid approach in Section 4. Finally, conclusions and future research directions are summarized in Section 5.

Assumptions
The assumptions underlying the manufacturing systems considered in this study are as follows: Fig. 1.Model of a manufacturing system.There are  machines and  − 1 buffers in a manufacturing system, as shown in Fig. 1 (Spinellis et al., 2009), in which the circles represent machines, the squares represent buffers, and the arrows represent the directions of part movement.Each buffer has a finite capacity that denotes the temporary storage area.Each machine has a service rate that denotes the working efficiency.The service rate of a machine obeys an exponential distribution.

Notation
The notations in this study are shown in Table 1.an integer that is less than or equal to .

Optimization problem
The goal of the study is to maximize the throughput and minimize the ABL of a manufacturing system subject to total buffer capacity and total service rate.This problem can be expressed as follows: Find  = ( , ⋯  ⋯ ,  ) and ̅ = ( , ⋯  ⋯ ,  ) to maximize  and minimize , (2)  ≥ 1,    , (3)  ≥ 0.1, the accuracy of μi is to one decimal place. (4) where Eq. ( 1) and Eq. ( 2) regulate the decision variables  and ̅ , respectively, and Eq.(3) and Eq. ( 4) represent bounds for these respective decision variables.

Solution methodology
Fig. 2 presents the framework of the proposed hybrid approach.The TS-NSGA-II is proposed to achieve a trade-off between solution quality and computation efficiency.The WOA-GBRT is proposed to build performance prediction models based on historical data generated by the generalized expansion method (GEM) (Kerbache & MacGregor Smith, 1988).The prediction models can rapidly predict the throughput and ABL of a manufacturing system with acceptable accuracy.

WOA-GBRT-based prediction model
The GBRT effectively solves both classification and regression problems by employing ensembles of regression trees as weak learners to generate a power learner (Ke et al., 2017).Consequently, the GBRT was used to build performance prediction models.Furthermore, the WOA in the previous literature (Mirjalili & Lewis, 2016) was used to optimize the hyperparameters in the GBRT.

Performance prediction models
Generate candidate solutions of SBSRA Evaluate the candidate solutions and select better ones The proposed

WOA-GBRT GEM
Provide fitting algorithm Provide historical data, including SBSRA solutions and their corresponding throughput and ABL

Generation of throughput and ABL data for manufacturing systems
The GEM was used to calculate the throughput and occupation probabilities of the buffers in a manufacturing system.The buffer occupation level ( ) was evaluated based on the occupation probabilities, and the ABL () was calculated according to literature (Ng et al., 2017).Four patterns were considered to generate data sets of the throughput and ABL based on a manufacturing system comprising five or eight machines.As shown in Table 2, the parameters in each sample included  and ̅ , as well as the corresponding  and .Note that  and ̅ were generated randomly subject to the total buffer capacity and total service rate.20000  2 (5, 20, 5)  3 (5, 16, 7.5)  4 (8, 28, 8)

WOA-GBRT
The WOA was utilized to optimize two critical hyperparameters involving the learning rate  and n_estimators .The default values were used for the other hyperparameters.For details on the hyperparameter selection process, readers can refer to the sklearn documentation (Sklearn, n.d.).For the details on the WOA, readers can refer to literature (Mirjalili & Lewis, 2016).Fig. 3 shows the flowchart of the WOA-GBRT.

TS-NSGA-II
In the TS-NSGA-II, the NSGA-II is responsible for global search, and a diversification strategy is utilized to intensify the diversity performance.As the TS has been proven to be effective in obtaining near-optimal solutions (Glover et al., 2021;Papadopoulos et al., 2013) and capable of utilization with an evolutionary algorithm (Su et al., 2017), it was used to intensify the local search and extend the search scope in every iteration.Fig. 4 shows the flowchart of the TS-NSGA-II.

Coding and initialization
The coding can be defined as ̅ = ( ,  ,  , ⋯ ,  ,  ),  = 1, ⋯ , .Each candidate solution to the SBSRAP corresponds to an individual in the TS-NSGA-II.The  initial individuals in the TS-NSGA-II can be generated randomly with the constraints given in Eqs.(1-4).

Individual ranking
The non-dominated sorting II method (Deb et al., 2002) is used to rank individuals by dividing them into different floor ranks.For individuals in the same floor rank, the crowding distance is used to rank their priorities.A larger crowding distance denotes better priority.In this study, the crowding distance was calculated using a method presented in the literature (Su et al., 2017).Assuming an individual and its two adjacent individuals are marked as ,  − 1, and  + 1, respectively, the crowding distance of individual  can be calculated by (5)  The number of neighborhood individuals is problem-dependent and can be expressed as follows: The move representation can be defined as random changes in the values of the buffer capacities and service rates of two locations in a candidate individual.Assuming that ̅ = (⋯ ,  , ⋯ ,  , ⋯ ,  , ⋯ ,  ⋯ ) is an individual.Here, 1, 2, 1, and 2 mark randomly selected buffers and machines, respectively.A neighborhood solution can be expressed as follows: (7) After each iteration of the TS, the individuals in  are chosen to generate neighborhood individuals in the next iteration if the stopping criterion is not satisfied.
(2) Moved attribute and tabu tenure According to previous studies (Gao, 2022;Glover et al., 2021), the buffer and machine locations describing the change in capacity and service rate [(2, 1), (2, 1)] are selected as the moved attribute.The tabu tenure employed is expressed as follows: If the tabu list becomes full, the moved attribute that first entered the tabu list is removed.
(3) Stopping criterion for the TS A maximum number of iterations is set as the stopping criterion of the TS, denoted by  .

Crossover operation
individuals are selected from  based on the proportion  , and  neighborhood individuals are generated by iteratively selecting each pair of the selected  individuals to generate a neighborhood individual.Thus,  can be calculated b The specific crossover process is shown in Fig. 5.The crossover operations of the buffers and service rates are respectively described by  =  , (0,1) < 0.5 where  and  are the capacity of buffer  and service rate of machine  , respectively, in the first selected individual;  and  are the capacity of buffer  and service rate of machine  , respectively, in the second selected individual.

Fig. 5. Flowchart of crossover operation 3.2.5 Mutation operation
Mutation operations are conducted on individuals selected from  with the proportion  to generate  neighborhood individuals.The specific mutation process is shown in Fig. 6.The mutation values for the buffer and service rate are respectively given by where

Diversification strategy
The diversification strategy employed in this study is described as follows: If the number of Pareto optimal solutions has not changed within a certain number of iterations, 10% of the new individuals in Π will be replaced by randomly generated individuals.The value of the iteration number is problem-dependent and set to one-fifth of  .

Stopping criterion for the NSGA-II
A maximum number of iterations was set as the stopping criterion of the NSGA-II, denoted by  .

Numerical examples
The numerical examples in this study comprised two experiments.Experiment 1 was conducted to test the effectiveness of the proposed WOA-GBRT-based performance prediction model by comparing it to the polynomial regression (PR), decision tree (DT), random forest (RF), and ANN prediction algorithms.To make the DT, RF, and WOA-GBRT more comparable, all hyperparameters for the decision trees in the three algorithms were the same.Furthermore, an ANN model with two hidden layers (Huang, 1999;Tsadiras et al., 2013) was used, and its number of neurons was determined based on a comparison of different ANN models.
Experiment 2 evaluated different manufacturing systems with five machines.The TS-NSGA-II was compared to the classical NSGA-II to determine the effectiveness of the TS-NSGA-II for manufacturing systems.To ensure that the results of the two algorithms were comparable, the coding, initialization, crossover operation, mutation operation, and stopping criterion used in the NSGA-II were the same as those used in the TS-NSGA-II.Moreover, the enumeration algorithm (EA) was used as the benchmark to obtain Pareto solutions and demonstrate the absolute solution quality of the TS-NSGA-II.In addition, three patterns were considered.Each pattern consists of 20 replicate calculations because of random initialization, and their mean value (MV) and median absolute deviations (MADs) are used to evaluate algorithm performance and robustness.
The WOA-GBRT, PR, DT, RF, ANN, TS-NSGA-II, NSGA-II, and EA were all written in Python 3.5.1 and executed for all experiments on a laptop computer with a 3.2 GHz Intel Core 4 CPU.In addition, the performance prediction models were trained using the sklearn library (version 0.19.0).

Comparison indicators in Experiment 1
The predictive values considered in this study were the system throughput  and ABL .The MSE, mean absolute error (MAE), and decision factor ( ) (Sklearn, n.d.) were determined for both of these values, and their averages served as the final evaluation measures.The MSE, MAE, and  values should all be in the 0-1 range; low MSE and MAE values and high  values indicate superior predictive performance.Readers can refer to the Sklearn document (Sklearn, n.d.) for details.

Comparison indicators in Experiment 2
The number of Pareto optimal solutions ( ), the number of Pareto solutions ( ^), average distance (ℎ ), and Zitlzler measure (Γ , ) were used to evaluate convergence performance.High  and  ^, and low ℎ indicate superior convergence.Furthermore, Γ , denotes the percentage of solutions obtained by algorithm 1 dominated by or equal to at least one solution obtained by algorithm 2.Algorithm 1 is considered to have better convergence than algorithm 2 if Γ , is lower than Γ , .Readers can refer to literature (Su et al., 2017) for the details of  ,  ^, and ℎ , and literature (Huang, 1999) for the details of Γ , .
The diversity metric (Δ) and limit distance () were used to evaluate the algorithm diversity performance.A smaller Δ represents a more uniform solution distribution, and  indicates the distribution of the solution scope.Readers can refer to literature (Su et al., 2017) for the details of Δ and .Additionally, the computation time  was evaluated for all methods to compare their computation efficiencies.

Setting
Four performance prediction models were constructed based on different total buffer capacities, total service rates, and numbers of machines, as shown in Table 2.For the initial parameters of the WOA-GBRT,  was set to 10 and  was set to 30 for all four patterns.

Results and discussion
Fig. 7 and Fig. 8 show the predicted throughput and ABL, respectively, obtained using the prediction models and the actual values in the testing samples.The closer the plots in each figure to the diagonal line passing through the origin, the smaller the difference between the predicted and actual values in the testing samples, indicating a higher predictive accuracy.• The WOA-GBRT achieves the best MAE, MSE, and  values among the models tested.The MAE and MSE of the WOA-GBRT are close to 0, whereas the  is close to 1, indicating that the WOA-GBRT provides the highest predictive accuracy.• The predictive performance of the DT is unsatisfactory.However, its regression performance can be significantly improved by using it as a weak learner ensembled in large numbers to realize a learner as strong as the RF and GBRT.The final prediction results of both the RF and GBRT are obtained by analyzing the results of several decision trees, which iteratively reduces the predictive error.Consequently, the predictive accuracy and robustness of these models can be improved.• The total buffer capacity, total service rate, and the number of machines increase from Patterns 1 to 4, enlarging the scope of solutions and making the training data increasingly complex.This increases the difficulty of building performance prediction models, decreasing the prediction accuracy of the WOA-GBRT.Furthermore, compared to the limited impact of the change in total buffer capacity or total service rate, an increase in the number of machines is found to severely decrease the predictive accuracy of the WOA-GBRT, as indicated by a comparison of Patterns 1-3 with Pattern 4. Consequently, it will be challenging to build prediction models with high predictive accuracy for large-scale manufacturing systems.

Setting
The manufacturing systems comprising five machines (Patterns 1-3) were evaluated in Experiment 2. Owing to the insufficient predictive accuracy of previously constructed prediction models at this stage, the manufacturing system involving eight machines (Pattern 4) was not tested.The input parameters of Patterns 1, 2, and 3 are listed in Table 2.The parameters of the NSGA-II and TS-NSGA-II are listed in Table 4.

Table 4
Parameter settings for the NSGA-II and TS-NSGA-II.

Results and discussion
Table 5 and Fig    The number of Pareto solutions obtained by the TS-NSGA-II is 340% greater than that obtained by the NSGA-II for Pattern 1.With increasing total buffer capacity and total service rate, the solution scope is extended, and the number of Pareto solutions obtained by both the TS-NSGA-II and NSGA-II increase.In addition, the numbers of Pareto solutions obtained by the TS-NSGA-II for Patterns 2 and 3 were respectively 307% and 279% greater than those obtained by the NSGA-II.This result demonstrates the improved convergence performance of the TS-NSGA-II and implies that it can still achieve superior convergence as the problem scale increases.
The average distances of the TS-NSGA-II are 49%, 16%, and 27% of those of the NSGA-II for Patterns 1, 2, and 3, respectively, indicating the superior convergence of the TS-NSGA-II.Consequently, the Pareto optimal solutions obtained by the TS-NSGA-II are shown in Fig. 9 to be closer to the Pareto front obtained by the EA than those obtained by the NSGA-II.
The Zitlzler measure of the TS-NSGA-II for each pattern is lower than that of the NSGA-II, which also demonstrates the superior convergence performance of the TS-NSGA-II.Furthermore, the difference in the Zitlzler measure of the NSGA-II and the TS-NSGA-II for Pattern 1 is determined to be 0.66 and increased by 12% as the solution scope increased for Pattern 2. The Zitlzler measure of the TS-NSGA-II and the NSGA-II for Pattern 1 is 0.07 and decreased by 42% in Pattern 3.This finding also proves that the convergence performance of the TS-NSGA-II is superior to that of the NSGA-II, especially for larger problems.
(2) Diversity analysis The diversity metrics of the NSGA-II are 89% of those of the TS-NSGA-II for Patterns 1 and 2, respectively, indicating the superior diversity performance of the NSGA-II.The TS searches for candidate individuals around local areas.However, the number of new individuals retained in each update iteration is limited, leading to some individuals being around local areas and thereby decreasing the diversity.Consequently, the diversity performance of the TS-NSGA-II becomes unsatisfactory.For Pattern 3, the increase in the total service rate and extension of the solution scope increases the diversity metric of the TS-NSGA-II to greater than that of the NSGA-II because the TS-NSGA-II can detect more Pareto optimal solutions in the larger-scale problem, and the effect of the local search on the diversity performance decreases accordingly.The limited distances of the TS-NSGA-II are 94%, 100%, and 93% of those of the NSGA-II for Patterns 1, 2, and 3, respectively.As the TS focuses on Pareto optimal solutions around local areas, the search scope of the TS-NSGA-II is limited, resulting in a lower limited distance performance for Patterns 1 and 3.
( (5) Insight The proposed hybrid approach integrates local and global search strategies, balancing exploration and exploitation and generating superior convergence performance and computation efficiency than the NSGA-II.
The diversity performance of the proposed approach decreases slightly compared to the NSGA-II.However, the difference between the TS-NSGA-II and NSGA-II is mainly situated in the low limit of the Pareto optimal front, which is hardly used in actual application.

Conclusions and future work
This article presented a data-driven hybrid approach integrating the TS-NSGA-II and WOA-GBRT to maximize the throughput and minimize the ABL of a manufacturing system.In the proposed hybrid approach, the WOA-GBRT-based performance prediction model can quickly calculate the throughput and ABL with acceptable accuracy; the TS-NSGA-II was used to effectively search for SBSRA solutions.To verify its effectiveness, two numerical examples were presented in which the results obtained using the proposed WOA-GBRT were compared to those obtained using the PR, DT, RF, and ANN, and the results obtained using the TS-NSGA-II were compared to those obtained using the NSGA-II and EA.These examples showed that the WOA-GBRT achieved high predictive accuracy for the throughput and ABL evaluation, and that the TS-NSGA-II efficiently solved the SBSRAP.
The primary objective of the study is to support efficient multi-objective optimization of manufacturing systems under multiple decision variables.Furthermore, the approach proposed to do so also establishes a foundation for efficient performance evaluation based on machine learning methods.
The WOA-GBRT was used to build high-accuracy performance prediction models for a manufacturing system with five machines in this study.One direction for future research is to construct performance prediction models for large-scale manufacturing systems that provide sufficient accuracy.

Fig. 7 .Fig. 8 .
Fig. 7. Comparison of predictive throughput and actual values of the PR, DT, RF, ANN, and WOA-GBRT Fig. 8.Comparison of the predictive ABL and actual values of the PR, DT, RF, ANN, and WOA-GBRT . 9 compare the performances of the NSGA-II and TS-NSGA-II.Fig.10shows the MVs of different evaluation indicators for NSGA-II and TS-NSGA-II, respectively.Fig.11presents the number of Pareto solutions obtained by the NSGA-II and TS-NSGA-II.

Fig. 11 .
Fig. 11.Number of Pareto solutions obtained by the NSGA-II and TS-NSGA-II ) Computation efficiency analysisThe computation times of the TS-NSGA-II are similar to those of the NSGA-II for Patterns 1 and 2, but as the solution scope increases for Pattern 3, the computation time of the TS-NSGA-II grows to 1.7 times that of the NSGA-II.However, according to Fig.11, the TS-NSGA-II obtains more Pareto solutions than the NSGA-II in the same computation time, although the NSGA-II can rapidly obtain solutions in the first few seconds.This results from the local search in the TS-NSGA-II, which improves the exploitation of candidate solutions in promising local areas.Moreover, the Pareto solutions obtained by the TS-NSGA-II, and NSGA-II may decrease with time.Because the number of new individuals retained in each update iteration is limited, some of the obtained Pareto solutions may be neglected in the individual ranking process owing to their low crowding distances.(4)Robustness analysis Although the MAD of  for the NSGA-II is lower than that for the TS-NSGA-II in Pattern 1, the TS-NSGA-II achieves lower MAD in other experiments, demonstrating the better robustness of the TS-NSGA-II.Because of the satisfactory global search ability of the TS-NSGA-II, it can reproduce a similar optimization process and obtain very close solutions.

Table 1
Notions in this study. is the best learning rate in the WOA-GBRT. number of trees, generally defined as n_estimators. is the best number of trees in the WOA-GBRT. mean squared error (MSE). is the best MSE in the WOA-GBRT.

Table 2
Throughput and ABL data generated by the GEM.

Table 3
MAE, MSE and  values of the PR, DT, RF, ANN, and WOA-GBRT

Table 3
summarizes the MAE, MSE, and  values of the PR, DT, RF, ANN, and WOA-GBRT.The following conclusions were drawn from this table:

Table 5
Performance comparison of the NSGA-II and TS-NSGA-II.