# Predictive model-based multi-objective optimization with life-long meta-learning for designing unreliable production systems

Ehsan Mahmoodi, Masood Fathi, Amos H.C. Ng, Alexandre Dolgui

► To cite this version:

Ehsan Mahmoodi, Masood Fathi, Amos H.C. Ng, Alexandre Dolgui. Predictive model-based multi-objective optimization with life-long meta-learning for designing unreliable production systems. Computers and Operations Research, 2025, 178, pp.107011. 10.1016/j.cor.2025.107011. hal-05081748

HAL Id: hal-05081748

https://nantes-universite.hal.science/hal-05081748v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Predictive model-based multi-objective optimization with life-long meta-learning for designing unreliable production systems

Ehsan Mahmoodi<sup>a, D</sup>, Masood Fathi<sup>a, b, D</sup>,*, Amos H.C. Ng<sup>a, b</sup>, Alexandre Dolgui<sup>c, D</sup>

## ARTICLE INFO

**Keywords:**

Multi-objective optimization

Simulation

Predictive model

Meta-learning

Buffer allocation

Resource allocation

## ABSTRACT

Owing to the realization of advanced manufacturing systems, manufacturers have more flexibility in improving their processes through design decisions. Design decisions in production lines primarily involve two complex problems: buffer and resource allocation (B&RA). The main aim of B&RA is to determine the best location and size of buffers in the production line and optimally allocate production resources, such as operators and machines, to workstations. Inspired by a real-world case from the marine engine production industry, this study addresses B&RA in high-mix, low-volume hybrid flow shops (HFSs) with feed-forward quality inspection. These HFSs can be characterized by uncertainties in demand, material handling, processing times, and quality control. In this study, the production environment is modeled via discrete-event simulation, which reflects the features of the actual system without requiring unreasonable or restrictive assumptions. To replace the expensive simulation runs, five widely used regressor machine learning algorithms in manufacturing are trained on data sampled from the simulation model, and the best-performing algorithm is selected as the predictive model. To obtain high-quality solutions, the predictive model is coupled with an enhanced non-dominated sorting genetic algorithm (En-NSGA-II) that incorporates lifelong meta-learning and features a customized representation and a variable neighborhood search. Additionally, a post-optimality analysis using a pattern-mining algorithm is performed to generate knowledge for allocating buffers and operators based on the optimization results, thus providing promising managerial insights.

## 1. Introduction

Hybrid flow shops (HFSs) comprise multiple sequentially organized production segments, where each segment contains single or multiple parallel machines with at least one segment containing parallel stations (Perez-Gonzalez and Framinan, 2024). Because of their flexibility, HFSs have been widely employed in industries involving customizable products, including the automotive, marine, aerospace, and electronics industries (Komaki et al., 2019; Shahvari and Logendran, 2018). Despite the advantages of HFSs, their performance and efficiency can be significantly affected by production flow disruptions caused by maintenance, machine failures, and process time fluctuations at workstations (Wu and Cao, 2022). These disruptions within a production system can result in severe throughput bottlenecks that adversely affect the overall efficiency and productivity of the manufacturing process. To mitigate these bottlenecks, the theory of constraints suggests that the strategic allocation of buffers and resources can improve throughput and thus the overall performance

of the system (Golmohammadi, 2015; Wu et al., 2020). Hence, the buffer allocation problem (BAP) (Puchkova et al., 2020) and resource allocation problem (RAP) (Li et al., 2023a) have garnered significant interest from both scholars and practitioners.

The BAP involves determining the appropriate storage space between two consecutive segments (Kassoul et al., 2022), whereas the RAP focuses on assigning appropriate resources such as machines, labor, and equipment to each workstation (Zhou et al., 2021). Over the past two decades, considerable research has been conducted to address BAP and RAP separately (Chang and Huang, 2022), with significant contributions from reconfigurable manufacturing systems (RMS) literature. RMS research has established sophisticated frameworks for resource allocation (Hammedi et al., 2024; Yelles-Chaouche et al., 2021), demonstrating capabilities to rapidly adjust both production capacity and functionality. While these advances provide valuable insights into system reconfiguration and resource optimization (Bahtat et al., 2023), few studies have simultaneously investigated the BAP

<footnote><sup>a</sup> Division of Intelligent Production Systems, School of Engineering Science, University of Skövde, Skövde, 54128, Sweden</footnote>

<footnote><sup>b</sup> Division of Industrial Engineering and Management, Department of Civil and Industrial Engineering, Uppsala University, Uppsala, Sweden</footnote>

<footnote><sup>c</sup> IMT Atlantique, LS2N-CNRS, La Chanterie, 4, rue Alfred Kastler, Nantes cedex 3, F-44307, France</footnote>

<footnote>* Corresponding author at: Division of Intelligent Production Systems, School of Engineering Science, University of Skövde, Skövde, 54128, Sweden. E-mail addresses: masood.fathi@his.se, fathi.masood@gmail.com (M. Fathi).</footnote>
and RAP (which can be named buffer and resource allocation, B&RA) to address HFSs (i.e., HFS-B&RA). Current HFS-B&RA studies primarily focus on machine-type selection, resulting in long-term, cost-intensive decisions. However, in advanced HFSs with mobile resources, the dynamic allocation of resources in accordance with production plans presents distinct optimization requirements that complement existing RMS frameworks while offering a more flexible and cost-effective approach for high-mix, low-volume production systems.

This study addresses a particular type of HFS featured with real-world characteristics inspired by a case in an engine manufacturing plant. The addressed HFS is characterized by complex and customizable products and benefits from mixed automation levels at various production stages. Owing to the difference in processing time, the investigated HFS is decoupled between stages in which the automation level changes, and material handling between these decoupled stages is performed by autonomous mobile robots (AMRs). Quality inspection is integral to this HFS as it increases the reliability of the products and reduces stochasticity in the processing time. In brief, this study investigates a complex HFS characterized by multiple distinguishing features. The system under consideration is multi-product, capable of processing diverse product variants. It exhibits machine unreliability, manifested through stochastic processing times and random failures, introducing significant operational variability. Furthermore, the HFS is unbalanced, with heterogeneous processing times across different workstations, potentially leading to bottlenecks and efficiency challenges. Notably, the system incorporates mixed levels of automation, presenting a hybrid environment where manual and automated processes coexist.

This study aims to maximize throughput and minimize the total number of buffers, operators, and AMRs in the above-explained HFS by simultaneously resolving the BAP and RAP (HFS-B&RA). Owing to the stochasticity and complexity of the HFS-B&RA, a discrete-event simulation (DES) model is used to estimate the throughput of the production line under alternative designs. Because executing and optimizing a simulation model for real-world problems is time-consuming, five commonly used machine learning (ML) algorithms in manufacturing industries (do Amaral et al., 2022), namely artificial neural networks (ANNs), gradient tree boosting (GBT), polynomial regression (PR), random forest regressor (RFR), and support vector machine (SVM), were trained on data sampled from the simulation model and compared against five performance measures. The best-performing regressor algorithm was then employed to construct a predictive model (PM) as an effective and agile replacement for the simulation model. The PM is coupled with an enhanced multi-objective optimization (MOO) algorithm to form a predictive model-based multi-objective optimization (PMB-MOO). The employed MOO algorithm is an enhanced version of the non-dominated sorting genetic algorithm (En-NSGA-II), featuring a variable neighborhood search (VNS) mechanism equipped with three local search strategies. Notably, it incorporates a Lifelong Meta-Learning (LLML) approach, which enables the self-tuning of algorithm parameters. This allows for a more intelligent exploration of the solution space.

The performance of En-NSGA-II is compared with those of the non-dominated sorting genetic algorithm (NSGA-II) and multi-objective particle swarm optimization (MOPSO) using a set of instances generated from a case study. The results indicate that the proposed PMB-MOO can obtain good solutions for HFS-B&RA in highly stochastic and complex environments within a reasonable computational timeframe. Additionally, a post-optimality analysis is performed on the optimization results, which provides important information and transferable rules concerning the allocation of resources and buffer sizes in an HFS.

In summary, this study contributes to the existing literature by proposing a PMB-MOO approach to address HFS-B&RA. This approach benefits from a high-performance PM and a customized optimization algorithm. The optimization algorithm is distinctively enriched with VNS and an LLML. Moreover, to provide practitioners with applicable rules for improving HFSs, the interconnections between the decision and

objective spaces are scrutinized and analyzed using a pattern-mining algorithm.

The remainder of this paper is organized as follows: Section 2 presents a comprehensive literature review pertaining to HFS-B&RA. Section 3 delineates the problem description, outlines the underlying assumptions, and introduces the real-world case study. Section 4 formally defines the optimization problem. Section 5 outlines the solution approach. The subsequent Sections 6, 7, and 8 elaborate on the key components of the solution approach: the simulation model, predictive model, and optimization algorithm, respectively. Section 9 details the computational experiments conducted, while Section 10 presents and analyzes the results. Section 11 discusses the managerial insights derived from these findings. Finally, Section 12 summarizes the key findings of the study and proposes directions for future research.

## 2. Literature review

This section reviews the most relevant literature on BAP and RAP in HFSs for a better understanding of the current state of the topic.

As a pioneering study, Dolgui et al. (2002) investigated a single-product BAP in an HFS to minimize the production cost by considering the buffer sizes and locations using the aggregation method and a genetic algorithm (GA). Computational studies demonstrated that the GA yields higher-quality results than the Monte Carlo method. Abu Qudeiri et al. (2008) addressed a problem similar to that of Dolgui et al. (2002) and proposed an optimal design for HFSs that considered buffer sizes in addition to the number of machines inside each workstation and the machine types. Abu Qudeiri et al. (2008) improved the performance of the GA using the gene family arrangement method. Colledani et al. (2008) investigated the BAP by considering a multi-product HFS and developing a hybrid aggregation-decomposition method to determine the maximum throughput of a production line. The proposed aggregation-decomposition approach coupled with an iterative algorithm provided accurate throughput calculation results but less precise average buffer level estimations.

Nahas et al. (2008) addressed redundant machines and buffers to optimize throughput subject to a total cost constraint using ant colony optimization (ACO). To evaluate the throughput, each pair of parallel machines was approximated as a single unreliable machine. An analytical decomposition approximation method was used throughout the optimization procedure to evaluate the production rate of each alternative configuration. In a subsequent study, Nahas et al. (2009) combined ACO with simulated annealing (SA) to extend their previous study. Massim et al. (2010) addressed the simultaneous allocation of buffers and the selection of machine types to obtain a desirable level of availability using an artificial immune algorithm (AIA) subject to the constraints of available machines. This algorithm selects system elements from a list of available machine types based on their cost and accessibility to satisfy customer demands. Luo et al. (2010) addressed the RAP for creating highly customized products with continually fluctuating orders. Luo et al. (2010) provided a simulation-based optimization approach to obtain the best RA strategy, guaranteeing that all tasks are completed within a specified time frame with the least resources and without buffer overflow. Specifically, simulation and a branch-and-bound (B&B) algorithm were combined to achieve the objective.

Amiri and Mohtashami (2012) developed a multi-objective formulation of the BAP in HFSs. The authors proposed a methodology for considering the general distribution functions of the parameters of a manufacturing line using a simulation-based approach. In their proposed methodology, a metamodel was proposed to estimate the throughput, and a hybrid evolutionary algorithm with a line-search method was used to provide near-optimal buffer allocations. Continuing the study by Amiri and Mohtashami (2012) and Mohtashami (2014) proposed a hybrid GA with a new dynamic mutation probability to solve simultaneous buffer sizing and machine allocation problems
in HFSs. The author proposed a solution that can suggest decisions related to purchasing new machines, selling old machines, and generating a series of non-dominated solutions. In another study, Dolgui et al. (2017) proposed a multi-objective BAP and employed aggregation methods to provide decision-makers with Pareto optimal solutions by considering throughput, buffer costs, and inventory cost. Two evolutionary algorithms, namely the simple evolutionary multi-objective optimizer (SEMO) and the Simple indicator-based evolutionary algorithm (SIBEA), were used. The proposed approach is applicable during the pre-construction phase of a manufacturing facility.

Recently, Chang and Huang (2022) proposed an efficient optimization framework for simultaneous machine selection and buffer allocation in HFSs that combined global exploration and local exploitation strategies by providing a problem-specific algorithm, namely, adaptive global particle and local hyperbox search (AGPLHS). The framework considers the selection of machine types and allocation of buffers and is applied to the resulting two-stage problem, in which machine types are first selected, and then buffers are allocated. Herps et al. (2022) studied the designing of an automated manufacturing system to optimize profit while simultaneously addressing the selection of automation concepts and buffer capacities. To address this problem, they employed a simulation-based approach.

