SENSITIVITY ANALYSIS OF UNCERTAIN DYNAMIC SYSTEMS USING SET-VALUED INTEGRATION∗

We present an extension of set-valued integration to enable efficient sensitivity analysis of parameter-dependent ordinary differential equation (ODE) systems, using both the forward and adjoint methods. The focus is on continuous-time set-valued integration, whereby auxiliary ODE systems are derived whose solutions describe high-order inclusions of the parametric trajectories in the form of polynomial models. The forward and adjoint auxiliary ODE systems treat the parameterization error of the original differential variables as a time-varying uncertainty, and propagate the sensitivity bounds forward and backward in time, respectively. This construction enables building on the sensitivity analysis capabilities of state-of-the-art solvers, such as CVODES in the SUNDIALS suite. Several numerical case studies are presented to assess the performance and accuracy of these set-valued sensitivity integrators.

1. Introduction. Set-valued integration methods for uncertain or parameterdependent dynamic systems find applications in many research areas, including reachability and invariance analysis for control systems, robust optimization and control, set-membership estimation, and global optimization; see, e.g., [5,9,18,19,22,26,36,48,50]. One possible classification for such methods distinguishes continuous-time and discretized approaches. The former derives an auxiliary differential system describing pointwise-in-time enclosures of the set of trajectories, which can be solved using a standard numerical ordinary differential equation (ODE) solver; the latter constructs parameterized tubes using time-series expansion and determines finite intervals where these tubes yield a guaranteed enclosure of the set of trajectories, following a predictor-validation scheme. Both methods can propagate a variety of set parameterizations, which describe (i) convex enclosures of the reachable set, in the form of interval boxes, ellipsoids, zonotopes, or polytopes [2,10,16,19,21,33,43,45]; (ii) nonconvex enclosures of the reachable set, in the form of Taylor or Chebyshev models [1,4,11,20,27,49]; and (iii) pairs of linear or convex/concave bounds on the variations of each dynamic state [41,42,44,46,47].
Methods and software for sensitivity analysis of dynamic models are also well developed [7,14,17,29]. They have reached a stage of maturity where both forward and adjoint sensitivities can be computed reliably and efficiently. These methods typ-

A3015
ically use automatic differentiation [15] and numerical integration techniques tailored to the special structure of the sensitivity system, such as staggered or simultaneous corrector methods for forward-in-time sensitivity integration, and state interpolation and check-pointing for backward-in-time adjoint integration. While forward sensitivity analysis is best suited to compute the sensitivities of a possibly large number of functions with respect to a small number of parameters, adjoint sensitivity analysis is best suited to the complementary situation of finding the sensitivity of a small number of functions with respect to a large number of parameters. Despite their widespread use in such areas as uncertainty quantification, experimental design, parameter estimation, optimization, and control, very little work has been devoted to set-valued integration for sensitivity analysis so far. In [39], the parametric ODEs and their first-order sensitivities are bounded as one combined set of ODEs. But more efficient staggered or simultaneous corrector methods have not been investigated in this context to date. Regarding adjoint sensitivity analysis, an approach to propagating bounds on the adjoint trajectories has been considered recently in [37] using interval analysis. However, more advanced bounding approaches based on ellipsoidal calculus are yet to be developed and analyzed.
Our main objective in this paper is an extension of set-valued integrators to enable sensitivity analysis for a class of switched nonlinear dynamic systems, using both the forward and adjoint methods. The main focus is on continuous-time setvalued integration, as it enables building on existing numerical solvers to efficiently propagate the sensitivity bounds, as well as high-order inclusion techniques in the form of polynomial models, which mitigate the wrapping effect and can stabilize the reachable set enclosures under certain conditions [20,49]. Of course, other set-valued integration approaches could be extended in a similar manner.
The outline of this paper is as follows. Section 2 gives a formal statement of the problem and reviews the corresponding forward and adjoint sensitivity systems. Section 3 presents an extension of continuous-time set-valued integration that relies on polynomial model inclusion techniques in order to handle both parametric and time-varying uncertainty. Two variants of this approach, one propagating interval remainders and the other ellipsoidal remainders, are described. Section 4 specializes this generic bounding capability to the forward and adjoint sensitivity systems, whereby the state parameterization error is treated as a time-varying uncertainty. A numerical implementation of the forward and adjoint sensitivity-bounding systems is also discussed, which enables building on top of existing numerical solvers. This implementation comes in the form of a C++ class as part of the library CRONOS, which relies on the numerical integrator CVODES in the SUNDIALS suite [17] and is available from https://github.com/omega-icl/cronos. Several numerical case studies are presented in section 5 to assess the performance and accuracy of the proposed set-valued sensitivity integrators. Finally, section 6 concludes the paper.
1.1. Notation. The set of compact subsets of R n is denoted by K n and the subset of convex subsets of K n by K n C . The Minkowski sum of two compact sets W, Z ∈ K n is given by The set of n-dimensional interval vectors is denoted by IR n [31]. The midpoint and radius of an interval vector P := [p, p] ∈ IR n are defined as mid(P ) := 1 2 (p + p) and rad(P ) := 1 2 (p − p), respectively, and we denote by diag rad(P ) ∈ R n×n the diagonal matrix whose elements are the components of rad(P ).
A set-valued function Ψ : IR n → K m is said to have convergence order k ≥ 1 on a set Z ∈ K n if there exists a constant C < ∞ such that (1) rad(Ψ(W )) ≤ C rad(W ) k for all W ∈ IR n with W ⊂ Z and sufficiently small rad(W ). The set of n-dimensional positive semidefinite symmetric matrices is denoted by S n + . An ellipsoid with center c ∈ R n and shape matrix Q ∈ S n + is denoted by or simply E(Q) for those ellipsoids centered at the origin. Given two n-dimensional symmetric matrices A and B, the property A − B ∈ S n + induces a partial order as A B. Moreover, tr B (A) denotes the trace of A scaled by B, given by (2) tr with B := max i=1...n B ii , and some relative tolerance RTOL > 0. Likewise, rad B (P ) denotes the scaled radius of an interval vector P ∈ IR n , given componentwise by with α κ ∈ R m ; and R q f ∈ K m is such that ∀p ∈ P , f (p) − P q f (p) ∈ R q f . Here, B κ : R n → R m is the basis function with multi-index κ ∈ N n in the selected polynomial representation. Taylor models [6,28,34] are a particular kind of polynomial model where the multivariate polynomial matches the qth-order Taylor expansion of f at some pointp ∈ P , provided that f is sufficiently smooth. In Chebyshev models [11,38], the multivariate polynomial approximates the qth-order Chebyshev expansion of f on P .
L n 2 denotes the set of square-integrable vector-valued functions of dimension n, and W n 1,2 denotes the set of n-dimensional weakly differentiable functions whose weak derivative is square-integrable. Unless otherwise stated, the domain of integrability is taken as R, and integration is understood with respect to the time variable. The abbreviation a.e. is used to indicate that a property holds almost everywhere.
2. Problem statement and background. We consider parametric dynamic systems in the form of switched nonlinear ODEs with fixed-time transitions t 0 < t 1 < · · · < t N ,ẋ with x(t 0 , p) = h 0 (p) (5) and

