StochasticBarrier.jl: A Toolbox for Stochastic Barrier Function Synthesis

StochasticBarrier .jl: A T oolbox for Stochastic Barrier Function Synthesis Rayan Mazouz rayan.mazouz@colorado.edu University of Colorado Boulder USA Frederik Baymler Mathiesen f.b.mathiesen@tudelft.nl Delft University of T e chnology The Netherlands Luca Laurenti l.laurenti@tudelft.nl Delft University of T e chnology The Netherlands Morteza Lahijanian morteza.lahijanian@colorado.edu University of Colorado Boulder USA Abstract W e present StochasticBarrier.jl , an open-source Julia-based toolbox for generating Stochastic Barrier Functions (SBFs) for safety verication of discr ete-time stochastic systems with additive Gauss- ian noise. StochasticBarrier.jl certies linear , polynomial, and piecewise ane (P W A) systems. The latter enables v erication for a wide range of system dynamics, including general nonlinear types. The toolbox implements a Sum-of-Squares (SOS) optimization ap- proach, as well as methods based on pie cewise constant (P W C) functions. For SOS-based SBFs, StochasticBarrier.jl leverages semi-denite programming solvers, while for PW C SBFs, it oers three engines: two using linear programming (LP) and one based on gradient descent (GD). Benchmarking StochasticBarrier.jl against the state-of-the-art shows that the to ol outp erforms existing tools in computation time, safety probability bounds, and scalabil- ity across over 30 case studies. Compared to its closest competitor , StochasticBarrier.jl is up to four orders of magnitude faster , achieves signicant safety probability improvements, and supports higher-dimensional systems. Ke y words Barrier Certicates, Sto chastic Systems, Probabilistic Safety , Convex Optimization, Formal V erication 1 Introduction Stochastic barrier functions (SBFs) are a class of tools used for the formal verication of stochastic systems [ 16 , 19 ]. These functions establish a lower bound on the probability of a system remaining within a safe set, using martingale theor y [ 11 , 19 ]. A key advantage of SBFs over other verication methods is their capacity to analyze the time evolution of a system through a set of static constraints, eliminating the need to propagate system uncertainty . Instead, SBF synthesis methods solve a functional optimization problem, for which various approaches have been developed based on, e.g., Sum- of-Squares (SOS) p olynomials [ 11 , 16 , 19 ] and piecewise (P W) func- tions [ 9 , 10 , 12 ]. Despite their p otential, however , there has been limited eort to develop r obust and ecient tools for SBF synthe- sis, restricting their broader application in practical elds such as robotics and autonomous vehicles. In this paper , we addr ess this gap by introducing an open-source Julia-based toolbox that generates SBFs for discrete-time stochastic systems with additive Gaussian noise - StochasticBarrier.jl 1 . The toolbox supports a wide range of dynamics, namely , linear , polynomial , and P W ane (P W A) inclusions , the latter enabling the verication of non-linear , non-polynomial systems. The estab- lished SOS optimization method is implemented, as well as more recent approaches base d on PW Constant (P WC) functions. For SOS optimization, existing semi-denite programming (SDP) solvers are leveraged. For PW C-based SBFs, three engines are included in StochasticBarrier.jl : tw o that utilize Linear Programming (LP) and one based on Gradient Descent (GD), as introduced in [10]. The toolbox is developed in Julia, a modern programming lan- guage designe d for the scientic community . Julia combines the interactivity and dynamic typing of scripting languages, facilitat- ing fast prototyping, with the high p erformance of compiled lan- guages [ 14 ]. Also, Julia’s support for parametric typing allows for customizable precision in numerical computations [ 8 ]. W e evalu- ate StochasticBarrier.jl on various benchmarks and compare its p erformance against two state-of-the-art to ols: PRoTECT [ 20 ] and StochasticBarrierFunctions [ 19 ], which, to the best of our knowledge, are the only available SBF tools for discrete-time sto- chastic systems. These existing tools e xclusively support SOS-based approaches. In addition, the three PWC-SBF methods are compared against the SOS-method. Benchmarks include over 30 case studies across eight dierent systems with linear , p olynomial, and P W A inclusion dynamics, ranging from one to six dimensions. The re- sults demonstrate that StochasticBarrier.jl outperforms ex- isting tools in three metrics: computation time, safety probabil- ity b ounds, and scalability . Compared to its closest competitor , StochasticBarrier.jl is up to 1000 × faster for successful syn- thesis, improves safety probability bounds from 0 to 1 . 0 in some cases, and handles systems with twice the dimensionality . 2 Theoretical Framework W e provide an overview of the SBF theoretical framework. 2.1 Problem Consider a stochastic process dened by dierence equation x 𝑘 + 1 = 𝑓 ( x 𝑘 ) + w 𝑘 , (1) 1 https://github.com/aria- systems- group/StochasticBarrier .jl. Mazouz et al. where x 𝑘 ∈ R 𝑛 is the state, 𝑓 : R 𝑛 → R 𝑛 is the vector eld represent- ing the one-step dynamics of System (1) , and w 𝑘 ∈ R 𝑛 is an inde- pendent and identically distributed random variable, representing noise. The probability distribution of w 𝑘 is assumed to b e zero-mean Gaussian, i.e., w 𝑘 ∼ 𝑝 w = N ( 0 , Σ ) , where Σ = diag  𝜎 2 1 , . . . , 𝜎 2 𝑛  . W e dene the stochastic kernel 𝑇 ( 𝑋 | 𝑥 ) : = ∫ R w 1 𝑋 ( 𝑓 ( 𝑥 ) + 𝑤 ) 𝑝 w ( 𝑑𝑤 ) . It follows from this denition of 𝑇 that, given an initial condi- tion x 0 = 𝑥 0 ∈ R 𝑛 , x 𝑘 is a Markov process [ 3 ]. Let 𝑋 s ⊂ R 𝑛 be a bounded set representing the safe set, 𝑋 0 ⊆ 𝑋 s be the initial set, and 𝑁 ∈ N be the time horizon. Then, probabilistic safety is dened as 𝑃 s ( 𝑋 s , 𝑋 0 , 𝑁 ) = inf 𝑥 0 ∈ 𝑋 0 Pr 𝑥 0 [∀ 𝑘 ∈ { 0 , . . ., 𝑁 } , x 𝑘 ∈ 𝑋 s ] . T o reason about 𝑃 s ( 𝑋 s , 𝑋 0 , 𝑁 ) , Stochastic Barrier Function (SBFs) can be used [ 11 , 19 ]. Function 𝐵 : R 𝑛 → R ≥ 0 is a SBF for System (1) if there exist scalars 𝜂, 𝛽 ≥ 0 , such that 𝐵 ( 𝑥 ) ≥ 1 ∀ 𝑥 ∈ R 𝑛 \ 𝑋 𝑠 , (2a) 𝐵 ( 𝑥 ) ≤ 𝜂 ∀ 𝑥 ∈ 𝑋 0 , (2b) E [ 𝐵 ( 𝑓 ( 𝑥 ) + w ) ) ] ≤ 𝐵 ( 𝑥 ) + 𝛽 ∀ 𝑥 ∈ 𝑋 𝑠 , (2c) which guarantees a lower-bound for the probabilistic safety 𝑃 s ( 𝑋 s , 𝑋 0 , 𝑁 ) ≥ 1 − ( 𝜂 + 𝑁 𝛽 ) . (3) The goal of the SBF synthesis problem is to nd a 𝐵 ( 𝑥 ) that maxi- mizes 𝑃 𝑠 , subject to Conditions (2a) - (2c) . StochasticBarrier.jl solves this problem by computing a 𝐵 along with a 𝑃 𝑠 lower bound. From a computational perspective, solving this functional opti- mization problem is non-trivial. Several methods hav e been devel- oped to ensure the convexity of the problem. One well-established approach assumes 𝐵 is a polynomial of a given degree, which allows the problem to be formulated as a SOS optimization, ultimately lead- ing to a semi-denite program [ 11 , 15 , 19 ]. A more recent approach, introduced in [ 10 ], assumes 𝐵 is a piecewise constant (P W C) func- tion, reducing the optimization problem to a linear program. Both methods are implemented in StochasticBarrier.jl . 2.2 Implemented SBF Metho ds StochasticBarrier.jl supports polynomial (SOS) and PWC class- es of SBFs for System (1) , wher e 𝑓 is ane , polynomial, or bounded by piecewise ane (PW A) functions. 2.2.1 Sum-of-Squares (SOS). A multivariate polynomial 𝜆 ( 𝑥 ) with 𝑥 ∈ R 𝑛 is an SOS polynomial if there exist polynomials 𝜆 𝑖 , 𝑖 = 1 , . . . , 𝑟 , for 𝑟 ∈ N , such that 𝜆 ( 𝑥 ) = Í 𝑟 𝑖 = 1 𝜆 2 𝑖 ( 𝑥 ) . If 𝜆 ( 𝑥 ) is SOS, then 𝜆 ( 𝑥 ) ≥ 0 for all 𝑥 ∈ R 𝑛 . The set of SOS polynomials is denote d by Λ ( 𝑥 ) , and a vector of SOS polynomials by L ( 𝑥 ) . Under the assumption that 𝑓 in System (1) is a polynomial of degree one or higher , and let 𝑋 0 , 𝑋 𝑠 be compact semi-algebraic sets, then, the following theorem [ 19 ] is implemented in StochasticBarrier.jl to reduce the SOS-SBF synthesis problem to a convex pr ogram. Theorem 2.1 (SOS-SBF for Pol ynomial Systems [ 19 ]). Let the safe set 𝑋 s = { 𝑥 ∈ R 𝑛 | ℎ s ( 𝑥 ) ≥ 0 } , initial set 𝑋 0 = { 𝑥 ∈ R 𝑛 | ℎ 0 ( 𝑥 ) ≥ 0 } , and unsafe set 𝑋 u = R 𝑛 \ 𝑋 s = { 𝑥 ∈ R 𝑛 | ℎ u ( 𝑥 ) ≥ 0 } , where ℎ 𝑠 , ℎ 0 , ℎ 𝑢 are vectors of polynomials, be given. Then, a SBF 𝐵 ( 𝑥 ) ∈ Λ ( 𝑥 ) for System (1) is obtained by solving the following SOS optimization problem min 𝛽 ,𝜂 ≥ 0 𝜂 + 𝑁 𝛽 subject to: (4a) 𝐵 ( 𝑥 ) ∈ Λ ( 𝑥 ) , (4b) − 𝐵 ( 𝑥 ) − L 𝑇 0 ( 𝑥 ) ℎ 0 ( 𝑥 ) + 𝜂 ∈ Λ ( 𝑥 ) , (4c) 𝐵 ( 𝑥 ) − L 𝑇 u ( 𝑥 ) ℎ u ( 𝑥 ) − 1 ∈ Λ ( 𝑥 ) , (4d) − E [ 𝐵 ( 𝑓 ( 𝑥 ) + w ) ] + 𝐵 ( 𝑥 ) + 𝛽 − L 𝑇 𝑠 ( 𝑥 ) ℎ 𝑠 ( 𝑥 ) ∈ Λ ( 𝑥 ) , (4e) where L 0 , L ℎ , and L 𝑠 are vectors of SOS polynomials of appropriate dimension. A valid solution to this optimization guarantees safety probability in Eq. (3) . PW A Inclusion Dynamics. Additionally to polynomial 𝑓 , non- polynomial 𝑓 are also supporte d in StochasticBarrier.jl by supplying piecewise ane (P W A) bounds of 𝑓 . In particular , given a partition 𝑄 = { 𝑞 1 , . . . , 𝑞 | 𝑄 | } of the safe set 𝑋 s for every 𝑞 ∈ 𝑄 , StochasticBarrier.jl supports 𝑓 expressed via bounds 𝑓 𝑞 ( 𝑥 ) = 𝐴 𝑞 𝑥 + 𝑏 𝑞 and 𝑓 𝑞 ( 𝑥 ) = 𝐴 𝑞 𝑥 + 𝑏 𝑞 (5) such that 𝑓 𝑞 ( 𝑥 ) ≤ 𝑓 ( 𝑥 ) ≤ 𝑓 𝑞 ( 𝑥 ) for every 𝑥 ∈ 𝑞 . The following theorem shows a convex optimization framework for synthesis of SOS-SBF for such dynamics. Theorem 2.2 (SOS-SBF for P W A Inclusion Systems [ 11 ]). Con- sider function 𝐵 ( 𝑥 ) and sets 𝑋 0 , 𝑋 𝑠 , 𝑋 𝑢 as in Theorem 2.1 and 𝑋 𝑠 - partition 𝑄 = { 𝑞 1 , . . . , 𝑞 | 𝑄 | } , where each 𝑞 ∈ 𝑄 is a semi-algebraic set 𝑞 = { 𝑥 ∈ R 𝑛 | ℎ 𝑞 ( 𝑥 ) ≥ 0 } . Given P W A functions 𝑓 𝑞 and 𝑓 𝑞 over 𝑄 as dened in Eq. (5) , an SOS SBF 𝐵 ( 𝑥 ) for System (1) with general 𝑓 ( 𝑥 ) ∈ [ 𝑓 𝑞 ( 𝑥 ) , 𝑓 𝑞 ( 𝑥 ) ] can b e obtained by solving the SOS optimiza- tion problem in (4) and replacing Constraint (4e) with − E [ 𝐵 ( 𝑦 + 𝑤 ) | 𝑥 ] + 𝐵 ( 𝑥 ) + 𝛽 − L 𝑇 𝑞,𝑥 ( 𝑥 ) ℎ 𝑞 ( 𝑥 ) − L 𝑇 𝑞, 𝑦 ( 𝑥 )  ( 𝑓 𝑞 ( 𝑥 ) − 𝑦 ) ⊙ ( 𝑦 − 𝑓 𝑞 ( 𝑥 ) )  ∈ Λ ( 𝑥 , 𝑦 ) , for all 𝑞 ∈ 𝑄 , where ⊙ is the Schur product. In StochasticBarrier.jl , the implementation of partition 𝑄 is based on a grid, wher e each 𝑞 ∈ 𝑄 is dened as a hyper-rectangle. Then, provided partition 𝑄 and PW A functions 𝑓 𝑞 and 𝑓 𝑞 , the tool performs the above optimization. 2.2.2 PWC Method. StochasticBarrier.jl also supports the class of P WC-SBFs. Consider a partition 𝑄 = { 𝑞 1 , . . . , 𝑞 | 𝑄 | } of 𝑋 𝑠 such that vector eld 𝑓 ( 𝑥 ) is continuous in each region 𝑞 ∈ 𝑄 . Dene PWC function 𝐵 : R 𝑛 → R as 𝐵 ( 𝑥 ) = 𝑏 𝑖 if 𝑥 ∈ 𝑞 𝑖 , other wise 𝐵 ( 𝑥 ) = 1 . (6) where 𝑏 𝑖 ∈ R . Further , let 𝑞 𝑢 = 𝑋 𝑢 and denote the bounds on the transition kernel 𝑇 from partition 𝑞 𝑖 ∈ 𝑄 by , for all 𝑗 ∈ { 1 , . . . , | 𝑄 | , 𝑢 } , 𝑝 𝑖 𝑗 = inf 𝑥 ∈ 𝑞 𝑖 𝑇 ( 𝑞 𝑗 | 𝑥 ) , and 𝑝 𝑖 𝑗 = sup 𝑥 ∈ 𝑞 𝑖 𝑇 ( 𝑞 𝑗 | 𝑥 ) . (7) Then, the set of feasible values for 𝑇 from 𝑞 𝑖 is given by P 𝑖 =  𝑝 𝑖 = ( 𝑝 𝑖 1 , . . . , 𝑝 𝑖 | 𝑄 | , 𝑝 𝑖 u ) ∈ [ 0 , 1 ] | 𝑄 | + 1 s.t. Í | 𝑄 | 𝑗 = 1 𝑝 𝑖 𝑗 + 𝑝 𝑖 u = 1 , 𝑝 𝑖 𝑗 ≤ 𝑝 𝑖 𝑗 ≤ 𝑝 𝑖 𝑗 ∀ 𝑗 ∈ { 1 , . . . , | 𝑄 | , u }  . The following theorem formalizes the optimization problem for PWC class of SBFs for System (1). Theorem 2.3 (PWC-SBF [ 10 ]). Consider System (1) with safe set 𝑋 𝑠 , initial set 𝑋 0 , and partition 𝑄 = { 𝑞 1 , . . . , 𝑞 | 𝑄 | } of 𝑋 𝑠 . Let B 𝑄 be the class of P WC functions in the form of Eq. (6) and P = P 1 × . . . × P | 𝑄 | , StochasticBarrier .jl: A T oolb ox for Stochastic Barrier Function Synthesis T able 1: Overview of SBF synthesis methods implemented in StochasticBarrier.jl for dierent classes of dynamics. Dynamics 𝑓 ( 𝑥 ) SBF 𝐵 ( 𝑥 ) Solver Method Linear SOS SDP PWC Dual LP, CEGIS LP , GD Polynomial SOS SDP PW A Inclusion SOS SDP (Eq. (5)) PWC Dual LP, CEGIS LP , GD where each P 𝑖 is the set of feasible transitions. Then, 𝐵 ∗ ∈ B 𝑄 is a PWC-SBF if 𝐵 ∗ is a solution to 𝐵 ∗ = arg min 𝐵 ∈ B 𝑄 max ( 𝑝 𝑖 ) | 𝑄 | 𝑖 = 1 ∈ P 𝜂 + 𝑁 𝛽 subject to, ∀ 𝑖 ∈ { 1 , . . . , 𝑄 } , 𝑏 𝑖 ≥ 0 , 𝑄  𝑗 = 1 𝑏 𝑗 · 𝑝 𝑖 𝑗 + 𝑝 𝑖 u ≤ 𝑏 𝑖 + 𝛽 𝑖 , 0 ≤ 𝛽 𝑖 ≤ 𝛽 , and for all 𝑖 : 𝑋 𝑖 ∩ 𝑋 0 ≠ ∅ , 𝑏 𝑖 ≤ 𝜂 . Then, 𝐵 ∗ guarantees the safety probability in Eq. (3) . W ork [ 10 ] intr oduces three methods for synthesizing P WC SBFs: (i) an exact duality-based approach formulate d as a zero-duality-gap Linear Program (LP), (ii) an e xact Counter-Example Guided Induc- tive Synthesis (CEGIS) method that alternates between generating candidate barriers and counter-examples via two LPs, and (iii) a gradient descent (GD) metho d oering improved scalability and memory eciency . StochasticBarrier.jl implements all three methods for System (1) with linear and P W A inclusion dynamics. 3 Overview of StochasticBarrier .jl StochasticBarrier.jl is a Julia toolbox for the synthesis of SBFs for System (1) . The tool includes multiple indep endent synthesis methods with the following features: • Synthesis of SOS-SBFs for systems in the form of (1) with linear , polynomial and P W A inclusion dynamics, • Computation of transition kernel bounds for System (1) with linear and PW A inclusion dynamics, • Synthesis of P WC-SBFs through LP (dual and CEGIS) and GD optimization methods given the kernel bounds. In T able 1, the SBF methods are outlined for dierent classes of dynamics 𝑓 . The solver method refers to the convex optimization class. T o instantiate an instance of SBF synthesis problem, the sys- tem model nee ds to be dened rst, which consists of the dynamics and the relevant sets (safe, unsafe and initial). In Section 3.1, we show how to setup these system related parameters. Then, the SBF synthesis method needs to be spe cied, as describe d in Section 3.2. 3.1 System Setup StochasticBarrier.jl dispatches and checks compatibility with the chosen synthesis algorithm based on the system type. W e give examples about how to construct a linear , polynomial, and P W A inclusion system, and specify the safe, initial, and unsafe sets. The examples are for 2D-systems, but the user can easily set up an 𝑛 -dimensional problem. 3.1.1 Linear D ynamics Example. With linear dynamics, both SOS and PWC classes of SBFs can be used. SOS-SBF. T o synthesize a SOS-SBF for System (1) with linear vector eld 𝑓 ( 𝑥 ) = 𝐴𝑥 + 𝑏 , we rst ne ed to specify matrix 𝐴 and bias vector 𝑏 . The noise distribution is assumed to b e zero-mean Gaussian 𝑝 w = N ( 0 , 𝜎 2 𝐼 ) characterized by standard deviation 𝜎 . An example setup is shown below for a 2D-linear system. # S y s t e m L i n e a r : S O S A = [ 0 . 9 5 0 . 0 ; 0 . 0 0 . 9 5 ] , b = [ 0 . 0 , 0 . 0 ] , 𝜎 = [ 0 . 1 , 0 . 1 ] s p a c e = H y p e r r e c t a n g l e ( l o w = [ - 1 . 0 , - 2 . 0 ] , h i g h = [ 1 . 0 , 2 . 0 ] ) s y s t e m = A d d i t i v e G a u s s i a n L i n e a r S y s t e m ( A , b , 𝜎 , s p a c e ) The Hyperrectangle struct of LazySets.jl , a Julia package for calculus with convex sets [ 4 ], is used to dene the state space. A struct of AdditiveGaussianLinearSystem is passed as a system object into StochasticBarrier.jl . PWC-SBF . Using P W C-SBF requires an additional parameter to dene a grid partitioning of the state space. The to ol uses this grid to dene the domain of P WC function in Eq. (6) and com- pute the transition kernel bounds in Eq. (7) . The grid is parameter- ized by 𝜖 ∈ R 𝑛 , where the 𝑖 -th entry is the grid cell’s half-width in the 𝑖 -th dimension. The generate_partitions() function takes the state space and 𝜖 as inputs, and returns the partitions. The system object AdditiveGaussianLinearSystem , along with the state partitions, are parsed into transition_probabilities() , which is a function in StochasticBarrier.jl that computes the transition kernel bounds. An example 2D-linear system setup is provided. # S y s t e m L i n e a r : P W C A = [ 0 . 5 0 0 . 0 ; 0 . 0 0 . 5 0 ] , b = [ 0 . 0 , 0 . 0 ] , 𝜎 = [ 0 . 1 , 0 . 1 ] s p a c e = H y p e r r e c t a n g l e ( l o w = [ - 1 . 0 , - 1 . 0 ] , h i g h = [ 0 . 5 , 0 . 5 ] ) # H a l f - w i d t h v e c t o r 𝜖 = [ 0 . 2 , 0 . 2 ] # C o m p u t e p r o b a b i l i t y d a t a s y s t e m = A d d i t i v e G a u s s i a n L i n e a r S y s t e m ( A , b , 𝜎 ) p a r t i t i o n s = g e n e r a t e _ p a r t i t i o n s ( s p a c e , 𝜖 ) p r o b a b i l i t i e s = t r a n s i t i o n _ p r o b a b i l i t i e s ( s y s t e m , p a r t i t i o n s ) In addition, the toolbox supports saving and loading transition probability bounds to the Network Common Data Form (NetCDF) [18]. See Appendix C for more details on this feature. 3.1.2 Polynomial D ynamics Example. For p olynomial systems, StochasticBarrier.jl only supp orts the SOS class of SBFs. T o set up, rst dene a vector of symbol variables 𝑥 , using the @polyvar macro of DynamicPolynomials.jl [ 6 ]. Then, the polynomial func- tion 𝑓 is dened symb olically . The programmatic construct of AdditiveGaussianPolySystem is passed as an object to the barrier synthesis method. An example is provided below for a 2D system. # S y s t e m P o l y n o m i a l : S O S @ p o l y v a r x [ 1 : 2 ] # d i m = 2 𝜏 = 0 . 1 f = [ x [ 1 ] + 𝜏 * x [ 2 ] , x [ 2 ] + 𝜏 * ( - x [ 1 ] + ( 1 - x [ 1 ] ) ^ 2 ) ] 𝜎 = [ 0 . 0 2 , 0 . 0 2 ] s p a c e = H y p e r r e c t a n g l e ( l o w = [ - 6 . 0 , - 6 . 0 ] , h i g h = [ 6 . 0 , 6 . 0 ] ) s y s t e m = A d d i t i v e G a u s s i a n P o l y S y s t e m ( f , 𝜎 , s p a c e ) 3.1.3 PW A Inclusion Dynamics Example. Setting up the system model for PW A inclusion dynamics requires PW A bounds on dy- namics 𝑓 of System (1) , as dened in Eq. (5) . These bounds can be generated using Linear Bound Propagation [ 21 ], and stored in an Mazouz et al. NetCDF le. These les can be loaded into the tool, which oers the options of using SOS- or PWC-SBF methods. SOS-SBF. Given a NetCDF le that contains the PW A bounds, the setup for the SOS-SBF option is as follows. First, the system ag and lename need to b e dened, which together specify the path to the PW A b ounds. Then, the load_dynamics() function r eads this data. The bounds dataset, along with the standard deviation 𝜎 , are passed into the AdditiveGaussianUncertainPWASystem struct. See the example below for a 2D pendulum model. # S y s t e m s y s t e m _ f l a g = " p e n d u l u m " f i l e n a m e = " p a t h / t o / $ s y s t e m _ f l a g / d a t a . n c " d a t a s e t = o p e n _ d a t a s e t ( f i l e n a m e ) d a t a = l o a d _ d y n a m i c s ( d a t a s e t ) 𝜎 = [ 0 . 0 1 , 0 . 0 1 ] s y s t e m = A d d i t i v e G a u s s i a n U n c e r t a i n P W A S y s t e m ( d a t a , 𝜎 ) PWC-SBF . For the P W C-SBF option, the le and system sp eci- cation are similar . In addition, the transition probability data must be computed using transition_probabilities() . # C o m p u t e p r o b a b i l i t y d a t a p r o b a b i l i t i e s = t r a n s i t i o n _ p r o b a b i l i t i e s ( s y s t e m ) 3.1.4 Specifying Initial and Unsafe Sets. T o complete the system model setup, the initial set and unsafe region(s) need to b e dened. This is done by using the Hyperrectangle struct. In the to olbox, the region R 𝑛 \ 𝑋 is considered unsafe by default. The user can further declare unsafe regions in the domain. If more than one unsafe region is dened, the UnionSet function of LazySets.jl can be used to take their union, as shown in the example below . i n i t i a l = H y p e r r e c t a n g l e ( l o w = [ - 0 . 2 , - 0 . 3 ] , h i g h = [ 0 . 0 , 0 . 1 ] ) u n s a f e 1 = H y p e r r e c t a n g l e ( l o w = [ - 0 . 5 5 , 0 . 4 ] , h i g h = [ - 0 . 4 5 , 0 . 6 ] ) u n s a f e 2 = H y p e r r e c t a n g l e ( l o w = [ 0 . 4 5 , 0 . 4 ] , h i g h = [ 0 . 5 5 , 0 . 6 ] ) u n s a f e _ r e g i o n s = U n i o n S e t ( u n s a f e 1 , u n s a f e 2 ) If domain 𝑋 does not contain any unsafe regions, the variable unsafe_regions must be declared empty: EmptySet(dim) , where dim is the dimension of the system. With that, the system model setup is concluded, and one can call on various SBF synthesis meth- ods, as explained below . 3.2 SBF Synthesis All methods for synthesizing SBFs are calle d under the function synthesize_barrier(alg, system, initial_region, unsafe_region; time_horizon=N) , which is designed with a unied interface of in- put parameters. Julia’s multiple dispatch will select the appr opriate method based on the alg parameter , and each algorithm returns a result obje ct, which contains, in addition to the barrier , the resulting safety certicate and computation time. The barrier type is deter- mined by the algorithm. Sp ecifying the time horizon is optional, with the default value set to 𝑁 = 1 . 3.2.1 SOS-SBF. T o synthesize an SOS-SBF, synthesize_barrier is called with input SumOfSquaresAlgorithm() as the rst parameter . The second parameter , the system model object dene d in Section 3.1, can be any of the three system types. N = 1 0 # O p t i m i z e : b a s e l i n e 1 ( s o s a p p r o a c h ) r e s = s y n t h e s i z e _ b a r r i e r ( S u m O f S q u a r e s A l g o r i t h m ( ) , s y s t e m , i n i t i a l _ r e g i o n , u n s a f e _ r e g i o n ; t i m e _ h o r i z o n = N ) # R e s u l t b a r r = b a r r i e r ( r e s ) 𝜂 , 𝛽 , p s = e t a ( r e s ) , b e t a ( r e s ) , p s a f e ( r e s , N ) c o m p u t a t i o n _ t i m e = t o t a l _ t i m e ( r e s ) The optimization returns a struct of values, which include the SOS-SBF 𝐵 , 𝜂 , and 𝛽 values, and the computation time. The al- gorithm SumOfSquaresAlgorithm() accepts keyword arguments to congure, e.g. the underlying SDP solver and the degree of the SOS barrier . Below we list the default values and ho w to construct the object for customization. a l g = S u m O f S q u a r e s A l g o r i t h m ( ; b a r r i e r _ d e g r e e = 4 l a g r a n g e _ d e g r e e = 2 s d p _ s o l v e r = M o s e k . O p t i m i z e r ) Setting SumOfSquaresAlgorithm(barrier_degree = 30) , for exam- ple, overrides the barrier degree. The default SDP solver for SOS optimization is Mosek [2]. 3.2.2 PWC-SBF . For PWC-SBFs, synthesize_barrier is calle d with DualAlgorithm() for the dual LP approach, CEGISAlgorithm() for the CEGIS approach, and GradientDescentAlgorithm() for GD ap- proach, respectively , as the rst parameter . The second parameter should be the result from transition_probabilities . N = 1 0 # O p t i m i z e : m e t h o d 1 ( d u a l a p p r o a c h ) d u a l = s y n t h e s i z e _ b a r r i e r ( D u a l A l g o r i t h m ( ) , p r o b a b i l i t i e s , i n i t i a l _ r e g i o n , o b s t a c l e _ r e g i o n , t i m e _ h o r i z o n = N ) # O p t i m i z e : m e t h o d 2 ( C E G I S a p p r o a c h ) c e g i s = s y n t h e s i z e _ b a r r i e r ( C E G I S A l g o r i t h m ( ) , p r o b a b i l i t i e s , i n i t i a l _ r e g i o n , o b s t a c l e _ r e g i o n ; t i m e _ h o r i z o n = N ) # O p t i m i z e : m e t h o d 3 ( G D a p p r o a c h ) g d = s y n t h e s i z e _ b a r r i e r ( G r a d i e n t D e s c e n t A l g o r i t h m ( ) , p r o b a b i l i t i e s , i n i t i a l _ r e g i o n , o b s t a c l e _ r e g i o n ; t i m e _ h o r i z o n = N ) Similar to SOS, each algorithm returns a result struct, which, in addition to the P WC-SBF 𝐵 , values for 𝜂 , 𝛽 and the computation time, also includes 𝛽 𝑖 , which is an upper b ound on Eq. (2c) in each region 𝑖 . Next, we highlight the core properties of each algorithm. Linear Programming. As outlined in Section 2.2.2, amongst the LP classes, there is the dual LP method, and the CEGIS method which uses inner and outer LP optimizations. The dual LP method only has the solver as a congurable parameter . a l g = D u a l A l g o r i t h m ( ; l i n e a r _ s o l v e r = H i G H S . O p t i m i z e r ) The default LP solver across the toolb ox is HiGHS solver [ 5 ], but is customizable to any LP solver supporte d by JuMP, Julia Mathemat- ical Programming [7]. The CEGISAlgorithm in addition be congured with the number of iterations and whether a l g = C E G I S A l g o r i t h m ( ; l i n e a r _ s o l v e r = H i G H S . O p t i m i z e r n u m _ i t e r a t i o n s = 1 0 d i s t r i b u t i o n _ g u i d e d = t r u e ) The default is the distribution guided approach, wher e in each iteration one LP generates a candidate P WC-SBF that optimizes the unsafety probability given a set of nite feasible distributions. Then, the second LP generates distribution witnesses that violate StochasticBarrier .jl: A T oolb ox for Stochastic Barrier Function Synthesis the safety probability guarantee of the candidate P WC-SBF . By de- fault, it runs for a xed number of iterations given by num_iterations (default = 10), or it can continue adaptively until no further coun- terexamples exist, guaranteeing nite-time convergence [10]. Projected Gradient Descent. The GD method solves the same problem as the LP approaches, utilizing the minimax formulation to dene a convex loss function. The projecte d gradient descent method is subject to hyperparameters. a l g = G r a d i e n t D e s c e n t A l g o r i t h m ( ; n u m _ i t e r a t i o n s = 1 0 0 0 0 i n i t i a l _ l r = 1 e - 2 d e c a y = 0 . 9 9 9 9 m o m e n t u m = 0 . 9 ) The struct GradientDescentAlgorithm denes key parameters for running GD. The num_iterations parameter controls the maximum number of iterations (default: 10,000). The initial_lr sets the ini- tial learning rate for the algorithm (default: 0.01), determining the size of each update step. The decay parameter (default: 0.9999) r e- duces the learning rate over time to ne-tune convergence. Finally , momentum (default: 0.9) helps to accelerate convergence by incor- porating a p ortion of the previous gradient step into the current. Importantly , at termination, the values for 𝜂 , 𝛽 are guarantee d to be sound with respect to the returned barrier , although not absolutely optimal in the sense of Theorem 2.3. StochasticBarrier.jl is also equipped with an independent interface to easily run benchmarks. See Appendix A for details. 4 Benchmark Studies In this section, we present benchmark results for various stochastic systems with linear , polynomial and P W A inclusion dynamics. W e performed two main comparison studies. First, the SOS-SBF algo- rithm in StochasticBarrier.jl is compared against the state-of- the-art SOS toolb oxes: StochasticBarrierFunctions [ 19 ] (MA T - LAB) and PRoTECT [ 20 ] (Python). The second evaluates all the StochasticBarrier.jl methods (SOS, DU AL-LP, CEGIS-LP , and GD) and compares their performance. All experiments are per- formed on a computer with a 3.9 GHz 8-Core CP U and 128 GB of memory . All methods use Mosek as the underlying SDP solver . The systems considered are summarized in T able 2 and details of the systems and their setups are provided in App endix B. For all benchmarks, time horizon 𝑁 = 10 is used. The tool is available at https://github.com/aria- systems- group/StochasticBarrier .jl. T able 2: Linear , polynomial and PW A inclusion stochastic dy- namic systems for benchmarking in StochasticBarrier.jl . Dynamics Name Dimension Linear Contraction Map 1 & 2 2 T wo T ank [17] 2 Room T emperature [20] 3 Quadrotor [10] 6 Polynomial Thermostat [13] 1 V an der Pol Oscillator [1] 2 PW A Inclusion Pendulum NNDM [11] 2 Unicycle [10] 4 Comparison Against Other SOS-SBF T ools. Fig. 1 shows StochasticBarrier.jl ’s performance of the SOS-SBF algorithm against StochasticBarrierFunc. [ 19 ] and PRoTECT [ 20 ]. A de- tailed table of results can b e found in App endix D. Note that as StochasticBarrierFunc. only supports linear systems, and PRo- TECT only linear and polynomial systems, we could only compare the tools with systems of this dynamics class. W e also note that the results from PRoTECT are random, and we pr esent average results over 10 runs. W e also note that 𝑃 𝑠 dened in Eq. (3) is computed dierently in PRoTECT ; we adjust it here to match our setting. This paper solely focuses on optimization-based SOS methods, i.e., dis- abling the objective is fundamentally dierent and not considered. Across all benchmarks, StochasticBarrier.jl outperforms both StochasticBarrierFunc. and PRoTECT in terms of computa- tion time and scalability , and tightness of safety probability . Stocha- sticBarrierFunc. fails to produce non-trivial results and cannot handle higher-degree SOS-SBFs (b eyond degree 12), often yielding a trivial 𝑃 𝑠 ≥ 0 or exceeding available memor y . PRoTECT performs better but remains signicantly slower ( by factors ranging from 5 × to 1000 × ) and often returns trivial safety bounds for complex systems. StochasticBarrier.jl , in contrast, scales eciently to high-degree SOS-SBFs (up to 30), achieves non-trivial or high safety probabilities, and maintains orders-of-magnitude faster computa- tion across all tested models (contraction map 1, two tank, room temperature, quadrotor , thermostat, and oscillator). Comparisons of StochasticBarrier.jl engines. Bench- mark results for the three P W C-SBF methods implemented in StochasticBarrier.jl are summarized in Fig.2, with detailed results presented in T able 4 of Appendix D. The results show that PWC-SBFs consistently outperform SOS-SBFs, particularly for the contraction map and hybrid nonlinear systems. For the linear con- traction map, the GD metho d achieves the fastest computation, especially as the number of partitions increases, while Dual LP and CEGIS yield tighter safety probability bounds that converge with higher partition counts. For P W A inclusion systems, PWC-SBFs are up to two orders of magnitude faster , and for the unicycle model, they provide non-trivial safety bounds where SOS fails. As expected, the Dual LP and CEGIS methods deliver the tightest b ounds at the cost of higher computation time. 5 Conclusion W e introduced StochasticBarrier.jl , a Julia toolbox for generat- ing barrier certicates and estimating safety probabilities of discrete- time stochastic dynamical systems. Benchmark results show that StochasticBarrier.jl is orders of magnitude faster and pr ovides tighter safety bounds than existing MA TLAB and Python implemen- tations, with sup erior scalability to higher-dimensional systems and higher-degree SOS barriers. The toolbox also supports piecewise constant barriers, which consistently outperform SOS formulations. Current limitations include its focus on additive Gaussian noise. Fu- ture extensions will address non-Gaussian and non-additive noise. References [1] Aba te, A., Blom, H. A., Ca uchi, N., Delicaris, J., Hartmanns, A., Khaled, M., La v aei, A., Pilch, C., Remke, A., Schupp, S., et al. Arch-comp20 categor y report: Stochastic models. In ARCH (2020). [2] ApS, M. Mosek optimization toolb ox for matlab. User’s Guide and Reference Manual, V ersion 4 , 1 (2019). Mazouz et al. 0 5 10 15 20 25 30 10 -2 10 0 10 2 10 4 (s) 0 5 10 15 20 25 30 Barrier Degree 0 0.2 0.4 0.6 0.8 1 P s StochasticBarrierFunc. PRoTECT StochasticBarrier.jl (a) Contraction Map 4 4.5 5 5.5 6 6.5 7 7.5 8 10 -1 10 0 10 1 10 2 10 3 (s) 4 4.5 5 5.5 6 6.5 7 7.5 8 Barrier Degree 0 0.2 0.4 0.6 0.8 1 P s StochasticBarrierFunc. PRoTECT StochasticBarrier.jl (b) T wo T ank 2 3 4 5 6 7 8 9 10 11 12 10 -2 10 0 10 2 10 4 (s) 2 3 4 5 6 7 8 9 10 11 12 Barrier Degree 0 0.2 0.4 0.6 0.8 1 P s StochasticBarrierFunc. PRoTECT StochasticBarrier.jl (c) Room T emp erature 2 2.5 3 3.5 4 4.5 5 5.5 6 10 -2 10 0 10 2 10 4 (s) 2 2.5 3 3.5 4 4.5 5 5.5 6 Barrier Degree 0 0.2 0.4 0.6 0.8 1 P s PRoTECT StochasticBarrier.jl (d) Quadrotor 4 4.5 5 5.5 6 6.5 7 7.5 8 10 -1 10 0 (s) 4 4.5 5 5.5 6 6.5 7 7.5 8 Barrier Degree 0 0.2 0.4 0.6 0.8 1 P s StochasticBarrier.jl PRoTECT (e) Thermostat 6 7 8 9 10 11 12 10 0 10 2 10 4 (s) 6 7 8 9 10 11 12 Barrier Degree 0 0.2 0.4 0.6 0.8 1 P s PRoTECT StochasticBarrier.jl (f ) Oscillator Figure 1: Polynomial systems benchmarking of StochasticBarrier.jl vs StochasticBarrierFunctions [19] and PRo TECT [20]. Plots depict optimization time 𝜏 ( 𝑠 ) and probability of safety 𝑃 𝑠 . If data points are not plotted beyond a certain barrier degree, it means the algorithm ran out of memory . 50 100 150 200 250 300 350 400 10 0 10 2 10 4 (s) 50 100 150 200 250 300 350 400 |Q| 0.2 0.4 0.6 0.8 1 P s SOS PWC Dual LP PWC CEGIS PWC GD (a) Contraction Map 100 150 200 250 300 350 400 450 500 10 -1 10 0 10 1 10 2 (s) 100 150 200 250 300 350 400 450 500 |Q| 0.88 0.9 0.92 0.94 0.96 0.98 P s SOS PWC Dual LP PWC CEGIS PWC GD (b) Pendulum 1200 1300 1400 1500 1600 1700 1800 10 1 10 2 10 3 (s) 1200 1300 1400 1500 1600 1700 1800 |Q| 0 0.2 0.4 0.6 0.8 P s SOS PWC Dual LP PWC CEGIS PWC GD (c) Unicycle Figure 2: Benchmark plots for the three P WC-SBF metho ds. For the pendulum and unicycle, the SOS results are also depicted for comparison. Plots depict total optimization time 𝜏 ( 𝑠 ) and probability of safety 𝑃 𝑠 . [3] Bertsekas, D. P., and Shreve, S. E. Stochastic Optimal Control: the Discrete-time Case , vol. 5. Athena Scientic, 2004. [4] Forets, M., and Schilling, C. LazySets.jl: Scalable Symbolic-Numeric Set Computations. Proceedings of the JuliaCon Conferences 1 , 1 (2021), 11. [5] Huangfu, Q., and Hall, J. J. Parallelizing the dual revised simplex method. Mathematical Programming Computation 10 , 1 (2018), 119–142. [6] Lega t, B. Multivariate polynomials in Julia. In JuliaCon ( July 2022). [7] Lubin, M., Dowson, O., Dias Garcia, J., Huchette, J., Lega t, B., and Vielma, J. P. JuMP 1.0: Recent improvements to a modeling language for mathematical optimization. Mathematical Programming Computation (2023). [8] Ma thiesen, F. B., Lahijanian, M., and La urenti, L. IntervalMDP. jl: Accel- erated Value Iteration for Interval Markov Decision Processes. arXiv preprint arXiv:2401.04068 (2024). [9] Ma thiesen, F. B., Romao, L., Calvert, S. C., La urenti, L., and Aba te, A. A data-driven approach for safety quantication of non-linear stochastic systems with unknown additive noise distribution. arXiv preprint arXiv:2410.06662 (2024). [10] Mazouz, R., Baymler Mathiesen, F., Laurenti, L., and Lahijanian, M. Pie ce- wise Stochastic Barrier Functions. arXiv preprint arXiv:2404.16986 (2024). [11] Mazouz, R., Muvv ala, K., Ratheesh, A., Laurenti, L., and Lahijanian, M. Safety Guarantees for Neural Network Dynamic Systems via Stochastic Barrier Functions. Advances in Neural Information Processing Systems (2022). [12] Mazouz, R., Skovbekk, J., Mathiesen, F. B., Frew, E., Laurenti, L., and Lahija- nian, M. Data-Driven Permissible Safe Control with Barrier Certicates. arXiv preprint arXiv:2404.16986 (2024). StochasticBarrier .jl: A T oolb ox for Stochastic Barrier Function Synthesis [13] Neja ti, A., Soudjani, S., and Zamani, M. Compositional abstraction-based synthesis for continuous-time stochastic hybrid systems. European Journal of Control 57 (2021), 82–94. [14] Perkel, J. M., et al. Julia: come for the syntax, stay for the speed. Nature 572 , 7767 (2019), 141–142. [15] Prajna, S. Barrier certicates for nonlinear model validation. Automatica 42 , 1 (2006), 117–126. [16] Prajna, S., Jadbabaie, A., and P appas, G. J. A framework for worst-case and stochastic safety verication using barrier certicates. IEEE Transactions on Automatic Control 52 , 8 (2007). [17] Ramos, J. A., and Dos Santos, P. L. Mathematical modeling, system identication, and controller design of a two tank system. In 2007 46th IEEE Conference on Decision and Control (2007), IEEE, pp. 2838–2843. [18] Rew, R., and Da vis, G. Netcdf: an interface for scientic data access. IEEE computer graphics and applications 10 , 4 (1990). [19] Santoy o, C., Dutreix, M., and Coogan, S. A barrier function approach to nite-time stochastic system verication and control. A utomatica 125 (2021). [20] W ooding, B., Horbanov, V ., and La v aei, A. Protect: Parallelized construction of safety barrier certicates for nonlinear polynomial systems. arXiv preprint arXiv:2404.14804 (2024). [21] Zhang, H., Weng, T .-W ., Chen, P.-Y ., Hsieh, C.-J., and Daniel, L. Ecient neural network robustness certication with general activation functions. Advances in neural information processing systems 31 (2018). A Benchmarking Setup StochasticBarrier.jl is also equipped with an independent in- terface to easily run benchmarks. Fig. 3 depicts the general structure of the interface. The input is a .yaml le, in which the dynamics, the sets (safe, unsafe and initial), and the SBF synthesis method are specied. The output includes the SBF 𝐵 , values for 𝜂 and 𝛽 , the lower bound on safety probability 𝑃 𝑠 , and the optimization time. Below , we show the .yaml input les for the systems discussed above for SOS-SBF and PWC-SBF . s y s t e m _ f l a g : " l i n e a r " d i m : 2 A : [ [ 0 . 9 5 , 0 . 0 ] , [ 0 . 0 , 0 . 9 5 ] ] b : [ 0 . 0 , 0 . 0 ] 𝜎 : [ 0 . 0 1 , 0 . 0 1 ] s t a t e _ s p a c e : l o w : [ - 1 . 0 , - 2 . 0 ] h i g h : [ 1 . 0 , 2 . 0 ] b a r r i e r _ s e t t i n g s : b a r r i e r _ t y p e : " S O S " t i m e _ h o r i z o n : 1 0 b a r r i e r _ d e g r e e : 1 2 s y s t e m _ f l a g : " p o l y n o m i a l " d i m : 2 f : [ " 0 . 1 * x [ 2 ] + x [ 1 ] " , " 0 . 1 + x [ 2 ] - 0 . 3 * x [ 1 ] + 0 . 1 * x [ 1 ] ^ 2 " ] 𝜎 : [ 0 . 0 2 , 0 . 0 2 ] s t a t e _ s p a c e : l o w : [ - 6 . 0 , - 6 . 0 ] h i g h : [ 6 . 0 , 6 . 0 ] i n i t i a l _ r e g i o n : l o w : [ - 5 . 0 , - 5 . 0 ] h i g h : [ 5 . 0 , 5 . 0 ] b a r r i e r _ s e t t i n g s : b a r r i e r _ t y p e : " S O S " t i m e _ h o r i z o n : 1 0 b a r r i e r _ d e g r e e : 6 Figure 3: Input-output overview of StochasticBarrier.jl . s y s t e m _ f l a g : " n o n l i n e a r " d y n a m i c s : " b o u n d s / p e n d u l u m / d y n a m i c s _ 1 2 0 . n c " p r o b a b i l i t i e s : " d a t a _ 4 8 0 _ s i g m a _ [ 0 . 0 1 , 0 . 0 1 ] . n c " d i m : 2 𝜎 : : [ 0 . 0 1 , 0 . 0 1 ] i n i t i a l _ r e g i o n : c : [ 0 . 0 , 0 . 0 ] r : [ 0 . 0 1 ] b a r r i e r _ s e t t i n g s : b a r r i e r _ t y p e : " P W C " o p t i m i z a t i o n _ t y p e : " D U A L _ A L G " t i m e _ h o r i z o n : 1 0 The input parameters can be change d according to the scope of the problem. For P WC-SBF , the options for the barrier optimiza- tion_type are: "DU AL_ALG" , "CEGIS_ALG" , or "GD_ALG" . This interface is primarily intended to serve as a means for utilizing the toolbox’s capabilities. However , the modular design ensures that StochasticBarrier.jl can be seamlessly integrated into larger frameworks or other toolboxes. Such integration enables researchers and dev elopers to leverage its features within broader systems. B Benchmark Systems B.1 Linear Systems Contraction Map. W e consider two settings of contractive linear stochastic system in R 2 . For the rst setting, use d for the SOS benchmarks, the dynamics (Contraction Map 1) are x ( 𝑘 + 1 ) = 0 . 95 𝐼 x ( 𝑘 ) + w , where 𝐼 is the identity matrix and w ∼ N ( 0 , 10 − 2 𝐼 ) . The domain and initial sets are 𝑋 = [ − 1 , 2 ] 2 and 𝑋 0 = [ − 0 . 05 , 0 . 05 ] 2 , respectively . For the second set of benchmarks (P WC vs SoS engines), the Contraction Map 2 dynamics are x ( 𝑘 + 1 ) = 0 . 50 𝐼 x ( 𝑘 ) + w , using the same distribution w ∼ N ( 0 , 10 − 2 𝐼 ) . The domain and initial sets are 𝑋 = [ − 1 , 0 . 5 ] 2 and 𝑋 0 = [ − 0 . 05 , 0 . 05 ] 2 , respectively . T wo T ank. Next, we consider a stochastic discrete-time two tank model from [17], with dynamics x 1 ( 𝑘 + 1 ) =  1 − 𝜏 𝛼 1 𝐴 1  x 1 ( 𝑘 ) + 𝜏 𝑞 1 𝐴 1 + w 1 , x 2 ( 𝑘 + 1 ) = 𝜏 𝛼 1 𝐴 2 x 1 ( 𝑘 ) +  1 − 𝜏 𝛼 2 𝐴 2  x 2 + 𝜏 𝑞 0 𝐴 2 + w 2 ( 𝑘 ) , where x is the uid heights in the two tanks and w ∼ N ( 0 , 10 − 4 𝐼 ) the additive noise. The system parameters are 𝜏 = 0 . 1 , 𝛼 1 𝐴 1 = 1 , 𝑞 1 𝐴 1 = 4 . 5 , 𝛼 1 𝐴 2 = 1 , 𝛼 2 𝐴 2 = 1 , and 𝑞 𝑜 𝐴 2 = − 3 . The domain and initial sets are 𝑋 = [ 1 , 9 ] 2 and 𝑋 0 = [ 2 . 75 , 3 . 25 ] 2 , respectively . Room T emperature model. The 3D room temperature model de- scribed in [20] has deterministic dynamics x 1 ( 𝑘 + 1 ) = ( 1 − 𝜏 ( 𝛼 + 𝛼 𝑒 ) ) x 1 ( 𝑘 ) + 𝜏 𝛼 x 2 ( 𝑘 ) + 𝜏 𝛼 𝑒 𝑇 𝑒 x 2 ( 𝑘 + 1 ) = ( 1 − 𝜏 ( 2 𝛼 + 𝛼 𝑒 ) ) x 2 ( 𝑘 ) + 𝜏 𝛼 ( x 1 ( 𝑘 ) + x 3 ( 𝑘 ) ) + 𝜏 𝛼 𝑒 𝑇 𝑒 x 3 ( 𝑘 + 1 ) = ( 1 − 𝜏 ( 𝛼 + 𝛼 𝑒 ) ) x 3 ( 𝑘 ) + 𝜏 𝛼 x 2 ( 𝑘 ) + 𝜏 𝛼 𝑒 𝑇 𝑒 to which noise w ∼ N ( 0 , 10 − 2 𝐼 ) is adde d. The states x 1 , x 2 , and x 3 dene the temperature of each room. The ambient temperature Mazouz et al. 𝑇 𝑒 = 10 ◦ C, while the heat coecients are 𝛼 𝑒 = 8 · 10 − 3 and 𝛼 = 6 . 2 · 10 − 3 . The sampling time 𝜏 = 5 minutes. Regions of interest are 𝑋 = [ 17 , 29 ] 3 and 𝑋 0 = [ 18 , 19 ] 3 . Quadrotor . W e consider the lateral quadrotor model from [10] Δ ¤ 𝑦 𝐸 = Δ 𝑣 , Δ ¤ 𝑣 = 𝑔 Δ 𝜙 , Δ ¤ 𝜙 = Δ 𝑝, Δ ¤ 𝑝 = 𝐼 − 1 𝑥 Δ 𝐿 𝑐 , where 𝑦 𝐸 ∈ [ − 0 . 5 , 2 . 0 ] is the 𝑦 − position, and 𝑣 ∈ [− 1 . 0 , 1 . 0 ] is the corresponding velocity , 𝑔 is the gravitational acceleration, 𝐼 𝑥 is the inertia about the 𝑥 − axis, 𝜙 ∈ [− 0 . 1 , 0 . 1 ] is the roll angle, 𝑝 ∈ [ − 0 . 1 , 0 . 1 ] is the roll rate, and 𝐿 𝑐 is the roll control moment. The lateral model is appended with guidance for vertical motion Δ ¤ 𝑧 𝐸 = Δ 𝑤 , Δ ¤ 𝑤 = 1 𝑚 Δ 𝑍 𝑐 , where 𝑚 denotes the quadrotor mass, 𝑧 𝐸 ∈ [ − 0 . 5 , 3 . 0 ] denotes the 𝑧 -p osition and 𝑤 ∈ [ − 0 . 5 , 1 . 5 ] is the yaw rate. The initial set is ℓ = [ − 0 . 01 , 0 . 01 , − 0 . 01 , 0 . 01 , − 0 . 01 , 0 . 01 ] , 𝑢 = [ 0 . 01 , 0 . 03 , 0 . 01 , 0 . 03 , 0 . 01 , 0 . 03 ] , 𝑋 0 = { 𝑥 ∈ R 6 | ℓ ≤ 𝑥 ≤ 𝑢 } . The control law follows dir e ctly from a stabilizing LQR contr oller at the nominal point. Using Euler’s method with Δ 𝑡 = 0 . 01 , a discrete- time model is obtained, to which w ∼ N ( 0 , 10 − 6 𝐼 ) is added. B.2 Polynomial Systems Thermostat Model. A polynomial room temp erature mo dule taken from [13]. The stochastic dierence equation is given by x 𝑘 + 1 = ( 1 − 𝛽 − 𝜃 𝜈 ) x 𝑘 + 𝜃 𝑇 ℎ 𝜈 + 𝛽𝑇 𝑒 + w , where x 𝑘 is the room temperature, the ambient temperature 𝑇 𝑒 = − 15 ◦ C, and the heater temperature 𝑇 ℎ = 45 ◦ C. Conduction factors are given as 𝜃 = 0 . 145 , 𝛽 = 0 . 6 , and 𝜈 = − 0 . 0120155 x 𝑘 + 0 . 8 . The system evolves o ver the state space 𝑋 = [ 1 , 50 ] with the initial set 𝑋 0 = [ 19 . 5 , 20 ] , and unsafe regions are 𝑋 𝑢 = [ 1 , 17 ] ∪ [ 23 , 50 ] . The noise parameter is w ∼ N ( 0 , · 10 − 2 𝐼 ) . V an der Pol Oscillator . The V an Der Pol oscillator is from [ 1 ] has the following dynamics x 1 ( 𝑘 + 1 ) = x 1 ( 𝑘 ) + 𝜏 x 2 ( 𝑘 ) + w 1 , x 2 ( 𝑘 + 1 ) = x 2 ( 𝑘 ) + 𝜏 ( − x 1 ( 𝑘 ) + ( 1 − x 1 ( 𝑘 ) 2 ) x 2 ( 𝑘 ) ) + w 2 . In this model, 𝜏 = 0 . 1 is the sampling time. The domain of operation is dened by 𝑋 = [ − 6 , 6 ] 2 and the initial set by 𝑋 0 = [ − 5 , 5 ] 2 . The noise parameter is w ∼ N ( 0 , 4 · 10 − 4 𝐼 ) . B.3 P W A Inclusion Systems Pendulum NNDM. In this case study , we consider the p endulum NNDM model from [ 11 ]. The dynamics ar e x = 𝑓 𝑁 𝑁 ( 𝑥 ) + w , wher e 𝑓 𝑁 𝑁 is trained on a Pendulum agent in a closed-loop with a NN controller . For this system, noise term w ∼ N ( 0 , 10 − 4 𝐼 ) . The control limit are 𝑢 ∈ [ − 1 , 1 ] and the dynamics of the p endulum are ¤ 𝜃 𝑘 + 1 = ¤ 𝜃 𝑘 + 3 𝑔 2 𝑙 sin ( 𝜃 𝑘 ) 𝛿 𝑡 2 + 3 𝑚𝑙 2 𝑢𝛿 𝑡 2 , 𝜃 𝑘 + 1 = 𝜃 𝑘 + ¤ 𝜃 𝑘 + 1 𝛿 𝑡 . The initial set for this sytem is 𝑋 0 = [ 0 . 0 , 0 . 01 ] 2 . Hybrid Unicycle. Finally , we considered the wheeled mobile robot example from [10] with the dynamics of a unicycle: ¤ 𝑥 = 𝑣 cos 𝜃 , ¤ 𝑦 = 𝑣 sin 𝜃 , ¤ 𝜃 = 𝜔 , ¤ 𝑣 = 𝑎, where 𝑥 ∈ [ − 1 . 0 , 0 . 5 ] and 𝑦 ∈ [ − 0 . 5 , 1 . 0 ] are the Cartesian positions. The angle 𝜃 ∈ [ − 1 . 75 , 0 . 5 ] is the orientation with respe ct to the 𝑥 − axis, and 𝑣 ∈ [ − 0 . 5 , 1 . 0 ] is the speed of the vehicle. The initial set is dened by a hyperrectangle, where 𝑥 ∈ [ − 0 . 01 , 0 . 01 ] 4 . The inputs are steering rate 𝜔 and acceleration 𝑎 . The combined feedback linearization and LQR stabilizer make this system hybrid. The Euler method with Δ 𝑡 = 0 . 01 denes the discrete-time dynamics, to which noise term w ∼ N ( 0 , 10 − 4 𝐼 ) is added. C Saving and loading of probability data Our toolbox supports storing transition probability bounds and PW A bounds for non-linear functions as NetCDF, which allows for ecient storage and r etrieving of multidimensional data in addition to axis annotation and attributes [ 18 ]. First, to save the transition data, one can run the following to rst create a NetCDF dataset and saving it. # S a v e t o a . n c f i l e f i l e n a m e = " p r o b a b i l i t y _ d a t a . n c " d a t a s e t = c r e a t e _ p r o b a b i l i t y _ d a t a s e t ( p r o b a b i l i t i e s ) # N e t C D F d a t a s e t s a v e d a t a s e t ( d a t a s e t , p a t h = f i l e n a m e ; d r i v e r = : n e t c d f , o v e r w r i t e = t r u e , c o m p r e s s = 1 ) The process in reverse is to rst load the NetCDF dataset followed by extracting and transforming data to the appropriate format via the load_probabilities function. f i l e n a m e = " s e t / p a t h / t o / f i l e . n c " d a t a s e t = o p e n _ d a t a s e t ( f i l e n a m e ) p r o b a b i l i t i e s = l o a d _ p r o b a b i l i t i e s ( d a t a s e t ) T o load PW A bounds, use the load_dynamics(dataset) function, which expects in the .nc le a num_regions property as well as 3 datasets: nominal_dynamics_A , nominal_dynamics_b , and regions . The dataset regions should have three axes: region , dir , and x of dimension 1 to num_regions , [ "lower" , "upper" ], and 1 to dim , respectively . nominal_dynamics_A and nominal_dynamics_b are similar with axes region , dir , y , x and region , dir , y , respec- tively . Ordering of axes is irrelevant; the tool will select the axis based on their name. Furthermore , the tool executes a number of consistency checks and will throw an informative err or ab out the inconsistency if any . D T able Results Detailed results for the SOS and P WC benchmarks are presented in T able 3 and 4, respectively . StochasticBarrier .jl: A T oolb ox for Stochastic Barrier Function Synthesis T able 3: Benchmark SOS methods: StochasticBarrier.jl vs MA TLAB and Python implementations for linear and polynomial systems. All verications are performed for N = 10 time-steps. OM denotes out-of-memor y , while FAILED indicates that the algorithm terminated without a valid barrier . The shaded areas indicate that the toolb ox cannot accommodate these systems. Class System Deg. StochasticBarrierFunc. [19] PRoTECT [20] StochasticBarrier.jl 𝜏 (s) 𝜂 𝛽 𝑃 𝑠 𝜏 (s) 𝜂 𝛽 𝑃 𝑠 𝜏 (s) 𝜂 𝛽 𝑃 𝑠 Linear 2D Contraction Map 1 2 23.56 1.00 0.00 0.00 0.58 1 . 7 𝑒 − 2 1 . 8 𝑒 − 2 0.80 0.01 1 . 7 𝑒 − 2 1 . 8 𝑒 − 2 0.80 4 74.86 1.00 0.00 0.00 1.51 3 . 7 𝑒 − 3 5 . 7 𝑒 − 3 0.94 0.02 3 . 7 𝑒 − 3 5 . 7 𝑒 − 3 0.94 8 543.29 1.00 0.00 0.00 12.85 1 . 1 𝑒 − 3 1 . 8 𝑒 − 3 0.98 0.05 1 . 1 𝑒 − 3 1 . 8 𝑒 − 3 0.98 12 2582.17 1.00 0.00 0.00 76.88 1 . 0 𝑒 − 3 1 . 3 𝑒 − 3 0.99 0.21 1 . 1 𝑒 − 3 1 . 3 𝑒 − 3 0.99 24 OM × × × 2442.43 8 . 3 𝑒 − 6 9 . 5 𝑒 − 8 0.99 18.08 9 . 9 𝑒 − 4 1 . 2 𝑒 − 3 0.99 30 OM × × × 9569.64 4 . 6 𝑒 − 7 7 . 8 𝑒 − 9 0.99 107.14 9 . 6 𝑒 − 4 1 . 2 𝑒 − 3 0.99 2D T wo Tank 4 90.87 1.00 0.00 0.00 2.37 1 . 00 3 . 6 𝑒 − 6 0.00 0.14 8 . 7 𝑒 − 2 4 . 1 𝑒 − 2 0.50 6 283.36 1.00 0.00 0.00 9.48 FAILED FAILED 0.00 0.10 2 . 9 𝑒 − 2 7 . 9 𝑒 − 5 0.97 8 884.79 1.00 0.00 0.00 53.60 FAILED FAILED 0.00 0.16 6 . 8 𝑒 − 4 4 . 2 𝑒 − 5 1.00 3D Room T emp. 2 109.74 FAILED FAILED 0.00 1.38 FAILED FAILED 0.00 0.12 FAILED FAILED 0.00 4 OM × × × 39.02 FAILED FAILED 0.00 0.27 9 . 7 𝑒 − 2 7 . 6 𝑒 − 3 0.83 6 OM × × × 2422.07 FAILED FAILED 0.00 11.08 9 . 0 𝑒 − 7 8 . 6 𝑒 − 8 1.00 8 OM × × × OM × × × 20.02 8 . 7 𝑒 − 6 2 . 1 𝑒 − 6 1.00 12 OM × × × OM × × × 3716.59 4 . 1 𝑒 − 8 9 . 7 𝑒 − 9 1.00 6D Quadrotor 2 OM × × × 9587.45 FAILED FAILED 0.00 0.14 1 . 0 𝑒 − 1 5 . 9 𝑒 − 4 0.89 4 OM × × × OM × × × 3.35 9 . 0 𝑒 − 4 1 . 0 𝑒 − 5 1.00 6 OM × × × OM × × × 148.65 8 . 8 𝑒 − 6 1 . 7 𝑒 − 7 1.00 Polynomial 1D Thermostat 4 0.59 FAILED FAILED 0.00 0.05 0.87 3 . 3 𝑒 − 3 0.09 6 0.79 FAILED FAILED 0.00 0.08 3 . 5 𝑒 − 11 3 . 3 𝑒 − 12 1.00 8 1.13 FAILED FAILED 0.00 0.16 4 . 2 𝑒 − 5 8 . 5 𝑒 − 6 1.00 2D Oscillator 6 21.21 FAILED FAILED 0.00 0.96 1.00 1 . 8 𝑒 − 8 0.00 8 168.79 FAILED FAILED 0.00 3.68 1.00 2 . 9 𝑒 − 8 0.00 12 5502.47 FAILED FAILED 0.00 77.03 1 . 2 𝑒 − 7 1 . 2 𝑒 − 8 1.00 T able 4: Benchmark: SOS vs PWC implementations for linear and PW A inclusion systems using StochasticBarrier.jl . All verications are performed for N = 10 time-steps. OM denotes out-of-memory . Class System SOS PWC DUAL LP PWC CEGIS PWC GD | Q | Deg. 𝜏 (s) 𝑃 𝑠 | Q | 𝜏 (s) 𝑃 𝑠 𝜏 (s) 𝑃 𝑠 𝜏 (s) 𝑃 𝑠 Linear 2D Contraction Map 2 NA 4 4.36 0.23 64 0.44 0.99 2.43 0.99 1.57 0.95 NA 8 13.33 0.43 225 105.26 0.99 19.73 0.99 10.08 0.97 NA 20 314.22 0.88 289 248.53 0.99 39.40 0.99 34.53 0.97 NA 30 2162.92 0.90 361 2056.58 0.99 50.14 0.99 49.18 0.98 PW A Inclusion 2D Pendulum 120 2 29.10 0.93 120 0.41 0.88 1.94 0.88 2.02 0.88 240 2 43.42 0.93 240 101.38 0.94 13.70 0.94 4.81 0.94 480 2 118.45 0.94 480 208.65 0.97 44.70 0.97 18.79 0.96 4D Unicycle 1250 4 1819.45 0.00 1250 1608.43 0.65 122.30 0.64 14.76 0.64 1800 4 2455.66 0.00 1800 2143.80 0.65 211.51 0.65 20.07 0.64