Lambert and Primer Vector Example

This topic illustrates the use of the LambertOrbitSolver to calculate delta-v vectors for maneuvers and the use of the a PropagationVector as a costate vector to target the directions of maneuvers with known delta-v magnitudes. Primer vector theory is applied to propagate the costate vector such that it points in the correct direction for both the initial and final maneuvers of a coplanar low Earth orbit (LEO) to geostationary Earth orbit (GEO) transfer.

Primer vectors were introduced by D. F. Lawden in "Optimal Trajectories for Space Navigation", published by Butterworths in 1963. Some of the equations and methodology used in this discussion are provided by J. E. Prussing in Chapter 2 "Primer Vector Theory and Applications" of "Spacecraft Trajectory Optimization", edited by B. A. Conway and published by Cambridge University Press in 2010.

Setup and Definition of Orbits

First, some preliminary setup needs to be completed. The reference frame, central body, central body radius, and central body gravitational parameter need to be defined to use the International Celestial Reference Frame (ICRF) and Earth-specific values. The initial date (or epoch) must be defined for the propagation definition to work, but its specific value is not particularly relevant to this application.

The LEO is defined as an equatorial circular orbit with an altitude of 300 km, and the GEO is defined using position and velocity vectors to begin 180 degrees away from the LEO starting conditions. The semi-major axis and time of flight of the geostationary transfer orbit (GTO) are defined to be consistent with an optimal Hohmann transfer from the LEO to the GEO.

Copy

// We want orbits with less than 1 complete revolution.
const int numberOfRevolutions = 0;

// Define central body parameters, reference frame, and initial date.
EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth;
ReferenceFrame icrfFrame = earth.InternationalCelestialReferenceFrame;
double gravitationalParameter = EarthGravitationalModel2008.GravitationalParameter; // 3.986004415e+14 m^3/s^2
double earthRadius = EarthGravitationalModel2008.SemimajorAxis; // 6378136.3 meters
JulianDate initialDate = new GregorianDate(2022, 3, 16).ToJulianDate();

// Parameters specific to initial LEO.
double initialSemimajorAxis = 300000.0 + earthRadius; // Circular orbit with 300 km altitude. 
double initialEccentricity = 0.0;
double initialInclination = 0.0;
double initialRightAscensionAscendingNode = 0.0;
double initialArgumentOfPerigee = 0.0;
double initialTrueAnomaly = 0.0;
Motion<Cartesian> initialLeo = new KeplerianElements(initialSemimajorAxis, initialEccentricity, initialInclination,
                                                     initialArgumentOfPerigee, initialRightAscensionAscendingNode,
                                                     initialTrueAnomaly, gravitationalParameter).ToCartesian();

// Parameters specific to final GEO.
Motion<Cartesian> finalGeo = new Motion<Cartesian>(new Cartesian(-42165000.0, 0.0, 0.0), // position vector, meters.
                                                   new Cartesian(0.0, -3074.6, 0.0)); // velocity vector, m/s.

// Time of flight of GTO transfer orbit.
double gtoSemimajorAxis = (initialSemimajorAxis + finalGeo.Value.Magnitude) / 2.0; // meters.
double gtoTimeOfFlight = 0.5 * OrbitalElements.ComputePeriod(1.0 / gtoSemimajorAxis, gravitationalParameter);

Solving for Initial Primer Vector Rate

The LambertOrbitSolver is used to compute the initial and final velocity vectors for the GTO. Because the GTO is a 180 degree transfer, there is an inherent ambiguity in the orbital plane of a transfer between the initial and final locations. The velocity vector of the initial LEO orbit is used to resolve that ambiguity and ensure a zero-inclination GTO.

The delta-v vector of the first maneuver is computed by subtracting the initial LEO velocity vector from the initial velocity vector of the Lambert solution. Similarly, the delta-v vector of the second maneuver is computed by subtracting the final velocity vector of the Lambert solution from the velocity vector of the desired GEO orbit.

