Box delivery is a complicated manual material handling task which needs to consider the box weight, delivering speed, stability, and location. This paper presents a subtask-based inverse dynamic optimization formulation for determining the two-dimensional (2D) symmetric optimal box delivery motion. For the subtask-based formulation, the delivery task is divided into five subtasks: lifting, the first transition step, carrying, the second transition step, and unloading. To render a complete delivering task, each subtask is formulated as a separate optimization problem with appropriate boundary conditions. For carrying and lifting subtasks, the cost function is the sum of joint torque squared. In contrast, for transition subtasks, the cost function is the combination of joint discomfort and joint torque squared. Joint angle profiles are validated through experimental results using Pearson’s correlation coefficient (r) and root-mean-square-error (RMSE). Results show that the subtask-based approach is computationally efficient for complex box delivery motion simulation. This research outcome provides a practical guidance to prevent injury risks in joint torque space for workers who deliver heavy objects in their daily jobs.