A Largest Empty Hypersphere Metaheuristic for Robust Optimisation with Implementation Uncertainty

We consider box-constrained robust optimisation problems with implementation uncertainty. In this setting, the solution that a decision maker wants to implement may become perturbed. The aim is to find a solution that optimises the worst possible performance over all possible perturbances. Previously, only few generic search methods have been developed for this setting. We introduce a new approach for a global search, based on placing a largest empty hypersphere. We do not assume any knowledge on the structure of the original objective function, making this approach also viable for simulation-optimisation settings. In computational experiments we demonstrate a strong performance of our approach in comparison with state-of-the-art methods, which makes it possible to solve even high-dimensional problems.


Introduction
The use of models to support informed decision making is ubiquitous. However, the size and nature of the decision variable solution space, and the model runtime, may make a comprehensive -exhaustive or simply extensive -evaluation of the problem space computationally infeasible. In such cases an efficient approach is required to search for global optima.
Mathematical programs are one form of model that are explicitly formulated as optimisation problems, where the model representation imposes assumptions on the structure of the decision variable space and objective function. Such models are well suited to efficient solution, and identification of global optima may be theoretically guaranteed when feasible solutions exist. However many real-world problems are not suited to expression as a mathematical program (e.g., a solution is evaluated by using a simulation tool). From an optimisation perspective models where no assumptions are made about the model structure can be thought of as a black-box, where decision variables values are input and outputs generated for interpretation as an objective. In this case optimisation search techniques such as metaheuristics are required, i.e., general rule-based search techniques that can be applied to any model.
An additional widespread feature of many real-world problems is the consideration of uncertainty which may impact on model outputs, and so on corresponding objective function values. One strategy is to simply ignore any uncertainty and perform a standard search, possibly assessing and reporting on the sensitivity of the optimum after it has been identified. However it has been established that optimal solutions which are sensitive to parameter variations within known bounds of uncertainty may substantially degrade the optimum objective function value, meaning that solutions sought without explicitly taking account of uncertainty are susceptible to significant sub-optimality, see [BTGN09,GS16]. In the face of uncertainty the focus of attention for an optimisation analysis shifts from the identification of a solution that just performs well in the expected case, to a solution that performs well over a range of scenarios.
In this paper we develop a new algorithm for box-constrained robust black-box global optimisation problems taking account of implementation uncertainty, i.e., the solution that a decision maker wants to implement may be slightly perturbed in practice, and the aim is to find a solution that performs best under the worst-case perturbation. Our method is based on an exploration technique that uses largest empty hyperspheres (LEHs) to identify regions that can still contain improving robust solutions. In a computational study we compare our method with a local search approach from the literature (see [BNT10b]) and a standard particle swarm approach. We find that our approach considerably outperforms these methods, especially for higher-dimensional problems.
Structure of this paper. We begin with a review of the literature on metaheuristics for robust optimisation in Section 2 before outlining the formal description of robust min max problems in Section 3. We also consider some of the details of the established local robust search technique due to [BNT10b]. In Section 4 we introduce a novel approach, an exploration-focused movement through the search space identifying areas that are free of previously identified poor points. We include a discussion and descriptions of the algorithms used to identify empty regions of the decision variable search space. The approach is then tested against alternative heuristics in Section 5, on test problems of varying dimension. The experimental set up is described and the results of this analysis presented. Finally we summarise and consider further extensions of this work in Section 6.
2. Literature review 2.1. Robust optimisation Different approaches to model uncertainty in decision making problems have been explored in the literature. Within robust optimisation, a frequent distinction is made between parameter uncertainty (where the problem data is not known exactly) and implementation uncertainty (where a decision cannot be put into practice with full accuracy). Implementation uncertainty is also known as decision uncertainty [BTGN09,Tal09,BNT10b].
A common approach to the incorporation of uncertainty for blackbox problems is stochastic optimisation. Here knowledge of the probability distributions of the uncertain parameters is assumed and some statistical measure of the fitness of a solution assessed, e.g. using Monte Carlo simulation to estimate the expected fitness. This may be the expected value, or a more elaborate model such as the variance in the fitness of a solution, or even a multi-objective optimisation setting, see [PBJ06,dMB14].
An alternative to a stochastic approach is robust optimisation, whose modern form was first developed in [KY97] and [BTN98]. Whereas with stochastic optimisation a knowledge of probability distributions over all possible scenarios is typically assumed, in robust optimisation it is only assumed that some set is identified containing all possible uncertainty scenarios (potentially infinite in number). A classic robust approach is then to find a solution across all scenarios that is always feasible (strictly robust) and optimises its performance in the worst case. This is known as min max. For a given point in the decision variable space there is an 'inner' objective to identify the maximal (worst case) function value in the local uncertainty neighbourhood, and an overall 'outer' objective to identify the minimum such maximal value.
The field of robust optimisation has been primarily aligned with mathematical programming approaches. There the methodology is based around the definition of reasonable uncertainty sets and the reformulation of computationally tractable mathematical programming problems. For specific forms of convex optimisation problems, the problem incorporating uncertainty can be re-formulated to another tractable, convex problem, see [BNT07,GS10]. To overcome concerns that the strictly robust worst case approach may be overly conservative, the concept of robustness can be expanded in terms of both the uncertainty set considered and the robustness measure [GS16]. On the assumption that it is overly pessimistic to assume that all implementation errors take their worst value simultaneously [BS04] consider an approach where the uncertainty set is reduced, and a robust model defined where the optimal solution is required to remain feasible for uncertainty applied to only a subset of the decision variables at any given time. Min max regret, see [ABV09], is an alternative to min max, seeking to minimise the maximum deviation between the value of the solution and the optimal value of a scenario, over all scenarios. [BTBB10] considers soft robustness, which utilises a nested family of uncertainty sets. The distributionally robust optimisation approach, see [GS10], attempts to bridge robust and stochastic techniques by utilizing uncertainty defined as a family of probability distributions, seeking optimal solutions for the worst case probability distribution. [CG16] use a biobjective approach to balance average and worst case performance by simultaneously optimising both.
Robust optimisation in a mathematical programming context has 4 been application-driven, so considerable work has been undertaken in applying robustness techniques to specific problems or formulations, see [BS07,GS16]. There has also been some cross-over into the application of specific heuristics, for example see [GLT97,VMdCM11]. However application to general problems has been less well addressed [GS16]. Furthermore robust approaches applied to black-box models are much less widely considered than approaches for mathematical programming problems, see [MWPL13,GS16,MWPL16]. Recently, robust optimisation with implementation uncertainty has also been extended to multi-objective optimisation, see [EKS17].