The initial and final values of the primer vectors are found by normalizing the delta-v vectors for the first and second maneuvers. In most trajectory optimization applications, the initial primer vector values would have to be guessed and corrected. However, this simple example uses a Hohmann transfer that is already known to be optimal.

The initial value of the first derivative of the primer vector (also known as the primer vector rate) can be calculated from the initial and final primer vectors, but this calculation requires the state transition matrix (STM) to be computed first. A separate propagator is configured that allows the computation of the STM from the initial state to the final state of the GTO. The STM is a 6x6 matrix because both the position and velocity vectors are included. However, only two 3x3 sub-matrices of the STM are needed to compute the initial primer vector rate. (See equations 2.36, 2.37, 2.46, and 2.49 of Prussing for details.)

Copy

// Find solution using Lambert solver.
LambertOrbitSolver lambertSolver = new LambertOrbitSolver(gravitationalParameter);
LambertResult lambertSolution = lambertSolver.SolveFixedDurationTransfer(initialLeo.Value,
                                                                         finalGeo.Value,
                                                                         Duration.FromSeconds(gtoTimeOfFlight),
                                                                         numberOfRevolutions,
                                                                         LambertPathType.ShortPath,
                                                                         OrbitDirectionType.Prograde,
                                                                         initialLeo.FirstDerivative);

// Calculate delta-v's required to complete Lambert solution.
Cartesian deltaV1 = lambertSolution.InitialPositionMotion.FirstDerivative - initialLeo.FirstDerivative;
Cartesian deltaV2 = finalGeo.FirstDerivative - lambertSolution.FinalPositionMotion.FirstDerivative;

double deltaV1Magnitude;
double deltaV2Magnitude;
UnitCartesian initialPrimerVector = deltaV1.Normalize(out deltaV1Magnitude);
UnitCartesian finalPrimerVector = deltaV2.Normalize(out deltaV2Magnitude);

// Setup preliminary propagator to calculate state transition matrix.
NumericalPropagatorDefinition stmPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = initialDate,
};

double initialMass = 1000.0; // kg.

string massElementId = "Mass";
PropagationScalar mass = new PropagationScalar(initialMass)
{
    Identification = massElementId,
    ScalarDerivative = 0.0,
};

// Define states using PropagationNewtonianPoint.
const string id = "PrimaryMotion";
PropagationNewtonianPoint point = new PropagationNewtonianPoint
{
    Identification = id,
    Mass = mass.IntegrationValue,
    InitialPosition = initialLeo.Value,
    InitialVelocity = initialLeo.FirstDerivative + deltaV1,
    IntegrationFrame = icrfFrame,
};

// Configure two-body gravity.
TwoBodyGravity twoBodyGravity = new TwoBodyGravity(point.IntegrationPoint, earth, gravitationalParameter);
point.AppliedForces.Add(twoBodyGravity);

stmPropagatorDefinition.IntegrationElements.Add(mass);
stmPropagatorDefinition.IntegrationElements.Add(point);

// Create state transition matrix and add it as an integration element.
const string stmId = "StateTransitionMatrix";
IPartialDifferentiable position = (IPartialDifferentiable)point.IntegrationPoint;
IPartialDifferentiable velocity = (IPartialDifferentiable)point.IntegrationPoint.CreateVectorVelocity(point.IntegrationFrame);
// point is the IPartialDifferentiable for the acceleration.

StateTransitionMatrix stateMatrix = new StateTransitionMatrix()
{
    Identification = stmId,
};
stateMatrix.AddStateParameter(position, velocity);
stateMatrix.AddStateParameter(velocity, point);
stmPropagatorDefinition.IntegrationElements.Add(stateMatrix);

// Configure numerical integrator.
RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator
{
    RelativeTolerance = Constants.Epsilon13,
    MinimumStepSize = 1.0,
    MaximumStepSize = 6.0 * 3600.0,
    InitialStepSize = 60,
    AbsoluteTolerance = Constants.Epsilon13,
    StepSizeBehavior = KindOfStepSize.Relative,
};
stmPropagatorDefinition.Integrator = integrator;

