SySCoRe: Synthesis via Stochastic Coupling Relations

TL;DR

SySCoRe is a MATLAB toolbox that synthesizes controllers for stochastic systems with complex temporal logic specifications, supporting nonlinear dynamics and providing formal robustness guarantees with efficient computation methods.

Contribution

It introduces a novel approach supporting nonlinear dynamics and infinite-horizon specifications, with probabilistic coupling and tensor-based computations for robust controller synthesis.

Findings

Provides non-trivial lower bounds for infinite-horizon properties

Supports complex co-safe temporal logic specifications

Efficiently computes transition probabilities using tensor representation

Abstract

We present SySCoRe, a MATLAB toolbox that synthesizes controllers for stochastic continuous-state systems to satisfy temporal logic specifications. Starting from a system description and a co-safe temporal logic specification, SySCoRe provides all necessary functions for synthesizing a robust controller and quantifying the associated formal robustness guarantees. It distinguishes itself from other available tools by supporting nonlinear dynamics, complex co-safe temporal logic specifications over infinite horizons and model-order reduction. To achieve this, SySCoRe generates a finite-state abstraction of the provided model and performs probabilistic model checking. Then, it establishes a probabilistic coupling to the original stochastic system encoded in an approximate simulation relation, based on which a lower bound on the satisfaction probability is computed. SySCoRe provides non-trivial lower bounds for infinite-horizon properties and unbounded disturbances since its computed error does not grow linearly in the horizon of the specification. It exploits a tensor representation to facilitate the efficient computation of transition probabilities. We showcase these features on several benchmarks and compare the performance of the tool with existing tools.

Tables

Table 1. TABLE I: Main functions of SySCoRe for steps (1)-(6), with optional steps (2a) and (2b).

Step	Function
(1) Translate the specification	TranslateSpec
(2) Finite-state abstraction	FSabstraction
(2a) Piecewise-affine approx.	PWAapproximation
(2b) Model-order reduction	ModelReduction
(3) Similarity quantification	QuantifySim
(4) Synthesize a controller	SynthesizeRobustController
(5) Control refinement	RefineController
(6) Deployment	ImplementController

Table 2. TABLE II: An overview of the different benchmarks and their total computation time in seconds (s) and memory usage in megabyte (MB). The details of the computation times for each step are reported in Table  III . Dim. and Comp. are abbreviations for Dimension and Computation , respectively. The size of the specification refers to the number of states of the DFA.

Benchmark	System		MOR	Specification			Comp. time (s)	Memory (MB)
	Dynamics	Dim.		Type	Time horizon	Size
Running example	Linear	2	No	Reach-avoid	Unbounded	3	7.94	27.53
Package delivery	Linear	2	No	Reach-avoid	Unbounded	3	11.02	133.4
Van der Pol oscillator	Nonlinear	2	No	Safety, reachability	Unbounded	3	3191.6	178.83
Building automation	Linear	7	Yes	Safety	Bounded	8	122.05	5365.6

Table 3. TABLE III: Computation times for the different steps (1)-(6) in seconds and as percentage of the total runtime. Steps (1)-(6) correspond to (1) translating the specification, (2) finite-state abstraction, (3) similarity quantification, (4) synthesizing a controller, (5) control refinement, and (6) deployment. Step (5) is almost instantaneous ( ≈ \approx 0.001 s), therefore, we have taken the computation times of steps (5) and (6) together.

	Step (1)	Step (2)	Step (3)	Step (4)	Step (5) and (6)	Total
Running example	0.259s (3.26%)	1.316s (16.57%)	5.590s (70.38%)	0.507s (6.39%)	0.204s (2.57%)	7.944s (100%)
Package delivery	0.284s (2.58%)	1.657s (15.04%)	6.193s (56.2%)	1.708s (15.5%)	0.702s (6.37 %)	11.02s (100%)
Van der Pol oscillator	0.590s (0.02%)	1440.1s (45.1%)	1748.6s (54.8%)	2.854s (0.09%)	1.417s (0.04%)	3191.6s (100%)
Building automation	0.361s (0.30%)	4.80s (3.94%)	67.92s (55.7%)	37.33s (30.6%)	9.19s (7.53%)	122.05s (100%)

Table 4. TABLE IV: Results of the benchmarks for different tools. Here, n.a. means that a tool is not applicable and n.s. means that the current version of the tool does not natively support the computations on the benchmark, but that we do not see fundamental limitations hindering such an extension. To compare the tools we exclude the deployment of the controller (step (6)), since this step is not performed by the other tools.

Tool	Run time (s)
AMYTISS	n.s.
FAUST	n.s.
SReachTools	n.a.
StocHy	n.s.
SySCoRe	10.319

Table 5. (a) Package delivery benchmark.

Tool	Run time (s)
AMYTISS	n.s.
FAUST	n.s.
SReachTools	n.a.
StocHy	n.s.
SySCoRe	10.319

Table 6. (b) Van der Pol benchmark.

Tool	Run time (s)
AMYTISS	n.s.
FAUST	n.a.
SReachTools	n.a.
StocHy	n.a.
SySCoRe	3190.2

Table 7. (c) Building automation benchmark.

Tool	Run time (s)	Max. reach probability
AMYTISS	312.14	$\approx 0.8$
FAUST	n.s.	n.s.
SReachTools	4.59	$\ge 0.99$
StocHy	$\ge$ 335.876	$\ge 0.8\pm 0.23$
SySCoRe	112.86	$\ge 0.9035$

Equations

