Introduction

Fatigue within joints during repetitive lifting tasks can substantially diminish the available muscle force capacity. It is a complex process arising from the interplay of multiple more basic physiological and neurological processes. Therefore, singling out the specific contribution of any individual phenomenon becomes notably intricate [1,2], necessitating the creation of algorithms capable of predicting subject-specific repetitive lifting strategies, navigating varying weights with each lift. In recent decades, researchers have developed different predictive models for simulating lifting motions based on both kinematics and physics [3]. Furthermore, considerable efforts have been made in the advancement of optimization-based algorithms to predict lifting motions [4][5][6][7][8][9][10]. These investigations utilize gradient-based sequential quadratic programing (SQP) methodologies to solve constrained nonlinear optimization problems for motion prediction. While these reported predictions can exhibit noteworthy precision for lifting without factoring in fatigue effects, a principal limitation emerges in accurately foreseeing fatigue progression within repetitive tasks. Incorporating the influence of fatigue emerges as a required step to enhance the accuracy of these simulated tasks.

Manual material handling (MMH) tasks constitute a critical aspect of industries ranging from manufacturing to healthcare. These tasks, however, bring with them a potential challenge of considering joint fatigue [1]. In workplaces, where weights and lifting scenarios vary widely, the imperative for a robust risk assessment tool is required. Such a tool must facilitate the routine determination of optimal motion, maximal weight capacities, and endurance times, with the objective of preventing work-related injuries. Though the development of physical prototypes remains essential for devising equipment and spatial arrangements, the realm of predictive biomechanical human simulations offers an efficient tool for interacting with new environments. The convergence of the optimization-based dynamic modeling and simulation method [6] with a muscle fatigue and recovery model [11,12] stands out as a promising approach for predicting fatigue in weight-lifting tasks, particularly involving the transportation of boxes.

While fatigue models have evolved through analytical or empirical roots [13], their applicability to lifting tasks remains a challenge [14][15][16]. Addressing the complexities intrinsic to such tasks calls for the inclusion of decay terms into existing analytical models, especially at the single muscle level [14]. A potential solution emerges using the motor unit models in the form of the three-and four-compartment controller fatigue models, encompassing resting, active, and fatigued compartments [11,12,17,18]. Noteworthy for their consistency with muscle physiology and joint biomechanics, these models offer a pathway to address the fatigue due to complex biomechanical tasks.

Pereira et al. [19] applied a physics-based fatigue model to a single elbow joint. Barman et al. [15] applied the three-compartment controller (3CC) fatigue model to lifting posture prediction to maximize the endurance time. Michaud et al. [16] used a 3CC fatigue model to optimize the repetitive elbow lifting motion. This study compares the repetitive lifting prediction between the 3CC and fourcompartment controller with augmented recovery (4CCr) muscle fatigue models through a 2D 10-degrees-of-freedom (DOF) human model. The central aim of this study is to elucidate an inverse dynamics optimization formulation to predict fatigue progression and the optimal lifting strategy through repetitive lifting with a given 20 kg box load. The repetitive lifting task is treated as a nonlinear programing (NLP) optimization problem, driven by the objective of minimizing dynamic effort and joint acceleration using the SQP algorithm. The novel contribution of this research lies in the development of an optimization-based motion prediction formulation and the subsequent integration of muscle fatigue models with the developed motion prediction formulation. The 4CCr model [17,18], an extension of the conventional 3CC fatigue model, accounts for both central and peripheral fatigue mechanisms. This advancement enhances the accuracy of fatigue predictions in joint space, thereby enriching the predictive capabilities of the overall model.

The 3CC muscle fatigue model has been validated only for isometric (static) tasks in literature although it can theoretically be applied to dynamic tasks [11]. In this study, we integrated the recently developed 4CCr fatigue model, which considers the effect of joint angular velocity on both central and peripheral fatigue into the motion prediction formulation. For comparison, we implemented both the 3CC and the 4CCr fatigue models to repetitive lifting tasks to validate their robustness for dynamic tasks. In the literature, muscle-based fatigue models have been reported but they can only be used for individual muscles or a model with a few muscles. However, both the 3CC and the 4CCr fatigue models can be applied to either individual muscles, motor units, or joints. If applied to muscles or motor units, the activation, resting, and fatigue compartments are defined as the respective fractions of total muscle forces, or of total motor unit forces. In contrast, if applied to joints, the activation, resting, and fatigue compartment sizes are defined as fractions of the net joint torque. Then, the fatigue governing equations are applied to the three compartments for the 3CC model and the four compartments for the 4CCr model. This study uses a simplified 2D whole-body model and applies fatigue to the major joints (in joint space) rather than to individual muscles throughout the body to reduce computational cost for dynamic motion prediction. A critical consideration for object lifting is stability, to address this, this work uses the ZMP criterion as a constraint. Examples of other stability metrics used in the applications of object lifting or load carriage include balanced state basin [20], ground projection of center of mass [21], and dynamic gait measure [22].

