Multi-objective optimization of mixed-model assembly lines incorporating musculoskeletal risks assessment using digital human modeling

the assembly lines for different product variants


Introduction
Notwithstanding the recent automation trends in a wide range of manufacturing industries owing to the Industry 4.0 paradigm, manual assembly lines are still prevalent because of the flexibility and capability of human labor.The human-centered assembly lines enable product customization while responding to the versatile demands of customers in the present competitive markets [1].Conventionally, companies have attempted to achieve higher productivity while configuring their assembly lines by distributing the tasks among the workstations as smoothly as possible, known as the assembly line balancing problem (ALBP) [2].Productivity can be improved by reducing inputs (e.g., a smaller number of workstations or less worker time) or increasing outputs (e.g., a higher throughput or delivery rate).
Over the past decades, assembly lines have been extended broadly to adapt to the production paradigms.One of the practical and effective extensions of the assembly lines that emerged with the mass customization trends is the production of several product models on the same assembly line.Accordingly, assembly lines have been categorized into single-model, multi-model, and mixed-model assembly lines [3].Mixed-model assembly lines (MMALs) are used widely in many industries to massively produce different models of products while responding to customers' demand for customized products.Notwithstanding the cost-effectiveness of MMALs, their efficient design and planning are challenging for manufacturers [4].One of the significant challenges in the (re)configuration of such assembly lines is the simultaneous balancing of workload for producing various models of products, known as the mixed-model assembly line balancing problem (MMALBP) while optimizing productivity measures.
On the other hand, human workers performing assembly line activities are subjected to various hazardous scenarios such as large weights, repetitive motions, and awkward postures.Impelling workers toward higher productivity without considering such pressures can significantly increase work-related musculoskeletal disorders (WMSDs), which can negatively affect the productivity of the assembly lines.Occupational disease studies in Europe have shown that using hand tools, exposure to mechanical pressure on the upper limbs, and awkward postures are among the most prevalent risk factors for WMSDs, particularly in the automotive industry [5].The statistics in the area of occupational safety and health (OSH) show that approximately 4% of the global gross domestic product is lost owing to accidents and unhealthy working conditions [6].In this regard, in line with the recent sustainable, human-centric, and resilient development goals initiated by Industry 5.0, many manufacturing industries have promoted better health and wellbeing for their workforce.Therefore, ergonomics is vital besides productivity measures while designing today's assembly lines.
The literature embraces appropriate studies in which human well-being aspects were incorporated in ALBP in terms of wellknown energy expenditure equations and/or rest allowance or fatigue [7,8] or ergonomic risk assessment methods [9,10].The latter studies included the National Institute for Occupational Safety and Health (NIOSH), Ovako Working posture Assessment System (OWAS), Occupational Repetitive Action (OCRA), rapid upper limb assessment (RULA), and rapid entire body analysis (REBA) as the most established risk assessment methods in practice due to their precise measurements of musculoskeletal risks [8,11].However, there are still empty gaps in the literature where the risk assessment methods (e.g., REBA) are integrated into MMALBP due to their numerous observational assessments [12].
On the other hand, recently the study of human well-being and ergonomic risk assessment has been enriched by the emergence of digital human modeling (DHM) and simulation modeling [13,14].Nevertheless, minimal research efforts have been undertaken to incorporate DHM in ergonomic assessment while addressing ALBP [15,16].Moreover, dealing with multiple objectives while addressing MMALBP with both ergonomic and economic aspects need tackling them separately rather than combining them into single-objective is another research gap as most of the studies have integrated the ergonomic and economic objectives [16].
Therefore, this study attempts to include musculoskeletal risk assessment through DHM while addressing MMALBP based on a multi-objective optimization approach.This study uses REBA for musculoskeletal risk assessment since (1) it is among the most widely used methods by ergonomists who assess postures of the upper limbs and the entire body [17], (2) only a few studies, particularly in MMALBP, have been conducted using REBA as the musculoskeletal risk assessment method, and (3) it is the most comparable to the specific ergonomic assessment method for measuring full-body posture and joint angles.
To find a trade-off between productivity and human well-being while addressing MMALBP, multiple objectives are considered, namely minimizing: (1) cycle time (CT ), (2) maximum ergonomic risk of workstations, and (3) the total ergonomic risk.The considered MMALBP is formulated as a mixed-integer linear programming (MILP) model considering the number of workstations, the precedence relationships among tasks for the mixed-model, and the musculoskeletal risks of the workstations based on REBA obtained by DHM.To deal with this multi-objective optimization problem, an enhanced non-dominated sorting genetic algorithm (NSGA-II) is proposed with customized mechanisms based on the problem characteristics.The enhanced NSGA-II (E-NSGA-II) has two embedded features: local search and multi-criteria decision-making (MCDM).The local search mechanism finds neighborhood solutions for the current feasible solutions.The MCDM mechanism ranks the neighborhood solutions using the technique for order of preference by similarity to the ideal solution (TOPSIS).The former feature enhances the diversification capability of the algorithm, while the latter improves the intensification capability of the algorithm.The E-NSGA-II is benchmarked against exact and multi-objective optimization algorithms while addressing a case study and test problems taken from the literature.The performance of the proposed algorithm is evaluated through statistical tests.In summary, this study aims to contribute to the MMALBP literature by: • Using DHM in evaluating the musculoskeletal risk of operators based on REBA, • Formulating the MMALBP with REBA as a MILP, • Developing an enhanced multi-objective optimization algorithm based on NSGA-II.
• Analysis of results based on a case study and standard test pro- blems.
The rest of this paper is organized as follows: Section 2 reviews the relevant literature.Section 3 presents the problem formulation, including the digital human modeling for measuring musculoskeletal risks and the proposed MILP.Section 4 describes the developed enhanced multi-objective optimization algorithm.In Section 5, the computational results show the comparison of the proposed algorithm with both exact and meta-heuristic approaches on a case study and test problems taken from the literature, followed by the analysis of results, and managerial insights.Finally, the conclusions and future research directions are outlined in Section 6.

