A Concrete Example Solved with Python
Introduction
A solar flare can disturb HF radio within minutes. A coronal mass ejection (CME) can trigger a geomagnetic storm a day or two later, degrading GNSS positioning, stressing power grids, and endangering satellites. Forecasters cannot watch everything all the time. Coronagraph time, magnetograph scans, ground magnetometer sampling, ionosonde sounding, GNSS-TEC processing, and L1 solar wind monitoring all have finite budgets.
This article poses a concrete allocation problem, derives its exact solution with Lagrange multipliers, and solves it in Python. Along the way we compare a generic solver with a fast semi-analytic algorithm and visualize everything in a single multi-panel figure that includes 3D graphs.
Problem Setting
We plan the next $T = 24$ hours in one-hour slots and have $I = 6$ observation resources:
| Index | Resource | Role |
|---|---|---|
| 1 | Coronagraph | CME detection and speed estimation |
| 2 | Magnetograph | Active-region magnetic-field evolution |
| 3 | Ground magnetometers | Geomagnetic disturbance monitoring |
| 4 | Ionosondes | Ionospheric state |
| 5 | GNSS-TEC receivers | Total electron content mapping |
| 6 | L1 solar wind monitor | In-situ solar wind measurement |
Let $x_{t,i} \ge 0$ be the observation effort given to resource $i$ in slot $t$. Each resource has:
- a value weight $w_i$,
- a saturation rate $a_i$ (how quickly extra effort stops adding information),
- a unit cost $c_i$.
Each slot has a risk level $r_t$, meaning how much an accurate observation is worth at that hour. We model one flare-related risk pulse around hour 6 and one CME-arrival pulse around hour 15:
$$
r_t = r_0 + s\left[A_1 \exp!\left(-\frac{(t-t_1)^2}{2\sigma_1^2}\right) + A_2 \exp!\left(-\frac{(t-t_2)^2}{2\sigma_2^2}\right)\right]
$$
with $r_0 = 0.15$, $s = 1$, $(A_1, t_1, \sigma_1) = (1.0, 6, 1.5)$, and $(A_2, t_2, \sigma_2) = (1.6, 15, 2.0)$.
Mathematical Formulation
The information gained from resource $i$ in slot $t$ saturates exponentially with effort. Maximizing the total risk-weighted information gain gives
$$
\max_{x}; U(x) = \sum_{t=1}^{T}\sum_{i=1}^{I} r_t, w_i \left(1 - e^{-a_i x_{t,i}}\right)
$$
subject to a total cost budget and per-slot effort caps:
$$
\sum_{t=1}^{T}\sum_{i=1}^{I} c_i, x_{t,i} \le B, \qquad 0 \le x_{t,i} \le x_{\max}.
$$
Each term of $U$ is concave in $x_{t,i}$ and all constraints are linear, so this is a convex optimization problem. The KKT conditions are therefore both necessary and sufficient.
Closed-form structure via the Lagrangian
Introduce a multiplier $\lambda \ge 0$ for the budget constraint:
$$
\mathcal{L}(x,\lambda) = \sum_{t,i} r_t w_i\left(1-e^{-a_i x_{t,i}}\right) - \lambda\left(\sum_{t,i} c_i x_{t,i} - B\right).
$$
Setting the derivative with respect to an interior $x_{t,i}$ to zero gives
$$
r_t, w_i, a_i, e^{-a_i x_{t,i}} = \lambda, c_i .
$$
In words, at the optimum the marginal information gain per unit cost is equalized across all active (slot, resource) pairs, and that common value is $\lambda$. Solving for $x$ and applying the bounds gives
A pair $(t,i)$ receives effort only if $r_t w_i a_i / c_i > \lambda$. The total cost $\sum c_i x^{*}_{t,i}(\lambda)$ is monotonically decreasing in $\lambda$, so the $\lambda$ that exactly exhausts the budget can be found by bisection. This is a water-filling type solution.
Source Code
1 | import time |
Detailed Code Walkthrough
1. Problem data. The arrays w, a, and c hold $w_i$, $a_i$, and $c_i$. The L1 monitor has the highest weight and saturation rate, since in-situ measurements are the most informative, but also the highest cost. The GNSS-TEC receivers are cheap but lower in value. risk_profile builds $r_t$ as a constant floor plus two Gaussian pulses, and the scale argument plays the role of $s$ so that storm intensity can be varied later.
2. Objective and gradient. utility evaluates $U(x)$ with fully vectorized broadcasting over the $T \times I$ matrix. gradient returns
$$
\frac{\partial U}{\partial x_{t,i}} = r_t, w_i, a_i, e^{-a_i x_{t,i}},
$$
which the generic solver needs for fast convergence.
3. Generic solver. solve_slsqp flattens the $24 \times 6$ matrix into 144 variables. np.tile(c, T) produces a cost vector that matches row-major flattening, so index $t \cdot I + i$ maps to $c_i$. The budget is passed as a linear inequality constraint with its exact Jacobian, and the bounds enforce $0 \le x \le x_{\max}$. After convergence, the solution is clipped to the bounds and rescaled if it exceeds the budget by rounding error, which guarantees feasibility.
4. Fast solver. allocation_for_lambda implements
in one vectorized expression. np.outer(r, w*a/c) builds the whole matrix of $r_t w_i a_i / c_i$, and np.maximum(arg, 1e-300) protects the logarithm. solve_fast then bisects on $\log \lambda$. The interval $[10^{-9}, 10^{3}]$ spans 12 orders of magnitude, and 100 halvings shrink it far below floating-point resolution. Because total cost is decreasing in $\lambda$, the invariant “lower end infeasible, upper end feasible” is maintained. If even $\lambda = 10^{-9}$ fits within the budget, every effort is already saturated and that allocation is returned directly.
5. Comparison and diagnostics. Two intuitive baselines are computed alongside the optimum: uniform allocation and allocation proportional to risk. Both spend exactly the budget $B$, so the comparison is fair. The KKT check evaluates the marginal gain per unit cost
$$
\frac{r_t w_i a_i e^{-a_i x_{t,i}}}{c_i}
$$
on all interior pairs and confirms that it equals $\lambda^{*}$ up to floating-point error.
6. Sensitivity surface. The fast solver is called for a $16 \times 16$ grid of budgets and storm scales, 256 full optimizations in total. This is only practical because each solve takes a fraction of a millisecond.
7. Figure. All six panels are drawn into one figure and shown with a single plt.show().
Making It Fast
The generic SLSQP solver treats the problem as a black-box nonlinear program with 144 variables. Each iteration solves a dense quadratic subproblem whose cost grows roughly cubically with the number of variables. It also needs many iterations, and its runtime grows quickly as $T$ or $I$ increases.
The Lagrangian approach exploits the structure of the problem. The coupling between all $T \times I$ variables happens only through one scalar constraint, so the problem collapses to a one-dimensional root-finding problem in $\lambda$. Each bisection step is a single vectorized array operation, and the complexity is $O(T \cdot I \cdot \text{iters})$, which is linear in problem size. This is why the 256-point sensitivity surface costs almost nothing. The same idea scales to thousands of slots and hundreds of instruments.
Execution Results
================================================================== Space Weather Observation Resource Allocation - Results ================================================================== Budget B : 60.00 Cost used (Lagrangian solution) : 60.0000 Optimal multiplier lambda* : 0.266670 Max KKT deviation (interior) : 2.220e-16 ------------------------------------------------------------------ Utility Uniform : 13.5547 Utility Risk-proportional : 19.9097 Utility SLSQP (generic) : 24.9295 Utility Lagrangian bisection : 24.9295 Gain over uniform : 83.92 % Gain over risk-proportional : 25.21 % ------------------------------------------------------------------ Time SLSQP : 0.4205 s Time Lagrangian bisection : 0.003741 s Speed-up : 112.4 x ------------------------------------------------------------------ Resource Total effort Cost share [%] Coronagraph 2.956 14.78 Magnetograph 1.326 5.53 Magnetometers 9.674 24.19 Ionosondes 6.213 12.43 GNSS-TEC 11.204 18.67 L1 Monitor 3.661 24.41 ==================================================================
Interpreting the Results
Console output. The cost used should match the budget $B$ almost exactly, which confirms that the budget constraint is active. The KKT deviation should be on the order of machine precision, showing that all active pairs share the same marginal value per unit cost $\lambda^{*}$. The utility list ranks the four strategies. The Lagrangian solution is the exact optimum of the convex problem, so SLSQP can match it but never beat it beyond numerical tolerance, and both must exceed the two heuristics. The timing lines quantify the benefit of exploiting problem structure. The per-resource table shows where the budget actually goes.
Panel 1 (3D allocation surface). Ridges appear at the two risk pulses, around hours 6 and 15. Along the resource axis the surface is highest for resources with a large ratio $w_i a_i / c_i$, because those cross the activation threshold 
first and are then pushed further by the logarithmic growth of $x^{*}$. Expensive resources appear mainly at the ridge peaks, and cheap resources also cover the quiet hours.
Panel 2 (heatmap). This is the same information from above. Bright vertical bands mark the flare and CME windows. The dark cells show pairs where the threshold is not met, meaning it is optimal to spend nothing there. This sparsity pattern is exactly the water-filling structure of the KKT solution.
Panel 3 (cost composition). The stacked areas show how the spending per hour is split among resources. The total height follows the risk profile, and the colored layers show that additional resources are “switched on” progressively as risk rises, rather than all scaling together.
Panel 4 (strategy comparison). Uniform allocation ignores when observations matter. Risk-proportional allocation ignores which resources are cost-effective and how quickly each saturates. The optimum accounts for both, which is why it clearly outperforms both baselines.
Panel 5 (3D sensitivity surface). Utility increases along both axes, but with different curvature. Along the budget axis the surface flattens, showing diminishing returns: the exponential saturation makes each extra unit of budget less valuable. Along the storm-scale axis, stronger events raise the achievable utility because more risk-weighted information is available to capture. The red marker shows the base case. The surface answers practical questions such as how much additional budget is needed to keep the same utility during a more active period.
Panel 6 (risk versus spending). The yellow curve is the risk profile and the teal bars are the total cost spent per hour. The bars follow the curve but are not proportional to it. Because of the concave utility, the optimal spending rises more gently than the risk itself, which avoids over-investing in a single hour.
Conclusion
The allocation of space weather observation resources can be written as a convex program with exponential-saturation utilities. Its KKT conditions yield a closed-form water-filling solution parametrized by a single Lagrange multiplier, which can be found by bisection in a fraction of a millisecond. Compared with uniform and risk-proportional allocation, the optimal plan concentrates effort where risk is high and where information is cheap, and the sensitivity surface shows how the optimal value responds to budget and storm intensity. The same framework extends naturally to additional instruments, finer time resolution, or richer risk forecasts.
