Solving the Optimal Launch to Orbit Two Point Boundary Value Problem

[Arrowstar]

Hi everyone,

I have a question regarding the solution technique for the two point boundary value problem that arises when computing the optimal launch to orbit of a rocket. The problem I'm posing is for a planar rocket launch with no atmospheric drag, constant thrust over mass (F/m), and uniform gravity (flat Earth). I currently have the following differential equations (having properly applied the Euler-Lagrange theorem and adjusted to remove the scale invariance problem):

x_dot = Vx
y_dot = Vy
Vx_dot = (F/m) * cos(alpha)
Vy_dot = (F/m) * sin(alpha)
alpha_dot = -beta*cos(alpha)
beta_dot = beta^2 * sin(alpha)

Here, x is the position of the rocket parallel to the surface of the flat Earth, y is the altitude of the rocket, Vx is the horizontal speed, Vy is the vertical speed, alpha is the steering angle of the rocket as measured from the horizontal (and is also a co-state that arises from removing the scale invariance problem) and beta is just another co-state.

In addition, I have the following boundary conditions:

x(0)=Xo
y(0)=Yo
Vx(0)=Vxo
Vy(0)=Vyo

y(t=tf) = r
Vx(t=tf) = Vc
Vy(t=tf) = 0

The initial conditions are known but are represented symbolically in order to maintain generality. The final conditions correspond to a circular orbit at altitude r (orbit radius would then be Rearth + r), where the final speed is purely tangential and equal to the circular speed, Vc, of the corresponding orbital radius. The final time computed is tf.

Now, here's my problem. I'd like to use a simple shooting method to solve this problem. To do so, I'll make some guesses for the initial state of alpha and beta and propagate forward in time. This will give me some final states of all my state variables. I can then create a residual vector, F, that is the difference between the propagated final states and the desired final states:

F1 = y(tf)-r
F2 = Vx(tf)-Vc
F3 = Vy(tf)-0

Here's where it gets weird. In the typical shooting method, you would compute the residual for a state variable and then use, say, secant method to compute a new initial guess for that same state variable. However, here, I have the interesting problem of having all of my final boundary conditions on variables which have defined initial conditions as well.

What I need, and what I'm hoping some folks here can provide insight into, is how I take the residual vector, F,and use that to compute new initial guesses for alpha and beta, the only two variables I have that don't have initial conditions.

Any help would be appreciated. Thanks!

 

[Notebook]

Intersting stuff Arrowstar, do you have any links to references I could look up?

Thanks, N.

 

[Jarod]

Thanks for working on that kind of stuff, well I mean if this is for Orbiter
I didn't work with that kind of maths for way too many years to be helpful, but I was thinking about autopilot for rockets and are you sure of your final boundary conditions ?
I mean that was my first thought too, and then I thought a two burn insertion would be more fuel efficient in some cases than a single one.
Or at least having inequalities, because if a solution more fuel efficient exists, right now you'll just fly through a less efficient trajectory just to end up at the apogee at your desired speed in one burn.
But this is still kind of blurry in my head so maybe I don't know enough.

edit: oups forget about some of this, I was thinking your tf variable was for maximum burn time possible.

 

[martins]

First of all, you can't really optimise your parameters individually, but you need to define an objective function that combines all data, for example

F(alpha,beta) = a dr^2 + b dVx^2 + c dVy^2

where dr, dVx and dVy are the errors in your model data, and a, b, c are suitable scaling parameters.

Then you can define a gradient vector that defines the rate of change of F as a function of the free model parameters:

g = [dF/dalpha, dF/dbeta]

Unfortunately, since the problem is fairly nonlinear, and since your forward model is Orbiter itself, there isn't an analytic way to get g, so you will have to calculate it numerically, by perturbing first alpha, then beta, and see what the result on F is:

dF/dalpha  ~ [F(alpha+dalpha,beta)-F(alpha,beta)]/dalpha
dF/dbeta ~ [F(alpha,beta+dbeta)-F(alpha,beta)]/dbeta

This means that you need to run 3 forward model runs to get g.

Once you have g, you can use a gradient descent method to get an update on alpha and beta:

[alpha,beta]_{n+1} = [alpha,beta]_n - gamma g

where gamma is another scaling value. Note that since the problem is nonlinear, the gradient is only correct in the vicinity of [alpha,beta]_n, so you cant make arbitrarily large steps (gamma needs to be "sufficiently" small to guarantee a downhill step). This could be implemented by a 1-dimensional line search, but that would be rather costly here, since it requires multiple evaluations of the forward model. Note that you have to recalculate g for each iteration.

There are better methods than steepest descent. Look for example for nonlinear conjugate gradients, or Levenberg-Marquardt.

NB: Since time t is also a free parameter in your model, you may have to add that to your parameter list, but maybe you can avoid that by simply running each forward model until it hits a minimum of F, from which you can extract an "optimum" t for the given [alpha,beta] set. This assumes that you don't hit local minima along the way, which may or may not be justified.

 

[Notebook]

What folk above said....

N.

 

[Arrowstar]

Thank you, Dr. Schweiger! I ended up creating another state for the final time and then using a normalized time, tau, to integrate along (where 0<tau<1). My equations of motion then have the following addition: t_dot = 0 (where the dot now signifies derivatives w.r.t. tau). Obviously since tau and t increase at the same rate, the derivative is 0. You can close out the equations of motion by then scaling the right-hand side of each equation by t.

I did end up using a root-finding technique very similar to what you suggested and it worked nicely, if a bit slow. I will look into the alternatives you noted, however.

Thank you, sir!

 

[Jarod]

Found these documents about a different way to approach the problem :
Rapid Ascent Trajectory Planning and Closed-Loop Guidance for Responsive Launch
http://www.responsivespace.com/Papers/RS7/SESSIONS/Session I/1004_Chavez/1004C.pdf
http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20030000638_2002153466.pdf
http://unica.academia.edu/MehranZar...guidance_law_for_an_orbital_injection_problem