Sampling-Based Planning and Model Predictive Control › Random Shooting and Cross-Entropy Methods

Random shooting is the most stripped-down planner you can write. You fix a planning horizon $H$, which is the number of future steps you intend to consider. You fix a control distribution, usually uniform over the admissible action set $\mathcal{U} = [-u_{\max}, u_{\max}]^{d_a}$. You draw $N$ independent action sequences,

roll each sequence through the model to obtain a predicted trajectory $\tau^{(i)}$, evaluate a cost $J(\tau^{(i)})$, and return the sequence with the lowest cost. An open-loop controller may execute the whole sequence. When random shooting is used inside MPC, only the first action or a short block is executed before replanning. That is the entire search algorithm.

Notice what is absent from this recipe. There is no gradient, no Hessian, no linearization, no convexity assumption, and no requirement that the model be differentiable at all. The planner treats the dynamics as an opaque function mapping a state and an action to a next state. It queries that function only in the forward direction. In optimization language this is zeroth-order or derivative-free optimization: it uses function values alone to search. This single property is what makes random shooting compatible with trained neural networks, with simulators that contain lookup tables or branch logic, and with physical systems that can only be probed by experiment. It is also why the method scales so well with parallel hardware: every one of the $N$ rollouts is independent, so the whole batch can be evaluated on a GPU at once.

The reason this works at all is that a model rollout may be cheap relative to acting in the real world. The limitation is the curse of dimensionality of action-sequence space: a regular grid with resolution $\epsilon$ along each coordinate contains $(1/\epsilon)^{H d_a}$ length-$H$ sequences for an action dimension $d_a$. A 10-step horizon with two action coordinates and 100 levels per coordinate would contain $100^{20}$ grid points. Random sampling avoids enumerating that grid but does not remove the dimensional dependence. For a regular grid with $N$ points in a $H d_a$-dimensional unit box, the per-coordinate spacing scales as $N^{-1/(H d_a)}$. Random samples do not form a regular grid, so this curve is a resolution heuristic rather than an exact nearest-neighbor law; it still exposes the exponential dependence on sequence-space dimension.

The shape of that heuristic is worth seeing. The cell below evaluates $N^{-1/(H d_a)}$ over five decades of point budget for three sequence-space dimensions and plots the result on log-log axes.

Multiplying the grid-point budget by $10{,}000$ reduces this spacing by a factor of 100 in two dimensions, but only by a factor of about 1.58 in twenty dimensions. The calculation is the argument: holding coordinate resolution fixed requires a budget exponential in dimension.

The practical consequence is that uniform sampling becomes inefficient as horizons and action dimensions grow. Whether it remains adequate depends on task structure, action parameterization, and compute budget.

CEM replaces the fixed proposal with a parametric distribution refitted to the lowest-cost samples in each current batch. In the diagonal-Gaussian form used here, every time-action coordinate has a mean $\mu$ and standard deviation $\sigma$ and is sampled independently before clipping. The refitted mean moves toward the current elites, while the spread may contract or expand coordinate by coordinate. In the version implemented below:

1. Initialize $\mu_0 = \mathbf{0}$ and $\sigma_0 = u_{\max} \mathbf{1}$.
2. For iteration $k = 1, \dots, K_{\text{iter}}$:
   • Draw $N$ samples $a^{(i)} \sim \mathcal{N}(\mu_{k-1}, \operatorname{diag}(\sigma_{k-1}^2))$, clipped to $[-u_{\max}, u_{\max}]$.
   • Evaluate each sample's cost $J_i = J(\text{rollout}(s, a^{(i)}))$.
   • Select the top $N_\text{elite}$ samples by lowest cost.
   • Set $\mu_k$ and $\sigma_k$ to the mean and standard deviation of those elites.
3. Return $\mu_{K_{\text{iter}}}$ as the plan.

The reason this works is the general elite-set idea of evolutionary methods. The elite samples identify a region with low observed cost, and refitting narrows the proposal distribution toward that region. Because the standard deviation usually shrinks, subsequent samples concentrate more finely on promising actions. The name "cross-entropy" comes from choosing a sampling distribution that is close in Kullback-Leibler divergence to an ideal rare-event or elite distribution, as developed in de Boer et al.'s tutorial. In the finite-sample Gaussian version used here, maximum-likelihood fitting to elites produces their sample mean and standard deviation. This is iterative distribution fitting; it should not be treated as literal projected gradient ascent on the original objective.

This iterative refitting distinguishes CEM from one larger random batch. Random shooting spends its budget on independent samples from a fixed proposal; CEM uses earlier evaluations to change later proposals. This can help when a representable low-cost region occupies little of the original proposal. It can hurt when elites span unrelated modes, because a single fitted Gaussian may average them into a poor region. A smaller elite fraction, explicit smoothing, covariance control, or a mixture proposal can change that behavior, but each must be tuned empirically.

Let us set up the running example. Our system is a 2D double integrator, a simplified model for a point mass commanded by acceleration in two Cartesian axes. The state is $s = (x, y, v_x, v_y)$ and the action $a = (a_x, a_y)$ is an acceleration command, bounded by $|a_i| \le u_{\max}$. The discrete-time dynamics are

