Spacecraft Control Toolbox To The Moon

Today, I will discuss two functions in release 2020.1 of the Spacecraft Control Toolbox (SCT) which can be used to get your spacecraft into a lunar orbit. They are LunarTargeting.m and LunarMissionControl.m. They are demonstrated together in the script LunarMission.m.

LunarTargeting.m produces a transfer orbit that starts at a Low Earth Orbit (LEO) altitude and ends up passing by the Moon with a specified perilune (periapsis of the Moon) and lunar orbital inclination. Its novel approach to the patched-conic-sections model of multibody orbital transfers uses the solution to Lambert’s problem to target a point on the gravitational boundary between the Earth and the Moon. Then it numerically optimizes over points on that surface until the initial velocity of the transfer is minimized. LunarTargeting.m requires the MATLAB optimization toolbox.

LunarMissionControl.m implements a control system which enables a spacecraft to propulsively enter lunar orbit. Like the other control systems implemented in the SCT, it stores its active state and degrees of freedom in a data structure, and accepts a list of commands as arguments. The commands we’ll see used here are ‘initialize,’ ‘lunar orbit insertion prepare,’ ‘align for lunar insertion,’ and ‘start main engine.’

LunarMission.m ties them both together and simulates a spacecraft, down to the attitude-control level. The simulation includes power and thermal models. The spacecraft can be controlled by reaction wheels or thrusters. Forces from the Sun, Earth, and Moon are included. The spacecraft starts on the trajectory returned by LunarTargeting.m, then acts in accordance to commands to LunarMissionControl.m. It takes the spacecraft 4.5 days to get to perilune, at which point it inserts itself into lunar orbit. Let’s take a look!

A figure produced by LunarMission.m. This figure shows the whole trajectory, from Low Earth altitude to a complete orbit around the Moon after insertion.

Take a look at the above figure. This is the entire mission trajectory in the Earth-Centered Inertial (ECI) frame. We can see the initial transfer orbit as the red line. Then it approaches the blue line (the Moon’s orbit), and begins corkscrewing around it after orbital insertion. Let’s look at that insertion in close-up:

One figure produced by LunarMission.m. This figure shows the lunar orbit insertion.

The above figure shows the final part of the trajectory in Moon-centered coordinates. The red line starts as the spacecraft passes the imaginary gravitational boundary between the Earth and the Moon. It falls closer to the Moon, and at its closest point, fires its engines to reduce its velocity. You can’t see it in this figure, but that process is actually resolved on a 2 second timescale. The spacecraft is commanded to point retrograde using a PID controller, waits until it has pointed correctly, then fires its engines for a prescribed duration. If you look closely, you will see that moon has a 3 dimension surface courtesy of the Clementine mission.

Let’s finish this post off with some technical details:

On the far left, you can see the reaction wheel rates. They stay at zero for 4.5 days, as the spacecraft coasts. Then, when the craft is commanded to point retrograde for its orbital insertion, you can see wheels 2 and 3 spin up. Wheel 1 stays near zero; its vertical scale is 10^-16. Then in the center, you can see fuel use. The only fuel use is the insertion burn, so fuel stays constant until 4.5 days in. Less than 2 kg of fuel is used for this example, as the spacecraft is a 6U cubesat. On the right, the components of the body quaternion are displayed. Again, they are constant until 4.5 days in, when the craft is commanded to point retrograde.

I hope you’ve enjoyed this demonstration of how to simulate a lunar mission with the SCT! For more information on our toolboxes check out our Spacecraft Control Toolbox for MATLAB. You can contact us directly by email if you have any questions.

Renewable Home Power, Backup, Heat and Air Conditioning

Princeton Satellite Systems has been in a leader in renewable energy with its SunStation home solar power system with battery backup. We introduced this product back in 2013. SunStation has lithium-ion phosphate batteries, the most stable and reliable batteries for home use. The core of the system is the Outback Inverter that seamlessly switches from grid power to internal power.

The solar system in the installation produces 7.3 kW of power, much more than the house needed for electric power including charging a Nissan Leaf and Toyota Prius Prime. The heating and air conditioning system was nearing its end-of-life so we decided to replace it with a geothermal heat pump. A heat pump is essentially an air conditioner that can both reject heat to a source and absorb heat from a source. The problem with both is that when the outside temperature is high, for rejecting heat, and low, for absorbing heat, the system loses efficiency. Modern air-source heat pumps are very efficient but do need backup resistance heating in some climates.

