NOTE: In the below, I'm implicitly discussing a ground-state, closed-shell wavefunction, where all occupied orbitals are doubly occupied. The discussion would be similar for open-shell wavefunctions, but there are complexities that I won't address here. Also, once one starts noodling at excited states, things get complicated pretty quickly, and (AFAIK) some manipulations that are valid for a ground state wavefunction are invalid for excited states.
In my understanding localized molecular orbitals (LMO) are equivalent to "standard" molecular orbitals, often called canonical orbitals (btw. why are they called canonical?). We can produce LMOs by linear combination of MOs and they represent the same physical state.
The "canonical" MOs (CMOs) are called such because they are in a certain sense the "default" MOs for a wavefunction. The sense in which they are the "default" is that they diagonalize the Fock matrix that defines the energy of the system and is used in the course of the self-consistent field (SCF) calculation of, e.g., the Hartree-Fock method for solving the approximate many-electron problem.
Occupied localized MOs (LMOs), when all considered together as a unified wavefunction, in total contain the exact same information as the wavefunction composed of the occupied CMOs. This is because the LMOs are obtained from a unitary transformation of the occupied CMOs, and the linear algebra indicates that the overall wavefunction is invariant to such a transformation.
If LMOs are equivalent and have the same energy then why do calculations always end up with canonical molecular orbitals and never with LMOs?
The italicized phrase of this question is false. LMOs individually are not equivalent to individual CMOs, and they do not have the same energy. Strictly speaking, the energies of individual LMOs are actually undefined. Only the total wavefunction composed of LMOs is equivalent to the total wavefunction composed of CMOs.
As noted above, the reason why we always get CMOs out of an SCF energy calculation is that they are the set of orbitals naturally associated with the solution of the SCF problem.
And wouldn't there be an unlimited amount of equivalent solutions?
Yes. Every possible unitary transformation among the occupied canonical orbitals, of which there is indeed an unlimited number, would provide a set of orbitals that gives an equivalent total wavefunction. Not all of these transformed orbital sets are especially useful, however. Specific examples of useful ones include the various localization transforms that are the topic of the question, the series of variously defined "natural orbitals", and (somewhat related) the Löwdin symmetric orthogonalization for atomic orbitals.
Is this just a mathematical gimmick or does this represent the reality? ... Or is this just a thing we do because it makes it easier to understand?
The transformed MOs do represent reality, and often provide a very powerful and/or convenient means for analysis and understanding of the overall wavefunction. You do have to know what is or is not a valid use of any particular set of transformed orbitals, however.
One specific application of transformed orbitals I know about that seems very powerful is the family of local pair natural orbital (LPNO) methods for accelerating correlated calculations. The v4.0 release of ORCA includes implementations of their advanced domain-based (DLPNO) version of the method for several calculation methods (MP2, CASSCF-NEVPT2, etc.; see various of the 2016/2017 Neese citations listed here). The DLPNO-CCSD method has already proven to be extraordinarily powerful in enabling correlated calculations of large systems—the below figure is a reproduction of Figure 5.5 from the ORCA v3.0.3 manual showing the enhanced performance of (D)LPNO-CCSD methods for a series of benchmark systems (click to enlarge):
Is it possible for a molecule to be in such a state where every electron is in such a localized orbital?
Yes, this is always true, as long as the transformation that is applied to the occupied orbitals includes no unoccupied orbitals. No matter what set of orbitals you use, whether canonical or transformed, every occupied orbital in the wavefunction is always fully (here, doubly) occupied. Transformation among the occupied orbitals doesn't change the occupancy.
Inclusion of some un-occupied orbitals (often called virtual orbitals) in the transformation, though, will lead to orbital occupations between zero and two. This is what happens in the classic natural bond orbital analysis of Weinhold & co., as can be seen in the sample natural population analysis output (pdf link; "Occupancy" column) provided on the GENNBO website .