Metaheuristic for robust optimisation
The min max approach has been tackled with standard metaheuristic techniques applied to both the inner maximisation and outer minimisation problems. In co-evolutionary approaches two populations (or swarms) evolve separately but are linked. The fitness of individuals in one group is informed by the performance of individuals in the other, see [CSZ09]. [Her99,Jen04] use such a two-population genetic algorithm (GA) approach, whilst [SK02,MKA11] consider two-swarm co-evolutionary particle swarm optimisation (PSO) techniques for min max problems. A brute force co-evolutionary approach is to employ complete inner maximisation searches to generate robust values for each individual in each generation of the outer minimisation, however this is expensive in terms of model runs (i.e., function evaluations), see [MWPL16]. More practical co-evolutionary approaches, for example using only small numbers of populations for the outer search and the inner (uncertainty) search which share information between populations from generation to generation, or following several generations, require the application of additional simplifications and assumptions, see [CSZ09,MKA11]. One general area of research is the use of emulation to reduce the potential burden of computational run times and the number of modelfunction evaluations, see [VDYL16]. [ZZ10] use a surrogate-assisted evolutionary algorithm to tackle the inner search for black-box min max problems. [MWPL13,MWPL16] employs Kriging meta-modelling coupled with an expected improvement (EI) metric, as well as a relaxation of the inner maximisation search. The EI metric is used to efficiently choose points in the decision variable space where nominal (expensive) function evaluation should be undertaken, see [JSW98], here with a view to most efficiently improving the estimate of the ro-bust global minimum. The relaxation involves iteratively performing the outer minimisation on a limited inner uncertainty neighbourhood followed by an extensive inner maximisation search in the region of the identified outer minimum. This continues whilst the inner search sufficiently deteriorates the outer solution, with the inner maximum point being added to the limited inner uncertainty set with each iteration.
A second approach due to [uRLvK14,uRL17] also uses Kriging and an EI metric, building on a meta-model of the expensive nominal problem by applying a robust analysis directly to the Kriging model and exploiting the fact that many more inexpensive function evaluations can be performed on this meta-model. A modified EI metric is calculated for the worst case cost function of the meta-model, to efficiently guide the search in the nominal expensive function space. In [uRL17] the approach is applied to a constrained non-convex 2 dimensional problem due [BNT10b,BNT10a], the unconstrained version of which is also considered here. The Kriging-based approach is shown to significantly outperform the approaches outlined here, in terms of the number of expensive function evaluations required to converge towards the robust optimum. In general we would expect the approach from [uRLvK14,uRL17] to outperform the approaches considered here, in terms of efficiency when applied to low dimensional non-convex problems. However the primary challenge with meta-model based approaches is their application to higher dimensional problems. The test cases considered in [MWPL13,MWPL16,uRLvK14,uRL17] have either been restricted to low dimensional non-convex problems, or simpler convex and convex-concave problems of up to 10 dimensions.
One local black-box min max approach is due to [BNT07,BNT10b,BNT10a]. Here a search is undertaken by iteratively moving along 'descent directions'. Uncertainty around individual points is assessed using local gradient ascents, based on which undesirable 'high cost points' (hcps) are identified. Steps are taken in directions which point away from these hcps, until no direction can be found.
Our approach is inspired by both elements of the descent directions technique and the concept of relaxation of the inner maximisation search. We extend the idea of locally moving away from identified hcps to a global perspective, seeking regions of the solution space currently empty of such undesirable points. Furthermore the nature of our outer approach enables the curtailing of an inner maximisation search if it is determined that the current point under consideration cannot improve on the current best robust global solution.

