Skip to content

Commit

Permalink
Eliminate mass costate entirely from PVG calcs
Browse files Browse the repository at this point in the history 
This is now likely never going to happen

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed Feb 27, 2024
1 parent 530ff86 commit b6ee029
Show file tree
Hide file tree
Showing 9 changed files with 26 additions and 105 deletions.
4 changes: 2 additions & 2 deletions MechJeb2/MechJebModulePVGGlueBall.cs
View file
Expand Up @@ -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)
Expand Down
7 changes: 2 additions & 5 deletions MechJebLib/PVG/ContinuityLayout.cs
View file
Expand Up @@ -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<double> other)
Expand All @@ -27,8 +26,7 @@ public void CopyTo(IList<double> 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<double> other)
Expand All @@ -38,8 +36,7 @@ public void CopyFrom(IList<double> 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<double> other)
Expand Down
9 changes: 2 additions & 7 deletions MechJebLib/PVG/InputLayout.cs
View file
Expand Up @@ -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;

Expand All @@ -34,8 +33,6 @@ public struct InputLayout

public double M;

public double Pm;

public double H0
{
get
Expand All @@ -55,7 +52,6 @@ public void CopyTo(IList<double> other)
PV.CopyTo(other, 6);
PR.CopyTo(other, 9);
other[M_INDEX] = M;
other[PM_INDEX] = Pm;
other[BT_INDEX] = Bt;
}

Expand All @@ -66,7 +62,6 @@ public void CopyFrom(IList<double> other)
PV.CopyFrom(other, 6);
PR.CopyFrom(other, 9);
M = other[M_INDEX];
Pm = other[PM_INDEX];
Bt = other[BT_INDEX];
}

Expand Down
3 changes: 1 addition & 2 deletions MechJebLib/PVG/Integrators/VacuumCoastAnalytic.cs
View file
Expand Up @@ -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;

Expand Down
2 changes: 0 additions & 2 deletions MechJebLib/PVG/Integrators/VacuumThrustAnalytic.cs
View file
Expand Up @@ -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);
}

Expand Down
4 changes: 0 additions & 4 deletions MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs
View file
Expand Up @@ -43,10 +43,6 @@ public void dydt(IList<double> yin, double x, IList<double> 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);
Expand Down
85 changes: 14 additions & 71 deletions MechJebLib/PVG/Optimizer.cs
View file
Expand Up @@ -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;
Expand All @@ -31,9 +31,9 @@ public partial class Optimizer : IDisposable
private readonly List<Vn> _terminal = new List<Vn>();
private readonly List<Vn> _residual = new List<Vn>();
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 }

Expand Down Expand Up @@ -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]);
Expand Down Expand Up @@ -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];
}
}

Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -398,8 +347,6 @@ private void Shooting(Solution? solution = null)

t0 += bt;
}

IntegrateMassCostate();
}

public Optimizer Bootstrap(V3 pv0, V3 pr0)
Expand Down Expand Up @@ -454,8 +401,6 @@ public Optimizer Bootstrap(V3 pv0, V3 pr0)
t0 += tf;
}

IntegrateMassCostate();

CalculateResiduals();

for (int p = 0; p <= lastPhase; p++)
Expand Down Expand Up @@ -544,8 +489,6 @@ public Optimizer Bootstrap(Solution solution)
t0 += tf;
}

IntegrateMassCostate();

CalculateResiduals();

for (int p = 0; p <= lastPhase; p++)
Expand Down
10 changes: 2 additions & 8 deletions MechJebLib/PVG/OutputLayout.cs
View file
Expand Up @@ -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;

Expand All @@ -28,8 +27,6 @@ public struct OutputLayout

public double M;

public double Pm;

public double H0
{
get
Expand All @@ -49,7 +46,6 @@ public OutputLayout(InputLayout other)
PV = other.PV;
PR = other.PR;
M = other.M;
Pm = other.Pm;
DV = 0;
}

Expand All @@ -60,7 +56,6 @@ public void CopyTo(IList<double> other)
PV.CopyTo(other, 6);
PR.CopyTo(other, 9);
other[M_INDEX] = M;
other[PM_INDEX] = Pm;
other[DV_INDEX] = DV;
}

Expand All @@ -71,7 +66,6 @@ public void CopyFrom(IList<double> other)
PV.CopyFrom(other, 6);
PR.CopyFrom(other, 9);
M = other[M_INDEX];
Pm = other[PM_INDEX];
DV = other[DV_INDEX];
}

Expand Down
7 changes: 3 additions & 4 deletions MechJebLibTest/PVGTests/AscentTests/BuggyTests.cs
View file
Expand Up @@ -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();

Expand All @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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);
}
}
}

0 comments on commit b6ee029

Please sign in to comment.