Introduction

Over the Summer of 2018, I spent several weeks working on a web-based Martian space simulator. This project was originally a university assignment and some skeleton code of the game was provided in C++, with the physics function that updates the position and velocity of the "Mars Lander". A secondary task was to use rudimentary control theory to program an autopilot capable of landing the spacecraft on the surface. As an extension to this assignment I reprogrammed the C++ into JavaScript to make use of WebGL graphics, which provides more modern graphics functions, and also enables the project to run in-browser, even on mobile. The three.js library was used to interface with WebGL, I'd never used any 3D graphics library before but found that it was very intuitive and straightforward to get started with, and very impressive results could be programmed quickly.

Terrain Generation

One of the first major differences from the OpenGL version that I incorporated into this project was how the terrain on the Martian surface was generated. The original terrain was simply the flat surface of a very large sphere, which was also the planet. To give a convincing terrain, I used a different approach. The Mars planet as seen from orbit is a textured and shaded sphere - but it is around 100x smaller than true scale. When the user passes close to the sphere, beneath a threshold altitude (9500 - 1000 m), the sphere is removed and replaced with procedurally generated terrain which is to scale. The game uses fog to partially shroud this transition.

It took many iterations to develop a satisfactory terrain algorithm and to position and generate the surface terrain properly. I used first a 2D Perlin noise algorithm, and then later a 2D Simplex algorithm (which is very similar) to give the feel of dusty, rolling, Martian hills and valleys. A THREE.PlaneGeometry object is used for each section of terrain, and each polygon is coloured based on altitude (plus some Math.random() noise), the polygons at lower altitudes are darker in colour whilst the higher polygons are lighter.

This initial implementation was good and gave convincing terrain, but what I really wanted was to have the terrain broken down into chunks and to "tile" together so that the player could fly in a straight line infinitely and have terrain procedurally generate, continuously. Adding in this functionality was considerably more complex but I designed a reasonably simple algorithm to handle this task. The algorithm is parameterised by two static constants, the minimum distance to the 'edge' of the world, CHUNK_RENDER_RADIUS, and the maximum distance, CHUNK_DISCARD_RADIUS.

The steps of the algorithm are as follows:

• An origin is defined as the point directly beneath the lander, at an altitude of zero metres, when it passes the threshold altitude. This origin is only computed once and stored, until the lander passes back into space.
• The integer values of which chunk, relative to the defined origin, the lander is currently flying over is calculated. If there's been no change since the last frame then return without updating any chunks.
• If we have changed chunk, we may need to generate new chunks, and may need to discard old ones. Using the chunk positions we enumerate all possible chunk positions in a square of length 2 * CHUNK_RENDER_RADIUS around the lander's position - and discard any with a cartesian distance greater than CHUNK_RENDER_RADIUS. This gives us a set of chunk positions that form a circle around the lander's current chunk - and all of which should ideally be rendered on the next frame. Let's refer to this set of positions as shouldRenderPositions.
• We then iterate through an object which stores all the rendered chunks. Let's refer to this set of positions as haveRenderedPositions.
  • The new chunks that should be generated are the positions that are in shouldRenderPositions but not haveRenderedPositions.
  • The previously generated chunks that should be discarded are the positions that are in haveRenderedPositions but not shouldRenderPositions.
  • Any chunks that are in haveRenderedPositions, but are a cartesian distance of greater than CHUNK_DISCARD_RADIUS, are discarded.

Note that CHUNK_DISCARD_RADIUS is larger than CHUNK_RENDER_RADIUS, so that if the lander very quickly passes back and forth over a chunk boundary, then new chunks won't need to be generated on each crossover.

The positions to be generated are then sent to a web worker running in a separate thread, which generates values from the Simplex algorithm, and computes the vertices of the THREE.PlaneGeometry object. Generating new chunks is CPU intensive, and without the web worker thread there is a noticeable framerate stutter when generating new chunks. This multi-threaded approach enables the asynchronous computation of new chunks without blocking the main render thread, leading to a smoother experience. A short clip showing a perspective of terrain generation not possible in the online demo illuminates this algorithm at work.

Each chunk corresponds to 512x512 metres of actual Martian terrain. There were several key considerations I made when selecting the optimal choice of chunk size. There were pros and cons of both small and large chunks.

For small chunks, chunk boundaries are crossed more often, so the worker must generate new chunks more often. However, generated chunks are smaller, so the overall effect is to reduce spikiness of CPU load and improve smoothness of framerate. Additionally, as the aforementioned algorithm only fires every time a chunk boundary is crossed, there is an “error” of ± half a chunk on the CHUNK_RENDER_RADIUS. That is, the user could be very close to the edge of a chunk (half a chunk away from the middle of the chunk) but until they cross into another chunk, new chunks will not be generated. For smaller chunks the absolute value of this half-chunk error decreases, meaning the CHUNK_RENDER_RADIUS can be lower, proportionally to the chunk size. However larger chunks can also be advantageous, as three.js performs better with fewer, larger chunks, than more, smaller ones. The final value chosen reflects a reasonable trade-off between these constraints.

