Planning a Collision-Avoidance Maneuver with Python
Every satellite operator eventually receives the message nobody wants: a piece of debris will pass within a few hundred meters of your spacecraft in a few hours. Do nothing and you accept a small but real chance of losing a very expensive asset. Maneuver and you burn fuel, which is your lifetime, and you also add new uncertainty to your own orbit.
So how do you decide how much to push and when? In this article we turn that question into a small optimization problem and solve it with Python. We will build a physical model, search for the risk-minimizing maneuver, verify the result with a Monte Carlo simulation, and visualize everything in one figure, including two 3D surfaces.
The Scenario
A satellite flies in a circular orbit at 500 km altitude, with an orbital period of about 94.6 minutes. A conjunction warning tells us:
- The predicted miss distance at the time of closest approach (TCA) is 40 m radial and 60 m along-track.
- The combined position uncertainty (1σ) in the encounter plane is 50 m radial and 200 m along-track.
- The combined hard-body radius of the two objects is $R = 10$ m.
- The warning arrives 4 orbits (about 6.3 hours) before TCA, so any burn must happen within that window.
We can perform one small along-track (prograde) burn of size $\Delta v$ at a lead time $t$ before TCA, followed later by an equal return burn to restore the orbit. Our decision variables are $(\Delta v,\ t)$.
Modeling the Physics
How a tiny burn becomes a large miss distance
For a near-circular orbit, the Clohessy–Wiltshire equations describe relative motion in the radial ($x$) and along-track ($y$) directions. For a tangential impulse $\Delta v$ applied $t$ seconds before TCA, the displacement at TCA is
$$
\Delta x = \frac{2}{n}\bigl(1-\cos nt\bigr),\Delta v, \qquad
\Delta y = \frac{4\sin nt - 3nt}{n},\Delta v
$$
where $n=\sqrt{\mu/a^{3}}$ is the mean motion. The secular term $-3nt$ is the key: the longer the lead time, the more a millimeter-per-second nudge is amplified. Small burns done early are far cheaper than big burns done late.
Execution error grows with lead time
A real thruster never delivers exactly the commanded impulse. We model the burn error as
$$
\sigma_{\Delta v}=\sqrt{\sigma_{0}^{2}+(k,\Delta v)^{2}}, \qquad \sigma_0 = 2\ \text{mm/s},\quad k=0.03
$$
and this error is amplified by the same geometry, so the covariance in the encounter plane becomes
$$
\sigma_x^{2}=\sigma_{x,0}^{2}+\bigl(c_x,\sigma_{\Delta v}\bigr)^{2},\qquad
\sigma_y^{2}=\sigma_{y,0}^{2}+\bigl(c_y,\sigma_{\Delta v}\bigr)^{2}
$$
Here $c_x$ and $c_y$ are the coefficients multiplying $\Delta v$ in the displacement equations above. This creates the central trade-off of the article: a long lead time reduces the required $\Delta v$, but it also inflates the uncertainty of where we end up.
Collision probability
With a small hard-body radius compared to the uncertainty, the collision probability is well approximated by
$$
P_c \approx \frac{R^{2}}{2\sigma_x\sigma_y}\exp!\left[-\frac12\left(\frac{m_x^{2}}{\sigma_x^{2}}+\frac{m_y^{2}}{\sigma_y^{2}}\right)\right]
$$
where $(m_x, m_y)$ is the miss vector after the maneuver.
The objective
We minimize the expected loss
$$
J(\Delta v, t) = C_{\text{col}},P_c(\Delta v, t) + 2,C_{\text{fuel}},\Delta v
$$
with $C_{\text{col}} = 5\times10^{8}$ USD (loss of the satellite and its mission) and $C_{\text{fuel}} = 5\times10^{6}$ USD per m/s. The factor 2 accounts for the return burn. In addition to the free optimum, we also solve the problem under a typical operational safety rule:
$$
\min_{\Delta v,,t}\ J(\Delta v, t)\quad \text{subject to}\quad P_c\le 10^{-5}
$$
and the pure minimum-fuel version of the same constraint,
$$
\min_{t}\ \Delta v_{\text{req}}(t), \qquad \Delta v_{\text{req}}(t)=\min{\Delta v : P_c(\Delta v’, t)\le 10^{-5}\ \ \forall \Delta v’\ge\Delta v}
$$
The Complete Source Code
Everything is in a single cell.
1 | import time |
Code Walkthrough
Section 1: Scenario constants
The mean motion N_MM and the orbital period T_ORB follow directly from the semi-major axis of a 500 km circular orbit. MISS, SIG_NOM, and R_HB are the conjunction data (predicted miss, uncertainty, and hard-body radius). SIG_DV0 and K_DV define the thruster execution error, and C_COL and C_FUEL are the economic weights of the objective. PC_LIMIT is the operational safety threshold of $10^{-5}$.
Section 2: The physical model
cw_coeffs returns the two coefficients $c_x$ and $c_y$ of the Clohessy–Wiltshire equations. The lead time is given in orbits, so the phase angle is $nt = 2\pi \times \text{orbits}$. Working in orbits keeps the plots easy to read.
encounter_state is the heart of the model. It shifts the nominal miss vector by the maneuver-induced displacement and inflates the covariance by the amplified execution error. The np.where(dv > 0, ..., 0) guard is important: when no burn is performed there is no execution error, so the do-nothing baseline is not unfairly penalized.
collision_probability implements the analytic $P_c$ formula, and expected_cost adds the fuel term to obtain $J$. Because everything is written with NumPy operations, these functions accept scalars and arrays of any shape, which is what makes the next section fast.
pc_monte_carlo is an independent check. Instead of sampling the Gaussian directly (which would need billions of samples to resolve a probability near $10^{-6}$), it uses importance sampling: points are drawn uniformly over the hard-body disk and weighted by the Gaussian density. The estimate is $\pi R^{2},\mathbb{E}[f(\mathbf{x})]$, which is accurate even for very rare events with only two million samples.
Section 3: Vectorized grid search
We evaluate $J$ and $P_c$ on a $301\times301$ grid of $(\Delta v, t)$ in a single call thanks to broadcasting. To show why this matters, the code also evaluates a $100\times100$ grid with a naive double for loop and compares it with the vectorized call. The assert np.allclose(...) line guarantees that both versions produce identical numbers. The vectorized version is roughly two orders of magnitude faster, which is what allows us to explore the design space interactively and extend it to finer grids or additional parameters.
Section 4: Optimization
The grid minimum gives a robust starting point in a landscape that has several local minima (the sine and cosine terms make the cost surface wavy). scipy.optimize.minimize with L-BFGS-B then polishes it to a continuous solution inside the bounds. The variables are scaled to mm/s so that both parameters have a similar magnitude, which helps the optimizer’s convergence.
Three answers are extracted:
- B (free optimum): the global minimum of $J$.
- C (constrained optimum): the minimum of $J$ among grid points with $P_c\le10^{-5}$, found by masking infeasible points with
np.inf. - D (minimum fuel): for every lead time, the smallest $\Delta v$ that keeps $P_c$ below the limit for all larger burns as well. This is computed with a reversed cumulative logical AND, which avoids being fooled by isolated dips in the wavy $P_c$ surface.
Section 5 and 6: Verification and reporting
The analytic result of every strategy is compared with the importance-sampling Monte Carlo estimate, and everything is printed in one table.
Section 7: A single combined figure
Six panels are drawn in one figure with layout="constrained" so that colorbars and 3D axes never overlap. Panels (a) and (b) are 3D surfaces, and (c) to (f) are 2D analyses. plt.show() is called only once.
Results
Console output
Orbital period : 94.62 min Grid evaluation (301x301) : 10.98 ms (vectorised) Loop vs vectorised (100x100): 231.2 ms vs 0.74 ms -> x314 faster Strategy dv [mm/s] Lead [orb] Pc (analytic) Pc (MC) Cost [k$] ---------------------------------------------------------------------------------------------- A: Do nothing 0.00 - 3.471e-03 3.464e-03 1735.5 B: Cost-optimal (free) 11.50 3.658 2.330e-05 2.346e-05 126.7 C: Cost-optimal (Pc <= 1e-5) 12.60 3.663 8.693e-06 8.762e-06 130.3 D: Min-fuel (Pc <= 1e-5) 12.60 3.575 9.688e-06 - 130.8 MC standard error (B): 2.55e-09
In our run, the vectorized evaluation was more than 100 times faster than the double loop. The table shows four strategies:
| Strategy | $\Delta v$ [mm/s] | Lead time [orbits] | $P_c$ | Expected cost [k$] |
|---|---|---|---|---|
| A: Do nothing | 0 | – | $3.5\times10^{-3}$ | 1735.5 |
| B: Cost-optimal | 11.50 | 3.658 | $2.3\times10^{-5}$ | 126.7 |
| C: Cost-optimal with $P_c\le10^{-5}$ | 12.60 | 3.663 | $8.7\times10^{-6}$ | 130.3 |
| D: Minimum fuel with $P_c\le10^{-5}$ | 12.60 | 3.575 | $9.7\times10^{-6}$ | 130.8 |
The most striking number is the drop from 1.7 million dollars of expected loss to about 127 thousand dollars, a reduction of more than 90 percent, achieved with a burn of only 1.15 cm/s. The Monte Carlo values agree with the analytic ones to within about one percent in every case, which validates the small-radius approximation.
The pure cost optimum (B) ends at $P_c=2.3\times10^{-5}$, slightly above the safety rule. Under the constraint, the answer moves to a somewhat larger burn of 12.6 mm/s (C), and it costs only about 3.6 thousand dollars more in expected loss. That is a cheap price for satisfying a hard operational threshold. The minimum-fuel solution (D) turns out to be almost identical to C, which tells us that in this scenario the fuel term dominates the constrained decision.
Combined figure
Reading the Graphs
(a) Expected cost surface (3D). The cost surface is high and wavy at the small-$\Delta v$ side, where the collision risk dominates, and it flattens into a broad valley toward larger $\Delta v$. There the fuel term takes over and slowly lifts the surface again. The red star marks the optimum at a lead time near 3.7 orbits. Notice the ridges along the lead-time axis: they are the fingerprints of the $\sin nt$ and $\cos nt$ terms, and they show that some lead times are much better than their neighbors.
(b) Collision probability surface (3D). On the logarithmic scale, $P_c$ falls by many orders of magnitude within only a few tens of mm/s. The staircase-like shelves reveal how the maneuver first moves the miss vector out of the dense core of the covariance ellipse, then out of the 1σ region, and finally into its tail. The floor is clipped at $10^{-12}$ for readability.
(c) $P_c$ map. This is the same surface seen from above. The white curve is the safety limit $P_c=10^{-5}$: everything to its right is acceptable. The curve is not a straight line, because it bends at roughly 1 and 2 orbits, where the burn geometry becomes inefficient. The three markers (B, C, and D) crowd together near the top of the map. This means the best strategy is to act as early as the warning allows.
(d) Encounter plane. Here we look at the geometry directly. The red ellipses (do nothing) sit right on top of the tiny black hard-body disk, which is exactly why $P_c$ is so high. After the maneuver, the blue and green ellipses have been moved about 700 m along-track. Their larger size, compared with the red ones, is the price of execution error amplified over almost four orbits, but the center of the distribution is now far from the target, so the probability mass overlapping the disk is negligible.
(e) Best cost per lead time. For every lead time we take the best $\Delta v$ and split the cost into its two parts. The collision-risk component (red) is always tiny compared with the fuel component (blue), which shows that once the optimizer has acted, most of the remaining cost is the fuel bill. The total (black) decreases with lead time overall, with oscillations caused by orbital geometry, and the gray line, the cost of doing nothing, is more than an order of magnitude above every point.
(f) Required $\Delta v$. The minimum burn needed to meet the safety limit falls from almost 60 mm/s with a quarter-orbit notice to about 12 mm/s with nearly four orbits. The local bump around one orbit is a well-known feature: after a full revolution the radial term of the burn returns to zero, so the maneuver is less effective at that phase. The orange square and the green circle sit at the bottom of the curve.
Conclusions
- A tiny along-track burn of about 1 cm/s, executed roughly four orbits before closest approach, lowers the collision probability from $3.5\times10^{-3}$ to below $10^{-5}$.
- Lead time is the most valuable resource. Every extra orbit of warning reduces the needed $\Delta v$ substantially, but the periodic geometry means the relationship is not monotonic, so a grid search followed by local refinement is safer than a blind gradient method.
- Enforcing a hard risk threshold costs only a few thousand dollars of expected loss here, which makes such a rule easy to justify.
- Vectorizing the model with NumPy broadcasting gave a speed-up of about two orders of magnitude and turned the whole analysis into a matter of milliseconds.
The same framework extends naturally to more realistic problems: full 3D encounter geometry, multiple burn options, non-Gaussian covariances, or several simultaneous conjunctions sharing one fuel budget.
