Jump to content

PEGAS - Powered Explicit Guidance Ascent System - devlog


Reddy

Recommended Posts

This PDF looks like it has the PEG algorithm from the JotAS article along with the iterative thrust integral calculation (which the JotAS paper leaves to the reader):

https://dspace.mit.edu/bitstream/handle/1721.1/35563/73828492-MIT.pdf;sequence=2

(also just interesting due to the application of sequential shutdown of engines to accomplish a lunar landing using unthrottled engines)

So what appear to be missing from this are the calculations of rthrust and vthrust, using Thomas Fill's method that I can't find anywhere.

 

Edited by Jim DiGriz
Link to comment
Share on other sites

I want to do something similar, but without using any existing approaches, made from scratch, and for the game with vanilla physics, parts and massive bodies. But, I have vast amounts of information I still have to find in the game code, and I can't seem to find any of it. Some of it, I know what it's based on, and I have good approximations, some of it I do not, but I want the exact code for all of it, retrieved myself.

Orbital stats of planets and moons, atmospheric density-altitude and pressure-altitude profiles (I think this depends on temperature too which depends on latitude and time of day), engine Isp-pressure profiles (each one has a unique curve, not linear, so I want the formula and constants behind those curves), the *full* drag system and code preferably (I have an idea for a simpler version, which assumes ship is always pointing prograde in the atmosphere, and assumes that cross-sectional area times drag coefficient is constant for each stage, and requires calculating that product experimentally for each stage in advance... it could, perhaps, save those constants along with a hash of the ship design for future reference.. I'd also have to verify that it really is mostly constant).

So, pretty much everything that isn't already obvious. I think that's it. I doubt I could persuade myself to do such a thing, if I didn't have all the data I need to do it perfectly.

Oh... it's also possible that I'll want to leverage the kind of computational power that's at my fingertips, which they didn't really have back in the early space shuttle days. kOS might not be fast enough for the precision I want, but I could probably interface it with Python and SciPy through a text file without too much hassle. I want to calculate truly optimal full ascents in the stock KSP universe, from scratch! This is the essential spirit of kOS after all, just taken to an extreme.

I am okay at examining decompiled C# games, good with code, physics, control theory, numerical analysis, all that stuff, buuut I've never made a mod for anything before. That might help to explain why I can't get the information I would need first.

Edited by Aru
Link to comment
Share on other sites

11 hours ago, Jim DiGriz said:

So what appear to be missing from this are the calculations of rthrust and vthrust, using Thomas Fill's method that I can't find anywhere.

It does not look that difficult to derive. Just fit Thrust/(m0 - kt) with a second-order polynomial, then integrate a polynomial fraction... Just need to get to my KSP notebook.

UPD: finding a "good" polynomial is a bit more difficult than it sounds. I think three conditions to find A, B and C must be correct acceleration value at t=t0, correct integral of a(t)dt (a.k.a. delta-V = -F/k*ln[1 - ktgo/m0] ) and correct integral of a(t)dsdt (because two similar integrals are calculated to get vthrust and rthrust). That gives A=F/m0 and a linear system to find B and C.

Edited by Pand5461
Link to comment
Share on other sites

I've tried to rewrite the algorithm and apply it for stock system. And the first iteration just doesn't converge. Also, vbias is insanely large.

Looking at output values, (dot lambda) magnitude on the first iteration is around 0.03 with tgo about 180 sec, so (dot lambda) * tgo has magnitude > 1.

Is the divergence of first iteration expected then because I'm outside the initial assumptions' validity range?

Link to comment
Share on other sites

On 8/18/2017 at 0:45 AM, Pand5461 said:

I've tried to rewrite the algorithm and apply it for stock system. And the first iteration just doesn't converge. Also, vbias is insanely large.

Looking at output values, (dot lambda) magnitude on the first iteration is around 0.03 with tgo about 180 sec, so (dot lambda) * tgo has magnitude > 1.

Is the divergence of first iteration expected then because I'm outside the initial assumptions' validity range?

You do have to initialize it something like this:

 

            // value of tgo from previous iteration
            double tgo_prev = 0;

            if (!initialized)
            {
                rbias = new Vector3d(0, 0, 0);
                // rd initialized to rdval-length vector 20 degrees downrange from r
                rd = QuaternionD.AngleAxis(20, -iy) * r.normalized * rdval;
                // vgo initialized to rdval-length vector perpendicular to rd, minus current v
                vgo = Vector3d.Cross(-iy, rd).normalized * vdval - v;
                tgo = 1;
            }
            else
            {
                tgo_prev = tgo;
                Vector3d dvsensed = v - vprev;
                vgo = vgo - dvsensed;
                vprev = v;
            }

although i'm debugging some wonkiness in my code with a consistent wiggle halfway through the ascent, and my tangent vector has an odd negative sign in it which I don't understand where its coming from, so i suspect i might be initializing something wrong.  you can find similar code in the MATLAB implemention.

Also here's an additional PDF (that appears to be after the UPFG PDF):  https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19760024151.pdf

As far as I can tell they moved some stuff around in the thrust integrals to reduce the operations but they're identical, the corrector also has some changes that algebraically work out to the same thing.  The rthrust and vthrust computations don't have that decomposition into phi and phidot and more closely follow the Journal of the Astronomical Sciences paper, but don't seem to change the burn a great deal.  The rgo calculation is both different from the UPFG PDF and from the JAS paper.  There's also some more details in the algorithm there -- somewhat frustratingly there are also missing details.

Link to comment
Share on other sites

@Jim DiGriz thanks for the PDF.

I figured what my problem was - I forgot to apply rgo,z correction (chapter 4.5 in 19740004402). After I included it, the convergence is good.

Examining this, I drew another conclusion.

Correction looks as follows (lambda denoted as L): rgoz = iz (S - (L, rgoxy)) / (L, iz). Then, substitute rgoxy = rgo - (iz, rgo) iz.

rgoz = iz ( S - (L, rgo) + (iz, rgo) (L, iz) ) / ( L, iz ) = iz ( S - ( L, rgo ) ) / ( L, iz ) + iz (iz, rgo).

then, the next correction (omitting the rbias):

rgo = rgoxy + rgoz = rgo - (iz, rgo) iz + iz ( S - ( L, rgo ) ) / ( L, iz ) + iz (iz, rgo) = rgo + iz ( S - ( L, rgo ) ) / ( L, iz )

meaning that a) separate calculation of rgoxy is redundant and b) iz doesn't even need to be normalized.

