The equations governing the motion of constrained multibody systems are not described solely by differential equations, but also by algebraic constraints arising from kinematic relations. The coupling of these two components leads to Differential-Algebraic Equations (DAEs), which are known to pose significant numerical challenges that do not arise for Ordinary Differential Equations (ODEs). While any mechanical system with holonomic constraints can, in principle, be reformulated as a set of Ordinary Differential Equations (ODEs) by adopting a minimal set of coordinates, identifying the analytical mapping between the full (redundant) set of coordinates and the minimal one is often challenging. This mapping and the resulting ODE in minimal coordinates is typically nonlinear and can become, in general, highly complex, which reduces its applicability for simulation purposes. In this work, we aim at learning solutions of constrained multibody systems, in order to support computationally efficient simulation. In doing so, we aim to ensure that the solutions of these surrogates not only provide accurate and efficient approximations but also remain consistent with the underlying physics, particularly, in satisfying the problem’s constraints exactly. Hereto, we adapt the index-aware learning framework, originally developed for electrical circuits, to the context of constrained multibody dynamics. By leveraging the concept of the dissection index, we are able to identify the truly differential variables associated with the system’s minimal coordinates. Subsequently, we aim to learn surrogates for the solutions of the unconstrained dynamics. This strategy enables us to learn the time- and parameter-dependent solutions of the system, while inherently respecting all kinematic constraints.