Click or drag to resize

Multivariable Function Solvers

Function solvers allow a system of multivariable equations to be driven to some set of values. DME Component Libraries includes one such solver type, the NewtonRaphsonMultivariableFunctionSolver. This solver will sequentially evaluate the objective function at the current variable vector, compute the derivative of the function, and iterate towards the solution. This process will be repeated until the difference between the computed function value and the constraint's target values are equal within tolerance. To assist with convergence the input variables have several settings such as, MaximumStep (get / set), VariableTolerance (get / set), and, Scaling (get / set). For more information on the available settings see the settings class.

Note Note

The functionality described in this topic requires a license for the Segment Propagation Library.

Newton Raphson Function Solver

Suppose you want to solve the following pair of equations for x and y:

x2 + y2 = 25

x + 3.5 sin(x) = y

First it would help to get both equations to be implicit equations. The first equation already has the constant on one side, but the second should be rewritten as:

x - y + 3.5 sin(x) = 0

Eventually, we will need to define a SolvableMultivariableFunction, but first we start by building up the variables and constraints. The following code sample creates two variables, one for x and one for y:

Java
SolverVariableSettings xVariable = new SolverVariableSettings(0.3); // maximum step
// initial value is left at the default 0
xVariable.setName("x");

SolverVariableSettings yVariable = new SolverVariableSettings(0.3);
yVariable.setName("y");

Next, we create two constraints, each corresponding to a value computed by the function. These define the criteria that the function will be solved for. In this case, the constraints are the answer for the individual equations, with the desired value corresponding to the answer of each of the equations.

Java
SolverConstraintSettings firstEquation = new SolverConstraintSettings(25.0, // desired value
                                                                      0.001); // tolerance
firstEquation.setName("firstEquation");

SolverConstraintSettings secondEquation = new SolverConstraintSettings(0.0, 0.001);
secondEquation.setName("secondEquation");

These variable and constraint settings are then added to the function solver, as seen in the code sample below. The order in which the variables and constraint settings are added to the function solver must match the order that the variables are used and the constraints are computed in the evaluate method of the function.

Java
NewtonRaphsonMultivariableFunctionSolver solver = new NewtonRaphsonMultivariableFunctionSolver();
solver.getConstraints().add(firstEquation);
solver.getConstraints().add(secondEquation);
solver.getVariables().add(xVariable);
solver.getVariables().add(yVariable);

With the variables and the constraints created, we now define the function. The function must derive from the SolvableMultivariableFunction abstract base class, and implement the evaluate method. This method will be given an array of the current values of the variables, as computed by the NewtonRaphsonMultivariableFunctionSolver. Then, it must compute each equation, and assign those values to the results, in the same order as the constraints on the NewtonRaphsonMultivariableFunctionSolver. The following code sample shows the implementation of the evaluate method:

Java
@Override
public SolvableMultivariableFunctionResults evaluate(double[] variables, ITrackCalculationProgress progressTracker) {
    double[] equationAnswers = new double[2];
    double x = variables[0]; // x is always the first variable in this type
    double y = variables[1]; // y is always the second variable in this type

    // the order of the constraint values must match the order of the constraints in this type
    equationAnswers[0] = x * x + y * y;
    equationAnswers[1] = x - y + 3.5 * Math.sin(x);

    return new SolvableMultivariableFunctionResults(variables, equationAnswers);
}

public final void setPerturbationValues(double xPerturbation, double yPerturbation) {
    setPerturbationValues(new double[] {
        xPerturbation,
        yPerturbation
    });
}

The following code shows how to create and run the NewtonRaphsonMultivariableFunctionSolver:

Java
// SimpleFunction is our example class derived from SolvableMultivariableFunction
SimpleFunction function = new SimpleFunction();
function.setPerturbationValues(0.1, 0.1);
solver.setFunction(function);

solver.findSolution(25);
MultivariableFunctionSolverResults<MultivariableFunctionSolverIterationResults> results = solver.getLastRunsResults();

double[] variablesUsed = results.getFinalIteration().getFunctionResult().getVariablesUsed();
double foundX = variablesUsed[0];
double foundY = variablesUsed[1];
// 1.3744 4.8073

You can also check the variables directly to see that x = 1.37448 and y = 4.8073. However, if you graph the functions, you will see that there are multiple solutions:

