Cyclically periodic structures, such as bladed disk assemblies in turbomachinery, are widely used in engineering systems. It is well known that small uncertainties exist among their substructures, which in certain situations may cause drastic change in the dynamic responses, a phenomenon known as vibration localization. Previous studies have suggested that the introduction of small, prespecified design modification, i.e., intentional mistuning, may alleviate vibration localization and reduce response variation. However, there has been no systematic methodology to facilitate the optimal design of intentional mistuning. The most significant challenge is the computational cost involved. The finite-element model of a bladed disk usually requires a very large number of degrees-of-freedom (DOFs). When uncertainties occur in a cyclically periodic structure, the response may no longer be considered as simple perturbation to that of the nominal structure. In this research, a suite of interrelated algorithms is proposed to enable the efficient design optimization of cyclically periodic structures toward alleviating their forced response variation. We first integrate model order reduction with a perturbation scheme to reduce the scale of analysis of a single run. Then, as the core of the new methodology, we incorporate Gaussian process (GP) emulation to conduct the rapid sampling-based evaluation of the design objective, which is a metric of response variation under uncertainties, in the parametric space. The optimal design modification can thus be directly identified to minimize the response variation. The efficiency and effectiveness of the proposed methodology are demonstrated by systematic case studies.