Real solar system mean anomaly

[dietmarwo]

Were did you get your mean anomaly data for the real solar system from?

Is meanAnomaly0 in https://github.com/Krafpy/KSP-MGA-Planner/blob/master/data/rss/bodies.yml referring to epoch0 at 1970/1/1 as defined in https://github.com/Krafpy/KSP-MGA-Planner/blob/master/data/rss/config.yml ?

Checked for Mercury and Venus comparing with https://in-the-sky.org/data/object.php?id=P1 and https://in-the-sky.org/data/object.php?id=P2 and also using https://esa.github.io/pykep/documentation/planets.html backpropagating 30 years (both use 2000/1/1 as epoch 0) and the difference is much larger than what is to be expected due to accuracy errors.

[Krafpy]

The meanAnomaly0 was provided by @TNQOYxNU in this pull request and it is referring to the mean anomaly on 1950/01/01 (epochOffset of -631152000 seconds in https://github.com/Krafpy/KSP-MGA-Planner/blob/master/data/rss/config.yml) as it was apparently validated in this comment by observation. Try backpropagating 50 years (to 1950/01/01). Do you observe issues in the planet positions on the solar system view in my tool compared to what pykep or other tools predict ?

[TNQOYxNU]

data is converted from RSS kopernicus config. Note it added an artificial 23.5 degree inclination to all planet orbit to simulate earth axial tilt.

If you need real world solar system data you can use https://ssd.jpl.nasa.gov/horizons/app.html#/

[dietmarwo]

Sorry, I missed the epochOffset parameter. Backpropagating 50 years (ESAs epoch 0 is at 01.01.2000) gives consistent mean anomalies:

Mecury 5.553920558995985 / 5.559496921012562 Venus 5.432267392125204 / 5.436181913693716

it added an artificial 23.5 degree inclination to all planet orbit to simulate earth axial tilt.

Doesn't this transformation require the ascending node to be equal for all planets? And additionally the "tilt" being consistent with these ascending nodes?

See https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node In your config https://github.com/Krafpy/KSP-MGA-Planner/blob/master/data/rss/bodies.yml this is ascNodeLongitude. They are similar but not equal:

Mercury: ascNodeLongitude: 10.86541167564728 Earth: ascNodeLongitude: 359.9965004168758 Pluto: ascNodeLongitude: 44.36099836994975

Do you observe issues in the planet positions on the solar system view in my tool compared to what pykep or other >tools predict ?

At least I observe different results when I try to compare with your results. Is it possible to dump the exact data for the computed optimal mission? I mean position + velocity at each planet before and after the GA?

https://esa.github.io/pykep/documentation/trajopt.html#pykep.trajopt.gym.cassini2 uses a very similar model (1DSMs, GAs and ESAs Lambert solver) so it could be used to verify your optimizer by comparing results. I wrote a tutorial how to solve these ESA problems in Python: https://github.com/dietmarwo/fast-cma-es/blob/master/tutorials/PYKEP.adoc .

In https://github.com/dietmarwo/fast-cma-es/blob/master/tutorials/MapElites.adoc I used Cassini2 to test a new Quality-Diversity optimization algorithm. QD-algorithms produce "heat maps" where you can identify "oportunity windows" for trajectory planning.

My idea is to complement your tool by a QD-"heat map"-generator (first implemented in Python, may be later to be ported to Typescript). But this only makes sense if the models/implementations produce at least similar results - which they should since both use 1DSMs + ESAs Lambert solver.

By the way: Do you have a reference to a scientific paper for your CR variation method ?

        const genCoeff = Math.pow(1 - genRatio, this.crInc);
        const cr = this.crMax + (this.crMin - this.crMax) * genCoeff;

Looks interesting.

[Krafpy]

Doesn't this transformation require the ascending node to be equal for all planets? And additionally the "tilt" being consistent with these ascending nodes?

I'll let @TNQOYxNU answer to this as the data come from them.

At least I observe different results when I try to compare with your results. Is it possible to dump the exact data for the computed optimal mission? I mean position + velocity at each planet before and after the GA?

I'm am actually in the process of adding a button that prints more details about the trajectory. See issue #29 The test branch is trajectory-to-text. I was only focusing on the maneuvers, but I can try to add more information about the trajectory itself. Can you be more precise about what information should add ? What exactly do you mean by "before" and "after" the genetic algorithm ? To avoid mixing subjects, maybe answer on the related issue, or open a new one.

https://esa.github.io/pykep/documentation/trajopt.html#pykep.trajopt.gym.cassini2 uses a very similar model (1DSMs, GAs and ESAs Lambert solver) so it could be used to verify your optimizer by comparing results. I wrote a tutorial how to solve these ESA problems in Python: https://github.com/dietmarwo/fast-cma-es/blob/master/tutorials/PYKEP.adoc .

That can be useful !

My idea is to complement your tool by a QD-"heat map"-generator (first implemented in Python, may be later to be ported to Typescript). But this only makes sense if the models/implementations produce at least similar results - which they should since both use 1DSMs + ESAs Lambert solver.

It can be an interesting idea, I'll have to look more into it. If you want to discuss it more in details maybe send me a private message on the KSP forum if you have an account there, or describe it in a new issue (but I indeed need to fix everything else first).

By the way: Do you have a reference to a scientific paper for your CR variation method ?

If I remember correctly I made up the calculation myself to simply add a power decay. It looked like it improved the results a bit, but I didn't test it quantitatively. It is possible that a similar idea was developed in a paper. I looked through several papers on differential evolution but the different methods they proposed didn't seem to improve the optimization. I can try to find the papers I looked through if you want. Again if you want to discuss this further maybe contact me directly on the KSP forum.

[TNQOYxNU]

What exactly do you mean by "before" and "after" the genetic algorithm ?

It's gravity assist.

Doesn't this transformation require the ascending node to be equal for all planets? And additionally the "tilt" being consistent with these ascending nodes?

You may already know it's a 3D rotation, but I didn't figure out how RSS mod developer do it either. Since this is intended for planning mission in KSP, I decide strictly following their data is better. So it's better to ask them about the detail of those data. If you need compare it with real world solver, I suggest use real world orbital data in this solver.

[dietmarwo]

Can you be more precise about what information should add ?

The most detailed approach to really see what is going on is to write for each time anything changes: (start, each DSM, each flyby, arrival) the time (or epoch), position, velocity of the spacecraft, and position and velocity of the planets - just to see if I get the same values. A machine readable format (csv-file?) would be nice, something like

type, planet id, time, posX, posY, posZ, velocityX, velocityY, velocityZ 

where type is "start", "DSM" "flyby", "arrival" for the spaceship or "planet" for planet data. For DSMs the planet id would be empty. For DSMs and flybys you should take the outgoing velocity (or spend two lines showing both incoming and outgoing velocity).

In fact for the application inside the game it is less important to have correct "real" data, but values which are consistent with the game. And the resulting trajectories should be compatible with the (simplified) model the game uses. It is possible to feed pykep with keplarian planets (the ones you defined in your yml files). I will try to reproduce trajectories computed by your tool using pykep as soon as I can dump more data.

By the way, optimizing 1DSM / Gravity assist trajectories is one of the hardest continuous optmization problems. If you leave the planet sequence open it becomes a mixed integer problem which is even harder. Which is the reason ESA created these reference probems https://www.esa.int/gsp/ACT/projects/gtop/ so that the optimization research community can compare results. In my experience you often need to use parallelization to solve these problems in reasonable time.

About differential evolution: I did some experiments myself for my own implementation (https://github.com/dietmarwo/fast-cma-es/blob/master/fcmaes/de.py is the Python code, I have a second one in C++). I also change CR (and F) but I prefer an oscillating change switching between two values. For trajectory optimization I combine it with CMA-ES inside a sophisticated parallel meta algorithm (which would be overkill for Typescript I fear). I think what you did is fine in the context of a Typescript application - I probably also would have chosen differential evolution in this context.

[Krafpy]

The most detailed approach to really see what is going on is to write for each time anything changes: (start, each DSM, each flyby, arrival) the time (or epoch), position, velocity of the spacecraft, and position and velocity of the planets - just to see if I get the same values. A machine readable format (csv-file?) would be nice

Alright I'll try to add a CSV download button.

By the way, optimizing 1DSM / Gravity assist trajectories is one of the hardest continuous optmization problems. If you leave the planet sequence open it becomes a mixed integer problem which is even harder. Which is the reason ESA created these reference probems https://www.esa.int/gsp/ACT/projects/gtop/ so that the optimization research community can compare results. In my experience you often need to use parallelization to solve these problems in reasonable time.

Thank you for the link ! That's indeed pretty useful to compare results.

I also change CR (and F) but I prefer an oscillating change switching between two values. For trajectory optimization I combine it with CMA-ES inside a sophisticated parallel meta algorithm (which would be overkill for Typescript I fear).

It's an interesting approach, I didn't think of an oscillating parameter. I will quickly experiment with it when I get time. My differential evolution algorithm is multithreaded using webworkers to distribute the trajectory evaluation of each individual in the population (it increases the speed a bit, but not by much). I don't know if that's the same aspect of the optimization process that you run in parallel.

Try not to mix several subjects in the same issue, open a new one for a specific feature/task. Don't hesitate to direct message me on the KSP forum if you want to discuss ideas that are not suited for a github issue.

[Krafpy]

I added a "Data" button which downloads the data of the trajectory as asked. The rows are:

type, bodyId, timeUT, posX, posY, posZ, velX, velY, velZ

where:

• type indicates the type of event on the trajectory: escape, dsm, flyby, arrival:
  • escape is the moment when the ship exits the departure body's sphere of influence;
  • dsm and flyby are clear;
  • arrival is the moment when the ship enters the destination body's sphere of influence.
• bodyId the identifier of the body related to the event. For dsm's, it is set to the body around which the ship orbits when it performs the maneuver.
• timeUT the time (UT) of the event.
• posXYZ the coordinates of the spaceship at that moment, relative to the attractor of the bodies in the flyby sequence. For flyby's, it's the position when exiting the body's sphere of influence.
• velXYZ the velocity of the spaceship. For dsm's, it's the velocity after the maneuver has been performed.

Tell me if there are missing information.

[dietmarwo]

My differential evolution algorithm is multithreaded using webworkers to distribute the trajectory evaluation of each >individual in the population (it increases the speed a bit, but not by much). I don't know if that's the same aspect of the >optimization process that you run in parallel.

In fact you may have parallelism on different levels:

• inside fitness, for instance in libraries used like BLAS
• fitness itself (what you do)
• parallelizing whole (sub-)optimisation runs.

Scaling increases largely from top to bottom because of the expensive thread/process creation/deletion operation. As larger the "chunk of work" you parallelize the better.

This is the reason the fcmaes meta algorithm breaks down the optimization process in smaller sub-optimizations. This way up to factor 14-18 can be achieved on a 16 core CPU with hyperthreading.

I added a "Data" button which downloads the data of the trajectory as asked.

Thks, will provide feedback as soon as I find time to continue here.

[dietmarwo]

By the way: https://optimize.esa.int/challenge/spoc-trappist-tour/p/trappist-tour contains verification / visualization code for an 1DSM tour using planets defined by their keplerian parameters: https://api.optimize.esa.int/media/problems/spoc-trappist-tour-trappist-tour-1648650812274.py It probably can easily be adapted to planets used here.