A3017
Here, the state trajectory x : R × P → R nx is regarded as a function of the uncertain parameter vector p ∈ P ∈ K np , and the notations t − k and t + k are used to make it clear when the state-and later other quantities too-is evaluated immediately before or after a transition time t k .
The problem addressed in this paper is the computation of bounds on the gradient of state-dependent functions expressed in the Bolza form The focus is on continuous-time set-propagation techniques that construct an auxiliary ODE bounding system based on both the forward and adjoint sensitivity methods, which are reviewed in subsections 2.1 and 2.2 below. Note that the quadrature terms in (7) could be reformulated as extra differential equations by augmenting (4). The main reason for not applying such a reformulation up front is that many numerical solvers can handle quadrature variables more efficiently in treating them separately from the other state variables and thereby reduce the size of the Jacobian matrix. For simplicity, we assume throughout the paper that the functions m k , q k , f k , and h k are at least continuously differentiable in all their arguments. Moreover, we make the assumption that a solution to the parametric initial value problem (4)- (6) exists for every parameter value p ∈ P .
2.1. Forward sensitivity method. Upon applying a chain rule of differentiation, the following set of expressions is obtained for the derivatives of the statedependent function (7) [29]: A set of differential equations and boundary conditions defining the state-sensitivity variable s j := ∂x ∂pj is also obtained by applying a chain rule of differentiation to the initial value problem (4)-(6), Since both the state and sensitivity problems are specified at the initial time t 0 , solving the sensitivity ODEs (9) simultaneously with the state ODEs (4), forward in time, removes the need to store the state trajectories. However, the overall number of ODEs with this approach scales as O(n x n p ), which may be computationally demanding when n p is large. Both staggered and simultaneous corrector methods [14,25] have been proposed to mitigate the computational burden for large-scale problems. In contrast, the computational effort remains largely unaffected by the number of state-dependent functions.
2.2. Adjoint sensitivity method. Following the derivations in [8,40], an alternative set of expressions for the derivatives (8) is Here, the adjoint (or costate) variables λ i [13] are defined by the following differential equations and boundary conditions: for each i ∈ {1, . . . , n F }, and A complication with the adjoint approach stems from the fact that the boundary conditions of the state and adjoint systems are specified at different points, with the former specified at initial time t 0 and the latter at terminal time t N . Therefore, these two systems may no longer be integrated simultaneously, as previously with the forward sensitivity approach. Nonetheless, since the state equations do not depend on the adjoint variables, a popular approach that preserves numerical stability proceeds in two steps: (i) integrate the state system (4)-(6) forward in time until t N , and store the computed state values obtained at (sufficiently many) grid points; then (ii) integrate the adjoint problem (13)-(15) backward in time until t 0 , by interpolating the state trajectories between grid points as necessary. Efficient implementations of this approach are available [7,17], which use accurate interpolation schemes alongside check-pointing to control the storage space requirement.
Another important observation about the adjoint problem (13)- (15) is that there is no explicit specification of the parameters p. This implies that the adjoint λ i associated to the function F i can be used to compute the gradient of that function with respect to any of the parameters. In other words, the overall number of ODEs required to compute the derivatives ∂F ∂p with this approach scales as O(n x n F ), making it more efficient than the forward sensitivity approach whenever n p n F .
3. Continuous-time set-valued integration with both parametric and time-varying uncertainty. A precursor to this work is the ongoing research into reachability analysis for uncertain nonlinear ODEs. This section presents an extension of existing continuous-time techniques based on polynomial model inclusions [9,49] in order to encompass time-varying uncertainty in addition to the parametric variations. This extension will prove necessary later in deriving auxiliary bounding systems for the forward and adjoint sensitivity problems (9)-(11) and (13)- (15), due to the presence of state parameterization errors that act as a time-varying uncertainty.
We consider the following generic uncertain dynamic system formulation: with y(t 0 , p, w) = η 0 (p) (17) and where w : R → R nw denotes the time-varying uncertainty, assumed to be in the set W := {w ∈ L nw 2 | ∀t ∈ [t 0 , t N ], w(t) ∈ W (t)} for a given set-valued function W : R → K nw . The objective is to describe a time-varying enclosure Y (t) ⊇ Y (t) of the reachable set of (16)- (18), defined as For the particular case of qth-order polynomial models, the reachable set enclosures come in the form (20) Y (t) := P q y (t, p) p ∈ P ⊕ R q y (t) . Based on the understanding that the polynomial models all share the same order q, we shall drop the superscript q in polynomial model notations in the rest of this section for simplicity.
Following [49], one natural approach to constructing the multivariate polynomial function P y : R × P → R ny is to propagate the O(n y (n p ) q ) monomial coefficients of the polynomial approximant through ODEs. After aggregating all these coefficients back into polynomial form, the following conditions hold: with P y (t 0 , p) = P η0 (p) (22) and where w ∈ W is a given uncertainty reference trajectory. The corresponding parameterization error r y : R × P × W → R ny , given by r y (t, p, w) := y(t, p, w) − P y (t, p), satisfies the initial value probleṁ r y (t, p, w) = ϕ k (t, P y (t, p) + r y (t, p, w), p, w(t)) − P ϕ k (t, P y (t, p), p, w(t)) (24) =: r ϕ k (t, r y (t, p, w), p, w(t)) a.e. t ∈ [t + k−1 , t − k ] , k = 1 . . . N, with r y (0, p, w) = η 0 (p) − P η0 (p) (25) and for every uncertainty scenario (p, w) ∈ P × W. Therefore, a reachable set enclosure in the form of (20) calls for an enclosure of the parameterization error, Two methods for computing such enclosures are described next, a simple method based on the classical theory of differential inequalities (subsection 3.1) and a secondorder method based on ellipsoidal calculus (subsection 3.2). The propagation of bounds for the quadrature variables is also discussed as a special case (subsection 3.3).

