A rigorous deterministic global optimization approach for the derivation of secondary information in digital maps

We derive a generic system that constructs an optimization model for an emergency stop scenario on the highway, based on map data from high definition maps that are used in Advanced Driver Assistance Systems (ADAS) and in Highly Automated Driving (HAD). New additional situative and scenario-based information is computed by applying a global maximization approach to the model. For this purpose, we develop two new rigorous and deterministic branch-and-bound algorithms that both determine the certified global optimal value up to a predefined tolerance. The underlying interval optimization algorithm, which uses first-order techniques, is enhanced by one of two second-order methods that are applied for specifically selected intervals. We investigate two approaches that either compute a concave overestimator for the objective function or approximate the function with a quadratic polynomial using Taylor expansion. We show the limits of interval arithmetic in our problem, especially for the interval versions of the derivatives, and present a local linearization of the curve data that improves the results significantly. The presented novel method for deriving secondary information is compared to state of the art methods on two exemplary and for the automotive context representative scenarios to show the advantages of our approach.


Introduction
Digital maps provide a range of important information for many applications in Highly Automated Driving (HAD). Since, nevertheless, the information in these maps is often not sufficient for some applications, we will present a method to compute additional locally and temporarily valid map data that hold for a specific vehicle in that situation. We will create a model based on the digital map data and on some vehicle information, and we develop a global optimization approach that derives those additional data from that model.

Technical background
Driver assistance functions and HAD are closely linked to future mobility. The resulting challenges in perception and recognition make large demands on connectivity and data management, and they open a broad field of research. It is important to provide high definition maps (Seif and Hu 2016;Liu et al. 2020; Schwab and Kolbe 2019) that consist of large amounts of street topology data, not only containing all lanes and boundaries, but also logical information about their connections. In addition, landmarks, e.g., street lamps and traffic signs, or other objects of interest, are included for further information. Apart from the given fixed map data, in the following referred to as primary information, we intend to derive additional scenario-based information, which we refer to as secondary information.
In this paper, we compute optimal emergency stops for next-generation driver assistance functions and automated driving functions. The suggested, novel approach for deriving secondary information from high definition maps is, however, designed to be generic so that it can be as well applied and adapted to many other use cases. We formulate a maximization problem with a nonlinear, nonconcave objective function that is explicitly defined by the emergency stop scenario. The global solutions to this problem describe new stopping points that are added to an additional secondary information layer of the digital map. We present two deterministic global optimization strategies, whereby one of the approaches is faster than all other deterministic algorithms that we have encountered.
The small highway section 1 in Fig. 1 indicates the large amount of map data that is handled by the optimization algorithm for a larger, more realistic map excerpt. The curves of different types in Fig. 1B are depicted in various colors. The visualized map data contribute to the objective function shown in Fig. 1A. Just the map elements that belong to a certain surrounding of the own vehicle's position, in the following referred to as ego position X 0 , are taken into account. The structure of the digital maps helps to filter for the relevant data, since they are usually divided into separate tiles, where only the closest tiles to X 0 are used. The geometric parametrization of the curves is Fig. 1 A The objective function is based on the digital map data below. The Cartesian Coordinates are given in meters. B The landmarks (two types, represented by red and green cylinders) and curves of different types in the digital orthophoto (©BKG 2021) define the above objective function different for the existing map formats. However, it is very beneficial for the derivation of secondary information to define the lanes by means of cubic B-spline curves (Piegl and Tiller 1997, chs. 2,3).
The following requirements on the solution method are relevant from the perspective of a car manufacturer or developer of automated driving functions: R.1 Three-dimensionality: Digital maps have evolved from two to three dimensions because of the importance of the height information for certain use cases. Therefore, the problem specifications explicitly require an optimization model that is not only defined along the map elements, but in the full three-dimensional space, to maintain the possibility of solutions that are off the available map elements. R.2 Genericity: Although we focus on a highway emergency stop here, the valuation of the map data and the developed algorithm must be able to cover general scenarios, e.g., urban scenarios. This allows to exploit synergy effects, as multiple functions can be built upon the same algorithmic foundation. As a further reason, the development of this kind of functionally safe optimization algorithm is very costly. Therefore, development costs can be shared across multiple modules. R.3 Determinism: The derivation of secondary information belongs to the category of applications for automated driving functions with high expectations to system safety. Thus, the optimization algorithm must be deterministic, according to the demanded Automotive Safety Integrity Level (ASIL) (International Organization for Standardization 2018). This functionally safe algorithm can also be used for assistant functions with less safety requirements. R.4 Real time: The real-time performance is essential, since decisions in a driving vehicle often must be made in a fraction of a second. The solutions are only temporarily valid, because a driving car changes its position continuously.
Although there exist practically efficient heuristic approaches to compute global optima (Hedar and Fukushima 2004;Hansen and Ostermeier 2001;Csendes et al. 2008), those algorithms neglect the required reproducibility of the solution stated in R.3. As a consequence, we develop a new hybrid optimization method that combines a deterministic global optimization method with local higher-order approaches. This requires enhancements in a number of aspects that use the structure of the considered class of optimization models to improve the results of the existing algorithms. The main focus lies on a fast algorithm (cf. R.4).