A summary of the problem features and solution methods is presented in Table 1, highlighting the differences between the current study and those investigated in the reviewed literature. Firstly, this work simultaneously considers buffer allocation, server allocation (machine or operator), and material handling units as decision variables, tackling a significant gap in the literature where these variables are often considered in isolation. Secondly, while most reviewed studies employed single-objective optimization, this work adopts a comprehensive multi-objective approach by considering the maximization of throughput and minimization of buffers, operators, and AMRs simultaneously. Thirdly, this study incorporates crucial problem characteristics such as variant-dependent processing times, quality inspection (intermediate, feed-forward), and material supply limits, as elaborated in Section 3.1. These aspects are essential in modern manufacturing but are overlooked in the reviewed literature that relies on simplified assumptions. By including these characteristics, this study enhances the applicability of the model to real-world scenarios. Finally, the solution approach employed in this work sets it apart from the literature by not only employing simulation but also investigating the performance of five ML algorithms to develop a PM. By selecting the best-performing algorithm, this study advances the literature beyond the previously used response surface method (RSM) as PM, in addressing HFS-B&RA.

### 3. Case study explanation

This section provides an overview of the case study.

#### 3.1. Case study

This study addresses the stochastic real-world HFS-B&RA in a marine engine production line, where three types of marine engines are assembled. These engines are manufactured according to over 300 different specifications (product variants), demonstrating the complexity and scope of the application. Each product variant has a unique processing time at each workstation due to the distinct set of operations required for its production. A schematic illustration of the production line is provided in Fig. 1. Engine Types 1 and 2 are manufactured in production line A, whereas engine Type 3 is manufactured in production line B. Assembling an engine in production line A involves 5 assembly areas (with 32 stations that perform 46 assembly operations), painting, leakage tests (cold tests), repair, final performance tests (hot tests), customization, and packaging.

The production process is configured as an HFS with nine segments ($S_1$ to $S_9$), two of which are arranged in parallel. Except for the paint

shop, which is completely automated without human intervention, all the other operations are performed either manually or partially by humans. The production line leverages a cross-trained workforce, enabling flexible operator deployment across segments.

Human involvement introduces a significant amount of variability into the processes, which implies frequent reworks. When a cold or hot test (Fig. 1) is performed, the engine might not pass the test. The dispatching logic of engines in both cold and hot test areas involves feed-forward quality inspections. In the feed-forward quality inspection, engines that require repair are sent to a “rework area” with a capacity of three engines, one of which is only reserved for emergencies. After being repaired, the engines are returned directly to the following buffer, i.e., Buffer 4, which is located before the fully automated paint shop. Occasionally, different engines complete their repair processes simultaneously or at short intervals, owing to the random nature of the repair process. Therefore, the feed-forward quality inspection changes the entering sequence (sequence scrambling) and inert-arrival rate of the engines to the paint shop.

The paint shop is equipped with two robots. Meanwhile, the engines are transported through the paint shop via overhead cranes to Assembly 4. Six manual workstations are installed in Assembly 4, and the line design does not allow the placement of buffers. Consequently, the fully automated area (paint shop), which operates in deterministic times, is followed by highly stochastic manual operations. This change in the automation level results in a significant difference between the variations in the paint shop and the subsequent assembly stations.

Furthermore, at three points (i.e., the paint shop, cooling-tunnel buffer, and Assembly 5), Line A intersects with Line B. Because of uncertainty in the line caused by manual operations, rework, intersection with engines transferred from Line B, and variability in material supply, including potential shortages and inconsistent delivery times, the size of buffers, assignment of operators to each segment, and number of AMRs available in the line significantly affect the throughput of the system. Under such circumstances, the manufacturing company typically aims to maximize the throughput while minimizing the total number of buffers, operators, and AMRs.

### 4. Problem definition

This section presents the characteristics of the production system under study and formally defines the optimization problem.

Consider an HFS with mixed automation levels, comprising *S* segments, operating in a high-mix, low-volume production environment, as depicted in Fig. 2. A segment represents a distinct grouping of workstations that execute multiple operations in either serial or parallel configurations. The production line operates asynchronously, allowing for both blockage and starvation. To smooth the production flow, buffers are strategically utilized between the segments.

In a high-mix production environment, different product variants may require distinct sets of processing tasks. Let *V* denote the set of product variants to be processed, and *N* denote the set of tasks required to complete each product. The processing requirements are represented by an *N* × *V* matrix $t_{ij}$, where $t_{ij}$ represents the processing time distribution for task *i* on variant *j* if the task is required, and 0 otherwise. The processing times follow stochastic, variant-dependent distributions, reflecting real-world variability in task execution.

Each manual or semi-automated segment is operated by a pool of operators, as depicted by the operator pools in Fig. 2. Let *O* denote the set of available operators. Each operator *o* ∈ *O* is qualified to work in a specific subset of segments $Q_o \subseteq S$, indicated by the binary parameter $\alpha_{o,s}$, where $\alpha_{o,s} = 1$ if operator *o* is qualified to work in segment *s*, and 0 otherwise. The operators are characterized by multiple skill categories (processing, setup, and repair capabilities), segment-specific qualifications defined by $\alpha_{o,s}$, and dynamic allocation based on availability and skill requirements.

The segments are interconnected by buffers and serviced by AMRs, denoted by *R*, for inter-segment material handling. Material handling follows a First-In-First-Out buffer management policy.
<table><caption>**Table 1**
Comparison of problem features and objectives pertaining to HFS-B&RA.</caption><thead><tr><th rowspan="2">Publications</th><th colspan="3">Decision variable</th><th colspan="4">Objectives</th><th colspan="3">Problem characteristics</th><th colspan="2">Solution approach</th></tr><tr><th>Buffer capacities</th><th>Server allocation</th><th>Material handling</th><th>Max(Throughput)</th><th>Min(Buffers)</th><th>Min(Operators)</th><th>Min(AMRs)</th><th>Variant-dependent processing time</th><th>Quality inspection</th><th>Material supply limit</th><th>Simulation</th><th>Predictive model</th><th>Optimization method</th></tr></thead><tbody><tr><td>Dolgui et al. (2002)</td><td>✓</td><td></td><td></td><td>✓</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>GA</td></tr><tr><td>Abu Qudeiri et al. (2008)</td><td>✓</td><td></td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>GA</td></tr><tr><td>Colledani et al. (2008)</td><td>✓</td><td></td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>GA</td></tr><tr><td>Nahas et al. (2008)</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>ACO</td></tr><tr><td>Nahas et al. (2009)</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>ACO+SA</td></tr><tr><td>Massim et al. (2010)</td><td>✓</td><td></td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>AIA</td></tr><tr><td>Luo et al. (2010)</td><td></td><td>✓</td><td></td><td></td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td>✓</td><td></td><td>B&B</td></tr><tr><td>Amiri and Mohtashami (2012)</td><td>✓</td><td></td><td></td><td>✓</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td>✓</td><td>GA</td></tr><tr><td>Mohtashami (2014)</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td>✓</td><td></td><td>Hybrid GA</td></tr><tr><td>Dolgui et al. (2017)</td><td>✓</td><td></td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>SEMOO, SIBEA</td></tr><tr><td>Chang and Huang (2022)</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td>✓</td><td></td><td>AGPLHS</td></tr><tr><td>Herps et al. (2022)</td><td>✓</td><td></td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td>✓</td><td></td><td>SA</td></tr><tr><td>This work</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>En-NSGA-II</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Process diagram of investigated marine engine production line.</figcaption></figure>

## 4.1. Objectives

The optimization problem seeks to maximize the expected throughput denoted as $E[TH(B_s, Op_s, AMR_r)]$ while minimizing the total number of buffers $\sum_{s \in S} B_s$, operators $\sum_{s \in S} Op_s$, and AMRs $\sum_{r \in R} AMR_r$. The throughput function $TH$ is evaluated through simulation to account for the system's stochastic nature and production line dynamics.

## 4.2. Constraints

The problem is subject to capacity constraints that limit the total number of buffers $\sum_{s \in S} B_s \le B$ and operators $\sum_{s \in S} Op_s \le OP$. Each segment's buffer size and operator count are bounded by their respective limits: $B'_s \le B_s \le B'_s$ and $Op'_s \le Op_s \le Op'_s$ for all $s \in S$.

The operator assignment is governed by three conditions: (i) the number of operators assigned to each segment equals the required

number $\sum_{o \in O} X_{o,s} = Op_s$ for all $s \in S$, (ii) each operator is assigned to at most one segment $\sum_{s \in S} X_{o,s} \le 1$ for all $o \in O$, and (iii) operators are only assigned to segments for which they are qualified $X_{o,s} \le \alpha_{o,s}$ for all $o \in O, s \in S$. The decision variables are defined as: $B_s, Op_s \in \mathbb{Z}^+$ for all $s \in S, AMR_r \in \{0,1\}$ for all $r \in R$, and $X_{o,s} \in \{0,1\}$ for all $o \in O, s \in S$.

## 5. Solution approach

As depicted in Fig. 3, the solution approach provided to manage HFS-B&RA comprises four components: the simulation model and data generation, a PM, an optimizer, and a knowledge extraction procedure. The suggested solution approach can manage the complexity and stochasticity of real-world production lines owing to the adoption of DES. However, performing and optimizing the simulation models of real-world HFSs require significant computational time.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Illustration of a segmented HFS.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Solution approach.</figcaption></figure>

In this study, five PMs were constructed by training five regressor ML algorithms – ANN, GBT, PR, RFR, and SVM – based on data sampled from the simulation model. The best-performing PM was then selected to replace the time-consuming simulation trials, thereby significantly reducing the computational time of the optimization process while maintaining high accuracy. The PM was coupled with the proposed optimization algorithm (i.e., En-NSGA-II) to form the PMB-MOO. In Sections 6 through 8, the PMB-MOO approach (simulation model, PM construction and evaluation, and optimization algorithm) is elaborated in depth. The knowledge extraction process, performed to extract insights from the optimization results, is further detailed in Section 10.3.

## 6. Simulation model and data generation

The first component of the solution approach is a DES model. Modeling a manufacturing system via simulation allows one to capture the mechanisms and dynamics of the system (Escudero-Ornelas et al., 2023). Additionally, the results of changes in operations and management policies can be analyzed using a simulation model, and different scenarios can be compared without manipulating the actual system. In this study, a simulation model was used to determine the throughput of a production line using alternative designs. Alternative designs were generated using different combinations of buffers and resources. The operational details, including the stochastic engine arrival rate, setup times, processing times, failure rates, and other operational constraints, were captured in the simulation model. To construct a valid simulation model suitable for further experimentation, optimization, and sampling, the simulation steps introduced by Law (2022) were

employed and performed rigorously. In this study, the simulation model was intended to provide the required input-output dataset for training the PM. The validated simulation model, combined with the sampling method, was later used to generate a larger dataset. This dataset comprises the design characteristics as input values and the throughput obtained from the simulation model as the output. This dataset was later used to construct the PM.

### 6.1. Validation of the simulation model

The required data were obtained using three data-acquisition techniques: (1) direct observation, (2) interviews with experts, and (3) data extracted from company databases and documents. After the data were obtained, a simulation model was constructed using the C++ based software package FACTS Analyzer (Ng and Bernedixen, 2018). The choice of FACTS Analyzer<fnref n="1" /> was motivated by its established use within the industrial partner's operations, allowing for the integration of simulation results into their processes. Additionally, the software's accessibility to the research team – including source code availability and collaboration with the developers – offered the adaptability needed to model the complexities of the production system while meeting the specific requirements of this research. Finally, access to the executable version of FACTS Analyzer, known as xsim-runner, enabled smooth integration of customized optimization algorithms.

In Fig. 4, the values generated for the hourly throughput using the simulation model are compared against the real throughput of

<footnote n="1">¹ EVOMA, https://evoma.se/.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Simulation results versus real data for hourly throughput in a 50-day period (excluding weekends).</figcaption></figure>

the production line. This comparison demonstrates that the simulation model can effectively capture the real system fluctuations. However, the validity of the model was evaluated via a three-step statistical analysis procedure, using the Kolmogorov-Smirnov (K-S) test for normality, Bartlett test for equality of variances, and an independent t-test for equality of means (Jeong et al., 2023). The K-S test results showed that both the simulation model data (test statistic: 0.1257, p-value: 0.3769) and real data (test statistic: 0.0973, p-value: 0.6942) did not significantly deviate from a normal distribution. The Bartlett test (test statistic: 0.0894, p-value: 0.7650) indicated no significant difference in variances between the two datasets. Finally, the independent t-test (test statistic: -0.1938, p-value: 0.8467) suggested no significant difference in mean throughput between the simulation model and real system. These results indicate that the simulation model data reasonably represents the real system, as they share similar statistical properties (mean, variance, and distribution).

Next, steady-state analysis was performed to determine the warm-up period. The throughput (TH) and work-in-process (WIP) were used to determine the warm-up period. The resulting output of the simulation model showed minimum warm-up periods of 10 and 15 days for the WIP and TH, respectively. Consequently, the warm-up period was set to 15 days.

## 7. Predictive model generation and evaluation

The second component of the solution approach is responsible for building the PM. In this study, five common regressor ML algorithms in the manufacturing industries (do Amaral et al., 2022), namely ANNs, GBT, PR, RFR, and SVM were trained on data sampled from the simulation model.

The use of simulation as a data generation tool for developing PMs of the complex systems is well-established in addressing design problems in production, particularly when evaluating alternative designs in the absence of direct analytical solutions (Dieleman et al., 2023; Altparmak et al., 2007; Bulgak, 2006; Tsadiras et al., 2013). However, it is important to acknowledge that simulated data inherently deviates from real data, potentially introducing errors in the learning process of PMs. To minimize this error to the greatest extent possible, careful selection of output measures with the least deviation from real values and following rigorous and systematic steps in building PMs are crucial.

