## 2 2D Skeletal Human Model  
Figure 1 shows a 2D skeletal model which has 14 degrees of freedom (DOFs). This skeletal model is adopted to simulate symmetric delivery motion. Three DOFs are used for global rotation (\(\beta\)) and translation (\(y\), \(z\)), and 11 DOFs (\(q_{1}\), \(q_{2}\), \(\ldots q_{11}\)) are used for modeling the human body joints. The global DOFs consist of one revolute joint and two prismatic joints. The left and right legs and arms  
are assumed to be symmetric. There are three body branches: two leg branches and one spine-arm branch. Each leg has a thigh, a shin, a hinder foot, and a forefoot. In the arm branch, two arms are combined as a single arm branch, because only 2D motion in sagittal plane is studied. The 50 th percentile anthropometric data from the experimental subjects are used for the skeletal model. Figure 1 shows how the DOFs are built in the Denavit-Hartenberg (DH) approach [28]. The DOF of each joint is defined in the local \(z\)-direction. In global Y-Z plane, the positive directions for the body rotation joints (\(q_{1}-q_{11}\)) are clockwise. The global rotation joint (\(\beta\)), spine joint (\(q_{1}\)), and two hip joints (\(q_{4}\) and \(q_{8}\)) coincide at the same location. A recursive Lagrangian dynamics formulation is utilized for this model, and the ground reaction force (GRF) is calculated from the gravity, joint kinematics, and external loads by a two-step inverse procedure [21]. The details are presented by [PERSON] et al. [29]. The forward joint kinematics are calculated as:  
\[\mathbf{A}_{i}=\mathbf{T}_{1}\mathbf{T}_{2}\mathbf{T}_{3}\cdots\mathbf{T}_{i }=\mathbf{A}_{i-1}\mathbf{T}_{i} \tag{1}\]  
\[\mathbf{B}_{i}=\dot{\mathbf{A}}_{i}=\mathbf{B}_{i-1}\mathbf{T}_{i}+\mathbf{A} _{i-1}\frac{\partial\mathbf{T}_{i}}{\partial q_{i}}\dot{q}_{i} \tag{2}\]  
\[\mathbf{C}_{i}=\dot{\mathbf{B}}_{i}=\ddot{\mathbf{A}}_{i}=\mathbf{C}_{i-1} \mathbf{T}_{i}+2\mathbf{B}_{i-1}\frac{\partial\mathbf{T}_{i}}{\partial q_{i}} \dot{q}_{i}+\mathbf{A}_{i-1}\frac{\partial^{2}\mathbf{T}_{i}}{\partial q_{i} ^{2}}\dot{q}_{i}^{2}+\mathbf{A}_{i-1}\frac{\partial\mathbf{T}_{i}}{\partial q _{i}}\dot{q}_{i} \tag{3}\]  
where \(q_{i}\) is the \(i\)th joint angle; \(\mathbf{T}_{i}\) is the link transformation matrix from the \((i-1)\)th link frame to the \(i\)th link frame; \(\mathbf{A}_{i}\), \(\mathbf{B}_{i}\), and \(\mathbf{C}_{i}\) are global position, velocity, and acceleration transformation matrices, respectively; and \(\mathbf{A}_{0}=[\mathbf{I}]\), \(\mathbf{B}_{0}=\mathbf{C}_{0}=[0]\). The backward joint torque equations are defined as follows:  
\[\tau_{i}=tr\bigg{(}\frac{\partial\mathbf{A}_{i}}{\partial q_{i}}\mathbf{D}_{ i}\bigg{)}-\mathbf{g}^{T}\frac{\partial\mathbf{A}_{i}}{\partial q_{i}} \mathbf{E}_{i}-\mathbf{f}_{k}^{T}\frac{\partial\mathbf{A}_{i}}{\partial q_{ i}}\mathbf{F}_{i}-\mathbf{G}_{i}^{T}\mathbf{A}_{i-1}\mathbf{z}_{0} \tag{4}\]  
\[\mathbf{D}_{i}=\mathbf{I}_{i}\mathbf{C}_{i}^{T}+\mathbf{T}_{i+1}\mathbf{D}_{ i+1} \tag{5}\]  
\[\mathbf{E}_{i}=m_{i}\mathbf{r}_{i}+\mathbf{T}_{i+1}\mathbf{E}_{i+1} \tag{6}\]  
\[\mathbf{F}_{i}=\mathbf{r}_{i}\mathbf{\delta}_{ik}+\mathbf{T}_{i+1}\mathbf{F}_{ i+1} \tag{7}\]  
Figure 1: The 2D skeletal model. \[\mathbf{G}_{i}=\mathbf{h}_{k}\delta_{ik}+\mathbf{G}_{i+1} \tag{8}\]  
where \(tr(\cdot)\) is the trace of a matrix, \(\mathbf{I}_{i}\) is the inertia matrix for link \(i\), \(\mathbf{D}_{i}\) is the matrix for inertia and Coriolis, \(\mathbf{F}_{i}\) is the external force torque vector, \(\mathbf{E}_{i}\) is the gravity torque vector, \(\mathbf{G}_{i}\) is the external moment torque vector, \(\mathbf{g}\) is the gravity vector, \(\mathbf{r}_{i}\) is the center of mass of link \(i\), \(m_{i}\) is the mass of link \(i\), \(\mathbf{h}_{k}=\left[ \begin{array}{cc}h_{x}&0&0&0\\ \end{array} \right]^{\text{T}}\) is the external moment applied on link \(k\), \(\mathbf{f}_{k}=\left[ \begin{array}{cc}0&f_{ky}&f_{kz}&0\\ \end{array}