Abstract: When simulating the dynamics of muti-rigid body systems, the Index-1 differential-algebraic equations (DAEs) only consider the acceleration constraint equations and ignore the position and velocity constraint equations. During the numerical integration procedure, the constraints violation at position and velocity levels increase, leading to unreliability of numerical results. Accordingly, constraint stabilization methods have to be introduced to reduce the violation of the constraints. In this work, the dynamic equations of muti-rigid body systems on SE(3) are offered based on the Hamilton’s principle. Then four constraint stabilization methods on SE(3) are introduced: the Baumgarte stabilization method, the penalty method, the augmented Lagrangian formulation and the constraint violation stabilization upgraded method. Compared with the traditional dynamic equations of muti-rigid body systems, the dynamic equations on SE(3) solved by integration scheme consider the coupling of rotations and transitions, leading more accurate configuration results and calculation of constraint error. The accurate calculation of constraint error improves the performance of the constraint stabilization methods on SE(3). The four kinds of constraint stabilization methods on SE(3) with dynamic equations of muti-rigid body systems are respectively simulated by RKMK(Runge-Kutta Munthe-Kass) integration scheme. Finally, two numerical examples, including a spatial double pendulum with spherical hinges and a crank slider mechanism with cylindrical hinges, are presented. The numerical results related to constraints violation at position and velocity levels, conservation of the total energy and simulation time are compared and analyzed. It is concluded that the four constraint stabilization methods on SE(3) are effective for preserving structure and energy conservation where RKMK is employed to simulate the dynamic equations of muti-rigid body systems. Compared with the other three constraint stabilization methods on SE(3), the augmented Lagrangian formulation method on SE(3) can provide better numerical accuracy and smaller constraints violation at both position and velocity levels. The Baumgarte stabilization method on SE(3) can balance the calculation efficiency and numerical accuracy well.