Nonlinear Energy Sinks (NESs) are a promising technique for passively reducing the amplitude of vibrations. Through nonlinear stiffiness properties, a NES is able to passively and irreversibly absorb energy. Unlike the traditional Tuned Mass Damper (TMD), NESs do not require a specific tuning and absorb energy from a wide range of frequencies. How- ever, they are only efficient over a limited range of excitations. In order to mitigate this limitation, this work studies the optimal design of several NESs configured in parallel. In addition, it is known that the efficiency of a NES is extremely sensitive to small perturbations in design parameters or initial conditions, with a near discontinuous behavior. In order to perform optimization while accounting for the high sensitivity to uncertainties, an optimization technique for the maximization of the mean efficiency is proposed. It is based on the identification, using a support vector machine (SVM) and clustering, of low and high efficiency regions as well as the corresponding two Kriging approximations. SVM and Kriging approximation are refined locally using the dedicated adaptive sampling scheme. This work will compare optimized designs for maximum mean efficiency with various numbers of parallel NESs. In all problems, the excitation of the main system is provided in the form of an initial velocity treated as an aleatory variable. The NES nonlinear stiffiness properties, considered random design variables, are optimized for cases with 1, 2, 3, 5, and 10 NESs in parallel.