In this study, the average hourly throughput (TH) over the simulation horizon (three-month period) is utilized, resulting in negligible error compared to real data. This metric is calculated as the average of the hourly throughput values, shown in Fig. 4. For a more granular analysis, the average hourly throughput (TH) is calculated using Eq. (1):

$$TH = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{1}{r} \sum_{j=1}^{r} \frac{P_{ij}}{H_s \cdot S} \right) \quad (1)$$

where:

* n is the number of days in the simulation horizon

* r is the number of replications

* $P_{ij}$ is the total parts produced on day i in replication j

* $H_s$ is the working hours per shift

* S is the number of shifts per day

TH provides a robust basis for optimization by smoothing out short-term variability and capturing overall system performance. To demonstrate the congruence between the throughput (TH) derived from empirical data and that obtained from simulated data, a comparative analysis was conducted using the data visualized in Fig. 4. Over a 50-day simulation horizon, the TH calculated from the real-world data is 0.8870, which closely aligns with the simulation model's TH of 0.8789. This yields a minimal cumulative error of approximately 0.91%. This negligible discrepancy serves as a robust indicator of the simulation model's fidelity and reliability in long-term performance evaluation.

The low errors observed in the simulation model over extended periods contribute directly to the PM's reliability. Given the small error margin between the simulation results and real-world data, the PMs built upon the simulated data perform effectively in generating reliable approximations of the TH (see Section 7.4). These models were trained to perform well over time horizons of three months, aligning with the operational needs of the company.

Another important factor in building reliable PMs in this study is the employment of a rigorous methodological approach. This approach ensures that the PMs are developed with a strong foundation in statistical and machine learning principles, further enhancing their accuracy and reliability. The rigorous methodological approach employed in building reliable PMs encompassed several key aspects. Diverse datasets were generated to represent a wide parameter space, ensuring the models could generalize well across different scenarios. Feature engineering and standardization techniques were applied to enhance model performance (Duboue, 2020). Multiple candidate models were trained and extensively evaluated to select the best-performing model. Through cross-validation and hyperparameter tuning, the selected PM achieved a close approximation of real-world throughput values. This comprehensive approach enabled the PMs to effectively capture the system's behavior, providing reliable inputs for the optimization algorithms. The integration of accurate simulation data and robust predictive modeling ensured that the overall solution approach delivered precise and dependable results. By combining these methodological elements, the developed PMs demonstrated high reliability and accuracy in representing the complex production system, thus forming a solid foundation for subsequent optimization processes.

To construct PMs, four steps, i.e., data sampling and preprocessing, training MLs, hyperparameter tuning, and performance evaluation, were performed as suggested in the literature (Malakouti et al., 2023; Norouzi et al., 2023) and are elaborate in Sections 7.1 to 7.4.
## 7.1. Data sampling and preprocessing

As the first step, input data were sampled from the simulation model. Because the problem described (Section 4) results in a complicated combinatorial optimization problem, the number of samples required to encompass the entire decision space must be sufficiently large. In this study, the automated hypercube sampling (Rajabi et al., 2015) method was used to generate 25 000 samples using the simulation model. The samples were generated at the minimum lower bound and the maximum upper bound. Hence, the PM remains valid for all inputs within the ranges provided in Table 3. Then, the dataset generated via hypercube sampling was preprocessed for duplicates and normalized for use in the subsequent learning step.

## 7.2. Model training

To train and test the performance of the PMs, the data was split into two subsets. 75% of the data was allocated for training the models, allowing them to learn patterns and relationships from a large portion of the available data. The remaining 25% of the data was reserved for testing and evaluating the trained models' performance on unseen data. Thus, for all the algorithms, the size of training and testing sets are 18 750 and 6250, respectively.

## 7.3. Hyperparameter tuning

Grid search with 10-fold cross-validation was utilized for tuning the hyperparameters of the PMs (Belete and Huchaiah, 2022). Grid search with 10-fold cross-validation is a process where, for each combination of hyperparameters, the model is trained and evaluated 10 times. The performance of each combination is averaged over the 10 folds. The hyperparameters that yield the best average performance are selected as the optimal set.

## 7.4. Performance evaluation of MLs

The performance of ML algorithms was evaluated using five performance measures, i.e., R-squared ($R^2$), mean absolute error (MAE), mean squared error (MSE), root mean squared error (RMSE), and mean absolute percentage error (MAPE), presented in Eqs. (2) to (6), respectively. In Eqs. (2) to (6), $y_i$ is the value generated by the simulation model, $\hat{y}_i$ is the value predicted by the RFR, $\bar{y}_i$ is the mean value of $y_i$, and $n$ is the sample size.

$$R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y}_i)^2} \qquad (2)$$

$$MAE = \frac{1}{n} \sum_{j=1}^{n} |y_j - \hat{y}_j| \qquad (3)$$

$$MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \qquad (4)$$

$$RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \qquad (5)$$

$$MAPE = \frac{1}{n} \sum_{i=1}^{n} \frac{|y_i - \hat{y}_i|}{|y_i|} \qquad (6)$$

The results of the evaluation process are reported in Table 2. The identical measures were computed for both the train and test datasets. Table 2 shows that Based on the performance evaluation of the machine learning algorithms for regression, RFR consistently outperforms the other algorithms across most measures. It achieves the highest R-squared value of 92.19%, indicating a strong fit to the data, and the lowest MAE of 0.27 and MAPE of 18.01%, suggesting minimal prediction errors. SVM also shows promising results, particularly in terms of MSE and RMSE, where it has the lowest values of 5.26 and

2.29, respectively. The GBT and Polynomial PR algorithms perform moderately well, while the ANN generally has the lowest performance across the measures.

Among the performance measures used to assess MLs in Table 2, the R-squared and MAPE are particularly crucial (Chicco et al., 2021; Prayudani et al., 2019). R-squared measures the model's fit to observed data, while MAPE quantifies the average error as a percentage. Experts recommend an R-squared above 50% and a MAPE below 20% for reliable models (Golmohammadi, 2015). Among the regressor ML algorithms evaluated in Table 2, RFR excelled with the highest R-squared and lowest MAPE, proving to be the most suitable choice for building PM.

The dominance of RFR is explainable by the fact that RFR is an appropriate approach for tackling HFS-B&RA due to its balance between predictive power and model interpretability, its robustness to outliers, its capability to handle missing values, and its resistance to overfitting due to its ensemble nature. Additionally, it has the ability to manage high-dimensional data and is effective in dealing with discrete decision spaces (El Mrabet et al., 2022). These considerations collectively establish the RFR as an excellent and appropriate choice for this study.

# 8. Optimization

In the proposed solution approach, the optimization component determines the optimal combination of decision variables (i.e., B&RA). For large search spaces, random search and metaheuristics are the preferred optimization approaches (Fares et al., 2022). An enhanced version of the NSGA-II, namely En-NSGA-II, is proposed in this study for optimization purposes. The pseudocode for the proposed En-NSGA-II, presented in Algorithm 1, features a customized solution representation, a mechanism for performing VNS, and an LLML as will be presented in Sections 8.1 to 8.3.

The notations used in Algorithm 1 are as follows: $N$ represents the population size; $G$ is the number of generations; $P_n$, $Q_n$, and $R_n$ are the set of parents, offspring, and combined population in the $n_{th}$ generation, respectively; and $P_c$, $\eta_c$, $P_m$, and $\eta_m$ are the crossover probability, distribution index of the crossover, mutation probability, and mutation parameter, respectively. $P_{vns}$ denotes the proportion of the population undergoing the VNS mechanism. $P_{Inc-Dec}$, $P_{FPR}$, and $P_{TSR}$ represent the probabilities of the three local search mechanisms within VNS. `LLML`_`Params` is a set of parameters for LLML, which will be introduced in Section 8.3.

## 8.1. Solution representation

Metaheuristic algorithms perform most effectively when the solution is effectively represented (Nourmohammadi et al., 2022). In this study, a solution representation of the defined problem was designed with a length of 2S+1. In this representation, the first S cells are related to the size of the buffers following each segment, whereas the second S cells refer to the number of operators assigned to each segment. The last cell determines the number of available AMRs and is shared between the segments. The values in the representation are integers. Fig. 5 illustrates the proposed solution representation.

## 8.2. Crossover and mutation

En-NSGA-II utilizes crossover and mutation operators to maintain diversity when generating a new population. This section presents a brief introduction and details regarding the pseudocode implementation for the crossover and mutation operators.

The applied crossover generates offspring by blending the parent solutions ($x_p^1$ and $x_p^2$) while maintaining the balance between exploration and exploitation. Additionally, it employs a distribution index ($\eta_c$) and
<table><caption>**Table 2**
Performance evaluation of PMs (regressor ML algorithms).</caption><thead><tr><th rowspan="2">Measure</th><th colspan="2">ANN</th><th colspan="2">GBT</th><th colspan="2">PR</th><th colspan="2">RFR</th><th colspan="2">SVM</th></tr><tr><th>Train</th><th>Test</th><th>Train</th><th>Test</th><th>Train</th><th>Test</th><th>Train</th><th>Test</th><th>Train</th><th>Test</th></tr></thead><tbody><tr><td>R<sup>2</sup> (%)</td><td>79.19</td><td>79.84</td><td>82.50</td><td>82.71</td><td>81.21</td><td>82.51</td><td>91.26</td><td><strong>92.19</strong></td><td>78.98</td><td>81.19</td></tr><tr><td>MAE</td><td>1.77</td><td>1.51</td><td>2.12</td><td>1.93</td><td>4.28</td><td>4.01</td><td>0.27</td><td><strong>0.27</strong></td><td>0.78</td><td>0.72</td></tr><tr><td>MSE</td><td>9.14</td><td>8.81</td><td>15.23</td><td>14.89</td><td>31.80</td><td>28.95</td><td>7.49</td><td>5.92</td><td>6.68</td><td><strong>5.26</strong></td></tr><tr><td>RMSE</td><td>4.59</td><td>4.45</td><td>3.90</td><td>3.85</td><td>5.64</td><td>5.37</td><td>2.74</td><td>2.43</td><td>2.58</td><td><strong>2.29</strong></td></tr><tr><td>MAPE (%)</td><td>29.61</td><td>25.31</td><td>23.82</td><td>23.01</td><td>28.23</td><td>27.01</td><td>18.03</td><td><strong>18.01</strong></td><td>36.05</td><td>35.09</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** The proposed solution representation.</figcaption></figure>

**Algorithm 1 En-NSGA-II Algorithm**

1: **Input:** N, G, $P_C$, $P_m$, $\eta_c$, $\eta_m$, $P_{vns}$, $P_{Inc-Dec}$, $P_{TPR}$, $P_{TSR}$, $lr_timeout$, **LLM L Params**
2: Generate population $P_0$ with N individuals
3: $P_0 \leftarrow Normalize(P_0)$
4: Evaluate each individual in $P_0$ using predictive model RFR($P_0$)
5: Perform non-dominated sorting on $P_n$, and assign rank $r_i$ to each $i$
6: Calculate crowding distance ($d_i$) for each $i$ in $P_n$
7: $n \leftarrow 0$
8: **while** $n < G$ **do**
9:   **while** $lr_time < lr_timeout$ **do**
10:    Perform LLM L($P_n$, LLM L Params)
11:    Update $P_C$, $P_m$, $P_{vns}$, $P_{Inc-Dec}$, $P_{TPR}$, $P_{TSR}$
12:    **end while**
13:    Select parents using binary tournament selection based on $r_i$ and $d_i$
14:    Perform crossover with $P_C$ and $\eta_c$, generate $Q_n$
15:    Perform mutation on $Q_n$ with $P_m$ and $\eta_m$
16:    Apply VNS to $P_{vns}$ portion of $Q_n$
17:    Evaluate $TH$, $\sum_{i} b_i$, $\sum_{i} Op_i$, and AMR for each individual in $Q_n$
18:    Combine $P_n$ and $Q_n$ to form $R_n$
19:    Perform non-dominated sorting on $R_n$, assign $r_i$ to each $i$
20:    Calculate $d_i$ for each $i$ in $R_n$
21:    Select top N individuals from $R_n$ based on $r_i$ and $d_i$, and form $P_{n+1}$
22:    $n \leftarrow n + 1$
23:    **end while**
24:    **return** non-dominated solutions from $P_G$ as Pareto-optimal solutions.

crossover probability ($P_C$) to control the offspring diversity, thereby yielding solutions similar to their parents yet sufficiently diverse for exploring the search space effectively. As shown in Algorithm 2, a random number is generated for each gene in the parent-solution representation to determine whether a crossover occurred. If the crossover probability is satisfied and the parent values differ, then the offspring values are calculated using a $\beta$ value based on $\eta_c$, a random number, and the difference between the parent values. If the parent values are the same or no crossover occurs, then the offspring values are set to be equal to the parent values. Finally, the algorithm returns two offspring solutions.

