Skip to content

Commit

Permalink
Ascent optimizer algorithm changes
Browse files Browse the repository at this point in the history 
It now always returns a fully numerical solution and has some other
tweaks

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed Feb 27, 2024
1 parent a75e950 commit eba4650
Show file tree
Hide file tree
Showing 6 changed files with 132 additions and 128 deletions.
5 changes: 1 addition & 4 deletions MechJeb2/Maneuver/TransferCalculator.cs
View file
Expand Up @@ -372,17 +372,14 @@ private void AdjustPeriapsis(ManeuverParameters maneuver, ref double utArrival)
//
alglib.minnlccreatef(VARS, x, DIFFSTEP, out alglib.minnlcstate state);
alglib.minnlcsetstpmax(state, 1e-3);
double rho = 250.0;
int outerits = 5;
alglib.minnlcsetalgoaul(state, rho, outerits);
alglib.minnlcsetalgoaul2(state, outerits);
//alglib.minnlcsetalgoslp(state);
//alglib.minnlcsetalgosqp(state);
alglib.minnlcsetcond(state, EPSX, MAXITS);

alglib.minnlcsetnlc(state, EQUALITYCONSTRAINTS, INEQUALITYCONSTRAINTS);

alglib.minnlcsetprecexactrobust(state, 0);

alglib.minnlcoptimize(state, PeriapsisObjective, null, null);
alglib.minnlcresults(state, out x, out alglib.minnlcreport rep);

Expand Down
119 changes: 91 additions & 28 deletions MechJebLib/PVG/Ascent.cs
View file
Expand Up @@ -68,6 +68,35 @@ public void Run()
}
}

public void Test(V3 pv0, V3 pr0)
{
(_smaT, _eccT) = Astro.SmaEccFromApsides(_peR, _apR);

using Optimizer.OptimizerBuilder builder = Optimizer.Builder()
.Initial(_r0, _v0, _u0, _t0, _mu, _rbody)
.TerminalConditions(_hT);

ForceNumericalIntegration();

if (_fixedBurnTime)
{
ApplyEnergy(builder);
}
else
{
if (_attachAltFlag || _eccT < 1e-4)
ApplyFPA(builder);
else
ApplyKepler(builder);
}

using Optimizer pvg = builder.Build(_phases);
pvg.Bootstrap(pv0, pr0);
pvg.Run();

_optimizer = pvg;
}

public Optimizer? GetOptimizer() => _optimizer;

