Space-weather monitoring satellites face a scheduling problem that looks deceptively simple: point the instrument at the most scientifically important target at every moment. In practice it is anything but simple. A single-instrument satellite (or a satellite sharing one pointing axis across several sensors) can only stare at one target at a time — a flaring active region on the Sun, the L1 point watching for an incoming CME, the radiation belts, the polar ionosphere, or the magnetopause boundary. Each target’s importance rises and falls over the day as forecasts update, and every time the satellite re-points, it pays a “slew” penalty in time and momentum-wheel fuel.
This is a textbook combinatorial optimization problem, and it fits neatly into a Mixed-Integer Linear Program (MILP). Below we formulate it, solve it in Google Colaboratory with scipy.optimize.milp, and visualize the result — including the slew path traced across the sky.
Problem Formulation
Let there be $T$ discrete observation slots (hours) and $N$ candidate targets. The binary decision variable $x_{t,i}$ equals $1$ if the satellite observes target $i$ during slot $t$:
$$x_{t,i}\in{0,1}, \quad t=0,\dots,T-1,\ \ i=0,\dots,N-1$$
To penalize re-pointing, we introduce a switching-indicator variable $y_{t,i,j}$ that turns on whenever the satellite moves from target $i$ (at slot $t-1$) to a different target $j$ (at slot $t$):
$$y_{t,i,j}\ge 0, \quad t=1,\dots,T-1,\ \ i\neq j$$
The objective maximizes accumulated scientific value minus total slew cost:
$$\max \ \sum_{t=0}^{T-1}\sum_{i=0}^{N-1} w_i(t),x_{t,i} ;-; \sum_{t=1}^{T-1}\sum_{i\neq j} c_{ij},y_{t,i,j}$$
subject to:
$$\sum_{i=0}^{N-1} x_{t,i} = 1 \qquad \forall t \quad \text{(observe exactly one target per slot)}$$
$$y_{t,i,j} \ge x_{t-1,i} + x_{t,j} - 1 \qquad \forall t\ge1,\ i\neq j \quad \text{(switch indicator linking)}$$
The time-varying scientific value of target $i$ is modeled as a baseline plus a Gaussian urgency bump centered on the forecast peak time $\tau_i$ (e.g., predicted flare or CME arrival):
$$w_i(t) = b_i + a_i \exp!\left(-\frac{(t-\tau_i)^2}{2\sigma_i^2}\right)$$
The slew cost between two targets is proportional to the great-circle angle between their boresight unit vectors $\hat{u}_i$:
$$c_{ij} = \kappa \arccos(\hat{u}_i\cdot \hat{u}_j)$$
Full Python Implementation (Google Colaboratory)
1 | import numpy as np |
Code Walkthrough
Section 1 — Problem setup. Six representative space-weather targets are defined with a baseline importance, an urgency amplitude, a forecast peak hour, and a spread. The value matrix w[t, i] is built in one vectorized NumPy expression rather than a nested loop, so it costs virtually nothing even for much larger $T \times N$ grids. Each target also gets a boresight direction on the unit sphere; the pairwise angular distance matrix angle (and thus switch_cost) is computed with a single matrix multiplication.
Section 2 — Variable indexing. Instead of a 3-D array of milp variables, everything is flattened into one long vector, since scipy.optimize.milp expects flat variable arrays. x_idx and y_idx are small helper closures that map $(t,i)$ and $(t,i,j)$ triples to flat positions.
Section 3 — Objective. The value terms go in with a negative sign because milp minimizes; the switching cost terms go in with a positive sign so that minimizing the total effectively maximizes value while penalizing slews.
Section 4 — Constraints. Constraint (a) forces exactly one target per slot. Constraint (b) is the standard linear-relaxation trick for an AND of two binaries: $y_{t,i,j}\ge x_{t-1,i}+x_{t,j}-1$. Because the cost coefficient on $y$ is strictly positive and the solver minimizes, $y$ is always pushed down to exactly $\max(0, x_{t-1,i}+x_{t,j}-1)$ at the optimum — which is automatically $0$ or $1$ whenever $x$ is binary. That means $y$ never needs to be declared integer, cutting the number of integer variables from 834 down to 144 in this example without losing any correctness.
Section 5 — Bounds & integrality. Only the x block is marked integer (integrality[:n_x] = 1); the y block stays continuous in $[0,1]$, relying on the argument above.
Section 6 — Solve & report. scipy.optimize.milp calls the HiGHS solver under the hood. The chosen target per slot is recovered via argmax over x_sol, and total value / number of slews / total slew cost are all recomputed directly from the schedule array — not from the y variables — so the reported numbers are correct regardless of any floating-point residue in y.
Section 7 — Visualization. A single Figure with four subplots is built in one plt.show() call: a value heatmap with the chosen schedule overlaid, a 3-D value surface with the schedule marked on top, a 3-D slew path traced on the unit sphere connecting the chosen targets in time order, and a Gantt-style bar showing the full day’s pointing plan.
============================================================ SATELLITE SPACE-WEATHER OBSERVATION SCHEDULE ============================================================ Solver status : Optimization terminated successfully. (HiGHS Status 7: Optimal) Solve time : 0.125 s Total science value : 18.203 Number of slews : 4 Total slew cost : 1.863 ------------------------------------------------------------ Hour 00:00 -> Polar Ionospheric TEC Scan Hour 01:00 -> Polar Ionospheric TEC Scan Hour 02:00 -> Polar Ionospheric TEC Scan Hour 03:00 -> AR3536 Flare Watch Hour 04:00 -> AR3536 Flare Watch Hour 05:00 -> AR3536 Flare Watch Hour 06:00 -> AR3536 Flare Watch Hour 07:00 -> L1 CME Arrival Monitor Hour 08:00 -> L1 CME Arrival Monitor Hour 09:00 -> L1 CME Arrival Monitor Hour 10:00 -> L1 CME Arrival Monitor Hour 11:00 -> L1 CME Arrival Monitor Hour 12:00 -> AR3541 Flare Watch Hour 13:00 -> AR3541 Flare Watch Hour 14:00 -> AR3541 Flare Watch Hour 15:00 -> AR3541 Flare Watch Hour 16:00 -> AR3541 Flare Watch Hour 17:00 -> GEO Radiation Belt Monitor Hour 18:00 -> GEO Radiation Belt Monitor Hour 19:00 -> GEO Radiation Belt Monitor Hour 20:00 -> GEO Radiation Belt Monitor Hour 21:00 -> GEO Radiation Belt Monitor Hour 22:00 -> GEO Radiation Belt Monitor Hour 23:00 -> GEO Radiation Belt Monitor ============================================================
Reading the Schedule
The console output lists, hour by hour, which target the satellite is pointed at, followed by the total accumulated scientific value, the number of slews performed over the day, and the cumulative slew cost. You should see the schedule cluster tightly around each target’s Gaussian peak hour $\tau_i$ — the satellite “camps” on a target through its urgency window rather than darting back and forth, because every switch away and back costs $c_{ij}$ twice. Where two targets’ urgency windows overlap, the one with the higher combined baseline-plus-amplitude value wins the contested slots, and the loser is picked up just before or after its own peak instead.
Visualizing the Optimized Schedule
The top-left heatmap shows brighter bands wherever a target’s Gaussian bump is active, with cyan dots marking exactly which slot the optimizer chose — these dots should sit inside or very near the brightest region of whichever target is active. The top-right 3-D surface is the same value landscape lifted into three dimensions, with red markers riding along the scheduled ridge line. The bottom-left sphere plot is the most physically intuitive view: each colored dot is a target’s boresight direction, and the cyan polyline is the literal path the satellite’s optics trace across the sky over 24 hours — long straight segments mean big, costly slews, while short segments mean the satellite stayed close to its previous pointing. The bottom-right Gantt bar gives the same schedule as a flat timeline, useful for quickly reading off exactly when each target is being watched.
Why This Formulation Scales Well
Two design choices keep this fast even as the constellation grows. First, the switching-indicator variables are relaxed to continuous rather than binary — a mathematically exact simplification given the positive-cost, greater-than-or-equal constraint structure — which removes the vast majority of integer variables from the branch-and-bound search. Second, the constraint matrix is built directly as a sparse csr_matrix from coordinate lists rather than as a dense array, so memory and construction time scale with the number of nonzero entries rather than $T^2N^2$. Together these mean the same code handles a week-long, twenty-target schedule about as comfortably as the one-day, six-target example shown here, without needing a separate “slow” and “fast” version.
Extensions
The same skeleton generalizes naturally: add a second and third satellite by indexing $x$ over satellites as well as time and target, add a downlink data-volume budget as an extra linear constraint, or replace the fixed Gaussian urgency bumps with a live Kp-index or flare-probability forecast feed. The optimization structure — one-target-per-slot assignment plus a linearized switching penalty — stays exactly the same.
