Sapo: Reachability Computation and Parameter
Synthesis of Polynomial Dynamical Systems
Tommaso Dreossi

arXiv:1607.02200v1 [cs.SY] 8 Jul 2016

University of California, Berkeley
Email: tommasodreossi@berkeley.edu

Abstract—Sapo is a C++ tool for the formal analysis of polynomial dynamical systems. Its main features are: 1) Reachability
computation, i.e., the calculation of the set of states reachable
from a set of initial conditions, and 2) Parameter synthesis,
i.e., the refinement of a set of parameters so that the system
satisfies a given specification. Sapo can represent reachable
sets as unions of boxes, parallelotopes, or parallelotope bundles
(symbolic representation of polytopes). Sets of parameters are
represented with polytopes while specifications are formalized as
Signal Temporal Logic (STL) formulas.

I. I NTRODUCTION
Formal verification involves the strict and exhaustive study
of systems using mathematically based techniques. The development and application of formal methods is motivated by
the fact that the formal study of a model implies the reliability
and robustness of the abstracted system.
Two important problems emerging from the formal analysis
of dynamical systems are the reachability computation and the
parameter synthesis problems. In the first case, for a given
set of initial conditions, it is asked to compute the set of
states reachable by the dynamical system. In the second case,
an initial set of conditions together with a set of parameters
and a specification are provided, and it is asked to refine
the set of parameters in such a way that the system satisfies
the specification. These two problems often recur in formal
analysis. For instance, the reachability computation can be
used to determine whether a systems reaches some undesired
states or remains within a particular subset of the state space,
while synthesis of parameters can be used to design dynamical
systems and tune their parameters so that to fit experimental
observations or satisfy requirements.
Despite the numerous tools available for the analysis of
linear dynamical systems (i.e., systems whose dynamics are
linear functions), not many tools for the study of nonlinear
systems are available. The main difficulties arising from
nonlinearity concern the transformation of sets with respect
to nonlinear functions. This operation, in general, does not
preserve nice properties of the sets (e.g., convexity), and
thus approximation techniques, that often suffer in terms of
scalability, are necessary.
In this work we present Sapo [1], a C++ tool that targets discrete-time polynomial dynamical systems (possibly
equipped with parameters) and that deals with both the
reachability computation and parameter synthesis problems.
Sapo can be used either to generate flowpipes that overapproximate reachable sets over bounded time horizons, or

to refine sets of parameters in such a way that the system
satisfies a Signal Temporal Logic (STL) [17] specification.
The tool allows the construction of flowpipes using boxes (i.e.,
hyperrectangles), parallelotopes (i.e., n-dimensional parallelograms), or parallelotope bundles (i.e., sets of parallelotopes
whose intersections symbolically represent polytopes), and the
refinement of parameter sets represented by polytopes. Sapo’s
underlying algorithms exploit some properties of Bernstein
coefficients of polynomials. Intuitively, the Bernstein coefficients can be used to bound a polynomial over the unit box
domain (i.e., [0, 1]n ). The core idea on which Sapo relies is the
transformation of unit boxes into generic sets. With this trick,
Sapo exploits Bernstein coefficients to maximize polynomials
and determine sets that over-approximate the trajectories of
the system. A similar approach is also used to reduce the
parameter synthesis problem into linear programs, where the
linear systems, constructed through parameterized Bernstein
coefficients, represent sets of valid parameters. The algorithms
implemented in Sapo are based on the theoretical results
exposed in [8], [11], [6], [7], [12], [10].
II. R EACHABILITY AND PARAMETER S YNTHESIS
Sapo deals with discrete-time parametric polynomial dynamical systems described by difference equations of the form:
xk+1 = f (xk , p)

(1)

where f : Rn ×Rm → Rn is a polynomial linear in p, xk ∈ Rn
is the state of the system at time k ∈ N, and p ∈ Rm are the
parameters. Given an initial condition x0 ∈ Rn and parameters
p ∈ Rm , the trajectory x0 , x1 , x2 , . . . describing the evolution
the system, can be obtained by iterating the function (1). The
problems addressed by Sapo are:
•

•

Reachability computation: Given a set of initial conditions X0 ⊂ Rn , a set of parameters P ⊂ Rm , and
a time instant T ∈ N, determine a sequence of tight
sets X0 , X1 , X2 , . . . XT , called flowpipe, such that all the
trajectories of length T with initial conditions in X0 and
parameters in P are included in the constructed flowpipe;
Parameter synthesis: Given a set of initial conditions
X0 ⊂ Rn , a set of parameters P ⊂ Rm , and a
specification ϕ, determine the largest set Pϕ ⊆ P such
that all the trajectories with initial conditions in X0 and
parameters in Pϕ satisfy ϕ.