Human Model and Fatigue Models

2.1 Two Dimensional Human Skeletal Model and Equations of Motion. First, the procedure for establishing a digital human model based on authentic anthropometric data is delineated. This process involves the development of a recursive approach for computing the kinematic and dynamic attributes of the human model. Then the motion prediction formulation incorporates the muscle fatigue model: 3CC or 4CCr model, facilitating an in-depth assessment of fatigue progression and enabling the anticipation of the repetitions of lifting prior to the onset of task failure.

To explore the progression of fatigue arising from repetitive lifting, a human model is characterized within joint space, comprising 10 DOFs in Fig. 1. The three initial DOFs-denoted as q 1 , q 2 , q 3 -encompass overall translation along the forward (z) and upward (y) axes, as well as global clockwise rotation (b). These translational and rotational motions correspond to the collective movements of the entire body around the hip joint, within the context of the global reference frame. The seven other DOFs encapsulate rotational actions specific to individual joints, encompassing spinal flexion (q 4 ), shoulder flexion (q 5 ), elbow extension (q 6 ), hip extension (q 7 ), knee flexion (q 8 ), ankle plantar flexion (q 9 ), and metatarsophalangeal flexion (q 10 ). The anthropometric data employed for the human model is drawn from GEBOD, a tool for anthropometric data regression [23]. Notably, the third (q 3 ), fourth (q 4 ), and seventh (q 7 ) DOFs coincide, indicating the alignment of these three joints. These DOFs are established along the local z-direction, with local coordinates established at the end of each joint, in adherence with the Denavit-Hartenberg (DH) methodology [24].

In this 2D model, the global coordinates constitute the parent branch, subsequently diverging into two child branches symbolizing the spine-arm and leg branches. This model is designed for an individual at a height of 1.7 m and weighing 67.3 kg. Note that mass, COM, and inertia of link 1 (L 1 ), are for the combined spine and head system, treated as a fused rigid body. Mass, COM, and inertia of link 4, L 4 are for the combined femur and pelvis system, also treated as a fused rigid body. Except for the spine link, other links combine the right and left limbs in the sagittal plane. In addition, the inertia terms Ixy ¼ Ixz ¼ Iyz ¼ 0. Also, the COM and inertia are measured in DH local coordinates from the end of each link.

The kinematic equations of the 2D model are deduced from the DH methodology, a technique enabling the conversion of local coordinates into global coordinates using transformation matrices. The DH parameters, denoted as h, d, a, and a, correspond, respectively, to the rotation angle around the preceding z-axis, the displacement along the preceding z-axis, the displacement along the current x-axis, and the rotation around the current x-axis. The DH transformation matrix is illustrated in Eq. (1), while the DH table specific to the 2D human model is given by Table 1 iÀ1 T i ¼

where the top-left 3 � 3 submatrix denotes the rotation matrix and the top-right 3 � 1 submatrix denotes the translation vector. This study employs a methodology encompassing recursive kinematics and Lagrangian dynamics for the human model including forward kinematics and backward dynamics. The former enables the propagation of motion from the origin to the end effectors, while the latter facilitates the transference of forces from (2)-(6), with s i denoting the torque exerted at the ith joint. Within the torque expression delineated by Eq. (2), the first component corresponds to the composite of inertia and Coriolis torque, the second element encompasses the torque attributable to gravitational loading, the third term signifies the torque arising from external forces, and the fourth term captures the torque resultant from external moments

I zz 2 I xy I xz m i x iCOM I xy I xx -I yy þ I zz 2 I yz m i y iCOM I xz I yz I xx þ I yy -I zz 2 m i z iCOM m i x iCOM m i y iCOM m i z iCOM m i 2 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 5 (7)