En-NSGA-II employs a mutation operator to avoid convergence to the local optima. The details of the applied mutation are presented in Algorithm 3, where $x_p$ is the parent solution; $x_c$ is the child solution; $x_l$ and $x_u$ are the lower and upper bounds of the solution variable, respectively, and $r$ is a random number uniformly distributed in the range [0, 1]. $P_m$ and $\eta_m$ are the mutation probability and parameter, respectively. As shown in Algorithm 3, a random value is generated in

the range [0, 1] and compared with $P_m$. If $r \le P_m$, then $x_p$ is sampled to the left-hand side, which is the region between $x_l$ and $x_p$; otherwise, it is sampled to the right-hand side, which is the region between $x_p$ and $x_u$. Based on the calculated $\delta_q$, a polynomial mutation can mutate $x_p$, even where $x_p$ is situated in close proximity to either of its boundaries. Finally, the value for the child ($x_c$) is calculated.

**Algorithm 2 Crossover operator**

1: Initialize parameters: $P_C$, $\eta_c$, $x_p^1$, $x_p^2$
2: Create empty arrays for $x_c^1$, $x_c^2$ with the same size as parents
3: **for** each gene $i$ in the representation of $x_p^1$, $x_p^2$ **do**
4:     Generate a random value $r$ in the range [0, 1]
5:     **if** $r \le P_C$ **then**
6:         $\delta = |x_p^1(i) - x_p^2(i)|$
7:         **if** $\delta > 0$ **then**
8:             Let $y_1 = \min(x_p^1(i), x_p^2(i))$, $y_2 = \max(x_p^1(i), x_p^2(i))$
9:             Generate a random value $u$ in the range [0, 1]
10:         **if** $u \le 0.5$ **then**
11:             $\beta = (2u)^{\frac{1}{n_c+1}}$
12:         **else**
13:             $\beta = (\frac{1}{2(1-u)})^{\frac{1}{n_c+1}}$
14:         **end if**
15:     $x_c^1(i) = 0.5((1+\beta)y_1 + (1-\beta)y_2)$
16:     $x_c^2(i) = 0.5((1-\beta)y_1 + (1+\beta)y_2)$
17:     **else**
18:         $x_c^1(i) = x_p^1(i)$
19:         $x_c^2(i) = x_p^2(i)$
20:     **end if**
21:     **else**
22:         No crossover for this gene
23:     **end if**
24:     Return $x_c^1(i)$ and $x_c^2(i)$
25: **end for**

## 8.3. Enhancement mechanism

In this research, two enhancement mechanisms were utilized to improve the NSGA-II algorithm's performance. The first mechanism employed a VNS structure consisting of three local search strategies, namely, Flexible Increment-Decrement (Flex-Inc-Dec) (Fig. 6(a)),
<figure><img src="image_8.png" /><figcaption>Fig. 6. Local search strategies.</figcaption></figure>

**Algorithm 3 Mutation operator**

1: Initialize parameters: $P_n, \eta_m$

2: Generate a random value $r$ in the range [0, 1]

3: Calculate $\delta_1 = (x_p - x_l) / (x_u - x_l)$ and $\delta_2 = (x_u - x_p) / (x_u - x_l)$

4: if ($r \le P_m$) then

5:     $\delta_q = [2r + (1 - 2r)(1 - \delta_1)^{\eta_m + 1}]^{\frac{1}{\eta_m + 1}} - 1$

6: else

7:     $\delta_q = 1 - [2(1 - r) + 2(r - 0.5)(1 - \delta_2)^{\eta_m + 1}]^{\frac{1}{\eta_m + 1}}$

8: end if

9: $x_c = x_p + \delta_q (x_u - x_l)$

Triple Point Rearrangement (TPR) (Fig. 6(b)), and Triple Section Rearrangement (TSR) (Fig. 6(c)).

As depicted in Fig. 6, the local search strategies were separately applied to the buffer and operator sections of the chromosome to ensure that the total number of buffers and operators remained constant throughout the search process. For a certain number of iterations, the Flex-Inc-Dec strategy randomly selects two locations and increases the size of the first location by a randomly generated value while decreasing the size of the second location by the same amount. TPR involves randomly selecting three adjacent genes at any position and exploring all possible permutations of these genes to find the most optimal solution. On the other hand, TSR selects two cut points on the chromosome to divide it into three segments and then permutes

these segments to generate five potential new permutations. Each permutation is evaluated to determine the best solution.

The second enhancement mechanism is an LLML strategy (Algorithm 4), which adaptively adjusts several key parameters of the algorithm, including the crossover probability, mutation probability, the overall probability of employing VNS, and the individual probabilities of selecting each local search strategy within the VNS framework. This LLML approach enables the algorithm to dynamically adapt to the problem landscape over time, optimizing its performance across various stages of the search process. The LLML algorithm starts with a set of input parameters: $P_n$ (the current population), $P_c$ (crossover probability), $P_m$ (mutation probability), $P_{vns}$ (VNS probability), $P_{Inc-Dec}$ (Flex-Inc-Dec probability), $P_{TPR}$ (TPR probability), $P_{TSR}$ (TSR probability), $lr_rate$ (learning rate), and $lr_iters$ (number of iterations of LLML loop).

Initially, the algorithm evaluates the fitness of each individual in the population $P_n$ and calculates the diversity of the population using crowding distance, a key component of NSGA-II. Crowding distance measures the density of solutions around each individual in the Pareto front, with higher values indicating better diversity. The crowding distance $d_l$ for an individual $i$ is calculated as $d_l = \sum_{k=1}^{m} \frac{f_{k,i+1} - f_{k,i-1}}{f_{k,max} - f_{k,min}}$ where $m$ is the number of objectives, $f_{k,i+1}$ and $f_{k,i-1}$ are the adjacent objective values, and $f_{k,max}$, $f_{k,min}$ are the maximum and minimum values for objective $k$. Low crowding distance indicates clustering, signaling the need for adaptive changes to parameters like crossover, mutation, or population size.

It then checks if there is a match in the memory (storage of past parameter sets and their outcomes) based on the diversity of the current
population. If a match is found, the algorithm updates the parameters $P_{vns}$, $P_{Inc-Dec}$, $P_{TPS}$, $P_{TSR}$, $P_c$, and $P_m$ with the values from memory, bypassing the need for further learning iterations in this step. If no match is found, the algorithm constructs an orthogonal array O from the parameters, and through a series of crossover, mutation, and VNS operations, it evaluates the standard deviation of the best fitness values obtained for each set of parameters in O. The set of parameters that yields the maximum standard deviation is considered the best, and its values are used to update the original parameters using a learning rate $lr$. The learning rate, $lr$, balances the exploration of new parameter spaces and exploitation of known good parameter sets, determining the extent to which new parameters are adopted. The updated parameters are stored in memory along with their corresponding diversity and standard deviation, allowing the algorithm to reference past successful parameter sets in future iterations. This process continually refines the parameters for optimizing the population's fitness over time, adapting to new challenges and environments in a lifelong learning context.

### Algorithm 4 Lifelong Meta-Learning Algorithm

1: **Input:** $P_n$, $P_s$, $P_c$, $P_m$, $P_{vns}$, $P_{Inc-Dec}$, $P_{TPR}$, $P_{TSR}$, $lr\_rate$, $lr\_iters$
2: fitness\_values ← fitness values for each individual in $P_n$
3: diversity ← Diversity($P_n$)
4: memory\_match ← False
5: if memory \* $\neq$ 0 then
6: for each (member, parameters) $\in$ memory do
7: if diversity $\approx$ diversity\_member then
8: $P_{vns}$, $P_{Inc-Dec}$, $P_{TPS}$, $P_{TSR}$, $P_c$, $P_m$ ← parameters
9: memory\_match ← True
10: break
11: end if
12: end for
13: end if
14: if not memory\_match then
15: O ← orthogonal array($P_{vns}$, $P_{Inc-Dec}$, $P_{TPR}$, $P_{TSR}$, $P_c$, $P_m$)
16: std\_devs ← []
17: for each $o \in O$ do
18: $P_{temp}$ ← $P_n$
19: for $i = 1$ to $lr\_iters$ do
20: $P_{temp}$ ← Crossover($P_{temp}$, $o$)
21: $P_{temp}$ ← Mutation($P_{temp}$, $o$)
22: $P_{temp}$ ← VNS($P_{temp}$, $o$)
23: Best\_fitness ← Append best value in $P_{temp}$ for each objective
24: end for
25: [std\_devs, $o$] ← Max(std\_Best\_fitness) and corresponding parameters $o$
26: end for
27: $\delta$ ← best\_orthogonal(std\_devs, O)
28: Update parameters using learning rate ($lr$):
29: $P_{vns}$ ← $lr \cdot \delta[1] + (1 - lr) \cdot P_{vns}$
30: $P_{Inc-Dec}$ ← $lr \cdot \delta[2] + (1 - lr) \cdot P_{Inc-Dec}$
31: $P_{TPR}$ ← $lr \cdot \delta[3] + (1 - lr) \cdot P_{TPR}$
32: $P_{TSR}$ ← $lr \cdot \delta[4] + (1 - lr) \cdot P_{TSR}$
33: $P_c$ ← $lr \cdot \delta[5] + (1 - lr) \cdot P_c$
34: $P_m$ ← $lr \cdot \delta[6] + (1 - lr) \cdot P_m$
35: memory ← Append the new parameters and their corresponding diversity and standard deviation
36: end if

## 9. Computational experiments

This segment presents the results of the computational experiments conducted to test the proposed PM-based MOO for addressing HFS-B&RA.

Considering that the PM model has already been validated, the performance of the solution approach heavily relies on the optimization algorithm efficiency, in this case En-NSGA-II. Initially, the study introduces five widely recognized performance metrics to evaluate the

algorithm's performance. Owing to the absence of established benchmark test instances for HFS-B&RA, the study generated test instances randomly, according to the characteristics of the industrial case. Subsequently, a series of parameter tuning experiments were performed to ascertain the optimal settings for the algorithms. The performance of En-NSGA-II was then benchmarked against the conventional NSGA-II and MOPSO algorithms.

NSGA-II was selected for enhancement due to its widespread use in solving multi-objective optimization problems, attributed to its elitist strategy and ability to maintain a diverse set of solutions (Neumann et al., 2024). The enhanced version, En-NSGA-II, aimed to further improve its performance in handling the complexities of the HFS-B&RA problem. On the other hand, MOPSO was chosen for comparison owing to its effectiveness in similar optimization problems (Shao et al., 2020; Chang and Huang, 2022). Thus, the selection of these two algorithms is justified by their popularity and broad use across a range of problems (Zhou et al., 2011).

### 9.1. Experimental settings

Due to the lack of test problems for HFS-B&RA in the literature, test instances are randomly generated. To maintain the essence of the real-world problem, the instances were generated following the case study observations and characteristics. The problem addressed in this study includes three sets of decision variables: the buffer sizes, the number of operators allocated to each segment, and the number of AMRs in the production line. The variable intervals are determined by engineers and system experts based on the availability of physical space, availability of resources, and compatibility with strategic and future expansion plans. The details of the generated instances can be found in Table 3. All algorithms were coded in Python, and the experiments were performed on a computer with a 12th Gen Intel(R) Core (TM) i9 operating at 2400 MHz and 128 GB of RAM. Several libraries were leveraged for data processing, visualization, ML operations, and optimization. Key libraries included NumPy (v1.24.3) and Pandas (v1.5.3) for numerical operations and data manipulation, while Matplotlib (v3.7.1) and Seaborn (v0.12.2) facilitated statistical visualizations. Machine learning tasks were implemented using Scikit-learn (v1.2.2), with models such as RFR, ANN, GBT, PR, and SVM. Interactive visualizations were created using Plotly (v5.13.1), and Joblib (v1.3.2) was used for efficient model serialization, ensuring scalability and reproducibility. For the optimization process, additional libraries and custom tools were utilized. These included random (Standard Library) for randomness, time (Standard Library) for timing operations, and multiprocessing (Standard Library) for parallelization. All functions essential for optimization were custom-coded in Python (v3.11.5), ensuring complete control and adaptability for the problem requirements.

All algorithms were compared based on five metrics: coverage metric (CM), hypervolume (HV), number of non-dominated solutions (NNDS), spacing metric (SM), and diversity metric (DM) (Li et al., 2023b; Rabet et al., 2024). Together, these metrics provide a robust framework for evaluating the performance of the En-NSGA-II algorithm in solving the HFS-B&RA problem, offering a comprehensive view of its effectiveness in terms of convergence, variety, distribution, and exploration of the solution space. Coverage metric (CM) is a measure used to evaluate the relative performance of two sets of solutions, denoted as P and Q, which are generated by different optimization algorithms. This metric assigns a value within the range of [0, 1] to the ordered pair (P, Q), as defined in Eq. (7):

$$CM(P, Q) = \frac{|\{q \in Q | \exists p \in A : p > q \text{ or } p = q\}|}{|Q|} \qquad (7)$$

