#newton's method

After some research, I found this dissertation by Hongsheng Sun which dealt pretty much with the problem at hand:

Closed-Loop Endo-Atmospheric Ascent Guidance by Hongsheng Sun

And after a failed attempt at my own approach, I tried to replicate his algorithm. Understanding the main parts of the algorithm used in this dissertation is the focus of this post.

The basic idea is to define initial and final conditions and then to get from one to the other in a smooth transition. We are not just talking about the smoothness of the flight path itself, but a smooth transition of all coordinates, including velocity.

Therefore, at any given time, our current set of coordinates has to be equal to our coordinates just slightly earlier, plus a slight perturbation:

$$y(t) = y(t-dt) + dt * \dot{y}(\frac{y(t-dt)+y(t)} {2}) $$

Where $\dot{y}$ is the time derivative of our coordinates and $dt$ is the time span we look into the past. For arbitrarily small $dt$ this is an exact solution. Since computers don’t work that well with infinities, we have to choose $dt$ small enough.

Since we can construct another set of coordinates based on the previous set, we can construct a solution which evolves smoothly from one point to the next, Dividing our smooth transition in many discrete time steps $y_i$.

If we now assume we have a path, connecting initial coordinates with final coordinates, we can determine how smooth such a transition is and apply some changes to make it even smoother. If that’s still not smooth enough, we can run the whole algorithm again and make the solution smoother still and keep repeating, until we are satisfied.

Calculating the changes we have to apply can be done by utilising Newton’s Method. Newton’s Method is a neat little trick to solve for zero:

$$x: f(x) = 0$$

Basically, Newtons method allows to get very close to any $x$ which would solve the problem. All we have to do now is to rewrite our expression so a smooth solution results in a zero:

$$0 {\overset{!}{=}} y_{i} - y_{i-1} - dt * f(\frac{y_{i-1}+y_{i}}{2}) = E_{i}$$

Where $E_{i}$ is simply the current actual value for this expression. So most probably, it’s not zero.

Now, through several iterations, we try to bring it as close to zero as we can! We do that by using Newton’s Method for vectors. In general, the next correction is given by

$$ \frac{\partial E}{\partial y} * \Delta y= -E(y)$$

Where $ \frac{\partial E}{\partial y} $ is the Jacobian of $E$ for all time steps of $y$.

Through this matrix multiplication, we get a set of linear equations we can solve for $ \Delta y $. Once $\Delta y$ is obtained, it can be added to our current $y$, improving our solution. Then we can do the same iteration again and again, until we have a $\Delta y$ sufficiently small. So thanks to our initial guess, we can optimise all time steps in just one operation! (In practice, we will find all solutions depending on the solutions of the previous time step.)

Now for the difficult part: Obtaining a time derivative and implementing initial and final conditions.

The time derivative of position is of course velocity, so that should be no problem. However, the time derivative of velocity is acceleration and therefore requires us to know the mass of our rocket. That’s difficult, because the rockets mass changes as it looses fuel during the ascent. So to know the mass at any given time, we have to first calculate the thrust we used previously.

To calculate the thrust we used ideally, we have to know by how much the velocity changed in previous time steps and then calculate the sum of all other forces acting on our rocket at the time. Through some mathematical juggling, the mass of time $t$ can be found and based on that any acceleration due to forces can be calculated. The thrust we used to calculate our weight is now actually the perfectly ideal thrust and should not be confused with the thrust we set and still have to optimise!

However, we want to fly most energy efficient. To get an energy efficient solution, we require some side conditions. I have to read up again on why these conditions are obtainable from the Hamiltonian, but that’s essentially what the dissertation does. It defines two side conditions, one for position and one for velocity, and defines them as partial derivative of the Hamiltonian. These side conditions have to evolve smoothly too, if we want an energy efficient solution.

Initial and final condition have to add up to the number of coordinates we want to optimise. If we use less, then we will find a solution which chooses one coordinate at random. If we use more, we will probably not find a solution at all, because there is none. (It would be equivalent to requiring our rocket to be in two places at once, just with nicer maths.)

Ideally, we specify half of the conditions as initial conditions, and half of the conditions as final conditions. That way we avoid strange solutions we didn’t expect and are nonsensical. For example, specifying only initial conditions would be equal to a simulation of a rocket left to itself, so most probably a falling rocket. Specifying only final conditions on the other hand would tell us to start with initial conditions to use as little energy as possible to get there.

Inital and final conditions again have to be implemented in a way that they can be “solved” with Newton’s Method. So instead of specifying our initial position as

$$ X(t=0) = X_init $$

we denote

$$ X(t=0) - X_init = 0 $$

Both statements are equivalent, but the second can be implemented seamlessly into our optimisation algorithm. Same goes for the final conditions. However, our velocity should allow us to reach a certain height, where we can then circularise to travel in a smooth circle in LEO. Therefore, what we do is to define an elliptic orbit according to Kepler’s laws with an apoapsis at the height we want to reach. From that condition, we can obtain:

$$ (GM(\frac{m}{L})^2 - \frac{1}{R_{Apo}})^2 = (GM(\frac{m}{L})^2-\frac{1}{R_{atmos}})^2 + (v_r \frac{m}{L})^2 $$

Where $R_{Apo}$ is the radius of our target orbit and $R_{atmos}$ is the radius at which the atmosphere runs out. (In Kerbal Space Program, that’s the specific height of 70km. Reality may differ.) The ellipsis condition can be transformed to:

$$ v_r^2 * \frac{R_{Apo}*R_{atmos}}{R_{Apo}-R_{atmos}} + v_\omega^2 \frac{R_{atmos}}{R_{Apo}}*(R_{Apo}+R_{atmos})-2GM=0$$

(Further final conditions are mentioned in the dissertation.)

So, that’s essentially it. We have initial and final conditions, have a function pointing from each position to the next one, and try to get from one point to the next with as little effort as possible. There are some tricks to determining the quality of the solution and simplifying the Jacobian, but they are beside the point. I hope you gained some insight into the topic.