Challenge Mars orbital transfer

MontBlanc2012

Member
Joined
Aug 13, 2017
Messages
136
Reaction score
19
Points
18
Here is a straightforward challenge designed to test skills in orbital mechanics:

A stock Delta Glider is in a circular orbit around Mars with an orbital radius of 9,000 km and coplanar with the plane of the ecliptic. For obscure operational reasons, the Delta Glider needs to transfer to a an orbit as follows:

Periapsis radius - 3,960 km
Apoapsis radius - 4,500 km
Inclination - 45 degrees
Longitude of ascending node - 135 degrees
Argument of periapsis - 65 degree

Ordinarily, this (somewhat arbitrary) orbit transfer could be constructed from a plane change manoeuvre; a retrograde burn to lower periapsis; and then another retrograde burn to lower apoapsis - so, three sequential burns in all.

However, the Delta Glider in question is fuel-constrained with just 2,030 m/s of dV for the main engines and 10 m/s in the RCS tanks. The goal of this challenge is, then, to find a more efficient way of making the orbital transfer - one that fits with this rather tight fuel budget.

The scenario for this challenge is given below. Just cut & paste into a txt file and save in the Orbiter Scenarios file as a .scn file.

Code:
BEGIN_DESC

END_DESC

BEGIN_ENVIRONMENT
  System Sol
  Date MJD 51982.0320109622
  Help CurrentState_img
END_ENVIRONMENT

BEGIN_FOCUS
  Ship GL-02
END_FOCUS

BEGIN_CAMERA
  TARGET GL-02
  MODE Extern
  POS 4.000000 -84.078820 -30.845422
  TRACKMODE TargetRelative
  FOV 50.00
END_CAMERA

BEGIN_SHIPS
GL-02:DeltaGlider
  STATUS Orbiting Mars
  RPOS -6436164.352 0.020 6290929.391
  RVEL -1524.8126 0.0000 -1560.0150
  AROT -11.011 -49.870 -19.751
  AFCMODE 7
  PRPLEVEL 0:0.044404 1:0.005
  NAVFREQ 586 466 0 0
  XPDR 0
  HOVERHOLD 0 1 0.0000e+000 0.0000e+000
  AAP 0:0 0:0 0:0
  SKIN BLUE
END
END_SHIPS


---------- Post added 02-18-18 at 12:58 AM ---------- Previous post was 02-17-18 at 06:29 AM ----------

SPOILER ALERT:

For those interested: one solution to the challenge is a two-burn solution which, in TransX coordinates, is given below:

Code:
Burn 1:
   MJD:                  51982.18271 
   Prograde:              -890.905 m/s
   Outward:                 74.329 m/s
   Plane:                 1105.666 m/s

Burn 2:
   MJD:                  51982.27468 
   Prograde:              -542.827 m/s
   Outward:                -32.921 m/s
   Plane:                 -261.894 m/s

The first burn is 1422 m/s and the second burn is 604 m/s making a total of around 2028 m/s. You can try this sequence out and you should end up with an insertion into the target orbit to within about 0.01 degrees. If anyone can do better than this, please report.
 

Ajaja

Active member
Joined
Apr 20, 2008
Messages
190
Reaction score
55
Points
28
It looks like the best possible total delta-V here is around 2004 m/s.

There are Delta Velocity programs for IMFD in P30 LVLH frame:
Burn 1: 51982.1825373219
dVf: -864.817 m/s
dVp: 1121.245 m/s

Burn 2: 51982.2734755480
dVf: -525.079 m/s
dVp: -263.381 m/s

Total: 1416 + 588 = 2004 m/s

A script for GMAT for solving this problem:
Code:
%General Mission Analysis Tool(GMAT) Script
%Created: 2018-08-23 14:13:29


%----------------------------------------
%---------- Spacecraft
%----------------------------------------

Create Spacecraft DefaultSC;
GMAT DefaultSC.DateFormat = UTCModJulian;
GMAT DefaultSC.Epoch = '21982.5320109622';
GMAT DefaultSC.CoordinateSystem = MarsMJ2000Ec;
GMAT DefaultSC.DisplayStateType = Cartesian;
GMAT DefaultSC.X = -6436164.352e-3;
GMAT DefaultSC.Y = 6290929.391e-3;
GMAT DefaultSC.Z = 0.020e-3;
GMAT DefaultSC.VX = -1524.8126e-3;
GMAT DefaultSC.VY = -1560.0150e-3;
GMAT DefaultSC.VZ = 0;


%----------------------------------------
%---------- ForceModels
%----------------------------------------

Create ForceModel DefaultProp_ForceModel;
GMAT DefaultProp_ForceModel.CentralBody = Mars;
GMAT DefaultProp_ForceModel.PointMasses = {Mars, Sun};

%----------------------------------------
%---------- Propagators
%----------------------------------------

Create Propagator DefaultProp;
GMAT DefaultProp.FM = DefaultProp_ForceModel;
GMAT DefaultProp.Type = RungeKutta89;

%----------------------------------------
%---------- Burns
%----------------------------------------

