Paper deep dive
Safety-Critical Contextual Control via Online Riemannian Optimization with World Models
Tongxin Li
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 4/26/2026, 11:42:39 PM
Summary
The paper introduces a framework for safety-critical contextual control using black-box world models. It addresses the challenge where a Planner must optimize a task objective without explicit knowledge of the system dynamics, relying instead on feasibility samples from a Simulator. The core contribution is a sample-based Penalized Predictive Control (PPC) framework that uses online Riemannian optimization. The Simulator compresses the time-varying feasibility manifold into a score-based density signal, which the Planner uses for gradient descent. A key theoretical result is the introduction of 'barrier curvature' (the minimum eigenvalue of the Fisher information), which replaces the traditional Lipschitz constant in safety bounds. The framework demonstrates that using a context signal to estimate a conditional density significantly improves safety and performance compared to marginal density models, especially during environmental shifts.
Entities (7)
Relation Signals (5)
Context Signal → conditions → Feasibility Density
confidence 100% · enabling conditional density estimation p̂(u|ξt)
Barrier Curvature → governs → Convergence Rate
confidence 100% · The barrier curvature κ(ξt)... governs both convergence rate and safety margin
Planner → optimizesobjectiveusing → Penalized Predictive Control (PPC)
confidence 100% · The Planner optimizes an objective over a feasibility set... using the score of the learned density
Simulator → providessamplesto → Planner
confidence 100% · the Simulator produces feasibility samples... which the Planner uses for gradient-based optimization
Penalized Predictive Control (PPC) → uses → Riemannian Optimization
confidence 100% · a sample-based Penalized Predictive Control (PPC) framework grounded in online Riemannian optimization
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Modern world models are becoming too complex to admit explicit dynamical descriptions. We study safety-critical contextual control, where a Planner must optimize a task objective using only feasibility samples from a black-box Simulator, conditioned on a context signal $\xi_t$. We develop a sample-based Penalized Predictive Control (PPC) framework grounded in online Riemannian optimization, in which the Simulator compresses the feasibility manifold into a score-based density $\hat{p}(u \mid \xi_t)$ that endows the action space with a Riemannian geometry guiding the Planner's gradient descent. The barrier curvature $\kappa(\xi_t)$, the minimum curvature of the conditional log-density $-\ln\hat{p}(\cdot\mid\xi_t)$, governs both convergence rate and safety margin, replacing the Lipschitz constant of the unknown dynamics. Our main result is a contextual safety bound showing that the distance from the true feasibility manifold is controlled by the score estimation error and a ratio that depends on $\kappa(\xi_t)$, both of which improve with richer context. Simulations on a dynamic navigation task confirm that contextual PPC substantially outperforms marginal and frozen density models, with the advantage growing after environment shifts.
Tags
Links
- Source: https://arxiv.org/abs/2604.19639v1
- Canonical: https://arxiv.org/abs/2604.19639v1
Trouble viewing inline? Open PDF directly →
Full Text
131,756 characters extracted from source content.
Expand or collapse full text
Safety-Critical Contextual Control via Online Riemannian Optimization with World Models Tongxin Li Abstract Modern world models are becoming too complex to admit explicit dynamical descriptions. We study safety-critical contextual control, where a Planner must optimize a task objective using only feasibility samples from a black-box Simulator, conditioned on a context signal ξt _t. We develop a sample-based Penalized Predictive Control (PPC) framework grounded in online Riemannian optimization, in which the Simulator compresses the feasibility manifold into a score-based density p^(u∣ξt) p(u _t) that endows the action space with a Riemannian geometry guiding the Planner’s gradient descent. The barrier curvature κ(ξt)κ( _t), the minimum curvature of the conditional log-density −lnp^(⋅∣ξt)- p(· _t), governs both convergence rate and safety margin, replacing the Lipschitz constant of the unknown dynamics. Our main result (Theorem 4) is a contextual safety bound showing that the distance from the true feasibility manifold is controlled by the score estimation error and a ratio that depends on κ(ξt)κ( _t), both of which improve with richer context. Simulations on a dynamic navigation task confirm that contextual PPC substantially outperforms marginal and frozen density models, with the advantage growing after environment shifts. I Introduction The increasing complexity of modern world models is fundamentally changing the interface between planning and control. Neural physics simulators [32], video-generation models [14], and foundation-model-based world models [17, 9] can predict the consequences of actions with remarkable fidelity, yet their internal representations (multi-modal embeddings, high-dimensional latent states, or raw image sequences) are far too complex to admit explicit and known dynamical descriptions of the form xt+1=ft(xt,ut)x_t+1=f_t(x_t,u_t). These black-box world models offer a compelling advantage over classical model-based control because they capture richer dynamics than any tractable closed-form approximation, generalize across environments without manual re-derivation, and directly leverage the rapid advances in foundation models and learned simulators. However, a controller that must guarantee safety while interacting with such a black-box world model faces a fundamental challenge. It cannot write down the dynamics, cannot form the constraint Jacobians that classical Model Predictive Control (MPC) requires, and cannot construct the explicit barrier functions used in safety-critical control (Figure˜1). Environmentxt+1=ft(xt,ut)x_t+1=f_t(x_t,u_t)ExplicitModel ftf_tJacobians∇uft _uf_tControllerminc(u) \,c(u)CBF-QPFilterxtx_t∇h∇ hutu_t(a)EnvironmentBlack-BoxSimulatorKDEp^t,s^t p_t,\, s_tPlannerminℱ(u) \,F(u)SafetyFilterxtx_tsampless^t,βt s_t,\, _tuppcu_ppcutu_tcontext ξt _t(b) Figure 1: Two closed-loop control paradigms. (a) Classical model-based control requires an explicit dynamics model ftf_t to form constraint Jacobians and a CBF-QP safety filter that projects the controller’s action onto the safe set. (b) The contextual control framework replaces the explicit model with a black-box Simulator that produces feasibility samples; a KDE compresses them into a score signal s^t s_t, which the Planner uses for gradient-based optimization. A safety filter (score-ascent retraction) projects the candidate uppcu_ppc onto ℳ^tα M_t^α if needed. The Planner also receives a context signal ξt _t from the environment (dashed arrow), enabling conditional density estimation p^(u∣ξt) p(u _t). See Figure˜3 for the detailed per-step information flow. We address this challenge by developing a framework for safe control through black-box world-model samples. The controller optimizes an objective over a feasibility set whose geometry is accessible only through feasibility samples drawn from the world model, rather than through an explicit dynamics model. The world model encapsulates all complexity of the physical environment, and the controller’s access to safety-relevant information is mediated entirely by a learned density signal that compresses the simulator’s knowledge of the feasibility manifold. This architecture separates what is feasible (determined by the black-box world model) from what is desirable (the controller’s objective), and the score function of the learned density provides a differentiable interface between the two. This separation is not merely an engineering convenience; it is an inevitability as world models become richer and more expressive. The architecture arises naturally in robotics with visual world models [17, 16], where a planner cannot extract closed-form obstacle constraints from a neural simulator’s latent state. It appears equally in multi-agent energy dispatch [23, 24], where an aggregator coordinates distributed resources whose local feasibility sets are private and time-varying, so the aggregator must plan using only sampled feasibility signals rather than explicit local models, and in autonomous driving [29], where a predictive model simulates traffic scenarios but its internal state is not available to the planner in closed form. In each scenario, the feasibility manifold ℳtM_t is a time-varying, non-convex subset of the decision space whose geometry is determined by the world model’s internal state. The controller must optimize over ℳtM_t without an explicit dynamics model, learning the manifold structure from online samples while guaranteeing safety at every step. This paper addresses the following question. Can a controller guarantee safety when its only access to the constraint geometry is through black-box world-model samples? We answer affirmatively. Our approach treats the world model as a Simulator that compresses the geometry of the time-varying feasibility manifold into a score-based density signal, which a Planner then uses for safe gradient-based optimization. When the Planner’s observation is a context ξt _t, a fixed-dimensional embedding of rich sensory data rather than the Simulator’s full state, we estimate a conditional feasibility density p^(u∣ξt) p(u _t) so that the learned manifold tracks changes in the environment that are reflected in ξt _t. This perspective builds on classical connections between control and probabilistic inference [35, 20] and extends recent score-based manifold optimization [21] to the online, safety-critical setting. I-1 Model-Free and Simulator-Based Control Classical MPC [15] and its robust variants [5, 25] require an explicit dynamics model to construct constraints and compute Jacobians. Data-driven MPC [6, 7, 30] learns surrogate models but still assumes the learned model is amenable to explicit constraint formulation. Recent work on predictive safety filters [38] and robust adaptive NMPC [10] further improves constraint handling, yet remains tied to an explicit or partially identified plant model. The Cross-Entropy Method (CEM) [31] is a popular sampling-based alternative that draws candidate actions from a parametric proposal, evaluates them against a model, and refits the proposal to an elite subset. CEM-MPC variants [13, 27] combine CEM with learned dynamics models, but provide no safety certificates and scale poorly with action dimension. The control-as-inference viewpoint [35, 20] and rate-distortion theory [34] provide a principled way to trade off performance and information cost, but existing instantiations assume known dynamics. Our framework makes no such assumption. The world model is accessed purely through feasibility samples, and the information-theoretic penalty emerges as the Lagrange multiplier of the sampling capacity, providing a mechanism for adaptive trust calibration unavailable in model-based or sampling-based approaches. I-2 Riemannian Optimization and Density Models Score matching [19, 37] estimates the gradient of the log-density without computing the normalizing constant. Classical Riemannian optimization [8, 1] requires explicit manifold descriptions, but recent work [21] shows that score functions of smoothed densities can recover the manifold projection and tangent-space operations needed for optimization in a static, offline regime. Normalizing Flows [28, 22] have been applied to trajectory forecasting [29] and policy parameterization [11], while the information-aggregation framework of [23, 24] encodes distributed constraint geometry into a feasibility density. None of these works addresses online score estimation with provable safety guarantees. Our framework fills this gap by extending density-based manifold operations to the online setting where the manifold is unknown and time-varying, and provides the first safety certificates for score-based feasibility learned from black-box samples. I-3 Data-Driven Safety Certificates Control Barrier Functions (CBFs) [33, 2, 3] parameterize safe sets as superlevel sets x:h(x)≥0\x:h(x)≥ 0\ of learned functions. Recent extensions address parametric uncertainty [26] and enable purely data-driven CBF synthesis from input-output measurements [18, 4]. These methods suffer from structural bias: maintaining tractability requires restricting h to specific function classes (polynomials, GPs [12]), which under-approximate complex, non-convex geometries. GP-based approaches additionally scale poorly with state dimension due to cubic kernel-matrix inversion. Our score-based approach can represent arbitrary manifold topologies without structural assumptions on h, scales favorably in the control dimension, and refines its safety certificate online as samples accumulate. Our contributions are three-fold. First, we formalize contextual control with black-box world models as a new problem class (Section˜I) in which the controller has no access to the dynamics ftf_t and receives only feasibility samples from a black-box Simulator, yet must certify safety at every step. To our knowledge, this is the first framework that unifies black-box simulator access, online density estimation, and provable safety guarantees under a single formulation. Second, we develop a sample-based instantiation of Penalized Predictive Control (PPC) [23] for the contextual black-box setting. We derive a rate-distortion formulation (Proposition 1) that yields the penalty weight β as the Lagrange multiplier of the sampling capacity, together with an online score-estimation pipeline that replaces the oracle density with a learned estimate. The resulting controller uses the score of the learned density as its sole interface to the constraint geometry, replacing the explicit Jacobians and barrier functions of classical MPC. A central novelty is that the barrier curvature κ, the minimum eigenvalue of the Fisher information ℐp^I_ p, plays the role that the dynamics Lipschitz constant plays in standard safety-critical control. Because κ is a property of the learned density rather than the unknown dynamics, it can be estimated from data and used to set β adaptively, which is impossible in classical barrier-function methods. Third, we establish a sequence of convergence and safety guarantees, each building on the preceding one. Theorem 1 shows that PPC iterates contract geometrically at a rate set by the barrier curvature κ, and Proposition 2 identifies the critical stiffness above which the equilibrium is guaranteed to lie inside the learned manifold. When the manifold drifts, Theorem 2 bounds the dynamic regret in terms of the score variation TV_T, revealing a dual role of κ, namely it speeds contraction and simultaneously dampens sensitivity to manifold drift. Theorem 3 then bounds the score estimation error at rate O(N−2/(m+4))O(N^-2/(m+4)), which feeds into Theorem 4 to yield a safety bound whose residual Gc/(βκ)G_c/(βκ) depends on the learned density’s curvature rather than a Lipschitz constant of the unknown dynamics. Finally, Proposition 3 derives, via a variational-inference argument, a mixture-Hessian identity showing that marginalizing out ξt _t subtracts a posterior-score covariance from the learned curvature, and Theorem 5 uses this to give a formal safety-guarantee gap of at least Gcσt2/(βtκ(ξt)κmarg)G_c\, _t^2/( _t\,κ( _t)\, _marg) between the Theorem 4 residuals of the context-aware and context-blind Planners whenever the current context is identifiable from the action, with σt2 _t^2 the minimum eigenvalue of the posterior conditional-score covariance. To our knowledge, this is the first formal characterization of when and by how much context improves safety in a black-box control setting. I Problem Formulation We consider an architecture composed of a Planner responsible for optimizing a global objective and a black-box Simulator (world model) that encapsulates the physical environment. The system evolves according to time-varying nonlinear dynamics: xt+1=ft(xt,ut),t=0,…,T,x_t+1=f_t(x_t,u_t), t=0,…,T, (1) where xt∈⊆ℝnx_t ^n is the system state and ut∈⊆ℝmu_t ^m is the control input. The dynamics ftf_t are embedded inside the Simulator and are not available to the Planner in closed form. The time-dependence of ftf_t captures exogenous factors (moving obstacles, changing network topology, or shifting workload patterns) that the Simulator can observe and simulate but that remain hidden from the Planner. Notation. Throughout, ∥⋅∥\|·\| denotes the Euclidean norm and A⪯BA B the Löwner order on symmetric matrices. We write dist(u,S)=infv∈S‖u−v‖dist(u,S)= _v∈ S\|u-v\| for point-to-set distance and dHd_H for Hausdorff distance. The following symbols recur across sections and are collected here for reference: T is the planning horizon; K is the number of inner gradient steps per PPC update; NtN_t is the cumulative number of feasibility samples available at time t; p^t p_t is the learned feasibility density and s^t=∇ulnp^t s_t= _u p_t its score function; ℳ^tα=u:p^t(u∣ξt)≥α M_t^α=\u: p_t(u _t)≥α\ is the learned α-superlevel set at time t; ℐp^(u)=−∇2lnp^(u)I_ p(u)=-∇^2 p(u) is the local Fisher information of the learned density, with eigenvalue bounds κ (lower, the barrier curvature) and Λ (upper) formalized in Definition 3 and Assumption 3. I-A The Unknown Feasibility Manifold Let safe⊂X_safe be a closed set of safe states. We define the feasibility manifold ℳt(xt)M_t(x_t) as the set of all inputs that maintain safety: ℳt(xt):=u∈∣ft(xt,u)∈safe.M_t(x_t):=\u f_t(x_t,u) _safe\. (2) Under smoothness of ftf_t, the set ℳtM_t constitutes a time-varying, potentially non-convex submanifold embedded in the control space U, whose topology is governed by the nonlinear dynamics inside the Simulator. The manifold ℳtM_t is unknown to the Planner, time-varying, and accessible only through feasibility samples from the Simulator. Assumption 1 (Compactness and Regularity). safeX_safe is compact and ftf_t is C2C^2-smooth for all t, so that ℳt(xt)M_t(x_t) is a compact embedded submanifold of U. Compactness of safeX_safe is standard in safety-critical control [3]. The C2C^2-smoothness of ftf_t ensures that the feasibility manifold ℳtM_t inherits a well-defined boundary, which is necessary for the score function to be well-defined in the interior. I-B The Simulator–Planner Interface The Simulator and Planner interact through a sample-based interface (Figure˜2). We formalize the black-box world model and the information flow between the two components. Definition 1 (Black-Box World Model). A black-box world model is a feasibility oracle t:×→0,1W_t:X×U→\0,1\ defined by t(x,u)=[ft(x,u)∈safe]W_t(x,u)= 1[f_t(x,u) _safe]. The dynamics ftf_t are embedded inside the oracle and are not accessible in closed form. Given the current state xtx_t and a proposal distribution q(u)q(u), the oracle returns i.i.d. feasibility samples111In practice, these samples can be obtained by any method that produces i.i.d. draws from the feasibility-conditioned distribution, e.g., rejection sampling from q (used in our experiments; see Section V-C), MCMC, or direct forward simulation of the world model. uii=1N\u_i\_i=1^N drawn from the conditional q(u∣t(xt,u)=1)q(u _t(x_t,u)=1). The Simulator has access to the full state xtx_t, the oracle tW_t, and the feasibility manifold ℳtM_t. It produces feasibility samples and transmits a learned density estimate p^t p_t (together with its score s^t=∇ulnp^t s_t= _u p_t) to the Planner. The Planner optimizes a stage cost c(u)c(u) but has no access to ftf_t, xtx_t, or the raw state; it receives only the score signal s^t s_t and the cost function c. Formally, the information sets are ℐtSim _t^Sim =xt,t(⋅),ℳt, =\x_t,\,W_t(·),\,M_t\, (3) ℐtPl _t^Pl =s^t(⋅),c(⋅). =\ s_t(·),\,c(·)\. (4) This information asymmetry (ft∉ℐtPlf_t _t^Pl) means the Planner cannot form constraint Jacobians or build explicit barrier functions; it must rely entirely on the Simulator’s score signal. I-C Score-Based Feasibility Signal Rather than transmitting the explicit manifold geometry to the Planner, the Simulator compresses ℳtM_t into a differentiable density whose score function encodes the manifold geometry. The ideal feasibility distribution and its score are from the existing information-aggregation literature [23, 24]. Definition 2 (Feasibility Distribution). The oracle feasibility distribution pt∗(⋅∣xt)p^*_t(· x_t) is a C2C^2 density supported on the closure of the single-step feasibility manifold ℳt(xt)M_t(x_t), with pt∗>0p^*_t>0 on intℳtint\,M_t, vanishing continuously on ∂ℳt _t, and having a nondegenerate gradient on ∂ℳt _t.222A canonical construction is a Gaussian-smoothed uniform density truncated to ℳtM_t, i.e., pt∗∝((0,σ2I)∗[⋅∈ℳt])⋅[⋅∈ℳt]/Zp^*_t (N(0,σ^2I)* 1[\,·\, _t])· 1[\,·\, _t]/Z, which recovers the manifold projection and tangent-space operations needed for Riemannian optimization [21, 19]. The smoothing scale σ plays the role of the KDE bandwidth h in Section IV-A. The associated oracle score is st∗(u):=∇ulnpt∗(u∣xt)s^*_t(u):= _u p^*_t(u x_t). We note that the maximum-entropy feasibility (MEF) density in [23, 24] assigns higher probability to actions that maximize future flexibility by weighting with the volume of reachable feasible sets over a prediction horizon. Our single-step density is the special case in which no look-ahead is used. Extending the theory to multi-step MEF densities is an interesting direction but beyond the scope of this work. The oracle score st∗s^*_t encodes the geometry of ℳtM_t in a form the Planner can use for gradient-based optimization, without ever accessing ftf_t. The critical challenge is that this score must be learned online from streaming feasibility samples, rather than from a pretrained model. Assumption 2 (Bounded Information Loss). The Simulator employs a learning algorithm (e.g., KDE or score matching) to produce an estimate p^t p_t of the feasibility density. The KL divergence between the oracle and the learned estimate is finite such that DKL(pt∗∥p^t)<∞.D_KL(p^*_t\| p_t)<∞. This assumption is satisfied by any consistent estimator [36]. We instantiate the estimator as a Gaussian KDE in Section IV-A, with convergence rates established in Section IV-D. SimulatorEnvironmentxt+1=ft(xt,ut)x_t+1=f_t(x_t,u_t)World Modelt(xt,u)W_t(x_t,u)KDE /Score ModelPlanneruppc≈argminℱu_ppc≈ *arg\,minFxtx_tuiu_is^t s_tuppcu_ppc Figure 2: The Simulator–Planner architecture. The Simulator contains the environment dynamics and the feasibility oracle tW_t; it draws feasibility samples ui\u_i\ and fits a density model whose score s^t s_t is transmitted to the Planner. The Planner minimizes the penalized objective c(u)−βlnp^t(u)c(u)-β p_t(u) via gradient descent on s^t s_t, without access to ftf_t. The Planner’s goal is to minimize a stage cost c(ut)c(u_t) subject to the unknown manifold constraint ut∈ℳtu_t _t, using the online score estimate s^t=∇ulnp^t s_t= _u p_t transmitted by the Simulator. The Planner does not observe the full state xtx_t but instead receives a context signal ξt∈ℝd _t ^d, a fixed-dimensional summary of the environment. For example, in the robot navigation task of Section V, xtx_t includes all obstacle positions and velocities, while ξt _t is a low-dimensional embedding of a top-down image produced by the Planner’s sensor. The feasibility manifold ℳtM_t depends on xtx_t through the obstacle geometry, but ξt _t correlates with xtx_t and allows the Planner to condition its density estimate as p^(u∣ξt) p(u _t), so that the learned manifold adapts when the environment changes. Without ξt _t, the Planner must rely on the marginal density p^(u) p(u), which averages over all states and cannot track manifold shifts. State-dependent objectives can be absorbed by conditioning the cost on ξt _t (e.g., c(u;ξt)=‖u−ugoal(ξt)‖2c(u; _t)=\|u-u_goal( _t)\|^2 for some context-dependent target ugoalu_goal); we suppress the ξt _t-dependence of c for notational clarity. We formalize the fundamental problem as follows. Problem P (Contextual Control): Given a black-box world model tW_t (Definition 1), a stage cost c:→ℝc:U , and a context signal ξt∈ℝd _t ^d, find a sequence of policies πtt=1T\ _t\_t=1^T that solves minπt∑t=1Tut∼πt[c(ut)]s.t.ut∈ℳt,∀t, _\ _t\\;\; _t=1^TE_u_t _t [c(u_t) ] .t. u_t _t,\;\;∀\,t, (5) using only feasibility samples from tW_t and the context ξt _t, without access to the dynamics ftf_t. Because the Planner observes ξt _t rather than xtx_t, the oracle density from Definition 2 induces two related distributions over actions. The conditional oracle density p∗(u∣ξt)≔xt∣ξt[pt∗(u∣xt)]p^*(u _t) _x_t _t[p^*_t(u x_t)] averages over the states consistent with the current context, while the marginal oracle density p∗(u)≔∫p∗(u∣ξ)p(ξ)ξp^*(u) p^*(u ξ)\,p(ξ)\,dξ averages over all contexts under the context-generating distribution p(ξ)p(ξ). Proposition 3 and Theorem 5 quantify the benefit of working with the former rather than the latter, via the mixture-Hessian identity (48) and the resulting safety-guarantee gap (52). Problem P cannot be solved directly because the Planner cannot evaluate the hard constraint u∈ℳtu _t, and the system dynamics couple the time steps through the state. Following [23, 24], we decouple the problem by optimizing each πt _t independently, treating the current manifold ℳtM_t as given. At each step t, the Simulator compresses the constraint into the learned density p^t(⋅∣ξt) p_t(· _t), and safety is enforced by requiring the policy to remain close to this density. This yields the per-step relaxation: Problem RDP_RD (Per-Step Rate-Distortion Relaxation): minπtu∼πt[c(u)]s.t.DKL(πt∥p^t(⋅∣ξt))≤Ic. _ _t\;\;E_u _t [c(u) ] .t. D_KL( _t\| p_t(· _t))≤ I_c. (6) Solving RDP_RD at each t and summing the costs recovers the structure of P. When p^t→pt∗ p_t→ p^*_t and Ic→0I_c→ 0, the KL constraint forces πt _t to concentrate on ℳtM_t, so per-step feasibility is recovered. The cumulative cost suboptimality of this online policy relative to the offline optimum of P is bounded by Theorem 2. For notational brevity, we often write p^t(u) p_t(u) for p^t(u∣ξt) p_t(u _t) when the context is clear; the full conditional notation is used whenever the distinction matters. I Control Synthesis via Free Energy Minimization We now solve Problem RDP_RD. Since RDP_RD is strictly convex in π, the KKT conditions yield the penalty weight β as the Lagrange multiplier of the capacity constraint IcI_c (not an arbitrary tuning parameter). The critical stiffness condition β>βcurv∗β>β^*_curv derived in Section˜IV ensures that the Planner adheres to the Simulator’s learned manifold geometry with strictness proportional to the barrier curvature κ. Proposition 1 (Thermodynamic Reformulation). The optimal policy πt∗(u)π^*_t(u) solving RDP_RD is the unique minimizer of the Free Energy functional: (π)=u∼π[c(u)]+βtDKL(π∥p^t). (π)=E_u π[c(u)]+ _tD_KL(π\| p_t). (7) The solution is given by the Gibbs-Boltzmann distribution: πt∗(u)=1Ztp^t(u)exp(−c(u)βt), π^*_t(u)= 1Z_t p_t(u) (- c(u) _t ), where βt>0 _t>0 is the Lagrange multiplier associated with the capacity constraint IcI_c, and Zt=∫p^t(v)exp(−c(v)βt)vZ_t= _U p_t(v) \! (- c(v) _t )\,dv is the normalizing constant (partition function). Proposition˜1 shows that the optimal policy balances economic performance against manifold adherence through a Gibbs-Boltzmann distribution. Here (π)G(π) is a functional over policies; the PPC objective in the next subsection is the pointwise free-energy function ℱ(u)=c(u)−βtlnp^t(u)F(u)=c(u)- _t p_t(u) obtained by taking the MAP of πt∗π^*_t. The stiffness βt _t controls how strongly the Planner trusts the Simulator’s feasibility signal. Proof. We form the Lagrangian functional ℒ(π,β,λ)L(π,β,λ): ℒ= = ∫π(u)c(u)u+β(∫π(u)lnπ(u)p^(u)du−Ic) π(u)c(u)du+β ( π(u) π(u) p(u)du-I_c ) +λ(∫π(u)u−1). +λ ( π(u)du-1 ). (8) Taking the Fréchet derivative with respect to π(u)π(u) and setting to zero, δℒδπ=c(u)+β(lnπ(u)−lnp^(u)+1)+λ=0. δπ=c(u)+β( π(u)- p(u)+1)+λ=0. Rearranging the terms we obtain lnπ(u)=lnp^(u)−c(u)/β−(λ/β+1) π(u)= p(u)-c(u)/β-(λ/β+1). Exponentiating yields π∗(u)∝p^(u)e−c(u)/βπ^*(u) p(u)e^-c(u)/β, matching the optimal controllers in Linearly Solvable MDPs [35] and path integral control [20]. ∎ I-A Deterministic Control Law The Planner implements a single action utu_t by taking the maximum a posteriori (MAP) estimate of πt∗π^*_t, yielding the Penalized Predictive Control (PPC) objective [23]. While PPC was originally derived for settings where the oracle density pt∗p^*_t is available in closed form, our contribution is to instantiate it in the black-box regime where pt∗p^*_t is replaced by the online estimate p^t p_t learned from simulator samples, and to provide convergence and safety guarantees for this sample-based variant. Corollary 1 (The PPC Objective). The MAP estimate utMAP=argmaxuπt∗(u)u_t^MAP= _uπ^*_t(u) is equivalent to: utMAP=argminu∈(c(u)−βtlnp^t(u)).u_t^MAP= _u (c(u)- _t p_t(u) ). (9) Proof. By definition, utMAP u_t^MAP =argmaxu1Ztp^t(u)e−c(u)/βt = _u\, 1Z_t\, p_t(u)\,e^-c(u)/ _t =argminu(c(u)/βt−lnp^t(u)). = _u\, (c(u)/ _t- p_t(u) ). Multiplying by βt _t yields the objective. Note that utMAPu_t^MAP is a conceptual unconstrained minimizer over U; we reserve the notation ut∗u_t^* for the restricted minimizer over ℳ^tα M_t^α introduced in Lemma 1(c), which is the object that PPC gradient descent from a warm start in ℳ^tα M_t^α converges to (Theorem 1). Proposition 2 further shows that ut∗u_t^* is itself an unconstrained critical point of ℱF (i.e., ∇ℱ(ut∗)=0 (u_t^*)=0) whenever βt>βcurv∗ _t>β^*_curv, although ℱF may admit additional critical points outside ℳ^tα M_t^α. ∎ The gradient of the PPC objective is ∇u(c(u)−βtlnp^t(u))=∇uc(u)−βts^t(u), _u (c(u)- _t p_t(u) )= _uc(u)- _t s_t(u), (10) where s^t(u)=∇ulnp^t(u) s_t(u)= _u p_t(u) is the Simulator’s learned score. Each gradient step thus combines the cost gradient ∇c∇ c (what is desirable) with the score s^t s_t (what is feasible), and the stiffness βt _t mediates between minimizing c and adhering to ℳtM_t. The Planner never accesses ftf_t; all constraint information enters through s^t s_t. The remainder of the paper answers Problem P by analyzing the per-step free energy ℱt(u)=c(u)−βtlnp^t(u)F_t(u)=c(u)- _t p_t(u). When the time index is clear from context, we write ℱ(u)F(u), p^(u) p(u), and β for brevity. IV Riemannian Score Descent In this section we fix a single time step t and context ξt _t, analyzing the PPC free energy ℱ(u)=c(u)−βlnp^(u∣ξt)F(u)=c(u)-β p(u _t) with p^≡p^t(⋅∣ξt) p≡ p_t(· _t) and β≡βtβ≡ _t held constant. Similarly, we simplify the learned α-superlevel set ℳ^tα=u:p^t(u∣ξt)≥α M_t^α=\u: p_t(u _t)≥α\ as ℳ^α M^α. The subscript k indexes the inner gradient iterations within this step; the extension to time-varying densities is in Section IV-C. IV-A Density Estimator We instantiate the Simulator’s estimator p^t p_t (Assumption 2) as a Gaussian kernel density estimator (KDE) [36], which admits a closed-form score and supports the convergence analysis of Section IV-D. At time t, the Simulator has collected NtN_t cumulative feasibility samples U1,…,UNt⊂ℝm\U_1,…,U_N_t\ ^m, and estimates the marginal feasibility density by p^N(u)=1N∑i=1NKh(u−Ui),Kh(z)=e−‖z‖2/(2h2)(2πh2)m/2, p_N(u)= 1N _i=1^NK_h(u-U_i), K_h(z)= e^-\|z\|^2/(2h^2)(2π h^2)^m/2, (11) with bandwidth h=Θ(N−1/(m+4))h= (N^-1/(m+4)). Setting p^t=p^Nt p_t= p_N_t gives the score s^t(u)=∇ulnp^t(u) s_t(u)= _u p_t(u) in closed form. When the Planner observes context ξt∈ℝd _t ^d, the Simulator stores joint pairs (Ui,ξi)∈ℝm×ℝd(U_i, _i) ^m×R^d and fits a product-kernel KDE, p^(u,ξ)=1N∑i=1NKhu(u−Ui)Khξ(ξ−ξi), p(u,ξ)= 1N _i=1^NK_h_u(u-U_i)\,K_h_ξ(ξ- _i), where KhuK_h_u and KhξK_h_ξ are Gaussian kernels with bandwidths huh_u and hξh_ξ for the action and context dimensions respectively. The conditional density p^(u∣ξt)=p^(u,ξt)/∫p^(v,ξt)v p(u _t)= p(u, _t)/ p(v, _t)\,dv and its score ∇ulnp^(u∣ξt) _u p(u _t) are used in place of p^t(u) p_t(u) and s^t(u) s_t(u) throughout. For the KDE (11), the KL divergence in Assumption 2 decreases with NtN_t and is bounded by the integrated squared score error (Theorem 3). IV-B The Riemannian Structure of PPC This subsection establishes the geometric structure that drives PPC convergence and safety. In particular, it links local density curvature to optimization curvature through the PPC Hessian. Recall from Corollary 1 that PPC minimizes the free energy ℱ(u)=c(u)−βlnp^(u)F(u)=c(u)-β p(u). The PPC gradient iterates are defined by uk+1=uk−η∇ℱ(uk),k=0,1,…,K−1,u_k+1=u_k-η\, (u_k), k=0,1,…,K-1, (12) with step size η>0η>0, inner step count K, and warm start u0=ut−1u_0=u_t-1. The gradient is ∇ℱ=∇c−βs =∇ c-β s, and the Hessian is Hℱ(u)=∇2c(u)+βℐp^(u),H_F(u)=∇^2c(u)+β\,I_ p(u), (13) where ℐp^(u):=−∇2lnp^(u)I_ p(u):=-∇^2 p(u) is the local Fisher information of the learned density. Thus each PPC gradient step implicitly uses a Riemannian metric shaped by the curvature of the log-density barrier. The action space U, equipped with the positive-definite matrix field ℐp^(u)I_ p(u), forms a Riemannian manifold whose metric is the Fisher-Rao metric of the learned density p^(⋅∣ξ) p(· ξ). Gradient descent on ℱF is then equivalent to Euclidean gradient descent with a geometry-aware preconditioner. Near the manifold boundary, ℐp^I_ p is large and the effective step size shrinks, preventing the iterate from leaving the feasible region; in the interior, ℐp^I_ p is small and the step expands. This geometry-awareness is a structural advantage over CBF-QP [3], which uses a flat Euclidean metric in its safety projection. Given the learned density p p and a threshold α>0α>0, let ℳ^α M^α denote the learned α-superlevel set. The threshold α controls the conservatism of the learned safe set and is chosen by the Simulator (e.g., as a percentile of the KDE density values; see Algorithm 2). We now formalize the central geometric quantity of the paper. Definition 3 (Barrier Curvature). Let p p be a learned feasibility density and ℳ^α=u:p^(u)≥α M^α=\u: p(u)≥α\ its α-superlevel set. The barrier curvature κ and the curvature upper bound Λ are the tightest constants satisfying κI⪯ℐp^(u)⪯ΛI∀u∈ℳ^α,κ\,I\; \;I_ p(u)\; \; \,I ∀\,u∈ M^α, (14) where ℐp^(u)=−∇2lnp^(u)I_ p(u)=-∇^2 p(u) is the local Fisher information defined in (13). When the density is conditioned on context ξt _t, the curvature bounds depend on ξt _t and we write κ(ξt):=minu∈ℳ^tαλmin(ℐp^(⋅|ξt)(u)).κ( _t)\;:=\; _u∈ M_t^α\, _ \! (I_ p(·| _t)(u) ). When the context is suppressed (as in the single-step analysis of this section), κ refers to the barrier curvature of whichever density is in use. The barrier curvature κ quantifies how sharply the log-density increases toward the interior of ℳ^α M^α; it is the analogue, in our score-based setting, of the barrier self-concordance parameter in interior-point methods. Its upper counterpart Λ controls the smoothness of the free energy ℱF. Assumption 3 (Local Log-Concavity). κ>0κ>0 (equivalently, −lnp^- p is strictly convex on ℳ^α M^α), Λ<∞ <∞, and ℳ^α M^α is a single connected component of the superlevel set u:p^(u)≥α\u: p(u)≥α\ that is convex (equivalently, α exceeds the density value at every saddle point of p p so that the component lies in a single basin). The cost c is convex, LcL_c-smooth, and has bounded gradient ‖∇c‖∞≤Gc\|∇ c\|_∞≤ G_c. The condition κ>0κ>0 requires −lnp^- p to be strictly convex on ℳ^α M^α. For the Gaussian KDE (11) used throughout this paper, p^(u)=1N∑i=1NKh(u−Ui) p(u)= 1N _i=1^NK_h(u-U_i) is a mixture of Gaussians, and −lnp^- p is smooth everywhere. The single-basin/convexity clause is not automatic for mixtures: a Gaussian mixture can have several modes, and if α is taken below the density at an inter-mode saddle, p^≥α\ p≥α\ becomes a disjoint union of convex components with the Hessian indefinite between them. Choosing α above every saddle value (equivalently, setting α as a high percentile of observed density values) isolates one basin; within that basin p p is log-concave, −lnp^- p is strictly convex, and the superlevel set is convex. The Hessian ℐp^(u)I_ p(u) is then positive definite throughout, and κ>0κ>0 holds with κ determined by the narrowest bottleneck of the set. For the Gaussian KDE, Λ=1/h2 =1/h^2 where h is the bandwidth, so Λ<∞ <∞ holds automatically. Both κ and Λ are computable from the KDE in closed form, and the single-basin condition can be checked online by locating the critical points of p p, so Assumption 3 is verifiable online by the Simulator at each step. The following lemma establishes the landscape properties of ℱF that underpin the convergence theory. Lemma 1 (Free Energy Landscape). Under Assumption 3, let μ:=βκμ:=βκ and L:=Lc+βΛL:=L_c+β . The PPC free energy ℱ(u)=c(u)−βlnp^(u)F(u)=c(u)-β p(u) satisfies: (a) μ-strong convexity: Hℱ(u)⪰μIH_F(u) μ\,I for all u∈ℳ^αu∈ M^α; (b) L-smoothness: Hℱ(u)⪯LIH_F(u) L\,I for all u∈ℳ^αu∈ M^α; (c) local uniqueness: ℱF has a unique minimizer on the compact convex set ℳ^α M^α, denoted u∗:=argminu∈ℳ^αℱ(u)u^*:= _u∈ M^αF(u), and at most one unconstrained critical point of ℱF lies in ℳ^α M^α; when such a point exists it equals u∗u^* (equivalently, u∗∈int(ℳ^α)u^* ( M^α) and ∇ℱ(u∗)=0 (u^*)=0); (d) Polyak–Łojasiewicz inequality: ‖∇ℱ(u)‖2≥2μ(ℱ(u)−ℱ(u∗))\| (u)\|^2≥ 2μ (F(u)-F(u^*) ) for all u∈ℳ^αu∈ M^α. The condition number χ:=L/μ=(Lc+βΛ)/(βκ)χ:=L/μ=(L_c+β )/(βκ) governs the convergence rate of gradient-based iterates on ℱF. Note that (c) asserts uniqueness only within ℳ^α M^α; since p p is a Gaussian mixture, ℱF need not be convex outside ℳ^α M^α and may admit additional critical points there. Proof. For every u∈ℳ^αu∈ M^α, the Hessian decomposes as Hℱ(u)=∇2c(u)+βℐp^(u)H_F(u)=∇^2c(u)+β\,I_ p(u) by (13). Part (a). Since c is convex, ∇2c(u)⪰0∇^2c(u) 0. By Assumption 3, ℐp^(u)⪰κII_ p(u) κ\,I. Hence Hℱ(u)⪰ 0+βκI=μI.H_F(u)\; \;0+βκ\,I\;=\;μ\,I. (15) Part (b). By LcL_c-smoothness of c, ∇2c(u)⪯LcI∇^2c(u) L_c\,I. By Assumption 3, ℐp^(u)⪯ΛII_ p(u) \,I. Hence Hℱ(u)⪯LcI+βΛI=LI.H_F(u)\; \;L_c\,I+β \,I\;=\;L\,I. Part (c). ℳ^α M^α is closed (as a superlevel set of the continuous function p p) and bounded (since −lnp^- p is coercive for the Gaussian KDE, the superlevel set p^≥α\ p≥α\ is contained in a ball of finite radius), hence compact, and convex by Assumption 3. ℱF is μ-strongly convex on this set by part (a); a continuous strongly convex function on a compact convex set attains its minimum at exactly one point, which we denote u∗u^*. Strong convexity further implies that any u′∈ℳ^αu ∈ M^α with ∇ℱ(u′)=0 (u )=0 satisfies ℱ(u)≥ℱ(u′)+μ2‖u−u′‖2F(u) (u )+ μ2\|u-u \|^2 for all u∈ℳ^αu∈ M^α, forcing u′=u∗u =u^*. Part (d). Fix u∈ℳ^αu∈ M^α. By Taylor expansion and (15), for any v∈ℳ^αv∈ M^α: ℱ(v)≥ℱ(u)+⟨∇ℱ(u),v−u⟩+μ2‖v−u‖2F(v) (u)+ (u),v-u + μ2\|v-u\|^2. Minimizing the right-hand side over v (unconstrained quadratic) yields the lower bound ℱ(u∗)≥ℱ(u)−12μ‖∇ℱ(u)‖2F(u^*) (u)- 12μ\| (u)\|^2. Rearranging gives the Polyak–Łojasiewicz inequality. ∎ The landscape established by Lemma 1 directly yields a geometric convergence guarantee. The strong convexity provided by the barrier curvature κ forces PPC iterates toward the unique minimizer at a rate controlled by the condition number χ. Theorem 1 (Score-Descent Contraction). Set μ:=βκμ:=βκ, L:=Lc+βΛL:=L_c+β , χ:=L/μχ:=L/μ. Let u∗:=argminu∈ℳ^αℱ(u)u^*:= _u∈ M^αF(u) be the restricted minimizer of Lemma 1(c). Suppose Assumption 3 holds, u∗∈int(ℳ^α)u^* ( M^α), and the warm start satisfies u0∈ℳ^α,‖u0−u∗‖≤dist(u∗,∂ℳ^α).u_0∈ M^α, \|u_0-u^*\| (u^*,∂ M^α). Then the PPC iterates uk\u_k\ of (12) with step size η≤1/Lη≤ 1/L satisfy uk∈ℳ^αu_k∈ M^α for all k, together with the following contractions for every K≥0K≥ 0: (i) Function-value contraction. ℱ(uK)−ℱ(u∗)≤(1−ημ)K[ℱ(u0)−ℱ(u∗)].F(u_K)-F(u^*)\;≤\;(1-ημ)^K [F(u_0)-F(u^*) ]. (16) (i) Iterate contraction. ‖uK−u∗‖≤χ(1−ημ)K/2‖u0−u∗‖.\|u_K-u^*\|\;≤\; χ\,(1-ημ)^K/2\,\|u_0-u^*\|. (17) Proof. Iterates remain in ℳ^α M^α. We first show that uk∈ℳ^αu_k∈ M^α for every k so that Lemma 1 applies at each iterate. By μ-strong convexity and L-smoothness of ℱF on ℳ^α M^α, for η≤1/Lη≤ 1/L the standard non-expansive bound ‖uk+1−u∗‖2≤(1−ημ)2‖uk−u∗‖2\|u_k+1-u^*\|^2≤(1-ημ)^2\|u_k-u^*\|^2 holds whenever uk∈ℳ^αu_k∈ M^α [8] (by the contractivity of I−η∇2ℱI-η∇^2F whose eigenvalues lie in [1−ηL,1−ημ]⊂[0,1−ημ][1-η L,1-ημ]⊂[0,1-ημ]), i.e., iterate distance to u∗u^* is monotonically non-increasing. Since u0∈ℳ^αu_0∈ M^α and ‖u0−u∗‖≤dist(u∗,∂ℳ^α)\|u_0-u^*\| (u^*,∂ M^α), every uku_k lies in the closed Euclidean ball B(u∗,‖u0−u∗‖)⊆ℳ^αB(u^*,\|u_0-u^*\|) M^α, using that ℳ^α M^α is convex (Assumption 3) and contains this ball by construction. A simple induction then confirms uk∈ℳ^αu_k∈ M^α for all k≥0k≥ 0. By the L-smoothness of ℱF (Lemma 1(b)), ℱ(uk+1)≤ℱ(uk)+⟨∇ℱ(uk),uk+1−uk⟩+L2‖uk+1−uk‖2.F(u_k+1)\;≤\;F(u_k)+ (u_k),\,u_k+1-u_k + L2\|u_k+1-u_k\|^2. Substituting the gradient step uk+1=uk−η∇ℱ(uk)u_k+1=u_k-η\, (u_k), ℱ(uk+1) (u_k+1) ≤ℱ(uk)−η‖∇ℱ(uk)‖2+η2L2‖∇ℱ(uk)‖2 (u_k)-η\| (u_k)\|^2+ η^2L2\| (u_k)\|^2 =ℱ(uk)−η(1−ηL2)‖∇ℱ(uk)‖2. =F(u_k)-η (1- η L2 )\| (u_k)\|^2. (18) For η≤1/Lη≤ 1/L, the factor 1−ηL/2≥1/2>01-η L/2≥ 1/2>0, so every gradient step strictly decreases ℱF. The Polyak–Łojasiewicz inequality from Lemma 1(d) gives ‖∇ℱ(uk)‖2≥ 2μ(ℱ(uk)−ℱ(u∗)).\| (u_k)\|^2\;≥\;2μ (F(u_k)-F(u^*) ). (19) Combining these, substituting (19) into (18) with 1−ηL/2≥1/21-η L/2≥ 1/2, ℱ(uk+1)−ℱ(u∗) (u_k+1)-F(u^*) ≤(ℱ(uk)−ℱ(u∗))−ημ(ℱ(uk)−ℱ(u∗)) ≤ (F(u_k)-F(u^*) )-ημ (F(u_k)-F(u^*) ) =(1−ημ)(ℱ(uk)−ℱ(u∗)). =(1-ημ) (F(u_k)-F(u^*) ). (20) Iterating (20) over K steps yields part (i). For the iterate distance, by μ-strong convexity and L-smoothness, μ2‖u−u∗‖2≤ℱ(u)−ℱ(u∗)≤L2‖u−u∗‖2. μ2\|u-u^*\|^2\;≤\;F(u)-F(u^*)\;≤\; L2\|u-u^*\|^2. (21) Combining (21) with part (i) gives μ2‖uK−u∗‖2≤(1−ημ)K⋅L2‖u0−u∗‖2. μ2\|u_K-u^*\|^2\;≤\;(1-ημ)^K· L2\|u_0-u^*\|^2. Dividing by μ/2μ/2 and taking square roots yields part (i) with χ=L/μχ=L/μ. ∎ Theorem 1 ensures fast convergence of PPC iterates (warm-started within ℳ^α M^α) to the restricted minimizer u∗u^* of Lemma 1(c), but takes the interior condition u∗∈int(ℳ^α)u^* ( M^α) as a hypothesis. Whether this condition actually holds is a separate question, since u∗u^* may sit safely interior to ℳ^α M^α or merely on its boundary. Because the KDE p p is a sum of Gaussians, ℱF is finite on all of U and may possess additional critical points outside ℳ^α M^α; when β is small, the cost gradient dominates the score barrier at the boundary, so the minimum of ℱF over ℳ^α M^α is pinned to ∂ℳ^α∂ M^α, and unconstrained gradient descent on ℱF drifts out of the high-density region. Proposition 2 below supplies a sufficient condition on the stiffness, β>βcurv∗β>β^*_curv, under which u∗u^* is strictly interior to ℳ^α M^α, or equivalently an unconstrained critical point of ℱF, and hence by Lemma 1(c) the unique such point within ℳ^α M^α. This closes the interior hypothesis of Theorem 1. Proposition 2 relies only on Lemma 1 and Assumption 3, not on Theorem 1, so there is no circularity. Under Assumption 3, the negative log-density g(u)=−lnp^(u)g(u)=- p(u) is κ-strongly convex on ℳ^α M^α, so it has a unique minimizer u¯:=argmaxu∈ℳ^αp^(u) u:= _u∈ M^α p(u), the point of highest learned density. Define the level-set radius rα=infu∈∂ℳ^α‖u−u¯‖>0r_α= _u∈∂ M^α\|u- u\|>0. Proposition 2 (Critical Stiffness from Manifold Curvature). Under Assumption 3, let u∗:=argminu∈ℳ^αℱ(u)u^*:= _u∈ M^αF(u) denote the restricted minimizer from Lemma 1(c). Whenever the stiffness satisfies β>βcurv∗:=Gcκ⋅rα,β\;>\;β^*_curv\;:=\; G_cκ· r_α, (22) the following hold (the threshold β=βcurv∗β=β^*_curv is a limiting case in which the interior margin degenerates to zero): (i) Interior placement. u∗∈int(ℳ^α)u^* ( M^α), with displacement and margin bounds ‖u∗−u¯‖≤Gcβκ,dist(u∗,∂ℳ^α)≥rα−Gcβκ> 0.\|u^*- u\|\;≤\; G_cβκ, (u^*,∂ M^α)\;≥\;r_α- G_cβκ\;>\;0. (i) Unconstrained stationarity. ∇ℱ(u∗)=0 (u^*)=0, so u∗u^* is also an unconstrained critical point of ℱF, and it is the unique critical point of ℱF inside ℳ^α M^α. Proof. Let g(u):=−lnp^(u)g(u):=- p(u). By Assumption 3, g is κ-strongly convex on ℳ^α M^α with u¯ u its unique minimizer on the set. Since rα>0r_α>0 by its definition, u¯∈int(ℳ^α) u ( M^α), so the unconstrained stationarity ∇g(u¯)=0∇ g( u)=0 holds. Strong convexity gives the monotonicity inequality ⟨∇g(u)−∇g(u¯),u−u¯⟩≥κ‖u−u¯‖2 ∇ g(u)-∇ g( u),\,u- u ≥κ\|u- u\|^2; combining with ∇g(u¯)=0∇ g( u)=0 and applying Cauchy–Schwarz to the left-hand side yields ‖∇g(u)‖≥κ‖u−u¯‖∀u∈ℳ^α,\|∇ g(u)\|\;≥\;κ\|u- u\| ∀\,u∈ M^α, (23) and hence ‖s^(u)‖=‖∇g(u)‖≥κ‖u−u¯‖\| s(u)\|=\|∇ g(u)\|≥κ\|u- u\| on ℳ^α M^α. Part (i): interior placement. By Lemma 1(c), u∗=argminu∈ℳ^αℱ(u)u^*= _u∈ M^αF(u) exists and is unique. We rule out the boundary case by contradiction. Suppose u∗∈∂ℳ^αu^*∈∂ M^α, so ‖u∗−u¯‖≥rα\|u^*- u\|≥ r_α. Since ℳ^α=u:g(u)≤−lnα M^α=\u:g(u)≤- α\ and g is smooth and convex, the outward unit normal to ℳ^α M^α at u∗u^* is n=∇g(u∗)/‖∇g(u∗)‖n=∇ g(u^*)/\|∇ g(u^*)\|. First-order optimality of the constrained minimizer u∗u^* requires that no feasible direction improves ℱF; equivalently, ⟨∇ℱ(u∗),n⟩≤0 (u^*),n ≤ 0. Expanding ∇ℱ=∇c+β∇g =∇ c+β∇ g, ⟨∇c(u∗),n⟩+β‖∇g(u∗)‖≤ 0. ∇ c(u^*),\,n +β\|∇ g(u^*)\|\;≤\;0. Using ⟨∇c(u∗),n⟩≥−‖∇c(u∗)‖≥−Gc ∇ c(u^*),n ≥-\|∇ c(u^*)\|≥-G_c and ‖∇g(u∗)‖≥κ‖u∗−u¯‖≥κrα\|∇ g(u^*)\|≥κ\|u^*- u\|≥κ r_α from (23), this forces βκrα≤Gcβκ r_α≤ G_c, i.e., β≤βcurv∗β≤β^*_curv, contradicting the hypothesis β>βcurv∗β>β^*_curv. Therefore u∗∈int(ℳ^α)u^* ( M^α) for β>βcurv∗β>β^*_curv, and unconstrained stationarity ∇ℱ(u∗)=0 (u^*)=0 gives ∇c(u∗)=βs^(u∗)∇ c(u^*)=β\, s(u^*). Taking norms and using (23), ‖u∗−u¯‖≤‖s^(u∗)‖κ=‖∇c(u∗)‖βκ≤Gcβκ.\|u^*- u\|\;≤\; \| s(u^*)\|κ\;=\; \|∇ c(u^*)\|βκ\;≤\; G_cβκ. (24) Because u¯∈int(ℳ^α) u ( M^α) with nearest boundary point at distance rαr_α, the triangle inequality gives dist(u∗,∂ℳ^α)≥rα−Gc/(βκ)>0dist(u^*,∂ M^α)≥ r_α-G_c/(βκ)>0 for β>βcurv∗β>β^*_curv. Part (i): unconstrained stationarity and uniqueness inside ℳ^α M^α. Interior placement of u∗u^* makes the constraint inactive, so ∇ℱ(u∗)=0 (u^*)=0 holds in the unconstrained sense. Any other point v∈ℳ^αv∈ M^α with ∇ℱ(v)=0 (v)=0 would, by Lemma 1(c), satisfy v=u∗v=u^*. Thus u∗u^* is the unique critical point of ℱF within ℳ^α M^α, although ℱF may possess additional critical points in ∖ℳ^αU M^α where log-concavity need not hold. ∎ The key mechanism is that strong convexity anchors u∗u^* near u¯ u, which is known to be interior to ℳ^α M^α with margin rαr_α. Above the threshold βcurv∗β^*_curv, the cost-driven displacement Gc/(βκ)G_c/(βκ) stays below rαr_α, so u∗u^* cannot reach ∂ℳ^α∂ M^α; equivalently, the inward pull of the score barrier dominates the outward pull of the cost gradient everywhere on ∂ℳ^α∂ M^α, ruling out boundary optima. The critical stiffness βcurv∗β^*_curv depends on the barrier curvature κ and the level-set radius rαr_α, both of which are properties of the learned density, rather than on a Lipschitz bound on the unknown dynamics. Both κ and rαr_α can be estimated online from the KDE (see Section˜V-C). IV-C Dynamic Regret under Manifold Drift When the manifold is time-varying, ℱtF_t changes at every step. We use the dynamic regret as our performance metric, which is defined as the cumulative suboptimality relative to the moving restricted minimizer ut∗:=argminu∈ℳ^tαℱt(u)u^*_t:= _u∈ M_t^αF_t(u). Under β>βcurv∗β>β^*_curv, Proposition 2 further identifies ut∗u^*_t with the unique unconstrained critical point of ℱtF_t in ℳ^tα M_t^α, so stationarity ∇ℱt(ut∗)=0 _t(u^*_t)=0 is available for the analysis. We first bound the sensitivity of ut∗u^*_t to score changes, then derive the regret bound. Lemma 2 (Comparator Sensitivity). Under the conditions of Lemma 1 and β>βcurv∗β>β^*_curv (so each ut∗∈int(ℳ^tα)u^*_t ( M_t^α) by Proposition 2), suppose ut−1∗∈ℳ^tαu^*_t-1∈ M_t^α for every t. Then the comparator path length T:=∑t=1T‖ut∗−ut−1∗‖C_T:= _t=1^T\|u^*_t-u^*_t-1\| satisfies T≤1κ∑t=1Tsupu‖s^t(u)−s^t−1(u)‖.C_T\;≤\; 1κ _t=1^T _u\| s_t(u)- s_t-1(u)\|. (25) Proof. Interior stationarity ∇ℱt(ut∗)=0 _t(u^*_t)=0 holds at each step by Proposition 2. From ∇ℱt−1(ut−1∗)=0 _t-1(u^*_t-1)=0, i.e., ∇c(ut−1∗)=βs^t−1(ut−1∗)∇ c(u^*_t-1)=β\, s_t-1(u^*_t-1), evaluate the gradient of ℱtF_t at the old minimizer: ‖∇ℱt(ut−1∗)‖ \| _t(u^*_t-1)\| =‖∇c(ut−1∗)−βs^t(ut−1∗)‖ =\|∇ c(u^*_t-1)-β\, s_t(u^*_t-1)\| =β‖s^t−1(ut−1∗)−s^t(ut−1∗)‖ =β\| s_t-1(u^*_t-1)- s_t(u^*_t-1)\| ≤βsupu‖s^t(u)−s^t−1(u)‖. ≤β _u\| s_t(u)- s_t-1(u)\|. (26) By μ-strong convexity of ℱtF_t (Lemma 1(a)) together with ∇ℱt(ut∗)=0 _t(u^*_t)=0 (interior stationarity at step t), ‖∇ℱt(ut−1∗)‖ \| _t(u^*_t-1)\| =‖∇ℱt(ut−1∗)−∇ℱt(ut∗)‖ =\| _t(u^*_t-1)- _t(u^*_t)\| ≥μ‖ut−1∗−ut∗‖=βκ‖ut−1∗−ut∗‖. ≥μ\,\|u^*_t-1-u^*_t\|\;=\;βκ\,\|u^*_t-1-u^*_t\|. (27) Combining (26) and (IV-C) yields ‖ut∗−ut−1∗‖≤supu‖s^t−s^t−1‖/κ\|u^*_t-u^*_t-1\|≤ _u\| s_t- s_t-1\|/κ. Summing over t gives (25). ∎ The hypothesis ut−1∗∈ℳ^tαu^*_t-1∈ M_t^α is mild. By Proposition 2, ut−1∗u^*_t-1 sits interior to ℳ^t−1α M_t-1^α with margin rα−Gc/(βκ)>0r_α-G_c/(βκ)>0, so it remains in ℳ^tα M_t^α whenever the feasible set drifts by less than that margin, the natural regime in which dynamic regret is informative. Theorem 1 treats a fixed manifold. In practice, the feasibility set drifts over time. The following result quantifies how the cumulative suboptimality depends on both the contraction rate and the speed of manifold drift, measured through the score variation T:=∑t=1Tsupu‖s^t(u)−s^t−1(u)‖2V_T:= _t=1^T _u\| s_t(u)- s_t-1(u)\|^2. Theorem 2 (Dynamic Regret under Manifold Drift). Let ut\u_t\ be the output of PPC with K inner gradient steps (step size η=1/Lη=1/L) warm-started from ut−1u_t-1, and let ut∗:=argminu∈ℳ^tαℱt(u)u_t^*:= _u∈ M_t^αF_t(u) be the restricted minimizer of ℱt(u):=c(u)−βlnp^t(u)F_t(u):=c(u)-β p_t(u). Assume, for every t: (i) Assumptions 1–3 hold for p^t p_t; (i) β>βcurv∗β>β^*_curv; (i) (Warm-start margin) Let R∗:=rα−Gc/(βκ)>0R^*:=r_α-G_c/(βκ)>0 denote the interior margin from Proposition 2. The initial distance and consecutive comparator drift are bounded by half this margin, i.e., ‖u0−u0∗‖≤R∗/2\|u_0-u^*_0\|≤ R^*/2 and ‖ut∗−ut−1∗‖≤R∗/2\|u^*_t-u^*_t-1\|≤ R^*/2 for every t; (iv) K is large enough that ρK:=χ(1−1/χ)K/2≤1/2 _K:= χ(1-1/χ)^K/2≤ 1/2 (suffices to take K≥χln(4χ)K≥χ (4χ), where χ:=L/μχ:=L/μ). Then the dynamic regret DRegT:=∑t=1T[ℱt(ut)−ℱt(ut∗)]DReg_T:= _t=1^T[F_t(u_t)-F_t(u^*_t)] satisfies DRegT≤LρK2D021−ρK2+LρK2T(1−ρK)2κ2,DReg_T\;≤\; L _K^2\,D_0^21- _K^2\;+\; L _K^2\,V_T(1- _K)^2\,κ^2, (28) where D0:=‖u0−u0∗‖D_0:=\|u_0-u^*_0\| and T:=∑t=1Tsupu‖s^t(u)−s^t−1(u)‖2V_T:= _t=1^T _u\| s_t(u)- s_t-1(u)\|^2. Assumptions (i), (i), and (iv) are standard, namely (i) and (i) are the per-step regularity and strict-stiffness conditions of Theorem 1 and Proposition 2, while (iv) is the usual choice of K making each inner PPC pass at least a factor-of-two contraction. The warm-start margin (i) is the natural slow-drift condition needed to keep the warm start inside the current attraction ball. Splitting the interior margin R∗R^* into a half-budget for the initial tracking error and a half-budget for the consecutive comparator shift is exactly what the induction at the start of the proof requires to conclude ‖ut−1−ut∗‖≤R∗≤Rt\|u_t-1-u^*_t\|≤ R^*≤ R_t at every step, so that Theorem 1(i) applies; by Lemma 2, the drift bound translates to supu‖s^t−s^t−1‖≤κR∗/2 _u\| s_t- s_t-1\|≤κ R^*/2, which is the regime in which dynamic regret is informative. The bound in (28) consists of two terms. The first, LρK2D02/(1−ρK2)L _K^2D_0^2/(1- _K^2), is a transient that captures the cost of starting away from the equilibrium and vanishes geometrically with the number of inner steps K through ρK _K. The second, LρK2T/[(1−ρK)2κ2]L _K^2V_T/[(1- _K)^2κ^2], is the steady-state cost of tracking a moving manifold, proportional to the cumulative score variation TV_T and inversely proportional to κ2κ^2. In classical online control, the analogous drift sensitivity depends on a Lipschitz constant of the known dynamics; here κ plays that role despite ftf_t being unknown, because all manifold information is mediated through the learned density. Proof. Define the tracking error et:=‖ut−ut∗‖e_t:=\|u_t-u^*_t\| and comparator shift dt:=‖ut∗−ut−1∗‖d_t:=\|u^*_t-u^*_t-1\|. By induction on t, hypotheses (i)–(iv) keep ‖ut−1−ut∗‖≤et−1+dt≤R∗/2+R∗/2=R∗≤Rt\|u_t-1-u^*_t\|≤ e_t-1+d_t≤ R^*/2+R^*/2=R^*≤ R_t (using Rt:=dist(ut∗,∂ℳ^tα)≥R∗R_t:=dist(u^*_t,∂ M_t^α)≥ R^* from Proposition 2), so ut−1u_t-1 lies in the attraction ball of Theorem 1 around ut∗u^*_t at every step. Applying Theorem 1(i) with K gradient steps on ℱtF_t starting from ut−1u_t-1 gives et≤ρK‖ut−1−ut∗‖≤ρK(et−1+dt),e_t\;≤\; _K\,\|u_t-1-u^*_t\|\;≤\; _K\,(e_t-1+d_t), (29) where the second inequality is the triangle inequality. Iterating (29) leads to (using summation index σ, with τ reserved for the manifold reach of Theorem 3) et≤ρKtD0+∑σ=1tρKt−σ+1dσ.e_t\;≤\; _K^t\,D_0\;+\; _σ=1^t _K^t-σ+1\,d_σ. (30) By the Cauchy–Schwarz inequality applied to the weighted sum (with weights ρKt−σ+1 _K^t-σ+1 and ∑σ=1tρKt−σ+1≤ρK/(1−ρK) _σ=1^t _K^t-σ+1≤ _K/(1- _K)), (∑σ=1tρKt−σ+1dσ)2≤ρK1−ρK∑σ=1tρKt−σ+1dσ2. ( _σ=1^t _K^t-σ+1d_σ )^2\;≤\; _K1- _K _σ=1^t _K^t-σ+1d_σ^2. (31) Summing over t and exchanging the order of summation, ∑t=1Tet2≤2ρK2D021−ρK2+2ρK2(1−ρK)2∑t=1Tdt2. _t=1^Te_t^2\;≤\; 2 _K^2D_0^21- _K^2\;+\; 2 _K^2(1- _K)^2 _t=1^Td_t^2. (32) Finally, by L-smoothness of ℱtF_t, ℱt(ut)−ℱt(ut∗)≤L2et2F_t(u_t)-F_t(u^*_t)≤ L2\,e_t^2. Substituting Lemma 2 (dt≤supu‖s^t−s^t−1‖/κd_t≤ _u\| s_t- s_t-1\|/κ) into (32) yields DRegT≤L2∑t=1Tet2≤LρK2D021−ρK2+LρK2T(1−ρK)2κ2.∎DReg_T≤ L2 _t=1^Te_t^2≤ L _K^2D_0^21- _K^2+ L _K^2\,V_T(1- _K)^2\,κ^2. (33) The regret is stated in terms of the free energy ℱtF_t because PPC descends on ℱtF_t, not on c directly. A cost regret bound follows since for any ut,ut∗∈ℳ^tαu_t,u^*_t∈ M_t^α the log-density difference is bounded, and hence ∑t[c(ut)−c(ut∗)] _t[c(u_t)-c(u^*_t)] differs from DRegTDReg_T by at most an additive term proportional to β. Safety is treated separately in Theorem 4. The barrier curvature κ plays a dual role in (28). It enters the contraction factor ρK _K through μ=βκμ=βκ, so larger κ accelerates convergence, and it appears in the drift sensitivity through 1/κ21/κ^2, so larger κ also reduces the impact of manifold motion. Both effects favor a density with sharp boundaries. In practice, κ is a property of the learned density and cannot be tuned directly. The controllable quantity is βt _t, which the adaptive schedule (58) sets to βcurv∗(1+C/Nt)β^*_curv(1+C/ N_t). When NtN_t is small the schedule inflates βt _t above the critical threshold, favoring safety over cost optimality. As NtN_t grows and the density estimate improves, βt _t decreases toward βcurv∗β^*_curv, gradually releasing the Planner to exploit more of the feasible space. IV-D Online Score Convergence This subsection provides the statistical bridge from the fixed-density analysis in Sections IV-A and IV-B to the online setting where the KDE estimate p^t p_t from (11) evolves with new simulator samples. Algorithms 1–2 present the controller using the marginal KDE (no context); the convergence guarantees below are stated for this case, and the contextual case follows analogously with m replaced by m+dm+d. Lemma 3 controls pointwise density and gradient errors, and Theorem 3 converts those controls into the integrated squared score error ϵt _t. Lemma 3 (KDE Uniform Convergence). Let p∗p^* be a C2C^2 density on a compact domain Ω⊂ℝm ^m satisfying infΩp∗>0 _ p^*>0 and ‖p∗‖C2<∞\|p^*\|_C^2<∞. Let p^N p_N be the Gaussian KDE with bandwidth h=Θ(N−1/(m+4))h= (N^-1/(m+4)) built from N i.i.d. samples from p∗p^*. Then with probability at least 1−δ1-δ: ‖p^N−p∗‖∞ \| p_N-p^*\|_∞ ≤C1N−2/(m+4)ln(N/δ), \;≤\;C_1\,N^-2/(m+4) (N/δ), (34) ‖∇p^N−∇p∗‖∞ \|∇ p_N-∇ p^*\|_∞ ≤C2N−1/(m+4)ln(N/δ), \;≤\;C_2\,N^-1/(m+4) (N/δ), (35) where C1,C2C_1,C_2 depend on ‖p∗‖C2\|p^*\|_C^2, m, and vol(Ω)vol( ). Proof. The KDE is p^N(u)=1N∑i=1NKh(u−Ui) p_N(u)= 1N _i=1^NK_h(u-U_i) with Kh(z)=h−mK(z/h)K_h(z)=h^-mK(z/h) and K the standard Gaussian kernel. By Taylor expansion of p∗p^* under convolution, the bias satisfies ‖[p^N]−p∗‖∞=O(h2)\|E[ p_N]-p^*\|_∞=O(h^2), and differentiating under the integral gives ‖∇[p^N]−∇p∗‖∞=O(h2)\| [ p_N]-∇ p^*\|_∞=O(h^2). The pointwise variance of p^N(u) p_N(u) is Var[p^N(u)]≤C/(Nhm)Var[ p_N(u)]≤ C/(Nh^m), and for the gradient Var[∇p^N(u)]≤C/(Nhm+2)Var[∇ p_N(u)]≤ C/(Nh^m+2), since ∇Kh∇ K_h introduces an extra 1/h1/h factor. A bounded-difference inequality over an ε -net of Ω (cardinality O((diam(Ω)/ε)m)O((diam( )/ )^m)) with a union bound yields, for any δ∈(0,1)δ∈(0,1) and with probability at least 1−δ1-δ, supu∈Ω|p^N(u)−[p^N(u)]|≤O(ln(N/δ)Nhm), _u∈ | p_N(u)-E[ p_N(u)]|≤ O\! ( (N/δ)Nh^m ), and the analogous bound for ∇p^N∇ p_N with hmh^m replaced by hm+2h^m+2. Setting h=Θ(N−1/(m+4))h= (N^-1/(m+4)) balances the density bias O(N−2/(m+4))O(N^-2/(m+4)) with its standard deviation O(lnN⋅N−2/(m+4))O( N· N^-2/(m+4)), yielding (34) (the logarithmic factor is absorbed into C1C_1). For the gradient, the bias is still O(N−2/(m+4))O(N^-2/(m+4)) but the standard deviation is O(lnN⋅N−1/(m+4))O( N· N^-1/(m+4)), which dominates [36], giving (35). ∎ Theorem 3 (Score Tracking Bound). Let p^t p_t be a Gaussian KDE with bandwidth h=Θ(Nt−1/(m+4))h= (N_t^-1/(m+4)) built from NtN_t cumulative feasibility samples on a manifold with reach τ>0τ>0. Then the integrated squared score error satisfies [∫ℳ^tα‖s^t(u)−st∗(u)‖2p^t(u)u]≤C(m,τ,α)Nt2/(m+4),E\! [ _ M_t^α\| s_t(u)-s^*_t(u)\|^2\, p_t(u)\,du ]\;≤\; C(m,τ,α)N_t^2/(m+4), (36) where C(m,τ,α)C(m,τ,α) depends polynomially on m, 1/τ1/τ, and 1/α1/α; the 1/α1/α dependence enters through the denominator p^t≥α p_t≥α on ℳ^tα M_t^α via pmin:=α/2p_ :=α/2. Proof. By the chain rule, s^(u)=∇ulnp^(u)=∇p^(u)/p^(u) s(u)= _u p(u)=∇ p(u)/ p(u) and similarly s∗(u)=∇p∗(u)/p∗(u)s^*(u)=∇ p^*(u)/p^*(u). Subtracting, we obtain s^(u)−s∗(u)=∇(p^−p∗)(u)p^(u)⏟=:A(u)−p^(u)−p∗(u)p^(u)s∗(u)⏟=:B(u), s(u)-s^*(u)\;=\; ∇( p-p^*)(u) p(u)_=:\,A(u)\;-\; p(u)-p^*(u) p(u)\,s^*(u)_=:\,B(u), (37) where the identity follows from a direct computation: A(u)−B(u) A(u)-B(u) =∇p^−∇p∗p^−(p^−p∗)s∗p = ∇ p-∇ p^* p- ( p-p^*)\,s^* p =∇p^⋅p∗−∇p∗⋅p^p^⋅p∗=s^−s∗. = ∇ p· p^*-∇ p^*· p p· p^*\;=\; s-s^*. On ℳ^tα M_t^α we have p^t≥α p_t≥α, and ‖p^t−pt∗‖∞=O(N−2/(m+4))\| p_t-p^*_t\|_∞=O(N^-2/(m+4)) by Lemma 3; therefore, for N large enough that ‖p^t−pt∗‖∞≤α/2\| p_t-p^*_t\|_∞≤α/2, we obtain pt∗≥α/2p^*_t≥α/2 on ℳ^tα M_t^α. Set pmin:=α/2>0p_ :=α/2>0; then p^t≥pmin p_t≥ p_ on ℳ^tα M_t^α as well (in fact p^t≥α=2pmin p_t≥α=2p_ ). Note that this is not an appeal to the oracle infℳtpt∗ _M_tp^*_t bound of Definition 2, which fails near ∂ℳt _t; the threshold α>0α>0 itself guarantees the bound on the learned set. The score s∗s^* is bounded on the compact set ℳ^tα M_t^α: ‖s∗‖∞≤Smax<∞\|s^*\|_∞≤ S_ <∞. Therefore: ‖A(u)‖ \|A(u)\| ≤2pmin‖∇(p^−p∗)‖∞≤2C2pminN−1/(m+4), ≤ 2p_ \|∇( p-p^*)\|_∞\;≤\; 2C_2p_ \,N^-1/(m+4), (38) ‖B(u)‖ \|B(u)\| ≤2Smaxpmin‖p^−p∗‖∞≤2C1SmaxpminN−2/(m+4). ≤ 2S_ p_ \| p-p^*\|_∞\;≤\; 2C_1S_ p_ \,N^-2/(m+4). (39) By (37) and ‖A−B‖2≤2‖A‖2+2‖B‖2\|A-B\|^2≤ 2\|A\|^2+2\|B\|^2: ∫ℳ^tα‖s^−s∗‖2p^u _ M_t^α\| s-s^*\|^2 p\,du ≤2∫ℳ^tα(‖A‖2+‖B‖2)p^u ≤ 2 _ M_t^α (\|A\|^2+\|B\|^2 ) p\,du ≤2(‖A‖∞2+‖B‖∞2)≤C(m,τ,α)N2/(m+4), ≤ 2 (\|A\|_∞^2+\|B\|_∞^2 )\;≤\; C(m,τ,α)N^2/(m+4), (40) where C(m,τ,α)C(m,τ,α) absorbs the constants from (38)–(39) and depends on pmin=α/2p_ =α/2, SmaxS_ , ‖p∗‖C2\|p^*\|_C^2, and vol(ℳ^tα)vol( M_t^α). These quantities depend polynomially on m, 1/τ1/τ (through the manifold geometry), and 1/α1/α (through the 1/pmin21/p_ ^2 prefactors of (38)–(39)). The bound (40) holds on the event of Lemma 3, which has probability at least 1−δ1-δ. The expectation in (36) follows by choosing δ=N−2/(m+4)δ=N^-2/(m+4) and noting that the complementary event contributes at most O(δ⋅Smax2)=O(N−2/(m+4))O(δ· S_ ^2)=O(N^-2/(m+4)), which is absorbed into C(m,τ,α)C(m,τ,α). The rate N−2/(m+4)N^-2/(m+4) is the minimax-optimal rate for nonparametric gradient estimation in m dimensions [36]. The plug-in score estimator inherits this rate because s^=∇p^/p s=∇ p/ p and the gradient error in the numerator dominates the density error in the denominator. For comparison, the density mean integrated squared error (MISE), [∫(p^N−p∗)2u]E[ ( p_N-p^*)^2\,du], converges at the faster rate N−4/(m+4)N^-4/(m+4), reflecting the additional difficulty of derivative estimation. ∎ Remark 1. For m=2m=2 (our experiments), the integrated squared score error decays as N−1/3N^-1/3, which is slower than the density MISE rate N−2/3N^-2/3 but still sufficient for the safety bound of Theorem 4 to shrink monotonically as N grows (Figure˜9). IV-E Safety Bound via Score Descent We now compose the contraction guarantee (Theorem 1) with the statistical bound (Theorem 3) to obtain the main safety result. The bound is analogous to Input-to-State Safety (ISSf) [2] in that the distance to the safe set is bounded by a function of the estimation error ϵt _t, but the sensitivity is governed by κ rather than a Lipschitz constant of the dynamics. A supporting lemma first controls the set approximation error. Lemma 4 (Level-Set Stability). Fix a time step t. Let ℳ^tα:=u:p^t(u∣ξt)≥α M_t^α:=\u: p_t(u _t)≥α\ and let ϵt _t be the score error of Theorem 3. Assume the oracle density pt∗p^*_t (Definition 2) has nondegenerate inward-normal derivative in a boundary tube. Specifically, there exist c∂,r∂>0c_∂,r_∂>0 with r∂≤τr_∂≤τ such that ⟨∇upt∗(u∣ξt),n^(u)⟩≥c∂ _up^*_t(u _t),\, n(u) \;≥\;c_∂ (41) for all u satisfying dist(u,∂ℳt)≤r∂dist(u, _t)≤ r_∂, where n^(u) n(u) is the inward unit normal at the nearest boundary point. If α and ϵt _t are small enough that ϵt<αand(α+ϵt)/c∂≤r∂, _t<α (α+ _t)/c_∂≤ r_∂, (42) then dH(ℳ^tα,ℳt)≤α+ϵtc∂.d_H( M_t^α,\,M_t)\;≤\; α+ _tc_∂. (43) The inward-normal nondegeneracy (41) is a standard transversality condition that rules out flat-approach boundaries where ∇pt∗∇ p^*_t vanishes; it is what lets density perturbations translate linearly into boundary displacements, and it holds generically for smooth densities whose support has a C2C^2 boundary. The smallness condition (42) simply requires the level-set threshold α and the sup-norm error ϵt _t to stay inside the boundary tube of width r∂r_∂, which is automatic for N large enough since ϵt→0 _t→ 0. Proof. Lemma 3 gives ‖p^t−pt∗‖∞2=O(N−4/(m+4)lnN)\| p_t-p^*_t\|_∞^2=O(N^-4/(m+4) N), which is eventually dominated by ϵt=O(N−2/(m+4)) _t=O(N^-2/(m+4)); we use the resulting sup-norm bound ‖p^t−pt∗‖∞≤ϵt\| p_t-p^*_t\|_∞≤ _t. Write sα:=(α+ϵt)/c∂s_α:=(α+ _t)/c_∂; by (42), sα≤r∂s_α≤ r_∂. Step 1: ℳ^tα⊆ℳt M_t^α _t. If u∈ℳ^tαu∈ M_t^α, then pt∗(u)≥p^t(u)−ϵt≥α−ϵt>0p^*_t(u)≥ p_t(u)- _t≥α- _t>0 by (42), so u∈ℳtu _t (the oracle density vanishes outside ℳtM_t). Hence supu∈ℳ^tαdist(u,ℳt)=0 _u∈ M_t^αdist(u,M_t)=0. Step 2: every v∈ℳtv _t is within sαs_α of ℳ^tα M_t^α. Let v∈ℳtv _t. If dist(v,∂ℳt)≥sαdist(v, _t)≥ s_α, let ub=π∂ℳt(v)u_b= _ _t(v) (well-defined since sα≤r∂≤τs_α≤ r_∂≤τ) and parametrize the inward normal ray γ(s)=ub+sn^(v)γ(s)=u_b+s n(v) for s∈[0,sα]s∈[0,s_α], which stays in the r∂r_∂-tube. The fundamental theorem of calculus and (41) give pt∗(γ(sα))=∫0sα⟨∇pt∗(γ(s)),n^(v)⟩s≥c∂sα=α+ϵt.p^*_t(γ(s_α))\;=\; _0^s_α ∇ p^*_t(γ(s)), n(v) \,ds\;≥\;c_∂s_α\;=\;α+ _t. (44) Because pt∗p^*_t is increasing along the inward ray in the tube, monotonicity gives pt∗(v)≥pt∗(γ(sα))≥α+ϵtp^*_t(v)≥ p^*_t(γ(s_α))≥α+ _t whenever dist(v,∂ℳt)≥sαdist(v, _t)≥ s_α, hence p^t(v)≥pt∗(v)−ϵt≥α p_t(v)≥ p^*_t(v)- _t≥α and v∈ℳ^tαv∈ M_t^α, i.e., dist(v,ℳ^tα)=0dist(v, M_t^α)=0. If instead dist(v,∂ℳt)<sαdist(v, _t)<s_α, let ub=π∂ℳt(v)u_b= _ _t(v) and set v′:=ub+sαn^(v)v :=u_b+s_α n(v). By the same FTC argument, v′∈ℳ^tαv ∈ M_t^α, and ‖v−v′‖≤sα\|v-v \|≤ s_α, so dist(v,ℳ^tα)≤sαdist(v, M_t^α)≤ s_α. Combining Steps 1 and 2, dH(ℳ^tα,ℳt)≤sα=(α+ϵt)/c∂d_H( M_t^α,M_t)≤ s_α=(α+ _t)/c_∂. ∎ Composing the contraction guarantee (Theorem 1) with the statistical bound (Theorem 3) yields the main safety result, a safety guarantee whose ultimate bound depends on manifold geometry rather than a Lipschitz constant of unknown dynamics. Let NtN_t denote the cumulative number of feasibility samples available at time t, and let ϵt:=C(m,τ,α)/Nt2/(m+4) _t:=C(m,τ,α)/N_t^2/(m+4) denote the integrated score-error rate from Theorem 3. Theorem 4 (Contextual Safety Bound). Let ut\u_t\ be the output of contextual PPC with stiffness βt _t and K inner gradient steps warm-started from ut−1u_t-1, and let κ(ξt)κ( _t) denote the barrier curvature of p^t(⋅∣ξt) p_t(· _t) (Definition 3). Assume, for every t: (i) Assumptions 1–3 hold for p^t(⋅∣ξt) p_t(· _t); (i) βt>βcurv∗ _t>β^*_curv; (i) the inward-normal nondegeneracy (41) and smallness conditions (42) of Lemma 4 hold; (iv) the warm-start margin of Theorem 2(i) holds, so that ut−1u_t-1 lies in the attraction ball of Theorem 1 around ut∗u^*_t; (v) K is large enough that ρKdiam(ℳ^tα) _K\,diam( M_t^α) is negligible relative to the other terms of (45) (quantified in the proof), where ρK:=χ(1−1/χ)K/2 _K:= χ\,(1-1/χ)^K/2 is the iterate contraction factor of Theorem 1(i) with condition number χ:=L/μχ:=L/μ, μ=βtκ(ξt)μ= _tκ( _t), and L=Lc+βtΛL=L_c+ _t . Then dist(ut,ℳt)≤ϵt+αc∂⏟set approx.+Gcβtκ(ξt)⏟cost–safety,dist(u_t,M_t)\;≤\; _t+αc_∂_set approx.\;+\; G_c _tκ( _t)_cost--safety, (45) where ϵt:=C(m,τ,α)/Nt2/(m+4) _t:=C(m,τ,α)/N_t^2/(m+4) is the integrated squared score error from Theorem 3. Remark 2 (Dependencies and asymptotic tightening). The bound (45) separates a statistical set-approximation term from a deterministic cost–safety term. The score error ϵt _t depends on the cumulative sample count NtN_t and the ambient dimension m (or m+dm+d for the conditional KDE), not on ξt _t directly; the context enters the bound through κ(ξt)κ( _t), which reflects the local geometry of the conditional density. As Nt→∞N_t→∞, ϵt→0 _t→ 0 and the bound tightens to α/c∂+Gc/(βtκ(ξt))α/c_∂+G_c/( _tκ( _t)). Proof. Let u¯t:=argmaxu∈ℳ^tαp^t(u∣ξt) u_t:= _u∈ M_t^α p_t(u _t) denote the maximum-density point of the learned conditional density, which is interior to ℳ^tα M_t^α with dist(u¯t,∂ℳ^tα)≥rαdist( u_t,∂ M_t^α)≥ r_α (Proposition 2, applied to p^t(⋅∣ξt) p_t(· _t)). Decompose the distance via u¯t u_t using the triangle inequality: dist(ut,ℳt)≤‖ut−ut∗‖+‖ut∗−u¯t‖ (u_t,M_t)≤\|u_t-u^*_t\|+\|u^*_t- u_t\| +dist(u¯t,ℳt) +dist( u_t,M_t) =‖ut−ut∗‖⏟planner residual+‖ut∗−u¯t‖⏟cost–safety offset = \|u_t-u^*_t\|_planner residual+ \|u^*_t- u_t\|_cost--safety offset +dist(u¯t,ℳt)⏟set approximation. + dist( u_t,M_t)_set approximation. (46) Cost–safety offset. Applying Proposition 2 to the conditional density p^t(⋅∣ξt) p_t(· _t), whose barrier curvature is κ(ξt)κ( _t) by Definition 3, yields ‖ut∗−u¯t‖≤Gcβtκ(ξt),\|u^*_t- u_t\|\;≤\; G_c _t\,κ( _t), (47) and ensures ut∗∈int(ℳ^tα)u^*_t ( M_t^α) with dist(ut∗,∂ℳ^tα)≥rα−Gc/(βtκ(ξt))>0dist(u^*_t,∂ M_t^α)≥ r_α-G_c/( _tκ( _t))>0 whenever βt>βcurv∗ _t>β^*_curv. This is where the contextual curvature enters the bound, since the Gc/(βκ)G_c/(βκ) term of Proposition 2 measures the offset of the cost-penalized minimizer from the density peak, not a residual of the gradient iterates. Set approximation. Since u¯t∈ℳ^tα u_t∈ M_t^α, Lemma 4 gives dist(u¯t,ℳt)≤dH(ℳ^tα,ℳt)≤ϵt+αc∂.dist( u_t,M_t)\;≤\;d_H( M_t^α,M_t)\;≤\; _t+αc_∂. Planner residual. By Theorem 1(i), K inner gradient steps warm-started from ut−1u_t-1 yield ‖ut−ut∗‖≤ρK‖ut−1−ut∗‖,ρK=χ(1−1/χ)K/2.\|u_t-u^*_t\|\;≤\; _K\,\|u_t-1-u^*_t\|, _K= χ\,(1-1/χ)^K/2. Since ρK _K decays geometrically and ‖ut−1−ut∗‖≤Dℳ,t:=diam(ℳ^tα)<∞\|u_t-1-u^*_t\|≤ D_M,t:=diam( M_t^α)<∞ (the learned superlevel set is compact as the closed α-superlevel set of a Gaussian KDE, which vanishes at infinity; Dℳ,tD_M,t is a diameter, distinct from the initial-distance symbol D0=‖u0−u0∗‖D_0=\|u_0-u^*_0\| of Theorem 2), any tolerance δK>0 _K>0 is reached with K≥2χln(χDℳ,t/δK)K≥ 2χ ( χ\,D_M,t/ _K). Taking δK _K smaller than the interior margin rα−Gc/(βtκ(ξt))r_α-G_c/( _tκ( _t)) also guarantees ut∈ℳ^tαu_t∈ M_t^α, so the score-ascent safety filter of Algorithm 2 is inactive in steady state. Summing the three terms of (46) and taking K large enough that the planner residual is negligible relative to the other two terms yields (45). For finite K, an additive remainder ρKdiam(ℳ^tα) _K\,diam( M_t^α) is added to the right-hand side. ∎ The bound is PPC-specific in that the residual Gc/(βκ)G_c/(βκ) depends on the barrier curvature κ rather than a Lipschitz constant of the unknown dynamics ftf_t. As NtN_t grows, ϵt→0 _t→ 0 at rate Nt−2/(m+4)N_t^-2/(m+4) (Theorem 3), so the bound tightens to Gc/(βκ)G_c/(βκ). The adaptive schedule (58) can then reduce β toward βcurv∗β^*_curv, allowing the Planner to exploit more of the feasible space as estimation improves. The value of context. The safety bound (45) depends on the context ξt _t through ϵt _t, c∂c_∂, and κ(ξt)κ( _t). The first two are properties of the oracle density pt∗(⋅∣ξt)p^*_t(· _t) (complexity for the score-estimation rate and inward-normal nondegeneracy for the set-approximation constant) that any Planner inherits from the underlying geometry. The only ξt _t-dependent quantity the Planner itself controls is the barrier curvature of the density it uses in its free energy, namely κ(ξt)κ( _t) for a context-aware Planner and κmarg:=infu∈ℳ^tαλmin(−∇u2lnp^t(u)) _marg:= _u∈ M_t^α _ \! (- _u^2 p_t(u) ) (i.e., Definition 3 applied to the marginal density) for a context-blind Planner that replaces p^t(⋅∣ξt) p_t(· _t) with the marginal p^t(u):=∫p^t(u∣ξ)π(ξ)ξ p_t(u):= p_t(u ξ)\,π(ξ)\,dξ, where π(ξ)π(ξ) denotes the context prior—the ξ-marginal of the joint product-kernel KDE of Section IV-A, i.e., π(ξ)=1Nt∑i=1NtKhξ(ξ−ξi)π(ξ)= 1N_t _i=1^N_tK_h_ξ(ξ- _i)—and the distribution underlying every expectation π[⋅]E_π[\,·\,] over contexts below (to be distinguished from the Planner’s policy πt _t, which is indexed by time). Since κ enters the cost–safety residual Gc/(βtκ)G_c/( _tκ) inversely, the value of context reduces to comparing these two curvatures. Proposition 3 records the variational-inference decomposition that underlies the comparison, and Theorem 5 turns it into a quantitative lower bound on the gap between the safety guarantees. Proposition 3 (Contextual Curvature Gain). With π(ξ)π(ξ) the context prior defined above, let w(ξ∣u)∝p^t(u∣ξ)π(ξ)w(ξ u) p_t(u ξ)\,π(ξ) be the posterior over contexts given action u, and s^t(u∣ξ):=∇ulnp^t(u∣ξ) s_t(u ξ):= _u p_t(u ξ) the conditional score. The marginal log-density Hessian ∇u2lnp^t(u) _u^2 p_t(u) equals w(ξ∣u)[∇u2lnp^t(u∣ξ)]+Covw(ξ∣u)(s^t(u∣ξ)).E_w(ξ u)\! [ _u^2 p_t(u ξ) ]\;+\;Cov_w(ξ u)\! ( s_t(u ξ) ). (48) Proof. Write q(ξ):=p^t(u∣ξ)q(ξ):= p_t(u ξ), so that p^t(u)=π[q(ξ)] p_t(u)=E_π[q(ξ)] and w(ξ∣u)=q(ξ)π(ξ)/p^t(u)w(ξ u)=q(ξ)π(ξ)/ p_t(u). Using ∇uq(ξ)=q(ξ)s^t(u∣ξ) _uq(ξ)=q(ξ)\, s_t(u ξ) and differentiating lnp^t(u)=lnπ[q(ξ)] p_t(u)= _π[q(ξ)] once gives the marginal-score identity ∇ulnp^t(u)=π[∇uq(ξ)]p^t(u)=w(ξ∣u)[s^t(u∣ξ)]. _u p_t(u)\;=\; E_π[ _uq(ξ)] p_t(u)\;=\;E_w(ξ u)\! [ s_t(u ξ) ]. (49) Differentiating lnw(ξ∣u)=lnq(ξ)+lnπ(ξ)−lnp^t(u) w(ξ u)= q(ξ)+ π(ξ)- p_t(u) in u and substituting (49) yields the standard variational-inference identity ∇ulnw(ξ∣u)=s^t(u∣ξ)−w(ξ∣u)[s^t(u∣ξ)], _u w(ξ u)\;=\; s_t(u ξ)-E_w(ξ u)\! [ s_t(u ξ) ], (50) i.e., the centered conditional score. Differentiating (49) once more and writing the expectation explicitly as ∇ulnp^t(u)=∫w(ξ∣u)s^t(u∣ξ)ξ _u p_t(u)= w(ξ u)\, s_t(u ξ)\,dξ, the product rule under the integral sign gives ∇u2lnp^t(u) _u^2 p_t(u) =∫[∇uw(ξ∣u)]s^t(u∣ξ)⊤ξ⏟=:T1 = [ _uw(ξ u)]\, s_t(u ξ)^\! \,dξ_=:\,T_1 +∫w(ξ∣u)∇us^t(u∣ξ)ξ⏟=:T2. + w(ξ u)\, _u s_t(u ξ)\,dξ_=:\,T_2. For T1T_1, the chain-rule identity ∇uw(ξ∣u)=w(ξ∣u)∇ulnw(ξ∣u) _uw(ξ u)=w(ξ u)\, _u w(ξ u) combined with the posterior log-derivative (50) yields ∇uw(ξ∣u)=w(ξ∣u)(s^t(u∣ξ)−w[s^t(u∣ξ)]) _uw(ξ u)=w(ξ u) ( s_t(u ξ)-E_w[ s_t(u ξ)] ), so T1 T_1 =w[(s^t(u∣ξ)−w[s^t(u∣ξ)])s^t(u∣ξ)⊤] =E_w\! [ ( s_t(u ξ)-E_w[ s_t(u ξ)] )\, s_t(u ξ)^\! ] =w[s^t(u∣ξ)s^t(u∣ξ)⊤]−w[s^t(u∣ξ)]w[s^t(u∣ξ)]⊤ =E_w\! [ s_t(u ξ)\, s_t(u ξ)^\! ]-E_w[ s_t(u ξ)]\,E_w[ s_t(u ξ)]^\! =Covw(s^t(u∣ξ)), =Cov_w\! ( s_t(u ξ) ), where the second line pulls the ξ-independent factor w[s^t(u∣ξ)]E_w[ s_t(u ξ)] out of the expectation. For T2T_2, the conditional score satisfies ∇us^t(u∣ξ)=∇u2lnp^t(u∣ξ) _u s_t(u ξ)= _u^2 p_t(u ξ) by definition, so T2=w[∇u2lnp^t(u∣ξ)]T_2=E_w\! [ _u^2 p_t(u ξ) ]. Adding T1T_1 and T2T_2 yields (48). ∎ Proposition 3 is the classical “total-Hessian” decomposition: the marginal log-density Hessian equals the posterior-averaged conditional Hessian plus the posterior score covariance. Flipping signs to the positive-semidefinite curvature −∇u2lnp^t⪰0- _u^2 p_t 0 used elsewhere in the paper, this says that marginalizing subtracts Covw⪰0Cov_w 0 from the expected conditional curvature, flattening the barrier pointwise in the Löwner order; the gap is exactly the posterior variance of the conditional score, which vanishes when the posterior concentrates on a single context and grows as different contexts push the conditional score in substantially different directions. Theorem 5 converts this PSD decomposition into a scalar safety-guarantee gap via Weyl’s inequality. Theorem 5 (Contextual Safety Gap). Let BctxB_ctx and BmargB_marg denote the Theorem 4 steady-state safety bounds for the context-aware and context-blind Planners, both sharing the statistical error ϵt _t, boundary constant c∂c_∂, and feasible set ℳ^tα M_t^α, so that κmarg:=infu∈ℳ^tαλmin(−∇u2lnp^t(u)) _marg:= _u∈ M_t^α _ (- _u^2 p_t(u)). Suppose the uniform identifiable-regime condition w(ξ∣u)[∇u2lnp^t(u∣ξ)]=∇u2lnp^t(u∣ξt)∀u∈ℳ^tαE_w(ξ u)\! [ _u^2 p_t(u ξ) ]\;=\; _u^2 p_t(u _t) ∀\,u∈ M_t^α (51) holds, and let σt2:=infu∈ℳ^tαλmin(Covw(ξ∣u)(s^t(u∣ξ)))≥0 _t^2:= _u∈ M_t^α _ \! (Cov_w(ξ u)( s_t(u ξ)) )≥ 0. Then Bmarg−Bctx≥Gcσt2βtκ(ξt)κmarg≥ 0,B_marg-B_ctx\;≥\; G_c\, _t^2 _t\,κ( _t)\, _marg\;≥\;0, (52) with strict inequality whenever σt2>0 _t^2>0. Proof. Theorem 4 applied separately to the two Planners gives Bctx=ϵt+αc∂+Gcβtκ(ξt),Bmarg=ϵt+αc∂+Gcβtκmarg,B_ctx= _t+αc_∂+ G_c _tκ( _t), B_marg= _t+αc_∂+ G_c _t _marg, and the common pre-residual cancels on subtraction: Bmarg−Bctx=Gc(κ(ξt)−κmarg)βtκmargκ(ξt).B_marg-B_ctx\;=\; G_c (κ( _t)- _marg ) _t\, _marg\,κ( _t). (53) It thus suffices to show κ(ξt)−κmarg≥σt2κ( _t)- _marg≥ _t^2. Fix u∈ℳ^tαu∈ M_t^α and write H(u):=−∇u2lnp^t(u)H(u):=- _u^2 p_t(u), Σ(u):=Covw(ξ∣u)(s^t(u∣ξ))⪰0 (u):=Cov_w(ξ u)( s_t(u ξ)) 0, and κ(ξt)(u):=λmin(−∇u2lnp^t(u∣ξt))κ( _t)(u):= _ (- _u^2 p_t(u _t)) for the pointwise context-aware curvature. Proposition 3 combined with the identifiable-regime condition (51) yields H(u)=−∇u2lnp^t(u∣ξt)−Σ(u).H(u)\;=\;- _u^2 p_t(u _t)\;-\; (u). Weyl’s inequality for Hermitian matrices states λmin(A+B)≤λmin(A)+λmax(B) _ (A+B)≤ _ (A)+ _ (B). Applying it with A=−∇u2lnp^t(u∣ξt)A=- _u^2 p_t(u _t) and B=−Σ(u)B=- (u), and using λmax(−Σ(u))=−λmin(Σ(u)) _ (- (u))=- _ ( (u)) together with λmin(Σ(u))≥σt2 _ ( (u))≥ _t^2, λmin(H(u))≤κ(ξt)(u)−λmin(Σ(u))≤κ(ξt)(u)−σt2. _ (H(u))\;≤\;κ( _t)(u)- _ ( (u))\;≤\;κ( _t)(u)- _t^2. (54) Taking the infimum over u∈ℳ^tαu∈ M_t^α of both sides of (54), and using κmarg=infuλmin(H(u)) _marg= _u _ (H(u)) and κ(ξt)=infuκ(ξt)(u)κ( _t)= _uκ( _t)(u), κmarg≤κ(ξt)−σt2, _marg\;≤\;κ( _t)- _t^2, i.e., κ(ξt)−κmarg≥σt2κ( _t)- _marg≥ _t^2. Substituting into (53) yields (52); strictness is immediate from σt2>0 _t^2>0. ∎ Theorem 5 turns the informal intuition behind Problem P into a quantitative lower bound on the safety-guarantee gap: the context-aware Planner’s Theorem 4 bound improves over the context-blind one by at least Gcσt2/(βtκ(ξt)κmarg)G_c\, _t^2/( _tκ( _t) _marg), where σt2 _t^2 is the uniform infimum (over the feasible set) of the smallest posterior score-covariance eigenvalue. This gap vanishes when the posterior w collapses to a single context, and is largest when different contexts place the Planner at substantially different conditional scores, which is exactly the regime in which the contextual advantage is most visible in Figure˜12. Environmentxt+1=ft(xt,ut)x_t+1=f_t(x_t,u_t)World ModeltW_tKDEp^t,s^t p_t, s_tPPCKK grad stepsSafetyFilterxtx_tui\u_i\s^t,βt s_t, _tuppcu_ppcutu_t Simulator (Alg. 2) Planner (Alg. 1) Figure 3: Per-step information flow (Algorithms 1–2). The Simulator draws N feasibility samples, fits a KDE to produce (s^t,βt)( s_t, _t), and after receiving the Planner’s candidate uppcu_ppc, applies a safety filter (score-ascent retraction if p^t(uppc)<α p_t(u_ppc)<α) before returning utu_t to the environment. IV-F Algorithms The control law is realized through two interacting components whose information flow is depicted in Figure˜3. The information sets ℐtSimI_t^Sim and ℐtPlI_t^Pl defined in (3)–(4) make the asymmetry concrete. The Simulator’s algorithm may read xtx_t, tW_t, and ℳtM_t, whereas the Planner’s algorithm may use only s^t s_t and c. At each time step t, the Simulator (Algorithm 2) draws N feasibility samples from the oracle tW_t, updates the density model p^t p_t, and transmits the score function s^t s_t and the stiffness parameter βt _t to the Planner. The Planner (Algorithm 1) then performs K gradient descent steps on the PPC free energy using s^t s_t, and returns the candidate action uppcu_ppc to the Simulator for safety projection. If the candidate lies outside ℳ^tα M_t^α, the Simulator retracts it onto the learned manifold via score ascent before applying it to the plant. Algorithm 1 Planner: Score-Based PPC 0: Score function s^t(⋅) s_t(·) from Simulator, cost c(⋅)c(·), stiffness βt _t, step size η, inner steps K, previous action ut−1u_t-1. 0: Candidate action uppcu_ppc. 1: u←ut−1u← u_t-1 // Warm start 2: for k=1,…,Kk=1,…,K do 3: u←u−η(∇uc(u)−βts^t(u))u← u-η ( _uc(u)- _t\, s_t(u) ) 4: end for 5: return uppc←u_ppc← u Algorithm 2 Simulator: Score Oracle and Safety Filter 0: State xtx_t, feasibility oracle tW_t, density model p^t−1 p_t-1, candidate uppcu_ppc from Planner, threshold α, retraction steps J, retraction rate ηr _r, stiffness-schedule constant C>0C>0 (Equation˜58). 0: Safe action ut∈ℳ^tαu_t∈ M_t^α, updated score s^t(⋅) s_t(·), stiffness βt _t. 1: Draw N feasibility samples ui∼t(xt,⋅)\u_i\ _t(x_t,·) 2: Update density model: p^t←KDE(ui∪buffer) p_t (\u_i\ ) 3: Compute score s^t(u)=∇ulnp^t(u) s_t(u)= _u p_t(u) 4: Set βt=βcurv∗(1+C/Nt) _t=β^*_curv (1+C/ N_t ) 5: Transmit s^t(⋅) s_t(·) and βt _t to Planner 6: Receive uppcu_ppc from Planner 7: if p^t(uppc)≥α p_t(u_ppc)≥α then 8: ut←uppcu_t← u_ppc 9: else 10: ut←uppcu_t← u_ppc 11: for j=1,…,Jj=1,…,J do 12: ut←ut+ηrs^t(ut)u_t← u_t+ _r\, s_t(u_t) // Score ascent 13: if p^t(ut)≥α p_t(u_t)≥α then 14: break 15: end if 16: end for 17: end if 18: Apply utu_t to plant IV-G Computational Properties IV-G1 Compositionality For composite manifolds ℳ=⋂iℳiM= _iM_i, the total barrier decomposes additively and s^total=∑is^i s_total= _i s_i, enabling modular fusion of heterogeneous constraints. IV-G2 Scalability Each PPC gradient step requires evaluating the KDE score s^t(u) s_t(u) at a single point, which costs O(Nm)O(Nm) where N is the sample count and m is the control dimension, plus an O(m)O(m) gradient update. In contrast, CBF-QP methods [3] must form and solve a quadratic program whose per-step cost is O(m3)O(m^3) (dominated by the dense linear system solve), with additional overhead when many constraints are active. The first-order structure of PPC avoids this cubic dependence on m, and the score evaluation is embarrassingly parallel across samples. IV-G3 Online Refinement Theorem 3 gives the integrated score error rate ϵt=O(Nt−2/(m+4)) _t=O(N_t^-2/(m+4)), so the set-approximation term of the safety bound in Theorem 4 shrinks as the cumulative sample count NtN_t grows. The decreasing statistical error then allows the Planner to drive the stiffness βt _t toward its critical value βcurv∗β^*_curv, approached from above so that the strict inequality required by Proposition 2 continues to hold. Lowering βt _t relaxes the score penalty and lets the Planner exploit a larger portion of the manifold interior, trading conservative safety margins for improved task performance once the data have certified the tighter bound. Figure 4: Contextual adaptation under obstacle reshuffle (T=1000T=1000, N=300N=300 samples/step, reshuffle at t=T/2t=T/2). Top row: Trajectories colored by phase (lighter before reshuffle, darker after); red × markers indicate safety violations. Left: PPC-Context (Ours) leverages the context signal ξt _t to rapidly adapt its density model after the manifold shift, achieving safety competitive with the oracle baseline using only black-box samples. Center: Offline DRGD uses a frozen score pretrained on 500500 samples at t=0t=0; after the reshuffle its stale model fails to track the new manifold, producing frequent violations. Right: CBF-QP (oracle) has direct access to exact obstacle positions and radii at every step, which is an information advantage unavailable in the black-box setting. Bottom row: (d) Rolling safety rate with reshuffle marked (magenta line); (e) cumulative violations; (f) rolling tracking cost. PPC-Context approaches oracle-level safety using only simulator samples, while offline pretraining without online updates collapses after the shift. V Experiments We validate the online score-based framework on a 2D robot navigation task with dynamic obstacles, designed to exercise the core challenges addressed by this paper, namely a non-convex, time-varying feasibility manifold that is unknown to the central planner and must be learned online from black-box interactions. Figure 5: Contextual observation pipeline at the most constraining warmed-up step of two structurally opposite obstacle modes (top: SW cluster; bottom: NE cluster), out of the four quadrant-clustered modes the experiment cycles through every 4040 steps. Column 1: three-channel raster sts_t (red = obstacles, green = robot, blue = goal). Column 2: embedding ξt=tanh(Rvec(st)) _t= (R\,vec(s_t)); the two rows carry visibly distinct patterns, so the context kernel can separate the modes. Column 3: conditional density p^(u∣ξt) p(u _t) over the admissible action disk, with feasibility boundary ∂ℳt _t overlaid in red dashed; annotations give the conditional feasibility-mass. Column 4: context signal p^(u∣ξt)−p^(u) p(u _t)- p(u)—red adds mass, blue suppresses it. The key quantity Δsafe=Pr(u∈ℳt∣p^(⋅∣ξt))−Pr(u∈ℳt∣p^(⋅)) _safe= (u\!∈\!M_t p(· _t))- (u\!∈\!M_t p(·)) reported above each Column 4 panel measures the directional value of context—the excess feasibility-mass the conditional concentrates inside ℳtM_t relative to the pooled marginal—and is the empirical counterpart of the strictly positive σt2 _t^2 that drives the safety gap in Theorem 5. V-A Problem Setup A robot navigates a bounded workspace [0,10]2[0,10]^2 to reach a goal position while avoiding 55 moving circular obstacles. The implementation uses single-integrator dynamics qt+1=qt+utq_t+1=q_t+u_t with ‖ut‖≤umax=1.0\|u_t\|≤ u_ =1.0, which yields a non-convex feasible action set at each step while keeping the planning problem two-dimensional; the same interface extends to higher-order dynamics by letting utu_t denote the first control block passed to the Simulator. Obstacle k has radius rk∈[0.4,0.8]r_k∈[0.4,0.8] and its center ok(t)∈ℝ2o_k(t) ^2 follows a Lissajous curve: ok(t)=ck+[Aksin(ωkt+ϕk)Bkcos(νkt+ψk)],o_k(t)=c_k+ bmatrixA_k ( _kt+ _k)\\ B_k ( _kt+ _k) bmatrix, (55) where ckc_k is a base position and (Ak,Bk,ωk,νk,ϕk,ψk)(A_k,B_k, _k, _k, _k, _k) are drawn randomly per seed. In the Simulator–Planner architecture of Section˜I-B, the Simulator (physics engine) has access to the obstacle states and dynamics, while the Planner’s task cost uses the robot position and goal. For the main experiments, the Planner’s density model is trained from the same black-box feasibility oracle as our method (no explicit obstacle parameters in the optimization). The Planner minimizes the tracking cost c(qt,u)=‖qt+u−gt‖2,c(q_t,u)= \|q_t+u-g_t \|^2, (56) where gtg_t is the current goal position (updated upon arrival). The feasibility manifold is ℳt(qt)=u∈:‖qt+u−ok(t+1)‖≥rk+dsafe,∀k,M_t(q_t)= \u :\|q_t+u-o_k(t+1)\|≥ r_k+d_safe,\;∀ k \, (57) with safety margin dsafe=0.3d_safe=0.3. This set is non-convex in u and its topology changes as obstacles move. V-B Baselines We compare our online PPC controller (Algorithm 1) against five baselines spanning offline generative models, data-driven safety methods, sampling-based control, and conservative approximation. V-B1 Offline DRGD [21] A KDE score model is fit to 500500 feasibility samples from the initial obstacle configuration (t=0t=0) only. At deployment, the density is fixed and used for manifold-constrained optimization without online updates. This baseline isolates the cost of not adapting the manifold model online. V-B2 CBF-QP [2] A Control Barrier Function hk(x)=‖q−ok(t)‖2−(rk+dsafe)2h_k(x)=\|q-o_k(t)\|^2-(r_k+d_safe)^2 with oracle access to obstacle positions projects the nominal input onto the safe set via a QP (class-K coefficient γ=0.5γ=0.5). This represents the strongest structural prior. V-B3 GP-CBF [12] A Gaussian Process models the constraint function h(x,u)=mink(‖qt+1(u)−ok(t+1)‖−rk−dsafe)h(x,u)= _k(\|q_t+1(u)-o_k(t+1)\|-r_k-d_safe) online; its posterior mean serves as a learned CBF within a QP filter (RBF kernel, hyperparameters refitted every 50 steps). V-B4 Sampling-Based MPC (CEM) The Cross-Entropy Method [31] draws 300300 candidate actions from a Gaussian proposal and selects the lowest-cost feasible candidate (KDE log-probability above α), mirroring CEM-MPC with a learned model [13, 27]. Like PPC, CEM collects N=300N=300 oracle feasibility samples per step (via rejection sampling) to update its KDE density model; feasibility of CEM’s own candidate actions is then assessed against this KDE rather than by re-querying the oracle. The proposal is updated via the elite set (top 10%, 5 iterations). V-B5 Static Conservative The reachable set is conservatively inner-approximated using the maximum swept obstacle radii rkmax=rk+dsafe+maxt‖ok(t)−ok(0)‖r_k =r_k+d_safe+ _t\|o_k(t)-o_k(0)\| computed from worst-case displacement. The robot treats these inflated obstacles as static and solves a QP at each step. Because the inflation depends on the randomly drawn Lissajous amplitudes, the feasible set size varies across seeds, leading to higher cost variance than the other baselines. V-C Our Method: Online PPC The Simulator maintains a kernel density estimate over feasibility samples. At each time step, it draws N=300N=300 feasible actions via rejection sampling against the black-box oracle for (57). The score ∇ulnp^(u) _u p(u) is evaluated in closed form from the KDE. The critical stiffness βcurv∗=Gc/(κrα)β^*_curv=G_c/(κ\,r_α) is estimated online from the KDE’s barrier curvature at the α-level-set boundary, and the stiffness schedule follows βt=βcurv∗⋅(1+C/Nt), _t=β^*_curv· (1+C/ N_t ), (58) where NtN_t is the cumulative number of feasibility samples and C>0C>0 is a tunable constant. The gradient step size is set to η=min(η0, 1/(Lc+βtΛ))η= ( _0,\,1/(L_c+ _t )) following Theorem 1, where η0=0.02 _0=0.02 is a base learning rate, Λ=1/h2 =1/h^2, and h is the KDE bandwidth. Contextual observation model. For the contextual experiment (Figure˜12), the Planner receives ξt∈ℝ12 _t ^12 obtained by flattening three top-down raster channels at resolution 16×1616×16 (obstacle occupancy, robot, goal) and applying a fixed random linear projection with tanh squashing. The Simulator fits p^(u∣ξt) p(u _t) with a product-kernel KDE over pairs (u,ξ)(u,ξ) from online samples. Figure˜5 illustrates the pipeline on two structurally opposite obstacle modes (SW and NE clusters) and quantifies the safety advantage of conditioning: at the most constraining warmed-up step in each mode, the conditional concentrates 90.7%90.7\% (SW) and 64.9%64.9\% (NE) of its mass inside the feasibility manifold ℳtM_t, against only 64.2%64.2\% and 61.8%61.8\% for the context-blind marginal; the resulting Δsafe _safe values of +26.5+26.5 p and +3.1+3.1 p are the empirical counterpart of the strictly positive σt2 _t^2 that drives the safety gap in Theorem 5. V-D Evaluation Metrics 1. Safety Rate: Fraction of time steps where ut∈ℳtu_t _t. By Theorem 4, dist(ut,ℳt)dist(u_t,M_t) is bounded; when that bound is zero (i.e., sufficient samples and βt>βcurv∗ _t>β^*_curv), the action lies inside ℳtM_t and the safety rate approaches 11. 2. Normalized Cost: Total tracking cost divided by the Oracle cost. The Oracle minimizes the single-step tracking cost c(u)c(u) subject to u∈ℳtu _t with full knowledge of the obstacle geometry at each t (no look-ahead). Values near 1.01.0 indicate near-optimal per-step performance. 3. Adaptation Speed: Number of steps to recover a safety rate ≥95%≥ 95\% after a sudden obstacle configuration change at t=500t=500 (all Lissajous parameters are re-randomized). 4. Learning Convergence: Empirical integrated squared score error ‖s^−s∗‖2\| s-s^*\|^2 vs. cumulative samples NtN_t, validating the (N−2/(m+4))O(N^-2/(m+4)) rate of Theorem 3 (=(N−1/3)=O(N^-1/3) for m=2m=2). Figure 6: Main comparison across all methods (T=1000T=1000, N=300N=300 feasibility samples per step, five seeds, mean ± std); see Table˜I. The trajectory figure (Figure˜4) uses the same per-step sample budget. (a) Safety rate. (b) Normalized cost (Oracle =1.0=1.0). (c) Wall-clock time per step. V-E Experimental Protocol V-E1 Experiment 1: Main Comparison All six methods including PPC and five baselines are run for T=1000T=1000 steps with N=300N=300 black-box feasibility samples per step (Offline DRGD pretrains on 500500 samples at t=0t=0 only), across five random seeds (random obstacle trajectories and initial positions). Table˜I lists mean ± standard deviation of safety rate and normalized tracking cost; Figure˜6 shows the same statistics as bar plots. Representative trajectories appear in Figure˜4. TABLE I: Main comparison (T=1000T=1000, N=300N=300 samples/step, five seeds, mean ± std). PPC achieves the best safety among black-box methods while maintaining the lowest cost. †The high cost variance of Static Conserv. reflects the sensitivity of episode-wide worst-case inflation to random obstacle motion amplitudes; Offline DRGD’s variance is similarly driven by the mismatch between its frozen model and the current manifold. Method Safety Norm. cost PPC (Ours) .964±.009.964±.009 0.65±.180.65±.18 Offline DRGD .676±.271.676±.271 0.70±.530.70±.53 CBF-QP .905±.043.905±.043 1.34±.301.34±.30 GP-CBF .762±.025.762±.025 1.12±.141.12±.14 CEM .810±.027.810±.027 1.13±.191.13±.19 Static Conserv. .993±.014.993±.014 1.57±.861.57±.86† Figure 7: Stiffness ablation validating Proposition 2. The critical stiffness βcurv∗=Gc/(κrα)β^*_curv=G_c/(κ\,r_α) is computed from the formula (median over seeds). Safety rate (blue, left axis) rises from ∼76% \!76\% to ∼96% \!96\% and crosses the 95%95\% threshold near β=βcurv∗β=β^*_curv (red dashed), confirming the predicted phase transition. Tracking cost (orange, right axis) is low for small β but rises sharply beyond βcurv∗β^*_curv, illustrating the cost-safety trade-off governed by barrier curvature. V-E2 Experiment 2: Stiffness Ablation The environment is fixed and β is varied from 0.1βcurv∗0.1β^*_curv to 10βcurv∗10β^*_curv. We plot safety rate and normalized cost as functions of β/βcurv∗β/β^*_curv, validating the critical stiffness condition of Proposition 2: safety should emerge at β≈βcurv∗β≈β^*_curv. Results are shown in Figure˜7. V-E3 Free Energy Landscape Figure˜8 complements the stiffness sweep with a spatial analysis of the PPC free energy at a representative time step (t=200t=200, N=200N=200 samples/step). Panel (a) shows the tracking cost c(u)c(u) with the true feasibility boundary (blue solid) and the learned level set ℳ^tα M_t^α (cyan dashed). Panel (b) overlays the free-energy contours with two key points: the PPC equilibrium u∗u^* (red triangle) and the density maximizer u¯=argmaxp u= p (cyan square). The equilibrium u∗u^* is pulled toward the cost minimum but anchored near u¯ u by the barrier penalty, exactly as Proposition 2 predicts. Panel (c) reports the geometric gaps: ‖u∗−u¯‖\|u^*- u\| is well within the theoretical bound Gc/(βκ)G_c/(βκ), and dist(u∗,∂ℳt)>0dist(u^*, _t)>0 confirms that u∗u^* lies safely inside the true manifold. Figure 8: Free-energy landscape at t=200t=200. (a) Cost c(u)c(u) with true manifold boundary (blue) and learned ℳ^tα M_t^α (cyan dashed). (b) Free energy ℱ(u)F(u) with the PPC equilibrium u∗u^* (red) and density maximizer u¯ u (cyan). (c) Geometric gaps: the empirical ‖u∗−u¯‖\|u^*- u\| is well within the theoretical bound Gc/(βκ)G_c/(βκ) from Proposition 2, and dist(u∗,∂ℳt)>0dist(u^*, _t)>0 confirms the equilibrium is safely interior. Figure 9: Effect of sample budget on safety and score estimation quality. Safety rate (blue, left) crosses 95%95\% between N=100N=100 and N=500N=500. The empirical squared score error (red, right, log scale) is measured against a reference KDE built from 10,00010,000 samples. The black dashed line is a least-squares fit of aNb^aN b in log-log space; the estimated exponent b^≈−0.96 b≈-0.96 is steeper than the worst-case −1/3-1/3 from Theorem 3. See Table˜I. Figure 10: Scalability with number of obstacles (Ko∈3,5,10,15,20K_o∈\3,5,10,15,20\, three seeds, T=300T=300); see Table˜I. (a) PPC degrades gracefully (0.92→0.850.92→ 0.85) while CBF-QP drops to 0.400.40 and CEM drops to 0.620.62. (b) Total tracking cost. (c) Wall-clock time per step. Figure 11: Dynamic regret experiment validating Theorem 2 (T=500T=500, five seeds). (a) Normalized tracking cost vs. obstacle speed multiplier. (b) Cost vs. manifold path length TP_T (proxy for TV_T) with a linear fit and Pearson correlation r; the positive trend is consistent with the O(T)O(V_T) scaling in (28). (c) Safety rate remains above 90%90\% across all speeds. V-E4 Experiment 3: Sample Budget The number of manifold samples per step is varied as N∈10,50,100,500,1000N∈\10,50,100,500,1000\. To validate the O(N−2/(m+4))O(N^-2/(m+4)) rate of Theorem 3, we measure the empirical integrated squared score error by building a test KDE from an increasing number of feasibility samples and comparing its score against a high-fidelity reference KDE (10,00010,000 samples, same bandwidth). We plot safety rate and the empirical ‖s^−s∗‖2\| s-s^*\|^2 vs. N; Table˜I lists the numerical results and Figure˜9 shows the data with a log-linear power-law fit. The fitted exponent (b^≈−0.96 b≈-0.96) is steeper than the worst-case O(N−1/3)O(N^-1/3) rate of Theorem 3 for m=2m=2, consistent with the bound being conservative for this smooth problem instance. TABLE I: Sample budget ablation (T=500T=500, five seeds). Safety rate rises monotonically with N; the empirical score error ‖s^−s∗‖2\| s-s^*\|^2 (computed against a reference KDE) decays consistently with the (N−2/(m+4))O(N^-2/(m+4)) rate of Theorem 3 for m=2m=2. N Safety rate ‖s^−s∗‖2\| s-s^*\|^2 10 0.762±0.0360.762± 0.036 10.51±5.1110.51± 5.11 50 0.874±0.0260.874± 0.026 2.59±0.852.59± 0.85 100 0.945±0.0100.945± 0.010 1.10±0.161.10± 0.16 500 0.976±0.0090.976± 0.009 0.29±0.060.29± 0.06 1000 0.999±0.0010.999± 0.001 0.13±0.050.13± 0.05 Figure 12: Contextual control ablation (five seeds; T=1000T=1000, recurring mode switches every 4040 steps across nmodes=4n_modes=4 structurally distinct obstacle layouts). (a) Average total safety-violation counts: PPC-Context incurs only ≈12≈12 violations per episode versus ≈69≈69 for PPC-Marginal and ≈94≈94 for Offline Contextual, a roughly 5×5× reduction attributable to the context kernel. (b) Per-seed scatter of steady-state (within-mode) versus post-switch (first 1010 steps after each mode change) safety: PPC-Context clusters near the (1,1)(1,1) corner, PPC-Marginal is shifted toward lower safety on both axes, and Offline Contextual is dispersed further down. (c) Cumulative violations over time: PPC-Context stays near zero, while PPC-Marginal and Offline Contextual accrue violations continuously as modes cycle. (d) Rolling safety (window=40=40): PPC-Context tracks each mode switch rapidly via the context kernel, whereas PPC-Marginal dips at every switch because its buffer now mixes feasibility samples from all past modes. The visible context–marginal gap (≈ 5.6%≈\,5.6\% on aggregate safety, ≈ 4.8%≈\,4.8\% on post-switch safety) is precisely the Gcσt2/(βtκ(ξt)κmarg)G_c\, _t^2/( _t\,κ( _t)\, _marg) residual predicted by Theorem 5: multiple mode layouts drive the posterior score-covariance σt2 _t^2 strictly positive, and this directly translates into a safety advantage for PPC-Context. V-E5 Experiment 4: Scalability The number of obstacles is varied as Ko∈3,5,10,15,20K_o∈\3,5,10,15,20\. We report safety rate, total cost, and wall-clock time per step for PPC, CBF-QP, and CEM (three seeds, T=300T=300). Table˜I lists the safety rates; Figure˜10 plots all three metrics. TABLE I: Scalability: safety rate (mean ± std, three seeds, T=300T=300) as the number of obstacles KoK_o grows. PPC degrades gracefully; CBF-QP drops sharply because the QP safety filter becomes increasingly infeasible; CEM degrades as the feasible volume shrinks and its sampling proposal becomes less efficient. KoK_o PPC (Ours) CBF-QP CEM 3 0.92±0.030.92± 0.03 0.99±0.010.99± 0.01 0.88±0.010.88± 0.01 5 0.92±0.030.92± 0.03 0.89±0.070.89± 0.07 0.84±0.040.84± 0.04 10 0.89±0.020.89± 0.02 0.66±0.040.66± 0.04 0.71±0.010.71± 0.01 15 0.87±0.070.87± 0.07 0.63±0.090.63± 0.09 0.63±0.040.63± 0.04 20 0.85±0.060.85± 0.06 0.40±0.200.40± 0.20 0.62±0.050.62± 0.05 V-E6 Experiment 5: Dynamic Regret under Manifold Drift We validate the qualitative prediction of Theorem 2 by varying the obstacle speed multiplier ωmult∈0.25,0.5,1.0,2.0,4.0,8.0 _mult∈\0.25,0.5,1.0,2.0,4.0,8.0\, which controls the manifold drift rate. For each speed, PPC is run for T=500T=500 steps over five seeds. We use the normalized tracking cost as an empirical proxy for the dynamic regret. Figure˜11 shows the results. (a) Normalized cost vs. the speed multiplier ωmult _mult. (b) Cost plotted against the manifold path length T=∑t∑k‖ok(t+1)−ok(t)‖P_T= _t _k\|o_k(t+1)-o_k(t)\|, a physical proxy for TV_T in (28); a positive trend is visible, consistent with the predicted O(T)O(V_T) scaling. (c) Safety remains above 90%90\% across all speeds. V-E7 Experiment 6: Contextual Feasibility Density We design the benchmark so that the posterior conditional-score covariance σt2 _t^2 in Theorem 5 is demonstrably positive and the value of context becomes empirically observable. Four obstacles are relocated to one of nmodes=4n_modes=4 pre-generated cluster layouts (obstacles concentrated in the SW, SE, NW, or NE quadrant of the workspace), with the active mode switched to a fresh random layout every 4040 steps throughout an episode of length T=1000T=1000. Because each mode layout recurs several times per episode, the KDE buffer contains feasibility samples from multiple modes simultaneously: this is precisely the regime in which the context kernel can disambiguate past data for the currently active mode, and in which the marginal KDE must average across modes. We compare three Planners under the identical schedule, with N=25N=25 samples per step and 55 random seeds. PPC-Context updates the conditional KDE p^(u∣ξt) p(u _t) online; PPC-Marginal uses the same samples but ignores ξt _t; Offline Contextual collects 200200 samples in each of the four layouts at t=0t=0 and then freezes the conditional KDE. Table˜IV reports aggregate statistics and Figure˜12 plots the runs. Aggregate safety is 98.8%98.8\% for PPC-Context versus 93.2%93.2\% for PPC-Marginal, a 5.6%5.6\% gap that concentrates at mode boundaries: in the first ten steps after each switch PPC-Context attains 99.0%99.0\% safety, while PPC-Marginal drops to 94.2%94.2\%. Frozen Offline Contextual achieves only 90.3%90.3\% because it cannot track within-mode geometry drift once deployed. Because raw trajectory cost would reward unsafe actions that cheat through obstacles, we report the mean per-step tracking cost restricted to safe steps, normalized by the oracle’s mean per-step cost; on this fair metric all three methods lie in a narrow band (0.720.72–0.840.84) while their safety rates differ by up to 8.5%8.5\%, confirming that the online context kernel delivers its advantage as safety rather than as an inflated tracking error. These numbers quantitatively instantiate the Gcσt2/(βtκ(ξt)κmarg)G_c\, _t^2/( _t\,κ( _t)\, _marg) term in Theorem 5: the mode-switching structure makes the per-mode conditional scores substantially different, yielding a strictly positive σt2 _t^2 and a safety advantage for context-aware estimation. TABLE IV: Contextual ablation (T=1000T=1000, five seeds, recurring mode switches every 4040 steps across nmodes=4n_modes=4 obstacle layouts, N=25N=25 samples/step). “Post-switch” = mean safety in the first 1010 steps after each mode change; “Steady” = mean safety in the remaining within-mode steps. “Cost” is the mean per-step tracking cost restricted to safe time steps, normalized by the oracle’s mean per-step cost, so that constraint violations cannot artificially lower the metric by letting the robot cut through obstacles. Method Safety Post-switch Steady Cost PPC-Context .988±.004.988±.004 .990±.003.990±.003 .987±.005.987±.005 .77±.09.77±.09 PPC-Marginal .932±.014.932±.014 .942±.011.942±.011 .929±.016.929±.016 .72±.09.72±.09 Offline Contextual .903±.021.903±.021 .864±.035.864±.035 .915±.026.915±.026 .84±.12.84±.12 VI Conclusion We have developed Penalized Predictive Control for safe online control when the Planner’s access to constraint geometry is mediated by a black-box Simulator. All guarantees flow from the PPC free energy ℱ(u)=c(u)−βlnp^(u)F(u)=c(u)-β p(u): Theorem 1 reveals that its Hessian contains the Fisher information of p p, yielding a contraction rate governed by the barrier curvature κ; Theorem 2 bounds dynamic regret in terms of score variation; Theorem 3 bounds score convergence; Theorem 4 composes the preceding results into a safety bound; Proposition 3 derives, via a variational-inference argument, a mixture-Hessian identity that subtracts a posterior-score covariance from the learned curvature under marginalization, and Theorem 5 promotes this into a formal safety-guarantee gap of at least Gcσt2/(βtκ(ξt)κmarg)G_c\, _t^2/( _t\,κ( _t)\, _marg) between the Theorem 4 residuals of the context-aware and context-blind Planners whenever ξt _t is identifiable from the action, with σt2 _t^2 the minimum eigenvalue of the posterior conditional-score covariance. Experiments validate each theoretical link. The stiffness ablation (Figure˜7) confirms the phase transition at βcurv∗β^*_curv; the sample-budget sweep (Figure˜9) confirms empirical score convergence consistent with Theorem 3; with N=300N=300 samples per step, the main comparison (Table˜I) reports 96.4%96.4\% mean safety for PPC with normalized cost 0.650.65, compared with CEM at 81.0%81.0\% / 1.131.13 and GP-CBF at 76.2%76.2\% / 1.121.12 under the same black-box interface; and the contextual experiment (Figure˜12) empirically confirms Theorem 5: under a benchmark with recurring mode switches across four structurally distinct obstacle layouts, PPC-Context attains 98.8%98.8\% safety while PPC-Marginal drops to 93.2%93.2\% and a frozen Offline Contextual baseline degrades to 90.3%90.3\%, with the 5.6%5.6\% context–marginal gap concentrating at mode boundaries where the posterior score-covariance σt2 _t^2 is largest. Contextual control. This work instantiates contextual control in a concrete Simulator–Planner split where the Planner never consumes raw obstacle parameters, yet online conditional density estimation preserves safety when the environment shifts. Natural extensions include replacing the hand-crafted raster with learned encoders for image or video observations, scaling certificates to high-dimensional actions, and analyzing non-stationarity when the Simulator’s internal state evolves faster than the Planner’s sample budget. References [1] P. Absil, R. Mahony, and R. Sepulchre (2008) Optimization algorithms on matrix manifolds. Princeton University Press. Cited by: §I-2. [2] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada (2019) Control barrier functions: theory and applications. In 2019 18th European control conference (ECC), p. 3420–3431. Cited by: §I-3, §IV-E, §V-B2. [3] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada (2016) Control barrier function based quadratic programs for safety critical systems. IEEE Transactions on Automatic Control 62 (8), p. 3861–3876. Cited by: §I-3, §I-A, §IV-B, §IV-G2. [4] M. Bajelani and K. van Heusden (2026) Data-driven input-output control barrier functions. IEEE Transactions on Automatic Control. Note: Early access Cited by: §I-3. [5] A. Bemporad and M. Morari (1999) Robust model predictive control: a survey. In Robustness in identification and control, p. 207–226. Cited by: §I-1. [6] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer (2020) Data-driven model predictive control with stability and robustness guarantees. IEEE Transactions on Automatic Control 66 (4), p. 1702–1717. Cited by: §I-1. [7] L. Bold, L. Grüne, M. Schaller, and K. Worthmann (2024) Data-driven mpc with stability guarantees using extended dynamic mode decomposition. IEEE Transactions on Automatic Control. Cited by: §I-1. [8] N. Boumal (2023) An introduction to optimization on smooth manifolds. Cambridge University Press. Cited by: §I-2, §IV-B. [9] J. Bruce, M. Dennis, A. Edwards, J. Parker-Holder, Y. Shi, E. Hughes, M. Lai, A. Mavalankar, R. Steiber, C. Rae, et al. (2024) Genie: generative interactive environments. In Proceedings of the 41st International Conference on Machine Learning, p. 4583–4612. Cited by: §I. [10] J. Buerger and M. Cannon (2026) Robust adaptive NMPC using ellipsoidal tubes. Automatica. Note: Submitted Cited by: §I-1. [11] C. Chao, C. Feng, W. Sun, C. Lee, S. See, and C. Lee (2024) Maximum entropy reinforcement learning via energy-based normalizing flow. Advances in Neural Information Processing Systems 37, p. 56136–56165. Cited by: §I-2. [12] R. Cheng, G. Orosz, R. M. Murray, and J. W. Burdick (2019) End-to-end safe reinforcement learning through barrier functions for safety-critical continuous control tasks. In Proceedings of the AAAI conference on artificial intelligence, Vol. 33, p. 3387–3395. Cited by: §I-3, §V-B3. [13] K. Chua, R. Calandra, R. McAllister, and S. Levine (2018) Deep reinforcement learning in a handful of trials using probabilistic dynamics models. In Advances in Neural Information Processing Systems, Vol. 31. Cited by: §I-1, §V-B4. [14] A. Escontrela, A. Adeniji, W. Tong, B. Mazoure, and P. Abbeel (2024) Video prediction models as rewards for reinforcement learning. In Advances in Neural Information Processing Systems, Vol. 36. Cited by: §I. [15] C. E. Garcia, D. M. Prett, and M. Morari (1989) Model predictive control: theory and practice—a survey. Automatica 25 (3), p. 335–348. Cited by: §I-1. [16] Y. Guo, L. Shi, J. Chen, and C. Finn (2026) Ctrl-world: a controllable generative world model for robot manipulation. In International Conference on Learning Representations (ICLR), Cited by: §I. [17] D. Ha and J. Schmidhuber (2018) Recurrent world models facilitate policy evolution. In Advances in Neural Information Processing Systems, Vol. 31. Cited by: §I, §I. [18] K. He, S. Shi, T. van den Boom, and B. De Schutter (2026) From learning to safety: a direct data-driven framework for constrained control. IEEE Transactions on Automatic Control. Note: Early access Cited by: §I-3. [19] A. Hyvärinen (2005) Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research 6, p. 695–709. Cited by: §I-2, footnote 2. [20] H. J. Kappen (2005) Linear theory for control of nonlinear stochastic systems. Physical Review Letters 95 (20), p. 200201. Cited by: §I-1, §I, §I. [21] A. Kharitenko, Z. Shen, R. De Santi, N. He, and F. Doerfler (2025) Landing with the score: riemannian optimization through denoising. In International Conference on Learning Representations (ICLR), Cited by: §I-2, §I, §V-B1, footnote 2. [22] I. Kobyzev, S. J. Prince, and M. A. Brubaker (2020) Normalizing flows: an introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence 43 (11), p. 3964–3979. Cited by: §I-2. [23] T. Li, Y. Chen, B. Sun, A. Wierman, and S. Low (2021) Information aggregation for constrained online control. ACM SIGMETRICS Performance Evaluation Review 49 (1), p. 7–8. Cited by: §I-2, §I-3, §I, §I-C, §I-C, §I-C, §I-A. [24] T. Li, B. Sun, Y. Chen, Z. Ye, S. H. Low, and A. Wierman (2021) Learning-based predictive control via real-time aggregate flexibility. IEEE Transactions on Smart Grid 12 (6), p. 4897–4913. Cited by: §I-2, §I, §I-C, §I-C, §I-C. [25] D. Q. Mayne, M. M. Seron, and S. Raković (2005) Robust model predictive control of constrained linear systems with bounded disturbances. Automatica 41 (2), p. 219–224. Cited by: §I-1. [26] T. Pati and S. Z. Yong (2025) Robust control barrier functions for uncertain parameter-varying control affine systems with set-membership parameter estimation. IEEE Transactions on Automatic Control. Cited by: §I-3. [27] C. Pinneri, S. Sawant, S. Blaes, J. Achterhold, J. Stueckler, M. Rolinek, and G. Martius (2021) Sample-efficient cross-entropy method for real-time planning. In Conference on Robot Learning, p. 1049–1065. Cited by: §I-1, §V-B4. [28] D. Rezende and S. Mohamed (2015) Variational inference with normalizing flows. In International Conference on Machine Learning, p. 1530–1538. Cited by: §I-2. [29] N. Rhinehart, R. McAllister, K. Kitani, and S. Levine (2019) Precog: prediction conditioned on goals in visual multi-agent settings. In Proceedings of the IEEE/CVF International Conference on Computer Vision, p. 2821–2830. Cited by: §I-2, §I. [30] U. Rosolia, X. Zhang, and F. Borrelli (2018) Data-driven predictive control for autonomous systems. Annual Review of Control, Robotics, and Autonomous Systems 1 (1), p. 259–286. Cited by: §I-1. [31] R. Y. Rubinstein (1999) The cross-entropy method for combinatorial and continuous optimization. Methodology and Computing in Applied Probability 1 (2), p. 127–190. Cited by: §I-1, §V-B4. [32] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia (2020) Learning to simulate complex physics with graph networks. In International Conference on Machine Learning, p. 8459–8468. Cited by: §I. [33] A. Taylor, A. Singletary, Y. Yue, and A. Ames (2020) Learning for safety-critical control with control barrier functions. In Learning for dynamics and control, p. 708–717. Cited by: §I-3. [34] N. Tishby, F. C. Pereira, and W. Bialek (2000) The information bottleneck method. Proceedings of the 37-th Annual Allerton Conference on Communication. Cited by: §I-1. [35] E. Todorov (2006) Linearly-solvable Markov decision problems. In Advances in Neural Information Processing Systems, Vol. 19. Cited by: §I-1, §I, §I. [36] A. B. Tsybakov (2009) Introduction to nonparametric estimation. Springer Series in Statistics, Springer. Cited by: §I-C, §IV-A, §IV-D, §IV-D. [37] P. Vincent (2011) A connection between score matching and denoising autoencoders. Neural Computation 23 (7), p. 1661–1674. Cited by: §I-2. [38] K. P. Wabersich and M. N. Zeilinger (2021) A predictive safety filter for learning-based control of constrained nonlinear dynamical systems. Automatica 129, p. 109597. Cited by: §I-1.