The emergence of cluster synchronization, as a universal phenomenon in complex systems, depends on the underlying networks and corresponding dynamics. Modifying network structures to achieve a desired cluster synchronization is a fundamental yet challenging task. Here, we address this problem by conducting a dimensional reduction in order to analyze the cluster stability of coupled phase oscillators. In particular, we exploit an optimization procedure to obtain a desired cluster synchronization pattern and analyze its stability and robustness based on Dirac delta perturbations. We find that altering the cluster synchronization pattern enhances the stability in comparison with the initial network configuration. Furthermore, we define the energy functions quantify the corresponding robustness, which rely on the phase differences between nodes in same clusters (intra-cluster) and nodes in different clusters (inter-cluster). The results show that, modifying cluster synchronization pattern has little impact on the intra-cluster energy variations, yet the inter-cluster energies monotonically increase with the number of clusters. Consequently, the global energy variations increase with the number of clusters, i.e., the higher the number of clusters, the more energy is spent to achieve the cluster synchronization pattern. Our approach to obtain a desired cluster synchronization pattern minimally modifies the network connections, while the corresponding system may not exhibit optimal stability. This work has implications for the development of optimal procedures to obtain a target cluster synchronization patterns with optimal stability.