In Eq. (7), CM(P, Q), quantifies the extent to which solutions in set Q are dominated by or equal to solutions in set P. If every solution in Q is dominated by at least one solution in P, then CM(P, Q) = 1. Conversely, if no solution in Q is dominated by any solution in P, then
<table><caption>**Table 3**
The experimental settings in the case study and test problems.</caption><thead><tr><th></th><th>Instance</th><th>B1</th><th>B2</th><th>B3</th><th>B4</th><th>B5</th><th>B6</th><th>B7</th><th>B8</th><th>B9</th><th>Op1</th><th>Op2</th><th>Op3</th><th>Op4</th><th>Op5</th><th>Op6</th><th>Op7</th><th>AMR</th><th>Size</th></tr></thead><tbody><tr><td rowspan="12">Small</td><td>1</td><td>3</td><td>9</td><td>2</td><td>8</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>5</td><td>2</td><td>8</td><td>4</td><td>3</td><td>3</td><td>3.76E+04</td></tr><tr><td>2</td><td>3</td><td>2</td><td>2</td><td>2</td><td>7</td><td>2</td><td>3</td><td>9</td><td>2</td><td>5</td><td>2</td><td>2</td><td>9</td><td>2</td><td>2</td><td>3</td><td>6</td><td>6.14E+04</td></tr><tr><td>3</td><td>3</td><td>7</td><td>2</td><td>3</td><td>2</td><td>2</td><td>2</td><td>5</td><td>2</td><td>5</td><td>6</td><td>4</td><td>2</td><td>4</td><td>2</td><td>2</td><td>6</td><td>8.64E+04</td></tr><tr><td>4</td><td>2</td><td>3</td><td>3</td><td>7</td><td>3</td><td>2</td><td>6</td><td>2</td><td>5</td><td>4</td><td>2</td><td>3</td><td>4</td><td>4</td><td>2</td><td>3</td><td>2</td><td>1.04E+05</td></tr><tr><td>5</td><td>2</td><td>9</td><td>4</td><td>7</td><td>10</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>3</td><td>5</td><td>10</td><td>2</td><td>2</td><td>2</td><td>3</td><td>1.87E+05</td></tr><tr><td>6</td><td>4</td><td>2</td><td>2</td><td>10</td><td>10</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>9</td><td>2</td><td>6</td><td>7</td><td>2</td><td>3</td><td>3</td><td>2.33E+05</td></tr><tr><td>7</td><td>2</td><td>7</td><td>2</td><td>4</td><td>3</td><td>2</td><td>2</td><td>4</td><td>5</td><td>2</td><td>2</td><td>10</td><td>2</td><td>10</td><td>8</td><td>2</td><td>2</td><td>2.45E+05</td></tr><tr><td>8</td><td>2</td><td>10</td><td>5</td><td>4</td><td>2</td><td>3</td><td>2</td><td>2</td><td>3</td><td>6</td><td>2</td><td>2</td><td>2</td><td>9</td><td>6</td><td>6</td><td>2</td><td>4.32E+05</td></tr><tr><td>10</td><td>3</td><td>7</td><td>2</td><td>4</td><td>8</td><td>3</td><td>2</td><td>2</td><td>2</td><td>6</td><td>3</td><td>4</td><td>5</td><td>4</td><td>3</td><td>3</td><td>2</td><td>7.26E+05</td></tr><tr><td>9</td><td>2</td><td>2</td><td>2</td><td>10</td><td>2</td><td>3</td><td>3</td><td>2</td><td>4</td><td>9</td><td>2</td><td>5</td><td>5</td><td>5</td><td>5</td><td>2</td><td>5</td><td>8.85E+05</td></tr><tr><td rowspan="12">Medium</td><td>11</td><td>8</td><td>5</td><td>2</td><td>6</td><td>2</td><td>6</td><td>7</td><td>4</td><td>2</td><td>5</td><td>2</td><td>2</td><td>4</td><td>9</td><td>2</td><td>2</td><td>3</td><td>2.42E+06</td></tr><tr><td>12</td><td>8</td><td>3</td><td>4</td><td>7</td><td>9</td><td>2</td><td>3</td><td>8</td><td>2</td><td>2</td><td>2</td><td>9</td><td>4</td><td>2</td><td>4</td><td>2</td><td>5</td><td>8.13E+06</td></tr><tr><td>13</td><td>7</td><td>2</td><td>2</td><td>3</td><td>9</td><td>2</td><td>6</td><td>5</td><td>4</td><td>3</td><td>2</td><td>2</td><td>7</td><td>10</td><td>5</td><td>4</td><td>2</td><td>2.24E+07</td></tr><tr><td>14</td><td>5</td><td>5</td><td>6</td><td>2</td><td>3</td><td>2</td><td>4</td><td>9</td><td>3</td><td>3</td><td>2</td><td>2</td><td>7</td><td>7</td><td>10</td><td>6</td><td>2</td><td>2.49E+07</td></tr><tr><td>15</td><td>2</td><td>10</td><td>9</td><td>8</td><td>2</td><td>6</td><td>6</td><td>9</td><td>3</td><td>2</td><td>6</td><td>2</td><td>10</td><td>2</td><td>6</td><td>2</td><td>3</td><td>9.07E+07</td></tr><tr><td>16</td><td>3</td><td>10</td><td>4</td><td>6</td><td>6</td><td>8</td><td>3</td><td>6</td><td>7</td><td>2</td><td>2</td><td>6</td><td>4</td><td>2</td><td>8</td><td>2</td><td>3</td><td>1.19E+08</td></tr><tr><td>17</td><td>3</td><td>10</td><td>4</td><td>3</td><td>2</td><td>4</td><td>7</td><td>13</td><td>2</td><td>2</td><td>6</td><td>9</td><td>2</td><td>9</td><td>4</td><td>5</td><td>4</td><td>3.58E+08</td></tr><tr><td>18</td><td>8</td><td>8</td><td>2</td><td>5</td><td>7</td><td>8</td><td>4</td><td>2</td><td>11</td><td>6</td><td>8</td><td>9</td><td>2</td><td>9</td><td>4</td><td>7</td><td>2</td><td>1.66E+09</td></tr><tr><td>19</td><td>5</td><td>2</td><td>2</td><td>7</td><td>9</td><td>7</td><td>7</td><td>3</td><td>2</td><td>8</td><td>6</td><td>8</td><td>9</td><td>10</td><td>3</td><td>2</td><td>6</td><td>2.44E+09</td></tr><tr><td>20<sup>a</sup></td><td>2</td><td>3</td><td>4</td><td>10</td><td>6</td><td>5</td><td>4</td><td>4</td><td>3</td><td>6</td><td>7</td><td>8</td><td>6</td><td>8</td><td>4</td><td>5</td><td>6</td><td>8.57E+09</td></tr><tr><td rowspan="12">Large</td><td>21</td><td>7</td><td>10</td><td>7</td><td>3</td><td>8</td><td>6</td><td>8</td><td>9</td><td>5</td><td>2</td><td>3</td><td>3</td><td>7</td><td>3</td><td>10</td><td>4</td><td>6</td><td>3.29E+10</td></tr><tr><td>22</td><td>8</td><td>2</td><td>3</td><td>6</td><td>9</td><td>8</td><td>6</td><td>11</td><td>8</td><td>10</td><td>9</td><td>7</td><td>2</td><td>2</td><td>10</td><td>6</td><td>3</td><td>5.33E+10</td></tr><tr><td>23</td><td>4</td><td>10</td><td>9</td><td>2</td><td>6</td><td>8</td><td>8</td><td>9</td><td>4</td><td>10</td><td>5</td><td>7</td><td>2</td><td>8</td><td>8</td><td>4</td><td>10</td><td>1.81E+11</td></tr><tr><td>24</td><td>5</td><td>10</td><td>10</td><td>2</td><td>3</td><td>7</td><td>4</td><td>14</td><td>5</td><td>6</td><td>7</td><td>3</td><td>2</td><td>8</td><td>8</td><td>6</td><td>4</td><td>2.85E+11</td></tr><tr><td>25</td><td>4</td><td>6</td><td>10</td><td>4</td><td>4</td><td>4</td><td>8</td><td>10</td><td>5</td><td>5</td><td>9</td><td>2</td><td>10</td><td>7</td><td>9</td><td>8</td><td>5</td><td>3.56E+11</td></tr><tr><td>26</td><td>8</td><td>8</td><td>9</td><td>8</td><td>7</td><td>5</td><td>5</td><td>8</td><td>7</td><td>5</td><td>9</td><td>9</td><td>10</td><td>2</td><td>2</td><td>7</td><td>4</td><td>4.59E+11</td></tr><tr><td>27</td><td>6</td><td>10</td><td>8</td><td>4</td><td>6</td><td>7</td><td>5</td><td>13</td><td>8</td><td>8</td><td>7</td><td>7</td><td>9</td><td>9</td><td>9</td><td>2</td><td>2</td><td>1.23E+12</td></tr><tr><td>28</td><td>5</td><td>10</td><td>9</td><td>10</td><td>5</td><td>8</td><td>8</td><td>14</td><td>7</td><td>3</td><td>10</td><td>10</td><td>6</td><td>6</td><td>7</td><td>9</td><td>4</td><td>2.57E+12</td></tr><tr><td>29</td><td>8</td><td>9</td><td>9</td><td>6</td><td>7</td><td>7</td><td>8</td><td>15</td><td>7</td><td>9</td><td>7</td><td>4</td><td>7</td><td>5</td><td>3</td><td>9</td><td>3</td><td>5.24E+12</td></tr><tr><td>30</td><td>8</td><td>6</td><td>9</td><td>5</td><td>6</td><td>6</td><td>5</td><td>11</td><td>8</td><td>4</td><td>6</td><td>5</td><td>5</td><td>10</td><td>10</td><td>10</td><td>5</td><td>5.49E+12</td></tr></tbody></table>

CM(P, Q) = 0. It is important to note that CM(P, Q) is not necessarily equal 1 – Coverage(Q, P), as there may exist solutions in P and Q that do not dominate each other.

Hypervolume (HV), defined in Eq. (8), measures the volume covered by the non-dominated solutions in the objective space.

$$HV(S) = \text{Volume}(\bigcup_{x \in S} [f(x), r]) \qquad (8)$$

where S is the set of non-dominated solutions, f(x) is the objective vector of a solution x, and r is a reference point.

Number of Non-dominated Solutions (NNDS), presented in Eq. (9), counts the solutions that are not dominated by any other solution.

$$NNDS(S) = |\{x \in S | \bar{y} \in S : y \text{ dominates } x\}| \qquad (9)$$

Spacing Metric (SM), defined in Eq. (10), assesses the distribution of solutions along the Pareto front.

$$SM(S) = \sqrt{\frac{1}{|S| - 1} \sum_{i=1}^{|S|} (\bar{d}_i - d_i)^2} \qquad (10)$$

where $d_i$ is the Euclidean distance between consecutive solutions, and $\bar{d}$ is the average of all $d_i$.

Diversity Metric (DM), defined in Eq. (11), Evaluates the extent of variation among the solutions.

$$DM(S) = \frac{\sum_{i=1}^{|S|-1} |d_i - \bar{d}|}{\max_{i,j} d(x_i, x_j)} \qquad (11)$$

where $d(x_i, x_j)$ is the Euclidean distance between two solutions $x_i$ and $x_j$, and $\max_{i,j}$ is the maximum distance between any two solutions in the set.

## 9.2. Parameter tuning

Parameter tuning is widely recognized as a critical factor influencing the performance of algorithms (Nourmohammadi et al., 2023; Mahmoodi et al., 2024). In this research, the gray-based Taguchi

method (Yang et al., 2024) was employed to identify the optimal parameter settings. For the NSGA-II algorithm, the parameters under consideration were population size (PS), crossover probability (Pc), mutation probability (Pm), and the number of generations (G). Each parameter was tested at four different levels. Consequently, an L₁₆(4⁴) orthogonal array from the Taguchi design was utilized to conduct the experiments for parameter tuning. The parameters for NSGA-II were fine-tuned to PS = 100, Pc = 0.90, Pm = 0.3, and G = 750. For the En-NSGA-II, the initial PS, Pc, and Pm were set to be identical to those of NSGA-II. To ensure a fair comparison in terms of the number of function evaluations (NFEs), the additional evaluations performed by the VNS within the En-NSGA-II algorithm must be considered. The VNS component of En-NSGA-II includes three different local searches, and in each generation, only one is applied to any given individual based on a predefined probability. Irrespective of the local search selected, five extra function evaluations will be carried out for the buffer section and five evaluations for the operator section. Consequently, each individual in En-NSGA-II undergoes 10 additional function evaluations in comparison to NSGA-II. Based on the aforementioned factors and keeping the population size constant as NSGA-II, taking into account the maximum VNS probability of 0.1 and the additional 10 evaluations per individual, the number of generations required for En-NSGA-II can be determined to be half of NSGA-II.

As for MOPSO, the population size and the number of generations for the MOPSO algorithm were set to be identical to those of NSGA-II to have the same NFEs. Additionally, the acceleration coefficients c1 and c2, as well as the inertia weight w, were independently tuned and found to be c1 = 1.25, c2 = 2, and w = 0.2.

## 10. Results and discussion

This section evaluates the performance of En-NSGA-II through extensive experiments and compares it with NSGA-II and MOPSO over 30 instances. The comparison begins with presenting the respective values generated by each algorithm for each performance measure,

