SMC-ES: Automated synthesis of formally verified control policies

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

SMC-ES: Automated synthesis of formally verified control policies

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.

SMC-ES: Automated synthesis of formally verified control policies

arXiv:2607.15003v1 [cs.AI] 16 Jul 2026

Riccardo Curcio, Toni Mancini, Enrico Tronci

Abstract—The deployment of autonomous cyber-physical systems in safety-critical environments requires closed-loop control strategies (i.e., policies) that are not only performant but also provably safe and robust. While learning-based methodologies such as Reinforcement Learning offer flexible and scalable approaches to automatically synthesize such controllers, they typically lack the formal guarantees necessary for safe deployment. To bridge this gap, we propose a novel simulation-based methodology to automatically synthesize policies with formal guarantees regarding performance, safety, and robustness specifications. Specifically, given a set of properties to verify, a confidence parameter δ and an allowable failure probability ε, our method guarantees that the synthesized policy comes with a certificate: with confidence at least 1 − δ, the probability of encountering a scenario where the given properties are violated is at most ε. We demonstrate the feasibility of our approach by developing SMC-ES, an algorithm that integrates Evolutionary Strategies with Statistical Model Checking-based verification. We evaluate SMC-ES on a suite of continuous control tasks using Gymnasium and Safety Gymnasium testbeds. Results show that, at the price of a sustainable increase in computational cost, our algorithm provides formal guarantees regarding performance, safety, and robustness specifications, while performing competitively against leading model-free Deep Reinforcement Learning (DRL) and Safe-DRL baselines. Index Terms—Reinforcement learning (RL), Neural Network (NN), Deep Reinforcement Learning (DRL), Safety-critical control, Evolutionary Strategy (ES), Statistical Model Checking (SMC)

I. I NTRODUCTION The ubiquity of autonomous systems in modern infrastructure (e.g., robotic manipulators, autonomous vehicles, UAVs) has highly intensified the demand for robust, closed-loop control strategies. Traditionally, control theory has relied on explicit model-based designs to provide theoretical guarantees of stability and safety. However, with system dynamics becoming progressively more difficult to characterize analytically, reliance on learning-based methodologies has naturally increased [1]. Despite the computational power available today, bridging the gap between the flexibility of learning-based control synthesis and the rigor of formal verification remains a central open problem. While standard model-free learning methods demonstrate impressive capabilities in simulation environments, the lack of strict formal guarantees regarding their worst-case behavior introduces risks for physical deployment that cannot be ignored [3]. This paper focuses on bridging this gap, enabling the automatic synthesis of control strategies that are formally verified to satisfy the system’s specifications. Riccardo Curcio is with the Department of Computer Science, Sapienza University of Rome, via Salaria 113, 00198 Rome, Italy Toni Mancini is with the Department of Computer Science, Sapienza University of Rome, via Salaria 113, 00198 Rome, Italy Enrico Tronci is with the Department of Computer Science, Sapienza University of Rome, via Salaria 113, 00198 Rome, Italy

A. Motivation Automated controller synthesis addresses the challenge of deriving executable control strategies (i.e., policies) directly from high-level specifications, either through analytical or learning-based methods. In this framework, the policy acts as the decision-making agent processing sensor data to drive the actuators, while the plant represents the physical hardware dynamics and the environmental disturbances acting upon them. The ultimate goal of controller synthesis methods is deployment in real-world applications. Over the past decades, learning-based methods, such as those based on Reinforcement learning (RL), have achieved remarkable success across diverse domains, e.g., transportation scheduling, recommender systems, and robotics (see [3] and citations therein). However, deploying these controllers in safety-critical systems requires strong guarantees regarding closed-loop behavior [4], [3]. Existing learning-based algorithms, particularly those scaling to complex systems, are generally unable to provide such guarantees. Indeed, extensive literature has studied and highlighted the inherent fragility of RL-based policies [5], [6], [7], [8]. This lack of robustness presents a substantial barrier to adoption in safety-critical domains, where the worst-case performance is often more significant than the average-case utility. In this work, we propose and evaluate a simulation-based methodology to automatically synthesize control policies with formal guarantees regarding performance, safety, and robustness specifications. B. Contribution Our main contributions are: 1) Automated synthesis of formally verified control policies. We present a methodology to automatically synthesize formally verified control policies from a simulation model of the plant and of its environment. The enforceable specifications include: • Performance: Ensure the system satisfies a lower bound on its operational metric. • Safety: Ensure the system avoids unsafe states throughout execution. • Robustness: Ensure the system maintains performance and safety guarantees even in presence of noise and/or faults in sensors and actuators. Consequently, the policy comes with a certificate: with confidence at least 1 − δ, the probability of encountering a scenario where the given specification is violated is at most ε. 2) Implementation: SMC-ES. We demonstrate the feasibility of our approach by introducing SMC-ES, an algorithm that couples Evolutionary Strategy (ES) with Statistical

