Manual material handling (MMH), particularly lifting, is one of the main reasons for work-related joint and back injuries, 1 which is the most common reason for seeking medical care for civilians 2 and the military. 3 Injuries associated with MMH play a significant role in the economy. A study on workers' compensation claims shows that about 32% of all claims and 36% of compensation costs are related to MMH. 4 In the United States, the economic impact of MMH related injuries such as low back pain and musculoskeletal disorder is more than $100 billion per year, considering the direct and indirect costs. 5 The direct costs of MMH injuries are over $13 billion in 2016 in the USA. 6 According to the U.S. Bureau of Labor Statistics, health care and social assistance, manufacturing, and retail trade are some of the most affected sectors in private industry, based on number of nonfatal occupational injuries and illnesses. 7 The lack of knowledge and training about proper lifting strategy, and awareness about the longterm consequences among the industrial workers make them the most vulnerable persons to work related musculoskeletal disorders. Therefore, the biomechanics of lifting is a critical issue in many industrial applications. It is necessary to determine subject-specific maximum lifting weight and explain why these lifting related injuries occur. However, it is challenging to determine the maximum lifting weight by experiments, as it is hazardous and risky for the participants.
In industrial settings, asymmetric lifting is more common than symmetric lifting. For symmetric lifting the center of gravity moves along the sagittal plane, whereas for asymmetric lifting the center of gravity moves along both the sagittal and frontal plane, and the lateral bending moment on spine reduces the material handling capability by 16%. 8 That is why there is a significant difference of the maximum weight predictions between symmetric and asymmetric lifting tasks.
Over the past few years, many researchers conducted lifting motion predictions. However, most of the studies focus on symmetric lifting. [9][10][11][12][13][14] Researches on asymmetric lifting prediction are very few, and most of them are based on static lifting. 15 Maximum lifting weight prediction requires the dynamic strength in joint space. [16][17][18][19] In muscle space, the crossing muscles' net moment generating capacity is based on muscles' activations, strength, and moment arms. Furthermore, muscles' strength surface (force-length-velocity) and moment arm properties are changing with the joint angles and angular velocities. Thus, the dynamic joint torque limit is a function of joint angle and angular velocity. Gu¨ndogdu et al. 9 studied 2D lifting prediction considering dynamic joint torque limits. The optimization predicted the maximum box weight, optimal lifting motion, and total time.
The goal of this work is to build and validate a predictive 3D asymmetric maximum weight lifting model. The maximum lifting weight, optimal lifting motion, and lifting time duration are all predicted based on the given box locations and the subject's anthropometric data and strength. The predicted results are validated with the lifting experimental data (mean and standard deviations).
A total of 23 healthy male participants (20-50 years old) were recruited for the experiments. We discarded those incomplete data during post-processing and 15 participants' data were valid. Because a back-lifting approach was used for three participants these three subjects were also discarded. The remaining 12 participants were used for this paper (age: 25.42 6 7.72 years; height: 182.2 6 3.6 cm; body mass: 84.16 6 10.16 kg) and all 12 participants used the squat-lifting strategy for heavy load. The anthropometries of these participants are shown in Table 1. The following criteria were used for participant selection: they should be mentally and physically sound; they should be able to perform the scripted task, and they should not be on any medication that might hamper their performance during the box-lifting task. In addition, these participants consisted of students and university staff that did not have explicit training in lifting or physical fitness but were otherwise healthy and capable of performing the boxlifting task. The Texas Tech University Institutional Review Board approved the lab experimental protocol and all participants signed their informed consent form.
One motion capture system with with cameras (Vicon Motion Systems, Oxford, UK) was used to collect 3D lifting motion data (100 Hz). Two force plates (Kistler, Winterthur, Switzerland) were used to collect ground reaction forces (GRFs) with 2000 Hz. The Vicon plug-in Gait model was used for marker protocol (markers with 9 mm, spherical) and an additional two markers were placed on the iliac crests so that a total of 42 markers were used. 20 The following anthropometric parameters were measured for each participant: leg length, ankle width, knee width, shoulder offset, elbow width, wrist width, hand thickness, waist circumference, inter-ASIS distance, height, and weight. 21,22 Each subject was instructed to psychophysically determine their maximum lifting load by gradually increasing the load on the box (65 cm 3 35 cm 3 15 cm) until the subject felt the load was too heavy to safely carry. The subjects were instructed to lift the box in the most comfortable and natural way and then to put the box on the table. Note that the real maximum lifting capacity was not used so that subjects could avoid injury during the experiment, i.e., each subject had a safety factor in determining his maximum weight. Therefore, the maximum weight value in this paper refers to the maximum weight a subject can safely carry. Once the maximum weight was determined, the lifting study was started. Because the box did not have any handle, it was initially placed in front of the participant on top of a weight disk that was 2.54 cm above the floor to allow for a better grasp on the box. While placing one foot on each force plate, the subject performed three trials of the lifting task that the participant lifted the box from the weight disk and set it down on a 1-m-tall table to their right, that is, an asymmetric lifting, shown in Figure 1. Between any two adjacent trials the participant took a 5 min break.
After data collection were done, data post processing was conducted in Vicon Ò software. Marker data were labelled, smoothed, and finally converted into a C3D file, which could be input into Visual 3D Ò software (C-Motion, Inc., Germantown, MD, USA). In Visual 3D, each participant's raw kinematic data and measured anthropometries were used to create hybrid models for 15 segments: CODA pelvis, trunk, thigh (bilateral), shank (bilateral) foot (bilateral), head, hand (bilateral), forearm (bilateral), and upper arm (bilateral). The kinematic and kinetic data were filtered using a Butterworth filter with cutoff frequencies of 6 and 25 Hz, respectively. Kinematic and kinetic data processing was conducted to extract the following variables: bilateral ankle flexion, bilateral knee flexion, bilateral hip flexion, spine flexion and rotation, bilateral elbow flexion, and bilateral vertical GRFs. ). The arms and legs are considered symmetric with respect to the sagittal plane of the spatial model. Each arm has three parts: upper arm, forearm, and hand. There are three DOFs for shoulder, two DOFs for elbow, and two DOFs for wrist. Each leg contains a rear foot, a forefoot, a shank, and a thigh. Each leg has seven DOFs: three for hip, one for knee, two for ankle, and one for the metatarsal joint at forefoot. 13 The anthropometric data are generated from Visual 3D software with experimentally measured height, weight, and stature data. The strength percentile is retrieved from symmetric maximum weight lifting. 24 The general equations of motion (EOM) of the spatial model can be expressed using the Recursive Lagrangian formulation in a matrix form which contains forward recursive kinematics and backward recursive dynamics. 25,26 Forward recursive kinematics:
where q i , _ q i , € q i are the joint angle, velocity, and acceleration, respectively, T i is the i th DH transformation matrix, 23 A i , B i , C i are the global recursive matrices for position, velocity, and acceleration, respectively, and
where tr Á ð Þ is the trace of a matrix, D i is the inertia and Coriolis matrix, E i is the gravity torque vector, F i is the external force torque vector, G i is the external moment torque vector, I i is the inertia matrix for link i, r i is the center of mass of link i, g is the gravity vector, m i is the mass of link i, f k = f kx 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 local frame k, h k = h x h y h z 0 ½ T is the external moment applied on link k,
½ are starting conditions for recursive matric and vectors.
The GRFs are calculated from a two-step active-passive algorithm 26,27 : in first step, the global forces and moments are calculated from given state variables without GRFs; in second step, the calculated global forces and moments are transferred to center of pressure and further applied to metatarsal joints as external forces and moments.
The maximum weight lifting task is formulated as a nonlinear programming (NPL) problem. In this formulation, the box dimensions, initial and final positions of the box, and initial, intermediate, and final key joint values and GRFs are obtained from the experiment.
Design variables. In the current optimization formulation, the design variables are the control points (c) of cubic B-spline interpolation of joint angle profiles for lifting motion, box weight W, and total time T as x = c T W T Â Ã T . The joint torques t t ð Þ are directly calculated from EOM using inverse dynamics, instead of integrating the differential equations.
Objective functions. The cost function J has large effect on the predicted motion. A multi-objective optimization (MOO) is used for the maximum weight lifting prediction by maximizing the box weight and minimizing the sum of joint torque squares. 28 The maximizing box weight cost function is transformed into a minimizing negative logarithmic function of box weight. There are two reasons for this transformation: one is that the optimizer can only handle the minimization type of problem; the other reason is that this transformation facilitates numerical convergence for optimization. Finally, the MOO cost function is defined as 28 :
where t L i is the lower joint torque limit and t U i is the upper joint torque limit, n is the number of DOF, w 1 and w 2 are coefficients for the two normalized cost functions where w 1 = 0:15 and w 2 = 0:85, 28 N Á ½ is the normalization function by dividing the function's maximum absolute value: for both negative logarithmic function of box weight and joint torque square function, their maximum absolute values are achieved by purely maximizing box weight at w 1 = 0 and w 2 = 1. 13,29,30 Constraints. The constraints imposed on the lifting motion can be divided into two types: time-dependent constraints and time-independent constraints. Timedependent constraints include dynamic joint torque limits, joint angle limits, dynamic stability, foot-contacting positions, box collision avoidance, hand distance, and box-ground parallel. These constraints are imposed throughout the lifting time interval T. Time independent constraints include initial and final static conditions, initial and final hand positions, initial, mid-time, and final key joint values from experiment, and initial, intermediate, and final GRF values from the experiment. These constraints are imposed only at specific time points of lifting motion.
(1) Joint angle limit represents the physical range of motion which is obtained from experiments:
where q L and q U represent the lower and upper limits on the joint angles respectively.
(2) Dynamic joint strength is imposed in this study. Dynamic joint torque limit is a function of strength percentile (z score ), joint angle (q), angular velocity (v), and time (t). The lower and upper dynamic joint strengths are:
Þrespectively. These two functions are regressed using logistic equations based on the dynamometer isometric and isokinetic strength data. [16][17][18][19]31
where e is the exponential function, c 1 ;c 8 are experimental regression coefficients, CV i U is the experimental upper coefficient covariance, and t i peak U is the upper peak torque value for the i th joint.
For the dynamic lower joint torque limit:
where e is the exponential function, d 1 ;d 8 are experimental regression coefficients, CV i L is the experimental lower coefficient covariance, and t i peak L is the lower peak torque value for the i th joint. The experimental dynamic strengths data for hip, knee, ankle, lower spine, shoulder, elbow, and wrist joints are obtained from the literature. [16][17][18][19]31 In equations ( 11) to ( 16), c 1 ;c 8 , d 1 ;d 8 , CV i U , and CV i L are regression values obtained from experiments. The subject's strength percentile z score is retrieved from an enumeration optimization process for symmetric maximum weight lifting based on experimental data. 24 The strength percentile is enumerated to increase the strength limits until the lifting optimization converges with the given box weight and find the optimal motion.
(3) The feet-contacting points are specified on the ground as follows:
where p feet are the calculated feet positions, p s are specified feet ground contacting positions.
(4) Dynamic stability is satisfied by constraining the zero-moment point (ZMP) inside the foot support region (FSR). 26,27,32 p ZMP (x, t) 2 FSR ð18Þ
(5) Box collision avoidance is considered in the optimization formulation to avoid penetration of the box into human body. The body thickness is represented by filling up the model with spheres on the hip, knee, ankle, thigh, shank, lower spine and higher spine. The distance (d) between the sphere center and box center is calculated at each time point for the box collision avoidance constraint.
where r is the radius of a sphere, and dep is the box depth. (6) The distance between the two hands in 3D space is a constant and equals to the width of the box. This constraint is expressed as,
where p right hand and p left hand are the right and left hand locations, respectively, and wid is the width of the box. (7) To keep the box parallel to the ground during the lifting process, the height of both hands should be same.
where h right hand and h left hand are the right and left hand heights, respectively.
Time-independent constraints.
(1) The initial and final hand locations are specified for lifting,
where p s box is the given experimental box location, and p hand is the calculated hand location.
(2) The human model is at rest at beginning and end of lifting motion.
(3) The initial, mid-time, and final key joint angles are specified. 14
where q E i is the given experimental joint angle including left and right ankle flexion, left and right knee flexion, left and right hip flexion, left and right elbow flexion, spine flexion and rotation, e = 10 degree at boundaries and e = 5 degree at mid-time point.
(4) Initial, intermediate, and final vertical GRFs are given from experimental data.
where GRF E left , GRF E right are the experimental vertical ground reaction force for right foot and left foot, respectively.
The experimental and predicted kinematics and kinetics were investigated in this study. The asymmetric lifting problem was solved by the optimizer SNOPT 33 using a sequential quadratic programming (SQP) method. It took about 7 CPU minutes to solve the problem on a laptop computer with Intel(R) Core(TM) i5-7200U 2.50 GHz processor and 8 GB RAM. The maximum lifted weight during experiment was 235.83 N (24.04 kg) for Subject #8. The predicted maximum lifted weight for Subject #8 on right hand and left hand were 128.39 and 128.39 N, respectively, that is, the predicted total lifted weight was 256.78 N that is 8.9% larger than the maximum weight of experiment. The optimal lifting time is 1.33 s. The strength (z score ) for the simulated model was retrieved from symmetric maximum weight lifting as 1.05. 14,24 Figure 3 depicts the predicted joint angle and vertical GRF profiles. The snapshots of the optimal asymmetric lifting motion are shown in Figure 4.
The predicted joint angles agree well with the experimental data. Although the predicted right ankle at the beginning of lifting (Figure 3(b)) and right knee at the last portion of lifting motion (Figure 3(d)) are outside of one standard deviation, the pattern and timing of phase change are consistent with the experimental data. The predicted hip flexions, elbow flexions, spine flexion, and spine rotation are all within one standard deviation of experimental joint angle profiles. However, for the simulation, the 3D human model started to straighten and rotate its spine earlier than the experimental subject did (Figure 3(g) and (h)). Figure 3(k) and (l) shows the comparisons of vertical GRFs on both feet during lifting. The predicted GRFs are within one standard deviation of experimental data, except the right GRF during 20%-35% of the task.
The predicted maximum lifting weight of Subject #8 is 8.9% higher than the experimental maximum lifting weight. As mentioned in the Experiments section, the maximum lifting weight determined during the experiment was a safe maximum weight. The true maximum lifting weight should be higher than the experimentally determined maximum lifting weight. The proposed MOO asymmetric maximum weight lifting prediction reveals this insight.
There are some minor discrepancies between the simulated and experimental joint angle profiles. One noticeable difference between the prediction and experimental data was the time lag of phase change for a small portion of the vertical GRF profiles (Figure 3(k) and (l)). The reason for this discrepancy might be early phase changes of spinal flexion and rotation compared to experimental data (Figure 3(g) and (h)). Early extension of the spine worked as a catalyst to give the model early upright standing stability and to start the rotation. That early upright standing is also the reason for the flat profile after 90% of the task for both GRFs. On the other hand, the experimental lifting strategy extended the spine later than the simulation model and rotated the body faster to place the box at the desired position on the right side. As a result, after the second peak (85% of the task) the subjects created higher GRF on the right side and lower GRF on the left side than the simulation did. Although there are some deviations of phase change for GRFs, the simulated profiles were almost within one standard deviation of the experimental data.
During the experiment, because of the heavy weight, the subjects tended to stand up straight first, and then rotate to the right. That is the reason of higher GRF on right side and lower GRF on left side after 90% of the task. The last-moment rotation of spine requires fast twisting of spine muscles with heavy weights on hand. Such kind of repetitive works can lead to chronic strain for spine muscles. Also, the higher GRF on right side and lower GRF on left side may also lead to high jointloads on right-side joints especially for lower extremity joints. Such kind of repetitive works may also lead to wear down of right hip and knee joints' cartilage, which is the cause of hip and knee osteoarthritis.
Although some experimental data are used in the optimization formulation (equations ( 24) and ( 25)) to guide the prediction, they are necessary to predict accurate results because of the complexity of the 3D asymmetric lifting motion. The previous studies showed that mid-time postures or key joint values could improve the accuracy of lifting prediction. 28,34 In this study, we tried to use minimal experimental data in the optimization formulation. It was found that three experimental intermediate GRF constraints (equation ( 25)) were necessary to capture the history of GRF profiles due to the fluctuating nature of asymmetric GRFs and the effects of spine flexion and rotation. Compared to regression models, 15 the proposed predictive model uses much less experimental data and has more powerful predictive capability. However, the model's predictive ability is compromised by the amount of experimental data used in the optimization formulation.
It is noted that we used symmetric maximum weight lifting strength for asymmetric maximum weight lifting prediction for subject #8. Here we assume that symmetric and asymmetric lifting strength percentiles are similar for this subject. This assumption is reasonably proved through these simulation and experiments in this study, as the simulation results predict accurate asymmetric lifting motion, box weight, and time duration compared to experimental data by using symmetric lifting strength. The determined strength percentile of a person is a global score considering interactions of all joints for a task. 24 The determined subject-specific strength value is critical to predict other strength related tasks to protect the subject from any injury risk in manual material handling. In case the symmetric and asymmetric strength percentiles are very different for a subject, the optimization-based enumeration retrieval approach 24 can be used to approximate the subject's strength percentile for the asymmetric lifting based on experimental data. This is the first study using MOO to predict 3D asymmetric maximum weight lifting motion considering the dynamic torque limits in the literature. Based on the comparisons with experimental data for both kinematics and kinetics, it is clear that, except for some minor discrepancies, the results of the predictive model demonstrated the ability to predict realistic 3D asymmetric lifting motion, accurate maximum lifting weight, and lifting time duration. This model also provides some insight view of 3D asymmetric maximum weight lifting by considering dynamic joint torque limits, which can be helpful when analyzing ergonomic safety problems involving lifting.
In this study, a 3D 40-DOF skeletal model was utilized to predict subject-specific asymmetric maximum lifting motion and box weight. The lifting task was formulated as a MOO problem by aggregating two cost functions: maximizing the box weight and minimizing the sum of joint torque squares. The lifting motion prediction problem was solved by an SQP based optimizer SNOPT considering dynamic joint strength as one of the constraints. The development of a predictive human model that can predict human kinematics and kinetics accurately is a big challenge. It is necessary to include a dynamic strength constraint to predict maximum lifting weight, optimal lifting motion, and lifting time duration. In addition, MOO can generate more accurate simulation than single objective optimization. 28 The validated dynamic-joint-strength-based 3D asymmetric lifting model will give researchers a robust tool to work on subject-specific motion analysis, which is helpful for designing workplace and ergonomic tools to avoid injury for lifting. In future work, we will (1) integrate the 3D skeletal motion prediction with a OpenSim model to study muscle activities 35 ; (2) integrate the model with the joint-space fatigue model 36 to study repetitive lifting; and (3) develop advanced ergonomic tools for the prevention of lifting injuries.
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Proc IMechE Part H: J Engineering in Medicine 235(4)