SET-VALUED INTERPOLATION, DIFFERENTIAL INCLUSIONS, AND SENSITIVITY IN OPTIMIZATION

. Set-valued interpolation and integration methods are introduced with special emphasis on error representations and error estimates with respect to Hausdor(cid:11) distance. The connection between order of convergence results and sensitivity properties of (cid:12)nite-dimensional convex optimization problems is discussed. The results are applied to the numerical approximation of reachable sets of linear control problems by quadrature formulae and interpolation techniques for set-valued mappings.


Introduction
Numerical methods for the solution of di erential inclusions follow three directions.
(i) Compute special solution trajectories with additional qualitative or quantitative properties: Such trajectories have to be computed by di erence methods with additional selection procedures choosing points from the set-valued right-hand side in an appropriate way. Common strategies result, e.g., in the discrete analogue of heavy solutions, or slow solutions, or selections with a certain kind of discrete minimal variation. These selection procedures require the solution of nite-dimensional optimization problems at every gridpoint. Qualitative and quantitative sensitivity properties of this family of optimization problems determine qualitative and quantitative properties of the approximated solution, cp. in this connection 23] and the survey 19].
(ii) Compute all solution trajectories, or at least all belonging to a given class of functions: This is theoretically and computationally an extremely di cult task. Applying the abstract framework of general discretization theory requires correct notions of stability and consistency. Conditions assuring order of convergence higher than 2 are not available until now. In principle, a proper calculus of higher order derivatives is required for set-valued mappings, guaranteeing Taylor expansions with valid error estimates with respect to Hausdor distance. Some results concerning Euler's method resp. Euler-Cauchy method and order of convergence equal to 1 resp. equal to 2 are available, cp. 33], 34]. Every solution belonging to an appropriate Sobolev space can be approximated in a theoretical sense by a higher order linear multistep method, where the relevant notion of consistency is related to stability properties of a family of perturbed optimization problems, cp. De nition 3.2 in 23].
(iii) Compute the reachable set of all solution trajectories at a prescribed point in time: The techniques mentioned in (ii) like Euler's method resp. Euler-Cauchy method yield, as a by-product, rst resp. second order discrete approximations of reachable sets of special classes of di erential inclusions. In 16] even higher order of convergence is proven for a method exploiting fully the structure of special linear di erential inclusions with polyhedral control region. In the sequel of papers 7], 6], 4], and in the thesis 5], the discrete approximation of reachable sets of linear di erential inclusions is totally reduced to the numerical integration of set-valued mappings. The basis of this approach consists in adaptations of quadrature formulae and extrapolation methods to the calculation of Aumann's integral for set-valued mappings. In principle, classical quadrature methods are applied to the support functional of the set-valued integrand. For every point in the integration interval and every unit vector in state space, the value of the support functional is determined by a convex optimization problem. Smoothness properties of this support functional as a function on the integration interval uniformly with respect to the unit ball in state space, thus strong stability and sensitivity properties of an in nite family of convex optimization problems, determine the order of the integration method and, consequently, the order of suitably de ned discrete approximations of reachable sets. In this framework, higher order discrete approximations to reachable sets can be de ned at least for special classes of linear di erential inclusions. Originally, only the use of quadrature formulae with nonnegative weights seemed to be reasonable, like some open or closed Newton-Cotes formulae, Gauss quadrature, or Romberg's extrapolation method with Romberg's stepsize sequence. But exploiting some ideas in 8], compare also 9], depending on the geometry of the set-valued integrand, even quadrature fomulae with negative weights could be applied, thus opening the way to all kinds of extrapolation methods, error estimates by inclusion, and stepsize control for set-valued integration.
As outlined above, there exists an intrinsic relationship between numerical methods for di erential inclusions and questions of sensitivity and stability analysis of nite dimensional optimization problems. The main objective of this paper is to describe this relationship. Hoping, that a numerical treatment of linear di erential inclusions in the very spirit of set-valued numerical analysis will also be of value for a more satisfactory numerical treatment of nonlinear di erential inclusions, we will concentrate on aspect (iii). Contrary to the thesis 5], where set-valued integration is the exclusive mathematical tool, we try to broaden the mathematical background to setvalued interpolation. The reader will easily recognize, that the techniques apply to set-valued mappings of several variables as well, thus opening the access to nite element methods for the discrete approximation of nonlinear di erential inclusions in the, hopefully, near future.