MMALBP and ergonomics
The study by Choi [18] was among the first studies in which ergonomics were included in MMALBP.The author proposed a goal programming model to optimize the time and physical workloads at workstations where the minimum deviation from the set goals could be reached using the mathematical modeling approach.The RULA assessment method was used by Barathwaj et al. [19] to calculate the accumulated risk posture.Then the weighted sum of the number of workstations (nw), workload variation, and accumulated risk posture was minimized using a genetic algorithm (GA) for MMALBP in a case study.The smoothness indexes in terms of time and energy expenditure were optimized by Sgarbossa et al. [20] using an exact heuristic algorithm to find a trade-off between these conflicting objectives while examining the effect of a virtual average model (VAM) for the mixed-model.The maximum and average absolute deviations of ergonomic risks as the primary and secondary objectives were optimized by Bautista et al. [21], respectively, given the workstations' cycle time, physical load, and space limitations.The effect of transportation modules and product sequencing modes on awkward postures and reduced recovery times of workers in MMAL was studied by Carrasquillo et al. [22] using a simulation approach.The MMALBP was addressed by Alghazi and Kurz [23] to minimize the number of workers in which the total weighted musculoskeletal risk of all tasks performed by each worker could not exceed a predefined threshold.A real case study of MMALBP was investigated by Tiacci and Mimmi [24] with musculoskeletal risk assessed by the OCRA method while assuming stochastic task times and equipment duplication.A GA was proposed to minimize the musculoskeletal risks besides the total cost of production measured in terms of economic and performance aspects.The MMALBP was addressed by Bautista and Alfaro [25] while reducing the ergonomic risk deviation between workstations and the maximum ergonomic risk level of the assembly line.The authors proposed a greedy randomized adaptive search procedure to solve the problem in an automotive factory.
In the study by Çil et al. [26] MMALBP was addressed with physical collaboration between human workers and robots by proposing a mixed-integer programming (MIP) model to minimize the CT of all models to solve small-size problems.A bee algorithm and an artificial bee colony algorithm were proposed to resolve large-size problems.Rest allowance and fatigue were introduced into MMALBP by Finco et al. [12] where CT was minimized in a real case using a mathematical model and heuristic approaches.Mokhtarzadeh et al. [9] used multi-criteria decision-making to classify the musculoskeletal risks obtained by different methods (e.g., NIOSH, OCRA, JSI, OWAS) into three groups in a mixed-model parallel U-shaped ALBP.The problem was formulated as a MIP, and constraint programming and heuristic approaches were proposed to minimize the nw first and then, the musculoskeletal risk level.Recently, Mura and Dini [27] considered the physical and environmental aspects of MMALBP while improving energy expenditure and noise.They developed a tool based on a GA that minimizes the weighted sum of the nw and the above ergonomic measures in an assembly line of electric scooters while considering job rotation and collaborative robots.