A ground source heat pump, or geothermal heat pump, uses the ground as the medium for absorbing or rejecting heat. The option we chose, due to land constraints, is to have two wells several hundred feet deep as the source. Alternatives are trenching, or a pond if you have one in your yard. The ground is always at around 50 deg F. The system was sized so that it rarely, if ever, needs resistance heating.

The geothermal system, which is made by WaterFurnace, was installed by Princeton Air. No changes to the SunStation were needed. The core geothermal system is shown below. The valves to the ground loops are in the foreground and the geothermal system is on the left.

The lines that run to the outside ground loops are shown below.

The system has a preheater for the (still gas) hot water heater. The gas water heater was less than a year old, so it didn’t make sense to replace it. The preheater is an electric hot water heater that does not have the heating coils connected.

The SunStation is shown below. The Outback inverter is on the bottom left. The boxes on top provide arc protection, which is now included in the inverter. The batteries on on the right and the battery management electronics between the inverter and the battery cabinet.

The well digging was quite a project. This picture shows the drilling rig.

This second picture shows the yard after the drilling was complete. Drilling took three days total.

The following system shows the SunStation with geothermal in operation. The Prius Prime is charging which is most of the load. The system is still sending considerable power to the grid. On average the house powers itself and two other houses.

Geothermal, with solar and battery backup is the ideal solution for new homes and for renovations to existing homes. There is no reason to even have a gas hookup anymore. Contact us at SunStation for more information

Three-Burn Solutions for Very Large Inclination Change

This is the second of two blog posts which introduce functionality from the Orbit Transfer Module of the Spacecraft Control Toolbox (SCT). The Orbit Transfer Module has lots of tools to allow you to numerically optimize the engine burns that a spacecraft needs to apply to go from some initial orbit to its target orbit, with a minimum of fuel used, time elapsed, or some other metric. You can do this to impulsive burns, continuous (low-thrust) burns, and a special model that splits large impulsive burns into many small impulsive burns.

In this blog post, I’d like to discuss an interesting thing that happens when a spacecraft aims to change its inclination by more than 37 degrees. In real life, of course, this would never happen. This is a truly impractical amount of inclination change. If you were truly faced with a mission that required coverage of orbits that were different by 37 degrees, you would simply launch two spacecraft.

But if you were to explore this hypothetical, you would find an interesting feature emerge. Here we will consider the case that a spacecraft starts in a circular, equatorial LEO and targets a circular LEO which is 49.5 degrees inclined to the equator. The well-known solution is to burn once, at the common node of the initial and target orbits, to change direction by the target inclination:

The 1-burn inclination change solution for 49.5 degrees. This requires a truly impractical amount of delta-V.

As we can see from the above figure, this inclination change requires 6.47 km/s of delta-V, almost the spacecraft’s entire orbital velocity!

But, as the Orbit Transfer Module’s OptimizeElementsImpulsiveSearch finds, there is actually another transfer which can marginally reduce the delta-V required for this inclination change. It finds that the following transfer saves delta-V:

  • Burn prograde to raise apoapsis (and a very slight inclination change)
  • Coast up to apoapsis, where inclination change is cheaper
  • Perform an inclination change
  • Coast down to periapsis
  • Burn retrograde to lower apoapsis (and a very slight inclination change)
The 3-burn inclination change solution for 49.5 degrees. Loops out to high altitude before changing inclination. This also requires an impractical amount of delta-V, but it is a smaller amount.

As we can see from the above figure, this transfer only expends 5.97 km/s of delta-V. While this is still impractical, it is interesting that there exists this other category of optimal inclination transfers which exists for high inclination changes.

It is interesting to note that, for 0 through 37 degree inclination changes, the direct approach is superior. Then, from 37 degrees to 60 degrees, this 3-burn solution produces a smaller required delta-V. Above 60 degrees, the intermediate, high-apoapsis orbit actually exceeds escape velocity so the transfer takes infinite time.

Hohmann Transfer, but Faster: Optimizing for fuel and elapsed time

In this blog post, we are going to introduce you to functionality found in the Orbit Transfer Module of our Spacecraft Control Toolbox. The Orbit Transfer Module has lots of tools to allow you to numerically optimize the engine burns that a spacecraft needs to apply to go from some initial orbit to its target orbit, with a minimum of fuel used, time elapsed, or some other metric. You can do this to impulsive burns, continuous (low-thrust) burns, and a special model that splits large impulsive burns into many small impulsive burns.