Create ImpulsiveBurn ImpulsiveBurn1;
GMAT ImpulsiveBurn1.CoordinateSystem = Local;
GMAT ImpulsiveBurn1.Origin = Mars;
GMAT ImpulsiveBurn1.Axes = LVLH;

Create ImpulsiveBurn ImpulsiveBurn2;
GMAT ImpulsiveBurn2.CoordinateSystem = Local;
GMAT ImpulsiveBurn2.Origin = Mars;
GMAT ImpulsiveBurn2.Axes = LVLH;

%----------------------------------------
%---------- Coordinate Systems
%----------------------------------------

Create CoordinateSystem MarsMJ2000Ec;
GMAT MarsMJ2000Ec.Origin = Mars;
GMAT MarsMJ2000Ec.Axes = MJ2000Ec;

%----------------------------------------
%---------- Solvers
%----------------------------------------

Create Yukon Yukon;

%----------------------------------------
%---------- Arrays, Variables, Strings
%----------------------------------------
Create Variable dV e t1 t2;
GMAT dV = 0;
GMAT e = 0;
GMAT t1 = 0;
GMAT t2 = 0;

%----------------------------------------
%---------- Mission Sequence
%----------------------------------------

BeginMissionSequence;

GMAT e = (4500-3960)/(4500+3960);

Optimize Yukon {SolveMode = Solve, ExitMode = SaveAndContinue, ShowProgressWindow = true};
   
   Vary Yukon(t1 = 13000, {MultiplicativeScaleFactor = 1e-6});
   Vary Yukon(t2 = 7860, {MultiplicativeScaleFactor = 1e-6});
   
   If t1 < 0
      Propagate BackProp DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t1};
   Else
      Propagate DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t1};
   EndIf;
   
   
   Vary Yukon(ImpulsiveBurn1.Element2 = -0.864817);
   Vary Yukon(ImpulsiveBurn1.Element3 = -1.121245);
   
   Maneuver ImpulsiveBurn1(DefaultSC);
   
   
   If t2 < 0
      Propagate BackProp DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t2};
   Else
      Propagate DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t2};
   EndIf;
   
   Vary Yukon(ImpulsiveBurn2.Element2 = -0.525078);
   Vary Yukon(ImpulsiveBurn2.Element3 =  0.263380);
   
   
   
   Maneuver ImpulsiveBurn2(DefaultSC);
   
   
   NonlinearConstraint Yukon(DefaultSC.Mars.ECC=e);
   NonlinearConstraint Yukon(DefaultSC.Mars.SMA=4230);
   NonlinearConstraint Yukon(DefaultSC.MarsMJ2000Ec.AOP=65);
   NonlinearConstraint Yukon(DefaultSC.MarsMJ2000Ec.RAAN=135);
   NonlinearConstraint Yukon(DefaultSC.MarsMJ2000Ec.INC=45);
   
   GMAT dV = sqrt(ImpulsiveBurn1.Element1^2+ImpulsiveBurn1.Element2^2+ImpulsiveBurn1.Element3^2)+sqrt(ImpulsiveBurn2.Element1^2+ImpulsiveBurn2.Element2^2+ImpulsiveBurn2.Element3^2);
   Minimize Yukon(dV);
EndOptimize;  % For optimizer Yukon
 

MontBlanc2012

Member
Joined
Aug 13, 2017
Messages
136
Reaction score
19
Points
18
Ajaja

Congrats! I wondered if anyone was making headway with this.

Although deceptively simple, this isn't an easy challenge: one needs to recognise that this orbital transfer problem can't be solved using simple analytical techniques; and the standard Orbiter suite of trajectory planning tools aren't quite up to the task either (although trial-and-error might work - eventually). So, as you've clearly worked out, some fidelity in using fairly sophisticated non-linear optimisation techniques is required to solve the problem

Whereas you've used GMAT as the optimisation engine, I used ESA's Pykep/Pygmo optimisation suite. It looks as if your run with GMAT produced a solution that is ~ 20 m/s lower than Pykep/Pygmo. I probably didn't drive the latter hard enough (or for long enough) in my optimisation run.

Thanks for posting.
 

Ajaja

Active member
Joined
Apr 20, 2008
Messages
190
Reaction score
55
Points
28
Thanks for the challenge. It was fun.
Honestly, I thought is's better to change plane as higher as possible.
So, I had started from 3 burns and then step by step GMAT eliminated the first burn and came to the similar 2-burn solution.

Unfortunately, I've found that the IMFD's autoburn is really not convenient for this scenario and makes big mistakes.

So, there is even better 2002 m/s solution for TransX:

Man. date 51982.1823268119
Prograde vel. -865.098224443
Ch. plane vel. 1119.6158848010
Outward vel. 44.3764138729

Man. date 51982.2757792
Prograde vel. -523.10120065
Ch. plane vel. -264.301774121
Outward vel. -21.1047051872

I reached in Orbiter the target orbit with 0.02 deg precision on the first run.