Musculoskeletal risk assessment in ALBP
While assessing the musculoskeletal risk, the intensity, frequency, and duration of abrupt or sustained exposure to force (vibration, repetitive motion, awkward posture) and environmental factors (temperature, noise, and light) are observed [28].Different methods have been used in the literature to measure musculoskeletal risks.The established methods are NIOSH, OWAS, OCRA, RULA, and REBA.The relevant studies focusing on these methods while addressing ALBP are reviewed below.
In the study by Otto and Scholl [11], the musculoskeletal risk assessment measures were incorporated into ALBP formulation and the authors proposed a hierarchical heuristic to solve the related problem.Rajabalipour et al. [29] included a simplified version of OWAS to examine the risk levels of the assembly tasks.The authors formulated the problem as a fuzzy goal programming model to optimize the weighted sum of CT , physical workload balance, and accumulated risk of postures.A GA was developed to solve certain test problems with randomly generated physical workloads.Mossa et al. [30] studied the job rotation scheduling among the assembly operators to balance human workloads and ergonomic risks evaluated by the OCRA method.The problem was formulated as a MIP to maximize the production rate and was successfully applied to solve an industrial case study.Baykasoglu et al. [31] considered the musculoskeletal risks evaluated by the OCRA method for both ALBP and layout optimization in which a rule-based constructive search algorithm was proposed to solve a real-world manufacturing case study.Bortolini et al. [15] evaluated the operators' musculoskeletal risk using the REBA method for component-picking activities in a single-model assembly line.The problem was formulated as an optimization model to minimize the takt time and musculoskeletal risks to solve a kitchen appliance case study.Moussavi et al. [32] combined an in-house risk assessment method based on video capturing with a MIP model to optimize job rotation scheduling and balancing workloads in an automotive assembly line.Akyol and Baykasoğlu [10] introduced human well-being in a variation of ALBP with worker assignment where the OCRA method was utilized for risk assessment.To minimize both CT and musculoskeletal risks, a preemptive goal programming approach was applied in a constructive rule-based random search algorithm.Zhang et al. [33] incorporated musculoskeletal risks evaluated by OCRA into a U-shaped single model ALBP with worker assignment.The problem with the objectives of CT and musculoskeletal risks was for- mulated as a MIP and a restarted iterated Pareto greedy algorithm was proposed to address it.

Ergonomics and digital human modeling in ALBP
Digital human modeling (DHM) belongs to the computerized tools to design and evaluate human-machine systems in a virtual environment.Using DHM and simulation software, digital mannequins can be created with the capability of mimicking most of the postures and critical joints of the human body with high precision [16].The utilization of DHM has been emphasized in different applications, such as human-centric manufacturing [34], vehicle interior design [35], and the development of ergonomic frameworks in highly automated vehicles [36].The relevant studies that implemented DHM in ALBP with an ergonomic perspective are reviewed below.
Longo and Mirabelli [37] proposed a multi-criteria approach based on virtual modeling and simulation software to find an efficient balance of time and ergonomic aspects in a heater production line.Ding and Hon [38] used DHM to simulate manual assembly lines by incorporating ergonomic constraints and finding a solution with minimum posture risks.DHM and motion capture technology were used by Wilhelm et al. [39] to measure RULA scores while balancing human well-being and time loads of workstations.An optimization framework with a holistic system perspective was proposed to minimize the integrated musculoskeletal risk and takt time.In another study by Ozdemir et al. [16], DHM and simulation software were used to measure the risk factors of tasks.A fuzzy mathematical model was proposed to minimize the weighted sum of CT and ergonomic risk imbalance in a real-world problem.

Research gap and motivation of the study
The literature review above shows that most studies, particularly in MMALBP, have used a weighted sum approach to combine ergonomic and economic objective functions into a single objective.In this regard, Sgarbossa et al. [20] developed a multi-objective optimization approach while dealing with time and energy aspects.However, the proposed method was applied to a small size problem.On the other hand, the literature review shows that fewer studies have used the REBA method in ALBP, while this method can measure the risk of musculoskeletal disorders for the entire body considering the weight and coupling of parts and components [28].In addition, the literature review in Section 2.3 shows that DHM and ergonomic simulation applications in ALBP have been limited to the singlemodel ALBP and not MMALBP.Thus, this study aims to fill the above gaps in the literature by (1) formulating an MMALBP with REBA as a MILP model, (2) applying DHM to evaluate the musculoskeletal risk of operators, and (3) developing an enhanced multi-objective optimization algorithm to address the considered MMALBP.

Problem specification
The problem in this study arises in many manufacturing industries with modular assembly processes where the aim is to satisfy the volatile demands of customers with different variants of the same product.However, from the production perspective, this yields a large variety in the assembly operations of the mixed-models.To ensure a smooth assembly of the product variants (mixed-models) in an existing line, it is necessary to balance the mixture of the assembly tasks among the workstations (MMALBP) while minimizing the CT for all variants.However, maximizing efficiency may produce undesirable ergonomic risks for manual assembly operators.Therefore, another objective to be optimized is the distribution of ergonomic loads among the operators at workstations.Thus, the MMALBP considering musculoskeletal risk can be defined as the assignment of a given number of tasks ( ) to a given number of workstations ( = … j nw 1, , ) for given task times for each model (t im ) while considering the musculoskeletal risk levels of tasks.The objectives are to minimize the CT and the imbalance of musculoskeletal risks among operators/workstations.It is assumed that each variant can have different tasks, precedence relationships, processing times, and musculoskeletal risks.With the above information, the problem can be represented by a combined precedence network as illustrated in Fig. 1 for a small 7-task problem with two variants.In this figure, the tasks are shown as nodes.The processing times of models and their musculoskeletal risks are presented above and below the task nodes, respectively, and the directed arrows among task nodes indicate the tasks' precedence relationships.