One issue that had to be overcome when tiling chunks together was the presence of anomalies along chunk borders. Due to the lander travelling through the world at a different scaling to the Simplex noise algorithm traversing through its own coordinate space (too low and the terrain would seem flat, too high and the hills are more like stalagmites), the edges didn't quite line up right on the chunk boundaries. An example of the kind of terrain anomalies is shown below.

To fix this, I used a simple linear interpolation between neighbouring pairs of values at edges of each chunk, and also to the adjoining corner values, which requires interpolation of four values instead of two. To my great surprise, this actually worked fantastically well, and gave a seamless 'stitch' between respective chunks.

Synchronising the physics and rendering

One noticeable improvement from the original C++ OpenGL software is that the physics simulation and OpenGL rendering are no longer tied together in the same thread. Originally, when the simulation speed was slowed down, for example to dt=0.1, (10 physics updates per second), the framerate was also reduced to 10 fps, leading to a very unfriendly user experience. In the javascript software, the physics engine always locks to a fixed framerate, which is separate to the rendering framerate. The main thread framerate and physics update framerate also ideally shouldn't be locked, or having a faster or a slower computer could affect the speed of the simulation.

The simulation computing the physics updates, using the Verlet Integrator, does run in the main thread, and so relies on a requestAnimationFrame call, but it's updates are scheduled as follows. If enough time has passed since the last animation frame for the physics simulation to be updated, then it is. If the rendering has caused the thread to slow and an entire frame has taken, say, 0.1 seconds to render (not ideal) then the simulation accounts for this and computes however many physics updates it missed. This way, the physics engine can always step with a constant value of dt in the Verlet integrator, for a given simulation speed. The real benefit of this system is that the simulation can be sped up, to say 100x real time factor, by simply adjusting dt; the physics engine will still update at the same frame rate, and getAnimationFrame will continue to render the Martian world to the user as fast as the user's device can - for the smoothest possible experience. The physics framerate is also adjusted for large real time factors (speeding time up in the game greater or equal to 200x), so as to prevent dt becoming too large, and precision being lost.

Put simply, the physics thread does this every time a new frame is rendered:

• How long has it been since the last frame, plus any left over time from last frame?
• How many physics updates can we do in that time, based on the physics framerate?
• Do those updates, and add any leftover time that wasn't enough to constitute a whole update to the left over for next frame.

This functionality is only enabled if the game is set to 'online' mode. In the 'offline' mode, the simulation is run without the WebGL render thread, and simply iterates with a fixed dt in a for loop. This offline mode is essential for running simulations as fast as possible, which is necessary when genetically training autopilots (see the next two sections...) and for verifying that the physics calculations compute as expected;

Which confirms the validity of the Verlet integration, as these values are known to be correct. Some slight discrepancies have been introduced by variable altitude terrain and crosswinds in the atmosphere.

The Autopilot

This second part of the project was a little more on topic - away from 3D graphics, though I was very much enjoying my first foray into the 3D world, and back to control theory. One of the main tasks of this assignment was to use control theory to program a simple one-dimensional autopilot that could land on the surface of Mars. I completed the assignment and later extended this autopilot to include a de-orbit burn, that lowers the periapsis into the atmosphere, and also the ability to deploy the parachute. A 'dead-time' engine delay, and an engine lag was also implemented (both 0.1 seconds), and realistic wind speeds of up to 36 km/h were added to the game - the autopilots can land comfortably at any altitude of terrain even with these handicaps (most of the time! Explosions do sometimes happen). Another autopilot uses infinite fuel to inject the lander from anywhere into an orbit with periapsis not below a set threshold - this autopilot works nearly all the time and only occasionally experiences malfunctions.

The configurable parameters for the autopilot are;

• gain - controller gain
• kh - constant to scale altitude value
• delta - offset to scale control output to throttle, as throttle can only be in the range [0, 1)
• parachuteAltitude - deploy parachute if at or below this altitude

I initially hand tuned the autopilots, and was able to land in the most basic scenarios without the added handicaps of engine delays and wind. However, I knew that to construct a high-performing autopilot I'd have to use more sophisticated techniques than tuning by hand…

Genetic Training Algorithm

To achieve finely tuned autopilots I turned to genetic algorithms, and implemented a simple one from scratch in javascript. This was designed to be run in NodeJS so that the programmed autopilot and simulation code could be reused, and so that the final autopilot would be compatible with the in-browser game.

To evolve a high quality autopilot the following steps are undertaken by the algorithm.