Model Checking (SMC)-based verification. Inspired by the method in [9] for low-dimensional feasibility design, we employ an iterative loop where an ES algorithm drives the optimization search, while SMC [10] is used to verify the resulting policy. If a candidate policy violates the specification(s), counterexample scenarios are generated and used to guide further refinement(s). This loop continues until the sought statistical guarantees are achieved. Crucially, to scale to high-dimensional problems, we implement our ES algorithm by extending the Basic Random Search (BRS) algorithm [11] to enforce constraints over multiple scenarios, leveraging massive parallelism to optimize complex Neural Network (NN)based policies. 3) Experimental evaluation. We evaluate SMC-ES on a suite of continuous control tasks using Gymnasium and Safety Gymnasium testbeds [12], [13]. We compare performance of our formally verified policies against leading model-free Deep Reinforcement Learning (DRL) and Safe-DRL algorithms. Results show that our policies are consistently competitive in terms of performance, despite facing much stricter requirements and evaluation standards [7]. Notably, we demonstrate this also on the Humanoid environment, widely considered one of the most challenging Gymnasium locomotion task. Importantly, our policy guarantees that the Humanoid agent never falls regardless of its initial state. Furthermore, it achieves a formally verified cumulative reward lower bound of 12 240, ranking among the highest in the literature. To address sim-to-real challenges [14], [15], we also assess SMC-ES on these benchmarks under synthetic noise conditions. Results confirm that SMC-ES successfully synthesizes formally guaranteed policies that maintain strong performance even in the presence of perturbations. Additionally, we empirically assess the robustness of official DRL policies sourced from RL Baselines3 Zoo [16], [17]. By rigorously testing these baselines against the safety properties implicitly encoded within their own reward functions, we demonstrate that they generally fail to satisfy any safety guarantees, echoing common findings in the field [5], [6], [7]. We conclude with a computational evaluation, acknowledging the increased, yet tractable, computational budget required to achieve formal guarantees relative to DRL algorithms. II. R ELATED WORK In the context of cyber-physical systems, controllers must guarantee consistent, safe behavior under a given set of scenarios. For instance, an autonomous mobile robot should never enter unsafe states and should strictly respect operational limits (e.g., maximum velocity, torque, or energy consumption). This practical requirement motivates two closely related research objectives: safe exploration and safe deployment [3]. Safe exploration targets failure avoidance during the learning process itself, whereas safe deployment ensures the final, synthesized policy adheres to system’s specifications before execution. Crucially, the choice of the learning paradigm dictates which objective is paramount: online learning requires rigorous safe

exploration [18], [19], [2], while simulation-based learning circumvents these physical risks by confining training hazards to a virtual environment [14]. Ultimately, both paradigms aim at guaranteeing trustworthy operation in the real world i.e., safe policies. Current methodologies for synthesizing safe policies can be generally categorized into learning-based and analytical methods (see [3] and references therein). 1) Learning-Based Methods: These methods integrate safety constraints directly into the optimization objective, often prioritizing performance and scalability over rigorous guarantees. • Safe RL: Safe RL algorithms aim to maximize reward while maintaining safety constraint violations below a specified threshold. Some of the most prominent SafeDRL algorithms include CPO, CUP, and FOCOPS (see [13] and citations therein). While effective for safe exploration, these methods are not designed to strictly enforce safety constraints (i.e., zero violations). • Gaussian Process: Methods using Gaussian Processes (GPs) leverage predictive uncertainty to guide exploration and improve sample efficiency [20]. While some GPbased algorithms provide rigorous probabilistic safety guarantees [21], they are computationally intensive and restricted to low-dimensional tuning tasks. Approaches that scale to higher dimensions sacrifice these rigorous guarantees for tractability. 2) Analytical methods: These methods focus on deriving mathematically rigorous safety certificates, but at the cost of scalability or flexibility. • Control Theory: Drawing on classical control, these methods enforce safety using certificates such as Lyapunov or Control Barrier Functions (CBFs) to guarantee stability and constraint satisfaction [22]. While they provide rigorous assurance, they typically assume explicit knowledge of system dynamics and require the manual, analytical construction of valid certificates, limiting their applicability to complex systems [23]. • Formal Methods: These approaches synthesize controllers that satisfy logic-based property specifications (e.g., LTL formulas) by construction. However, classical synthesis algorithms often suffer from state-space explosion and assume discrete states or full knowledge of the environment. These restrictive assumptions often hinder their direct applicability to high-dimensional, continuous realworld tasks [24]. Despite substantial progress, a critical gap remains: safe deployment of automatically synthesized policies for complex systems lacks rigorous solutions. Given the unavailability of mathematical guarantees, rigorous empirical testing becomes paramount. In this context, continuous control testbeds such as Gymnasium and Safety Gymnasium [12], [13] have become instrumental in allowing researchers to assess controllers’ performance empirically and standardize comparisons. III. BACKGROUND In this manuscript, we denote with Z, R and R+ the set of integer, real and strictly positive real numbers, respectively.

