Alberto Casagrande
email: acasagrande@units.it

Thao Dang
email: thao.dang@univ-grenoble-alpes.fr

Luca Dorigo
email: dorigo.luca001@spes.uniud.it

Tommaso Dreossi
email: tommasodreossi@gmail.com

Carla Piazza
email: carla.piazza@uniud.it

Eleonora Pippia
email: pippia.eleonora@gmail.com

Parameter Synthesis of Polynomial Dynamical Systems

Keywords: Parametric Dynamical Systems, Bernstein Coefficients, Signal Temporal Logic, Neural Network Verification

Parametric dynamical systems emerge as a natural formalism for modeling natural and engineered systems ranging from biology, epidemiology, and medicine to cyber-physics. Parameter tuning is a complex task which, in many cases, is performed exploiting heavy simulations that have high computational complexity and do not ensure the correctness of the synthesized systems. In this manuscript, we consider the problem of parameter synthesis for discrete-time polynomial systems. We propose a formal method based on Bernstein coefficients that allows refining the set of parameters according to a temporal specification defined as a Signal Temporal Logic formula. We prove that the synthesized system is correct with respect to the specification and we demonstrate the scalability of the approach by implementing it in the C++ library Sapo, also available within a stand-alone application and as a web application. Finally, we illustrate the tool usage and the interface through a simple yet realistic epidemiological model and consider an intriguing application enhancing the accuracy of the verification of neural networks.

Introduction

Dynamical systems are standard mathematical models used to describe the temporal evolution of systems (e.g., see [START_REF] Brin | Introduction to Dynamical Systems[END_REF][START_REF] Galor | Discrete dynamical systems[END_REF]). Parameters allow specifying and analysing the behaviour of a possibly infinite set of models at one time (e.g., see [START_REF] Paul | Modeling, identification and simulation of dynamical systems[END_REF]). In the case of partial knowledge over the observed system parameters, they allow drawing general conclusions despite the lack of information. They are also helpful in abstracting stochastic noise when this has no dramatic effects on global behaviour. During the design phases of engineered systems, parameters have to be refined to guarantee the desired result. In the context of parametric dynamical systems, we are interested in defining algorithmic techniques for solving both reachability analysis and parameter synthesis. Reachability analysis aims to determine the values assumed by the system variables during an evolution that starts from an initial set of states. Instead, parameter synthesis refines the set of parameters to ensure that the reachable states satisfy a specification.

A dynamical system is either discrete-time or continuous-time depending on whether the changes occur in discrete time instants or in continuous time, respectively. The typical formalisms used to describe discrete-time and continuoustime dynamical systems are difference and differential equations, respectively. In this work, we consider discrete-time dynamical systems since, in most contexts, the system under analysis can be observed only at fixed time instants. In some cases, discrete numerical integration methods (e.g., Euler, Runge-Kutta) are applied in the study of continuous-time models, see, e.g. [START_REF] Althoff | Combining zonotopes and support functions for efficient reachability analysis of linear systems[END_REF][START_REF] Bak | Simulation-equivalent reachability of large linear systems with inputs[END_REF][START_REF] Althoff | Set propagation techniques for reachability analysis[END_REF] where discretization is used to compute discrete step dynamics and, then, the continuoustime flowpipe is over-approximated.

We focus, in particular, on parametric discrete-time dynamical systems in which the dynamics are polynomials. Many works exist on the analysis of linear systems, whose properties, such as convexity preservation, allow to define efficient algorithms, e.g., for computing sets of reachable points. Some attempts for tackling the reachability problem in the case of nonlinearity used approximation and symbolic computation techniques. A recent and detailed review of the state-of-the-art in the context of reachability analysis can be found in [START_REF] Althoff | Set propagation techniques for reachability analysis[END_REF].

Polynomial systems allow us to exploit the properties of Bernstein coefficients [START_REF] Rida | The Bernstein polynomial basis: A centennial retrospective[END_REF] to over-approximate the reachability problem. By imposing linear dependence on the parameters, we can rely on linear optimization algorithms to efficiently solve our problem.

As far as parameter synthesis is concerned, we use temporal logic for specifying the properties of interest. Temporal logics were introduced in formal verification of hardware and software systems for defining specifications involving runtime acceptable computations [START_REF] Pnueli | The temporal logic of programs[END_REF]. Their success in that field is not only due to their balance between expressive power and ease-of-use, but also to the development of efficient algorithms for checking whether a system satisfies a temporal logic specification (e.g., see [START_REF] Clarke | Model checking[END_REF]).

Temporal logic has found applications outside formal verification, for instance, in monitoring. Within this context, the system is treated as a black box whose observable behaviours can be analysed by evaluating the satisfaction of the desired temporal specification. Signal Temporal Logic (STL [START_REF] Maler | Monitoring temporal properties of continuous signals[END_REF][START_REF] Maler | Checking temporal properties of discrete, timed and continuous behaviors[END_REF]) has been introduced for specifying properties of dense-time real-valued signals. It is particularly suitable for monitoring both industrial case studies (e.g., see [START_REF] Jones | Analog property checkers: a DDR2 case study[END_REF][START_REF] Jin | Mining requirements from closed-loop control models[END_REF][START_REF] Jyotirmoy | Robust online monitoring of signal temporal logic[END_REF]) and biological systems (e.g., see [START_REF] Donzé | Robustness analysis and behavior discrimination in enzymatic reaction networks[END_REF][START_REF] Stoma | STL-based analysis of TRAIL-induced apoptosis challenges the notion of type I/type II cell line classification[END_REF][START_REF] Dang | Set-based analysis for biological modeling[END_REF]). It has also been used in the study of parametric systems (e.g., see [START_REF] Asarin | Parametric identification of temporal properties[END_REF][START_REF] Hoxha | Mining parametric temporal logic properties in model-based design for cyber-physical systems[END_REF]) where parametric disturbance rejection properties are formalized in STL and then verified. One of its exciting aspects is its semantics: in addition to the classical semantics, where the result of the evaluation of a formula is a truth value, STL offers a quantitative semantics that gives the idea of "how robustly" a property is satisfied [START_REF] Fainekos | Robustness of temporal logic specifications for continuous-time signals[END_REF][START_REF] Bartocci | On the robustness of temporal properties for stochastic models[END_REF][START_REF] Haghighi | Control from signal temporal logic specifications with smooth cumulative quantitative semantics[END_REF].

In our setting, we exploit Bernstein coefficients' properties again and we automatically infer new linear constraints on the parameters. Such constraints ensure that the specification is satisfied.

Recently several methods and tools have been developed for reachability analysis and parameter synthesis. Some examples are HyTech [START_REF] Henzinger | HyTech: A model checker for hybrid systems[END_REF], d/dt [START_REF] Asarin | The d/dt tool for verification of hybrid systems[END_REF], Breach [START_REF] Donzé | Breach, a toolbox for verification and parameter synthesis of hybrid systems[END_REF], SpaceEx [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF], Ariadne [START_REF] Benvenuti | Ariadne: Dominance checking of nonlinear hybrid automata using reachability analysis[END_REF], Flow* [START_REF] Chen | Flow*: An analyzer for non-linear hybrid systems[END_REF], pyHybridAnalysis [START_REF] Casagrande | pyHybridAnalysis: A Package for Semantics Analysis of Hybrid Systems[END_REF][START_REF] Casagrande | ϵ-semantics computations on biological systems[END_REF], dReach [START_REF] Kong | dReach: δ-Reachability Analysis for Hybrid Systems[END_REF], CORA [START_REF] Althoff | Implementation of interval arithmetic in CORA 2016[END_REF], DynIbex [START_REF] Mullier | Validated computation of the local truncation error of Runge-Kutta methods with automatic differentiation[END_REF], JuliaReach [START_REF] Bogomolov | Juliareach: a toolbox for set-based reachability[END_REF], Kaa [START_REF] Kim | Kaa: A Python implementation of reachable set computation using bernstein polynomials[END_REF], and HyPro [START_REF] Schupp | HyPro: A C++ library of state set representations for hybrid systems reachability analysis[END_REF].

Similarly to what we present here, Breach performs parameter synthesis. However, it exploits linearizations to compute the flow over a finite set of points. On the other hand, Flow* and HyPro deal with non-linear systems by exploiting Taylor model arithmetic techniques. They also handle hybrid systems, but not parametric ones. More detailed comparison can be found in [START_REF] Nuzzo | A platform-based design methodology with contracts and related tools for the design of cyber-physical systems[END_REF] where a methodology that integrates different tools for the design of cyber-physical systems is described. ARCH-COMP (International Workshop on Applied Verification of Continuous and Hybrid Systems -Competition) reports the state of the art in the field, see, e.g., [START_REF] Geretti | ARCH-COMP20 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics[END_REF][START_REF] Geretti | ARCH-COMP21 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics[END_REF]. This manuscript extends [START_REF] Dang | Parameter synthesis through temporal logic specifications[END_REF] and [START_REF] Dreossi | Sapo: Reachability computation and parameter synthesis of polynomial dynamical systems[END_REF]. On the one hand, we present the algorithms described in [START_REF] Dang | Parameter synthesis through temporal logic specifications[END_REF] trying to keep the technicalities as simple as possible and give more intuitions through simple examples. On the other hand, we introduce some improvements to obtain more precise results. As far as the implementation is concerned, we extended the C++ library Sapo introducing the input language SIL and a web application, significantly improving usability. Finally, we briefly show the usefulness of Sapo in the contexts of epidemiological models and neural networks analysis.

The work is organized as follows. Section 1 formally introduces our setting and the problems we consider. Section 2 gives more details on the dynamics and on the set of points/parameters we manipulate. Such details are necessary for defining our reachability algorithm. Section 3 focus on the parameter synthesis problem. Sapo library, stand-alone tool, and the web application are described in Section 4. Section 5 presents applications in two different realms, while Section 6 concludes the work and presents some possible future developments.

Preliminaries

In this section, we describe the formal notions at the basis of our work. We are interested in dynamical systems, i.e., phenomena evolving in time, and we focus on models of such systems in which time evolves through discrete steps. Specifically, we consider polynomial discrete-time dynamical systems whose parameters are used to tune the model to match given observations. We also present a formal language for specifying the expected evolutions of the system called Signal Time Temporal Logic (STL) [START_REF] Maler | Monitoring temporal properties of continuous signals[END_REF] that was originally introduced in the context of monitoring.

Parametric Dynamical Systems

Let us start with basic notions [START_REF] Kalman | Topics in mathematical system theory[END_REF] necessary to define dynamical systems. The state of a system is a description that is sufficient to predict the system's future. In this work, we deal with memoryless systems, i.e., systems for which at time t it is possible to predict the future states without recurring to states prior to t. The space of all the possible system states is called the state space of the dynamical system.

The nature of a dynamical system is related to the structure of time it relies on. If the time of a system ranges on non-negative real values, then the system is called continuous, while if the time is described by naturals, then the system is said to be discrete. The evolution of the system over time is a continuous trajectory (in the continuous-time case) or a sequence (in the discrete-time case) of states through the state space.

The rules that allow us to determine the state of the systems are called dynamics or laws of evolution. Typically the dynamics of dynamical systems are described by differential equations or difference equations depending on whether they are continuous-time or discrete-time, respectively. Finally, the initial condition is the state at an initial time from which an evolution starts.

Parametric discrete-time dynamical systems are dynamical systems that evolve in discrete time and include parameters in their definition.

Definition 1 (Parametric Discrete-Time Dynamical System). A parametric discrete-time dynamical system is a tuple D = (X , P, f ) where:

• X ⊆ R n is the state space;

• P ⊆ R m is the parameter space;

• f : X × P → X is a function.

The evolutions of parametric discrete-time dynamical systems are governed by difference equation of the form:

x t+1 = f (x t , p) (1)
where x ∈ X is the state of the system, p ∈ P are the parameters, and t ∈ N is a discrete-time variable. During the evolution of the dynamical system the state variables can change their values, while the parameters remain constant.

Definition 2 (Trajectory of Parametric Discrete-Time Dynamical System). A trajectory of a parametric discrete-time dynamical system D = (X , P, f ) starting from an initial state x ∈ X with parameters p ∈ P is a function ξ p x : N → X such that ξ p x is the solution of the difference equation x t+1 = f (x t , p) with initial condition x, that is:

ξ p x (0) = x and ∀t ∈ N >0 , ξ p x (t + 1) = f (ξ p x (t), p). (2)
Let D = (X , P, f ) be a parametric discrete-time dynamical system, X 0 ⊆ X be a set of initial conditions, and P ⊆ P be a set of parameters. The set of all possible trajectories of D with initial conditions in X 0 and parameters in P is defined as: Ξ(X 0 , P ) = {ξ p x0 | x 0 ∈ X 0 , p ∈ P, and ξ p x0 is a trajectory of D}.

Given two states x, x ′ ∈ X , we say that x ′ is reachable from x in time 0 ≤ t < +∞ if there are a parameter p ∈ P and a trajectory ξ p x of S starting in x such that x ′ = ξ p x (t). The set of all the states reached by the system from x 0 ∈ X with parameter p ∈ P is defined as:

Reach p (x 0 ) = {x ′ | x ′ = ξ p x0 (t), t ∈ N}. (4)
We can extend the notion of reachability to sets, that is, given a set of initial conditions X 0 ⊆ X and a parameter set P ⊆ P, the reachable set is the set of all the states reachable by the system:

Reach P (X 0 ) = x0∈X p∈P Reach p (x 0 ). ( 5
)
The definition of reachable set reflects the behaviour of the considered dynamical system for an infinite amount of time. However, we might be interested in studying a model for a bounded time horizon. In particular, let I = [a, b] be a closed bounded interval of N. The set of states reachable within the time interval I is defined as:

Reach p I (x 0 ) = {x ′ | x ′ = ξ p x0 (t), t ∈ I} (6)
Reach P I (X 0 ) = x0∈X0 p∈P Reach p I (x 0 ). ( 7
)
Example 1. Let us introduce a simple yet significant example that will be used throughout the paper to help the reader. SIR is a classic 3D epidemic model [START_REF] Kermack | A contribution to the mathematical theory of epidemics[END_REF] that describes the spread of an epidemic disease in the simplest possible form.

In the SIR model the population is partitioned into 3 classes: S susceptible, I infected, and R removed. Two parameters characterize the disease: β represents the transmission rate, i.e., the probability for a susceptible to become infected once there is a contact with a sick person; α models the recovery rate, i.e., the inverse of the time required for moving from infected to removed. The susceptible individuals at time t + 1 are the susceptible individuals at time t decremented by the individuals that at time t + 1 become infected, i.e., decremented by β * S t * I t . Similarly, the infected at time t + 1 gain the new infected and lose the removed, i.e., α * I t . Finally, the removed are incremented with the new removed. Every removed individual is considered immune to the disease in this simple model. In this work, we will refer to a normalized population of size 1, i.e., the system's variables will range from 0 to 1.

With our notation an SIR model is defined as D = (X , P, f ) where:

• X = [0, 1] 3 is the state space described by the variables S, I, and R;

• P = [0, 1] 2 is the parameter space described by the parameters α and β;

• f : X × P → X is defined as

S t+1 = S t -βS t I t I t+1 = I t + βS t I t -αI t R t+1 = R t + αI t
Let us consider the case of β = 0.5 and α = 1/5. If the initial point x 0 is (S 0 , I 0 , R 0 ) where S 0 = 0.99, I 0 = 0.01, and R 0 = 0, we have: (x 0 ), where x 0 = (0.99, 0.01, 0), described in Example 1. The points are labelled as x i instead of ξ p x 0 (i) to simplify the plot.

Reach (0,5,1/5) [0,6]
As soon as one of the parameters or one of the variables is estimated through an interval, we obtain an infinite set of points.

Synthesis Language: Signal Temporal Logic

Temporal logics [START_REF] Pnueli | The temporal logic of programs[END_REF] are well-known formalisms for specifying desirable properties of systems under investigation. In formal verification, temporal logic properties are used to identify system's bugs. Model-checking algorithms take as input a model and a temporal formula and return a counter-example if the model does not satisfy it, i.e., a proof of the existence of unwanted behaviours. In our setting, we use temporal logic formulas to synthesize a set of parameters that guarantees the temporal property.

Example 2. Let us consider the SIR model described in Example 1 again and imagine that we do not know the transmission rate β. We may try to estimate it by using some statistics. For instance, if we know that the number of infected ranged is in the interval [0.01, 0.02] at time 0, while it roses to [0.015, 0.03] after 6 days, and if we also know that the recovery time is 4 days, i.e., α = 1/4, then we can intuitively deduce that the transmission rate cannot be too close to 1.

Signal Temporal Logic (STL [START_REF] Maler | Monitoring temporal properties of continuous signals[END_REF][START_REF] Maler | Checking temporal properties of discrete, timed and continuous behaviors[END_REF][START_REF] Donzé | Efficient robust monitoring for STL[END_REF]) formalizes properties of dense-time real-valued signals, i.e., functions defined on dense intervals. In our context, the trajectories of a dynamical system are the signals over which STL formulas are evaluated. We specifically focus on STL formulas in positive normal form over a bounded time horizon.

Let Σ = {σ 1 , σ 2 , . . . , σ k } be a finite set of predicates mapping R n into Boolean values. For a given j ∈ {1, 2, . . . , k}, the predicate σ j is of the form σ j = s j (x 1 , x 2 , . . . , x n ) ≤ 0 where s j : R n → R is a function over the variables x 1 , x 2 , . . . , x n .

A Signal Temporal Logic formula is generated by the following grammar:

φ := σ | ¬φ | φ ∧ φ | φ ∨ φ | φU I φ (8)
where σ ∈ Σ is a predicate and I = [a, b] is a bounded closed interval of N. For t ∈ N, the shifted interval t + I is the set {t + t ′ | t ′ ∈ I}. A STL formula is in positive normal form when does not use the negation operator ¬.

There are two elements that distinguish STL from other logics:

• the predicates σ are evaluated on real-values, that in our case are the states of the dynamical system;

• the temporal operators φU I φ are decorated with intervals that determine the temporal windows on which the operators are defined.

Other standard temporal operators can be derived as follows:

• true ⊤ ≡ 0 ≤ 0;

• false ⊥ ≡ 1 ≤ 0;

• ¬⊥ ≡ ⊤ and ¬⊤ ≡ ⊥;

• release φ 1 R [a,b] φ 2 ≡ (φ 2 U [a,b] (φ 1 ∧ φ 2 )) ∨ (φ 2 U [b+1,b+1] ⊤);
• eventually in the future F I φ ≡ ⊤U I φ;

• always in the future

G I φ ≡ ⊤U [a,a] (⊥R [0,b-a] φ).
Any STL formula ϕ can be rewritten as an equivalent STL formula in positive normal form φ such that φ ≡ ϕ (e.g., see [START_REF] Sadraddini | Robust temporal logic model predictive control[END_REF]). Because of this result, in the remaining part of this manuscript, we exclusively deal with STL formulas in positive normal form.

Example 3. The property informally described in Example 2 can be stated as: x 0 (i) to simplify the plot. The formula does not hold because ξ p x 0 (6) does not satisfy the formula 0.015 ≤ I ≤ 0.03.

F [6,6] (0.015 ≤ I ≤ 0.03). x 0 x 1 x 2 x 3 x 4 x 5 x 6
An interesting aspect of STL is its semantics. Two semantics for STL formulas can be defined: qualitative semantics, also known as Boolean semantics, and quantitative semantics. Intuitively, the former establishes the truth value of a formula over a trace, telling us whether a formula holds or not; the latter provides additional information on how robustly a trace satisfies (or not) a formula. This article refers to the qualitative semantics of STL which is defined over trajectories as follows.

Definition 3 (Qualitative Semantics [START_REF] Maler | Monitoring temporal properties of continuous signals[END_REF]). Let D = (X , P, f ) be a parametric discrete-time dynamical system. Let ξ p x be a trajectory of D starting from x ∈ X with parameter p ∈ P. Let t ∈ R ≥0 be a time instant, and φ be an STL formula. The qualitative semantics of φ at time t over ξ p

x is given by the following inductive definition: In order to simplify the plot, the points are labelled as x i meaning ξ p x 0 (i).

ξ p x , t |= σ iff σ(ξ p x (t)) is true ξ p x , t |= φ 1 ∧ φ 2 iff ξ p x , t
(I ≤ 0.02)U [0,3] (R ≥ 0.01), since I ≤ 0.02 becomes false in ξ p x0 (3), while R ≥ 0.01 is not yet true.

Parametric Reachability and Synthesis Problems

Within the above described context, i.e., models defined in terms of parametric discrete-time dynamical systems and properties specified as STL formulas, we are ready to formalize two problems of interest: reachability computation and parameter synthesis. Definition 4 ((Bounded) Parametric Reachability Computation). Let D = (X , P, f ) be a dynamical system, X 0 ⊆ X be a set of initial conditions, P ⊆ P be a set of parameters, and I be a closed interval of N .

The Parametric Reachability Computation problem asks to compute the set of points that can be reached with a trajectory of D starting from x ∈ X 0 and p ∈ P , i.e., it asks to compute the set Reach P (X 0 ).

The Bounded Parametric Reachability Computation problem asks to compute the set of points that can be reached with a trajectory of D starting from x ∈ X 0 and p ∈ P within the time interval I, i.e., it asks to compute the set Reach P I (X 0 ).