private Optimizer ConvergedOptimization(Optimizer.OptimizerBuilder builder, Solution solution)
Expand Down Expand Up @@ -101,7 +130,7 @@ private Optimizer ConvergedOptimization(Optimizer.OptimizerBuilder builder, Solu
(_, _, _, _, _, double tanof, _) =
Astro.KeplerianFromStateVectors(_mu, rf, vf);

if (_attachAltFlag || Abs(ClampPi(tanof)) < PI / 2.0)
if (_attachAltFlag || _fixedBurnTime || Abs(ClampPi(tanof)) < PI / 2.0)
return pvg;

ApplyFPA(builder);
Expand Down Expand Up @@ -205,6 +234,10 @@ private Optimizer InitialBootstrappingFixed(Optimizer.OptimizerBuilder builder)

private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder builder)
{
/*
* Analytic initial bootstrapping with infinite upper stage and no coast, forced FPA attachment
*/

ApplyFPA(builder);

// guess the initial launch direction
Expand All @@ -229,13 +262,18 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build
bootphases[bootphases.Count - 1].Unguided = false;
bootphases[bootphases.Count - 1].OptimizeTime = true;

Print("*** PHASE 1: DOING INITIAL INFINITE UPPER STAGE ***");
using Optimizer pvg = builder.Build(bootphases);
pvg.Bootstrap(pvGuess, _r0.normalized);
pvg.Run();

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

/*
* Analytic bootstrapping with finite upper stage and coast, forced FPA attachment
*/

using Solution solution = pvg.GetSolution();
ApplyOldBurnTimesToPhases(solution);

Expand All @@ -250,46 +288,74 @@ private Optimizer InitialBootstrappingOptimized(Optimizer.OptimizerBuilder build
_phases[_optimizedCoastPhase].bt = total;
}

Optimizer pvg2 = builder.Build(_phases);
Print("*** PHASE 2: ADDING COAST AND FINITE UPPER STAGE ***");
using Optimizer pvg2 = builder.Build(_phases);
pvg2.Bootstrap(solution);
pvg2.Run();

/*
* Numerical re-bootstrapping with finite upper stage and coast, forced FPA attachment
*/

ForceNumericalIntegration();
Optimizer pvg3 = builder.Build(_phases);
if (!pvg2.Success())
{
// 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");
Print("*** PHASE 3: REDOING WITH NUMERICAL INTEGRATION ***");
using Solution solution2 = pvg2.GetSolution();
pvg3.Bootstrap(solution2);
pvg3.Run();

if (!pvg3.Success())
throw new Exception("Target unreachable after analytic convergence");
}
else
{
Print("*** PHASE 3: ATTEMPTING NUMERICAL INTEGRATION AFTER ANALYTICAL FAILURE ***");
pvg3.Bootstrap(solution);
pvg3.Run();

if (!pvg3.Success())
throw new Exception("Target unreachable");
}

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

// we have a periapsis attachment solution, redo with free attachment
using Solution solution2 = pvg2.GetSolution();
/*
* Numerical re-bootstrapping with finite upper stage and coast, free attachment
*/

using Solution solution3 = pvg3.GetSolution();

ApplyKepler(builder);
ApplyOldBurnTimesToPhases(solution2);
ApplyOldBurnTimesToPhases(solution3);

Optimizer pvg3 = builder.Build(_phases);
pvg3.Bootstrap(solution2);
pvg3.Run();
Print("*** PHASE 4: RELAXING TO FREE ATTACHMENT ***");
Optimizer pvg4 = builder.Build(_phases);
pvg4.Bootstrap(solution3);
pvg4.Run();

if (!pvg3.Success())
return pvg2;
if (!pvg4.Success())
{
Print("*** FREE ATTACHMENT FAILED, FALLING BACK TO PERIAPSIS ***");
return pvg3;
}

// sanity check to force back near-apoapsis attatchment back to periapsis attachment
using Solution solution3 = pvg3.GetSolution();
/*
* Sanity check to ensure free attachment did not attach to the apoapsis
*/

(V3 rf, V3 vf) = solution3.TerminalStateVectors();
using Solution solution4 = pvg4.GetSolution();

(_, _, _, _, _, double tanof, _) =
Astro.KeplerianFromStateVectors(_mu, rf, vf);
if (solution3.Vgo(solution3.T0) < solution4.Vgo(solution4.T0))
{
// this probably means apoapsis attachment or wrong side of the periapsis
Print($"*** PERIAPSIS ATTACHMENT IS MORE OPTIMAL ({solution3.Vgo(solution3.T0)} < {solution4.Vgo(solution4.T0)}) THAN FREE ATTACHMENT SOLN ***");
return pvg3;
}

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

private void ForceNumericalIntegration()
Expand All @@ -313,10 +379,7 @@ private void ApplyOldBurnTimesToPhases(Solution oldSolution)
if (oldSolution.CoastPhase(j) != _phases[i].Coast)
continue;

if (!_phases[i].Coast)
_phases[i].bt = Min(oldSolution.Bt(j, _t0), _phases[i].bt);
else
_phases[i].bt = oldSolution.Bt(j, _t0);
_phases[i].bt = Min(oldSolution.Bt(j, _t0), _phases[i].tau * 0.99);
}
}
}
Expand Down
5 changes: 3 additions & 2 deletions MechJebLib/PVG/Phase.cs
View file
Expand Up @@ -22,6 +22,7 @@ public class Phase
public double maxt;
public double mint;
public double ve;
public double c => ve; // alias for normalized ve
public readonly double a0;
public double tau;
public double mdot;
Expand Down Expand Up @@ -151,11 +152,11 @@ public Phase Rescale(Scale scale)
return phase;
}

public static Phase NewFixedCoast(double m0, double ct, int kspStage, int mjPhase) => new Phase(m0, 0, 0, m0, ct, kspStage, mjPhase);
public static Phase NewFixedCoast(double m0, double ct, int kspStage, int mjPhase) => new Phase(m0, 0, double.PositiveInfinity, m0, ct, kspStage, mjPhase);

public static Phase NewOptimizedCoast(double m0, double mint, double maxt, int kspStage, int mjPhase)
{
var phase = new Phase(m0, 0, 0, m0, mint, kspStage, mjPhase) { OptimizeTime = true, mint = mint, maxt = maxt };
var phase = new Phase(m0, 0, double.PositiveInfinity, m0, mint, kspStage, mjPhase) { OptimizeTime = true, mint = mint, maxt = maxt };
return phase;
}