// Propagate and extract results.
NumericalPropagator stmPropagator = stmPropagatorDefinition.CreatePropagator();
NumericalPropagationStateHistory stmPropagatorResults = stmPropagator.Propagate(Duration.FromSeconds(gtoTimeOfFlight), 1);
Matrix finalStm = stmPropagatorResults.GetDateMotionCollection<Matrix>(stmId).Values.Last();
// Extract dPosition/dPosition0 submatrix.
Matrix finalStmSubmatrix1 = finalStm.GetMatrix(0, 3, 0, 3);
Matrix3By3 convertedSubmatrix1 = new Matrix3By3(finalStmSubmatrix1[0, 0], finalStmSubmatrix1[0, 1], finalStmSubmatrix1[0, 2],
                                                finalStmSubmatrix1[1, 0], finalStmSubmatrix1[1, 1], finalStmSubmatrix1[1, 2],
                                                finalStmSubmatrix1[2, 0], finalStmSubmatrix1[2, 1], finalStmSubmatrix1[2, 2]);
// Extract dPosition/dVelocity0 submatrix.
Matrix finalStmSubmatrix2 = finalStm.GetMatrix(0, 3, 3, 3);
Matrix3By3 convertedSubmatrix2 = new Matrix3By3(finalStmSubmatrix2[0, 0], finalStmSubmatrix2[0, 1], finalStmSubmatrix2[0, 2],
                                                finalStmSubmatrix2[1, 0], finalStmSubmatrix2[1, 1], finalStmSubmatrix2[1, 2],
                                                finalStmSubmatrix2[2, 0], finalStmSubmatrix2[2, 1], finalStmSubmatrix2[2, 2]);

// Solve for initial primer vector rate.
Cartesian rightHandSide = finalPrimerVector.Subtract(convertedSubmatrix1 * initialPrimerVector);
Cartesian initialPrimerVectorRate = convertedSubmatrix2.Invert() * rightHandSide;

Main propagator and segment definitions

The primer vector is defined as a costate PropagationVector that uses the pre-computed initial primer vector and initial primer vector rates as its initial conditions. The VectorDerivative, which is the second derivative of the primer vector in this example, is defined by multiplying the instantaneous value of the gravity-gradient matrix of the Earth's two-body gravity by the instantaneous value of primer vector. Scalar arithmetic is used to describe these computations in order to delay the execution of these operations until the propagator actually runs.

Next, the main propagator is configured with the constant mass, spacecraft point, and primer vector as PropagationStateElements. The magnitudes of the maneuvers were calculated earlier in the code from the Lambert solutions. The directions of the maneuvers are defined by the primer vector values that occur immediately before each maneuver.

Finally, a SegmentList is configured that uses an initial state segment, a first maneuver segment, a coast segment, and a second maneuver segment.

Copy

// Now set-up main propagator that includes primer-vector propagation and maneuvers.
point.InitialVelocity = initialLeo.FirstDerivative; // Needs to be reset to no longer include deltaV1.

// Define PrimerVector as a costate PropagationVector
string primerVectorId = "PrimerVector";
PropagationVector primerVector = new PropagationVector(initialPrimerVector, initialPrimerVectorRate)
{
    Identification = primerVectorId,
    PropagationAxes = icrfFrame.Axes,
    // The VectorDerivative is the second derivative and is set in the following code.
};

// Get partial derivatives of the gravity.
Scalar radius = point.IntegrationPoint.GetScalarElement(CartesianElement.Magnitude, point.IntegrationFrame);
Scalar radiusX = point.IntegrationPoint.GetScalarElement(CartesianElement.X, point.IntegrationFrame);
Scalar radiusY = point.IntegrationPoint.GetScalarElement(CartesianElement.Y, point.IntegrationFrame);
Scalar radiusZ = point.IntegrationPoint.GetScalarElement(CartesianElement.Z, point.IntegrationFrame);
Scalar radiusSquared = radius * radius;
Scalar factor = gravitationalParameter / (radiusSquared * radius);