III. M AIN F EATURES
A. Bernstein Coefficients
At the core of Sapo for both reachability computation
and parameter synthesis there are Bernstein coefficients of
polynomials. Intuitively, Bernstein coefficients can be used
to represent a polynomial in Bernstein form, i.e., as the
linear combination of Bernstein basis and coefficients. Given
a polynomial π : Rn → R:
X
ai xi
(2)
π(x) =
i∈I π

where i = (i1 , . . . , in ) ∈ Nn is a multi-index, xi = xi11 . . . xinn
is a monomial, ai ∈ R is a coefficient, and I π ⊂ Nn is the
multi-index set of π, the i-th Bernstein coefficient is:

X ji
(3)
bi =
 aj
d
j≤i

j


i
where i ≤ j if ik ≤ jk for k =
1, . . . , n,
j is the product of


the binomial coefficients ji11 . . . jinn , and d is the smallest
multi-index such that i ≤ d for all i ∈ I π .
One of the interesting properties of Bernstein coefficients
is the range enclosure property [4] stating that minπ bi ≤
i∈I

π(x) ≤ max
bi for x ∈ [0, 1]n . This implies that the maximum
i∈I π
Bernstein coefficient is an upper bound of the maximum of
π over the unit box domain. Thus, to optimize a polynomial,
instead of using nonlinear/nonconvex optimization techniques,
one can compute the Bernstein coefficients and extract their
maximum. The main drawback of this approach is that the
enclosure property holds only on the unit box domain.
In [10] we extended the range enclosure property to parametric polynomials of the form π : Rn × Rm → R,
showing that minπ min bi (p) ≤ π(x, p) ≤ max
max bi (p) for
π
i∈I

p∈P

i∈I

p∈P

all x ∈ [0, 1]n and parameters p ∈ P . This property, with a
slight adaptation of the treated polynomial, is the key element
of Sapo’s reachability and parameter synthesis algorithms.
B. Reachability Computation
The flowpipe X0 , X1 , . . . , XT that includes the reachable
set can be obtained as a sequence of set image transformations
Xk+1 = f (Xk , P ). Unluckily, this operation becomes difficult
when f is nonlinear, since many properties of sets (e.g.,
convexity) can be lost. A common method to circumvent this
problem consists in over-approximating the image of sets with
simpler objects such as polytopes. A polytope X ⊂ Rn is a
convex set that can be represented by the solutions of a system
Dx ≤ c, where D ∈ Rm×n and c ∈ Rm . The matrix D
and vector c are called template and offset, and the polytope
generated by D and c is denoted by hD, ci.
Given a template D, the offset c such that the set hD, ci
over-approximates the set f (Xk , P ) can be determined by
solving optimization problems of the form:
ci ≥

max

x∈Xk ,p∈P

Di f (x, p)

(4)

for i = 1, . . . , m. An upper-bound of the maximum of
Di f (x, p) over the unit box can be found computing its
Bernstein coefficients and maximizing them over the set P .
However, in order to apply this technique on generic domains,
we need a slight adjustment of the treated polynomial.
In [11] we introduced a method that combines Bernstein
coefficients and parameter synthesis on generic box domains.
Intuitively, for a given box Xk , we can compute a map v :
[0, 1]n → Xk that transforms the unit box into Xk . It is easy
to see that f (Xk , P ) = f (v([0, 1]n ), P ), which implies that
an upper bound of Di f (x, p), with x ∈ Xk and p ∈ P , can be
established determining the maximum Bernstein coefficient of
Di f (v(x), p). Repeating this operation for all the directions
of the template D, we obtain a new offset c that leads to an
over-approximating box hD, ci ⊇ f (Xk , P ).
Changing the map v we can alter the shape of the overapproximation set and define reachability algorithms based
on sets different from boxes. For instance, in [6] we shown
how the map v can be defined to transform unit boxes into
parallelotopes. In doing so, we defined a parallelotope-based
set image approximation technique that is more flexible in the
choice of the initial set and preciser than the box-based one.
The transformation of a unit box into a generic polytope
is in general difficult. However, in [12] we introduced a
new way or representing polytopes as finite intersections of
polytopes. We called these sets parallelotope bundles. Once
that a polytope is decomposed in a collection of parallelotopes,
we can reason separately on each of its parallelotopes and,
using the Bernstein coefficients, we can determine a new
parallelotope bundle that over-approximates the image of the
starting polytope. This method sensibly increases the accuracy
of the over-approximation flowpipes at the cost of a higher
number of optimizations.
In summary, for reachability computation, Sapo supports
the construction of flowpipes based on boxes, parallelotopes,
and parallelotope bundles. All these objects are represented
by Sapo as both solutions of linear systems (constraint representation) and affine transformations of unit boxes (generator
representation). Sapo uses either of these representations depending on the operation to perform on the set.
C. Parameter Synthesis
Given a set of initial conditions X0 , a set of parameters
P , and a specification ϕ, we want to determine the largest
set Pϕ ⊆ P such that the reachable set starting in f (X0 , Pϕ )
satisfies ϕ.
Sapo allows the user to formalize the specifications in terms
of STL formulas [17] in positive normal form, i.e., formulas
generated by the following grammar:
ϕ := σ | ϕ ∧ ϕ | ϕ ∨ ϕ | ϕUI ϕ
n

