Paper deep dive
Safe Continuous-time Multi-Agent Reinforcement Learning via Epigraph Form
Xuefeng Wang, Lei Zhang, Henglin Pu, Husheng Li, Ahmed H. Qureshi
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 7/21/2026, 12:32:13 AM
Summary
This paper proposes a novel framework for Safe Continuous-time Multi-Agent Reinforcement Learning (CT-MARL) using an epigraph-based reformulation. It addresses the limitations of discrete-time MDPs and existing continuous-time methods that struggle with safety constraints due to value discontinuities. The authors introduce a Continuous-Time Constrained MDP (CT-CMDP) formulation and transform it into a continuous epigraph form using an auxiliary state variable. This allows for the use of Physics-Informed Neural Networks (PINNs) within an actor-critic architecture to solve the resulting Hamilton-Jacobi-Bellman (HJB) partial differential equations. The method integrates inner and outer optimization to compute optimal auxiliary states, ensuring stable training and improved performance on continuous-time safe multi-particle environments and MuJoCo benchmarks.
Entities (8)
Relation Signals (7)
Epigraph Reformulation → transforms → CT-CMDP
confidence 96% · transforms discrete MDPs into CT-CMDPs via an epigraph-based reformulation.
CT-MARL → addresses → discrete-time limitations
confidence 95% · most algorithms still rely on a discrete-time Markov Decision Process (MDP)... leading to degraded performance and motivating the development of continuous-time MARL
Epigraph Reformulation → enables → PINN-based learning
confidence 94% · transforms discontinuous constrained values into a continuous form suitable for PDE-based learning... proposing a novel physics-informed neural network (PINN)-based actor–critic method
PINN → solves → HJB PDE
confidence 93% · PINNs have emerged as a common approach to approximate HJB solutions
Auxiliary State z → usedin → Epigraph Reformulation
confidence 92% · We introduce an auxiliary state variable z(t) to reformulate Eq. 2 using the epigraph forms.
Proposed Method → evaluatedon → MuJoCo
confidence 90% · and safe multi-agent MuJoCo benchmarks.
Proposed Method → evaluatedon → MPE
confidence 90% · We evaluate our approach on continuous-time safe multi-particle environments (MPE)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Multi-agent reinforcement learning (MARL) has made significant progress in recent years, but most algorithms still rely on a discrete-time Markov Decision Process (MDP) with fixed decision intervals. This formulation is often ill-suited for complex multi-agent dynamics, particularly in high-frequency or irregular time-interval settings, leading to degraded performance and motivating the development of continuous-time MARL (CT-MARL). Existing CT-MARL methods are mainly built on Hamilton-Jacobi-Bellman (HJB) equations. However, they rarely account for safety constraints such as collision penalties, since these introduce discontinuities that make HJB-based learning difficult. To address this challenge, we propose a continuous-time constrained MDP (CT-CMDP) formulation and a novel MARL framework that transforms discrete MDPs into CT-CMDPs via an epigraph-based reformulation. We then solve this by proposing a novel physics-informed neural network (PINN)-based actor-critic method that enables stable and efficient optimization in continuous time. We evaluate our approach on continuous-time safe multi-particle environments (MPE) and safe multi-agent MuJoCo benchmarks. Results demonstrate smoother value approximations, more stable training, and improved performance over safe MARL baselines, validating the effectiveness and robustness of our method.
Tags
Links
- Source: https://arxiv.org/abs/2602.17078v1
- Canonical: https://arxiv.org/abs/2602.17078v1
Trouble viewing inline? Open PDF directly →
Full Text
93,246 characters extracted from source content.
Expand or collapse full text
Safe Continuous-time Multi-Agent Reinforcement Learning via Epigraph Form Xuefeng Wang1, Lei Zhang111footnotemark: 1, Henglin Pu1, Husheng Li1, Ahmed H. Qureshi122footnotemark: 2 1Purdue University Equal contribution.Corresponding author. Abstract Multi-agent reinforcement learning (MARL) has made significant progress in recent years, but most algorithms still rely on a discrete-time Markov Decision Process (MDP) with fixed decision intervals. This formulation is often ill-suited for complex multi-agent dynamics, particularly in high-frequency or irregular time-interval settings, leading to degraded performance and motivating the development of continuous-time MARL (CT-MARL). Existing CT-MARL methods are mainly built on Hamilton–Jacobi–Bellman (HJB) equations. However, they rarely account for safety constraints such as collision penalties, since these introduce discontinuities that make HJB-based learning difficult. To address this challenge, we propose a continuous-time constrained MDP (CT-CMDP) formulation and a novel MARL framework that transforms discrete MDPs into CT-CMDPs via an epigraph-based reformulation. We then solve this by proposing a novel physics-informed neural network (PINN)-based actor–critic method that enables stable and efficient optimization in continuous time. We evaluate our approach on continuous-time safe multi-particle environments (MPE) and safe multi-agent MuJoCo benchmarks. Results demonstrate smoother value approximations, more stable training, and improved performance over safe MARL baselines, validating the effectiveness and robustness of our method. Code is available at this link. 1 Introduction MARL has achieved remarkable success in diverse domains, ranging from strategic games (Samvelyan et al., 2019; Vinyals et al., 2019), multi-robot coordination (Haydari and Yılmaz, 2020; Kuyer et al., 2008), and wireless communication (Wang et al., 2023). These advances demonstrate the potential of MARL as a powerful framework for solving complex cooperative and competitive decision-making problems. Despite these achievements, most existing MARL algorithms are formulated in discrete time and fundamentally rely on the Bellman equation (Bellman, 1966). This formulation often assumes fixed time intervals between decision steps, which is adequate in settings where the decisions naturally occur at uniform time intervals. However, this assumption is not well-suited for complex high-frequency domains such as autonomous driving (Kiran et al., 2021; Chen et al., 2021), financial trading (Shavandi and Khedmati, 2022), where decision-making requires continuous-time control. In such cases, discrete-time RL often struggles to learn accurate policy (Doya, 2000; Mukherjee and Liu, ), as fixed-step discretization fails to represent non-uniform temporal dynamics, resulting in degraded performance and unstable learning (Tallec et al., 2019; Park et al., 2021; De Asis and Sutton, 2024). These limitations highlight the necessity of developing an alternative framework beyond discrete-time Bellman equations, which is compatible with CT-MARL. Recent studies (Wang et al., 2025) have explored the HJB equations to solve CT-MARL problems. The HJB can be viewed as the continuous-time analogue of the Bellman recursion, where the value function is characterized as the viscosity solution of a nonlinear Partial Differential Equation (PDE) (Shilova et al., 2024). In practice, PINNs have emerged as a common approach to approximate HJB solutions: they train neural networks to minimize HJB PDE residuals and leverage gradient-consistent signals for policy improvement (Mukherjee and Liu, ; Meng et al., 2024). This formulation eliminates the need for fixed time discretization and enables MARL to operate in continuous-time domains. However, in safety CT-MARL settings, state constraints (e.g., when they are treated as collision penalties) introduce value discontinuities, making it difficult for HJB-based PINNs to approximate the value functions accurately (Zhang et al., 2024). To address these challenges, we first cast safe CT-MARL as a CT-CMDP with explicit state constraints. We then introduce a revised epigraph reformulation that augments the system with an auxiliary state z, transforming the discontinuous constrained values into a continuous form suitable for PDE-based learning. On top of this reformulation, we adopt an actor–critic framework to learn values and policies under continuous-time state constraints. Specifically, we improve epigraph-based training by integrating the inner and outer optimization into a unified scheme. At each rollout, we compute the optimal auxiliary state z∗z^* and uses it directly for training, while keeping all networks z-independent. This design avoids the noise of random z sampling, yields more accurate policy updates, and eliminates costly root-finding at execution. Our main contributions are summarized as follows. (1) To the best of our knowledge, this is the first work to explicitly incorporate state constraints into the formulation of CT-MARL. We introduce an epigraph-based reformulation to bounds discounted cumulative cost and state constraints within a unified objective, effectively transforming discontinuous values into continuous ones. (2) We design an improved epigraph training scheme that integrates inner and outer optimization, providing more stable learning signals and removing the need for costly root-finding algorithms. (3) We prove the existence and uniqueness of viscosity solutions for epigraph-based HJB PDEs, providing theoretical support for our method. Extensive experiments on adapted continuous-time safe MPE and multi-agent MuJoCo benchmarks further demonstrate that our approach consistently outperforms current safe MARL methods. 2 Related Work 2.1 Continuous-Time Reinforcement Learning Discrete-time reinforcement learning (DTRL) often performs poorly in continuous-time environments, particularly when decision intervals are irregular (Tallec et al., 2019; Park et al., 2021; De Asis and Sutton, 2024). Consequently, continuous-time reinforcement learning (CTRL) has received growing attention as a more suitable framework for such problems (Doya, 2000; Yildiz et al., 2021; Wang et al., 2020; Bradtke and Duff, 1994; Jia and Zhou, 2022b; a). Most existing studies focus on the single-agent setting, proposing various approaches for value function approximation (Mukherjee and Liu, ; Wallace and Si, 2023; Lee and Sutton, 2021). For example, Mukherjee and Liu employ PINNs to approximate the value function and guide a PPO-based policy update, while Jia and Zhou (2022a) address stochastic dynamics through a Martingale loss designed for stochastic differential equations. In contrast, research on CT-MARL remains limited. Prior works (Luviano and Yu, 2017; Jiang et al., 2023) have considered multi-agent problems in continuous time, but largely in application-specific contexts rather than as general-purpose algorithms. The study in Wang et al. (2025) represents the first systematic attempt to design CT-MARL methods, combining PINNs with value gradient iteration to improve value approximation and performance. However, these approaches still inherit the limitations of PINNs that they can only approximate smooth value functions and therefore neglect safety constraints. 2.2 Multi-agent systems with Safety Concerns Multi-agent scenarios often raise critical safety concerns, and directly learning under combined reward and safety signals poses significant challenges. A number of studies have explored safe MARL frameworks to address these issues (Gu et al., 2023b; ElSayed-Aly et al., 2021; Gu et al., 2024; Shalev-Shwartz et al., 2016). For instance, Chow et al. (2018) employ primal–dual methods to enforce safety constraints, while Althoff et al. (2019) adopt a trust-region approach. Gu et al. (2021) introduce MACPO and MAPPO-Lagrange, which provide theoretical guarantees for both monotonic reward improvement and safety constraint satisfaction. In addition, Zhang et al. (2025b) leverage epigraph forms to formulate multi-agent safe optimal control problems, improving stability during training. However, these approaches are primarily developed in discrete-time settings, which limits their ability to capture continuous-time dynamics. Some efforts have incorporated safety into continuous-time multi-agent systems (e.g., Tayal et al. (2025)), but they assume fully known system dynamics and rely on optimal control algorithms, significantly restricting applicability. In more realistic scenarios, where dynamics are only partially known or highly complex, such methods fail to provide practical solutions. Existing methods remain limited in handling discontinuities and safety constraints in CT-MARL. Discrete-time safe MARL algorithms provide theoretical guarantees but do not naturally extend to continuous dynamics, while continuous-time approaches struggle with discontinuous value functions. To address these challenges, we propose an epigraph-based reformulation that unifies safety constraints and standard cost functions within a single objective, enabling principled and stable learning in CT-MARL. 3 Methodology In this section, we present our epigraph-based PINN actor–critic iteration (EPI) for solving CT-MARL with state constraints. 1) We first formalize the learning problem as CT-CMDP. Secondly, 2) we reformulate the CT-CMDP using an epigraph form. By introducing an auxiliary state z to augment system states, this reformulation converts discontinuous value functions into continuous ones. Building on this reformulation, 3) we develop an actor-critic learning architecture that aligns with the epigraph inner-outer optimization scheme. Specifically, the outer optimization computes the optimal auxiliary state z∗z^* along the rollout, ensuring that the critic captures the tightest feasible trade-off between return and safety. Based on this, the inner optimization trains the critic using PINNs, which jointly update the return and constraint networks together with z∗z^* to approximate the epigraph-based value function. This stabilized critic then serves as the foundation for actor training: we derive an advantage function consistent with the epigraph-based HJB PDEs, which provides the key learning signal for policy improvement. 3.1 Problem Formulation 3.1.1 Continuous-time Constrained Markov Decision Process We consider a CT-CMDP problem, formally defined by the tuple ℳ=⟨,ii=1N,N,f,lii=1N,c,tkk≥0,γ⟩,M= ,\U_i\_i=1^N,N,f,\l_i\_i=1^N,c,\t_k\_k≥ 0,γ , (1) where ⊆ℝnX ^n is the global state space, and =1×⋯×N⊆ℝmU=U_1×…×U_N ^m is the joint control space for N agents. The system evolves according to time-invariant nonlinear dynamics x˙(t)=f(x(t),u(t)) x(t)=f(x(t),u(t)) with x(0)=x0x(0)=x_0, where f:×→f:X×U . Each agent i applies a decentralized policy πi:×[0,∞)→i _i:X×[0,∞) _i, and the joint policy is denoted as π=(π1,…,πN)π=( _1,…, _N). All agents share the non-negative cost function l=∑i=1Nlil= _i=1^Nl_i, where li:×i→ℝl_i:X×U_i is the independent cost function of agent i. The system is further subject to state-dependent safety constraints specified by a function c:→ℝc:X , with the feasible set defined as ℱ=x∈|c(x)≤0F=\x \;|\;c(x)≤ 0\. Control actions are updated at irregular decision times tkk≥0\t_k\_k≥ 0, with strictly positive intervals τk=tk+1−tk _k=t_k+1-t_k. γ∈(0,1]γ∈(0,1] is the discount factor. Throughout the paper, we assume that iU_i is compact and convex, f and c are Lipschitz continuous, and lil_i is Lipschitz continuous and bounded. The joint objective is to minimize the cumulative cost under joint control input u=(u1,…,uN)u=(u_1,…,u_N) subject to state constraints c(x)c(x): v(x) v(x) =minu∈∫t∞γτ−tl(x(τ),u(τ))τ = _u _t^∞γ^τ-t\,l(x(τ),u(τ))\,dτ (2) s.t.c(x(τ))≤0,∀τ≥t. .t. c(x(τ))≤ 0, ∀τ≥ t. 3.1.2 Epigraph Reformulation The value becomes discontinuous (Altarovici et al., 2013) when state constraints are violated in Eq. 2, which hinders the convergence of HJB-based PINN training. To address this, we leverage an epigraph reformulation that converts value in Eq. 2 into a continuous representation. Definition 1 (Epigraph Reformulation). We introduce an auxiliary state variable z(t)∈ℝz(t) to reformulate Eq. 2 using the epigraph forms. Here, z follows the dynamic z˙(t)=−l(x(t),u(t))−lnγ⋅z(t) z(t)=-l(x(t),u(t))- γ· z(t). Therefore, the auxiliary value function is defined as V(x,z)=minu∈maxmaxτ∈[t,∞]c(x(τ)),∫t∞γτ−tl(x(τ),u(τ))τ−z,V(x,z)= _u \ _τ∈[t,∞]c(x(τ)), _t^∞γ^τ-tl(x(τ),u(τ))dτ-z \, (3) Lemma 3.1 (Value Equivalence). Suppose the assumptions in Sec. 3.1.1 hold. For all (t,x,z)∈[0,∞)×ℝ(t,x,z)∈[0,∞)×X×R, the constrained value v and auxiliary value V are related by v(x)=minz∈ℝ|V(x,z)≤0. v(x)= \z ~|~V(x,z)≤ 0\. (4) Here, the sub-zero level set of auxiliary value V becomes the epigraph of the constrained value v. The proof is listed in Appendix A.1. Lemma 3.2 (Optimality Condition). For all (t,x,z)∈[0,∞)×ℝ(t,x,z)∈[0,∞)×X×R, consider a small enough h>0h>0, the auxiliary value function V satisfies V(x,z)=minu∈maxmaxτ∈[t,t+h]c(x(τ)),γhV(x(t+h),z(t+h)).V(x,z)= _u \ _τ∈[t,t+h]c(x(τ)),γ^hV(x(t+h),z(t+h)) \. (5) The proof is listed in Appendix A.2. Theorem 3.3 (Epigraph-based HJB PDE). Let V:×ℝ→ℝV:X×R\!→\!R be the auxiliary value function defined in Eq. 3. Then V is the unique viscosity solution of the following HJB PDE for all (t,x,z)∈[0,∞)×ℝ(t,x,z)∈[0,∞)×X×R. maxmaxτ∈[t,∞]c(x)−V(x,z),minu∈ℋ(x,z,∇xV,∂zV)=0, \ _τ∈[t,∞]c(x)-V(x,z),\; _u H(x,z, _xV, _zV) \=0, (6) where ℋ(x,z,∇xV,∂zV)H(x,z, _xV, _zV) is Hamiltonian and satisfies ℋ=∇xV⋅f(x,u)−∂zV⋅l(x,u)+lnγ⋅VH= _xV· f(x,u)- _zV· l(x,u)+ γ· V and optimal control u∗=argminu∈ℋu^*= *arg\,min_u H. The derivation proof is provided in Appendix A.3. 3.2 Epigraph Learning Framework Figure 1: Overview of the proposed epigraph-based CT-MARL framework. The pipeline begins with data collection, where individual agent rollouts are aggregated into a centralized rollout RX_R for the training; the outer optimization computes optimal z∗z^* to balance discounted cumulative cost and safety constraints; the inner optimization corresponds to critic learning, where return networks Vψret(x)V^ret_ψ(x) and constraint value networks Vϕcons(x)V^cons_φ(x) are optimized jointly with the optimal auxiliary state z∗z^*; and actor learning leverages the advantage function to improve policies. ‘ As illustrated in Fig. 1, our framework integrates the epigraph-based inner-outer optimization (Zhang et al., 2025b) into the actor-critic paradigm. The outer loop updates z∗z^* along the rollout by solving Eq. 7, ensuring that the critic is trained with the minimal z that simultaneously satisfies both costs and safety constraints. minz∈ℝzs.t.minπmaxsupτ≥tc(x(τ)),∫t∞γτ−tl(x(τ),π(τ))τ−z≤0. _z \;z .t. _π \ _τ≥ tc(x(τ)),\; _t^∞γ^τ-t\,l(x(τ),π(τ))\,dτ-z \≤ 0. (7) In the inner loop, the critic is trained as follows: the return and constraint value networks (Vψret(x)V^ret_ψ(x) and Vϕcons(x)V^cons_φ(x)) are optimized using z∗z^* to approximate the auxiliary value function V~(x,z∗) V(x,z^*). This stabilized critic subsequently supplies the learning signals for decentralized actors, which map local observations to continuous-time policies under the standard centralized training decentralized execution setup (Foerster et al., 2018; Lowe et al., 2017). We next describe the revised outer optimization in detail, focusing on solving the optimal auxiliary state z∗z^* that trades off discounted cost against safety violations without costly root-finding algorithms (So and Fan, 2023; So et al., 2024; Zhang et al., 2025b). 3.2.1 Revised Outer Optimization We seek the minimal z such that the epigraph-based value V remains non-positive, as defined in Eq. 4. Using the return and constraint value network learned by the critic, the optimal auxiliary state z∗z^* can be found by solving for the minimal feasible solution: z∗=minz∈ℝ|maxVϕcons(x),Vψret(x)−z≤0,z^*= \z ~|~ \V^cons_φ(x),\,V^ret_ψ(x)-z\≤ 0 \, (8) where return value network Vψret(x)V^ret_ψ(x) that approximates the discounted cumulative cost ∫t∞γτ−tl(x(τ),π(τ))τ _t^∞γ^τ-t\,l(x(τ),π(τ))\,dτ, and constraint value network Vϕcons(x)V^cons_φ(x) represents the violation for worst-case future constraints supτ≥tc(x(τ)) _τ≥ tc(x(τ)). In previous epigraph formulations (Tayal et al., 2025; Zhang et al., 2025b), the outer problem is solved during the execution phase: z is sampled along the rollouts during training, and z∗z^* is computed at execution time via root-finding (Stoer et al., 1980). This design has two drawbacks in CT-MARL: (1) the random sampling of z introduces nonstationary noise that destabilizes the updates of actor and critic and further leads to poor convergence; (2) at execution, root-finding must be performed at every step, which is computationally expensive and often incompatible with real-time requirements. In contrast, we design the return and constraint value networks as functions of the states x solely. We then integrate the outer optimization into actor-critic training: for each episode, z∗z^* is computed using the current learned value V~ V along the predicted rollout. The actor is then trained against a z-independent critic, yielding a z-independent policy π(x)π(x). This design ensures stable actor training, and enables real-time deployment by eliminating the need for root-finding during execution. Since the critic’s value networks are z-independent, the outer optimization is simplified to a scalar search for z∗z^*, which adds negligible cost to model training. 3.2.2 Inner Optimization with Critic Learning The inner optimization is responsible for updating the PINN-based critic networks. Given a task-dependent range [zmin,zmax][z_ ,z_ ], the outer optimization computes z∗z^*, which is then clipped to this range (i.e., z∗←minmaxz∗,zmin,zmaxz^*← \ \z^*,z_ \,z_ \) before being used to train the critic module. The critic consists of two value networks: a return value network Vψret(x)V^ret_ψ(x), and a constraint value network Vϕcons(x)V^cons_φ(x). Together with the computed z∗z^* from Eq. 8, these define the composite epigraph-based value function: V~(x,z∗)=maxVϕcons(x),Vψret(x)−z∗. V(x,z^*)\;=\; \V^cons_φ(x),\;\;V^ret_ψ(x)-z^* \. (9) To ensure stable and accurate training, we employ three complementary losses: (i) Residual Loss. We use PINN architecture (Mukherjee and Liu, ) to approximate the value function governed by epigraph-based HJB PDEs, and introduce a residual loss that penalizes violations of the corresponding PDEs: ℒResidual=(maxc(x)−V~,minu∈[∇xV~⋅f(x,u)−∂zV~⋅l(x,u)+lnγ⋅V~])2.L_Residual= ( \c(x)- V,\; _u [ _x V· f(x,u)- _z V· l(x,u)+ γ· V ] \ )^2. (10) (i) Target Loss. In standard PINNs, a boundary loss is combined with the PDE residual to approximate PDE solutions (Cai et al., 2021; Raissi et al., 2019). In the infinite-horizon setting, however, no boundary condition is available, and training the critic only on residuals is insufficient: optimization may converge, but to incorrect PDE solutions (Wang et al., 2022). To address this, we add a rollout-based target loss that measures the discrepancy between the epigraph-based value approximation with a numerical target defined by Eq. 3. For each episode, the current value V~ V generates a closed-loop trajectory x(τ),u(τ)τ=t∞\x(τ),u(τ)\_τ=t^∞; from this trajectory we construct the target Vtgt(x,z)=maxmaxτ∈[t,∞]c(x(τ)),∫t∞γτ−tl(x(τ),u(τ))τ−z∗V_tgt(x,z)= \ _τ∈[t,∞]c(x(τ)),\;\; _t^∞γ^τ-t\,l(x(τ),u(τ))\,dτ-z^* \ and minimize the squared error: ℒTarget=(Vtgt(x,z∗)−maxVϕcons(x),Vψret(x)−z∗)2.L_Target\;=\; (V_tgt(x,z^*)\;-\; \\,V^cons_φ(x),\;V^ret_ψ(x)-z^*\,\ )^2. (11) (i) Value Gradient Iterations. Standard PINN training in multi-agent settings often struggles to approximate accurate value functions, primarily because the learned value gradients are inaccurate or unstable (Wang et al., 2025; Zhang et al., 2024). The VGI techniques (Eberhard et al., 2025; Wang et al., 2025) are designed to enhance the quality of learned value gradients. In our framework, accurate gradients ∇xV(x) _xV(x) are crucial for precise value approximations, which in turn affect actor learning and ultimately determine the quality of the resulting policies. To establish the theoretical basis of this module, we follow Theorem 3.4 in Bokanowski et al. (2021) and Theorem 2 in Hermosilla and Zidani (2023): ∇xtV~(xt)=∇xt(χ(xt)l(xt,ut)+(1−χ(xt))c(xt))Δt+γΔt∇xt+ΔtV~(xt+Δt)⋅∇xtf(xt,ut), _x_t V(x_t)= _x_t(χ(x_t)l(x_t,u_t)+(1-χ(x_t))c(x_t)) t+γ t _x_t+ t V(x_t+ t)· _x_tf(x_t,u_t), (12) where the χ(xt):= 1Vψret(xt)−zt≥Vϕcons(xt)χ(x_t)\;:=\;1\\,V^ret_ψ(x_t)-z_t\;≥\;V^cons_φ(x_t)\. As shown in Eq. 12, the value gradient satisfies a recursive relation coupling the local cost gradient with the backpropagated dynamics term. The overall critic objective is a weighted sum of the three losses as: ℒCritic=λresℒResidual+λtgtℒTarget+λvgiℒVGI,L_Critic= _resL_Residual+ _tgtL_Target+ _vgiL_VGI, (13) where the weights (λres,λtgt,λvgi)( _res, _tgt, _vgi) are selected to keep the losses on comparable scales and are determined via grid search. 3.2.3 Actor Learning After introducing the inner-outer optimization for critic learning, we turn to the actor learning. We first define the epigraph-based Q-function, which is used for deriving policy update rules. Definition 2 (Epigraph-based Q-function). Following the definition in (So and Fan, 2023), for any state-action pair (xt,ut)(x_t,u_t) and auxiliary state ztz_t, the epigraph-based Q-function is defined Q(xt,zt∗,ut)=maxc(xt),γhV(xt+h,zt+h∗).Q(x_t,z_t^*,u_t)= \c(x_t)\,,\,γ^hV(x_t+h,z_t+h^*) \. (14) where xt+hx_t+h and zt+h∗z_t+h^* are the states and optimal auxiliary state at t+ht+h, respectively. h is a short time interval. Lemma 3.4 (Epigraph-based advantage function). The epigraph-based advantage function A(xt,zt∗,ut)=Q(xt,zt∗,ut)−V(xt,zt∗)A(x_t,z_t^*,u_t)=Q(x_t,z_t^*,u_t)-V(x_t,z_t^*) (15) is equivalent to epigraph-based HJB PDE when h→0h→ 0 A(xt,zt∗,ut)=maxc(xt)−V(xt,zt∗),∇xtV⋅f(xt,ut)−∂ztV⋅l(xt,ut)+lnγ⋅V.A(x_t,z_t^*,u_t)= \c(x_t)-V(x_t,z_t^*), _x_tV· f(x_t,u_t)- _z_tV· l(x_t,u_t)+ γ· V\. (16) In practice, evaluating the epigraph-based advantage in Eq. 16 requires knowledge of the true dynamics f(x,u)f(x,u) and cost function l(x,u)l(x,u). Since these quantities are generally unknown in model-free reinforcement learning, we replace them with neural networks that are jointly trained alongside the actor. The derivation of the epigraph-based advantage function is listed at Appendix A.4. Dynamics and Cost Networks. To assist with the policy training, we employ two neural networks: a dynamics network fξ(x,u,Δt)f_ξ(x,u, _t) that predicts the next state x′x given the current state–action pair, and a cost network lϕ(x,u,Δt)l_φ(x,u, _t) that estimates the instantaneous stage cost. Both models are trained via supervised regression using observed transitions (x,u,x′,l)(x,u,x ,l) from the environment. Specifically, the training losses are ℒdyn(ξ)=‖fξ(x,u,Δt)−x′‖2,ℒrew(ϕ)=(lϕ(x,u,Δt)−l(x,u))2,L_dyn(ξ)= \|f_ξ(x,u, _t)-x \|^2, _rew(φ)= (l_φ(x,u, _t)-l(x,u) )^2, (17) where x′x is the observed next state and l(x,u)l(x,u) is the empirical cost signal. Equivalently, the dynamics learning can be interpreted as approximating the continuous-time derivative dynamics (fξ(x,u,Δt)−x)/Δt(f_ξ(x,u, _t)-x)/ t. Actor Update with Learned Models. By substituting V~(x,z∗) V(x,z^*), fξf_ξ and lϕl_φ into the epigraph advantage expression Eq. 16, we obtain a differentiable surrogate Aθ(x,z∗,u)=maxc(x)−V~,∇xV~⋅fξ−∂zV~⋅lϕ+lnγ⋅V~.A_θ(x,z^*,u)= \c(x)- V, _x V· f_ξ- _z V· l_φ+ γ· V\. (18) The actor πθ(u∣x) _θ(u x) is updated by minimizing the expected surrogate advantage: ℒactor(θ)=x∼R,u∼πθ(⋅∣x)[Aθ(x,z∗,u)],L_actor(θ)=E_x _R,u _θ(· x) [\,A_θ(x,z^*,u)\, ], (19) where RX_R is the sampled data along the rollout. Specifically, we adopt a centralized-training decentralized-execution structure: each agent’s actor πi(oi,Δt) _i(o_i, t) takes its local observation oio_i as input, while the training signal is derived from the state x. The overall training pipeline is summarized in Algorithm 1 in Appendix B. 4 Experimental Results Figure 2: Overall results for adapted MPE environments. We organize our empirical study around the following research questions: Q1. How well does our method balance discounted cumulative cost and constraint satisfaction compared to state-of-the-art baselines? Q2. How does the different loss component in critic learning contribute to stable training and accurate value approximations? Q3. How does performance change when training with versus without the epigraph reformulation? Q4. How sensitive is the epigraph formulation to the choice of the auxiliary variable z during training? Q5. How robust is the method under stochastic disturbances, and how does performance degrade under model-mismatch noise? Q6. How does the performance change under different discretization resolutions Δt t? 4.1 Benchmarks and baselines. To evaluate our approach under continuous-time environments with safety constraints, we consider two adapted benchmarks: the safe continuous-time MPE (Lowe et al., 2017; Wang et al., 2025) and continuous-time Safe MA-MuJoCo (Gu et al., 2023a; Wang et al., 2025). In MPE, we design several scenarios including Corridor, Formation, Line, Target, Simple Spread, and Cooperative Predator–Prey. These tasks typically place agents in environments with obstacles and require them to avoid both collisions with obstacles and collisions with other agents while navigating or pursuing their objectives. In MuJoCo, we adapt several scenarios such as Half Cheetah and Ant into continuous-time versions and introduce randomly placed walls as obstacles. The agents must coordinate to move forward efficiently while avoiding crashing into walls, ensuring that the learned policies account for both locomotion and safety considerations. Figure 3: Performance of constraints and cost over MPE settings. Lastly, we design a didactic example based on a constrained coupled oscillator, which admits an analytical ground-truth solution for both value functions and actions. This example provides a transparent testbed to directly validate the correctness of our learned critics against exact solutions. Full details of the agent setups, metrics, state and action spaces, and cost specifications are provided in the Appendix C. Figure 4: Overall results for adapted multi-agent MuJoCo environments. We compare our approach EPI with MACPO (Gu et al., 2021), MAPPO-Lag (Gu et al., 2021), SAC-Lag (Haarnoja et al., 2018), EPPO (Zhang et al., 2025b) and CBF (Zhang et al., 2025a). The first three represent the most widely used families of safe MARL algorithms: trust-region based methods (MACPO) and Lagrangian based methods (MAPPO-Lag, SAC-Lag), covering both on-policy and off-policy learning. We also include EPPO as an epigraph-based baseline that follows the traditional epigraph optimization framework. We additionally include a control barrier function (CBF) baseline, which enforces safety through model-based barrier certificates and is commonly used in safe multi-agent control. Although these algorithms were originally developed in the discrete-time setting, we adapt them to continuous time by equipping their critics with the same PDE residual loss used in our method. Since the performance gap between discrete-time and continuous-time algorithms has already been well studied (Tallec et al., 2019; De Asis and Sutton, 2024), our baselines focus only on isolating the effect of different safety mechanisms (trust-region, Lagrangian, or epigraph). Figure 5: Ablation study of different loss terms in critic network over MPE. 4.2 Results Analysis In this section, we present a systematic analysis of the results, addressing each research question in turn. Q1. Our method consistently outperforms all baselines across both adapted MPE and MuJoCo environments in Fig. 2 and Fig. 4. We adopt the same reward design commonly used in prior safe MARL works such as MACPO (Gu et al., 2021). Specifically, the reward is the combination of the task cost provided by the environment (e.g., distance to the target in MPE) and the safety penalty provided by the environment (e.g., collision penalty between agents or with obstacles), as detailed in Appendix C, which directly reflects performance under both objectives. In Fig. 3 and 6, each point corresponds to the average performance of one algorithm, with horizontal and vertical bars denoting standard deviations. Since the goal is to minimize both cost and constraint violations, the lower-left corner of each panel represents the desirable region. These results show that our algorithm EPI achieves nearly the lowest cost and constraint violation in every scenarios. Specifically, EPPO often remains stuck at suboptimal solutions because it randomly samples the auxiliary state z instead of using z∗z^* for model training, introducing noise that disrupts policy updates and prevents stable convergence. Figure 6: Performance of constraints and cost over MuJoCo settings. MACPO enforces constraints through a hard trust-region style update, which yields strong violation rejection but tends to be overly conservative. SAC-Lag and MAPPO-Lag rely on Lagrangian relaxation, which is known to suffer from instability when balancing objectives under tight safety requirements (Zhang et al., 2025b). CBF achieves reasonable constraint-violation levels but tends to be conservative. The CBF condition relies on the gradient of a learned barrier function ∇B(x)∇ B(x), approximation errors in this component can distort the effective safe set and degrade the overall performance. Q2. The ablation results in Fig. 5 clearly demonstrate the importance of each loss component in critic learning. It presents the cumulative reward performance of our full method compared with its ablation variants across representative continuous-time MPE tasks. Figure 7: Weighted loss performance. Removing the target loss or the VGI loss significantly degrades performance, whereas removing the residual loss has only a minor effect. This difference stems from the fact that, unlike existing HJ-based PINN methods (Zhang et al., 2024; Tayal et al., 2025; Cai et al., 2021) that address finite-horizon problems with boundary conditions, our framework targets the infinite-horizon setting where no such boundary conditions are available. In this case, the target loss serves as an anchor to stabilize value approximations, ensuring that value function V(x)V(x) does not drift arbitrarily, while the VGI loss enforces consistency of the learned value gradients, which are crucial for both accurate value approximations and policy improvement. In contrast, the HJB residual loss mainly regularizes the PDE structure, but its role becomes less critical once the value gradients are optimized by VGI. As a result, the removal of VGI has a severe impact, since inaccurate value gradients directly harm both critic accuracy and actor updates, while the residual loss contributes less critically to overall training stability. The grouped bars in Fig. 7 report the average distance to the target (lower is better) for three MPE tasks (Formation, Line, and Target) under different loss weightings. The balanced setting (EPI) attains the smallest distance in all tasks and shows the tightest variability. Over-emphasizing any single component degrades performance: increasing the target loss weight is particularly harmful on Line (large increase in distance), while overweighting (×20× 20) the residual or the VGI loss also worsens results relative to EPI, though to a lesser extent. These ablations support using the balanced weighting adopted by EPI. Q3. We generate one trajectory using EPI and collect the visited states. On these same states, we compare the value and policy from three methods: EPI, Ground Truth via the LQR method (details in the Appendix C), and an ablation without the epigraph reformulation, where the state constraint is treated as a collision penalty added to the cost function l, making the value function discontinuous. Figure 8: Performance with and without epigraph reformulation. For Ground Truth, the value is computed as the discounted cumulative cost. While for the EPI and ablation without the epigraph form, the value is predicted through the trained value network. EPI closely tracks the Ground Truth in both value and actions for both agents, indicating accurate value approximation and stable control policies. In contrast, the ablation without the epigraph form exhibits severely mis-scaled value predictions (we plot it after a ×120×\! 120 scaling to share the same y-axis) and noticeably unstable actions, which in practice are more likely to violate constraints because the discontinuous value function is not addressed by the epigraph form. The poor performance of the ablation without epigraph stems from the discontinuity of the value function when state constraints are directly encoded as hard penalties. Such discontinuities are notoriously difficult to approximate with neural networks, leading to severely mis-scaled value predictions and unstable gradients for policy updates. By contrast, the epigraph reformulation converts the discontinuous penalty into a continuous and smooth upper-bound optimization, which stabilizes critic learning and yields reliable policies. Q4. To better understand how model performance depends on z, we test two MPE tasks (Formation and Line) under different values of z. Specifically, we train the EPI model with z∈z∗−0.5zmax,z∗−0.2zmax,z∗,z∗+0.2zmax,z∗+0.5zmaxz∈\z^*-0.5z_ ,z^*-0.2z_ ,z^*,z^*+0.2z_ ,z^*+0.5z_ \. Fig. 9 reports the results, where the x-axis indicates cost and the y-axis denotes the constraint violation rate. Figure 9: Sensitivity test of different z choices. Compared with the optimal auxiliary state z∗z^*, using a suboptimal z shifts the trade-off between cost and constraint satisfaction, often resulting in either much higher violation rates or larger costs. Specifically, a smaller z (e.g., z∗−0.2zmaxz^*-0.2z_ , z∗−0.5zmaxz^*-0.5z_ ) significantly increases the violation rate while only slightly reducing cost. Getting back to the epigraph form maxVϕcons(x),Vψret(x,z)−z \V^cons_φ(x),\;V^ret_ψ(x,z)-z\, a smaller z makes Vret(x,z)−zV^ret(x,z)-z lager than Vϕcons(x)V^cons_φ(x), so the return term dominates in the epigraph form. As a result, the optimization prioritizes reward improvement while neglecting constraint satisfaction, leading to frequent violations. In constrast, when z is larger than z∗z^* (e.g., z∗+0.2zmaxz^*+0.2z_ , z∗+0.5zmaxz^*+0.5z_ ), the term Vret(x,z)−zV^ret(x,z)-z becomes smaller than Vϕcons(x)V^cons_φ(x), making constraint value dominate in the epigraph form. This forces the critic and actor to emphasize constraint satisfaction, which reduces violations but increases cost. 5 Conclusion In this paper, we propose an epigraph-based framework for CT-MARL that addresses the challenges of balancing reward maximization with constraint satisfaction. By reformulating the problem through the epigraph forms, we introduced an inner–outer optimization procedure that enables stable critic learning and effective policy updates. Our design further integrates different losses in critic learning, including target, residual, and VGI losses, to anchor value approximations and improve gradient accuracy in the infinite-horizon setting. Through extensive experiments in both adapted MPE and MuJoCo benchmarks, we demonstrated that our method consistently outperforms state-of-the-art baselines in terms of both cost reduction and constraint satisfaction. Acknowledgement This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-24-1-0233. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the United States Air Force. Ethics Statement This work focuses on decision-making for the continuous-time constrained MDP problems. All experiments are conducted entirely in simulation and do not involve human subjects or personal data. References A. Altarovici, O. Bokanowski, and H. Zidani (2013) A general Hamilton-Jacobi framework for non-linear state-constrained control problems. ESAIM: Control, Optimisation and Calculus of Variations 19 (2), p. 337–357. Cited by: §3.1.2. M. Althoff, A. Giusti, S. B. Liu, and A. Pereira (2019) Effortless creation of safe robots from modules through self-programming and self-verification. Science Robotics 4 (31), p. eaaw1924. Cited by: §2.2. R. Bellman (1966) Dynamic programming. science 153 (3731), p. 34–37. Cited by: §1. O. Bokanowski, A. Désilles, and H. Zidani (2021) Relationship between maximum principle and dynamic programming in presence of intermediate and final state constraints. ESAIM: Control, Optimisation and Calculus of Variations 27, p. 91. Cited by: §3.2.2. S. Bradtke and M. Duff (1994) Reinforcement learning methods for continuous-time markov decision problems. Advances in neural information processing systems 7. Cited by: §2.1. S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis (2021) Physics-informed neural networks (PINNs) for fluid mechanics: A review. Acta Mechanica Sinica 37 (12), p. 1727–1738. Cited by: §3.2.2, §4.2. J. Chen, S. E. Li, and M. Tomizuka (2021) Interpretable end-to-end urban autonomous driving with latent deep reinforcement learning. IEEE Transactions on Intelligent Transportation Systems 23 (6), p. 5068–5078. Cited by: §1. Y. Chow, M. Ghavamzadeh, L. Janson, and M. Pavone (2018) Risk-constrained reinforcement learning with percentile risk criteria. Journal of Machine Learning Research 18 (167), p. 1–51. Cited by: §2.2. K. De Asis and R. S. Sutton (2024) An idiosyncrasy of time-discretization in reinforcement learning. arXiv preprint arXiv:2406.14951. Cited by: §1, §2.1, §4.1. K. Doya (2000) Reinforcement learning in continuous time and space. Neural computation 12 (1), p. 219–245. Cited by: §1, §2.1. O. Eberhard, C. Vernade, and M. Muehlebach (2025) A Pontryagin perspective on reinforcement learning. In 7th Annual Learning for Dynamics &\& Control Conference, p. 233–244. Cited by: §3.2.2. I. ElSayed-Aly, S. Bharadwaj, C. Amato, R. Ehlers, U. Topcu, and L. Feng (2021) Safe multi-agent reinforcement learning via shielding. arXiv preprint arXiv:2101.11196. Cited by: §2.2. L. C. Evans (2022) Partial differential equations. Vol. 19, American mathematical society. Cited by: §A.2, §A.3, §A.3. J. Foerster, G. Farquhar, T. Afouras, N. Nardelli, and S. Whiteson (2018) Counterfactual multi-agent policy gradients. In Proceedings of the AAAI conference on artificial intelligence, Vol. 32. Cited by: §3.2. S. Gu, J. G. Kuba, Y. Chen, Y. Du, L. Yang, A. Knoll, and Y. Yang (2023a) Safe multi-agent reinforcement learning for multi-robot control. Artificial Intelligence, p. 103905. Cited by: §4.1. S. Gu, J. G. Kuba, Y. Chen, Y. Du, L. Yang, A. Knoll, and Y. Yang (2023b) Safe multi-agent reinforcement learning for multi-robot control. Artificial Intelligence 319, p. 103905. Cited by: §2.2. S. Gu, J. G. Kuba, M. Wen, R. Chen, Z. Wang, Z. Tian, J. Wang, A. Knoll, and Y. Yang (2021) Multi-agent constrained policy optimisation. arXiv preprint arXiv:2110.02793. Cited by: Appendix C, §2.2, §4.1, §4.2. S. Gu, L. Yang, Y. Du, G. Chen, F. Walter, J. Wang, and A. Knoll (2024) A review of safe reinforcement learning: methods, theories and applications. IEEE Transactions on Pattern Analysis and Machine Intelligence. Cited by: §2.2. T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine (2018) Soft actor-critic: off-policy maximum entropy deep reinforcement learning with a stochastic actor. In International conference on machine learning, p. 1861–1870. Cited by: §4.1. A. Haydari and Y. Yılmaz (2020) Deep reinforcement learning for intelligent transportation systems: a survey. IEEE Transactions on Intelligent Transportation Systems 23 (1), p. 11–32. Cited by: §1. C. Hermosilla and H. Zidani (2023) Relationship between the maximum principle and dynamic programming for minimax problems. Applied Mathematics & Optimization 87 (2), p. 34. Cited by: §3.2.2. Y. Jia and X. Y. Zhou (2022a) Policy evaluation and temporal-difference learning in continuous time and space: a martingale approach. Journal of Machine Learning Research 23 (154), p. 1–55. Cited by: §2.1. Y. Jia and X. Y. Zhou (2022b) Policy gradient and actor-critic learning in continuous time and space: theory and algorithms. Journal of Machine Learning Research 23 (275), p. 1–50. Cited by: §2.1. Y. Jiang, W. Gao, J. Wu, T. Chai, and F. L. Lewis (2023) Reinforcement learning and cooperative H∞ output regulation of linear continuous-time multi-agent systems. Automatica 148, p. 110768. Cited by: §2.1. B. R. Kiran, I. Sobh, V. Talpaert, P. Mannion, A. A. Al Sallab, S. Yogamani, and P. Pérez (2021) Deep reinforcement learning for autonomous driving: a survey. IEEE transactions on intelligent transportation systems 23 (6), p. 4909–4926. Cited by: §1. L. Kuyer, S. Whiteson, B. Bakker, and N. Vlassis (2008) Multiagent reinforcement learning for urban traffic control using coordination graphs. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2008, Antwerp, Belgium, September 15-19, 2008, Proceedings, Part I 19, p. 656–671. Cited by: §1. D. Lee (2022) Safety-guaranteed autonomy under uncertainty. Ph.D. Thesis, University of California, Berkeley. Cited by: §A.1, §A.2, §A.3. J. Lee and R. S. Sutton (2021) Policy iterations for reinforcement learning problems in continuous time and space—fundamental theory and methods. Automatica 126, p. 109421. Cited by: §2.1. R. Lowe, Y. I. Wu, A. Tamar, J. Harb, O. Pieter Abbeel, and I. Mordatch (2017) Multi-agent actor-critic for mixed cooperative-competitive environments. Advances in neural information processing systems 30. Cited by: §3.2, §4.1. D. Luviano and W. Yu (2017) Continuous-time path planning for multi-agents with fuzzy reinforcement learning. Journal of Intelligent & Fuzzy Systems 33 (1), p. 491–501. Cited by: §2.1. Y. Meng, R. Zhou, A. Mukherjee, M. Fitzsimmons, C. Song, and J. Liu (2024) Physics-informed neural network policy iteration: algorithms, convergence, and verification. arXiv preprint arXiv:2402.10119. Cited by: §1. [32] A. Mukherjee and J. Liu Bridging Physics-Informed Neural Networks with Reinforcement Learning: Hamilton-Jacobi-Bellman Proximal Policy Optimization (HJBPPO). In ICML Workshop on New Frontiers in Learning, Control, and Dynamical Systems, Cited by: §1, §1, §2.1, §3.2.2. S. Park, J. Kim, and G. Kim (2021) Time discretization-invariant safe action repetition for policy gradient methods. Advances in Neural Information Processing Systems 34, p. 267–279. Cited by: §1, §2.1. M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, p. 686–707. Cited by: §3.2.2. M. Samvelyan, T. Rashid, C. S. De Witt, G. Farquhar, N. Nardelli, T. G. Rudner, C. Hung, P. H. Torr, J. Foerster, and S. Whiteson (2019) The starcraft multi-agent challenge. arXiv preprint arXiv:1902.04043. Cited by: §1. S. Shalev-Shwartz, S. Shammah, and A. Shashua (2016) Safe, multi-agent, reinforcement learning for autonomous driving. arXiv preprint arXiv:1610.03295. Cited by: §2.2. A. Shavandi and M. Khedmati (2022) A multi-agent deep reinforcement learning framework for algorithmic trading in financial markets. Expert Systems with Applications 208, p. 118124. Cited by: §1. A. Shilova, T. Delliaux, P. Preux, and B. Raffin (2024) Learning HJB viscosity solutions with PINNs for continuous-time reinforcement learning. In ICML 2024 Workshop: Foundations of Reinforcement Learning and Control–Connections and Perspectives, Cited by: §1. O. So and C. Fan (2023) Solving stabilize-avoid optimal control via epigraph form and deep reinforcement learning. arXiv preprint arXiv:2305.14154. Cited by: §3.2, Definition 2. O. So, C. Ge, and C. Fan (2024) Solving minimum-cost reach avoid using reinforcement learning. Advances in Neural Information Processing Systems 37, p. 30951–30984. Cited by: §3.2. J. Stoer, R. Bulirsch, R. Bartels, W. Gautschi, and C. Witzgall (1980) Introduction to numerical analysis. Vol. 1993, Springer. Cited by: §3.2.1. C. Tallec, L. Blier, and Y. Ollivier (2019) Making deep q-learning methods robust to time discretization. In International Conference on Machine Learning, p. 6096–6104. Cited by: §1, §2.1, §4.1. M. Tayal, A. Singh, S. Kolathaya, and S. Bansal (2025) A physics-informed machine learning framework for safe and optimal control of autonomous systems. arXiv preprint arXiv:2502.11057. Cited by: §2.2, §3.2.1, §4.2. O. Vinyals, I. Babuschkin, W. M. Czarnecki, M. Mathieu, A. Dudzik, J. Chung, D. H. Choi, R. Powell, T. Ewalds, P. Georgiev, et al. (2019) Grandmaster level in StarCraft I using multi-agent reinforcement learning. nature 575 (7782), p. 350–354. Cited by: §1. B. A. Wallace and J. Si (2023) Continuous-time reinforcement learning control: a review of theoretical results, insights on performance, and needs for new designs. IEEE Transactions on Neural Networks and Learning Systems. Cited by: §2.1. H. Wang, T. Zariphopoulou, and X. Y. Zhou (2020) Reinforcement learning in continuous time and space: a stochastic control approach. Journal of Machine Learning Research 21 (198), p. 1–34. Cited by: §2.1. S. Wang, X. Yu, and P. Perdikaris (2022) When and why PINNs fail to train: A neural tangent kernel perspective. Journal of Computational Physics 449, p. 110768. Cited by: §3.2.2. X. Wang, X. Li, J. Shao, and J. Zhang (2023) AC2C: adaptively controlled two-hop communication for multi-agent reinforcement learning. arXiv preprint arXiv:2302.12515. Cited by: §1. X. Wang, L. Zhang, H. Pu, A. H. Qureshi, and H. Li (2025) Continuous-time value iteration for multi-agent reinforcement learning. arXiv preprint arXiv:2509.09135. Cited by: §1, §2.1, §3.2.2, §4.1. C. Yildiz, M. Heinonen, and H. Lähdesmäki (2021) Continuous-time model-based reinforcement learning. In International Conference on Machine Learning, p. 12009–12018. Cited by: §2.1. L. Zhang, M. Ghimire, W. Zhang, Z. Xu, and Y. Ren (2024) Value approximation for two-player general-sum differential games with state constraints. IEEE Transactions on Robotics 40, p. 4631–4649. Cited by: §A.1, §A.2, §A.3, §1, §3.2.2, §4.2. S. Zhang, O. So, M. Black, and C. Fan (2025a) Discrete gcbf proximal policy optimization for multi-agent safe optimal control. arXiv preprint arXiv:2502.03640. Cited by: §D.6, §4.1. S. Zhang, O. So, M. Black, Z. Serlin, and C. Fan (2025b) Solving multi-agent safe optimal control with distributed epigraph form marl. arXiv preprint arXiv:2504.15425. Cited by: §2.2, §3.2.1, §3.2, §3.2, §4.1, §4.2. Appendix A Mathematical Proof A.1 Lemma 3.1: Equivalence of two value functions Proof. Following proofs in (Lee, 2022; Zhang et al., 2024)), Eq. 4 implies the following equivalence v(x)−z≤0⟺V(x,z)≤0v(x)-z≤ 0 V(x,z)≤ 0\\ To prove the above relation, we first start from v(x)−z≤0v(x)-z≤ 0, which implies that there exists a joint control input u∈u such that ∫t∞γτ−tl(x(τ),u(τ))τ−z≤0, _t^∞γ^τ-tl(x(τ),u(τ))dτ-z≤ 0, with c(x(τ))≤0c(x(τ))≤ 0 for ∀τ≥t∀τ≥ t. Thus, there will exist a joint control u such that V(x,z)≤0V(x,z)≤ 0. Second, when V(x,z)≤0V(x,z)≤ 0 and c(x(τ))≤0c(x(τ))≤ 0 for ∀τ≥t∀τ≥ t hold, it implies that there exists u∈u such that ∫t∞γτ−tl(x(τ),u(τ))τ−z≤0, _t^∞γ^τ-tl(x(τ),u(τ))dτ-z≤ 0, which concludes v(x)−z≤0v(x)-z≤ 0. Therefore, the Lemma 3.1 is proved. ∎ A.2 Lemma 5: Optimality Condition Proof. Following proofs in (Lee, 2022; Zhang et al., 2024; Evans, 2022), given all (t,x,z)∈[0,∞)×ℝ(t,x,z)∈[0,∞)×X×R and select a enough small h>0h>0. There exist two different joint control inputs (u1(⋅),u2(⋅))∈(u_1(·),u_2(·)) such that u(τ)=u1(τ),τ∈[t,t+h],u2(τ),τ∈(t+h,∞).u(τ)= casesu_1(τ),&τ∈[t,t+h],\\ u_2(τ),&τ∈(t+h,∞). cases Then we have the following transformation for Eq. 3 V(x,z)= V(x,z)= minu1∈,u2∈maxmaxτ∈[t,t+h]c(x(τ)),maxτ∈[t+h,∞)c(x(τ)), _u_1 ,u_2 \ _τ∈[t,t+h]c(x(τ)),\, _τ∈[t+h,∞)c(x(τ)), ∫t∞γτ−tl(x(τ),u(τ))dτ−z(t) _t^∞γ^τ-tl(x(τ),u(τ))dτ-z(t) \ = = minu1∈maxmaxτ∈[t,t+h]c(x(τ)),minu2∈maxmaxτ∈[t+h,∞)c(x(τ)),∫t+hγτ−tl(x(τ),u(τ))dτ _u_1 \ _τ∈[t,t+h]c (x(τ) ), _u_2 \ _τ∈[t+h,∞)c (x(τ) ), _t^t+hγ^τ-tl(x(τ),u(τ))dτ +∫t+h∞γτ−tl(x(τ),u(τ))dτ−(z(t+h)+∫t+hγτ−tl(x(τ),u(τ))dτ) + _t+h^∞γ^τ-tl(x(τ),u(τ))dτ- (z(t+h)+ _t^t+hγ^τ-tl(x(τ),u(τ))dτ ) \ \ = = minu1∈maxmaxτ∈[t,t+h]c(x(τ)),minu2∈maxmaxτ∈[t+h,∞)c(x(τ)), _u_1 \ _τ∈[t,t+h]c(x(τ)), _u_2 \ _τ∈[t+h,∞)c (x(τ) ), ∫t+h∞γτ−tl(x(τ),u(τ))dτ−z(t+h) _t+h^∞γ^τ-tl(x(τ),u(τ))\,dτ-z(t+h) \ \ ≈ ≈ minu1∈maxmaxτ∈[t,t+h]c(x(τ)),minu2∈maxmaxτ∈[t+h,∞)c(x(τ)), _u_1 \! \ _τ∈[t,t+h]c(x(τ)), _u_2 \ _τ∈[t+h,∞)c (x(τ) ), γh(∫t+h∞γτ−(t+h)l(x(τ),u(τ))dτ−z(t+h)) γ^h ( _t+h^∞γ^τ-(t+h)l(x(τ),u(τ))\,dτ-z(t+h) ) \ \ = = minu1∈maxmaxτ∈[t,t+h]c(x(τ)),γhV(x(t+h),z(t+h)) _u_1 \ _τ∈[t,t+h]c(x(τ)),\,γ^hV(x(t+h),z(t+h)) \ = = minu∈maxmaxτ∈[t,t+h]c(x(τ)),γhV(x(t+h),z(t+h)) _u \ _τ∈[t,t+h]c(x(τ)),\,γ^hV(x(t+h),z(t+h)) \ Therefore, the Lemma 5 is proved. ∎ A.3 Theorem 6: Epigraph-based HJB PDE Proof. Following proofs in (Lee, 2022; Zhang et al., 2024; Evans, 2022), given all all (t,x,z)∈[0,∞)×ℝ(t,x,z)∈[0,∞)×X×R with a small horizon Δt>0 t>0, we apply Lemma 5 and Taylor expansion to derive the epigraph-based HJB PDE as follows V(x,z) V(x,z) =minu∈maxmaxτ∈[t,t+Δt]c(x(τ)),γhV(x(t+Δt),z(t+Δt)) = _u \ _τ∈[t,t+ t]c(x(τ)),\,γ^hV(x(t+ t),z(t+ t)) \ ≈minu∈maxc(x),(1+lnγΔt)(V(x,z)+∇xV⋅f(x,u)Δt−∂zV⋅l(x,u)Δt+o(Δt)) ≈ _u \c(x),(1+ γ t)(V(x,z)+ _xV· f(x,u) t- _zV· l(x,u) t+o( t)) \ =maxc(x),(1+lnγΔt)minu∈(V(x,z)+∇xV⋅f(x,u)Δt−∂zV⋅l(x,u)Δt+o(Δt)) = \c(x),(1+ γ t) _u (V(x,z)+ _xV· f(x,u) t- _zV· l(x,u) t+o( t)) \ Subtracting V(x,z)V(x,z) from both sides of above equality, dividing by Δt t, and letting Δt→0 t→ 0 yields the following HJB PDE, where V(x,z)V(x,z) is the optimal solution to such PDE. maxc(x)−V(x,z),minu∈[∇xV⋅f(x,u)−∂zV⋅l(x,u)+lnγ⋅V]=0. \! \c(x)-V(x,z), _u [ _xV· f(x,u)- _zV· l(x,u)+ γ· V ] \=0. Here ℋ=∇xV⋅f(x,u)−∂zV⋅l(x,u)+lnγ⋅VH= _xV· f(x,u)- _zV· l(x,u)+ γ· V is Hamiltonian and optimal control u∗=argminu∈ℋu^*= *arg\,min_u H. Next we prove that V(x,z)V(x,z) is the unique viscosity solution to the epigraph-based HJB PDE using the contradiction technique. First, for U∈C∞(×ℝ)U∈ C^∞(X×R) such that V−UV-U has local maximum at (x0,z0)∈×ℝ(x_0,z_0) ×R and (V−U)(x0,z0)=0(V-U)(x_0,z_0)=0, we will prove maxc(x0)−U(x0,z0),minu∈[∇xU(x0,z0)⋅f(x0,u)−∂zU(x0,z0)⋅l(x0,u)+lnγ⋅U(x0,z0)]≥0. \c(x_0)-U(x_0,z_0),\; _u [ _xU(x_0,z_0)· f(x_0,u)- _zU(x_0,z_0)· l(x_0,u)+ γ· U(x_0,z_0) ] \≥ 0. Suppose the above inequality is not correct. We consider that there exists θ>0θ>0 and u~∈ u such that c(x)−U(x0,z0)≤−θ, c(x)-U(x_0,z_0)≤-θ, ∇xU⋅f(x,u~)−∂zU⋅l(x,u~)+lnγ⋅U≤−θ. _xU· f(x, u)- _zU· l(x, u)+ γ· U≤-θ. for all points (x,z)(x,z) sufficiently close to (x0,z0)(x_0,z_0): ‖x(s)−x0‖+|z(s)−z0|<h\|x(s)-x_0\|+|z(s)-z_0|<h for small enough h>0h>0, where s∈[t0,t0+h]s∈[t_0,t_0+h]. Under the assumptions in Sec. 3.1.1, and given state trajectories x and z evolved from the initial conditions x=x0x=x_0 and z=z0z=z_0 according to the corresponding dynamics, the following inequality holds c(x(s))−U(x0,z0)≤−θ, c(x(s))-U(x_0,z_0)≤-θ, ∇xU(x(s),z(s))⋅f(x(s),u~)−∂zU(x(s),z(s))⋅l(x(s),u~)+lnγ⋅U(x(s),z(s))≤−θ. _xU(x(s),z(s))· f(x(s), u)- _zU(x(s),z(s))· l(x(s), u)+ γ· U(x(s),z(s))≤-θ. Since V−UV-U has a local maximum at (x0,z0)(x_0,z_0), we can have that minu∈[γhV(x(t0+h),z(t0+h))−V(x0,z0)] _u [γ^hV(x(t_0+h),z(t_0+h))-V(x_0,z_0) ] ≤ ≤ minu∈[γhU(x(t0+h),z(t0+h))−U(x0,z0)] _u [γ^hU(x(t_0+h),z(t_0+h))-U(x_0,z_0) ] = = minu∈[(∇xU(x(t0),z(t0))⋅f(x(t0),u)−∂zU(x(t0),z(t0))⋅l(x(t0),u)+lnγ⋅U(x(t0),z(t0)))h] _u [( _xU(x(t_0),z(t_0))· f(x(t_0),u)- _zU(x(t_0),z(t_0))· l(x(t_0),u)+ γ· U(x(t_0),z(t_0)))h ] ≤ ≤ −θh -θ h We know that Lemma 2 implies V(x0,z0)=minu∈maxmaxsτ∈[t0,t0+h]c(x(s)),γhV(x(t0+h),z(t0+h)).V(x_0,z_0)= _u \ _sτ∈[t_0,t_0+h]c(x(s)),\;γ^hV(x(t_0+h),z(t_0+h)) \. By subtracting U(x0,z0)U(x_0,z_0) on both side, we have (V−U)(x0,z0)=minu∈maxc(x(s))−U(x0,z0),γhV(x(t0+h),z(t0+h))−U(x0,z0).(V-U)(x_0,z_0)= _u \c(x(s))-U(x_0,z_0),\;γ^hV(x(t_0+h),z(t_0+h))-U(x_0,z_0) \. Since (V−U)(x0,z0)=0(V-U)(x_0,z_0)=0 holds such that V(x0,z0)=U(x0,z0)V(x_0,z_0)=U(x_0,z_0), then we will have that minu∈maxc(x(s))−V(x0,z0),γhV(x(t0+h),z(t0+h))−V(x0,z0)=minu∈maxθ,θh>0, _u \c(x(s))-V(x_0,z_0),\;γ^hV(x(t_0+h),z(t_0+h))-V(x_0,z_0) \= _u \θ,θ h\>0, which has a contradiction with (V−U)(x0,z0)=0(V-U)(x_0,z_0)=0. Thus we prove that maxc(x0)−U(x0,z0),minu∈[∇xU(x0,z0)⋅f(x0,u0)−∂zU(x0,z0)⋅l(x0,u0)+lnγ⋅U(x0,z0)]≥0. \c(x_0)-U(x_0,z_0),\; _u [ _xU(x_0,z_0)· f(x_0,u_0)- _zU(x_0,z_0)· l(x_0,u_0)+ γ· U(x_0,z_0) ] \≥ 0. Second, for U∈C∞(×ℝ)U∈ C^∞(X×R) such that V−UV-U has local minimum at (x0,z0)∈×ℝ(x_0,z_0) ×R and (V−U)(x0,z0)=0(V-U)(x_0,z_0)=0, we will prove maxc(x0)−U(x0,z0),minu∈[∇xU(x0,z0)⋅f(x0,u0)−∂zU(x0,z0)⋅l(x0,u0)+lnγ⋅U(x0,z0)]≤0. \c(x_0)-U(x_0,z_0),\; _u [ _xU(x_0,z_0)· f(x_0,u_0)- _zU(x_0,z_0)· l(x_0,u_0)+ γ· U(x_0,z_0) ] \≤ 0. The definition of auxiliary value V(x,z)V(x,z) shows that V(x,z) V(x,z) =minu∈maxmaxτ∈[t,∞]c(x(τ)),∫t∞γτ−tl(x(τ),u(τ))τ−z = _u \ _τ∈[t,∞]c(x(τ)), _t^∞γ^τ-tl(x(τ),u(τ))dτ-z \ ≥minu∈maxc(x0),∫t∞γτ−tl(x(τ),u(τ))τ−z ≥ _u \c(x_0), _t^∞γ^τ-tl(x(τ),u(τ))dτ-z \ for all u∈u . By subtracting U(x0,z0)U(x_0,z_0) on both sides, we have 0=(V−U)(x0,z0)≥maxc(x0)−U(x0,z0),∫t∞γτ−tl(x,u)τ−z0−U(x0,z0).0=(V-U)(x_0,z_0)≥ \c(x_0)-U(x_0,z_0), _t^∞γ^τ-tl(x,u)dτ-z_0-U(x_0,z_0)\. The rest of the proof is to show minu∈[∇xU(x0,z0)⋅f(x0,u)−∂zU(x0,z0)⋅l(x0,u)+lnγ⋅U(x0,z0)]≤0. _u [ _xU(x_0,z_0)· f(x_0,u)- _zU(x_0,z_0)· l(x_0,u)+ γ· U(x_0,z_0) ]≤ 0. Suppose the above inequality is not correct. We consider that there exists θ>0θ>0 such that minu∈[∇xU(x,z)⋅f(x,u)−∂zU(x,z)⋅l(x,u)+lnγ⋅U(x,z)]≥θ, _u [ _xU(x,z)· f(x,u)- _zU(x,z)· l(x,u)+ γ· U(x,z) ]≥θ, for all points (x,z)(x,z) sufficiently close to (x0,z0)(x_0,z_0): ‖x−x0‖+|z−z0|<h\|x-x_0\|+|z-z_0|<h for small enough h>0h>0, where s∈[t0,t0+h]s∈[t_0,t_0+h]. Given state trajectories x and z that evolve from the initial conditions x=x0x=x_0 and z=z0z=z_0 under the corresponding dynamics with any control u~∈ u , where u~(s)= u(s)= argminu~∈∇xU(x(s),z(s))⋅f(x(s),u~)−∂zU(x(s),z(s))⋅l(x(s),u~) _ u \ _xU(x(s),z(s))· f(x(s), u)- _zU(x(s),z(s))· l(x(s), u) +lnγ⋅U(x(s),z(s)). + γ· U(x(s),z(s)) \. Then we have the following condition that holds ∇xU(x(s),z(s))⋅f(x(s),u~)−∂zU(x(s),z(s))⋅l(x(s),u~)+lnγ⋅U(x(s),z(s))≥θ. _xU(x(s),z(s))· f(x(s), u)- _zU(x(s),z(s))· l(x(s), u)+ γ· U(x(s),z(s))≥θ. Consider V−UV-U has a local minimum at (x0,z0)(x_0,z_0), we will have that minu~∈[γhV(x(t0+h),z(t0+h))−V(x0,z0)] _ u [γ^hV(x(t_0+h),z(t_0+h))-V(x_0,z_0) ] ≥ ≥ minu~∈[γhU(x(t0+h),z(t0+h))−U(x0,z0)] _ u [γ^hU(x(t_0+h),z(t_0+h))-U(x_0,z_0) ] = = minu~∈[(∇xU(x(t0),z(t0))⋅f(x(t0),u~)−∂zU(x(t0),z(t0))⋅l(x(t0),u~)+lnγ⋅U(x(t0),z(t0)))h] _ u [( _xU(x(t_0),z(t_0))· f(x(t_0), u)- _zU(x(t_0),z(t_0))· l(x(t_0), u)+ γ· U(x(t_0),z(t_0)))h ] ≥ ≥ θh θ h Based on this derivation, we finally have that minu~∈γhV(x(t0+h),z(t0+h))≥V(x0,z0)+θh>V(x0,z0). _ u γ^hV(x(t_0+h),z(t_0+h))≥ V(x_0,z_0)+θ h>V(x_0,z_0). However, we know that Lemma 5 implies that minu~∈γhV(x(t0+h),z(t0+h))≤V(x0,z0), _ u γ^hV(x(t_0+h),z(t_0+h))≤ V(x_0,z_0), which is a contradiction. Thus, we prove that maxc(x0)−U(x0,z0),minu∈[∇xU(x0,z0)⋅f(x0,u0)−∂zU(x0,z0)⋅l(x0,u0)+lnγ⋅U(x0,z0)]≤0. \c(x_0)-U(x_0,z_0),\; _u [ _xU(x_0,z_0)· f(x_0,u_0)- _zU(x_0,z_0)· l(x_0,u_0)+ γ· U(x_0,z_0) ] \≤ 0. Hence, we prove that V(x,z)V(x,z) is the viscosity solution to the epigraph-based HJB PDE. The uniqueness follows Theorem 1 of Chapter 10 in Evans (2022). ∎ A.4 Advantage Function We define the Q(xt,zt,ut)=maxc(xt),rhV(xt+h,zt+hQ(x_t,z_t,u_t)= \c(x_t),r^hV(x_t+h,z_t+h\ over a short time interval h>0h>0 and compute Q(xt,zt,ut)−V(xt,zt)= Q(x_t,z_t,u_t)-V(x_t,z_t)= maxc(xt),rhV(xt+h,zt+h)−V(xt,zt) \c(x_t),r^hV(x_t+h,z_t+h)\-V(x_t,z_t) = = maxc(xt)−V(xt,zt),(1+lnγh)(V(xt,zt)+∇xtV⋅f(xt,ut)h \c(x_t)-V(x_t,z_t),(1+ γ h)(V(x_t,z_t)+ _x_tV· f(x_t,u_t)h −∂ztV⋅l(xt,ut)h−V(xt,zt)+o(h) - _z_tV· l(x_t,u_t)h-V(x_t,z_t)+o(h)\ = = maxc(xt)−V(xt,zt),(∇xtV⋅f(xt,ut)−∂ztV⋅l(xt,ut)+lnγ⋅V)h \c(x_t)-V(x_t,z_t),( _x_tV· f(x_t,u_t)- _z_tV· l(x_t,u_t)+ γ· V)h\ We divide h on both sides of the above equation and let h→0h→ 0 to compute the advantage function as A(xt,zt,ut) A(x_t,z_t,u_t) =limh→0Q(xt,zt,ut)−V(xt,zt)h = _h→ 0 Q(x_t,z_t,u_t)-V(x_t,z_t)h =maxc(xt)−V(xt,zt),∇xtV⋅f(xt,ut)−∂ztV⋅l(xt,ut)+lnγ⋅V = \c(x_t)-V(x_t,z_t), _x_tV· f(x_t,u_t)- _z_tV· l(x_t,u_t)+ γ· V\ A.5 Convergence of Epigraph Value Function Consider the augmented state (x,z)(x,z) with state constraint c(x)c(x) and non-negative cost l(x,u)l(x,u). Define the discounted epigraph-Bellman operator over a short step Δt>0 t>0: (V)(xt,zt):=(1−γΔt)c(xt)+γΔtminu∈maxc(xt),V(xt+Δt,zt+Δt),(TV)(x_t,z_t):=(1-γ t)c(x_t)+γ t _u \ \c(x_t),\,V(x_t+ t,z_t+ t) \ \, for V:×ℝ→ℝV:X×R bounded. Then the value iteration Vk+1=VkV_k+1=TV_k converges uniformly to the unique fixed point of T. Proof. (i) Contraction. For any c(xt)c(x_t) and bounded functions V,WV,W, we have the following condition satisfying the contraction. |maxc(xt),V(xt+Δt,zt+Δt)−maxc(xt),W(xt+Δt,zt+Δt)| | \c(x_t),V(x_t+ t,z_t+ t)\- \c(x_t),W(x_t+ t,z_t+ t)\| ≤ ≤ |V(xt+Δt,zt+Δt)−W(xt+Δt,zt+Δt)| |V(x_t+ t,z_t+ t)-W(x_t+ t,z_t+ t)| ≤ ≤ ‖V−W‖∞ \|V-W\|_∞ (i) Existence and uniqueness. By Banach’s fixed-point theorem, T admits a unique fixed point V, and for value iteration Vk+1=VkV_k+1=TV_k we have that ‖Vk−V‖∞≤γk‖V0−V‖∞→0,\|V_k-V\|_∞≤γ^k\|V_0-V\|_∞→ 0, (i) Approximate evaluation. If each iteration uses an approximate operator ~ T satisfying ‖~V−V‖∞≤ε\| TV-TV\|_∞≤ , then lim supk→∞‖Vk−V‖∞≤ε1−γΔt. _k→∞\|V_k-V\|_∞≤ 1-γ t. ∎ Appendix B Training Algorithms In this part, we provide additional details on the overall algorithmic pipeline and clarify the key implementation choices. Algorithm 1 Epigraph-Based Continuous-Time MARL 1:Initialize actor πθ _θ, return critic VψretV^ret_ψ, constraint critic VϕconsV^cons_φ, dynamics network fξf_ξ, reward network lφl_ , and local rollout ℛR. 2:for l=1,…,Tl=1,…,T do 3: ⊳ Collect one rollout: 4: x←env.reset()x .reset() 5: for k=1,…,Kk=1,…,K do 6: sample arbitrary decision time t∼t 7: for each agent i=1,…,Ni=1,…,N do 8: ui∼πθi(ui∣x)u_i _ _i(u_i x) 9: end for 10: set joint action u=(u1,…,uN)u=(u_1,…,u_N) 11: (x′,r)←env.step(u)(x ,r) .step(u) 12: append (x,u,r,x′)(x,u,r,x ) to local rollout ℛR 13: x←x′x← x 14: end for 15: ⊳ Outer optimization: epigraph update 16: find z∗=infz∈ℝ:maxVϕcons(x),Vψret(x,z)−z≤0z = \z : \V^cons_φ(x),\,V^ret_ψ(x,z)-z\≤ 0\ 17: ⊳ Dynamics and Cost Model learning on ℛR 18: update ξ,φξ, as per the Eq. 17. 19: ⊳ Inner optimization given z∗z : Critic update on RX_R 20: update ψ,ϕψ,φ by losses ℒcons,ℒret,ℒHJBL_cons,L_ret,L_HJB and ℒVGIL_VGI as per the Eq. 11, Eq. 10 and Eq. 12. 21: ⊳ Actor update for each agent 22: for i=1,…,Ni=1,…,N do 23: compute A(x,u,z∗)A(x,u,z ) for all (x,u,z∗)∈R(x,u,z ) _R and update the θ as the Eq. 18. 24: end for 25:end for Appendix C Environmental Settings We provide detailed descriptions of all benchmark environments used in our experiments. For each scenario, we list the number of agents, the number of obstacles, the safety constraints imposed, and the specific task objective with metrics. Metrics. We report two primary metrics—one reward-style training score that aggregates task cost and constraint penalty, and one violation rate measured over held-out rollouts. (1) Cumulative penalty / reward-style training score. In many standard environments (e.g., MPE and multi-agent MuJoCo), the task reward often consists of two independent components: (i) a task term such as distance-to-target or velocity tracking, and (i) a safety penalty that is activated only when constraint-relevant events occur (e.g., collisions or proximity violations). This design is also used in prior safe MARL methods such as MACPO and Lagrangian baselines (Gu et al., 2021). For clarity of notation, we write the task cost as ℓt≥0 _t≥ 0 (derived from the negative reward of the task term) and denote the constraint penalty as κt≥0 _t≥ 0. The environment therefore provides a composite instantaneous cost ψt:=ℓt+κt, _t\;:=\; _t+ _t, which simply aggregates the task objective and the constraint penalty already defined in the environment. For a trajectory τ with horizon T(τ)T(τ), we define the total episode cost as J(τ):=∑t=0T(τ)−1ψt,S(τ):=−J(τ),J(τ):= _t=0^T(τ)-1 _t, S(τ):=-J(τ), where S(τ)S(τ) is the cumulative reward used for performance plots. (2) Violation rate (evaluation). Given NevalN_eval episodes (we use Neval=100N_eval=100 by default), define the episode-level violation indicator v(τ):= 1∃ts.t.κt>0,v(τ)\;:=\;1\! \\,∃\,t\ s.t.\ _t>0\ \, i.e., an episode is counted as violating if it ever incurs a positive state-constraint penalty.111If κt _t is an indicator of hard violations, this coincides with “any violation.” If κt _t is a continuous hinge, we use the same criterion κt>0 _t>0. The violation rate is then Viol. Rate=1Neval∑i=1Nevalv(τi).Viol.\ Rate\;=\; 1N_eval _i=1^N_evalv( _i). C.1 Safe MPE In the MPE, we setup the details as follows: Action. Continuous 2-D acceleration for x and y axis. Reward and costs. Each agent is assigned a per-agent target gig_i. The dense goal reward is rigoal(t)=−‖xi(t)−gi‖2.r^goal_i(t)=-\|x_i(t)-g_i\|_2. A discrete collision cost with obstacles or other agents applies: cidisc(t)=10,if agent–obstacle overlap0,otherwise.c^disc_i(t)= cases10,&if agent--obstacle overlap\\ 0,&otherwise. cases We also record a continuous proximity/penetration cost (not added into the dense goal reward): cicont(t)=12∑o∈ϕ((ri+ro)−‖xi−xo‖),ϕ(δ)=20δ,δ>0(overlap)0.5δ,δ≤0c^cont_i(t)\;=\; 12\!\! _o \!\!φ ((r_i+r_o)-\|x_i-x_o\| ), φ(δ)= cases20\,δ,&δ>0\ \ (overlap)\\ 0.5\,δ,&δ≤ 0 cases where ri,ror_i,r_o are radius (sizes). Difference from the original discrete-time MPE. The standard MPE environment uses a fixed and discrete integration step Δt t, where each simulation step updates the agent states according to pt+1=pt+vtΔtp_t+1=p_t+v_t t and vt+1=vt+ftΔtv_t+1=v_t+f_t t with a fixed time increment. In contrast, our continuous-time MPE adapts the physical integration step to an arbitrary Δt t provided by the learning algorithm. The state evolution follows p˙(t)=v(t),v˙(t)=f(t)m−damping⋅v(t), p(t)=v(t), v(t)= f(t)m-damping· v(t), and is numerically integrated via p←p+v⋅Δt,v←v+fmΔt,p← p+v· t, v← v+ fm t, using the user-specified Δt t. For clarity, the update used in the original environment is: step (F): (F): p p =p+v⋅0.1(fixed as 0.1), =p+v· 1 (fixed as 1), v v =v+Fm⋅0.1(fixed as 0.1). =v+ Fm· 1 (fixed as 1). Our continuous-time version introduces: step_continuous (F,Δt): (F, t): p p =p+v⋅Δt(depend on the input Δt), =p+v· t (depend on the input t), v v =v+Fm⋅Δt(depend on the input Δt). =v+ Fm· t (depend on the input t). so that the state update depends directly on the argument Δt t rather than a fixed constant. Corridor. This scenario contains 3 agents with 2 large corridor walls. Agents must avoid collisions with the corridor walls and with each other while navigating from their starting positions to reach the assigned target locations on the opposite side. Formation. This scenario also involves 3 agents and 2 obstacles. The agents are required to bypass obstacles and then coordinate to form a triangular formation at the designated region, under the constraint of avoiding collisions with both obstacles and other agents. Line. In this task, 3 agents operate in an environment with 2 obstacles. After avoiding the obstacles, the agents must position themselves to form a straight line. The safety constraints enforce that no agent collides with obstacles or with other agents during navigation. Target. This scenario uses 2 agents with 1 obstacle placed in the environment. Each agent is assigned a fixed target position, and the agents must navigate to their respective goals while avoiding collisions with the obstacle and with each other. Cooperative Navigation. This is a cooperative navigation task with 3 agents and no obstacles. The agents must spread out to cover multiple target landmarks while avoiding collisions among themselves. Specifically, the agents’ goals are the one closest to them rather than fixed ones. Cooperative Predator–Prey. This task includes 3 controllable predator agents and 1 prey that moves randomly. There are no obstacles, but predators must avoid colliding with each other. The predators’ objective is to coordinate their movements to capture the prey. C.2 Safe Multi-agent MuJoCo Half Cheetah. We adapt the Half Cheetah environment into three multi-agent variants: Half Cheetah-2x3, Half Cheetah-3x2, and Half Cheetah-6x1. In each case, the body is partitioned into joints agents with different grouping configurations. For example, Half Cheetah-3x2 is three agents with 2 moving joints for each agent. Randomly placed walls are introduced into the environment, requiring the agents not only to coordinate efficient forward locomotion but also to avoid collisions with obstacles. Reward. r=rrun=xt+1−xtΔt.r=r_run= x_t+1-x_t t. Safety cost. A binary proximity cost to the wall: ct=|xwall−xagent|<9∈0,1.c_t=1\\,|x_wall-x_agent|<9\,\∈\0,1\. Observation augments the usual state with wall velocity, wall force proxy, and clipped distance to the wall; the environment also recolors the wall when unsafe. Difference from the original MuJoCo environment. In standard MuJoCo control tasks, the simulation uses a fixed micro time step 0.010.01 (each frame takes 0.010.01), and each environment step corresponds to a fixed number of internal physics frames (e.g., frame_skip =5=5), resulting in a fixed control interval Δt=0.05 t=0.05. Our continuous-time MuJoCo variant removes this fixed control interval. For any desired Δt t, we execute do_simulation(a,N),N=Δt0.01, do\_simulation(a,\;N), N= t0.01, i.e., the number of internal physics frames is chosen dynamically according to the requested integration step. Thus the effective control interval is no longer fixed but fully determined by Δt t, enabling variable-resolution continuous-time rollouts. The reward terms (forward velocity, control cost, contact cost) are normalized by the actual Δt t, ensuring consistency across different temporal resolutions. The original update is: step(u):N=5,do_simulate(u,N). step(u): N=5, do\_simulate(u,N). Our continuous-time version becomes: step_continuous(u,Δt):N=Δt/0.01,do_simulate(u,N). step\_continuous(u, t): N= t/0.01, do\_simulate(u,N). Ant. We construct four multi-agent variants of the Ant: Ant-2x4, Ant-4x2, Ant-8x1, and Ant-2x4d. In all cases, the body is controlled by joints agents arranged in different groupings across the legs. As with Half Cheetah, walls are introduced as obstacles, and the agents must coordinate locomotion while ensuring safety by avoiding collisions with these obstacles. The reward is set same as the Half Cheetah Safety shaping. Identical piecewise-slant corridor: compute yoffy_off from (x,y)(x,y) and define ctobj=|yoff|<1.8.c^obj_t=1\|y_off|<1.8\. C.3 Constrained Coupled Oscillator Environment We consider a two–agent coupled spring–damper system. The state and control are x=[x1v1x2v2]⊤,u=[u1u2]⊤.x\!=\! bmatrixx_1&v_1&x_2&v_2 bmatrix , u\!=\! bmatrixu_1&u_2 bmatrix . Each agent i∈1,2i∈\1,2\ controls one mass with continuous–time dynamics x˙i x_i =vi, =v_i, v˙i v_i =−kxi−bvi+ui, =-k\,x_i-b\,v_i+u_i, with spring constant k=1.0k=1.0 and damping coefficient b=0.5b=0.5. Stacking the states gives x˙=Ax+Bu x=Ax+Bu with A=[0100−k−b00000100−k−b],B=[00100001].A= bmatrix0&1&0&0\\ -k&-b&0&0\\ 0&0&0&1\\ 0&0&-k&-b bmatrix, B= bmatrix0&0\\ 1&0\\ 0&0\\ 0&1 bmatrix. Control limits and discretization. Actions are normalized u~i∈[−1,1] u_i∈[-1,1] and mapped to physical inputs by ui=umaxu~iu_i=u_ u_i with umax=10u_ =10 (component–wise box constraint). vit+1 v_i^t+1 =vit+(−kxit−bvit+uit)Δt, =v_i^t+ (-k\,x_i^t-b\,v_i^t+u_i^t ) t, xit+1 x_i^t+1 =xit+vit+1Δt, =x_i^t+v_i^t+1 t, for a horizon of N=30N=30 steps. Stage cost. The per–step quadratic cost is ℓ(x,u)=x12+x22+λc(x1−x2)2+β(u12+u22), (x,u)=x_1^2+x_2^2+ _c\,(x_1-x_2)^2+β\,(u_1^2+u_2^2), with coupling strength λc=2.0 _c=2.0 and control penalty β=0.01β=0.01. Equivalently, ℓ(x,u)=x⊤Qx+u⊤Ru (x,u)=x Qx+u Ru where Q=[1+λc0−λc00000−λc01+λc00000],R=βI2.Q= bmatrix1+ _c&0&- _c&0\\ 0&0&0&0\\ - _c&0&1+ _c&0\\ 0&0&0&0 bmatrix, R=β I_2. For training we use a shaped reward rt=−130ℓ(xt,ut).r_t=- 130\, (x_t,u_t). Hard state constraint. We impose an ordering constraint between the two positions, x1≤x2+0.02,x_1\;≤\;x_2+0.02, and record an additional penalty pt=−10⋅x1,t>x2,t+0.02,p_t=-10·1\! \\,x_1,t>x_2,t+0.02\, \, returned alongside rtr_t. Smooth violation signal (for logging). We also log a smooth surrogate of the constraint violation, ϕ(x)= 2σ(s(x1−x2+0.02))−1,σ(z)=11+e−z,s=20,φ(x)\;=\;2\,σ\! (s\,(x_1-x_2+0.02) )-1, σ(z)= 11+e^-z,\;s=20, which maps to (−1,1)(-1,1) and grows monotonically with the amount of violation. Unconstrained LQR. The continuous-time algebraic Riccati equation (CARE) A⊤P+PA−PBR−1B⊤P+Q=0A P+PA-PBR^-1B P+Q=0 is solved for the unique positive semidefinite matrix P. The unconstrained optimal linear feedback is K=R−1B⊤P,uLQR(x)=−Kx.K\;=\;R^-1B P, u_LQR(x)\;=\;-Kx. Hard state constraint and CBF condition. We impose the safety constraint x1−x2−0.02≤ 0⟺h(x):= 0.02−(x1−x2)≥ 0.x_1-x_2-0.02\;≤\;0 h(x)\;:=\;0.02-(x_1-x_2)\;≥\;0. Let ∇h(x)=[−1010]⊤∇ h(x)= bmatrix-1&0&1&0 bmatrix^\! . A (first-order) control barrier function (CBF) condition enforces forward invariance of the safe set =x:h(x)≥0C=\x:h(x)≥ 0\ by requiring h˙(x,u)=∇h(x)⊤(Ax+Bu)≥−αh(x), h(x,u)\;=\;∇ h(x) (Ax+Bu)\;≥\;-α\,h(x), with a user-chosen class-K parameter α>0α>0. Defining a(x):=∇h(x)⊤B∈ℝ2,b(x):=−∇h(x)⊤Ax−αh(x)∈ℝ,a(x)\;:=\;∇ h(x) B ^2, b(x)\;:=\;-\,∇ h(x) Ax\;-\;α\,h(x) , the CBF condition Eq. C.3 is the single affine-in-u half-space constraint a(x)⊤u≥b(x).a(x) u\;≥\;b(x). Closed-form safety projection. To obtain a safe control with minimal distortion from uLQRu_LQR, we solve the weighted projection minu∈ℝ212(u−uLQR)⊤W(u−uLQR)s.t.a(x)⊤u≥b(x), _u ^2\; 12\,(u-u_LQR) W\,(u-u_LQR) .t. a(x) u≥ b(x), with W=RW=R (“R-metric”; Euclidean W=IW=I is also possible). Because Eq. C.3 has a single linear constraint, it admits a closed form: u⋆(x)=uLQR(x),if a⊤uLQR≥b,uLQR(x)+τW−1a,otherwise, with τ=b−a⊤uLQRa⊤W−1a.u (x)\;=\; casesu_LQR(x),&if a u_LQR≥ b,\\[2.0pt] u_LQR(x)+τ\,W^-1a,&otherwise, with τ\;=\; b-a u_LQR\,a W^-1a\,. cases Finally we saturate to the actuator limits umax>0u_ >0: uGT(x)=clip(u⋆(x),−umax,umax).u_GT(x)\;=\;clip\! (u (x),\,-u_ ,\,u_ ). Appendix D Additional Environmental Results D.1 Visiual Trajectories (a) EPI-1 (b) EPI-2 (c) EPI-3 (d) EPI-4 (e) EPI-5 (f) EPPO-1 (g) EPPO-2 (h) EPPO-3 (i) EPPO-4 (j) EPPO-5 (k) MACPO-1 (l) MACPO-2 (m) MACPO-3 (n) MACPO-4 (o) MACPO-5 Figure 10: Trajectory demonstrations (key frames) across methods in Formation. Row 1: EPI results, Row 2: EPPO results, Row 3: MACPO results. The trajectory demonstrations in Fig. 10 highlight clear behavioral differences across algorithms in Formation scenario. Our proposed method EPI learns smooth trajectories that avoid obstacles while consistently reaching the target, demonstrating both constraint satisfaction and goal achievement. In contrast, EPPO occasionally captures the avoidance behavior but often gets stuck at suboptimal solutions. This is because during training, its randomized sampling of the auxiliary state z prevents stable policy convergence in continuous-time settings; even if outer optimization is applied at execution, the learned policy lacks accurate control signals. On the other hand, MACPO, which enforces hard constraints via a trust-region style update, tends to overestimate the obstacle region. As a result, agents often exhibit overly conservative behaviors—such as retreating toward corners to avoid violations—rather than efficiently pursuing their targets. Together, these comparisons confirm that EPI achieves the most balanced and effective behavior among the three approaches. D.2 Performance Under Stochastic Settings Figure 11: Performance under Different Noise Levels. To evaluate robustness under stochastic dynamics, we perturb the continuous-time transition model as xt+Δt=f(xt,ut)Δt+εt,εt∼(0,σ2I).x_t+ t=f(x_t,u_t)\, t+ _t, _t (0,σ^2I). in Fig.11. We consider three noise magnitudes: Low Noise: σ2=0.1σ^2=0.1 Mid Noise: σ2=0.5σ^2=0.5 and High Noise: σ2=1.0σ^2=1.0. We observe that No Noise and Low Noise yield similar identical cost and constraint-violation behavior across all three tasks. Because the PINN-based value approximation are inherently robust to small local perturbations, as long as the injected disturbance is within a moderate range, the learned dynamics model, cost model, and value gradients remain accurate. In contrast, Mid Noise and High Noise introduce much larger deviations in the state propagation. These disturbances accumulate over time, causing the PINN to receive significantly deviated training signals. Since our method does not incorporate explicit uncertainty modeling or stochastic HJB formulations, the serious noise directly degrades the learned critic and value gradients, eventually leading to unstable or even failed policies. Figure 12: Average Distance to the Target under Different Δt t. D.3 Effect of the Discretization Interval. Figure 12 evaluates how the choice of discretization interval Δt t affects the performance of EPI. For each fixed Δt t, we roll out complete trajectories using the learned policy and measure the average distance to the target over the entire trajectory. Across all three scenarios, we observe a consistent trend: the average distance to the target increases as Δt t becomes larger. This behavior is expected in continuous-time control. When Δt t is small, the temporal resolution is high and the policy is updated frequently, allowing the learned value gradients to provide fine-grained control corrections. In contrast, larger Δt t leads to coarser control updates, reducing the precision of the policy’s response to the evolving system dynamics. Moreover, both the HJB residual and the VGI update rely on local differential information. As Δt t grows, the mismatch between the continuous-time formulation and the discrete rollout increases, which in turn amplifies approximation errors in the learned value gradients. These errors accumulate along the trajectory and result in the observed degradation in task accuracy. D.4 Trajectory of z∗z^* through the Training Figure 13: z∗z^* Trajectory through the Training in Target. Figure 13 illustrates the evolution of the optimal epigraph variable zt∗z_t and the active branch (return vs. constraint) at three representative stages of training. In early training (Episode 1000), the policy frequently visits infeasible states, causing Vcons(xt)>0V_cons(x_t)>0 and forcing the epigraph to select the constraint branch; consequently zt∗z_t remains at the clipped upper bound zmaxz_ . By mid training (Episode 12000), the critic starts to maintain Vcons(xt)≤0V_cons(x_t)≤ 0 for part of the trajectory, producing intermittent switching and a decreasing zt∗z_t . In late training (Episode 23000), the trajectory remains feasible, the return branch is consistently selected, and zt∗z_t decreases smoothly along the rollout. These behaviors align with the expected epigraph semantics: infeasible states produce zmaxz_ , while improved policies yield stable return-dominated gradually decreasing zt∗z_t . D.5 Compare EPI with Traditional Epigraph Method Figure 14: Performance of EPI and Traditional Epigraph under MPE settings. Figure 14 compares our z-independent epigraph formulation (EPI) with the traditional z-dependent epigraph used in EPPO-like methods on the Formation, Line, and Target tasks. In the traditional design, a scalar z is randomly sampled at the initial state of each episode and then propagated through its auxiliary dynamics, so that both critic and actor are conditioned on this randomly chosen epigraph level. As shown in Fig. 14, converges to a lower cumulative reward, and exhibits substantially larger variance across seeds. In contrast, EPI learns z-independent critics (Vcons(x),Vret(x)) (V^cons(x),V^ret(x) ) and computes z∗z via a one-dimensional search during training, while the actor depends only on the physical state x. This removes the nonstationary noise introduced by random z sampling: for a fixed x, the policy gradient under EPI is unique, whereas in the traditional epigraph it fluctuates with the sampled z even when the critic has already converged. In continuous-time settings this issue is amplified, since small changes in z shift the switching time between the constraint and return branches and thereby alter the entire rollout. D.6 Comparison between EPI and Discrete-time Baselines Figure 15: Performance of EPI and Discrete-time Baselines under MPE settings. To validate the performance of traditional discrete-time based methods in continuous-time settings, the Fig 15 compares EPI with two discrete-time baselines (DT-CBF and DT-MAPPO-LAG) on the Line and Target tasks in the continuous-time MPE environment. All baselines are adapted to the discrete-time setting by removing their residual-loss components. Apart from this modification, all implementation details follow their original published versions (Zhang et al., 2025a). Across both tasks, EPI consistently achieves lower mean distance to the target and smaller variance, demonstrating the performance gain from the modules that designed for the continuous-time settings. Table 1: Hyperparameter settings used. Parameter Value Episode length for MPE 50 Episode length for MuJoCo 100 Episode length for Didactic 50 Total number of episode for MPE 30000 Total number of episode for MuJoCo 30000 Total number of episode for Didactic 3000 z range for MPE 0-10 z range for MuJoCo 0-5 z range for Didactic 0-2 Discount factor γ 0.99 Actor learning rate 0.0001 Critic (Return) learning rate 0.001 Critic (Constraint)learning rate 0.001 Dynamics model learning rate 0.001 Reward model learning rate 0.001 Exploration steps 1000 Model save interval 1000 Random seed 113-120 Appendix E Hyperparameters and Neural Network Structures Experiments were conducted on hardware comprising an Intel(R) Xeon(R) Gold 6254 CPU @ 3.10GHz, four NVIDIA A5000 GPUs and eight NVIDIA A6000 GPUs. This setup ensures the computational efficiency and precision required for the demanding simulations involved in multi-agent reinforcement learning and safety evaluations. Table LABEL:tab:hp lists the defaults used in all experiments. Episode lengths are chosen so that a single rollout covers a full interaction cycle (50 steps for MPE and the didactic environment, 100 for MuJoCo). We train for 30000 episodes in MPE and MuJoCo and for 3000 episodes in the didactic setting, reflecting simulator cost and convergence speed. The z range controls epigraph sampling for the VGI updates and is set wider in MPE (0–10), moderate in MuJoCo (0–5), and narrow in the didactic task (0–2). The actor uses a conservative learning rate (1e-4) for stable policy updates; the critics and the dynamics/reward models use 1e-3 to accelerate value/model fitting. Training is warm-started with 1000 exploration steps, checkpoints are saved every 1000 episodes, and reported results are averaged over seeds 113–120. Table 2: Summary of neural network architectures used in our framework. Network Input Dimension Architecture and Activation Return Value Network State (d) FC(128) → FC(128) → FC(1), ReLU or Tanh Constraint Value Network State (d) FC(128) → FC(128) → FC(1), ReLU or Tanh Dynamics Network State + Joint Action (d+nad+na) FC(128) → FC(128) → FC(d), ReLU Reward Network State + Joint Action (d+nad+na) FC(128) → FC(128) → FC(1), ReLU PolicyNet Observation + Time Interval (o+1o+1) FC(128) → FC(128) → FC(64) → FC(a), ReLU Table 2 summarizes the five multilayer perceptrons used in our framework. Two scalar critics map the state x∈ℝdx ^d to the return value and the constraint value, each with two hidden layers of width 128 and ReLU or Tanh activations. The dynamics and reward models take the concatenated state–action input (x,u)∈ℝd+na(x,u) ^d+na and output, respectively, a d-dimensional state derivative/increment and a scalar reward; both use two 128-width hidden layers with ReLU. The policy network consumes the observation o∈ℝo ^\,o augmented with a scalar time-interval feature Δt t to condition actions on continuous-time step size, and produces an a-dimensional action through a 128–128–64 hidden stack with ReLU. Notation: d = state dimension, o = observation dimension, a = per-agent action dimension, n = number of agents, so the joint action has dimension nana. The value heads output scalars; the dynamics head outputs ℝdR^d; the policy head outputs ℝaR^a. Action squashing or clipping to environment bounds (if used) is applied after the final linear layer. Appendix F The Use of Large Language Models (LLMs) We employed LLMs as a writing assistant to polish the paper.