If your tangent vector has the wrong sign the first thing to triple-check is cross-products. Also, I'd use projection of r into target plane in the initialization instead of r itself.

As for computation of rthrust and vthrust, that's the same decomposition into phi written differently.

The main difference in rgo computation  from 19740004402 is whether rbias is added before or after downrange correction. I don't think it changes too much.

Edited by Pand5461
Link to comment
Share on other sites

Sorry for joining the party late. Thanks @Jim DiGriz for the research and confirmation that the solution I developed is sound :)

Regarding the difference between correctors in 19790048206 and 19740004402 (and hence my code) - these are actually the same equations. To see this, let's first look at 19790048206, which does:

vgo = vgo - 1.0 * vmiss

where vmiss = vp - vd = v + vthrust + vgrav - vd, so this gives:

vgo = vgo - v - vthrust - vgrav + vd.

Which means we update our vgo by (vd - v - vthrust - vgrav) in each iteration. Now look at what 19740004402 does:

vgo = vd - v - vgrav + vbias (I reduced the vgop and dvgo terms)

but vbias = vgo - vthrust (block 5) so we obtain:

vgo = vd - v - vgrav + vgo - vthrust.

Let's reorder:

vgo = vgo - vthrust + vd - v - vgrav

and what we obtain is an update of vgo by (reordering further) vd - v - vthrust - vgrav. Which is exactly the same thing :)

 

Jim, can you elaborate on the difference in how they do rthrust and vthrust?

Edited by Reddy
Link to comment
Share on other sites