Set-Valued Interpolation
In the following, we introduce set-valued interpolation as a mathematical tool to approximate set-valued mappings by simpler set-valued mappings. Deliberately, we avoid the technique of embedding spaces of convex sets into normed linear spaces, cp. the papers 28], 21], 30], 10], and 18]. This technique leaves the question unanswered how to interpret the results in the original spaces. Instead, we stay completely in the framework of setvalued mappings. Naturally, the problem arises how to de ne di erences of sets in an appropriate way. This is done by a method already used in 8] for the proof of error estimates for set-valued quadrature formulae with negative weights, and in 5] for the derivation of inclusions of set-valued integrals by extrapolation methods. Only for simplicity we restrict ourselves to interpolation by set-valued polynomials, extensions to other function classes and even to interpolation of set-valued mappings of several variables by set-valued nite elements being rather obvious.
2.1. Interpolation Problem. Let I = a; b] with a < b and F : I =) R n be a set-valued mapping with non-empty, convex and compact values. Choose N 2 N and a grid a t 0 < t 1 < : : : < t N b ; and compute for every l 2 R n the polynomial p N (l; ) of degree N with p N (l; t j ) = ? (l; F(t j )) (j = 0; : : : ; N) : which can be computed as follows.
By de nition, cp. e.g. 29], we have p ? N (z; t) = sup l2R n z ? l p N (l; t)] = 8 < : 0; if z ? l p N (l; t) for all l 2 R n ; 1; if z ? l > p N (l; t) for at least one l 2 R n : Hence, p ? N (z; t) is the indicator function of the set P N (t) = fz 2 R n : z ? l p N (l; t) for all l 2 R n g ; (2.2) and therefore p ?? N ( ; t) = ? ( ; P N (t)) is the support functional of P N (t) for every t 2 I. 2.2. Lemma. On the standard assumptions of Interpolation Problem 2.1, the set P N (t) is closed, convex and bounded for every t 2 I. Proof. According to (2.2) the set P N (t) is the intersection of closed half spaces in R n , therefore P N (t) is convex and closed. Moreover, (2.1) shows that p N (l; t) is bounded uniformly for all l 2 R n with klk 2 = 1, p N (l; t) c(t) (klk 2 = 1) ; (2.3) this implies for z 2 P N (t) kzk 2 2 p N (z; t) and hence, for kzk 2  ; t c(t) (t 2 I) : Since p ?? N (l; t j ) = ??? (l; F(t j )) = ? (l; F(t j )) and F(t j ) is closed and convex, P N (t j ) = F(t j ) (j = 0; : : : ; N) : Therefore, in a very natural way, we can de ne the set-valued interpolation \polynomial" which solves Interpolation Problem 2.1.
2.3. De nition. For every l 2 R n let p N (l; ) be the interpolation polynomial which solves Interpolation Problem 2.1. Then the set-valued mapping P N : I =) R n ; de ned by P N (t) = fz 2 R n : z ? l p N (l; t) for all l 2 R n g (t 2 I) ; is called the set-valued solution of Interpolation Problem 2.1.
At this point, we should add a warning: Neither is P N (t) in general polynomial with respect to t, nor is P N (t) necessarily non-empty for all t 2 I. Hence, it is crucial to give conditions which guarantee P N (t) 6 = ; for all t 2 I. In addition, these conditions should allow the proof of error estimates with respect to Hausdor distance between F(t) and P N (t) which are analogous to error estimates between the scalar functions ? (l; F(t)) and p N (l; t). For this purpose, we use the following result which was already exploited in 8] for the proof of error estimates for set-valued quadrature formulae with negative weights.
2.4. Lemma. Consider a xed t 2 I where p N ( ; t) is not itself a support functional. Assume moreover, that there exists a ball B(m(t); r(t)) = fz 2 R n : kz m(t)k 2 r(t)g with center m(t) 2 R n and radius r(t) > 0, which is contained entirely in P N (t), B(m(t); r(t)) P N (t) : De ne, as in (2. Then the following error estimate holds haus (F (t); P N (t)) 2c(t) r(t) sup klk 2 =1 j ? (l; F(t)) p N (l; t)j : Here, haus( ; ) denotes Hausdor distance with respect to Euclidean norm k k 2 . The proof is contained in 8] and 5]. More convenient in applications is the following condition on F(t) itself.
2.5. Corollary. Consider again a xed t 2 I where p N ( ; t) is not itself a support functional. Assume moreover that the ball B(m(t); r(t)) with center m(t) 2 R n and radius r(t) > 0 is contained entirely in F(t).
If for a xed t 2 I the interpolating function is itself a support functional, which is clear for all grid points, and for linear interpolation or other interpolation techniques with non-negative basis functions, then the error estimate does not depend any longer on the geometry of the set-valued mapping F( ). Then the following estimate, cp. 11], 21], holds.
2.6. Lemma. Consider a xed t 2 I where p N ( ; t) is itself a support functional of a non-empty convex and compact set P N (t). Then haus (F (t); P N (t)) = sup The last representation of Hausdor distance is extremely useful for the direct proof of error estimates for set-valued quadrature formulae with nonnegative weights without recourse to set-valued interpolation, cp. 7], 6], 4], 5], and Section 5.
By Lemma 2.4, Corollary 2.5 and Lemma 2.6, the error between F(t) and P N (t) with respect to Hausdor distance is totally reduced to the classical error between ? (l; F(t)) and p N (l; t) and, eventually, some upper bounds for c(t) and positive lower bounds for r(t) which depend on the geometry of P N (t) resp. F(t). As we will see in Section 3, continuity and di erentiability properties of ? (l; F(t)) with respect to t 2 I uniformly for all l 2 R n with klk 2 = 1 play a crucial role for the classical error.
But, we want to stress that such regularity properties of ? (l; F( )) can only be expected to hold for special classes of set-valued mappings F( ), cp. Section 4. In any case, all subsequent error representations and error estimates have to be done very cautiously to exploit at least some absolute continuity properties for reasonably large classes of problems.