(5)

with σ = g(x1 , . . . , xn ) ≤ 0, where g : R → R, and I ⊂ N
is an interval.
Similarly to the reachability computation, also the parameter synthesis algorithm implemented in Sapo exploits the
Bernstein coefficients. Given the sets Xk , P and a predicate

σ = g(x1 , . . . , xn ) ≤ 0, the set P is a valid parameter
set if g(f (x, p)) ≤ 0 holds for all x ∈ Xk and p ∈ P .
This constraint can be verified by determining a map v :
[0, 1]n → Xk , computing the Bernstein coefficients of the
function g(f (v(x), p)), and checking that bi (p) ≤ 0 for all
bi (p) ∈ I g◦f ◦v and p ∈ P . In particular, a refinement Pσ of
the set P with respect to the predicate σ can be obtained by
adding the linear constraints bi (p) ≤ 0 to the linear system
representing P . The solutions of the new linear system are
the valid parameters. This condition can be verified solving a
single linear program.
Reasoning by structural induction on the given specification,
it is possible to reduce the parameter synthesis problem into
a collection of refinements over predicates. The partial results
are then combined accordingly with the treated formula. The
refinement of the until formula requires also the computation
of the evolution of the system, task that can be achieved by
the reachability methods previously introduced. The parameter
synthesis algorithm, its correctness, and its complexity are
presented in detail in [7], [10].
In summary, for the parameter synthesis, Sapo receives in
input a set of initial conditions representable with a box or
a parallelotope, a polytopic set of parameters specified as a
linear system, an STL formula, and it produces a collection of
polytopes representing the set of parameter set under which
the system satisfies the specification.
IV. S TRUCTURE AND U SAGE
A. Tool Architecture
Sapo is implemented in C++. It relies on the libraries
GiNaC1 for symbolically manipulating polynomials and on
GLPK2 for solving linear programs.
Figure 1 summarizes the structure of the tool. The user
can specify a dynamical system and a possible specification
using the Model and STL modules, respectively. The set
of initial conditions and possible parameters can be specified using the modules Bundle and LinearSystemSet,
respectively (a single box or parallelotope can be seen as
a bundle with a single template). With these elements, the
user can analyze the dynamical system invoking the main
module SapoCore. This block implements the reachability
and parameter synthesis algorithms and returns the computed
results. An important module with which SapoCore interacts
is the BaseConverter. This block computes Bernstein
coefficients of polynomial functions heavily exploited by our
algorithms. The BaseConverter symbolically computes
the Bernstein coefficients of a polynomial only once. The
symbolic coefficients are stored in a data structure from
which they are fetched and numerically instantiated whenever
needed. This trick allows the tool to save computations and
sensibly speeds up the analysis. Bernstein coefficients are
computed by default using our improved matrix method [11].
To visualize the results, Sapo gives the possibility to generate a
1
2

http://www.ginac.de
https://www.gnu.org/software/glpk/

(a) Box-based

(b) Bundle-based

Fig. 2: Reachable sets of SIR model (300 steps)

Matlab script that displays 2/3-dimensional sets or projections
of higher dimensional sets.

