Skip to content

Commit

Permalink
Fix some PVG optimizer edge conditions
Browse files Browse the repository at this point in the history 
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 <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed Aug 11, 2023
1 parent 7f3cd29 commit 83323fd
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 40 deletions.
5 changes: 3 additions & 2 deletions MechJeb2/MechJebLib/PVG/InputLayout.cs
View file
Expand Up @@ -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;

Expand Down Expand Up @@ -48,7 +49,7 @@ public void CopyTo(IList<double> 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;
}
Expand All @@ -59,7 +60,7 @@ public void CopyFrom(IList<double> 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];
}
Expand Down
6 changes: 3 additions & 3 deletions MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs
View file
Expand Up @@ -60,8 +60,8 @@ public void dydt(IList<double> yin, double x, IList<double> 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;
Expand All @@ -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;
Expand Down
20 changes: 8 additions & 12 deletions MechJeb2/MechJebLib/PVG/Optimizer.cs
View file
Expand Up @@ -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()
{
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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()
{
Expand Down
28 changes: 5 additions & 23 deletions MechJeb2/MechJebModulePVGGlueBall.cs
View file
Expand Up @@ -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()
Expand Down Expand Up @@ -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)
{
Expand Down
44 changes: 44 additions & 0 deletions MechJebLibTest/PVG/AscentTests/BuggyTests.cs
View file
Expand Up @@ -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();
}
}
}

0 comments on commit 83323fd

Please sign in to comment.