SUMMARYDynamics of mechanical systems during motion usually involves reaction forces and moments due to the interaction with external objects or constraints from the environment. The problems of determining the external reaction loads have often been addressed in the literature in terms of forward dynamics solution methods for non-redundant systems. In this article, we propose a numerical formulation of distributing the external reaction loads of redundant systems concurrently with optimal motion planning and simulation. For dynamic equilibrium, the resultant reaction loads that include the effects of inertia, gravity, and general applied loads are distributed to each contact point. Unknown reactions are determined along with the system configuration at each time step using iterative non-linear optimization algorithm. The required actuator torques as well as the motion trajectories are obtained while satisfying the given holonomic constraints. The formulation is applied with two types of discretization methods for the reactions, and several example motions of multi-rigid-body systems, such as a simple 3-degree-of-freedom mechanical manipulator and a highly articulated whole-body human mechanism, are demonstrated. The example results are compared with the cases where the reactions are pre-assigned. The proposed numerical formulation demonstrates the realistic distribution of external reaction loads and the associated goal-oriented motions that are dynamically consistent, and provides alternate schemes for hybrid motion/force control of redundant robots.