B. Experimental Evaluation
As a demonstration, we apply Sapo to the SIR epidemic
model, a well-known 3d nonlinear system that describes the
progress of a disease in a population. The model considers
three groups of individuals: s the susceptible healthy, i the
infected, and r the removed from the system (e.g., the recovered ones). Two parameters regulate the system’s evolution: β
the contraction rate and γ, where 1/γ is the mean infective
period.
We used Sapo to compute the bounded time reachable set
of the SIR model with different configurations. At first, with
parameters β = 0.34, γ = 0.05, and set of initial conditions
s0 ∈ [0.79, 0.80], i0 ∈ [0.19, 0.20], and r0 ∈ [0.00, 0.00], we
computed in 0.12s the reachable set for 300 steps using a
single box (see Figure 2a). Then, by adding two surrounding
parallelotopes to the set of initial conditions, we obtained in
2.83s a bundle-based flowpipe that better over-approximates
the reachable set (see Figure 2b). Notice how, for the same
initial set, the bundle-based flowpipe is more accurate than the
box-based one.
Sapo can also be used to synthesize sets of valid parameters.
For instance, considering the set of initial conditions s0 ∈
[0.79, 0.80], i0 ∈ [0.19, 0.20], and r0 = [0.00, 0.00], and the
initial set of parameters β ∈ [0.18, 0.20] and γ ∈ [0.05, 0.06],
we can ask Sapo to refine the parameter set so that the
STL formula G[50,100] (i ≤ 0.44) holds. Figure 3a shows the
original parameter set (in white) and the synthesized one (in
gray) computed in 0.83s, while Figure 3b depicts the evolution
of the system under the synthesized parameter set.
We evaluated both Sapo’s reachability and parameter synthesis algorithms on several dynamical systems [11], [6], [7],
[10], [12]. We computed the reachable sets of biological
models, such as SIR (3d) or a honeybees nest choice model
(5d), or of a quadrotor drone (17d) (see Appendix A), and
synthesized parameters of epidemic systems, such as SIR (3d,
2params), influenza (4d, 2params), Ebola (5d, 4params), and
SARS (6d, 4params) (see Appendix B).

Sapo
Sapo core

STL formula

Model

Base
converter

Bundle

Linear
system set

Parallelotope

Linear
system

Fig. 1: Architecture of Sapo.

scalability. Finally, a comparison with the available tools for
nonlinear reachability is on our schedule.
R EFERENCES

(a) Synthesized parameters

(b) Constrained evolution

Fig. 3: Parameter synthesis of SIR model (G[50,100] (i ≤ 0.44))

V. C ONCLUSION
A. Related Work
In the last two decades numerous tools have been developed for the reachability analysis of linear dynamical system and hybrid automata. Some examples are SpaceEx [14],
HyTech [15], or d/dt [3]. Unlike for the linear case, not
many tools for the nonlinear analysis are available. Some
tools, such as Breach [9], S-TaLiRo [2], or C2E2 [13] address
nonlinearity using a finite number of simulations. Tools that
directly deal with nonlinear reachability are dReach [16] and
Flow*[5]. The latter is probably the closest tool to Sapo in
terms of both problem formulation and provided results. The
main difference between Flow* and Sapo is that the first
produces flowpipes grouping collections of Taylor models,
while the second produces a collection of polytopes computed
using Bernestein coefficients. Another difference is that Flow*
considers hybrid automata, while Sapo, at the moment, focuses
only on dynamical systems.
B. Future Developments
Sapo is still under development. As future work, we will
extend the reachability computation to polynomial hybrid
automata, where invariants and guards need to be considered.
An interesting aspect of operations on parallelotope bundles
is that they can be easily parallelized. It could be interesting
to implement a parallel version of the tool and investigate its