For any n-dimensional vector x, we denote its i-th element as xi , hence x = [x1 , . . . , xn ]. When comparing two vectors a, b ∈ Rn , the notation a ≤ b denotes a component-wise inequality, meaning that ai ≤ bi holds for all i ∈ [1, n]. Additionally, we use 0(n) to represent the zero n-dimensional vector. Given a random variable X, the expression x ∼ X indicates that x is a specific realization of X. Furthermore, we denote with Pr(E) the probability that a logical predicate or event E holds true, and with ∆(X) the space of all probability distributions over a given support set X. A. Reinforcement learning Reinforcement learning (RL) models sequential decisionmaking as an interaction between an agent and an environment, formalized as a Markov Decision Process (MDP). An MDP is defined as a tuple (X, U, R, P, γ, I), comprising a set of states X ⊆ Rn (n > 0), a set of actions U ⊆ Ra (a > 0), a reward function R : X × U × X → R, a stochastic transition dynamics model P : X × U → ∆(X), a discount factor γ ∈ [0, 1) and an initial state distribution I ∈ ∆(X). Given an initial state x 0 ∼ I, at each time step t the agent observes its state x t ∈ X, selects an action u t ∈ U governed by a policy π (which may be deterministic π : X → U , or stochastic π : X → ∆(U )), transitions to the next state x t+1 ∼ P (x t , u t ) and receives a scalar reward rt = R(x t , u t , x t+1 ). The agent’s objective is to find a policy π that maximizes P∞ the expected discounted cumulative reward J(π) = E[ t=0 γ t rt ]. Deep Reinforcement Learning (DRL) extends this framework by integrating MDP formulations with Deep Neural Networks (DNN). This allows agents to learn complex representations from high-dimensional inputs, enabling them to scale to tasks that were previously computationally infeasible. Policy gradient algorithms [25] form the backbone of DRL. Research has since progressed to address key limitations: improving stability and scalability via on-policy updates (i.e., learning exclusively from current experience e.g., TRPO [26], PPO [27], A2C [28]), and maximizing sample efficiency through off-policy mechanisms (i.e., reusing collected past experiences e.g., TD3 [29], SAC [30], TQC [31]). More recently, performance has been further improved by integrating diffusion models (DSAC-T [32], DACER [33]) and advanced representation learning (TD7 [34]). However, despite these empirical successes, DRL methods fundamentally lack mechanisms to guarantee consistent, safe policy behavior under a given set of scenarios. B. Evolutionary strategy Evolutionary Strategy (ES) is a class of Black Box Optimization (BBO) algorithms [35] that are heuristic search procedures inspired by natural evolution. In contrast to RL, ES offers a population-based paradigm: instead of relying on backpropagated gradients, this approach evolves a distribution of potential solutions and uses them to empirically estimate a search direction. The most widely known member of the ES class is Covariance Matrix Adaptation Evolution

Strategy (CMA-ES) [36], which represents the population by a full-covariance multivariate Gaussian. CMA-ES has been extremely successful in solving optimization problems in low to medium dimensions. A core appeal of ES lies in their minimal reliance on differentiability or an explicit MDP formulation. Since ES require only a scalar fitness (i.e., objective function) value, they are applicable to discrete, non-differentiable, or simulator-based problems. From a systems perspective, ES are straightforward to parallelize with negligible communication overhead, a property effectively leveraged by OpenAI-ES [37] to scale optimization to high-dimensional tasks. Furthermore, recent ES algorithms (ARS [11]) have demonstrated competitive performance relative to standard DRL baselines by extending a simple Basic Random Search (BRS) algorithm. Borrowing from the MDP notation established in Section III-A, in the BRS setting the optimization problem is formulated as argmax Eξ [r(πθ , ξ)]. Here, πθ : X → U θ∈Rd

denotes a deterministic policy parameterized by parameters θ (formally defined later in Definition 3), the random variable ξ encodes the randomness of the environment (i.e., random initial states and stochastic transitions) and r(πθ , ξ) is the (undiscounted) cumulative reward achieved on one trajectory generated from the system. To solve this problem, BRS shifts exploration to the parameter space: it optimizes a smoothed objective function argmax Eζ Eξ [r(πθ+σζ , ξ)] by introducing a θ∈Rd

zero-mean Gaussian perturbation vector ζ ∼ N (0(d) , I(d×d) ) (with N denoting a Gaussian distribution and I the identity matrix), and a positive exploration noise scale (i.e., standard deviation) σ > 0. These perturbations allow the algorithm to estimate the search direction via antithetic sampling [38], effectively guiding optimization. Notably, the smoothed objective function provides an accurate approximation of the original objective function as the scaling parameter σ becomes small [11]. However, ES are not without limitations. They are sample-inefficient compared to DRL algorithms, requiring extensive function evaluations to estimate the search direction [39]. Consequently, ES is generally computationally expensive and, under ineffective recombination schemes, remain susceptible to getting trapped in local optima. C. Black-Box Optimization A Black Box Optimization (BBO) problem is a tuple (J, C) where: d • J : R → R (objective function) and d • C : R → Rm (constraints violation function), with d, m ≥ 0. A solution (if any) to the BBO problem (J, C) is a real vector θ∗ ∈ Rd such that: 1) C(θ∗ ) ≤ 0 (all constraints are satisfied) 2) ∀θ ∈ Rd [(C(θ) ≤ 0) → (J(θ∗ ) ≥ J(θ))] (the objective function is maximized). 1) Augmented Lagrangian: Augmented Lagrangian methods are constraint-handling approaches that combine penalty functions with the Karush-Kuhn-Tucker (KKT) necessary con-

ditions for optimality [40]. Formally, we define the Augmented Lagrangian function h in terms of θ as follows: m X µi 2 h(θ; λ, µ) = J(θ) − λi C̃i (θ) + C̃i (θ) , 2 i=1 where C̃i (θ) = max{Ci (θ), −λi /µi }, and λ ∈ Rm and µ ∈ Rm + are the Lagrange multipliers and positive penalty coefficients, respectively. Augmented Lagrangian approaches have been widely adopted in the context of ES [40], [41], [42]. In this work, we adopt the specific adaptive scheme proposed in [40] and implemented in the pycma [43] python package. To formalize this constraint-handling mechanism, we introduce the operator LES . This operator acts directly upon the objective value j = J(θ) and constraint values c = C(θ). Consequently, we (re)-define the scalar fitness function for a given realization (j, c) as: m X µi 2 h(j, c; λ, µ) = j − λi c̃i + c̃i 2 i=1 m

