Advanced Question Multi Stage PEG Algorithm...

[fred18]

Hi all,

I am developing a launcher family in dll, every is quite finished on the first two-stage rocket. I actually drive the first stage on a simple gravity turn and the second stage using the correct PEG algorithm that is doing just fine, very precise and nice.

In the next rockets of the family anyway I will have to develop the multistage algorithm. I have been trying on the two stage one for months, with hundreds of thousands of trials but up to now i did not succeded in having a multi stage algorithm working...

It seems that the estimated values at staging do not work properly, even if all the steering integrals and procedures are exactly as per the documents availables around (like orbiterwiki or the PEG addons already published on OH).
It is not a simple mistake of typos or anything silly because I double and triple and multiple checked everything thousands of times... so it looks like I need someone who is willing to help me, explaining the procedure correctly!

please help

 

[indy91]

I recently looked a little bit into this, too. At the Orbiterwiki page about the PEG algorithm I noticed a strange difference (mistake?) between the deltaV calculation of the single and multi stage guidance. In the single stage guidance the term v_e*T is added and in the multi stage guidance it is multiplicated to another term. I never got the algorithm to work, but also didn't look too much into this. I have no idea if this is helpful though.

 

[Enjo]

Hi Fred,

I know exactly what you feel, as I've been trying to implement MultiPEG two years ago, and since I got married in meantime, the effort is on hold for now until we sort our lives here. I remember very well the inability to tell anything about what's wrong whatsoever. I thought that I could draw some output data from what is seemingly working: the Java implementation, but couldn't make myself to do it.

Here's what I've been able to code so far, but as I said, it's not working and I don't have time to develop a methodological framework to figure out what's wrong:
http://sourceforge.net/p/enjomitchsorbit/code/HEAD/tree/launchmfd/PEG/PEGMulti.cpp
I know that woo wanted to work on this subject. I don't know how far we went.

Anyway, wishing you luck! If you ever succeed, please be polite enough to fulfill the GPL license of the Java implementation, that you are probably basing your work on, and share your code, so we can all benefit. :tiphat:

 

[Urwumpe]

I also plan working on a multistage PEG soon, but with the modification, that the first stage is supposed to be a coast phase. Not sure yet how to implement this coast properly.

 

[Enjo]

What is coast phase?

 

[Urwumpe]

What is coast phase?

No thrust, ballistic flight. Just predicting the ignition of the second stage to reach orbit.

 

[Donamy]

No thrust, ballistic flight. Just predicting the ignition of the second stage to reach orbit.

Also put's the 2nd stage into proper position.

 

[Urwumpe]

Also put's the 2nd stage into proper position.

Yes, that as well, that is what PEG also solves. Though it can for the start be simplified by just looking at the vertical speed and pointing horizontal.

 

[BruceJohnJennerLawso]

Also put's the 2nd stage into proper position.

And giving a bit of time to safely clear the last stage. It shouldn't be too long though, time is definitely a lot of money in this particular situation.

 

[fred18]

Thank you very much for your support guys!

At the moment the only solution that I can find is to take the original Nasa paper, take the assumptions and redo the calculation by myself... If I will find how to work this out I will share with all of you of course!

Another thing that I was thinking about but I still did not have time to try: single stage PEG works perfectly, and inside PEG the only stage specific parameter is Ve (or isp in orbiter because of the moltiplication by G). Now let's assume that I calculate my multistage rocket total dV, and then I imagine that my rocket is a SSTO with the last stage empty mass, the total fuel mass and I calculate a Ve that will grant the same dV as the multistage. Then I use this "average" Ve in the single stage PEG... I have to try this anyway, because i have no idea if it could work, I will let you know in case

 

[fred18]

Another important mistake on orbiter wiki:

The formula:

is not for h[T1], it is for delta H from instant now to instant T1!

so h[T1] will be equal to H[0]+delta H.

step by step we will see the end of it :yes:

 

[fred18]

All right gentlemen, as promised. I made it! :sweet::woohoo:

My multistage peg alghoritm works both with two and three stages rockets, and to me seems to do its job really good!

I am not a programmer, and the code is still a lot messy and not clear since I did a lot of trials and error with it so first of all I will post here the alghoritm I use, so you can implement it in your own way. For the code I will need a lot of refinery to be acceptable to you programming master guys :lol:

Let's start:

first of all, since the normal peg alghoritm worked perfectly, and I didn't want to loose its functionality I keep multistage and regular peg separate, so I use the multistage peg to all the stages that are not the last stage, and I then use the regular peg for the last stage and orbit insertion, simply initializing it with the outputs of the multistage peg.


I also have an important semplification that helped me a lot: I always point either to the apogee or to the perigee, never in the middle.

With my rocket the most comfortable boundary for choosing is around 260 km, so if the user inputs an apogee above 260 km the guidance will target the perigee. Of course if both target apogee and perigee are above 260 km it's ok, but targeting the perigee save us a bit of energy (and fuel).

then the procedure:

steering integrals: they cannot be kept defined as they are now, but they must become stage specific, so b0 must become function of T, isp and tau, bn=f(T,n,isp,tau); c0=f(T,isp,tau); cn=f(T,n,isp,tau)

so when you call them for a future stage you are sure that you are calculating those with the correct values. This is crucial, you miss one correct tau or isp here, and you're done, no multipeg...