$$M:\{\begin{array}{l}{x}_{t+1}=f(x_t,u_t)+B_w{w}_t\\ y_t=C{x}_t,\forall t\in \left\{0,1,2,\dots \right\},\end{array}$$

$$\{\begin{array}{l}{x}_{t+1}=A_i{x}_t+B_i{u}_t+a_i+B_{w,i}{w}_t+{\kappa }_t\text{ for }{x}_t\in P_i\\ y_t=C{x}_t,\end{array}$$

$$\{\begin{array}{l}{x}_{t+1}=A{x}_t+B{u}_t+B_w{w}_t\\ y_t=C{x}_t,\end{array}$$

$${\varphi }_{park}=\neg p_2\mathsf{U}{p}_1,$$

$$\mathbb{P}(M\times C\models \varphi )\ge \rho .$$

$$u_{r,t}={\hat{u}}_t$$

$$u_{r,t}={\hat{u}}_t+K(x_{r,t}-{\hat{x}}_t)$$

$$u_t=u_{r,t}+Q{x}_{r,t}+K_{MOR}(x_t-P{x}_{r,t}).$$

$$M_r:\{\begin{array}{l}{x}_{r,t+1}=A_r{x}_{r,t}+B_r{u}_t+B_{rw}{w}_{r,t}\\ y_{r,t}=C_r{x}_{r,t},\end{array}$$

$$\text{SwapAboveDisplaySkip}{\mathcal{R}}_{MOR}:=\left\{(x_r,x)\in {\mathbb{X}}_r\times \mathbb{X}\mid ||x-P{x}_r|{|}_{D_r}\le {ϵ}_r\right\}$$

$$\mathcal{R}:=\left\{(\hat{x},x_r)\in \hat{\mathbb{X}}\times {\mathbb{X}}_r\mid ||x_r-\hat{x}|{|}_D\le ϵ\right\},$$

$${\mathbf{T}}^{\hat{u}}(V)(\hat{x},q):=\mathbf{L}\left({\mathbb{E}}_{\hat{u}}\left(\underset{q^{+}\in Q^{+}}{\min }\max \{1_{Q_f}(q^{+}),V({\hat{x}}^{+},q^{+})\}\right)-\delta \right),$$

$$Q^{+}(q,{\hat{y}}^{+}):=\left\{{\tau }_{{\mathcal{A}}_{\varphi }}(q,L(y^{+}))\mid ||y^{+}-{\hat{y}}^{+}||\le ϵ\right\},$$

$$\text{SwapAboveDisplaySkip}{x}_{1,t+1}$$

$$x_{2,t+1}$$

Full text

I Introduction

The design of provably correct controllers is crucial for the development of safety-critical systems such as autonomous vehicles and smart energy grids [25, 6]. To this end, methods for synthesizing controllers for dynamical systems that are guaranteed to satisfy temporal logic specifications have gained an increasing amount of attention in the control community [8, 37, 9, 24]. Besides establishing the theory underlying these methods, it is equally important to develop tools that facilitate their application. For stochastic systems, a collection of tools that can perform formal controller synthesis is already available. A subset of these tools include in alphabetical order: AMYTISS [23], FAUST [36], hpnmg [17], HYPEG [30], Mascot-SDS [27], the Modest Toolset [15], ProbReach [34], SReachTools [41], and StocHy [11]. A complete list of these tools with their descriptions and capabilities can be found in the ARCH Competition Report (stochastic category) [1]. These tools perform the computations either using analytical methods or employing statistical model checking. The approaches in the analytical methods can further be divided into abstraction-based [36, 11, 23, 27] and abstraction-free techniques [41, 19]. Abstraction-free techniques are generally less prone to suffering from the curse of dimensionality, however, they are often limited to simple invariance and reachability specifications. In contrast, abstraction-based tools can be applied to a breath of systems and specifications. A survey on formal verification and control synthesis of stochastic systems is given in [24].

SySCoRe contributes to the category of tools that employ analytical abstraction-based methods. It is a MATLAB toolbox applicable to stochastic nonlinear systems with a possibly unbounded disturbance. Furthermore, it can perform the controller synthesis to satisfy arbitrary co-safe specifications that can have unbounded time horizons. To this end, it uses the $(\epsilon,\delta)$-approximate simulation relation provided in [14], that explicitly designs the coupling between the continuous-state model and its (reduced) finite-state abstraction [39]. Hence, SySCoRe extends the capabilities of the current tools by considering properties that are unbounded in time and by considering systems with an unbounded disturbance.

SySCoRe is a comprehensive toolbox for temporal logic control of stochastic continuous-state systems, implementing all necessary steps in the control synthesis process. Moreover, it supports model-order reduction in the abstraction process with formal error quantification quarantees, which makes it applicable to a larger classes of systems. To increase its computational efficiency, SySCoRe performs computations based on tensors and sparse matrices. Furthermore, computations based on efficient convex optimizations for polytopic sets are implemented where possible. The tool is developed with a focus on ease of use and extensibility, such that it can easily be adapted to suit individual research purposes. The development of SySCoRe is a step towards solving the tooling need for temporal logic control of stochastic systems as it expands both the class of models and the class of specifications for which abstraction-based methods can provide controllers with formal guarantees.

This tool paper is organized as follows. We discuss in Section II the temporal logic control problem and the set-up in SySCoRe. We then give an overview of SySCoRe in Section III by introducing the associated functions and classes. Section IV discusses multiple benchmarks that show the capabilities of SySCoRe and how it compares to existing tools. We end the paper with a summary and a discussion of possible extensions. Throughout, we give the core functions of SySCoRe in framed white boxes and example code in gray boxes.

II Temporal logic control

The main purpose of SySCoRe is to perform the complete control synthesis procedure in abstraction-based temporal logic control. It is applicable to discrete-time models with a possibly unbounded stochastic disturbance and synthesizes a controller for satisfying co-safe linear temporal logic specifications that may have an unbounded time horizon. The computational approach is based on the theory of approximate simulation relations [14], the coupling between models [14, 39] and robust dynamic programming mappings [13]. In this section, we introduce the class of models and specifications handled by SySCoRe, and show how to set up the problem. Furthermore, we provide a high-level description of the theory underlying the implementations in SySCoRe.

II-A Problem parameters

Model. We consider discrete-time systems described by stochastic difference equations

[TABLE]

with state $x_{t}\in\mathbb{X}$, input $u_{t}\in\mathbb{U}$, (unbounded) stochastic disturbance $w_{t}\in\mathbb{W}$, measurable function $f:\mathbb{X}\times\mathbb{U}\rightarrow\mathbb{X}$, and matrices $B_{w}$ and $C$ of appropriate sizes.

To handle nonlinear systems of the form (1) we perform a piecewise-affine (PWA) approximation that yields a system described by

[TABLE]

with $P_{i}$ a partition of $\mathbb{X}$ and $\kappa_{t}\in\mathscr{K}_{i}$ the error introduced by performing the PWA approximation. For ease of notation, we denote the state-dependent error $\kappa_{x_{t}}$ as $\kappa_{t}$. Furthermore, $A_{i},B_{i},B_{w,i}$ and $a_{i}$ are matrices of appropriate sizes. Details of temporal logic control for nonlinear stochastic systems via piecewise-affine abstractions can be found in [40]. Besides nonlinear systems, we also consider the special case of linear time-invariant (LTI) systems:

[TABLE]

with $A$ and $B$ matrices of appropriate sizes.

Remark 1

This first release of SySCoRe assumes the disturbance $w_{t}$ has unbounded Gaussian distribution $w_{t}\sim\mathcal{N}(0,I)$. The implementation for other classes of distributions is under way and will be included in the future release of the tool. Note that the assumption of standard Gaussian distribution with zero mean and identity covariance matrix is without loss of generality since any system (1)-(3) with disturbance $w\sim\mathcal{N}(\mu,\Sigma)$ can be rewritten to a system in the same class with an additional affine term [5].

To specify the model, that is a nonlinear system (1), a PWA system (2) or an LTI system (3), we have developed the classes NonLinModel, PWAmodel, and LinModel, respectively. The state space, input space, and the sets needed for defining the specification should be defined in these class descriptions.

Running example

Consider a two-dimensional (2D) case study of parking a car with dynamics of the form (3) with $A=0.9I_{2}$, $B=0.7I_{2}$, and $B_{w}=C=I_{2}$. Furthermore, we have state space $\mathbb{X}=[-10,10]^{2}$, input space $\mathbb{U}=[-1,1]^{2}$, and disturbance $w\sim\mathcal{N}(0,I_{2})$. After specifying the matrices $A,B,C,B_{w}$, and setting the values for the disturbance $w$ with mean mu and covariance matrix sigma equal to zero and identity respectively, we can initialize a model in SySCoRe as follows:

\MakeFramed\FrameRestore

\endMakeFramed

The state and input spaces are defined using Polyhedron from the multi-parametric toolbox (MPT3) [16] as follows.

\MakeFramed\FrameRestore

⬇

% Define bounded input spacesLTI.U = Polyhedron(combvec([-1,1],[-1,1])');

\endMakeFramed

Specifications. In SySCoRe, we consider formal specification written using co-safe linear temporal logic (scLTL) [9, 20], which consists of atomic proposition (AP) $AP=\left\{p_{1},p_{2},\dots,p_{N}\right\}$ that are either true or false. To connect the system and the specification, we label the output space of the system, such that we can relate the trajectories of the system $\boldsymbol{y}=y_{0},y_{1},y_{2},\dots$ to the atomic propositions of the specification $\phi$.

Running example cont’d

For the 2D car park, we consider reach-avoid specification $\phi_{park}$ with the region to reach $P_{1}$ and with avoid region $P_{2}$. First, we define the regions

\MakeFramed\FrameRestore

⬇

= Polyhedron([4, 0; 4, 4; 10, 4; 10 0]);

\endMakeFramed

and add them to the system object:

\MakeFramed\FrameRestore

⬇

% Propositions corresponding to the regionssLTI.AP = {'p1', 'p2'};

\endMakeFramed

Implicitly, this means that states inside regions P1 and P2 are labeled using the corresponding atomic propositions 'p1' and 'p2', respectively. Now, we can write the scLTL specification

[TABLE]

using the syntax from [12] as follows

\MakeFramed\FrameRestore

\endMakeFramed

Denote the system $M$ under the controller $C$ by $M\times C$ as in [37]. The goal is to synthesize a controller $C$, such that the controlled system satisfies an scLTL specification $\phi$, denoted as $M\times C\models\phi$. Since we consider stochastic systems, we compute the satisfaction probability denoted as $\mathbb{P}(M\times C\models\phi)$. This goal is formulated mathematically next.

Problem statement. Given model $M$, scLTL specification $\phi$, and probability threshold $\rho\in(0,1)$, design controller $C$ such that

[TABLE]

SySCoRe automatically synthesizes a controller by maximizing the right-hand side of (5) on a simplified abstract model and makes the computations robust with respect to the abstraction errors. It provides a robust lower bound on the satisfaction probability, which then can be used by the user to compare with the probability threshold $\rho$.

II-B Stochastic coupling relations for control synthesis

To solve the above problem, we use an abstraction-based approach and the dynamic programming mappings from [13], which allows us to consider infinite-horizon properties. More specifically, the abstraction-based temporal logic control implemented in SySCoRe has six main steps, namely (1) translating the specification to an automaton, (2) constructing a (reduced) finite-state abstraction, (3) quantifying the similarity, (4) synthesizing a controller, (5) control refinement, and (6) deployment.

As visualized in Figure 1, we start from a temporal logic specification that expresses the desired behavior of the controlled system and translate it to an automaton (see top layer). A finite abstract model $\hat{M}$ of the system is also constructed (step 2). For this abstract model $\hat{M}$, its bounded deviation from the original model can be quantified using simulation relations (step 3) [14, 39]. Computing these bounds is based on an efficient invariant set computation formulated as an optimization problem constrained by a set of parameterized matrix inequalities [39]. Based on the automaton, an abstract controller over the abstract model can be synthesized. In step 4, we synthesize an abstract controller $\hat{C}$ and compute the robust satisfaction probability. The robust satisfaction probability takes the deviation bounds computed in step 3 into account and gives a lower bound on the actual satisfaction probability.

To compute the robust satisfaction probability and to synthesize an abstract controller $\hat{C}$, SySCoRe solves a reachability problem over the abstract system combined with the automaton corresponding to the specification. This reachability problem is then solved as a dynamic programming problem. It is shown in [13] that leveraging the deviation bounds from step 3, the controller for the abstract model can be refined to the original continuous-state model while preserving the guarantees. To construct this controller $C$, SySCoRe refines the abstract controller in step 5. The resulting controller $C$ is a policy that can be represented with finite memory. Finally, SySCoRe deploys the controller on the model (step 6). It is important to note that the abstraction step (step 2 in Figure 1) can additionally contain model-order reduction or piecewise-affine approximation, which shows the comprehensiveness of SySCoRe enabled by establishing coupled simulation relations.

The next section gives a complete overview of the toolbox and specifies how each of the steps from Figure 1 is implemented.

III Toolbox overview

After setting-up the problem by specifying the system using the classes NonLinModel, PWAmodel or LinModel, and the specification as an scLTL formula, we continue with the steps illustrated in Figure 1. Each step corresponds to a specific function as in Table I. Note that the abstraction step may have multiple (formal) approximation stages depending on the type of the model or its dimension.

III-A Translating the specification

For control synthesis, the scLTL specification is written as a deterministic finite-state automaton (DFA) [9]. Examples of such DFAs are given in Figure 2. We use the tool LTL2BA111Tool available at http://www.lsv.fr/~gastin/ltl2ba/index.php to translate an scLTL specification, which constructs a non-deterministic Büchi automaton for a general LTL specification [12]. Additionally, we check whether the given formula is written using scLTL (instead of full LTL) and then (if possible) rewrite the non-deterministic Büchi automaton to a DFA. This step is based on powerset conversion [31] that is used to convert a nondeterministic finite-state automaton to a DFA. The complete translation from an scLTL specification to a DFA is implemented in the function TranslateSpec.

⬇

% Translate an scLTL formula to a DFAA = TranslateSpec(formula, AP);

The input formula is given using the syntax of LTL2BA in [12].

Running example cont’d

For the 2D car park, we consider the reach-avoid specification $\phi_{park}$ in (4), which we translate to a DFA using TranslateSpec with AP and formula given respectively in code lines 13 and 15.

\MakeFramed\FrameRestore

\endMakeFramed

Besides reach-avoid specifications it is also possible to describe many other types of specifications, such as more complex reach-avoid specification, e.g. $\phi_{PD}=\lozenge(p_{1}\land(\neg p_{2}\mathbin{\sf U}p_{3}))$, or time-bounded and unbounded safety specifications, e.g. $\phi_{BAS}=\bigwedge_{i=0}^{5}\bigcirc^{i}p_{1}$ and $\phi_{vdPol}=p_{1}\mathbin{\sf U}p_{2}$. These specifications are written in SySCoRe as

[TABLE]

Note that it is also possible to directly pass a DFA as an input to SySCoRe instead of giving the specification as an scLTL formula. SySCoRe is able to natively handle both acyclic and cyclic DFAs (see Figure 2), in contrast to many other tools [36, 32, 11, 23] that do not natively support DFAs but often rely on external tools such as PRISM [21] to compute the controller.

III-B Abstraction

SySCoRe includes two possible abstraction methods, namely finite-state abstraction for continuous-state systems in (1)-(3) and model-order reduction for continuous-state LTI systems (3). However, in order to create a finite-state abstraction of a nonlinear system (1) we require an additional approximation step before constructing a piecewise-affine finite-state abstraction. Note that the piecewise affine approximation itself is considered as an integral part of the finite-state abstraction method.

Piecewise affine approximation. To approximate a nonlinear system (1) by a PWA system (2), we partition the state space and use a standard first-order Taylor expansion to approximate the nonlinear dynamics in each partition by affine dynamics. Additionally, we compute the error introduced by this approximation. In SySCoRe, this is performed by the function PWAapproximation.

⬇

% Perform piecewise-affine approximationsPWA = PWAapproximation(sysNonLin, Np);

Here, the nonlinear system (1) is given by sysNonLin and the number of partitions in each direction is given by Np. The result is a PWA system (2) sysPWA.

Interface function. SySCoRe can construct a reduced-order abstract model $M_{r}$ and a finite-state abstract model $\hat{M}$ of the original model $M$. Let us denote the control inputs of these models respectively by $u_{r}$ and $\hat{u}$. The abstract control inputs $u_{r}$ and $\hat{u}$ need to be refined to a control input $u$ for $M$ as illustrated in Figure 3. The input refinement is performed by one or multiple interface functions, namely

[TABLE]

To refine the input $\hat{u}$ of a finite-state model to the input $u_{r}$ of a continuous-state reduced-order model, we implemented two different interface functions in the format of (6ga) and (6gb). For many cases the default interface function (6ga) should work fine, however, the option (6gb) gives more influence on the refined controller by including a feedback term. When the interface function (6gb) is used, we have to take this into account when constructing the finite-state abstraction to avoid the input bounds being violated, therefore, the interface function must be chosen before constructing the finite-state abstraction. We further use the interface function (6gc) to refine the input $u_{r}$ of a reduced-order model to the input $u$ of the full-order model. It should be noted that if only a finite-state abstraction is performed without using model-order reduction (MOR), we have $P=I,Q=0$ and $x_{t}=x_{r,t}$, hence we obtain interface functions (6ga) and (6gb) with $u_{t}=u_{r,t}$ and $x_{t}=x_{r,t}$.

Running example cont’d

It is required to select an interface function for the input refinement before starting with the temporal logic control steps. For this running example only use the default interface function (6ga) without model-order reduction, that is $u_{t}=\hat{u}_{t}$. However, if desired, the user can select the option (6gb) by setting int_f = 1 and passing this to the functions.

In the remainder of this section, we discuss how to obtain the reduced-order and finite-state abstract models.

Model-order reduction.  It is essential to include model-order reduction for high-dimensional models. For LTI systems (3) this yields a reduced-order model $M_{r}$ of the form

[TABLE]

with $x_{r}\in\mathbb{X}_{r},u\in\mathbb{U}$, $y\in\mathbb{Y}$, $w_{r}\in\mathbb{W}$, and matrices $A_{r},B_{r},B_{rw}$ and $C_{r}$ of appropriate sizes.

In SySCoRe, the function ModelReduction constructs a reduced-order model sysLTIr of dimension dimr based on the original model sysLTI by using balanced truncations on a closed loop system with a feedback matrix $F$. This feedback matrix is computed by solving discrete-time algebraic Riccati equations that can be tuned using constant f [29]. The syntax of ModelReduction is

⬇

% Construct reduced-order modelysLTIr, F] = ModelReduction(sysLTI, dimr, f)

We couple the inputs $u,u_{r}$ from $M$ (3) and $M_{r}$ (6h) using the interface function (6gc) as illustrated in Figure 3b. This is based on the theoretical results presented in [14, 39]. To compute matrices $P$ and $Q$ for the interface function, we have the function ComputeProjection that adds the matrices automatically to the object sysLTIr.

⬇

% Compute matrices P and QsLTIr = ComputeProjection(sysLTI, sysLTIr);

Finite-state abstraction. We grid the state space to construct a finite-state abstraction $\hat{M}$ of the continuous-state models (2), (3) or (6h). More specifically, we compute the abstract state space $\mathbb{\hat{X}}$ as the set consisting of the centers of the grid cells. Next, the dynamics of the abstract model is defined by using the operator $\Pi:\mathbb{X}\rightarrow\mathbb{\hat{X}}$ that maps states $x$ to the center of the grid cell it is in. Details on how to construct a finite-state abstraction of a nonlinear system or an LTI system can be found in [40, Section III], and [13, Section IV] or [39, Section IV] respectively.

In SySCoRe, the construction of the finite-state abstraction is implemented in the functions GridInputSpace and FSabstraction. The function GridInputSpace constructs the abstract input space uhat by selecting a finite number of inputs from the input space sysLTI.U.

⬇

% Construct abstract input spacehat, InputSpace] = GridInputSpace(lu, sys.U, options);

Here, lu is the number of abstract inputs in each direction and options are used to select an interface function from (6g). If interface function (6gb) or (6gc) is chosen, GridInputSpace also divides the continuous input space into a part for actuation and for feedback, and returns these spaces as output InputSpace. This is done to make sure that the input bounds $u\in\mathbb{U}$ of the original model are satisfied. Next, we use FSabstraction to compute a probability matrix that contains the transition probabilities between states for all possible inputs in uhat.

⬇

% Construct abstract modelsAbs = FSabstraction(sys, uhat, l, tol, DFA, options);

Here, sys is the continuous-state system, uhat is the abstract input space $\mathbb{\hat{U}}$, l is the number of grid cells in each direction and tol is the tolerance for truncating to zero. This means that if a probability is smaller than the value set by tol, then we set it to zero to increase sparsity and hence decrease computation time. Via efficient tensor computations, we split the computation of the probability matrix into two parts: one for the deterministic part of the transitions computed as a sparse matrix, and one for the stochastic part of the transitions. This reduces the required memory allocation and computation time drastically. For development purposes options can be used to select whether or not to use this efficient tensor computation. The complete probability matrix can then be obtained by using a tensor multiplication, however, we do not store the complete probability matrix and compute it when necessary in order to save memory.

Running example cont’d

To construct a finite-state abstraction of the car park model sysLTI (defined in code lines 1-13), we compute the abstract input space uhat:

\MakeFramed\FrameRestore

⬇

t = GridInputSpace(lu, sysLTI.U);

\endMakeFramed

and construct the abstract model sysAbs using the DFA constructed in code line 17 as follows:

\MakeFramed\FrameRestore

⬇

= 10^-6;

sysAbs = FSabstraction(sysLTI, uhat, l, tol, DFA, 'TensorComputation', true);

\endMakeFramed

III-C Similarity quantification

To quantify the similarity between the model and its abstraction (either reduced order or finite state), we compute $\epsilon$ and $\delta$ such that they satisfy the $(\epsilon,\delta)$-stochastic simulation relation as defined in [39, Definition 4]. Here, $\epsilon$ and $\delta$ represent bounds on the output and probability deviations, respectively. This simulation relation allows us to consider scLTL specifications with unbounded time properties [13].

When using model-order reduction, we construct two simulation relations, one relation $\mathscr{R}_{MOR}$ between the original model $M$ (3) and reduced-order model $M_{r}$ (6h), and one relation $\mathscr{R}$ between $M_{r}$ and the finite-state model $\hat{M}$. The simulation relations are of the form

[TABLE]

with $||x||_{D}=\sqrt{x^{\top}Dx}$ the weighted two-norm, where $D=D^{\top}\succeq 0$ is positive semi-definite. Following [39], these simulation relations can be combined into one total simulation relation between $M$ and $\hat{M}$. Following [13, Section IV.A], we can now compute the initial state of the reduced-order model as the state $x_{r,0}$ that minimizes $||x_{0}-Px_{r,0}||_{D_{r}}$, that is $x_{r,0}:=(P^{\top}D_{r}P)^{-1}P^{\top}D_{r}x_{0}$.

The computation of the simulation relation relies heavily on the coupling of the inputs $u,\hat{u}$ and disturbances $w,\hat{w}$ of the two models. The inputs are coupled through an interface function and the disturbances via a coupling kernel. This is illustrated in Figure 3 and is based on the method developed in [39]. More specifically, the underlying computation is based on finding an invariant set for the error dynamics $x_{r,t+1}-\hat{x}_{t+1}$. To this end, an optimization problem constrained by parameterized linear matrix inequalities is used to find a value for $\delta$ that corresponds with the given value of $\epsilon$ [39]. To solve this optimization problem, we use the multi-parametric toolbox (MPT3) [16] with YALMIP [26] and with either solver SeDuMi [22] or MOSEK [7].

In SySCoRe, similarity quantification is implemented in the function QuantifySim.

⬇

% Quantify similarityimRel, interface] = QuantifySim(sys, sysAbs, epsilon, options)

This function quantifies the similarity between the models sys and sysAbs, with sysAbs either a reduced-order or a finite-state approximation of sys, hence in terms of behavior $\texttt{sysAbs}\preceq\texttt{sys}$. QuantifySim yields a simulation relation simRel of the form (6i) that is stored in the object SimRel. This object includes a method to check whether two states belong to the simulation relation and a method to combine the two simulation relations from (6i) if necessary. Besides that, the function QuantifySim also returns the feedback-matrix of the interface function, when interface (6gb) or (6gc) is chosen through the options.

Running example cont’d

Next, we quantify the similarity between the model of the car stored in sysLTI and its finite-state abstraction sysAbs constructed in code line 24 by choosing a suitable value for $\epsilon$ and using the function QuantifySim.

\MakeFramed\FrameRestore

⬇

% Quantify similaritymRel = QuantifySim(sysLTI, sysAbs, epsilon);

\endMakeFramed

Piecewise affine systems. The function QuantifySim can handle both PWA (2) and LTI models (3). However, for PWA systems the probability deviation is a PWA function $\boldsymbol{\delta}(\hat{x})$ that depends on the partition of the abstract state [40].

III-D Synthesizing a robust controller

We synthesize a robust (finite-state) controller based on the dynamic programming approach described in [13], which is robust in the sense that it takes the deviation bounds $\epsilon$ and $\delta$ into account to compute a lower bound on the actual satisfaction probability. Furthermore, it is proven in [13, Theorem 4] that the resulting control policy synthesized for the abstract model can always be refined to a control policy for the actual model.

More specifically, we implicitly construct a product composition of the finite-state model $\hat{M}$ with the DFA such that computing the satisfaction probability becomes a reachability problem over this product composition. This can in turn be solved using dynamic programming by associating a robust dynamics programming operator that allows for an iterative computation of the lower bound on the satisfaction probability. Denote the state of the DFA by $q$, then the probability that a trajectory starting at $(\hat{x},q)$ reaches the set of accepting states by applying policy $\boldsymbol{\mu}$ within horizon $[1,2,\dots N]$ is denoted as $V_{N}^{\boldsymbol{\mu}}(\hat{x},q)$. This is equivalent to the probability of satisfying the specification $\phi$ over this time horizon. The probability $V$ is computed iteratively by defining the operator

[TABLE]

where $\hat{x}^{+}$and $q^{+}$ are resp. the next state of the abstract model and of the DFA, $\mathbb{E}$ is expectation with respect to the probabilistic transitions in the abstract model, $1_{Q_{f}}(q)$ is an indicator function that is equal to 1 if $q$ is inside the set of accepting states $Q_{f}$ of the DFA and is [math] otherwise, $\boldsymbol{L}:\mathbb{R}\rightarrow[0,1]$ is a truncation function, and with

[TABLE]

where $\tau_{\mathcal{A}_{\phi}}$ is the transition function of the DFA and $L(y^{+})$ is the label of the next output. This operator is robust in the sense that the probability gets reduced by $\delta$ at every time step and the worst case transition of the DFA is considered with respect to $\epsilon$. The derivation of this operator for Markov decision processes can be found in [13].

Synthesis of an abstract control strategy pol and the computation of the robust satisfaction probability satProb is performed by the function SynthesizeRobustController and it is based on the abstract model sysAbs, the specification as a DFA and the simulation relation simRel.

⬇

% Compute satisfaction probabilty and policyatProb, pol] = SynthesizeRobustController(...

sysAbs, DFA, simRel, thold, options)

We include the possibility to set the threshold thold that stops the value iteration when the difference between two iterations is smaller than this threshold. The default value is set to $1\cdot 10^{-12}$. This choice is justified by the fact that the operator in (6j) is contractive and will always converge monotonically to a fixed-point. Additionally, we include the options to compute the value function only for the initial DFA state and to compute an upper bound on the satisfaction probability. Internally, the dynamic programming algorithm computes the product between large-scale matrices (one of which is the probability matrix as mentioned in Section III-B on finite-state abstractions). By performing these computations using a tensor product [28], we gain superior computational efficiency.

The resulting control policy pol is a mapping $\mu:\mathbb{\hat{X}}\times Q\rightarrow\mathbb{\hat{U}}$ from the pair of abstract and DFA states to the abstract input space. The abstract controller can now be written as $\hat{C}:\hat{u}=\mu(\hat{x},q)$.

Running example cont’d

After specifying the desired threshold for convergence thold, we synthesize a robust control policy pol based on the finite-state abstract model sysAbs, the specification as a DFA and the simulation relation simRel constructed in code line 28. In this case, we are only interested in the satisfaction probability satProb of the initial DFA state, hence we set the options to true.

\MakeFramed\FrameRestore

⬇

% Synthesize an abstract robust controlleratProb, pol] = SynthesizeRobustController(...

sysAbs, DFA, simRel, thold, true);

