Organic photovoltaic devices are promising candidates for efficient energy harvesting from sunlight. Designing new dye molecules suitable for such devices is a challenging task restricted by the rapid increase of computational cost with system size. Solar cell material properties are closely related to the electronic structure of the dye, and an effective molecular orbital energy screening method for a family of dyes is therefore desired. In this work, a machine learning approach is used to sort through the chemical space of peripheral double-substituted boron-Subphthalocyanine dyes. A database of 12,102 PM6 optimized structures was built and for each of the structures time-dependent density functional theory (LC-omega HPBE/6-31+G(d)) calculations were performed. We investigated the changes of the molecular orbital energies of the molecular orbitals related to reduction and oxidation of the compounds. With the Electrotopological-state index moleculear representation all the tested algorithms, Support Vector Machine, Random Forest Regression, Neural Network, and Simple Linear Regression, captured the calculated frontier orbital energies with a prediction root-mean-square-error in the order of 0.05 eV. Finally, frontier orbital energies were predicted for more than 40,000 new structures by the trained Support Vector Machine algorithm. Compared to the parent boron-Subphthalocyanine structure, 237 and 132 functionalized dyes were predicted to have upshifted molecular orbital energies using the Electrotopological-state index and OneHot encoding feature vector, respectively. Out of 27 investigated donor and acceptor ligands, the acetamide and hydroxyl ligands gave rise to the desired increase in frontier molecular orbital energy.