[1] Sapo. Available at http://www.eecs.berkeley.edu/∼tommasodreossi/sapo/.
[2] Y. Annpureddy, C. Liu, G. Fainekos, and S. Sankaranarayanan. S-taliro:
A tool for temporal logic falsification for hybrid systems. In Tools and
Algorithms for the Construction and Analysis of Systems, TACAS, pages
254–257, 2011.
[3] E. Asarin, T. Dang, and O. Maler. The d/dt tool for verification of hybrid
systems. In Computer Aided Verification, CAV, pages 365–370, 2002.
[4] G.T. Cargo and O. Shisha. The Bernstein form of a polynomial. In
Journal of Research, pages 79–81, 1966.
[5] X. Chen, E. Ábrahám, and S. Sankaranarayanan. Flow*: An analyzer for
non-linear hybrid systems. In Computer Aided Verification, CAV, pages
258–263, 2013.
[6] T. Dang, T. Dreossi, and C. Piazza. Parameter synthesis using parallelotopic enclosure and applications to epidemic models. In Hybrid Systems
and Biology, HSB, pages 67–82, 2014.
[7] T. Dang, T. Dreossi, and C. Piazza. Parameter synthesis through temporal
logic specifications. In International Synmposium on Formal Methods,
FM, pages 213–230, 2015.
[8] T. Dang and R. Testylier. Reachability analysis for polynomial dynamical
systems using the Bernstein expansion. Reliable Computing, 17(2):128–
152, 2012.
[9] A. Donzé. Breach, a toolbox for verification and parameter synthesis of
hybrid systems. In Computer Aided Verification, CAV, pages 167–170,
2010.
[10] T. Dreossi. Reachability Computation and Parameter Synthesis for
Polynomial Dynamical Systems. PhD thesis, University of Udine and
Joseph Fourier University, 2016.
[11] T. Dreossi and T. Dang. Parameter synthesis for polynomial biological
models. In Hybrid Systems: Computation and Control, HSCC, pages
233–242, 2014.
[12] T. Dreossi, T. Dang, and C. Piazza. Parallelotope bundles for polynomial
reachability. In Hybrid Systems: Computation and Control, HSCC, pages
297–306, 2016.
[13] P. S. Duggirala, M. Potok, S. Mitra, and M. Viswanathan. C2E2: a tool
for verifying annotated hybrid systems. In Hybrid Systems: Computation
and Control, HSCC, pages 307–308, 2015.
[14] G. Frehse, C. Le Guernic, A. Donzé, S. Cottont, R. Ray, O. Lebeltel, R.
Ripado, A. Girard, T. Dang, and O. Maler. SpaceEx: Scalable verification
of hybrid systems SpaceEx: Scalable verification of hybrid systems In
Computer Aided Verification, CAV, pages 379–395, 2011.
[15] T. A. Henzinger, P. Ho, and H. Wong-Toi. HyTech: A model checker for
hybrid systems. In Computer Aided Verification, CAV, pages 460–463,
1997.
[16] S. Kong, S. Gao, W. Chen, and E. Clarke. dreach: δ-reachability analysis
for hybrid systems. In Tools and Algorithms for the Construction and
Analysis of Systems, TACAS, pages 200–205, 2015.

Model
SIR [5]

Vars
2

Nest choice [1]

5

Quadcopter [2]

17

Dir/Temps
4/2
3/1
6/4
5/1
7/3
17/1
18/2

Steps
25
60
60
1500
1500
300
300

Time
0.24
0.12
2.83
26.90
81.27
17.74
39.07

TABLE I: Reachability computation benchmark
Model
SIR [5]

Vars
3

Params
2

Influenza [4]
Ebola [3]

4
5

2
4

