Click or drag to resize

Multivariable Function Solvers

Function solvers allow sets of multivariable equations to be driven to some set of values. STK Components includes the NewtonRaphsonMultivariableFunctionSolver type. The solver will evaluate the function at the current value, numerically compute the derivative of the function, and assume a linear function to guess the solution. Using that predicted solution, it will start again until the difference between the computed value and the constraint's desired values are equal within tolerance. Because the function being evaluated may not be linear, the variables have a MaximumStep (get / set) property limiting them from stepping too far.

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:

x^2 + y^2 = 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 2 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, // 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 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. As a result, if you change your initial guess, it will find different answers:

Java
// re-create 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);

// re-create 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) -> {
            ImpulsiveManeuverInformation info = configuration.get(maneuverInfo);
            info.setX(info.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 the 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
PropagationNewtonianPoint propagationPoint = createPropagatedNewtonianPoint(SatelliteMotionIdentification);

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, creating the variables themselves is similar across the different types of segment.

In the above example, the evaluate method of the function was created such that the variables were simply the X and Y used in the function. However, that can cause problems with function solvers in more complicated functions where the variables are coupled or applied sequentially. This is the case with the TargetedSegmentListFunction. The way that the function deals with it is that instead of solving for the variable directly, the change of the variable is solved for.

To illustrate how the sequential application 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). So take the case where the inclination variable is applied first and it happens to change the right ascension as well as the inclination. If the right ascension variable then would set itself to what the function solver said it should be, it will clobber part of the change that the inclination made. This will confuse the function solver because the intended values for the inclination and right ascension were not actually used in the function evaluation and it will often not converge or converge improperly.

However, that clobbering of previous changes is not an issue if we have the function solver variables actually be the delta of the relevant value in the function. The TargetedSegmentListFunction is designed so that the variable values it receives from the function solver are actually the total change of relevant value. Since the values in the state that is getting edited by the variables are now incremented instead of set absolutely, later variables can't undo the changes made by previous variables.

So the way that the inclination and right ascension example actually works: in the case where the function solver computes such a delta to the inclination that it will need to flip some of the other orbital elements. The inclination variable is applied by adding its delta to the old inclination, and it flips the other orbital elements. Then when the right ascension variable is applied, it has the new state that the inclination variable created, and it adds its delta to the value of the right ascension that the inclination computed.

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((variableDelta, 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() + variableDelta;
            if (inclination > Math.PI) {
                // this 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 RightAscension.
                inclination = 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 = 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((variableDelta, configuration) -> {
            Motion1<Cartesian> currentMotion = configuration.getMotion(SatelliteMotionIdentification);
            ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter);
            ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(
                    oldElements.getRadiusOfPeriapsis(), 
                    oldElements.getInverseSemimajorAxis(),
                    oldElements.getInclination() + variableDelta,  // this will not work properly
                    oldElements.getRightAscensionOfAscendingNode() + variableDelta,  // 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((variableDelta, configuration) -> {
            ImpulsiveManeuverInformation info = configuration.get(firstBurnInfo);
            info.setX(info.getX() + variableDelta);
        }));

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.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
ScalarAtEndOfSegmentConstraint periodConstraint = new ScalarAtEndOfSegmentConstraint(
        dv2, 
        TimeConstants.SecondsPerDay, // desired value, seconds
        0.0001);                     // tolerance, seconds
ScalarModifiedKeplerianElement periodScalar = new ScalarModifiedKeplerianElement(
        gravitationalParameter, 
        new ParameterizedOnStatePoint(
                periodConstraint.getParameter(), 
                frame, 
                motionID),
        KeplerianElement.PERIOD, 
        frame);
periodConstraint.setScalar(periodScalar);
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 full example using a nested TargetedSegmentList:

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

NumericalPropagatorHelperMethods.jplData().useForCentralBodyPositions(CentralBodiesFacet.getFromContext());
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
final 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.getManeuvers().add(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);



final 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.getManeuvers().add(dv2Info);

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


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


final 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((variableDelta, configuration) -> {
            ImpulsiveManeuverInformation info = configuration.get(dv1Info);
            info.setX(info.getX() + variableDelta);
        }));
dv1VelocityXVariable.setName("ImpulsiveManeuver_ThrustVector_X");

// constraint 1, dv1's

ScalarAtEndOfSegmentConstraint radiusOfApoapsisPostDv1Constraint = new ScalarAtEndOfSegmentConstraint(
        transferOrbit, 
        0, // desired value, meters, leaving as null since a outer variable will set this
        0.0001); // tolerance, meters


ScalarModifiedKeplerianElement apoapsisScalar = new ScalarModifiedKeplerianElement(
        gravitationalParameter, 
        new ParameterizedOnStatePoint(
                radiusOfApoapsisPostDv1Constraint.getParameter(),
                frame, 
                motionID),
        KeplerianElement.RADIUS_OF_APOAPSIS, 
        frame);

radiusOfApoapsisPostDv1Constraint.setScalar(apoapsisScalar);
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((variableDelta, configuration) -> {
            ImpulsiveManeuverInformation info = configuration.get(dv2.getManeuvers().get(0));
            info.setX(info.getX() + variableDelta);
        }));
dv2VelocityXVariable.setName("ImpulsiveManeuver_ThrustVector_X");

// constraint 2, eccentricity constraint after the second maneuver
ScalarAtEndOfSegmentConstraint eccentricityPostDv2Constraint = new ScalarAtEndOfSegmentConstraint(
        dv2, 
        0.1,       // desired value, unitless
        0.000001); // tolerance, unitless

ScalarModifiedKeplerianElement eccentricityScalar = new ScalarModifiedKeplerianElement(
        gravitationalParameter, 
        new ParameterizedOnStatePoint(
                eccentricityPostDv2Constraint.getParameter(),
                frame, 
                motionID), 
        KeplerianElement.ECCENTRICITY, 
        frame);
eccentricityPostDv2Constraint.setScalar(eccentricityScalar);
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
ScalarAtEndOfSegmentConstraint periodConstraint = new ScalarAtEndOfSegmentConstraint(
        dv2, 
        TimeConstants.SecondsPerDay, // desired value, seconds
        0.0001);                     // tolerance, seconds
ScalarModifiedKeplerianElement periodScalar = new ScalarModifiedKeplerianElement(
        gravitationalParameter, 
        new ParameterizedOnStatePoint(
                periodConstraint.getParameter(), 
                frame, 
                motionID),
        KeplerianElement.PERIOD, 
        frame);
periodConstraint.setScalar(periodScalar);
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();