Multivariable Function Plot

The initial guess will influence the results that the solver finds. Because of this fact, if you change your initial guess, it will find different answers:

Java
// recreate the variables with new initial values
xVariable = new SolverVariableSettings(0.3, 5.0, 0.001);
yVariable = new SolverVariableSettings(0.3, 5.0, 0.001);

function = new SimpleFunction();
function.setPerturbationValues(0.1, 0.1);

// recreate the solver
solver = new NewtonRaphsonMultivariableFunctionSolver();
solver.getVariables().add(xVariable);
solver.getVariables().add(yVariable);
solver.getConstraints().add(firstEquation);
solver.getConstraints().add(secondEquation);
solver.setFunction(function);
solver.findSolution(25);

results = solver.getLastRunsResults();
variablesUsed = results.getFinalIteration().getFunctionResult().getVariablesUsed();

foundX = variablesUsed[0];
foundY = variablesUsed[1];
// 4.8178 1.3372

One limitation of the Newton Raphson method is that it assumes a linear function. Consider if the slope of the function was near 0 when one of the variable is perturbed. The variable would be asked to take a large step that, if taken, could take it well beyond the parts of the function you are interested in. The MaximumStep (get / set) property on the variables prevents those large jumps.

Segmented Propagation and Function Solvers

The function you are solving can be very abstract. For example, you can set up a list of segments to compute a spacecraft's ephemeris, compute some value based on that ephemeris (such as the eccentricity and inclination at a point in time), and set up a variable that modifies the delta-v of some maneuver. In effect, you have the following equation, which can be run through the TargetedSegmentListDifferentialCorrector:

f(delta-v-x, delta-v-z) = [inclination, eccentricity]

Targeted Segment List Differential Corrector

TargetedSegmentListDifferentialCorrector requires that at least one SegmentPropagatorConstraint and one SegmentPropagatorVariable are configured on it, as shown below:

Java
// create the differential corrector
TargetedSegmentListDifferentialCorrector differentialCorrector = new TargetedSegmentListDifferentialCorrector();
differentialCorrector.setName("Corrector");

// add the differential corrector to the TargetedSegmentList
targetedSegmentList.getOperators().add(differentialCorrector);