Problem description
We consider a general optimisation problem of the form where x x x = (x 1 , x 2 , . . . , x n ) T denotes the n-dimensional vector of decision variables, f : R n → R is the objective function, and X ⊆ R n is the set of feasible solutions. We write [n] := {1, . . . , n}. In this paper, we assume box constraints X = i∈[n] [l i , u i ]. Any other potential feasibility constraints are assumed to be ensured through a penalty in the objective.
In implementation uncertainty, we assume that a desired solution x x x might not be possible to put into practice with full accuracy. Instead, a "close" solutionx x x = x x x + ∆x x x may be realised. The aim is to find a robust x x x such that for any such solutionx x x from the neighbourhood of x x x, the worst case performance is optimised.
More formally, we follow the setting of [BNT10b] and consider the so-called uncertainty set where Γ > 0 defines the magnitude of the uncertainty, and · refers to the Euclidean norm. The worst case costs of a solution x x x ∈ X are then given as and so the robust optimisation problem is given by: We therefore have an inner maximisation and outer minimisation problem, such that the identification of the robust global optimum is based on finding the (outer) minimum worst case cost objective function value in the decision variable space, and that objective is determined by the (inner) maximisation of the nominal objective function in the uncertainty neighbourhood around each point in the decision variable space. This type of problem is also known as min max.
Note that x x x + ∆x x x may not be in X , for which reason we assume that the definition of f is not limited to X . However, if it is desired that x x x + ∆x x x ∈ X for all ∆x x x ∈ U, then this can be ensured by reducing the size of the feasible search space by Γ.
As an example for our problem setting, consider the 2-dimensional polynomial function due to [BNT10b]: f (x, y) =2x 6 − 12.2x 5 + 21.2x 4 + 6.2x − 6.4x 3 − 4.7x 2 − y 6 −11y 5 + 43.3y 4 − 10y − 74.8y 3 + 56.9y 2 − 4.1xy −0.1y 2 x 2 + 0.4y 2 x + 0.4x 2 y (poly2D) For a feasible solution space within bounds [−1, 4] in each dimension, and uncertainty defined by a Γ-radius value of 0.5, the nominal and worst case plots for (poly2D) are shown in Figure 1. In min max the problem is one of finding the global minimum for the worst case cost function. If uncertainty is ignored the problem is just one of finding the global minimum of the (nominal) objective as shown in Figure 1a, whereas including uncertainty the problem becomes one of finding the (worst case cost) objective as shown in Figure 1b. In both cases the search proceeds based on generating nominal objective values but for the worst case cost we must further undertake some assessment of the impact of uncertainty on those objective outputs.  Here the global optimum value for the nominal problem is -20.8 at (2.8, 4.0). The worst case plot is estimated by randomly sampling large numbers of points in the Γ-uncertainty neighbourhood around each plotted point. The worst case cost at each point is then approximated as the highest value of f (x) identified within each Γ-uncertainty neighbourhood. The global optimum for the worst case problem is approximately 4.3 at (-0.18, 0.29). The significant shift in the nominal versus robust optima, both in terms of its location and the optimum objective, emphasises the potential impact of considering uncertainty in decision variable values. The difference between the nominal and robust optimal objective function values is the 'price of robustness', see [BS04].

Local robust search using descent directions
We briefly summarise the local search approach for (ROP) that was developed in [BNT10b]. Here, (ROP) is solved using a local robust optimisation heuristic illustrated by Figure 2. An initial decision variable vectorx x x is randomly sampled. Then a series of gradient ascent searches are undertaken within the Γ-uncertainty neighbourhood of this candidate solution to identify hcps, see Figure 2a. This approximates the inner maximisation problem max ∆x x x f (x x x+∆x x x). Using a threshold value that is dynamically adjusted during the algorithm, a subset H(x x x) of all evaluated points is identified, see Figure 2b.
In the next step, a descent direction is identified that points away from the set H(x x x), see Figure 2c. To this end, a mathematical programming approach is used, minimising the angle between the hcps and the candidate solution. This leads to the following second order cone problem.
Here, d d d is the descent direction, which is normalised by Constraint (2). Constraints (3) ensure that β is the maximum angle between d d d and all high cost points h h h. Through Constraint (4), we require a feasible descent direction to point away from all points in H(x x x). When an optimal direction cannot be found, the algorithm stops -a robust minimum has been reached. Next the size of the step to be taken is calculated, see Figure 2d. A step size just large enough to ensure that all of the hcps are outside x (a) Candidate point x x x (centre), and points evaluated for the inner maximisation problem (blue). points.
x (c) A descent direction is identified by solving a second order cone problem.
x (d) The step size is determined. Figure 2: Description of the descent direction robust local search approach [BNT10b].
of the Γ-uncertainty neighbourhood of the next candidate solution is used. Using the identified descent direction and step size the algorithm moves to a new candidate point, and so the heuristic repeats iteratively until a robust minimum has been identified.