where tr � ð Þ is the trace of a matrix, A i , C i are, respectively, the global position and acceleration transformation matrices, I i is the inertia matrix for link i given by Eq. ( 7), D i is the recursive inertia and Coriolis matrix, E i is the recursive vector for gravity torque calculation, F i is the recursive vector for external force-torque calculation, G i is the recursive vector for external moment-torque calculation, g is the gravity vector, m i is the mass of link i, r i is the COM of link i in the ith local frame, f k ¼ 0 f ky f kz 0 � � T is the external force applied on link k, r k is the position of the external force in the kth local frame, h k ¼ h x 0 0 0 � � T is the external moment applied on link k, z 0 ¼ 0 0 1 0 � � T for a revolute joint, z 0 ¼ 0 0 0 0 � � T for a prismatic joint, and finally, d ik is the Kronecker delta given by Eq. ( 8).The starting conditions are

For a comprehensive derivation of the equations of motion and their sensitivity concerning state variables, refer to Refs. [25] and [26].

The computation of ground reaction forces (GRF) employs a force balance model associated with foot-ground interaction. The active forces on the body are inertia, Coriolis forces, gravitational influence, and the external load of the box. The principle of active-passive balance, introduced in earlier studies [5,27] is used in this model: in the first step, the global forces and moments are calculated from given state variables without the GRFs; and in the second step, the calculated global forces and moments are transferred to center of pressure and further applied to the metatarsal joints as external forces and moments. Through the implementation of this contact model, the resultant global forces and torques will converge to zero, given the fulfillment of the balance condition during the optimization process.

Three-Compartment Controller Fatigue

Model. The 3CC fatigue model comprises three fundamental components: the resting, active, and fatigued compartments, where each joint motor unit is assigned to one of the three compartments. Illustrated in Fig. 2, the system's schematic demonstrates the interplay between these compartments, which is governed mathematically by three first-order differential equations as shown in the following equations:

where M Ai t ð Þ denotes the active joint torque for the ith joint at time instant t, M Fi t ð Þ represents the fatigued joint torque for the ith joint at time instant t and M Ri t ð Þ denotes the resting joint torque for the ith joint at time instant t. Additionally, R i stands for the recovery coefficient for the ith joint, F i denotes the fatigue coefficient for the ith joint, and

ð Þ is a bidirectional, time-varying torque activation-deactivation drive for the ith joint at time instant t, acting as a link between M Ai t ð Þ and M Ri t ð Þ. The value and direction of C i t ð Þ is determined by the value of TL i t ð Þ which represents the target load for each joint at time instant t which is calculated using the following equation:

, for positive torque direction

where s U i denotes the maximum allowable torque and s L i denotes the lower limit for joint torque, for the ith physical joint of interest. If the active compartment size fails to meet target load, the activation/

Fig. 2 Schematic of the 3CC fatigue model depicting the flow of motor units between the compartments

Table 1 Denavit-Hartenberg table for the 2D human model

Hind foot, L 6 0 q 10 0 0 Fore foot, L 7 0 deactivation drive will recruit more motor units from the resting compartment to the active compartment so to meet the target load.

Conversely, if more motor units are assigned to the active compartment than the target load demands, the drive will reassign some of the motor units from the active compartments back to the resting compartment. The values of C i t ð Þ adhere to the conditions detailed in the following equations:

where L Di and L Ri denote the force development and the relaxation factors, respectively. Both are merely tracking parameters with minimal effect on endurance time, so arbitrary values close to literature [28] are assigned to these parameters, i.e., L Di ¼ 5 and L Ri ¼ 5. Notably, the values of the fatigue and recovery coefficients for the five pertinent physical joints are obtained from literature [28], focusing solely on the study of fatigue progression in these joints.

The coefficients used for the hip joint are the recommended generalized joint coefficients from the same reference.

Four-Compartment Controller With Augmented Recovery Fatigue

Model. The 4CCr model represents a sophisticated advancement in modeling neuromuscular fatigue dynamics at the level of individual joints [17,18]. In essence, this model extends the conventional 3CC fatigue model by introducing a new fatigued compartment allowing for a distinction to be made between central and peripheral fatigue mechanisms. Fundamentally, the 4CCr model divides the fatigued compartment of the traditional 3CC model into two discrete components: peripherally fatigued M Fp t ð Þ, and centrally fatigued M Fc t ð Þ compartments. This partitioning allows for a more detailed analysis of fatigue processes, distinguishing between central mechanisms originating in the central nervous system (CNS) and peripheral mechanisms, which encompass a broad array of factors occurring distally to the neuromuscular joint. These factors may include but are not limited to, processes such as muscle action potential propagation, excitation-contraction coupling (ECC), and chemical alterations within the contractile elements of muscle tissue.