Representation of the Interpolation Error
There are several methods, to prove estimates for the interpolation error R N (l; t) = ? (l; F(t)) p N (l; t) : One could follow classical lines, cp. e.g. 32], which usually requires a little bit too strong smoothness assumptions on ? (l; F( )). One could also follow an approach exploiting systematically moduli of smoothness of ? (l; F( )), cp. 31]. This approach yields the weakest estimates for problems in one variable. Instead, we present an approach which leads to weak error estimates for an especially important class of problems, and which can easily be extended to interpolation problems in several variables, cp. 15] and 14].
For simplicity, in this section we use the abbreviation f(t) = ? (l; F(t)) and suppress the explicit indication of l whenever possible.
If in addition f (N 1) ( ) is absolutely continuous, then f (N ) ( ) exists almost everywhere and is integrable on I, partial integration is justi ed and gives f(t j ) = f(t) + f 0 (t)(t j t) + : : Hence, for almost all t 2 I, we have the representation, which holds for Consider t as a xed parameter, then the polynomial of degree at most equal to N p(z) = f(t) + f 0 (t)(z t) + : : satis es d dz p(z) j z=t = f ( ) (t) ( = 0; : : : ; N) (3.2) and, for j = 0; : : : ; N, Therefore, it coincides with the Lagrange interpolation polynomial of degree at most equal to N which attains the same values at the nodes t j , we get from (3.2) the following 3.1. Error Representation. Let f (N 1) ( ) be absolutely continuous on I. Then for = 0; : : : ; N and almost all t 2 I the following representation This error representation clearly shows that the variation of f (N ) ( ) on I plays a crucial rôle, where this variation has to be de ned in an appropriate way, since f (N ) ( ) is only integrable. Fortunately, f (N ) ( ) appears only in integrated form. Hence, the following de nition is su cient for our purposes, cp. 5], p. 15, g( ) : I ! R n is integrable and g(t) = f (N ) (t) for almost all t 2 I ; (3.3) where var( ) denotes the usual variation of a vector valued function with respect to Euclidean norm. In the rest of this paper, the variation of integrable functions is to be understood in the sense of (3.3).
Assuming f( ) to be absolutely continuous and N = 1; = 0, we get for almost all t 2 I as a very special case. This example is special in another sense as well. Since the basis functions t t 1 t 0 t 1 ; t t 0 t 1 t 0 are nonnegative on I = t 0 ; t 1 ], is itself a support functional of the non-empty convex and compact set hence Lemma 2.6 applies, and (3.4) results directly in the following error estimate for the Hausdor distance between F(t) and P 1 (t).