Interval remainder bounds.
We consider the case of a time-varying un- . It is convenient to rescale this time-varying uncertainty within [−1, 1] nw , by applying the linear transformation w(t) := c w (t) + ∆ w (t) (t), with c w (t) := mid(W (t)) and ∆ w (t) := diag rad(W (t)). Moreover, a natural choice for the uncertainty reference trajectory w is the midpoint c w . The rescaled problem (24)-(26) reads aṡ and The classical theory of differential inequalities [24,51] can be used to construct an interval enclosure of the parameterization error in the form R y (t) := [r y (t), r y (t)]. A sufficient condition is the availability of piecewise Lipschitz-continuous sub-and superfunctions r y (t), r y (t) : R → R ny satisfying [37] a.e. t for each i = 1 . . . n y , along with the initial conditions Overall, this approach entails the solution of O(n y (n p ) q ) auxiliary ODEs to describe a reachable set enclosure Y (t) for (16)-(18) on P × W. In particular, an interval box enclosure is obtained with q = 0.

Ellipsoidal remainder bounds.
The use of axis-aligned interval boxes for bounding the parameterization error can result in a large overestimation for some problems-even for asymptotically stable dynamic systems-due to the well-known wrapping effect [20]. This subsection presents a new approach to enclosing the parameterization error within an ellipsoid, as a means to mitigate wrapping. This is essentially a unification between (i) existing algorithms to bound ODEs with time-varying uncertainty in [19,23] and (ii) high-order inclusion algorithms based on polynomial models with ellipsoidal remainders to bound parametric ODEs in [49].
We consider the case of a time-varying uncertainty tube W with ellipsoidal cross-  (38) and Our approach to propagating ellipsoidal bounds R y (t) := E(Q y (t)) for the parameterization error entails a decomposition of (37) and (39) into linear and nonlinear parts asṙ Here, the matrices A ϕ k (t), B ϕ k (t), A η k , and B η k may result from a linearization of the functions ϕ k and η k , e.g., using Taylor or Chebyshev expansion. The nonlinear term N ϕ k (resp., N η k ) is so chosen that (37) and (40) (resp., (39) and (41)) are equivalent. Moreover, we shall assume for now that nonlinearity bounders Ω ϕ k (t), Ω η0 , Ω η k ∈ IR ny can be constructed such that A computational procedure for automating such bounders is outlined in subsection 4.2.
A direct consequence of Proposition 5 in [49] and Theorem 4.1 in [19] is that an ellipsoidal enclosure E(Q y (t)) of the parameterization error may be obtained for any piecewise Lipschitz-continuous function Q y : along with the initial conditions and the transition conditions for some functions ν : R → R ++ and κ : Particular choices for all these functions and scalars are discussed later in subsection 4.2.
Overall, the ellipsoidal bounding approach requires the solution of O(n y (n p ) q + (n y ) 2 ) auxiliary ODEs to describe a reachable set enclosure for (16)-(18) on P × W.
3.3. Quadrature bounds. The focus thus far has been on propagating the reachable set for systems of nonlinear ODEs under uncertainty, without consideration of quadrature terms. Clearly, any quadrature term of the form with initial condition q(t 0 , p, w) = 0 and transition conditions q(t + k , p, w) = q(t − k , p, w), k = 1 . . . N − 1. Hence, a straightforward approach to dealing with such quadratures proceeds by appending one extra ODE (48) for each quadrature term to the original problem (16) and then bounding the extended system using a similar method as previously. In an ellipsoidal remainder bounding approach, handling the quadrature equations (49) in the same way as the other ODEs (16) could artificially inflate the state parameterization error bound R y (t) since it does not recognize that the righthand side of the extended ODE system (16), (49) is independent of the quadrature variable q(t) itself. Therefore, this simple approach is not used here.
Instead, each quadrature term can be bounded separately by using the available polynomial model inclusion (20) of the reachable set and treating the state parameterization error in R y (t) as a time-varying uncertainty. Specifically, a polynomial model inclusion in the form {P q (t, p) | p ∈ P } ⊕ [r q (t), r q (t)] may be obtained for each quadrature variable, where the polynomial part relies on the propagation of O((n p ) q ) extra monomial coefficients through ODEs in order to satisfy with P q (t 0 , p) = 0 and P q (t + k , p) = P q (t − k , p), k = 1 . . . N − 1; and the quadrature parameterization error bounds are any pair of Lipschitz-continuous functions along with the initial and transition conditions Interestingly, this approach will also allow exploiting certain features of state-of-theart numerical ODE solvers to handle quadrature terms more efficiently.
4. Application to forward and adjoint sensitivity bounding. Enclosures for the derivatives of the state-dependent function (7), in the form can be computed by combining the sensitivity analysis techniques in subsections 2.1 and 2.2 with the generic continuous-time set-propagation techniques in section 3. Table 1 summarizes how the variables and functions in the generic uncertain dynamic system (16)- (18) specialize to the case of the state system (4)-(6), the forward sensitivity system (9)- (11), and the adjoint sensitivity system (13)- (15). Based on this specialization, the derivative bounds (55) can be obtained by propagating either polynomial models with interval remainders based on the sufficient conditions in (21)- (23) and (31)- (36), or polynomial models with ellipsoidal remainders based on the sufficient conditions in (21)- (23) and (45)- (47). At this point, it is worth reiterating that the polynomial approximant P q x (t, p) and the corresponding parameterization error R q x (t) in a reachable set inclusion (20) act as a parametric uncertainty and a time-varying disturbance, respectively, in both the forward and adjoint sensitivity bounding systems. Table 1 Correspondence between the variables and functions in the generic uncertain dynamic system (16)- (18) and those in the state system (4)-(6), the forward sensitivity system (9)- (11), and the adjoint sensitivity system (13)- (15).
In the forward sensitivity approach, polynomial models of both the state variables and their sensitivities may be jointly propagated forward in time since their respective initial conditions are all specified at t = t 0 . That is, the required state polynomial model inclusions (P q x (t, p), R q x (t)) are available in order to compute the sensitivity polynomial model inclusions (P q sj (t, p), R q sj (t)) at each time t. Algorithmic advances for reducing the computational burden in classical forward sensitivity analysis, such as staggered or simultaneous corrector methods [14,25], could be used in the scope of set-valued sensitivity integration. In doing so, however, it is important to define analytical sensitivity right-hand sides, and not let a solver create or evaluate these equations using finite differences or automatic differentiation which would result in incorrect sensitivity bounds. Another possible caveat is that the state and forward sensitivity bounding systems no longer share the same Jacobian matrices, due to different sources of overestimation in the state and forward sensitivity parameterization error inclusions. Those implementations of the staggered corrector method, whereby the Jacobian matrix from the state system is reused subsequently for the sensitivity system for improved efficiency, could therefore result in incorrect sensitivity bounds. In contrast, simultaneous corrector methods are safe to use even though the Jacobian matrices of the state and sensitivity systems could be different since the right-hand side of the sensitivity equations needs reevaluating at each time step.
(ii) In the adjoint sensitivity approach, joint forward propagation in time of the state and adjoint polynomial models is no longer possible since the adjoint initial values are specified at t = t N . In order to enable the propagation of adjoint polynomial model inclusions (P q λi (t, p), R q λi (t)) backward in time, one may use polynomial model inclusions (P q x (t, p), R q x (t)) interpolated from a number of mesh points stored during a preliminary state bounding step. This two-step approach is essentially identical to the one used in classical adjoint sensitivity analysis [8], and it may therefore benefit from algorithmic advances thereof, such as check-pointing in order to control the forward storage requirement.
In both approaches, the main practical difficulty in using state-of-the-art numerical solvers for propagating the sensitivity bounds is constructing and evaluating the right-hand side, initial value, and transition value functions of the auxiliary sensitivity bounding systems. Subsections 4.1 and 4.2 provide details for such constructions. We assume that the functions f k , h k , m k , and q k in the problem definition (4)-(7) are factorable-namely defined by a finite number of bivariate sums, products, and univariate compositions [30]-and sufficiently smooth, so that automatic differentiation and set-valued computer algebra can be applied.

