Undulatory locomotion is common to nematodes as well as to limbless vertebrates, but its control is not understood in spite of the identification of hundred of genes involved in Caenorhabditis elegans locomotion. To reveal the mechanisms of nematode undulatory locomotion, we quantitatively analysed the movement of C. elegans with genetic perturbations to neurons, muscles, and skeleton (cuticle). We also compared locomotion of different Caenorhabditis species. We constructed a theoretical model that combines mechanics and biophysics, and that is constrained by the observations of propulsion and muscular velocities, as well as wavelength and amplitude of undulations. We find that normalized wavelength is a conserved quantity among wild-type C. elegans individuals, across mutants, and across different species. The velocity of forward propulsion scales linearly with the velocity of the muscular wave and the corresponding slope is also a conserved quantity and almost optimal; the exceptions are in some mutants affecting cuticle structure. In theoretical terms, the optimality of the slope is equivalent to the exact balance between muscular and visco-elastic body reaction bending moments. We find that the amplitude and frequency of undulations are inversely correlated and provide a theoretical explanation for this fact. These experimental results are valid both for young adults and for all larval stages of wild-type C. elegans. In particular, during development, the amplitude scales linearly with the wavelength, consistent with our theory. We also investigated the influence of substrate firmness on motion parameters, and found that it does not affect the above invariants. In general, our biomechanical model can explain the observed robustness of the mechanisms controlling nematode undulatory locomotion.