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, VariableTolerance, and, Scaling. For more information on the available settings see the settings class.
 Note |
|---|
The functionality described in this topic requires a license for the Segment Propagation Library. |
Newton Raphson Function SolverSuppose 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:
SolverVariableSettings xVariable = new SolverVariableSettings(0.3); // maximum step // initial value is left at the default 0 xVariable.Name = "x"; SolverVariableSettings yVariable = new SolverVariableSettings(0.3); yVariable.Name = "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.
SolverConstraintSettings firstEquation = new SolverConstraintSettings(25, // desired value 0.001); // tolerance firstEquation.Name = "firstEquation"; SolverConstraintSettings secondEquation = new SolverConstraintSettings(0, 0.001); secondEquation.Name = "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.
NewtonRaphsonMultivariableFunctionSolver solver = new NewtonRaphsonMultivariableFunctionSolver();
solver.Constraints.Add(firstEquation);
solver.Constraints.Add(secondEquation);
solver.Variables.Add(xVariable);
solver.Variables.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:
public override 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 void SetPerturbationValues(double xPerturbation, double yPerturbation) { PerturbationValues = new[] { xPerturbation, yPerturbation }; }
The following code shows how to create and run the NewtonRaphsonMultivariableFunctionSolver:
// SimpleFunction is our example class derived from SolvableMultivariableFunction SimpleFunction function = new SimpleFunction(); function.SetPerturbationValues(0.1, 0.1); solver.Function = function; solver.FindSolution(25); var results = solver.LastRunsResults; double[] variablesUsed = results.FinalIteration.FunctionResult.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:
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:
// 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.Variables.Add(xVariable); solver.Variables.Add(yVariable); solver.Constraints.Add(firstEquation); solver.Constraints.Add(secondEquation); solver.Function = function; solver.FindSolution(25); results = solver.LastRunsResults; variablesUsed = results.FinalIteration.FunctionResult.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 property on the variables prevents those large jumps.
Segmented Propagation and Function SolversThe 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 CorrectorTargetedSegmentListDifferentialCorrector requires that at least one SegmentPropagatorConstraint and one SegmentPropagatorVariable are configured on it, as shown below:
// create the differential corrector TargetedSegmentListDifferentialCorrector differentialCorrector = new TargetedSegmentListDifferentialCorrector(); differentialCorrector.Name = "Corrector"; // add the differential corrector to the TargetedSegmentList targetedSegmentList.Operators.Add(differentialCorrector); // create a variable DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable = impulsiveManeuverSegment.CreateVariable(200.0, // maximum step, meters/second 1.0, // perturbation, meters/second (variable, configuration) => { configuration.Maneuver.X += variable; }); velocityXVariable.Name = "Velocity_X_Variable"; // add the variable to the differential corrector differentialCorrector.Variables.Add(velocityXVariable); // create a constraint DelegateBasedConstraint periodConstraint = new DelegateBasedConstraint( segmentResults => { Motion<Cartesian> objectsMotion = segmentResults.StateForNextSegment.GetMotion<Cartesian>(SatelliteMotionIdentification); return new ModifiedKeplerianElements(objectsMotion, WorldGeodeticSystem1984.GravitationalParameter).ComputePeriod(); }, impulsiveManeuverSegment, TimeConstants.SecondsPerDay, // desired value, seconds 0.01); // tolerance, seconds periodConstraint.Name = "Period_Constraint"; // add the constraint to the differential corrector differentialCorrector.Constraints.Add(periodConstraint);
Segment Propagation ConstraintsThere 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:
var constraint = new DelegateBasedConstraint( segmentResults => { Motion<Cartesian> objectsMotion = segmentResults.StateForNextSegment.GetMotion<Cartesian>(satelliteID); return new ModifiedKeplerianElements(objectsMotion, WorldGeodeticSystem1984.GravitationalParameter).ComputePeriod(); }, segment, TimeConstants.SecondsPerDay, // a desired value of 1 day in seconds 0.01); // tolerance, seconds solver.Constraints.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:
var constraint = new ScalarAtEndOfSegmentConstraint(segment, 500000.0, // desired value of 500 km 1.0); // 1 meter tolerance var altitudeScalar = new ScalarCartographicElement(earth, new ParameterizedOnStatePoint(constraint.Parameter, earth.FixedFrame, SatelliteMotionIdentification), CartographicElement.Height); constraint.Scalar = altitudeScalar; solver.Constraints.Add(constraint);
Segment Propagation VariablesMany 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 VariableDelegateBasedVariable<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:
DelegateBasedVariable<InitialStateSegmentConfiguration> inclinationVariable = initialStateSegment.CreateVariable(
0.2, // maximum step, radians
0.01, // perturbation, radians
(variable, configuration) =>
{
Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(SatelliteMotionIdentification);
ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter);
double rightAscensionCorrection = 0;
double argumentOfPeriapsisCorrection = 0;
double inclination = oldElements.Inclination + 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.ArgumentOfPeriapsis + argumentOfPeriapsisCorrection);
double newRightAscension = Trig.ZeroToTwoPi(oldElements.RightAscensionOfAscendingNode + rightAscensionCorrection);
ModifiedKeplerianElements newElements =
new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis,
oldElements.InverseSemimajorAxis,
inclination,
newArgumentOfPeriapsis,
newRightAscension,
oldElements.TrueAnomaly,
oldElements.GravitationalParameter);
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:
// DO NOT DO THIS DelegateBasedVariable<InitialStateSegmentConfiguration> incorrectVariable = initialStateSegment.CreateVariable( 0.2, // maximum step, radians 0.01, // perturbation, radians (variable, configuration) => { Motion<Cartesian> currentMotion = configuration.GetMotion<Cartesian>(SatelliteMotionIdentification); ModifiedKeplerianElements oldElements = new ModifiedKeplerianElements(currentMotion, gravitationalParameter); ModifiedKeplerianElements newElements = new ModifiedKeplerianElements(oldElements.RadiusOfPeriapsis, oldElements.InverseSemimajorAxis, oldElements.Inclination + variable, // This will not work properly oldElements.RightAscensionOfAscendingNode + variable, // This will not work properly oldElements.ArgumentOfPeriapsis, oldElements.TrueAnomaly, 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. Varying a delta-v is a common use case:
DelegateBasedVariable<ImpulsiveManeuverSegmentConfiguration> velocityXVariable =
firstBurn.CreateVariable(100, // maximum step, meters/second
0.1, // perturbation, meters/second
(variable, configuration) => { configuration.Maneuver.X += variable; });Parameterized VariablesThroughout 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 property, in place of the Scalar, as shown below:
// What drag coefficient will give my satellite an altitude of 250 km after 1 day of propagation? // common items EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth; 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.PropagatorDefinition = propagatorDefinition; propagatorDefinition.Epoch = TimeConstants.J2000; RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator(); integrator.RelativeTolerance = Constants.Epsilon13; integrator.MinimumStepSize = 1.0; integrator.MaximumStepSize = 86400.0; integrator.InitialStepSize = 300.0; integrator.StepSizeBehavior = KindOfStepSize.Relative; propagatorDefinition.Integrator = integrator; propagationSegment.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromDays(1.0))); // create the numerical propagator point PropagationNewtonianPoint satellitePoint = new PropagationNewtonianPoint(); propagatorDefinition.IntegrationElements.Add(satellitePoint); satellitePoint.Identification = SatelliteMotionIdentification; satellitePoint.Mass = 500.0; //Values of an orbit starting at 255 km satellitePoint.InitialPosition = new Cartesian(6633140.0, 0.0, 0.0); satellitePoint.InitialVelocity = new Cartesian(0.0, 6812.52, 3689.9); // add gravity satellitePoint.AppliedForces.Add(new TwoBodyGravity(satellitePoint.IntegrationPoint, earth, gravityParameter)); // make the differential corrector TargetedSegmentListDifferentialCorrector differentialCorrector = new TargetedSegmentListDifferentialCorrector(); differentialCorrector.Name = "Differential_Corrector"; differentialCorrector.Solver.Multithreaded = false; differentialCorrector.MaximumIterations = 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.Name = "Drag Coefficient"; dragCoefficientVariable.Settings.VariableTolerance = 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.Parameter, satellitePoint.IntegrationFrame, satellitePoint.Identification), CartographicElement.Height); altitudeConstraint.Scalar = altitudeScalar; // add the variable and the constraint differentialCorrector.Variables.Add(dragCoefficientVariable); differentialCorrector.Constraints.Add(altitudeConstraint); // add the differential corrector targetedSegment.Operators.Add(differentialCorrector); // create the drag force ScalarAtmosphericDensity density = new ScalarDensityMsis90(satellitePoint.IntegrationPoint, solarData); AtmosphericDragForce dragForce = new AtmosphericDragForce(density, dragCoefficientVariable.Value, // here we use the parameterized scalar 50.0); // surface area, meters squared // add the force satellitePoint.AppliedForces.Add(dragForce); // add the segment targetedSegment.Segments.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.OperatorResults[0]; SolvableMultivariableFunctionResults finalFunctionResults = correctorResults.FunctionSolverResults.FinalIteration.FunctionResult; 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 of a constraint in one differential corrector can be varied by an outer differential corrector using a ParameterizedDoubleVariable, as shown below:
TargetedSegmentList outerTargetedSegmentList = new TargetedSegmentList(); outerTargetedSegmentList.Name = "Outer_Target_List"; TargetedSegmentListDifferentialCorrector solverFor24HourPeriod = new TargetedSegmentListDifferentialCorrector(); solverFor24HourPeriod.Name = "Solving_For_24_Hour_Period"; // 24 hour period constraint ScalarModifiedKeplerianElement periodScalar = new ScalarModifiedKeplerianElement(gravitationalParameter, integrationPoint.IntegrationPoint, KeplerianElement.Period, frame); ScalarAtEndOfNumericalSegmentConstraint periodConstraint = new ScalarAtEndOfNumericalSegmentConstraint(periodScalar, dv2, TimeConstants.SecondsPerDay, // desired value, seconds 0.0001); // tolerance, seconds periodConstraint.Name = "Orbital_Period"; // vary the desired value of the nested Radius of Apoapsis ParameterizedDoubleVariable nestedDesiredRadiusOfApoapsisVariable = new ParameterizedDoubleVariable(42000000, // initial value, meters 5000000.0, // maximum step, meters 500000.0, // perturbation, meters transferOrbit); nestedDesiredRadiusOfApoapsisVariable.Name = "Desired_Radius_Of_Apoapsis"; radiusOfApoapsisPostDv1Constraint.DesiredValue = nestedDesiredRadiusOfApoapsisVariable.Value; // add the variable and constraint solverFor24HourPeriod.Constraints.Add(periodConstraint); solverFor24HourPeriod.Variables.Add(nestedDesiredRadiusOfApoapsisVariable); // add the differentialCorrector to the segment outerTargetedSegmentList.Operators.Add(solverFor24HourPeriod);
Below is a complete example using a nested TargetedSegmentList, demonstrating many of the features described above:
const string motionID = "Satellite"; const string fuelMassName = "Fuel_Mass"; const string dryMassName = "Dry_Mass"; EarthCentralBody earth = CentralBodiesFacet.GetFromContext().Earth; ReferenceFrame frame = earth.InertialFrame; double gravitationalParameter = WorldGeodeticSystem1984.GravitationalParameter; JulianDate epoch = TimeConstants.J2000.ToTimeStandard(TimeStandard.CoordinatedUniversalTime); Motion<Cartesian> initialConditions = new Motion<Cartesian>(new Cartesian(8000000.0, 0.0, 0.0), new Cartesian(1500.0, 7500.0, 0.0)); NumericalPropagatorDefinition numericalPropagatorDefinition = new NumericalPropagatorDefinition(); numericalPropagatorDefinition.Epoch = epoch; PropagationNewtonianPoint integrationPoint = new PropagationNewtonianPoint(motionID, frame, initialConditions.Value, initialConditions.FirstDerivative); // simple gravity TwoBodyGravity gravity = new TwoBodyGravity(integrationPoint.IntegrationPoint, earth, gravitationalParameter); integrationPoint.AppliedForces.Add(gravity); // fuel mass PropagationScalar fuelMass = new PropagationScalar(900.0); fuelMass.Identification = fuelMassName; fuelMass.ScalarDerivative = new ScalarFixed(0.0); // dry mass AuxiliaryStateScalar dryMass = new AuxiliaryStateScalar(100.0); dryMass.Identification = dryMassName; // total mass integrationPoint.Mass = fuelMass.IntegrationValue + dryMass.AuxiliaryScalar; // configure the propagator RungeKuttaFehlberg78Integrator integrator = new RungeKuttaFehlberg78Integrator(); integrator.RelativeTolerance = Constants.Epsilon13; integrator.MinimumStepSize = 1.0; integrator.MaximumStepSize = 86400.0; integrator.InitialStepSize = 300.0; integrator.StepSizeBehavior = KindOfStepSize.Relative; numericalPropagatorDefinition.Integrator = integrator; numericalPropagatorDefinition.IntegrationElements.Add(integrationPoint); numericalPropagatorDefinition.IntegrationElements.Add(fuelMass); numericalPropagatorDefinition.AuxiliaryElements.Add(dryMass); NumericalInitialStateSegment initial = new NumericalInitialStateSegment(); initial.Name = "Initial_State_Segment"; initial.PropagatorDefinition = numericalPropagatorDefinition; NumericalPropagatorSegment propagateToPerigeeSegment = new NumericalPropagatorSegment(); propagateToPerigeeSegment.Name = "Propagate_To_Perigee"; propagateToPerigeeSegment.PropagatorDefinition = numericalPropagatorDefinition; // perigee stopping condition ScalarStoppingCondition perStop = new ScalarStoppingCondition(new ScalarModifiedKeplerianElement(gravitationalParameter, integrationPoint.IntegrationPoint, KeplerianElement.TrueAnomaly, frame), 0.0, // threshold in radians 0.0000000001, // value tolerance, radians StopType.AnyThreshold); perStop.AngularSetting = CircularRange.NegativePiToPi; // to avoid the branch cut perStop.Name = "Periapsis"; propagateToPerigeeSegment.StoppingConditions.Add(perStop); TargetedSegmentList innerTargetSegment = new TargetedSegmentList(); innerTargetSegment.Name = "Inner_Target_List"; ImpulsiveManeuverSegment dv1 = new ImpulsiveManeuverSegment(); dv1.Name = "First_Maneuver"; // maneuver info ImpulsiveManeuverInformation dv1Info = new ImpulsiveManeuverInformation(motionID, new Cartesian(200.0, 0.0, 0.0), // first guess fuelMassName, dryMassName, 4500.0, // exhaust velocity, meters/second InvalidFuelStateBehavior.ThrowException); dv1Info.Orientation = new AxesVelocityOrbitNormal(dv1Info.PropagationPoint, earth); dv1.Maneuver = dv1Info; NumericalPropagatorSegment transferOrbit = new NumericalPropagatorSegment(); transferOrbit.Name = "Transfer_Orbit"; transferOrbit.PropagatorDefinition = numericalPropagatorDefinition; // apogee stopping condition ScalarStoppingCondition apoapsisStoppingCondition = new ScalarStoppingCondition(new ScalarModifiedKeplerianElement(gravitationalParameter, integrationPoint.IntegrationPoint, KeplerianElement.TrueAnomaly, frame), Math.PI, // threshold, radians 0.00000001, // value tolerance, radians StopType.AnyThreshold); apoapsisStoppingCondition.AngularSetting = CircularRange.ZeroToTwoPi; // let the system know to avoid the branch cut apoapsisStoppingCondition.Name = "Apoapsis_Stopping_Condition"; transferOrbit.StoppingConditions.Add(apoapsisStoppingCondition); ImpulsiveManeuverSegment dv2 = new ImpulsiveManeuverSegment(); dv2.Name = "Second_Maneuver"; ImpulsiveManeuverInformation dv2Info = new ImpulsiveManeuverInformation(motionID, new Cartesian(500.0, 0.0, 0.0), fuelMassName, dryMassName, 4500.0, // exhaust velocity, meters/second InvalidFuelStateBehavior.DoFullDeltaV); dv2.Maneuver = dv2Info; dv2Info.Orientation = new AxesVelocityOrbitNormal(dv2Info.PropagationPoint, earth); innerTargetSegment.Segments.Add(dv1); innerTargetSegment.Segments.Add(transferOrbit); innerTargetSegment.Segments.Add(dv2); TargetedSegmentListDifferentialCorrector hohmannTransferDifferentialCorrector = new TargetedSegmentListDifferentialCorrector(); hohmannTransferDifferentialCorrector.Name = "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 (variable, configuration) => { configuration.Maneuver.X += variable; }); dv1VelocityXVariable.Name = "ImpulsiveManeuver1_ThrustVector_X"; // constraint 1, dv1's ScalarModifiedKeplerianElement apoapsisScalar = new ScalarModifiedKeplerianElement(gravitationalParameter, integrationPoint.IntegrationPoint, KeplerianElement.RadiusOfApoapsis, frame); ScalarAtEndOfNumericalSegmentConstraint radiusOfApoapsisPostDv1Constraint = new ScalarAtEndOfNumericalSegmentConstraint(apoapsisScalar, transferOrbit, 0, // desired value, meters, leaving as null since a outer variable will set this 0.0001); // tolerance, meters radiusOfApoapsisPostDv1Constraint.Name = "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 (variable, configuration) => { configuration.Maneuver.X += variable; }); dv2VelocityXVariable.Name = "ImpulsiveManeuver2_ThrustVector_X"; // constraint 2, eccentricity constraint after the second maneuver ScalarModifiedKeplerianElement eccentricityScalar = new ScalarModifiedKeplerianElement(gravitationalParameter, integrationPoint.IntegrationPoint, KeplerianElement.Eccentricity, frame); ScalarAtEndOfNumericalSegmentConstraint eccentricityPostDv2Constraint = new ScalarAtEndOfNumericalSegmentConstraint(eccentricityScalar, dv2, 0.1, // desired value, unitless 0.000001); // tolerance, unitless eccentricityPostDv2Constraint.Name = "Eccentricity"; // add the variables hohmannTransferDifferentialCorrector.Variables.Add(dv1VelocityXVariable); hohmannTransferDifferentialCorrector.Variables.Add(dv2VelocityXVariable); // add the constraints hohmannTransferDifferentialCorrector.Constraints.Add(radiusOfApoapsisPostDv1Constraint); hohmannTransferDifferentialCorrector.Constraints.Add(eccentricityPostDv2Constraint); // add the corrector to the targeted segment innerTargetSegment.Operators.Add(hohmannTransferDifferentialCorrector); TargetedSegmentList outerTargetedSegmentList = new TargetedSegmentList(); outerTargetedSegmentList.Name = "Outer_Target_List"; TargetedSegmentListDifferentialCorrector solverFor24HourPeriod = new TargetedSegmentListDifferentialCorrector(); solverFor24HourPeriod.Name = "Solving_For_24_Hour_Period"; // 24 hour period constraint ScalarModifiedKeplerianElement periodScalar = new ScalarModifiedKeplerianElement(gravitationalParameter, integrationPoint.IntegrationPoint, KeplerianElement.Period, frame); ScalarAtEndOfNumericalSegmentConstraint periodConstraint = new ScalarAtEndOfNumericalSegmentConstraint(periodScalar, dv2, TimeConstants.SecondsPerDay, // desired value, seconds 0.0001); // tolerance, seconds periodConstraint.Name = "Orbital_Period"; // vary the desired value of the nested Radius of Apoapsis ParameterizedDoubleVariable nestedDesiredRadiusOfApoapsisVariable = new ParameterizedDoubleVariable(42000000, // initial value, meters 5000000.0, // maximum step, meters 500000.0, // perturbation, meters transferOrbit); nestedDesiredRadiusOfApoapsisVariable.Name = "Desired_Radius_Of_Apoapsis"; radiusOfApoapsisPostDv1Constraint.DesiredValue = nestedDesiredRadiusOfApoapsisVariable.Value; // add the variable and constraint solverFor24HourPeriod.Constraints.Add(periodConstraint); solverFor24HourPeriod.Variables.Add(nestedDesiredRadiusOfApoapsisVariable); // add the differentialCorrector to the segment outerTargetedSegmentList.Operators.Add(solverFor24HourPeriod); NumericalPropagatorSegment finalPropagator = new NumericalPropagatorSegment(); finalPropagator.Name = "Final_Orbit"; finalPropagator.PropagatorDefinition = numericalPropagatorDefinition; finalPropagator.StoppingConditions.Add(new DurationStoppingCondition(Duration.FromDays(1.0))); finalPropagator.StoppingConditions[0].Name = "Duration"; SegmentList overallList = new SegmentList(); overallList.Name = "Overall_List"; overallList.Segments.Add(initial); overallList.Segments.Add(propagateToPerigeeSegment); // add the inner targeted segment list to the outer targeted segment list outerTargetedSegmentList.Segments.Add(innerTargetSegment); // and add the outer list overallList.Segments.Add(outerTargetedSegmentList); overallList.Segments.Add(finalPropagator); SegmentPropagator propagator = overallList.GetSegmentPropagator(new EvaluatorGroup()); SegmentListResults propagatorResults = (SegmentListResults)propagator.Propagate();