A new largest empty hypersphere approach 4.1. Algorithm overview
Building on the notion of a search that progresses locally by moving away from already identified poor (high cost) points, we develop a global approach that iteratively moves to the region of the decision variable solution space furthest away from recognised hcps. This is an exploration-focused approach, although rather than concentrating on examining unvisited regions the intention here is to identify and visit regions devoid of hcps. Assuming uncertainty as considered previously in terms of a single value Γ that defines a radius of uncertainty in all decision variables, we associate the idea of the largest empty region (empty of hcps) with the idea of the largest empty hypersphere (LEH), or largest empty circle in 2D. The approach is then to locate the next point in the search at the centre of the identified LEH, and to iteratively repeat this as more regions are visited and hcps identified. The approach is described in Figure 3. We start by randomly sampling one or more points and evaluating the objective function f at each. From these start points a candidate point is selected and an inner analysis undertaken in the candidate's Γ-uncertainty neighbourhood with a view to identifying the local maximum, Figure 3a. This local worst case cost for the candidate is the first estimate of a robust global minimum, that is a global min max, and is located at the candidate point. The aim is now to move to a point whose uncertainty neighbourhood has a lower worst case cost than the current global value. We seek to achieve this by identifying the largest hypersphere of radius at least Γ within the defined feasibility bounds which is completely empty of hcps, and moving to the centre of that LEH, see Figures 3b -3c.
All points evaluated are recorded in a history set, a subset of which forms the high cost set. The high cost set contains a record of all points evaluated so far with an objective function value greater or equal to a high cost threshold, and here the high cost threshold is set as the current estimate of the robust global minimum. Both the history set and the high cost set are updated as more points are visited and the high cost threshold reduces, see Figures 3d -3e. On performing all inner searches after the first candidate, a candidate's robust value may be no better than the current estimate of the robust global minimum (and therefore the current high cost threshold), in which case at least one point will be added to the high cost set. Alternatively if a candidate's robust value is better than the current estimate of the robust global minimum, this current recorded optimum is overwritten and the high cost threshold reduced accordingly. Again this introduces at least one high cost point to the high cost set, but the reducing threshold may also introduce additional points from the history set; this is suggested in Figure 3e.
The search stops when no LEH of radius greater than Γ exists or some pre-defined resource limit has been reached. Then the candidate point around which the current estimate of the robust global minimum has been determined is deemed the robust global minimum. Otherwise the search repeats, performing analysis in the Γ-uncertainty neighbourhood around candidates to estimate the local (inner) max, updating the global minimum worst case cost if appropriate, and moving to the next identified LEH, Figure 3f. (c) Identify the largest empty hypersphere, the centre of which is the next candidate point.
(d) Inner search around the new candidate. The robust value here is less than the current global minimum.
(e) The current high cost set, including more previously evaluated points due to the reduced high cost threshold.
(f) Identify the largest empty hypersphere, the centre of which is the next candidate point. The critical feature of such an approach is the identification of regions of the solution space that are currently empty of, and furthest from, the undesirable hcps. As defined here this corresponds to identifying the largest hypersphere devoid of all hcps.
Given a discrete history set H of all points evaluated so far, high cost points are those members of H with objective value which is at least the current high cost threshold τ , i.e., The identification of a point p p p ∈ X which is furthest from all N τ high cost points in H τ is a max min problem: where d(p p p, q q q) is the Euclidean distance between two points p p p and q q q, see [OS97].
In the following, we specify this general LEH approach by considering two aspects in more detail: The outer search is concerned with placing the next candidate point x x x by solving (LEHP). The inner search then evaluates this candidate by calculating g(x x x) approximately.

Outer search methods
Here we will introduce different approaches to identifying the largest empty hypersphere, given a set of high cost points H τ . It should be noted that none of these approaches requires additional function evaluations, which is usually considered the limiting resource in black-box settings.

Randomly sampled LEH algorithm
A very simple approach is to generate potential candidates randomly within the feasible region, then determine whether they are more than Γ away from all hcps. If so they are a valid candidate, if not resample up to some defined maximum number of times beyond which it is assumed that no such candidate point can be found and the solution has converged on a robust global minimum. Rather than being a largest empty hypersphere approach this is just a valid empty hypersphere approach, and the size of the identified empty hypersphere might vary considerably from one candidate to the next.

