Computer-implemented method for generating a dynamics model

Abstract

Provided is a computer-implemented method for generating a dynamics model of a system. The method comprises generating the dynamics model, optionally including an IMC, by determining unknown dynamics of the system using a predetermined kernel-based method designed for regression and / or classification, optionally BTGP regression. Generating the dynamics model comprises computing probability bounds for each transition or multiple transitions of the dynamics model by obtaining transition probability distributions from the kernel-based method that are amenable to a predefined model checking algorithm, and determining an abstraction error of the model based on the computed probability bounds.

Description

[0001]Tautz & Schuhmacher Law MPI1117P12WO June 3, 20251 / 39 COMPUTER-IMPLEMENTED METHOD FOR GENERATING A DYNAMICS MODEL Provided is a computer-implemented method for generating a dynamics model of a system, involving generating the model using kernel-based regression and classification methods using a kernel that leads to models that are easily amenable to model checking, a data processing device, a computer program, and a computer- readable medium related to this method. Developing the control software for safety-critical cyber-physical systems such as power grids, autonomous vehicles, medical implants, and robotic systems has raised the need for formal approaches to guarantee the systems will behave as expected, and not lead to catastrophe. Designing controllers for such systems is especially challenging for their stochasticity and heterogeneous nature, i.e., the components of the system evolve over both discrete and continuous spaces. Previous approacheswere mainly based on physics-based models of the underlying systems. To meet theincreasing need for more scalable and uncertainty-robust techniques, data-driven approaches govern the contemporary landscape of active research. One powerful offspring that is being explored is Gaussian process (GP) regression (Rasmussen and Williams, 2018). GPs use a Bayesian paradigm to learning unknown functions from observations, that is particularly valuable for safety-critical applications as the statistical approximation error is naturally quantified. In the following a brief review of related approaches on data-driven formal verification and synthesis of uncertain stochastic systems is given. The works by Badings et al.(2023) and Lavaei et al. (2022) address epistemic uncertainty in stochastic systems byabstracting the system to an interval Markov decision process (IMDP). While prior is limited to finite horizons and linear systems, latter requires solving a scenario optimization problem. Jackson et al. (2020) and Jiang et al. (2022) have used GPs instead of scenario optimization. In a related kernel-based setting, Thorpe et al. (2022) use conditional mean embeddings (CMEs) to embed the conditional distribution of a non-Markovian random trajectory into a reproducing kernel Hilbert space (RKHS). Temporal logic control of systems captured by stochastic neural network dynamic models (NNDM) is addressed by Adams et al. (2022) through translating the NNDMsTautz & Schuhmacher Law MPI1117P12WO June 3, 20252 / 39 to IMDPs. The starting point, however, is the NNDM and no guarantees on thecorrectness of the learned model with respect to the data-generating system areconsidered, rendering the approach only robust from a model-based perspective. Capturing epistemic uncertainty using interval MDPs relies on the assumption of state- wise independent epistemic uncertainty, leading to inconsistent and overly conservative results. In contrast, Schön et al. (2023) capture epistemic uncertainty explicitly using parametric MDPs and probabilistic coupling relations to the latent true system. Furthermore, their results are applicable to a wider range of specifications and distributions. The neural abstractions studied by Abate et al. (2022) rely on a counterexample-producing scheme via SMT solvers. A contraction-guided adaptive partitioning approach is introduced by Harapanahalli et al. (2023a).In figure 1 a method 11 for formal verification according to the state-of-the-art is shown.The method 11 comprises determining a model 13 of an actual system 12, the model13 being a system description of the actual system 12, wherein the model 13 isobtained by abstraction of data describing the actual system 12 (or its behavior). Themodel 13 is input into an algorithm 14 together with a temporal logics specification 16of a property 15 that the actual system 12 should satisfy, e.g., a target state in whichthe actual system should be kept in. The algorithm 14 is configured to check or(formally) very if the actual system 12 is able to satisfy the temporal logics specification16 based on the model 13 and the temporal logics specification 16. This allows a checkfor probabilistic satisfaction, a discretization error compensation, and a data- driven / statistical error compensation. In the specific nonlimiting example shown infigure 1, the actual system 12 is a driver assistance function for a car, namely a lanekeeping assistant configured to automatically take steps to ensure the car or vehiclestays in its lane. In other words, the lane keeping assistant acts as a controller for thecar. The property 15 the driver assistance function has to satisfy is to keep the car inits lane (i.e., “Stay in lane”). The model 13 represents the driver assistance function orthe behavior of the driver assistance function and is determined by abstraction of datarelating to said behavior. The algorithm 14 is configured to check or (formally) verify ifthe driver assistance function is configured to keep the car in the lane as required bythe temporal logics specification 16 based on the model 13 of the driver assistanceTautz & Schuhmacher Law MPI1117P12WO June 3, 20253 / 39 function. One drawback is that the computational power needed to generate and verifythe model 13 is very high.A modified state-of-the-art approach is shown in figure 2. This approach is – similar tothe one described with respect to figure 1 – an abstraction-based approach. Themethod 21 for formal verification shown in figure 2 comprises – in difference to themethod 11 shown in figure 1 – determining a model 23 or system description of theactual system 22 in a two-step process. The model 23 comprises two models, anintermediate data or first model 231 generated in a first step and a second model 232 generated based on the first model 231 in a second step. The first model 231 is generated from the data describing the system 22. The second model 232 is generatedbased on an abstraction of the intermediate data or first model 231. The intermediatedata model 231 is generated based on the data describing the (behavior of the) actualsystem 22 using Gaussian Processes (GP). The second model 232 is represented byan Interval Markov Chain (IMC) or an Interval Markov decision process (IMDP). Themodel 23 is input into an algorithm 24 together with a (temporal logics, optionallyreachability) specification(s) 26 of the property 25 (e.g., “Stay in lane”) that the actualsystem 22 should satisfy. The algorithm 24 is configured to check or (formally) very ifthe actual system 22 is able to satisfy the specification 26 based on the model 23 andthe specification 26 using value iteration. The output of the value iteration are boundson satisfaction probability. However, generating the second model 232 based on theintermediate data or first model 231 still needs a lot of computational power.The steps for generating the second (IMC / IMDP) model 232 via the intermediate data(GP) model 231 and for the value iteration using algorithm 24 are shown in more detailin figure 3. The GP model 231 learns a mean function and a covariance function fromthe data describing the actual system 22. This provides a probabilistic estimate oftransition probabilities rather than fixed values. Confidence bounds can be determinedbased on the learned covariance function and mean function. The confidence boundsrepresent an uncertainty in the predictions. The IMC / IMCP model 232 is thengenerated based on the trained GP model 231 by abstraction of the GP model 231.More specifically, the second model 232 defines bounds ^, ^ on the transitionprobability, i.e., an upper bound ^ and a lower bound ^ for the probability that the actualTautz & Schuhmacher Law MPI1117P12WO June 3, 20254 / 39system 22 transitions from one (discrete) state ^^ into another (discrete) state ^^^^(wherein t is the time), as follows: ^(^^^^|^^) < ^(^^^^|^^) < ^(^^^^|^^).The bounds ^, ^ are then used together with the specification 26 to perform valueiteration in order to determine bounds on the satisfaction probability. However, thebounds ^, ^ are computed via integration for every state pair ^^, ^^^^ by optimizing overcontinuous cells. Therefore, a lot of computing power is needed to compute the bounds^, ^ on the transition probability.Against the background of this state of the art, one object of the present disclosure is to specify a device and / or a method, each of which is suitable for enriching the state of the art. One optional concrete object of the disclosure may be formulated as to reduce a computing power needed to generate the model from the intermediate model, i.e., to reduce the computing power needed to compute the upper and lower bounds of the transition probability. The object is solved by the features of the independent claims. The dependent claims each relate to optional further embodiments of the disclosure. Provided is a computer-implemented method for generating a dynamics model of a system, the method comprising generating the dynamics model, optionally includingan interval Markov chain (IMC) or interval Markov decision process (IMDP), bydetermining unknown dynamics of the system using a (optionally predetermined) kernel-based method designed for regression and / or classification, optionally binary- tree Gaussian process (BTGP) regression. The method further includes computing probability bounds for each transition (or optionally multiple transitions) of the generated dynamics model by obtaining transition probability distributions from the kernel-based method that are amenable to a predefined model checking algorithm, and determining an abstraction error of the model based on the computed probability bounds.Tautz & Schuhmacher Law MPI1117P12WO June 3, 20255 / 39 That is, kernel-based methods for regression or classification equipped with a kernel yielding an RKHS of functions efficiently amenable to formal verification / model checking can be provided. The system can be a continuous-time or discrete-time, heterogeneous-space including discrete and / or continuous variables, stochastic or non-stochastic system. The obtained transition probability distributions can be piecewise-constant or of other piecewise-predefined, periodic, or otherwise hierarchically or structurally correlated shape in order to be amenable to the predefined model checking algorithm. The abstraction error can characterize a deviation of a behavior of the generated model from the system. It is possible that only finitely many observations of the system are available for generating the dynamics model. It is possible that not more than two integrations are performed for each of the transitions to identify the worst-case probabilities for the respective transition. A kernel of the kernel-based method can be discontinuous and defined on a finite space of binary or non-binary strings. The generating of the model can include directly learning the model from the observations of the system without performing explicit steps of discretization. The system can represent a control system, optionally designed to be used in safety- critical environments. A data processing device comprising means for carrying out the above-described method can be provided.Tautz & Schuhmacher Law MPI1117P12WO June 3, 20256 / 39 A computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the above-described method can be provided. A computer-readable medium comprising instructions which, when executed by a computer, cause the computer to carry out the above-described method can be provided.A dynamics model can be defined as a computer-generated representation thatcaptures the evolution (or “dynamics”) of a system over time. It may be derived fromdata (optionally via a (machine) learning method such as Gaussian processregression) and quantify both the expected behavior and (potentially) the uncertainty(through computed probability bounds) of the system’s transitions. The term “system” refers to the entity whose behavior is to be modeled and verified. It can be any process or control mechanism whose state evolves over time. In this context, the system may be continuous-time or discrete-time, with a heterogeneous state space (i.e., containing both discrete and continuous variables) and may exhibit stochastic (random) or non-stochastic behavior.A Kernel-based method designed for regression and / or classification may be definedas a data-driven approach that employs reproducing kernels to estimate unknownfunctions based on observed data. As such, kernel-based methods can be used forboth regression (predicting continuous outputs) and classification tasks. For example, a function estimator, such as a Gaussian process, can be constructed by utilizing a kernel that maps input data into a reproducing kernel Hilbert space (RKHS) to capture complex, possibly non-linear relationships.A binary-tree Gaussian process (BTGP) can be defined as a specialized Gaussianprocess that uses a binary-tree (BT) kernel. The BT kernel can be defined on a finitespace of binary strings, which induces an RKHS of piecewise-constant functions. As a result, the posterior mean and covariance functions of the BTGP are inherentlyTautz & Schuhmacher Law MPI1117P12WO June 3, 20257 / 39 piecewise-constant, naturally yielding a finite (discrete) abstraction of the underlying continuous system that is easier to verify.Probability bounds can be defined as the computed upper and lower limits on thelikelihood of the system transitioning from one state to another. These bounds quantify the uncertainty in the predictions made by the model (for instance, the GP or BTGP) and are used in constructing the abstract model (such as an interval Markov chain) for formal verification. In the context of the abstracted dynamics model, a transition denotes the change of the system’s state from one discrete (abstracted) state to another. Each transition is associated with a probability that reflects the likelihood of moving between those states according to the learned model.Transition probability distributions (that are amenable to a predefined model checkingalgorithm) can be defined as the probability distributions over possible next states(transitions) as derived e.g. via a kernel-based method. In an easily amenable form,they are structured, optionally with a piecewise-predefined (e.g., piecewise-constant)form, so that they can be efficiently integrated and processed by a model checking algorithm during formal verification.A model checking algorithm can be defined as a formal verification tool thatsystematically explores the abstract model (for example, an interval Markov chain)against a given (temporal logic) specification. It computes whether (and to whatprobability) the modeled system satisfies certain safety or performance properties by processing the transition probability bounds.The Abstraction error may quantify the deviation between the behavior of theabstracted dynamics model (e.g., generated via a data-driven method) and the true, underlying system. It accounts for errors introduced by approximations, finite observations, and the inherent limitations of the chosen kernel and modeling approach. A continuous-time or discrete-time, heterogeneous-space including discrete and / or continuous variables, stochastic or non-stochastic system may define the broad classTautz & Schuhmacher Law MPI1117P12WO June 3, 20258 / 39 of systems that the method is applicable to. It covers systems whose state changes either continuously over time or in discrete time steps; whose state space may include both discrete and continuous elements; and which can be inherently stochastic (with randomness) or non-stochastic (deterministic).Worst-case probabilities for the respective transition can be defined as the extreme(upper and lower) probability values computed for a particular state transition. Theycan be determined by optimizing the computed probability bounds, optionally viaintegration, to identify the worst-case scenarios.A discontinuous kernel can be defined as a reproducing kernel that is not continuousand (typically) leads to discontinuous function estimates. Here, the binary-tree kernelis an example of a discontinuous kernel since it is defined on a finite space of binary(or non-binary) strings and induces e.g. piecewise-constant (and hence discontinuous) functions in the resulting Gaussian process model.Explicit steps of discretization may refer to the conventional procedures for partitioninga continuous state space into discrete regions (e.g., grid-based methods) prior to modeling. In contrast, the disclosed method learns the dynamics model directly from continuous observations without performing such separate, explicit discretization stepssince the discretization is inherent in the structure of the e.g. BT kernel.Safety-critical environments can be defined as application areas where system failurescan have severe consequences (such as significant property damage, or criticalservice disruptions etc.). Examples include autonomous vehicles, medical devices,power grids, and robotic systems. In these settings, ensuring the correctness and robustness of the system through formal verification is essential. BRIEF DESCRIPTION OF THE DRAWINGS In the drawings.FIG. 1 schematically shows a flowchart of a method for formal verificationaccording to the state of the art;Tautz & Schuhmacher Law MPI1117P12WO June 3, 20259 / 39FIG. 2 schematically shows a flowchart of another method for formal verificationaccording to the state of the art; FIG.3 shows some steps of the method of figure 2 in more detail;FIG. 4 schematically shows a flowchart of a method for formal verificationaccording to the disclosure;FIG. 5 shows an output of a GP used in the method of figure 4;Fig. 6 visualizing a technical effect induced by the disclosure;Fig. 7 depicts a piecewise-constant posterior mean ^̂^ and double standarddeviation ^̂^ ± 2^^^ as functions of state of a GP utilizing a BT kernel ofprecision ^ = 2;Fig. 8 shows a lower bound on the infinite-horizon reachability probability as afunction of the initial state. DESCRIPTION In the following, details are set forth to provide a more thorough explanation of the disclosure. However, it will be apparent to those skilled in the art that the disclosure can be practiced without these specific details. In other instances, well-known structures and devices are shown in block diagram form or in a schematic view rather than in detail in order to avoid obscuring the disclosure. In addition, features describedhereinafter can be combined with each other, unless specifically noted otherwise. It isalso to be understood that other features can be utilized and structural or logical changes can be made without departing from the scope defined by the claims. The following detailed description, therefore, is not to be taken in a limiting sense. Further, equivalent or like elements or elements with equivalent or like functionality are denoted in the following description with the same or like reference numerals. As theTautz & Schuhmacher Law MPI1117P12WO June 3, 202510 / 39 same or functionally equivalent elements are given the same or like reference numbers in the figures, a repeated description for elements provided with the same or like reference numbers can be omitted. Hence, descriptions provided for elements having the same or like reference numbers are mutually exchangeable as long as not indicated otherwise. Directional terminology, such as “top”, “bottom”, “above”, “below”, “front”, “back”, “behind”, “leading”, “trailing”, “over”, “under”, etc., can be used with reference to the orientation of the figures and / or elements being described. Because the features can be positioned in a number of different orientations, the directional terminology is used for purposes of illustration and is in no way limiting. In some instances, directional terminology can be exchanged with equivalent directional terminology based on the orientation of a feature so long as the general directional relationships between elements, and the general purpose thereof, is maintained. In the present disclosure, expressions including ordinal numbers, such as “first”, “second”, and / or the like, can modify various elements. However, such elements are not limited by the above expressions. For example, the above expressions do not limitthe sequence and / or importance of the elements. The above expressions are usedmerely for the purpose of distinguishing an element from the other elements. For example, a first feature and a second feature indicate different features, although both are features. For further example, a first element could be termed a second element, and similarly, a second element could also be termed a first element without departing from the scope of the present disclosure. Neither the drawings nor the following descriptions thereof are in anyway intended to limit the scope of the disclosure and are provided for exemplary purposes only. In figure 4 a flowchart of a computer-implemented method 41 for generating adynamics model 43 of a system 42 is shown.The method 41 comprises a step of generating the dynamics model 43, wherein thestep of generating the dynamics model 42 in turn comprises two sub steps 431, 432.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202511 / 39 The generated dynamics model 43 comprises a GP model. The GP model is learnedin the first sub step 431 directly, i.e., without discretization, from data describing thesystem 42 and / or its behavior. That is, the method comprises directly learning themodel 43 from the observations of the system 42 without performing explicit steps ofdiscretization. The system 42 is a continuous-time or discrete-time, heterogeneous- space including discrete and / or continuous variables, stochastic or non-stochasticsystem 42. The system 42 may represent a control system, optionally designed to beused in safety critical environments. Furthermore, only finitely many observations ofthe system 42, i.e., a limited amount of data or data points, are available for generatingthe dynamics model 43. The GP model is learned or generated by determiningunknown dynamics of the system 42 using a kernel-based method designed for regression and / or classification. The kernel used in the GP of the first sub step 431dictates the reproducing kernel Hilbert space (RKHS) of GP mean functions. Ingeneral, a kernel may be used that induces an RKHS of a piecewise constant or linearshape. Therefore, a kernel of the kernel-based method 431 can be discontinuous anddefined on a finite space of binary or non-binary strings. Here BTGP regression is used,i.e., a GP with a kernel inducing an RKHS of piecewise-constant functions is used. TheBT kernel is defined as follows: where ^ ∈ ℕ is the precision of the kernel, ^^, … , ^^ ∈ ℝ are weights, and Ψ^(^) maps^ to a binary string of length ^.Using the BT kernel leads to a piecewise constant shape of the mean and of thecovariance, as can be gathered from figure 5. That is, the mean function and the covariance function are both piecewise-constant.Fig. 6 visualizes the technical effect induced by the disclosure, i.e. not more than twointegrations are performed / needed for each of the transitions to identify the worst-case probabilities for the respective transition. That is, due to the piecewise constant shapeof the BTGP predictions (right in figure 6), only one integration per IMC / IMDP bound isnecessary. In contrast, for regular GPs (left in figure 6), the arbitrary shape brings theTautz & Schuhmacher Law MPI1117P12WO June 3, 202512 / 39 need to optimise for the worst-case transition probabilities over the continuous cellss_0, s_1, … Therefore, less computing power needed is needed when using themethod according to the disclosure. For the confidence bound, the piecewise-constant shape gives rise to additional error terms epsilon_2 and epsilon_3. This shows that it is a tradeoff between efficiency andconservatism, but also that the error is considered, and it is still formally sound.In the second sub step 432 probability bounds for each transition or multiple transitionsof the dynamics model 43 are computed by obtaining transition probability distributionsfrom the kernel-based method 431 that are amenable to a predefined model checkingalgorithm. As a result of the BT kernel, the obtained transition probability distributionsare of a piecewise-predefined shape, here piecewise-constant, in order to beamenable to the predefined model checking algorithm. In other words, because the obtained transition probability distributions are of a piecewise-predefined shape, herepiecewise-constant, they are amenable to the predefined model checking algorithm.The obtained transition probability distributions are used for an IMC model formingtogether with the GP model the dynamics model 43. In a last step of the method 41 an abstraction error of the generated dynamics model 43 is determined based on the computed probability bounds, wherein the abstraction error characterizes a deviation of a behavior of the generated dynamics model 43 fromthe system 42.In the following a specific example of implementing the disclosure is described in detail,the following description should not be construed as limiting the disclosure:Autonomous vehicles, medical devices, energy and communication grids, and robotic systems are just some products of modern engineering that embody machines operating in safety-critical environments. Designing suitable control software for such systems challenges engineers to create tools that can certify that the system will behave as expected once deployed based on its training performance. In doing so, such tools need to deal with the fact that neither the behaviour of the system can be tested on every possible operating scenario, nor do systems in practise behave in aTautz & Schuhmacher Law MPI1117P12WO June 3, 202513 / 39 deterministic way. Here, a novel approach to obtaining models of dynamical systems for which the abstraction error can be quantified and that can be directly amenable to symbolic formal verification and synthesis, is provided. Said abstract models can be directly learned from observations of the system without requiring explicit steps of discretization. This can be achieved by constraining the scope of candidate models to inherently discrete and hence naturally amenable representations. By establishing results on the relation between the generated abstraction and the actual data- generating system, probabilistic certifications of the well-behavedness of the actual system can be computed and transferred even when the dynamics of the true system are stochastic and unknown, and only finitely many observations are available. More specifically, that is, GPs based on the so-called binary-tree (BT) kernel (Cohen et al., 2022b) are used. The BT kernel is a discontinuous kernel that has been shown to reduce the computational burden of GP regression significantly, increasing its scalability to higher dimensional systems and larger data sets. In comparison toalternative approaches to constructing scalable GPs (see, e.g., Quinonero-Candelaand Rasmussen (2005)), the BT kernel achieves GP regression in linear time – apolynomial improvement. The resulting binary-tree Gaussian process (BTGP) is characteristic for its piecewise-constant posterior mean and covariance functions, resulting in a natural partitioning of its input space. The main idea of this disclosure can be to exploit this direct learning of an abstract representation of the unknown data- generating system for formal verification. Formal verification of known stochastic systems against temporal logic properties has been well studied (Baier and Katoen, 2008). For settings where the system dynamicsare uncertain – termed epistemic uncertainty — and specifications have infinitehorizons, most available abstraction-based approaches can be divided into two categories. Approaches based on statistical relations replace the original system with an estimated surrogate model that does not capture the epistemic uncertainty, generally a variant of a Markov chain (MC). Results obtained on this model hence rely on quantifying the probabilistic deviation to the true system to transfer guarantees (Abate et al., 2008; Haesaert et al., 2018). For example, probabilistic coupling relationscan be used to transfer guarantees obtained on a parameterized surrogate model ofTautz & Schuhmacher Law MPI1117P12WO June 3, 202514 / 39 the unknown system back to the original latent system by quantifying the expected parametric uncertainty (Schön et al., 2023). Verification is conducted using a robust version of dynamic programming (Haesaert and Soudjani, 2020). Establishing suchcoupling relations, however, requires the computation of local error parametersquantifying the error between the true system and its surrogate model for every discrete partition, deeming their application to high-dimensional systems computationally excruciating. In an effort to avoid the construction of surrogate models requiring such types of relations, approaches based on interval Markov chains (IMCs) incorporate both stochasticity and epistemic uncertainty (Badings et al., 2023; Lavaei et al., 2022; Jackson et al., 2020). The IMCs constructed based on upper and lower bounds on the transition probabilities, learned, e.g., using GP regression (Jackson et al., 2020), can be subsequently verified via the interval iteration algorithm (Haddad and Monmege, 2018). Whilst this alleviates the need for constructing additional coupling relations, the computations required to construct and verify the IMC are substantially more involving than the corresponding computations on the MC in the coupling-based approach. This compromise can be overcome using the solution according to the disclosure. BTGP regression can be used to learn the unknown dynamics of a discrete-time continuous-space stochastic system and show that the probability bounds of the resulting IMC are considerably easier to compute by obtaining transition probability distributions from the BTGP that are naturally piecewise-constant. Whilst prior IMC- based approaches require the optimization of both probability bounds to identify the worst-case probabilities for every transition, this change allows to perform only twointegrations for each transition – one for either bound – and obtain equivalentrobustness results. Hence, a speedup for the construction of the abstraction that is polynomial w.r.t. the size of the abstraction can be achieved. Thus, the disclosure provides the following advantages:- A data-driven approach is disclosed to verify uncertain stochastic systems definedover continuous spaces against finite-horizon and / or infinite-horizon specifications based on an IMC generated via BTGP regression. In comparison to existing approaches, the IMC can be constructed without the need for the identification of theTautz & Schuhmacher Law MPI1117P12WO June 3, 202515 / 39 worst-case bounds for every partition. This yields a computational speedup over IMC abstraction that is polynomial in the number of partitions.- Statistical guarantees on the accuracy of the BTGP model can be obtained, thestatistical guarantees combining available measure-theoretic results with an error inflicted by the approximation error of the BTGP even if the true dynamics function does not live in the BT kernels reproducing kernel Hilbert space (RKHS).- The discontinuous kernel can be used to generate a finite abstraction of the latentmodel and this abstraction can be efficiently verified.- Closed-form expressions for the function space generated by the BT kernel can begiven.In the following one optional aspect of the disclosure will be described in detail. Thefollowing description is not intended to limit the disclosure to any specific aspectsdescribed therein. In this respect reference is made to Schön et. al. (Oliver Schön,Shammakh Naseer, Ben Wooding, Sadegh Soudjani: Data-Driven Abstractions viaBinary-Tree Gaussian Processes for Formal Verification, IFAC-PapersOnLine, Volume58, Issue 11, 2024, Pages 115-122, ISSN 2405-8963,https: / / doi.org / 10.1016 / j.ifacol.2024.07.434) which is incorporated herein in its entirety.The remainder of this detailed description is structured as follows. After providingpreliminary definitions and the problem statement system identification using BTGP regression will be described. Here, results examining the BTGP and its associatedstatistical error bound will be provided. Then it is described how the BTGP allows togenerate an IMC model efficiently and provide a suitable verification algorithm. Thenthe presented approach is showcased on a nonlinear case study. Due to spacerestrictions, the proofs of statements are relegated to the appendix. Preliminaries and Problem Statement:Let ℝ, ℝ^^, ℕ, and ^ be the sets of reals, non-negative reals, positive integers, andbinary numbers. Let s ∈ ^^ be a binary string of length |s| : = q. We denote its substringprefix of length l ∈ ℕ by [s]^^. Let I^ be the N × N identity matrix. The transpose of avector or matrix A is indicated by A^. For some topological space ^, let X^Tautz & Schuhmacher Law MPI1117P12WO June 3, 202516 / 39: = [x^, … , x^]^ be a column vector with x^ ∈ ^. To save space, we may use a compactnotation and write the column vector as X^ : = [x^]^^^^. We denote the element-wiseevaluation of a function f: ^ → ^ on X^ by f(X^) : = [f(x^)] Similarly, we may writeA = [a^^]^^,^^^to denote a matrix with its elements. The trace of such a matrix is denotedby tr(A) ≔ ∑^ a^^ . For b ∈ ⊤, ⊥, the indicator function is given by 1(b) : = 1, if b =⊤, and 0 otherwise.Probability theory: A probability space is a tuple (^, ℬ(^), p) equipped with a samplespace ^, a Borel σ-algebra ℬ(^) defined over ^, i.e., the smallest σ-algebra containingopen subsets of ^, and a probability measure p, which has realizations x ∼ p( ⋅ ). Inthis work, we restrict our attention to Polish sample spaces (Bogachev, 2007). Aprobability measure p on a measurable space ^^, ℬ(^)^ is a map p: ℬ(^) → [0,1] suchthat for all countable collections {A^^}^^^ of pairwise disjoint sets in ℬ(^) it holds thatp(⋃ A^) = ∑^ p (A^) and p(^) = 1. For two measurable spaces ^^, ℬ(^)^ and ^^, ℬ(^)^,a probability kernel is a mapping p: ^ × ℬ(^) → [0,1] such that p(x, ⋅ ): ℬ(^) → [0,1] isa probability measure for all x ∈ ^, and p( ⋅ , Y): ^ → [0,1] is measurable for all Y ∈ ℬ(^).A probability kernel associates to each point x ∈ ^ a measure, also denoted by p( ⋅ |x).Given n ∈ ℕ, the (Gaussian) normal probability measure with mean μ ∈ ℝ^ andcovariance matrix Σ ∈ ℝ^×^ is denoted as ^( ⋅ |μ, Σ).RKHS theory: A positive-definite, symmetric function k: ^ × ^ → ℝ is called a kernel(note the distinction from probability kernels) if for all N ∈ ℕ, a^, … , a^ ∈ ℝ, andx^, … , x^ ∈ ^ we have ∑^ ^^^ ∑^ ^^^a^ a^k^x^, x^^ ≥ 0. Note that a positive-definite kernelgives rise to a positive-definite Gram matrix K : = [k(x , x^)] Given a kernel k onnon-empty set ^, there exists a corresponding unique reproducing kernel Hilbert space (RKHS) ℋ^of functions. ℋ^is equipped with an inner product ^ ⋅ , ⋅ ^ℋ^with thereproducing property such that for any function f: ^ → ℝ, f ∈ ℋ^ and x ∈ ^ we havef(x) = ^f, k( ⋅ , x)^ where k( ⋅ , x): ^ → ℋ is a real-valued function fo ^ℋ^ ^ r which k(x, x ) =^k( ⋅ , x), k( ⋅ , x^)^ℋ^ for all x, x^ ∈ ^. The inner product induces the norm ||f||ℋ^: = ^^f, f^ℋ^. A kernel k is called translation-invariant if it can be expressed as a functionof distance between points x, x^ ∈ ^, i.e., there exists some function Δ: ^ → ℝ such thatk(x, x^) = Δ(x − x^).Tautz & Schuhmacher Law MPI1117P12WO June 3, 202517 / 39Dynamical systems: Consider a system M with a state space ^ ⊂ ℝ^ of dimension n ∈ℕ and dynamicsM: x^^^ = f(x^) + v^, v^ ∼ ^( ⋅ |0, σ^^I^), (1)where x^ ∈ ^ denotes the system state at the time instance t ∈ ℕ ∪ 0. The independent,identically distributed (i.i.d.) process noise v^is zero-mean Gaussian with covariance matrix σ^^I^. Note that this implies that v^is dimension-wise independent. The transitionfunction f: ^ → ^ is unknown and we assume that it lives in the RKHS of sometranslation-invariant kernel k, i.e., f ∈ ℋ^. The covariance matrix of the process noiseis known.Reachability specifications: Given a target set ^ ⊂ ^, the system satisfies thereachability specification with probability at least ∈ [0,1] if the probability that thetrajectories remain in ^ until eventually reaching the target set ^ is greater than orequal to p^. This is denoted by ℙ(M ⊨ ψ^) ≥ p^. For the purpose of this work, we limitourselves to reachability specifications, but the provided verification algorithm can be applied to any linear temporal logic specification by adding a preprocessing step on the graph structure of the constructed abstract model (Baier and Katoen, 2008; Haddad and Monmege, 2018).Problem statement: Let ^ : = {x , y ^^ ^ ^}^^^ ⊂ ^ × ^ be a set of N ∈ ℕ data samples fromthe latent true system in (1) such that y^ : = f(x^) + v^ are output samples with i.i.d. noise We may write ^^ = (X^, Y^), where X^ : = [x^, … , x^] and Y^: = [y^, … , y^]^.Problem 1. Let data ^^ from the system in (1) with unknown function f ∈ ℋ^ be given.For an infinite-horizon reachability specification ψ^, verify that the system satisfies with probability at least p^ ∈ [0,1].To address this problem, we construct a function estimator ^ f: ^ → ^ based on ^^ thatapproximates the unknown function f ∈ ℋ^ in (1). In particular, we specify a kernel k^characterizing the function class of ^ f ∈ ℋ^^. To be able to address a wealth of functionsf and perform the function approximation in a function space ℋ^^ ⊅ ℋ^, we raise thefollowing assumption.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202518 / 39Assumption 1. Let the unknown function f ∈ ℋ^ and the function estimator^ f ∈ ℋ^^, withknown kernels k, k^: ^ × ^ → ℝ and complexity bound B ≥ ||f||ℋ^ ≥ 0 .Moving forward, we use ‘^’ to denote symbols referring to the estimated model.Remark 1. We raise Assumption 1 to allow for greater flexibility in choosing the kernelk^ used for estimation. In particular, we will use a discontinuous kernel. To still be ableto address more arbitrary functions, such as continuous functions living in the RKHS of the squared exponential kernel, we need to quantify the error introduced by constructing an estimator in a different RKHS. System Identification Using Binary-Tree Gaussian Processes: In this section, we give a brief overview of the classical GP regression framework, followed by an examination of GPs based on the binary-tree kernel. For the latter, we characterize the underlying RKHS and conclude with a statistical error bound on its predictions. Gaussian Process Regression: A Gaussian process (GP) is a parameter-free regression model that allows the user tolearn input-output mappings of the form^ f^: ^ → ℝ, d ∈ 1, … , n, from empirical data ^^ ={x^, y^} ^^^^ , y^ : = [y^,^, … , y^,^]^, such that the prediction errors ∑^^^^ | |y^,^ − ^ f^(x^)|| areminimized for all dimensions d (Rasmussen and Williams, 2018). Recall that, asspecified in (1), we assume that the process noise on the samples y^is dimension-wise independent, which allows us to use individual GPs for each output dimension. This is without loss of generality, and the results of this paper are applicable to non-diagonalcovariance matrices by applying a linear transformation on the data gathered from thesystem. In the following, we may drop the subscript d indicating the dimension forsimplicity. As an infinite-dimensional generalization of the multivariate normal distribution, GPs are particularly popular for their ability to capture statistical information about the latentTautz & Schuhmacher Law MPI1117P12WO June 3, 202519 / 39mapping. The following definition is analogous to Definition 2.2 by Kanagawa et al.(2018).Definition 1 (Gaussian process (GP)). Let k^: ^ × ^ → ℝ be a (symmetric, positive-definite) kernel or covariance function and m^ ^: ^ → ℝ. A random function ^ f^: ^ → ℝ isa GP, denoted by ^ f^ ∼ ^^^m^ ^, k^^, if for any finite X^ : = [x^, … , x^], x^ ∈ ^, N ∈ ℕ, thevector ^ f : = [^f ^^(x^)]^^^ has the distribution ^ f^,^ ∼ ^,^, K^^^, with ^,^: = [m^ ^(x^)] ^^^^ and Gram matrix K^^ : = [k^(x^, x^)]^^,^^^.The goal of GP regression is to generate predictions of the output at test points x ∈ ^that are not captured by the data ^^. Even though we do not assume that the unknowndata-generating function f is a sample from the GP with kernel k^ (see Assumption 1),we can construct a so-called posterior GP and obtain a probabilistic bound on itsprediction error (e.g., via Theorem 2 by Chowdhury and Gopalan (2017) or the resultsby Fiedler et al. (2021)).Definition 2 (Posterior GP). Let ^^ = (X^, Y^) from the system (1). Then, the posteriordistribution of f^(x) given the data ^^ is ^ f^(x) ∼ ^ ^ ⋅  characterized by the posterior mean and covariance functions Based on Definition 2, predictions of the next state of the system (1) for the currentstate x ∈ ^ can be generated via y^ (x) ∼ ^^ We establish probabilistic bound on the prediction error in the following after introducing a specific type of GP in the upcoming subsection. Binary-Tree Gaussian Processes: We focus on GPs based on the binary-tree (BT) kernel, which is defined on a finite space of binary strings and wields a range of practical advantages that renders itTautz & Schuhmacher Law MPI1117P12WO June 3, 202520 / 39 particularly interesting for formal approaches (Cohen et al., 2022b). Let us define amap Ψ^: ^ → ^^ encoding every state x ∈ ^ from the continuous space ^ as a binarystring s ∈ ^^ in the finite space of binary strings of length q. Note that this implies apartitioning of the space ^ into 2^ disjoint sets, thus q ∈ ℕ is called the bit-depth orprecision.In the following, we denote the set of states represented by a string s ∈ ^^ as Ξ^(s) ⊂^, i.e., we have Ψ^(x) = s for all x ∈ Ξ^(s). It follows that the state space is the unionof the partitions, i.e., ^ =∪^∈^^ Ξ^(s) for any precision q ∈ ℕ. For a set A ⊂ ^, we maywrite Ψ^(A) ⊂ ^^to denote its projection onto ^^. To simplify the notation, we maywrite Ψ^(x) : = [Ψ^(x)]^^, x ∈ ^, i ∈ ℕ, i.e., the operator that takes the first i elements ofthe binary string Ψ^(x). With this, we define the BT kernel formally.Definition 3 (Binary-tree (BT) kernel). Given a map Ψ^: ^ → ^^, the BT kernelk^ ^: ^ × ^ → ℝ of bit-depth or precision q ∈ ℕ is defined as with weight coefficients w^ ∈ ℝ such that w^ = 1.Intuitively, the BT kernel k^ ^ assigns its covariance mass based on whether pointsx, x^ ∈ ^ fall into the same partitions. The final covariance value is a weighted sum ofthis coincidence on all discretization levels up to precision q. It stands out that, in contrast to well-known alternatives such as the squared-exponential or Matérn kernel, the BT kernel is discontinuous since it is translation-variant. Hence, it producesdiscontinuous functional mappings when used in the context of GP regression. Fig. 7depicts the piecewise-constant posterior mean μ^ ^ and double standard deviation μ^^ ±2σ^^ as functions of state of a GP utilizing a BT kernel of precision q = 2, i.e. the BTGPcontains functions that are piecewise constant over the partitions of ^, wherein theposterior mean μ^ ^ and double standard deviation μ^ ^ ± 2σ^^ are displayed as apiecewise-constant graph and a shaded area plot, respectively. For this example withchosen precision q = 2, the posterior mean takes values μ^^ ∈ {y^^ , … , y^^}.. Movingforward, we will call such a GP a binary-tree Gaussian process (BTGP), which yieldsTautz & Schuhmacher Law MPI1117P12WO June 3, 202521 / 39 inherently discrete abstractions / representations of the underlying system.Furthermore, we may decide to write them as explicit functions of s ∈ ^^, i.e., μ^^(s)and σ^^^ (s). We note that although function samples from the BTGP (or any GP in general) do not live in the RKHS of the associated kernel almost surely, the resulting posterior mean function does (Kanagawa et al., 2018, Sec.4). Remark 2. The BT kernel violates the conditions for a continuous RKHS (Berlinet and Thomas-Agnan, 2011, Theorem 17, condition a)). Despite being discontinuous, the BTkernel is still a valid positive-definite kernel as shown in Proposition 1 by Cohen et al.(2022b) and therefore by the Moore-Aronszajn theorem acts as a reproducing kernel for a unique RKHS (Berlinet and Thomas-Agnan, 2011, Theorem 3), which contains only discontinuous functions.Remark 3. An implication of using the BT kernel of precision q in GP regression is thatthe resulting Gram matrix K^ ^ ^^ will be capped at a size of at most 2 × 2 . This is sinceM ∈ ℕ data samples {x^, y^} ^^^^ can be replaced by a modified dataset {x^, y^}^∈^^, wherex^ ∈ Ξ^(s) is a representative state and y^ y with A(s) : = (averaging overy^’s with x^being in the same partition set). The BTGP contains functions that are piecewise constant over the partitions of ^. Theposterior mean ^̂^ and double standard deviation ^̂^ ± 2^^^ are displayed as apiecewise-constant graph and an area plot in light blue, respectively. For this examplewith chosen precision ^ = 2, the posterior mean takes values ^^̂ ∈ ^^^, … , ^^^.The following theorem establishes a closed-form representation of the functions in the RKHS of the BT kernel.Theorem 1 (BT function space). Let k^ ^: ^ × ^ → ℝ be the BT kernel in (4). Then, itsassociated RKHS contains only functions f ∈ that can be represented for some y ∈ ℝ. The resulting finite dimensional RKHS (of dimensioncharacterized by the piecewise-constant feature map ϕ^ (x) : = [ w^1(x ∈ Ξ^(s))]^∈^^,^∈^,^,…,^. Tautz & Schuhmacher Law MPI1117P12WO June 3, 202522 / 39 A measure for the complexity of the functions in a given RKHS is the associated RKHS norm. For the BT kernel, a bound on the RKHS norm can be computed as follows.Corollary 1 (BT RKHS norm bound). Let k^ ^: ^ × ^ → ℝ be the BT kernel in (4) andℋ^^^ its associated RKHS. For all functions f ∈ ℋ^^^ of the form (5), we have that Statistical Error Bound To use GP regression in the context of formal approaches, we must quantify the errorbetween the unknown data-generating function f ∈ ℋ^ in (1) and the GP posteriormean μ^^ ∈ ℋ^^^ in (2). To be able to address a wealth of functions f beyond thosediscontinuous functions living in the RKHS ℋ^^^associated with the BT kernel, weleverage Assumption 1 to quantify the additional error introduced by the erroneousfunction space. The following theorem takes inspiration from Theorem 2.1 by Hsu etal. (2012) and Proposition 2 by Fiedler et al. (2021), where latter is itself an adaptationof Theorem 2 by Chowdhury and Gopalan (2017). As mentioned before, the outputdimensions — indexed by d — are modeled by one GP each.Theorem 2 (Piecewise-constant error bound). Consider the BT kernel k^ ^ in (4) withprecision q ∈ ℕ and the associated RKHS ℋ^^^. Let the unknown function f ∈ ℋ^ witha known translation-invariant kernel k that has a constant c satisfying c^ ≥ k(x, x) ≥ 0for all x ∈ ^, and a complexity bound B ≥ ||f||ℋ^ ≥ 0 (Assumption 1). Given data ^^ =(X^, Y^) and a constant δ ∈ (0,1), we have for each dimension d ∈ 1, … , n for thecorresponding BTGP posterior mean μ^^,^(s) ∈ ℋ^^^ in (2) that with the error ε^(s) : = ε^,^(s) + ε^,^(s) + ε^,^(s) given byε^,^(s) Tautz & Schuhmacher Law MPI1117P12WO June 3, 202523 / 39 where C^^ : = ^K^^ + σ^^I^^, k^ ^(s) : = [k(x^, Ψ^(s))]^^^^, the posterior covariance σ^^^,^ (s) in(3), Σ : = K^ ^C^^^^, k^^^ (s) : = sup^∈^(^)k(x^, x) for some fixed set of representative points Remark 4 (Approximating a continuous kernel). The performance benefits of GP regression via the BT kernel stem from the fact that it constitutes a series of subsequent partitionings of the input space such that in every discretization level the precision of the previous partitioning is refined (or kept the same). For each level, functional mass represented by the corresponding weight coefficient is appliedwhenever two points fall into the same partition. As shown in Theorem 1 by Cohen etal. (2022b), this allows for an efficient computation of the Gram matrix and its inverse(computations for the inversion have a complexity of O(N) for every precision level i: = 1, … , q instead of O(iN), which overall leads to O(qN)).Remark 5. The theorem above establishes a finite-sample bound on the error between the posterior mean of the BTGP and the actual unknown dynamics f. The complexitybound B on ||f||ℋ^ enables us to provide guarantees even when f does not live in theRKHS of the BTGP, ℋ^^^, e.g., if f is continuous. This is called the ‘misspecified case’.In comparison to Fiedler et al. (2021), we quantify the introduced error ε^without computing a global upper bound on the error between the two kernels, which would beat its maximum in our case as k^ ^ does generally not approximate k. The first term inthe minimization of ε^,^is equivalent to prior works and captures the statistical error of the properly specified case, which is generally very conservative. We provide an alternative second term in an effort to reduce this error by exploiting the low- dimensional BT kernel Gram matrices. Nevertheless, note that whilst the error bounds^ε^,^ , ε^,^^ can be tightened by increasing the precision q, this is not the case for ε^,^.Determining the error bound ε^,^ requires us to compute Gram matrices via both k andk^ ^, reducing the scalability benefits of the BTGP. It is interesting to study thequantification of ε^without the need for constructing K^.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202524 / 39One might be tempted to use the BT kernel in Definition 3, or an analogously definedalternative kernel bearing the same properties, to approximate a continuous kernelsuch as the squared-exponential kernel for the limit case q → ∞, allowing for asimplification of Theorem 2. However, such a kernel does not exist. In fact, let k^ ^(x, x^)be such a kernel with precision q ∈ ℕ and consider x ∈ ^ to be fixed and x^ ∈ ^variable. Then, for q < q^ ∈ ℕ the kernels k^ ^ and k^ ^^ define two discrete distancemetrics where the exponential increase in the level of precision of k^ ^^ compared to k^ ^is focused solely on the neighborhood of x. Hence, escalating the level of precisionmerely increases the smoothness of the covariance in the neighborhood of x and hasno effect on the quantification elsewhere. In contrast to methods such as probabilistic coupling relations (Schön et al., 2023) where the error between the true system and the discrete abstraction is captured bycomputing local parameters bounding the local deviations, the error bound ε^(s) inTheorem 2 summarizes the error in a formula that similarly yields local deviations whenevaluated for individual partitions s ∈ ^^.Verification Approach: The verification approach we consider is based on the construction of a finite-stateabstraction of the dynamical system in (1) via the BTGP introduced in Sec. 3.2. For thelearned BTGPs, let ^^be a finite set of binary strings representing the states of thisabstraction. Recall that a BTGP defines a map Ψ ^^: ^ → ^ that maps every state x ∈ ^from the continuous state space ^ to a string s ∈ ^^ associated with a partition of ^and a representative state x^ ∈ ^. To capture the probabilistic ambiguity w.r.t.transitions between discrete partitions, recent approaches for verifying infinite-horizon specifications are based on interval Markov chains (IMCs), which can be learned from data via GP regression (Badings et al., 2023; Lavaei et al., 2022; Jackson et al., 2020). In this manuscript, we follow a similar approach, however, we perform GP regressionusing the BT kernel introduced in Sec. 3.2, that allows us to reduce the computationaleffort in constructing an IMC from data significantly.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202525 / 39Definition 4 (Interval Markov chain (IMC)). An IMC is a tuple M^ : = ^^, sinit , t,‾ t^, with afinite state space ^ of discrete states s ∈ ^; initial state sinit ∈ ^; and probability boundst,‾ t: ^ × ^ → [0,1] describing lower and upper bounds on the transition probability.Computing the Abstract IMC via BTGP:The abstract IMC M^ : = ^^, sinit , t,‾ t^ is constructed with the state space ^ : = ^^, and theinitial state sinit : = Ψ^(xinit) with xinit being the initial state of the system M. Recall thatwe use one BTGP to model the dynamics for each output dimension. We hence obtain use the same precision q ∈ ℕ togenerate a common partitioning ^ = Ψ^(^). Computing the discrete transition kernelst,‾ t from the constructed BTGPs is comparatively straightforward due to the piecewise-constant nature of the posterior mean and covariance functions of the BTGP. In otherwords, for every partition s ∈ ^, the consecutive state (for each dimension d ∈ 1, … , n)is distributed according to a Gaussian distribution whose mean is uncertain in the formof μ^^,^(s) + ^^(s) with error ^^(s) ∈ E^(s) : = [−ε^(s), ε^(s)] obtained via thecorresponding stochastic learning errors (Theorem 2). The Gaussian distribution hasthe covariance (see Fig. 6, right). We obtain lower and upper bounds on theprobability of transitioning from s to any s^ ∈ ^ ast(s^|s) =^m∈^i(n^)g(s, ^)   and   ‾t(s^|s) =^m∈^a(x^)g(s, ^), (8) with E(s) : = ∏^ ^^^ E^ (s), ^ = [^^] ^^^^ , and x^ = [x^^]^^^^. If the partitions of the BTGPsare selected to be rectangular, the integration in (9) and subsequent optimizations in(8) can be done separately for each dimension. In comparison, for a kernel that is notpiecewise-constant, e.g., the squared-exponential kernel in Fig. 6 (top), obtaining thebounds requires optimizing with respect to ^ the function^ g(s, ^) = optim(iz)e^ ^^ ^dx^ ^^^μ^^,^(x) + ^^(x), σ^^,^ (x)^ , (10)^∈^ ^^^∈^(^^) ^^^with optimize ∈ max, min for each pair (s, s^). This shows that the BTGP eliminates theadditional inner optimization, thus significantly simplifying the computational burden associated with constructing the IMC to one integration per bound per transition.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202526 / 39 Interval Iteration Algorithm for Verification: Once the abstract IMC is constructed via (8)-(9), the interval iteration algorithm from the work by Haddad and Monmege (2018) can be employed to compute lower and upper bounds on the reachability probability. The algorithm is based on the iterative update of two value functions forming an envelope or interval around the trueprobability of satisfying the specification starting from some initial partition. We deferthe details to Appendix E. With the next theorem, we provide a robust satisfaction probability for the latent truesystem in (1) satisfying a reachability specification as follows.Theorem 3 (Robust satisfaction probability). Consider the setup of Theorem 2. Let^ f^ σ^ ^^,^ ^, d ∈ 1, … , n be the posterior BTGPs based on data ^^ from thesystem M in (1). Let M^ be the corresponding IMC constructed as outlined in (8). Then,we have for an infinite-horizon reachability specification and M initialized at anystate xinit ∈ ^ that Vmaxare the lower and upper bounds on the value function obtained via the intervaliteration algorithm (Algorithm 1).Case Study: Consider the discrete-time nonlinear system evolving on the continuous state space ^ = [−10,10]^ ⊂ ℝ^ with τ^ = 0.5 and v^ ∼^( ⋅ |0, σ^^I^). The goal is to compute a lower bound on the probability of the systemreaching the target set ^ = [−3,3]^ ⊂ ^ based on N = 5000 observations ^^ from thesystem generated with σ^ = 3.16. We construct two BTGPs — one for each outputdimension — with a precision of q = 12. Training and evaluation of the BTGPs at therepresentative states takes 13 seconds on an Apple MacBook Pro M1. To compute anerror bound via Theorem 2, we choose the confidence bound δ : = 0.2, squared-exponential kernels k ^^(x, x ) : = c^^exp(−(x − x^)^M ^^ ^^ (x − x )) for d ∈ 1,2 withTautz & Schuhmacher Law MPI1117P12WO June 3, 202527 / 39hyperparameters c^ = 12, c^ = 7, and M^ = diag^l ^^^ where l^ = [4000,2500] and l^ =[500,2000], and dimension-wise complexity bounds B^ = 0.015, B^ = 0.006. For{x^}^∈^^, we choose the center points of the partitions. Due to its piecewise-constantnature, generating the IMC via (8) from the BTGPs and error bounds takes less than 3seconds. In comparison, abstraction of a comparable continuous GP via (10) takesaround 8 hours. The IMC is subsequently verified via the interval iteration algorithm,which converges after 9 seconds (ν = 10^^). The resulting lower bound on theprobability of reaching ^ is depicted as a function of the initial state in Fig. 8.Conclusion: In this non-limiting specific example, we showcased how to use binary-tree Gaussian processes (BTGPs) for learning model representations from data that are naturally ideal for finite-state abstraction and verification due to the piecewise-constant natureof the employed binary tree kernel. We also formulated the error bounds in the learning,which were used to generate the abstraction as an interval Markov chain. We are currently working on improving the computation of the error bounds and benchmarking the proposed approach against alternative data-driven methods. References: Abate, A., D’Innocenzo, A., Di Benedetto, M.D., and Sastry, S.S. (2008). Markov set- chains as abstractions of stochastic hybrid systems. In International Workshop on Hybrid Systems: Computation and Control, 1–15. Springer. Abate, A., Edwards, A., and Giacobbe, M. (2022). Neural abstractions. Advances in Neural Information Processing Systems, 35, 26432–26447. Adams, S., Lahijanian, M., and Laurenti, L. (2022). Formal control synthesis for stochastic neural network dynamic models. IEEE Control Systems Letters, 6, 2858– 2863.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202528 / 39 Badings, T., Romao, L., Abate, A., and Jansen, N. (2023). Probabilities are not enough: Formal controller synthesis for stochastic dynamical models with epistemic uncertainty. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, 14701– 14710. Baier, C. and Katoen, J.P. (2008). Principles Of Model Checking, volume 950. MIT Press. Berlinet, A. and Thomas-Agnan, C. (2011). Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media. Bogachev, V.I. (2007). Measure theory. Springer Science & Business Media. Cao, M.E., Bloch, M., and Coogan, S. (2022). Efficient learning of hyperrectangularinvariant sets using Gaussian processes. IEEE Open Journal of Control Systems, 1,223–236. Chowdhury, S.R. and Gopalan, A. (2017). On kernelized multi-armed bandits. In International Conference on Machine Learning, 844–853. PMLR. Cohen, M.H., Belta, C., and Tron, R. (2022a). Robust control barrier functions for nonlinear control systems with uncertainty: A duality-based approach. In 2022 IEEE 61st Conference on Decision and Control (CDC), 174–179. IEEE. Cohen, M.K., Daulton, S., and Osborne, M.A. (2022b). Log-linear-time Gaussian processes using binary tree kernels. Advances in Neural Information Processing Systems, 35, 8118–8129. Fiedler, C., Scherer, C.W., and Trimpe, S. (2021). Practical and rigorous uncertainty bounds for Gaussian process regression. In Proceedings of the AAAI conference on artificial intelligence, volume 35, 7439–7447. Haddad, S. and Monmege, B. (2018). Interval iteration algorithm for MDPs and IMDPs. Theoretical Computer Science, 735, 111–131.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202529 / 39 Haesaert, S., Nilsson, P., Vasile, C., Thakker, R., Aghamohammadi, A., Ames, A., and Murray, R. (2018). Temporal logic control of POMDPs via label-based stochastic simulation relations. ADHS. Haesaert, S. and Soudjani, S. (2020). Robust dynamic programming for temporal logic control of stochastic systems. IEEE Transactions on Automatic Control, 66(6), 2496– 2511. Harapanahalli, A., Jafarpour, S., and Coogan, S. (2023a). Contraction-guided adaptive partitioning for reachability analysis of neural network controlled systems. In 202362nd IEEE Conference on Decision and Control (CDC), 6044–6051. IEEE. Harapanahalli, A., Jafarpour, S., and Coogan, S. (2023b). Forward invariance in neural network controlled systems. IEEE Control Systems Letters. Hasanbeig, M., Abate, A., and Kroening, D. (2019). Logically-constrained neural fitted Q-iteration. In AAMAS. Hsu, D., Kakade, S.M., and Zhang, T. (2012). A tail inequality for quadratic forms ofsubgaussian random vectors. Electronic Communications in Probability, 17, 1–6.Jackson, J., Laurenti, L., Frew, E., and Lahijanian, M. (2020). Safety verification ofunknown dynamical systems via Gaussian process regression. In 202059thIEEE Conference on Decision and Control (CDC), 860–866. IEEE.Jafarpour, S., Harapanahalli, A., and Coogan, S. (2023). Interval reachability ofnonlinear dynamical systems with neural network controllers. In Learning for Dynamicsand Control Conference, 12–25. PMLR.Jagtap, P., Soudjani, S., and Zamani, M. (2020). Formal synthesis of stochasticsystems via control barrier certificates. IEEE Transactions on Automatic Control, 66(7), 3097–3110.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202530 / 39 Jiang, J., Zhao, Y., and Coogan, S. (2022). Safe learning for uncertainty-awareplanning via interval MDP abstraction. IEEE Control Systems Letters, 6, 2641–2646.Kanagawa, M., Hennig, P., Sejdinovic, D., and Sriperumbudur, B.K. (2018). Gaussian processes and kernel methods: A review on connections and equivalences. arXiv:1807.02582.Kazemi, M., Majumdar, R., Salamati, M., Soudjani, S., and Wooding, B. (2024a). Data-driven abstraction-based control synthesis. Nonlinear Analysis: Hybrid Systems, 52,101467.Kazemi, M., Perez, M., Somenzi, F., Soudjani, S., Trivedi, A., and Velasquez, A.(2024b). Assume-guarantee reinforcement learning. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, 21223–21231.Kazemi, M. and Soudjani, S. (2020). Formal policy synthesis for continuous-statesystems via reinforcement learning. In International Conference on Integrated FormalMethods, 3–21.Lavaei, A., Soudjani, S., Frazzoli, E., and Zamani, M. (2022). Constructing MDPabstractions using data with formal guarantees. IEEE Control Systems Letters, 7, 460–465. Lew, T. and Pavone, M. (2021). Sampling-based reachability analysis: A random settheory approach with adversarial sampling. In Conference on robot learning, 2055–2070. PMLR.Lopez, B.T. and Slotine, J.J.E. (2023). Unmatched control barrier functions: Certaintyequivalence adaptive safety. In 2023 American Control Conference (ACC), 3662–3668. IEEE.Tautz & Schuhmacher Law MPI1117P12WO June 3, 202531 / 39Quinonero-Candela, J. and Rasmussen, C.E. (2005). A unifying view of sparseapproximate Gaussian process regression. The Journal of Machine LearningResearch, 6, 1939–1959.Rasmussen, C.E. and Williams, C.K.I. (2018). Gaussian Processes for MachineLearning. The MIT Press.Salamati, A., Lavaei, A., Soudjani, S., and Zamani, M. (2024). Data-driven verificationand synthesis of stochastic systems via barrier certificates. Automatica, 159, 111323.Schön, O., van Huijgevoort, B., Haesaert, S., and Soudjani, S. (2023). Bayesian formalsynthesis of unknown systems via robust simulation relations. arXiv:2304.07428.Schön, O., Zhong, Z., and Soudjani, S. (2024). Data-driven distributionally robust safety verification using barrier certificates and conditional mean embeddings. In 2024American Control Conference (ACC), arXiv:2403.10497.Thorpe, A., Lew, T., Oishi, M., and Pavone, M. (2022). Data-driven chance constrainedcontrol using kernel distribution embeddings. In Learning for Dynamics and ControlConference, 790–802. PMLR.Wajid, R., Awan, A.U., and Zamani, M. (2022). Formal synthesis of safety controllersfor unknown stochastic control systems using Gaussian process learning. In Learningfor Dynamics and Control Conference, 624–636. PMLR.Zhang, Z., Ma, C., Soudijani, S., and Soudjani, S. (2024). Formal verification ofunknown stochastic systems via non-parametric estimation. International Conferenceon Artificial Intelligence and Statistics (AISTATS).Appendix A. Additional Figures on the Case Study: (Omitted)Appendix B. Proof of Theorem 1:We can write the BT kernel in (4) as the sum Tautz & Schuhmacher Law MPI1117P12WO June 3, 202532 / 39of kernels feature maps where we expanded each ith summand of k^ ^ w.r.t. all possible binary substrings oflength i. Note that the binary substrings s ∈ ^^ in κ^: ^ × ^ → [0,1] are unique and κ^ isclearly positive-definite. The associated RKHS ℋ^^is of dimension 2^. Since finite dimensional spaces are automatically complete, ℋ^^is simply the span of the functionsκ^(x, ⋅ ), ∀x ∈ ^. All functions f^ ∈ ℋ^^ are hence linear combinations of the form for arbitrary m ∈ ℕ, α^ ∈ ℝ, and x^ ∈ ^ (Berlinet and Thomas-Agnan, 2011). Note thatϕ^(x^) only places functional weight on the partition s = Ψ^(x^). Therefore, f^ ∈ ℋ^^ isalways of the form f^(x) = y^^ϕ^(x),with y^ ∈ ℝ^^. Since the BT kernel is a weighted sum of kernels κ^as shown in (B.1), we obtain f(x) = y^ϕ^ ^(x),with y^ : = [y^, … , y^] and feature map ϕ^ ( ^ ^ ^^ ^ ^ x) : = [ϕ^(x) , … , ϕ^(x) ] via the kernelsum rule (Berlinet and Thomas-Agnan, 2011, Theorem 5), which is equivalent to (5) .Note that k^ ^(x, x^) = ^^^(x), ϕ^^(x)^ℋ is (∑^2)-Appendix C. Proof of Corollary 1:We can write the BT kernel in (4) as the sum of (positive-definite) kernels Following similar reasoning as in the proof of Theorem 1, we find We compute the norm of f^ ∈ ℋ^^ given κ^ via Tautz & Schuhmacher Law MPI1117P12WO June 3, 202533 / 39where we used κ^(x^, x^) = 1. From Aronszajn’s sums of kernels theorem (Berlinet andThomas-Agnan, 2011, Theorem 5) we have ∀f ∈ concluding the proof.Appendix D. Proof of Theorem 2:Similar to the proof of Theorem 2 by Chowdhury and Gopalan (2017), we start bysplitting the error bound into one term related to the noise and one to the statisticalapproximation error. For the data ^^ = (X^, Y^), let Y^ : = [y^,^, … , y^,^]^and recognizethat we have Via Cauchy-Schwarz we obtain where we abbreviate C^^^^ : = [K^^ + σ^^I^]^^. We bound the noise term following similarsteps to Proposition 2 by Fiedler et al. (2021) to get^ 2^Nlog1 ^ + 1 ^+ 2log, δ δ which yields the first term in the minimization of ε^,^(s). As an alternative bound, wecan follow the same steps as in the proof of Theorem 2 by Chowdhury and Gopalan(2017) to derive that From Theorem 2.1 by Hsu et al. (2012) we have that^ ℙ^||ΣV ||^ ≤ tr(Σ) + 2(tr(Σ^)t) ^ + 2||Σ| ^^^ |t^ ≥ 1 − e ,where Σ : = k^ C^ ^^. Setting t : = l^ ^^ og ^^^ yields the second term in the minimization ofε^,^(s). Depending on the setting, one or the other bound will provide a less conservative error for ε^,^(s).Next, we select some fixed x^ ∈ ^ for all s ∈ ^ such that‾ f^(s) : = f^(x^). Expanding thesecond term in (D.1) with ‾ f^(s) and Cauchy-Schwarz yields for every s ∈ ^ that for allx ∈ Ξ(s)Tautz & Schuhmacher Law MPI1117P12WO June 3, 202534 / 39 The latter term is bounded from above via Let Φ^ : = [ϕ(x^), … , ϕ(x^)]^. With this and Cauchy-Schwarz, the prior term in (D.2) yields which reduces to ε^,^(s), concluding the proof. Appendix E. Interval iteration algorithm: In this section, we provide the details of the interval iteration algorithm by Haddad and Monmege (2018) (Algorithm 1). Inter alia, they establish the following result.Proposition 1 (Value convergence). For a given IMC M^ : = (^, sinit , t,‾ t) and an infinite-horizon reachability specification ψ∗^, let V: ^ → [0,1] be the optimal value functionsuch that ℙ(M^ ⊨ ψ^) = V∗(sinit). Then, the interval iteration in Algorithm 1 convergesto V∗ as i → ∞. Furthermore, for every iteration i ∈ ℕ ∪ {0} we have V ∗^(sinit) ≤ V (sinit) ≤V‾^(sinit). Algorithm 1. Interval iteration algorithmInput: IMC M^ = (^, sinit , t, t), specification ψ^, stopping threshold ν > 0for s ∈ ^ do end forTautz & Schuhmacher Law MPI1117P12WO June 3, 202535 / 39i ← 0, V^ ← 1, V^ ← 0while ^ V^ − V^^^ > ν dor, t ← solve min LP in (E.1)for s ∈ ^ do end for i← i + 1, V^ ← V^^^, V^ ← V^^^end whileVmin ← V^ (Take underapproximation)i ← 0, V^ ← 1, V^ ← 0while ^ V^ − V^^^ > ν dor, t ← solve max LP in (E.1)for s ∈ ^ do end whileVmax ← V^ (Take overapproximation)Output: Vmin, VmaxBefore giving a summary of Algorithm 1, let us define the following linear program (LP),where optimize ∈ {max, min}:Tautz & Schuhmacher Law MPI1117P12WO June 3, 202536 / 39 In lines 2 to 5 of Algorithm 1, we compute functional lower and upper bounds (r, r‾) onthe reward – the probability of reaching the projected target region Ψ^(^) ⊂ ^ in onetime-step starting from the state s ∈ ^. Analogously, (‾l, l) represent functional boundson the loss – the probability of leaving the state space ^. In line 7, theoverapproximation V^ and underapproximation V^ of the value function – the probabilityof eventually reaching the target starting from partition s – are initialized. Untilconvergence, we repeat the following two steps. First, for every s, s^ ∈ ^, we minimizethe LP in (E.1) in line 9 to find feasible rewards r and transition probabilities t(s^|s) thatcapture the worst-case behavior of the system. Then, we update V^and V^via theBellman equation in lines 11 and 12. The while-loop converges when the maximaldifference between the over- and under-approximation is less than the predefinedcloseness bound ν. After convergence, the under-approximation V^is taken as the lower bound Vminon the value function in line 16. The same procedure is repeated toobtain an overapproximation Vmax in lines 17 to 26, where the LP is maximized in line19. From Proposition 1 we get that the true latent value V∗ lies between the obtainedupper and lower bounds.Following from Theorem 2 and Proposition 1, Theorem 3 gives bounds on the robustsatisfaction probability.

Claims

Tautz & Schuhmacher Law MPI1117P12WO June 3, 202537 / 39 Claims1. Computer-implemented method (41) for generating a dynamics model (43) ofa system (42), the method (41) comprising:- generating the dynamics model (43), optionally including an IMC orIMDP, by determining unknown dynamics of the system (42) using a kernel-based method (431) designed for regression and / or classification,optionally BTGP regression, characterized in that generating the dynamics model (43) further includes:- computing (432) probability bounds for each transition or multipletransitions of the dynamics model (43) by obtaining transition probability distributions from the kernel-based method (431) that are amenable to a predefined model checking algorithm, and -determining (44) an abstraction error of the generated dynamics model(43) based on the computed probability bounds.

2. Computer-implemented method (41) according to claim 1, characterized inthat the system (42) is a continuous-time or discrete-time, potentiallyheterogeneous-space including discrete and / or continuous variables, stochastic or non-stochastic system (42).

3. Computer-implemented method (41) according to claim 1 or 2, characterizedin that the obtained transition probability distributions are of a piecewise-predefined, optionally piecewise-constant, periodic, or otherwise hierarchically or structurally correlated shape in order to be amenable to the predefinedmodel checking algorithm.

4. Computer-implemented method (41) according to any of claims 1 to 3,characterized in that the abstraction error characterizes a deviation of a behavior of the generated dynamics model (43) from the system (42).Tautz & Schuhmacher Law MPI1117P12WO June 3, 202538 / 395. Computer-implemented method (41) according to any of claims 1 to 4,characterized in that only finitely many observations of the system (42) are available for generating the dynamics model (43).

6. Computer-implemented method (41) according to any of claims 1 to 5,characterized in that not more than two integrations are performed for eachof the transitions to identify the worst-case probabilities for the respective transition.

7. Computer-implemented method (41) according to any of claims 1 to 6,characterized in that a kernel of the kernel-based method (431) isdiscontinuous and defined on a finite space of binary or non-binary strings.

8. Computer-implemented method (41) according to any of claims 1 to 7,characterized in that the generating of the model (43) includes directlylearning the model (43) from the observations of the system (42) withoutperforming explicit steps of discretization.

9. Computer-implemented method (41) according to any of claims 1 to 8,characterized in that the system (42) represents a control system, optionallydesigned to be used in safety critical environments.

10. A data processing device comprising means for carrying out the method (41)according to any of claims 1 to 9.

11. A computer program comprising instructions which, when the program isexecuted by a computer, cause the computer to carry out the method (41)according to any of claims 1 to 9.

12. A computer-readable medium comprising instructions which, when executedby a computer, cause the computer to carry out the method (41) according toany of claims 1 to 9.

Images