• Firstly, a population of autopilots is created with randomly selected parameters for the PD controller that constitutes the core of the autopilot. There is no prior knowledge used when selecting parameter values - each of the control theory parameters gain, kh, and delta, are initialised as 5 * Math.random().
• Then we iterate through each autopilot agent and run offline simulations (using the offline mode mentioned in the simulation section) of landing scenarios, until either the lander reaches the Martian ground, or a simulation-time time limit is reached.
• We can then score each autopilot based on how well it performed in the simulations. A fitness score is computed based on the final state of the lander. The fitness function can be arbitrarily chosen depending on the property to be optimised.
• Then - the next generation is created. Fitter autopilots are weighted such that they are more likely to be selected for reproduction and to pass their parameters (analogous to chromosomes of DNA) into the next generation. Each child is generated by selecting a value for each parameter from a random choice of one of two parents. There is a chance (I used 10%, somewhat arbitrarily) that any given parameter will mutate, which is realised by scaling the value by a random scalar in the range [0.5, 2.0).
• We then repeat this process for as many generations as it take the population to converge its optimum fitness.

Over time the fitness of the population increases, as the fitter agents mate more often and produce fitter offspring, and the weaker are gradually filtered out of the population. The fittest individual in the last generation is then chosen as an autopilot with superior parameters to that that might be hand-tuned.

The true power of the genetic algorithm comes from being able to arbitrarily specify the fitness function. We expect that as we iterate through generations, the value returned from the fitness function increases as the population grows fitter. For my first fitness function I chose to minimise the impact velocity, whilst also minimising the amount of left over fuel. The actual fitness function programmed was as follows:

Taking the Genetic Algorithm Further

Adapting the autopilot population to cope with variable altitude terrain required training with simulated variable altitude terrain (samples from a gaussian distribution, scaled accordingly) - otherwise the autopilots learn to use the atmosphere to slow down, but are overfitted to the density of atmosphere as approaching 'sea level'. If they are trained to always land at 0 metres, that is their radial distance from the centre of Mars if exactly simulation's value for the the radius of Mars, and the ground is above sea level then the atmosphere will be thinner as they approach the ground than they were trained on, and so the lander will be travelling too fast on impact.

The autopilot initially coped very well, but engine delay and lag that were added as an additional challenge for the autopilot induced oscillations in the throttle upon landing. These oscillations can be seen in a clip of the lander coming in to land at ground level.

Due to these oscillations I experimented with switching from a PD controller seeking to bring altitude to zero from a PID controller seeking to take descent velocity to 0.5 m/s, rather than minimising altitude error. This didn't seem to fix the oscillations problem, and also lowered performance of autopilot and made it more difficult to train using the genetic algorithm (taking longer to converge). I reverted to the original strategy. With a lower dead time and lag the autopilot could cope fine, but above 0.2 seconds for either parameter it struggled. The fitness of the population as it tried to converge with my experimental PID controller is shown below:

Even after 300 generations of a 150 strong population the fittest candidates still wasted far more fuel than the fuel optimised PD controller. I only tried this reasonably briefly so it may be that my implementation of the PID controller was flawed - as a full PID controller should be easily capable of preventing these oscillations.

I did however find that the versatility of the fitness function paid dividends when changing the terrain, or optimisation constraints, in any way. For example, after adding wind into the simulation I didn't need to retune the parameters by hand, I just ran my evolution script again and let a new optimum autopilot converge. Even more useful was when I tried optimising the genetic algorithm to not try and minimise fuel usage, but to try and minimise either total descent time, or peak acceleration value. The results produced by evolving an autopilot with descent time or peak acceleration included in the fitness function were fantastic. The algorithm very clearly produced a different optimum solution that fitted the constraints provided. The following descent profiles were produced on a 10 km descent test (with parachute):

These profiles are exactly as we would expect - the minimum time optimised autopilot has the steepest descent down, the minimum acceleration autopilot has the shallowest, and as expected the fuel optimised autopilot outperforms the hand-tuned one. To really appreciate just how effectively these autopilots converge we need to look at the fuel, landing times, and peak acceleration values for each autopilot:

Each autopilot trained to optimise a given specification outperformed every other autopilot, on that specification. The final values of parameters were as follows:

It is interesting to note some differences in the final parameter values based on the specification of the autopilot. Unsurprisingly, the minimum time autopilot has a much lower parachute deploy altitude, and the minimum peak acceleration autopilot has a much higher altitude. It is also worth noting that whilst minimising time or fuel leads to a very similar value for kh, minimising acceleration seems to lead to a lower value. This fits with the earlier results plotted for varying values of kh - we would expect that a shallower descent profile favours lower peak acceleration. The simplicity of the genetic algorithm is wonderful - to change the parameter the algorithm was optimising, all that was done was to change the line:

To minimise peak acceleration, to the following:

Or to minimise descent time, to the following:

Note that for the acceleration fitness, up to 400 'fitness points' can be 'earned' if the lander doesn't accelerate too fast. It is important to use this method of giving a reward of a capped size, that is reduced if the autopilot performs poorly, rather than a negative feedback 'punishment' where the maximum acceleration is subtracted from the fitness score. The time optimised fitness scoring uses a similar logic with the value cap at 1000 due to the different expected range of values. The reason for this is because the fitness increase from optimising fuel, acceleration, or time, is only available as a bonus to the autopilot, once it has successfully landed, as discussed in the section above. If a subtractive method is used the autopilot learns to never earn its bonus points, because the bonus points are always a punishment, and so the autopilot learns to land very slowly - but not so slowly that it will avoid exploding on impact.

Other Improvements

In addition to the improved terrain generation, simulation thread and timing handling, and genetically evolved autopilots, there were many small tweaks and modifications I added into this project, for general aesthetics and realism.

Rotation of camera, terrain, and lander

Getting everything in the right rotation was possibly the hardest part of this entire project. I'd never used 3D graphics before so Euler angles and quaternions were new to me - and proved challenging to get to grips with at first. I made extensive use of three.js's ability to group meshes together, with a parent-child structure. It took a long time to work out the required structure, but after a lot of trial and error, I worked it out. The key points are as follows:

• landerGroup is the parent group containing all lander objects, and includes the THREE.PerspectiveCamera object that views the scene. When the position of the lander is updated, it is landerGroup.position that is changed.
• The camera object is in its own group. This camera group is added to the landerGroup. This camera group is rotated when transitioning between terrain and space mode, so that the 'up' orientation is always the surface normal vector.
• In space, the camera group always takes the same orientation - the up axis is the axis of rotation of Mars (which happens to the y axis, in the scene constructed).
• To rotate the lander using the WASD keys, the landerMesh which is added directly to the landerGroup has its quaternion updated based on quaternion rotations of unit vectors, and its angular velocity.
• The particle exhaust trail group is added to the landerMesh, whereas the explosion animation group and parachute group (parachute mesh and cords) is added to the landerGroup. This is because the exhaust trail should always point out the bottom of the lander, but the parachute and explosion's quaternion rotations are not dependant on the landerMesh.
• The terrain itself is rotated to a completely different orientation depending on which direction the lander lands from. This is so it always appears from the right direction and angle, that is if the lander travels straight down its own position vector, no matter where over Mars they are, they will always be travelling normal to the plane that the terrain is oriented in. To compute the distance travelled, which needs to be known when generating new chunks (see earlier), two unit vectors are created that initially corresponded to the X and Z axis in the terrain, but have had the same rotation quaternion as the terrain, applied to them. The position on each of the North/South East/West axis is then the 3D position vector of the lander, relative to the terrain origin (see earlier again), dotted with each of these unit vectors.

Exhaust trail and explosion animations

I wrote these animations from Scratch using three.js meshes and various Math functions in javascript. The exhaust trail is rather simple, just a trail of octahedral meshes with the right opacity / velocity / rotation / colour in the right place, that turns on and off at the right times and scales with the simulation's real time factor, but I'm rather pleased with the end result. The explosion animation uses essentially the same behaviour but simplified with no starting offsets or particle re-generation at the end of lifecycle, and with velocities generated differently.

Other Miscellaneous Points

• The background colour changes and the starry skybox fades out when approaching and travelling through exosphere. Fog mostly shrouds the transition between the “Space mode” and the “Terrain mode”.
• A directional light source emanates from what would be the Sun, which I didn't get round to texturing in, but does light up Mars, beautifully catching light on craters at sunrise and sunset. Information for this shading behaviour comes from the topological satellite scan of Mars set as a bump map. I am very grateful to Planet Pixel Emporium for their fantastic textures and bump maps of Mars.
• A lightweight 1D Perlin noise implementation is used to generate the random gusts of wind, plus a base magnitude that is always present.
• Although Mars does graphically rotate there is no rotation considered in the physics engine. This is an inaccuracy that given enough time I'm sure I could fix, but the rotation of Mars whilst on the planet surface would introduce new challenges to overcome - as the terrain Mesh would have to all rotate, along with the camera group. I decided I'd probably spent enough time working on this project and wanted to get onto some other projects, so left this problem unfixed!

Overall Summary

All in all this was a project that took rather longer and went into more depth than I originally intended. I initially set out to see if it would be possible to improve on the OpenGL software by writing it in WebGL, and started with the hardest parts of the project that I thought would be the limiting factors. It took me a good few days to sink my teeth into these problems, but by the time I'd done that I decided I was already invested in working it all the through.

I am certainly pleased with the end result, and am glad to have built on my experience in front-end development and Javascript in general, and I look forward to using this knowledge in future projects…