Digital human modeling
In this study, an embedded DHM tool known as IPS IMMA [40] is applied to assess the musculoskeletal risk of the assembly tasks.Within this DHM tool, three steps are conducted to evaluate the ergonomic risks of the assembly tasks.Fig. 2(a)-(c), illustrates the three steps for calculating the musculoskeletal risks of the tasks which are explained below: • Step 1: The product is modeled in CAD and imported into the DHM tool.Then the related assembly line including all the workstations and resources (e.g., material and equipment) are modeled and imported into the DHM tool as shown for a sample station in Fig. 2(a).Then the manikins representing the average human workers are added to the simulation model.Next, the sequence editor controls the manikins to perform the assembly tasks.Finally, the IPS IMMA engine virtually simulates the assembly tasks.
• Step 2: The IPS IMMA digitally captures the manikins' body postures and movements while performing the assembly tasks in terms of some inputs (e.g., postures position, direction, and angles).Then these data are used to calculate the ergonomic scores, using the original REBA method [41].Fig. 2(b), illustrates the procedure for obtaining the REBA scores for a sample assembly task in which the operator is mounting the upper portion of a workpiece on the lower component.First, the operator's body movements in the lower limbs and neck, including the trunk, legs, and neck, are measured in terms of the joints' angles, as shown by T°, L°, and N° in this figure.The respective angles are then converted into standard scores for each body part provided by the REBA method.For instance, considering the range of neck angles in this sample task (i.e., 30°), the neck score is 2. Accordingly, the related scores for the trunk and leg parts are obtained based on the standard ranges provided by the REBA method.Then a scoring table, as shown in Table (A) in Fig. 2(b), is used to combine all these scores (i.e., neck, trunk, and leg scores) into a posture A score which ranges from 1 to 9. Then posture A score is added by a force/load score, for example, if the task includes some weight/force handling.The resulting score is known as score A. Similarly, the position of the upper limbs is determined based on the joint angles of the upper arm, lower arm, and wrist, as shown by UA°, LA°, and W°, respectively, in Fig. 2(b).These angles are then converted into standard scores for these upper limb parts using the standard positions defined by the original REBA method.For example, the upper arm joint angle of 40° in this example equals a standard score of 2 based on the standard rates provided by the REBA method for this body position.After obtaining the remaining scores for the lower arm and wrist, another scoring table shown as Table (B) for the upper limbs is used to combine all these scores into a posture B score, similar to the posture score A, which ranges from 1 to 9. Then posture score B is added by a coupling score in cases where a handling/coupling is needed while performing the task, particularly, e.g., when handling or gripping an object.The resulting score is known as score B. Having calculated scores A and B, then a combined score C is obtained by matching scores A and score B together using another standard table provided by the REBA method.Finally, the REBA score of the assembly task is calculated by adding an activity score to score C through which the intensity ranges of body postures while performing the assembly task are considered.The obtained REBA score ranges from 1 to 15, representing the risk level of musculoskeletal disorder of an assembly task interpreted by a conversion table ranging from negligible to very high risks, as shown in Fig. 2(b).For the mounting task in this example, the above steps result in a REBA score of 5, which belongs to the medium ergonomic risk level.
• Step 3: the musculoskeletal risks for each task are visualized with further details regarding the body postures and their associated risks.At the end of this step, the REBA scores for each assembly task are calculated as shown for a sample task in

Mathematical model
The considered MMALBP is now formulated as a MILP model.The notations used in this formulation are presented in Table 1.
The objective functions of MMALBP are to minimize (1) the CT , (2) the maximum REBA of all workstations (SE max ), and (3) the sum of REBA over all workstations (SE sum ).In Eqs. ( 1)-( 3), the values of the above objective functions are measured, respectively.

Min OF CT
1