// create a variable
DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable = impulsiveManeuverSegment.createVariable(
        200.0, // maximum step, meters/second
        1.0,   // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));
velocityXVariable.setName("Velocity_X_Variable");

// add the variable to the differential corrector
differentialCorrector.getVariables().add(velocityXVariable);

// create a constraint
DelegateBasedConstraint periodConstraint = new DelegateBasedConstraint(
        DelegateBasedConstraintCallback.of(segmentResults -> {
            Motion1<Cartesian> objectsMotion = segmentResults.getStateForNextSegment().getMotion(SatelliteMotionIdentification);
            return new ModifiedKeplerianElements(
                    objectsMotion,
                    WorldGeodeticSystem1984.GravitationalParameter).computePeriod();
        }),
        impulsiveManeuverSegment,
        TimeConstants.SecondsPerDay, // desired value, seconds
        0.01); // tolerance, seconds

periodConstraint.setName("Period_Constraint");

// add the constraint to the differential corrector
differentialCorrector.getConstraints().add(periodConstraint);

Segment Propagation Constraints

There are several constraints included in the Segment Propagation Library. The first is DelegateBasedConstraint, which allows you to specify a delegate that returns some computed value from the ephemeris of the segments propagated. For example, if you wanted to constrain an orbit such that its period is 24 hours, you could do so as shown below:

Java
DelegateBasedConstraint constraint = new DelegateBasedConstraint(
        DelegateBasedConstraintCallback.of(segmentResults -> {
            Motion1<Cartesian> objectsMotion = segmentResults.getStateForNextSegment().getMotion(satelliteID);
            return new ModifiedKeplerianElements(objectsMotion, WorldGeodeticSystem1984.GravitationalParameter).computePeriod();
        }),
        segment,
        TimeConstants.SecondsPerDay, // a desired value of 1 day in seconds
        0.01); // tolerance, seconds
solver.getConstraints().add(constraint);

The other constraints generally consist of evaluating a Scalar of some sort. For example, the ScalarAtEndOfSegmentConstraint will evaluate a Scalar at the end of a segment in the TargetedSegmentList. For example, if you want your segment final state to be at a specific altitude, you could specify that as shown below:

Java
ScalarAtEndOfSegmentConstraint constraint =
        new ScalarAtEndOfSegmentConstraint(segment,
                                           500000.0, // desired value of 500 km
                                           1.0); // 1 meter tolerance
ScalarCartographicElement altitudeScalar =
        new ScalarCartographicElement(earth,
                                      new ParameterizedOnStatePoint(constraint.getParameter(),
                                                                    earth.getFixedFrame(),
                                                                    SatelliteMotionIdentification),
                                      CartographicElement.HEIGHT);
constraint.setScalar(altitudeScalar);

solver.getConstraints().add(constraint);

Segment Propagation Variables

Many segments have values in them that can be used as variables. Although each segment has different properties that can be used as variables, the process used to create variables is similar across the different types of segments.

In the first example, where we showed solving a simple system of multivariable equations, the evaluate method of the function was created such that the variables were the familiar X and Y used in mathematics literature. In more complicated functions present within the Segment Propagation Library, the variables can be parameterized primitives or elements within the propagation state. These variables may be trivially altered, require auxiliary calculation, or even require that they are applied sequentially. To handle the more complex behavior involved in modeling changes to the propagation elements and segment primitives, we include a class that handles most of the heavy lifting- TargetedSegmentListFunction.

To illustrate how the complex behavior of variables could be a problem: consider an example where the initial state of propagation has two variables being solved for, the inclination and the right ascension of the ascending node. As we will explore shortly, the inclination may also need to change the right ascension (if the function solver tries to push the inclination to be less than 0 or more than 180 degrees, other orbital elements must flip otherwise the actual position of the satellite will be on the other side of the central body). This complex behavior requires the use of DelegateBasedVariable<T>. See the example under the appropriate section for the details.

Now let's dive a little further into the specifics of the underlying mechanics in SegmentPropagatorVariables. For the delegate type variable one typically modifies a value directly present in the state or performs a complex calculation which then modifies the state in some manner. If you look at the below example variable in the impulsive maneuver segment, you'll notice that when we call our variable setter we add the input value to the configuration's maneuver element. Each time the TargetedSegmentListFunction is evaluated in the solver it passes a clean configuration containing any initial state either passed-in or set when constructing the propagation sequence. Thus the passed in configuration functions as an initial value for the variable that the delegate alters. The solver itself maintains a record of the variable's current value with no knowledge of the initial value supplied by the configuration.

This behavior is similar for the other types of SegmentPropagatorVariables. For each of these types the initial value set on the definitional object itself (the variable) is used in propagation/function evaluation but is not present in iteration output. What this means is that when querying the values of the variables the iteration results will all be relative to the initial value on the definitional object or found within the clean configuration. If inspecting state elements directly the initial value will be incorporated and thus be the true value used in the computations and propagation.

Delegate Based Variable

DelegateBasedVariable<T> works with a setter delegate. The setter will take in the segments configuration and the actual value of the variable, and this delegate must modify the appropriate value in the configuration.

InitialStateSegment<T> can have any value in its state edited. If, for example, you don't know what inclination your satellite should start at, you can create such variables by calling the createVariable method on the segment:

Java
DelegateBasedVariable<InitialStateSegmentConfiguration> inclinationVariable = initialStateSegment.createVariable(
        0.2,  // maximum step, radians
        0.01, // perturbation, radians
        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(SatelliteMotionIdentification);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter);
            double rightAscensionCorrection = 0.0;
            double argumentOfPeriapsisCorrection = 0.0;
            double inclination = oldElements.getInclination() + variable;
            if (inclination > Math.PI) {
                // these checks are in here because it is possible that the differential corrector could
                // push the inclination to be greater than PI, and if it does,
                // we need to flip the inclination and right ascension.
                inclination -= Math.PI;
                rightAscensionCorrection = -Math.PI;
                argumentOfPeriapsisCorrection = -Math.PI;
            } else if (inclination < 0) {
                // This is correct. If the inclination becomes less than 0,
                // to avoid discontinuities in position and velocity,
                // the right ascension and AOP need to flip too.
                inclination += Math.PI;
                rightAscensionCorrection = Math.PI;
                argumentOfPeriapsisCorrection = Math.PI;
            }

            double newArgumentOfPeriapsis = Trig.zeroToTwoPi(oldElements.getArgumentOfPeriapsis() + argumentOfPeriapsisCorrection);
            double newRightAscension = Trig.zeroToTwoPi(oldElements.getRightAscensionOfAscendingNode() + rightAscensionCorrection);

            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(
                    oldElements.getRadiusOfPeriapsis(),
                    oldElements.getInverseSemimajorAxis(),
                    inclination,
                    newArgumentOfPeriapsis,
                    newRightAscension,
                    oldElements.getTrueAnomaly(),
                    oldElements.getGravitationalParameter());

            configuration.modifyMotion(SatelliteMotionIdentification, newElements.toCartesian());
        }));