Rm +

where c̃i = max{ci , −λi /µi }, and λ ∈ R and µ ∈ are the Lagrange multipliers and positive penalty coefficients, respectively. The LES operator is initially configured with parameters λ and µ. Then, at each call, it receives as input a population of η pairs {(jk , ck )}ηk=1 and performs two operations: 1) Parameter Update: It updates the internal multipliers λ and penalties µ based on the population statistics (i.e., the average of objective values and constraints’ values), following the adaptive scheme detailed in [40]. 2) Fitness Evaluation: It returns the (vector) of penalized fitness values computed via the evaluation of h(jk , ck ; λ, µ) for each element of the population. D. Statistical Model Checking Statistical Model Checking (SMC) (see, e.g., [44] for an overview) is a set of methods and algorithms to verify that, with a given statistical confidence, a system satisfies given requirements. We will use Monte Carlo–based SMC to compute an upper bound to the probability that a given function is not identically 0. To that end, we employ the same algorithm following from the theorem in [9]. Theorem 1 (from [9]). Let X be an n-dimensional real-valued random variable, f : Rn → {0, 1}, and ε, δ ∈ (0, 1]. Then, an algorithm exists such that: 1) it terminates; 2) it returns a finite set (of counterexamples) Ω ⊂ Rn such that ∀ω ∈ Ω f (ω) = 1; 3) if Ω = ∅, then, with probability at least (1 − δ) we have that Pr f (X) = 1 ≤ ε. We denote by verify(f, ε, δ) the algorithm from [9]. Its guarantee can be stated succinctly: for any predicate f and any chosen ε, δ ∈ (0, 1], if verify(f, ε, δ) terminates without producing a counterexample then, with confidence at least 1 − δ, the probability that f (X) = 1 is at most ε.

IV. P ROBLEM S TATEMENT In the following we deal with discrete-time dynamical systems. Let T = N be our (discrete) time set. Given a set of values A, we denote with a ∈ AT a function associating to each time point t ∈ T a value at ∈ A (time function). Definition 1. A discrete-time dynamical system is a tuple S = (X, U, W, ϕ) where n • X ⊆ R is a set of states. a • U ⊆ R is a set of control inputs (actuations). q • W ⊆ R is a set of uncontrollable inputs (e.g., disturbances from the environment). • ϕ : X × U × W → X is a state-transition function. For any time point t ∈ T, the state at time point t + 1 is defined as x t+1 = ϕ(x t , u t , w t ) where x t ∈ X is the state at time t, u t ∈ U is the control input at time t, and w t ∈ W is the disturbance at time t. The sequence of such states defines a time function x ∈ X T (i.e., system trajectory). Thus, a dynamical system S = (X, U, W, ϕ) as in Definition 1 takes as (uncontrollable) inputs an initial state x 0 ∈ X and a disturbance time function w ∈ W T . Definition 2. An environment model for a system S (Definition 1) is a pair EnvS = (I, W), where: • I ∈ ∆(X) is an initial state distribution (a probability distribution over the state space X of S) • W is a discrete-time real-valued stochastic process, taking values in the set of uncontrollable inputs W . Any realization (x 0 , w) ∼ EnvS is an actual operational scenario for S, defining both the initial state x 0 and a specific time function of disturbances w under which S operates. We denote by (x 0 , w) ∈ EnvS a scenario entailed by EnvS , and with (x 0 , w) ∼ EnvS a sample scenario drawn from EnvS according to its underlying distributions. Next, we define a control strategy (i.e., a policy) that generates a (control) input time function u ∈ U T for a system S (Definition 1) as follows. Definition 3. A deterministic policy for a system S (Definition 1) is a function: πθ : X → U parametrized by θ ∈ Rd , which assigns a control action πθ (x ) ∈ U to each state x ∈ X. Lastly, we define transition indicators to characterize rewards and constraint violations associated with each statetransition. Definition 4. Given a discrete-time dynamical system S as in Definition 1, we define Transition indicators (r, c) for S, consiting of: • r : X ×U ×X → R, representing the transition’s reward.

c : X × U × X → Rm , representing violation costs for m distinct safety constraints.

(x 0 , w) ∈ EnvS . Therefore, we formalize our optimal safe policy design problem as follows.

The combination of a discrete-time dynamical system S, a control policy πθ and an environment model EnvS defines the following closed-loop system.

Definition 7. An optimal safe policy design problem is a tuple Γ = (S, EnvS , (r, c), d, H), where: • S is a discrete-time dynamical system. • EnvS is an environment model for S. • (r, c) are Transition indicators for S. • d is the number of dimensions of the policy parameter space Rd • H ∈ T is a time horizon This tuple formulates the following constrained optimization problem:

•

Definition 5. A Closed Loop System (CLS) is a tuple S cl = (S, EnvS , πθ ), where • S is a discrete-time dynamical system. • EnvS is an environment model for S. • πθ is a policy for S. Given a scenario (x 0 , w) ∼ EnvS , the closed-loop system S cl evolves at each time point t ∈ T as follows: x t+1 = ϕ x t , πθ (x t ), w t for all t ≥ 0.

(1)

