Skip to content

Commit

Permalink
Merge pull request #1784 from MuMech/lcg/extreme-pvg-convergence
Browse files Browse the repository at this point in the history
  • Loading branch information
@lamont-granquist
lamont-granquist committed Oct 20, 2023
2 parents 73d4653 + 9f4d60a commit cba3f84
Show file tree
Hide file tree
Showing 5 changed files with 156 additions and 27 deletions.
32 changes: 23 additions & 9 deletions MechJeb2/MechJebLib/PVG/Ascent.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@ private Optimizer InitialBootstrappingFixed(Optimizer.OptimizerBuilder builder)
bootphases2[_optimizedCoastPhase].bt = total;
}

using Optimizer pvg2 = builder.Build(bootphases2);
Optimizer pvg2 = builder.Build(bootphases2);
pvg2.Bootstrap(solution);
pvg2.Run();

Expand All @@ -216,21 +216,27 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build

List<Phase> bootphases = DupPhases(_phases);

bootphases[bootphases.Count - 1].Infinite = true;
bootphases[bootphases.Count - 1].Unguided = false;

// set the coast phase to fixed time of zero
if (_optimizedCoastPhase > -1)
{
bootphases[_optimizedCoastPhase].OptimizeTime = false;
bootphases[_optimizedCoastPhase].bt = 0;
}

// switch the optimized phase to the top stage of the rocket
bootphases[_optimizedPhase].OptimizeTime = false;

// set the top stage to infinite + optimized + unguided
bootphases[bootphases.Count - 1].Infinite = true;
bootphases[bootphases.Count - 1].Unguided = false;
bootphases[bootphases.Count - 1].OptimizeTime = true;

using Optimizer pvg = builder.Build(bootphases);
pvg.Bootstrap(pvGuess, _r0.normalized);
pvg.Run();

if (!pvg.Success())
throw new Exception("Target unreachable (bootstrapping)");
throw new Exception("Target unreachable (infinite ISP)");

using Solution solution = pvg.GetSolution();
ApplyOldBurnTimesToPhases(solution);
Expand All @@ -246,14 +252,22 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build
_phases[_optimizedCoastPhase].bt = total;
}

using Optimizer pvg2 = builder.Build(_phases);
Optimizer pvg2 = builder.Build(_phases);
pvg2.Bootstrap(solution);
pvg2.Run();

if (!pvg2.Success())
throw new Exception("Target unreachable");
{
// when analytic thrust integrals fail, the numerical thrust integrals may succeed.
ForceNumericalIntegration();
pvg2 = builder.Build(_phases);
pvg2.Bootstrap(solution);
pvg2.Run();
if (!pvg2.Success())
throw new Exception("Target unreachable");
}

if (_attachAltFlag)
if (_attachAltFlag || _eccT < 1e-4)
return pvg2;

// we have a periapsis attachment solution, redo with free attachment
Expand All @@ -262,7 +276,7 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build
ApplyKepler(builder);
ApplyOldBurnTimesToPhases(solution2);

using Optimizer pvg3 = builder.Build(_phases);
Optimizer pvg3 = builder.Build(_phases);
pvg3.Bootstrap(solution2);
pvg3.Run();