When making the delegates for variables, you must remember that the Newton Raphson solver does not know anything about the physics of the problem. It doesn't know that the inclination can only be between 0 and 180 degrees. When it computes how much the variable must change in order to meet its constraints, it may try to push the inclination beyond its bounds. It is up to the variable to handle that properly, hence the extra checks in this setter delegate.

Notice how by modifying the inclination you may also need to modify other values too. This is perfectly acceptable as long as conceptually one variable is being modified.

However, if you want to solve for both inclination and for the right ascension of the ascending node, you would need a second variable. Do not be tempted to try to apply two separate variables in one delegate. Do not do this:

Java
// DO NOT DO THIS

DelegateBasedVariable<InitialStateSegmentConfiguration> incorrectVariable = initialStateSegment.createVariable(
        0.2,  // maximum step, radians
        0.01, // perturbation, radians
        SetVariableCallback.of((variable, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(SatelliteMotionIdentification);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter);
            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(
                    oldElements.getRadiusOfPeriapsis(),
                    oldElements.getInverseSemimajorAxis(),
                    oldElements.getInclination() + variable,  // this will not work properly
                    oldElements.getRightAscensionOfAscendingNode() + variable,  // this will not work properly
                    oldElements.getArgumentOfPeriapsis(),
                    oldElements.getTrueAnomaly(),
                    gravitationalParameter);
            configuration.modifyMotion(SatelliteMotionIdentification, newElements.toCartesian());
        }));

The function solver will try to solve it and it may even converge, but the results you get will not be what you expect.

Any of the values in an ImpulsiveManeuverInformation can be directly edited, except for the Orientation (get / set). Varying a delta-v is a common use case:

Java
DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable = firstBurn.createVariable(
        100.0, // maximum step, meters/second
        0.1,   // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));

Parameterized Variables

Throughout the library, there are many settings that you might want to vary in a function solver. For those, consider using ParameterizedScalarVariable, ParameterizedDoubleVariable, ParameterizedDurationVariable, or ParameterizedDateVariable. These variables all work the same way, by providing a parameterized placeholder for the value in question.

For example, any Scalar in any segment can be a variable by using the ParameterizedScalarVariable. This variable requires that you use its Value (get) property, in place of the Scalar, as shown below:

Java
// What drag coefficient will give my satellite an altitude of 250 km after 1 day of propagation?

// common items
EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
double gravityParameter = WorldGeodeticSystem1984.GravitationalParameter;

SolarGeophysicalData solarData = new ConstantSolarGeophysicalData(150.0, 3.0);

// make the segment list
TargetedSegmentList targetedSegment = new TargetedSegmentList();

// the propagate segment
NumericalPropagatorSegment propagationSegment = new NumericalPropagatorSegment();

// configure the numerical propagator
NumericalPropagatorDefinition propagatorDefinition = new NumericalPropagatorDefinition();
propagationSegment.setPropagatorDefinition(propagatorDefinition);
propagatorDefinition.setEpoch(TimeConstants.J2000);
RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
integrator.setRelativeTolerance(Constants.Epsilon13);
integrator.setMinimumStepSize(1.0);
integrator.setMaximumStepSize(86400.0);
integrator.setInitialStepSize(300.0);
integrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);
propagatorDefinition.setIntegrator(integrator);
propagationSegment.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(1.0)));

// create the numerical propagator point
PropagationNewtonianPoint satellitePoint = new PropagationNewtonianPoint();
propagatorDefinition.getIntegrationElements().add(satellitePoint);
satellitePoint.setIdentification(SatelliteMotionIdentification);
satellitePoint.setMass(Scalar.toScalar(500.0));