Genetic Algorithm for LEH
The solution of (LEHP) is an optimisation problem. Furthermore, given a point p p p which is a potential candidate for the centre of the largest empty hypersphere, the inner minimisation calculation in (LEHP) involves just an enumeration over the N τ Euclidean distance calculations between each hcp and p p p to identify the minimum distance d(p p p, h h h k ), where h h h k is the closest hcp. Therefore the focus for the solution of (LEHP) is the outer maximisation, for which we may consider an approximate heuristic approach. We employ a genetic algorithm (GA), a commonly cited evolutionary algorithm (EA) [Tal09]. Here each individual represents a point p p p in the decision variable space, and the objective function f LEH (p p p) := min h h h∈Hτ d(p p p, h h h) is the minimum distance between a given point p p p and all hcps in H τ . We seek to maximise this minimal distance by evolving a population of points starting from randomly selected feasible points in the decision variable space X . The best point generated by the GA is the next candidate point -that is estimated centre of the LEH, for the current H, τ and H τ .

Voronoi based LEH
Within the literature a widely referenced approach for tackling low dimensional LEH problems is due to [Tou83], and is based on the geometric Voronoi diagram approach, see [Cha93,OS97]. The Voronoi approach partitions a space into regions (cells). For a given set of points each cell corresponds to a single point such that no point in the cell is closer to any other point in the set. Points on the edges between cells are equidistant between the set points which lie on either side of that edge. For our LEH problem the set of points is H τ , and the Voronoi diagram approach corresponds to segmenting the feasible space X into N τ separate cells, one for each hcp. The (Voronoi) vertices that lie at the intersection of these cell (Voronoi) edges maximise the minimum distance to the nearby set points, see [Cha93,OS97]. So for a given H τ if we can determine the Voronoi diagram we can use the identified Voronoi vertices as potential candidate points p p p. The solution of (LEHP) is then simply a matter of enumeration, for each p p p calculating the (inner) minimum Euclidean distance to all hcps, and then selecting the (outer) maximum such minimal distance.
The original approach due to [Tou83] includes the identification of vertices (candidate centres of LEHs) that can be sited outside of defined boundaries, in infeasible regions. This is not exactly as required here.
To deal with this edges that cross feasibility boundaries are identified and the associated vertices which are outside of X are relocated to an appropriate point on the boundary of X . Here any coordinate i ∈ [n] of such an external vertex that is either less than l i or greater than u i is re-set to l i or u i as appropriate.
However the Voronoi approach has exponential dependence on n, as constructing the Voronoi diagram of N τ points requires O(N τ logN τ + N n/2 τ ) time [Cha93]. This suggests that such an approach in not computationally viable for anything other than low dimensional problems. On the basis that a Voronoi diagram based approach is the primary recognised heuristic for identifying the largest empty circle we will consider a Voronoi based robust LEH heuristic here only in the context that for 2D problems in our experimental analysis this approach will serve as a good direct comparator for our other robust LEH heuristics.

Inner search methods
Discussions of the LEH approach have so far focussed on the outer minimisation search, assuming some form of inner search that provides the inner robust maximum for each candidate point in the minimisation search. In [BNT10b] a two-stage gradient ascent search is recommended for each inner search around a candidate point. This assumes gradient information is available and proposes (n + 1) individual twostage gradient ascents for each candidate. For a 100-dimensional problem this would require several thousand function evaluations around each candidate point. In practical terms both the number of function evaluations required to undertake a global search and the requirement for gradient information may make such extensive inner searches prohibitive. Given, for example, budgetary restrictions on the number of function evaluations, some trade-off must be achieved between the extent of each inner Γ-radius uncertainty neighbourhood search and globally exploring the search space. But this trade-off between robustness in terms of the extent of the inner searches, and performance in terms of the outer global search, is complex, see [MLM15,DHX17]. For example the determination of an appropriate inner approach -type of search, extent of search and parameter settings -may be both instance (problem and dimension) dependent and dependent on the outer approach.
Here we do not propose to recommend a definitive inner search approach. From a theoretical point of view we assume the information is provided by some oracle. From an experimental point of view in the algorithm testing and comparisons below we assume the same basic inner Γ-radius uncertainty neighbourhood analysis for all heuristics, to ensure a consistency when comparing results for alternative search approaches.
There is, however, an aspect of our LEH approach that enables an additional feature, the forcing of an early end to an inner search. The LEH approach is exploration-led, the objective being to locate and move to the candidate point in the decision variable space furthest from all hcps. Hcps are designated based on the determination of a high cost threshold τ , set here as the current estimate of the robust global minimum (min max) value. The nature of this approach enables (inner) uncertainty neighbourhood searches around each candidate point to be restricted when appropriate. If an inner search identifies a local point with objective function value above τ the inner search can be immediately curtailed on the basis that the candidate is not distant from hcps. This equates to the recognition that the candidate point is not an improvement on the current estimated robust optima. Such regulating of inner searches has the potential to significantly reduce the number of function evaluations expended on local neighbourhood analysis. In the case of budgetary limitations on numbers of function evaluations this further enables more exploration of the decision variable space.