of course when I am talking about ISP, i'm referring to Ve, but in orbiter ISP is already considered isp*g;


Initialization:

all variables should be initialized to promote the converging of the system, I have function that initialize all of them that is called whenever anything goes wrong so it resets the system

I introduce here two numbers that I will used later on that are important

a10 is the acceleration of the last drop of fuel of the current stage, so basically thrust of current stage / empty mass of current stage (inclusive full mass of all next stages and payload of course)

a20 initial acceleration of next stage: is simply the thrust of next stage / sum of all masses that will be on the rocket just immediatly after the separation


then usual alghoritm:


Multi-navigate:

velocity vector, radius vector,angular momentum vector, directions of all three of them, and direction of theta, as per normal alghoritm then angular speed and everything as per reguar peg except tau definition:
they shall be defined here as

tau of current stage=isp of current stage / current acceleration.

for acceleration I used the one along the velocity vector, the ACC value in surface mfd, there is the code to obtain it some where in the forum

tau of next stage=isp of next stage / a20


Multi-Estimate

here things are a bit more complicated:

first of all we have a series of new variables:

T1= time remaining to staging, i made a function that calculate this, basically is remaining burning time = current ISP * current Fuel mass / current thrust

then r[T], radius at cutoff is either radius earth + target apogee or radius earth + target perigee depending on the input (as said above)

the condition at staging:

r dot [T1] = current rdot + b0 (T1,isp[current stage],tau[current stage])*A+bn[T1,1,isp[current stage],tau[current stage])*B

r [T1] = (current r + r dot * T1) + c0(T1,isp[current stage], tau [current stage])*A + cn [T1,1,isp[current stage],tau[current stage])*B

then the usual procedure bu referred to T1, so

fr = A + (mu / (r*r) - omega*omega*r) / current acceleration
fr [T1] = A + B*T1 + ((mu / (r[T1] * r[T1]) - (omega[T1]*omega[T1]*r[T1])) / a10

frdot=(fr[T1]-fr)/T1
ftheta=1-(fr*fr)/2
f dot theta = -(fr*frdot)
f dot dot theta = -(frdot*frdot)/2

then one thing very important that I already stated above, the next formula was wrong in orbiter wiki and in some documents I found around!

in some documents now there was the definition of H[T1], but is actually delta_H! so

delta_H=(r+r[T])/2 * (ftheta*b0(T1,isp[current stage],tau[current stage])+f dot theta * bn(T1,1,isp[current stage], tau[current stage]) + f dot dot theta * bn(T1,2,isp[current stage], tau[current stage]))

then

h[T1] = current h + delta_H

vtheta[T1] = h[T1] / r[T1]
omega[T1]=vtheta[T1]/r[T1]

then for delta A and delta B I post the formulas of orbiter wiki since they are correct, take for a1[T1] and a2[0] what I called above a10 and a20, and of course ve,1 is isp[current stage] and ve,2 is isp[next stage]

then final conditions

h[T] = r[T] * final orbital velocity (calculated apart from target apogee or perigee)


next stage delta h = h[T] - h[T1]
next stage rbar = (r[T] + r[T1])/2
omega[T]= final velocity / r[T]

then I have here another semplifications, but actually works perfectly for me, so ...

next stage deltaV=next stage delta h / next stage rbar

then next stage burn time

T2 = tau[next stage] * (1-exp(next stage deltaV/isp[next stage])

this last passages are for two stages rockets, if i have three or more it changes in this way

I calculate apart deltaV for all the intermediate stages and burning time for all the intermediate stages

then

instat T2 (cutoff) is

T2=sum of burning time of intermediate stages + tau[last stage]*(1-exp(deltaV - sum of dVs of intermediate stages)/isp[last stage])


and now the final set

MultiGuide
to solve the matrix for obtaining A and B the coefficients are

a = b0(T1,isp[current stage],tau[current stage])+b0(T2,isp[next stage],tau[next stage]

b = bn(T1,1,isp[current stage],tau[current stage])+bn(T2,1,isp[next stage],tau[next stage])+b0(T2,isp[next stage],tau[next stage])*T1

c= c0(T1,isp[current stage],tau[current stage])+c0(T2,isp[next stage],tau[next stage]) + b0(T1,isp[current stage],tau[current stage])*T2

d= cn(T1,1,isp[current stage],tau[current stage])+bn(T1,1,isp[current stage],tau[current stage])*T2+c0(T2,isp[next stage], tau[next stage])*T1+cn(T2,1,isp[next stage],tau[next stage])

then the two solution terms:

y1= r dot T- rodt - b0(T2,isp[next stage],tau[next stage])*deltaA-bn(T2,1,isp[next stage],tau[next stage])*deltaB

y2=r[T]- r - rdot * (T1+T2) - c0(T2,isp[next stage],tau[next stage])*deltaA-cn[T2,1,isp[next stage],tau[next stage])*deltaB

anyway for this part orbiter wiki is perfect so here are the formulas, that can be cross checked with the "next stage, current stage" references of what i wrote above:

and then to solve of course

and then solve for A and B as usual and find

pitch = A + (mu/(r*r) - (omega*omega*r)) / current accelearation


that is most of it...

I am sure there will be some typo or something somewhere, so forgive me in advance.

I hope that this would be helpful to those who are willing to implement this difficult system not present anywhere around in orbiter

Cheers to every one!