We use Bayesian optimization algorithms in combination with a nonequilibrium Green's function transport model to increase the maximum operating temperature of terahertz quantum cascade lasers. This procedure lead to the recent temperature record of 210 K in terahertz quantum cascade lasers, and here we provide even-further-improved structures. The Bayesian optimization algorithm, which takes into account all the available history of the optimization, converges much faster and more securely than the commonly used genetic algorithm. Designs based on two and three wells per period are considered, and using the large amount of data generated, we systematically evaluate the studied schemes in terms of optimal extraction energy and relevance of electron-electron correlations. This analysis shows that the two-well scheme is superior for reaching high operating temperatures, while the three-well scheme is more robust to variations in layer thicknesses. Furthermore, we study the sensitivity to flux-rate fluctuations during growth and simulation-model inaccuracies, showing the period thickness needs to be controlled to within a few percent, which is challenging but achievable with present-day molecular-beam epitaxy. These limits to the growth accuracy can be a guiding principle for experimentalists, along with the suggestion to fabricate devices across the wafer radius so as to find the optimal period thickness.