\endMakeFramed

The robust satisfaction probability is computed for all $x_{0}\in\mathbb{X}$. For initial states $x_{0}=[-4,-5]^{\top}$, $x_{0}=[-8,2]^{\top}$, and $x_{0}=[4,8]^{\top}$, it equals respectively $0.60,0.52$, and $0.42$.

III-E Control refinement

To refine an abstract finite-state controller to a controller $C$ that can be implemented on the original continuous-state system (see step 5 in Figure 1) we use one or multiple interface functions from (6g) as illustrated in Figure 3. In SySCoRe, control refinement is included in the class RefineController, where it is possible to select an interface function using the options.

⬇

% Refine abstract controllerntroller = RefineController(satProb, pol, sysAbs, simRel, sys, DFA, options);

This class not only refines the finite-state input to the actual input, but also determines the state of the finite-state model based on the state of the original model.

Running example cont’d

To construct a controller $C$ that can be implemented on the original model $M$ based on the abstract control policy pol computed in code line 32, we use the following.

\MakeFramed\FrameRestore

\endMakeFramed

III-F Deployment

The final step is to deploy the controller on the model and perform simulations using ImplementController.

⬇

% Implement the controller on the modelim = ImplementController(x0, N, Controller, option);

Here, N is the desired time horizon for the simulation and option is used to supply the number of trajectories and/or additional model-order reduction inputs.