<footnote>a Case study.</footnote>
<table><caption>**Table 4**
Coverage metric comparison.</caption><thead><tr><th>Problem size</th><th>Test instance</th><th>CM(En-NSGA-II, NSGA-II)</th><th>CM(NSGA-II, En-NSGA-II)</th><th>CM(En-NSGA-II, MOPSO)</th><th>CM(MOPSO, En-NSGA-II)</th></tr></thead><tbody><tr><td rowspan="10">Small</td><td>1</td><td><strong>0.92</strong></td><td>0.05</td><td><strong>0.87</strong></td><td>0.15</td></tr><tr><td>2</td><td><strong>0.78</strong></td><td>0.11</td><td><strong>0.89</strong></td><td>0.17</td></tr><tr><td>3</td><td><strong>0.88</strong></td><td>0.01</td><td><strong>0.90</strong></td><td>0.12</td></tr><tr><td>4</td><td><strong>0.86</strong></td><td>0.00</td><td><strong>0.82</strong></td><td>0.05</td></tr><tr><td>5</td><td><strong>0.97</strong></td><td>0.07</td><td><strong>0.90</strong></td><td>0.06</td></tr><tr><td>6</td><td><strong>0.89</strong></td><td>0.02</td><td><strong>0.98</strong></td><td>0.08</td></tr><tr><td>7</td><td><strong>0.82</strong></td><td>0.11</td><td><strong>0.85</strong></td><td>0.08</td></tr><tr><td>8</td><td><strong>0.87</strong></td><td>0.04</td><td><strong>0.95</strong></td><td>0.07</td></tr><tr><td>9</td><td><strong>0.93</strong></td><td>0.02</td><td><strong>0.85</strong></td><td>0.08</td></tr><tr><td>10</td><td><strong>0.92</strong></td><td>0.06</td><td><strong>0.89</strong></td><td>0.14</td></tr><tr><td rowspan="10">Medium</td><td>11</td><td><strong>0.81</strong></td><td>0.07</td><td><strong>0.96</strong></td><td>0.08</td></tr><tr><td>12</td><td><strong>0.97</strong></td><td>0.03</td><td><strong>0.87</strong></td><td>0.16</td></tr><tr><td>13</td><td><strong>0.92</strong></td><td>0.03</td><td><strong>0.82</strong></td><td>0.12</td></tr><tr><td>14</td><td><strong>0.82</strong></td><td>0.06</td><td><strong>0.77</strong></td><td>0.07</td></tr><tr><td>15</td><td><strong>0.87</strong></td><td>0.08</td><td><strong>0.79</strong></td><td>0.17</td></tr><tr><td>16</td><td><strong>0.86</strong></td><td>0.03</td><td><strong>0.96</strong></td><td>0.17</td></tr><tr><td>17</td><td><strong>0.95</strong></td><td>0.03</td><td><strong>0.78</strong></td><td>0.10</td></tr><tr><td>18</td><td><strong>0.84</strong></td><td>0.04</td><td><strong>0.87</strong></td><td>0.13</td></tr><tr><td>19</td><td><strong>0.94</strong></td><td>0.05</td><td><strong>0.98</strong></td><td>0.07</td></tr><tr><td>20</td><td><strong>0.81</strong></td><td>0.08</td><td><strong>0.97</strong></td><td>0.16</td></tr><tr><td rowspan="10">Large</td><td>21</td><td><strong>0.87</strong></td><td>0.00</td><td><strong>0.83</strong></td><td>0.12</td></tr><tr><td>22</td><td><strong>0.80</strong></td><td>0.00</td><td><strong>0.91</strong></td><td>0.08</td></tr><tr><td>23</td><td><strong>0.97</strong></td><td>0.04</td><td><strong>0.82</strong></td><td>0.06</td></tr><tr><td>24</td><td><strong>0.75</strong></td><td>0.05</td><td><strong>0.91</strong></td><td>0.10</td></tr><tr><td>25</td><td><strong>0.85</strong></td><td>0.01</td><td><strong>0.81</strong></td><td>0.15</td></tr><tr><td>26</td><td><strong>0.92</strong></td><td>0.06</td><td><strong>0.77</strong></td><td>0.07</td></tr><tr><td>27</td><td><strong>0.88</strong></td><td>0.08</td><td><strong>0.83</strong></td><td>0.14</td></tr><tr><td>28</td><td><strong>0.90</strong></td><td>0.01</td><td><strong>0.94</strong></td><td>0.11</td></tr><tr><td>29</td><td><strong>0.79</strong></td><td>0.02</td><td><strong>0.89</strong></td><td>0.12</td></tr><tr><td>30</td><td><strong>0.82</strong></td><td>0.02</td><td><strong>0.97</strong></td><td>0.17</td></tr><tr><td>Hit rate</td><td></td><td>30/30</td><td></td><td>30/30</td><td></td></tr></tbody></table>

followed by the relevant statistical analysis. Next, a comprehensive runtime analysis is performed to underscore the advantages of the proposed PMB-MOO over traditional simulation-based multi-objective optimization (SB-MOO) and to assess the computational efficiency of the algorithms under varying conditions. Subsequently, to evaluate the scalability of the algorithms, their performance is compared across small, medium, and large problems.

Table 4 presents the coverage metric of En-NSGA-II with NSGA-II and MOPSO. Bold values in Table 4 indicate superior performance, while the hit rate tallies the instances where an algorithm outperformed its counterparts. Table 4 shows the clear dominance of En-NSGA-II across all test instances. This superiority highlights En-NSGA-II's enhanced ability to converge on higher-quality solutions compared to NSGA-II and MOPSO.

Figs. 7(a) to 7(d) present the performance comparison of the algorithms over NNDS, HV, DM, and SM, across 30 test instances. Based on the t-test results summarized in Table 5, the performance of the algorithms across various measures can be analyzed as follows. When comparing the NNDSs found by each algorithm, En-NSGA-II significantly outperformed both NSGA-II and MOPSO. Regarding the DM measure, En-NSGA-II significantly outperformed NSGA-II, with a less significant difference between En-NSGA-II and MOPSO. For the SM measure, where a lower value is better, En-NSGA-II significantly outperformed both NSGA-II and MOPSO. In summary, the t-test results indicate that En-NSGA-II consistently outperformed NSGA-II and, in most cases, showed significant improvements over MOPSO. These results highlight the overall superior performance of En-NSGA-II in addressing HFS-B&RA, across the given test instances.

## 10.1. Runtime analysis

To highlight the advantages of the PMB-MOO approach over traditional SB-MOO, their runtime performances are compared. Additionally, the runtime of the algorithms is analyzed across 30 instances.

<table><caption>**Table 5**
Pairwise comparison of algorithms based on t-test results, across 30 test instances.</caption><thead><tr><th>Measure</th><th></th><th></th><th>Mean difference</th><th>p-value</th></tr></thead><tbody><tr><td rowspan="2">NNDS</td><td>En-NSGA-II</td><td>NSGA-II</td><td>6.15</td><td>0.00</td></tr><tr><td>En-NSGA-II</td><td>MOPSO</td><td>3.08</td><td>0.00</td></tr><tr><td rowspan="2">HV</td><td>En-NSGA-II</td><td>NSGA-II</td><td>2.00</td><td>0.05</td></tr><tr><td>En-NSGA-II</td><td>MOPSO</td><td>0.87</td><td>0.39</td></tr><tr><td rowspan="2">DM</td><td>En-NSGA-II</td><td>NSGA-II</td><td>2.33</td><td>0.02</td></tr><tr><td>En-NSGA-II</td><td>MOPSO</td><td>1.49</td><td>0.14</td></tr><tr><td rowspan="2">SM</td><td>En-NSGA-II</td><td>NSGA-II</td><td>-3.64</td><td>0.00</td></tr><tr><td>En-NSGA-II</td><td>MOPSO</td><td>-2.45</td><td>0.02</td></tr></tbody></table>

In the SB-MOO approach, the high-fidelity DES model is applied directly within the optimization loop, resulting in substantial computational time due to the complexity of the simulations. For instance, running the simulation model (of the case study) for 100 replications requires approximately 150 ± 5 min. In contrast, the PMB-MOO approach replaces these costly simulation runs with a fast-running PM. The PM completes the same 100 evaluations in three milliseconds. This significant decrease in computational time enables rapid iterations within the optimization process, facilitating the exploration of larger solution spaces and the attainment of optimal results more efficiently. It is worth mentioning that SB-MOO was performed only on the case study due to its extensive computational time.

Fig. 8 provides the comparison of runtimes for 30 instances using EN-NSGA-II, NSGA-II, and MOPSO. According to Fig. 8 as the instance size increases, the runtimes of all three algorithms exhibit a predictable rise due to growing computational complexity. Among the algorithms, EN-NSGA-II exhibits marginally extended runtimes, ranging from 3 to 150 s (14% on average) longer than NSGA-II across various instances. However, this modest computational overhead is justified by the algorithm's superior exploration capabilities, as evidenced by its enhanced performance across all quality metrics (NNDS, HV, DM, and SM), as illustrated in Fig. 7.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Performance of algorithms in 30 test instances.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Runtime comparison of EN-NSGA-II, NSGA-II, and MOPSO.</figcaption></figure>

This trade-off between computational cost and solution quality aligns with the fundamental principle in evolutionary computation that more thorough exploration of the objective space often requires additional computational resources (Deb and Jain, 2014). The observed time-quality trade-off is particularly favorable, as the relatively small increase in computational overhead yields significant improvements in the diversity and convergence of the Pareto front. This phenomenon is consistent with previous findings in multi-objective optimization literature, where enhanced exploration mechanisms have demonstrated that marginal increases in computational time can lead to substantial improvements in solution quality.

It is worth noting that the increased runtime in EN-NSGA-II can be attributed to two distinct factors, each contributing differently to the overall computational overhead. The first factor contributing to

the runtime increase in EN-NSGA-II is the integration of an additional learning loop. However, this factor has a minimal impact on the overall runtime due to the following considerations. The meta-learning loop is specifically designed for efficiency, operating under a time constraint controlled by the parameter *lr_timeout* (see Algorithm 1), which ensures that the computational overhead remains manageable. Furthermore, the total iterations of the learning loop are constrained by the product of the size of orthogonal arrays and the number of learning iterations (which remains less than 100 iterations).

On the other hand, the second factor, which arises from EN-NSGA-II identifying a larger number of non-dominated solutions compared to NSGA-II and MOPSO, has a more significant impact on runtime. This increased computational overhead is due to the additional time required to evaluate and process the larger set of non-dominated solutions (while
performing non-dominated sorting and crowding distance procedures), contributing substantially to the overall runtime increase.

## 10.2. Scalability of the algorithms

Fig. 9 presents box-plots that analyze the scalability of En-NSGA-II, NSGA-II, and MOPSO across different problem sizes based on various performance measures, complementing the pairwise comparison in Table 5. For the NNDS measure (Figs. 9(a)–9(c)), En-NSGA-II demonstrates superior scalability, consistently achieving higher mean values across medium and large size problems when compared to NSGA-II and MOPSO. Similarly, En-NSGA-II outperforms the other algorithms in terms of the HV measure (Figs. 9(d)–9(f)), achieving the highest mean HV values across all problem sizes. This indicates better scalability for En-NSGA-II in terms of hypervolume. Regarding the DM measure (Figs. 9(g)–9(i)), En-NSGA-II consistently achieves the highest mean values across all problem sizes, suggesting better performance and scalability compared to NSGA-II and MOPSO. For the SM measure (Figs. 9(j)–9(l)), En-NSGA-II performs better by achieving lower mean values across all problem sizes compared to NSGA-II and MOPSO, indicating more favorable performance and scalability.

Overall, En-NSGA-II demonstrates superior scalability in terms of NNDS, HV, and DM, achieving higher values across different problem sizes. MOPSO generally performs moderately across all measures. These insights provide a comprehensive understanding of the scalability and performance of the algorithms across various measures and problem sizes.

## 10.3. Post-optimality analysis

The results yielded by MOO provide many near-optimal choices (Ramírez et al., 2019). Consequently, an element of ambiguity is created that renders it difficult for decision-makers to reach their final decision. To avoid this ambiguity, decision-makers should be provided with simple and easy-to-implement rules. In this study, a data mining algorithm known as flexible pattern mining (FPM) (Smedberg and Bandaru, 2022) was applied to extract applicable and transferable knowledge from MOO results. Fig. 10(a) depicts the parallel coordinates for the objective (first four coordinates from the left) and decision space (fifth coordinate from the left to the last).

To extract practical knowledge from MOO results (Fig. 10(a)), the PM was optimized using the En-NSGA-II algorithm, which demonstrated superior performance compared to NSGA-II and MOPSO. Then, solutions with the highest throughput were selected as the desirable solutions. Next, rules for the decision variables ($b_s$, $Op_s$, AMR) that result in desirable solutions, i.e., solutions with the highest throughput, were identified. The association rules generated via the FPM exhibit two important characteristics: significance (sig) and unselected significance (unsig). Sig is the frequency of generated association rules for decision variables that result in the highest throughput, whereas unsig is the frequency of the same rules in solutions with a lower throughput. An FPM rule with both a high sig and a low unsing is desirable.