Min OF SE
2 max (2) (5)  ) .( 1) .( 1 (14) Eq. ( 4) ensures that each task is assigned to only one workstation.The precedence relationships among tasks are satisfied using Eq. ( 5).At each workstation and for each model, if there is a precedence relationship among tasks, then the completion times of those tasks are determined by Eq. ( 6).Otherwise, Eqs.(7) or (8) ensure that the tasks with no direct or indirect precedence relationships are scheduled without overlapping each other.Eq. ( 9) determines the CT as the highest completion time of all the tasks, while Eq. ( 10) sets the lower bound of the completion time of each task of each model to be at least equal to its processing time.The musculoskeletal risks of the workstations are defined using Eqs.( 11)- (13).Eq. ( 11) calculates the E jm as the time-weighted average REBA (TWAR) score of workstation j for model m in that the longer the processing time of a task, the larger its weight on the TWAR of the workstation.Using Eq. ( 12), the ergonomic risk of a workstation (SE j ) is set to the maximum ergonomic risk of all models, while Eq. ( 13) determines the maximum musculoskeletal risk of all workstations (SE max ).Finally, Eq. ( 14) defines the decision variables' domain as the binary and real numbers.
The Eq. ( 11 , can be introduced so that Eq. ( 11) can be transformed into a linear form using constraints ( 15)- (19), where = R 15 is the upper bound of E jm .

=
) 1); , , The proposed E-NSGA-II for MMALBP Owing to the successful application of the NSGA-II in different multi-objective optimization problems, in this study, an enhanced NSGA-II (E-NSGA-II) has been proposed to address the MMALBP.The original NSGA-II is based on the Pareto-dominance concept.It was embedded with an elitism-preserving strategy and a niching operator with no parameter that provide a good spread of solutions and converge to the Pareto-optimal solutions [42].The main steps of NSGA-II are as follows: First, a set of solutions are generated randomly as the initial population (P 0 ) with a size of popsize, and their fitness values are evaluated.Then, the solutions are ranked into different Pareto fronts using the non-domination sorting procedure, where front 1 includes all non-dominated solutions (NDS), front 2 includes another NDS without considering front 1, and so on.In addition, the crowding distance sorting procedure orders each solution in comparison with other solutions lying on the same front in terms of a distance measure, i.e., crowding distance (CD), where a higher CD indicates a better dispersion of solutions.Next, a parallel Decision variables: 1 if task i is assigned to workstation j; 0 otherwise population Q G is produced iteratively for each population (P G ) at generation G using the binary tournament selection, crossover, and mutation operators.Finally, P G and Q G are combined, and then the non-domination sorting and CD sorting procedures are reused to select the Pareto optimal solutions identified in a lower front/rank (i.e., 1, 2, …) and/or a higher CD (if fronts/ranks are equal) so that the exact popsize of solutions are transferred to the next generation + P G 1 .The main loop of NSGA-II is repeated until a specific stopping condition (e.g., a maximum number of generations, G max ) is attained.To improve the performance of the NSGA-II for solving MMALBP, in this study, two improvement mechanisms are incorporated into this algorithm, namely; (1) local search and (2) multi-criteria decision-making.The local search mechanism aims to find a set of neighborhood solutions for each generated feasible solution to enhance the diversification capability of the algorithm.On the other hand, the multi-criteria decision-making mechanism attempts to rank the obtained neighborhood solutions by the local search mechanism based on the technique for order preference by similarity to the ideal solution (TOPSIS).The latter feature can improve the intensification capability of the algorithm while finding the NDS.In the following subsections, the elements of the E-NSGA-II are discussed in detail.

Encoding and decoding
The structure of the solutions in the algorithm space is indicated by a string of numbers where the length of the string equals the number of tasks (nt).The values in this string are permutations of 1 to nt.Fig. 3 shows an example of a solution representation for a 7- task example.Here, the number in the i-th task index shows the relative priority of task i to be ordered in the sequence of tasks.
A two-step procedure known as encoding and decoding is followed to generate a feasible solution.The encoding step produces a possible ordering of tasks, or task sequence (TS), based on the given permutation of numbers.In this step, the tasks are sequenced based on their priorities and precedence relationships.The pseudo-code of the encoding procedure is presented as follows: Encoding procedure: 1.
Input: a given permutation of numbers, relationships 2.
Select the task i with the highest priority and free precedence rela- tionship; 5.
Put the selected task i into TS and remove it from TS ; = TS i i , TS = TS -i ; 6. End 7.
Output: TS Next, decoding is performed on the obtained TS from the encoding step.The decoding step attempts to generate a feasible assignment of tasks between a fixed number of workstations (nw) based on the given TS.In the decoding procedure, the tasks in the TS are assigned sequentially to workstations so that the work/time load of the tasks in the first nw-1 workstations does not exceed the theoretical cycle time (TCT ), as presented in Table 1.However, this limit is relaxed for the last workstation, where all remaining tasks are assigned to the nw-th workstation.The pseudo-code of the de- coding procedure is presented as follows: Decoding procedure: 1.
Open the first workstation j = 1, task counter i = 1, workstation time wt When a feasible solution for MMALBP is generated, its fitness function can be evaluated based on the three objectives, namely (1) CT , (2) SE max , and (3) SE sum , as presented in Table 1.