Running example cont’d

To simulate the controlled system with the Controller constructed in code line 34, we use ImplementController to obtain the state trajectory starting at $x_{0}$. Trajectories of the controlled system with three initial states are illustrated in Figure 4.

\MakeFramed\FrameRestore

⬇

40; % time horizon

% Simulate controlled system

xsim = ImplementController(x0, N, Controller);

\endMakeFramed

IV Benchmarks

To show the capabilities of SySCoRe, we included multiple benchmarks, of which some are discussed here. The package delivery has a complex specification with a cyclic DFA, the building automation system includes model-order reduction and the Van der Pol oscillator is nonlinear. We evaluate the run time and memory usage of the benchmarks, and compare SySCoRe to some existing tools.

IV-A Package delivery

With the package delivery benchmark [4], we show the capability of SySCoRe to handle complex scLTL specifications beyond basic reach-avoid scenarios, i.e., cyclic DFAs. Consider an agent traversing in a 2D space, whose dynamics can be described by an LTI system (3) with $A:=0.9I_{2}$, $B:=I_{2}$, $B_{w}:=\sqrt{0.2}I_{2}$, $C:=I_{2}$, and disturbance $w_{k}\sim\mathcal{N}(0,I_{2})$. We initialize the system using LinModel.

Define the state space $\mathbb{X}=[-6,6]^{2}$, input space $\mathbb{U}=[-1,1]$, output space $\mathbb{Y}=\mathbb{X}$, and regions $p_{1}$, $p_{2}$ and $p_{3}$ as follows: $p_{1}:=[5,6]\times[-1,1]$, $p_{2}:=[0,1]\times[-5,1]$ and $p_{3}:=[-4,-2]\times[-4,-3]$. The agent can pick up a package at $p_{1}$ and must deliver it to $p_{3}$. If the agent visits $p_{2}$ while carrying a package, it loses the package and has to pick up a new package at $p_{1}$. This corresponds to the scLTL specification $\lozenge(p_{1}\wedge(\neg p_{2}\mathbin{\sf U}p_{3}))$ implemented as in (6a). We generate the corresponding DFA using TranslateSpec.

