Maximally-localized Wannier functions are quantum wavefunctions resembling atomic orbitals that are used to describe electrons in condensed matter1. Since their introduction in 19972, these functions have become ubiquitous in ab initio materials simulations, including applications in linear-scaling methods3, strongly correlated electron systems4, quantum transport5, electron-phonon interactions6, and topological materials7. Despite their widespread adoption in a vast software ecosystem8, Wannier functions have not yet attained their fullest potential in the presence of entangled bands, as their optimization remains challenging and labor-intensive. Here, we introduce a universal meta-optimization method that leverages workflow abstraction and machine-learning techniques like differential evolution and Bayesian optimization to generate globally optimized Wannier functions without human intervention. We demonstrate this approach through three applications: (i) autonomous interpolation of entangled band structures with millielectronvolt accuracy starting from coarse Brillouin zone grids, (ii) thousand-fold acceleration of fully ab initio Boltzmann transport calculations via the use of minimal coarse Brillouin zone grids, and (iii) ultrafast high-throughput calculations of high-precision Wannier functions for large materials libraries. This work brings calculations that previously required supercomputers within the reach of personal computers.