The following is an example of a case in which the impulsive burn optimization functions can reveal new and interesting classes of orbit transfer trajectories.

The spacecraft starts in a circular, equatorial low Earth orbit (LEO) and targets a circular, geosynchronous equatorial orbit (GEO). The fuel-optimal trajectory is well known: a Hohmann Transfer, wherein the craft burns once to raise its apoapsis to the target altitude, then at apoapsis it burns again to raise its periapsis to the target altitude. It looks like this:

A Hohmann Transfer from LEO to GEO, which is also the fuel-optimal transfer. Burn once to raise the apoapsis to the target altitude, then again to raise the periapsis.

The craft starts in the blue LEO orbit, then at the big black arrow it burns its engine until 2.42 km/s of delta-V has been applied, then coasts for 5.28 hours on the yellow transfer orbit until it reaches the altitude of the purple dashed, target orbit. Then a final insertion burn until 1.47 km/s has been expended, and GEO has been achieved! Final tally: 3.89 km/s of delta-V has been expended and it took 5.28 hours to get to the requested orbit.

This Hohmann transfer results from optimizing for the minimum delta-V expended, but what happens if you optimize both fuel and coast time? If you care just as much about decreasing the elapsed time by one hour as expending 1 km/s of fuel? The function OptimizeElementsImpulsiveSearch can solve this problem also. Take a look:

In these plots, we’ve told the optimizer that we care about both fuel and time. In the plot on the left, we’ve specified a tradeoff of 1 hour of elapsed time to 1 km/s of delta-V. If the optimizer can decrease the coast time by more than 1 hour at the cost of 1 km/s delta-V, it will do it. And we can see that the optimizer has in fact returned a faster, more delta-V intensive trajectory. Final tally: 5.06 km/s of delta-V has been expended and it took 3.29 hours to get to GEO.

In the plot on the right, we’ve specified a tradeoff of 1 hour of elapsed time to 2 km/s of delta-V, indicating that we care much more about elapsed time than the previous case. And again, we can see that the optimizer returns an even faster, even more delta-V intensive trajectory. Final tally: 6.23 km/s of delta-V has been expended and it took 2.55 hours to get to GEO.

Time tradeoff:
1 hour is worth how much delta-V?
Delta-V expendedTime spent coasting
0 km/s (I don’t care about time)3.89 km/s5.28 hours
1 km/s (I care about time a little)5.06 km/s3.29 hours
2 km/s (I care about time a lot)6.23 km/s2.55 hours
Numerically optimized LEO GEO transfer trajectories, given different priorities between delta-V and time. These 3 points are part of a family of trajectories which are called fuel-time optimal.

We can visually verify the fingerprints of this faster transfer on the plots. The transfer orbit in the Hohmann transfer ends right at the target altitude, but for 1km/s = 1 hour, the transfer orbit extends out beyond this altitude, and for 2 km/s = 1 hour, the transfer orbit extends farther still. This is a hallmark of being a faster transfer orbit. The spacecraft doesn’t actually loop all the way out there; it performs a braking and insertion burn when it gets to the required altitude.

There are real-world scenarios where you’d want to do this. If your spacecraft is sensitive to radiation and you want to minimize the time spent in the Van Allen belts, for example. Or taking a more science-fiction approach, what if you need to urgently resupply or rescue a stranded crew?


Yes, that is the name of the very famous computer game.

In our last blog post we talked about optical navigation for lunar missions. Optical navigation is also very valuable for asteroid missions. It is being used today on OSIRIS-REx. The system we developed can be used anywhere in the solar system without, hopefully, too much attention from the ground.

Now that you’ve decided to go to an asteroid, where are they? As it turns out there is a fabulous website with a downloadable database of asteroids You can download a file with thousands of asteroid.

In Version 2020.1 of the Spacecraft Control Toolbox, we’ve added a function to plot the asteroids. It reads the file and returns the orbital elements, name, number and epoch. Here are two examples of asteroids in motion.

Ceres is the most famous asteroid. You may find one named after you! If you want the bigger picture, here is a second plot with 1000 asteroids. The circles are the orbits of Earth and Jupiter. The function will propagate all the asteroids you select if you don’t request any outputs.

