Wannier functions provide a localized representation of spectral subspaces of periodic Hamiltonians, and play an important role for interpreting and accelerating Hartree-Fock and Kohn-Sham density functional theory calculations in quantum physics and chemistry. For systems with isolated band structure, the existence of exponentially localized Wannier functions and numerical algorithms for finding them are well studied. In contrast, for systems with entangled band structure, Wannier functions must be generalized to span a subspace larger than the spectral subspace of interest to achieve favorable spatial locality. In this setting, little is known about the theoretical properties of these Wannier functions, and few algorithms can find them robustly. We develop a variational formulation to compute these generalized maximally localized Wannier functions. When paired with an initial guess based on the selected columns of the density matrix (SCDM) method, our method can robustly find Wannier functions for systems with entangled band structure. We formulate the problem as a constrained nonlinear optimization problem, and show how the widely used disentanglement procedure can be interpreted as a splitting method to approximately solve this problem. We demonstrate the performance of our method using real materials including silicon, copper, and aluminum. To examine more precisely the localization properties of Wannier functions, we study the free electron gas in one and two dimensions, where we show that the maximally-localized Wannier functions only decay algebraically. We also explain using a one dimensional example how to modify them to obtain super-algebraic decay.