Next, we construct a finite-state abstraction using GridInputSpace and FSabstraction. For this study, we choose a comparatively fine state abstraction $l=[400,400]$, which allows us to generate a simulation relation using QuantifySim with an epsilon of just 0.075. Note that the partition size $l$ is a tuning parameter which is determined empirically. We synthesize a robust controller for the discrete abstraction using \MakeFramed\FrameRestore

⬇

tProb, pol] = SynthesizeRobustController( ...

sysAbs, DFA, rel, thold, false);

\endMakeFramed

Since the resulting control policy is conditional on both the current system state and the DFA state, we set the 5th argument to false. By doing so, we synthesize a controller for all DFA states instead of only the initial one. The obtained robust satisfaction probability satProb over different initial states $x_{0}$ is displayed in Figure 5a and can be obtained by running \MakeFramed\FrameRestore

⬇

% Plot satisfaction probabilityotSatProb(satProb, sysAbs, 'initial', DFA);

\endMakeFramed

The peak satisfaction probability is 0.663.

Finally, we refine the controller using RefineController. To demonstrate the performance of the obtained controller, we simulate the controlled system using ImplementController for $N=60$ time steps and an initial state of $x_{0}=[-5,-5]^{T}$. Note that $N$ is an empirical parameter and should be set high enough for the DFA to terminate. As expected, the agent moves to region $p_{1}$ to pick up a package, and delivers it to $p_{3}$ whilst avoiding $p_{2}$. To plot trajectories we included the function plotTrajectories.

