Skip to content

Commit

Permalink
Misc cleanup
Browse files Browse the repository at this point in the history 
- rename input/output wrapper
- was playing around with BS3 method, i don't think it'll help PVG at
  this point

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>
  • Loading branch information
@lamont-granquist
lamont-granquist committed Jun 8, 2023
1 parent 26de45d commit 9dfc8a8
Show file tree
Hide file tree
Showing 22 changed files with 113 additions and 105 deletions.
4 changes: 2 additions & 2 deletions MechJeb2/MechJeb2.csproj
View file
Expand Up @@ -121,8 +121,8 @@
<Compile Include="MechJebLib\Maneuvers\ReturnFromMoon.cs" />
<Compile Include="MechJebLib\Primitives\Vn.cs" />
<Compile Include="MechJebLib\Primitives\Scale.cs" />
<Compile Include="MechJebLib\PVG\InputWrapper.cs" />
<Compile Include="MechJebLib\PVG\OutputWrapper.cs" />
<Compile Include="MechJebLib\PVG\InputLayout.cs" />
<Compile Include="MechJebLib\PVG\OutputLayout.cs" />
<Compile Include="MechJebLib\PVG\AscentBuilder.cs" />
<Compile Include="MechJebLib\PVG\Integrators\IPVGIntegrator.cs" />
<Compile Include="MechJebLib\PVG\Integrators\VacuumCoastAnalytic.cs" />
Expand Down
8 changes: 6 additions & 2 deletions MechJeb2/MechJebLib/Core/ODE/BS3.cs
View file
Expand Up @@ -14,6 +14,10 @@ namespace MechJebLib.Core.ODE
{
using IVPFunc = Action<IList<double>, double, IList<double>>;

/*
* BS3 requires a smaller Rtol+Atol of 1e-8 and higher MaxIterations of 20,000 (10x DP5) or so.
* - Without fixing the minstepsize to always take at least a minimum floating point step, BS3 may hang.
*/
public class BS3 : AbstractRungeKutta
{
protected override int Order => 3;
Expand Down Expand Up @@ -53,11 +57,11 @@ protected override void RKStep(IVPFunc f, Vn err)
f(Ynew, T + C3 * h, K[3]);

for (int i = 0; i < N; i++)
Ynew[i] = Y[i] + h * (A41 * K[1][i] + A42 * K[2][i] + A43 * K[3][i]);
Ynew[i] = Y[i] + h * (A41 * Dy[i] + A42 * K[2][i] + A43 * K[3][i]);
f(Ynew, T + h, K[4]);

for (int i = 0; i < N; i++)
err[i] = K[1][i] * E1 + K[2][i] * E2 + K[3][i] * E3 + K[4][i] * E4;
err[i] = Dy[i] * E1 + K[2][i] * E2 + K[3][i] * E3 + K[4][i] * E4;

K[4].CopyTo(Dynew);
}
Expand Down
21 changes: 8 additions & 13 deletions MechJeb2/MechJebLib/PVG/InputWrapper.cs → MechJeb2/MechJebLib/PVG/InputLayout.cs
View file
Expand Up @@ -4,9 +4,10 @@

namespace MechJebLib.PVG
{
public struct InputWrapper
public struct InputLayout
{
public const int INPUT_WRAPPER_LEN = 15;
public const int INPUT_LAYOUT_LEN = 15;
public const int BT_INDEX = 14;

public V3 R;

Expand Down Expand Up @@ -42,7 +43,7 @@ public void CopyTo(IList<double> other)
PR.CopyTo(other, 9);
other[12] = M;
other[13] = Pm;
other[14] = Bt;
other[BT_INDEX] = Bt;
}

public void CopyFrom(IList<double> other)
Expand All @@ -53,20 +54,14 @@ public void CopyFrom(IList<double> other)
PR.CopyFrom(other, 9);
M = other[12];
Pm = other[13];
Bt = other[14];
Bt = other[BT_INDEX];
}

public static InputWrapper CreateFrom(IList<double> other)
public static InputLayout CreateFrom(IList<double> other)
{
var a = new InputWrapper();
var a = new InputLayout();

a.R.CopyFrom(other);
a.V.CopyFrom(other, 3);
a.PV.CopyFrom(other, 6);
a.PR.CopyFrom(other, 9);
a.M = other[12];
a.Pm = other[13];
a.Bt = other[14];
a.CopyFrom(other);

return a;
}
Expand Down
6 changes: 3 additions & 3 deletions MechJeb2/MechJebLib/PVG/Integrators/VacuumCoastAnalytic.cs
View file
Expand Up @@ -14,8 +14,8 @@ public class VacuumCoastAnalytic : IPVGIntegrator
{
public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf)
{
var y0 = OutputWrapper.CreateFrom(yin);
var yf = new OutputWrapper();
var y0 = OutputLayout.CreateFrom(yin);
var yf = new OutputLayout();

(V3 rf, V3 vf, M3 stm00, M3 stm01, M3 stm10, M3 stm11) = Shepperd.Solve2(1.0, tf - t0, y0.R, y0.V);

Expand All @@ -35,7 +35,7 @@ public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf)

public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf, Solution solution)
{
var interpolant = Hn.Get(OutputWrapper.OUTPUT_WRAPPER_LEN);
var interpolant = Hn.Get(OutputLayout.OUTPUT_LAYOUT_LEN);
interpolant.Add(t0, yin);
for (int i = 1; i < 21; i++)
{
Expand Down
6 changes: 3 additions & 3 deletions MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustAnalytic.cs
View file
Expand Up @@ -18,8 +18,8 @@ public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf)
{
Check.True(phase.Normalized);

var y0 = OutputWrapper.CreateFrom(yin);
var yf = new OutputWrapper();
var y0 = OutputLayout.CreateFrom(yin);
var yf = new OutputLayout();

double rm = y0.R.magnitude;

Expand Down Expand Up @@ -113,7 +113,7 @@ public void Integrate(Vn yin, Vn yfout, Phase phase, double t0, double tf)
public void Integrate(Vn y0in, Vn yfout, Phase phase, double t0, double tf, Solution solution)
{
// kinda janky way to compute interpolants with intermediate points
var interpolant = Hn.Get(OutputWrapper.OUTPUT_WRAPPER_LEN);
var interpolant = Hn.Get(OutputLayout.OUTPUT_LAYOUT_LEN);
interpolant.Add(t0, y0in);
const int SEGMENTS = 20;

Expand Down
16 changes: 9 additions & 7 deletions MechJeb2/MechJebLib/PVG/Integrators/VacuumThrustIntegrator.cs
View file
Expand Up @@ -18,16 +18,16 @@ public class VacuumThrustIntegrator : IPVGIntegrator
{
private class VacuumThrustKernel
{
public static int N => OutputWrapper.OUTPUT_WRAPPER_LEN;
public static int N => OutputLayout.OUTPUT_LAYOUT_LEN;

public Phase Phase = null!;

public void dydt(IList<double> yin, double x, IList<double> dyout)
{
Check.True(Phase.Normalized);

var y = OutputWrapper.CreateFrom(yin);
var dy = new OutputWrapper();
var y = OutputLayout.CreateFrom(yin);
var dy = new OutputLayout();

double at = Phase.thrust / y.M;
if (Phase.Infinite)
Expand Down Expand Up @@ -60,17 +60,19 @@ 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 = false;
_solver.Rtol = 0;
_solver.ThrowOnMaxIter = true;
_solver.Maxiter = 20000;
_solver.Rtol = 1e-9;
_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 = false;
_solver.Rtol = 0;
_solver.ThrowOnMaxIter = true;
_solver.Maxiter = 2000;
_solver.Rtol = 1e-9;
_solver.Atol = 1e-9;
_ode.Phase = phase;
var interpolant = Hn.Get(VacuumThrustKernel.N);
Expand Down
69 changes: 41 additions & 28 deletions MechJeb2/MechJebLib/PVG/Optimizer.cs
View file
Expand Up @@ -48,24 +48,24 @@ private Optimizer(Problem problem, IEnumerable<Phase> phases)
private void ExpandArrays()
{
while (_initial.Count < _phases.Count)
_initial.Add(Vn.Rent(OutputWrapper.OUTPUT_WRAPPER_LEN));
_initial.Add(Vn.Rent(OutputLayout.OUTPUT_LAYOUT_LEN));
while (_terminal.Count < _phases.Count)
_terminal.Add(Vn.Rent(OutputWrapper.OUTPUT_WRAPPER_LEN));
_terminal.Add(Vn.Rent(OutputLayout.OUTPUT_LAYOUT_LEN));
while (_residual.Count < _phases.Count)
_residual.Add(Vn.Rent(ResidualWrapper.RESIDUAL_WRAPPER_LEN));
}

private void CopyToInitial(double[] yin)
{
for (int i = 0; i < yin.Length; i++)
_initial[i / OutputWrapper.OUTPUT_WRAPPER_LEN][i % OutputWrapper.OUTPUT_WRAPPER_LEN] = yin[i];
_initial[i / OutputLayout.OUTPUT_LAYOUT_LEN][i % OutputLayout.OUTPUT_LAYOUT_LEN] = yin[i];
}

private double CalcBTConstraint(int p)
{
var yfp = OutputWrapper.CreateFrom(_terminal[p]);
var y0p = InputWrapper.CreateFrom(_initial[p]);
var yf = OutputWrapper.CreateFrom(_terminal[lastPhase]);
var yfp = OutputLayout.CreateFrom(_terminal[p]);
var y0p = InputLayout.CreateFrom(_initial[p]);
var yf = OutputLayout.CreateFrom(_terminal[lastPhase]);

// handle coasts
if (_phases[p].Coast && _phases[p].OptimizeTime)
Expand All @@ -92,7 +92,7 @@ private double CalcBTConstraint(int p)

if (_phases[p].DropMass > 0 && p < lastPhase)
{
var y0p1 = OutputWrapper.CreateFrom(_initial[p + 1]);
var y0p1 = OutputLayout.CreateFrom(_initial[p + 1]);
return H(yfp, p) - H(y0p1, p + 1);
}

Expand All @@ -103,15 +103,15 @@ private double CalcBTConstraint(int p)
return y0p.Bt - _phases[p].bt;
}

private double H(OutputWrapper y, int p)
private double H(OutputLayout y, int p)
{
return y.H0 + _phases[p].thrust * y.PV.magnitude / y.M - y.Pm * _phases[p].mdot;
}

private void BaseResiduals()
{
var y0 = InputWrapper.CreateFrom(_initial[0]);
var yf = OutputWrapper.CreateFrom(_terminal[lastPhase]);
var y0 = InputLayout.CreateFrom(_initial[0]);
var yf = OutputLayout.CreateFrom(_terminal[lastPhase]);
using var z = ResidualWrapper.Rent(_residual[0]);

z.R = y0.R - _problem.R0Bar;
Expand All @@ -126,8 +126,8 @@ private void ContinuityConditions()
{
for (int p = 1; p < _phases.Count; p++)
{
var y0 = InputWrapper.CreateFrom(_initial[p]);
var yf = OutputWrapper.CreateFrom(_terminal[p - 1]);
var y0 = InputLayout.CreateFrom(_initial[p]);
var yf = OutputLayout.CreateFrom(_terminal[p - 1]);
using var z = ResidualWrapper.Rent(_residual[p]);

z.RContinuity = yf.R - y0.R;
Expand Down Expand Up @@ -218,14 +218,27 @@ private void UnSafeRun()
AnalyzePhases();
ExpandArrays();

double[] yGuess = new double[_phases.Count * OutputWrapper.OUTPUT_WRAPPER_LEN];
double[] yNew = new double[_phases.Count * OutputWrapper.OUTPUT_WRAPPER_LEN];
double[] yGuess = new double[_phases.Count * InputLayout.INPUT_LAYOUT_LEN];
double[] yNew = new double[_phases.Count * InputLayout.INPUT_LAYOUT_LEN];
double[] z = new double[_phases.Count * ResidualWrapper.RESIDUAL_WRAPPER_LEN];
double[] bndu = new double[_phases.Count * InputLayout.INPUT_LAYOUT_LEN];
double[] bndl = new double[_phases.Count * InputLayout.INPUT_LAYOUT_LEN];

for (int i = 0; i < bndu.Length; i++)
{
bndu[i] = double.PositiveInfinity;
bndl[i] = double.NegativeInfinity;
}

for (int i = 0; i < yGuess.Length; i++)
yGuess[i] = _initial[i / OutputWrapper.OUTPUT_WRAPPER_LEN][i % OutputWrapper.OUTPUT_WRAPPER_LEN];
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;

alglib.minlmcreatev(OutputWrapper.OUTPUT_WRAPPER_LEN * _phases.Count, yGuess, LmDiffStep, out _state);
alglib.minlmcreatev(ResidualWrapper.RESIDUAL_WRAPPER_LEN * _phases.Count, yGuess, LmDiffStep, out _state);
alglib.minlmsetbc(_state, bndl, bndu);
alglib.minlmsetcond(_state, LmEpsx, MaxIter);
alglib.minlmoptimize(_state, _residualHandle, null, null);
alglib.minlmresultsbuf(_state, ref yNew, _rep);
Expand Down Expand Up @@ -286,9 +299,9 @@ public Optimizer Run()

private void Shooting(Solution? solution = null)
{
using var integArray = Vn.Rent(OutputWrapper.OUTPUT_WRAPPER_LEN);
using var integArray = Vn.Rent(OutputLayout.OUTPUT_LAYOUT_LEN);

var y0 = new InputWrapper();
var y0 = new InputLayout();

double t0 = 0;
double lastDv = 0;
Expand All @@ -312,7 +325,7 @@ private void Shooting(Solution? solution = null)
double tf = t0 + bt;
phase.u0 = GetIntertialHeading(p, y0.PV);

var y0p = new OutputWrapper(y0);
var y0p = new OutputLayout(y0);
y0p.DV = lastDv;
y0p.CopyTo(integArray);

Expand All @@ -321,7 +334,7 @@ private void Shooting(Solution? solution = null)
else
phase.Integrate(integArray, _terminal[p], t0, tf);

var yf = OutputWrapper.CreateFrom(_terminal[p]);
var yf = OutputLayout.CreateFrom(_terminal[p]);
lastDv = yf.DV;

t0 += bt;
Expand All @@ -336,7 +349,7 @@ public Optimizer Bootstrap(V3 pv0, V3 pr0)

ExpandArrays();

using var integArray = Vn.Rent(OutputWrapper.OUTPUT_WRAPPER_LEN);
using var integArray = Vn.Rent(OutputLayout.OUTPUT_LAYOUT_LEN);

double t0 = 0;
double lastDv = 0;
Expand All @@ -345,7 +358,7 @@ public Optimizer Bootstrap(V3 pv0, V3 pr0)
{
Phase phase = _phases[p];

var y0 = new InputWrapper();
var y0 = new InputLayout();

if (p == 0)
{
Expand All @@ -369,13 +382,13 @@ public Optimizer Bootstrap(V3 pv0, V3 pr0)
double tf = t0 + y0.Bt;
phase.u0 = GetIntertialHeading(p, y0.PV);

var y0p = new OutputWrapper(y0);
var y0p = new OutputLayout(y0);
y0p.DV = lastDv;
y0p.CopyTo(integArray);

phase.Integrate(integArray, _terminal[p], t0, tf);

var yf = OutputWrapper.CreateFrom(_terminal[p]);
var yf = OutputLayout.CreateFrom(_terminal[p]);
lastDv = yf.DV;

t0 += tf;
Expand Down Expand Up @@ -419,7 +432,7 @@ public Optimizer Bootstrap(Solution solution)

ExpandArrays();

using var integArray = Vn.Rent(OutputWrapper.OUTPUT_WRAPPER_LEN);
using var integArray = Vn.Rent(OutputLayout.OUTPUT_LAYOUT_LEN);

//double tbar = solution.Tbar(_problem.t0);

Expand All @@ -430,7 +443,7 @@ public Optimizer Bootstrap(Solution solution)
{
Phase phase = _phases[p];

var y0 = new InputWrapper();
var y0 = new InputLayout();

if (p == 0)
{
Expand All @@ -455,13 +468,13 @@ public Optimizer Bootstrap(Solution solution)

phase.u0 = GetIntertialHeading(p, y0.PV);

var y0p = new OutputWrapper(y0);
var y0p = new OutputLayout(y0);
y0p.DV = lastDv;
y0p.CopyTo(integArray);

phase.Integrate(integArray, _terminal[p], t0, tf);

var yf = OutputWrapper.CreateFrom(_terminal[p]);
var yf = OutputLayout.CreateFrom(_terminal[p]);
lastDv = yf.DV;

t0 += tf;
Expand Down

0 comments on commit 9dfc8a8

Please sign in to comment.