Algorithm summary
Given one of our three approaches to identifying the LEH devoid of hcps, random, GA or Voronoi, the overarching algorithm for the robust exploratory LEH heuristic is given in Algorithm 1. Here one of these three approaches to the outer search is applied in line 16 as LEH Calculator(H τ ), for a defined high cost set H τ . It is assumed that this routine will return a candidate point x x x LEH and an associated radius r LEH , that is the minimal distance between x x x LEH and all points in H τ . The heuristic will halt if r LEH is not greater than Γ.
For a defined number of initialisation points, random points in X are selected and the function f evaluated at these points. The points and their function evaluations are recorded in history sets H and F H , lines 1 -6. Having randomly selected a candidate point x x x c from H we perform an inner maximisation in the Γ-uncertainty neighbourhood around x x x c , see line 10. The description of the inner maximisation is given below as Algorithm 2. If this is the first candidate point, or the local robust value for this candidateg(x x x c ) is less than the current best solution τ , this minimum is updated and the associated global minimum point x x x Op replaced by x x x c , see lines 11 -14.
Next the high cost set H τ is established as all members of H with corresponding function values in F H that are greater than or equal to the current high cost threshold τ , see line 15. Based on H τ , the next candidate point is identified via one of the outer search approaches, see line 16. If the heuristic is halted at this stage due to an inability to identify a valid LEH or at any stage due to the budget being exceeded, the extant estimate for the robust global minimum x x x Op is returned. Here the point to be evaluated is determined by random sampling in the Γ-radius hypersphere centred on x x x c , line 11. Under other inner maximisation rules this would be determined by some explicit maximisation search heuristic. As the function is evaluated at the inner search points the local robust value (inner maximum) Local Robust is updated as appropriate, line 18. If Local Robust exceeds the high cost threshold τ the inner maximisation is immediately terminated, lines 19 -21. Algorithm 2 ends by returning an estimate for the worst case cost value at x x x c ,g(x x x c ) into Algorithm 1.

Example LEH application
In order to give some indication of the nature of our LEH search we have applied it to the 2-dimensional problem (poly2D) and plotted the points evaluated and associated search path of the current estimate of the robust global minimum in Figures 4e and 4f. Here the LEH Voronoi algorithm is used. For comparison we have also plotted corresponding results for two alternative heuristics, a robust Particle Swarm Optimisation (PSO) approach shown in Figures 4a and 4b, and the local descent directions approach from Section 3.2 shown in Figures 4c and 4d. Here the robust PSO is used as a proxy to a brute force or co-evolutionary approach. The basic global PSO formulations have been used, as described in [SE98]. The descent directions approach has been extended by using random re-starts, as a proxy to extending it to a global approach. In all cases inner random sampling in a hypersphere of 100 Γ-uncertainty neighbourhood points is used, and a maximum budget of 10,000 function evaluations employed.
The plots shown in Figure 4 are for only a single run of each heuristic, and as such should only be seen as exemplars intended to give some indication of the different natures of these outer search approaches. It can be seen that whilst the robust PSO explores the decision variable space somewhat, and the re-starting descent directions follows (exploits) a series of local paths, the LEH approach features both considerable exploration globally and more intense analysis of promising points. It is clear that the curtailing of the inner searches in the LEH approach enables much wider exploration for fewer function evaluations. In this example less than 1,000 function evaluations have been required before the LEH heuristic has stopped because an LEH of radius greater than Γ cannot be found, but for larger (dimensional) problems such stopping prior to reaching the budgetary limit will not apply. One striking feature of Figure 4e is how many of the inner searches stop immediately on the evaluation of a candidate point. This is because the objective value at these candidate points exceeds the current threshold τ .
The Voronoi based search exemplified by Figures 4e and 4f is a good indicator of the nature of the searches due to all three LEH approaches, random, GA and Voronoi. However the radii of the LEH identified for each candidate will vary with the use of each of these algorithms. Figure 9 in Appendix B gives some indication of how the radii of the hyperspheres generated by each of these LEH heuristics progress as the exploration proceeds.

Set up
In order to assess the effectiveness of the LEH approach the heuristic has been applied to eight test problems, and results compared against the two alternative search heuristics described in Section 5.2. Experiments have been performed on 2D, 4D, 7D, 10D and 100D instances of these test problems; results have also been generated for (poly2D). Both the genetic algorithm and random forms of the LEH heuristic have been assessed for all instances. The LEH Voronoi has additionally been applied to the 2D instances, with the intention of giving some indication of the differences due to a 'best' LEH identifier algorithm (Voronoi) versus the alternatives. All LEH approaches are initialised by randomly sampling a single point in X . Assuming that for most real-world problems the optimisation analysis will be limited by resources, a fixed budget of 10,000 function evaluations (model runs) is assumed. The same inner approach is employed for all heuristics. A simple random sampling in a hypersphere of 100 points in a point's local Γ-uncertainty neighbourhood is used for all instances, and the local robust maximum is estimated as the maximum due to this sampling. For the LEH approaches this inner sampling is curtailed if a point is identified in the uncertainty neighbourhood that has objective value exceeding the current high cost threshold τ .
All experiments have have been performed using Java, on an HP q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  Pavilion 15 Notebook laptop computer, with 64 bit operating system, an Intel Core i3-5010U, 2.10GHz processor, and 8GB RAM. Each heuristic search has been applied to each test problem-dimension instance 50 times to reduce variability. For the solution of the Second Order Cone Problem as part of the descent directions algorithm [BNT10b], the IBM ILOG CPLEX Optimization Studio V12.6.3 package is called from Java.