Remark 1. We point out that in this work we adopt discretetime dynamical systems (Definition 1), rather than the classical MDP formulation used in RL. This choice allows for the explicit characterization of trajectories (including disturbances), which is essential to record failures. Note, however, that this formulation retains the full expressive power of an MDP: any stochastic transition function can be equivalently represented by a deterministic state-transition function driven by a stochastic disturbance process, exactly as captured by our environment model (Definition 2) [45]. We define the following trajectory evaluation to value the behavior of a CLS S cl under a given scenario (x 0 , w) ∼ EnvS . Definition 6 (Policy trajectory evaluation). Let S cl = (S, EnvS , πθ ) be a CLS and let (r, c) be transition indicators for S. Given a scenario (x 0 , w) ∼ EnvS and a time horizon H ∈ T, the evaluation of policy πθ over the resulting trajectory x ∈ X T (Equation (1)) is the pair J(πθ , x 0 , w, H), C(πθ , x 0 , w, H) , where: J(πθ , x 0 , w, H) = C(πθ , x 0 , w, H) =

H−1 X t=0 H−1 X

r x t , πθ (x t ), x t+1

(2)

n o max 0(m) , c x t , πθ (x t ), x t+1 .

t=0

(3) J(πθ , x 0 , w, H) and C(πθ , x 0 , w, H) are, respectively, the cumulative reward and the (vectorial) cumulative constraints’ violation along x (in the formula for C, max is evaluated element-wise and results in an m-dimensional vector). Remark 2. Equation (2) in Definition 6 aligns with the standard cumulative return used in RL. However, rather than embedding constraints within the reward formulation as commonly done in RL practices, we model them explicitly. This decouples the problem definition from the solution method, enabling algorithms to handle constraints more efficiently. Ideally, we would like to find an optimal policy that satisfies constraints for every scenario (uncertainty realization)

maximize γ subject to J(πθ , x 0 , w, H)≥γ C(πθ , x 0 , w, H)=0(m)

∀(x 0 , w)∈EnvS

(4) (5)

∀(x 0 , w)∈EnvS

(6)

θ∈R , γ ∈R.

(7)