State of the art
Rigorous optimization algorithms (Floudas 2005, chs. 3-5, 11-13, Kearfott 1996, chs. 1, 5, Pintér 2002Rios and Sahinidis 2009;Neumaier 2004) identify global solutions with absolute certainty up to a predefined tolerance. They often are based on branchand-bound type methods (Mitten 1970) that divide the initial search space into smaller subsets (usually boxes) step-by-step (branching) and discard the irrelevant ones that do not contain a global optimum (bounding). Branch-and-bound methods often use convex underestimators or concave overestimators, such as those obtained by the αBB technique (Adjiman et al. 1998;Adjiman and Floudas 1996) or by bounding schemes for multilinear and other functions (Bao et al. 2015;Rikun 1997;Ryoo and Sahinidis 2001). They can be further enhanced, e.g., by range reduction techniques that generate tightened valid inequalities (Ryoo and Sahinidis 1996;. One of the most popular methods is the DIRECT (DIviding RECtangles) algorithm by Jones et al. (1993) which has been developed to overcome convergence speed problems of Shubert's algorithm (Shubert 1972), a sequential method using a global Lipschitz constant. Instead of applying this Lipschitz constant globally to estimate upper bounds for the function value on each box (in the maximization case), a locally valid positive constant is used to find potentially optimal boxes. Since then, there have been various attempts to improve or modify the DIRECT algorithm for specific types of objective functions (for instance, Gablonsky and Kelley 2001;Chiter 2006;Liu et al. 2015;di Serafino et al. 2011). Other rigorous methods are the Multilevel Coordinate Search (MCS) (Huyer and Neumaier 1999), which is regarded as a branch-without-bound type method that, such as DIRECT, evaluates a single point per box. Besides the box size and the subdivision level, the branching rules of MCS include the computation of an expected gain by a quadratic model that is based on a few point evaluations. Furthermore, the MULTK algorithm (Sergeyev and Kvasov 2017, ch. 4) uses locally approximated Lipschitz gradients for differentiable objective functions.
Moreover, global interval optimization algorithms introduced by Hansen (1980) and further investigated in Gau and Schrage (2004), (Hansen and Walster 2004, chs. 7-8, 12-13), Neumaier (2015) use interval arithmetic (Hickey et al. 2001;Hansen 1975), (Moore et al. 2009, chs. 2, 4, 5), so that the bounds on the function value ranges are used to make optimality statements for whole boxes at once. Based on the interval values, one can discard those subsets that cannot contain a global optimum. There are various refinements and extensions of the original interval algorithm with new bounding strategies, e.g., the Krawczyk Newton approach (Krawczyk and Neumaier 1986), interval constraint propagation (Kjøller et al. 2007), McCormick-based relaxation (Mitsos et al. 2009), and relatively recent implementations such as GOP (Pál and Csendes 2009).
Apart from these deterministic and rigorous approaches, we want to mention two stochastic algorithms. GLOBAL (Csendes et al. 2008) is a multistart clustering algorithm that uses the best locally minimizing samples to gain information about the regions of attraction, whereas CMA-ES (Hansen and Ostermeier 2001) is an evolution strategy that adapts the covariance matrix for the sampling of new points iteratively. We apply these methods later on to compare them to some deterministic algorithms and to get an impression of how good our solutions compete with stochastic approaches.

Key contributions of the paper
In this paper, we construct a system to derive secondary information for the given map data, i.e., we design a new type of model that values the various map elements and we develop optimization approaches that exploit the structure of this model as much as possible. Thereby, we build two hybrid branch-and-bound optimization algorithms that use interval arithmetic in the first place, based on Hansen's considerations (Hansen 1980). We illustrate and also verify the advantages of the interval-valued evaluation in comparison to simple point evaluations of other global methods. Then, second-order bounding methods that are only applied to certain boxes in the interval algorithm extend both algorithms to generate better function value bounds. Several modifications in both algorithms are needed due to the fact that the objective function is just piecewise differentiable.
The αBB approach locally overestimates the objective function by a concave function whose maximum can easily be computed. The techniques by Gerschgorin (1931) or Hertz (2009) require good estimates for the corresponding interval Hessian matrices to compute concave overestimators. The interval Hessian and the overestimation term have been modified by Meyer and Floudas (2005), Skjäl et al. (2012), Skjäl and Westerlund (2014) to construct tighter overestimators. In our concrete problem, the interval-valued outcome for the function evaluation and its derivatives are significantly overestimated, and, therefore, these interval matrix approaches do not provide meaningful overestimators. It is a key point of this paper to obtain sufficiently tight  2 The pipeline presents the various steps to derive secondary information from primary map data. We also refer to the corresponding sections that cover the details estimates for the parameter α, which defines the concave overestimator, without the need of computing interval matrices as in the mentioned approaches. To reach this goal, we exploit the form of the model and estimate the α-values for all map element functions separately. For the curves, in particular, we present a novel method that locally linearizes the curves in the digital map and computes the maximum curvature of this locally very accurate approximation of the part of the objective function corresponding to this curve (which, for short, we call the curve function).
In the second approach, we build local quadratic Taylor models of the objective function to benefit from existing theories about quadratic programming with box constraints (Burer and Letchford 2009;De Angelis et al. 1997). It is essential to bound the Taylor approximation error for a proper interval bound, e.g., presented by Bompadre et al. (2013). We again apply interval arithmetic to the map element functions individually, including the local curve linearizations for the curve functions, to calculate a tighter bound for the overall Taylor error than, e.g., in (Mitsos et al. 2009). Thereby, the distances between geometric primitives, e.g., the branch-andbound boxes and the lines as approximations of the curves, are used to speed up the function evaluation and to improve the result significantly for an interval-valued evaluation.

Structure of the paper
The proposed pipeline for deriving secondary information by means of a mathematical optimization is shown in Fig. 2. It depicts the particular subtasks with references to the corresponding sections. This paper is structured as follows: In Sect. 2, we derive the underlying optimization problem for the derivation of secondary information from digital map data and show how to adapt the generic formulation to the specific use case of deriving suitable emergency stops. Section 3 deals with suitable approaches for solving the optimization problem derived in Sect. 2. We introduce an interval algorithm as well as the computation of locally valid concave overestimators and quadratic approximations, that are combined and extended to construct a suitable algorithm. Additionally, we describe a one-dimensional optimization method that provides good initial values to reduce the computational effort of the interval algorithm significantly. In Sect. 4, we focus on aspects that are relevant for an efficient function evaluation, in particular exploiting the B-spline form of the curves. In Sect. 5, we evaluate our prototypical implementation for a highly relevant practical scenario of deriving emergency stops for advanced driver assistance and automated driving functions. Furthermore, we point out how the alterations of the optimization algorithms significantly improve the suitability and performance, and give guidance for further optimizations. In the last section, we summarize our results, discuss them, and give an outlook.

Mathematical problem formulation
The input data must be converted into a three-dimensional mathematical model in order to compute the optimal emergency stops. We define a suitable objective function that is based upon the map elements, the scenario, and the ego position. For every type of map element, we construct a generic function model into which we insert the elements of the specific data type. Then, the position, the scenario-dependent weight of each specific map element and the ego position define the map element functions that are all summed up to define the objective function of the corresponding maximization problem.

Representation of the input data
Let the disjoint union X = X M∪ X C denote the set containing all map elements, where X M is the set of all sorts of landmarks, and X C is the set of all lanes in the digital map, i.e., X M ={y ∈ R 3 ; y is a landmark}, (2.1) We assume that the cardinality of X is finite and write X = {χ i ; i ∈ I }, where |I | = |X | < ∞ and (χ i ) i∈I enumerates the elements of X . In our setting, landmarks in X M are points, whereas lanes are C 2 -curves, indicated by the subscript C. We would like to stress that our model and our optimization method can be extended to the case where also curves are allowed as landmarks, but for the remainder of the paper, we only consider the case of points as landmarks.
As part of the scenario S, which means performing an emergency stop on the highway in our case, each element χ i is assigned a constant weight w(χ i , S) ∈ R which models the strength of the contribution of the particular map element to the objective function, i.e., χ i ∈ X → w(χ i , S) ∈ R. In practice, this choice is made based on the categories of the map elements which arise from the element properties in the digital map. They are exemplarily shown in Table 1. According to those categories, X C is partitioned into L ∈ N disjoint subsets, each containing lanes of equal weight analogous to the weighting on the basis of the map data properties, i.e., X C =˙ l=1,...,L C l with w(χ i , S) = w l for all χ i ∈ C l , l = 1, . . . , L. (2.2) Furthermore, the function w X 0 ,S : R 3 → R weights the distance to the ego position X 0 by w X 0 ,S (x) = ω S ( x − X 0 2 ) with a function ω S : R ≥0 → R.

Formulation of the objective function
We are targeting an objective function that models only the interaction between the lanes and the landmarks, i.e., curves of different categories should not influence each other. Accordingly, every function evaluation can be assigned to a specific category C l , l = 1, . . . , L from Equation (2.2). Since a category weight represents the maximum contribution of that category to the objective function, the objective function obtained by aggregating all curve evaluations of the same category shall not exceed this weight. Therefore, we define the objective function F S : R 3 → R by with the category functions f l : R 3 → R and the landmark function f M : R 3 → R that are given by The cut-off w l > 0 bounds the function value of f C l , which sums up all curve functions of category l. We point out that f C l (x) > w l only occurs near points where different curves are very close to each other. F S is composed of the leading curve category value by f l and the values for the landmarks by f M . Every function f χ i : R 3 → R models the contribution of a map element χ i ∈ X to F S . Since the f χ i are at least C 2 , the objective function is piecewise continuously differentiable depending on the active category and whether the cut-off w l in f l is active. In that case, the local derivative of f l vanishes, and the derivative of F S equals the derivative of f M . This piecewise differentiability is exploited in Sect. 3.2.3 for the optimization algorithm.
Since we perform the maximization in three dimensions to generate solutions that do not coincide with the given map data, it is required to develop a generic function model that extends weightings w(χ i , S) from map elements χ i to the map element functions f χ i . Even in the case where an optimization would be performed only along the curves contained in X C , such an extension would be required to obtain a global interplay between the different map elements in forming a global objective function. One idea could be to extend w(χ i , S) to any point x ∈ R 3 based on the distance between x and χ i . However, computing this distance for a fixed point x requires a global optimization along χ i , and, if χ i is a curve that is not a line segment, the distance function would in general be nonsmooth with respect to x. Due to these drawbacks, a radial function φ( We convolve the resulting weighted kernel function w X 0 ,S w(χ i , S)φ(χ i , S) with a Dirac-type measure (Brokate and Kersting 2015, chs. 3-5), (Kanwal 1998, ch. 1) supported on χ i to obtain the desired extension to the three-dimensional space. For general d ∈ N and p ∈ R d , the Dirac δ-function over (2.5) Note, that χ i is an abbreviation for χ i ([0, 1]) if χ i ∈ X C in this formula. We define the test function ξ x : χ i → R as the weighted kernel centered at a point x ∈ R 3 , i.e., ξ x (y) = w X 0 ,S (x)w(χ i , S)φ(χ i , S)(y − x). Therefore, the corresponding map element function looks as follows in the case of a landmark χ i = y ∈ X M : The · -placeholder marks the variable on which the function in the second entry of ·, · depends. Analogously, Equation (2.5) yields the following function for curves (2.7) The kernel function φ(χ i , S) must be monotonically decreasing to model diminishing influence of χ i with growing distance between x and χ i . The norm · in the integral expression depends on the structure of the kernel function and will be discussed later. Due to the convolution, the curve function f γ fails to equal its target value w(γ , S) near the ends of γ . Fortunately, this behavior cancels out when two curves with coinciding weights, or rather of the same category, are connected via a C 2 -transition. From the mathematical point of view, the sum of such two curve functions and a curve function of a longer merged curve are equal.
The construction of the map functions f χ i requires that each kernel function takes its maximum at the origin and vanishes with increasing distance. The probability density function of the multivariate Gaussian distribution (Vinga 2004), also referred to as the Gaussian function, of the form μ, , fulfills these requirements under the assumption that tiny function values numerically equal zero. We set μ = 0 and choose the covariance matrix The standard deviation σ S,k indicates the size of the support in x k -direction, k = 1, 2, 3, depending on χ i . In the following, we will abbreviate σ S,k by σ k . The global coordinate system and the ellipsoidal shape of the support of the Gaussian function are depicted in Fig. 5 in Sect. 4. φ(χ i , S) is supposed to be invariant under rotations of the coordinate system in the x 1 -x 2 -plane with equal values for σ 1 and σ 2 , so that the objective function does not depend on the orientation of the map. σ 3 is small, since the vehicle motion in x 3 -direction is restricted by the elevation profile of the street.
With a normalization factor c χ i ensuring that f χ i restricted on χ i approximately equals w(χ i , S) (disregarding the factor w X 0 ,S for the ego position X 0 ), the kernel functions are defined in the following way with ϕ(r ) = exp(−r 2 /2) for each x ∈ R 3 : (2.8) We use (2π) 3 | | to compensate the factor of the Gaussian function. The norm · −1 is induced by the scalar product x, y −1 := x T −1 y. For a landmark y ∈ X M , this constant simply equals 1, because of f y (y) = c y w X 0 ,S (y)w(γ , S). In the case of a curve γ ∈ X C , we choose c γ = 1/ √ 2π based on the following considerations: We specify the norm · in Equation (2.7) with · −1 to remove the dependencies from the curve direction in such a way that the actual function value does not depend on the orientation of the curve, but mainly on the distance to the curve. If we approximate the curve γ locally with a straight lineγ = vt + y, t ∈ (−∞, ∞), y ∈ , cγ = c γ , then (2.9) Hence, we are able to state the contribution of the landmarks and the curves to the objective function.

Optimization approaches for solving for the problem formulation
The maxima of the derived objective function F S in Equations (2.3) and (2.4) represent points of interest for the given scenario S, i.e., the optimal locations for an emergency stop in our example. We determine them by solving the maximization problem where X ⊂ R 3 is a closed box (Cartesian product of closed intervals) that contains all relevant points. We target at methods that provide certificates of global optimality within a specifiable tolerance, as already mentioned in Sect. 1. Since this problem is nonsmooth, we will reformulate it to be able to compute higher-order derivatives. We need those derivatives to formulate optimality criteria and to apply higher-order optimization methods, which reduce the number of function evaluations significantly. The maximization problem is decomposed into L subproblems of the form: Here, X l ⊂ X is a large enough closed box containing all relevant map elements that are necessary for the definition of F l and in particular all maximizers of F l . We resolve the nonsmoothness of min{w l , ·} in Problem (3.2) by reformulating it as an equivalent constrained optimization problem with an auxiliary variable u l . For l = 1, . . . , L, we obtain max x∈X l ,u l ∈R Hence, the global solution set X * of Problem (3.1) is a subset of l=1,...,L X * l , where X * l ⊂ R 3 is the solution set of the l-th subproblem (3.3). X * consists of all x that fulfill F S (x) ≤ F S (y) for all y ∈ l=1,...,L X * l . The subproblems are independent of each other, i.e., they can be solved in parallel to reduce the total run time. The number of subproblems is induced by the number of categories C l , i.e., based on the road model complexity.
As we have seen in Sect. 2.2, Problem (3.3) is nonlinear and nonconvex, but sufficiently smooth to apply second-order optimization methods (Bertsekas 1999, chs. 1, 4), since we use smooth kernel functions and C 2 -curves. Due to the nonconvexity of these L subproblems, ascent methods (Boyd and Vandenberghe 2004, ch. 9) are in general not able to determine the global optima, rather they provide only local solutions that depend on the initial point. In fact, we need a measure for approximate global optimality to obtain certified solutions that contribute to the rigorous algorithm. We develop a hybrid optimization method specially tailored to the modeled problem that combines a rigorous deterministic branch-and-bound global maximization method with local approaches for selected boxes. For each of the L categories, we globally solve Problem (3.2), or Problem (3.3), to a prescribed accuracy. The following theory is used to achieve a rigorous branch-and-bound algorithm.

Interval arithmetic
Standard interval arithmetic (Hickey et al. 2001;Hansen 1975), (Moore et al. 2009, chs. 2, 4, 5), (Neumaier 1991, ch. 1) is used to obtain a box y that is an enclosure, i.e., y ⊃ g(x), for the image of a box x under a map g.
The principal goal of interval arithmetic is to evaluate a function on a whole interval x at once without having to evaluate the function at every single point x ∈ x. This is done by computing an interval (the tighter the better) that contains the image of x under this function. Hence, we build an interval extension of our objective function and, where required, of its derivatives. However, it must be noted that this approach in general leads to an overestimation of the interval-based function value since possible dependencies between intervals are usually lost.
Given a function g : refers to the output of interval arithmetic when we evaluate g for the box x, i.e., [g] is an interval extension for g. For instance, interval extensions of + and · yield for x = [−1, 1]: In comparison, the optimal enclosure (or interval hull, i.e., the smallest containing interval) for g(x) = · +(·) 2 (x) would be . We see that, in general, we can only expect promising results from interval arithmetic if the sequence of interval calculations is arranged and processed in an attentive way. Due to the form of the objective function, we compute several interval extensions of all map element functions f χ i in order to derive an expression for [F l ]. In particular, defining suitable inclusion functions is a non-trivial task because of the convolution of the map elements with the Gaussian kernel function, and, therefore, The monotonicity of ϕ μ,σ and ϕ μ,σ , which help to derive tight bounds for the interval-based evaluation, are deduced from ϕ μ, Therefore, we are able to calculate optimal enclosures for φ( , and for the interval evaluation of ∇φ(χ i , S). With the constant η 1 := exp(−1/2)/( √ 2πσ 2 ), the one-dimensional interval-valued Gaussian function and its derivative have the following form for z = [z L , z U ]: (3.5) The exact evaluation of the interval Gaussian function also provides optimal enclosures for the interval extensions of the landmark functions f y , y ∈ X M . Unfortunately, we encounter widely overestimating interval extensions [ f γ ] for curves γ ∈ X C due to the following reason: The integration along a curve γ in Equation (2.7) is performed with a quadrature formula. We apply the trapezoidal rule, because higher-order Newton-Cotes formulas (Davis and Rabinowitz 1984) do not provide better results for our purposes. Hence, we define a set Q : between the summands gets lost, and the interval sum is expansive. In general, the larger N Q , i.e., the more quadrature points are used to cover γ , the larger is the size of the output interval. Upper-bounding [ f C l ](x) by w l avoids inaccurate upper limits for the interval enclosure, since for The cut-off with w l also removes small oscillations of f C l along its ridges (cf. Fig. 1A) and at the same time greatly improves the quality of interval extension of min{ f C l , w l } compared to the extension of f C l .

INT: Interval optimization approach
With this interval arithmetic knowledge, we will specify our rigorous branch-andbound maximization algorithm based on the considerations for interval optimization (Hansen 1980), (Hansen and Walster 2004, chs. 7-8, 12-13), convex overestimation (Adjiman et al. 1998;Adjiman and Floudas 1996), and quadratic approximation (Burer and Letchford 2009;De Angelis et al. 1997). The basic methodology of the branch-andbound algorithm and the advanced bounding methods will be adapted to the emergency stop scenario by using the Problem (3.2) as well as the smoother reformulation in Problem (3.3).
After an initialization step, the branching and bounding strategies are repeatedly applied until a stopping criterion is fulfilled. The general pattern for the interval algorithm INT, with various refinements that are described later, has the following form:

Initialization
Let x init = X l be the initial interval box containing the feasible set of Problem (3.2) for the l-th category, thus . Instead of computing a reference value assumed by F l on x init , we perform a simplified search restricted to X in Sect. 3.4 that approximates the maximum of F S on X before investigating Problem (3.2) for all categories. Letf best,l denote the current best value for category l that is initialized with the result of this lower-dimensional maximization on X . We define B := {x init } to be the set of active interval boxes which remain to be investigated.

Branching
In every iteration, we select a box Then, x current is subdivided by multisection (Csallner et al. 2000) into four subboxes {x j sub } j=1,...,4 through its midpoint along the two directions with the longest side lengths so that and reference values sub for all j = 1, . . . , 4, and update the current best function valuef best,l iff x j sub >f best,l for some j. For efficiency reasons, we compute the reference value only when the box is used for further considerations, i.e., it holds x current = x j sub .

Bounding
is smaller thanf best,l for every point x ∈ x, the box x clearly cannot contain a global optimum. Apart from this derivative-free bounding, we also apply a gradient check (Hansen and Walster 2004, sec. 12.4). The interval gradient is used to identify whether F l is monotone along at least one coordinate direction in x, so that there is no local solution in the interior of the box and we can remove it from B. Since F l is not differentiable everywhere in x init due to the upper cut-off, special attention must be paid to the construction of the interval gradient [∇ F l ](x current ). For holds. Therefore, the gradient of F l has the following form for all x not lying on the boundary of x init or on the boundary of the cutting area R l ( For x ∈ ∂ R l (x current ), F l is in general only directionally differentiable. In order to combine the gradient criterion with the different expressions for ∇ F l , the interval gradient check proceeds with (3.8) The gradient criterion must be fulfilled for both interval extensions [∇ F l ] to be applicable without explicit knowledge of R l (x current ). We compute those gradient interval enclosures in the branching step for every new box resulting from the subdivision.
Since the boxes in B always cover the set of global optima X * l , the union of the boxes in B converges monotonically decreasing (in terms of set inclusion) to X * l when it is ensured that the maximum size of the boxes, i.e., the length of the longest edge of all boxes, tends to zero. The gap between the best function valuef best,l and the upper bound of its interval function value defines a certificate for global optimality for a box. In the algorithm, we determine the global optimal value only up to a guaranteed accuracy f , i.e., we stop to branch if this gap is smaller than f . Furthermore, we also stop branching if the box size is smaller than a predefined minimum box size >0 , because the limited accuracies for the vehicle localization and the map data make further box examinations needless. Finally, the category l * withf best := f best,l * ≥f best,l for all l = 1, . . . , L, provides the global optimal valuef best and the global solution box set.

Two local second-order methods
We consider the set of active boxes B after several steps of the interval optimization algorithm. Since the interval extensions [ f γ ] of the curve functions f γ significantly overestimate the interval enclosure for large boxes x, the upper limit of [ f l ](x) often matches w l , i.e., the upper bounds do not improve by further branching, and the number of removed boxes decreases. Therefore, we will investigate two different second-order methods, αBB and a quadratic approximation approach which we apply for boxes that fulfill specific subsequently introduced criteria. Hence, we derive the algorithms INT-αBB (INTerval αBB) and INT-QUAP (INTerval QUadratic APproximation).

INT-˛BB: An˛−based Branch-and-Bound Method
The idea of αBB (Adjiman et al. 1998;Adjiman and Floudas 1996;Meyer and Floudas 2005) is to compute locally valid concave overestimators of F l using second-order curvature information to generate better upper bounds than the interval extension [F l ] that we use in Sect. 3.2. The αBB algorithm is summarized in Algorithm 3.1. F l is regularized by a nonnegative quadratic term, weighted by α ≥ 0, to construct a local If α is large enough, then the Hessian H L x (x) = H F l (x) − α I is negative semi-definite and, therefore, L x is concave. The maximum value max x∈x L x can be exactly computed by local searches for the concave function, e.g., with SQP methods (Boggs and Tolle 1995). This maximum value, just as f U x , is an upper bound for F l (x). In our case, F l is not smooth enough to compute such an α. For this purpose, we present a method that combines αBB with the cut-off for f C l to derive a concave overestimator of F l in x.
Instead of a single α-value, we compute a vector (α 1 , . . . , α p ) ∈ R p ≥0 of α-values with p = |X M ∪ C l | to obtain individual concave overestimators L x χ i for the functions f χ i in order to define a concave overestimator L x l of F l and to compute the corresponding maximum value. We define (3.9) (cf. Maranas and Floudas 1994) An example of αBB is visualized in Fig. 3. The function f χ i is overestimated by the concave function L x (3.10) We use the maximum difference L x l − F l to decide whether it is promising to solve this problem. Since , and, therefore, with Equation (3.9) the maximum gap between F l and L x l in x is bounded by   (3.9)), which are symmetric matrices with interval enclosures in the entries. However, the interval Hessians [H f γ ](x) for the curves γ ∈ X C are usually widely overestimated for the same reasons as explained for the interval-valued function evaluation in Sect. 3.1. As a remedy, we replace w X 0 ,S by its maximum value w max X 0 ,S (x) := max x∈x w X 0 ,S (x) on x to reduce the complexity of the function. This modified curve function overestimates f γ slightly by the function since the slope of w X 0 ,S in x is negligible in small boxes that qualify for αBB. Furthermore, we can expect the size of those boxes to be smaller than the size of the kernel function support, which is roughly comparable to the width of a lane. In combination with our highway scenario, we assume that, in this case, the curve γ is approximately a straight lineγ . We derive the curvature of the corresponding function fγ in Appendix A. It depends only on the distance toγ , and, therefore, it is sufficient to determine [· − proj · −1 (·,γ ) −1 ](x), which contains all −1 -distances between the points in x and their projections onto the lineγ . If we insert this enclosure into the Algorithm 3.1 αBB: Local concave overestimation of the objective function Input: box x ∈ B, function F l : R 3 → R 1 determine the values (α 1 , . . . , α p ) for all map elements and define the concave overestimator L x l : x → R using Equations (3.9), (3.12) and Appendix A 2 compute a global maximum x * = arg max x∈x L x l (x) with the SQP method 3 Check box status with possible optimal value update forf best,l (Section 3.2.1): 4 if L x l (x * ) <f best,l then 5 remove x from B 6 else 7 if F l (x * ) >f best,l then 8f best,l := F l (x * ), x best,l := x * 9 end 10 end result of Appendix A, we get an interval extension for the maximum eigenvalue of Hfγ in x, i.e., an upper bound for the eigenvalues of Hfγ (x). In fact, an optimal interval enclosure for the second derivative of the Gaussian function, analogous to Equation (3.5), in Equation (A.7) provides a tight bound for αγ = max u 2 =1, x∈x u T Hfγ (x)u. However, this is only a tight bound for an overestimatorfγ while, strictly speaking, we are looking for a bound α γ for f γ .

INT-QUAP: local quadratic approximation approach
In the following, we will describe an alternative method beside the αBB theory to extend the interval algorithm. The local linearization of a curve γ is used again, but the objective function is locally approximated instead overestimated. We approximate the objective function F l locally in a box x by a quadratic function Q x l : x → R. The maximum value of Q x l is used to improve the upper bound provided by [F l ](x), although we have to add the approximation error R x l to this value. In the following, we will construct the quadratic approximation, determine the remainder term and present the computation of their upper bounds.
The Taylor expansion at a point x 0 = (x 0,1 , x 0,2 , x 0,3 ) T ∈ x provides quadratic approximations Q x for some θ = θ(x, x 0 ) ∈ [0, 1] and a multi-index ν = (ν 1 , ν 2 , ν 3 ), |ν| = ν 1 + ν 2 + ν 3 , ν! = ν 1 !ν 2 !ν 3 !. Then, we can write F l (x) = min{Q x l, (3.15) Each argument in this min-representation of F l contains a quadratic function and a remainder term, and, therefore, an upper-bounding approach for each of these two arguments provides an upper bound for F l in x by In practice, we compute the exact maximum value of the quadratic function Q x l,m in this formula, but we determine only interval enclosures of all the remainders R x χ i (x; x 0 ) to get an upper bound for the remainder R x l,m . However, we can use the cut-off min{w l , [ f C l ](x)} =: [a, b] to reduce the computational effort. If the upper bound w l for f C l is not active anywhere in x, i.e., b < w l , On the other hand, if the upper bound w l is active for every x ∈ x, i.e., a = b = w l , then f l = w l in x, and If only b equals w l , then we cannot determine whether or where the cut-off is active in x. Therefore, we compute the right-hand side of Equation (3.16) for both m = 1 and m = 2, and we choose the smaller value as an upper bound for F l (x). An example for this case is visualized in Fig. 4, which shows the partial cut-off in x and the resulting versions of the Taylor approximation. It should be noted that this figure gives just a qualitative description, since the approximation error here is too large for real application.
This quadratic approximation approach is also only performed on sufficiently small boxes (cf. Sect. 3.3.1), i.e., if the approximation error in Equation (3.16) is sufficiently small. Analogously to the αBB method, we simplify our problem by using the locally valid functionfγ (cf. Equation (A.3)) with the linearized curveγ , since widely overestimates the true enclosure in general. If we apply ∂ from the first line of Equation (A.4) and we recap v −1 = 1, then the third direc- (3.17) Fig. 4 A The function F l is bounded above by f C l + f M . They only differ inside the red surrounded area, where the cut-off with w l is active. B The functions w l + f M and f C l + f M are locally approximated by the quadratic Taylor expansions Q x l,1 and Q x l,2 We maximize these quadratic functions and add the approximation error to get an upper bound for F l in x Inserting these formulas into Equation (3.14) provides an integral-free interval extension [R x γ ](x, x 0 ) with θ = [0, 1] and an estimate for [ · −proj · −1 (·,γ ) −1 ](x). Then, the second term in the objective function of Problem (3.16) is bounded by where this choice of the interval extension [R x l ](x) for the Taylor remainder depends on the explicit form of the approximation that we discussed in Equation (3.15). It remains to show how to maximize the quadratic function Q x l . In general, the quadratic optimization problems with box constraints of the form  (Banchoff 1990, ch. 4). For a box x = [x L , x U ] ⊂ R 3 , this subdivision is given by (3.20)

Algorithm 3.2 Local quadratic approximation of the objective function
Input: box x ∈ B, function F l : R 3 → R 1 examine F l at the midpoint x 0 ∈ x, define quadratic function Q x l with Equations (3.13) and (3.15) 2 compute global maximum x * = arg max x∈x Q x l (x) with Equation (3.20) 3 compute an enclosure [R x l ](x) of the approximation error by Equation (3.18) 4 Check box status with possible optimal value update forf best,l (Section 3.2.1): Even though this approach is inefficient for higher dimensions, it is computationally feasible for d = 3. This completes our description of an upper-bounding approach for the objective function F l on x based on maximizing quadratic Taylor polynomials and bounding the remainder terms. The resulting method is summarized in Algorithm 3.2. In the worst case, we must apply this optimization method two times per box to maximize the quadratic functions Q x l,1 and Q x l,2 . Analogous to the concave overestimator in αBB, we could merge these quadratic problems and obtain a quadratically constrained quadratic program of the form max x∈R 3 ,u l ∈R (3.21) However, this problem cannot be solved with the previous box-constrained approach. Polynomial optimization (Lasserre 2001) was also considered, but it was significantly slower than our approach for this application.

1D optimization along the map data
Large values for the current best function valuef best are essential for the bounding step in the interval algorithm. Therefore, we perform a one-dimensional maximization along all curves γ ∈ X C to obtain an initial a priori reference value (i.e., a lower bound) Algorithm 3.3 One-dimensional maximization along the curves in X C to determine an initial value for the interval algorithm.
Input : map data X , scenario S Output: optimal point x best ∈ X C with function valuef best 1 foreach γ ∈ X C do 2 define function f 1D γ and the Problem (3.22) 3 apply extended TGO with optimal solution x best γ and valuef best,γ = f 1D γ (x best γ ) 4 end 5 determine γ * = arg max γ ∈X Cf best,γ and set for the interval algorithm. For every curve γ ∈ X C , we solve the one-dimensional problem ( 3.22) where f 1D γ approximates the function values along γ . Since the solutions of these problems are only initial values for the interval algorithm, they do not need to be solved rigorously. We use the idea of topographical global optimization (TGO) (Endres et al. 2018) with equidistant sampling, which compares the function values of neighboring samples to identify subintervals in [0, 1] where local optima must occur. We will extend this procedure by also computing the derivative d/dt f 1D γ (t) = w(γ , S)∇w X 0 ,S (γ (t)) T γ (t) + ∇ f M (γ (t)) T γ (t) at the samples.
Let {t 1 , . . . , t N ; t 1 = 0, t N = 1, t i < t j for i < j} be a set of N ∈ N sampling points. We choose a subset of the quadrature knots that are used for the evaluation of the curve functions. We define the set . . , N } for every γ ∈ X C . Note, that d f 1 and d f N are only one-sided derivatives. Our goal is to determine for every interval [t j , t j+1 ], j = 1, . . . , N − 1, whether it contains a local maximum. Therefore, we will state three criteria that prove the existence of local optima due to the differentiability (and the continuity) of f 1D γ : 1. The value at the right bound is at least as large as the value at the left bound and the function is non-increasing at the right bound, thus ( f j ≤ f j+1 ) ∧ (d f j+1 ≤ 0). 2. The value at the right bound is not larger than the value at the left bound and the function is non-decreasing at the left bound, thus ( f j ≥ f j+1 ) ∧ (d f j ≥ 0). 3. The function is non-decreasing at the left bound and non-increasing at the right bound, thus d f j ≥ 0 ∧ (d f j+1 ≤ 0).
On the boundary intervals [t 1 , t 2 ] and [t N −1 , t N ], we check the derivative-free criteria f 1 ≥ f 2 or f N −1 ≤ f N , respectively, to ensure the existence of local optima. For the intervals that demonstrably contain local solutions, we apply the algorithm by (Brent 1973, chs. 3-5). It performs a combination of the golden-section search (Kiefer 1953) and parabolic interpolation near the interval bounds.
However, it is necessary to compare f 1D γ (t) with F S (γ (t)) for all local solutions to overcome discrepancies between f γ and its approximation f 1D γ . The maximum of these optimal values and the values F S (y) for all landmarks y ∈ X M is chosen as an initial reference value for the full-dimensional interval approach.

Overall algorithms INT-˛BB and INT-QUAP
The steps of our global optimization algorithm are shown in Algorithm 3.4. We perform the one-dimensional maximization along the curve data X C in Line 5 before starting the interval algorithm according to Sect. 3.2. We either apply INT-αBB that utilizes concave overestimation, or INT-QUAP that incorporates the computation of quadratic approximations to extend the interval algorithm INT with second-order methods. Note, that we use two different rules with d αBB in Line 17 to decide whether we apply INT-αBB or INT-QUAP in x at all.
The design of the algorithm enables the solution to be computed by independent optimization problems for the different categories. Here, we do not exploit this parallelism, since this requires further effort for the synchronization to keep determinism. We use the optimal valuef best,l of category l as an initial value for the category l + 1 problem in Line 13, even though this creates dependencies between the categories. The solutions x best,l are compared in Line 34 to obtain the global optimal solution x best .
Interval Newton (Krawczyk and Neumaier 1986), interval constraint propagation (Kjøller et al. 2007), and non-diagonal shift matrices for αBB (Skjäl and Westerlund 2014;Skjäl et al. 2012) were also considered for the sake of completeness, but these approaches did not improve the performance of the overall algorithm for our problem.

Function evaluation improvement by curve linearization
An efficient evaluation of the objective function is decisive for a fast execution of the optimization. It strongly depends on the number of point and curve functions (cf. Eqs. (2.6) and (2.7)) that are derived for each map element. We use geometrical observations for the evaluation of these map element functions to decide in advance whether we can set the function value validly without an explicit evaluation. These geometrical checks reduce the computation time significantly for the curve evaluations in particular

Underestimating the distance to a curve
The functions f χ i in Equation (2.4) become infinitesimally small with increasing distance to the map elements χ i ∈ X due to the convolution with the Gaussian kernel function φ(χ i , S). With Equation (2.8), it can be shown that φ( In other words, if the distance between x ∈ R and y is larger than r φ,χ i for all y ∈ χ i , then f χ i (x) has no Algorithm 3.4 The extended interval algorithm. The αBB method for INT-αBB and the quadratic approximation for INT-QUAP are shown in Algorithms 3.1 and 3.2.
Input : map data X , scenario S (Section 2.1) Output: x best ∈ X * ⊂ R 3 global optimal point for scenario S 1 Modeling of Problem (3.3) for all categories based on scenario S (Section 2): 2 build C l with weights w l for all categories l = 1, . . . , L according to Equation (2.2) 3 determine the initial boxes X l containing the feasible sets of the subproblems in Equation (3.2) 4 define thresholds d , f > 0, x ∈ R 3 >0 (Sections 3.3.1 and 3.2.3) 5 call Algorithm 3.3: 1D optimization along all curves in X C (Section 3.4) 6 return reference point x best ∈ X C with function valuef best = F S (x best ) 7 if there exists a landmark y ∈ X M with F S (y) >f best then 8 determine y * = arg max y∈X M F S (y) 9 set x best = y * ,f best = F S (y * ) 10 end 11 Apply interval algorithm for every category: 12 for l = 1 to L do 13 set initial box set B = {x init } with X l ⊂ x init , x best,l = x best ,f best,l =f best and compute upper bound impact on the objective function and we can set the value to zero. It is straightforward to verify this condition for the map points. In the following, we present a method how to perform this relevance check for the map curves γ by using their B-spline representation.
We want to provide a simple lower bound for min t∈[τ j ,τ j+1 ] x − γ (t) −1 for every B-spline knot interval [τ j , τ j+1 ], j = 4, . . . , n − 3, to compare it with r φ,γ in order to decide whether or not to compute f γ (x) explicitly. Therefore, we approximate γ with a polygonal chain (or polyline) s γ using a set of vertices by (4.1) s γ interpolates γ and, hence, provides a better approximation than the polyline induced by the B-spline control points P γ of γ . Since This lower-bounding property is visualized in Fig. 5.
Here, the minimum distance d j s γ (x) between x and the line segment [S j , S j+1 ] := {s γ (t); t ∈ [τ j , τ j+1 ]} equals x − proj · −1 (x, [S j , S j+1 ]) −1 , given that the projection onto the line segment [S j , S j+1 ] has the form The second term on the righthand side of Equation (4.2) is bounded above by d j max , which is given in Equation (B.3) in Appendix B. In total, we decide to set and r φ,γ are independent of x and can be computed a priori.
Analogously, we need a lower bound for d j s γ (x) to apply that rule for a box x. An efficient computation of the exact minimum distance inf(d j s γ (x)) between a box and a line segment is demonstrated in (Schneider and Eberly 2002, sec. 10.9.4). Then, we set In particular, this result transmits to subsets of x, i.e., we can set f γ also zero for points in x and subboxes of x, which is a huge advantage in the branch-and-bound algorithm to reduce the computational effort. For every box, we store the active curve segments that contribute to the objective function.

Intersection between a box and a curve
Analogously, we can use geometrical tools to identify whether f γ takes its maximum value inside x. If the curve γ intersects the box x, then we set S) according to the line approximation in Equation (2.9). u x γ is an upper bound for max x∈x fγ (x), whereγ is a line approximation of γ in x. The difference between max x∈x f γ (x) and u x γ depends on the slope of w X 0 ,S , the curvature of γ in x, and to a largest extent on the fact whether the intersection of γ and x is only close to one of the ends of the curve, because f γ fails to approximate w(γ , S) there as discussed below Equation (2.7). We want to stress that we use this curve intersection approach especially for large boxes to avoid the numerical integration, whose computational effort is high due to the number of required quadrature points.
Since we determine only the intersection of the polyline s γ and x in practice (this is much easier than for B-spline curves), we define a subbox x sub ⊂ x so that γ intersects x if s γ intersects x sub . Consider again the line segment [S j , S j+1 ] and the curve segment γ | [τ j ,τ j+1 ] and assume thatx ∈ [S j , S j+1 ]∩x sub . Forȳ = proj · −1 (x, γ | [τ j ,τ j+1 ] ), it holds that |x k −ȳ k | ≤ max t k ∈[τ j ,τ j+1 ] |q jk (t k )| =: e k for every component k = 1, 2, 3 according to the considerations for Equation (B.3). Therefore, we define x sub = [x L + e, x U − e] with e = (e 1 , e 2 , e 3 ) T and x = [x L , x U ]. The side lengths of x must be larger than 2e so that x sub is well-defined and this box intersection criterion can be applied. This distance computation of the previous section with zero distance detects these intersections.

Prototypical implementation and evaluation
The optimization program has been implemented in MATLAB R2018b that is running on an Ubuntu 16.04 virtual machine with a four core processor and 16 GB RAM. The underlying hardware specifications are an Intel(R) Core(TM) i7-8850 H 2.60GHz processor with six cores and 32 GB RAM in total. Apart from MATLAB built-in functionalities, we use the CORA software package (Althoff and Grebenyuk 2017) that provides an implementation of interval arithmetic in MATLAB. In practice, basic arithmetic operations (addition, subtraction, multiplication, division) and elementary functions (cos, sin, exp,...) are overloaded such that they take intervals as inputs and map them to the smallest interval that contains the image of the input set (or, if this is too costly, to a reasonably tight upper estimate). Furthermore, we include C ++ libraries to MATLAB by using MEX files. The Geometric Tools Engine (Eberly 2014, chs. 5-7) helps to compute the distances between the boxes and the polyline segments. We also use the interval arithmetic functionality of the C-XSC library (Hofschuster and Krämer 2004) to accelerate the expensive interval Gaussian described in Sect. 3.2 by means of C ++ code. A single implementation of interval arithmetic is advised for production-ready code.
The fminbnd function of the MATLAB Optimization Toolbox performs the onedimensional maximization described in Sect. 3.4. Furthermore, we use the fmincon function to apply the SQP-algorithm for the concave overestimator. We limit the maximum number of steps to 20 due to the fast convergence. We choose the knot parameters in τ γ (cf. Sect. 2.1) for the quadrature to be piecewise equidistant in such a way that y m − y m+1 −1 ≈ 10 −1 for two neighboring quadrature points y m = γ (t m ), y m+1 = γ (t m+1 ) (cf. Sect. 3.1). Those parameters must be determined a priori for the continuity of f γ .

Framework for the emergency stop example
We define two examples for an emergency stop while driving on the highway to evaluate the developed algorithm and compare the results. For this scenario, we observe a short highway section of around five kilometers length on the Autobahn A9 near Ingolstadt, Bavaria, Germany. This highway section reflects the usual complexity of a german highway, since there are long monotonous passages, but it also contains complex areas, e.g., highway entries and exits, as well as a service station 1 . The map excerpt and the map elements are visualized in Fig. 6. The evaluation is applied for two distinct ego vehicle positions X 1 0 in northern direction, this example is referred to as Ego 1, and X 2 0 in southern direction of the highway 2 , this example is referred to as Ego 2. Apart from the lanes along the route of the vehicle, we also reward emergency call boxes and kilometer posts since they help to call for the breakdown service or to localize the emergency stop. In total, there are 76 lane objects, which are split into 11 categories and 21 landmarks for our example. We perform a further selection of the relevant map data based on the ego position, since, e.g., opposite lanes or road sections behind the vehicle are not considered. Therefore, we reduce the map data to 12 lanes (8 categories) and six landmarks corresponding for Ego 1 and 23 lane objects (10 categories) and six landmarks for Ego 2.
For an explicit formulation of the objective function, we need to define a weighting for the curve categories and the landmark types, as well as for the variances (σ 2 1 , σ 2 2 , σ 2 3 ) of the Gaussian functions φ(χ i , S). Here, we have tried to find a plausible weighting visualized in Table 1, but explicit rules for this weighting need to be established for a real application.
We determine the ego position factor w X 0 ,S (x) = w δ ( x − X 0 2 ) by the decreasing radial function w δ : R ≥0 → R that is defined by the positive, scenario-dependent constant δ = δ(S) > 0: The smaller δ is, the faster the contribution of the map data to the objective function decays with increasing distance from the ego position. We set δ X 1 0 = 5000 for a larger relevant area around X 1 0 and δ X 2 0 = 1000 for a smaller surrounding of X 2 0 . Furthermore, we choose f = 10 −3 for the maximum error of the guaranteed optimal function value, x = 0.05 for the minimum box side length and d = 10 −1 (αBB), respectively d = 5 · 10 −1 (quadratic approximation) for the threshold that activates the second-order method. Table 2 show the results of INT-αBB and INT-QUAP for the first example Ego 1. For every category, the weight w l and the number of active curves |C l | are shown, as well as the number of iterations and boxes. The tables A and B also depict the number of boxes that fulfill the box intersection criterion of Sect. 4.2, the gradient monotonicity criterion and the second-order method criterion to apply concave overestimation or quadratic approximation in Algorithm 3.4. Furthermore, the number of real-valued and 1 The UTM 32U coordinates in meters (x 1 -and x 2 -component) and heights in meters above mean sea level (x 3 -component) of the boundaries to the relevant box are [(680848, 5408861, 424)   We start the algorithm with the category of the highest weight w l , because this category is most likely to provide the global optimal solution. Since the result for the first category already returns the global optimum here, the interval algorithms for the remaining categories stop after a few steps. INT-αBB performs better than INT-QUAP, since a closer examination reveals that αBB always removes irrelevant boxes, whereas there remain active boxes by quadratic approximation that do not contain a global optimum. Since the higher-order approaches are expensive, those unsuccessful quadratic approximations affect the run time negatively.

Example Ego 1
The maximization for the leading category is visualized in Fig. 7. We observe all the boxes of the interval algorithm with αBB in Fig. 7A. The curves that belong to the specific category are marked in pink. In Fig. 7B are all boxes colored in green for which the αBB method is applied. In practice, these boxes must be sufficiently small. The optimal solution is close to an emergency lane, with a slight offset towards the emergency call box as the weighting of the map elements has already suggested. As the reader can imagine, this is a valid emergency stop.

Example Ego 2
For the second example Ego 2 with the vehicle position X 2 0 , the computation of the global maximum results in an emergency stop that is close to the next kilometer sign in the driving direction. Note that δ X 2 0 has been reduced compared to the first example. We use this example to compare the performance of our algorithm to existing global optimization methods. The results are shown in Table 3.
Apart from our algorithms INT-αBB and INT-QUAP, we examine INT without any second-order methods (third method) and the αBB algorithm which computes concave estimators for every observed box (fourth method). We must only adjust the parameter d in our implementation to enforce these borderline cases. Furthermore, we apply several global optimization algorithms implemented in MATLAB by handing over the objective function as a black box. These include MCS (Huyer and Neumaier 1999) and DIRECT (Finkel 2003) (cf. Sect. 1.2), as well as the GOP algorithm (Pál and Csendes 2009), which uses the INTLAB software package (Rump 1999). We apply this advanced interval algorithm to get a comparison for the performance of a state of the art implementation without handing over specific properties of our objective function. Furthermore, we compare the deterministic algorithms to the stochastic algorithms GLOBAL (Csendes et al. 2008) and CMA-ES (Hansen and Ostermeier 2001). Table 3 presents the number of function evaluations (as detailed as in Table 2) and the run time of the different methods, as well as the returned optimal value (the average values as well as the appearing optimal value ranges). The plots in Fig. 8 depict the ranges for the number of total function evaluations and for the resulting optimal values, which are especially interesting for the stochastic algorithms. Every dot corresponds to one of the 100 runs that we perform for every of the compared algorithms (only five runs of GOP). We can see that INT-αBB is the overall winner. This hybrid algorithm is almost twice as fast as our interval algorithm INT without αBB enhancement. The large number of pointwise function evaluations in pure αBB results from the maximization of the overestimators. In this example, INT-QUAP is even a bit worse than INT, which does not use second-order approaches. It creates less boxes but takes quite some effort in the quadratic approximation part. However, we also have investigated this example with f = 10 −2 (rather than f = 10 −3 ). Then, without reporting details here, INT-αBB and INT-QUAP process approximately the same number of boxes and, in these cases, INT-QUAP is slightly faster, since solving the quadratic problem with box constraints is cheaper than maximizing the concave overestimator. Although computing more boxes, INT required roughly the same run time as INT-αBB. GOP fails here with regard to the run time, although we even provided two function handles, one for the pointwise and one for the interval-based evaluation. Hence, we see the benefit of our customized implementation for this problem. The runtime of MCS is comparable to our quadratic approximation interval approach, but the optimal value is a bit worse and achieving a better value would require a significantly higher runtime of MCS. We used X 2 0 as initial point. Applying MCS to various other examples showed that in some situations MCS does not find a sufficiently good approximate global solution within an acceptable run time. We made similar observations for the DIRECT algorithm. Both methods can have difficulties to find the narrow ridges around the map elements where good candidates for global solutions are located. Also, they do not offer a certificate for approximate global optimality and thus the stopping criteria are based on other conditions that do not guarantee that always a sufficiently good solutions is returned.
The stochastic methods GLOBAL and CMA-ES are suitable alternatives to the interval αBB method in terms of the computation time, but they are not deterministic and, therefore, the range of the returned optimal values is significantly larger. The plots in Fig. 8 show their ranges of the number of function evaluations and of the optimal values. CMA-ES performs fastest and it is also able to return the global optimal value to the predefined accuracy. However, this stochastic algorithm is not able to provide the true optimal value in all of the test runs (cf. last row in Fig. 8). Furthermore, this result has only been reached for an explicitly chosen covariance matrix . When we use this choice for Example Ego 1, then CMA-ES misses the global optimal value by more than 0.5. In total, those algorithms that only use pointwise evaluations to find the global optima have limited success here. The support of the objective function is small compared to the search space and the solutions are close to the map elements. The large gradients of the curve functions f γ impede the success of algorithms that principally use single samples to infer the behavior of the function, e.g., MCS or DIRECT.

Conclusion and outlook
We developed a novel optimization-based approach for deriving secondary information from high definition digital maps based on a given scenario S of lanes, landmarks and weights. Thereby, new points of interest are generated by solving a global maximization problem derived from the scenario. A new rigorous deterministic global maximization algorithm was developed. It combines and, in several aspects, extends state-of-the-art approaches, enhancing them by exploiting specific properties of our problem. The resulting algorithm achieves a performance that compares favorably with other global optimization methods.
We proposed two variants of our method that both build on an interval algorithm. INT-αBB outperforms other deterministic algorithms and shows a performance that is competitive with non-deterministic approaches (cf . Table 3). However, the latter do not provide certificates of global optimality and they would not meet the safety requirements demanded for autonomous driving applications. Although INT-QUAP cannot fully compete with INT-αBB here, we have observed scenarios with an increased tolerance f for the optimal function value, where both algorithms examine roughly the same number of boxes. In that case, the advantages of INT-QUAP by solving boxconstrained quadratic problems result in less run times. For the future, higher-order polynomials could reduce the approximation error, but need more elaborate methods for the optimization with box constraints.
The local linearization of the curves has been fundamental to apply second-order techniques. This approach is valid for the highway scenario, since sufficiently small curvatures are guaranteed there. We plan to extend our results to urban scenarios, which are expected to be more complex and challenging due to the increased information density. On the other hand, this can be partly alleviated as the foresight could be significantly reduced. Therefore, we aim for investigating the validity of the local curve linearization for urban maps. As the suggested optimization approach can be highly parallelized, we plan to exploit this in future work. We also will focus on specifically tailored data structures for the management of the quadrature points.