Computing polynomial models with interval remainder bounds.
A high-order inclusion of the reachable set (19), in the form of a qth-order polynomial model with interval reminder (P q y (t, ·), [r q y (t), r q y (t)]), can be computed based on the sufficient conditions in (21)-(23) and (31)-(36), with equality signs in all of the inequalities. At a given time t ∈ [t + k−1 , t − k ], k = 1 . . . N , and for each component i = 1, . . . , n y , we compute a qth-order polynomial model on P of the composite function and R w (t) ∈ IR nw an interval enclosure of the time-varying uncertainty W at t. By construction, the coefficients of the multivariate polynomial P q ϕ k,i provide the time derivatives of the polynomial coefficients in (21), whereas the lower bound (min R q ϕ k,i ) provides the right-hand side of (31). The right-hand side of (32) is constructed likewise. At the initial and transition times t k , k = 0 . . . N , the polynomial parts in (22) and (23) and the interval remainder bounds in (33)-(36) are directly obtained from polynomial models of the initial and transition value functions η k on P .

Computing polynomial models with ellipsoidal remainder bounds.
A high-order inclusion of the reachable set (19), in the form of a qth-order polynomial model with ellipsoidal remainder (P q y (t, ·), E(Q q y (t))), can be computed based on the sufficient conditions in (21)- (23) and (45)-(47), with equality signs in all of the inequalities.
Of the alternatives for constructing the right-hand sides of (21) and (45) at a given time t ∈ [t + k−1 , t − k ], k = 1 . . . N , we use the following two-step procedure, based on a modification of the mean-value approach described in [49]: (i) Compute a qth-order polynomial model with an interval remainder of the composite function ϕ k (t, P q y (t, ·), ·, c w (t)) on P , such that ∀p ∈ P , ϕ k (t, P q y (t, p), p, c w (t)) ∈ with α ϕ k ,γ ∈ R ny and R 0 ϕ k ∈ IR ny . (ii) Compute q th-order polynomial models with interval remainders and q ≤ q of the Jacobian matrices ∂ϕ k ∂y (t, P q y (t, ·) + R q y (t), ·, R w (t)) and ∂ϕ k ∂w (t, P q y (t, ·) + R q y (t), ·, R w (t)) on P , such that ∂ϕ k ∂y (t, P q y (t, p) + r, p, w) ∈ C ∂yϕ k ⊕ R ∂yϕ k and ∂ϕ k ∂w (t, P q y (t, p) + r, p, w) ∈ {C ∂wϕ k } ⊕ R ∂wϕ k , with C ∂yϕ k ∈ R ny×ny , C ∂wϕ k ∈ R ny×nw , R ∂yϕ k ∈ IR ny×ny , and R ∂wϕ k ∈ IR ny×nw . Here, (P q y (t, ·), R q y (t)) denotes a q th-order polynomial model with an interval re- mainder of the reachable set (19), as derived from (P q y (t, ·), E(Q q y (t))) by appending all the terms of order greater than q to the remainder term; and the interval vector R w (t) ∈ IR nw associated with the time-varying uncertainty W is such that R w (t) ⊇ E(Q w (t)).
By construction, the time derivatives of the polynomial coefficients in (21) can be taken as the coefficients α ϕ k ,γ in (56); and the matrices A ϕ k (t) and B ϕ k (t) in (40) and the nonlinearity bounder Ω ϕ k (t) in (42) can be taken as . One simple choice for the ellipsoid parameterization trajectories κ(t) and ν(t) in (45) that minimize tr Q q tr (Q q y (t)) + for some small regularization > 0. Where the magnitude of the parameterization error term R q y (t) varies by multiple orders in different directions, it is often advantageous to minimize the scaled trace tr Q q for some finite tolerance RTOL 1 in (2) and (3). The same two-step procedure may be applied for constructing the right-hand sides of the transition conditions (23) and (47), merely by substituting ϕ k for η k . Here, one particular choice of the ellipsoid parameterization scalars α k , β k , and γ k in (47) that minimize tr(Q q y (t + k )), for each k = 1 . . . N , is . Like previously, a scaled version may be used instead, e.g., with the scalars α k , β k , and γ k minimizing tr Q q ). Finally, the right-hand sides of the initial conditions (22) and (46) may be readily obtained from a qth-order polynomial model (P q η0 , R q η0 ) of η 0 on P , by setting Ω η0 := R q η0 and choosing the ellipsoid parameterization scalars γ 0 in (46) as γ 0 := rad(Ω η0 ) tr (diag rad(Ω η k )) + .

