We perform a general optimization of the parameters in the Multilevel Monte Carlo (MLMC) discretization hierarchy based on uniform discretization methods with general approximation orders and computational costs. Moreover, we discuss extensions to non-uniform discretizations based on a priori refinements and the effect of imposing constraints on the largest and/or smallest mesh sizes. We optimize geometric and non-geometric hierarchies and compare them to each other, concluding that the geometric hierarchies, when optimized, are nearly optimal and have the same asymptotic computational complexity. We discuss how enforcing domain constraints on parameters of MLMC hierarchies affects the optimality of these hierarchies. These domain constraints include an upper and lower bound on the mesh size or enforcing that the number of samples and the number of discretization elements are integers. We also discuss the optimal tolerance splitting between the bias and the statistical error contributions and its asymptotic behavior. To provide numerical grounds for our theoretical results, we apply these optimized hierarchies together with the Continuation MLMC Algorithm that we recently developed, to several examples. These include the approximation of three-dimensional elliptic partial differential equations with random inputs based on FEM with either direct or iterative solvers and It\^o stochastic differential equations based on the Milstein scheme.