This thesis presents an optimization architecture for on-orbit servicing mission design in the long-range rendezvous phase. We develop a methodology to generate Pareto Optimal trajectories for long-range rendezvous of a servicing satellite with a moving target. The methodology employs a multi-impulse shape-based trajectory planning algorithm for in-plane orbit transfer, based on the two-body problem. We first derive the necessary and sufficient conditions that determine the set of smooth impulsive trajectories connecting the servicing satellite to the orbiting target. The Pareto Optimal trajectories from this set are then obtained using a constrained multiobjective optimization algorithm developed based on the Non-dominated Sorting Genetic Algorithm-II (NSGA-II). Transfer time and control effort are the two Pareto cost functions that are considered in the multi-objective optimization. To reduce the risk of collision in populated orbits and to remain in an orbital regime, we include restrictions on orbital elements as part of the constraints. Further, a maximum available impulse is considered as an upper-bound for velocity changes in an impulsive trajectory. The number of impulses along with the location of the first impulse in the parking orbit and the orbital parameters of the intermediate orbits form the set of design variables. We demonstrate the superiority of the developed trajectory planner by comparing its results with those obtained from another multi-objective evolutionary algorithm called the Multi-Objective Genetic Algorithm and an optimal Lambert approach. In an on-orbit servicing mission, a solution from the Pareto frontier set of optimal trajectories may be selected based on the mission requirements.To robustly follow the generated reference trajectories in a J 2 -perturbed orbital environment, we propose a Nonlinear Model Predictive Control (NMPC) scheme. The control signals are velocity increments at the time of applying each impulse, and the variable horizon is considered to be the time difference between every two impulses in the reference trajectory. To avoid singularities, the equinoctial orbital elements are used in the process model that includes the secular and first-order long-periodic effects iii