@kerbinorbiter I assume it literally said "No boot file loaded! Crashing..." - if this is the case, it means PEGAS hasn't found one of the variables you are supposed to create. Take another look at the tutorial (I linked to a section of interest), make sure all 4 variables (vehicle, controls, sequence, mission) exist and are named correctly.

Link to comment
Share on other sites

@Reddy do you have an idea, by chance, how to deal with coasting phases?

In the current form, 5 orbital elements out of 6 must be specified (except for argument of periapsis). While this is good in RSS, it's next to unusable in stock.

So, to include a coast, we need two additional constraints (to determine when to start coasting and when to start circularizing). But if circularization dV is fixed (i.e. Ap and Pe of suborbital trajectory already defined), then only one constraint is needed to decide when to start coasting. The constraint is obviously the minimal tgo and I don't have the idea how to find the final altitude from this.

Link to comment
Share on other sites

@Pand5461 What a coincidence, I was about to reply to you :) First, I agree with your finding that rgoxy wouldn't have to be calculated separately, and this optimization is something I would expect to find in assembly code for the actual guidance computer. But for the purpose of explaining the logic of the solution, I think they left it in a longer but clearer form deliberately. Same with normalization of iz, it might just make it easier to explain the maths in the simplest way possible.

Regarding the usage of projection of r instead of r in the initialization, I was surprised that you noticed this as I was certain that my final code already does it. I am 100% sure I had it in place at some point, since I spent several hours on trying to come up with a good initialization method - none of the documents I read contained anything on that. Simplest answer would be that I temporarily removed this at one of the debugging sessions and forgotten to put it back in :D

Now, onto your question, but please clear a few things for me first. Why is the current state "next to unusable in stock", and how would inclusion of coasting fix it? I see a potential source of confusion in the fact that you're mistaken about the current constraints: actually only 4 of 6 orbital elements must be specified: semi-major axis, eccentricity, inclination and LAN (argument of periapsis and true anomaly are free).

Edited by Reddy
wording
Link to comment
Share on other sites

@Reddy  what you've done is AMAZING WORK! Your work hasn't got the recognition it truly deserves .

PEGAS can be a life saver for anyone who wants to take RO seriously. I haven't run it yet but I can already see the endless possibilities .I tried doing something way simpler and after months of trying I finally gave up , But thanks to people like you I may still have a chance at enjoying a super-accurate , closer to real life KSP without getting overwhelmed by the sheer complexity of things .

Go PEGAS !

 

 

Link to comment
Share on other sites

22 minutes ago, Reddy said:

Why is the current state "next to unusable in stock", and how would inclusion of coasting fix it?

Continuous burn to orbit, though possible, is suboptimal in stock. Either upper stages have reasonable TWR and need to pitch wildly, or upper stages have extremely low TWR inducing gravity losses. The real Soyuz launchers suffered this problem because they had somewhat overpowered upper stages that could not be reignited. We all know the solution - burn to a suborbital trajectory at full thrust, then circularize at full thrust. It is mathematically proven that if some maneuver can be done either by a single throttled burn or two separated burns at full thrust, the latter case is more optimal.

And there are 5 elements constrained. Inclination and LAN are defined explicitly. Altitude, speed (scalar) and angle of velocity to local horizon are three parameters, so they define SMA (from total energy), eccentricity (from rotational momentum) and anomaly (from altitude, considering SMA and ECC are known).

29 minutes ago, Reddy said:

Regarding the usage of projection of r instead of r in the initialization, I was surprised that you noticed this as I was certain that my final code already does it.

Erm... I was sure I saw the projection in your code. Actually, you didn't even need to go that far - just projecting current location on the target plane and scaling it accordingly works just as well.

What I haven't seen is the account of rotating reference frame below 100 km. Again, this is more stock-related because of faster rotation and switch to inertial frame well above atmosphere (and above my typical parking orbit).

Link to comment
Share on other sites

2 hours ago, Pand5461 said:

@Reddy do you have an idea, by chance, how to deal with coasting phases?

