I need someone help me do some math for launch optimization

[alterbaron]

AFAIK, the error in this method will be on-par with regular velocity-verlet.

There's a document which mentions this method here.

Still looking into PSOPT stuff, but mostly working on reading up on the math.

EDIT: A simpler and more practical problem would be to optimize on the parameters given to the mechjeb ascent guidance system. Results from this sort of investigation would be immediately useful to players as well.

(Though the problem itself isn't as cool . . .)

Edited August 23, 2013 by alterbaron

[tavert]

Edited away talk of Lagrange polynomials and LGL points? Was getting interesting.

I just remembered I think I still have a copy of GPOPS-II (http://vdol.mae.ufl.edu/SubmittedJournalPublications/TOMS-GPOPS-II-February-2013.pdf) somewhere, since one of the authors was asking me for help a few months ago with the Matlab interface to Ipopt. That might be easier to use than PSOPT since it's in Matlab, though it's much more restrictively licensed (only free if you work for the federal government or the University of Florida...).

[alterbaron]

Just a small update on this stuff:

PSOPT is really hard to install on windows, so instead, I threw together a little pseudospectral optimal control solver in MATLAB (using fmincon from the optimization toolbox). I've tested it on the brachistochrone problem with good results: http://imgur.com/LjYLYdM.

I should be able to modify the code to solve for an optimal ascent trajectory (fingers crossed).

[tavert]

Oh dear. Good luck, but I anticipate this problem to be more numerically challenging. ICLOCS maybe? http://www.ee.ic.ac.uk/ICLOCS/

I finished writing a direct RK4 optimization formulation, but wasn't providing a good enough initial guess to get reasonable convergence.

Edit: I also spent too much time today reading and figuring out this paper that closette also linked to way back in the old max-altitude challenge: http://arxiv.org/pdf/math/0703911v1.pdf

Revisited Pontryagin principle for the first time in a while, I now have a Mathematica notebook where I worked out the switching function and singular-arc throttle for KSP's drag and atmospheric-Isp model (in 2D, assuming equatorial launch). Unfortunately NDSolve didn't want to solve the boundary value problem with a piecewise-defined throttle function, and the symbolic solutions are way too long and messy to try porting them to some other language.

Interestingly, in 2 or 3 dimensions, the optimal throttle direction isn't exactly aligned with velocity, rather it's aligned with the velocity costate (adjoint) vector. The costates have slightly different time evolution, but I imagine it shouldn't be too far away from the velocity direction for decent TWR.

Edited August 28, 2013 by tavert

[alterbaron]

Thanks for all the useful info! ICLOCS looks like a nice solution -- I'll give it a shot.

I've been able to get my homebrew solver to generate pretty accurate solutions to the (purely vertical) Goddard maximum ascent problem as defined here. (My Result)

It's probably too slow to solve more complicated problems, but it was a real learning experience to write.

In any case, I have a problem formulation (for a single-stage rocket) ready to plug into whichever solver I can get to work first:

EDIT: Problem Formulation

Edited August 28, 2013 by alterbaron

[K^2]

Have you checked your solution against the analytic solutions for the simplified case? In other words, with simple exponentially dropping atmosphere, do you get optimal ascent at terminal velocity for most of the trip?

[ThreeKelvin]

EDIT: Problem Formulation

Looks nice. You might want some further constraints:

T(t_f) = 0 (No thrusting at the last instant in order to cheat the constraint \dot{r}(t_f) = 0)

r(t) >= r_{planet} (No shortcut through the planet - It might not be needed though.)

[alterbaron]

Have you checked your solution against the analytic solutions for the simplified case? In other words, with simple exponentially dropping atmosphere, do you get optimal ascent at terminal velocity for most of the trip?

I ran a test using KSP's dynamics and a test craft with a high TWR.

The objective was to maximize the final altitude for a vertical ascent.

Using my solver, this is what I got:

Velocity as a Function of Altitude

Ascent Overview

Using a better solver should allow these results to be refined much further.

Looks nice. You might want some further constraints:

T(t_f) = 0 (No thrusting at the last instant in order to cheat the constraint \dot{r}(t_f) = 0)

r(t) >= r_{planet} (No shortcut through the planet - It might not be needed though.)

Nice catch!

[tavert]

Thanks for all the useful info! ICLOCS looks like a nice solution -- I'll give it a shot.

I've been able to get my homebrew solver to generate pretty accurate solutions to the (purely vertical) Goddard maximum ascent problem as defined here. (My Result)

It's probably too slow to solve more complicated problems, but it was a real learning experience to write.

In any case, I have a problem formulation (for a single-stage rocket) ready to plug into whichever solver I can get to work first:

EDIT: Problem Formulation

Yeah, writing nontrivial stuff like this is a good way to learn, but a bad way to get something that runs quickly and reliably enough for challenging problems.

Your formulation has some errors:

- Drag should be defined in terms of surface-relative velocity, since the atmosphere rotates with Kerbin

- rdot^2 and thetadot^2 are different units, can't add them - the magnitude of orbital velocity is sqrt(rdot^2 + r^2 thetadot^2)

- Similarly, missing a factor of r in the angular component of drag

- Wrong sign in the exponent for pressure

- Why a factor of 10 in the denominator of mdot? If you use kN for all forces and tons for all masses, then units will be fine

- Shouldn't need to constrain theta, pitching below the horizon should be allowed (probably won't be optimal, but who knows)

- Isp® not defined (more an omission than an error)

Looks nice. You might want some further constraints:

T(t_f) = 0 (No thrusting at the last instant in order to cheat the constraint \dot{r}(t_f) = 0)

I disagree with this one, rdot is a state so last-step thrusting can't make a difference.

Edited August 29, 2013 by tavert

[alterbaron]

Your formulation has some errors:

<snip>

Thanks for finding those problems. I'll have to make a couple of changes.

[tavert]

In case anyone's interested, here's the Mathematica code I was playing with that shows my dynamic equations and derivation of the switching function and singular-arc throttle. I also was working in polar, but instead of theta as a state I was using horizontal speed.

Pressure = AtmSLPressure*Exp[-Altitude[t]/AtmScaleHeight];
(* no cutoff included yet *)

HSpeedSurf = HSpeedOrb[t] - 2*Pi*(Altitude[t] + RBody)/RotPeriod;

SpeedSurfMagnitude = Sqrt[VSpeed[t]^2 + HSpeedSurf^2];

MassDot = -Sqrt[VThrust[t]^2 + HThrust[t]^2]/
   (g0*(IspVac + (IspAtm - IspVac)*Pressure));
(* no Isp saturation for > 1 atm included yet *)

AltitudeDot = VSpeed[t];

VSpeedDot = VThrust[t]/Mass[t] -
   MuBody/(Altitude[t] + RBody)^2 + 
   HSpeedOrb[t]^2/(Altitude[t] + RBody) -
   (1/2)*DragCoefficient*Density1atm*Pressure*
   VSpeed[t]*SpeedSurfMagnitude;

HSpeedOrbDot = HThrust[t]/Mass[t] -
   VSpeed[t]*HSpeedOrb[t]/(Altitude[t] + RBody) -
   (1/2)*DragCoefficient*Density1atm*Pressure*
   HSpeedSurf*SpeedSurfMagnitude;

Hamiltonian = CostateMass[t]*MassDot + CostateAltitude[t]*AltitudeDot +
   CostateVSpeed[t]*VSpeedDot + CostateHSpeedOrb[t]*HSpeedOrbDot;

ThrustTerms = CostateMass[t]*MassDot +
   CostateVSpeed[t]*VThrust[t]/Mass[t] +
   CostateHSpeedOrb[t]*HThrust[t]/Mass[t];
(* terms from the Hamiltonian that depend on thrust *)

VThrust[t] = MaxThrust*Throttle*CostateVSpeed[t]/
    Sqrt[CostateVSpeed[t]^2 + CostateHSpeedOrb[t]^2];

HThrust[t] = MaxThrust*Throttle*CostateHSpeedOrb[t]/
    Sqrt[CostateVSpeed[t]^2 + CostateHSpeedOrb[t]^2];

SwitchingFunction = Assuming[{MaxThrust >= 0, Throttle >= 0},
   FullSimplify[ThrustTerms]/(MaxThrust*Throttle)];

MassDotNew = Assuming[{MaxThrust >= 0, Throttle >= 0},
   Simplify[MassDot]];

CostateMassDot = Assuming[{MaxThrust >= 0, Throttle >= 0},
   Simplify[D[-Hamiltonian, Mass[t]]]];

CostateAltitudeDot = Assuming[{MaxThrust >= 0, Throttle >= 0},
   Simplify[D[-Hamiltonian, Altitude[t]]]];

CostateVSpeedDot = Assuming[{MaxThrust >= 0, Throttle >= 0},
   Simplify[D[-Hamiltonian, VSpeed[t]]]];

CostateHSpeedOrbDot = Assuming[{MaxThrust >= 0, Throttle >= 0},
   Simplify[D[-Hamiltonian, HSpeedOrb[t]]]];

SwitchingFunctionDot = Simplify[D[SwitchingFunction, t] /.
    {Mass'[t] -> MassDotNew,
     Altitude'[t] -> AltitudeDot,
     VSpeed'[t] -> VSpeedDot,
     HSpeedOrb'[t] -> HSpeedOrbDot,
     CostateMass'[t] -> CostateMassDot,
     CostateAltitude'[t] -> CostateAltitudeDot,
     CostateVSpeed'[t] -> CostateVSpeedDot,
     CostateHSpeedOrb'[t] -> CostateHSpeedOrbDot}];

SwitchingFunctionDotDot = D[SwitchingFunctionDot, t] /.
   {Mass'[t] -> MassDotNew,
    Altitude'[t] -> AltitudeDot,
    VSpeed'[t] -> VSpeedDot,
    HSpeedOrb'[t] -> HSpeedOrbDot,
    CostateMass'[t] -> CostateMassDot,
    CostateAltitude'[t] -> CostateAltitudeDot,
    CostateVSpeed'[t] -> CostateVSpeedDot,
    CostateHSpeedOrb'[t] -> CostateHSpeedOrbDot};

SingularSolution = Solve[SwitchingFunctionDotDot == 0, Throttle];

Length[SingularSolution]
(* this should return 1 *)

SingularThrottle = Throttle /. SingularSolution[[1]];

Throttle = Piecewise[{{0, SwitchingFunction < 0},
    {1, SwitchingFunction > 0},
    {SingularThrottle, True}}];

Then setting numerical values for all the constant parameters and boundary conditions and attempting to solve the boundary value problem with both state and costate dynamics, which hasn't worked yet. Tried a simple shooting method but no idea what appropriate values for costates are going to be at time 0.

[alterbaron]

Alright, so I finally got PSOPT working.

As a test, I tried to optimize for the case of the purely vertical ascent of a simple craft:

Command Pod Mk1 + FL-T400 Fuel Tank + LV-909 Liquid Fuel Engine

The goal was to maximize the penultimate height.

Here's some plots from the run:

Here's the output from PSOPT:

PSOPT results summary
=====================

Problem: 					KSP Maximum Ascent
CPU time (seconds):				4.790000e+00
NLP solver used: 				IPOPT
Optimal (unscaled) cost function value: 	-8.153408e+04
Phase 1 endpoint cost function value:		-8.153408e+04
Phase 1 integrated part of the cost: 		0.000000e+00
Phase 1 initial time:				0.000000e+00
Phase 1 final time:				2.525323e+02
Phase 1 maximum relative local error:		2.604555e-03
NLP solver reports: 				The problem solved!

It looks like the craft couldn't quite manage terminal velocity, but its velocity follows the curve quite closely in the densest part of the atmosphere.

I'll try to replicate these results in-game to see how accurate they are.

[tavert]

Nice, looks pretty good. About 135 seconds of burn time is what you'd expect for 2 tons of fuel at an average Isp around halfway between atmospheric and vacuum values for the LV-909. Final altitude may be a bit off what you'd get from an in-game test if you didn't include the surface rotation in the simplified problem. A higher-TWR test will also be interesting to look at, naturally.

What part of installing PSOPT gave you the most trouble, was it the Fortran-related stuff? Remember before you get too far that if you downloaded any of the HSL linear solvers for Ipopt, you can't redistribute binaries. I should've offered some help there, a lot of the libraries are actually easier to build if you use MinGW, you can use the full GCC compilers including for Fortran and run configure scripts almost the same way as if you're using Linux.

Edited August 29, 2013 by tavert

[alterbaron]

Final altitude may be a bit off what you'd get from an in-game test if you didn't include the surface rotation in the simplified problem.

That's a good point. My model for this problem neglected rotation, so there will be a bit of a discrepancy due to that.

A quick test in-game gave a final altitude of about 85Km, so there's some other source of error at work here as well.

(I'd expect that neglecting rotation would make my result too high rather than too low.)

Edit: This might be in part due to the noise in the control function.

What part of installing PSOPT gave you the most trouble?

The biggest problem I had initially was that I was trying to build PSOPT and its dependencies using visual studio. It got messy pretty fast.

I ended up just installing it on my Linux VM, though the installation was still a bit tricky. Had to add some includes to a few of the source files to get it to build, along with editing some Makefiles.

Also had to make some small changes to the gnuplot interface code in PSOPT, since my system's version of gnuplot used a slightly different syntax than the one expected by PSOPT.

Everything else went relatively smoothly. PSOPT's build system is just a bit brittle.

Edited August 29, 2013 by alterbaron

[tavert]

That's a good point. My model for this problem neglected rotation, so there will be a bit of a discrepancy due to that.

A quick test in-game gave a final altitude of about 85Km, so there's some other source of error at work here as well.

(I'd expect that neglecting rotation would make my result too high rather than too low.)

Why are you thinking this direction? The centrifugal acceleration from Kerbin's surface rotation slightly reduces the effective gravity, so in-game with rotation effective gravity is slightly lower than 9.81 m/s^2, whereas in your simulation it's exactly 9.81 m/s^2 (at the surface, mu/r^2 elsewhere). You also start from a bit higher altitude than sea level on the launch pad, though from the resolution of your plot you may have already included that.

The biggest problem I had initially was that I was trying to build PSOPT and its dependencies using visual studio. It got messy pretty fast.

I ended up just installing it on my Linux VM, though the installation was still a bit tricky. Had to add some includes to a few of the source files to get it to build, along with editing some Makefiles.

Also had to make some small changes to the gnuplot interface code in PSOPT, since my system's version of gnuplot used a slightly different syntax than the one expected by PSOPT.

Everything else went relatively smoothly. PSOPT's build system is just a bit brittle.

Ah, interesting. Doesn't look like PSOPT has been updated for some time. I'm more familiar with the dependency libraries than PSOPT itself. Build systems are tricky and often require active maintenance by someone who knows what they're doing, unfortunately.

[K^2]

Keep in mind that the error isn't necessarily on your end. KSP might be sloppy about integrating the forces. A naive Euler integration of drag while rocket is accelerating would result in underestimating drag, and that might be responsible for KSP rocket going a few km higher.

[alterbaron]

The centrifugal acceleration from Kerbin's surface rotation slightly reduces the effective gravity, so in-game with rotation effective gravity is slightly lower than 9.81 m/s^2, whereas in your simulation it's exactly 9.81 m/s^2 (at the surface, mu/r^2 elsewhere).

Good point! I often make stupid mistakes when working with polar coordinates (like the problem of adding r_dot and theta_dot in the earlier problem formulation). In this instance I'm quite happy, though, since it might mean the results are less far off than I had thought.

[tavert]

And did you include the hard atmosphere cutoff at scaleheight*ln(10^6)? Probably best to do this as a problem phase, with an unknown time at the cutoff altitude. Drag up there should be almost negligible, but might integrate out to a couple m/s total?

[alterbaron]

I neglected the hard cutoff in this first run. I wouldn't imagine it would have any measurable effect (other than perhaps making the optimization problem a bit more difficult to solve). Drag acceleration at the cutoff point should be on the order of 10^-3 m/s^2 from a quick calculation.

A naive Euler integration of drag while rocket is accelerating would result in underestimating drag, and that might be responsible for KSP rocket going a few km higher.

I wonder how significant this effect actually is.

Probably best to do this as a problem phase, with an unknown time at the cutoff altitude.

Correct me if I'm wrong, but couldn't you just zero out the drag term in the differential constraints if r>r_cutoff?

Edited August 29, 2013 by alterbaron

[tavert]

Correct me if I'm wrong, but couldn't you just zero out the drag term in the differential constraints if r>r_cutoff?

You could, but this would technically be a discontinuity (though a small one) and might lead to slight numerical troubles. I wonder how/if the auto-differentiation code handles conditional statements, as well. Try both and see if there's any noticeable difference?

Regarding the type of integration in Unity, check this out (scroll down a bit):

http://answers.unity3d.com/questions/48179/rigidbody-acceleration.html

If the way Unity works is you calculate forces at each time step then have Unity update to the next step, then I imagine it must be doing Euler on the velocity, or maybe Stormer-Verlet on the position with a first order velocity estimate (either forward or backward difference)? It could hypothetically be using some kind of multistep method, but those are quite a bit more complicated when timestep isn't constant.

[K^2]

It's almost certainly Verlet. That's practically a standard in gaming industry. Mostly because Verlet is conservative (at least on average) under harmonic potential, meaning you don't run into problems with collisions.

[ThreeKelvin]

I disagree with this one, rdot is a state so last-step thrusting can't make a difference.

Well, yes, \dot{r} is a state, but the problem formulation is continuous, not discrete. There is no last step.

Even if the problem formulation was discrete, e.g. x[n+1] = f(x[n], u[n]), you'd need to ensure that there'd be no last step thrusting, i.e., T[N_f - 1] = 0.

[alterbaron]

Alright, another quick update:

I ran the vertical ascent solver on the same craft, but doubled the TWR.

The results are in agreement with the analytical solution of terminal velocity ascent.

This is quite reassuring. I'll try to solve the full ascent into orbit problem next.

[tavert]

Well, yes, \dot{r} is a state, but the problem formulation is continuous, not discrete. There is no last step.

Even if the problem formulation was discrete, e.g. x[n+1] = f(x[n], u[n]), you'd need to ensure that there'd be no last step thrusting, i.e., T[N_f - 1] = 0.

Why? Your first post seemed to imply that some instantaneous thrust at the final time could satisfy the zero-vertical-speed terminal constraint in an unphysical way? That's not possible, since thrust enters into acceleration, it needs to be integrated over a finite time to have any effect on the vertical speed. For insertion into orbit, you actually usually do want to be thrusting at the end of the time horizon, in order to circularize.

[SaturnV]

Hey guys, good discussion here

K^2 's derivation is still too complicated for me, so I did some numeric simulation, the result shows that stick to terminal velocity is practically optimal

Since I only care the ascend, profiles are compared with delta V used by resistance when they hit the same altitude, the less unnecessary consumption by resistance, the better a profile is.

This way we can avoid simulating the whole journey since we can tell which is better earlier.

Simulation parameters

0.001 // step 1ms
999  // record every 1 second(1000ms)
0  // stop after fuel runs out
0 // start altitude 0m
10000  // stop if hit altitude 10km
5.29158E22 // kerbin mass, kg
600000 // kerbin radius, m
5000 // kerbin scale height, m
1.223 // kerbin sea level atmosphere density, kg/m^3
0.2 // rocket drag coefficient. cross-section area is not included now, because currently the area is calculated by 0.008*mass
15E3 // start mass, kg (mainsail 6t + X200-16 9t)
7E3 // dry mass, kg
544 // fuel flow rate, kg/s (sea level, consider it's constant while it's not, but I only concern the low altitude optimization, so let's make it simple)
1500E3 // thrust, N

Flight report:

Ascend 0, thrust is always maximized

Time: 14.7060 s
Max. Height: 5434.0120 m   // max height this rocket hit
Max. Velocity: 661.8994 m/s
Max. Thrust Acceleration: 214.2710 m/s^2
Max. Drag Acceleration: 144.5729 m/s^2
Max. Net Acceleration: 90.6002 m/s^2
Start ÃŽâ€v:2101.4891m/s
ÃŽâ€v Left(5434 m): 0.0000m/s
ÃŽâ€v used by resi.(5434 m): 1439.5897m/s  // delta v used by gravity and air drag
Fuel Left:0.0000kg(0.000%)

Ascend 1, thrust controlled to make velocity close to terminal velocity by a simple degeneration controller (proportional factor only)

Time: 63.9200 s
Max. Height: 10000.0295 m
Max. Velocity: 267.8369 m/s
Max. Thrust Acceleration: 103.7592 m/s^2
Max. Drag Acceleration: 15.3129 m/s^2
Max. Net Acceleration: 90.6002 m/s^2
Start ÃŽâ€v:2101.4891m/s
ÃŽâ€v Left(10000 m): 599.5070m/s
ÃŽâ€v used by resi.(10000 m): 1234.1452m/s
Fuel Left:1700.0723kg(21.251%)

Ascend 2, thrust controlled to make velocity close to terminal velocity by a advanced prediction module controller

Time: 64.1270 s
Max. Height: 10000.1547 m
Max. Velocity: 267.7240 m/s
Max. Thrust Acceleration: 103.7592 m/s^2
Max. Drag Acceleration: 9.8077 m/s^2
Max. Net Acceleration: 90.6002 m/s^2
Start ÃŽâ€v:2101.4891m/s
ÃŽâ€v Left(10000 m): 600.1021m/s
ÃŽâ€v used by resi.(10000 m): 1233.6630m/s
Fuel Left:1701.9503kg(21.274%)

Unnecessary delta V consumption:

Ascend 0: 1439.5897m/s

Ascend 1: 1234.1452m/s

Ascend 2: 1233.6630m/s

You must noticed that if you simply max your thrust all the way, you even can not make it to 10km, however small fluctuation won't make a difference, if you get to terminal speed close enough high performance is achieved.

Here I only shows you 3 picture, but I also tried the profile of close-to-terminal-velocity-minus-10, close-to-terminal-velocity-plus-10, high-twr-ascend, low-twr-ascend and so on, result is, as expected, the terminal velocity ascend is always the best.

Edited August 31, 2013 by SaturnV