For every feasible solution (θ, γ) to the problem above, tuple (S, EnvS , πθ ) defines a CLS (Definition 5) whose policy is safe, i.e., always satisfies all safety constraints (cumulative violation cost equal to 0(m) under all possible operational scenarios), and, at the same time, guarantees a minimum cumulative reward of γ. We denote with (θ∗∗ , γ ∗∗ ) the optimal solution, where θ=θ∗∗ constitutes an optimal safe policy πθ and γ ∗∗ is its associated (maximum) objective value. V. A LGORITHM DESIGN For systems with general nonlinear dynamics and continuous uncertainty spaces, even the feasibility problem behind that of Definition 7 (with a given value for γ used as a threshold) is undecidable, let alone the maximization problem. In the following, we propose methods and an algorithm aimed at computing a solution to our optimal safe policy problem which is quality-guaranteed with high-enough statistical confidence. Specifically, given a confidence parameter δ ∈ (0, 1] and an error tolerance ε ∈ (0, 1], our algorithm aims to find an (ε, δ)-solution (θ∗ , γ ∗ ) to the given optimal safe policy design problem. This solution guarantees, with confidence at least 1 − δ, that the probability of encountering a scenario where πθ∗ violates safety constraints or yields a cumulative reward below γ ∗ (the minimum guaranteed cumulative reward declared by the algorithm for πθ∗ ) is at most ε. We formally define an (ε, δ)−solution as follows. Definition 8 ((ε, δ)-solution). Let Γ = (S, EnvS , (r, c), d, H) be an optimal safe policy design problem (Definition 7). Given ε, δ ∈ (0, 1], a pair (θ, γ) is an (ε, δ)−solution to Γ if, with probability at least 1 − δ, the following holds:   Performance violation }| { z    J(πθ , x 0 , w, H) < γ      ∨ Pr  ≤ ε.    x ,w ∼ Env ( 0 ) S  C(πθ , x 0 , w, H) > 0(m)   | {z }  Safety violation

An (ε, δ)−solution (θ, γ) (Definition 8) is a solution satisfying safety constraints over fewer scenarios compared to the original optimal safe design problem Γ (Definition 7). Consequently, an optimal (ε, δ)−solution (θ∗ , γ ∗ ) serves as an upper bound (super-optimal, [46]) on the true global optimum γ ∗∗ for Γ (i.e., γ ∗ ≥ γ ∗∗ ). While this makes γ ∗ an optimistic estimate of the performance lower bound, the performance violation property ensures that the risk of γ ∗ being too optimistic can be arbitrarly reduced acting on parameters ε, δ. In what follows, we slightly abuse notation by referring to a solution (θ, γ) simply as (πθ , γ), directly referencing the parameterized policy rather than its parameters. A. Main algorithm and policy verification Algorithm 1 shows our main algorithm, which, in the spirit of [9], repeatedly alternates a solving phase and a SMC–based verification phase. Algorithm 1: (ε, δ)-solving an optimal safe policy design problem input Γ = (S, EnvS , (r, c), d, H); // the problem input ε, δ ∈ (0, 1]; // error and confidence thresholds 3 Σ ← a non-empty set of scenarios (x 0 , w) ∼ EnvS ; 4 while computation budget not exhausted do 5 πθ∗ , γ ∗ ← solve(ΓΣ ); 6 if πθ∗ =⊥ then break; 7 cntrex ← verify(πθ∗ , γ ∗ , ε, δ); // return counterex’s 8 if cntrex = ∅ then return (πθ∗ , γ ∗ ); // (ε, δ)-solution 9 else Σ ← Σ ∪ cntrex; // include counterexamples 10 return ⊥; // no solution found 1

Proof. Function verify(πθ∗ , γ ∗ , ε, δ) in Algorithm 1 implements the algorithm described in the proof of [9, Theorem 1], verify(f, ε, δ), with (implicit) predicate f being: J(πθ∗ , x 0 , w, H) < γ ∗ ∨ C(πθ∗ , x 0 , w, H) > 0(m) . Therefore, upon termination (see Line 8 in Algorithm 1), the solution (πθ∗ , γ ∗ ) satisfies the property as in Definition 8. We point out that forcing termination of Algorithm 1 when a given computational budget is exhausted is a necessary precaution. This is because, due to the undecidability of the underlying problem, it is theoretically possible that function verify() keeps finding counterexamples for infinitely many iterations. B. Effective policy finding through ES We aim to synthesize highly parameterized NN-based policies (each one defined by parameters θ ∈ Rd ). Hence, function solve() used in Algorithm 1 must scale to highly-dimensional parameter spaces Rd . Furthermore, since our sequence of (finitized) optimization problems ΓΣ are expected to define constrains over many scenarios, the optimizer requires robustness against noisy performance evaluations. To this end, we build upon the Basic Random Search algorithm described in (Section III-C) and extend it with a few standard practical devices. Algorithm 2: Function solve(ΓΣ ) based on ES

param α ∈ R+ ; // learning rate θ ← initial policy (NN) parameters; ∗ ∗ 3 (θ , γ ) ← (⊥, −∞); // current optimum 4 while computation budget not exhausted do 5 γ ← min(x 0 ,w)∈Σ J(πθ , x 0 , w, H) ; The algorithm starts with a finite collection Σ of operational 6 C ← max(x 0 ,w)∈Σ C(πθ , x 0 , w, H) ; scenarios (any finite collection would work in theory, but // possibly update current optimum practical devices to ensure effectiveness are discussed in Sec- 7 if C = 0(m) ∧ γ > γ ∗ then (θ∗ , γ ∗ ) ← (θ, γ); tion V-B), and attempts to solve a finitization ΓΣ of problem Γ 8 ĝ ← estimate dir(θ, ΓΣ ); // estimate search direction where constraints J(πθ , x 0 , w, H) ≥ γ and C(πθ , x 0 , w, H) = 9 θ ← θ + α × ĝ; // perform move ∗ ∗ 0(m) are defined only for scenarios (x 0 , w) ∈ Σ. 10 return (θ , γ ); If the solving phase fails, the algorithm terminates. Otherwise, the solution found (πθ∗ , γ ∗ ) is verified through Statistical First, since we account for multiple scenarios Σ, we must Model Checking, using the algorithm described in [9, Theoevaluate and aggregate the policy πθ performance across the ∗ rem 1]. Specifically, function verify(πθ∗ , γ , ε, δ) repeatedly entire set of scenarios. The most natural aggregation choices samples scenarios in EnvS . For each scenario, it checks are either the average (average objective function and average whether the candidate policy πθ∗ satisfies the constraints and constraints’ violation) or the worst-case (minimum objective ∗ yields a cumulative reward of at least γ (Definition 8). function and maximum constraints’ violation) (see line 10 In case the function finds counterexamples (set cntrex), from Algorithm 3). However, estimating a search direction these are added to Σ, before a new iteration of the algorithm (see Section III-B) using the full set of scenarios Σ quickly is performed on the enlarged problem. Conversely, whenever becomes computationally prohibitive as the set grows. To function verify() reaches its termination condition without mitigate this, we use a lightweight heuristic: at the start of each finding counterexamples, the algorithm terminates, since the iteration t, given a policy πθ , we identify the subset Σ̃ ⊆ Σ found policy πθ∗ is an (ε, δ)-solution to the problem, as stated by selecting the Pareto-worst [47] scenarios among the current by Proposition 1. set. Specifically, this subset comprises the scenarios exhibiting Proposition 1 (Correctness of Algorithm 1). A pair (πθ∗ , γ ∗ ) the most adverse combinations of minimum (worst) objective returned by Algorithm 1 is an (ε, δ)-solution to the given values j and maximum (worst) constraints’ violations c. We optimal safe policy design problem Γ. then estimate the search direction using only that subset (see 1 2

Algorithm 3: Function estimate dir(θ, ΓΣ ) 1 param mode ∈ avg, worst ; // dir. estimation mode 2 param σ ∈ R+ ; // noise standard deviation 3 param n ∈ N+ and even; // nb. of policy perturbations 4 Σ̃ ← Pareto-worst scenarios in Σ for πθ ; 5 J̃ ← empty n-dim. real vector; 6 C̃ ← empty n × m-dim. real matrix; n 7 for i from 1 to 2 do 8 ζi ← draw sample from N (0(d) , I(d×d) ); 9 (θ2i−1 , θ2i ) ← θ ± ζi σ; 10 for j ∈ {2i − 1, 2i} do /* aggregate cumulative rewards and constraint violations across scenarios Σ̃ using avg if mode = ‘avg’ and min / max otherwise */ 11 J̃j ← avg | min (x 0 ,w)∈Σ̃ J(πθj , x 0 , w, H) ; 12 C̃j ← avg | max (x 0 ,w)∈Σ̃ C(πθj , x 0 , w, H) ; 13

h ← LES (J̃, C̃); // compute (vector) of scalar fitnesses (see Section III-C1) P n2 return n1 i=1 (h2i − h2i−1 )ζi

line 10 from Algorithm 3). This concentrates computation on the most informative/critical scenarios while keeping computational cost manageable. Lastly, we convert the constrained ΓΣ̃ problem into an unconstrained one using the augmented-lagrangian operator LES described in Section III-C1 (see Line 13 from Algorithm 3). We refer to Algorithm 1 when instantiated with the solver from Algorithm 2 as SMC-ES. VI. E XPERIMENTS Our goal is to assess whether our algorithm can learn formally verified policies while achieving competitive performance with respect to DRL baselines. To evaluate our methodology, we conduct experiments on three distinct types of MuJoCo environments. First, we use the subset of standard MuJoCo tasks [12] that feature explicit unhealthy (i.e., unsafe) states. Second, we test on the SafetyVelocity MuJoCo suite [13], which imposes a strict velocity limit. This setting is particularly significant because the reward function incentivizes speed maximization, creating a direct conflict between task performance and safety compliance. Lastly, to enforce robustness specifications and tackle sim-to-real challenges [14], we also conduct experiments on all these tasks augmented with synthetic noise (i.e., disturbances W ). We adopt default environment configurations, including default initial state distribution I (Definition 2). However, we exclude the healthy reward term (when present) from the transition reward function (Definition 4), a bonus reward given by default in MuJoCo for reaching a safe state. We reframe this as a strict requirement and therefore as a safety constraint. To ensure fair comparison with baseline DRL algorithms, we reintegrate this reward component back when reporting performance metrics.