In the current form, 5 orbital elements out of 6 must be specified (except for argument of periapsis). While this is good in RSS, it's next to unusable in stock.

So, to include a coast, we need two additional constraints (to determine when to start coasting and when to start circularizing). But if circularization dV is fixed (i.e. Ap and Pe of suborbital trajectory already defined), then only one constraint is needed to decide when to start coasting. The constraint is obviously the minimal tgo and I don't have the idea how to find the final altitude from this.

In the PDFs there's code to manually add coast phases, there's no facility in PEG to optimize the time of the coast.  

There are approaches which are not closed-loop guidance that NASA used - OPGUID and SWITCH:  https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19710026292.pdf

For closed-loop guidance that optimizes coast phases (and can even fade-in the atmosphere) there is more recent work towards that:

https://arc.aiaa.org/doi/abs/10.2514/2.5045?journalCode=jgcd

https://arc.aiaa.org/doi/abs/10.2514/1.36084

Link to comment
Share on other sites

4 hours ago, Reddy said:

Jim, can you elaborate on the difference in how they do rthrust and vthrust?

The JAS paper uses:

rgo = rd - r - v tgo - rgrav + rbias

tlambda = J / L

vthrust = lambda [ L - lambdadot2 ( H - J tlambda ) / 2 ]

rthrust = lambda [ S - lambdadot2 ( P - 2Qtlambda + Stlambda2 ) ] + ( Q - S tlambda ) lambdadot

With lambda and lambdadot the same.

Your code has the decomposition of rgo into rgoxy and rgoz (which it seems is also in the 1976 PDF, although it seems to be computed differently) and then it uses the expression involving phi and phidot which the journal article in the 1976 PDF have dropped.

And yeah, after some algebraic substitution I figured out the correctors were identical, just rearranged.

 

2 hours ago, Pand5461 said:

What I haven't seen is the account of rotating reference frame below 100 km. Again, this is more stock-related because of faster rotation and switch to inertial frame well above atmosphere (and above my typical parking orbit).

Ugh, you may have just solved my rocket wiggle problem.  Inverse rotation.  That would explain why it very suddenly adjusts at one spot halfway through the ascent.

Also:  thanks for all the other info, but I'm still drinking coffee and I wake up slowly...

Edited by Jim DiGriz
Link to comment
Share on other sites

34 minutes ago, Jim DiGriz said:

Confirmed, my launch wiggle is at exactly 145,000m on Earth which exactly matches the inverseRotThresholdAltitude

Ugh, now how to extract this from kOS (I'm starting to suspect it's not necessarily 100 km)?

The only way I can think of is to compare SolarPrimeVector between the iterations, like this:

// On every PEG loop initialization
  wait 0.
  local dt to missiontime - peg_prev["t"].
  local SPV to SolarPrimeVector.
  if peg_prev["SPr"] <> SPV {
    local omega to body:angularvel.
    local rot to
    {
      parameter vec.
      return vec - vcrs(omega, vec)*dt.
    }.
  set peg_now["Rbias"] to rot(peg_now["Rbias"]).
  //rotate everything else you need
  }
  set peg_prev["SPr"] to SPV.
  set peg_prev["t"] to peg_prev["t"] + dt.

Jim, Thanks for the references, trying to dig into.

Edited by Pand5461
Link to comment
Share on other sites

1 minute ago, Pand5461 said:

Ugh, now how to extract this from kOS (I'm starting to suspect it's not necessarily 100 km)?

Correct, its 145km on Earth, probably 100km on Kerbin?  Its a method that can get called on the CelestialBody.  I think if the body is airless there's no inv rot threshold?  Not sure, its kind of deep KSP voodoo that isn't easy to google for.

What I'm thinking about trying is getting position and velocity out of the Orbit class instead of off of the Vessel.  Probably more expensive, but I would hope KSP is somewhat optimized for UT == now case.  Gotta eat some food first tho...

Link to comment
Share on other sites

@NathanKell gave me this formula on IRC

r = vesselState.CoM - vessel.mainBody.position;
Vector3d rot = vessel.orbit.GetRotFrameVelAtPos(mainBody, r.xzy).xzy;
v = rot + vessel.velocityD;