Contact us for more information!

Optical Navigation to the Moon

There is great interest in lunar missions. The U.S. plans to land astronauts on the moon in this decade. Several commercial companies are working on landers. Many other national space programs are working on their own landers and rovers.

As with all spacecraft, you need to know where you are. The traditional way is to use a communications link from the ground to get range and range rate (via Doppler shift.) This works well but you need to collect data over a long period of time to get an 3-dimensional fix.

Recently, NASA discovered via the Magnetospheric Multiscale Mission (MMS), that GPS could be used at altitudes higher than the GPS altitude, and maybe all the way to the moon!

Princeton Satellite Systems developed an Optical Navigation System for NASA as part of an SBIR project. It uses two cameras. One is a navigation camera mounted on a 2-axis gimbal that looks for nearby planets or asteroids. A second camera points out of the orbit plane just looking at the star field. The navigation camera has sufficient dynamic range so that it can image a planet or moon and still see stars. The whole system works as a sextant, that has been used by mariners for hundreds of years. The Apollo astronauts used a sextant for backup navigation. Our system automates the process so that it is fully autonomous. Our simulations for the SBIR were for deep space missions, like simulated NASA Messenger and New Horizons spacecraft, and for communication satellites.

Recently we’ve customized the system to lunar missions. The following images are from a MATLAB script that simulates a transfer from the Earth to the Moon. The first is from the star camera pointing out of the orbit plane. The numbers are the star id from a list and don’t correspond to numbers in the Hipparcos Catalog which is used for the star simulation.

The second is the navigation camera, that points at the moon. The star moves as the relative angle changes. The circle around the moon is its disk.

The third shows the trajectory. It starts from behind the Earth. A lunar orbit insertion is done when the spacecraft reaches the moon.

Here is the spacecraft in lunar orbit.

The simulation runs until the spacecraft reaches the moon. The following video shows the simulation in action.

The simulation uses functions coming out in the 2020.1 release of the Spacecraft Control Toolbox. Contact us for more information!

Fun with Wavelet Analysis

Not all the new functions in 2020.1 are specific to spacecraft. We have also been hard at work adding new functionality to the core toolbox. Here, I’d like to give an example of one of our new functions for performing a Wavelet analysis.

But what is a Wavelet analysis? Well, you plot the Wavelet transform of a signal when you want to visualize how the frequency spectrum changes in time. The Wavelet transform is a lot like a Fourier transform that you perform at every possible starting point, with an appropriate window function multiplied in so that you’re only looking at a portion of the signal.

But there’s one added wrinkle, because the frequency spectrum at a specific frequency at a specific time doesn’t technically exist. It’s not technically possible to know what the component of a signal at 100 kHz is at 0.5 seconds in, because the frequency spectrum depends on the entire signal. There has to be some trade-off between time resolution and signal resolution. If we look at a very long chunk of the signal, we can nail down its frequency components very well but we can’t see them change quickly. If we look at a very short chunk of the signal, we know precisely when the frequency changes but we can’t tell the difference between two similar frequencies. It’s a trade-off.

Now let’s get to the examples! The new function in 2020.1 is called WaveletMorlet because the specific window function we use is called the Morlet wavelet (A Wavelet transform using a Morlet wavelet is also called a Gabor transform). Here’s the signal that we’ll be analyzing:

A test signal. We want to visualize how the frequency spectrum changes in time.

We already know what we’re going to expect in this example. It looks like there’s a persistent, low-frequency component, then a higher-frequency component whose frequency goes up, peaks around 0.25 seconds, then goes down, and bottoms around 0.75 seconds. Here’s what the Wavelet transform looks like:

A visualization of the wavelet transform of the test signal. Wavelet width parameter: 10

Great! Exactly what we expected. This was a simple case, but you can imagine how this analysis would be useful if there were a greater spread in frequencies, a longer signal, or both.

Now, let’s explore that trade-off that I mentioned earlier. What does the signal look like when you choose a different value on that trade-off? For the above analysis, I kept the default wavelet width parameter of 10. Here’s what it looks like when we prioritize time resolution over frequency resolution by choosing 5, then frequency resolution over time by choosing 25:

For a wavelet width parameter of 5, all that happens is that the signal gets broader in the frequency direction. For 25, what’s happening here? Sure, at 0.25 seconds it appears that the visualization is able to nail down the frequency to a tighter band, but what’s happening to the rest of the image? The answer is that the frequency is changing too fast for this chosen time resolution. The signal doesn’t spend long enough at any given frequency for the algorithm to identify a significant component there.