Expand Down
8 changes: 4 additions & 4 deletions MechJeb2/MechJebLib/PVG/Optimizer.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@ public partial class Optimizer : IDisposable
{
public double ZnormTerminationLevel = 1e-9;
public double Znorm;
public int MaxIter { get; set; } = 20000; // this could maybe be pushed down
public double LmEpsx { get; set; } = 1e-10; // we terminate manually at 1e-9 so could lower this?
public double LmDiffStep { get; set; } = 1e-9;
public int MaxIter { get; set; } = 200000; // rely more on the optimizertimeout instead of iterations
public double LmEpsx { get; set; } = EPS; // rely more on manual termination at znorm=1e-9
public double LmDiffStep { get; set; } = 1e-10;
public int OptimizerTimeout { get; set; } = 5000; // milliseconds
public int LmStatus;
public int LmIterations;
Expand Down Expand Up @@ -236,7 +236,7 @@ private void UnSafeRun()
{
// pin the maximum time of any finite burn phase to below the tau value of the stage
if (!_phases[i].Coast && !_phases[i].Infinite && _phases[i].OptimizeTime)
bndu[i * InputLayout.INPUT_LAYOUT_LEN + InputLayout.BT_INDEX] = _phases[i].tau * 0.999;
bndu[i * InputLayout.INPUT_LAYOUT_LEN + InputLayout.BT_INDEX] = _phases[i].tau;

// pin the time of any phase which isn't allowed to be optimized
if (!_phases[i].OptimizeTime)
Expand Down
6 changes: 6 additions & 0 deletions MechJeb2/MechJebLib/Utils/BackgroundJob.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,12 @@ protected virtual void OnTaskFaulted(Task<T> task)
private readonly Action<Task<T>> _onTaskFaultedDelegate;
private readonly Action<Task<T>> _onTaskCancelledDelegate;

public void RunSync(object? o)
{
Result = _executeDelegate(o);
ResultReady = true;
}

public void StartJob(object? o)
{
ResultReady = false;
Expand Down
51 changes: 51 additions & 0 deletions MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,57 @@ public void SubOrbitalThor()
ClampPi(tanof).ShouldBeZero(1e-7);
}

[Fact]
private void BuggyDelta4ExtremelyHeavy()
{
// this is a buggy rocket, not an actual delta4 heavy, but it is being used to test that
// some of the worst possible targets still converge and to test the fallback from failure
// with the analytic thrust integrals to doing numerical integration.
var r0 = new V3(5591854.96465599, 126079.439022067, 3050616.55737457);
var v0 = new V3(-9.1936944030452, 407.764494724287, 0.000353003400966649);
var u0 = new V3(0.877712545724556, 0.0197197822130759, 0.478781640529633);
double t0 = 81786.024077168;
double mu = 398600435436096;
double rbody = 6371000;

double decmass = 2403; // can be decreased to add more payload

Ascent ascent = Ascent.Builder()
.Initial(r0, v0, u0, t0, mu, rbody)
.SetTarget(6571000, 6571000, 6571000, 0.558505360638185, 0, 0.349065850398866, true, false)
.AddStageUsingFinalMass(731995.326793174 - decmass, 328918.196876464 - decmass, 408.091742854086, 159.080051443926, 4, 4)
.AddStageUsingFinalMass(268744.821741949 - decmass, 168530.570801559 - decmass, 408.091702159863, 118.652885602609, 2, 2)
.AddStageUsingFinalMass(40763.8839289468 - decmass, 13796.4420050909 - decmass, 465.500076480742, 1118.13132581707, 1, 1, true)
.FixedBurnTime(false)
.Build();

ascent.Run();

Optimizer pvg = ascent.GetOptimizer() ?? throw new Exception("null optimzer");

using Solution solution = pvg.GetSolution();

pvg.Znorm.ShouldBeZero(1e-9);
Assert.Equal(8, pvg.LmStatus);

// re-run fully integrated instead of the analytic bootstrap
Ascent ascent2 = Ascent.Builder()
.Initial(r0, v0, u0, t0, mu, rbody)
.SetTarget(6571000, 6571000, 6571000, 0.558505360638185, 0, 0.349065850398866, true, false)
.AddStageUsingFinalMass(731995.326793174 - decmass, 328918.196876464 - decmass, 408.091742854086, 159.080051443926, 4, 4)
.AddStageUsingFinalMass(268744.821741949 - decmass, 168530.570801559 - decmass, 408.091702159863, 118.652885602609, 2, 2)
.AddStageUsingFinalMass(40763.8839289468 - decmass, 13796.4420050909 - decmass, 465.500076480742, 1118.13132581707, 1, 1, true)
.FixedBurnTime(false)
.OldSolution(solution)
.Build();

ascent2.Run();

Optimizer pvg2 = ascent2.GetOptimizer() ?? throw new Exception("null optimzer");

using Solution solution2 = pvg2.GetSolution();
}

[Fact]
private void BiggerEarlyRocketMaybe()
{
Expand Down
86 changes: 72 additions & 14 deletions MechJebLibTest/PVGTests/AscentTests/Titan2Tests.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,16 +115,72 @@ public void Kepler3()
ClampPi(tanof).ShouldEqual(0.090826664012324976, 1e-7);
}

// target a 185x186 without getting apopasis attachment. unfortuantely, i don't have a use case that
// exercizes the logic to re-pin to the periapsis because it grabs the apopasis incorrectly.
// extreme 1010 x 1012 and still get free attachment
[Fact]
public void Kepler3NoApoapsisAttachment()
public void Kepler3ExtremeElliptical()
{
var r0 = new V3(5593203.65707947, 0, 3050526.81522927);
var v0 = new V3(0, 407.862893197274, 0);
double t0 = 0;
double PeR = 6.371e+6 + 185e+3;
double ApR = 6.371e+6 + 186e+3;
double PeR = 6.371e+6 + 1.010e+6;
double ApR = 6.371e+6 + 1.012e+6;
double rbody = 6.371e+6;
double incT = Deg2Rad(28.608);
double mu = 3.986004418e+14;

Ascent ascent = Ascent.Builder()
.AddStageUsingThrust(157355.487476332, 2340000, 301.817977905273, 148.102380138703, 4, 4)
.AddStageUsingThrust(32758.6353093992, 456100.006103516, 315.000112652779, 178.63040653022, 3, 3, true)
.Initial(r0, v0, r0.normalized, t0, mu, rbody)
.SetTarget(PeR, ApR, PeR, incT, Deg2Rad(270), 0, false, false)
.Build();

ascent.Run();

Optimizer pvg = ascent.GetOptimizer() ?? throw new Exception("null optimzer");

Solution solution = pvg.GetSolution();

(V3 rf, V3 vf) = solution.TerminalStateVectors();

(double smaf, double eccf, double incf, double lanf, double argpf, double tanof, _) =
Maths.KeplerianFromStateVectors(mu, rf, vf);

solution.R(0).ShouldEqual(r0);
solution.V(0).ShouldEqual(v0);
solution.M(0).ShouldEqual(157355.487476332);

// the lower epsilon values below here probably indicate that the free attachment solution is breaking down
solution.Constant(0).ShouldBeZero(1e-3);
solution.Constant(solution.Tmax).ShouldBeZero(1e-3);

solution.Vgo(0).ShouldEqual(18968.729888080918, 1e-7);
solution.Pv(0).ShouldEqual(new V3(0.39420604314896607, 0.024122181832433174, 0.21558209101216347), 1e-7);

double aprf = Maths.ApoapsisFromKeplerian(smaf, eccf);
double perf = Maths.PeriapsisFromKeplerian(smaf, eccf);

perf.ShouldEqual(PeR, 1e-9);
aprf.ShouldEqual(ApR, 1e-9);

pvg.Znorm.ShouldBeZero(1e-5);
smaf.ShouldEqual(7381999.9997271635, 1e-5);
eccf.ShouldEqual(0.00013546395462809413, 1e-5);
incf.ShouldEqual(incT, 1e-5);
lanf.ShouldEqual(Deg2Rad(270), 1e-2);
argpf.ShouldEqual(1.5561731710314275, 1e-5);
tanof.ShouldEqual(0.036325746839258599, 1e-5);
}

// extreme 1130 x 1132 which breaks free attachment and falls back to periapsis
[Fact]
public void Kepler3MoreExtremerElliptical()
{
var r0 = new V3(5593203.65707947, 0, 3050526.81522927);
var v0 = new V3(0, 407.862893197274, 0);
double t0 = 0;
double PeR = 6.371e+6 + 1.130e+6;
double ApR = 6.371e+6 + 1.132e+6;
double rbody = 6.371e+6;
double incT = Deg2Rad(28.608);
double mu = 3.986004418e+14;
Expand Down Expand Up @@ -153,22 +209,24 @@ public void Kepler3NoApoapsisAttachment()

solution.Constant(0).ShouldBeZero(1e-7);
solution.Constant(solution.Tmax).ShouldBeZero(1e-7);
solution.Vgo(0).ShouldEqual(8498.9352440057573, 1e-7);
solution.Pv(0).ShouldEqual(new V3(0.24009395607850137, 0.21526467329187604, 0.13094696426599889), 1e-7);

// this is 0.12 seconds before tau with an 18 kg final mass
solution.Vgo(0).ShouldEqual(27198.792190399912, 1e-7);
solution.Pv(0).ShouldEqual(new V3(0.42013374453123059, 0.013974700407340215, 0.22914045811295167), 1e-7);

double aprf = Maths.ApoapsisFromKeplerian(smaf, eccf);
double perf = Maths.PeriapsisFromKeplerian(smaf, eccf);

perf.ShouldEqual(PeR, 1e-9);
aprf.ShouldEqual(ApR, 1e-9);

pvg.Znorm.ShouldBeZero(1e-9);
smaf.ShouldEqual(6556500.0000104867, 1e-7);
eccf.ShouldEqual(7.6259809968559891E-05, 1e-7);
pvg.Znorm.ShouldBeZero(1e-7);
smaf.ShouldEqual(7502000.0002534958, 1e-7);
eccf.ShouldEqual(0.0001332978218926981, 1e-7);
incf.ShouldEqual(incT, 1e-7);
lanf.ShouldEqual(Deg2Rad(270), 1e-7);
argpf.ShouldEqual(1.6494150858682159, 1e-7);
tanof.ShouldEqual(0.09094016908398217, 1e-7);
argpf.ShouldEqual(1.6027749993714844, 1e-7);
tanof.ShouldEqual(0.00011182292330591537, 1e-7);
}

[Fact]
Expand Down Expand Up @@ -207,8 +265,8 @@ public void CircularTest()

solution.Constant(0).ShouldBeZero(1e-4);
solution.Constant(solution.Tmax).ShouldBeZero(1e-4);
solution.Vgo(0).ShouldEqual(8498.6610944186486, 1e-7);
solution.Pv(0).ShouldEqual(new V3(0.24011527786624715, 0.21525367023286618, 0.13089596115507932), 1e-7);
solution.Vgo(0).ShouldEqual(8498.6702629355023, 1e-7);
solution.Pv(0).ShouldEqual(new V3(0.24008912999871768, 0.21525232653820714, 0.1309443327389036), 1e-7);

double aprf = Maths.ApoapsisFromKeplerian(smaf, eccf);
double perf = Maths.PeriapsisFromKeplerian(smaf, eccf);
Expand Down

0 comments on commit cba3f84

Please sign in to comment.