//Values of an orbit starting at 255 km
satellitePoint.setInitialPosition(new Cartesian(6633140.0, 0.0, 0.0));
satellitePoint.setInitialVelocity(new Cartesian(0.0, 6812.52, 3689.9));

// add gravity
satellitePoint.getAppliedForces().add(new TwoBodyGravity(satellitePoint.getIntegrationPoint(),
                                                         earth,
                                                         gravityParameter));

// make the differential corrector
TargetedSegmentListDifferentialCorrector differentialCorrector =
    new TargetedSegmentListDifferentialCorrector();
differentialCorrector.setName("Differential_Corrector");
differentialCorrector.getSolver().setMultithreaded(false);
differentialCorrector.setMaximumIterations(25);

// the variable for the differential corrector
ParameterizedScalarVariable dragCoefficientVariable =
        new ParameterizedScalarVariable(0.1, // initial value of the drag coefficient
                                        0.3, // maximum step
                                        0.001, // perturbation
                                        propagationSegment);
dragCoefficientVariable.setName("Drag Coefficient");
dragCoefficientVariable.getSettings().setVariableTolerance(0.0001);

// the constraint
ScalarAtEndOfSegmentConstraint altitudeConstraint =
        new ScalarAtEndOfSegmentConstraint(propagationSegment,
                                           250000.0, // desired value, meters
                                           0.1); // tolerance, meters

ScalarCartographicElement altitudeScalar =
        new ScalarCartographicElement(earth,
                                      new ParameterizedOnStatePoint(altitudeConstraint.getParameter(),
                                                                    satellitePoint.getIntegrationFrame(),
                                                                    satellitePoint.getIdentification()),
                                      CartographicElement.HEIGHT);
altitudeConstraint.setScalar(altitudeScalar);

// add the variable and the constraint
differentialCorrector.getVariables().add(dragCoefficientVariable);
differentialCorrector.getConstraints().add(altitudeConstraint);

// add the differential corrector
targetedSegment.getOperators().add(differentialCorrector);

// create the drag force
ScalarAtmosphericDensity density =
        new ScalarDensityMsis90(satellitePoint.getIntegrationPoint(), solarData);
AtmosphericDragForce dragForce =
        new AtmosphericDragForce(density,
                                 dragCoefficientVariable.getValue(), // here we use the parameterized scalar
                                 Scalar.toScalar(50.0)); // surface area, meters squared

// add the force
satellitePoint.getAppliedForces().add(dragForce);

// add the segment
targetedSegment.getSegments().add(propagationSegment);

// propagate
SegmentListPropagator propagator = targetedSegment.getSegmentListPropagator();

// since we know we were propagating a TargetedSegmentList, we can cast the results safely
TargetedSegmentListResults results = (TargetedSegmentListResults)propagator.propagateSegmentList();

TargetedSegmentListDifferentialCorrectorResults correctorResults =
    (TargetedSegmentListDifferentialCorrectorResults)results.getOperatorResults().get(0);

SolvableMultivariableFunctionResults finalFunctionResults =
    correctorResults.getFunctionSolverResults().getFinalIteration().getFunctionResult();

double finalAltitude = finalFunctionResults.getConstraintValues()[0];
double dragCoefficient = finalFunctionResults.getVariablesUsed()[0];
// final altitude will be 250000.0 meters + or - 0.1 meters
// the drag coefficient will be 0.03803

In fact, the DesiredValue (get / set) of a constraint in one differential corrector can be varied by an outer differential corrector using a ParameterizedDoubleVariable, as shown below:

Java
TargetedSegmentList outerTargetedSegmentList = new TargetedSegmentList();
outerTargetedSegmentList.setName("Outer_Target_List");

TargetedSegmentListDifferentialCorrector solverFor24HourPeriod = new TargetedSegmentListDifferentialCorrector();
solverFor24HourPeriod.setName("Solving_For_24_Hour_Period");

// 24 hour period constraint
ScalarModifiedKeplerianElement periodScalar =
        new ScalarModifiedKeplerianElement(
                gravitationalParameter,
                integrationPoint.getIntegrationPoint(),
                KeplerianElement.PERIOD,
                frame);
ScalarAtEndOfNumericalSegmentConstraint periodConstraint =
        new ScalarAtEndOfNumericalSegmentConstraint(
                periodScalar,
                dv2,
                TimeConstants.SecondsPerDay, // desired value, seconds
                0.0001);// tolerance, seconds