// Two-body gravity has simple partial derivatives. Even then, these operations only work because
// all of the states and co-states are in the same reference frame.
Scalar dGravityX_dX = factor * (3.0 * radiusX * radiusX / radiusSquared - 1);
Scalar dGravityY_dY = factor * (3.0 * radiusY * radiusY / radiusSquared - 1);
Scalar dGravityZ_dZ = factor * (3.0 * radiusZ * radiusZ / radiusSquared - 1);
Scalar dGravityX_dY = factor * 3.0 * radiusX * radiusY / radiusSquared; // Equal to dGravityY_dX
Scalar dGravityX_dZ = factor * 3.0 * radiusX * radiusZ / radiusSquared; // Equal to dGravityZ_dX
Scalar dGravityY_dZ = factor * 3.0 * radiusY * radiusZ / radiusSquared; // Equal to dGravityZ_dY

// Primer vector is multiplied by partials matrix to get derivatives.
Scalar primerVectorX = primerVector.IntegrationValue.GetScalarElement(CartesianElement.X, icrfFrame.Axes);
Scalar primerVectorY = primerVector.IntegrationValue.GetScalarElement(CartesianElement.Y, icrfFrame.Axes);
Scalar primerVectorZ = primerVector.IntegrationValue.GetScalarElement(CartesianElement.Z, icrfFrame.Axes);

Scalar primerSecondDerivativeX = primerVectorX * dGravityX_dX
                                 + primerVectorY * dGravityX_dY
                                 + primerVectorZ * dGravityX_dZ;
Scalar primerSecondDerivativeY = primerVectorX * dGravityX_dY
                                 + primerVectorY * dGravityY_dY
                                 + primerVectorZ * dGravityY_dZ;
Scalar primerSecondDerivativeZ = primerVectorX * dGravityX_dZ
                                 + primerVectorY * dGravityY_dZ
                                 + primerVectorZ * dGravityZ_dZ;
Vector primerSecondDerivative = primerSecondDerivativeX * VectorFixed.GetUnitXVector(icrfFrame.Axes)
                                + primerSecondDerivativeY * VectorFixed.GetUnitYVector(icrfFrame.Axes)
                                + primerSecondDerivativeZ * VectorFixed.GetUnitZVector(icrfFrame.Axes);

primerVector.VectorDerivative = primerSecondDerivative;

// Set-up main propagator.
NumericalPropagatorDefinition mainPropagatorDefinition = new NumericalPropagatorDefinition
{
    Epoch = initialDate,
    Integrator = integrator,
};
mainPropagatorDefinition.IntegrationElements.Add(mass);
mainPropagatorDefinition.IntegrationElements.Add(point);
mainPropagatorDefinition.IntegrationElements.Add(primerVector);

// Set-up segment propagation.
NumericalInitialStateSegment initial = new NumericalInitialStateSegment
{
    Name = "Initial State",
    PropagatorDefinition = mainPropagatorDefinition,
};

// To illustrate the use of primer vectors, the delta-v will be purely in the x-direction of the maneuver's orientation axes.
ImpulsiveManeuverSegment maneuver1 = new ImpulsiveManeuverSegment();
maneuver1.Maneuver = new ImpulsiveManeuverInformation(id, new Cartesian(deltaV1Magnitude, 0.0, 0.0));
maneuver1.Maneuver.CostateElementIdentification = primerVectorId;

// The maneuver's orientation axes are defined using the costate vector of the maneuver, which is the primer vector.
AxesAlignedConstrained axes1
    = new AxesAlignedConstrained(maneuver1.Maneuver.CostateVector,
                                 new VectorTrueDisplacement(earth.CenterOfMassPoint,
                                                            maneuver1.Maneuver.PropagationPoint));
maneuver1.Maneuver.Orientation = axes1;