Linear
Interpolation. Let ? (l; F( )) be absolutely continuous, and let d dt ? (l; F( )) be of bounded variation in I uniformly for all l 2 R n with klk 2 = 1.
Then, for linear set-valued interpolation, the following error estimate holds haus(F (t); P 1 (t)) = sup Naturally, such error representations suggest the use of piecewise polynomial interpolation of set-valued mappings to get error estimates in terms of stepsize. Piecewise linear interpolation leads to corresponding error estimates for the composite trapezoidal rule for set-valued mappings which is the basis for extrapolation methods for set-valued integration, cp. Section 5 and 7], 4], 5].
For later use, we add another special case, set-valued interpolation by polynomials of second degree. Assuming now d dt f( ) to be absolutely continuous and N = 2; = 0, we get Now, clearly, for xed t 2 I is positively homogeneous with respect to l, but not any longer necessarily convex, since the Lagrangean elementary polynomials generally have different signs. Therefore, applying Corollary 2.5, we get the following error representation for set-valued quadratic interpolation. ? (l; F( )) of bounded variation in I uniformly for all l 2 R n with klk 2 = 1. Assume moreover that for t 2 I the ball B(m(t); r(t)) with center m(t) 2 R n and radius r(t) > 0 is contained in F(t), and that (t) = sup klk 2 =1 j ? (l; F(t)) p 2 (l; t)j is small enough, i.e. 0 (t) < r(t). Let c(t) = sup klk 2 =1 p 2 (l; t) .
All Here, the perturbation parameter is t 2 I, the vector l 2 R n is considered to be xed, and continuity and di erentiability properties of the corresponding value function ? (l; F(t)) with respect to t 2 I uniformly for all l 2 R n with klk 2 = 1 play a crucial rôle. Naturally, such additional regularity properties of ? (l; F( )) can only be expected to hold for special classes of set-valued mappings F( ) or for some concrete problems. Then, ? (l; F( )) is Lipschitz continuous, and d dt ? (l; F( )) of bounded variation uniformly for all l 2 R n with klk 2 = 1.
(b) Let U R m be compact and non-empty, and let the n m-matrix function A( ) describe the following parametrization of F,

F(t) = A(t)U (t 2 I) :
Let A( ) be absolutely continuous, and d dt A( ) of bounded variation. Then, ? (l; F( )) is Lipschitz continuous, and d dt ? (l; F( )) of bounded variation uniformly for all l 2 R n with klk 2 = 1.
(iv) Consider a real function : I ! R and a nonempty subset U R n . (v) Let ( ) be an n m-matrix function and B(m(t); r(t)) = fz 2 R m : kz m(t)k 2 r(t)g (t 2 I) a varying ball in R m with center m(t) 2 R m and radius r(t) 0. De ne F(t) = (t)B(m(t); r(t)) (t 2 I) : Then ? (l; F(t)) = l ? (t)m(t) + r(t)k ? (t)lk 2 (t 2 I) : Hence, as long as ? (t)l 6 = 0 R m (which is, e.g., the case for all l 2 R n with klk 2 = 1 if the rows of (t) are linearly independent for all t 2 I), di erentiability properties of m( ); r( ), and ( ) are inherited by ? (l; F( )).
The situation is much worse with a varying ball in R m with respect to in nity norm, B 1 (m(t); r(t)) = fz 2 R m : kz m(t)k 1 r(t)g (t 2 I) :
Summarizing, we want to stress that Theorem 4.2 just su ces to justify the error estimate for (piecewise) linear set-valued interpolation for broader classes of parametrized set-valued mappings, whereas Examples 4.3 justify even higher order set-valued interpolation by (piecewise) polynomials for more restricted classes of set-valued mappings. In the following sections, these results are used for the derivation of error estimates for set-valued integration and discrete approximations of attainable sets.

Set-Valued Integration
Set-valued integration can be introduced in di erent ways either following 10] and 18] exploiting abstract embedding theorems for spaces of convex sets 28], 21], or in a direct way for quadrature formulae with nonnegative weights, cp. 7], 4], 6], 5], resp. for quadrature formulae with negative weights, cp One drawback of the above approach is, that one needs additional geometric assumptions to guarantee that the interpolatory set-valued polynomial has non-empty values. For quadrature formulae with nonnegative weights, there is a direct approach avoiding this drawback and with a wider range of applicability even to set-valued Gau quadrature formulae and extrapolation methods, cp. 7], 4], 6], 5]. In the following, we give a brief sketch of this direct approach.
It is well-known that, under all assumptions of Theorem 5. jR(l; F)j : Obviously, the order of the error is determined again by regularity properties of the scalarized integrand ? (l; F(t)) (t 2 I) with respect to t uniformly with respect to all l 2 R n with klk 2 = 1.
To be more speci c, we cite only two special cases, the set-valued analogues of composite trapezoidal rule and composite Simpson's rule. In fact, these two methods form the rst two stages of a set-valued analogue of Romberg's extrapolation method, cp. 4], 6], and 5] for more details and complete proofs. 5.4. Theorem. Let F : I ) R n be a measurable and integrably bounded set-valued mapping with non-empty, convex, and compact values, and let ? (l; F( )) be absolutely continuous with rst derivative with respect to t of bounded variation uniformly for all l 2 R n with klk 2 = 1.  If, in addition, ? (l; F( )) has an absolutely continuous second derivative and if its third derivative is of bounded variation uniformly for all l 2 R n with klk 2 = 1, then the error of composite Simpson's rule h 3 F(t 0 )+4F (t 1 )+2F (t 2 )+4F (t 3 ) + : : : +2F (t N 2 )+4F (t N 1 )+F (t N )] on the same grid for even N is of order 4 in h with respect to Hausdor distance.
Under additional smoothness assumptions Romberg extrapolation gives approximations of even higher order with respect to Hausdor distance, compare the examples and numerical tests in 4], 6], 5].
We want to emphasize that, due to Theorem 4.2, we get second order convergence for a remarkable big class of parametrized set-valued mappings. Due to Examples 4.3, we get even higher order of convergence for more speci c problem classes. Exploiting inclusion techniques, based on the asymptotic expansion of composite trapezoidal rule, for such problem classes, one can even get inner and outer approximations to the exact Aumann integral which converge of the same order with respect to Hausdor distance as the underlying Romberg extrapolation method, cp. 5].

