Skip to content

Commit

Permalink
PVG: Use Analytic for bootstrap and Integrator for converged
Browse files Browse the repository at this point in the history 
Makes bootstrapping fast and then makes it accurate.

Also can find issues with the integrator by having it be the
actual second pass through.  It should be very definitive bug
reports as long as people describe it properly or show a video.

Had to initialize the Vn's to zeros because otherwise weird things
happen due to some kind of reuse bug.  Someting is getting +='d or
something like that without being zero'd.  Really that's hard
enough to track down that zero'ing everything makes a lot of sense.

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed May 31, 2023
1 parent a32c207 commit 1547d2b
Show file tree
Hide file tree
Showing 6 changed files with 168 additions and 106 deletions.
45 changes: 18 additions & 27 deletions MechJeb2/MechJebLib/PVG/Ascent.cs
View file
Expand Up @@ -7,8 +7,8 @@

using System;
using System.Collections.Generic;
using MechJebLib.Core;
using MechJebLib.Primitives;
using MechJebLib.PVG.Integrators;
using static MechJebLib.Utils.Statics;

namespace MechJebLib.PVG
Expand Down Expand Up @@ -53,24 +53,7 @@ private Ascent(AscentBuilder builder)

public void Run()
{
(_smaT, _eccT) = Core.Maths.SmaEccFromApsides(_peR, _apR);

foreach (Phase phase in _phases)
{
// FIXME: the analytic coast integrator is definitely buggy so the Shepperd solver must be buggy
if (phase.Coast)
//phase.Integrator = new VacuumThrustIntegrator();
phase.Integrator = new VacuumCoastAnalytic();
else
// FIXME: make a debug setting to flip between these somewhere
phase.Integrator = new VacuumThrustAnalytic();
//phase.Integrator = new VacuumThrustIntegrator();
}

//_phases[0].Integrator = new VacuumThrustIntegrator();
//_phases[1].Integrator = new VacuumThrustIntegrator();

//_phases[3].Integrator = new VacuumThrustIntegrator();
(_smaT, _eccT) = Maths.SmaEccFromApsides(_peR, _apR);

using Optimizer.OptimizerBuilder builder = Optimizer.Builder()
.Initial(_r0, _v0, _u0, _t0, _mu, _rbody)
Expand All @@ -93,6 +76,8 @@ public void Run()

private Optimizer ConvergedOptimization(Optimizer.OptimizerBuilder builder, Solution solution)
{
ForceNumericalIntegration();

if (_fixedBurnTime)
{
ApplyEnergy(builder);
Expand All @@ -115,7 +100,7 @@ private Optimizer ConvergedOptimization(Optimizer.OptimizerBuilder builder, Solu
(V3 rf, V3 vf) = solution2.TerminalStateVectors();

(_, _, _, _, _, double tanof, _) =
Core.Maths.KeplerianFromStateVectors(_mu, rf, vf);
Maths.KeplerianFromStateVectors(_mu, rf, vf);

if (_attachAltFlag || Math.Abs(ClampPi(tanof)) < PI / 2.0)
return pvg;
Expand All @@ -138,7 +123,7 @@ private void ApplyFPA(Optimizer.OptimizerBuilder builder)
// just "SetTarget" that fixes this correctly there.
double attR = _attachAltFlag ? _attR : _peR;

(_vT, _gammaT) = Core.Maths.ConvertApsidesTargetToFPA(_peR, _apR, attR, _mu);
(_vT, _gammaT) = Maths.ConvertApsidesTargetToFPA(_peR, _apR, attR, _mu);
if (_lanflag)
builder.TerminalFPA5(attR, _vT, _gammaT, _incT, _lanT);
else
Expand All @@ -147,7 +132,7 @@ private void ApplyFPA(Optimizer.OptimizerBuilder builder)

private void ApplyKepler(Optimizer.OptimizerBuilder builder)
{
(_smaT, _eccT) = Core.Maths.SmaEccFromApsides(_peR, _apR);
(_smaT, _eccT) = Maths.SmaEccFromApsides(_peR, _apR);

if (_lanflag)
builder.TerminalKepler4(_smaT, _eccT, _incT, _lanT);
Expand All @@ -168,9 +153,9 @@ private Optimizer InitialBootstrappingFixed(Optimizer.OptimizerBuilder builder)
ApplyEnergy(builder);

// guess the initial launch direction
V3 enu = Core.Maths.ENUHeadingForInclination(_incT, _r0);
V3 enu = Maths.ENUHeadingForInclination(_incT, _r0);
enu.z = 1.0; // add 45 degrees up
V3 pvGuess = Core.Maths.ENUToECI(_r0, enu).normalized;
V3 pvGuess = Maths.ENUToECI(_r0, enu).normalized;

List<Phase> bootphases = DupPhases(_phases);

Expand Down Expand Up @@ -224,9 +209,9 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build
ApplyFPA(builder);

// guess the initial launch direction
V3 enu = Core.Maths.ENUHeadingForInclination(_incT, _r0);
V3 enu = Maths.ENUHeadingForInclination(_incT, _r0);
enu.z = 1.0; // add 45 degrees up
V3 pvGuess = Core.Maths.ENUToECI(_r0, enu).normalized;
V3 pvGuess = Maths.ENUToECI(_r0, enu).normalized;

List<Phase> bootphases = DupPhases(_phases);

Expand Down Expand Up @@ -289,11 +274,17 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build
(V3 rf, V3 vf) = solution3.TerminalStateVectors();

(_, _, _, _, _, double tanof, _) =
Core.Maths.KeplerianFromStateVectors(_mu, rf, vf);
Maths.KeplerianFromStateVectors(_mu, rf, vf);

return Math.Abs(ClampPi(tanof)) > PI / 2.0 ? pvg2 : pvg3;
}

private void ForceNumericalIntegration()
{
for (int i = 0; i < _phases.Count; i++)
_phases[i].Analytic = false;
}

private void ApplyOldBurnTimesToPhases(Solution oldSolution)
{
for (int i = 0; i < _phases.Count; i++)
Expand Down
11 changes: 8 additions & 3 deletions MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs
View file
Expand Up @@ -30,7 +30,8 @@ public void dydt(IList<double> yin, double x, IList<double> dyout)
using var dy = ArrayWrapper.Rent(dyout);

double at = Phase.thrust / y.M;
if (Phase.Infinite) at *= 2;
if (Phase.Infinite)
at *= 2;

double r2 = y.R.sqrMagnitude;
double r = Math.Sqrt(r2);
Expand All @@ -57,14 +58,18 @@ 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.ThrowOnMaxIter = false;
_solver.Rtol = 0;
_solver.Atol = 1e-9;
_ode.Phase = phase;
_solver.Solve(_ode.dydt, y0, yf, t0, tf);
}

public void Integrate(Vn y0, Vn yf, Phase phase, double t0, double tf, Solution solution)
{
_solver.ThrowOnMaxIter = true;
_solver.ThrowOnMaxIter = false;
_solver.Rtol = 0;
_solver.Atol = 1e-9;
_ode.Phase = phase;
var interpolant = Hn.Get(VacuumThrustKernel.N);
_solver.Solve(_ode.dydt, y0, yf, t0, tf, interpolant);
Expand Down
109 changes: 68 additions & 41 deletions MechJeb2/MechJebLib/PVG/Phase.cs
View file
Expand Up @@ -16,58 +16,83 @@ namespace MechJebLib.PVG
{
public class Phase
{
public double m0;
public double thrust;
public double isp;
public double mf;
public double bt;
public double maxt;
public double mint;
public double ve;
public double a0;
public double tau;
public double mdot;
public double dv;
public double DropMass = 0.0; // FIXME: unused
public bool OptimizeTime;
public bool Infinite = false;
public bool Unguided;
public bool MassContinuity = false;
public bool LastFreeBurn = false;
public bool Normalized = false;
public int KSPStage;
public IPVGIntegrator Integrator;
public bool Coast => thrust == 0;
public V3 u0;
public double m0;
public double thrust;
public double isp;
public double mf;
public double bt;
public double maxt;
public double mint;
public double ve;
public double a0;
public double tau;
public double mdot;
public double dv;
public double DropMass; // FIXME: unused
public bool OptimizeTime;
public bool Infinite = false;

private bool _analytic = true;

public bool Analytic
{
get => _analytic;
set
{
_ipvgIntegrator = null;
_analytic = value;
}
}

public bool Unguided;
public bool MassContinuity = false;
public bool LastFreeBurn = false;
public bool Normalized;
public int KSPStage;

private IPVGIntegrator? _ipvgIntegrator;

private IPVGIntegrator _integrator => _ipvgIntegrator ??= GetIntegrator();

private IPVGIntegrator GetIntegrator()
{
if (Coast)
return new VacuumCoastAnalytic();

return Analytic ? (IPVGIntegrator)new VacuumThrustAnalytic() : new VacuumThrustIntegrator();
}

public bool Coast => thrust == 0;
public V3 u0;

public Phase DeepCopy()
{
var newphase = (Phase) this.MemberwiseClone();
var newphase = (Phase)MemberwiseClone();

return newphase;
}

private Phase(double m0, double thrust, double isp, double mf, double bt, int kspStage)
{
this.KSPStage = kspStage;
this.m0 = m0;
this.thrust = thrust;
this.isp = isp;
this.mf = mf;
this.bt = bt;
ve = isp * G0;
a0 = thrust / m0;
tau = thrust == 0 ? double.PositiveInfinity : ve / a0;
mdot = ve == 0 ? 0 : thrust / ve;
dv = thrust == 0 ? 0 : -ve * Math.Log(1 - bt / tau);
OptimizeTime = false;
KSPStage = kspStage;
this.m0 = m0;
this.thrust = thrust;
this.isp = isp;
this.mf = mf;
this.bt = bt;
ve = isp * G0;
a0 = thrust / m0;
tau = thrust == 0 ? double.PositiveInfinity : ve / a0;
mdot = ve == 0 ? 0 : thrust / ve;
dv = thrust == 0 ? 0 : -ve * Math.Log(1 - bt / tau);
OptimizeTime = false;
}

public Phase Rescale(Scale scale)
{
Check.False(Normalized);

var phase = (Phase)this.MemberwiseClone();
var phase = (Phase)MemberwiseClone();

phase.ve = ve / scale.VelocityScale;
phase.tau = tau / scale.TimeScale;
Expand All @@ -86,15 +111,16 @@ public Phase Rescale(Scale scale)

public void Integrate(Vn y0, Vn yf, double t0, double tf)
{
Integrator.Integrate(y0, yf, this, t0, tf);
_integrator.Integrate(y0, yf, this, t0, tf);
}

public void Integrate(Vn y0, Vn yf, double t0, double tf, Solution solution)
{
Integrator.Integrate(y0, yf, this, t0, tf, solution);
_integrator.Integrate(y0, yf, this, t0, tf, solution);
}

public static Phase NewStageUsingFinalMass(double m0, double mf, double isp, double bt, int kspStage, bool optimizeTime = false, bool unguided = false)
public static Phase NewStageUsingFinalMass(double m0, double mf, double isp, double bt, int kspStage, bool optimizeTime = false,
bool unguided = false)
{
Check.PositiveFinite(m0);
Check.PositiveFinite(mf);
Expand All @@ -112,7 +138,8 @@ public static Phase NewStageUsingFinalMass(double m0, double mf, double isp, dou
return phase;
}

public static Phase NewStageUsingThrust(double m0, double thrust, double isp, double bt, int kspStage, bool optimizeTime = false, bool unguided = false)
public static Phase NewStageUsingThrust(double m0, double thrust, double isp, double bt, int kspStage, bool optimizeTime = false,
bool unguided = false)
{
Check.PositiveFinite(m0);
Check.PositiveFinite(thrust);
Expand Down
3 changes: 2 additions & 1 deletion MechJeb2/MechJebLib/Primitives/Vn.cs
View file
Expand Up @@ -191,7 +191,8 @@ public static Vn Rent(double[] other)
private static void Clear(Vn obj)
{
// do not resize the array to zero to avoid creating garbage
// we will return nonsense in the array values (not zeros) for performance
for (int i = 0; i < obj.Count; i++)
obj[i] = 0;
}

public static void Return(Vn obj)
Expand Down
44 changes: 43 additions & 1 deletion MechJebLibTest/PVG/AscentTests/TheStandardTests.cs
View file
Expand Up @@ -11,7 +11,7 @@ public class TheStandardTests
{
// This is fairly sensitive to initial bootstrapping and will often compute a negative coast instead of the optimal positive coast
[Fact]
public void OptmizedBoosterWithCoastElliptical()
public void OptimizedBoosterWithCoastElliptical()
{
var r0 = new V3(-521765.111703417, -5568874.59934707, 3050608.87783524);
var v0 = new V3(406.088016257895, -38.0495807832894, 0.000701038889818476);
Expand Down Expand Up @@ -63,6 +63,48 @@ public void OptmizedBoosterWithCoastElliptical()
lanf.ShouldEqual(3.0481642123046941, 1e-7);
argpf.ShouldEqual(1.8684926416804641, 1e-7);
tanof.ShouldEqual(0.091089077094440363, 1e-7);

// re-run fully integrated instead of the analytic bootstrap
Ascent ascent2 = Ascent.Builder()
.Initial(r0, v0, u0, t0, mu, rbody)
.SetTarget(PeR, ApR, PeR, incT, 0, 0, false, false)
.AddStageUsingFinalMass(49119.7842689869, 7114.2513992454, 288.000034332275, 170.308460385726, 3)
.AddStageUsingFinalMass(2848.62586760223, 1363.71123994759, 270.15767003304, 116.391834883409, 1, true)
.AddOptimizedCoast(678.290157913434, 0, 450, 1)
.AddStageUsingFinalMass(678.290157913434, 177.582604389742, 230.039271734103, 53.0805126571005, 0, false, true)
.OldSolution(solution)
.Build();

ascent2.Run();

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

using Solution solution2 = pvg2.GetSolution();

solution.Tgo(solution2.T0, 0).ShouldBePositive();
solution.Tgo(solution2.T0, 1).ShouldBePositive();
solution.Tgo(solution2.T0, 2).ShouldBePositive();
solution.Tgo(solution2.T0, 3).ShouldBePositive();

pvg2.Znorm.ShouldBeZero(1e-9);

(V3 rf2, V3 vf2) = solution2.TerminalStateVectors();

(double smaf2, double eccf2, double incf2, double lanf2, double argpf2, double tanof2, _) =
MechJebLib.Core.Maths.KeplerianFromStateVectors(mu, rf2, vf2);

solution2.R(t0).ShouldEqual(r0);
solution2.V(t0).ShouldEqual(v0);
solution2.M(t0).ShouldEqual(49119.7842689869);
solution2.Vgo(t0).ShouldEqual(9677.8444697307314, 1e-7);
solution2.Pv(t0).ShouldEqual(new V3(0.4936213839655178, -0.37228178807752599, 0.17701920733179485), 1e-7);

smaf2.ShouldEqual(11463499.98898875, 1e-7);
eccf2.ShouldEqual(0.4280978753095433, 1e-7);
incf2.ShouldEqual(incT, 1e-7);
lanf2.ShouldEqual(3.0481648247809443, 1e-7);
argpf2.ShouldEqual(1.8659359635671597, 1e-7);
tanof2.ShouldEqual(0.091785663682881768, 1e-7);
}
}
}

0 comments on commit 1547d2b

Please sign in to comment.