From 83323fd6cb6b5b1ce4cbe95816670f6068eefe13 Mon Sep 17 00:00:00 2001 From: Lamont Granquist Date: Fri, 11 Aug 2023 13:17:19 -0700 Subject: [PATCH] Fix some PVG optimizer edge conditions This pins the initial mass to m0 of the phase so that the optimizer doesn't play around with it. That is necessary to combine with the max burntime constraint in order to avoid the rocket burning past the tau restriction and off into infinity. If m0 is free then the tau constraint becomes dynamic based on the selected value of m0 and no longer a box constraint (so this may need to change for e.g. stage and a half or coasts in the middle of a stage). Also though this flips the boolean so that the integrator shouldn't blow up when it goes past max interations, which seems to have been a mistake (although it surfaced this issue pretty nicely). Signed-off-by: Lamont Granquist --- MechJeb2/MechJebLib/PVG/InputLayout.cs | 5 ++- .../PVG/Integrators/VacuumThrustIntegrator.cs | 6 +-- MechJeb2/MechJebLib/PVG/Optimizer.cs | 20 ++++----- MechJeb2/MechJebModulePVGGlueBall.cs | 28 +++--------- MechJebLibTest/PVG/AscentTests/BuggyTests.cs | 44 +++++++++++++++++++ 5 files changed, 63 insertions(+), 40 deletions(-) diff --git a/MechJeb2/MechJebLib/PVG/InputLayout.cs b/MechJeb2/MechJebLib/PVG/InputLayout.cs index 1da4c6a83..51492181c 100644 --- a/MechJeb2/MechJebLib/PVG/InputLayout.cs +++ b/MechJeb2/MechJebLib/PVG/InputLayout.cs @@ -15,6 +15,7 @@ public struct InputLayout { public const int INPUT_LAYOUT_LEN = 15; public const int BT_INDEX = 14; + public const int M_INDEX = 12; public V3 R; @@ -48,7 +49,7 @@ public void CopyTo(IList other) V.CopyTo(other, 3); PV.CopyTo(other, 6); PR.CopyTo(other, 9); - other[12] = M; + other[M_INDEX] = M; other[13] = Pm; other[BT_INDEX] = Bt; } @@ -59,7 +60,7 @@ public void CopyFrom(IList other) V.CopyFrom(other, 3); PV.CopyFrom(other, 6); PR.CopyFrom(other, 9); - M = other[12]; + M = other[M_INDEX]; Pm = other[13]; Bt = other[BT_INDEX]; } diff --git a/MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs b/MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs index 3592a6a2b..87a6400a8 100644 --- a/MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs +++ b/MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs @@ -60,8 +60,8 @@ public void dydt(IList yin, double x, IList dyout) public void Integrate(Vn y0, Vn yf, Phase phase, double t0, double tf) { - _solver.ThrowOnMaxIter = true; - _solver.Maxiter = 20000; + _solver.ThrowOnMaxIter = false; + _solver.Maxiter = 2000; _solver.Rtol = 1e-9; _solver.Atol = 1e-9; _ode.Phase = phase; @@ -70,7 +70,7 @@ public void Integrate(Vn y0, Vn yf, Phase phase, double t0, double tf) public void Integrate(Vn y0, Vn yf, Phase phase, double t0, double tf, Solution solution) { - _solver.ThrowOnMaxIter = true; + _solver.ThrowOnMaxIter = false; _solver.Maxiter = 2000; _solver.Rtol = 1e-9; _solver.Atol = 1e-9; diff --git a/MechJeb2/MechJebLib/PVG/Optimizer.cs b/MechJeb2/MechJebLib/PVG/Optimizer.cs index 3a2380fc5..7c3c2aec8 100644 --- a/MechJeb2/MechJebLib/PVG/Optimizer.cs +++ b/MechJeb2/MechJebLib/PVG/Optimizer.cs @@ -103,10 +103,7 @@ private double CalcBTConstraint(int p) return y0p.Bt - _phases[p].bt; } - private double H(OutputLayout y, int p) - { - return y.H0 + _phases[p].thrust * y.PV.magnitude / y.M - y.Pm * _phases[p].mdot; - } + private double H(OutputLayout y, int p) => y.H0 + _phases[p].thrust * y.PV.magnitude / y.M - y.Pm * _phases[p].mdot; private void BaseResiduals() { @@ -235,8 +232,12 @@ private void UnSafeRun() yGuess[i] = _initial[i / InputLayout.INPUT_LAYOUT_LEN][i % InputLayout.INPUT_LAYOUT_LEN]; for (int i = 0; i < _phases.Count; i++) + { if (!_phases[i].Coast && !_phases[i].Infinite) bndu[i * InputLayout.INPUT_LAYOUT_LEN + InputLayout.BT_INDEX] = _phases[i].tau * 0.999; + bndl[i * InputLayout.INPUT_LAYOUT_LEN + InputLayout.M_INDEX] = _phases[i].m0; + bndu[i * InputLayout.INPUT_LAYOUT_LEN + InputLayout.M_INDEX] = _phases[i].m0; + } alglib.minlmcreatev(ResidualLayout.RESIDUAL_LAYOUT_LEN * _phases.Count, yGuess, LmDiffStep, out _state); alglib.minlmsetbc(_state, bndl, bndu); @@ -536,16 +537,11 @@ public Solution GetSolution() return solution; } - public bool Success() - { + public bool Success() => // even if we didn't terminate successfully, we're close enough to a zero to use the solution - return Znorm < 1e-5; - } + Znorm < 1e-5; - public static OptimizerBuilder Builder() - { - return new OptimizerBuilder(); - } + public static OptimizerBuilder Builder() => new OptimizerBuilder(); public void Dispose() { diff --git a/MechJeb2/MechJebModulePVGGlueBall.cs b/MechJeb2/MechJebModulePVGGlueBall.cs index 4cb27b20d..3b1d55e3d 100644 --- a/MechJeb2/MechJebModulePVGGlueBall.cs +++ b/MechJeb2/MechJebModulePVGGlueBall.cs @@ -53,30 +53,15 @@ protected override void OnModuleDisabled() { } - public override void OnFixedUpdate() - { - Core.StageStats.RequestUpdate(); - } + public override void OnFixedUpdate() => Core.StageStats.RequestUpdate(); - public override void OnStart(PartModule.StartState state) - { - GameEvents.onStageActivate.Add(HandleStageEvent); - } + public override void OnStart(PartModule.StartState state) => GameEvents.onStageActivate.Add(HandleStageEvent); - public override void OnDestroy() - { - GameEvents.onStageActivate.Remove(HandleStageEvent); - } + public override void OnDestroy() => GameEvents.onStageActivate.Remove(HandleStageEvent); - private void HandleStageEvent(int data) - { - _blockOptimizerUntilTime = VesselState.time + _ascentSettings.OptimizerPauseTime; - } + private void HandleStageEvent(int data) => _blockOptimizerUntilTime = VesselState.time + _ascentSettings.OptimizerPauseTime; - private bool IsUnguided(int s) - { - return _ascentSettings.UnguidedStages.Contains(s); - } + private bool IsUnguided(int s) => _ascentSettings.UnguidedStages.Contains(s); // FIXME: maybe this could be a callback to the Task? private void HandleDoneTask() @@ -134,10 +119,7 @@ public void SetTarget(double peR, double apR, double attR, double inclination, d HandleDoneTask(); if (_task is { IsCompleted: false }) - { - Debug.Log("active PVG task not completed"); return; - } if (_ascentSettings.OptimizeStageFlag) { diff --git a/MechJebLibTest/PVG/AscentTests/BuggyTests.cs b/MechJebLibTest/PVG/AscentTests/BuggyTests.cs index 917cd78e9..d9493f48c 100644 --- a/MechJebLibTest/PVG/AscentTests/BuggyTests.cs +++ b/MechJebLibTest/PVG/AscentTests/BuggyTests.cs @@ -61,5 +61,49 @@ public void SubOrbitalThor() argpf.ShouldEqual(1.7723635725563422, 1e-7); ClampPi(tanof).ShouldBeZero(1e-7); } + + [Fact] + private void BiggerEarlyRocketMaybe() + { + var r0 = new V3(-4230937.57027061, -3658393.88789034, 3050613.04457008); + var v0 = new V3(266.772640873606, -308.526291373473, 0.00117499917444357); + var u0 = new V3(-0.664193346276844, -0.574214958863673, 0.478669494390488); + double t0 = 134391.70903729; + double mu = 398600435436096; + double rbody = 6371000; + + Ascent ascent = Ascent.Builder() + .Initial(r0, v0, u0, t0, mu, rbody) + .SetTarget(6621000, 6623000, 6621000, 0.558505360638185, 0, 0.349065850398866, false, false) + .AddStageUsingFinalMass(51814.6957082437, 12381.5182945184, 277.252333393375, 139.843831752395, 2, 2) + .AddStageUsingFinalMass(9232.92402212159, 887.216491554419, 252.019434034823, 529.378964006269, 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(6621000, 6623000, 6621000, 0.558505360638185, 0, 0.349065850398866, false, false) + .AddStageUsingFinalMass(51814.6957082437, 12381.5182945184, 277.252333393375, 139.843831752395, 2, 2) + .AddStageUsingFinalMass(9232.92402212159, 887.216491554419, 252.019434034823, 529.378964006269, 1, 1, true) + .FixedBurnTime(false) + .OldSolution(solution) + .Build(); + + ascent2.Run(); + + Optimizer pvg2 = ascent2.GetOptimizer() ?? throw new Exception("null optimzer"); + + using Solution solution2 = pvg2.GetSolution(); + } } }