Jump to content

PEGAS - Powered Explicit Guidance Ascent System - devlog


Reddy

Recommended Posts

So the problem is not that bad.  The typical values vessel.obt_velocity, vessel.body.position, and vessel.CoM and such are in normal inertal non-rotating worldspace.  As long as you translate from worldspace to body-centered coordinates it all works out and can be run through PEG.

Where it gets weird seems to be isolated to the Orbit class.  In addition to getting xzy "swizzled" the velocity values there need to be offset by the frame rotation.  I think the issue here is that internally KSP is doing things in the rotating reference frame and it is handing those values off to the Orbit class and expecting to get them back in the same rotating frame.  But most of the values that are exposed on Vessel and CelestialBody are fine.

Where I ran into this was that I cheated with my CSE routine by using the Orbit class in KSP directly (which you can't do in KoS at all).  Swapping in the CSESimple routine fixed the 145km wobble in my ascents.  Where I'm winding up now looks something like this:

        private void CSEKSP(Vector3d r0, Vector3d v0, double t, out Vector3d rf, out Vector3d vf)
        {
            Vector3d rot;

            if ( vessel.altitude < mainBody.inverseRotThresholdAltitude )
                rot = mainBody.getRFrmVel(r0 + mainBody.position);
            else
                rot = new Vector3d(0, 0, 0);

            CSEorbit.UpdateFromStateVectors(r0.xzy, (v0 + rot).xzy, mainBody, vesselState.time);
            CSEorbit.GetOrbitalStateVectorsAtUT(vesselState.time + t, out rf, out vf);

            if ( vessel.altitude < mainBody.inverseRotThresholdAltitude )
                rot = mainBody.getRFrmVel(rf + mainBody.position);
            else
                rot = new Vector3d(0, 0, 0);

            rf = rf.xzy;
            vf = vf.xzy - rot;
        }

I've still got a little bit of a wiggle now, but its of the order of 5 degrees and not 30.

The only person I know of that really knows this stuff though is @taniwha

EDIT: so what i suspect is missing here is that the Earth needs to rotate around the ship and probably rf and mainBody.position are off.

Edited by Jim DiGriz
Link to comment
Share on other sites

Ha! Looks like my "corrections" increased the rotation instead of compensating it.

This is what happens when you adapt formulas written in right-handed inertial frame for left-handed rotating one. (Alleged violation of 2.2g)

I'll just go brute-force and recompute coordinates like this:

function chFrame {
  parameter Vec, oldIx, newIx.
//changes vector from frame with ix = oldIx, iy = V(0,1,0) to frame with ix = newIx, iy = V(0,1,0)
//this works under the assumption that all planetary axes are V(0,1,0) in KSP
//Assumed that oldIx and newIx are Solar Prime Vectors at different moments of time
  local oldIz to V(-oldIx:z, 0, oldIx:x).
  local newIz to V(-newIx:z, 0, newIx:x). 
  return vdot(Vec, oldIx)*newIx + vdot(Vec, oldIz)*newIz + Vec:y * V(0,1,0).
}

 

Link to comment
Share on other sites

@Jim DiGriz I clearly remember that even without state prediction, dVorbit/dt was not mu/R^2. The offset was exactly what's expected from the rotating coordinate axes.

I'm trying to rewrite PEG for maneuver execution now. Instead of burning the specified deltaV, it takes the trajectory after node as the reference and tries to put the ship onto it. The unbound orbital element in this case is the anomaly, i.e. where on that trajectory ship completes the burn. Compare with the unbound AoP for ascent, i.e. how the SMA is aligned relative to DN-AN  line.

So far the result does not seem as much improvement over burning specified dV vector, but the main goal is to test my implementation of the core functions.

@Reddy can you explain your logic of control loop convergence? Do you iterate until both Vgo and tgo converged and then update the steering command? I saw in the russian source that new steering vector can be issued right after tgo has converged and Vgo convergence is the condition to start next loop. But that produces side effects when tgo "converges" to the total fuel burnout time. On the other hand, setting convergence time threshold as absolute value (say, 0.1s) leads to too many iterations initially, while the actual steering direction does not change much between them.

On the fun side: tried to understand why my numerical integration-based CSE is off by a few km on 30-second timescale. Realized I just skipped the last integration step.

Link to comment
Share on other sites

7 hours ago, Pand5461 said:

@Jim DiGriz I clearly remember that even without state prediction, dVorbit/dt was not mu/R^2. The offset was exactly what's expected from the rotating coordinate axes.

Well, I've got a very good replication case of the inverse rotation issue.  I've also debugged nearly all the relevant vectors that I'm using and the only wonkiness I see is isolated to the Orbit class.  The only place I was using that was in the CSE routine.  And this makes sense because if the vessel position and velocity or mainbody position, etc were wonky then there would be substantially more issues around players not being able write plane autopilots, hoverslam scripts and other endoatmospheric routines.

Link to comment
Share on other sites

When it comes to node execution there's a characteristic time of tlambda = J/L -- the "burn centroid time" according to 19760024151 -- which i suspect is (close?) to what you should lead the node by to start the burn.  I haven't gotten there yet though.

Link to comment
Share on other sites

@Jim DiGriz Setting the lead time to J/L actually improved the convergence a lot. Now ship does not try somersaulting from the very beginning. Still wants to misbehave by the end of the burn, but I just restrict magnitude of (lambda dot) * J/L.

The accuracy of apses for a stock 6-minute 850 m/s prograde burn is now Ap 10320 km (target 10217), Pe 94 km (target 86.75). SMA is rotated by a few degrees unfortunately. I'm trying to find out whether I need to restrict steering a bit more or allow more aggressive steering.

UPD: restricted magnitude of (lambda dot) * J/L to 0.15, that improved accuracy of maneuver. Test case: cheat craft to 686750m orbit (min available from cheat menu), maneuver node 850m/s with equal NRM and PRG components. Initial TWR 0.23, "nominal" burn time 5m35s. Result: plane difference 0.15 degrees, AoP difference about 0.5 degrees, SMA difference about 0.5% (didn't check eccentricity). Ap accuracy is about 50 km, much better than burning straight along node prograde. However, the maneuver took extra 65-70 m/s.

I've noticed also a small error in kOS update procedure. dVsensed is the velocity gain measured by accelerometers, so it does not include effect of gravity IRL. So, setting dVsensed = Vcurrent - Vprevious introduces a bias to the initial guess of Vgo at the start of PEG loop. This gets corrected after one or two iterations of course but convergence is faster if gravity effects are subtracted at loop update routine.

Edited by Pand5461
Link to comment
Share on other sites

@Reddy is this a bug in your CSE routine?  https://github.com/Noiredd/PEGAS-MATLAB/blob/master/MATLAB/CSEroutine.m#L98-L99

It seems like the if should be dropped, and the if prevents the while loop from ever executing...  (EDIT: and the conditional in the while loop needs reversing)

Also in USS the loop is supposed to always execute once, but the k = k + 2 line needs to happen after the check, so I started with your code:

https://github.com/Noiredd/PEGAS-MATLAB/blob/master/MATLAB/CSEroutine.m#L185-L192

But mine wound up looking like this in c#:

        static double USS(double xarg, double a, int kmax)
        {
            double du1 = xarg/4;
            double u1 = du1;
            double f7 = -a*du1*du1;
            double k=3;
            while( true )
            {
                du1 = f7*du1 / (k*(k-1));
                double u1old = u1;
                u1 = u1+du1;
                if ( u1 == u1old )
                    break;
                if ( k >= kmax )
                    break;
                k=k+2;
            }
            return k;
        }

 

Edited by Jim DiGriz
Link to comment
Share on other sites

8 hours ago, Pand5461 said:

@Jim DiGriz Setting the lead time to J/L actually improved the convergence a lot. Now ship does not try somersaulting from the very beginning. Still wants to misbehave by the end of the burn, but I just restrict magnitude of (lambda dot) * J/L.

Excellent.  I'm still working on getting the CSE routine working fast + in the atmosphere and getting ascents working.

Link to comment
Share on other sites

10 hours ago, Pand5461 said:

@Jim DiGriz Still wants to misbehave by the end of the burn

You need to freeze guidance.

lambda and lambdaDot should be frozen and you should stop iterating PEG something like 10 seconds before the end of the burn.

Then you can do something like this:

dt = vesselState.time - last_PEG

iF = ( lambda + lambdaDot * ( dt - J/L ) ).normalized

You can also use that to reduce the frequency you run PEG to reconverge to 1 Hz.  You can also update tgo and vgo, and can cut the burn when tgo < 0, but practically I find its better to use something like the orbital angular momentum to cut the burn, since bugs in staging/dV/thrust/isp calculations can throw that off.

Dunno what to do yet if the initial time of the burn is less than 10 seconds...

Link to comment
Share on other sites

18 hours ago, Jim DiGriz said:

You need to freeze guidance.

lambda and lambdaDot should be frozen and you should stop iterating PEG something like 10 seconds before the end of the burn.

That's what I did. Thing is, lambdaDot goes crazy way before 10 seconds to finish. At tgo = 120 sec it has magnitude of over unity, so I need to scale it down. The constraint for rd being on a very specific trajectory is too harsh, I guess. Here's about as close as I can get:

ax3nS0x.png

18 hours ago, Jim DiGriz said:

Dunno what to do yet if the initial time of the burn is less than 10 seconds...

Just burn maneuver prograde / Vgo vector, it must be just fine.

Link to comment
Share on other sites

I hate to barge in on all this heady math speak, but I've got a bit of a question. Well, two actually, which are somewhat related.

 

Okay. So First is, I know that I have to have inclination and LAN specified to make things work. Great. I know what inclination is, that's easy. I have a general understanding of the concept of LAN once in orbit, but that as a function of when to launch I'm a bit confused about. In short, sitting stationary on the ground, I have no idea where I am in relation to a given LAN. So if I want to launch from the Cape at say 0500Z, how do I translate that into a LAN that will launch into that orbit?

 

Or, how do I tell PEGAS to launch into an inclination upon activation of the script? In other words, "I don't care what my LAN is going to be, let's light this candle."

Link to comment
Share on other sites

21 minutes ago, gsarducci said:

Or, how do I tell PEGAS to launch into an inclination upon activation of the script? In other words, "I don't care what my LAN is going to be, let's light this candle."

That's easier. Let's say you have the inclination i. Then, the trajectory crosses the meridian at the angle xi = arcsin( cos(i)/cos(latitude) ) Then, if you are instantaneously at that orbit, its "virtual" LAN is given by the following function:

function azlan {
  parameter xi.
  
  local Vvirt to heading(xi, 0):vector.
  local nvirt to vcrs(body:position, Vvirt).
  local ANvirt to vcrs(nvirt, V(0,1,0)).
  local vlan to arctan2( vdot( V(0,1,0), vcrs(ANvirt, solarprimevector) ), vdot(ANvirt, solarprimevector) ).
  if vlan<0 set vlan to vlan+360.
  return vlan.
}

As @Reddy mentioned, he planned launches 8 minutes before orbit passes over launch site, which translates into Earth rotation by 2 degrees. Therefore, adding 2 degrees (maybe 2.5, to be on the safe side) to the calculated value is more or less what you want.

EDIT: if you want to launch at other time than "now", just calculate angle of rotation as 360*(ETA to launch)/body:rotationperiod and add that to previous calculation.

Edited by Pand5461
Link to comment
Share on other sites

I think I missed something here. How can I see the value of "vlan" or practically apply it to my solution? Running that program that you gave to me doesn't display a solution, it just simply runs and ends. I assume the value went somewhere, just not into my brain. :) Is this something intended to drop into my mission.ks script in lieu of an actual LAN entry?

 

EDIT: Also, when applying your equation to solve for value xi, how then do I apply that solution to the program cited? Simply replacing xi with an integer value breaks the program. My assumption is that the program takes the solution for xi and outputs the LAN for an instantaneous launch, thus I would need to add some "slop" to it considering rate of rotation for a launch time using the equation you put in your edit?

 

Thanks for your patience with this Pand5461. I am learning a lot with this!

Edited by gsarducci
Link to comment
Share on other sites

This should print vlan out:

set xi to arcsin( cos(incl)/cos(latitude) ).
set vlan to azlan(xi).
set launch_now_lan to vlan + 2.5.
print launch_now_lan.

Inclination must be higher than latitude of course and azlan function must be loaded to work of course.

Link to comment
Share on other sites

Excellent. Thanks for that. I think I can see how this all fits together now. Hopefully I can give it a shot tonight when I get home.

 

Since I've got your ear I have two unrelated questions. First, I am curious if there is a way to impart a spin on a stage prior to stage/ignition. I am launching a Thor-Delta with a Telstar type satellite. The third solid fueled stage is spin stabilized but I don't see a way to command the program to "spin up" the second stage prior to release.

Second is a roll program. Given that the spacecraft is pointing directly up, I believe the best solution for a roll while heading vertical would be a directional roll without regard to heading, so say roll left 24 degrees and stabilize, followed by a specific roll lock of "0" (or whatever is appropriate) once pitchover is initiated and a horizon is established.

 

Both of these questions aren't so much "how do I do it", but more "how or can I do it within the confines of the PEGAS routine?"

 

Also, I'm gonna need to know where to send your bottle of 14 year old scotch when you're done educating me!

Edited by gsarducci
Link to comment
Share on other sites

Challenge:  can anyone get this craft file into orbit with PEGAS?

https://gist.github.com/cb281f07dd7b4d4da0ac970d69ed13e7

You'll need RSS/RO and RealScaleBoosters.  That is the RSB Atlas V HLV with enough payload on it to push the SLT of the first stage down to 1.05 so that it only crawls off the launchpad.  I get PEG convergence, but it never makes orbit, once it gets to the upper stage it needs a pitch of 40 degrees and the vgo starts continually increasing and it falls back into the atmosphere.  Wondering if its my implementation of PEG, or if I've simply built a rocket that will never make orbit, or if it takes a particular launch profile to make it, or if the PEG algorithm itself has limitations.  Squinting at it, perhaps the 1.05 SLT costs so much in gravity losses that there's not enough dV in the higher TWR stages to get it 'lofted' enough for the 0.60 TWR upper stage to ever make orbit.

And it should have a 1.13 starting vac TWR, 1.05 starting SLT 2814 dV first stage, a 0.98 starting TWR 2261 dV "sustainer", and 0.61 starting TWR 8763 dV upper stage for 13838 vac DV total.

(and i understand this is a faintly ridiculous rocket, but i'm trying to figure out how to determine more analytically why it is ridiculous...)

Edited by Jim DiGriz
Link to comment
Share on other sites

On 8/26/2017 at 2:03 PM, Jim DiGriz said:

But mine wound up looking like this in c#:


        static double USS(double xarg, double a, int kmax)
        {
            double du1 = xarg/4;
            double u1 = du1;
            double f7 = -a*du1*du1;
            double k=3;
            while( true )
            {
                du1 = f7*du1 / (k*(k-1));
                double u1old = u1;
                u1 = u1+du1;
                if ( u1 == u1old )
                    break;
                if ( k >= kmax )
                    break;
                k=k+2;
            }
            return k;
        }

I'm guessing this was just in the forum version, but I figured I'd point out that you probably want to return u1, rather than k.

 

Link to comment
Share on other sites

16 hours ago, Jim DiGriz said:

Challenge:  can anyone get this craft file into orbit with PEGAS?

https://gist.github.com/cb281f07dd7b4d4da0ac970d69ed13e7

You'll need RSS/RO and RealScaleBoosters.  That is the RSB Atlas V HLV with enough payload on it to push the SLT of the first stage down to 1.05 so that it only crawls off the launchpad.  I get PEG convergence, but it never makes orbit, once it gets to the upper stage it needs a pitch of 40 degrees and the vgo starts continually increasing and it falls back into the atmosphere.  Wondering if its my implementation of PEG, or if I've simply built a rocket that will never make orbit, or if it takes a particular launch profile to make it, or if the PEG algorithm itself has limitations.  Squinting at it, perhaps the 1.05 SLT costs so much in gravity losses that there's not enough dV in the higher TWR stages to get it 'lofted' enough for the 0.60 TWR upper stage to ever make orbit.

And it should have a 1.13 starting vac TWR, 1.05 starting SLT 2814 dV first stage, a 0.98 starting TWR 2261 dV "sustainer", and 0.61 starting TWR 8763 dV upper stage for 13838 vac DV total.

(and i understand this is a faintly ridiculous rocket, but i'm trying to figure out how to determine more analytically why it is ridiculous...)

Not going to happen, not on that rocket and you've already touched on the reasons why. Your payload of 170 tons means that you're only going to get about 3000m/s DV  before the Centaur stage starts up and whether it goes shallow or steep, most of that is eaten up by gravity losses. The Atlas V Heavy only had a hypothetical payload capacity of 30 tons. 170 is just not going to happen, not even to LEO

Edited by Starwaster
Link to comment
Share on other sites

12 hours ago, ubik2 said:

I'm guessing this was just in the forum version, but I figured I'd point out that you probably want to return u1, rather than k.

 

nope, it wasn't, although i figured that out as well.

11 hours ago, Starwaster said:

Not going to happen, not on that rocket and you've already touched on the reasons why. Your payload of 170 tons means that you're only going to get about 3000m/s DV  before the Centaur stage starts up and whether it goes shallow or steep, most of that is eaten up by gravity losses. The Atlas V Heavy only had a hypothetical payload capacity of 30 tons. 170 is just not going to happen, not even to LEO

Yeah, although if I back it off to 22t to LEO it wigs out completely as well so its both that i've got bugs and that rocket never would have made it to orbit.

Link to comment
Share on other sites

On 8/27/2017 at 10:43 AM, Pand5461 said:

That's what I did. Thing is, lambdaDot goes crazy way before 10 seconds to finish. At tgo = 120 sec it has magnitude of over unity, so I need to scale it down. The constraint for rd being on a very specific trajectory is too harsh, I guess. Here's about as close as I can get:

I think I'm seeing something like this on ascents as well.  With that Atlas V HLV with the payload tweaked to need a few hundred dV out of the upper centaur stage it has all kinds of difficulties, and behaves similarly to what you describe with lambdaDot freaking out, and the craft starts to spin.  Trying to constrain lambdaDot doesn't seem to help either.  I'm wondering if the PEG predictor needs to be replaced by a proper 4th order RK integrator with proper transversality constraints between the stages.

Link to comment
Share on other sites

Okay, I think I've got PEGAS more or less working in MechJeb now, although, its a mix of PEGAS and PEG-4 code (PEGAS steering and corrector, with PEG-4 thrust integrals and rthrust/vthrust integration).

I'm still seeing odd effects in the case of having burnout use a bit of a low-TWR "TLI" stage after a high-TWR upper stage.  I suspect the problem is the rthrust/vthrust expansion and that higher order terms need to be kept.

There's an updated approach to PEG in these references which has some good background:

https://arc.aiaa.org/doi/10.2514/3.25350

http://www.sciencedirect.com/science/article/pii/0094576589900969

Most of the launch I'm seeing the tangent angle at 45 degrees which produces a lot of errors and as a result the arc is too high and it wants to burn to too high of a periapsis and then pitch down and burn at high thrust, and the error created by the high TWR booster produces a bit crazy need to burn straight down in the low-TWR centaur to make up the difference.  And at 45 degrees of tangent angle the error estimate in the pitch angle is about 30% according to those authors, with a lot of improvement from keeping the 2nd or 3rd order terms in the expansion.  I haven't validated if PEGAS's kOS or MATLAB code sees this same behavior.

Also the best paper I've found so far on PEG's mathematical background is:

https://arc.aiaa.org/doi/abs/10.2514/6.1977-1051

That covers a lot of ground that the Journal of the Astronautical Sciences paper glosses over a bit (https://ntrs.nasa.gov/search.jsp?R=19790048206)

Link to comment
Share on other sites

  • 1 month later...

I am sorry I was absent for a while. As some of you might know, I am a PhD student - there are periods where I simply don't have time to pay attention to PEGAS and this thread.

I've just pushed some updates to the main repository, here's the list of new features:

  • roll control - you can define initial roll (vehicle will roll to that angle during the initial pitchover maneuver) via controls["initialRoll"], and then control the vehicle's attitude using events,
  • better throttle handling - RO does not treat the throttle slider in an absolute way: 50% slider does not correspond to 50% throttle, but 50% of throttle range that an engine allows. PEGAS can now handle this (fixing issues with g-limited stages), at the expense of users having to input an additional key ("minThrottle") along the "gLim" in their stage definitions,
  • new event type: jettison - allows user to inform the system of the amount of mass lost in a staging event, which is also important in acceleration-limited stages (eg. if a stage's dry mass included fairings, and then these were jettisoned, the vehicle will weigh less than PEGAS planned for - as a consequence, throttle will be lowered at a slower than expected rate, consuming fuel slower, causing the engine shutdown later than PEGAS scheduled the next stage ignition: as an effect, the constant-acceleration stage could smash into the next stage),
  • several minor tweaks: PEGAS will change the control part to the currently running terminal on boot, lack of "message" key in user event will no longer cause it to crash, etc.

We're almost halfway to the 1.1 release, which will include also the following features:

  • launch with no LAN constraint (something like a "launch now" function - LAN will still be constrained, but it will be hidden from user and dynamically calculated),
  • launch opportunities (ability to pick northerly or southerly launch, together with the option "nearest from now"),
  • better atmospheric ascent control (I'm thinking of something akin to aim points for optional, additional control over pitch during early ascent).

I'll read up on that thread later and try to reply you guys :)

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...