From b6ee0296e25747c8d23e5e2a19c7613d3b3036cc Mon Sep 17 00:00:00 2001 From: Lamont Granquist Date: Tue, 27 Feb 2024 11:22:16 -0800 Subject: [PATCH] Eliminate mass costate entirely from PVG calcs This is now likely never going to happen Signed-off-by: Lamont Granquist --- MechJeb2/MechJebModulePVGGlueBall.cs | 4 +- MechJebLib/PVG/ContinuityLayout.cs | 7 +- MechJebLib/PVG/InputLayout.cs | 9 +- .../PVG/Integrators/VacuumCoastAnalytic.cs | 3 +- .../PVG/Integrators/VacuumThrustAnalytic.cs | 2 - .../PVG/Integrators/VacuumThrustIntegrator.cs | 4 - MechJebLib/PVG/Optimizer.cs | 85 +++---------------- MechJebLib/PVG/OutputLayout.cs | 10 +-- .../PVGTests/AscentTests/BuggyTests.cs | 7 +- 9 files changed, 26 insertions(+), 105 deletions(-) diff --git a/MechJeb2/MechJebModulePVGGlueBall.cs b/MechJeb2/MechJebModulePVGGlueBall.cs index e804a0a23..6c515fec8 100644 --- a/MechJeb2/MechJebModulePVGGlueBall.cs +++ b/MechJeb2/MechJebModulePVGGlueBall.cs @@ -86,8 +86,8 @@ private void HandleDoneTask() Debug.Log("failed guidance, znorm: " + pvg.Znorm); } - LastLmStatus = pvg.LmStatus; - LastLmIterations = pvg.LmIterations; + LastLmStatus = pvg.TerminationType; + LastLmIterations = pvg.Iterations; LastZnorm = pvg.Znorm; if (LastLmIterations > MaxLmIterations) diff --git a/MechJebLib/PVG/ContinuityLayout.cs b/MechJebLib/PVG/ContinuityLayout.cs index a8e0dc16f..f1801830f 100644 --- a/MechJebLib/PVG/ContinuityLayout.cs +++ b/MechJebLib/PVG/ContinuityLayout.cs @@ -17,7 +17,6 @@ public class ContinuityLayout public V3 Pv; public V3 Pr; public double M; - public double Pm; public double Bt; public void CopyTo(IList other) @@ -27,8 +26,7 @@ public void CopyTo(IList other) Pv.CopyTo(other, 6); Pr.CopyTo(other, 9); other[12] = M; - other[13] = Pm; - other[14] = Bt; + other[13] = Bt; } public void CopyFrom(IList other) @@ -38,8 +36,7 @@ public void CopyFrom(IList other) Pv.CopyFrom(other, 6); Pr.CopyFrom(other, 9); M = other[12]; - Pm = other[13]; - Bt = other[14]; + Bt = other[13]; } public static ContinuityLayout CreateFrom(IList other) diff --git a/MechJebLib/PVG/InputLayout.cs b/MechJebLib/PVG/InputLayout.cs index 293fce8b7..682bf70da 100644 --- a/MechJebLib/PVG/InputLayout.cs +++ b/MechJebLib/PVG/InputLayout.cs @@ -18,9 +18,8 @@ public struct InputLayout public const int VY_INDEX = 4; public const int VZ_INDEX = 5; public const int M_INDEX = 12; - public const int PM_INDEX = 13; - public const int BT_INDEX = 14; - public const int INPUT_LAYOUT_LEN = 15; + public const int BT_INDEX = 13; + public const int INPUT_LAYOUT_LEN = 14; public V3 R; @@ -34,8 +33,6 @@ public struct InputLayout public double M; - public double Pm; - public double H0 { get @@ -55,7 +52,6 @@ public void CopyTo(IList other) PV.CopyTo(other, 6); PR.CopyTo(other, 9); other[M_INDEX] = M; - other[PM_INDEX] = Pm; other[BT_INDEX] = Bt; } @@ -66,7 +62,6 @@ public void CopyFrom(IList other) PV.CopyFrom(other, 6); PR.CopyFrom(other, 9); M = other[M_INDEX]; - Pm = other[PM_INDEX]; Bt = other[BT_INDEX]; } diff --git a/MechJebLib/PVG/Integrators/VacuumCoastAnalytic.cs b/MechJebLib/PVG/Integrators/VacuumCoastAnalytic.cs index 30dd2e80e..0aa8b291e 100644 --- a/MechJebLib/PVG/Integrators/VacuumCoastAnalytic.cs +++ b/MechJebLib/PVG/Integrators/VacuumCoastAnalytic.cs @@ -23,8 +23,7 @@ public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf) yf.PV = stm00 * y0.PV - stm01 * y0.PR; yf.PR = -stm10 * y0.PV + stm11 * y0.PR; - yf.M = y0.M; - yf.Pm = y0.Pm; + yf.M = y0.M; yf.DV = y0.DV; diff --git a/MechJebLib/PVG/Integrators/VacuumThrustAnalytic.cs b/MechJebLib/PVG/Integrators/VacuumThrustAnalytic.cs index 44cb2739f..5df3c8d46 100644 --- a/MechJebLib/PVG/Integrators/VacuumThrustAnalytic.cs +++ b/MechJebLib/PVG/Integrators/VacuumThrustAnalytic.cs @@ -103,8 +103,6 @@ public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf) else yf.DV = y0.DV + phase.thrust / phase.m0 * dt; - yf.Pm = y0.Pm; // FIXME: this is certainly wrong - yf.CopyTo(yfout); } diff --git a/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs b/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs index e491d5d65..02678a851 100644 --- a/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs +++ b/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs @@ -43,10 +43,6 @@ public void dydt(IList yin, double x, IList dyout) dy.PV = -y.PR; dy.PR = y.PV / r3 - 3 / r5 * V3.Dot(y.R, y.PV) * y.R; dy.M = Phase.thrust == 0 || Phase.Infinite ? 0 : -Phase.mdot; - dy.Pm = 0; - /*dy.Pm = _phase.FinalMassProblem && !_phase.Infinite - ? _phase.ThrustBar * V3.Dot(y.PV, u) / (y.M * y.M) - : 0; */ dy.DV = at; dy.CopyTo(dyout); diff --git a/MechJebLib/PVG/Optimizer.cs b/MechJebLib/PVG/Optimizer.cs index 6d900f681..e06eb9a7b 100644 --- a/MechJebLib/PVG/Optimizer.cs +++ b/MechJebLib/PVG/Optimizer.cs @@ -21,8 +21,8 @@ public partial class Optimizer : IDisposable public double DIFFSTEP { get; set; } = 1e-9; public double STPMAX = 1e-4; public int OptimizerTimeout { get; set; } = 5000; // milliseconds - public int LmStatus; - public int LmIterations; + public int TerminationType; + public int Iterations; public OptimStatus Status; private readonly Problem _problem; @@ -31,9 +31,9 @@ public partial class Optimizer : IDisposable private readonly List _terminal = new List(); private readonly List _residual = new List(); private int lastPhase => _phases.Count - 1; - private readonly alglib.minnlcreport _rep = new alglib.minnlcreport(); + private readonly alglib.minnlcreport _rep = new alglib.minnlcreport(); private readonly alglib.ndimensional_fvec _residualHandle; - private alglib.minnlcstate _state = new alglib.minnlcstate(); + private alglib.minnlcstate _state = new alglib.minnlcstate(); public enum OptimStatus { CREATED, BOOTSTRAPPED, SUCCESS, FAILED } @@ -67,44 +67,17 @@ private double CalcBTConstraint(int p) var y0p = InputLayout.CreateFrom(_initial[p]); var yf = OutputLayout.CreateFrom(_terminal[lastPhase]); - // handle coasts - if (_phases[p].Coast && _phases[p].OptimizeTime) - { - return yfp.H0; // coast after jettison or an initial first coast - } - if (_phases[p].OptimizeTime) { if (_phases[p].Coast) - { - if (p == 0) - return yfp.H0; // initial first coast - return yfp.H0; // coast after jettison - // FIXME: coasts during stages - } - - // handle the optimized burntime that gives rise to the free final time constraint - if (_phases[p].LastFreeBurn) - { - //if (_phases[p].FinalMassProblem) return H(yf[phases.Count - 1], phases.Count - 1); - return yf.CostateMagnitude - 1; - } - - if (_phases[p].DropMass > 0 && p < lastPhase) - { - var y0p1 = OutputLayout.CreateFrom(_initial[p + 1]); - return H(yfp, p) - H(y0p1, p + 1); - } + return yfp.H0; - // any other optimized burntimes - return yfp.H0 - y0p.H0; + return yf.CostateMagnitude - 1; } return y0p.Bt - _phases[p].bt; } - 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() { var y0 = InputLayout.CreateFrom(_initial[0]); @@ -153,9 +126,9 @@ private void CalculateResiduals() private void CopyToZ(double[] z) { z[0] = 0; - for (int i = 0; i < z.Length-1; i++) + for (int i = 0; i < z.Length - 1; i++) { - z[i+1] = _residual[i / ResidualLayout.RESIDUAL_LAYOUT_LEN][i % ResidualLayout.RESIDUAL_LAYOUT_LEN]; + z[i + 1] = _residual[i / ResidualLayout.RESIDUAL_LAYOUT_LEN][i % ResidualLayout.RESIDUAL_LAYOUT_LEN]; } } @@ -259,8 +232,8 @@ private void UnSafeRun() alglib.minnlccreatef(InputLayout.INPUT_LAYOUT_LEN * _phases.Count, yGuess, DIFFSTEP, out _state); alglib.minnlcsetbc(_state, bndl, bndu); alglib.minnlcsetstpmax(_state, STPMAX); - //alglib.minnlcsetalgoslp(_state); - alglib.minnlcsetalgosqp(_state); + alglib.minnlcsetalgoslp(_state); + //alglib.minnlcsetalgosqp(_state); alglib.minnlcsetcond(_state, EPSX, MAXITS); alglib.minnlcsetnlc(_state, ResidualLayout.RESIDUAL_LAYOUT_LEN * _phases.Count, 0); alglib.minnlcoptimize(_state, _residualHandle, null, null); @@ -274,11 +247,11 @@ private void UnSafeRun() alglib.minlmresultsbuf(_state, ref yNew, _rep); */ - LmStatus = _rep.terminationtype; - LmIterations = _rep.iterationscount; + TerminationType = _rep.terminationtype; + Iterations = _rep.iterationscount; - Print("lmstatus: " + LmStatus); - Print("lmiterations: " + LmIterations); + Print("terminationtype: " + TerminationType); + Print("iterations: " + Iterations); if (_rep.terminationtype != 8) ResidualFunction(yNew, z, null); @@ -331,30 +304,6 @@ public Optimizer Run() return this; } - private void IntegrateMassCostate() - { - double pm = 1; - - for (int p = lastPhase; p >= 0; p--) - { - _terminal[p][OutputLayout.PM_INDEX] = pm; - - // FIXME: okay this is harder than I thought, so here's what needs to get done I think: - // - need dense output high fidelity interpolant from the integrator for just the Pv 3-vector. - // - need to write a real simple 1-dimensional integration kernel for Pm that uses that - // interpolant to integrate the Pm backwards. - // - then need to figure out the jump condition due to mass loss between stages. - // - alternatively it might be possible to use analyic expressions for Pv to integrate backwards - // in closed-form without any interpolant or using RK here. - // - obviously, when using analytic thrust integrals it makes sense to do exactly that, so maybe - // that is really the first step? - - _initial[p][InputLayout.PM_INDEX] = pm; - - pm += 0; // replace with jump condition between stages - } - } - private void Shooting(Solution? solution = null) { using var integArray = Vn.Rent(OutputLayout.OUTPUT_LAYOUT_LEN); @@ -398,8 +347,6 @@ private void Shooting(Solution? solution = null) t0 += bt; } - - IntegrateMassCostate(); } public Optimizer Bootstrap(V3 pv0, V3 pr0) @@ -454,8 +401,6 @@ public Optimizer Bootstrap(V3 pv0, V3 pr0) t0 += tf; } - IntegrateMassCostate(); - CalculateResiduals(); for (int p = 0; p <= lastPhase; p++) @@ -544,8 +489,6 @@ public Optimizer Bootstrap(Solution solution) t0 += tf; } - IntegrateMassCostate(); - CalculateResiduals(); for (int p = 0; p <= lastPhase; p++) diff --git a/MechJebLib/PVG/OutputLayout.cs b/MechJebLib/PVG/OutputLayout.cs index 7139cd89b..54e3c2ba2 100644 --- a/MechJebLib/PVG/OutputLayout.cs +++ b/MechJebLib/PVG/OutputLayout.cs @@ -12,9 +12,8 @@ namespace MechJebLib.PVG public struct OutputLayout { public const int M_INDEX = 12; - public const int PM_INDEX = 13; - public const int DV_INDEX = 14; - public const int OUTPUT_LAYOUT_LEN = 15; + public const int DV_INDEX = 13; + public const int OUTPUT_LAYOUT_LEN = 14; public V3 R; @@ -28,8 +27,6 @@ public struct OutputLayout public double M; - public double Pm; - public double H0 { get @@ -49,7 +46,6 @@ public OutputLayout(InputLayout other) PV = other.PV; PR = other.PR; M = other.M; - Pm = other.Pm; DV = 0; } @@ -60,7 +56,6 @@ public void CopyTo(IList other) PV.CopyTo(other, 6); PR.CopyTo(other, 9); other[M_INDEX] = M; - other[PM_INDEX] = Pm; other[DV_INDEX] = DV; } @@ -71,7 +66,6 @@ public void CopyFrom(IList other) PV.CopyFrom(other, 6); PR.CopyFrom(other, 9); M = other[M_INDEX]; - Pm = other[PM_INDEX]; DV = other[DV_INDEX]; } diff --git a/MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs b/MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs index bd77ccf42..fcbd8d2a9 100644 --- a/MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs +++ b/MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs @@ -51,7 +51,7 @@ public void SubOrbitalThor() using Solution solution = pvg.GetSolution(); pvg.Znorm.ShouldBeZero(1e-9); - Assert.Equal(8, pvg.LmStatus); + Assert.Equal(8, pvg.TerminationType); (V3 rf, V3 vf) = solution.TerminalStateVectors(); @@ -72,7 +72,6 @@ public void SubOrbitalThor() ClampPi(tanof).ShouldBeZero(1e-7); } - [Fact] // 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 @@ -106,7 +105,7 @@ private void BuggyDelta4ExtremelyHeavy() using Solution solution = pvg.GetSolution(); pvg.Znorm.ShouldBeZero(1e-9); - Assert.Equal(8, pvg.LmStatus); + Assert.Equal(8, pvg.TerminationType); } [Fact] @@ -135,7 +134,7 @@ private void BiggerEarlyRocketMaybe() using Solution solution = pvg.GetSolution(); pvg.Znorm.ShouldBeZero(1e-9); - Assert.Equal(8, pvg.LmStatus); + Assert.Equal(8, pvg.TerminationType); } } }