As one can imagine in the general case both reachability and bounded reachability are not computable (see, e.g., [START_REF] Henzinger | What's decidable about hybrid automata?[END_REF][START_REF] Brihaye | On reachability for hybrid automata over bounded time[END_REF]). However, in Section 2 we will discuss how bounded reachability can be over-approximated in the case of polynomial systems.

Definition 5 (Parameter Synthesis). Let D = (X , P, f ) be a dynamical system, X 0 ⊆ X be a set of initial conditions, P ⊆ P be a set of parameters, and φ be an STL specification. The Parameter Synthesis Problem asks to compute the largest subset P φ ⊆ P such that ∀x 0 ∈ X 0 , ∀p ∈ P φ , ξ p x0 |= φ, where ξ p x0 is a trajectory of D. Now that we have the definition of the problem it is natural to observe that we had to refer the notion of satisfiability to time t = 1, since at time t = 0 the set X 0 is given. Hence, it makes no sense to consider parameter refinement at time t = 0, since they have not been used yet.

Again, while the problem is in general not computable, in Section 3 we will present a method for under-approximating the solution over polynomial systems.

Parametric Reachability Computation

Set Based Reachability

In the context of formal verification, the interest is not only in simulating one specific trajectory of a dynamical system starting from a given initial condition, but also in computing a flowpipe that contains all the states reachable by the system from a possibly infinite set of initial conditions.

A standard methodology for computing flowpipes in the case of parametric discrete-time dynamical systems consists of the following steps:

• Choice of a domain DX providing finite representations for the set of points to be manipulated. The domain should allow any set of points to be both over and under-approximated. The most used domains are ellipsoids, hyper-rectangles, zonotopes.

• Choice of a domain DP providing finite representations for the set of parameters to be considered.

• Choice of a family F of admissible functions to be used in the dynamical systems. For instance, linear systems or polynomial systems.

• Definition of an algorithm reachStep that for each representable set of points X ∈ DX, each representable set of parameters P ∈ DP , and each admissible function f ∈ F , approximates the set of points f (X, P ) = {f (x, p) | x ∈ X and p ∈ P } with a representable set of points X ′ ∈ DP .

Notice that in the general case the set of points f (X, P ) could be not exactly representable in the chosen domain of sets of points DX. For instance, if one considers hyper-rectangles and quadratic functions, the image of a hyper-rectangle is not a hyper-rectangle. Hence, the algorithm has to approximate such image with a new hyper-rectangle. Given a set X 0 ∈ DX, a set of parameters P ∈ DP , a function f ∈ F , and a time interval I = [a, b], Algorithm 1 can be used to compute an approximated solution of the bounded reachability problem, i.e., to approximate Reach P I (X 0 ), by iterating reachStep computations. In the next section, we will focus on our choice of the domains DX and DP and on the family F of polynomial functions and describe our reachStep algorithm based on Bernstein coefficients.

Algorithm 1 Reachability

Sets of Points and Sets of Parameters: Polytopes

The domain that we consider for sets of points is that of polytopes.

Definition 6. A polytope Q ⊂ R n is a closed, compact, bounded subset of R n such that there is a finite set of half-spaces H = {h 1 , . . . , h m } whose intersection is Q, that is: Q = m i=1 h i ( 9
)
where a half-space h = {x | dx T ≤ c} is a set characterized by a non-null normal vector d ∈ R n and an offset c ∈ R.

The linear constraints that generate the half-spaces can be organized in a matrix D ∈ R m×n called template and a vector c ∈ R m called offset vector, in short offsets. The polytope generated by the template D and the offset vector c is denoted by ⟨D, c⟩. Notice that not all the pairs ⟨D, c⟩ define a nonempty polytope.

Example 5. Let us consider the SIR model of Example 1 again. Imagine that we start observing the system after some days from the initial spread of the disease. It can be the case that we do not precisely know the number of susceptible, infected, and removed, but we only have some estimations saying that the number of removed is at most the 2% of the population. This means that we can consider the following polytope:

           S ≥ 0 I ≥ 0 R ≥ 0 R ≤ 0.02 S + I ≥ 0.98
The corresponding template matrix and offset vector are

D =       -1 0 0 0 -1 0 0 0 -1 0 0 1 1 1 0       c =       0 0 0 0.02 0.98      
Our reachability algorithm fixes a template matrix and computes at each step a new offset. Boxes and parallelotopes are particular classes of polytopes and are helpful in our context to provide canonical representations. Boxes are n-dimensional generalizations of rectangles.

Definition 7 (Box).

A set B ⊂ R n is a box if and only if it can be expressed as the Cartesian product of n intervals, that is:

B = [x 1 , x 1 ] × • • • × [x n , x n ] = n i=1 [x i , x i ], ( 10
)
where x i , x i ∈ R, for i ∈ {1, . . . , n}. It is easy to see that a box B = [x 1 , x 1 ] × • • • × [x n ,
x n ] is a polytope since it can be represented by a template D ∈ R 2n×n and offsets c ∈ R 2n where:

D =           1 . . . 0 . . . . . . . . . 0 . . . 1 -1 . . . 0 . . . . . . . . . 0 . . . -1           c =           x 1 . . . x n -x 1 . . . -x n           (11)
The unit box is the box

[0, 1] × [0, 1] × • • • × [0, 1].
A parallelotope is a centrally symmetric convex polytope whose opposite facets are parallel. It can be represented as a collection of linear constraints. Definition 8 (Parallelotope). Let Λ ∈ R 2n×n be a template matrix such that, for each i ∈ {1, 2, . . . , n}, Λ i = -Λ i+n , and let c ∈ R 2n . The parallelotope generated by Λ and c is the polytope:

⟨Λ, c⟩ = {x | Λx ≤ c}. (12)
Example 6. We consider again the SIR model. A parallelotope over the variables of the system is given by:

Λ =         -1 0 0 -1 -1 0 0 0 -1 1 0 0 1 1 0 0 0 1         c =         -0.8 -0.95 0 0.85 1 0        
According to the definition, each line in Λ and c defines a linear inequality on S, I, and R which corresponds to a half-space in S × I × R. For instance, the second line defines the inequality -S -I ≤ -0.95; the associated half-space consists of all the points that do not lay below the green line in Figure 4. Since all the inequalities belong to the same system defined by ⟨Λ, c⟩, all of them must hold at the same time. Thus, the parallelotope is the intersection of all the corresponding half-spaces.

In this example, ⟨Λ, c⟩ constrains the variable S to lay in [0.8, 0.85], the variable R to be 0, and the sum of S + I to belong to the interval [0.95, 1]. Figure 4 depicts a projection of ⟨Λ, c⟩ itself on the plane R = 0.

The above representation is also called constraint representation. Another way to characterize parallelotopes, similar to the one adopted for zonotopes [START_REF] Scott | Regular polytopes[END_REF], is to fix a point of origin and use vectors to define the parallelotope. From Equation [START_REF] Bogomolov | Juliareach: a toolbox for set-based reachability[END_REF], it is immediate that boxes are parallelotopes. be a vector of coefficients. We consider the set G of generators defined as

G = {g 1 = β 1 u 1 , g 2 = β 2 u 2 , . . . , g n = β n u n }.
Moreover, let q be a point in R n . The parallelotope generated by G and q is:

P gen (G, q) =    q + n j=1 λ j g j | (λ 1 , . . . , λ n ) ∈ [0, 1] n    . ( 13
)
This representation is called generator representation. The vectors g j are called generators of the parallelotope and q is called base vertex. Given a set of generators G = {g 1 , g 2 , . . . , g n } and a base vertex q, we also represent the parallelotope generated by G and q with the notation:

P gen (G, q) = {γ (q,G) (λ 1 , . . . , λ n ) | (λ 1 , . . . , λ n ) ∈ [0, 1] n }, (14)
where γ (q,G) (λ 1 , . . . , λ n ) is the linear function defined as:

γ (q,G) (λ 1 , . . . , λ n ) = q + n j=1 λ j g j . ( 15
)
The generator representation of parallelotopes resembles the way in which zonotopes are usually defined, i.e., as the sets {c+

m j=1 λ j g j | (λ 1 , . . . λ m ) ∈ [-1, 1]
m } where c is the center of the zonotope and g 1 , . . . , g m are its generators. This similarity is not surprising as the parallelotopes in R n are special kinds of zonotopes having exactly n linearly independent generators. However, it is worth noticing that our representation of parallelotopes uses base vertices as opposed to the zonotope definition which refers to the zonotope centers. This choice affects the domain of the coefficient λ 1 , . . . λ n , which is reduced from [-1, 1] n to [0, 1] n in parallelotope representation, and emphasizes that any parallelotope can be seen as the affine transformation of the unit box [0, 1] n . The reachability algorithm described in the remaining part of this section extensively uses this property.

Standard algebraic procedures can be used to switch from one representation to the other and back (e.g., see [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF]).

Example 7. Let us consider the parallelotope described in Example 6. In order to compute a parallelotope generator representation for it, we can proceed as follows. We consider as base vertex the point of the parallelotope obtained by minimizing all the variables, i.e., q = (0.8, 0.15, 0). We compute the vectors that go from q to the other vertexes and normalize them. So we obtain U = {(0.7071, -0.7071, 0), (0, 1, 0), (0, 0, 1)} and β = (0.07071, 0.05, 0). This means that the parallelotope is defined as the set of points whose coordinates satisfy:

q + λ 1 g 1 + λ 2 g 2
where q = (0.8, 0.15, 0), g 1 = (0.05, -0.05, 0), g 2 = (0, 0.05, 0), and λ 1 , λ 2 ∈ [0, 1]. In this case, g 3 and λ 3 disappear since β 3 = 0. This is because R has a single value instead of ranging in a proper interval in Example 6.

S I

Figure 5: The projection on the plane R = 0 of the paralletope ⟨Λ, c⟩ described in Example 6, its base vertex q = (0.8, 0.15, 0), and the generators g 1 = (0.05, -0.05, 0) and g 2 = (0, 0.05, 0).

Figure 5 depicts the projection of the plane R = 0 of the parallelotope de-scribed in Example 6 together with the base vertex q = (0.8, 0.15, 0) and its generators.

Parallelotopes can be used to decompose polytopes through parallelotope bundles.

Definition 10 (Parallelotope Bundle). A parallelotope bundle B = {P 1 , . . . , P b } is a finite set of parallelotopes whose intersection, denoted by

I(B) = ∩ b i=1 P i , generates a polytope.
The polytope I(B) represented by a parallelotope bundle B corresponds to the set of all the templates and offsets of the parallelotopes that constitute B itself. The following Lemma can be easily proved (see [START_REF] Dreossi | Reachability computation for polynomial dynamical systems[END_REF]) and establishes an upper bound on the number of parallelotopes necessary to represent a polytope.

Lemma 1 (Polytope Decomposition). Any polytope Q ⊂ R n defined by m constraints can be represented by a bundle involving ⌈m/n⌉ parallelotopes.

Example 8. Intuitively the polytope described in Example 5 can be decomposed in the two parallelotopes defined as follows:

• in the parallelotope P 1 the variable S ranges in [0, 0.98], the variable I ranges in [0, 0.98], and the variable R ranges in [0, 0.02];

• in the parallelotope P 2 the sum S + I ranges in [0, 0.98], the variable I ranges in [0, 0.98], and the variable R ranges in [0, 0.02].

As far as the domain for sets of parameters is concerned, we refer again to polytopes: P ⊆ P is going to be a polytope. While the described representations for sets of states in terms of boxes and parallelotopes are at the basis of our algorithm, the polytopes representing sets of parameters will be handled as systems of linear inequalities with variables in P.

Admissible Functions: Polynomials

The family of functions that we consider for defining dynamical systems is based on polynomials.

Definition 11 (Parametric Polynomial System).

A parametric polynomial system D = (X , P, f ) is a parametric discrete-time dynamical system in which f : X × P → X is such that f (x, p) = (f 1 (x, p), . . . , f n (x, p)) and each of the f i is a polynomial and it is linear over p.

Bernstein coefficients and their properties are valuable tools for bounding polynomials. In [START_REF] Dreossi | Reachability computation for polynomial dynamical systems[END_REF], we investigated the use of Bernstein coefficients for the over-approximation of bounded reachability in the case of polynomial systems without parameters. Since this work deals with parametric dynamical systems, we want to bound parametric polynomials of the kind f i (x, p) : R n × R m → R. Thus, we have to extend the definition of Bernstein basis to the parametric case. ⟨D, c⟩ can be decomposed in the two parallelotope P 1 and P 2 and, as a consequence, in the bundle B = {P 1 , P 2 }. The intersection of the bundle parallelotopes equals the polytope, i.e.,

I(B) = P 1 ∩ P 2 = ⟨D, c⟩
We will see that all the presented properties related to the Bernstein coefficients also hold when parameters are involved.

A multi-index is a vector i = (i 1 , i 2 , . . . , i n ) where each i j ∈ N. Given two multi-indices i and d of the same dimension, we write i ≤ d (d dominates i) if for all j ∈ {1, 2, . . . , n}, i j ≤ d j . We also write i/d for the multi-index (i 1 /d 1 , i 2 /d 2 , . . . , i n /d n ) and d i for the product of the binomial coefficients d1 i1 d2 i2 . . . dn in . Let us now focus on each of the components of the function f . To avoid using further indexes we denote it as π(x, p) : R n × R m → R.

A parametric polynomial π(x, p) : R n × R m → R is usually represented in the power basis as:

π(x, p) = i∈I π a i (p)x i , (16)
where I π is the set of multi-indexes of the polynomial π over the variables x and the coefficients a i (p) : R m → R, for i ∈ I f , are functions defined over the parameters p ∈ R m . In the case of parametric polynomial systems we further required a linear dependence over P, i.e., the coefficients a i (p) are linear in p.

Example 9. Consider the parametric polynomial π(

x 1 , x 2 , p 1 ) = 1/3 p 1 x 2 1 - 1/2 p 1 x 2 + 1/4 x 1 x 2 + 1/2. Its non-null parametric coefficients are a (2,0) (p 1 ) = p 1 /3, a (0,1) (p 1 ) = -p 1 /2, a (1,1) (p 1 ) = 1/4, and a (0,0) (p 1 ) = 1/2.
Bernstein basis polynomials of degree d are basis for the space of polynomials of degree at most d over R n . For x = (x 1 , x 2 , ..., x n ) ∈ R n , the i-th Bernstein polynomial of degree d is defined as:

B (d,i) (x) = β d1,i1 (x 1 )β d2,i2 (x 2 ) . . . β dn,in (x n ) ( 17
)
where, for any real number x ∈ R and for any d, i ∈ N,

β d,i (x) = d i x i (1 -x) d-i . ( 18
)
Example 10. It is easy to verify that β 0,0

(x) = 1, β 1,0 (x) = 1 -x, β 1,1 (x) = x, β 2,0 (x) = (1 -x) 2 , β 2,1 (x) = 2x(1 -x), and β 2,2 (x) = x 2 . Thus, B ((2,1),(0,1)) (x) = β 2,0 (x 1 )β 1,1 (x 2 ) = (1 -x 1 ) 2 x 2 .
Analogously, B ((2,1),(0,0)

) (x) = (1 -x 1 ) 2 (1 -x 2 ), B ((2,1),(1,0)) (x) = 2x 1 (1 - x 1 )(1 -x 2 ), B ((2,1),(1,1)) (x) = 2x 1 (1 -x 1 )x 2 , B ((2,1),(2,1)) (x) = x 2 1 x 2 , and B ((2,1),(2,0)) (x) = x 2 1 (1 -x 2 ).
Any parametric polynomial π(x, p) : R n × R m → R can be represented using Bernstein basis as:

π(x, p) = i∈I π b i (p)B (d,i) (x) (19)
where, for each i ∈ I π , the parametric Bernstein coefficient, is defined as:

b i (p) = j≤i i j d j a j (p). ( 20
)
The parametric Bernstein coefficients are functions of the form b i (p) : R m → R. However, the b i (p) have a linear dependence on p.

Example 11. Consider again the parametric polynomial presented in Example 9, i.e., π(

x 1 , x 2 , p 1 ) = 1/3p 1 x 2 1 -1/2p 1 x 2 + 1/4x 1 x 2 + 1/2.
For illustrative purposes, we study only the multi-index (1, 1), the corresponding parametric Bernstein coefficient is:

b (1,1) (p 1 ) = (1,1) (1,1) (2,1) (1,1) a (1,1) (p 1 ) + (1,1) (0,1) (2,1) (0,1) a (0,1) (p 1 ) + (1,1) (1,0) (2,1) (1,0) a (1,0) (p 1 ) + (1,1) (0,0) (2,1) (0,0) a (0,0) (p 1 ) = 1 1 1 1 2 1 1 1 a (1,1) (p 1 ) + 1 0 1 1 2 0 1 1 a (0,1) (p 1 ) + 1 1 1 0 2 1 1 0 a (1,0) (p 1 ) + 1 0 1 0 2 0 1 0 a (0,0) (p 1 ) = 1 1 2 1 a (1,1) (p 1 ) + 1 0 2 0 a (0,1) (p 1 ) + 1 1 2 1 a (1,0) (p 1 ) + 1 0 2 0 a (0,0) (p 1 ). Since n 0 = 1 and n 1 = n for any n ∈ N, b (1,1) (p 1 ) = 1 2 a (1,1) (p 1 ) + 1a (0,1) (p 1 ) + 1 2 a (1,0) (p 1 ) + 1a (0,0) (p 1 ).
As noticed in Example 9, the non-null parametric coefficients of the investigated polynomial are a (2,0) (p 1 ) = p 1 /3, a (0,1) (p 1 ) = -p 1 /2, a (1,1) (p 1 ) = 1/4, and a (0,0) (p 1 ) = 1/2. Thus,

b (1,1) (p 1 ) = 1 2 1 4 + 1(-p 1 /2) + 1 2 0 + 1 1 2 = 1 8 + 4 8 -p 1 /2 = 5/8 -p 1 /2.
Analogously, we obtain the parametric Bernstein coefficients b (0,0

) (p 1 ) = 1/2, b (0,1) (p 1 ) = (1 -p 1 )/2, b (1,0) (p 1 ) = 1/2, b (2,0) (p 1 ) = p 1 /3 + 1/2, and b (2,1) (p 1 ) = 3/4 -p 1 /6.
Since the parametric polynomial π(x 1 , x 2 , p 1 ) has the degree d = (2, 1), its Bernstein representation is:

π(x 1 , x 2 , p 1 ) =b (0,0) (p 1 )B (d,(0,0)) (x) + b (0,1) (p 1 )B (d,(0,1)) (x) + b (1,0) (p 1 )B (d,(1,0)) (x) + b (1,1) (p 1 )B (d,(1,1)) (x) + b (2,0) (p 1 )B (d,(2,0)) (x) + b (2,1) (p 1 )B (d,(2,1)) (x) = 1 2 (1 -x 1 ) 2 (1 -x 2 ) + 1 -p 1 2 (1 -x 1 ) 2 x 2 + 1 2 (2x 1 (1 -x 1 )(1 -x 2 )) + 5 8 - p 1 2 (2x 1 (1 -x 1 )x 2 ) + p 1 3 + 1 2 (x 2 1 (1 -x 2 )) + 3 4 - p 1 6 (x 2 1 x 2 ).
We now extend the range enclosing property of Bernstein coefficients to the parametric case. Let P ⊂ R m be a polytope.

Lemma 2 (Parametric Range Enclosing [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF]).

min i∈I π min p∈P b i (p) ≤ π(x, p) ≤ max i∈I π max p∈P b i (p), (21)
for all x ∈ [0, 1] n and p ∈ P , where b i (p), for i ∈ I π , are the parametric Bernstein coefficients of π.

Proof. Since we let p linearly range over P , this is an immediate consequence of the range enclosing property of Bernstein coefficients that states that the minimum and maximum Bernstein coefficients of a non-parametric polynomial π are a lower bound and an upper bound of the image of π over the unit box domain, respectively [START_REF] Cargo | The Bernstein form of a polynomial[END_REF].

Single Step Reachability Algorithm

In [START_REF] Dreossi | Reachability computation for polynomial dynamical systems[END_REF], we presented an algorithm for single-step reachability in the case of polytopes and polynomial functions without parameters. In particular, we exploited Bernstein coefficients to compute a polytope that over-approximates the image of a given polytope produced by a polynomial transformation. The idea can be generalized to the case of parametric polynomial functions provided that one can determine the minimum and the maximum of a parametric Bernstein coefficient over the set P of parameters. We recall here the main ingredients of the procedure and refer the reader to [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF] for all the details.

Let us start from the trivial case of the unit box [0, 1] n . So let X = [0, 1] n , P be a polytope, f = (f 1 , . . . , f n ) be polynomial, and proceed as follows:

1. for each j = 1, . . . , n compute the Bernstein coefficients b i (p) of the polynomial f j , with i ∈ I fj ; 2. for each Bernstein coefficient b i (p) of f j determine the minimum min j (i) = min p∈P b i (p) and the maximum max j (i) = max p∈P b i (p). This is possible since we are bounding a linear function over a polytope; 3. for each j = 1, . . . , n consider the minimum min j = min i∈I f j min j (i) and the maximum max j = max i∈I f j max j (i).

The output box, which over-approximates f

([0, 1] n , P ), is [min 1 , max 1 ] × • • • × [min n , max n ].
The reader should be aware that the parameter selection that minimizes or maximizes one of the Bernstein coefficients may differ from that that minimizes or maximizes another coefficient. So, while the returned box certainly over-approximates the set reachable in one step, the larger the parameter set, the coarser the over-approximation. In both the cases, if we want to compute a further reachability step, we cannot reuse the above-described method because the initial set is no more a unit box. The next paragraph presents a generalization of the method that deals with any kind of parallelotopes.

If we are instead interested in bounding the image of a parallelotope ⟨Λ, c⟩ we have to first rely on its generator representation, e.g., following these steps:

1. compute the generator representation of ⟨Λ, c⟩. Let γ (q,G) (λ 1 , . . . , λ n ) be the linear function over the unit box that occurs in such representation; 2. for each k = 1, . . . , 2n let Λ k be k-th row of the template matrix Λ.

Consider the function Λ k (f (γ (q,G) (λ 1 , . . . , λ n ), p)) over [0, 1] n × P ; 3. determine the Bernstein coefficients of Λ k (f (γ (q,G) (λ 1 , . . . , λ n ), p)). Identify the maximum c ′ k of such coefficients over P . Since both γ (q,G) and Λ k are linear, in the function Λ k (f (γ (q,G) (λ 1 , . . . , λ n ), p)) the dependence on p is still linear and the maximum can be found exploiting linear optimization techniques. Moreover, the property of Bernstein coefficients ensures that Λ k (f (γ (q,G) (λ 1 , . . . , λ n ), p)) ≤ c ′ k over [0, 1] n × P , and hence, Λ k (f (x, p)) ≤ c ′ k over ⟨Λ, c⟩ × P . The output parallelotope which over-approximates f (⟨Λ, c⟩, P ) is ⟨Λ, c ′ ⟩, where c ′ = (c ′ 1 , . . . , c ′ 2n ). The procedure is sketched in Algorithm 2.

Algorithm 2 Parallelotope-based reachability step

1: function reachStep(X = ⟨Λ, c⟩, P, f ) 2:
for k ∈ {1, . . . , 2n} do 3:

γ (q,G) (λ 1 , . . . , λ n ) ←con2gen(X) ▷ Compute generator function 4: c ′ k ←maxBernCoeff(Λ k (f (γ (q,G) (λ 1 , . . . , λ n ), p)), P ) ▷ Compute maximum coefficient
Λ =         -1 0 0 0 -1 0 0 0 -1 1 0 0 0 1 0 0 0 1         c =         0 0 0 1 0.7 1.6         .
This means that q = (0, 0, 0) and γ (q,G) (λ 1 , λ 2 , λ 3 ) = (λ 1 , 0.7λ 2 , 1.6λ 3 ). The Bernstein coefficients of the polynomials in the variables λ 1 , λ 2 , λ 3 are: If the initial set is (S, In the case of polytopes, the bundle decomposition allows us to efficiently solve the problem of defining a function from [0, 1] n to a generic polytope Q. However, this is crucial only for transforming the input, while the directions for computing the new offset vector can be changed or combined differently in the bundles during the computation. Consequently, in [START_REF] Dreossi | Reachability computation for polynomial dynamical systems[END_REF], we described two different strategies for computing the new bundle. One of them is faster to compute, while the other is the more precise of the two. The above considerations, that allow us to generalize the methods presented in [START_REF] Dreossi | Reachability computation for polynomial dynamical systems[END_REF] to the parametric case, also apply to such strategies [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF].

ρ S = λ 1 -0.7βλ 1 λ 2 ρ R = 0.7λ 2 + 0.7βλ 1 λ 2 -0.7αλ 2 ρ I = 1.6λ 3 + 0.7αλ 2 The Bernstein coefficients of ρ S are in {0, 1, 1 -0.7β}, those of ρ R are in {0, 0.7(1 -α), 0.7(1 + β -α)},
I, R) ∈ [0, 1] × [0, 0.65] × [0, 1.55], γ (q,G) (λ 1 , λ 2 , λ 3 ) is (λ 1 , 0.65λ 2 , 1.
The presented algorithms correctly compute an over-approximation of the reachable set, as proved by Lemma 2. It is worth noticing that the proposed strategy is memory-less and the evaluated evolution from a point does not consider the parameter values that let that point be reached. While this approach simplifies the parametric reachability problem and, in some sense, makes it solvable, it produces an over-approximation of the solution that gets coarser and coarser as the evolution proceeds because, at each step, we admit any value in the parameter set as a valid choice for the parameter values. (2) = (0.9786718, 0.0167382, 0.00459) and ξ (2) = (0.98162695, 0.00481305, 0.01356).

If we over-approximate the reachable set by avoiding to select α once for all and preserving the choice between α = 1/5 and α = 2/5 at each evolution step, we produce an infinite number of "spurious" trajectories and, in particular, in two steps we can reach ξ (4/5,1/2) x1

(1) = (0.9786718, 0.0089682, 0.01236) where

x 1 = ξ (1/5,1/2) x0
(1). While the proposed algorithm exclusively deals with dense sets of parameters, it applies the just-mentioned approach to over-approximate the flowpipe. Mixing α

x 1 x 2 x 0 = x ′ 0 x ′ 1 x ′ 2 ξ p ′ x1 (1)
ξ p x ′ 1 (1
Figure 8: Each step of the reachability procedure admits any possibile value in the parameter set. Thus, when the parameter set is not a singleton, the algorithm follows also "spurious" trajectories along which the parameters changes. The trajectory ξ p x 0 described in Example 1 for p = (1/5, 1/2) and p ′ = (4/5, 1/2). The points ξ p x 0 (i) and ξ p ′ x 0 (i) are labelled as x i and x ′ i , respectively. The triangles depict two points that are included in the over-approximation of the reachable set exclusively because the algorithm selected two different values for α during the two-time step horizon.

As far as computational complexity is concerned, the most expensive part of the computation is the Bernstein coefficients' computation. However, since the template parallelotopes do not change during the evolution of the reach set, the coefficients can be symbolically computed once for all and instantiated at each step. This approach significantly reduces the complexity (see [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF]).

Parameter Synthesis

The parameter synthesis problem requires handling flows of parametric traces that, in most cases, cannot be exactly computed, as we saw in the previous section. Therefore, we need to find a compromise between precision and tractability of the problem. For this reason, we introduce the notion of approximated sets of trajectories which allows us to deal with the synthesis problem.

Let D = (X , P, f ) be a discrete-time dynamical system, X 0 ⊆ X , and P ⊆ P. The set of trajectories Ξ(X 0 , P ) of D can be over-approximated by the flowpipe

W (X0,P ) = {X t | t ∈ N and for t > 0 it holds X t = f (X t-1 , P )}.
The over-approximating flowpipe W (X0,P ) can be computed with the setintegration algorithms presented in Section 2. However, it is worth noticing that this is an over-approximation of Ξ(X 0 , P ) since the relationship between a single trace and its corresponding parameter is not kept.

We can define a semantics on flowpipes that reflects the parameter synthesis problem. In particular, we define the following synthesis semantics.

Definition 12 (Synthesis Semantics). Let D = (X , P, f ) be a discrete-time dynamical system, W (X0,P ) be a flowpipe with X 0 ⊆ X and P ⊆ P, t ∈ N be a time instant, and φ be an STL formula in positive normal form. The synthesis semantics Φ(φ, W (X0,P ) , t) of φ at time t over the flowpipe W (X0,P ) is given by the following inductive definition: Φ(σ, W (X0,P ) , t) = P σ , where P σ ⊆ P is the largest subset s.t. ∀x 0 ∈ W (X0,P ) (t), ∀p ∈ P σ , σ(f (x 0 , p)) is true Φ(φ 1 ∧ φ 2 , W (X0,P ) , t) = Φ(φ 1 , W (X0,P ) , t) ∩ Φ(φ 2 , W (X0,P ) , t) Φ(φ 1 ∨ φ 2 , W (X0,P ) , t) = Φ(φ 1 , W (X0,P ) , t) ∪ Φ(φ 2 , W (X0,P ) , t)

Φ(φ 1 U I φ 2 , W (X0,P ) , t) = t ′ ∈t+I (Φ(φ 2 , W (X0,P ) , t ′ )∩ t ′′ ∈[t,t ′ ] Φ(φ 1 , W (X0,P ) , t ′′ )) (22)
From an intuitive point of view, Φ(φ, W (X0,P ) , t) returns a subset P φ ⊆ P of parameters that ensures that φ is satisfied at time t starting from any point in X 0 and assigning to the parameters any value in P φ . Again, the synthesis semantics at time t returns a set P φ that steers the system from time t + 1 on. This is slightly counterintuitive since the semantics usually returns evaluations referring to the time instant in which they are applied. However, it is consistent with our choice on the semantics of STL.

We say that a flowpipe W (X0,P ) satisfies a formula φ, denoted with W (X0,P ) |= φ if and only if Φ(φ, W (X0,P ) , t) = P . By using Definition 12 and structural induction on φ, it is easy to see that Φ(φ, W (X0,P ) , t) is idempotent on P , i.e., Φ(φ, W (X0,Pφ) , t) = P φ where P φ = Φ(φ, W (X0,P ) , t). Thus, W (X0,Pφ) |= φ. This statement, which is formally reported in the following theorem, asserts that the synthesis semantics provides a refined set of parameters that satisfies the formula φ.

Theorem 1. If Φ(φ, W (X0,P ) , 0) = P φ , then W (X0,Pφ) |= φ.

Since W (X0,P ) over-approximates Ξ(X 0 , P φ ) with P φ ⊆ P , if φ is satisfied by W (X0,P ) , then it is satisfied also by Ξ(X 0 , P φ ). The following lemma shows that the parameters generated by the synthesis semantics are also valid for the exact trajectories of a dynamical system. Lemma 3. If Φ(φ, W (X,P ) , t) = P φ , then for each x ∈ X and for each p ∈ P φ it holds that ξ p x , t |= φ.

Proof. By structural induction on φ.

• (σ) Let p ∈ P σ , x ∈ X, x t-1 = ξ p x (t -1), and X t-1 = W (X,P ) (t -1). Since P σ ⊆ P we have that x t-1 ∈ X t-1 and x t = f (x t-1 , p) ∈ f (X t-1 , P σ ). Hence, by definition of P σ we have that σ(x t ) is true, i.e., the thesis;

• (φ 1 ∧ φ 2 ) Let Φ(φ 1 , W (X0,P ) , t) = P φ1 and Φ(φ 2 , W (X0,P ) , t) = P φ2 . We have that P φ = P φ1 ∩ P φ2 . Let p ∈ P φ and x ∈ X, since p belongs to both P φ1 and P φ2 , by inductive hypothesis we have that ξ p x (t) |= φ 1 and ξ p x (t) |= φ 2 . Hence, we get the thesis; • (φ 1 ∨ φ 2 ) Similar to the conjunction;

• (φ 1 U I φ 2 ) By definition, if p ∈ P φ , there exists t ′ ∈ t + I such that p ∈ P t ′ φ2 ∩ t ′′ ∈[t,t ′ ] P t ′′ φ1 , where P t ′ φ2 = Φ(φ 2 , W (X0,P ) , t ′ ) and P t ′′ φ1 = Φ(φ 1 , W (X0,P ) , t ′′ ). Hence, by inductive hypothesis on φ 1 and φ 2 , we get the thesis.

As a consequence, we are correctly under-approximating the original parameter synthesis problem, as stated by the following theorem.

Theorem 2. If Φ(φ, W (X0,P ) , 0) = P φ , then Ξ(X 0 , P φ ) |= φ.

Proof. This is an immediate consequence of Lemma 3.

In summary, we proved that all the trajectories starting in X 0 and having parameters in P φ , where P φ is the result of the synthesis semantics of φ on the flowpipe W (X0,P ) , satisfy the formula φ. Thus, the synthesis semantics can be used to produce under-approximations of the parameter synthesis problem (see Definition 4).

In the rest of this section, we will see how to synthesize a set P ′ such that Ξ(X 0 , P ′ ) |= φ in the case of the domains and systems defined in Section 2. Our goal will be achieved by exploiting the synthesis semantics and, in particular, Φ(φ, W (X0,P ) , 0). Even though this last set would be the largest subset of P satisfying the specification, it will be enough to under-approximate it, rather than exactly compute it, because, if Ξ(X 0 , P ) |= φ, then Ξ(X 0 , P ′ ) |= φ too for any subset P ′ of P .

We will proceed by induction on the structure of φ and rely on the main procedure described in Algorithm 3 which takes in input a set of states X 0 ⊆ X , a set of parameters P ⊆ P, a function f describing the evolution of the system, and an STL formula φ in positive normal form, and produces a refinement P φ ⊆ P .

Algorithm 3 Parameter synthesis.

1: function ParaSynth(X 0 , P, f , φ)

2: if φ = σ then ▷ Predicate 3:
return RefPredicate(X 0 , P, f , σ)

4:
end if

5: if φ = φ 1 ∧ φ 2 then ▷ Conjunction 6:
return ParaSynth(X 0 , P, f , φ 1 ) ∩ ParaSynth(X 0 , P, f , φ 2 )

7:
end if

8: if φ = φ 1 ∨ φ 2 then ▷ Disjunction 9:
return ParaSynth(X 0 , P, f , φ 1 ) ∪ ParaSynth(X 0 , P, f , φ

end if

11: if φ = φ 1 U I φ 2 then ▷ Until 12:
return UntilSynth(X 0 , P, f , φ 1 U I φ 2 )

13:

end if 14: end function

Basic Refinement Step φ ≡ σ

Let us consider the case of sets of points in the domain of parallelotopes, set of parameters in the domain of polytopes and polynomial functions with a linear dependence on the parameters. Let us also assume that σ ≡ s(x) ≤ 0, where s is a linear function. In this case, similarly to what has been done in Section 2.4, we can compute the Bernstein coefficients that are going to be linear in the parameters. By imposing the parametric Bernstein coefficients to be nonpositive, we obtain a set of new linear constraints that we add to the definition of the set of parameters. Specifically, if X 0 = ⟨Λ, c⟩, the RefPredicate procedure operates as follows:

1. compute the generator representation of ⟨Λ, c⟩. Let γ (q,G) (λ 1 , . . . , λ n ) be the linear function over the unit box that occurs in such representation; 2. consider the function s(f (γ (q,G) (λ 1 , . . . , λ n ), p)). The image of this function over [0, 1] n × P coincides with the image of s(f (x, p)) over ⟨Λ, c⟩ × P ; 3. determine the Bernstein coefficients b i (p) of s(f (γ (q,G) (λ 1 , . . . , λ n ), p)); 4. add to P all the constraints of the form b i (p) ≤ 0.

We now present two examples, dealing with different formulas and two sets of initial conditions, to clarify the RefPredicate procedure. First of all, the generator for [0, 1] × [0, 0.7] × [0, 1.6], i.e., γ (q,G) (λ 1 , λ 2 , λ 3 ) = (λ 1 , 0.7λ 2 , 1.6λ 3 ), is found (see the first part of Example 14). Then, the algorithm computes s(f (γ (q,G) (λ 1 , . . . , λ n ), p)), where s(x, p) ≤ 0 is the constraint to be satisfied. In this example, s(f Example 14). Finally, the algorithm finds s(f (γ (q,G) (λ 1 , . . . , λ n ), p)) Bernstein coefficients, i.e., {-0.365, 0.335 -0.7α, 0.335 + 0.7β -0.7α}, and refines the parameter set by using the constraints b i (p) ≤ 0, for each coefficient b i (p). Thus, the synthesized parameter set is such that α ∈ [0.5, 0.6], β ∈ [0.1, 0.2], -0.365 ≤ 0, 0.335 -0.7α ≤ 0, and 0.335 + 0.7β -0.7α ≤ 0, which is equivalent to α ≤ 0.6, β ≥ 0.1, and α -β ≥ 67/140.

(γ (q,G) (λ 1 , . . . , λ n ), p)) is 0.7(λ 2 + βλ 1 λ 2 -αλ 2 ) -0.365 (see the first part of
If, instead, the initial conditions are S ∈ [0, 1], The generator representation γ (q,G) (λ 1 , λ 2 , λ 3 ) equals (λ 1 , 0.7λ 2 , 1.6λ 3 ), and s(f (γ (q,G) (λ 1 , . . . , λ n ), p)) is 1.6λ 3 +0.7αλ 2 -2 (see the first part of Example 14). The Bernstein coefficients of s(f (γ (q,G) (λ 1 , . . . , λ n ), p)) are -2,-0.4, 0.7α -2, and 0.7α -0.4 and the algorithm returns the set such that α ∈ [0.5, 4/7] and β ∈ [0.1, 0.2].

I ∈ [0, 0.65], R ∈ [0, 1.55], β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55], γ (q,G) (λ 1 , λ 2 , λ 3 ) = (λ 1 , 0.65λ 2 , 1.55λ 3 ) (see Example 14) and s(f (γ (q,G) (λ 1 , . . . , λ n ), p)) is 0.65(λ 2 + βλ 1 λ 2 -αλ 2 ) -0.
If, instead, the initial conditions are

S ∈ [0, 1], I ∈ [0, 0.65], R ∈ [0, 1.55], β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55], γ (q,G) (λ 1 , λ 2 , λ 3 ) = (λ 1 , 0.65λ 2 , 1.55λ 3 ) and s(f (γ (q,G) (λ 1 , . . . , λ n ), p)) is 1.55λ 3 + 0.

Conjunction and Disjunction Formulas

In the case of formulas of the form φ 1 ∧ φ 2 , if ParaSynth(X 0 , P, f , φ 1 ) and ParaSynth(X 0 , P, f , φ 2 ) return two polytopes P φ1 and P φ2 defined as a set of linear constraints, it is sufficient to put together all the constraints and we obtain the set P φ1∧φ2 .

On the other hand, in the case of formulas of the form φ 1 ∨ φ 2 , we have to accumulate the outputs of the two calls. A similar issue will occur in the case of the until operator described in the next subsection. 3. a = 0 and b = 0: the interval coincides with the single time instant 0.

Intuitively, the function UntilSynth recursively transforms the cases 1 and 2 into the base case 3. Before defining the algorithm, it is worth to point out that a single until formula φ 1 U [a,b] φ 2 may require several basic refinements. Consider, for instance, the case where φ 1 always holds and φ 2 holds at several time instants inside [a, b]. Here, the number of basic refinements that φ 1 U [a,b] φ 2 requires equals the number of time instants in which φ 2 holds. This means that an until can admit several valid refinements. Some of the refined parameter sets might be included in others, but since we do not know it in advance, we need to compute and accumulate all the possible solutions.

Let us now analyse the three cases constituting the structure of UntilSynth (Algorithm 4):

1. a > 0 and b > 0: the until formula is satisfied if φ 1 holds until φ 2 is true inside the interval [a, b]. We first refine the parameters at time 0 over φ 1 , obtaining the subset P φ1 (Line 3). Then the algorithm performs a reachability step using the valid parameter set P φ1 to produce the new set X ′ (Line 4). Now the algorithm proceeds with the recursive call (Line 5). This can be seen as a step towards the interval [a, b], except that instead of restoring the synthesis from time 1, we shift the interval backwards; 2. a = 0 and b > 0: there are two ways to satisfy the until formula:

(a) φ 2 is satisfied right now at time 0; (b) φ 1 holds until φ 2 is satisfied before the time instant b.

In the first case, we need to refine the parameter set with respect to φ 2 .

In the second case, the algorithm refines with respect to φ 1 , then the procedure performs a reachability step under the refined parameter set P φ1 , obtaining the new set X ′ . Similarly to the previous case, we execute Algorithm 4 Until synthesis.

1: function UntilSynth(X 0 , P, f , φ 1 U [a,b] φ 2 ) 2:
if a > 0 and b > 0 then ▷ Outside interval

3: P φ1 ← ParaSynth(X 0 , P, f , φ 1 ) ▷ Check φ 1 4: X ′ ← ReachStep(X 0 , P φ1 , f ) 5: return UntilSynth(X ′ , P φ1 , f , φ 1 U [a-1,b-1] φ 2 ) 6: end if 7:
if a = 0 and b > 0 then ▷ In interval

8: P φ1 ← ParaSynth(X 0 , P, f , φ 1 ) ▷ Check φ 1 9: P φ2 ← ParaSynth(X 0 , P, f , φ 2 ) ▷ Check φ 2 10: X ′ ← ReachStep(X 0 , P φ1 , f ) 11: return P φ2 ∪ UntilSynth(X ′ , P φ1 , f , φ 1 U [a,b-1] φ 2 ) 12:
end if

13:
if a = 0 and b = 0 then ▷ Base 14:

return ParaSynth(X 0 , P, f , φ

end if 16: end function a step forward by shortening the interval by one. The procedure then returns the union of P φ2 and the result provided by the recursive call;

F [1,1] (I ≤ 0.365 ∧ R ≤ 2)
does hold. First of all, this formula is translated into the corresponding Ubased formula, ⊤U [START_REF] Althoff | Combining zonotopes and support functions for efficient reachability analysis of linear systems[END_REF][START_REF] Althoff | Combining zonotopes and support functions for efficient reachability analysis of linear systems[END_REF] (I ≤ 0.365 ∧ R ≤ 2), as explained in Section 1.2. Then, Algorithm 4 performs the following steps:

1. synthesize the parameters for the formula ⊤; 2. compute an over-approximation of the set reachable by using the evaluated parameters; 3. further restrict α and β by synthesizing their admissible values for the formula ⊤U [0,0] (I ≤ 0.365 ∧ R ≤ 2) by using as initial set that computed during step 2.

The formula ⊤ always holds by definition, thus, step 1 does not change the parameter set. So, step 2 is equivalent to that reported in the first part of Example 13 and it computes the parallelotope (S,

I, R) ∈ [0, 1] × [0, 0.7] × [0, 1.6].
Then, Algorithm 4 reduces step 3 to synthesizing a set for the formula I ≤ 0.365 ∧ R ≤ 2 by using the just-mentioned conditions as initial region. We know that the resulting set is empty from the first part of Example 18. The correctness of the above described algorithms in synthesizing a set of parameters that guarantees the STL specification is ensured by the correctness of the reachability algorithm that at each step over-approximates the set of reachable points together with Theorem 2.

As far as the computational complexity is concerned, in [START_REF] Dang | Parameter synthesis through temporal logic specifications[END_REF] it has been shown that the synthesis could require a number of reachability steps, basic refinements, unions and intersections which is exponential with respect to the dimension of the STL formula. The dimension of a formula includes both the number of operators and the size of the observed time horizon. Some scalability tests are available in [START_REF] Dang | Parameter synthesis through temporal logic specifications[END_REF][START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF], where time intervals having length up to 120 and formulas with increasing nesting of until operators up to 20 have been considered on an Ebola system involving 5 variables and 6 parameters (see Section 6.2.4 and Tables 6.1 and 6.2 in [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF]).

Parameter Splits

The described method guarantees the synthesis of a subset of parameters satisfying the specification over all the possible starting points. In many cases, we can obtain an empty set of parameters because of the over-approximations of the reached points accumulated during the computation. In particular, we do not store, for each point x in the reach set, the parameters that have allowed the model to flow to that very point and we propagate x in the future steps by using all the parameter values that are considered valid up to that step. This approach is hazardous when the initial set of parameters is large and produces an extremely conservative over-approximation of the real reach set.

Moreover, our procedure refines only the parameters directly involved in the dynamics of the variables occurring in the specification. Thus, when a parameter does not appear in the dynamic of a variable v occurring in the specification, but it rules the evolutions of other variables involved in the dynamic of v itself, we cannot refine the parameter set and we return an empty answer.

We take care of both issues by introducing an automatic split procedure over the parameters set. If we have not identified any parameter satisfying the specification, sapo resets the computation to time 0 and autonomously splits the set of parameters into two subsets over each of the hyperplanes defining it. For instance, we get four sub-rectangles in the case of two parameters ranging over a rectangle. Then, each subset is refined independently. The splitting step is iterated until either a non-empty parameter set is synthesized or the maximum number of iterations allowed by the user has been reached. Since each split exponentially increases the computational requirements, users should carefully choose a trade-off between efficiency and precision.

sapo and Sapo Input Language (SIL)

Originally, Sapo was meant to be a library rather than a stand-alone tool. Users were forced to represent models using the C++ programming language and write the code for calling the aimed analysis methods. Moreover, a new compilation was required every time either the model or the reachability/synthesis problem changed. In order to simplify the usage of Sapo and broaden its user base, we developed a stand-alone application, named sapo, that avoids recompilations by accepting in input problems described by the Sapo Input Language (SIL for short). This tool is an efficient and practical alternative to the Sapo library for non-expert users; it also embeds for free some nice features -such as multi-threading computations-that, during the writing of this article, were not fully integrated by the Sapo C++ API.

sapo takes in input SIL file which consist of four main sections:

Header which, among the other options, specifies whether a synthesis or a reachability analysis must be performed;

Symbols definition which declares the variables, the parameters, the constants, the dynamics, and, whenever needed, the STL specification;

Matrices definition which defines the template matrices and the initial set;

Footer which contains some Sapo-specific options to fine-tune its behaviour.

Pairs of rows in the parallelotope template matrices are linearly dependent and, in these pairs, each line can be obtained from the other one by multiplying the latter by -1 (e.g., see the template matrix Λ of Example 14).

In order to reduce specification redundancy and decrease user effort, SIL decouples the template matrices of all the bundles used during the computation into two components: the vector of the planes/directions and the template vector whose values specify the directions of each parallelotope. Under this convention, the initial set is declared by specifying an interval of admitted values for each direction.

Section 5.1 presents a case study together with some SIL examples. Please, refer to the SIL page of the project wiki for a complete overview of the syntax.

Thanks to SIL, the sapo application avoids the need for a recompilation at each new analysis. The class Driver interprets the input provided in SIL format and returns an InputData object to represent the model description and, possibly, its specification. As soon as it is executed, the sapo tool creates an instance of the Driver class and, through it, parses the input. The resulting InputData object is used to fill a Model object and create a sapo opt object to represent the computation options reported in the input. Then, the application creates an instance of the class Sapo and performs the reachability or the synthesis analysis by calling either the Sapo::reach or the Sapo::synthesize, respectively. Finally, it prints the results: in the case of reachability, users obtain a temporal sequence of sets representing the flowpipe, while, when the synthesis analysis is requested, the tool also shows the set of synthesized parameters. In both cases, the sets are represented as the union of the solution sets of some linear inequalities.

Web Application

With the main aim of relieving the user from installing Sapo and from updating it any time a new version is released, we developed a web interface for the sapo tool which, among the other features, can plot flowpipes and sets of synthesized parameters.

The web interface named webSapo, allows the user to input variables, system dynamics, and parameters. According to the user's choices for the domains of variables and parameters, it establishes which matrices must be input by users and automatically creates the others. If the addressed problem is parameter synthesis, it also offers a section for the STL specification.

Once the investigated problem has been analysed, the results are locally saved and they can be depicted in both 2D or 3D plots. The step-by-step evolution of the reached set can be also presented as a sequence of images.

Among all its features, webSapo also supports off-line computations: any sapo output file in JSON format can be loaded by webSapo and the analysis results can be plotted as they were just computed by webSapo itself.

The webSapo source code can be downloaded from https://github.com/ LucaDorigo/webSapo and a manual for it is available on the project wiki page. A free-to-use online webSapo service is offered at http://encase.uniud.it: 3001.

Case Studies

In this section we consider two case studies drawn from real scenarios meant to provide a better insight on the Sapo potentials. Here, we do not present any scalability test because, on the one hand, designing adequate tests having increasing complexity and avoiding empty parameter set scenarios is not an easy task, on the other hand, completely factitious models may miss to provide realistic results. Some scalability tests performed on Sapo are instead presented in [START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF].

COVID-19 Simple Model

The SIR example allows us to keep the description simple, show the main functionalities of Sapo, and, at the same time, make some consideration on the current COVID-19 outbreak.

Many models of COVID-19 can be found in the literature (e.g., see [START_REF] Arino | A simple model for COVID-19[END_REF][START_REF] Giordano | Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy[END_REF][START_REF] Engbert | Sequential data assimilation of the stochastic SEIR epidemic model for regional COVID-19 dynamics[END_REF]). Here we are not interested in depth investigating these models. Some of them are very complex, for instance, involving partial differential equations or stochasticity, and a full paper could be devoted to their implementation and analysis with Sapo. Instead, we aim to present our framework using fascinating yet simple examples.

Let us consider a more realistic SIR model with two extra parameters: Vaccination rate µ and loss of immunity rate γ. The system dynamics are:

S t+1 = S t -βS t I t -µS t + γR t I t+1 = I t + βS t I t -αI t R t+1 = R t + µS t -γR t + αI t
To capture the evolution of COVID-19 pandemic, let us assume that:

• the transmission rate β ranges in [0.055, 0.1]. This is equivalent to an R 0 rate ranging in [0, 77, 1.4];

• the recovery rate α ranges in [0.05, 0.07], i.e., between 14 days and 20 days are necessary for recovery;

• the vaccination rate µ ranges in [0.00001, 0.001], i.e., we assume that we can vaccinate at most 0.1% of the susceptible individuals each day. We consider that only a low fraction of the population can be vaccinated each day, i.e., we are assuming that we want to rely on the existing vaccination centers and not on the hub built during the emergency. Notice that to be more precise µ represents the rate of susceptible individuals that each day are vaccinated and over which the vaccine has efficacy;

• the loss of immunity rate γ ranges in [0.0027, 0.0055], representing that the immunity duration lays between 180 and 360 days.

In order to model the situation in Italy in October 2021 at time 0, we also assume that:

• R ranges in [0.7, 0.8]. This accounts of all the individuals that are either vaccinated or recovered from the disease;

• I ranges in [0.001, 0.1]. This range is a pessimistic estimation in which we assume a very high number of asymptomatic individuals;

• S ranges in [0.2, 0.3]. This value is due to the other choices.

One could be initially interested in forecasting the pandemic situation within 100 days, i.e., computing the set of points reachable within 100 steps. The model together with the just-mentioned problem is described by the following SIL file. The last three lines in the SIL file define the directions vector and the intervals vector representing the initial parallelotope. The expression between the reserved words direction and in could be any expression linear in s, i, and r; for instance, 3 * s/4 + r or s -2 * i + r/3.

Since our example uses the canonical directions, i.e., our parallelotopes are boxes, the last three lines can be omitted by specifying the box boundaries during the variable declarations as follows: It is easy to see that this number surges after around 80 days. Thus, one could be interested in synthesizing the parameters that ensure that the number of infected stays below 0.1 within the next 100 days. Since this last property corresponds to the STL formula G [0,100] i < 0.1, our goal can be achieved by replacing the word reachability with synthesis in the first line of the SIL file and, at the same time, by adding the line spec: G[0,100] i < 0.1; at the end of it.

Figure 13 depicts a graphical representation, produced by webSapo, of the synthesized values for α and β. The meaning of the result is quite intuitive: since we cannot revaccinate a high fraction of the population, we have at least to introduce restrictions on the circulation that reduce the transmission rate β below 0.9. Moreover, β has to be further reduced by introducing more stringent measures if we cannot provide treatments, such as antiviral drugs, that increase the recovery rate α. As a matter of fact, a low recovery rate increases the probability for an infected individual to transmit the disease to a susceptible one. Figure 14 depicts the evolution of the infected with respect to the synthesized parameters.

Finally, one could adjust the vaccination rate µ and loss of immunity rate γ. Let us assume that we can allow the transmission rate β to grow up to 0.15 (i.e., R 0 = 2.1), i.e., we would prefer to avoid restrictions. On the other hand, we are prepared for investing in massive revaccination, i.e., µ can grow up to 0.1 of the susceptibles. This approach is not realistic when most of the population is susceptible, but it becomes feasible after the first round of vaccination when the susceptibles are a low fraction of the population. Moreover, this is one of the parameters we are interested in synthesizing, so we let it range in a large interval and check whether it gets reduced after the computation. We are now interested in ensuring that the infected are less than 15% until the susceptible are less than 10% within the next 100 days, i.e., we are considering the specification (I < 0.15)U [0,100] (S < 0.1).

In this case, the parameter splitting procedure is necessary to obtain a nonempty result. Since sapo does not split parameter set during the synthesis by default, we added the line: max_parameter_splits: 3; in the header of the SIL file. This syntax sets the maximum number of admitted parameter splits; in this case, three. The first synthesis round considers the whole initial parameter set. If the synthesis produces a non-empty set, sapo returns the result. Otherwise, the parameter set is split in 2 n sub-sets and the process is repeated until either the maximum number of splits has been reached or sapo synthesizes a non-empty set of parameters. It is worth noticing that a non-empty set may be returned during the first round. However, in the worst case scenario, the synthesis takes exponential time with respect the maximal number of admitted splits. Because of this, we decided to let users tune this option according to their wishes.

Figure 15 represents the 2D projection of the synthesis over γ and µ and it proves that the parameter set was split during the computation. As expected, the picture shows that if the loss of immunity rate γ is high, i.e., the immunity period is short, we need to invest in revaccination drastically, i.e., µ has to be at least 0.075. Let us notice again that the splits are the key ingredient to getting a non-empty result.

The above described examples have been run on a MacBook Pro M1 2020, with 16GB of RAM, taking less than a second, a part from the last one that took around 34 seconds with multi-threading and 118 with single-threaded computation.

Enhancing Accuracy of Neural Networks Analysis

In this section we illustrate another potential application of Sapo, which is formal analysis of neural networks (NN). This problem has become crucial in the recent development of autonomous applications such as robots, self-driving vehicles, and medical devices. Indeed, neural networks are increasingly used in such applications thanks to their approximation power and ability to "learn" from data, demonstrated by their high performance in image classification, natural language processing and speech recognition (see for example [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Kamath | Deep Learning for NLP and Speech Recognition[END_REF]). However, assuring the safety and robustness of their usage, in particular under the impact of adversarial inputs or perturbations, is still challenging. Beyond classification purposes, NN can be used to approximate and implement complex control laws, and verifying a closed-loop system with such neural network controllers requires the ability to approximate the image of an input set by a neural network within a given error bound. Various set propagation and reachability computation from formal verification and abstract interpretation of dynamical systems have been applied to address this problem (see, e.g. [START_REF] Xiang | Output reachable set estimation and verification for multilayer neural networks[END_REF][START_REF] Dutta | Sherlock -a tool for verification of neural network feedback systems: Demo abstract[END_REF][START_REF] Huang | Reachnn: Reachability analysis of neural-network controlled systems[END_REF][START_REF] Ivanov | Verisig: verifying safety properties of hybrid systems with neural network controllers[END_REF][START_REF] Tran | NNV: the neural network verification tool for deep neural networks and learning-enabled cyber-physical systems[END_REF]). While the application of these approaches does not pose serious theoretical problems, computational complexity is an issue mainly because of the non-linearities of activation functions and the large size of practical neural networks. Piecewise-linear approximations, which are commonly used in these approaches, provide good accuracy for ReLU activation functions, but may lead to overly conservative results for sigmoid activation functions. Polynomial/rational approximations and the Bernstein techniques can thus be combined to enhance the accuracy for parts of the network with high sensitivity. This is what we want to demonstrate in this section. In particular, we focus on the usage of Sapo for reachability computation in the non-parametric case, while some considerations on possible extensions to the parameter synthesis functionalities conclude this section. Therefore, in [START_REF] Pippia | Image approximation for feed forward neural nets[END_REF][START_REF] Pippia | Optimization and Modeling Techniques for Food Service Appliances[END_REF], we propose an algorithm for accurately computing the image of activation functions, by first approximating them using polynomial and rational functions and then applying the Bernstein techniques on the resulting functions.

Let us now demonstrate the idea with feed-forward neural networks. Each layer l of such a network N has the form:

h l (x) = s(ℓ l (x)),
where ℓ l is a linear function, while s is an activation function. Given an initial set of values for x, the image of the neural networks is the result of the application of all the layers, denoted by N (x). First, our algorithm [START_REF] Pippia | Image approximation for feed forward neural nets[END_REF][START_REF] Pippia | Optimization and Modeling Techniques for Food Service Appliances[END_REF] approximates the functions h l by polynomials. Then it over-approximates the outputs of the activation functions h l using Sapo, parametric reachability computation. In [START_REF] Pippia | Image approximation for feed forward neural nets[END_REF], we show how hyperbolic tangent activation functions are approximated using polynomials. In particular, given a desired degree we can sample the hyperbolic tangent in the confident interval and we can tune the coefficients of the polynomial function in order to minimize the maximum absolute error, using, for instance, the polyfit function available in MATLAB ® . Considering degree 5, we obtained the following polynomial: tanh(x) ≈ 0.9568x -0.2107x 3 + 0.02354x 5 with maximum error below 0.01.

An implementation of this algorithm is then integrated in the tool ERAN (ETH Robustness Analyzer for Neural Networks software) [START_REF] Singh | Boosting robustness certification of neural networks[END_REF], which is a welldeveloped platform implemented in Python. It has been used for proving NN robustness against adversarial perturbations based on changes in pixel intensity and geometric transformations of images. What makes this tool also powerful is that it supports convolutional, feedforward, and residual networks with different types of activation functions. The tool deals with hyperbolic tangent activation functions via piecewise-linear approximations which, as mentioned earlier, can lead to imprecise results.Combining ERAN and the Bernstein techniques implemented in Sapo enables benefitting from the advantages of both, namely the accuracy of image approximations for activation functions of Sapo, and the optimized propagation of sets through NN layers of ERAN. Indeed, ERAN propagates the dependency between the neurons through the layers by grouping k activation functions [START_REF] Singh | Beyond the single neuron convex barrier for neural network certification[END_REF]. Note that propagating the full dependency between the neurons is prohibitively expensive, and the group size k is a user-defined parameter for fine-tuning the tradeoff between accuracy and computation cost, by retaining the dependency only between the outputs of the activation functions of the same group. ERAN over-approximates the inputs to these groups by sets of fixed form (also a way to reduce the geometric complexity of the sets generated through the propagation) which can be represented using the paral-lelotope bundles of Sapo. While ERAN uses a piecewise-linear approximation to deal with hyperbolic tangent activation functions, our Sapo-based algorithm uses the Bernstein expansion.

We illustrate now the advantage of this combination with two examples. One involves two neural network controllers for a robotic arm, and the other involves the neural network benchmarks provided by ERAN. The second example aims at showing the accuracy advantage of combining Sapo with ERAN by a comparison between this combination and ERAN alone. Their implementations are available at https://github.com/PippiaEleonora/SapoForNN.

Neural Network Controllers. The neural networks here are trained using data generated by a nominal discrete-time feedback PID controller for a robotic arm provided by MathWorks ® . The aim is to control the trajectory of this system while the value of the input reference signal r(•) randomly changes in [-0.5, 0.5] every 10 seconds. The control u k at time k is the result of the neural network on 10 inputs: the reference signals r k , r k-1 , r k-2 , r k-3 ; the outputs y k , y k-1 , y k-2 , y k-3 ; the controls u k-1 , u k-2 . Note that since PID controllers have memory and in order to capture it using a feedforward neural network, previous values of the variables should be fed as inputs to the network.

Two neural networks trained using different data sets are considered: one "bad" using data with diverging behaviour and one "good" with correct behaviour. We train two different neural nets with 2 hidden layers with 30 neurons each. We call the first neural network N N bad since it has a diverging behaviour as shown in Figure 17 for the constant reference with value -0.2. The second is called N N ok since it has a correct behaviour shown in Figure 18. We define a noise factor ϵ = |r| • p where p has three values of p = {0.01%, 0.1%, 0.5%}. We run our algorithm where the first 4 inputs r k , . . . , r k-3 take constant value -0.2, the next 4 inputs have the range [-0.2 -ϵ, -0.2 + ϵ] and the last two have the range u = 10y. We obtain the approximations in Table 1.

Here we use the above-described algorithm combining Sapo and ERAN to compute the set of reachable states under uncertain inputs for the two networks. The computation using our algorithm confirms the expected results, i.e., the "bad" network diverges, while the "good" one does not. For instance, when the perturbation p is at most 1% the "bad" network reaches the interval [-7.812, 1.049], while the "good" one remains in [-4.146 that the output interval of NN bad is 2.4 times the size of the output interval of NN ok , and this result is coherent with the behaviour of the two nets.

Sapo-ERAN versus ERAN. Our goal here is to evaluate if our algorithm using Sapo is more precise compared to ERAN which uses piecewise-linear approximations. In the experiments, we use a neural network with hyperbolic tangent activation functions for image recognition provided by ERAN (named mnist_tanh_3_50.tf). We describe in the following the results obtained for different sizes of the input box. The group size parameter k is 3, meaning that the activation functions of each layer are considered in groups of size 3.

We run ERAN and Sapo-ERAN over different input boxes (e.g., [-6, 6], [0, 1], and [START_REF] Althoff | Implementation of interval arithmetic in CORA 2016[END_REF][START_REF] Arino | A simple model for COVID-19[END_REF]) obtaining that Sapo-ERAN always returned a region included in the one returned by ERAN and, in particular, the rate between the volumes of the regions produced by the former and that of produced by the latter is always about 0.008. In Table 2 In the examples on neural networks we exploited the reachability functionalities of Sapo over systems without parameters. In the current implementation of Sapo the dynamics can be parametric with a linear dependence on the parameters. Such linear dependence allows the synthesis algorithms to rely on linear optimization subroutines. It would be very interesting if we could keep the linear functions defining the neural network parametric and synthesize the parameters that guarantee the expected behaviours. Unfortunately, the dependence on the parameters does not remain linear when the activation function is applied. We can easily synthesize linear parameters added at the level of the activation function, while it is more natural to have them in the linear part. The extension of Sapo with non linear optimization is currently under investigation, since it is of interest also in other applicative domains.

Conclusion and Future Work

We described algorithms for reachability and synthesis of parametric polynomial dynamical systems based on Bernstein coefficients that are implemented in the Sapo C++ library. We also presented a stand-alone application sapo that integrates these algorithms and avoids the need for problem recompilations. This stand-alone tool implements all the features offered by the Sapo library, presents them in a simplified interface, and offers some free extra bonus such as multi-threading. We are aware that libraries are usually more flexible than stand-alone applications and ease integrations in other tools and, because of this reason, sapo is meant to be an alternative way of using the Sapo library rather than a replacement for it.

The webSapo application further extends the audience of Sapo to any possible user on the internet. We presented two examples: the first from epidemiological models using a simple yet realistic COVID setting; the second from neural networks reachability analysis.

Sapo is a well-established tool, but we still consider to add features and expand its applicability domain. In the near future, it will support continuous dynamic laws and hybrid models. In particular, it will be able to over-approximate a flowpipe within a time interval by explicitly adding a time variable ranging in the interval itself. We also aim to support disturbances, control inputs, and we intend to integrate a model checker.

As far as the implementation is concerned, while Sapo already partially supports multi-threading and parallelizes the computation after parameter splits, we plan to extensively use this technology in linear algebra and symbolic operations too. Moreover, we would like to implement status dump, so to save the status of an incomplete analysis at wish and restart it afterward, and MPI multi-processing to distribute the computation on a cluster.

Another intriguing technical development consists of computing Bernstein coefficients and, more in general, each bundle's single evolution step by using GPU. In this case, the aimed efficiency boost will be related to convex components' complexity rather than their number. As far as the synthesis is concerned, the latter approach is symmetrical to the thread-based one: the more Sapo can refine a set of parameters without nullifying it, the less parameter splittings are performed.

Finally, we plan to explore the use of Sapo in other applicative contexts such as security verification [START_REF] Rossi | Model checking adaptive multilevel service compositions[END_REF], performances analysis [START_REF] Alzetta | Lumpingbased equivalences in markovian automata: Algorithms and applications to product-form analyses[END_REF], and blockchain [START_REF] Malakhov | Workload prediction in btc blockchain and application to the confirmation time estimation[END_REF].

(x 0

 0 Figure 1: The set Reach (0,5,1/5) [0,6]

Figure 2 :

 2 Figure 2: The evaluation of the STL formulas F [6,6] (0.015 ≤ I ≤ 0.03) on the trajectory ξ p x 0 described in Example 1. The points are labelled as x i instead of ξ px 0 (i) to simplify the plot. The formula does not hold because ξ p x 0 (6) does not satisfy the formula 0.015 ≤ I ≤ 0.03.

  formula G [0,3] I ≤ 0.025 since x0, x1, x2, and x3 satisfy the formula I ≤ 0.025. Since x0, x1, x2, and x3 do not satisfy R ≥ 0.01 and I ≤ 0.02 doesn't hold on x3, ξ p x 0 does not satisfy (I ≤ 0.02)U [0,3] (R ≥ 0.01).

Figure 3 :

 3 Figure 3: The evaluation of the two STL formulas G [0,3] I ≤ 0.025 and (I ≤ 0.02)U [0,3] (R ≥ 0.01) on the trajectory ξ px 0 described in Example 1. In order to simplify the plot, the points are labelled as x i meaning ξ p x 0 (i).

1 :

 1 function reach(X 0 , P, f , a, b)

Definition 9 ( 8 S 5 S + I ≤ 1 SIFigure 4 :

 98514 Figure 4: The projection on the plane R = 0 of the paralletope ⟨Λ, c⟩ described in Example 6. Each line in Λ and c defines a linear inequality on S, I, and R which corresponds to a halfspace in S × I × R. The paralletope is the intersection of all the half-spaces. The gray colored region depicted in this figure is the projection of ⟨Λ, c⟩.

2 Figure 6 :

 26 Figure 6: The projection on the plane R = 0 of the polytope ⟨D, c⟩ described by Example 5.⟨D, c⟩ can be decomposed in the two parallelotope P 1 and P 2 and, as a consequence, in the bundle B = {P 1 , P 2 }. The intersection of the bundle parallelotopes equals the polytope, i.e., I(B) = P 1 ∩ P 2 = ⟨D, c⟩

Example 12 .

 12 Let us consider once more the parametric polynomial described by Example 9 and assume that p 1 ranges in P = [0, 10]. Since p 1 is linear in b i (p 1 ), either b i (0) = min p1∈P (b i (p 1 )) and b i (10) = max p1∈P (b i (p 1 )) or b i (10) = min p1∈P (b i (p 1 )) and b i (0) = max p1∈P (b i (p 1 )). By Lemma 2, the polynomial π(x 1 , x 2 , p 1 ) ranges in the interval between min i∈I π {b i (0), b i (10)} and max i∈I π {b i (0), b i (10)} when x ∈ [0, 1] 2 . Moreover, because of the values of the Bernstein coefficients evaluated in Example 11, min i∈I π {b i (0), b i (10)} = b (0,1) (10) = -9/2 and max i∈I π {b i (0), b i (10)} = b (2,0) (10) = 23/6. Thus, the image of x ∈ [0, 1] 2 through the polynomial π(x 1 , x 2 , p 1 ) is contained by [-9/2, 23/6] (see Fig. 7).

z 1 zFigure 7 :

 17 Figure 7: The projection of π(x 1 , x 2 , p 1 ) = p 1 * x 2 1 /3 -p 1 * x 2 /2 + x 1 * x 2 /4 + 1/2 on the plane x 2 = 0 for x ∈ [0, 1]2 and p 1 ∈ [0, 10]. Under these conditions, π(x 1 , x 2 , p 1 ) lays in the interval [-9/2, 23/6] for all x ∈ [0, 1] 2 . See Example 9, 10, 11, and 12 for all the details.

Example 13 .

 13 Let us consider the SIR model of Example 1. Let S, I, and R range in [0, 1], while β ∈ [0.1, 0.2] and α ∈ [0.5, 0.6]. The Bernstein coefficients for the dynamic law of S are in {0, 1, 1 -β}, while those for I are in {0, 1α, 1+β -α}, and those for R are in {0, 1, α, 1+α}. Notice that some coefficients are repeated more than once for each function, e.g., 0. However, since we are only interested in minimum and maximum values, the repetitions are irrelevant. The minimum coefficient for S is 0, while the maximum is 1. The minimum value for I is 0 and the maximum is 0.7 (1 + β -α = 0.7 when β = 0.2 and α = 0.5). Finally, the minimum value for R is 0 and the maximum is 1.6 (1 + α = 1.6 when α = 0.6). So, the set reached after one step is a subset of [0, 1] × [0, 0.7] × [0,1.6]. It is worth noticing that the values for α that maximize I and R -i.e., 0.5 and 0.6, respectively-are different.If the parameters set is β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55], then the symbolic Bernstein coefficient sets are the same we computed above. Thus, the minima for S, I, and R and the maxima for S do not change, whereas the maxima for I and R are 0.65 (1 + β -α = 0.65 when β = 0.15 and α = 0.5) and 1.55 (1 + α = 1.55 when α = 0.55), respectively. So, when β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55], the parallelotope [0, 1] × [0, 0.65] × [0, 1.55] is an over-approximation of the reachable set in one step from the unit box.

6 :

 6 return X ′ = ⟨Λ, c ′ ⟩ 7: end function Example 14. Let us now compute another step from the set reached in Example 13. We are considering the set [0, 1]×[0, 0.7]×[0, 1.6] that can be represented by the parallelotope {x | Λx T ≤ c} where

  and those of ρ I are in {0, 1.6, 0.7α, 1.6 + 0.7α}. By minimizing and maximizing them over β ∈ [0.1, 0.2] and α ∈ [0.5, 0.6] we get the new intervals for S, I, and R, i.e., S ∈ [0, 1], I ∈ [0, 0.49], and R ∈ [0, 2.02].

55λ 3 )

 3 and the Bernstein coefficients sets for the dynamic laws of S, I, and R are {0, 1, 1 -0.65β}, {0, 0.65(1 -α), 0.65(1 + β -α)}, and {0, 1.55, 0.65α, 1.55 + 0.65α}, respectively. Thus, if β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55], (S, I, R) ranges in [0, 1] × [0, 0.4225] × [0, 1.9075].

Example 15 .( 1

 151 Let us consider the classical SIR model (see Example 1) and its evolution from x 0 = (0.99, 0.01, 0.0) when α ∈ {1/5, 4/5} and β = 1/2.Since α and β are parameters, their values should be fixed once for all, thus, ξ are the only valid trajectories from x 0 and, for instance, ξ (1/5,1/2) x0

Example 16 .

 16 Let us consider the SIR model described by Example 1 and let us synthesize a set of parameters satisfying the STL formula I ≤ 0.365 when S ∈ [0, 1], I ∈ [0, 0.7], R ∈ [0, 1.6], β ∈ [0.1, 0.2] and α ∈ [0.5, 0.6].

  365. Thus, the corresponding Bernstein coefficients are -0.365, 0.285 -0.65α, and 0.285 + 0.65β -0.65α and the algorithm returns the set such that α ∈ [0.5, 0.55], β ∈ [0.1, 0.15], -0.365 ≤ 0, 0.285 -0.65α ≤ 0, and 0.285 + 0.65β -0.65α ≤ 0, which are equivalent to α ≤ 0.55, β ≥ 0.1, and α -β ≥ 57/130. α -β ≥ 67/140 α ∈ [0.5, 0.6] β ∈ [0.1, 0.2] (a) The synthesized set when S ∈ [0, 1], I ∈ [0, 0.7], R ∈ [0, 1.6], β ∈ [0.1, 0.2] and α ∈ [α -β ≥ 57/130 α ∈ [0.5, 0.55] β ∈ [0.1, 0.15] (b) The synthesized set when S ∈ [0, 1], I ∈ [0, 0.65], R ∈ [0, 1.55], β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55].

Figure 9 :

 9 Figure 9: The synthesized parameter set for the SIR model presented in Example 1 satisfying the formula I ≤ 0.365. See Example 16 for the details.

Example 17 .

 17 Let us consider the SIR model described by Example 1 and let us synthesize a set of parameters satisfying the STL formula R ≤ 2 when S ∈ [0, 1], I ∈ [0, 0.7], R ∈ [0, 1.6], β ∈ [0.1, 0.2] and α ∈ [0.5, 0.6].

  65αλ 2 -2 (see the last part of Example 14). Thus, the corresponding Bernstein coefficients are -2, -0.45, 0.65α -2, and 0.65α -0.45 and the original parameter set, α ∈ [0.5, 0.55] and β ∈ [0.1, 0.15], is returned.

Example 18 .

 18 Let us synthesize the parameter set satisfying the STL formula I ≤ 0.365 ∧ R ≤ 2 on the SIR model described by Example 1 when S ∈ [0, 1], I ∈ [0, 0.7], R ∈ [0, 1.6], β ∈ [0.1, 0.2] and α ∈ [0.5, 0.6]. Under these conditions, the synthesized set for I ≤ 0.365 is such that α ≤ 0.6, β ≥ 0.1, and α -β ≥ 67/140 (see the first part of Example 16), while that for R ≤ 2 is such that α ∈ [0.5, 4/7] and β ∈ [0.1, 0.2] (see the first part of Example 17). Their intersection is empty because α -β ≥ 67/140 is not even satisfied by the maximum of α, i.e., 4/7, and the minimum for β, i.e., 0.1. If, instead, the initial conditions are S ∈ [0, 1], I ∈ [0, 0.65], R ∈ [0, 1.55], β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55], then, the synthesized set for I ≤ 0.365 is such that α ≤ 0.55, β ≥ 0.1, and α -β ≥ 57/130 (see the last part of Example 16), while that for R ≤ 2 is such that α ∈ [0.5, 0.55] and β ∈ [0.1, 0.15] (see the last part of Example 17). Their intersection equals the former. 3.3. Until Formulas The function UntilSynth(X, P, φ 1 U [a,b] φ 2 ) (Algorithm 4) refines the set P with respect to an until formula φ 1 U [a,b] φ 2 . It is structured in three main blocks, depending on the values a, b that define the interval of the until formula. The cases are the following: 1. a > 0 and b > 0: the interval is far from time 0; 2. a = 0 and b > 0: the interval starts at time 0 and ends somewhere else; Synthesized for I ≤ 0.365 Synthesized for R ≤ 2 (a) The synthesized set when S ∈ [0, 1], I ∈ [0, 0.7], R ∈ [0, 1.6], β ∈ [0.1, 0.2] and α ∈ [0.5, 0.6]. Synthesized for I ≤ 0.365 Synthesized for R ≤ 2 (b) The synthesized set when S ∈ [0, 1], I ∈ [0, 0.65], R ∈ [0, 1.55], β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55].

Figure 10 :

 10 Figure 10: The synthesized parameter set for the SIR model presented in Example 1 satisfying the formula I ≤ 0.365 ∧ R ≤ 2. See Example 18 for the details.

3 .Example 19 .

 319 a = 0 and b = 0: this is the base case of the recursive calls. It suffices to refine P with respect to φ 2 and return P φ2 . Let us consider once more the original SIR model as presented by Example 1 when the initial values for S, I, and R range in [0, 1], and the parameters α and β belong to [0.1, 0.2] and [0.5, 0.6], respectively. Let us synthesize the parameters α and β such that the STL formula

  Instead if the initial parameter set is α ∈ [0.5, 0.55] and β ∈ [0.1, 0.15], step 2 computes the set (S, I, R) ∈ [0, 1] × [0, 0.65] × [0, 1.55] as detailed in the last part of Example 13. Thus, step 3 synthesizes a parameter set satisfying the formula I ≤ 0.365 ∧ R ≤ 2 when the initial conditions are those just-mentioned and β ∈ [0.1, 0.15] and α ∈ [0.5, 0.55]. As seeing in the last part of Example 18, the algorithm returns the non-empty set α ≤ 0.55, β ≥ 0.1, and α -β ≥ 57/130.

  beta in [0.055, 0.1]; param mu in [0.00001, 0.001]; param gamma in [0.0027, 0.0055]; param alpha in [0.05, 0.07]; dynamic(s) = s -beta*s*i -mu*s + gamma*r; dynamic(i) = i + beta*s*i -alpha*i; dynamic(r) = r + mu*s -gamma*r + alpha*i; direction s in [0.2, 0.3]; direction i in [0.001, 0.1]; direction r in [0.7, 0.8];

Figure 12 ,

 12 Figure 12, which was produced by webSapo, depicts the number of infected individuals over time as computed by sapo.

Figure 12 :

 12 Figure 12: COVID-19 infected individuals over time. This image was produced by webSapo and then modified for editorial needs.

Figure 13 :

 13 Figure 13: Synthesized values of β and α with respect to G [0,100] (i < 0.1). This image was produced by webSapo and then modified for editorial needs.

Figure 14 :

 14 Figure14: COVID-19 infected individuals over time when using parameters satisfying G [0,100] (i < 0.1). This image was produced by webSapo and then modified for editorial needs.

Figure 15 :

 15 Figure 15: 2D projection of the synthesis over γ and µ. This image was produced by webSapo and then adapted to the editorial style of this journal.

Figure 17 :

 17 Figure 17: Comparison of the control for N N bad and a standard PID controller.

Figure 18 :

 18 Figure 18: Comparison of the control output for N N ok and a standard PID controller.

  |= φ 1 U I φ 2 iff ∃t ′ ∈ t + I s.t. ξ p x , t ′ |= φ 2 and for all t ′′ ∈ [t, t ′ ], ξ p x , t ′′ |= φ 1

	|= φ 1 and ξ p x , t |= φ 2 x , t |= φ 1 or ξ p x , t |= φ 1 ∨ φ 2 iff ξ p ξ p x , t |= φ 2 ξ p x , t We say that a trajectory ξ p x satisfies φ, denoted by ξ p x |= φ, if and only if ξ p x , 1 |= φ. Notice that the semantics is given with respect to t = 1, while
	usually t = 0 is considered. In other terms is like if we are considering only
	formulas starting with a next operator. It will become clear soon that this is a
	reasonable requirement for talking about synthesis.
	Example 4. Let us consider again the bounded trajectory described in Exam-
	ple 1. It satisfies the formula:
	G [0,3] (I ≤ 0.025),
	since the maximum value reached by I among ξ p x0 (0), ξ p x0 (1), ξ p x0 (2), and ξ p x0 (3)
	is 0.0215812. Instead, it does not satisfy the formula:

Table 1 :

 1 , 0.528]. We can see Interval approximation of NN bad and NN ok adding the noise factor ϵ.

	p	NN bad	NN ok
	0.01%	[-2.344, -1.684]	[-2.133, -1.862]
	0.1%	[-5.342, 1.313]	[-3.368, -0.627]
	0.5%	[-18.673, 14.600]	[-8.850, 4.858]

  , we report the outputs for the input box [-1, 1].

	Variable	ERAN	Sapo-ERAN
	n. 1	[-0.4212406, 0.5770473]	[0.0163105, 0.5726539]
	n. 2	[-0.3896767, 0.0377295]	[-0.3207002, -0.0517504]
	n. 3	[-0.3349207, 0.2658570]	[-0.1952389, 0.2020677]
	n. 4	[-0.3848480, 0.0126727]	[-0.3804094, -0.0099185]
	n. 5	[-0.2343316, 0.3914449]	[0.0450511, 0.3863766]
	n. 6	[0.0583854, 0.4344926]	[0.1958679, 0.3695503]
	n. 7	[0.1373136, 0.4805078]	[0.1809298, 0.3411388]
	n. 8	[0.2759886, 0.8369578]	[0.4734756, 0.8258935]
	n. 9	[0.0675532, 0.8206135]	[0.1222727, 0.5918586]
	n. 10	[-0.0174709, 0.5043359]	[-0.0100259, 0.4139967]

Table 2 :

 2 ERAN and Sapo-ERAN outputs for the input box [-1,1].

Acknowledgments Alberto Casagrande is partially supported by the INdAM-GNCS project "InSANE: Investigating Sparse Algorithms in the post von Neumann Era". Thao Dang is partially supported by the ANR-JST CREST project CyPhAI funded by the French National Research Agency ANR and Japan Science and Technology Agency JST, and the project Modylam (Plan Cancer, INSERM). Carla Piazza is partially supported by the Italian PRIN project "Noninterference and Reversibility Analysis in Private Blockchains (NiRvAna)" -20202FCJMH and by the INdAM/GNCS project "LESLIE: LogichE non-claSsiche per tooL Intelligenti ed Explainable".

We thank Alessandro Grisafi and Gianluca Ermacora for their contribution on the first version of the GUI for sapo.

The reachability and parameter synthesis algorithms described in the previous sections have been implemented in Sapo [START_REF] Dreossi | Sapo: Reachability computation and parameter synthesis of polynomial dynamical systems[END_REF][START_REF] Dreossi | Reachability Computation and Parameter Synthesis for Polynomial Dynamical Systems[END_REF]: an open source tool for the formal analysis of discrete-time polynomial dynamical systems freely available at https://github.com/dreossi/sapo.

Figure 11 sketches Sapo architecture. The Model and STL modules implement dynamical systems and temporal logic properties, respectively. The initial set of points, X 0 , and the parameter set, P , can be defined using the modules Bundles and LinearSystemSet. The main module Sapo performs the reachability and parameter synthesis computation. The BaseConverter module symbolically computes the Bernstein coefficients of polynomials. Thanks to the method described in [START_REF] Dreossi | Parameter synthesis for polynomial biological models[END_REF][START_REF] Dang | Parameter synthesis using parallelotopic enclosure and applications to epidemic models[END_REF], Sapo symbolically evaluates these coefficients only once and, then, fetches and numerically instantiates them during the computation.

r i s

Figure 16: The reachable set after synthesis. This image was produced by webSapo and then modified for editorial needs.