Spec
G[10,30] (i ≤ 0.68)
G[20,50] (i ≤ 0.50)
G[0,50] (i ≤ 0.43)
(i ≤ 200)U[6,10 (q ≤ 20)
(q ≤ 40)U[10,15] (i ≤ 270)
(q ≤ 50)U[5,15] (e > 100 ∨ q > 25)

TABLE II: Parameter synthesis benchmark (scalability in
system’s dimension)

[17] O. Maler and D. Nickovic. Monitoring temporal properties of continuous
signals. In Formal Techniques, Modelling and Analysis of Timed and
Fault-Tolerant Systems, pages 152–166, 2004.

A PPENDIX A
R EACHABILITY
Table I reports the evaluation of Sapo’s reachability computation algorithm. The tool has been evaluated on three
nonlinear systems: the SIR epidemic model [5], a honeybees
site choice model [1] describing the mechanism adopted by
a swarm of honeybees to choose one among two different
nest-sites, and a quadrotor drone model [2] characterizing
both the drone and its controller dynamics. Each system has
been tested on different combinations of number of directions
and templates. The table reports the dimension of the system
(i.e., the number of variables), the number of directions and
templates used for the parallelotope bundle, the total number
of reachability steps, and the computation times expressed in
seconds.
A PPENDIX B
PARAMETER S YNTHESIS
The scalability of the parameter synthesis algorithm has
been studied in terms of system dimension (Section B-A) and
specification’s length (Section B-B).
A. Benchmark 1
Table II reports the evaluation considering models of increasing dimension. We considered three nonlinear dynamical
systems: the SIR model [5], a simplification of the influenza
model presented in [4], and a variation of the model describing
the Ebola outbreak in Congo 1995 and Uganda 2000 exposed
in [3]. For each model we considered different specifications.
The table reports the dimension of the system, the numbers of
synthesized parameters, the considered STL specification, and
the computation time expressed in seconds.

Time
0.10
0.65
6.27
0.10
0.14
0.14

a
5
5
5
5
5
5
15
20
30
50
100

b
15
20
30
50
100
125
20
30
50
100
200

φ1
0.20 (11)
0.35 (16)
0.82 (26)
2.63 (46)
0.29 (6)
0.64 (11)
1.88 (21)
13.00 (51)
-

φ2
0.15 (7)
0.24 (11)
0.55 (21)
1.65 (41)
10.95 (91)
0.22 (4)
0.45 (11)
1.27 (21)
6.89 (51)
-

φ3
0.14 (4)
0.21 (4)
0.36 (4)
0.80 (4)
1.69 (4)
1.69 (4)
0.19 (0)
0.35 (0)
0.78 (0)
1.72 (0)
1.66 (0)

(a) Increasing until interval.
N
1
2
3
4
5
6
7
8
9
10
15
16
17

φN
1
0.11 (5)
0.26 (9)
0.48 (13)
0.74 (17)
1.10 (21)
1.50 (25)
1.97 (29)
2.59 (33)
3.23 (37)
4.98 (42)
10.33 (61)
11.75 (65)
-

φN
2
0.14 (2)
0.36 (7)
0.69 (12)
1.09 (17)
1.61 (22)
2.28 (27)
3.05 (32)
3.97 (37)
5.06 (42)
6.70 (47)
11.80 (62)
15.98 (67)
-

φN
3
0.13 (4)
0.29 (3)
0.50 (3)
0.70 (3)
1.01 (3)
1.20 (3)
1.65 (3)
1.81 (3)
2.16 (3)
2.56 (3)
4.85 (3)
5.43 (3)
5.97 (3)

(b) Nesting until.

TABLE III: Parameter synthesis benchmark (scalability in
specification’s size). Times and computed polytopes per refinement. φ1 ≡ (i ≤ 200)U[a,b] (q ≤ 20), φ2 ≡ (q ≤ 40)U[a,b] (i ≤
270), φ3 ≡ (q ≤ 50)U[a,b] (e > 100 ∨ q > 25).

B. Benchmark 2
For the second evaluation of the parameter synthesis algorithm we considered the six dimensional Ebola model [3]
and we synthesize four of its controllable parameters over
different specifications. Let φ1 ≡ (i ≤ 200)U[a,b] (q ≤ 20),
φ2 ≡ (q ≤ 40)U[a,b] (i ≤ 270), and φ3 ≡ (q ≤ 50)U[a,b] (e >
100 ∨ q > 25) be three STL specifications. Table IIIa reports
the evaluations when stretching the time intervals [a, b] of
the specifications φ1 , φ2 , and φ3 . As second evaluation, we
nested several until on the (most critical) right hand side.
For instance, the double nesting of φ1 is φN
1 ≡ (I(t) ≤
200)U[6, 10]((I(t) ≤ 200)U[6, 10](Q(t) ≤ 20)) with N = 2.
Table IIIb reports the running times and the number of
obtained parameter sets.
R EFERENCES
[1] N. F. Britton, N. R. Franks, S.C. Pratt, and T. D.Seeley. Deciding on
a new home: how do honeybees agree? In Royal Society of London B:
Biological Sciences, pages 1383–1388, 2002.
[2] A. E. Carrilho da Cunha. Benchmark: Quadrotor Attitude Control. In
Applied Verification for Continuous and Hybrid Systems, ARCH, 2015.
[3] G. Chowell, N. W. Hengartner, C. Castillo-Chavez, P. Fenimore, and J.
M. Hyman. The basic reproductive number of Ebola and the effects of
public health measures: the cases of Congo and Uganda. In Journal of
Theoretical Biology, 2004.

[4] P. A. González-Parra, S. Lee, L. Velazquez, and C. Castillo-Chavez. A
note on the use of optimal control on a discrete time model of influenza
dynamics. In Mathematical Biosciences and Engineering, pages 183–197,
2011.
[5] W. Kermack, O. William, and A. G. McKendrick. A contribution to
the mathematical theory of epidemics. In Royal Society of London A:
Mathematical, Physical and Engineering Sciences, pages 700–721, 1927.