Local search mechanism
The above decoding procedure can generate only one feasible MMALBP solution based on the calculated TCT .To enhance the search quality of NSGA-II and the resulting NDS, a local search mechanism was introduced to find a set of feasible solutions for each solution generated by the decoding procedure.This local search enhances the exploration capability of the multi-objective optimization algorithm by allowing a set of neighborhood searches over the solution space.Using this local search, the decoding step above is repeated for a specific number of iterations (Max IT ) in which a set of feasible solutions can be found.Therefore, instead of having a single , where is a random number generated in the interval (1,2).The above local search stops when a maximum number of iterations (Max IT ) is reached.The pseudo-code of the proposed local search procedure is presented as follows: Local search procedure: 1.
For IT = 1: Perform the decoding procedure to find a feasible assignment of tasks to workstations 5.
Evaluate the fitness function (CT , SEmax , and SEsum) for the obtained feasible solution 6. End 7.
Output: A set of feasible solutions

Multi-criteria decision-making
To ensure that the best-found solutions by the above local search are maintained during the search, another mechanism was embedded in the algorithm.This mechanism aims to evaluate all the generated solutions by the local search step and preserve the one with the highest quality based on the considered objectives.Therefore, a multi-criteria decision-making (MCDM) approach based on TOPSIS by Hwang et al. [43] was implemented in the NSGA-II to rank the set of feasible solutions obtained by the local search mechanism.The TOPSIS assumes the set of feasible solutions by the local search as the alternatives, with the size of max IT , and the values of their fitness function based on CT , SE sum and SE max as the decision criteria.It is worth mentioning that the original TOPSIS method sorts the feasible solutions based on a similarity measure calculated based on the distances to the calculated negative and positive ideal solutions.The higher this similarity measure, ranging between (0,1), the better the solution.Next, the solution with the highest similarity measure is selected as the best-found feasible solution among all the solutions obtained by the local search.

Generation of new solutions
To generate new solutions, the tournament selection, crossover, and mutation operators were utilized.It is worth mentioning that the non-domination sorting and CD sorting steps are not presented here for brevity purposes, and interested readers are referred to Deb et al. [42] for further details.The tournament selection ensures that the solutions with lower ranks or higher CDs are more likely to be selected as the parents for the recombination process.First, two solutions are chosen randomly from the current population.Then, the chosen solutions are compared in terms of their ranks.If these have equal ranks, their CDs are compared with each other.The dominating solution with a lower rank or higher CD (if ranks are similar) is then selected as the first parent (Par1).The above procedure is repeated to determine the second parent (Par2).Fig. 4 illustrates a sample tournament selection process.
Then, the recombination process of the selected parents is performed by the crossover operator.Based on the solution representation, the parent matching crossover shown in Fig. 5 was applied to recombine two selected parents [44].
After the crossover, the mutation operator transmutes the task priorities in a few offspring to maintain a higher diversity among the obtained solutions.Fig. 6 illustrates the applied mutation operator.
Based on the above elements, the flowchart of the E-NSGA-II proposed for MMALBP is presented in Fig. 7.

Computational results
In this section, the performance of the proposed E-NSGA-II is benchmarked against the MILP and other solution methods.First, the MILP and the E-NSGA-II are compared based on some small-sized test problems, where the Pareto front solutions by the original epsilon-constraint method in Carvalho and Ribeiro [45] are compared with and the E-NSGA-II.Second, the performance of the E-NSGA-II is evaluated against the conventional version of NSGA-II (without enhancement) by Deb et al. [42] and the adapted version of MOGA by Nearchou [46], while addressing a case study and a few test problems taken from the literature.
The case study information and the data were taken from a major automotive company.Regarding the test problems, the basic information was based on the benchmark problems accessible at [47].The rest of the data that were not included in the original data set, e.g., task times for additional variants (mixed-model) and REBA scores, were randomly generated.All the information about the considered problems is available in the supplementary data.
The MILP and epsilon-constraint method were coded in the "General Algebraic Modeling System" software and solved by the CPLEX solver.The multi-objective optimization algorithms were coded in MATLAB R2023a and ran on a Core i9 PC with a 3.10 GHz processor and 64 GB of RAM.
A statistical analysis of the results was performed to test the hypothesis about the superiority of the algorithms in terms of performance metrics.Three following metrics were calculated to measure the quality of NDS obtained by the multi-objective optimization algorithms [48]: • Number of non-dominated solutions ( NDS | |): This indicator mea- sures the quantity of the non-dominated solutions in the obtained Pareto front.The higher the |NDS|, the better the approach.
• Hyper volume (HV): This indicator measures the quality of the NDS in terms diversity and uniformity when the real Pareto front is unknown.The higher the HV, the better the approach.
• Spacing metric (SM): This indicator measures the non-uniformity of the distribution of NDS over the Pareto front.The lower the SM, the more uniform the solutions are dispersed over the Pareto front.The parameters of the algorithms were calibrated by the Taguchi Method [49][50][51] using the orthogonal array L 9 (3 4 ) and the Signal-to- Noise (S/N) where = S N NDS HV / 10 log(| |/ ).As a result, the levels of parameters for the E-NSGA-II and NSGA-II were set as popsize = 100, crossover rate = 0.99, mutation rate = 0.1, and G max = 1000.Like- wise, the MOGA's parameters were set to popsize = 100, crossover rate = 0.8, mutation rate = 0.2, elitism rate = 0.2, and G max = 1000.