Convergence and stability properties of the sensitivity bounds.
An important property in comparing the tightness of the derivative enclosures ∂ p F in (55) is how quickly they converge to the actual set of derivatives when the size of the parameter set P vanishes. This convergence rate is dictated by how quickly the reachable set inclusions for the state and forward/adjoint sensitivity systems themselves converge to the exact reachable sets.
We start by recalling that, for any factorable and (q +1)-times continuously differentiable function Ψ : R n → R m , the remainder bound R q Ψ in a qth-order polynomial model (P q Ψ , R q Ψ ) of Ψ on Z ∈ IR n has convergence order (at least) q + 1 on Z in the sense of (1) when it is constructed using either Taylor or Chebyshev model arithmetic [6,38]. Moreover, it has been established in [49] that the error parameterization bounds R q x (t) in polynomial model inclusions of the state system (4)-(6), as obtained using the approaches in subsections 4.1 and 4.2 with Taylor or Chebyshev model arithmetic, have convergence order q + 1 on P whenever the functions f k , h 0 , and h k are (q + 1)-times continuously differentiable in x and p. We argue here that the error parameterization bounds R q sj (t) (resp., R q λi (t)) in polynomial model inclusions of the forward sensitivity system (9)-(11) (resp., adjoint sensitivity system (13)-(15)) also have convergence order q + 1 on P when the participating functions ϕ k , η 0 , and η k are sufficiently smooth. Since the time-varying disturbance thereof is enclosed by R q x (t), which is known to have convergence order q + 1, the same method of proof as in [49] can be used to establish this convergence result.
As far as stability is concerned, it has been established in [20] that stable reachable set inclusions on infinite time horizons can be computed for nonlinear dynamic systems in the vicinity of their locally asymptotically stable equilibrium points, or their stable periodic orbits insofar as the cycle time is independent of the parameters, when polynomial models with ellipsoidal remainders are propagated by the set-valued integrator. The reason behind this stabilization is that the magnification of the parameterization error due to the wrapping effect is overpowered by the natural contraction of the trajectories for such systems. Since both the forward and adjoint sensitivity systems inherit the asymptotic stability properties of the state system [7], we argue that forward and adjoint system bound propagation will enjoy stability properties similar to state bound propagation.