Comparator heuristics
Our experiments have been conducted on LEH, a re-starting descent directions, and robust PSO metaheuristics. We have applied parameter tuning to 3 of the 5 comparator heuristics -LEH Voronoi and LEH Random do no have tunable parameters -employing an evolutionary tuning approach using a genetic algorithm to generate a single set of parameters for each heuristic, for all test problems. For each of the 3 tuned heuristics the same subset of the test instances was used, running each member of an evolving population on each of these instances multiple times to generate mean result for each member of a population on each test instance. The performance of each individual in a population was ranked separately for each test instance, across the members of the population, leading to mean overall ranks which were used as the utility measure in tournament selection; see e.g. [ES12]. The effectiveness of the local descent directions approach [BNT10b] suggests that extending this to a global search by using random restarts will provide a reasonable comparator. A local descent directions search is undertaken from a random start point, and when this is complete it is repeated from another random start point. This is repeated until the function evaluations budget is reached. In descent directions a set of high cost points leads to the identification of an optimal stepping direction and step size, if a valid direction exists. However the algorithm includes a number of dynamically changing parameters which adapt the high cost set and enforce a minimum step size. Here we have tuned 5 parameters relating to these stages of the heuristic; see [BNT10b] for further information. Labelled 'd.d. Re' in the results section.
As a proxy to a brute force or co-evolutionary approach an outer particle swarm search is considered. The basic formulations for the global PSO approach have been used as described in [SE98] and 5 parameters have been tuned: swarm size, number of iterations, and for the velocity equation the C 1 and C 2 acceleration parameters and inertia weight parameter ω. The combined swarm size times number of iterations was limited to 100 in order to align with the budget of 10,000 function evaluations and the level of inner sampling. Labelled 'PSO' in the results section.
Our robust LEH metaheuristic is considered for the three alternative ways of identifying the largest hypersphere that is empty of hcps: • Randomly sampled valid empty hypersphere, see Section 4.2.1. This includes re-sampling up to 1,000 potential candidates in an attempt to identify a valid empty hypersphere, otherwise it is assumed that a valid point cannot be found and a robust global minimum has been reached. Labelled 'LEH Rnd' in the results section.
• Genetic algorithm LEH, see Section 4.2.2. Here we have tuned 6 parameters: the size of the population, number of generations, number of elites, tournament size, and mutation probability and size; we have fixed the use of tournament selection and the choice of midpoint crossover. The combined population size times number of generations was limited to 100, which is somewhat based on runtime considerations associated with the large value of N τ , the number of candidate points visited with a budget of 10,000 function evaluations. Labelled 'LEH GA' in the results section.
• Voronoi based [Tou83] LEH, see Section 4.2.3. Here the construction of the Voronoi diagram for the input points H τ is performed using the Java library due to [Nah17]. This generates geometric data, Voronoi vertices and edges, which are used to determine a set of potential candidate points -Voronoi vertices, including those originally outside of X relocated to the boundary of X -for the centre of the LEH. Labelled 'LEH Vor' in the results section.

Test functions
A large number of test functions are available for benchmarking optimisation algorithms, and posing a variety of difficulties, see [Kru12,JY13].
The full description of these eight test functions is given in Appendix A. To give some indication of the nature of these functions contour plots of the 2D instances are shown in Figure 5, for both the nominal and worst cases.