// The GTO transfer orbit is defined next.
NumericalPropagatorSegment gtoSegment = new NumericalPropagatorSegment()
{
    Name = "GTO coast segment",
    PropagatorDefinition = mainPropagatorDefinition,
    MaximumDuration = Duration.FromDays(1.0),
};

// Stop after the GTO is finished.
DurationStoppingCondition durationStop = new DurationStoppingCondition(Duration.FromSeconds(gtoTimeOfFlight));
gtoSegment.StoppingConditions.Add(durationStop);

// To illustrate the use of primer vectors, the delta-v will be purely in the x-direction of the maneuver's orientation axes.
ImpulsiveManeuverSegment maneuver2 = new ImpulsiveManeuverSegment();
maneuver2.Maneuver = new ImpulsiveManeuverInformation(id, new Cartesian(deltaV2Magnitude, 0.0, 0.0));
maneuver2.Maneuver.CostateElementIdentification = primerVectorId;

// The maneuver's orientation axes are defined using the costate vector of the maneuver, which is the primer vector.
AxesAlignedConstrained axes2
    = new AxesAlignedConstrained(maneuver2.Maneuver.CostateVector,
                                 new VectorTrueDisplacement(earth.CenterOfMassPoint,
                                                            maneuver2.Maneuver.PropagationPoint));
maneuver2.Maneuver.Orientation = axes2;

// Make a segment list, populate the segment list, and propagate the segment list.
SegmentList mainSegmentList = new SegmentList();
mainSegmentList.Segments.Add(initial);
mainSegmentList.Segments.Add(maneuver1);
mainSegmentList.Segments.Add(gtoSegment);
mainSegmentList.Segments.Add(maneuver2);

Propagation and Results

A SegmentListPropagator is generated from the SegmentList and then propagated to get results. The results for each of the four inner segments are extracted from the SegmentListResults.

The initial and final values of the position, velocity, and costate primer vectors of the spacecraft are extracted from the results of each relevant segment.

Copy

// Define propagator and propagate it.
SegmentListPropagator mainPropagator = mainSegmentList.GetSegmentListPropagator();
SegmentListResults mainResults = mainPropagator.PropagateSegmentList();

// Extract results for each segment inside segment list.
SegmentResults initialSegmentResults = mainResults.GetResultsOfSegment(initial);
ImpulsiveManeuverSegmentResults maneuver1Results
    = (ImpulsiveManeuverSegmentResults)mainResults.GetResultsOfSegment(maneuver1);
SegmentResults gtoSegmentResults = mainResults.GetResultsOfSegment(gtoSegment);
ImpulsiveManeuverSegmentResults maneuver2Results
    = (ImpulsiveManeuverSegmentResults)mainResults.GetResultsOfSegment(maneuver2);

Motion<Cartesian> propagationInitialState = initialSegmentResults.FirstPropagatedState.GetMotion<Cartesian>(id);

Motion<Cartesian> maneuver1InitialState = maneuver1Results.FirstPropagatedState.GetMotion<Cartesian>(id);
Motion<Cartesian> maneuver1FinalState = maneuver1Results.FinalPropagatedState.GetMotion<Cartesian>(id);

Motion<Cartesian> gtoInitialState = gtoSegmentResults.FirstPropagatedState.GetMotion<Cartesian>(id);
Motion<Cartesian> gtoFinalState = gtoSegmentResults.FinalPropagatedState.GetMotion<Cartesian>(id);

Motion<Cartesian> gtoInitialCostate = gtoSegmentResults.FirstPropagatedState.GetMotion<Cartesian>(primerVectorId);
Motion<Cartesian> gtoFinalCostate = gtoSegmentResults.FinalPropagatedState.GetMotion<Cartesian>(primerVectorId);

Motion<Cartesian> maneuver2InitialState = maneuver2Results.FirstPropagatedState.GetMotion<Cartesian>(id);
Motion<Cartesian> maneuver2FinalState = maneuver2Results.FinalPropagatedState.GetMotion<Cartesian>(id);