This paper proposes an approach for the optimal design of nonlinear energy sinks (NES) installed in structures subjected to stochastic excitation. The response of these nonlinear systems subjected to stationary or non stationary stochastic excitations has been characterized by the time-history of its covariance matrix. However, time-invariant equivalent linearization methods of NES equipped structures will not always converge due to its essential nonlinearity. A solution strategy is proposed using equivalent linearization at each time step of the governing matrix differential equation, obtained through a sequence of Lyapunov equations. The adequacy of this method for optimization problems is evaluated by comparison with Monte Carlo simulations of two representative NES examples. Moreover, the influence of the Gaussian assumption of the equivalent system response is studied in these examples. An optimization formulation to design systems equipped with NES devices is proposed, in which the performance function is defined in terms of the maximum variance over the time history. This result is obtained by using the proposed method with a Gaussian assumption of the equivalent linear system. To illustrate the proposed methodology, parametric optimization of an NES installed at the top of a shear building structure is performed.