Summary: The geometrical structure of constraints leads to the important conclusion that components of velocities and accelerations in the normal space of the constraint hypersurface are totally determined by the constraints. Orthogonal base vectors of the constraint’s normal and tangent spaces can be obtained by the QR decomposition with column permutation. A numerical method of closed-loop multibody systems with singular configurations is presented, in which the traditional hypothesis on independence of constraints is abandoned. Instead of correcting constraint violations at the end of each integration step, positions and velocities are modified to satisfy their constraint equations before they are used to form equations of motion. Such an approach can collaborate with any standard ODE solver. In order to systematically generate constraint equations and obtain the corresponding joint’s reaction forces, rotational and translational constraints of a closed-loop are standardized as part of six explicit equations, and a clear relationship between Lagrange multipliers and joint’s reaction forces is derived based on the principle of virtual power equivalence. The proposed method can dynamically identify independent constraints and modify the equations of motion accordingly. The numerical examples validated its effectiveness.