1. Introduction and Related Works
Many scientific applications in physical fields, such as mechanics, robotics, chemistry and electronics, require solving differential equations. This particular type of equation appears, for example, when the location is required; however, only the velocity and/or the acceleration are available when modeling a system. In the general case, these differential equations cannot be formally integrated—that is to say, closed form solutions are not available, and a numerical integration scheme is used to approximate the state of the system. Let us consider the most classical continuous differential equations, the Ordinary Differential Equations (ODEs):
One equation as above has a unique solution
if certain hypotheses occur:
f is Lipschitz, and the associated vector field is smooth. The most classical approach is to use a Runge–Kutta scheme—carefully chosen with respect to the problem, desired accuracy, etc.—to simulate the system behavior.
Historically, the first method for numerical solutions of differential equations was proposed by Euler in
Institutiones Calculi Integralis [
1]. His main idea is based on a simple principle: if a particle is located at
at time
and if its velocity at this time is known to be equal to
, then, at time
, the particle will be at approximately position
under the condition that
is sufficiently close to
(that is, after a short time), and thus the velocity does not change “too much” over
. Based on this principle, around the year 1900, C. Runge and M. W. Kutta developed a family of iterative methods, now called Runge–Kutta methods. While many such methods have been proposed since then, a unified formalism and a deep analysis was first proposed by John Butcher in the 1960s [
2].
The requirement of an unified and elegant formalism led to what we now call the Butcher series, or B-series (in respect to the seminal paper of John Butcher [
3]). Indeed, from the works of Cayley [
4] and Merson [
5], it has been shown that the expressions appearing in the derivatives of the solution of Equation (
1),
,
, are similar to rooted trees [
6]. Therefore, considering the formal series, such as
is natural and leads to express the solution of Equation (
1) under this form with the series coefficients
,
a(
) = 1/2,
a(
) = 1/3, etc. These series, with arbitrary coefficients, are named B-series [
7]. B-series are then the more general formalism to express computable solutions of ODEs. It is then a natural choice to define our set-based formalism.
Models of physical systems defined with the help of ODEs are, in general, simplified and, thus, approximate the actual behavior of such systems. It is, then, common to consider the uncertainties on the parameters to improve the representativeness of the model. These uncertainties can be defined by probabilistic methods or enclosed in sets. If the first technique is the most common—with Markov Chains [
8], Monte-Carlo approaches [
9], Graph theory [
10], etc.—the latter is more interesting in a verification procedure due to the guarantees provided (for safety, for example). Uncertainties, if bounded, can be enclosed in sets (for example, intervals), such that the models obtained are guaranteed to contain the actual ones [
11]. Other sets have been considered, such as polytopes, zonotopes (particular polytopes) and ellipsoids.
In the special case of ODEs, when parameters are given under the form of intervals, the computation of solutions is called reachability, i.e., the goal is then to compute the sets that enclose the solution of an ODE for all time and all possible values of the parameters.
Set-membership methods, associated with interval analysis, lead one to producing validated numerical integrations of ODEs. Validated numerical integration is an appealing approach to produce rigorous bounds of the solution of ODEs. Most of the rigorous techniques defined thus far, since Moore’s first algorithm based on Euler’s method [
12], are based on a Taylor series approach [
13]; see, for example, [
14,
15,
16,
17] and the references therein. Nevertheless, it is unlikely that only one kind of methods is adapted to all various classes of ODEs. Alternative recent works [
18,
19,
20,
21,
22,
23] deal with the adaptation of Runge–Kutta methods for ODEs and Differential-Algebraic Equations (DAEs) [
24].
Currently, important efforts are being made regarding the implementation of set-based methods for differential equation simulation. This leads to the releases of several tools. The most efficient, promising and innovative tools participate each year in the ARCH competition [
25,
26]: Ariadne, which adopts an approach with compositionality [
27]; CORA, which proposes various set representations [
28,
29]; DynIbex using Runge–Kutta methods and constraint programming [
23]; JuliaReach, which defines sets using Taylor models [
30]; Kaa based on parallelotope bundles [
31]; and KeYmaera X, a theorem prover able to consider differential Equations [
32].
In this article, we propose a unification of these validated numerical integration methods using B-series. B-series were introduced by E. Hairer and G. Wanner in [
7] after the seminal work of J. Butcher [
2]. As said before, they follow the results from Cayley [
4] and Merson [
5] on the expressions arising when studying the derivatives of the solution of an ODE. Since then, B-series have been intensely studied (see, e.g., [
6,
33] or [
34]); however, no papers refer to set-based B-series.
Outline.
We review the theory of trees, B-series and Runge–Kutta methods, in
Section 2. In
Section 3, we present the set-based framework proposed in this paper. In
Section 4, we show how validate set-based B-series using Runge–Kutta schemes. A small experiment is presented in
Section 5 to exhibit a result. In
Section 6, we summarize the main contributions of the paper.
2. Theories behind B-Series
As said in the introduction, B-series are strongly connected to rooted trees, and thus they are recalled in this section, followed by B-series and a Runge–Kutta overview.
2.1. Rooted Trees
A rooted tree [
35] is a graph with a distinguished node, called the root, in which every other node is connected to the root by a unique path. If the root of a tree
is removed and the nodes connected to the original node by a single bond are taken as new roots, the tree
breaks up into rooted trees
. Reversing this process, a new tree
can be constructed by joining the roots of the trees to a new common root. The number of nodes in a tree is denoted by
(order).
The number of equivalence classes of heap-orderings on a particular tree is denoted by
and can be computed using the Butcher’s formula:
where
denotes the symmetry group of
, and the tree factorial is defined recursively by
with the tree factorial of an isolated root defined to be 1
We denote, by , the symmetry of the tree .
2.2. Initial Value Problems
An initial value problem (IVP) consists on the merge of an ODE and a value, while an ODE has a solution (a function of time) only when a formal expression such that exists, an IVP has always a solution (if some hypotheses on the flow are respected).
An IVP can be expressed as follows:
Definition 1 (An IVP with ODE [
7]).
Starting from , we want to compute w.r.t. ODE ; . An exact solution is . As this expression makes the function (what we are looking for) appear, this solution can not be computed explicitly. This is why an approximation is often computed using numerical integration.
2.3. B-Series
Rooted trees allow one to build the formal series given in Equation (2), solution of the ODE
. The elementary differentials
are expressed regarding trees such that
with
,
,
, … } the set of trees and
the result of
k trees grafting.
A B-series [
36] is a formal series of the form
with
the elementary differentials as described above.
Lemma 1 (B-series and IVP [
36]).
The exact solution of an IVP is given by a B-series Lemma 2 (Taylor expansion).
The terms in a B-series correspond to the terms in the Taylor expansion of a map [37]. Operations on B-series:
Two important operations can be done on B-series, the composition and the substitution.
Composition law:
For , the composition is also a B-series, where and .
Substitution law:
For , the substitution is also a B-series, where and .
Definition of ordered subtrees S and partition P.
For the composition law, the notion of ordered subtrees (or admissible cuts) is required. An ordered subtree of a tree is a subset s of the set of all vertices connected and containing the root of . All the subtrees of are gathered in the set . In the previous paragraph, is then a collection of rooted trees (a forest) obtained by removing the subtree s to the tree , and is the rooted tree obtained after the cuts.
For the substitution law, the notion of partitions is needed. A partition p of a tree is a subset of the edges of the tree. All the partitions are gathered in the set . In the previous paragraph, is the forest remaining when edges of p are removed in , and is the tree built from all trees from (also called the skeleton).
2.4. Runge–Kutta Methods
With the previous notations and definitions, a Runge–Kutta scheme can be defined as a B-series.
Lemma 3 (Runge–Kutta and B-series [
36]).
A Runge–Kutta scheme approximates the solution of ODE (1) with a B-series:if the Runge–Kutta scheme is defined as followswith coefficients given in a Butcher tableau [2] of the form In terms of the form of the matrix A, consisting of the coefficients , a Runge–Kutta method can be:
Explicit—for example, as in the classical Runge–Kutta method of order 4 given in
Figure 1a. It means that the computation of the intermediate steps
only depends on the previous steps
for
.
Diagonally implicit—for example, as in the diagonally implicit fourth-order method given in
Figure 1b. In this case, the computation of an intermediate step
involves the value
itself; therefore, non-linear systems in
must be solved (with a Newton algorithm, for example). A method is
singly diagonally implicit if the coefficients on the diagonal are all equal.
Fully implicit—for example, the Runge–Kutta fourth-order method with a Lobatto quadrature formula given in
Figure 1c. In this last case, the computation of intermediate steps involves the solution of a non-linear system of equations in all the values
for
.
The order of a Runge–Kutta method is
p if and only if the
local truncation error, in other words, the distance between the exact solution
and the numerical solution
is such that:
In term of B-series, it means that a Runge–Kutta method is of order
p if
which means that
.
The main idea behind Runge–Kutta methods.
As a B-series (or a Taylor expansion) is an infinite and noncomputable, sum of terms that gives the solution of an ODE. By cutting it at a given order, one can obtain a sufficient approximation. Runge–Kutta methods provide these approximations at a given order if the previous condition on coefficients holds without derivative computation.
3. B-Series and Sets
In this section, a formalism of set-based B-series is presented. First, different kinds of ODEs with initial conditions and parameters defined as sets are discussed. Indeed, considering uncertainties bounded in sets, an IVP with the ODE (
1) can be:
a parametrized ODE , with and ;
a non precise initial value problem , with ; or
a non precise initial value problem with a parametrized ODE , with and .
Case 1 and 3 have to be discussed. Case 1 can be seen as the differential inclusion ; however, the nature of the parameter is important as it could lead to errors in some contributions for reachability computations.
In terms of the literature, one paper studied the computation of rigorous reachable set in the particular case of differential inclusion [
38]. The authors proposed some conditions on parameters and a novel algorithm based on [
13].
3.1. Discussion on the Nature of Parameters
Two kind of parameters could appear in a model of physical systems: (i) a physical parameter, unknown but constant or (ii) a varying parameter, bounded but with an unknown derivative.
Physical Parameters
Considering a physical parameter in a model, such as a length, mass, etc, marred by uncertainties during the modeling process, the actual value is, of course, not perfectly known but is unique, constant (even if the mass can vary with respect to the fuel consumption in some applications, in these cases, the mass without fuel consideration is the parameter that interests us in this discussion) and bounded. This leads to the fact that the parameter derivative is equal to zero, and then the ODE can be extended as follows:
with
(and then
). The problem is then similar to a non precise initial value problem with parameter added to states. If we split the parameter set
such that
, then the reachable set corresponding to
is the set-union of the reachable sets associated with
and
due to the monotonicity of sets. It leads to the fact that algorithms based on bisection will be correct. Derivatives of
p will not appear in the high order terms in the B-series; thus, the order of the B-series is respected.
3.2. Not Fixed Parameters
The other kind of parameter is time varying—for example, an input of a controlled system, a noise or some perturbations. These parameters (or variables) are unknown but bounded. Their derivatives are also unknown and different than zero. This leads to the following extended ODE:
with
. This ODE cannot be integrated because it does not fulfill the requirements (non smooth and non Lipschitz). Moreover, if we split the parameter set
such that
, then the reachable set corresponding to
is not the set-union of the reachable sets associated with
and
, because, at each instant, all the values for
p are possible, i.e., taken in
or in
. Therefore, bisection-based approaches must be avoided in this case (some papers have proposed bisection-based algorithms—for example, for control synthesis, which leads to mistakes in reachable sets and then in safety verification).
Derivatives of p will appear in the high order terms in the B-series; thus, the order of the B-series is not respected. In this objective, parameters have to be considered as constant during an integration step (that can be small regarding a control frequency, for example).
To conclude this discussion, cases 1, 2 and 3 can be seen as differential inclusions, under the condition that parameters are considered as piecewise constant. This strong condition is finally due to a lack of knowledge in the model that we attempt to simulate. In order to define a more precise model, one can add parameters in the states and define these derivatives as parameters (at least one of its derivatives has to be piecewise constant).
3.3. Set-Based B-Series
We propose to define set-based B-series to express the solution of differential inclusions (with the assumption of piecewise constant parameters as discussed above), such as
From the B-series definition, the following can be proposed:
Proposition 1 (Set-based B-series for differential inclusions).
A set-based B-series (with ) is the exact reachable set of the differential inclusion (6) at instant h. Set-based B-series being defined as set of punctual B-series, we state what follows:
Proposition 2 (B-series define an element of set-based B-series). With the B-series , with and , the following inclusion holds due to the monotonic inclusion of sets: .
Considering the pointwise union (or collecting), the following holds
This latter is equivalent to the principle behind the Monte-Carlo approach with uniform distributions for parameters and initial condition and an infinite number of simulations. The main advantage of set-based approaches is that one reachability computation is equivalent to an infinite number of pointwise computations.
Proposition 3 (Composition law is applicable to set-based B-series).
Proof. A set-based B-series describes the reachable set at instant h such that , and then choosing leads to the composition law. □
Proposition 4 (Substitution law is applicable to set-based B-series).
Proof. A set-based B-series describes a set inclusion that can be seen as a differential inclusion, which leads to the substitution law. □
4. Validated B-Series
In this section, a generalization of validated integration (reachability analysis or guaranteed simulation) is proposed using set-based B-series. In these methods, computation must be guaranteed by considering numerical errors from the integration scheme (i.e., the truncation of the B-series) and from the floating point arithmetic (this point is discussed in
Section 4.5).
4.1. Validated Truncated B-Series
B-series are truncated in order to be computable (with Runge–Kutta methods or Taylor approximations). First of all, by considering the truncation error, one can validate a finite B-series at order
p with the following formula:
with
the unknown Lagrangian remainder point,
and the functions
and
This latter formula comes directly from the Lagrangian remainder of Taylor expansion [
39]. The term
is called the local truncation error (LTE). The term
can be computed formally or with automatic differentiation for Taylor methods, or with Runge–Kutta methods (see [
40] for details on Runge–Kutta coefficients for validated integration) while the term
can be computed formally, with automatic differentiation or rooted trees (more details in [
23,
41]).
Remark 1. The addition between a B-series and a set-based B-series is a Minkowski sum. This depends on the representation chosen for the sets (see Section 4.5), and this addition is never computed in practice as it is preferable to directly add two set-based B-series. 4.2. With ODEs
The main issue in the above formula is to find the Lagrangian term to compute the LTE. By the help of set-based B-series, an enclosure is computable. Indeed, the exact solution of an ODE, such as (
1), can be enclosed in a validated B-series with the help of the following formula:
The second term on the right part is a set-based B-series, evaluated on the set
, which is proven by construction (see
Section 4.4) to contain the Lagrangian remainder point, i.e.,
.
Runge–Kutta methods:
Let us consider a Runge–Kutta scheme at order
p and defined by the Butcher tableau
, then a validated B-series can be constructed:
with
, and
or
The expression
comes from the LTE definition, obtained with the difference of the two Lagrangian remainders at order
of two B-series (the one from the exact solution and the one from the numerical approximation), see [
2] for more details. An important remark is that, if
, then the LTE of the Runge–Kutta scheme will be smaller than the Taylor’s one, i.e., the approximation is then better with the Runge–Kutta scheme than with the order equivalent Taylor series.
Remark 2. With a zero order method, i.e., , we obtain the Euler method.
4.3. With Differential Inclusions
One can compute the reachable set of a differential inclusion, such as (
6), with a validated truncated set-based B-series with the following formula:
Runge–Kutta methods can be used to compute these B-series, as described above.
4.4. Picard–Lindelöf
To validate a B-series, a Lagrangian remainder has to be computed in an unknown point, see Equation (
7). In Equation (
12), the proposed solution consists of enclosing this unknown point in a set
. In order to obtain this latter, the well-known Picard–Lindelöf theorem can be used. Furthermore, called the Cauchy–Lipschitz theorem, or existence and uniqueness theorem, it gives a set of conditions, under which, an initial value problem accepts a unique solution. Details on this theorem and proof can be found in several books, such as in [
42].
A Picard–Lindelöf operator has been defined, already based on sets and B-series:
Definition 2 (Picard–Lindelöf operator with B-series [
3]).
If there exists a set such thatthen, due to Brouwer theorem (also proposed/studied by Picard, Poincaré and Hadamard) and Picard–Lindelöf theorem, the initial value problem has a unique solution in the set . As a consequence, this proves that (existence and uniqueness simultaneously). Using this operator in a simple inflate and contract algorithm, one can obtain the set
, such that
, and can then replace
by
to compute the LTE in Equation (
12). The infinite B-series
can be replaced by a validated truncated B-series such that
Picard–Lindelöf operator for Runge–Kutta methods:
Let us consider a Runge–Kutta scheme at order
p and defined by the Butcher tableau
, then a validated B-series can be constructed:
In addition, the latter can be itself used as a Picard–Lindelöf operator:
Definition 3 (Picard–Lindelöf operator with validated Runge–Kutta)
If there exists a set such thatthen, the initial value problem has a unique solution in the set . With differential inclusions:
Due to the monotonic inclusion, the previous results are directly applied to set-based B-series. Therefore, one can compute the reachable set of a differential inclusion, such as (
6) with a validated truncated set-based B-series with the following formula:
Furthermore, the same result is obtained with validated truncated set-based B-series using Runge–Kutta (by replacing by and ).
Remark 3. With a zero order method, i.e., , we obtain the validated Euler method, also called the rectangle rule, as proposed by Moore [12]. 4.5. Remarks on Set Representation
The abstraction used to represent sets in a computer must have several properties in order to implement the presented techniques. First, as said in the section on validated B-series, the chosen abstraction must handle numerical errors due to floating numbers. Second, it has to come with an arithmetic, such as a Minkowski sum, as well as all the guaranteed arithmetic operators, even the nonlinear ones (for non linear ODES, e.g., with the cosine or square root). Finally, set operators, at least the inclusion test, are needed for the Picard–Lindelöf operator.
Some abstractions are well-known and mainly used for reachability, such as intervals [
12], zonotopes [
43] or Taylor models [
44]. Other approaches, less famous but efficient in certain cases/applications, have been tested, such as ellipsoids [
45], parallelotopes [
46], polytopes [
47], etc.
A trick to compute with polytopes:
As said above, polytopes can be used in set-based B-series, even if no associated arithmetic is available. Indeed, a polytope can be described by the intersection of a finite number of zonotopes [
47]. Let us define an initial condition
enclosed by a polytope
. If
, then the following holds:
4.6. A Basic Algorithm to Compute Reachability
An algorithm can be built gathering previous definitions, able to compute reachability of differential inclusions. An overview of this algorithm, based on the well-known Lohner two-step method [
13], is given in the following. Many steps are not described; for more details on the algorithm and its implementation, refer to [
23].
- 0.
Starting with an initial condition , a differential inclusion , a stepsize h, a guess .
- 1.
Initialize reachable set .
- 2.
While , inflate .
- 3.
Let .
- 4.
Compute reachable set .
- 5.
Go back to step 2.
This algorithm can be strongly improved [
23]—for example, with some tricks for the guess
, for the stepsize controller, with additive contraction on
, the techniques used to compute the B-series terms, set operations, etc.