Figure 3 shows the schematic of this system and the flow between the compartments which is governed mathematically by the four first-order differential equations as in the following equations:

where M Ai t ð Þ is the active compartment size for the ith joint at time instant t, M Fpi t ð Þ is the peripherally fatigued compartment size for the ith joint at time instant t, M Fci t ð Þ is the centrally fatigued compartment size for the ith joint at time instant t, M Ri t ð Þ is the resting compartment size for the ith joint at time instant t, R pi is the recovery coefficient for the peripherally fatigued compartment of the ith joint, R ci is the recovery coefficient for the centrally fatigued compartment of the ith joint, F pi t ð Þ is the fatigue coefficient due to the contribution of the peripherally fatigued compartment of the ith joint at time instant t, F ci t ð Þ is the fatigue coefficient due to the contribution of the centrally fatigued compartment of the ith joint at time instant t, and C i t ð Þ is a bidirectional, time-varying torque activation-de-activation drive for the ith joint at time instant t which relates M Ai t ð Þ and M Ri t ð Þ. The fatigue and recovery coefficient values for all five physical joints are reported in Tables 2 and 3 [17,18] and can be described by the Eqs. ( 20)- (24). Due to availability of experimental data for only shoulder, hip, and knee joint coefficients for the 4CCr model, coefficients from the 3CC model [28] are used for the spine and elbow as F ci0 and R ci0 . The values of C i t ð Þ satisfy the conditions given by Eq. ( 20):

where L denotes the force development/relaxation factor (value set to 40) that ensures the model reacts to changes in target load quickly. Table 3 Fatigue and recovery coefficients (For TL i (t) > 0)

Spine (q 4 ) 0.0 0.0 0.00755 0.00075 0.0 Shoulder (q 5 ) 0.0072488 0.0000082 0.01820 0.00168 0.0233 Elbow (q 6 ) 0.0 0.0 0.00912 0.00094 0.0 Hip (q 7 ) 0.017664 0.000032 0.00970 0.00091 0.0319 Knee (q 8 ) 0.004906 0.000022 0.01500 0.00149 0.0147

Table 2 Fatigue and recovery coefficients (For TL i (t) 5 0)

Spine (q 4 ) 0.0 0.0 0.00755 0.00075 0.0 Shoulder (q 5 ) 0.014532 0.000012 0.01820 0.00168 0.0086 Elbow (q 6 ) 0.0 0.0 0.00912 0.00094 0.0 Hip (q 7 ) 0.00848 0.00001 0.00970 0.00091 0.0254 Knee (q 8 ) 0.00927 0.00001 0.01500 0.00149 0.0182 r is a scaling factor that helps better prediction during dynamic tasks (i.e., box lifting), and its value is set to 15.

Repetitive Lifting Motion Prediction Formulation With Muscle Fatigue

3.1 Optimization Formulation With the Three-Compartment Controller Fatigue Model. In this section, the constrained optimization formulation for a single cycle of repetitive tasks is first presented. First, the elements of the formulation, including the design variables, their initial guesses, the objective functions, and the constraints, are discussed. It is important to note that each cycle of the repetitive task consists of two phases: the lifting phase, where the box is lifted to a designated location, and the setting down phase, where the box is returned to its initial position. The time spent in each phase is experimentally observed to be equal, representing half of the total lifting time. Additional repetitive lifting constraints are outlined and added for the optimization formulation to consider multiple cycles.

Considering that the optimization simulation unfolds within the joint space, the joint angle profiles, denoted as q i t ð Þ, are the first design variables. Due to the progression of fatigue within a repetitive task, the design variables also encompass additional profiles originating from the 3CC fatigue model. These encompass the resting compartment size profiles, M Ri t ð Þ, the active compartment size profiles, M Ai t ð Þ, and the fatigued compartment size profiles, M Fi t ð Þ, all examined as part of the design variables. The continuous time profiles of the design variables can be transformed into discrete control points using B-spline interpolation, given by Eq. (25). With this representation, the control points become the design variables for the parameterized optimization problem. Quartic B-splines represent a numerical interpolation technique replete with critical attributes including continuity, differentiability, and localized control. These attributes, particularly differentiability and localized control, render B-splines an apt choice for capturing trajectories of joint angles, resting compartment size, active compartment size, and fatigued compartment size. These trajectories necessitate both smoothness and adaptability, qualities which are aptly provided by B-splines. The attainment of interpolation continuity, spanning both first and second derivatives within time intervals and at interpolation nodes, is achieved through spline interpolation