periodConstraint.setName("Orbital_Period");

// vary the desired value of the nested Radius of Apoapsis
ParameterizedDoubleVariable nestedDesiredRadiusOfApoapsisVariable = new ParameterizedDoubleVariable(
        42000000.0, // initial value, meters
        5000000.0,  // maximum step, meters
        500000.0,   // perturbation, meters
        transferOrbit);
nestedDesiredRadiusOfApoapsisVariable.setName("Desired_Radius_Of_Apoapsis");
radiusOfApoapsisPostDv1Constraint.setDesiredValue(nestedDesiredRadiusOfApoapsisVariable.getValue());

// add the variable and constraint
solverFor24HourPeriod.getConstraints().add(periodConstraint);
solverFor24HourPeriod.getVariables().add(nestedDesiredRadiusOfApoapsisVariable);

// add the differentialCorrector to the segment
outerTargetedSegmentList.getOperators().add(solverFor24HourPeriod);

Below is a complete example using a nested TargetedSegmentList, demonstrating many of the features described above:

Java
final String motionID = "Satellite";
final String fuelMassName = "Fuel_Mass";
final String dryMassName = "Dry_Mass";

EarthCentralBody earth = CentralBodiesFacet.getFromContext().getEarth();
ReferenceFrame frame = earth.getInertialFrame();
double gravitationalParameter = WorldGeodeticSystem1984.GravitationalParameter;
JulianDate epoch = TimeConstants.J2000.toTimeStandard(TimeStandard.getCoordinatedUniversalTime());
Motion1<Cartesian> initialConditions = new Motion1<Cartesian>(
    new Cartesian(8000000.0, 0.0, 0.0),
    new Cartesian(1500.0, 7500.0, 0.0));

NumericalPropagatorDefinition numericalPropagatorDefinition = new NumericalPropagatorDefinition();
numericalPropagatorDefinition.setEpoch(epoch);

PropagationNewtonianPoint integrationPoint = new PropagationNewtonianPoint(
    motionID,
    frame,
    initialConditions.getValue(),
    initialConditions.getFirstDerivative());

// simple gravity
TwoBodyGravity gravity = new TwoBodyGravity(
    integrationPoint.getIntegrationPoint(),
    earth,
    gravitationalParameter);
integrationPoint.getAppliedForces().add(gravity);

// fuel mass
PropagationScalar fuelMass = new PropagationScalar(900.0);
fuelMass.setIdentification(fuelMassName);
fuelMass.setScalarDerivative(new ScalarFixed(0.0));

// dry mass
AuxiliaryStateScalar dryMass = new AuxiliaryStateScalar(Scalar.toScalar(100.0));
dryMass.setIdentification(dryMassName);
integrationPoint.setMass(Scalar.add(fuelMass.getIntegrationValue(), dryMass.getAuxiliaryScalar()));

// configure the propagator
RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator();
integrator.setRelativeTolerance(Constants.Epsilon13);
integrator.setMinimumStepSize(1.0);
integrator.setMaximumStepSize(86400.0);
integrator.setInitialStepSize(300.0);
integrator.setStepSizeBehavior(KindOfStepSize.RELATIVE);

numericalPropagatorDefinition.setIntegrator(integrator);
numericalPropagatorDefinition.getIntegrationElements().add(integrationPoint);
numericalPropagatorDefinition.getIntegrationElements().add(fuelMass);
numericalPropagatorDefinition.getAuxiliaryElements().add(dryMass);

NumericalInitialStateSegment initial = new NumericalInitialStateSegment();
initial.setName("Initial_State_Segment");
initial.setPropagatorDefinition(numericalPropagatorDefinition);

NumericalPropagatorSegment propagateToPerigeeSegment = new NumericalPropagatorSegment();
propagateToPerigeeSegment.setName("Propagate_To_Perigee");
propagateToPerigeeSegment.setPropagatorDefinition(numericalPropagatorDefinition);

// perigee stopping condition

ScalarStoppingCondition perStop = new ScalarStoppingCondition(
        new ScalarModifiedKeplerianElement(
                gravitationalParameter,
                integrationPoint.getIntegrationPoint(),
                KeplerianElement.TRUE_ANOMALY,
                frame),
        0.0,          // threshold in radians
        0.0000000001, // value tolerance, radians
        StopType.ANY_THRESHOLD);
