Many low earth satellites use GNSS for orbit determination, both for operational purposes and post-facto scientific orbits. Using GNSS systems for Lunar navigation proves difficult, however, due to the signals being weak and having limited coverage. This thesis investigates, through numerical simulation, the usage of a constellation of Polar or Equatorial Medium Earth Orbit (PECMEO) satellites for Lunar navigation. The system is proposed as alternative for orbit determination using ground-based observations, making navigation of the spacecraft more autonomous, without limitations on the number of receivers. A navigation constellation of 9 satellites at an orbital radius of 14 000 km is designed by minimizing the Geometric Dilution of Precision (GDOP). The resulting Pareto front of optimum solutions consists of two distinct groups of constellations. Firstly, constellation with lower mean GDOP and higher maximum GDOP have two orbits approximately orthogonal to the ecliptic plane, while the third orbit is in the ecliptic plane. Meanwhile, a second group of solutions with slightly higher mean GDOP, but significantly lower maximum GDOP are found for constellation with one plane approximately orthogonal to the ecliptic plane, while the remaining two planes have an inclination of near 45 degrees. One of the latter group of solutions has been chosen for the system for its lower maximum GDOP. The effect of omitting observations passing through the ionosphere is investigated to determine if there is a need to obtain observations at multiple frequencies to allow for ionospheric corrections. Omitting those observations results in a slightly increased GDOP right before and after a satellite passes behind Earth, of up to 56.4%. However, only 13.9% of the samples are effected, making the overall effect rather small. The positioning accuracy of the system is assessed by simulating pseudorange and carrier phase observations, and solving these with point position and kinematic least squares methods. Two main error contributions on the solutions are caused by observation noise and navigation satellite ephemerides errors. A solution for observations without ephemeris error, a mean 3-dimensional position error of 3.3m is obtained. A solution for noiseless observations yields a mean 3-dimensional position error of 26.8m. Finally, the solution obtained from observations with all simulated errors results in a mean 3-dimensional position error of 27.0m, to which the ephemeris errors are thus the main contributors, despite optimistic assumptions on their magnitude. The observed accuracy shows that the system is a viable navigation method for Lunar spacecraft.