Skip to content

Commit

Permalink
ReturnFromMoon tweaks
Browse files Browse the repository at this point in the history 
- fixes a bug in the initial guess generator
- splits the infeasibility in the guess over both SOIs by averaging
  the vsoi at the interface.
- tweaks some parameters trying to target faster convergence.

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed Aug 29, 2023
1 parent 4287b09 commit 56d7513
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 14 deletions.
27 changes: 15 additions & 12 deletions MechJeb2/MechJebLib/Maneuvers/ReturnFromMoon.cs
View file
Expand Up @@ -297,7 +297,7 @@ private static double[] GenerateInitialGuess(Args args)

// calculate the true anomaly of the SOI point on the transfer orbit (using the planet scale)
double moonSOI2 = moonSOI / moonToPlanetScale.LengthScale;
double nuSOI = TAU-SafeAcos((sma * (1 - ecc * ecc) / moonSOI2 - 1) / ecc); // FIXME: probably assumes we're going to a periapsis?
double nuSOI = TAU - SafeAcos((sma * (1 - ecc * ecc) / moonSOI2 - 1) / ecc); // FIXME: probably assumes we're going to a periapsis?

if (Maths.EccFromStateVectors(1.0, r0, v0) < 1 && Maths.ApoapsisFromStateVectors(1.0, r0, v0) < moonSOI)
{
Expand All @@ -306,18 +306,18 @@ private static double[] GenerateInitialGuess(Args args)

// now compute the lunar hyperbolic burn.
V3 vneg, vpos, rburn;
(vneg, vpos, rburn, dt) = Maths.SingleImpulseHyperbolicBurn(1.0, r0, v0, vsoiPlanet * moonToPlanetScale.VelocityScale);
(vneg, vpos, rburn, dt) = Maths.SingleImpulseHyperbolicBurn(1.0, r0, v0, (vsoiPlanet - moonV0) * moonToPlanetScale.VelocityScale);

dv = vpos - vneg;
tt1 = Maths.TimeToNextRadius(1.0, rburn, vpos, moonSOI);
dv = vpos - vneg;
tt1 = Maths.TimeToNextRadius(1.0, rburn, vpos, moonSOI);
(rsoi, vsoi) = Shepperd.Solve(1.0, tt1, rburn, vpos);
}
else
{
// if we're already on an escape trajectory, try an immediate burn and hope we're reasonably close
dt = 0;
dv = V3.zero;
tt1 = Maths.TimeToNextRadius(1.0, r0, v0, moonSOI);
dt = 0;
dv = V3.zero;
tt1 = Maths.TimeToNextRadius(1.0, r0, v0, moonSOI);
(rsoi, vsoi) = Shepperd.Solve(1.0, tt1, r0, v0);
}

Expand All @@ -334,7 +334,10 @@ private static double[] GenerateInitialGuess(Args args)
(rf, vf) = Shepperd.Solve(1.0, tt2, rsoiPlanet, vsoiPlanet3);

V3 r2Sph = rsoi.cart2sph;
V3 v2Sph = vsoi.cart2sph;
//V3 v2Sph = vsoi.cart2sph;
// taking the average of the mismatch in vsoi and spreading the infeasibility over both SOIs seems to actually
// produce better convergence properties.
V3 v2Sph = (((vsoiPlanet3 - moonV2) * moonToPlanetScale.VelocityScale + vsoi) / 2).cart2sph;

x[0] = dv.x; // maneuver x
x[1] = dv.y; // maneuver y
Expand Down Expand Up @@ -364,8 +367,8 @@ private static double[] GenerateInitialGuess(Args args)
}

private const double DIFFSTEP = 1e-9;
private const double EPSX = 1e-8;
private const int MAXITS = 1000;
private const double EPSX = 1e-4;
private const int MAXITS = 10000;

private const int NVARIABLES = 20;
private const int NEQUALITYCONSTRAINTS = 17;
Expand Down Expand Up @@ -397,7 +400,7 @@ private static (V3 V, double dt) ManeuverScaled(Args args, double dtmin = double
alglib.minnlccreate(NVARIABLES, x, out alglib.minnlcstate state);

alglib.minnlcsetbc(state, bndl, bndu);
alglib.minnlcsetstpmax(state, 1e-4);
alglib.minnlcsetstpmax(state, 1e-2);
//double rho = 1000.0;
//int outerits = 5;
//alglib.minnlcsetalgoaul(state, rho, outerits);
Expand All @@ -420,7 +423,7 @@ private static (V3 V, double dt) ManeuverScaled(Args args, double dtmin = double
alglib.minnlcresults(state, out double[] x2, out alglib.minnlcreport rep);

sw.Stop();
Print($"optimization took {sw.ElapsedMilliseconds}ms");
Print($"optimization took {sw.ElapsedMilliseconds}ms: {rep.iterationscount} iter, {rep.nfev} fev");

/*
alglib.minnlcoptguardresults(state, out alglib.optguardreport ogrep);
Expand Down
16 changes: 14 additions & 2 deletions MechJebLibTest/Maneuvers/ReturnFromMoonTests.cs
View file
Expand Up @@ -24,15 +24,15 @@ private void NextManeuverToReturnFromMoonTest()
// this forces JIT compilation(?) of the SQL solver which takes ~250ms
ChangeOrbitalElement.ChangePeriapsis(1.0, new V3(1, 0, 0), new V3(0, 1.0, 0), 1.0);

//Logger.Register(o => _testOutputHelper.WriteLine((string)o));
// Logger.Register(o => _testOutputHelper.WriteLine((string)o));
double centralMu = 398600435436096;
double moonMu = 4902800066163.8;
var moonR0 = new V3(325420116.073166, -166367503.579338, -138858150.96145);
var moonV0 = new V3(577.012296778094, 761.848508254181, 297.464594270612);
double moonSOI = 66167158.6569544;
var r0 = new V3(4198676.73768844, 5187520.71497923, -3.29371833446352);
var v0 = new V3(-666.230112925872, 539.234048888927, 0.000277598267012666);
double peR = 6.3781e6 + 60000; // 60Mm
double peR = 6.3781e6 + 60000; // 60km

// this test periodically just fails, although its about 1-in-150 right now
for (int i = 0; i < 1; i++)
Expand All @@ -48,6 +48,18 @@ private void NextManeuverToReturnFromMoonTest()
continue;

/*
r0 = new V3(3033960.04435434, -538512.74449519, -1569171.01639252);
v0 = new V3(344.550626047212, -973.108221764261, 907.813925253141);
r0 = new V3(-23744.6019871556, -935848.195057236, -2970820.75826717);
v0 = new V3(-511.683969531061, -141.692448007731, 975.982327934346);
r0 = new V3(-1358340.16497667, -2896749.43027493, -2706207.90479207);
v0 = new V3(774.176044284448, -468.922061598358, -642.571722922182);
r0 = new V3(-2370135.82005065, -78770.3300753334, 504857.052794559);
v0 = new V3(-735.735041897621, -590.266316007015, -137.364944041411);
r0 = new V3(1922287.76973033, 1688857.60199125, -1128186.04080648);
v0 = new V3(-644.272529354446, 90.1035308326145, -74.2516010414118);
Expand Down

0 comments on commit 56d7513

Please sign in to comment.