One of the most interesting and challenging aspects of formation guidance algorithm design is the coupling of the orbit design and the science return. The effectiveness of the formation as a science instrument is intimately coupled with the relative geometry and evolution of the collection of spacecraft. Therefore, the science return can be maximized by optimizing the orbit design according to a performance metric relevant to the science mission goals. In this work, we present a general method for optimal formation guidance that is applicable to missions whose performance metric, requirements, and constraints can be cast as functions that are explicitly dependent upon the orbit states and spacecraft relative positions and velocities. The approach is fully nonlinear, accommodates orbital perturbations, and is applicable to multiple flight regimes including highly eccentric, hyperbolic, interplanetary, or libration point orbits. Furthermore, the method is applicable to small formations, as well as large formations and constellations. We present a general form for the cost and constraint functions, and derive their semi-analytic gradients with respect to the formation initial conditions. The gradients are broken down into two types. The first type are gradients of the mission-specific performance metric with respect to formation geometry. The second type are derivatives of the formation geometry with respect to the orbit initial conditions. The fact that these two types of derivatives appear separately allows us to derive and implement a general framework that requires minimal modification to be applied to different missions or mission phases. To illustrate the applicability of the approach, we conclude with applications to two missions: the Magnetospheric Multiscale Mission, and the Laser Interferometer Space Antenna. Nomenclature A = upper left 3 3 partition of B = upper right 3 3 partition of C = lower left 3 3 partition of c k = quadrature constant at point k D = lower right 3 3 partition of m = no. of unique sides in the formation N = no. of path constraints n = no. of spacecraft n k = no. of quadrature points r = position vector s ik = vector defining side i at quadrature point k _ s ik = rate vector of side i at quadrature point k s ik = length of side i at quadrature point k _ s ik = rate of change in length of side i at quadrature point k v = velocity vector = orbit state transition matrix Subscripts C = indicates association with constraints i = side index J = indicates association with cost function j = spacecraft index k = quadrature point index ' = spacecraft index p = dummy index