and identically in $y$. The cost encodes our intent. To reach a goal $g$ we penalize squared distance from $g$ over the trajectory plus a small penalty on control effort:

A few remarks on why this example is worth the space. Position-only feedback does not directly damp the double integrator's velocity, so reaching a target smoothly requires accelerating and then decelerating. The displayed update is the exact one-step solution of continuous constant-acceleration dynamics when the action is held fixed over each interval. Because the nominal dynamics are available in closed form, we can use them as the planner's model and separately perturb the simulated plant. That isolates the effect of replanning under a controlled model-plant mismatch. The effort term with weight $\alpha$ discourages unnecessarily large or oscillatory controls; it does more than merely break exact ties.

Let us implement this system and the two planners.

Two small comments on the helpers. First, we batch everything. A planner that evaluates thousands of candidate sequences per control tick benefits strongly from vectorized rollout, so we operate on arrays of shape (B, T, 2). A Python loop over candidate plans would add interpreter overhead and would make the later timing comparison measure implementation style as much as planner structure. Second, the cost is defined on the equally spaced trajectory returned by batched_rollout. In a learned system, this function could call a neural transition model instead of the analytic equations. A replacement model can use the same planner interface when it matches the expected batch shapes, state semantics, time step, and output convention.

Now the two planners. Random shooting is a one-liner wrapped in a function.

Read the function from top to bottom and you see the whole algorithm in four lines: sample, roll out, evaluate, return the argmin. The returned cost array supports diagnostics, but its shape needs context. A flat high-cost batch suggests the proposal found no useful distinction; a flat low-cost batch may mean that many plans work. A single low outlier can be a better plan or a noisy or model-exploiting anomaly, so it should be checked rather than trusted automatically.

The cross-entropy method shares the same evaluation loop but keeps the sample distribution as explicit state. The keyword mu_init lets a caller provide a warm-start mean; we will use it later for MPC.

There are three implementation details worth dwelling on because they are easy to get wrong. The first is the cost_fn indirection. Passing the cost as a callable, rather than hard-coding the goal-distance cost, is what lets us later bolt on obstacle penalties and other constraint terms without touching the planner. The planner's job is to minimize whatever scalar the cost function returns; it has no opinion about what that scalar means.

The second is the np.clip inside the sample loop. Without it, the elite set could contain inadmissible actions. Clipping keeps evaluated controls inside the box but creates point mass at its boundaries: the elite moments are fitted to clipped samples while the next batch comes from an unclipped normal. When substantial proposal mass lies outside the action box, a truncated distribution or bounded reparameterization avoids that mismatch.

The third is the tiny constant 1e-4 added to the standard deviation. It prevents the distribution from collapsing to a point, which would make every subsequent sample identical and freeze the search, and it also prevents division-by-zero type pathologies in any downstream code that uses sigma. The value is a floor on exploration; if the planner needs to keep searching aggressively, it should be larger, and if the planner should commit to the plan it has found, it should be smaller.

The concentration that the mean and standard deviation updates produce is worth seeing directly. The cell below runs a short CEM with a fixed seed and records the first-step action samples at the first and fourth iterations, together with the elite mean fitted at each. Ties inside the elite set are broken by lowest cost, then lowest sample index.

Let us now compare the two planners on a single starting state. The start is the origin, the goal is $(3, 2)$, and the horizon is 25 steps, which at $\Delta t = 0.1$ s is 2.5 seconds of simulated time.

In this fixed-seed run, CEM achieves a lower objective value than random shooting. The accounting matters: random shooting evaluates 500 sequences, while CEM evaluates $15 \times 800 = 12{,}000$, so this is a mechanism demonstration rather than an efficiency benchmark. A fair efficiency comparison would equalize wall-clock time or model evaluations and repeat across seeds. Here the paths show the plans selected under the stated, unequal budgets; the objective values establish only which plan scores lower under this chapter's cost function.

The figure records a cost advantage for this run, and the right panel shows the lowest sampled cost in each CEM iteration. A steep early drop followed by a plateau is consistent with the proposal locating and refining one low-cost region. A curve still falling at the budget boundary suggests that extra iterations are worth testing. A flat, high curve is less diagnostic: poor initialization, inadequate proposal variance, too few samples, an unsuitable elite fraction, or a flat objective can all produce it.

The elite fraction is one knob behind that diagnostic. The sweep below holds the start state, horizon, sample budget, iteration count, and pseudorandom seed fixed, changing the number of elites. Within this deterministic implementation, the configured independent variable is the fraction; the result is still one optimizer seed on one task. Ties are ordered by cost and then sample index.

Two practical notes before we move on. First, the split between samples and CEM iterations matters, not just their product. Iterations let the proposal migrate, whereas samples improve each empirical elite estimate; too few samples can make that estimate unstable. Which split is better therefore depends on dimension, elite fraction, objective, and overhead. Second, in a deterministic open-loop setting, both planners produce a plan calibrated to one starting state. That plan is not a feedback controller: it does not react to disturbances, model error, or unexpected observations. To close the loop, we need MPC.