The selected set shown in Fig. 10(b) represents the solutions with the highest throughputs and the lowest number of operators. The FPM was later applied to the selected solutions to generate association rules with Sig values above 95% and unsig values below 5%. The generated association rules are listed in Table 6. As shown in the second column of Table 6, the rules were generated at three to five levels (a rule with n levels was constructed using n variables). Rule extraction was performed separately for the buffers and operators.

Table 6 lists some rules for buffers and operators that repeatedly appeared while FPM was performed on desirable solutions (higher throughput and lower number of buffers, operators, and AMRs). All high-Sig and low-unsig rules included mutual rules for buffer capacities after segments 3 and 4 ($B3 \le 3 \& B4 = 10$) and the number of operators

in segments 4 and 5 ($Op4 \ge 3 \& Op5 = 7$), thus highlighting $B3$, $B4$, $Op4$, and $Op5$ for further investigation.

To obtain deeper knowledge and assist decision-makers with practical and easy-to-comprehend rules, FPM was performed on the results of the 30 solved instances. The rules generated for $B3$ and $Op4$ for 30 instances were inconclusive, thus resulting in diverse conclusions. Nevertheless, the rules generated for $Op5$ and $B4$ were conclusive. Thus, the rules generated for $Op5$ and $B4$ are discussed in detail. In all instances, the generated rules were in the form of “equal to” or “greater than and equal to”. For example, in instance 8, where the lower bound was equal to 1 and the upper bound was equal to 9, the generated rule for $Op5$ was $Op5 \ge 7$. Similarly, in instance 20, where the lower bound was equal to 1 and the upper bound was equal to 10, the rule for $B4$ was $B4 \ge 8$. The analysis of the rules generated over 30 random instances resulted in the following general conclusions:

1. In all instances of the problem, wherein the upper bound for $Op5$ exceeded 4, the best results were achieved when the $Op5$ value was equal to or greater than 70% of the permitted upper bound. Therefore, the number of operators allocated after the paint shop should be greater than 70% of its upper bound when the upper bound exceeds 4.

2. In all instances of the problem, wherein the upper bound for $B4$ exceeded 4, the best results were achieved when the $B4$ value was equal to or greater than 67% of the permitted upper bound. When the upper bound for $B4$ exceeded 6, the best results were achieved when the $B4$ value was equal to or greater than 80% of the permitted upper bound. Therefore, the size of $B4$ (located before the paint shop) should be greater than 67% of its upper bound when the upper bound is less than or equal to 6. When the upper bound is greater than or equal to 6, the size of $B4$ should be greater than 80% of the upper bound.

The above rules can be summarized in Eq. (12):

