Accurate computation of singlet-triplet energy gaps of diradicals remains a challenging problem in density-functional theory (DFT). In this work, we propose a variational extension of our previous work [D. H. Ess, E. R. Johnson, X. Q. Hu, and W. T. Yang, J. Phys. Chem. A 115, 76 (2011)], which applied fractional-spin density-functional theory (FS-DFT) to diradicals. The original FS-DFT approach assumed equal spin-orbital occupancies of 0.5 α-spin and 0.5 β-spin for the two degenerate, or nearly degenerate, frontier orbitals. In contrast, the variational approach (VFS-DFT) optimizes the total energy of a singlet diradical with respect to the frontier-orbital occupation numbers, based on a full configuration-interaction picture. It is found that the optimal occupation numbers are exactly 0.5 α-spin and 0.5 β-spin for diradicals such as O(2), where the frontier orbitals belong to the same multidimensional irreducible representation, and VFS-DFT reduces to FS-DFT for these cases. However, for diradicals where the frontier orbitals do not belong to the same irreducible representation, the optimal occupation numbers can vary between 0 and 1. Furthermore, analysis of CH(2) by VFS-DFT and FS-DFT captures the (1)A(1) and (1)B(1) states, respectively. Finally, because of the static correlation error in commonly used density functional approximations, both VFS-DFT and FS-DFT calculations significantly overestimate the singlet-triplet energy gaps for disjoint diradicals, such as cyclobutadiene, in which the frontier orbitals are confined to separate atomic centers.