Results
Results of the 50 samples runs for each heuristic applied to each test problem-dimension instance are presented here. In each run the best solution as identified by the heuristic is used. However the points in the decision variable space that have been identified as best have robust values generated using the simple inner random sampling approach, with a budget of up to 100 sample points. To better approximate the true robust values at these points their robust values have been re-estimated based on randomly sampling a large number of points (nominally 1,000,000) in the Γ-uncertainty neighbourhood of the identified robust point. This is a post processing exercise and does not affect the min max search.
Mean results due to each set of 50 sample runs are shown in Tables 1 and 2. We have applied the Wilcoxon rank-sum test with 95% confidence to identify the statistically best approaches. Results highlighted in bold indicate the approaches that are statistically equivalent to the best one observed, for a given problem-dimension instance. Corresponding box plots, giving some indication of how the results are distributed across the 50 samples, are shown in Figures 6, 7 and 8. Additional results, the standard deviations due to each set of 50 sample runs, the average number of candidate points visited and average number of function evaluations undertaken, are given in Appendix C.
From Table 2 we see that for 100D instances the LEH GA approach is best for all test problems, and in several cases the mean LEH GA result is substantially better than all of the alternative heuristics. From Tables 1 and 2 the LEH approach is among the best in at least 6 of the instances for all other dimensions.
For 2D instances the LEH Voronoi approach is among the best results for 7 of the 9 problems, whilst LEH GA and LEH Rnd are each amongst        the best results for 3 and 2 problems respectively. It should also be noted that in 5 of the 7 instances where LEH Voronoi is among the best, LEH GA is either statistically equivalent or the mean value is second best. For the 2D Sphere instance d.d. Re is marginally better than LEH Voronoi, whilst d.d. Re and LEH heuristics are statistically equivalent for the (poly2D) and 2D Volcano and Rosenbrock instances. The robust PSO approach is statistically equivalent to LEH heuristics for the 2D Sawtooth instance. For the 4D -10D instances d.d. Re is statistically equivalent to LEH GA in the 4D Volcano problem and the 7D and 10D instances of the Rastrigin problem, and better than LEH GA for the Rosenbrock and Sphere problems. For the 4D Rosenbrock and Sphere problems the differences between d.d. Re and LEH GA are reasonably small, however in the 7D and 10D instances d.d. Re is substantially better. Considering the shape of the Rosenbrock and Sphere functions it can be expected that a local search will perform particularly well for these problems.
LEH GA is better than LEH Rnd for all instances excluding (poly2D). In a number of instances LEH GA is substantially better than LEH Rnd. The number of candidate points that LEH can visit is substantially increased by the early stopping of inner searches as soon as the high cost threshold is exceeded, see Tables 3 and 5 in Appendix C. Although this feature must unquestionably play a role in the success of the LEH GA approach, the fact that LEH Rnd visits a comparable number of candidate points indicates that the additional pro active seeking of the largest hypersphere devoid of high cost points is also a significant factor in the success of LEH GA.

Conclusions and further work
We have introduced a new metaheuristic for box-constrained robust optimisation problems with implementation uncertainty. We do not assume any knowledge on the structure of the original objective function, making the approach applicable to black-box and simulationoptimisation problems. We do assume that the solution is affected by uncertainty, and the aim is to find a solution that optimises the worst possible performance in this setting. This is the min max problem. Previously, few generic search methods have been developed for this setting.
We introduce a new approach for a global search based on distinguishing undesirable high cost -high objective value -points (hcps), identifying the largest hypersphere in the decision variable space that is completely devoid of hcps, and exploring the decision variable space by stepping between the centres of these largest empty hyperspheres.
We demonstrated the effectiveness of the approach using a series of test problems, considering instances of varying dimension, and comparing our LEH approach against one metaheuristic that employs an outer particle swarm optimisation and one from the literature that uses multiple re-starts of the local descent directions approach. For low and moderate dimensional instances the approach shows competitive performance; for high-dimensional problems the LEH approach significantly outperforms the comparator heuristics for all problems.
There are several ways in which this work can be developed. Further consideration can be given to the inner maximisation search approach in order to better understand the trade-off between expending function evaluations on the local Γ-radius uncertainty neighbourhood search versus globally exploring the search space, in the context of our LEH approach.
The repeated calculation of large numbers of Euclidean distances each time a new LEH needs to be identified within the LEH GA heuristic is computationally expensive. Rather than only calculating a single next candidate point each time the GA is performed, identifying multiple points could speed up computation or alternatively enable the use of larger population-generation sizes to improve the estimation of the largest empty hypersphere.
Results of the mid-dimension experiments on the Rosenbrock and Sphere test problems suggest that an exploitation based approach works well in these instances, indicating a direction for extending our exploration focussed LEH approach.
It is clear that within the LEH algorithm the early stopping of the inner searches when it is established that the current robust global value cannot be improved upon has significant advantages. It is worth considering whether alternative search approaches could take advantage of this feature.

A. Test functions
Functions used to assess the effectiveness of the Largest Empty Hypersphere robust metaheuristics taken from [Kru12,JY13].

B. Radii due to alternative LEH algorithms
Whilst the Voronoi based search exemplified by Figures 4e and 4f in Section 4.5 is a good indicator of the nature of the searches due to all three alternative LEH approaches, random, GA and Voronoi, the radii of the LEH identified for each candidate will vary across these approaches. Here Figure 9 gives some indication of how the radii of the hyperspheres generated by each these three LEH heuristics progress as the exploration proceeds. The three curves represent separate runs of the LEH algorithm when applied to (poly2D),and should be considered indicative.
As would be expected the general nature of the size of the radius of the LEH steadily decreases with increasing numbers of candidate points evaluated. However superimposed on this overall decrease are  Figure 9: Alternative LEH approaches applied to the same problem: variation in empty hypersphere radius with numbers of candidates evaluated for robustness.
the indicative patterns due to the alternative heuristics. For the random algorithm the size of the LEH is quite variable, whilst for Voronoi the curve is smooth. The GA algorithm sits somewhere between the two.

C. Additional results
The standard deviations due to each set of 50 sample runs, the average number of candidate points visited and average number of function evaluations undertaken, are shown in Tables 3, 4, 5 and 6 here respectively. Labelling of comparator heuristics in the tables is as follows: • PSO: particle swarm optimisation.
• LEH GA: LEH using a genetic algorithm.