scheme. Autonomous agents deployed in safety-critical environments must operate with a clear awareness of safety constraints embedded in their dynamics. Barrier certificates (BCs) are a powerful formalism for verifying such safety properties in continuous dynamical systems. However, synthesizing BCs remains a computationally dificult problem. In this work, we propose a novel approach for constructing BCs using the Augmented Lagrangian framework, combined with Bayesian optimization over Gaussian Process models to eficiently guide the search for valid certificates. We demonstrate our algorithm on several test cases, illustrating the success of our proposed Barrier certificate, Dynamical systems, Gaussian process, Safety verification Autonomous agents operating in complex environments must not only achieve task performance but also remain continuously aware of safety constraints that govern their permissible behaviors. Dynamical systems theory, which studies the long-term behavior of evolving systems [1], provides a formal foundation for modeling such agents. This is especially relevant in hybrid systems, where continuous and discrete dynamics interact to define the agent's trajectory in the physical world, as commonly encountered in cyber-physical systems.
https://ge-xinyu.github.io/ (X. Ge); https://hycodev.com/ssoudjani (S. Soudjani); https://pzuliani.github.io/ (P. Zuliani)
CEUR Workshop
ISSN1613-0073 time-delay systems is proposed in [17]. Recent application of barrier certificates for checking properties of quantum systems is addressed in [18].
With the advances in learning methods, various data-driven approaches have emerged to address safety verification and synthesis of autonomous agents, showcasing the potential integration of machine learning with traditional formal methods. Data-driven abstraction with formal guarantees is studied in [19, 20] for satisfying temporal specifications. Data-driven reachable set computation with formal guarantees is studied in [21]. Learning algorithms for computing BCs that are in the form of neural networks are also proposed [22, 23]. Safety verification of autonomous agents with unknown dynamics using noisy data from observations, Gaussian process regression, and abstracting the system into a ifnite Markov model has been studied in [ 24].
Gaussian process (GP) regression is an eficient and accurate data-driven approach for approximating computationally complex functions. An advantage of GPs is their ability to ofer an eficient analytical approximation over the entire domain of the target function. Furthermore, other approaches such as deep neural networks typically need vast amounts of training data to ensure convergence of the learning phase. This contrasts with GPs, which in many cases can learn good approximations even with relatively small quantities of training data. GPs are used in optimization through a Bayesian approach known as Gaussian process optimization (GPO) [25], which is often more sample-eficient than traditional optimization methods.
We concentrate on using GPO for the safety verification of a given autonomous agent through the computation of BCs. Building upon insights from [29] and leveraging the augmented Lagrangian framework as described in [30], we propose a framework that transforms the BC synthesis problem into a constrained optimization problem using the augmented Lagrangian method. Our approach computes the parameters of a BC from a fixed template through GPO for handling constraints. Given that the constraints are in the form of parameterized optimizations, we fit a GP to estimate their values from a finite number of evaluations. This gives a probabilistic interpretation of the constraints, which are transformed to chance-constraints and subsequently incorporated into a reliability-based design optimization (RBDO) [31]. We further validate the feasibility and efectiveness of our method through implementation in three case studies. The main contributions of this article are threefold: 1. We introduce an optimization approach that synthesizes a BC from a parameterized set of functions for autonomous agents evolving continuously in time. We reconfigure the problem into a robust parameter optimization task and use augmented Lagrangian.
the functions appearing as maximum over a continuous domain and satisfy the related constraints.
safety verification, BCs, and Gaussian process regression is provided in Section 2. Our proposed solution approach is presented along with the algorithm in Section 3 containing the details of the augmented Lagrangian, RBDO, and verification of the BC using SMT Solvers. Case studies are provided in Section
2. Preliminaries and Problem Statement In the following, ℝ denotes the set of real numbers; arg min(⋅)returns an argument that minimizes the input function, max(⋅)represents the largest value taken by the input function; Φ(⋅)is the cumulative distribution function of the standard normal distribution, and (⋅) is the probability density function of the standard normal distribution. 2.1. Autonomous Agents and the Safety Problem An autonomous agent modeled as a -dimensional dynamical system over the continuous state space ⊂ ℝ is described by
̇ = ( ), where = [ 1, … , ] is a column vector, ̇ denotes the derivative of with respect to time, ( ) = [ 1( ), … , ( ) ] is a Lipschitz-continuous vector field, and ⊆ ℝ is an open set defining the state space of the system. The region for the initial state of the system 0 is denoted by 0; the unsafe region of the system is denoted by .
Definition 1 (Safety verification problem) . Let Ψ( 0, ) be the solution of the state equation (
0) ∩ = ∅,
then the system (
0) ∩ , then is reachable, which indicates the system is unsafe.
(a) A safe system (b) An unsafe system states 0 (in green), the reachable set ( a trajectory of the system. Panel (1a): if ( 0)(in yellow) and the unsafe region (in red). The arrow represents
0) does not intersect , then is not reachable from 0, and therefore the system is safe. Panel (1b): if there exists a trajectory that reaches from 0, then the system is unsafe.
Due to the dificulty of computing the reachable set (
0), the concept of barrier certificate (BC) [4]
has emerged as an efective way to prove the safety of a system. The idea of a BC is to find a so-called
barrier function that separates the reachable set and the unsafe set of the system, which means there
is no trajectory of the system that reaches from 0. Thus, as a suficient condition, the system is
considered safe. A BC is a function of state that satisfies a set of inequalities on both the function itself
and its time derivative along the flow of the system [ 4]. We recall that for a continuously diferentiable
function ( ), the gradient of ( ) is defined by the column vector
Definition 2 (Lie derivative). Given a vector field ( ), the Lie derivative of a continuously diferentiable
function ( ) is defined as the inner product of ∇( ) and ( ):
∇( ) =
( )
= [
( )
1
, … ,
( )
(
Input: A dynamical system ̇ = ( ), initial region 0, unsafe region and state space .
for safety verification. A BC represents a formal proof of safety for the system. For a system that is safe, the existence of multiple BCs is plausible: multiplying a BC with a positive constant will give another BC, and the set of BCs is a convex set meaning that any convex combination of BCs is again a
Gaussian process regression (GPR) is an approach for non-parametric Bayesian inference on functions over continuous domains [32]. Theoretical and practical developments over the last decade have made GPR a suitable approach for machine learning applications, in particular when training data are limited or costly. This non-parametric approach assumes that a function ∶ → ℝ is a Gaussian process (GP) meaning that its values at any ∈ is a random variable that is distributed according to a Gaussian distribution. Therefore, the full probabilistic information on the function is characterized by a mean function ( ) and a covariance function ( , ′), defined as follows:
( ) = [( )] and ( , ′) = [(( ) − ( ))(( ′) − ( ′))].
( ) ∼ ((
), ( , ′)).
When a set of observations = [ 1, 2, … , ] on the values of the function (⋅) at data points [ 1, … , ] is available, the posterior distribution conditioned on the training data (̂ |) , is again a GP with mean ̂ and variance ̂: =̂ ( ) −1
and =̂ ( , ) − ( ) −1 ( )
A frequently used covariance function is the squared exponential function with the covariance matrix ∶= [( , )] ,=1 and the vector ( ) ∶= [( , 1), … , ( , )].
( , ) = 2 exp(−|| − ||2 )
likelihood estimation.
with prior variance 2 and scaling matrix . The scaling matrix , such as diag(1/12, … , 1/ 2), assigns
individual scaling factors to each component of the vector . The vector of length scales 1, … , , and
2 are called hyperparameters, which are typically tuned through methodologies such as maximum
(
∈ 1( ) ∶= max ( , ), 2( ) ∶= max −( , ), and 3( , ) ∶= max (̇ , ), ∈Δ( ) ( ) = −∏ | 1 − 2 |.
=1 ⊆ Δ( ), For simplicity, the constraints in (14) are denoted collectively as ( , ), = 1, 2, 3 , even though 1 and 2 do not depend on .
The construction of barrier functions is often not simple. For barrier functions that are polynomials with unknown coeficients, the verification problem is transformed into obtaining the value of said coeficients. In this section, as a first step, we transform the problem into an optimization problem on the coeficients and state. Then, we propose a computational method and illustrate the algorithm.
Let ( , ) denote a polynomial with fixed structure, where is the vector of parameters in , and
= [ 1, … , ] are the polynomial variables; = [ 11, 12, … , 1 , 2 ] is a 2 vector; Δ( ) is a hyperrectangle with edges that are parallel to the parameter axes defined as
Δ( ) = { ∈ | 1 ≤ ≤ 2 , = 1, … , }.
the system’s safety on the largest possible space. Reworking the general constraints described by (
> 0 being penalty parameters that penalize the positive values of ( , ).
The functions ( , ) in (14) do not admit a closed-form solution in general since they are in the form of parameterized maximization over a continuous domain. Therefore, we move towards using Gaussian process regression that can build a likelihood for the values of such functions when a finite number of evaluations of these functions is available. Denote by ̂ ( , ) the GP model of ( , ) and define the following objective function based on the expectation of the assigned GPs: 3 =1 3 =1
3 1 2 =1
3 1 2 =1 (17) (18) (19) (20) (21) for some small admissible failure probability . To handle such chance constraints, so-called reliability-based design optimization (RBDO) [31] approaches are employed that are capable of considering constraints on the failure probability. Numerical methods for RBDO are studied previously, e.g., in [31]. The core idea is to estimate the failure probability (i.e., the complement of the right-hand side of (20)) using Monte Carlo estimation ( , ) ≈
∑ ( ̂ ( , ) > 0), Then the diference ( , ) −
where (⋅) is the indicator function and ̂ , = 1, 2, … , , are the samples drawn from the GP model ̂ ( , ).
has to be negative and is considered as part of the optimization.
The constructed optimization approach is presented in Algorithm 1. The details of the algorithm and its subroutines are given as follows with a flowchart of the algorithm presented in Figure 3. Initialization The initialization consists of: 1) selecting values for the Lagrange multipliers = [ 1, 2, 3] and the penalty factor ; 2) randomly sampling ( , ) over their domain and forming the initial dataset ; 3) computing the functions ( , ) on and storing the values in a set . where the GP model ̂ ( , ) has the mean ̂ ( , ) and variance ̂ ( , ), and we have from [30] that [ ̂
( , )] = ̂( , ), [ max(0, ̂( , ))2] = ̂ 2( , )[(1 + ̂ ( , ))2Φ( ̂ ( , ) ̂ ( , ) ̂ ( , ) ) + ( ̂ ( , ) ̂ ( , ) )].
By interpreting the values of the functions in the constraint as Gaussian random variables, the constraints can no longer be hold, but need to be interpreted probabilistically and hold with a high probability. Therefore, the constraint ( , ) ≤ 0 can be replaced with ( , , , ) = ( ) +∑ [ ̂ ( , )] +
∑ [ max(0, ̂( , ))2], Prob( ̂( , ) ≤ 0) ≥ 1 − factor when the constraints are violated in new data points. The iterations are performed until an ∗ is found that satisfies ( 16) and (13).
Loop 2 (lines 6-13) The inner loop solves the maximization in (18) on and , given the most recent candidates for and from the outer loop. The optimization (line 10) is done using Bayesian optimization, during which GP models of ( , )are learned. Then ′ is included in the dataset (line 14). With updated ′( , , , ) (line 15), the iterations are performed until an ideal ∗ is found or ′ does not improve over with current and . reasonable time constraints.
In order to verify the BC computed by the optimization algorithm, Satisfiability Modulo Theories
(SMT) solvers [33, 34] are employed for formally checking the validity of the BC. An SMT solver is a
tool for deciding correctly, i.e., without numerical errors, the satisfiability of a logical formula over a
mathematical theory, e.g., linear real arithmetic. In our case, each negated constraint can be seen as a
formula that we aim to refute (which means the constraint is satisfied). Define the logical formula
(∃ ∈ 0, −( ) < 0) ∨ (∃ ∈ , ( ) ≤ 0) ∨ (∃ ∈ , − (̇ ) < 0),
(22)
which is the negation of constraints in (
In this section, we present four case studies that demonstrate the validity of our approach. All experiments are performed with Matlab R2025a with an Apple M3, 24 GB RAM Macbook Air. The estimation The defined objective function
!(, , , )
Choose a new set of data " = (", ") Add ′ into the data set and observe c# ′
Update !" ", ", , No
Does !" ", ", , have no improvement over !?
Loop 1
Loop 2 Update ,
No
Yes
An ideal (′, ) Yes
Does c! ′ meet (13) and (16)? 1 Select initial values for the Lagrange multipliers = [ 1, 2, 3] and the penalty factor ; 2 Select a dataset by randomly choosing ( , ); 3 Compute the functions ( , ) using (14) on and store the values in a set ; 4 Compute the augmented Lagrangian ( , , , ) on the dataset ; 5 if ( , ) meet (13) and Δ( ) meets (16) for some ( , ) ∈ then 7 end 8 repeat ∗ ← ( , ) repeat 6 9 10 11 12 13 14 15 16 17 18
Select a new ′ = ( ′, ′) by maximizing (18); if ( ′) meet (13) and Δ( ′) meets (16) then ∗ ← ′ ; end ← ∪
′, ← ∪ update the Lagrangian ′
( ′, ′) ; until ≤ ′ or ∗ ≠ ∅; ← max(0, + ( ′, ′) ) ;
Set ← 2 1 if ( ′, ′) > 0 19 until ∗ ≠ ∅;
Output: ∗ = ( ∗, ∗) satisfying (13) and (16). 20 Verify that (13) is satisfied with ∗ = ( ∗, ∗) using SMT solvers; ( , , , ) with the GPs on the new dataset ; of hyperparameters for the GPR is executed using the fitrgp function in MATLAB with its default settings.
Case 4.1. Consider the two-dimensional system the BC is polynomial as ( , ) = 1 12 + 2 2 + 3. with initial set 0 = [−1, 1] × [−1, 1], unsafe set = [−1, 1] × [3, 5], and = [−5, 5] × [−5, 5] . Assume
value of the ’s are ∈ {−1, −12 , 0, 12 , 1} for = 1, 2, 3 , thereby resulting in 53 = 125 data points for the initial dataset . However, the algorithm terminates at the 27th data point, returning a valid result well component of . The ∗ column reports [ ∗, ∗] returned by Algorithm 1.
Barrier Certificate Synthesis using Algorithm 1. The column shows the initial dataset for each Row
Case Study 1 2 3 4 5
Case 1 Case 1 Case 2 {−1, −12 , 0, 1 , 1}
2 {−1, 0, 1} {−12 , 0, 1 }
2
{−1, 0, 1}
{−1, 1}
∗ = [ ∗, ∗] = [⏟⏟−⏟⏟1⏟⏟.⏟3⏟9⏟⏟9⏟3⏟⏟,⏟2⏟.⏟4⏟⏟4⏟9⏟⏟2⏟,⏟−⏟⏟⏟2⏟.⏟9⏟6⏟⏟46, ⏟−⏟⏟⏟9⏟.⏟3⏟9⏟⏟7⏟7⏟⏟,⏟5⏟.⏟2⏟⏟6⏟8⏟⏟1⏟,⏟−⏟⏟9⏟⏟.4⏟⏟5⏟9⏟⏟0⏟,⏟7⏟⏟.1⏟⏟6⏟3⏟0],
∗
∗
∗
system, and the system is safe.
which is listed in Table 1 Row 3. Therefore, the BC ( ) = −0.2465 12 − 0.6794 1 2 − 0.8787 1 −
1.1328 2 − 0.8086 satisfies the condition (13) for Δ( ) = [−4.9977, 4.8801] × [−4.9696, 4.7986]that covers
75 , 2] × [−75 , 1 ]. Hence, (
1
which means that the candidate BC ( ) = −1.3993 12 + 2.4492 2 − 2.9646 meets (13) for Δ( ) =
[−9.3977, 5.2681] × [−9.4590, 7.1630], which covers = [−5, 5] × [−5, 5] . Hence, (
Case 4.2. Consider the two-dimensional system with initial set 0 = {( 1 − 1.5)2 + 22 ≤ 0.25}, unsafe set = {( 1 + 1)2 + ( 2 + 1)2 ≤ 0.16}, and 57 , 2] × [−75 , 12 ]. Assume the BC is a polynomial as ( ) = 1 12 + 2 1 2 + 3 1 + 4 2 + 5.
Due to the conservativeness of the convex condition, the method based on the condition (
∗ = [ ∗, ∗] = [⏟0⏟.⏟7⏟⏟6⏟5⏟⏟4⏟,⏟−⏟⏟2⏟⏟.⏟7⏟8⏟5⏟⏟8⏟,⏟⏟1⏟.⏟8⏟0⏟⏟12, ⏟−⏟⏟⏟1⏟.⏟8⏟2⏟⏟1⏟0⏟⏟,⏟2⏟.⏟2⏟⏟3⏟8⏟⏟3⏟,⏟−⏟⏟0⏟⏟.0⏟⏟9⏟4⏟⏟5⏟,⏟1⏟⏟.4⏟⏟8⏟3⏟3] afirming that ( ) = 0.7654 1 − 2.7858 2 + 1.8012 fulfills condition (13) for Δ( ) = [−1.8210, 2.2383] × [−0.0945, 1.4833], which covers = [0, 1 2
] × [13 , 1]. Consequently, (
Case 4.4. The Moore-Greitzer jet engine model is described as Assume the BC is a polynomial as ( ) = 1 12 + 2 1 2 + 3 22 + 4. with initial set 0 = [0.1, 0.5] × [0.1, 0.5], unsafe set = [0.7, 1] × [0.7, 1], and = [0.1, 1] × [0.1, 1] .
∗ = [ ∗, ∗] = [⏟1⏟⏟,⏟1⏟⏟,⏟1⏟⏟,⏟−⏟ 1, 0⏟⏟.⏟1⏟,⏟⏟1⏟,⏟⏟0⏟.⏟1⏟⏟,1],
∗
∗
Therefore, the BC ( ) = 2 12 + 2 1 2 + 2 22 − 2 satisfies the condition (13) for Δ( ) = [0.1, 1] × [0.1, 1,]
which covers . Hence, (
Our Work
of 2 in Row 5. Consequently, the initial dataset comprises 8 data points instead of 27. Notably, despite the reduction in the size of the dataset, the algorithm consistently returns valid results, of course at the expense of a higher number of iterations and thus higher CPU time. In Figure 4 (bottom panels) we display the BC found for Row 5 of Table 1.
The barrier certificates reported in Table 1 are verified by the SMT solvers Z3 [ 34] and dReal [33]. Specifically, according to the applicability, for polynomial-based Cases 4.1,4.2 and 4.4, the Z3 solver [34] was employed. For Case 4.3, which has transcendental elements, the dReal solver [33] was utilized for verification. Note that the unsatisfiable
answer by dReal is formally correct.
both the efectiveness of our approach. An insightful observation from the provided table highlights that, even for the same system, distinct initializations yield diferent results. Compared with computation time in the order of hours in [37] in the case study 4.3, our method successfully computes a BC over a bounded state space, requiring only a small dataset and computation times in the order of minutes on a standard laptop.
In this paper, we proposed an optimization method for computing barrier certificates that give guarantees on safety of dynamical systems. The optimization utilizes the augmented Lagrangian framework and Gaussian process regression to eficiently represent black-box functions appearing in the constraints. The dataset needed for training the Gaussian process are sequentially generated based on Bayesian optimization and the augmented Lagrangian. The performance of the proposed optimization on the case studies shows that the method returns results with small-size samples, which are generated automatically, thus reducing the computational time from hours to minutes on the tested case studies.
Theory. Our approach currently assumes the system dynamics are known and replaces the parameterized optimizations in the barrier certificate conditions with black-box functions that can be evaluated a finite number of times. In the future, we plan to relax this assumption by treating the system as a black-box model and couple the current optimization approach directly to the data gathered from the system.
P. Zuliani is supported by the SERICS project (PE00000014) under the Italian MUR National Recovery and Resilience Plan funded by the European Union - NextGenerationEU. X. Ge is supported by the Independent Research Fund Denmark 10.46540/3120-00041B. S. Soudjani is supported by the grants EIC 101070802 and ERC 101089047.
During the preparation of this work, the authors did not use Generative AI for preparation of our work. [12] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, P. Tabuada, Control barrier functions: Theory and applications, in: European control conference, IEEE, 2019, pp. 3420–3431. [13] C. Li, Z. Zhang, N. Ahmed, Q. Liu, F. Liu, M. Buss, Safe feedback motion planning in unknown environments: An instantaneous local control barrier function approach, Journal of Intelligent & Robotic Systems 109 (2023) 40. [14] L. Dai, T. Gan, B. Xia, N. Zhan, Barrier certificates revisited, Journal of Symbolic Computation 80 (2017) 62–86. SI: Program Verification. [15] C. Sloth, G. J. Pappas, R. Wisniewski, Compositional safety analysis using barrier certificates, in:
[16] P. Jagtap, S. Soudjani, M. Zamani, Formal synthesis of stochastic systems via control barrier certificates, IEEE Transactions on Automatic Control 66 (2020) 3097–3110. [17] W. Liu, Y. Bai, L. Jiao, N. Zhan, Safety guarantee for time-delay systems with disturbances, Science
[18] M. Lewis, P. Zuliani, S. Soudjani, Verification of quantum systems using barrier certificates, in:
International Conference on Quantitative Evaluation of Systems, Springer, 2023, pp. 346–362. [19] M. Kazemi, R. Majumdar, M. Salamati, S. Soudjani, B. Wooding, Data-driven abstraction-based control synthesis, Nonlinear Analysis: Hybrid Systems 52 (2024) 101467. [20] L. Romao, A. R. Hota, A. Abate, Distributionally robust optimal and safe control of stochastic systems via kernel conditional mean embedding, in: CDC’23, 2023, pp. 2016–2021. [21] S. Bogomolov, J. Fitzgerald, S. Soudjani, P. Stankaitis, Data-driven reachability analysis of digital twin FMI models, in: International Symposium on Leveraging Applications of Formal Methods, Springer, 2022, pp. 139–158. [22] A. Perufo, D. Ahmed, A. Abate, Automated and formal synthesis of neural barrier certificates for dynamical models, in: International Conference on Tools and Algorithms for the Construction and Analysis of Systems, Springer, 2021, pp. 370–388. [23] H. Zhao, N. Qi, L. Dehbi, X. Zeng, Z. Yang, Formal synthesis of neural barrier certificates for continuous systems via counterexample guided learning, ACM Trans. Embed. Comput. Syst. 22 (2023). [24] J. Jackson, L. Laurenti, E. Frew, M. Lahijanian, Safety verification of unknown dynamical systems via Gaussian process regression, in: IEEE Conference on Decision and Control, 2020, pp. 860–866. [25] M. A. Osborne, R. Garnett, S. J. Roberts, Gaussian processes for global optimization (2009). [26] P. Jagtap, G. J. Pappas, M. Zamani, Control barrier functions for unknown nonlinear systems using
[27] J. Jiang, Y. Zhao, S. Coogan, Safe learning for uncertainty-aware planning via interval MDP abstraction, IEEE Control Systems Letters 6 (2022) 2641–2646. [28] R. Reed, L. Laurenti, M. Lahijanian, Promises of deep kernel learning for control synthesis, IEEE
Control Systems Letters 7 (2023) 3986–3991. [29] J. Stecher, L. Kiltz, K. Graichen, Semi-infinite programming using Gaussian process regression for robust design optimization, in: 2022 European Control Conference (ECC), 2022, pp. 52–59. [30] R. B. Gramacy, G. A. Gray, S. Le Digabel, H. K. Lee, P. Ranjan, G. Wells, S. M. Wild, Modeling an augmented Lagrangian for blackbox constrained optimization, Technometrics 58 (2016) 1–11. [31] C. A. Aoues Younes, Benchmark study of numerical methods for reliability-based design optimization, Structural and Multidisciplinary Optimization 41 (2010) 277–294. [32] C. E. Rasmussen, C. K. I. Williams, Gaussian Process for Machine Learning, MIT Press, 2006. [33] S. Gao, S. Kong, E. M. Clarke, Dreal: An SMT solver for nonlinear theories over the reals, in:
ADE’13, Springer-Verlag, Berlin, Heidelberg, 2013, pp. 208–214. [34] L. De Moura, N. Bjørner, Z3: An eficient SMT solver, in: TACAS’08, Springer, 2008, pp. 337–340. [35] H. Kong, F. He, X. Song, W. N. N. Hung, M. Gu, Exponential-condition-based barrier certificate generation for safety verification of hybrid systems, in: CAV, Springer, 2013, pp. 242–257. [36] H. R. Wilson, J. D. Cowan, Excitatory and inhibitory interactions in localized populations of model neurons, Biophysical Journal 12 (1972). [37] J. F. Ingham, Y. Wang, P. Zuliani, S. Soudjani, Barrier certificates for a computational model of epileptic seizures, in: IEEE Int. Conf. on Systems, Man, and Cybernetics, 2023, pp. 4728–4733. [38] B. Wooding, V. Horbanov, A. Lavaei, PRoTECT: Parallel construction of barrier certificates for safety verification of polynomial systems, in: Proceedings of the ACM/IEEE ICCPS, ACM, 2025.