IV-B Van der Pol oscillator

In this benchmark, we show how SySCoRe can be applied to nonlinear stochastic systems. For this, consider the discrete-time dynamics of the Van der Pol oscillator [4], given by

[TABLE]

where the sampling time $\tau$ is set to $0.1s$, $w_{t}\sim\mathcal{N}(0,0.2I_{2})$, and $y_{t}=x_{t}$. We define the state space $\mathbb{X}=[-4,4]^{2}$, input space $\mathbb{U}=[-1,1]$, and output space $\mathbb{Y}=\mathbb{X}$. For the Van der Pol oscillator, we are looking at an unbounded safety specification (cf. (6c)), where the objective is to synthesize a controller such that the system remains in the region $p_{1}:=\mathbb{X}$ until reaching region $p_{2}:=[-1.4,-0.7]\times[-2.9,-2]$, corresponding to the scLTL specification $p_{1}\mathbin{\sf U}p_{2}$. First, we construct a DFA for the formula (6c) using TranslateSpec.

Since the dynamics of the oscillator (sysNonLin) in (6l) are nonlinear, the abstraction process is split into two parts as outlined in Section III-B. First, we construct a PWA approximation as follows: \MakeFramed\FrameRestore

⬇

% Number of grid points in each direction= [41 41];

% Perform PWA approximation

