In this paper we present a self-consistent ensemble optimization (SCEO) theory for efficient conformational search, which we have applied to predicting the effects of mutations on protein thermostability. This approach takes advantage of a statistical mechanical self-consistency condition to home in iteratively on the global minimum structure. We employ a fast potential of mean-force approximation to cut computation time to a few minutes for a typical protein mutation, with only linear time-dependence on the size of the prediction problem. Rather than seeking a single, static structure of minimum energy, the new method optimizes an ensemble of many conformations, seeking to predict the most likely ensemble for the native state at a desired temperature. Testing this approach with a simple physical model focusing entirely on steric interactions and side-chain rearrangement, we obtain robustly convergent prediction of core side-chain conformation, and of hydrophobic core mutations' effects on protein stability. Self-consistent ensemble optimization is superior to simulated annealing in its speed and convergence to the global minimum, and insensitive to starting conformation. In calculations on lambda repressor protein, structural predictions for an eight-residue molten-zone had side-chain r.m.s. error of 0.49 A for the wild-type protein. Evaluation of the method's mutant structure predictions should become possible, as structures of these mutant repressors are solved. Predicted energies for a series of nine hydrophobic core mutants correlated with measured free energies of unfolding with a coefficient of 0.82.