5.
Numerical case studies. The main objective of the numerical case studies below is to illustrate various properties of the proposed set-valued ODE sensitivity methods and make comparisons between several options or variants. Our implementation of these methods comes in the form of a C++ class called ODEBNDS as part of the library CRONOS [9], which is available from https://github.com/omega-icl/cronos. We make use of the library MC++ (https://github.com/omega-icl/mcpp) for creating and manipulating expression trees, and for bounding these expressions using Taylor or Chebyshev models, in combination with verified interval libraries such as PROFIL (http://www.ti3.tu-harburg.de/). For numerical ODE integration, we use the code CVODES in SUNDIALS [17] (version 2.7), which comes with both forward and adjoint sensitivity analysis capabilities.
Unless otherwise noted, the following case studies use the BDF scheme with a maximum order of 5, the simultaneous corrector method with analytical sensitivity right-hand sides for forward sensitivity analysis, and a cubic Hermite interpolation without check-pointing for adjoint sensitivity analysis. The option of an approximate diagonal Jacobian formed by way of a difference quotient is used for simplicity. The relative and absolute numerical integration tolerances for the state, sensitivity, and adjoint systems are set to 10 −8 and 10 −10 , respectively, for all the variables, apart from the shape matrix entries in the ellipsoidal remainders whose absolute tolerances are set to 10 −20 . 1 The reported CPU times are for an Intel(R) Core(TM) i7-3770