A. Implementation details Training uses Multi Layer Perceptron (MLP) policies implemented in PyTorch [48], with a tanh activation on the output layer. Adopting the standard practices in [37], we apply column-normalized Gaussian initialization for the NN parameters θ. Furthermore, we normalize the NN inputs using the running average and standard deviation of the sensor readings (i.e., agent’s observations). We continuously update these statistics by randomly incorporating trajectories seen during the search direction estimation (see Algorithm 3) with a small probability p. To allow deployment on physical hardware, we use PyTorch’s native export functionality to convert the trained policy into an Open Neural Network Exchange (ONNX) model [49]. This representation enables efficient inference across a wide range of hardware platforms, either through direct execution with ONNX Runtime or via compilation into portable C/C++ code for embedded targets using converters such as onnx2c [50] or full-stack compilers like TVM [51]. At the start of our synthesis process we initialize the set of scenarios Σ (see Algorithm 1 line 3) by sampling j candidates from the initial state distribution I and applying the K-Medoids algorithm [52] (implemented in scikit [53]) to select the k most representative medoids. Then we start our optimization algorithm. We use Algorithm 2 and additionally employ centered rank-based fitness on the objective function [54] to mitigate influence of outliers. Moreover, we use the AdamW optimizer [55] (implemented in pytorch [48]) to automatically regulate the update step size α via adaptive moment estimation and decoupled weight decay (Algorithm 2 line 9). Exploiting the inherent parallelization capabilities of ES, we distribute policy evaluations (i.e., independent simulations) across multiple machines. To optimize parallel efficiency and mitigate synchronization latency caused by variable episode durations, we employ an adaptive simulation horizon hsteps . This horizon progressively expands together with policy improvement, eventually reaching the environment’s limit H [37]. Each call to the solver (see line 5 of Algorithm 1) is constrained by both a maximum iteration limit and a patience mechanism, which terminates optimization if no improvement is observed within a set number of steps. Furthermore, we adapt the optimization strategy across sequential solver calls. The first call uses average aggregation (see Algorithm 3) and a higher learning rate α for faster convergence. Upon obtaining a candidate solution (πθ , γ) from Algorithm 2, we verify if it is indeed an (ε, δ)−solution using the algorithm verify(πθ∗ , γ̂, ε, δ) (with ε, δ = 1%), but setting the target performance lower bound to γ̂ = γ(1 − τ ) with τ ∈ [0, 1) and retaining a specified maximum number of counterexamples as done in [9] (if any). Here, τ acts as a slight tolerance to accommodate minor performance variations. If verification fails, subsequent solver calls retain the LES multipliers λ, µ, but switch to worst-case aggregation (see Algorithm 3) with a reduced learning rate α′ ≤ α to facilitate refinement(s). These subsequent solver calls also employ early stopping to

TABLE I: Standard MuJoCo evaluation: SMC verification of pre-trained official RL Baselines3 Zoo [16] policies. Columns show whether verification of the safety property succeeded (✓) or failed (✗). “-” indicates unavailable policies from RL Baselines3 Zoo [16]. Environment Humanoid Ant Walker2d Hopper

A2C ✗ ✗ ✗ ✗

PPO ✗ ✗ ✓

SAC ✗ ✗ ✗ ✗

TD3 ✗ ✗ ✓ ✓

ARS ✓ ✗ ✓

TQC ✗ ✗ ✗ ✓

TABLE II: Standard MuJoCo evaluation: comparison of SMCES performances with respect to DRL baselines. Multiple values denote diverging results reported in the literature. Complete details regarding DRL algorithms can be found in [34], [33], [32].

TRPO

Algorithm

