## 1 Introduction  
The box delivery process includes multi-task jobs (e.g., that involves box lifting, carrying, and lowering it). In multi-task case, workers combine subtasks with transition phases which affect the pace, posture, motion, and the resulting biomedical loads of the workers [1]. Biomechanical modeling provides a useful tool to investigate the cause and effect without worrying about injury risks. However, it proves that simulating a human delivering motion is a challenging task. The objective of this study is to explore the subtask-based inverse dynamic optimization formulation to efficiently analyze and predict dynamic human delivering motion in ergonomic applications. The box delivery task can be divided into four main subtasks: lifting, carrying, transition, and unloading. The transition task includes the motion from lifting to carrying and carrying to unloading with varied boundary conditions. Each subtask can be formulated as an optimization problem by associating boundary conditions between the two adjacent subtasks. The task decomposition approach was used to simulate walking motion with double support and single support phases [2], but only few studies have been conducted on the delivery task (multi-task jobs). Here is a summary of the literature for carrying, lifting, and transition tasks. There are two different methods for optimization-based lifting planning problems: static lifting optimization [3; 4] and dynamic lifting optimization [5; 6; 7; 8; 9]. Various cost functions were used in the optimization formulation, such as muscle intensity [3], joint torque squared [7], compressive force on spine [10], and smoothness of the lifting motion [6; 8]. [PERSON] and [PERSON] [3] investigated various cost functions for static lifting using a simple geometrical model. They found that muscle intensity was a more efficient cost function for lifting. Three objective functions for optimal lifting postures were studied by [PERSON] and [PERSON] [4]. They concluded that minimizing the overall effort generated the most accurate motion with subject validation. In addition, the hand's relative height location was an important factor for the accuracy of the prediction. [PERSON] and [PERSON] [10] calculated the stresses on major joints using a 2D biomechanical lifting model. The effects of lifting technique on the spine was investigated. [PERSON] and [PERSON] [6] studied three cost functions using dynamic lifting simulation. The three cost functions were: to minimize the sudden kinematic variance of the center of gravity, to maximize the lifting motion's smoothness, and to minimize the sum of joint torque squared. They found that the peak spine loads were minimal when the load took the smoothest path. [PERSON] et al. [11] simulated a tossing motion with three phases. It was found that mobilization and torque utilization were important cost functions in the final phase of tossing. [PERSON] et al. [7] aggregated dynamic effort and stability objective functions to predict lifting motion using a multi-objective optimization method. The proposed approach successfully found the optimal lifting motion which was validated with experimental data. The lifting height and weight correlate with the shear and compressive forces on the spine. Many task-analysis tools (e.g., the revised NIOSH lifting equation) are rooted in a task-based analysis approach [12]. These tools divide the multi-task jobs into single subtasks and analyze each individual subtask separately, and synthesize the results into the entire job. Several studies reported differences in the metabolic rates, perceived exertion, heart rates, hand forces, compression, and shear forces of workers when they performed these multi-task jobs compared to single-task jobs [13; 14; 15; 16]. It has been concluded that the risk of injury during multi-task jobs could not be accurately assessed using tools that analyze single tasks [13; 14; 17].