Expand Down
37 changes: 4 additions & 33 deletions MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs
View file
Expand Up @@ -2,6 +2,7 @@
using MechJebLib.Functions;
using MechJebLib.Primitives;
using MechJebLib.PVG;
using MechJebLib.Utils;
using Xunit;
using Xunit.Abstractions;
using static MechJebLib.Utils.Statics;
Expand All @@ -22,6 +23,7 @@ public BuggyTests(ITestOutputHelper testOutputHelper)
[Fact]
public void SubOrbitalThor()
{
Logger.Register(o => _testOutputHelper.WriteLine((string)o));
var r0 = new V3(1470817.12150277, -5396417.72296515, 3050610.07863681);
var v0 = new V3(393.513790634158, 107.250777770179, 0.00161788515083675);
var u0 = new V3(0.23076841882078, -0.846631607518583, 0.47954249382019);
Expand Down Expand Up @@ -73,6 +75,7 @@ public void SubOrbitalThor()
[Fact]
private void BuggyDelta4ExtremelyHeavy()
{
Logger.Register(o => _testOutputHelper.WriteLine((string)o));
// 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.
Expand Down Expand Up @@ -102,28 +105,12 @@ private void BuggyDelta4ExtremelyHeavy()

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()
{
Logger.Register(o => _testOutputHelper.WriteLine((string)o));
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);
Expand All @@ -147,22 +134,6 @@ private void BiggerEarlyRocketMaybe()

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();
}
}
}
60 changes: 14 additions & 46 deletions MechJebLibTest/PVGTests/AscentTests/TheStandardTests.cs
View file
Expand Up @@ -2,17 +2,27 @@
using MechJebLib.Functions;
using MechJebLib.Primitives;
using MechJebLib.PVG;
using MechJebLib.Utils;
using Xunit;
using Xunit.Abstractions;
using static MechJebLib.Utils.Statics;

namespace MechJebLibTest.PVGTests.AscentTests
{
public class TheStandardTests
{
private readonly ITestOutputHelper _testOutputHelper;

public TheStandardTests(ITestOutputHelper testOutputHelper)
{
_testOutputHelper = testOutputHelper;
}

// This is fairly sensitive to initial bootstrapping and will often compute a negative coast instead of the optimal positive coast
[Fact]
public void OptimizedBoosterWithCoastElliptical()
{
Logger.Register(o => _testOutputHelper.WriteLine((string)o));
var r0 = new V3(-521765.111703417, -5568874.59934707, 3050608.87783524);
var v0 = new V3(406.088016257895, -38.0495807832894, 0.000701038889818476);
var u0 = new V3(-0.0820737379089317, -0.874094973679233, 0.478771328926086);
Expand Down Expand Up @@ -54,57 +64,15 @@ public void OptimizedBoosterWithCoastElliptical()
solution.R(t0).ShouldEqual(r0);
solution.V(t0).ShouldEqual(v0);
solution.M(t0).ShouldEqual(49119.7842689869);
solution.Vgo(t0).ShouldEqual(9595.3503336684062, 1e-7);
solution.Pv(t0).ShouldEqual(new V3(0.4907116486773232, -0.35249571092720933, 0.16642316543642413), 1e-7);
solution.Vgo(t0).ShouldEqual(9673.8952125677752, 1e-7);
solution.Pv(t0).normalized.ShouldEqual(new V3(0.76795016838086438, -0.57846262219206013, 0.27501551521775636), 1e-7);

smaf.ShouldEqual(11463499.98898875, 1e-7);
eccf.ShouldEqual(0.4280978753095433, 1e-7);
incf.ShouldEqual(incT, 1e-7);
lanf.ShouldEqual(3.0481642123046941, 1e-7);
argpf.ShouldEqual(1.8684926416804641, 1e-7);
lanf.ShouldEqual(3.0481683004886317, 1e-7);
argpf.ShouldEqual(1.9887149934489514, 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, 3)
.AddStageUsingFinalMass(2848.62586760223, 1363.71123994759, 270.15767003304, 116.391834883409, 1, 1, true)
.AddOptimizedCoast(678.290157913434, 0, 450, 1, 1)
.AddStageUsingFinalMass(678.290157913434, 177.582604389742, 230.039271734103, 53.0805126571005, 0, 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, _) =
Astro.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 eba4650

Please sign in to comment.