Small-sized test problem
The epsilon-constraint method was used to find the Pareto-optimal solutions in some small-sized test problems from the literature.The problems' specifications are shown in the first five columns of Table 2.While using the epsilon-constraint method, the CT was considered as the main objective whereas the other objectives, namely SE max and SE sum were dealt as the additional constraints delimited by the related epsilon bounds.Some pilot studies were conducted to select three epsilon bounds for the latter two objectives based on their minimum, average, and maximum values.Accordingly, the number of Pareto-optimal solutions found by the epsilon-constraint method for each instance can be at most 9 NDS.In Table 2, the resulting performance metrics of Pareto-optimal solutions obtained by the epsilon-constraint and E-NSGA-II in terms of |NDS|, HV, and SM are presented.
Considering |NDS| and SM, the epsilon-constraint is dominated by E-NSGA-II in all instances.In terms of HV, E-NSGA-II dominates epsilon-constraint, except in instance 8. Therefore, is the dominant approach due to the quantity and quality of NDS considering the performance metrics.Fig. 8 shows the box plots of objective values of NDS found by epsilon-constraint and E-NSGA-II over all considered instances.According to this figure, the lower bounds of CT , SE max , and SE sum obtained by E-NSGA-II are either less than (SE max and SE sum ) or equal (CT ) to the values obtained by the epsilonconstraint.Moreover, the connecting lines the means of each objective by each method show that the means of objectives by E-NSGA-II is mostly less than the related means obtained by the epsilon-constraint method.Thus, it can be concluded that the solutions obtained by the E-NSGA-II lie close to the exact NDS obtained by the epsilon-constraint method while having better performance metrics.

Case study and test problems
The efficiency and effectiveness of the proposed E-NSGA-II were also evaluated by solving a case study and some test problems.
The case study has two models of a modular product from the automotive industry assembled on a straight line.There are 29 tasks with different processing times for each model.Given their precedence relationships, the tasks can be performed at either three, four, or five workstations.The information about the assembly process of these variants (including the precedence relationships  and processing times) is provided by company experts.Moreover, the musculoskeletal risks of the assembly tasks are digitally calculated based on a DHM tool using the REBA method, as discussed in Section 3.2.The input data for the case study is presented in the supplementary data.Additionally, a set of test problems taken from the literature were solved.The specifications of the considered instances for both case study and test problems are shown in the first five columns of Table 3 in terms of the problem, nt, nm, and nw, repectively.In the rest of columns, the performance of E-NSGA-II is compared with the other two algorithms, namely NSGA-II and MOGA.The results yielded by the multi-objective optimization algorithm are reported in this table in terms of |NDS|, HV, and SM.The CPU times of the algorithms ranged between 160 and 1600 s Considering the results in Table 3, the following observations can be made.In the case study, which includes three instances, E-NSGA-II dominates MOGA and NSGA-II, considering |NDS| and SM.In the HV metric, E-NSGA-II overcomes other algorithms except in instance 3, where NSGA-II was better.In the Thomopoulos problem with 9 instances, E-NSGA-II is superior in terms of both |NDS| and HV with only one exception in instance 8 where NSGA-II is better.Considering SM, in instances 4, 7, 8, 10, and 11, E-NSGA-II was better, while in instances 5, 6, 9, and 12, NSGA-II was superior.Considering the Kim test problem with 9 instances, E-NSGA-II was superior to  To better demonstrate the results reported in Table 3, Fig. 9 graphically illustrates the performance metrics of all algorithms over all 30 instances using line diagrams in logarithmic scale.According to this figure, the |NDS| obtained by E-NSGA-II is always better than the other algorithms.Similarly, in the HV metric, in 27 instances, E-NSGA-II was superior, while only in three instances (3,8 and 23) NSGA-II was better.Finally, regarding the SM, in 22 instances, E-NSGA-II was dominant, while in the remaining 8 instances, NSGA-II was better.The above results indicate that the E-NSGA-II outperforms MOGA and NSGA-II, particularly in terms of the NDS | | and HV over the considered instances.Considering the SM, E-NSGA-II is superior to other algorithms in the majority of the instances followed by NSGA-II.
In addition, some statistical tests were performed to verify the above observations about the superiority of E-NSGA-II on the performance metrics.The one-way ANOVA test was applied to the 30 instances as the sample data while assuming the algorithms as the groups.The null hypothesis tested for each metric was whether the  Furthermore, the analysis of the above results shows that among the assembly tasks in the case study, task 12 (tightening screw in the cap) in Model1 and task 13 (mounting steering wheel) in Model2 have the highest musculoskeletal risks measured in terms of REBA.Therefore, further studies should be conducted from the ergonomic perspectives on these tasks to find how their musculoskeletal risks could be moderated, e.g., by improving the operator's body postures, the workpiece repositioning, etc.