sysPWA = PWAapproximation(sysNonLin, N);

\endMakeFramed

In the second part of the abstraction step, a finite-state abstraction (sysAbs) of the PWA approximation (sysPWA) is constructed using GridInputSpace and FSAbstraction with l=[600,600] grid cells. To generate a simulation relation between this abstraction and the original model, we set $\epsilon=0.1$ and compute a suitable weighting matrix $D$ for the simulation relation on $(\hat{x},x)$, as described in Section III-C. To reduce computation time, we only use a finite number of states to compute this weighting matrix. Details on why we need this global weighting matrix can be found in [40]. \MakeFramed\FrameRestore

⬇

% Compute weighting matrix D for the simulation relation based on the following statesates = [1/8x1l, 6/10x2u; 5/7x1u, 5/17x2u; 2/13x1u, 5/9x2l; 3/4x1l, 1/7x2l; 0, 0]';

[D, $\sim$] = ComputeD(epsilon, sysPWA, sysAbs, 'interface', int_f, 'states', States);

% Quantify similarity

[rel, sysPWA] = QuantifySim(sysPWA, sysAbs, epsilon, 'interface', int_f, 'weighting', D);

\endMakeFramed

Note that QuantifySim returns sysPWA instead of the usual interface, because each piecewise-affine system gets its own interface function and we store this directly in sysPWA.

Next, we use SynthesizeRobustController to synthesize a robust controller for sysAbs and show the satisfaction probability (displayed in Figure 5b) using plotSatProb. Finally, we refine the controller as follows: \MakeFramed\FrameRestore

⬇

troller = RefineController(satProb, pol, sysAbs, rel, sysPWA, DFA, int_f);

\endMakeFramed

As before, ImplementController is used to simulate the system.

IV-C Building automation system

In the last benchmark, we address a large-scale system showcasing the model-order reduction capabilities of SySCoRe. We consider a 7D affine stochastic system of a building automation system, regulating the temperature in two zones influenced by a 6D disturbance. A detailed description including the system dynamics can be found in [3, 10]. The goal is to synthesize a controller maintaining the temperature in zone one at $20^{\circ}C$ with a maximum permissible deviation of $\pm 0.5^{\circ}C$ for 6 consecutive time steps. We translate the specification (6b) to a DFA using TranslateSpec.

The dynamics of this building automation system are not of the form (3), since it is influenced by a Gaussian disturbance with mean $\mu\neq 0$ and variance $\Sigma\neq I$. Furthermore, it is not an LTI system, but affine, which cannot be handled by our current implementation of model-order reduction. To deal with the disturbance, we first transform the system to a system with Gaussian disturbance $w\sim\mathcal{N}(0,I)$ using the following: \MakeFramed\FrameRestore

⬇

% Transform the modelysLTI, a] = NormalizeDisturbance(sysLTI,a);

\endMakeFramed

To deal with the affine dynamics, we perform a steady-state shift and simulate the steady-state system that has LTI dynamics. After performing the control synthesis steps, we compensate for this steady-state shift again to obtain the dynamics of the actual system.

Now, we can start with the synthesis steps. First, we reduce the 7D model to a 2D reduced-order model (see Eq. (6h)) using function ModelReduction with f = 0.098 and dimr=2. \MakeFramed\FrameRestore

⬇

% Perform model-order reductionysLTIr, $\sim$] = ModelReduction(sysLTI, dimr, f);

\endMakeFramed

As mentioned in Section III-B, we use an interface function of the form (6gc), which is selected using int_f = 1 and compute the matrices $P$ and $Q$ using ComputeProjection. Next, we define the state and input spaces, and the output regions and APs for the reduced-order model as before.

To construct the finite-state abstraction of the reduced-order model, we first grid the input space with lu = 3. \MakeFramed\FrameRestore

⬇

% Construct abstract input spacehat,sysLTIr.U] = GridInputSpace(lu, sysLTIr.U, 'interface', int_f, 0.6, 0.175);

\endMakeFramed

Here, we have chosen to use $60\%$ of the input space for actuation and $17.5\%$ for feedback. This leaves $22.5\%$ for the $Qx_{r,t}$ part of the interface function, which is currently not guaranteed to be satisfied.

Before constructing a finite-state abstraction of the reduced-order model, we reduce the state space to increase the computational speed. This step is currently only available for invariance specifications and is performed by ReduceX, which performs a number of backwards iterations on the safety region $P_{1}$ to determine a good guess of the invariant set. This set is then used as the reduced state space. The construction of the finite-state abstraction of the reduced-order model is as before, except that we give the total number of grid cells as input l, instead of the number of grid cells in each direction. \MakeFramed\FrameRestore

⬇

% Reduce the state space to speed up computationsysLTIr, $\sim$] = ReduceX(sysLTIr, sysLTIr.U{2}, P1, 'invariance', 5);

% Construct finite-state abstraction

l = [3000*3000]; % Total number of grid cells

tol=10^-6;

sysAbs = FSabstraction(sysLTIr, uhat, l, tol, DFA, 'TensorComputation', true);

\endMakeFramed

To relate the reduced-order finite-state model sysAbs to the original model sysLTI, we construct two simulation relations: relation rel_1 with $\epsilon_{1}=0.2413$ between sysLTI and sysLTIr, and relation rel_2 with $\epsilon_{2}=0.1087$ between sysLTIr and sysAbs: \MakeFramed\FrameRestore

