This work presents a novel cluster based optimization procedure for estimating parameter values that yield stable, periodic responses with desired amplitude in nonlinear vibrating systems. The parameter values obtained by conventional nonlinear optimization schemes, with minimization of amplitude as the objective, may not furnish periodic and stable responses. Moreover, global optimization strategies may converge to isolated optima that are sensitive to parametric perturbations. In practical engineering systems, unstable or isolated optimal orbits are not practically realizable. To overcome these limitations, the proposed method tries to converge to a cluster in the r-dimensional parameter space in which the design specifications including the specified optimality, periodicity, stability and robustness are satisfied. Thus, it eliminates the need for computationally expensive bifurcation studies to locate stable, periodic parameter regimes before optimization. The present method is based on a hybrid scheme which involves the algebraic form of the governing equations in screening phase and its differential form in the selection phase. In the screening phase, force and energy balance conditions are used to rephrase the nonlinear governing equations in terms of the design parameter vector. These rephrased equations are reduced to algebraic form using a harmonic balance procedure which also specifies the desired amplitude and frequency of the response. An error norm based on this algebraic form is defined and is used to contract the search bounds in the parameter space leading to convergence to a cluster. The selection phase of the algorithm uses shooting method coupled with evaluation of Floquet multipliers to retain only those vectors in the arrived cluster yielding stable periodic solutions. The method is validated with Den Hartog's vibration absorbers and is then applied to vibration absorbers with material nonlinearity and vibration isolators with geometric nonlinearity. In both the cases, the converged cluster is shown to yield stable, periodic responses satisfying the amplitude condition. Parametric perturbation studies are conducted on the cluster to illustrate its robustness. The use of algebraic form of governing equations in the screening phase drastically reduces the computational time needed to converge to the cluster. The fact that the present method converges to a cluster in the parameter space rather than to a single parameter value offers the designer more freedom to choose the design vector from inside the cluster. It also ensures that the design is robust to small changes in parameter values.