The structural optimization of the components in multibody systems is performed using a fully coupled optimization method. The design's predicted response is obtained from a flexible multibody system simulation under various service conditions. In this way, the resulting optimization process enhances most existing studies which are limited to weakly coupled (quasi-) static or frequency domain loading conditions. A level set description of the component geometry is used to formulate a generalized shape optimization problem which is solved via efficient gradient-based optimization methods. Gradients of cost and constraint functions are obtained from a sensitivity analysis which is revisited in order to facilitate its implementation and retain its computational efficiency. The optimizations of a slidercrank mechanism and a 2-dof robot are provided to exemplify the procedure.