## 4 Subtask-Based Optimization Formulation
### Lifting and Unloading Subtasks  
For the lifting subtask, the total lifting time, initial and final postures, and lifting weight are specified in the simulation from the experimental mean values. The design variables are joint angle profile cubic B-spline control points (\(\mathbf{P}_{l}\)). The cost function \(J_{1}\) is the sum of the normalized joint torque squared as in Equation (9) for the lifting period \(T_{l}\).  
\[J_{1}(\mathbf{P}_{l})=\sum_{i=4}^{n}\int_{t=0}^{T_{l}}\left(\frac{\tau_{i}( \mathbf{P}_{l})}{\tau_{i}^{U}-\tau_{i}^{L}}\right)^{2}\mathrm{d}t \tag{9}\]  
where \(n\) is the number of DOF, \(\tau_{i}^{U}\) is the upper limit for the \(i\)th joint torque, and \(\tau_{i}^{L}\) is the lower limit for the \(i\)th joint torque.  
The lifting subtask is subject to the following constraints:  
Joint angle limits,  
\[\mathbf{q}^{L}\leq\mathbf{q}(t)\leq\mathbf{q}^{U} \tag{10}\]  
where \(\mathbf{q}^{L}\) is the lower joint angle limit and \(\mathbf{q}^{U}\) is the upper joint angle limit.  
Joint torque limits,  
\[\boldsymbol{\tau}^{L}\leq\boldsymbol{\tau}(t)\leq\boldsymbol{\tau}^{U} \tag{11}\]  
In addition, the feet contacting positions are given from experiments,  
\[p_{foot}(\mathbf{q},\ t)=p_{foot}^{E}(t) \tag{12}\]  
where \(p_{foot}\) is the calculated foot position, and \(p_{foot}^{E}\) is the measured Cartesian coordinates of footstep position during delivering motion based on the mean experimental step lengths.  
The foot is contacting the ground, the height of the foot-contacting points is zero. The other points' heights should be greater than zero.  
\[\begin{array}{l}y_{i}(t)=0;\ i\in\ contact\\ y_{j}(t)\geq 0;\ j\
otin\ contact\end{array} \tag{13}\]  
where \(i\) is the point index on foot which is contacting the ground [21].  
Balance is considered during the box delivering process using the zero-moment point (ZMP) constraint,  
\[p_{ZMP}(\mathbf{q},\ t)\in\text{FSR} \tag{14}\]  
where FSR is the foot support region, and \(p_{ZMP}\) is the ZMP location. This constraint is defined by constraining the sign of cross product of two vectors: one vector is composed of ZMP and FSR vertex, the other vector is the adjacent FSR boundary. This process is repeated until all FSR boundaries are satisfied. Details are given in reference [32].  
Collision avoidance between box and body,  
\[d(\mathbf{q},\ t)\geq r_{1}+r_{2} \tag{15}\]  
where \(r_{2}\) is the radius of the circle filled in the body segment, \(r_{1}\) is half of the box width, and \(d\) is the distance between each circle center and the hand. There are five circles filled in the lower body at hip, knee, ankle, mid of shank, and mid of tibia, respectively. In addition, there are two circles filled on spine to represent low spine and high spine segments.  
The differences of joint angle between the experiment and the skeletal model are constrained in a small range \(\varepsilon=10\) degrees at initial and final lifting/unloading postures, where \(\mathbf{q}^{E}\) is the experimental joint angle.  
\[|\mathbf{q}(t)-\mathbf{q}^{E}(t)|\leq\varepsilon;\ t=0,T_{l} \tag{16}\]It is noted that the unloading subtask starts with holding the box at a higher position and ending at a lower position during the time period \(T_{ul}\). The optimization formulation for unloading subtask is similar to the lifting subtask.