Ensemble Optimal Control

• Ensemble optimal control problems design a universal input to steer a continuum of systems with parameter variations while minimizing an averaged cost functional.
• An iterative linearization and sweep algorithm, incorporating forward simulations and backward Riccati equations, efficiently addresses the inherent bilinear and infinite-dimensional challenges.
• Robust implementations in quantum control and neuroscience demonstrate the method's ability to synchronize heterogeneous systems with energy-efficient, convergent strategies.

Ensemble optimal control problems concern the design of a universal control law that simultaneously steers a continuum (or a large collection) of parameter-varying dynamical systems—often bilinear, linear, or nonlinear stochastic ODEs—indexed by a compact set Θ of parameter values. Each system receives the same control input, and the performance objective is typically an averaged cost across the ensemble, possibly incorporating state-dependent terminal terms and quadratic penalties on the common open-loop control. These problems are central to robust quantum control, neuroscience, molecular imaging, and any discipline requiring simultaneous actuation of multiple uncertain or structurally heterogeneous systems. The intrinsic challenge is the non-convex, infinite-dimensional, and often bilinear nature of the parameterized system family, which precludes naive decoupling or direct application of classical optimal control approaches.

1. Problem Formulation and Fundamental Structure

An ensemble optimal control problem is typically specified as follows:

• Ensemble member dynamics: Each system’s state $x(t,θ)\in\mathbb{R}^n$ evolves by $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$ where $θ \in Θ$, a compact parameter set. The bilinear term $N_i(θ) x u_i$ introduces nonlinearity and parameter dependency in the control-channel interactions, with all $A,B,N_i,f_0$ continuous in $θ$.
• Control signal: $u(t)\in\mathbb{R}^m$ is broadcast identically to all ensemble members.
• Cost functional:

$$J(u) = \int_{0}^{T} [x(t,θ)^{\top} Q(θ) x(t,θ) + u(t)^{\top} R u(t)]\,dt + x(T,θ)^{\top} G(θ) x(T,θ)$$ with $Q(θ)$ (possibly zero), $G(θ)$ positive semidefinite, $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$0, and typically integrated (weighted) over $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$1.

• Free-endpoint penalty: Rather than imposing hard terminal state constraints $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$2 (which may be infeasible for all $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$3 with a single $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$4), the deviation from a desired profile is penalized in the terminal cost.

In the stochastic extension, each $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$5 may be driven by system-specific independent noise terms (additive Wiener or Poisson processes).

2. Iterative Linearization and the Sweep Algorithm

The bilinearity of the ensemble system renders direct solution intractable. The principal innovation is a Picard–type sweep (fixed-point) algorithm:

1. Freeze cross-terms: Given the iterate $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$6, linearize the bilinear term:

$$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$7

2. Form linear time-varying ensemble:

$$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$8

1. Pontryagin’s Principle and Riccati PDEs:
   • The optimal open-loop $$\frac{\partial}{\partial t}x(t,θ) = A(θ)x(t,θ) + B(θ)u(t) + \sum_{i=1}^m N_i(θ) x(t,θ) u_i(t) + f_0(θ)$$9 is synthesized via stationarity in the Hamiltonian:

   $θ \in Θ$0

• The costate satisfies

  $θ \in Θ$1

  where $θ \in Θ$2 evolve backward via Riccati PDEs (or ODEs if $θ \in Θ$3 discrete): \begin{align*} \frac{\partial K}{\partial t} &= -K A - A^{\top} K + KΛR^{-1}Λ^{\top}K \ \frac{\partial s}{\partial t} &= -[A^{\top} - KΛR^{-1}Λ^{\top}] s - Kf_0 \end{align*} with terminal conditions $θ \in Θ$4, $θ \in Θ$5.

1. Forward–backward integration: Propagate $θ \in Θ$6 forward in time with updated controls, then iterate.

This process is repeated until convergence.

3. Convergence Theory and Optimality Characterization

Contraction and optimality analysis rest on quantifying the mapping from $θ \in Θ$7 to $θ \in Θ$8 in sup–norms over $θ \in Θ$9:

• Contraction condition: If $N_i(θ) x u_i$0 (control penalty) is chosen sufficiently large relative to the magnitude of the bilinear terms and $N_i(θ) x u_i$1 is uniformly bounded,

$N_i(θ) x u_i$2

• Limit behavior: Iteration converges to $N_i(θ) x u_i$3, yielding an optimal control $N_i(θ) x u_i$4.

• Optimality in the sense of HJB: With mild regularity (the value function $N_i(θ) x u_i$5 is $N_i(θ) x u_i$6 and Lipschitz), the solution satisfies the Hamilton–Jacobi–Bellman equation for the original bilinear problem, and $N_i(θ) x u_i$7 is a global optimizer due to strict convexity in $N_i(θ) x u_i$8 (since $N_i(θ) x u_i$9).

4. Stochastic Extensions

If the ensemble is driven by additive noise, each trajectory $A,B,N_i,f_0$0 evolves by

$A,B,N_i,f_0$1

or includes Poisson jumps. The same ensemble iterative method applies, but the control is designed to guide the mean state according to

$A,B,N_i,f_0$2

with a cost penalizing the expected terminal error. Thus, the structure of the control synthesis and the convergence guarantees are preserved for the mean trajectory under independent noise realizations.

5. Computational Implementation and Applications

Efficient practical implementation is achieved by discretizing the parameter domain $A,B,N_i,f_0$3 and employing forward-backward sweeps and numerical Riccati integration. Notable applications:

• Spiking control for an ensemble of integrate-and-fire neurons:

  • Example: $A,B,N_i,f_0$4 with varying $A,B,N_i,f_0$5 or $A,B,N_i,f_0$6, discretized into $A,B,N_i,f_0$7 samples.
  • Converges in 10–20 iterations.
  • Quantitatively achieves robust, energy-efficient, and synchronized firing across neuronal heterogeneity.
• Broadband excitation in quantum spin ensembles:
  • Example: Bloch equations for two-level systems with Larmor frequency dispersion, steering from $A,B,N_i,f_0$8 to $A,B,N_i,f_0$9.
  • Ensemble control pulse designed using $θ$0 iterations, achieving ≥99% fidelity across the frequency band.

Computationally, the dominant costs arise from discretization (integration over $θ$1), the numerical solution of coupled Riccati equations, and propagation of the linearized ensemble dynamics. The approach is robust and efficient for moderate parameter-space dimensionality; scalability remains a topic for further refinement (e.g., sparse grid quadrature or randomized SVD for very high–dimensional parameter spaces).

6. Algorithmic Workflow and Theoretical Implications

The overall algorithm executes the following key steps:

1. Initialization: Discretize $θ$2, initialize $θ$3, $θ$4, $θ$5.
2. Iterative linearization: Freeze the bilinear term at each iterate.
3. Solve time-varying ensemble LQR: Integrate Riccati equations and compute optimal $θ$6.
4. Forward simulation: Update $θ$7.
5. Check convergence: Terminate if the change is below threshold; else, repeat.

This approach leverages the structure of bilinear ensemble control to obtain global, robust, and energy-efficient solutions. The results extend to stochastic systems and are directly validated in real-world applications relevant to quantum control and neuroscience. Theoretical guarantees—contraction with large $θ$8 and satisfaction of the ensemble HJB equation—solidify its position as a foundational method in ensemble control theory.

7. Broader Significance and Extensions

This framework demonstrates the crucial role of ensemble-based optimization for system families with structural (parametric) inhomogeneity, especially where hard endpoint constraints are infeasible. The method provides tractable, convergent algorithms for challenging infinite-dimensional optimal control tasks with strong applicability in experimental domains. Extensions to constraints on the control, non-quadratic costs, and integration with uncertainty quantification (e.g., stochastic parameter distributions) are well-supported by the bilinear iterative paradigm, underpinning a broad class of ensemble control strategies within engineering and physical sciences (Wang et al., 2016).