Managerial insights
The proposed E-NSGA-II provides the decision-makers with a wide range of opportunities to improve the line efficiency and wellbeing of operators in human-centric MMALBP.However, further studies should be conducted to assess which solution is more advantageous regarding the company's finances.In this regard, the ergonomic significance of different solutions can be evaluated in terms of the employees' well-being (e.g., absence, sickness) and their economic effects in the short and long term.Then, the most efficient alternative, where a trade-off between ergonomic risks and line efficiency occurs, can be achieved from the company's perspective by conducting a cost-and-benefit analysis of the results.
The application of the DHM tool in this study could enrich the risk assessment by digitalizing the manual steps of calculating musculoskeletal risks in the body postures of the operator while performing the assembly tasks.In addition, the anthropometric properties of operators, age, gender, etc., can also be considered while using DHM in measuring ergonomic risks.In this study, the musculoskeletal risks for average male operators are calculated by the DHM tool in the case study.Thus, it is recommended to conduct comparable studies to investigate how gender might affect ergonomic risks and the results of the MMALBP.

Conclusions and future research
The emergence of Industry 5.0 has emphasized the necessity of human-centric workplaces given today's aging labor resources.Alternatively, mass customization trends have urged many industries to maintain higher agility and more efficiency in their manufacturing processes.Consequently, it is necessary to consider the ergonomic aspects while focusing on production efficiency.This study concentrates simultaneously on ergonomic and efficiency perspectives while balancing mixed-model assembly lines known as MMALBP.From the ergonomic perspective, this study attempted to implement the DHM while measuring the musculoskeletal risks of operators through a risk assessment method known as REBA.This was accomplished by: (1) simulating the assembly tasks using a DHM tool, (2) calculating the musculoskeletal risks digitally, and (3) visualizing the REBA scores for the whole body.Three objectives were considered to optimize the time and ergonomic loads among workstations: (1) cycle time, (2) maximum ergonomic risk, and (3) the sum of ergonomic risks.The resulting problem was formulated as a mixed-integer linear programming model.A multi-objective optimization algorithm based on NSGA-II was developed to address the problem.The proposed approach (E-NSGA-II) incorporates a local search mechanism that generates neighborhood solutions to improve the search capability of the proposed algorithm.In addition, an MCDM approach based on TOPSIS was embedded in the algorithm to ensure that quality solutions are selected as the population members.The computational results indicated that the E-NSGA-II can find Pareto front solutions close to the exact front solution found by the epsilon-constraint method while solving small-sized problems.The proposed E-NSGA-II was benchmarked against other algorithms, namely MOGA, and NSGA-II, in terms of the quantity and quality of the obtained NDS while solving a case study from the automotive industry as well as test problems from the literature.The significance of E-NSGA-II results, while compared with MOGA and NSGA-II, was tested through one-way ANOVA.The analysis of results showed that when all three considered objectives are simultaneously considered in evaluating the obtained Pareto front solutions, a smooth distribution of both time and ergonomic measures in MMALBP can be attained in the case study.
For future research direction, finding a trade-off between cycle time, ergonomic risks, and the environmental aspects while addressing MMALBP is deemed an important future research direction.This can be further studied by expanding the current problem while considering the sustainability and energy-efficiency aspects.Moreover, developing efficient multi-objective optimization approaches to deal with such optimization problems is another exciting research topic.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Three steps of using DHM to calculate ergonomic scores.

Fig. 2 (
Fig.2(c).In this step, the REBA values are ranged from negligible to very high risk with different highlight colors for each body part of the operator.

TCT
, a set of cycle times at different iterations of the local search step is defined, as shown by TCT IT .The initial value of TCT IT at the first iteration (TCT 1 ) is set to TCT .Then, the values of TCT IT at iteration IT change dynamically using the equation:

Fig. 8 .
Fig. 8. Box plots of objectives values in NDS found by epsilon-constraint and E-NSGA-II.

Fig.
Fig. Line diagrams for performance metrics of the considered algorithms in logarithmic scale.

Fig. 12 .
Fig. 12. Selected solutions for each scenario in the case study.
) above is non-linear because E jm and X ij are multiplied as decision variables, where ij.Therefore, a new auxiliary variable, i.e., + Z ijm

Table 1
Notations for the formulation of MMALBP.

Table 2
Comparison of epsilon-constraint and E-NSGA-II in small-sized test problems.

Table 3
Comparison of multi-objective algorithms in the case study and test problems.

Table 4
ANOVA one-way test on the performance metrics' means of different algorithms.