Solar-wind monitors such as SOHO, ACE, Wind and DSCOVR sit near the Sun–Earth L1 point, about 1.5 million km sunward of Earth. Solar wind at 400–500 km/s takes roughly an hour to reach Earth from there, so an instrument at L1 gives a short but valuable warning of geomagnetic storms. This post designs the orbit of such a satellite as a concrete optimization problem and solves it in Python.
1. The example problem
A spacecraft flies a halo orbit around L1 in the Sun–Earth system. Two things must be chosen:
- $A_z$, the z-amplitude of the halo orbit.
- $\tau$, the interval between station-keeping maneuvers.
The mission requirements are:
- Solar radio interference. The Sun–Earth–spacecraft (SEV) angle $\alpha$ must never fall below $4.5^\circ$. Below that, solar radio noise swamps the downlink.
- Antenna pointing. The angle must never exceed $28^\circ$.
The cost to minimize is the total $\Delta V$ over the mission:
$$
\min_{A_z,;\tau}; J = \Delta V_{\rm ins}^{+}(A_z) + N_{\rm yr},\Delta V_{\rm sk}(A_z,\tau)
$$
$$
\text{s.t.}\quad \alpha_{\min}(A_z)\ge 4.5^\circ,\quad \alpha_{\max}(A_z)\le 28^\circ,\quad 10^5\le A_z\le 4\times10^5\ \text{km},\quad 10\le\tau\le 90\ \text{day}
$$
| Symbol | Meaning | Value |
|---|---|---|
| $\mu$ | Sun–Earth mass ratio | $3.00348\times10^{-6}$ |
| $N_{\rm yr}$ | Mission lifetime | 5 yr |
| $\sigma_v$ | Dispersion growth per maneuver cycle | 0.02 m/s |
| $\delta_v$ | Fixed cost per maneuver | 0.05 m/s |
The coefficients $\sigma_v$ and $\delta_v$ are illustrative values for this example.
2. Mathematical model
2.1 Circular restricted three-body problem
In nondimensional units, one length unit is 1 AU and one time unit is $1/n$. The frame rotates with the Sun–Earth line, with the Sun at $(-\mu,0,0)$ and the Earth at $(1-\mu,0,0)$:
$$
\ddot x-2\dot y=\Omega_x,\qquad \ddot y+2\dot x=\Omega_y,\qquad \ddot z=\Omega_z
$$
$$
\Omega=\frac{x^2+y^2}{2}+\frac{1-\mu}{r_1}+\frac{\mu}{r_2}+\frac{\mu(1-\mu)}{2},\quad
r_1=\sqrt{(x+\mu)^2+y^2+z^2},\quad r_2=\sqrt{(x-1+\mu)^2+y^2+z^2}
$$
The Jacobi constant is conserved along any trajectory. This makes it a good check on numerical accuracy:
$$
C = 2\Omega - (\dot x^2+\dot y^2+\dot z^2)
$$
The collinear point L1 is the root of $\Omega_x(x,0,0)=0$ between the Sun and the Earth.
2.2 State transition matrix
Write the state as $\mathbf s=[x,y,z,\dot x,\dot y,\dot z]^\top$. Its sensitivity to the initial state obeys
$$
\dot\Phi = A(t),\Phi,\qquad
A=\begin{bmatrix}0 & I\ \Omega_{\mathbf{rr}} & \Xi\end{bmatrix},\qquad
\Xi=\begin{bmatrix}0&2&0\-2&0&0\0&0&0\end{bmatrix},\qquad \Phi(0)=I
$$
where $\Omega_{\mathbf{rr}}$ is the Hessian of $\Omega$ with respect to position.
2.3 Halo orbit by differential correction
A halo orbit is symmetric about the $xz$-plane. Start from
$$
\mathbf s_0=[x_0,,0,,z_0,,0,,\dot y_0,,0]^\top,\qquad z_0=A_z
$$
and require a perpendicular crossing of the $y=0$ plane at the half period:
$$
\dot x(T/2)=0,\qquad \dot z(T/2)=0
$$
The unknowns are $(x_0,\dot y_0)$. Newton’s method uses the STM $\Phi$ evaluated at the crossing. The crossing time itself moves when the initial state changes, which adds a correction term:
$$
\begin{bmatrix}\delta x_0\ \delta\dot y_0\end{bmatrix}
=-J^{-1}\begin{bmatrix}\dot x_f\ \dot z_f\end{bmatrix},\qquad
J=\begin{bmatrix}
\Phi_{41}-\dfrac{\ddot x_f}{\dot y_f}\Phi_{21} & \Phi_{45}-\dfrac{\ddot x_f}{\dot y_f}\Phi_{25}\[3mm]
\Phi_{61}-\dfrac{\ddot z_f}{\dot y_f}\Phi_{21} & \Phi_{65}-\dfrac{\ddot z_f}{\dot y_f}\Phi_{25}
\end{bmatrix}
$$
2.4 Instability of the orbit
A halo orbit around L1 is unstable. The monodromy matrix $M=\Phi(T)$ has one real eigenvalue $\lambda_u\gg1$. Its growth rate is
$$
s=\frac{\ln\lambda_u}{T}
$$
and the e-folding time of an orbit error is $1/s$.
2.5 The SEV angle constraint
The Earth-to-spacecraft vector is $\boldsymbol\rho=(x-1+\mu,;y,;z)$. The Earth-to-Sun direction is $(-1,0,0)$, so
$$
\alpha=\arccos!\left(\frac{-(x-1+\mu)}{|\boldsymbol\rho|}\right)
$$
2.6 Cost model
Insertion cost. Arriving with the L1 energy level $C_{L1}$ and moving onto a halo orbit requires raising the energy by $\Delta C=C_{L1}-C_{\rm halo}$. A burn at the highest orbital speed $v_{\max}$ changes $C$ by $-2v,\delta v$, which gives a first-order estimate. The result is expressed as an increment over the smallest orbit in the search range, $A_{z,\rm lo}=10^5$ km. A constant offset does not change the optimum:
$$
\Delta V_{\rm ins}^{+}(A_z)=\left[\frac{\Delta C(A_z)}{2v_{\max}(A_z)}-\frac{\Delta C(A_{z,\rm lo})}{2v_{\max}(A_{z,\rm lo})}\right]V_U
$$
Station keeping. An error in the unstable mode grows as $e^{s\tau}$ between maneuvers. The per-year cost is
$$
\Delta V_{\rm sk}(A_z,\tau)=\frac{365.25}{\tau}\Big(\sigma_v,e^{s(A_z),\tau}+\delta_v\Big)
$$
Short intervals waste fixed maneuver costs. Long intervals let the instability amplify errors. Setting $\partial\Delta V_{\rm sk}/\partial\tau=0$ gives the optimal interval in closed form:
$$
\sigma_v,e^{s\tau},(s\tau-1)=\delta_v
$$
2.7 Optimality conditions
Form the Lagrangian with multipliers $\nu_1,\nu_2\ge0$ for the two angle constraints:
$$
\mathcal L=J-\nu_1\big(\alpha_{\min}-4.5^\circ\big)-\nu_2\big(28^\circ-\alpha_{\max}\big)
$$
At an optimum with only the first constraint active:
$$
\frac{\partial J}{\partial\tau}=0,\qquad
\frac{\partial J}{\partial A_z}=\nu_1\frac{\partial\alpha_{\min}}{\partial A_z}
$$
The multiplier $\nu_1$ has a direct meaning. It is the extra $\Delta V$ paid for each additional degree of exclusion angle.
3. Source code
1 | import time |
4. Code walkthrough
Section 1: constants and L1
All dynamics run in nondimensional units, with 1 AU as the length unit and $1/n$ as the time unit. Then $\mu$ is the only physical parameter, and one revolution of the rotating frame takes $2\pi$ time units. The constants TU_DAY and VU_MS convert results back to days and m/s.
lagrange_l1() finds L1 by Brent’s method on the axial force balance. pseudo_potential() is $\Omega$, and C_L1 is the Jacobi constant of a spacecraft at rest at L1. It serves as the energy reference for the insertion cost.
Section 2: dynamics and differential correction
rhs_with_stm. The right-hand side packs the 6 state variables and the flattened $6\times6$ STM into one 42-element vector. It uses scalar math.sqrt rather than NumPy ufuncs for the small operations, because NumPy call overhead dominates at this scale. This function is called hundreds of thousands of times, so the saving is real.
cross_y_down. This event function detects the half-period plane crossing. With direction = -1 it fires only when $y$ passes through zero downward. The initial state also has $y=0$, but $\dot y_0>0$ there, so it is ignored automatically and no manual “skip the first step” logic is needed.
correct_halo. This is the Newton iteration from Section 2.3. Because the analytic Jacobian comes from the STM, each iteration needs only one integration. A finite-difference Jacobian would need three. It also converges quadratically, typically in a handful of iterations.
linear_guess. Around L1 the linearized dynamics have an in-plane frequency $\omega_p$ and an amplitude ratio $k$. Together they give a starting guess for $x_0$ and $\dot y_0$ from a chosen in-plane amplitude. The seed loop tries several amplitudes so that one failed guess cannot stop the script.
evaluate_orbit. This integrates one full period with dense output at 3001 points and extracts everything the optimizer needs:
- the monodromy eigenvalue $\lambda_u$ and growth rate $s$,
- the SEV angle history $\alpha(t)$,
- the Jacobi constant history, whose drift $\Delta C$ verifies the integration,
- the insertion-cost estimate.
Section 3: continuation
Solving each $A_z$ from scratch would be wasteful. The family is computed by continuation: the first orbit uses the linear guess, and every later orbit uses a secant predictor,
$$
\mathbf p_{i+1}^{\rm guess}=2\mathbf p_i-\mathbf p_{i-1}
$$
Newton then converges in one or two iterations. After the loop, the tabulated quantities $\alpha_{\min}$, $\alpha_{\max}$, $s$, $\Delta V_{\rm ins}$, $x_0$ and $\dot y_0$ are turned into cubic splines over $A_z$.
Section 4: cost model
dv_keep_per_year and total_cost are pure NumPy expressions with no ODE solves inside. tau_optimal solves $\sigma_v e^{x}(x-1)=\delta_v$ by a vectorized Newton iteration, where $x=s\tau$. The function is convex and increasing for $x>1$, so starting at $x=2$ converges monotonically.
Section 5: optimization
The optimization runs in two stages:
- Brute-force grid. A $301\times301$ grid over $(A_z,\tau)$ is evaluated in a single vectorized call. Infeasible points are set to $+\infty$, and the minimum gives a starting point.
- SLSQP refinement. Variables are scaled to order one with
SCALE, which matters for SLSQP conditioning. The two SEV constraints enter as inequality constraints.
Section 6: verification
Splines are approximations, so the answer is verified against the full nonlinear model. true_orbit() corrects a halo at any $A_z$ and uses the splines of $x_0$ and $\dot y_0$ as initial guesses, so it converges almost immediately. Brent’s method then finds the exact $A_z$ at which $\alpha_{\min}=4.5^\circ$ and $\alpha_{\max}=28^\circ$. The final cost is recomputed with the true growth rate at that orbit.
Speed-up strategy
A direct approach would call the nonlinear halo solver at every point of a $301\times301$ grid. At about 0.04 s per orbit, that is roughly an hour. Four measures keep the whole script at a few seconds:
- Analytic Newton Jacobian from the STM, so one integration per iteration.
- Continuation with a secant predictor.
- Splines of the expensive quantities, so the optimizer never calls an ODE solver.
- Vectorized Newton for $\tau^*$ and a single vectorized grid evaluation.
Section 8: the combined figure
A single GridSpec of 2×3 panels holds three 3-D plots and three 2-D plots, saved as one PNG. Panel (a) draws the orbit with gray wall projections. Panel (c) uses an azimuthal equidistant projection: the radius is the SEV angle $\alpha$ and the polar angle is $\arctan(z/y)$. This makes the $4.5^\circ$ exclusion zone an exact circle.
5. Execution results
Console output
========================================================================
Sun-Earth L1 halo orbit design for a space-weather satellite
========================================================================
L1 distance from Earth : 1491550.9 km
Time unit / velocity unit : 58.1324 day / 29784.74 m/s
Family computed (n = 31) : 7.41 s
------------------------------------------------------------------------
Az [km] T [day] lambda a_min[deg] a_max[deg] dV_ins+[m/s]
100000 177.880 1744.3 3.441 25.621 0.000
150000 177.856 1718.0 5.156 25.945 3.788
200000 177.821 1681.9 6.865 26.392 8.909
250000 177.777 1636.7 8.567 26.953 15.203
300000 177.721 1583.3 10.259 27.621 22.496
350000 177.654 1522.5 11.940 28.389 30.604
400000 177.575 1455.5 13.610 29.246 39.354
------------------------------------------------------------------------
Feasible Az range (exact, nonlinear model)
lower bound (a_min = 4.5 deg) : 130848.1 km
upper bound (a_max = 28.0 deg) : 325444.1 km
------------------------------------------------------------------------
Optimisation results
grid search : Az = 131000.0 km, tau = 36.67 day, J = 9.3064 m/s
SLSQP : Az = 130848.1 km, tau = 36.68 day, J = 9.2944 m/s (success = True)
verified : Az = 130848.1 km, tau = 36.68 day, J = 9.2944 m/s
------------------------------------------------------------------------
Optimal orbit (nonlinear model)
period : 177.8660 day
unstable eigenvalue : 1729.26 (e-folding time 23.86 day)
SEV angle range : 4.5000 - 25.8063 deg
initial state (x0, vy0): 0.988883314, 0.008921282
closure error : 1.393e-11
Jacobi constant drift : 1.776e-15
insertion increment : 2.172 m/s
station keeping : 1.425 m/s/yr x 5 yr = 7.123 m/s
total cost J : 9.294 m/s
------------------------------------------------------------------------
Total computation time : 9.48 s
========================================================================
Result image
6. Reading the results
The optimum
The optimizer finds:
- $A_z^*\approx 130{,}848$ km
- $\tau^*\approx 36.7$ days
- $J^*\approx 9.29$ m/s, made up of a 2.17 m/s insertion increment and 1.42 m/s per year of station keeping (7.12 m/s over five years)
The optimal orbit has a period of 177.87 days. Its unstable eigenvalue is about 1729, which corresponds to an error e-folding time of about 23.9 days. The closure error is of order $10^{-11}$ and the Jacobi constant drift is of order $10^{-15}$. The orbit is therefore a genuine periodic halo orbit, and the integration is accurate.
The grid search lands within one grid cell of this answer, at 131,000 km and 36.67 days. SLSQP and the independent nonlinear verification agree to all printed digits.
Panel (a): the optimal orbit in 3-D
The halo orbit is a large tilted loop around L1. Relative to L1 it reaches about −171 to +245 thousand km in $x$ (toward the Sun and toward the Earth), ±667 thousand km in $y$, and −105 to +131 thousand km in $z$. The gray curves are its projections onto the three walls. The dashed line is the Sun–Earth line, and the spacecraft distance from Earth stays between roughly 1.25 and 1.67 million km. The loop is far larger in $y$ than in $z$, which is why the orbit looks flattened.
Panel (b): the halo family
The colored curves show the family for $A_z$ from 100 to 400 thousand km, with the optimal orbit in red. Larger $A_z$ tilts the loop out of the ecliptic plane and enlarges it slightly. The red orbit sits at the small-amplitude end of the family. This is what the constraint analysis below predicts.
Panel (c): sky view from Earth
This is the view of the orbit as seen from Earth, with the Sun at the origin. The pink disk is the $4.5^\circ$ exclusion zone. The orange dashed orbit ($A_z=100{,}000$ km) dips into the zone, so it is infeasible. The red orbit at the optimum just grazes the top edge of the circle. The green orbit ($A_z=300{,}000$ km) keeps a wide margin, but it pays for that margin in insertion cost.
The bound is easy to check by hand. The tightest approach to the Sun occurs where the spacecraft is at the Sun-side extreme of the orbit with $z=A_z$. That point is about 171,000 km sunward of L1, so it is about 1,662,600 km from Earth. Then
$$
A_z \approx 1{,}662{,}600\ \text{km}\times\tan 4.5^\circ\approx 130{,}850\ \text{km}
$$
which matches the numerical result.
Panel (d): SEV angle over one revolution
The angle oscillates twice per revolution. The deeper minimum occurs at the Sun-side crossing, where $z$ is at its maximum, and this is what sets the lower bound on $A_z$. The second minimum, at the Earth-side crossing where $z=-105{,}000$ km, is slightly higher at about $4.8^\circ$. The orange curve ($A_z=100{,}000$ km) dips into the lower red band, with a minimum of $3.44^\circ$. The purple curve ($A_z=400{,}000$ km) rises to $29.2^\circ$ and pokes into the upper red band, violating the antenna limit. The optimal red curve touches the $4.5^\circ$ line exactly. The feasible range of $A_z$ is therefore about 130.8 to 325.4 thousand km.
Panel (e): the trade-off
The blue curve shows that the insertion increment grows rapidly and convexly with $A_z$: from 0 at 100,000 km to about 39 m/s at 400,000 km. The green curve shows that the five-year station-keeping cost is almost flat, falling only from about 7.13 to 6.97 m/s. Larger orbits are slightly less unstable, with $\lambda_u$ dropping from about 1744 to 1456, but the benefit is small. The black total cost is therefore dominated by the insertion term and increases with $A_z$. Without the SEV constraint, the best orbit would be the smallest one. The exclusion angle is what pushes the optimum up to the boundary, and the optimum sits exactly at the edge of the pink infeasible band.
The multiplier from Section 2.7 can be read off numerically. Near the optimum, $\partial J/\partial A_z\approx0.079$ m/s per 1000 km and $\partial\alpha_{\min}/\partial A_z\approx0.034^\circ$ per 1000 km. Their ratio gives $\nu_1\approx2.3$ m/s per degree. Each extra degree of solar exclusion angle costs about 2.3 m/s over the mission.
Panel (f): cost surface
This 3-D surface shows $J(A_z,\tau)$ over the feasible region. Along the $\tau$ direction it forms a clear valley. At $A_z^*$, the cost at $\tau=10,,20,,30,,36.7,,50,,70,,90$ days is about 16.9, 11.0, 9.5, 9.29, 9.9, 13.3, 20.8 m/s. Short intervals pay the fixed maneuver cost too often. Long intervals let the instability amplify errors exponentially. The black curve on the surface is the analytic $\tau^*(A_z)$, which stays close to 37 days across the whole family because $s$ varies only slightly. Along the $A_z$ direction the surface rises steadily, which is the insertion term at work. The red star sits on the boundary of the feasible region, exactly where the constraint is active.
7. Takeaways
- The optimal design is a compromise between two effects: radio interference forces the halo to be large, and insertion energy forces it to be small.
- The station-keeping interval has a genuine interior optimum of about 37 days, set by balancing fixed maneuver cost against exponential error growth with an e-folding time near 24 days.
- Every closed-form check agrees with the numerics: the exclusion-angle geometry for $A_z^*$, and the stationarity condition for $\tau^*$.
- A spline surrogate built from a handful of nonlinear orbit solutions turns a problem that would take about an hour by brute force into a few seconds, without giving up verification against the full model.