perStop.setAngularSetting(CircularRange.NEGATIVE_PI_TO_PI); // to avoid the branch cut
perStop.setName("Periapsis");

propagateToPerigeeSegment.getStoppingConditions().add(perStop);

TargetedSegmentList innerTargetSegment = new TargetedSegmentList();
innerTargetSegment.setName("Inner_Target_List");

ImpulsiveManeuverSegment dv1 = new ImpulsiveManeuverSegment();
dv1.setName("First_Maneuver");

// maneuver info
ImpulsiveManeuverInformation dv1Info = new ImpulsiveManeuverInformation(
        motionID,
        new Cartesian(200.0, 0.0, 0.0), // first guess
        fuelMassName,
        dryMassName,
        Scalar.toScalar(4500.0), // exhaust velocity, meters/second
        InvalidFuelStateBehavior.THROW_EXCEPTION);

dv1Info.setOrientation(new AxesVelocityOrbitNormal(
        dv1Info.getPropagationPoint(),
        earth));

dv1.setManeuver(dv1Info);

NumericalPropagatorSegment transferOrbit = new NumericalPropagatorSegment();
transferOrbit.setName("Transfer_Orbit");
transferOrbit.setPropagatorDefinition(numericalPropagatorDefinition);

// apogee stopping condition
ScalarStoppingCondition apoapsisStoppingCondition = new ScalarStoppingCondition(
        new ScalarModifiedKeplerianElement(
                gravitationalParameter,
                integrationPoint.getIntegrationPoint(),
                KeplerianElement.TRUE_ANOMALY,
                frame),
        Math.PI,    // threshold, radians
        0.00000001, // value tolerance, radians
        StopType.ANY_THRESHOLD);
apoapsisStoppingCondition.setAngularSetting(CircularRange.ZERO_TO_TWO_PI); // let the system know to avoid the branch cut
apoapsisStoppingCondition.setName("Apoapsis_Stopping_Condition");
transferOrbit.getStoppingConditions().add(apoapsisStoppingCondition);

ImpulsiveManeuverSegment dv2 = new ImpulsiveManeuverSegment();
dv2.setName("Second_Maneuver");
ImpulsiveManeuverInformation dv2Info = new ImpulsiveManeuverInformation(
        motionID,
        new Cartesian(500.0, 0.0, 0.0),
        fuelMassName,
        dryMassName,
        Scalar.toScalar(4500.0),
        InvalidFuelStateBehavior.DO_FULL_DELTA_V);
dv2.setManeuver(dv2Info);

dv2Info.setOrientation(new AxesVelocityOrbitNormal(
        dv2Info.getPropagationPoint(),
        earth));

innerTargetSegment.getSegments().add(dv1);
innerTargetSegment.getSegments().add(transferOrbit);
innerTargetSegment.getSegments().add(dv2);

TargetedSegmentListDifferentialCorrector hohmannTransferDifferentialCorrector =
    new TargetedSegmentListDifferentialCorrector();
hohmannTransferDifferentialCorrector.setName("Solve_For_Hohmann_Transfer");