where m is the total number of control points, s ¼ fs 0 ,…., s l g is the knot vector, P ¼ fP 1 ,……, P m g is the control point vector and N i s, t ð Þ is the basis function. The joint angle profiles and their derivatives, q i t ð Þ, _ q i t ð Þ, and € q i t ð Þ, are functions of t and P; therefore torque s ¼ s t, P ð Þ is an explicit function of the knot vector and control points from the equations of motion. Thus, the derivatives of torque with respect to control points can be computed using the chain rule as given by the following equation:

In this investigation, multiplicity is incorporated at the termini within the knot vector. In the context of the quartic B-spline, this multiplicity characteristic ensures that the initial and terminal values of joint angles and compartment sizes for a DOF precisely match those associated with the initial and terminal control point values.

The count of knots is contingent upon both the number of control points and the degree of the B-spline. A greater number of control points results in more flexible curves, but it concurrently escalates computational demands due to their role as design variables. To ascertain the optimal count of control points, practical numerical assessments are required. In this study, we opted for 10 control points for each quartic B-spline, thereby leading to 13 knots [10 (number of control points) þ 4 (quartic degree) -1]. Consequently, a total of 280 design variables are used in the optimization [10 (number of DOF) * 10 (number of control points) þ 3 (number of fatigue model compartments) * 6 (number of physical joints of interest) * 10 (number of control points)]. These design variable ensembles can be systematically aggregated and expressed through the following equation:

Initial guesses for all the joint angles and active and fatigued compartment sizes are set to 0, but the resting compartment sizes are set to 1, reflecting the optimization starting with the joint at rest. The optimization problem is solved in SNOPT [29], which employs a sequential quadratic programing (SQP) solver. The next subsection discusses the mathematical formulation of the objective function to study the fatigue progression during each cycle of the repetitive task.

The goal of this study is to minimize the dynamic effort and joint angular accelerations of all the physical joints, mathematically given by the expression in the following equation:

where i denotes the index of the physical joints considered in the motion study, s i denotes the i-th joint torque, s U i and s L i are upper and lower limits of the joint torques, respectively, w s denotes the weight associated with the joint torques, and w a denotes the weight associated with local angular accelerations. Note that all physical joints are present in the objective function expression-not just the five physical joints for which fatigue is considered. This ensures that the lifting and setting down motions are predicted correctly. The weight associated with the joint torques is set to 0.5 and that associated with the local angular accelerations is set to 1.0 for all the physical joints. Some form of numerical testing for Pareto optimization [30] is required to find these joint weights and these ensure the lifting and setting down motions are smooth. The optimization formulation for each cycle of the repetitive task encompasses both time-dependent and time-independent constraints: time-dependent constraints are computed sequentially throughout the optimization process at each discretized time point. Conversely, the optimization process evaluates time-independent constraints at specific designated times. The constraints imposed upon the design variables are discussed in Table 4, where the designated times for the time-independent constraints are also recorded. For the joint angle time-independent constraints, the experimental joint angle values are specified in Table 5.

The optimization formulation for a single cycle of the repetitive task discussed above is repeated until SNOPT fails to find an optimal solution for the current cycle. The additional time-independent constraints are given by Eqs. ( 29)-( 35) are imposed to ensure a smooth transition between consecutive cycles. Equation ( 29) is imposed to ensure that there is no fatigue present for any of the five physical joints of interest at the beginning of the first lifting cycle

While simulating the optimization process for the repetitive task, starting with the second cycle, the time-independent constraints given in Eqs. ( 30)-(35) are imposed. For example, the final compartment sizes from the optimal solution of the first cycle become the initial compartment sizes for the second cycle, the final compartment sizes from the optimal solution of the second cycle become the initial compartment sizes for the third cycle, and so on, as shown in the following equations:

where

resting, active and fatigued compartment sizes at initial and final times of the kth cycle, respectively, with k ¼ 2, 3, 4, …, n iter : A small tolerance of 60:02 is considered for the constraints in Eqs. ( 30) and ( 31) and a small tolerance of 60:001 is considered for the constraint in Eq. (32).

Also, to ensure a general increase in fatigue compartment size across all five physical joints of interest, the derivatives of the final compartment sizes from the optimal result of the previous cycle are also set equal to the derivatives of the initial compartment sizes of the current cycle as shown in the following equations:

where

the derivatives of the resting, active and fatigued compartment sizes at the initial and final time instants of the kth cycle, respectively. Note that a small tolerance of 60:02 is considered for the constraints in Eqs. ( 33)-(35) to ensure numerical convergence. Also, additional constraints are given by Eqs. (36)

Table 4 Constraints considered for the optimization formulation (single cycle)

Constraint name Equation Explanation

Physical constraints

Joint angle limits

limit s U i upper limit Foot contacting position P feet x, t ð Þ ¼ P specified feet heel (z ¼ -0:05 m), ankle (z ¼ 0:0 m), and toe (z ¼ 0:235 m) Zero-moment point (ZMP) P ZMP x, t ð Þ � FSP ZMP stays inside foot support polygon (FSP) Box forward position z wrist x, t ð Þz knee x, t ð Þ � 0 z wrist x, t ð Þ wrist location, z knee x, t ð Þ the knee location Collision avoidance d human x, t ð Þ � r human þ dep 2 d human distance between hand and circle center filling on human body, r human radius of the circle, dep depth of the box Monotonically changing wrist position

Initial and final box locations

Initial, midtime and final joint velocity conditions

Initial, quarter-time, midtime, three-quarter-time, and final joint angles

Fatigue-based constraints for 3CC fatigue model

Governing equations of the fatigue model

Three governing equations of the 3CC fatigue model set up as equality constraints Unit summation

Any motor units assigned to the compartments of a joint must add up to 1

Residual capacity

Motor units available from the resting and active compartment sizes must meet the target load for sustained task

Table 5 Experimental joint angles (in degrees) at quarter cycle intervals Time 0 T 4 T 2 3T 4 T DOF Lifting phase (deg) Setting down phase (deg) Spine (q 4 ) 76.29 70.19 17.09 65.16 81.86 Shoulder (q 5 ) -73.51 -60.78 -4.83 -58.75 -78.70 Elbow (q 6 ) -40.95 -36.38 -78.52 -44.87 -30.65 Hip (q 7 ) -62.75 -30.75 -7.79 -38.42 -55.47 Knee (q 8 ) 110.08 42.94 5.56 55.43 93.86 Ankle (q 9 ) -42.01 -10.97 -5.83 -14.75 -32.25

and (37) ensure the continuity between joint angles and angular velocities from one cycle to the next

where q i k 0 ð Þ, q i k T ð Þ denote the joint angles at the initial and final time instants of the kth cycle, respectively, with k ¼ 2, 3, 4, …, n iter : A small tolerance of 60:001 is considered for the constraints in Eq. (36). _ q i k 0 ð Þ and _ q i k T ð Þ denote the joint velocities at the initial and final time instants of the kth cycle. A small tolerance of 60:02 is considered for the constraints in Eq. (37).

Table 6 Fatigue-based constraints for 4CCr fatigue model Fatigue-based constraints for the 4CCr fatigue model Governing equations of the fatigue model

Four governing equations of the 4CCr fatigue model set up as equality constraints Unit summation

Any motor units assigned to the compartments for a joint must add up to 1

Residual capacity

Motor units available from the resting and active compartment sizes must meet or exceed the target load for a sustained task Fig. 5 Comparison of joint angle profiles predicted by the 3CC fatigue model to experimental data through 11 lifting cycles Fig. 6 Predicted optimal trajectory for each cycle predicted by the 3CC fatigue model: (a) t 5 0 s, (b) t 5 T/4 s, (c) t 5 T/2 s (d) t 5 3T/4 s, and (e) t 5 T s 4. Note that the optimization formulation with the 4CCr fatigue model has additional design variables, fatigue-based constraints, and continuity conditions.

3.2.1 Design Variables.

The additional design variables include the peripheral and central fatigue compartment profiles.

Alternative Fatigue Constraints.

The optimization formulation with the 4CCr fatigue model uses the same physics-based constraints as the formulation with the 3CC fatigue model as in Table 4 and has additional alternative fatigue-based constraints as in Table 6.

Alternative Continuity

Constraints. The boundary condition constraints for fatigue compartment sizes are given by the following equations:

The boundary condition constraints for fatigue compartment size velocities are given by the following equations:

The initial conditions for the fatigue compartment sizes in the first lifting cycle are given by the following equations: Fig. 8 Compartment size profiles of the physical joints of interest predicted by the 4CCr fatigue model through 13 lifting cycles Fig. 9 Optimal trajectory for each cycle predicted by the 4CCr fatigue model: : (a) t 5 0 s, (b) t 5 T/4 s, (c) t 5 T/2 s (d) t 5 3T/4 s, and (e) t 5 T s 061005-10 / Vol. 147, JUNE 2025 Transactions of the ASME

The continuity conditions for the joint angle and joint angular velocity profiles are the same as in the 3CC fatigue model optimization formulation.

Experiments

The sole participant of this study was a 33-year-old healthy adult male, with a height of 1.7 m and a weight of 67.3 kg. An Xsens Awinda inertial motion capture system (Movella Inc., Henderson, NV) running at 60 Hz was used to record the participant's motion. The box weight is 20 kg for this study.

For the data collection, the subject was instructed to stand with their feet spread apart and positioned symmetrically about the sagittal plane. The box was then placed centrally in front of them at a distance they self-determined to be comfortable for lifting from, and this position on the ground was marked prominently with painter's tape. They were then instructed to adopt a squatlifting strategy to lift the box up to pelvis height and lower it to the marked original position without touching the ground and then repeat the cycle till exhaustion. An experimenter provided continuous guidance to help the subject achieve consistent lift heights.

Collected data was processed in Xsens MVN Analyze 2023.0 (Movella Inc.) and Visual3D v6 (C-Motion Inc., Boyds, MD) to calculate the joint angle profiles.

Results

Using the optimization formulations described in Sec. 3, the optimal motion and fatigue progression during repetitive lifting tasks with a 20 kg box have been obtained. It is worth noting that the box weight falls within the maximum recommended weight limit for box lifting set by the U.S. Department of Labor, Occupational Safety, and Health Administration. The simulations are performed on a laptop equipped with an 11th Gen Intel Core i9-11900H processor with a base speed of 2.50 GHz featuring eight cores and 16 logical processors, and 32 GB RAM. Throughout the optimization process, various constraints are active, including joint limits, foot contacting position, initial and final box locations, residual capacity, collision avoidance, compartment size continuity, and unit summation constraints. Additionally, in the figures depicting joint angle profiles, joint torque profiles, and joint fatigue profiles, the first three vertical gray dotted lines denote the first three cycles as shown in Fig. 4. These lines are positioned at regular intervals of T seconds (average cycle time) to distinguish each cycle.

Three-Compartment Controller Fatigue Model

Results. For the 20 kg box case, the simulation predicts a maximum of 11 cycles with an average cycle time of 2.8 s before task failure. In contrast, the subject lifted 13 times in the experiment. The primary reason for task failure after 11 cycles is that the residual capacity drops below the target load for the hip joint. Figure 4 illustrates the fatigue progression in the five physical joints of interest, indicating a general trend of increasing fatigued compartment size over time. The active compartment size profiles closely follow the shape of the target load profiles, which are determined by optimal joint torques. Conversely, the resting compartment size profiles exhibit the opposite trend, as dictated by the governing equations of the 3CC fatigue model, ensuring the active compartment is adequately supported throughout the repetitive task. For joint angle validation in Fig. 5, compared to experimental data, similar trends are observed in the joint angle profiles for the six physical joints during the lifting and setting down phases of each cycle. The simulations generally match the experimental joint angle values, with the two differing by less than 12 deg. These similarities are also confirmed by the motion snapshots (Fig. 6). The joint torque profiles (Fig. 7) also demonstrate the periodic nature of the lifting motion.

Four-Compartment Controller With Augmented

Recovery Fatigue Model Results. For repetitive lifting motion prediction using 4CCr fatigue model, it successfully predicts 13 lifting cycles, which is the same as the experiment. The fatigued compartment sizes predicted are shown in Fig. 8. The optimal snapshots for one lifting cycle are depicted in Fig. 9. The predicted joint angle profile are shown in Fig. 10. Finally, the predicted joint torques are shown in Fig. 11.

Discussion

The RMSE comparison between the 3CC and 4CCr fatigue models across different joints provides valuable insights into their respective performances. Examining Table 7, we observe that the 4CCr model exhibits similar RMSE values compared to the 3CC model, indicating both fatigue models have similar level of accuracy for joint angle profiles predictions.