Cubic oscillator. Consider the dynamic systeṁ
with uncertain initial conditions x 1 (0) = p 1 ∈ [1.5, 2.5] and x 2 (0) = p 2 ∈ [−0.5, 0.5], describing an oscillator of mass m = 0.1. Our focus is on the cumulative squareddeviation of the oscillator with respect to the origin, measured as over a time horizon of T = 100. Bounds derived from 3rd-and 5th-order Chebyshev models with ellipsoidal remainders are shown in Figure 1 for selected state, sensitivity, and adjoint trajectories. These bounds are found to be tight throughout the oscillations for 5th-order Chebyshev models, whereas the sensitivity and adjoint bounds are progressively diverging in their respective directions of integration with 3rd-order models. At this point, we would like to recall that the cubic oscillator problem (57), (58) presents a stable periodic orbit, in the vicinity of which the state bounds can be stabilized [20]. The results in Figure 1 confirm that the sensitivity and adjoint bounds may enjoy a similar stability for this problem, but the size of uncertainty or the expansion order for which such stabilization occurs is different for the state, forward sensitivity, and adjoint sensitivity problems.
Approximate inner bounds, using Monte Carlo sampling, for the derivatives of (59) are given by In the case of 4th-order Chebyshev models, the bounds remain rather inaccurate, with the remainder bounds having the same magnitude as the polynomial approximant In the case of 8th-order Chebyshev models finally, the remainder bounds are in the order of 10 −3 , so the polynomial approximant may be assimilated to the actual reachable sets, Notice the residual gap between the derivative enclosures and sampled ranges, despite the remainder bounds becoming small. This conservatism is due to the bounding of the multivariate polynomial approximant, which is not exact here [38]. These results are also confirmed by the plots of various reachable tube enclosures in Figure 2, both for the sensitivity variables at time t = T and the adjoint variables at t = 0. Finally, we note that the forward sensitivity bounds are slightly less conservative than their adjoint counterparts, although the latter are computationally cheaper since the problem consists of a single function but two parameters.
The problem involves bounding the derivatives of the least-squares function where the measurements S m 1 (t k ), S m 2 (t k ), and C m (t k ) are every 6 hours over a 3-day horizon, as given in Table 2.
Bounds derived from 2nd-order Chebyshev models with both interval and ellipsoidal remainders are shown in Figure 3 for selected state, sensitivity, and adjoint trajectories. For such a low expansion order, the forward and adjoint sensitivity bounds  Fig. 3. Projections onto S 2 (left), ∂S 2 /∂X 2 (0) (middle), and λ S 2 (right) of the reachable tube enclosures for the state, sensitivity, and adjoint trajectories. The bounds derived from 2nd-order Chebyshev models with interval and ellipsoidal remainders correspond to the blue and red lines, respectively. Inner approximations for the actual reachable sets are represented in shaded gray. Color is available in the online version only.  appear to be diverging when interval remainders are considered, whereas those corresponding to ellipsoidal remainders are significantly less conservative. Notice also the discontinuous adjoint trajectories at the times t k in (71).
A more systematic comparison of the conservatism associated with Chebyshev models of 2nd-, 3rd-, and 5th-order is presented in Figure 4 for various sensitivity bounding approaches. The overapproximation thereof is measured in terms of the radius of a Chebyshev model's remainder, the magnitude of which should be compared with the following approximate inner bounds (using Monte Carlo sampling) for the derivatives of (71): A general trend in Figure 4 is the tighter derivative bounds that are obtained by propagating (i) higher-order Chebyshev models and (ii) ellipsoidal remainders rather than interval remainders-two strategies known for mitigating the wrapping effect and improving the stability of the bounds. Moreover, Chebyshev models of order 5 or higher produce very little conservatism in this problem.
Similar to the previous case study, forward sensitivity analysis appears to provide more accurate bounds than adjoint sensitivity analysis. However, this is not without additional computational overhead. The computational comparison in Table 3   with the adjoint sensitivity approach (considering the combined state-adjoint system). These differences are consistent with the fact that the problem has a single function and three parameters. Also, the reason for not observing a three-fold improvement with the adjoint sensitivity approach is linked to the larger number of integration steps taken by the discontinuous adjoint problems.