⬇

% Compute MOR simulation relationel_1, K, kernel] = QuantifySim(sysLTI, sysLTIr, epsilon_1, 'MOR', sysAbs);

% Compute finite-state simulation relation

[rel_2] = QuantifySim(sysLTIr, sysAbs, epsilon_2);

% Combine simulation relations

rel = CombineSimRel(rel_1, rel_2, sysLTIr, sysAbs);

\endMakeFramed

For model-order reduction we have to explicitly define the coupling kernel matrix $F$, that is later used to compute the disturbance of the reduced-order model as $w_{r}\!=\!w+F(x\!-\!Px_{r})$. For details see [39].

Synthesizing and refining the controller are done as before and the satisfaction probability of the reduced-order model is shown in Figure 5c (obtained through plotSatProb). We simulate the controlled system $N_{s}=6$ times, making sure the output is shifted with respect to the steady-state solution. \MakeFramed\FrameRestore

⬇

% Simulate controlled system Ns timess = 6;

xsim = ImplementController(x0, N, Controller, Ns, 'MOR', sysLTIr, kernel);

\endMakeFramed

The resulting trajectories can be evaluated using plotTrajectories.

IV-D Performance evaluation

The performance of SySCoRe is evaluated on the benchmarks mentioned above. The details of the benchmarks and their total run time and memory usage are reported in Table II. The computation times per step are reported in Table III. The data has been obtained on a computer with a 2,3 GHz Quad-Core Intel Core i5 processor and 16 GB 2133 MHz memory by taking the average over $5$ computations. Here, we observed a maximum $6\%$ standard deviation.

Table II can be used to compare the different benchmarks with respect to the computations performed by SySCoRe. The main difference between the running example and the package delivery benchmark is the DFA. The DFA of the package delivery benchmark requires more memory, however, the increase in computation time is small. Due to the simple DFA of the running example, we only compute the satisfaction probability for the initial DFA state. This will not suffice for the package delivery benchmark, which is the reason that more computation time is spent on steps (4)-(6) compared to the running example (see Table III). The computation time for the nonlinear benchmark is large, however, the memory usage remains reasonable. The increase in computation time is mainly due to the fine gridding. We can also see in Table III that the similarity quantification takes a considerable amount of time. This is because we perform this step for each partition separately (1600 times in this case). For higher-dimensional systems that require model-order reduction (building automation system benchmark), the computation time and memory usage increase substantially, mainly due to the the fact that the similarity quantification has to be performed multiple times. However, we also see from Table III a large increase in the computation time for the controller synthesis.

Table III shows that the similarity quantification of step (3) requires the most computation time, followed by the finite-state abstraction of step (2). The large computation time of the similarity quantification is due to solving an optimization problem constrained by parameterized matrix inequalities that could be non-convex. For most abstraction-based approaches in the literature, the main bottleneck is the finite-state abstraction. This shows the efficiency of our tensor-based implementations. It should also be noted that the efficiency of the tensor computations is also exploited in the control synthesis step.

IV-E Comparison to existing tools

A comparison of the results on the benchmarks obtained by SySCoRe and current tools is given in Table IV. The package delivery benchmark has a complex DFA and cannot be handled natively by tools AMYTISS, FAUST, and StocHy (see Table IVa). SReachTools can only handle safety specifications and is not applicable to this benchmark. The Van der Pol oscillator benchmark poses significant challenges for the tools due to its nonlinear dynamics, as reported in Table IVb. Only AMYTISS can solve a benchmark that resembles this one as considered in [1] with multiplicative noise instead of additive noise. AMYTISS can only handle systems with a bounded disturbance, hence it cannot directly solve the benchmark as presented here.

The benchmark on the building automation system can be solved by AMYTISS, SReachTools, and StocHy without being able to use model-order reduction. This benchmark considers a stochastic safety problem and the performance of multiple tools is compared in [2, 1]. Table IVc reports the results of SySCoRe together with the results from running the repeatability packages of [2, 1] on a computer with a 2,3 GHz Quad-Core Intel Core i5 processor and 16 GB 2133 MHz memory. For StocHy, there was no repeatability package available, however, since the computational power of the CPU used for the results in [1] was more than our computer, we included the results of [1] as a lower bound on the computation time required by StocHy. Note that this benchmark belongs to the class of partially degenerate systems [35]. The formulation of the abstraction error for this class is available but the current version of FAUST does not natively support partially degenerate systems. With respect to the computation time, SReachTools performs best, and AMYTISS and StocHy require a longer computation time. Though from the results in [1], we see that AMYTISS could be faster than the current implementation of SySCoRe when parallel execution within CPUs is available (this parallel computation will be exploited in future versions of SySCoRe). With respect to accuracy, both AMYTISS and StocHy obtain a maximum reachability probability smaller than SySCoRe, while SReachTools still outperforms SySCoRe. This shows that SReachTools is the best option for this benchmark, which is expected since it is developed exactly for linear systems and stochastic reach-avoid problems with small disturbances.

V Summary and extensions

This paper described the first release of SySCoRe, a tool that excels at control synthesis problems for systems with a large (unbounded) stochastic disturbances and temporal specifications with possibly unbounded time horizon. It combines reduced-order models and finite abstractions with formal guarantees obtained by coupled stochastic simulation relations. SySCoRe substantially extends the class of models and temporal specifications that current tools can handle for control synthesis. Furthermore, the modular development of SySCoRe allows ease of use and facilitates future extensions. The efficient implementation of tensor computations in SySCoRe allows for fast computations, which can be exploited further by including more parallel computations as done in AMYTISS.

An important direction for future releases is to extend the current implementation of model-order reduction to piecewise-affine systems, such that it can also be applied to nonlinear systems. Currently, only Gaussian disturbances are implemented in SySCoRe, however, extensions to other distributions are under way and require deriving new inequality constraints for the optimization problem solved in the similarity quantification. The computation time of the similarity quantification is large due to solving optimization problems constrained by parameterized matrix inequalities that could be non-convex. We are working on improving the efficiency of solving this optimization.

The modular implementation of SySCoRe can be utilized to integrate model-order reduction with discretization-free approaches such as the kernel method of SReachTools [41, 38] and the barrier certificates [18], or to perform synthesis for stochastic systems with parametric uncertainty [33]. To get non-trivial lower bounds, SySCoRe currently requires fine-tuning the hyper parameters (e.g., the grid size and the output deviation). It is of interest to automatically design these parameters or to provide guidelines to the user on the appropriate values depending on the case study.