Future ambitious solar system exploration missions are likely to require ever larger propulsion capabilities and involve innovative interplanetary trajectories in order to accommodate the increasingly complex mission scenarios. Two recent advances in trajectory design can be exploited to meet those new requirements: the use of low-thrust propulsion which enables larger cumulative momentum exchange relative to chemical propulsion; and the consideration of low-energy transfers relying on full multi-body dynamics. Yet the resulting optimal control problems are hypersensitive, time-consuming and extremely difficult to tackle with current optimization tools. Therefore, the goal of the thesis is to develop a methodology that facilitates and simplifies the solution finding process of low-thrust optimization problems in multi-body environments. Emphasis is placed on robust techniques to produce good solutions for a wide range of cases despite the strong nonlinearities of the problems. The complete trajectory is broken down into different component phases, which facilitates the modeling of the effects of multiple bodies and makes the process less sensitive to the initial guess. A unified optimization framework is created to solve the resulting multi-phase optimal control problems. Interfaces to state-of-the-art solvers SNOPT and IPOPT are included. In addition, a new, robust Hybrid Differential Dynamic Programming (HDDP) algorithm is developed. HDDP is based on differential dynamic programming, a proven robust second-order technique that relies on Bellman's Principle of Optimality and successive minimization of quadratic approximations. HDDP also incorporates nonlinear mathematical programming techniques to increase efficiency, and decouples the optimization from the dynamics using first- and second-order state transition matrices. Crucial to this optimization procedure is the generation of the sensitivities with respect to the variables of the system. In the context of trajectory optimization, these derivatives are often tedious and cumbersome to estimate analytically, especially when complex multi-body dynamics are considered. To produce a solution with minimal effort, an new approach is derived that computes automatically first- and high-order derivatives via multicomplex numbers. Another important aspect of the methodology is the representation of low-thrust trajectories by different dynamical models with varying degrees of fidelity. Emphasis is given on analytical expressions to speed up the optimization process. In particular, one novelty of the framework is the derivation and implementation of analytic expressions for motion subjected to Newtonian gravitation plus an additional constant inertial force. Example applications include low-thrust asteroid tour design, multiple flyby trajectories, and planetary inter-moon transfers. In the latter case, we generate good initial guesses using dynamical systems theory to exploit the chaotic nature of these multi-body systems. The developed optimization framework is then used to generate low-energy, inter-moon trajectories with multiple resonant gravity assists.