If the energy of a system is actually constrained (as it generally is), then we should be applying a second constraint, besides the one that allows us to normalize our probabilities. $$\mathcal{L} = -k_B\sum_iP_i\ln P_i + \alpha k_B\left(1-\sum_i P_i\right) + \beta k_B \left(U - \sum_i P_i E_i\right)$$ where $\alpha$ and $\beta$ are the two Lagrange multipliers. We want to maximize this, so we set its derivatives to zero: $$\frac{\partial\mathcal{L}}{\partial P_i} = 0 $$ $$= -k_B\left(\ln P_i + 1\right) - k_B\alpha - \beta k_B E_i$$ $$\ln P_i = -1 -\alpha - \beta E_i$$ At this point, it is convenient to invoke the normalization constraint… $$\sum_i P_i = 1$$ $$1= \sum_i e^{-1-\alpha-\beta E_i}$$ $$1= e^{-1-\alpha}\sum_i e^{-\beta E_i}$$ $$e^{1+\alpha} = \sum_i e^{-\beta E_i}$$ $$Z \equiv \sum_i^\text{all states} e^{-\beta \epsilon_i}$$ $$P_i = \frac{e^{-\beta \epsilon_i}}{Z}$$ $$P_i = \frac{Boltzmann factor}{partition function}$$ At this point, we haven't yet solved for $\alpha$, and to do so, we'd need to invoke the internal energy constraint: $$U = \sum_i E_i P_i$$ $$U = \frac{\sum_i E_i e^{-\beta E_i}}{Z}\label{eq:UfromZ}$$