With the pursuit of increasingly innovative and complex space missions, the focus of the space industry has been turning towards electric propulsion systems. Due to their high specific impulse - about ten times that of a chemical engine - they provide large savings in propellant mass, decreasing the overall cost of the mission. This proves to be essential for small low cost missions, such as interplanetary CubeSats, and more ambitious endeavours such as asteroid retrieval or crewed missions to Mars. Designing a low-thrust trajectory is a more complex task than doing so for a high-thrust one, since computing the thrust sequence that minimizes the fuel spent requires a search over a huge and complex design space. Setting up the optimal control problem generally requires a good first-guess solution, a fine tuning of the parameters involved, and the definition of feasible bounds for the trajectory. In order to converge to a solution, the problem settings are simplified as much as possible. This includes the dynamical framework used, which often may not be sensitive enough to describe the low-energy trajectory regime necessary for some of the mission examples mentioned above. This abstract proposes a new set of equations of motion to solve the optimal control problem. These are derived from the disturbing function of the previously studied Keplerian Map, formulated from the Hamiltonian of the CR3BP. Its motion corresponds to the propagation of Gauss’s planetary equations with both the disturbing potential of the CR3BP, and the accelerations of the electric engine. The novelty of this formulation is that it describes a third-body motion in terms of the orbital elements that define the osculating orbit of the spacecraft, in a barycentric coordinate system. This is advantageous in several respects: first, low-thrust sub-optimal control laws can be easily generated and explored to find a first guess solution near global optima. Second, bounds for the optimal control problem, as well as the boundary values, can be easily defined, which allows for a much faster convergence. This dynamical framework is accurate until very close to the sphere of influence of the perturbing body, and thus can be efficiently used to target low-energy hyperbolic invariant manifold structures associated with periodic orbits near it. The paper presents the methodology as well as a full retrieval trajectory for asteroid 2018 AV2, a small co-orbital asteroid that could be retrieved during its next Earth encounter in 2037.