Computational modeling of disordered crystal structures is essential for the study of composition–structure–property relations for many families of functional materials. Efficient and reliable solutions are of great interest due to the promising reduction of labor-intensive calculations through data-driven prediction of target properties. One of the modern problems directly related to this topic is the enhancement of phase stability of functional materials with the required properties. In this work, we address the effects of Cd and Zn substitutions on the Pb sites on the structural stability of inorganic lead halide perovskite CsPbI3 — a semiconductor promising for optoelectronic devices. At room temperature, this material undergoes an undesirable phase transition from the direct band gap “black” to indirect “yellow” (ca. 15 meV/atom more favored) phase. For evaluations of effects caused by the substitutions, we combined the density functional theory (DFT) calculations and graph neural network (GNN) models. At low Cd and Zn contents, complete composition/configuration spaces comprising ca. 200 structures were set and studied using DFT. The GNN models trained on this data were used to predict energetics of ca. 74k structures with higher contents of substitutions and evaluate the energy differences between the phases. Because of the small amount of DFT data, an impact of data subsampling on the quality of predictions was comprehensively studied. Additionally, transfer learning routines were tested using available collections of DFT-derived properties. For the examined models, the root-mean-square errors of ca. 0.7 meV/atom were achieved in predicting formation energies. The data treatment workflows proposed and tested in this study open up perspectives for data-driven research of disordered materials of particular interest with low computational costs.