Multibody dynamics simulations are currently widely accepted as valuable means for dynamic performance analysis of mechanical systems. The evolution of theoretical and computational aspects of the multibody dynamics discipline makes it conducive these days for other types of applications, in addition to pure simulations. One very important such application is design optimization for multibody systems. In this paper, we focus on gradient-based optimization in order to find local minima. Gradients are calculated efficiently via adjoint sensitivity analysis techniques. Current approaches have limitations in terms of efficiently performing sensitivity analysis for complex systems with respect to multiple design parameters. To improve the state of the art, the adjoint sensitivity approach of multibody systems in the context of the penalty formulation is developed in this study. The new theory developed is then demonstrated on one academic case study, a five-bar mechanism, and on one real-life system, a 14 degree of freedom (DOF) vehicle model. The five-bar mechanism is used to validate the sensitivity approach derived in this paper. The full vehicle model is used to demonstrate the capability of the new approach developed to perform sensitivity analysis and optimization for large and complex multibody systems with respect to multiple design parameters with high efficiency.