For the spine, the 4CCr model shows an RMSE of 10.4935, slightly higher than the 3CC model's 8.9914. The discrepancy is more pronounced in the shoulder, where the RMSE for the 4CCr model is 14.6236, compared to 9.8194 for the 3CC model. Conversely, for the elbow, the 4CCr model demonstrates better performance with an RMSE of 12.2949, as opposed to the 3CC model's higher value of 16.4669. The hip results again favor the 3CC model, with RMSE values of 9.1381 compared to 11.6677 for the 4CCr model. This trend continues with the knee, where the 4CCr model has an RMSE of 16.1694, significantly higher than the 3CC model's 13.0435. Finally, for the ankle, the 4CCr model has a slightly lower RMSE of 8.4887 compared to the 3CC model's 9.1383. This indicates that joint angle predictions alone may not provide a sufficient basis to favor one model over another for a similar simulation. This is likely because the joint angle predictions are mainly controlled by physics-based constraints which are the same in both fatigue models' optimization formulations. 061005-12 / Vol. 147, JUNE 2025 Transactions of the ASME

In this study, we used multi-objective cost functions which include a dynamic effort function and a joint angular acceleration function. The dynamic effort function is related to the subject's mechanical energy during repetitive lifting which is to be minimized. The joint angular acceleration is minimized to reduce the oscillations in the predicted lifting motion. Note that in Figs. 7 and 11, there are sudden changes in joint torques between lifting cycles. This is because no continuity conditions are imposed on the joint accelerations between the cycles in the optimization formulation. Joint torques are computed from EOM which requires joint acceleration. Due to the discontinuity in the joint acceleration, there are sudden changes in the predicted joint torques.

The critical finding is that the 3CC fatigue model predicts two fewer lifting cycles (11 cycles) than the experiment (13 cycles). In contrast, the 4CCr fatigue model predicts exactly 13 cycles as observed in the experiment. The discrepancy in predicted cycles may primarily be attributed to the inherent nature of the fatigue models. The 4CCr fatigue model is able to predict repetitive lifting cycles more accurately than the 3CC fatigue model as the former has been designed and validated for dynamic tasks, while the latter has only been validated for static tasks. It must also be considered that validation with a larger number of subjects will be required for future investigation. There remain, however, important limitations to both fatigue models. Fatigue coupling between adjacent joints was not considered in the current optimization formulation and is not inherently provided for by either model. Model constants representing fatigue and recovery coefficients were also obtained from an aggregate population, even though they were used for a single individual in this simulation. Furthermore, the subject was unable to achieve perfectly continuous lifts with zero delay between the end of one lift and the start of another. This short delay may have allowed for partial recovery, resulting in less fatigue than predicted.

For future improvements, the formulation's robustness could be enhanced by exploring other global optimization algorithms, such as the genetic algorithm to solve the optimization process. This study provides valuable insights into fatigue progression and lifting strategies during repetitive tasks, with the potential for practical applications in understanding and mitigating work-related injuries. The generalized nature of the formulation allows its application in various studies concerning repetitive motions and lays the foundation for further research in biomechanics and ergonomics to improve workplace safety and human performance [31][32][33].

Conclusions

This study presents a comprehensive examination of fatigue progression and optimal motion trajectory during a repetitive lifting task by comparing two fatigue models. The 4CCr fatigue model considers the effect of joint angular velocity on the fatigue parameters for both centrally and peripherally fatigued compartments, in contrast with the 3CC model which is unaffected by joint angular velocity, and does not differentiate between fatigue mechanisms. The results obtained from the simulations are subsequently validated through experiments conducted under similar conditions. The simulation results demonstrate periodic patterns in joint angle and torque profiles, showcasing the distinctive phases of the repetitive lifting task. For a subject with the same anthropometry as the participant of this study lifting a 20 kg box, the simulations predict a maximum of 11 cycles for 3CC fatigue model, compared to the 13 cycles for 4CCr fatigue model and experimental result. The results indicate that the 4CCr muscle fatigue model provides enhanced accuracy over the 3CC model for predicting task duration (number of cycles) of repetitive lifting. The fatigue progression across all five physical joints displays a general trend of increasing compartment size over time, with active and resting compartment profiles following expected patterns. These predictive methodologies have significant implications for ergonomic design and workplace injury prevention.

Table 7 RMSE comparison of 3CC and 4CCr fatigue models DOF 4CCr (deg) 3CC (deg) Spine (q 4 ) 10.4935 8.9914 Shoulder (q 5 ) 14.6236 9.8194 Elbow (q 6 ) 12.2949 16.4669 Hip (q 7 ) 11.6677 9.1381 Knee (q 8 ) 16.1694 13.0435 Ankle (q 9 ) 8.4887 9.1383