Approximating Reachable Sets by Set-Valued Integration and Interpolation
We want to outline in the following how set-valued integration and interpolation could be applied in a combined form to the discrete approximation of reachable sets of linear di erential inclusions. To be more speci c, we analyse the following standard 6.1. Linear Control Problem. Let the n n-matrix function A( ) and the n m-matrix function B( ) be (at least) integrable on I = a; b], and the control region U R m be non-empty, compact, and convex. Let Y 0 R n be a non-empty, compact, and convex starting set.
Find an absolutely continuous function y : I ! R n with y 0 (t) 2 A(t)y(t) + B(t)U for almost all t 2 I ; y(a) 2 Y 0 : The can only be computed up to an error of order q with respect to some chosen stepsize h according to procedures a) and b), one never will get the exact set-valued interpolatory polynomial P q ( ), but only some set-valued approximationP q ( ). Assuming, that the error haus P q (t);P q (t) The grid data for interpolation result from set-valued numerical integration with nonnegative weights, and we recommend the dual approach by supporting hyperplanes discussed thoroughly in 5] including an analysis of the additional error due to the fact that only nitely many hyperplanes can be computed. For higher dimensional problems, actually support points of these hyperplanes should be plotted, not the hyperplanes themselves. The dual representation of^ (t; a)P q (t) between gridpoints is more di cult, since in the representation (2.2) ofP q (t) even non-supporting hyperplanes will occur, i.e. one has to develop an appropriate device to get rid of such hyperplanes and to restrict this representation to actually supporting hyperplanes resp. corresponding support points. The complexity of this task is dimension dependent. In the following numerical tests, we restrict ourselves to state space dimension 2.
The combined procedures a), b), c), and d) are visualized for the fol- Due to Theorems 4.2, 5.4, and 6.2, the sets Y (1); Y (1:5); resp: Y (2) can be computed by a combined method of order 2 very precisely. The results are plotted in Figures 1, 2, resp. 3. Since the control region is polyhedral, the methods in 16], 20] could have been used equally well.   According to De nition 2.3, this set-valued mapping is related to a scalar quadratic interpolation polynomialp 2 (l; ) in the following way, P 2 (t) = n z 2 R 2 : l ? z p 2 (l; t) for all l 2 R 2 o : (6.2) This opens the way, at least in R 2 , for a dual representation ofP 2 (t) as intersection of halfspaces. In Figure 4, the image of this representation under the linear transformation^ (t; 1) is visualized for t = 1:1. Observe, that not all of these halfspaces are supporting ones due to the non-convexity ofp 2 (l; 1:1) with respect to l, cp. the magni cation around the upper left corner in Figure 5.  Observe, that, in principle, it would have been su cient to compute the data sets for set-valued interpolation, cp. Figures 1, 2, 3, within an error of order 2 with respect to the stepsize for interpolation to get an overall error of order 2. In practical computations, one has to restrict oneself in (6.2) to a nite collection of vectors l 2 R 2 with klk 2 = 1. Therefore, to retain order of convergence equal to 2 for the actually computed discrete approximations in Figures 6,7,8,9, one has to choose the vectors l from an appropriately adjusted grid on the unit sphere in R 2 .
Acknowledgement. I appreciate the help of Robert Baier in preparing the plots. They are essentially based on his program package on the approximation of reachable sets of linear di erential inclusions by set-valued integration methods.