// variable 1, dv1's x velocity
DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> dv1VelocityXVariable = dv1.createVariable(
        200.0, // maximum step, meters/second
        0.1,   // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));
dv1VelocityXVariable.setName("ImpulsiveManeuver1_ThrustVector_X");

// constraint 1, dv1's
ScalarModifiedKeplerianElement apoapsisScalar =
        new ScalarModifiedKeplerianElement(
                gravitationalParameter,
                integrationPoint.getIntegrationPoint(),
                KeplerianElement.RADIUS_OF_APOAPSIS,
                frame);
ScalarAtEndOfNumericalSegmentConstraint radiusOfApoapsisPostDv1Constraint =
        new ScalarAtEndOfNumericalSegmentConstraint(
                apoapsisScalar,
                transferOrbit,
                0.0, // desired value, meters, leaving as null since a outer variable will set this
                0.0001); // tolerance, meters
radiusOfApoapsisPostDv1Constraint.setName("Radius_Of_Apoapsis");

// variable 2, dv2's x velocity
DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> dv2VelocityXVariable = dv2.createVariable(
        200.0, // maximum step, meters/second
        0.1,   // perturbation, meters/second
        SetVariableCallback.of((variable, configuration) -> {
            configuration.getManeuver().setX(configuration.getManeuver().getX() + variable);
        }));
dv2VelocityXVariable.setName("ImpulsiveManeuver2_ThrustVector_X");

// constraint 2, eccentricity constraint after the second maneuver
ScalarModifiedKeplerianElement eccentricityScalar =
        new ScalarModifiedKeplerianElement(
                gravitationalParameter,
                integrationPoint.getIntegrationPoint(),
                KeplerianElement.ECCENTRICITY,
                frame);
ScalarAtEndOfNumericalSegmentConstraint eccentricityPostDv2Constraint =
        new ScalarAtEndOfNumericalSegmentConstraint(
                eccentricityScalar,
                dv2,
                0.1, // desired value, unitless
                0.000001); // tolerance, unitless
eccentricityPostDv2Constraint.setName("Eccentricity");

// add the variables
hohmannTransferDifferentialCorrector.getVariables().add(dv1VelocityXVariable);
hohmannTransferDifferentialCorrector.getVariables().add(dv2VelocityXVariable);

// add the constraints
hohmannTransferDifferentialCorrector.getConstraints().add(radiusOfApoapsisPostDv1Constraint);
hohmannTransferDifferentialCorrector.getConstraints().add(eccentricityPostDv2Constraint);

// add the corrector to the targeted segment
innerTargetSegment.getOperators().add(hohmannTransferDifferentialCorrector);


TargetedSegmentList outerTargetedSegmentList = new TargetedSegmentList();
outerTargetedSegmentList.setName("Outer_Target_List");

TargetedSegmentListDifferentialCorrector solverFor24HourPeriod = new TargetedSegmentListDifferentialCorrector();
solverFor24HourPeriod.setName("Solving_For_24_Hour_Period");

// 24 hour period constraint
ScalarModifiedKeplerianElement periodScalar =
        new ScalarModifiedKeplerianElement(
                gravitationalParameter,
                integrationPoint.getIntegrationPoint(),
                KeplerianElement.PERIOD,
                frame);
ScalarAtEndOfNumericalSegmentConstraint periodConstraint =
        new ScalarAtEndOfNumericalSegmentConstraint(
                periodScalar,
                dv2,
                TimeConstants.SecondsPerDay, // desired value, seconds
                0.0001);// tolerance, seconds
periodConstraint.setName("Orbital_Period");

// vary the desired value of the nested Radius of Apoapsis
ParameterizedDoubleVariable nestedDesiredRadiusOfApoapsisVariable = new ParameterizedDoubleVariable(
        42000000.0, // initial value, meters
        5000000.0,  // maximum step, meters
        500000.0,   // perturbation, meters
        transferOrbit);
nestedDesiredRadiusOfApoapsisVariable.setName("Desired_Radius_Of_Apoapsis");
radiusOfApoapsisPostDv1Constraint.setDesiredValue(nestedDesiredRadiusOfApoapsisVariable.getValue());

// add the variable and constraint
solverFor24HourPeriod.getConstraints().add(periodConstraint);
solverFor24HourPeriod.getVariables().add(nestedDesiredRadiusOfApoapsisVariable);

// add the differentialCorrector to the segment
outerTargetedSegmentList.getOperators().add(solverFor24HourPeriod);


NumericalPropagatorSegment finalPropagator = new NumericalPropagatorSegment();
finalPropagator.setName("Final_Orbit");
finalPropagator.setPropagatorDefinition(numericalPropagatorDefinition);
finalPropagator.getStoppingConditions().add(new DurationStoppingCondition(Duration.fromDays(1.0)));
finalPropagator.getStoppingConditions().get(0).setName("Duration");

SegmentList overallList = new SegmentList();
overallList.setName("Overall_List");
overallList.getSegments().add(initial);
overallList.getSegments().add(propagateToPerigeeSegment);

// add the inner targeted segment list to the outer targeted segment list
outerTargetedSegmentList.getSegments().add(innerTargetSegment);

// and add the outer list
overallList.getSegments().add(outerTargetedSegmentList);
overallList.getSegments().add(finalPropagator);

SegmentPropagator propagator = overallList.getSegmentPropagator(new EvaluatorGroup());
SegmentListResults propagatorResults = (SegmentListResults) propagator.propagate();