Exothermic batch reactor.
We consider the model of a batch reactor with a cooling jacket, where a first-order exothermic reaction A → B is taking place. The ODEs describing the overall conversion X and the temperature T arė with initial values X(0) = 0 and T (0) ∈ [350, 370] K. The cooling jacket temperature profile T j is discretized as a piecewise constant control over N stages, . Values for all the other parameters are the same as in [12]. The problem involves bounding the derivatives of the final conversion X(t N ) at t N = 900 s with respect to the uncertain parameters p = (T (0), T  Table 4. The overestimation present in the computed derivative bounds with the forward and adjoint sensitivity bounding methods described in section 4 is 1-2 orders of magnitude smaller for 2nd-order Chebyshev models, and 2-3 orders of magnitude smaller with 3rd-order Chebyshev models, than the derivative ranges represented in shaded gray on Figure 5. A general trend is that interval remainder propagation provides tighter bounds than ellipsoidal remainder propagation,  an indication that the dynamics (72), (73) are not causing too much wrapping here. Nonetheless, interval remainder propagation typically entails a higher computational burden than its ellipsoidal counterpart, due to a larger number of integration steps.
Quite expectedly when the number of parameters increases from n p = 3 to 9, adjoint sensitivity bounding becomes much cheaper than forward sensitivity bounding. A difference with real-valued adjoint integration, however, is that the computational cost of propagating adjoint bounds increases significantly with more parameters, all else being the same, due to a larger number of auxiliary bounding ODEs. Another general observation is that the computed adjoint sensitivity bounds are tighter than their forward sensitivity counterparts here. Table 4 and Figure 5 also compare the forward and adjoint sensitivity bounding methods described in section 4 with alternative forward sensitivity methods that use existing set-valued integration to bound an extended state-sensitivity system directly; that is, these methods do not exploit the built-in forward sensitivity capabilities of CVODES, and in particular the remainder bound equations are different from those presented in section 4. Two variants thereof are whether one single state-sensitivity system with n x (n p + 1) is considered for all the parameters (all-at-once approach) or separate state-sensitivity systems with 2n x are considered for each parameter (one-byone approach). The results in Table 4 and Figure 5 clearly show that neither of these variants is competitive with the proposed methods in terms of computational effort. The corresponding sensitivity bounds are also much weaker when ellipsoidal remainders are used, due to the need for propagating these ellipsoids in higher-dimensional spaces. 5.4. Belousov-Zhabotinskii reaction system. We consider the following reaction-diffusion system [32], given in dimensionless form by the partial differential equations (PDEs) with boundary conditions for u 1 and u 2 as This problem presents two uncertain kinetic parameters r ∈ [9.8, 10.2] and b ∈ [1.225, 1.275], and it involves bounding the derivatives of the average value of u 1 (T, ·) with respect to these parameters at a given time T > 0.
Following [32], we discretize this PDE system into an ODE system by applying the method of lines over the finite spatial domain −30 ≤ x ≤ 20:        Moreover, we simulate a moving front by defining the following initial conditions: here assuming N is even. The average value of the concentration u 1 (T, ·) over the spatial domain is approximated by the function with the time horizon set to T = 20, and the uncertain parameter vector p = (r, b). Approximate inner bounds (using Monte Carlo sampling) for the derivatives of (81) take values in the order of 10 −3 for discretization levels N in the range of 10-50. For N = 50, the derivative bounds computed with 5th-order Chebyshev models with ellipsoidal remainders are accurate, having remainders in the order of 10 −7 and 10 −6 with the forward and adjoint sensitivity bounding approaches, respectively. These tight derivative bounds are reflected in Figure 6 showing projected enclosures for selected sensitivity and adjoint trajectories, which are indeed tight across the discretized spatial domain. An interesting observation here is that the tightness of the derivative bounds increases as the mesh size is refined. For a coarser discretization level of N = 20, for instance, remainders an order of magnitude larger are obtained with both the forward and adjoint sensitivity bounding approaches. We also note that the propagation of interval remainders does not produce meaningful derivative bounds in this case study. This discretized ODE system (77)-(78) provides a means to investigate how the proposed sensitivity-bounding approaches scale with an increasing number of states by refining the spatial discretization. Figure 7 compares the performance of the state, forward sensitivity, and adjoint sensitivity set-valued integrators, with either 3rd-or 5th-order Chebyshev models with ellipsoidal remainders, as a function of the discretization size. Generally, it is the discretization size N , and not the polynomial expansion order q, that is responsible for increasing the computational burden here. Computational comparison of the state, sensitivity, and adjoint bounding problems with 3rdand 5th-order Chebyshev models in terms of the overall CPU time (left) and the time per righthand side evaluation of the bounding ODEs (right). The red, blue, and green lines correspond to the state, sensitivity, and adjoint bounding systems, respectively. The solid lines are for 3rd-order Chebyshev models and the dashed lines for 5th-order Chebyshev models. Color is available in the online version only.
In fact, 5th-order expansions even turn out to be faster than 3rd-order expansions, due to their taking significantly fewer integration steps, which balances the larger computational effort for computing the right-hand side of the auxiliary bounding ODEs. This is because the size of the auxiliary bounding system in the ellipsoidal remainder approach is dominated by (n x ) 2 rather than 2 q n x for low expansion orders q in this two-parameter problem. Finally, Figure 7 reveals that adjoint sensitivity bounding is consistently faster than forward sensitivity bounding, for the problem contains two parameters and a single function.
6. Conclusions. The main contribution of this paper is an extension of continuous-time set-valued integration methods to enable sensitivity analysis of parametric ODEs. An approach to propagating polynomial model inclusions for either of the forward and adjoint sensitivity systems has been described, which entails a specialization of generic continuous-time set-valued integration for ODEs with both parametric and time-varying uncertainties. Similar to classical sensitivity analysis, bounds for the forward sensitivity systems are propagated forward in time, alongside the state bounds, whereas bounds for the adjoint sensitivity systems call for a backward pass and interpolation between grid points of the state bounds computed during a forward pass. Two variants have been described, whereby the remainder terms in the forward/adjoint sensitivity polynomial model inclusions, which enclose the forward or adjoint sensitivity parameterization errors, are interval boxes or ellipsoids. A key aspect of the proposed approach is that the state parameterization error acts as a time-varying uncertainty in either of the sensitivity bounding approaches. Other setvalued integration approaches could be extended likewise in order to enable sensitivity bounding. The paper has also described an implementation of the set-valued sensitivity integrators that builds upon, and takes advantage of, state-of-the-art numerical integrators with sensitivity capability. From a practical standpoint, the main task thereby involves evaluating the right-hand side of the auxiliary dynamic bounding system and of the corresponding transition conditions. An open source implementation is made available as a C++ class in the library CRONOS, based on the numerical integrator CVODES in the SUNDIALS suite. The latter enables the use of simultaneous correcc 2017 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license tor methods-and potentially staggered corrector methods too-in forward sensitivity analysis and check-pointing techniques in adjoint sensitivity analysis.
Finally, a number of numerical case studies have been presented, which compare the performance of the forward and adjoint sensitivity bounding variants. For problems having more than just a few parameters, adjoint sensitivity is computationally advantageous compared with forward sensitivity in bounding the derivatives of a single or a few state-dependent functions, a situation similar to real-valued integration. The derivative bounds computed by propagating polynomial model inclusions with ellipsoidal remainders are typically tighter than with interval remainders, due to a lower prevalence of the wrapping effect, but which of the derivative bonds computed with forward or adjoint sensitivity analysis are tighter appears to be problem dependent. It is also shown that the forward and adjoint sensitivity bounds may remain stable when the state bounds are themselves stable-e.g., in the vicinity of asymptotically stable equilibrium points or periodic orbits of the dynamic system at hand-by propagating polynomial model inclusions with ellipsoidal remainders.
As part of future work, this generic sensitivity bounding capability could be extended to compute bounds on second-order directional derivatives [35] via the combination of forward and adjoint sensitivity bounding. One promising application is in the area of complete-search methods for global optimization and constraint satisfaction with dynamic systems embedded in order to expedite convergence.