Paper deep dive
Policy Synthesis for Finite Populations of MDP Agents under Aggregate Reach-Avoid Chance Constraints
Jie Fu, Anamika Dubey
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 10/9/2026, 7:16:22 AM
Summary
This paper addresses the synthesis of policies for finite populations of Markov Decision Process (MDP) agents subject to aggregate reach-avoid chance constraints. Unlike standard mean-field methods that enforce constraints only in expectation, this work accounts for stochastic fluctuations at finite fleet sizes by propagating the second-order moment (variance) of the empirical density alongside the mean-field trajectory using a discrete-time Lyapunov recursion. The authors apply the Cantelli inequality to convert probabilistic chance constraints into deterministic conditions on the moments, which are then incorporated into a gradient-based sequential convex approximation procedure. The method includes moment-error bounds to provide rigorous finite-N certificates and is validated on gridworld and EV-charging aggregation problems.
Entities (10)
Relation Signals (7)
Jie Fu → affiliatedwith → University of Florida
confidence 100% · Jie Fu is with the Department of Electrical and Computer Engineering, University of Florida
Anamika Dubey → affiliatedwith → Washington State University
confidence 100% · Anamika Dubey is with the School of Electrical Engineering and Computer Science, Washington State University
Cantelli Inequality → usedfor → Reach-Avoid Chance Constraints
confidence 95% · apply the Cantelli inequality to convert chance constraints into tractable deterministic conditions
Sequential Convex Approximation → usedfor → Policy Synthesis
confidence 90% · develop a model-based Sequential Convex Approximation (SCA) algorithm that updates the density-feedback policy parameters
Gridworld → usedinevaluation → Proposed Method
confidence 90% · The method is evaluated on a gridworld environment
EV-Charging Aggregation → usedinevaluation → Proposed Method
confidence 90% · evaluated on ... a power-system EV-charging aggregation problem
Mean-Field Control → limitation → Stochastic Fluctuations
confidence 85% · standard mean-field methods enforce these constraints only in expectation, which fails to account for stochastic fluctuations at finite fleet size N
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Consider a finite population of agents with decoupled Markov transition dynamics and empirical-density feedback, subject to the following constraints: with probability at least $1-\delta_r$, at least a fraction $\alpha_r$ of agents must reach a target region at some time $t^*$, while, at each time up to $t^*$, the unsafe population fraction must remain below $\beta_u$ with probability at least $1-\delta_u$. However, standard mean-field methods enforce these constraints only in expectation, which fails to account for stochastic fluctuations at finite fleet size $N$. To address this control problem, we propagate the second-order moment (variance) of the empirical density alongside the mean-field trajectory via a discrete-time Lyapunov recursion, and apply the Cantelli inequality to convert chance constraints into tractable deterministic conditions on the moments of the empirical density. We then incorporate these moment-based surrogate constraints into a gradient-based sequential convex approximation procedure for density-feedback policy synthesis. We further introduce additional moment-error bounds to construct a rigorous finite-$N$ certificate. The method is evaluated on a gridworld environment and a power-system EV-charging aggregation problem and compared with a standard deterministic population-level LP baseline.
Tags
Links
- Source: https://arxiv.org/abs/2610.12028v1
- Canonical: https://arxiv.org/abs/2610.12028v1
Trouble viewing inline? Open PDF directly →
Full Text
56,274 characters extracted from source content.
Expand or collapse full text
Policy Synthesis for Finite Populations of MDP Agents under Aggregate Reach-Avoid Chance Constraints Jie Fu Affiliation: Jie Fu is with the Department of Electrical and Computer Engineering, University of Florida, Gainesville, FL 32611, USA fujie@ufl.edu. Anamika Dubey Affiliation: Anamika Dubey is with the School of Electrical Engineering and Computer Science, Washington State University, Pullman, WA 99164, USA anamika.dubey@wsu.edu. Abstract Consider a finite population of agents with decoupled Markov transition dynamics and empirical-density feedback, subject to the following constraints: with probability at least 1−δr1- _r, at least a fraction αr _r of agents must reach a target region at some time t∗t^*, while, at each time up to t∗t^*, the unsafe population fraction must remain below βu _u with probability at least 1−δu1- _u. However, standard mean-field methods enforce these constraints only in expectation, which fails to account for stochastic fluctuations at finite fleet size N. To address this control problem, we propagate the second-order moment (variance) of the empirical density alongside the mean-field trajectory via a discrete-time Lyapunov recursion, and apply the Cantelli inequality to convert chance constraints into tractable deterministic conditions on the moments of the empirical density. We then incorporate these moment-based surrogate constraints into a gradient-based sequential convex approximation procedure for density-feedback policy synthesis. We further introduce additional moment-error bounds to construct a rigorous finite-N certificate. The method is evaluated on a gridworld environment and a power-system EV-charging aggregation problem and compared with a standard deterministic population-level LP baseline. Index Terms: chance constraints, mean-field, Markov decision process, covariance propagation. I Introduction Controlling large populations of autonomous agents subject to safety and mission requirements is a fundamental challenge in swarm robotics [4, 11, 2], multi-agent systems [3], and smart-grid applications [26]. Tracking the state of every agent individually becomes intractable as N grows because the resulting joint state space scales exponentially with N; mean-field methods instead describe the population through an aggregate distribution whose dimension is independent of N [1, 12, 13]. For finite N, however, the empirical population density exhibits O(1/N)O(1/ N) stochastic fluctuations around its mean-field limit, so a policy whose limiting trajectory satisfies a constraint may still violate that constraint with non-negligible probability. We therefore study policy synthesis for a finite-N population subject to aggregate reach–avoid chance constraints: with probability at least 1−δr1- _r, at least a fraction αr _r of the agents must occupy the target region at some time t∗t^*, while, at every time up to t∗t^*, the unsafe-region occupancy must not exceed βu _u with probability at least 1−δu1- _u. Related work. Mean-field games study equilibria among noncooperative agents, whereas mean-field control (MFC) considers a cooperative population whose aggregate objective is optimized by a centralized planner. Foundational continuous-time mean-field game theory characterizes population equilibria through coupled backward optimality and forward population-dynamics equations [12, 13]. For discrete-time mean-field MDPs, Bäuerle establishes the existence of optimal policies and convergence of finite-population control problems to their mean-field limits [7]. Related mean-field reinforcement learning (MFRL) studies develop and analyze learning algorithms for both mean-field games and cooperative mean-field control [14, 15, 16, 17]. For the class of reach-avoid chance constraints, stochastic reach–avoid methods has been developed to ensure probabilistic specifications for individual stochastic systems [6]. Yet, the method does not extend directly to enforce constraints on the empirical density of a finite multi-agent population when individual agents are stochastic systems. For MFC, safety constraints have been imposed on the deterministic population distribution or the system’s mean-field limit [1, 3, 4, 5]. Safe MFRL method [23] accounts for uncertainty in unknown transition dynamics when learning policies subject to global population constraints, but it does not account for stochastic fluctuations of the realized empirical density at finite N. Relating the proposed problem to other work focusing on stochastic variability in control design. Risk-sensitive mean-field games [20] use an exponential objective to account for variability in an individual agent’s accumulated return, rather than fluctuations in the aggregate occupancy of a finite population. Covariance steering uses the mean and covariance of a stochastic system’s state to enforce trajectory-level chance constraints [9], while distribution- and density-steering methods control state distributions under non-Gaussian disturbances or in robotic swarms [21, 22]. These approaches characterize either noise intrinsic to an individual system’s dynamics, or spatial dispersion of an continuum-population density. In contrast, the covariance we aim to constrain is the O(1/N)O(1/N) sampling covariance of the discrete empirical occupancy vector, which arises only because the population is finite and vanishes as N→∞N→∞. Related finite-population work [19] formulates mean-field-type multi-agent control as a measure-valued MDP and develops discretization-based approximations with near-optimality guarantees. This measure-valued state representation avoids the exponentially large joint state space, but the planning method discretizes the measure-valued state space into a finite set of quantized empirical distributions, whose number grows combinatorially with the population size and the number of quantization bins. Moreover, the resulting approximate policies carry no chance-constraint guarantee for the finite-N empirical density. For uncontrolled stochastic population models, [8] derives finite-size corrections to the mean-field limit at large but finite N. This is done by solving a Lyapunov equation which yields a bias-correction term of order 1/N1/N that sharpens the mean-field estimate of a performance metric’s expected value. This work does not characterize the variance of the finite-N population distribution, nor does it address chance-constrained policy synthesis. Contribution. We study a finite population of homogeneous MDP agents with population-independent transition kernels and density-feedback policies. The goal is to synthesize a policy that, when each agent follows this policy, their collective behavior optimizes an aggregate objective subject to reach–avoid chance constraints on the empirical density. We derive first-order recursions for the mean and covariance of the finite-population empirical density under density-feedback policies and use Cantelli’s inequality to obtain deterministic sufficient conditions for the reach–avoid chance constraints. Because the mean and covariance recursions are first-order approximations, we introduce bounds on the moment-approximation errors to the reach and avoid margins to derive rigorous finite-N certificates. Using these constraint surrogates, we develop a model-based Sequential Convex Approximation (SCA) algorithm that updates the density-feedback policy parameters through a sequence of second-order cone programs, using automatic differentiation to linearize the reach–avoid constraints locally. The method is validated using two experiments. I Problem Formulation I-A Swarm MDP and Empirical Density A finite-horizon MDP is M=(,,P,r,μ0)M=(S,A,P,r, _0), where S and A are finite state and action spaces, P(s′∣s,a)P(s s,a) is the transition kernel, r:×→ℝr:S×A is the reward, and μ0∈Δ() _0∈ (S) is the mean initial population density. A swarm of N homogeneous agents evolves on M. The empirical density is μtN(s)=1N∑i=1N[Xti=s],s∈, _t^N(s)= 1N _i=1^N1[X_t^i=s], s , (1) where Xti∈X_t^i is agent i’s state at time t. We allow μ0N _0^N to be random with known mean μ0=[μ0N] _0=E[ _0^N] and covariance Σ0 _0; deterministic initialization is recovered by setting Σ0=0 _0=0. Each agent selects actions via the time-varying density-feedback policy πθt(a∣s,μtN) _ _t(a s, _t^N), parameterized by θ=(θ0,…,θT−1)∈ℝT×||×||×nfeatθ=( _0,…, _T-1) ^T×|S|×|A|× n_ feat, where θt∈ℝ||×||×nfeat _t ^|S|×|A|× n_ feat is the slice of parameters used at time t and nfeatn_ feat is the number of features. Specifically, we use the softmax-linear policy πθt(a∣s,μ)=exp(θt[s,a,:]⊤ϕ(s,μ))∑a′∈exp(θt[s,a′,:]⊤ϕ(s,μ)), _ _t(a s,μ)= \! ( _t[s,a,:] φ(s,μ) ) _a \! ( _t[s,a ,:] φ(s,μ) ), (2) where ϕ(s,μ)∈ℝnfeatφ(s,μ) ^n_ feat is a twice continuously differentiable feature vector that may include state and population-density information. We write πθ _θ for brevity when the time index is clear from context. Because the transition kernel and reward are population-independent, the density-feedback policy couples the agents through the controlled transition dynamics. Remark 1 (Density-Dependent Reward: General Case). The reward r(s,a)r(s,a) in Problem 1 depends only on an agent’s state-action pair and not on the population density. More general density-dependent rewards arise in congestion models and coverage control but are outside the present scope. I-B Reach-Avoid Chance Constraints Let ,⊆T,U be disjoint target and unsafe regions, and define the free region as ℱ=∖(∪)F=S (T ). Define the aggregate density in region ℛR as μtN(ℛ)=∑s∈ℛμtN(s) _t^N(R)= _s _t^N(s). Let [T]:=1,…,T[T]:=\1,…,T\ denote the set of time steps. Problem 1 (Chance-Constrained Reach-Avoid). Given αr,βu,δr,δu∈(0,1) _r, _u, _r, _u∈(0,1), find policy parameters θ∗θ^* and a reach time t∗∈[T]t^*∈[T] solving maxθ,t∗∈[T] _θ,\,t^*∈[T] J(θ)=πθ[∑t=0T−1Rt] J(θ)=E _θ\! [ _t=0^T-1R_t ] (3a) s.t. ℙπθ(μt∗N()≥αr)≥1−δr, _θ\! ( _t^*^N(T)≥ _r )≥ 1- _r, (3b) ℙπθ(μtN()≤βu)≥1−δu,∀t≤t∗, _θ\! ( _t^N(U)≤ _u )≥ 1- _u, ∀ t≤ t^*, (3c) where Rt=∑s,ar(s,a)μtN(s)πθt(a∣s,μtN)R_t= _s,ar(s,a) _t^N(s) _ _t(a s, _t^N) is the population-average one-step reward conditional on μtN _t^N. The constraints in (3b)–(3c) are bound the reach and avoid violation probabilities separately, with the avoid guarantee imposed at each time t≤t∗t≤ t^*. I Mean-Field Dynamics and Covariance Propagation I-A Mean-Field Trajectory For a fixed finite horizon and under the regularity and initial-concentration assumptions stated below, μtN _t^N concentrates around the deterministic mean-field trajectory μ¯t μ_t as N→∞N→∞. This trajectory is defined by the mean-field flow equation [1]: μ¯t+1(s′,θ) μ_t+1(s ;θ) =∑s,aP(s′∣s,a)πθt(a∣s,μ¯t)μ¯t(s), = _s,aP(s s,a)\, _ _t(a s, μ_t)\, μ_t(s), (4) μ¯0 μ_0 =μ0. = _0. (5) We distinguish μ¯t μ_t from the true finite-N expectation mtN:=[μtN]m_t^N:=E[ _t^N]. Define the conditional-mean map gθt(μ)(s′):=∑s,aP(s′∣s,a)πθt(a∣s,μ)μ(s).g_ _t(μ)(s ):= _s,aP(s s,a)\, _ _t(a s,μ)μ(s). (6) As shown below, the finite-population process satisfies [μt+1N∣μtN]=gθt(μtN)E[ _t+1^N _t^N]=g_ _t( _t^N) exactly, whereas the deterministic trajectory satisfies μ¯t+1=gθt(μ¯t) μ_t+1=g_ _t( μ_t). By Jensen’s inequality, mt+1N=[gθt(μtN)]≠gθt(mtN)m_t+1^N=E[g_ _t( _t^N)]≠ g_ _t(m_t^N) (7) in general when density feedback makes gθtg_ _t nonlinear. So setting μ¯t=mtN μ_t=m_t^N on the right-hand side does not reproduce mt+1Nm_t+1^N on the left. Accordingly, μ¯t μ_t is a nominal deterministic population trajectory and a leading-order approximation to the finite-N mean, not the exact expectation. Define the time-t occupation measure by zt(s,a,θ):=πθt(a∣s,μ¯t(θ))μ¯t(s,θ).z_t(s,a;θ):= _ _t(a s, μ_t(θ)) μ_t(s;θ). (8) For compactness, define ft(θ,μ)(s,a):=πθt(a∣s,μ)μ(s)f_t(θ,μ)(s,a):= _ _t(a s,μ)μ(s) and t(θ)=ft(θ,μ¯t(θ))z_t(θ)=f_t(θ, μ_t(θ)). The flow matrix F∈ℝ||×||||F ^|S|×|S||A| is defined by F[s′,sa]=P(s′∣s,a)F[s ,sa]=P(s s,a). The mean-field flow can then be written compactly as μ¯t+1(θ)=Ft(θ)=Fft(θ,μ¯t(θ)),μ¯0=μ0. μ_t+1(θ)=F\,z_t(θ)=Ff_t(θ, μ_t(θ)), μ_0= _0. (9) I-B Exact Conditional Moments and First-Order Covariance Propagation The empirical density itself is a finite-state Markov process under a density-feedback policy. For an arbitrary density μ, define the effective one-agent transition matrix Qθt(μ)[s′,s]:=∑a∈πθt(a∣s,μ)P(s′∣s,a).Q_ _t(μ)[s ,s]:= _a _ _t(a s,μ)P(s s,a). (10) The conditional-mean map in (6) then satisfies gθt(μ)=Qθt(μ)μ=Fft(θ,μ)g_ _t(μ)=Q_ _t(μ)μ=Ff_t(θ,μ). Conditional on μtN=μ _t^N=μ for an attainable empirical density μ∈ΔN():=ν∈Δ():Nν(s)∈ℤ≥0∀sμ∈ _N(S):=\ν∈ (S):Nν(s) _≥ 0\ ∀ s\, there are Nμ(s)Nμ(s) agents in state s. Given this count vector, agents randomize their actions independently and their physical transitions are conditionally independent. Consequently, [μt+1N∣μtN=μ]=gθt(μ)E[ _t+1^N _t^N=μ]=g_ _t(μ) (11) exactly. Define qs(μ) q_s(μ) :=Qθt(μ)[⋅,s], :=Q_ _t(μ)[·,s], (12) Λθt(μ) _ _t(μ) :=∑s∈μ(s)[Diag(qs(μ))−qs(μ)qs(μ)⊤]. := _s μ(s) [Diag(q_s(μ))-q_s(μ)q_s(μ) ]. (13) Summing the covariances of the conditionally independent categorical next-state indicators gives the exact identity Cov(μt+1N∣μtN=μ)=1NΛθt(μ).Cov( _t+1^N _t^N=μ)= 1N _ _t(μ). (14) For density feedback, fluctuations of μtN _t^N propagate through the Jacobian of gθtg_ _t. Let Ht(θ):=∂gθt(μ)∂μ|μ¯t=FJtμ(θ),H_t(θ):= . ∂ g_ _t(μ)∂μ |_ μ_t=FJ_t^μ(θ), (15) where, coordinatewise, Jtμ(s,a)[s′]=πθt(a∣s,μ¯t)δss′+μ¯t(s)∂μ(s′)πθt(a∣s,μ)|μ¯t.J_t^μ(s,a)[s ]= _ _t(a s, μ_t) _s + μ_t(s)\\ . _μ(s ) _ _t(a s,μ) |_ μ_t. (16) When the policy is density independent, the second term vanishes and Ht(θ)=Qθt(μ¯t)H_t(θ)=Q_ _t( μ_t). We say that the initial empirical density satisfies the initial concentration condition if μ0N−μ0=Op(N−1/2). _0^N- _0=O_p(N^-1/2). (17) This condition holds for deterministic initialization and for independently sampled initial states; in the latter case, [‖μ0N−μ0‖2]=1−‖μ0‖22N.E[\| _0^N- _0\|^2]= 1-\| _0\|_2^2N. (18) Proposition 1 (First-Order Covariance Approximation). Assume a fixed horizon, the initial concentration condition, and that gθtg_ _t is twice continuously differentiable in a neighborhood of the nominal trajectory. Linearizing the conditional-mean map around the nominal trajectory gives the following first-order approximation Σ^t _t to Σt:=Cov(μtN) _t:=Cov( _t^N): Σ^t+1=HtΣ^tHt⊤+1NΛθt(μ¯t). _t+1=H_t _tH_t + 1N _ _t( μ_t). (19) The recursion is initialized with Σ^0=Σ0 _0= _0, including Σ^0=0 _0=0 for deterministic initialization. Proof sketch. The law of total covariance and (11)–(14) give Σt+1 _t+1 =1N[Λθt(μtN)]+Cov(gθt(μtN)). = 1NE[ _ _t( _t^N)]+Cov(g_ _t( _t^N)). (20) Linearizing the conditional-mean map around μ¯t μ_t yields gθt(μtN)=gθt(μ¯t)+Ht(μtN−μ¯t)+higher-order terms.g_ _t( _t^N)=g_ _t( μ_t)+H_t( _t^N- μ_t)+higher-order terms. Under this linearization, HtΣ^tHt⊤H_t _tH_t describes how existing population fluctuations propagate to the next time step. Because μtN _t^N fluctuates around μ¯t μ_t at rate Op(N−1/2)O_p(N^-1/2), the variation of Λθt(μtN) _ _t( _t^N) contributes only a smaller-order term after multiplication by 1/N1/N; we therefore evaluate Λθt _ _t at μ¯t μ_t. This gives (19). ∎ Remark 2 (Aggregate action constraints). For constraints on actions aggregated across the population, an analogous first-order variance approximation can be derived for the empirical occupation measure; see Appendix VII. IV Cantelli Chance-Constraint Certification The standard MF approach enforces only the nominal constraints μ¯t()≥αr μ_t(T)≥ _r and μ¯t()≤βu μ_t(U)≤ _u, but these mean-field conditions alone do not imply the finite-N chance constraints (3b)–(3c). Write mtN:=[μtN]m_t^N:=E[ _t^N] for the exact finite-population mean and mtN(ℛ):=cℛ⊤mtNm_t^N(R):=c_R m_t^N for the expected fraction of agents in ℛR. Here, cℛ(s)=[s∈ℛ]c_R(s)=1[s ] is the indicator vector for ℛR. We use Cantelli’s inequality to derive deterministic sufficient conditions based on the first two moments, following a construction similar to that used in chance-constrained MPC [10]. Proposition 2 (Cantelli Bounds). For any t and policy π, the reach bound below holds when mtN()≥αrm_t^N(T)≥ _r, and the avoid bound holds when mtN()≤βum_t^N(U)≤ _u: ℙπ(μtN()≥αr) ^π\! ( _t^N(T)≥ _r ) ≥1−σt2()σt2()+(mtN()−αr)2, ≥ 1- _t^2(T) _t^2(T)+(m_t^N(T)- _r)^2, (21) ℙπ(μtN()≤βu) ^π\! ( _t^N(U)≤ _u ) ≥1−σt2()σt2()+(βu−mtN())2, ≥ 1- _t^2(U) _t^2(U)+( _u-m_t^N(U))^2, (22) where σt2(ℛ)=cℛ⊤Σtcℛ _t^2(R)=c_R _tc_R denotes the exact finite-population variance. Specifically, when the relevant variance is zero, the corresponding probability bound equals one, and the Cantelli condition reduces to mtN()≥αrm_t^N(T)≥ _r for reach or mtN()≤βum_t^N(U)≤ _u for avoid. Proof. For the reach bound, let X=μtN()X= _t^N(T) and dr=mtN()−αr≥0d_r=m_t^N(T)- _r≥ 0. If σt2()>0 _t^2(T)>0, Cantelli’s inequality gives ℙ(X−mtN()≤−dr)≤σt2()σt2()+dr2.P (X-m_t^N(T)≤-d_r )≤ _t^2(T) _t^2(T)+d_r^2. Since X<αr⊆X−mtN()≤−dr\X< _r\ \X-m_t^N(T)≤-d_r\, taking complements yields (21). If the variance is zero, X is constant almost surely and the conclusion follows directly. The bound on the avoid condition is derived analogously. ∎ Lemma 1 (Sufficient Cantelli Conditions). Setting κr=(1−δr)/δr _r= (1- _r)/ _r and κu=(1−δu)/δu _u= (1- _u)/ _u, the chance constraints (3b)–(3c) are implied by: mt∗N()−αr m_t^*^N(T)- _r ≥κrσt∗(), ≥ _r\, _t^*(T), (23) βu−mtN() _u-m_t^N(U) ≥κuσt(),∀t≤t∗. ≥ _u\, _t(U), ∀ t≤ t^*. (24) Proof. For the reach constraint, let dr:=mt∗N()−αrd_r:=m_t^*^N(T)- _r. Condition (23) ensures dr≥0d_r≥ 0 and dr2≥1−δrδrσt∗2().d_r^2≥ 1- _r _r\, _t^*^2(T). Therefore, σt∗2()σt∗2()+dr2≤δr, _t^*^2(T) _t^*^2(T)+d_r^2≤ _r, and substituting this inequality into (21) gives the required probability bound. For the avoid constraint, set du:=βu−mtN()d_u:= _u-m_t^N(U). Condition (24) gives du≥0d_u≥ 0 and du2≥((1−δu)/δu)σt2()d_u^2≥((1- _u)/ _u) _t^2(U), so the same argument applied to (22) yields the avoid bound. If either variance is zero, the corresponding random variable is constant almost surely and the conclusion follows directly. ∎ Corollary 1 (Robust Finite-N Certificate). Let Σ^t _t be the covariance approximation produced by Proposition 1. Suppose validated deterministic bounds bt,ℛ,et,ℛ≥0b_t,R,e_t,R≥ 0 are available for the mean and variance approximation errors, respectively, such that |cℛ⊤(mtN−μ¯t)|≤bt,ℛ,cℛ⊤Σtcℛ≤cℛ⊤Σ^tcℛ+et,ℛ. |c_R (m_t^N- μ_t) |≤ b_t,R, c_R _tc_R≤ c_R _tc_R+e_t,R. (25) Then the exact finite-N chance constraints are implied by c⊤μ¯t∗−bt∗,−αr c_T μ_t^*-b_t^*,T- _r ≥κrc⊤Σ^t∗c+et∗,, ≥ _r c_T _t^*c_T+e_t^*,T, (26) βu−c⊤μ¯t−bt, _u-c_U μ_t-b_t,U ≥κuc⊤Σ^tc+et,,∀t≤t∗. ≥ _u c_U _tc_U+e_t,U, ∀ t≤ t^*. (27) Proof. The first enclosure lower-bounds the exact reach mean and upper-bounds the exact unsafe mean; the second upper-bounds the relevant exact variance. Thus (26)–(27) imply the conditions of Lemma 1. ∎ Replacing the moments in (23)–(24) with cℛ⊤μ¯t(θ)c_R μ_t(θ) and Σ^t(θ) _t(θ), for ℛ∈,R∈\T,U\, gives surrogate constraints. These constraints are nonlinear in θ through (9) and (19). V Chance-Constrained Mean-Field Policy Synthesis Using the mean-field and first-order covariance approximations, we formulate the following approximate synthesis problem: maxθ,t∗ _θ,t J(θ)=∑t=0T−1∑s,ar(s,a)μ¯t(s,θ)πθ(a|s,μ¯t(θ)) J(θ)= _t=0^T-1 _s,ar(s,a) μ_t(s;θ) _θ(a|s, μ_t(θ)) =∑t=0T−1∑s,ar(s,a)zt(s,a,θ) = _t=0^T-1 _s,ar(s,a)z_t(s,a;θ) s.t. (27),(26),(9),(19). eq:robust_suff_avoid, eq:robust_suff_reach, eq:mf_recurrence, eq:lyapunov. Recall that ztz_t is the occupancy measure. Because the problem is nonconvex, we solve it using the sequential convex approximation (SCA) method. At each SCA iteration, we linearize the mean-field trajectory and solve a trust-region SOCP that takes one policy-improvement step subject to the Cantelli reach-avoid constraints. We first derive the required Jacobians, then formulate the SOCP subproblem, and finally describe the trust-region update and the search over candidate reach times. V-A Linearization of the objective function Given the policy parameter θ(k)θ^(k) at iteration k, the objective function has the following first-order Taylor expansion: J(θ(k)+Δθ) J(θ^(k)+ θ) = = ∑t=0T−1∑s,ar(s,a)zt(s,a,θ(k)+Δθ) _t=0^T-1 _s,ar(s,a)z_t(s,a;θ^(k)+ θ) ≈ ≈ ∑t=0T−1∑s,ar(s,a)[zt(s,a,θ(k))+∇θzt(s,a,θ(k))Δθ] _t=0^T-1 _s,ar(s,a) [z_t(s,a;θ^(k))+ _θz_t(s,a;θ^(k)) θ ] Recall that zt(s,a,θ)=μ¯t(s,θ)πθt(a∣s,μ¯t(θ))z_t(s,a;θ)= μ_t(s;θ) _ _t(a s, μ_t(θ)), which is a function of θ through both the parameterized policy and μ¯t(θ) μ_t(θ). Let t(θ)=ft(θ,μ¯t(θ))z_t(θ)=f_t(θ, μ_t(θ)), and Dt(k):=∂μ¯t/∂θ|θ(k)∈ℝ||×nθD_t^(k):=∂ μ_t/∂θ |_θ^(k) ^|S|× n_θ denote the Jacobian of the mean-field density with respect to the policy parameters, with D0(k)=0D_0^(k)=0 because μ¯0=μ0 μ_0= _0 is fixed. Here nθ=T||||nfeatn_θ=T|S||A|n_ feat is the total number of policy parameters. Proposition 3. The Jacobian of ztz_t at iterate θ(k)θ^(k) is ∂t(θ)∂θ|θ=θ(k):=Gt(k)=Jtθ+JtμDt(k), . _t(θ)∂θ |_θ=θ^(k):=G_t^(k)=J_t^θ+J_t^μD_t^(k), (29) where Jtθ:=∂θft|(θ(k),μ¯t(k))J_t^θ:= _θf_t |_(θ^(k), μ_t^(k)) and Jtμ:=∂μft|(θ(k),μ¯t(k))J_t^μ:= _μf_t |_(θ^(k), μ_t^(k)), and Dt(k)D_t^(k) propagates forward via Dt+1(k)=FGt(k),D0(k)=0.D_t+1^(k)=F\,G_t^(k), D_0^(k)=0. (30) Consequently, the linearized occupation measure and mean trajectory are tlin(Δθ)=t(k)+Gt(k)Δθz_t lin( θ)=z_t^(k)+G_t^(k) θ and μ¯tlin(Δθ)=μ¯t(k)+Dt(k)Δθ μ_t lin( θ)= μ_t^(k)+D_t^(k) θ. Proof. By (9), μ¯t+1(θ)=Ft(θ)=Fft(θ,μ¯t(θ)) μ_t+1(θ)=Fz_t(θ)=Ff_t(θ, μ_t(θ)). Since ftf_t depends on θ both directly and through μ¯t(θ) μ_t(θ), the chain rule with ∂θft=Jtθ _θf_t=J_t^θ, ∂μft=Jtμ _μf_t=J_t^μ, dμ¯t/dθ=Dt(k)d μ_t/dθ=D_t^(k) gives exactly (29): Gt(k)=dztdθ|θ(k)=Jtθ+JtμDt(k).G_t^(k)= dz_tdθ |_θ^(k)=J_t^θ+J_t^μD_t^(k). Here, Jtθ(s,a)=μ¯t(k)(s)∇θπθt(a∣s,μ¯t(k))J_t^θ(s,a)= μ_t^(k)(s) _θ _ _t(a s, μ_t^(k)) and Jtμ(s,a)[s′]=π(k)(a∣s,μ¯t(k))δss′+μ¯t(k)(s)∂μ(s′)πθ(a∣s,μ)|μ¯t(k).J_t^μ(s,a)[s ]=π^(k)(a s, μ_t^(k))\, _s +\\ μ_t^(k)(s)\, _μ(s ) _θ(a s,μ) |_ μ_t^(k). The flow matrix F is constant, θ-independent, so differentiating μ¯t+1(θ)=Fzt(θ) μ_t+1(θ)=Fz_t(θ) directly gives Dt+1(k)=Fdzt/dθ=FGt(k)D_t+1^(k)=F\,dz_t/dθ=FG_t^(k), i.e. (30). Substituting the expression for Gt(k)G_t^(k) above into this combines the pair (29)–(30) into a single linear recursion in Dt(k)D_t^(k): Dt+1(k)=FJtμDt(k)+FJtθ,D0(k)=0.D_t+1^(k)=FJ_t^μD_t^(k)+FJ_t^θ, D_0^(k)=0. (31) The first term propagates the density sensitivity, while the second captures the direct effect of the policy parameters. Gt(k)G_t^(k) is then recovered from Dt(k)D_t^(k) via (29). For θ=θ(k)+Δθ=θ^(k)+ θ, zt(θ(k)+Δθ) z_t(θ^(k)+ θ) ≈zt(θ(k))+Gt(k)Δθ=t(k)+Gt(k)Δθ, ≈ z_t(θ^(k))+G_t^(k) θ=z_t^(k)+G_t^(k) θ, μ¯t(θ(k)+Δθ) μ_t(θ^(k)+ θ) ≈μ¯t(θ(k))+Dt(k)Δθ=μ¯t(k)+Dt(k)Δθ, ≈ μ_t(θ^(k))+D_t^(k) θ= μ_t^(k)+D_t^(k) θ, ∎ V-B SOCP Subproblem For a finite-N certificate, let bt,ℛ(k)b_t,R^(k) and et,ℛ(k)e_t,R^(k) be validated upper bounds satisfying (25) throughout the current trust region. Linearising the mean and freezing these bounds gives the k-th SCA subproblem: maxΔθ,ξ≥0 _ θ,\,ξ≥ 0\; ∑t=0T−1⊤(t(k)+Gt(k)Δθ)−λ∑tξt _t=0^T-1r (z_t^(k)+G_t^(k) θ)-λ _t _t (32a) s.t. c⊤μ¯t∗lin(Δθ)≥αr+bt∗,(k)+κrσ^t∗,(k), c_T μ_t^* lin( θ)≥ _r+b_t^*,T^(k)+ _r σ_t^*^T,(k), (32b) c⊤μ¯tlin(Δθ)−ξt≤βu−bt,(k)−κuσ^t,(k),∀t≤t∗, c_U μ_t lin( θ)- _t≤ _u-b_t,U^(k)- _u σ_t^U,(k),\;∀ t≤ t^*, (32c) ‖Δθ‖2≤ρ(k), \| θ\|_2≤ρ^(k), (32d) where ξt _t are slack variables, r is the reward function written in the vector form, t∗t^* is fixed (see below), ρ(k)>0ρ^(k)>0 is the trust-region radius, μ¯tlin(Δθ):=μ¯t(k)+Dt(k)Δθ μ_t lin( θ):= μ_t^(k)+D_t^(k) θ (33) is the linearized mean trajectory of Proposition 3, and σ^tℛ,(k):=max(cℛ⊤Σ^t(k)cℛ+et,ℛ(k), 0),ℛ∈,, σ_t^R,(k):= (c_R _t^(k)c_R+e_t,R^(k),\,0 ), ∈\T,U\, (34) is a certified upper bound on the Cantelli standard deviation of region ℛR. Here Σ^t(k) _t^(k) is the first-order covariance recursion evaluated at the current iterate. The resulting constraints are the robust conditions (26)–(27) with a linearized mean. Problem (32) is an SOCP (the objective and constraints are affine in Δθ θ; the trust region is a second-order cone constraint), solvable in polynomial time. Fixed-covariance surrogate The mean trajectory is linearized as μ¯tlin(Δθ) μ_t lin( θ), whereas the covariance surrogate is held at Σ^t(k) _t^(k). This is because Σ^t(θ) _t(θ) depends on θ through the nonlinear Lyapunov recursion (19), which contains the term Ht(θ)Σ^t(θ)Ht(θ)⊤H_t(θ) _t(θ)H_t(θ) . Freezing σ^tℛ,(k) σ_t^R,(k) keeps the constraints affine in Δθ θ and preserves the SOCP form of the subproblem. The trust region keeps this approximation local, and the covariance is recomputed after every accepted step. By Proposition 1, the change in covariance over a gradient step is typically O(‖Δθ‖/N)O(\| θ\|/N). Such a small error can be negligible or absorbed into the error et,ℛe_t,R of Corollary 1. Trust-Region Update For each SCA run, t∗t^* is fixed to one candidate reach time. Since a step that improves J can violate the robust constraints after linearization, we accept or reject it using an L1L_1-penalty merit function evaluated at the nonlinear robust margins, M(θ):=J(θ)−λM[ M(θ):=J(θ)- _M [ max(0,−greach(θ)) (0,-g_ reach(θ) ) (35) +max(0,−gavoid(θ))], + (0,-g_ avoid(θ) ) ], where we define the robust standard-deviation margin qt,ℛ(θ):=cℛ⊤Σ^t(θ)cℛ+et,ℛ(θ)q_t,R(θ):= c_R _t(θ)c_R+e_t,R(θ). Then greach(θ) g_ reach(θ) =c⊤μ¯t∗(θ)−αr−bt∗,(θ)−κrqt∗,(θ), =c_T μ_t^*(θ)- _r-b_t^*,T(θ)- _rq_t^*,T(θ), gavoid(θ) g_ avoid(θ) =mint≤t∗[βu−c⊤μ¯t(θ)−bt,(θ) = _t≤ t^* [ _u-c_U μ_t(θ)-b_t,U(θ) −κuqt,(θ)]. - _uq_t,U(θ) ]. The bounds bt,ℛ(θ)b_t,R(θ) and et,ℛ(θ)e_t,R(θ) are recomputed or validated at the candidate policy; λM>0 _M>0 penalizes constraint violation. If the SOCP is feasible, set θ′=θ(k)+Δθ(k)θ =θ^(k)+ θ^(k) and evaluate the candidate using the nonlinear merit function. The trust-region radius is then updated as: ρ(k+1)=min(βincρ(k),ρmax)if M(θ′)≥M(θ(k))−ε,βdecρ(k)otherwise.ρ^(k+1)= cases ( _ incρ^(k), _ max)&if M(θ )≥ M(θ^(k))- ,\\ _ decρ^(k)&otherwise. cases (36) where βdec∈(0,1) _ dec∈(0,1), βinc>1 _ inc>1, and ε≥0 ≥ 0 is the acceptance tolerance. Candidates that fail the merit-function acceptance test are rejected and the trust-region radius is reduced. If instead the SOCP (32) is infeasible, the trust region may be too small to reach the required constraint shifts. We therefore increase ρ(k)ρ^(k) up to ρmax _ max. If infeasibility persists at ρmax _ max, we terminate the current SCA run and proceed to the next candidate t∗t^*. Choose candidate reach time The existential reach constraint requires only that (32b) holds at some t∗t^*, not a specific terminal time T. Given an initial warm-start policy, we enumerate the top-Kt∗K_t^* candidate time steps by their initial reach margins c⊤μ¯t(0)−αr−bt,(0)−κrσ^t,(0)c_T μ_t^(0)- _r-b_t,T^(0)- _r σ_t^T,(0) (descending). We run SCA for each candidate in turn and stop when both margins in the merit function (35) satisfy rm:=greach≥0r_m:=g_ reach≥ 0 and am:=gavoid≥0a_m:=g_ avoid≥ 0. If no candidate satisfies both margins, we return the infeasible result with the largest M(θ)M(θ) across all Kt∗K_t^* candidates. Because the search stops at the first candidate satisfying both margins, rather than comparing all Kt∗K_t^* candidates, the returned policy may be suboptimal in reward. Algorithm 1 SCA for Reach-Avoid Swarm Control Input: θ(0)θ^(0), N, K, Kt∗K_t^*, ρ0 _0, βinc _ inc, βdec _ dec, ρmax _ max, λ, λM _M, ε Output: θ∗θ^*, μ¯∗ μ^*, Σ^∗ ^*, rm∗r_m^*, am∗a_m^* 1 Compute (μ¯t(0),Σ^t(0))( μ_t^(0), _t^(0)) via forward sim. and Lyapunov recursion at θ(0)θ^(0); 2 Compute or validate bt,ℛ(0),et,ℛ(0)b_t,R^(0),e_t,R^(0) for the initial trust-region radius ρ0 _0; 3 Rank t∈[T]t∈[T] by c⊤μ¯t(0)−αr−bt,(0)−κrσ^t,(0)c_T μ_t^(0)- _r-b_t,T^(0)- _r σ_t^T,(0) descending; set ∗←T^*← top-Kt∗K_t^* candidates; 4 foreach t∗∈∗t^* ^* do 5 θ←θ(0)θ←θ^(0); 6 ρ←ρ0ρ← _0; 7 for k=0,…,K−1k=0,…,K-1 do 8 Compute (μ¯t(k),Σ^t(k))( μ_t^(k), _t^(k)); 9 Compute or validate bt,ℛ(k),et,ℛ(k)b_t,R^(k),e_t,R^(k) for the current trust region; 10 Compute Gt(k),Dt(k)G_t^(k),D_t^(k) via Prop. 3; 11 Solve SOCP (32) →Δθ(k)→ θ^(k); 12 if infeasible then 13 if ρ=ρmaxρ= _ max then 14 break; // proceed to the next t∗t^* 15 ρ←min(βincρ,ρmax)ρ← ( _ incρ, _ max); continue; 16 θ′←θ+Δθ(k)θ ←θ+ θ^(k); 17 if M(θ′)≥M(θ)−εM(θ )≥ M(θ)- then 18 θ←θ′θ←θ ; 19 ρ←min(βincρ,ρmax)ρ← ( _ incρ, _ max); 20 else 21 ρ←βdecρ← _ decρ; 22 Evaluate rm,am,M(θ)r_m,a_m,M(θ); 23 if rm≥0r_m≥ 0 and am≥0a_m≥ 0 then 24 Store the candidate; 25 break; // first candidate satisfying both margins 26 Store if M(θ)M(θ) is the largest infeasible value so far; 27 return stored (θ∗,μ¯∗,Σ^∗,rm∗,am∗)(θ^*, μ^*, ^*,r_m^*,a_m^*) VI Experimental Validation We validate the algorithm using two case studies.11 1 Code is available at https://github.com/jiefu2017/mfmdp-swarm-control. As a baseline, we use a mean-field LP that solves the mean-field relaxation of Problem 1 directly over the occupation measure t(s,a)≥0z_t(s,a)≥ 0, rather than over policy parameters θ: max0,…,T−1≥0 _z_0,…,z_T-1≥ 0\; ∑t=0T−1⊤t _t=0^T-1r z_t (37a) s.t. μ¯t+1=Ft,μ¯0=μ0, μ_t+1=Fz_t, μ_0= _0, (37b) c⊤μ¯t∗≥αr,c⊤μ¯t≤βu∀t≤t∗, c_T μ_t^*≥ _r, c_U μ_t≤ _u\;\;∀ t≤ t^*, (37c) enforcing only the mean constraints. Since t∗t^* is existential, (37) is solved once for each candidate t∗∈[T]t^*∈[T] and the best feasible reward is kept. Because the objective and constraints are linear in tz_t, (37) is an LP and is solved exactly. The policy is recovered as πt(a∣s)=zt∗(s,a)/μ¯t∗(s) _t(a s)=z_t^*(s,a)/ μ_t^*(s) (uniform where μ¯t∗(s)=0 μ_t^*(s)=0). VI-A Gridworld Environment Our first case study is a 6×66× 6 discrete gridworld (Fig. 1, ||=36|S|=36, ||=5|A|=5: up/down/left/right/stay, T=20T=20). Transitions: intended action with probability 0.90.9; otherwise a uniformly random other action executes instead (probability 0.1/40.1/4 each). Moves into a wall cell, or off the grid, leave the agent in place. Start: (5,0)(5,0). A wall at x=2x=2 separates the swarm from three reach targets =(0,0),(0,3),(0,4)T=\(0,0),(0,3),(0,4)\ (reward 0.50.5 each) and one reward-only target (0,5)(0,5) (reward 1.01.0, not in T): higher reward, but not reach-designated, mirroring the reward/safety tension of the reach-avoid setting considered here. The wall has three physically impassable cells (2,0),(2,2),(2,4)(2,0),(2,2),(2,4) and three openings (2,1),(2,3),(2,5)(2,1),(2,3),(2,5); every route from the start to any target must pass through one of these openings. Fig. 1: The 6×66× 6 gridworld environment. Start (red, S) at (5,0)(5,0); wall at x=2x=2 with three impassable cells (black) and three open gap cells s1,s2,s3s_1,s_2,s_3 (orange) that every route must cross; three reach targets (green, T) and one reward-only target (blue, R), both past the wall. Rather than one avoid constraint, we consider three independent per-cell instances of the avoid clause in (37) (and its Cantelli counterpart (32c)), tracking each cell’s instantaneous occupancy (congestion): for opening i, ℙ(μtN(si)>0.30)≤0.40P( _t^N(s_i)>0.30)≤ 0.40, i.e. βu,i=0.30 _u,i=0.30, δu,i=0.40 _u,i=0.40 (κu,i=(1−δu,i)/δu,i=1.5≈1.22 _u,i= (1- _u,i)/ _u,i= 1.5≈ 1.22), which means the probability that more than 30%30\% of agents occupy the same opening simultaneously is at most 40%40\%. Reach parameters: αr=0.3 _r=0.3, δr=0.1 _r=0.1 (κr=3 _r=3), which means with probability at least 90%90\%, at least 30%30\% of the agents reach the target. Table I reports the surrogate Cantelli margins rm,amr_m,a_m at the t∗t^* selected by each method, together with Monte Carlo estimates over 500 trials (fleet sizes N∈10,20,50,100N∈\10,20,50,100\), using Algorithm 1 with the LP-based warm start. P^v P_ v reports the maximum over the three per-cell violation rates. The target values are P^reach≥0.90 P_ reach≥ 0.90 and P^v≤0.40 P_ v≤ 0.40. TABLE I: Gridworld results. rmr_m, ama_m: surrogate Cantelli margins (≥0≥ 0 indicates that the surrogate condition is satisfied). J J: empirical mean total reward. P^r P_ r: empirical reach probability. P^v P_ v: maximum empirical per-cell avoid violation rate. Emp. met? indicates whether both empirical reach and avoid requirements are met. MF-LP (baseline) SCA (proposed) N J J P^r P_ r P^v P_ v emp. met? J J P^r P_ r P^v P_ v rmr_m/ama_m 10 9.39 0.91✓ 0.38✓ ✓ 9.38 1.00✓ 0.40✓ −.108×/−.194×-.108^×/-.194^× 20 9.38 0.95✓ 0.39✓ ✓ 8.97 1.00✓ 0.12✓ +.029/+.000+.029/+.000 50 9.37 1.00✓ 0.44× × 9.23 1.00✓ 0.08✓ +.000/+.000+.000/+.000 100 9.38 1.00✓ 0.47× × 9.28 1.00✓ 0.12✓ +.090/+.000+.090/+.000 Key findings. (1) SCA achieves positive Cantelli margins and stays under budget (P^v≤0.12 P_ v≤ 0.12) for N≥20N≥ 20. At N=10N=10, the empirical reach rate is 100%100\% and the empirical per-cell violation rate is also 40%40\%. However, rm=−0.108r_m=-0.108 and am=−0.194a_m=-0.194 mean that the surrogate Cantelli conditions are not satisfied. This reflects the conservatism of the distribution-free Cantelli bound: the empirical trials meet both requirements even though the surrogate test does not. In comparison, the MF-LP violation rate does not decrease with N, rising from 0.380.38 at N=10N=10 to 0.470.47 at N=100N=100. This is because the MF-LP’s mean trajectory sits exactly at the boundary, μ¯t(si)=βu,i=0.30 μ_t(s_i)= _u,i=0.30. By the central limit theorem, the occupancy μtN(si) _t^N(s_i) is approximately Gaussian for large N, and because a Gaussian is symmetric about its mean, ℙ(μtN(si)>βu,i)≈50%P( _t^N(s_i)> _u,i)≈ 50\%. Lastly, for N≥20N≥ 20, the SCA reward is within 5%5\% of the MF-LP reward: the three reach targets give SCA freedom to redistribute reward-collection time without sacrificing much of the safety margin. VI-B EV-Charging Aggregator Our second case study is an EV-charging application in which the avoid constraint is a transformer capacity limit on the aggregate action, i.e. on the occupation measure tNz_t^N itself. Environment. States =0,…,10S=\0,…,10\ are state-of-charge (SOC) bins in 10%10\% increments; the target set is =8,9,10T=\8,9,10\ (SOC ≥80%≥ 80\%). Three actions: fast charging (reward 33; advances SOC by 22 bins w.p. 0.60.6, else 11 bin), slow charging (reward 11; advances by 11 bin w.p. 0.90.9, else unchanged), and idle/no-charging (reward 00; SOC unchanged deterministically). SOC 1010 is absorbing. All agents start at SOC 00. The policy feature vector is ϕ(s,μ)=[1,s/10,μ()]φ(s,μ)=[1,\,s/10,\,μ(T)]: a bias, normalized SOC, and a global density-feedback feature reporting how much of the fleet has already reached the target. Rather than bounding how many agents occupy an unsafe state, the operational constraint here is that, with probability at least 1−δu1- _u, no more than a fraction βu _u of the fleet may fast-charge simultaneously, at any time. Define the empirical occupation measure as ztN(s,a):=N−1∑i=1N[Xti=s,Ati=a]z_t^N(s,a):=N^-1 _i=1^N1[X_t^i=s,A_t^i=a], the fraction of agents in state s taking action a at time t. With cF(s,a)=[a=fast]c_F(s,a)=1[a= fast], ℙπ(cF⊤tN≤βu)≥1−δu,∀t.P^π\! (c_F z_t^N≤ _u )≥ 1- _u, ∀ t. (38) For the aggregate-action constraint (38), we use the occupation-measure variance approximation in Appendix VII, with c=cFc=c_F. The surrogate Cantelli condition and SCA subproblem then use Gt(k)G_t^(k), the Jacobian of tz_t, for this constrained quantity. The following parameters are used: T=12T=12, αr=0.8 _r=0.8, δr=0.1 _r=0.1 (κr=3 _r=3), βu=0.5 _u=0.5, δu=0.3 _u=0.3 (κu=1.53 _u=1.53). Unlike in the gridworld case study, the reach time is fixed at t∗=T=12t^*=T=12 because the requirement is that at least 80%80\% of the fleet attain an SOC of at least 80%80\% by the terminal time. The avoid/congestion constraint is then active over the entire charging window. TABLE I: EV-charging results. rmr_m, ama_m: surrogate Cantelli margins (≥0≥ 0 indicates that the surrogate condition is satisfied), with ama_m computed using (39). J J: empirical mean total reward. P^r P_ r: empirical reach probability. P^v P_ v: maximum empirical per-timestep marginal avoid violation rate. Emp. met? indicates whether both empirical requirements (P^r≥0.90 P_ r≥ 0.90 and P^v≤0.30 P_ v≤ 0.30) are met. MF-LP (baseline) SCA (proposed) N J J P^r P_ r P^v P_ v emp. met? J J P^r P_ r P^v P_ v rmr_m/ama_m 20 19.11 1.00✓ 0.45× × 17.62 1.00✓ 0.27✓ +.138✓/−.134×+.138 /-.134^× 50 19.09 1.00✓ 0.49× × 17.56 1.00✓ 0.06✓ +.196✓/+.000✓+.196 /+.000 100 19.09 1.00✓ 0.52× × 17.99 1.00✓ 0.08✓ +.199✓/+.000✓+.199 /+.000 200 19.09 1.00✓ 0.51× × 18.31 1.00✓ 0.07✓ +.200✓/+.000✓+.200 /+.000 Key findings. (1) The MF-LP maximum per-time-step violation rate remains near the CLT-predicted 50%50\% level (0.450.45–0.520.52) across all tested fleet sizes, which is a property of the LP’s mean-only formulation. (2) SCA uses the occupation-measure variance approximation in (39). It satisfies the surrogate conditions for N≥50N≥ 50 (am≈0a_m≈ 0, P^v≤0.08 P_ v≤ 0.08 at N=50,100,200N=50,100,200). As N grows, the aggregate charging fraction concentrates around its mean, reducing the uncertainty margin required by the surrogate constraint. For N≥50N≥ 50, the observed violation rate remains at or below 0.080.08. Figure 2 illustrates the corresponding mean charging profiles: the MF-LP fast-charging fraction sits at the βu=0.5 _u=0.5 limit, whereas the SCA fraction remains around 0.420.42. We also observe that SCA’s slow-charging fraction is higher than MF-LP’s, especially at earlier timesteps: because SCA’s average fast-charging fraction is lower, it must rely more on slow-charging throughout to still satisfy the reach constraint. Lastly, as an ablation study, we tested replacing the scalar feature μ()μ(T) with the full density vector μt∈ℝ|| _t ^|S| (raising nfeatn_ feat from 33 to 2+||=132+|S|=13, and θ from 1,1881,188 to 5,1485,148 parameters). The two feature designs produced similar results, so we omit the full-density results. Computational performance. On an Apple M1 Pro with 16 GB RAM, the gridworld policy has 18,000 parameters and uses at most 25 SCA iterations for each of up to three candidate reach times. Including SCA and the 500-trial evaluations of both methods, the runtime per fleet size was 53–174 s. The EV policy has 1,188 parameters and uses one fixed reach time; the corresponding runtime was 12–15 s. The experiments used Python 3.14, JAX 0.11, and CVXPY 1.9. Fig. 2: Aggregate fast-charging, slow-charging, and idle fractions over time (summing to 1 at every t), MF-LP’s optimal solution (top) versus the SCA solution at N=100N=100 (bottom). MF-LP saturates the βu=0.5 _u=0.5 cap at every timestep and avoids slow-charging until the departure deadline forces it (t=9t=9–1111); SCA holds a margin below the cap throughout. VII Conclusion We presented a policy synthesis framework for chance-constrained reach-avoid control of finite-N swarms on discrete MDPs. The proposed method optimizes a density-feedback policy via Sequential Convex Approximation over a sequence of SOCPs for constrained policy gradient updates. In both case studies, the MF-LP baseline’s mean-only solution lies on the active constraint boundary, causing its violation probability to approach approximately 50%50\% rather than vanish as N increases. The standard-deviation penalties in the SCA constraints decrease at rate O(1/N)O(1/ N), reducing conservatism for larger populations. Future work includes extending this approach to true mean-field-type MDPs, where the transition kernel and reward function are coupled to the population’s empirical density, and to chance-constrained finite-N MDPs with common noise. [Variance of Aggregate Occupation Measures] For a linear statistic Yt=c⊤tNY_t=c z_t^N of the empirical occupation measure, where c∈ℝ||||c ^|S||A|, define hc(μ) h_c(μ) :=∑s,ac(s,a)πθt(a∣s,μ)μ(s), := _s,ac(s,a) _ _t(a s,μ)μ(s), vt(s,c) v_t(s,c) :=Vara∼πθt(⋅∣s,μ¯t)[c(s,a)]. := Var_a _ _t(· s, μ_t)[c(s,a)]. Conditioned on μtN _t^N, agents at each state choose actions independently, so their action counts have multinomial distributions. Applying the law of total variance and linearizing the conditional mean hch_c around μ¯t μ_t gives Var[Yt]≈1N∑sμ¯t(s)vt(s,c)+∇hc(μ¯t)⊤Σ^t∇hc(μ¯t), Var[Y_t]≈ 1N _s μ_t(s)v_t(s,c)+∇ h_c( μ_t) _t∇ h_c( μ_t), (39) where Σ^t _t is the state-density covariance approximation from Proposition 1. The first term is the conditional action-sampling variance; the second propagates state-density fluctuations through the density-dependent conditional mean. Define σ^t(c) σ_t(c) as the square root of the right-hand side of (39). For the EV-charging example, c=cFc=c_F selects the fraction of agents taking the fast-charging action. References [1] B. Açıkmeşe and D. S. Bayard, “A Markov chain approach to probabilistic swarm guidance,” in Proc. ACC, 2012, p. 6300–6307. [2] N. Demir, U. Eren, and B. Açıkmeşe, “Decentralized probabilistic density control of autonomous swarms with safety constraints,” Autonomous Robots, vol. 39, no. 4, p. 537–554, 2015. [3] M. El Chamie, Y. Yu, and B. Açıkmeşe, “Convex synthesis of randomized policies for controlled Markov chains with density safety upper bound constraints,” in Proc. ACC, 2016, p. 6290–6295. [4] X. Yang, Y. Hu, H. Gao, K. Ding, Z. Li, P. Zhu, Y. Sun, and C. Liu, “Risk-aware non-myopic motion planner for large-scale robotic swarm using CVaR constraints,” in Proc. IROS, 2024, p. 5784–5790. [5] S. Wu Fung and L. Nurbekyan, “Mean-field control barrier functions: A framework for real-time swarm control,” in Proc. ACC, 2025. [6] Đ. Žikelić, M. Lechner, T. A. Henzinger, and K. Chatterjee, “Learning control policies for stochastic systems with reach-avoid guarantees,” in Proc. AAAI, 2023, p. 11926–11935. [7] N. Bäuerle, “Mean field Markov decision processes,” Appl. Math. Optim., vol. 88, art. 12, 2023. [8] N. Gast and B. Van Houdt, “A refined mean field approximation,” Proc. ACM Meas. Anal. Comput. Syst., vol. 1, no. 2, art. 33, 2017. [9] K. Okamoto, M. Goldshtein, and P. Tsiotras, “Optimal covariance control for stochastic systems under chance constraints,” IEEE Control Syst. Lett., vol. 2, no. 2, p. 266–271, 2018. [10] M. Farina and R. Scattolini, “Model predictive control of linear systems with multiplicative unbounded uncertainty and chance constraints,” Automatica, vol. 70, p. 258–265, 2016. [11] B. Doerr and R. Linares, “Control of large swarms via random finite set theory,” in Proc. ACC, 2018, p. 2904–2909. [12] J.-M. Lasry and P.-L. Lions, “Mean field games,” Japanese J. Math., vol. 2, p. 229–260, 2007. [13] M. Huang, R. P. Malhamé, and P. E. Caines, “Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle,” Commun. Inf. Syst., vol. 6, no. 3, p. 221–252, 2006. [14] Y. Yang, R. Luo, M. Li, M. Zhou, W. Zhang, and J. Wang, “Mean field multi-agent reinforcement learning,” in Proc. ICML, 2018, p. 5571–5580. [15] J.-P. Fouque, M. Laurière, and M. Zhang, “Convergence of actor-critic learning for mean field games and mean field control in continuous spaces,” arXiv:2511.06812, 2025. [16] A. Angiuli, J.-P. Fouque, and M. Laurière, “Unified reinforcement Q-learning for mean field game and control problems,” Math. Control Signals Syst., vol. 34, p. 217–271, 2022. [17] R. Carmona, M. Laurière, and Z. Tan, “Linear-quadratic mean-field reinforcement learning: convergence of policy gradient methods,” arXiv:1910.04295, 2019. [18] K. Cui, C. Fabian, A. Tahir, and H. Koeppl, “Major-minor mean field multi-agent reinforcement learning,” in Proc. 41st Int. Conf. Machine Learning, vol. 235, 2024, p. 9603–9632. [19] E. Bayraktar, N. Bäuerle, and A. D. Kara, “Finite approximations for mean field type multi-agent control and their near optimality,” Appl. Math. Optim., vol. 92, art. 7, 2025. [20] N. Saldi, T. Başar, and M. Raginsky, “Approximate Markov-Nash equilibria for discrete-time risk-sensitive mean-field games,” Math. Oper. Res., vol. 45, no. 4, p. 1596–1620, 2020. [21] V. Sivaramakrishnan, J. Pilipovsky, M. Oishi, and P. Tsiotras, “Distribution steering for discrete-time linear systems with general disturbances using characteristic functions,” in Proc. ACC, 2022, p. 4183–4190. [22] C. Sinigaglia, A. Manzoni, F. Braghin, and S. Berman, “Robust optimal density control of robotic swarms,” Automatica, vol. 176, art. 112218, 2025. [23] M. Jusup, B. Pásztor, T. Janik, K. Zhang, F. Corman, A. Krause, and I. Bogunovic, “Safe model-based multi-agent mean-field reinforcement learning,” in Proc. 23rd Int. Conf. Autonomous Agents and Multiagent Systems (AAMAS), 2024, p. 973–982. [24] P. Mestres, B. Werner, R. K. Cosner, and A. D. Ames, “Probabilistic control barrier functions: safety in probability for discrete-time stochastic systems,” in Proc. ACC, 2026, p. 110–117. [25] G. Rapakoulias, A. R. Pedram, and P. Tsiotras, “Steering large agent populations using mean-field Schrödinger bridges with Gaussian mixture models,” IEEE Control Syst. Lett., vol. 9, p. 1760–1765, 2025. [26] R. Couillet, S. M. Perlaza, H. Tembine, and M. Debbah, “A mean field game analysis of electric vehicles in the smart grid,” in Proc. IEEE INFOCOM Workshops, 2012, p. 79–84.