<p>Maximally-localized Wannier functions are quantum wavefunctions resembling atomic orbitals that are used to describe electrons in condensed matter. Since their introduction in 1997, these functions have become ubiquitous in ab initio materials simulations, including applications in linear-scaling methods, strongly correlated electron systems, quantum transport, electron-phonon interactions, and topological materials. Despite their widespread adoption in a vast software ecosystem, 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.<br><br>The data generated in this work and reported in the manuscript is attached in the Data.tar.gz with full description in the Readme.md. </p>