The GMAT script:
Code:
%----------------------------------------
%---------- Spacecraft
%----------------------------------------

Create Spacecraft DefaultSC;
GMAT DefaultSC.DateFormat = UTCModJulian;
GMAT DefaultSC.Epoch = '21982.5320109622';
GMAT DefaultSC.CoordinateSystem = MarsMJ2000Ec;
GMAT DefaultSC.DisplayStateType = Cartesian;
GMAT DefaultSC.X = -6436.16435199976;
GMAT DefaultSC.Y = 6290.929390995172;
GMAT DefaultSC.Z = 2.000600670726271e-005;
GMAT DefaultSC.VX = -1.524812600000001;
GMAT DefaultSC.VY = -1.560015;
GMAT DefaultSC.VZ = 1.110223024625157e-016;


%----------------------------------------
%---------- ForceModels
%----------------------------------------

Create ForceModel DefaultProp_ForceModel;
GMAT DefaultProp_ForceModel.CentralBody = Mars;
GMAT DefaultProp_ForceModel.PointMasses = {Mars, Sun};

%----------------------------------------
%---------- Propagators
%----------------------------------------

Create Propagator DefaultProp;
GMAT DefaultProp.FM = DefaultProp_ForceModel;
GMAT DefaultProp.Type = RungeKutta89;

%----------------------------------------
%---------- Burns
%----------------------------------------

Create ImpulsiveBurn ImpulsiveBurn1;
GMAT ImpulsiveBurn1.CoordinateSystem = Local;
GMAT ImpulsiveBurn1.Origin = Mars;
GMAT ImpulsiveBurn1.Axes = VNB;

Create ImpulsiveBurn ImpulsiveBurn2;
GMAT ImpulsiveBurn2.CoordinateSystem = Local;
GMAT ImpulsiveBurn2.Origin = Mars;
GMAT ImpulsiveBurn2.Axes = VNB;

%----------------------------------------
%---------- Coordinate Systems
%----------------------------------------

Create CoordinateSystem MarsMJ2000Ec;
GMAT MarsMJ2000Ec.Origin = Mars;
GMAT MarsMJ2000Ec.Axes = MJ2000Ec;

%----------------------------------------
%---------- Solvers
%----------------------------------------

Create Yukon Yukon;

%----------------------------------------
%---------- Arrays, Variables, Strings
%----------------------------------------
Create Variable dV e t1 t2;
GMAT dV = 0;
GMAT e = 0;
GMAT t1 = 0;
GMAT t2 = 0;

%----------------------------------------
%---------- Mission Sequence
%----------------------------------------

BeginMissionSequence;

GMAT e = (4500-3960)/(4500+3960);

Optimize Yukon {SolveMode = Solve, ExitMode = SaveAndContinue, ShowProgressWindow = true};
   
   Vary Yukon(t1 = 12987.228081366, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1e-6});
   Vary Yukon(t2 = 8069.655123212838, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1e-6});
   
   If t1 < 0
      Propagate BackProp DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t1};
   Else
      Propagate DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t1};
   EndIf;
   
   
   Vary Yukon(ImpulsiveBurn1.Element1 = -0.8652682609766254, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
   Vary Yukon(ImpulsiveBurn1.Element2 = -1.119858871110692, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
   Vary Yukon(ImpulsiveBurn1.Element3 = 0.04327717637027173, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
   
   Maneuver ImpulsiveBurn1(DefaultSC);
   
   
   If t2 < 0
      Propagate BackProp DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t2};
   Else
      Propagate DefaultProp(DefaultSC) {DefaultSC.ElapsedSecs = t2};
   EndIf;
   
   Vary Yukon(ImpulsiveBurn2.Element1 = -0.5230668164953295, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
   Vary Yukon(ImpulsiveBurn2.Element2 = 0.2637472641019023, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
   Vary Yukon(ImpulsiveBurn2.Element3 = -0.02189236242268765, {Perturbation = 0.001, MaxStep = 9.999999e300, AdditiveScaleFactor = 0.0, MultiplicativeScaleFactor = 1.0});
   
   
   
   Maneuver ImpulsiveBurn2(DefaultSC);
   
   
   NonlinearConstraint Yukon(DefaultSC.Mars.ECC=e);
   NonlinearConstraint Yukon(DefaultSC.Mars.SMA=4230);
   NonlinearConstraint Yukon(DefaultSC.MarsMJ2000Ec.AOP=65);
   NonlinearConstraint Yukon(DefaultSC.MarsMJ2000Ec.RAAN=135);
   NonlinearConstraint Yukon(DefaultSC.MarsMJ2000Ec.INC=45);
   
   GMAT dV = sqrt(ImpulsiveBurn1.Element1^2+ImpulsiveBurn1.Element2^2+ImpulsiveBurn1.Element3^2)+sqrt(ImpulsiveBurn2.Element1^2+ImpulsiveBurn2.Element2^2+ImpulsiveBurn2.Element3^2);
   Minimize Yukon(dV);
EndOptimize;  % For optimizer Yukon
 
Last edited:
Top