Humanoid

Ant

Walker2d

Hopper

✗ ✗ ✗

SMC-ES ✓

12 240

6505

6482

3763

10 829 11 888 10 281 8951 9335 / 6555 8361 6869 5631 / 5433 965 5291

7086 9108 10 133 8547 6427 / 4615 6329 6156 6184 / 5589 6203 4549

6424 6701 7397 6379 6200 / 5681 6137 4831 5237 / 5078 5502 4095

3688 / 3660 4104 4075 3423 2483 / 3167 3462 2647 / 1679 3569 / 3682 3474 / 3138 2644 / 2933

enhance efficiency, terminating optimization immediately once the objective function reaches at least γ̂ (which then remains the target lower bound for verification).

DSAC-T ✗[32] DACER ✗[33] TD7 ✗[34] TD3-OFE ✗[34] SAC ✗ TQC ✗ PPO ✗ TD3 ✗ TRPO ✗ DDPG ✗

B. Policy quality evaluation In this section we use the symbols (✓) and (✗) to denote if a policy is formally verified (through SMC with ε, δ = 1%) or not, respectively. Moreover, when comparing with DRL baselines we highlight the 1st , 2nd , and 3rd best results. 1) Standard MuJoCo environments: Here we empirically test DRL policies and evaluate SMC-ES on all the standard MuJoCo environments where the agent is able to reach an unsafe state (failure). a) Baselines: We evaluate pre-trained DRL-based policies using the same SMC verification procedure verify() as in our experiments (setting ε, δ = 1%), but exclusively checking for the safety violation (thereby omitting our stricter performance lower bound guarantees, see Definition 8). This metric is well-suited for standard DRL agents, as their reward function naturally encodes a safety specification through healthy state reward bonuses. Specifically, we test all the official (MuJoCo) policies available from RL Baselines3 Zoo [16], a training framework designed for Stable Baselines3 [17]. Table I empirically demonstrates that these policies generally fail to provide robust behavioral assurances (in accordance to studies regarding RL instability [5], [6], [7]). b) Results: Table II benchmarks our results against DRL baselines sourced from recent state of the art RL literature [34], [33], [32]. In cases where discrepancies arise between sources, due to environment versions mismatch or stochastic training outcomes, we preserve all reported values to ensure a transparent comparison. It is crucial to contextualize our results against the evaluation standards used for the baseline methods. Standard DRL algorithms typically report empirical averages derived from peak performance over the final 10% of training, evaluated across 5 to 10 random seeds (i.e., scenarios). In contrast, our reported metric is a formally verified cumulative reward lower bound, a rigorous guarantee that requires succesful validation across well over 1000 consecutive random scenario evaluations. Despite being subject to these more conservative evaluation criteria, SMC-ES achieves highly competitive performance. Specifically, our algorithm achieves the top-1 rank in the (MuJoCo most challenging) Humanoid task and top-3 positions in Hopper and Walker2D, while remaining competitive in the Ant environment. We exclude Inverted Pendulum and

Inverted Double Pendulum, as they are absent from the results reported by TD7 [34] and are generally regarded as trivial tasks. However, for the sake of completeness, we report that SMC-ES successfully guarantees verified performance lower bounds of 1000.00 and 9358.94 for these two environments, respectively. It is worth noting that, by directly comparing our verified performance lower bounds against the self-reported metrics of baseline papers, we evaluate our method against the strongest possible representations of prior work [11]. 2) SafetyVelocity MuJoCo environments: In this section we evaluate SMC-ES on all the SafetyVelocity MuJoCo environments [13]. These tasks impose a strict velocity limit that directly conflicts with the standard speed-maximizing reward definition. This creates a rigorous test of safety compliance under opposing incentives. Additionally, for agents capable of entering unsafe states, we enforce a dual safety constraint: the policy must simultaneously satisfy the velocity threshold and maintain agent safeness at each time step. a) Results: Table III compares our formally verified policies against prominent safe-DRL baselines. Baseline results are sourced from [13]. Quoting the source literature, these values are computed as average values over 10 assessment iterations across multiple random seeds (i.e., scenarios). Although the considered safe-DRL algorithms are not designed to strictly enforce safety constraints (i.e., zero violations), we include them as a benchmark to assess relative performance. Results show that SMC-ES is able to consistently provide formally verified safety guarantees (i.e., zero violations), while achieving highly competitive performance with respect to the lowest-violation baselines. 3) MuJoCo environments augmented with noise: While simulations provide safe and data-rich training environments, modeling discrepancies often prevent policies from transferring successfully to the real world. To assess the efficacy of our method against these sim-to-real challenges [14], we evaluate SMC-ES on the same MuJoCo tasks used in Section VI-B1 and Section VI-B2, augmented with synthetic noise (i.e., disturbances W ). Specifically, we apply Gaussian perturbations w ∼ N (0, ν 2 ) at each time step: multiplicatively to sensors (i.e., NN inputs) by a factor of (1 + w ), and additively to actuations (i.e., NN outputs, scaled by the actuation range and

TABLE III: SafetyVelocity MuJoCo evaluation: comparison of SMC-ES performance with respect to Safe-DRL baselines results sourced from [13]. Top 3 best results are prioritized first by lower constraint violations and second by higher cumulative reward (i.e., Performance). Complete details regarding the Safe-DRL algorithms can be found in [13]. Algorithm

Environment

Violations

Performance

SMC-ES ✓