That should correct for the inverse rotation and krakensbane velocities.  The magnitudes seem sane, but breaks my code (because I think I need to audit my cross products again, hopefully making them look less screwy...)

Link to comment
Share on other sites

@Alpha_Mike_741 Thank you for that comment! Glad to hear such words of support :)

@Pand5461 It looks right... INC and LAN explicitly, SMA, ECC and anomaly from the remaining 3 - although in a not so very transparent way. I agree. Confusion was on my end, probably due to my lacking understanding of anomaly, which is now more thorough thanks to you pointing that out.

Regarding projections - I use that in my position corrector, yes, but not in initialization. I choose the initial rd by rotating the current r, instead of projection of r on target plane. But I suppose it would save 1 iteration of convergence at most, since right after that first iteration it is constrained to the plane by said corrector.

You and @Jim DiGriz mention an important thing: accounting of rotating reference frame. This is of course merely a KSP quirk, since the SolarPrimeVector constantly changes for some reason and the axes of the reference frame with it. However, both our current state (r, v) and desired state (rd, vd) change in the same way - shouldn't that make fixes unnecessary? Or do the intermediate vectors, like biases, thrust integrals etc. have to be explicitly corrected? I spent a while pondering on that when coding the kOS version, but decided to try running without any corrections first, see how it works and fix if it's wrong - but it doesn't look wrong. This conversation would be nice to have on github issues, too, as it has potential of improving PEGAS :wink:

I will also let you know that I don't feel strong enough in control theory to do anything regarding coast phases optimization. I tried reading a few papers months back, but this level of math is far beyond me.

Link to comment
Share on other sites

@Reddy looking at your video - I don't think rotating frame is really much of an issue, at least for RSS. Kerbin's rotation is 4 times faster, and the extra apparent acceleration due to rotating frame is 0.6 m/s2 at LKO, very much noticeable.

As far as I understand it, below certain altitude KSP internally uses rotating frame with the good old centrifugal and Coriolis forces. In that frame, cartesian coordinates of launch site (or any fixed geoposition) remain fixed for any observer resting on ground. All constant vectors in inertial frame (of which I only know the Solar Prime vector), on the contrary, appear rotating. Now here's the trick: orbital velocity components exposed in kOS are the velocity in inertial frame projected onto instantaneous position of rotating axes. This makes orbital velocity at low altitudes inconsistent with position changes ("inertial force"). At high altitudes, vorbit = dR/dt, at low vsurface = dR/dt.

What is even more confusing is that at low altitudes it appears as though dvorbit / dt <> -body:mu*R / |R|3 on coasting trajectory. Again, this is because we can only compare this component-wise in kOS. And vx = (v, ix), so dvx / dt = (dv/dt, ix) + (v, dix / dt) =  gx + (v, [Omega x ix ]), where gx = -body:mu*Rx / |R|3 is the projection of gravity vector onto X axis and Omega is rotational speed. Overall, we get dv' / dt = g + [ v' x Omega ] where v' is the set of coordinates for vorbit in the rotating frame.

This appears as quite a drastic change if you log body:position and ship:velocity:orbit along orbit with Ap above critical altitude and Pe below it.

Now, about consequences of this. Both current state (r, v) and desired state (rd, vd) change in the same way, yes, but they are rotating, not translating. That means that the difference between them is not constant but also rotates.

To illustrate, let's assume you initialize r at launchpad with coordinates V(600, 0, 0) and at the end of loop you compute rd 100 km directly above it at V(700, 0, 0). Let's say a computation loop takes one second. In the inertial frame, you'd expect that V(600, 0, 0) does not point at launchpad anymore because the launchpad has rotated. However, V(600, 0, 0) still points exactly at launchpad in the rotating frame. Therefore, rd = V(700, 0, 0) in the new coordinate frame is again directly above KSC and not a bit westward as you probably want. The solution is to rotate all computed positions, velocities and directions a little bit every iteration.

Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

×
×
  • Create New...