$$ \left\{ \begin{array}{ll} Op5^* \ge 0.7 \times \max(Op5), & \text{where } \max(Op5) > 4 \\ B4^* \ge 0.67 \times \max(B4), & \text{where } 4 < \max(B4) \le 6 \\ B4^* \ge 0.8 \times \max(B4), & \text{where } \max(B4) > 6 \end{array} \right. \qquad (12) $$

where * shows the proposed values.

## 11. Managerial insight

These analyses indicate that the following three essential factors can assist line managers in handling HFSs. First, the buffer allocated after the feed-forward quality inspection significantly affects the performance of the production line. In other words, the size of the buffer in this area must be allocated such that it can absorb the fluctuations and variations that occur because of the quality inspection. These fluctuations can lead to sequence scrambling, where the order of products on the production line is disrupted due to the need to remove, repair, and reintegrate defective products. An adequately sized buffer after the quality inspection station is crucial for maintaining the flow of production, as it provides a temporary holding area for products awaiting repair or reintegration. This buffer acts as a shock absorber, preventing the immediate effects of sequence scrambling from propagating downstream and turning subsequent stations into bottlenecks. Sequence scrambling poses a significant challenge, as it can disrupt the planned order of operations and lead to inefficiencies, such as increased waiting times and reduced throughput. The insights from the FPM analysis suggest that addressing sequence scrambling requires a holistic approach that includes both buffer management and operator allocation. By ensuring that buffers are sufficiently large to handle the variations introduced by quality inspections and that manual sections are adequately staffed, production lines can mitigate the adverse effects of sequence scrambling. This approach not only helps maintain the efficiency and productivity of the production line but also ensures the
<figure><img src="image_11.png" /><figcaption>Fig. 9. Box-plots for NNDS, HV, DM, and SM over small, medium, and large size problems.</figcaption></figure>

<table><caption>**Table 6**
Rules generated via FPM.</caption><thead><tr><th>Variable name</th><th>Rule levels<sup>a</sup></th><th>Generated rules</th><th>Sig</th><th>Unsig</th></tr></thead><tbody><tr><td rowspan="3">Buffers</td><td>3-level</td><td>B2 == 3 &amp; B3 ≤ 3 &amp; B4 == 10</td><td>97.5</td><td>0.27</td></tr><tr><td>4-level</td><td>B1 == 1 &amp; B3 ≤ 3 &amp; B4 == 10 &amp; B9 == 2</td><td>96.2</td><td>1.31</td></tr><tr><td>5-level</td><td>B1 == 1 &amp; B2 == 3 &amp; B3 ≤ 3 &amp; B4 == 10 &amp; B6 ≥ 2</td><td>95.1</td><td>2.67</td></tr><tr><td rowspan="3">Operators</td><td>3-level</td><td>Op1 == 5 &amp; Op4 ≥ 3 &amp; Op5 == 7</td><td>98.2</td><td>1.64</td></tr><tr><td>4-level</td><td>Op1 == 5 &amp; Op2 ≥ 5 &amp; Op4 ≥ 3 &amp; Op5 == 7</td><td>97.8</td><td>3.47</td></tr><tr><td>5-level</td><td>Op1 == 5 &amp; Op2 ≥ 5 &amp; Op3 == 7 &amp; Op4 ≥ 3 &amp; Op5 == 7</td><td>95.5</td><td>4.25</td></tr></tbody></table>

<footnote>a Number of variables included in generated rules.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Objective and decision space (a) filtered by solutions with the highest throughput (b).</figcaption></figure>

high quality of the products by allowing for the repair and reintegration of defective items without compromising the overall production flow.

Second, Designers of the production line must prioritize examining the area of the HFS with mixed automation levels, in which automation changes among manual, semi-automatic, and fully automatic. This is because of the inherent stochasticity in the processing time, which changes significantly in the abovementioned areas along the line. Therefore, sufficiently large buffers must be used to absorb the fluctuations between the manual and automatic sections of the line. Additionally, a sufficient number of operators must be assigned to the manual sections of the line that appear after a fully automated section.

Last but not least, combining simulation, optimization, and artificial intelligence, including FPM, as part of the Industry 4.0 revolution (Mahmoodi et al., 2022), has been a significant development, emphasizing the need for a holistic approach to enhance the production line performance and manufacturing processes. This integration is crucial for leveraging the vast amounts of data generated in modern manufacturing environments. This comprehensive strategy not only streamlines operations but also significantly reduces downtime and waste, leading to enhanced productivity and global market competitiveness. Given these advancements, managers should increasingly adopt these technologies to harness their full potential. This proactive adoption is essential for realizing the benefits of reduced operational costs and increased production efficiency, ultimately contributing to the organization's long-term success.

## 12. Conclusions

Inspired by a real-world problem in the marine engine production industry, this study addressed a highly stochastic and complex problem of simultaneous buffers and resources (B&RA) problem in high-mix, low-volume hybrid flow shops (HFSs) with feedforward quality inspection. This problem is known as HFS-B&RA and aims to optimize the allocation of operators and the sizing of the autonomous mobile robots (AMRs) fleet to maximize the throughput and minimize the number of buffers, operators, and AMRs.

A discrete-event simulation (DES) was employed to manage the stochasticity inherent in the production environment. Because of the long runtime of the simulation model, five widely used regressor

machine learning (ML) algorithms in manufacturing industries were compared, and the best-performing one, i.e., random forest regressor, was used to construct a predictive model (PM) to replace the time-consuming simulation model. To solve the HFS-B&RA, an enhanced non-dominated sorting genetic algorithm (En-NSGA-II) equipped with a variable neighborhood search (VNS) and a lifelong meta-learning (LLML) mechanism was presented and coupled with the PM to form a PM-based multi-objective optimization (PMB-MOO) approach. The performance of the proposed En-NSGA-II was compared with those of the NSGA-II and multi-objective particle swarm optimization (MOPSO) using a set of instances, including a real case. The comprehensive analysis of the performance measures demonstrated the superiority of the En-NSGA-II algorithm. The results of the application of En-NSGA-II to the case study yielded significant improvements in system performance. Notably, the implementation of En-NSGA-II resulted in an approximate 14% enhancement in system throughput. Concurrently, the algorithm achieved a reduction in the total number of buffers from 41 to 35 and a decrease in the total number of operators from 44 to 36 without necessitating an increase in the number of AMRs.

Additionally, the rules generated using knowledge extraction over all instances (Eq. (12)) revealed that in a hybrid flow shop (HFS) with feedforward quality inspection, the size of buffers located immediately after the inspection process plays a crucial role in the system's performance. Another important aspect of HFSs is the location where the automation level differs from one segment to another, specifically if two consecutive segments of the HFS are joined rigidly without an appropriate location for buffers. In this scenario, synchronizing the production speed of the segments in which the automation level differs (e.g., from automatic to manual or vice versa) becomes challenging. Based on the optimization results and patterns generated in the computational study, it was observed that the difference between variations in two consecutive segments adversely affected the general performance of the production line.

Future research pertaining to B&RA presents immense potential for optimizing production processes and enhancing overall efficiency. A promising area of investigation is the application of B&RA strategies in diverse production environments, including closed-loop systems and job shops. Integrating logistics with the BAP and RAP would result in a novel set of combinatorial optimization problems wherein the supply to
the production line cannot be presumed to be unlimited. Furthermore,
determining the appropriate dimensions for sideline inventory buffers
(the buffers required for maintaining inventory), as well as the buffer
spaces between machines, introduces a new level of complexity. In
addition, identifying the optimal placement of sideline buffers and
effectively allocating both buffers and resources may yield promising
results. Another recent trend for further research is the combination
of simulation and reinforcement learning (sim-learning) in mobile re-
source allocation (operators and AMRs), which represents a promising
shift towards more responsive and intelligent manufacturing systems.
Sim-learning, which merges simulation with the learning algorithms,
enables systems to learn and adapt from historical data, optimizing
resource distribution in real-time. Sim-learning method is particularly
advantageous in complex production environments characterized by
variability and uncertainty, such as those with feed-forward quality
inspection and sequence scrambling.

CRediT authorship contribution statement

**Ehsan Mahmoodi:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Masood Fathi:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. **Amos H.C. Ng:** Writing – review & editing, Supervision, Project administration, Funding acquisition. **Alexandre Dolgui:** Writing – review & editing, Supervision.

Acknowledgments

The authors gratefully acknowledge funding from the
Sweden Knowledge Foundation (KKS) through the ACCURATE 4.0
project (grant agreement No. 20200181) and extend their gratitude to
Volvo Penta of Sweden for their collaborative support throughout this
study.

Data availability

Data will be made available on request.

References

Abu Qudeiri, J., Yamamoto, H., Ramli, R., Jamali, A., 2008. Genetic algorithm for buffer size and work station capacity in serial-parallel production lines. Artif. Life Robot. 12 (1), 102–106. http://dx.doi.org/10.1007/s10015-007-0449-5

Altparmak, F., Dengiz, B., Bulgak, A.A., 2007. Buffer allocation and performance modeling in asynchronous assembly system operations: An artificial neural network metamodeling approach. Appl. Soft Comput. 7 (3), 946–956. http://dx.doi.org/10.1016/j.asoc.2006.06.002

Amiri, M., Mohtashami, A., 2012. Buffer allocation in unreliable production lines based on design of experiments, simulation, and genetic algorithm. Int. J. Adv. Manuf. Technol. 62 (1–4), 371–383. http://dx.doi.org/10.1007/s00170-011-3802-8

Bahtat, C., El Barkany, A., Abdelouahhab, J., 2023. Reconfigurable manufacturing systems: From automation through industry 4.0. Int. J. Ind. Eng. Prod. Res. ..

Belete, D.M., Huchaiah, M.D., 2022. Grid search in hyperparameter optimization of machine learning models for prediction of HIV/AIDS test results. Int. J. Comput. Appl. 44 (9), 875–886. http://dx.doi.org/10.1080/1206212X.2021.1974663

Bulgak, A.A., 2006. Analysis and design of split and merge unpaced assembly systems by metamodelling and stochastic search. Int. J. Prod. Res. 44 (18–19), 4067–4080. http://dx.doi.org/10.1080/00207540600564625

Chang, K.-H., Huang, C.-C., 2022. An interactive two-stage framework for simultaneous machine selection and buffer allocation. IEEE Trans. Autom. Sci. Eng. 1–15. http://dx.doi.org/10.1109/TASE.2022.3173641

Chicco, D., Warrens, M.J., Jurman, G., 2021. The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation. PeerJ Comput. Sci. 7, e623. http://dx.doi.org/10.7717/peerj-cs.623

Colledani, M., Gandola, F., Matta, A., Tolio, T., 2008. Performance evaluation of linear and non-linear multi-product multi-stage lines with unreliable machines and finite homogeneous buffers. IIE Trans. 40 (6), 612–626. http://dx.doi.org/10.1080/07408170701745345

Deb, K., Jain, H., 2014. An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints. IEEE Trans. Evol. Comput. 18 (4), 577–601. http://dx.doi.org/10.1109/TEVC.2013.2281535.

Dieleman, N.A., Berkhout, J., Heidergott, B., 2023. A neural network approach to performance analysis of tandem lines: The value of analytical knowledge. Comput. Oper. Res. 106124. http://dx.doi.org/10.1016/j.cor.2022.106124.

do Amaral, J.V.S., de Carvalho Miranda, R., Montevechi, J.A.B., dos Santos, C.H., Gabriel, G.T., 2022. Metamodeling-based simulation optimization in manufacturing problems: a comparative study. Int. J. Adv. Manuf. Technol. 120 (7), 5205–5224. http://dx.doi.org/10.1007/s00170-022-09072-9

Dolgui, A., Eremeev, A., Kolokolov, A., Sigaev, V., 2002. A genetic algorithm for the allocation of buffer storage capacities in a production line with unreliable machines. J. Math. Model. Algorithms 1 (2), 89–104. http://dx.doi.org/10.1023/A:1016560109076, 2002 1: 2.

Dolgui, A.B., Eremeev, A.V., Sigaev, V.S., 2017. Analysis of a multicriterion buffer capacity optimization problem for a production line. Autom. Remote. Control. 78 (7), 1276–1289. http://dx.doi.org/10.1134/S0005117917070098, 2017 78: 7.

Duboue, P., 2020. The Art of Feature Engineering: Essentials for Machine Learning. Cambridge University Press, Cambridge, http://dx.doi.org/10.1017/9781108671682.

El Mrabet, Z., Sugunaraj, N., Ranganathan, P., Abhyankar, S., 2022. Random forest regressor-based approach for detecting fault location and duration in power systems. Sensors 22 (2), 458. http://dx.doi.org/10.3390/s22020458.

Escudero-Ornelas, I.O., Tiwari, D., Farnsworth, M., Tiwari, A., 2023. Modelling interdependencies in an electric motor manufacturing process using discrete event simulation. J. Simul. 1–22. http://dx.doi.org/10.1080/17477778.2023.2202338.

Fares, D., Fathi, M., Mekhilef, S., 2022. Performance evaluation of metaheuristic techniques for optimal sizing of a stand-alone hybrid PV/Wind/Battery system. Appl. Energy 305, 117823. http://dx.doi.org/10.1016/j.apenergy.2021.117823.

Golmohammadi, D., 2015. A study of scheduling under the theory of constraints. Int. J. Prod. Econ. 165, 38–50. http://dx.doi.org/10.1016/j.ijpe.2015.03.015.

Hammedi, S., Elmelliani, J., Nabli, L., 2024. Optimizing production in reconfigurable manufacturing systems with artificial intelligence and Petri nets. Int. J. Adv. Comput. Sci. Appl. 15 (10), 413–420.

Herps, K., Dang, Q.-V., Martagan, T., Adan, I., 2022. A simulation-based approach to design an automated high-mix low-volume manufacturing system. J. Manuf. Syst. 64, 1–18. http://dx.doi.org/10.1016/j.jmsy.2022.05.013.

Jeong, H.Y., Kim, Y., Lee, S., Moreland, J., Zhou, C., 2023. Disruption propagation and repair response in interdependent system: Network model and simulation approach. Simul. Model. Pr. Theory 124, 102730. http://dx.doi.org/10.1016/j.simpat.2023.102730.

Kassoul, K., Cheikhrouhou, N., Zufferey, N., 2022. Buffer allocation design for unreliable production lines using genetic algorithm and finite perturbation analysis. Int. J. Prod. Res. 60 (10), 3001–3017. http://dx.doi.org/10.1080/00207543.2021.1909169.

Komaki, G.M., Sheikh, S., Malakooti, B., 2019. Flow shop scheduling problems with assembly operations: a review and new trends. Int. J. Prod. Res. 57 (10), 2926–2955. http://dx.doi.org/10.1080/00207543.2018.1550269.

Law, A.M., 2022. How to build valid and credible simulation models. In: 2022 Winter Simulation Conference. WSC, pp. 1283–1295. http://dx.doi.org/10.1109/WSC57314.2022.10015411.

Li, L., Firouz, M., Ahmed, A., Delen, D., 2023a. On the Egalitarian–Utilitarian spectrum in stochastic capacitated resource allocation problems. Int. J. Prod. Econ. 262, 108900. http://dx.doi.org/10.1016/j.ijpe.2023.108900.

Li, P., Xue, Q., Zhang, Z., Chen, J., Zhou, D., 2023b. Multi-objective energy-efficient hybrid flow shop scheduling using Q-learning and GVNS driven NSGA-II. Comput. Oper. Res. 159, 106360. http://dx.doi.org/10.1016/j.cor.2023.106360.

Luo, X., Li, W., Tu, Y., Xue, D., Tang, J., 2010. Optimal resource allocation for hybrid flow shop in one-of-a-kind production. Int. J. Comput. Integr. Manuf. 23 (2), 146–154. http://dx.doi.org/10.1080/09511920903440339.

Mahmoodi, E., Fathi, M., Ghobakhloo, M., 2022. The impact of industry 4.0 on bottleneck analysis in production and manufacturing: Current trends and future perspectives. Comput. Ind. Eng. 174, 108801. http://dx.doi.org/10.1016/j.cie.2022.108801.

Mahmoodi, E., Fathi, M., Tavana, M., Ghobakhloo, M., Ng, A.H.C., 2024. Data-driven simulation-based decision support system for resource allocation in industry 4.0 and smart manufacturing. J. Manuf. Syst. 72, 287–307. http://dx.doi.org/10.1016/j.jmsy.2023.11.019.

Malakouti, S.M., Menhaj, M.B., Suratgar, A.A., 2023. The usage of 10-fold cross-validation and grid search to enhance ML methods performance in solar farm power generation prediction. Clean. Eng. Technol. 15, 100664. http://dx.doi.org/10.1016/j.clet.2023.100664.

Massim, Y., Yalaoui, F., Amodeo, L., Chatelet, E., Zeblah, A., 2010. Efficient combined immune-decomposition algorithm for optimal buffer allocation in production lines for throughput and profit maximization. Comput. Oper. Res. 37 (4), 611–620. http://dx.doi.org/10.1016/j.cor.2009.06.016.
Mohtashami, A., 2014. A new hybrid method for buffer sizing and machine allocation in unreliable production and assembly lines with general distribution time-dependent parameters. Int. J. Adv. Manuf. Technol. 74 (9–12), 1577–1593. http://dx.doi.org/10.1007/s00170-014-6098-7.

Nahas, N., Nourelfath, M., Ait-Kadi, D., 2008. Ant colonies for structure optimisation in a failure prone series-parallel production system. J. Qual. Maint. Eng. 14 (1), 7–33. http://dx.doi.org/10.1108/13552510810861914/FULL/PDF.

Nahas, N., Nourelfath, M., Ait-Kadi, D., 2009. Selecting machines and buffers in unreliable series-parallel production lines. Int. J. Prod. Res. 47 (14), 3741–3774. http://dx.doi.org/10.1080/00207540701806883.

Neumann, A., Hajji, A., Rekik, M., Pellerin, R., 2024. Genetic algorithms for planning and scheduling engineer-to-order production: a systematic review. Int. J. Prod. Res. http://dx.doi.org/10.1080/00207543.2023.2237122.

Ng, A.H.C., Bernedixen, J., 2018. Production systems analysis and optimization using facts analyzer. In: Proceedings of the 2018 Winter Simulation Conference. In: WSC '18, IEEE Press, Gothenburg, Sweden, p. 4258.

Norouzi, A., Heidarifar, H., Borhan, H., Shahbakhti, M., Koch, C.R., 2023. Integrating machine learning and model predictive control for automotive applications: A review and future directions. Eng. Appl. Artif. Intell. 120, 105878. http://dx.doi.org/10.1016/j.engappai.2023.105878.

Nourmohammadi, A., Fathi, M., Ng, A.H., 2022. Balancing and scheduling assembly lines with human-robot collaboration tasks. Comput. Oper. Res. 140, 105674. http://dx.doi.org/10.1016/j.cor.2021.105674.

Nourmohammadi, A., Ng, A.H.C., Fathi, M., Vollebregt, J., Hanson, L., 2023. Multi-objective optimization of mixed-model assembly lines incorporating musculoskeletal risks assessment using digital human modeling. CIRP J. Manuf. Sci. Technol. 47, 71–85. http://dx.doi.org/10.1016/j.cirpj.2023.09.002.

Perez-Gonzalez, P., Framinan, J.M., 2024. A review and classification on distributed permutation flowshop scheduling problems. European J. Oper. Res. 312 (1), 1–21. http://dx.doi.org/10.1016/j.ejor.2023.02.001.

Prayudani, S., Hizriadi, A., Lase, Y.Y., Fatmi, Y., Al-Khowarizmi, 2019. Analysis accuracy of forecasting measurement technique on random K-nearest neighbor (RKNN) using MAPE and MSE. J. Phys.: Conf. Ser. 1361 (1), 012089. http://dx.doi.org/10.1088/1742-6596/1361/1/012089.

Puchkova, A., McFarlane, D., Srinivasan, R., Thorne, A., 2020. Resilient planning strategies to support disruption-tolerant production operations. Int. J. Prod. Econ. 226, 107614. http://dx.doi.org/10.1016/j.ijpe.2020.107614.

Rabet, R., Ganji, M., Fathi, M., 2024. A simheuristic approach towards supply chain scheduling: Integrating production, maintenance and distribution. Appl. Soft Comput. 153, 111264. http://dx.doi.org/10.1016/j.asoc.2024.111264.

Rajabi, M.M., Ataie-Ashtiani, B., Janssen, H., 2015. Efficiency enhancement of optimized Latin hypercube sampling strategies: Application to Monte Carlo uncertainty analysis and meta-modeling. Adv. Water Resour. 76, 127–139. http://dx.doi.org/10.1016/j.ADVWATRES.2014.12.008.

Ramírez, A., Romero, J.R., Ventura, S., 2019. A survey of many-objective optimisation in search-based software engineering. J. Syst. Softw. 149, 382–395. http://dx.doi.org/10.1016/j.JSS.2018.12.015.

Shahvari, O., Logendran, R., 2018. A comparison of two stage-based hybrid algorithms for a batch scheduling problem in hybrid flow shop with learning effect. Int. J. Prod. Econ. 195, 227–248. http://dx.doi.org/10.1016/j.ijpe.2017.10.015.

Shao, H., Moroni, G., Li, A., Liu, X., Xu, L., 2020. Simultaneously solving the transfer line balancing and buffer allocation problems with a multi-objective approach. J. Manuf. Syst. 57, 254–273. http://dx.doi.org/10.1016/j.jmsy.2020.09.009.

Smedberg, H., Bandaru, S., 2022. Interactive knowledge discovery and knowledge visualization for decision support in multi-objective optimization. European J. Oper. Res. http://dx.doi.org/10.1016/j.ejor.2022.09.008.

Tsadiras, A.K., Papadopoulos, C.T., O’Kelly, M.E., 2013. An artificial neural network based decision support system for solving the buffer allocation problem in reliable production lines. Comput. Ind. Eng. 66 (4), 1150–1162. http://dx.doi.org/10.1016/j.CIE.2013.07.024.

Wu, X., Cao, Z., 2022. An improved multi-objective evolutionary algorithm based on decomposition for solving re-entrant hybrid flow shop scheduling problem with batch processing machines. Comput. Ind. Eng. 169, 108236. http://dx.doi.org/10.1016/j.cie.2022.108236.

Wu, K., Zheng, M., Shen, Y., 2020. A generalization of the theory of constraints: Choosing the optimal improvement option with consideration of variability and costs. IISE Trans. 52 (3), 276–287. http://dx.doi.org/10.1080/24725854.2019.1632503.

Yang, P., Wang, L., Qing, T., Liu, Y., 2024. Multi-objective optimization of geometrical parameters of laterally perforated on silt fin by gray relation analysis based on Taguchi method. Int. J. Therm. Sci. 196, 108705. http://dx.doi.org/10.1016/j.ijthermalsci.2023.108705.

Yelles-Chaouche, A.R., Gurevsky, E., Brahimi, N., Dolgui, A., 2021. Reconfigurable manufacturing systems from an optimisation perspective: a focused review of literature. Int. J. Prod. Res. 59 (21), 6400–6418. http://dx.doi.org/10.1080/00207543.2020.1813913.

Zhou, B., Bao, J., Li, J., Lu, Y., Liu, T., Zhang, Q., 2021. A novel knowledge graph-based optimization approach for resource allocation in discrete manufacturing workshops. Robot. Comput.-Integr. Manuf. 71, 102160. http://dx.doi.org/10.1016/j.rcim.2021.102160.

Zhou, A., Qu, B.-Y., Li, H., Zhao, S.-Z., Sughanthan, P.N., Zhang, Q., 2011. Multiobjective evolutionary algorithms: A survey of the state of the art. Swarm Evol. Comput. 1 (1), 32–49. http://dx.doi.org/10.1016/j.swevo.2011.03.001.