Thanks, all, for tuning in to this update from PSS, and thanks for this opportunity to get into the nitty gritty of one of our new mathematics functions!

Rendezvous Made Simple

In the “good old days” the only people worried about rendezvous were those designed space missions with crews. ISS established the need for robotic rendezvous and docking on a regular basis.

Orbit dynamics can be complex. If you are looking at rendezvous with any other spacecraft in a very different orbit you can start with Lambert’s Time-of-Flight algorithm. Given initial velocity and position vectors, and desired final position and velocity vectors and time of flight, it will give you the initial impulse velocity change needed to rendezvous. There are numerous formulation as it is complex math problem.

If you happen to be close to your target you can formulate your orbits as a relative orbit problem with Hill’s equations, shown below in state space form.

n is the orbit rate of the target spacecraft. x, y and z are the position of the “chase” spacecraft in the Hill’s frame. a is the control acceleration. You want to reduce the positions and velocities to zero. This can be done with a Proportional Derivative (PD) Controller, or with a Linear Quadratic (LQ) Controller. If your chase and target spacecraft have GPS it is relatively easy to find this state vector in the above equation. A PD controller will ignore the coupling in the above equations while the LQ will accommodate the coupling.

It is interesting to look at the gain matrices for the two cases and the corresponding eigenvalues. We tweaked the PD to make its position gains close to that for the LQ. The PD is designed for a damping ratio of 1. The eigenvalues are identical. The cross-axis gains are small, but non-zero.

Gain Matrix LQ
0.0032 -0.0001 -0.0000 0.0796 0.0000 -0.0000
0.0001 0.0032 -0.0000 0.0000 0.0796 -0.0000
0.0000 -0.0000 0.0032 0.0000 -0.0000 0.0796

LQ eigenvalues
-0.0398 + 0.0409i
-0.0398 – 0.0409i
-0.0398 + 0.0386i
-0.0398 – 0.0386i
-0.0398 + 0.0397i
-0.0398 – 0.0397i

Gain Matrix PD
0.0036 0 0 0.1200 0 0
0 0.0036 0 0 0.1200 0
0 0 0.0036 0 0 0.1200

PD eigenvalues
-0.0398 + 0.0409i
-0.0398 – 0.0409i
-0.0398 + 0.0386i
-0.0398 – 0.0386i
-0.0398 + 0.0397i
-0.0398 – 0.0397i

The simulation results for the LQ are:

And for the PD are:

The results are very close. The PD has no overshoot, as expected. The LQ is slightly faster but has some overshoot. Both get the chase spacecraft to the target in a few minutes, assuming, of course, that you have the acceleration capability shown in the plots.

Both are linear controllers. You can approximate a linear controller with thrusters by using pulse width modulation. An issue will be the minimum impulse bit of the thrusters, that will lead to a minimum velocity and position error that can be achieved.

This script is included in the Spacecraft Control Toolbox 2020.1 coming soon!

Flying Near the ISS

Many small satellites are launched from the ISS or near the ISS by one of the transfer vehicles. A new function in the Spacecraft Control Toolbox 2020.1, coming in May, allows you to visualize your spacecraft near the ISS. Here is an example. You can also display a trajectory.

Easy Low Thrust Orbit Raising

It is sometimes necessary to change your orbit semi-major axis, ascending node and inclination with a low-thrust engine. It is easy to do, as long as you can point your engine along orbit normal and tangential to the orbit.

It is easiest to see how this is done by looking at the Gauss’ Variational equations, simplified for small eccentricity.

I is inclination, a is semi-major axis, mu is the gravitational parameter, omega is argument of perigee, nu is true anomaly, Omega is ascending node. a_T and the orbit tangential acceleration and a_Nis the orbit normal acceleration.

The resulting simulation is shown below. Mode 0 is semi-major axis change, Mode 1 is ascending node change, Mode 2 is inclination change and Mode 3 is off. It is best to change inclination and ascending node at the highest semi-major axis. You should change ascending node at the lowest inclination. The burns are done where the rate of changes are higher. Some change in inclination and ascending. node will happen when the other is being corrected.

The script for this simulation with the controller is part of the Spacecraft Control Toolbox Release 2020.1 coming in May.