Electric low-thrust propulsion has nowadays found wide application in space dynamics as it entails considerable savings in spacecraft propellant mass, thanks to the very high specific impulse that this kind of engine is able to generate. However, continuous thrust opens new extensive sets of feasible trajectories, and optimization algorithms are needed to mine the admissible search space and find optimal transfers. Practical methods to solve complex optimal control problems as low-thrust trajectory optimization typically involve high computational times, the major bottleneck of these techniques. The objective of this work is the development of a novel multiple-shooting optimization tool, employing a variational approach for quick derivative computation, and the assessment of its performance against a variety of test cases. Indeed, after a first analysis of practical optimization methods, the derivative estimation by finite-difference approximations has been found as the major contributor to the computational burden, and the propagation of the variational dynamics has been selected as an accurate approach to speed-up their computation. The theory of variational dynamics for multiple-shooting application has been analyzed in detail, and further developed for what concerns the second-order equations. After practical considerations on the method implementation (scaling procedures, sparsity patterns, et cetera) and its interface with WORHP, the selected non-linear programming solver, the tool has been applied to a broad range of test cases, spanning from elementary problems to practical applications. For what concerns the latter cases, two complex problems were analyzed and optimized: a CubeSat rendezvous departing from Earth-Moon L2 and arriving at asteroid 2000SG334, resulting in a propellant mass convenient trajectory suitable for asteroid reconnaissance well before a proposed NASA manned mission in 2069; The Kessler Run, i.e. the 9th Global Trajectory Optimization Competition (GTOC9), in which the developed tool, employed as last step of the optimization cascade of Strathclyde++ team, managed to make the solution constraint-feasible, valid and further mass-optimal.