Global solution of constrained min-max problems with inflationary differential evolution

This paper proposes a method for the solution of constrained min-max problems. The method is tested on a benchmark of representative problems presenting different structures for the objective function and the constraints. The particular min-max problem addressed in this paper finds application in optimisation under uncertainty when the constraints need to be satisfied for all possible realisations of the uncertain quantities. Hence, the algorithm proposed in this paper search for solutions that minimise the worst possible outcome for the objective function due to the uncertainty while satisfying the constraint functions in all possible scenarios. A constraint relaxation and a scalarisation procedure are also introduced to trade-off between objective optimality and constraint satisfaction when no feasible solutions can be found.


Introduction
A min-max optimisation problem aims at minimising, with respect to a vector d defined in some space D, the maximum value of a given cost function with respect to a vector u, different from d. In this paper we are interested in a class of constrained min-max problems where the constraints have to be always satisfied for all values of u in a given set U .
Min-max problems appear in a wide range of applications. They can be found in the optimal sequence of moves for a computer player in famous board games: the Egyptian ancient game Seega (Abdelbar et al. 2003), the game of checkers (Hughes 2005) and chess (Shannon 1950). They are used to formulate the design for robustness of systems affected by uncertainty: design of electric circuits (Agnew 1981), design of on-line controllers (Sebald and Schlenzig 1994), of aerodynamic shapes (Ong et al. 2006), of the location of sensors and pursuit-evasion games between missile and target (Kim and Tahk 2001).
Due to their wide applicability, they can appear in different forms. They can be seen indeed as a special case of bi-level optimisation problems where two agents quantify their utility with opposite fitness functions. From a game theoretic point of view minmax problems are two-players zero-sum games where the optimal strategy brings the two agents to a Nash equilibrium. From an engineering point of view min-max problems can be seen as deterministic approaches to make decisions under uncertainty. In this case one can understand the engineering problem as a zero-sum game where the two antagonist players are respectively the designer, handling the decision variables, and Nature, handling the uncertainty variables.
Mathematical Programming (Agnew 1981;Chaney 1982) has been widely used to solve min-max problems. An overview of classical discrete approaches can be found in Aissi et al. (2008). The use of Mathematical Programming, however, requires, strong assumptions on the nature of the problem and tends to be problem specific.
On the other hand heuristic approaches (Cramer et al. 2009;Lung and Dumitrescu 2011;Cavalier et al. 2007) and in particular evolutionary based algorithms (Cramer et al. 2009;Lung and Dumitrescu 2011;Laskari et al. 2002;Zhou and Zhang 2010) appear to be more flexible. For example one can find a form of Genetic Algorithm (GA) in Cramer et al. (2009), a special version of the Differential evolution (DE) called Crowding Differential Evolution (CrDE) combined with the Nash ascendancy relation in Lung and Dumitrescu (2011) and a form of Particle Swarm Optimisation (PSO) in Laskari et al. (2002).
As stated in Zhou and Zhang (2010) we can divide the class of evolutionary approaches for min-max problems in three branches. The first one includes problems with a discrete set of possible scenarios U . In Laskari et al. (2002) U is by definition discrete and small, in Cavalier et al. (2007), initially continuous, it is instead discretised by a uniform grid while in Cramer et al. (2009) it is reduced to be discrete by a random selection. The second branch considers a continuous space U and directly solves a sequence of nested problems (Agnew 1981;Sebald and Schlenzig 1994). This approach could be computationally intractable if no strategy is implemented to alleviate the cost. Finally, a third category uses a form of co-evolution optimisation (Kim and Tahk 2001;Cramer et al. 2009;Tahk and Sun 2000;Shi and Krohling 2002) where two separate populations evolve simultaneously with a predator-prey interaction in D and U separately each one being competitor and environment to the other. As demonstrated in Jensen (2001) however most of the co-evolutionary algorithms for min-max have to satisfy a symmetric assumption known as Issac's condition (Tamer Basar 1982). A class of problems that can be properly modelled under this assumption exists (Tahk and Sun 2000;Barbosa 1999) and for them the co-evolution approach is proved to converge to the global solution (Power 1998). However this is not the general case and new methods to overcome this limitation are presented in Branke and Rosenbusch (2008).
In the interest of completeness, we include also the Best Replay approach cited by Marzat et al. (2013) where two optimisation problems are alternated until convergence, one minimising over D and the other maximising over U . It has been proven, however, that this approach often does not converge or it cycles through wrong design candidates, problem known as the Red Queen Effect.
An other research line, applicable to any previously presented method, deals with the reduction of the computational complexity. It makes use of surrogates which approximate the fitness function (Lung and Dumitrescu 2011;Zhou and Zhang 2010;Marzat et al. 2013). In particular, Zhou and Zhang (2010), Marzat et al. (2013) consider response surfaces within an Efficient Global Optimisation (EGO) where they alternate the minimisation in D of the surrogate of the maxima in U and the maximisation in U of the real function for the obtained optimal solution in D.
When it comes to constraint handling, in the existing literature, only few papers could be found that have explored how to deal with constraints in min-max optimisation. Most of them need to start from some strong assumptions on the nature of constraints and cost functions and have been developed for constrained bi-level problems and not specifically for the treatment of min-max problems (Sinha et al. 2020(Sinha et al. , 2018. In Kim and Tahk (2001) the duality between primal constrained minmax problem and dual unconstrained min-max problem with the Lagrange multipliers is demonstrated for both separable and non-separable constraints under the Issac's condition. The constrained problem is then translated in the unconstrained one and a co-evolutionary algorithm is finally used to converge to the saddle point. Other works instead, still dealing with both min-max problem and constraint handling, go in the opposite direction. They indeed translate a single-level constrained problem in an unconstrained min-max problem and use the co-evolution approach to solve it (Tahk and Sun 2000;Shi and Krohling 2002).
This paper proposes a novel algorithm for the solution of a class of constrained min-max problems. The proposed algorithm is inspired by the procedure suggested by Shimizu and Aiyoshi (1980) and further elaborated in Zhou and Zhang (2010), Marzat et al. (2013). The idea is to alternate a minimisation and a restoration (or maximisation) process. As the alternation of these two processes progresses, we incrementally build a discrete representation of the space of the maxima in U . This representation is then used in the minimisation process to converge to an approximation of the desired min-max solution. This mechanism takes memory of the previous solutions avoiding the optimiser to follow the same path again and again. For this reason the proposed algorithm does not suffer from the Red queen Effect. The main novelty is the extension of this logic to the handling of constraints that need to be satisfied for all maxima in a given set. The paper proposes also a constraint relaxation and a scalarisation techniques that allow convergence to a solution in the case the given set U is too restrictive. A multi-population inflationary differential evolution algorithm is proposed to address both the constrained minimisation and the maximisation.
The paper is structured as follows. First we introduce the constrained min-max problems of interest. We then present the general algorithm together with a complexity analysis. The paper then introduces a benchmark of synthetic functions with increasing complexity. Finally the proposed method is tested on a real application case of space systems engineering under epistemic uncertainty.

Problem statement
This paper is concerned with the following class of single-objective constrained minmax problems: where f is the objective function and c i is the i-th constraint function. Both f and c i are defined on the space D × U and depend on a vector of design (or decision) variables d ∈ D ⊂ R n and a vector of uncertain (or environmental) variables u ∈ U ⊂ R m . The solution d opt of Eq.
(1) has two properties: it satisfies all the constraint functions c i over the whole uncertain domain U and minimises the worst realisation of the objective function f over U . Furthermore, we assume that both f and c i are locally C 2 . This assumption can be relaxed if the local search can handle non-differentiable or discontinuous problems.

A memetic single objective constrained min-max approach
The proposed algorithm is a bi-level optimisation procedure based on the alternation of a minimisation and a restoration step. The minimisation step searches for a global solution to the constrained min-max problem: (2) While the restoration step searches for a solution to the following two global maximisation problems, given the solutiond coming from problem (2): The two archivesĀ u f andĀ uc contain respectively the solutions of problem (3) and (4). By iteratively alternating the minimisation and restoration steps, one fills the two archives with the maxima of problems (3) and (4). Thus we can say that problem (2) searches for an optimal d over a discrete representation of the space of the maxima of problems (3) and (4). This mechanism was first proposed for unconstrained problems by Shimizu and Aiyoshi (1980) and Marzat et al. (2013). A solution approach using a memetic algorithm was then introduced in Vasile (2014) and the constrained version was first formulated in . In this paper we provide five main contributions: (i) a detailed algorithmic presentation of a memetic solution approach, (ii) a constrained relaxation strategy, (iii) a scalarisation strategy, (iv) an algorithmic complexity analysis and (v) an extensive benchmark on which our proposed solver of constrained min-max problems can be tested.
The solution approach is summarised in Algorithms 1-8 and explained in the following subsections with the help of the flow diagrams in Figs. 1 and 2 and the examples in Figs. 3, 5 and 6.

Initialisation
The algorithm is initialised either by selecting a random valued or with a given first guess and then by searching for a first solution to problems (3) and (4). These first solutions are then saved into the two archives A u f and A uc . Problems (3) and (4) are global optimisation problems in general. In the following we will propose the use of a memetic global optimiser which combines Differential Evolution with a multi-restart mechanism and local search. Given the global nature of the search for a solution, one needs to set the amount of computational resources allocated to the solution of each of the minimisation and maximisation problems. We quantify this resource in terms of objective function evaluations. More specifically, with reference to problem (3) the maximum number of function evaluations n feval,max is the number of calls to the objective function f (d, u) while for problem (4) we count the function calls to max i c i (d, u). Once the maximum number of function evaluations per sub- Fig. 1 Flow diagram of the constrained min-max alternative approaches described in Sects. 3.1, 3.2 and 3.4 -3.7 and summarised in Algorithms 1-8. In particular, this diagram describes the sequence of optimisation problems applied for the standard approach and its alternative strategies: the constraint relaxation and the scalarisation. The first optimisation problem is defined in block A and it refers to the general constrained min-max in Eq. (1). If a solution exists for this problem, no further analysis is required (link 1). If differently, two alternatives are given. The first (link 2) brings to the relaxation strategy (block B). Here an unconstrained min-max problem is solved to find . The relaxed constrained min-max problem is then solved (Sect. 3.4). Link 3 instead brings to the scalarisation strategy (Sect. 3.5). Block C is first solved to find the reference points in the Pareto front (ideal and nadir). Finally, the scalarisation procedure is activated (block D) problem are defined one needs to define the maximum number of iterations n loop,max of minimisation/restoration.
The steps for the initialisation are summarised in Algorithm 1. Line 1 presents how to define the problem in order to be solved by Algorithm 2: the required information on objective function f (d, u), set of constraint functions c(d, u), dimensions n D and n U of the design and uncertain vectors d ∈ D ⊂ R n D and u ∈ U ⊂ R n U and the bounds of each variable in d and u. Line 2 defines the constraints on the computational resources for the single sub-problems in Eqs. (2) -(4). They are the maximum number of function evaluations n outer feval,max for problem (2), n inner feval,max for problem (3) and n inner,c feval,max for problem (4). Line 3 defines the threshold σ stop below which Algorithm 2 is considered to have converged. Line 4 defines the limit on the computational cost for the whole procedure (1) with its maximum number of function evaluations n fval,max and maximum number of iterations n loop,max between minimisation and restoration.  Fig. 1. The relaxation and scalarisation alternatives (blocks B,C and D) follow however the same logic. Algorithms 2, 1, 4, 7 and 8 and Sects. 3.1, 3.2, 3.6 and 3.7 explain in detail all the represented blocks. Each of them lists the main operations that are performed: the optimisation problems and the archives updating. The algorithm is initialised. Then there is the iteration between blocks 'Restoration' and 'Minimisation'. Finally the 'Cross Check' is performed and the solution is chosen Line 5 defines the optimisation algorithms and its parameter settings that are used to solve the sub-problems in Eq. (2) -(4). On line 6 the initial design vector is defined by the user or initialised randomly with a hyper-cube sampling procedure. In the case some important uncertain scenarios in U or design configurations in D are already known, the initial archives can be seeded with these scenarios, otherwise they are initialised as empty sets (line 7). Similarly, it is done for the archives relative to functions f and c (line 8). Lines 9 and 10 initialise the counters of the number of function evaluation n feval and the iteration n loop between minimisation and restoration processes. The initial accepted constraint violation (line 11) is initialised to zero.

Minimisation-restoration loop
The main algorithm is summarised in Algorithm 2. Using the first design guessd from the initialisation, Eqs. (3) and (4) are solved in parallel (lines 4 and 5) within the first restoration, or inner, loop (lines from 2 -8). Withd fixed, the former equation evaluates the feasible worst case condition of f while the latter determines the worst constraint violation. In the following we will adopt a multi-population version of Infla- The main loop is then started (line 9 of Algorithm 2 and Eq. (2)) with an alternation of the minimisation, or outer, step (lines 10, 11 and 12) and restoration, or inner, step (lines from 13 to 17). The latter has been already described. A further check is, however, performed in lines 1 to 8 of Algorithm 4 before the updating procedures of the archives. If the condition in line 2 holds, indeed, the inner loop in Eq. (3) has failed. Then the solution generated in this loop is discarded and the one of the previous outer loop is maintained. This condition follows the same criterion used in Algorithm 3 by giving priority to the constraint satisfaction. Similarly the condition in line 6 for Eq. (4).
The outer step addresses the solution of the constrained minimisation in Eq. (4). Also for this minimisation loop we propose the use of a multi-population version of Inflationary Differential Evolution . For each new design vector generated by the optimiser the cost function and constraints are evaluated for all the u vectors inĀ u f andĀ uc and the worst cost function and constraint violation values are retained (line 11). Note that in some cases it is desirable to run a local search every time a vector in the two archivesĀ u f andĀ uc is evaluated (see for example Vasile (2014)). This added local search significantly increases the computational cost of each single evaluation of the outer loop but in some cases improves convergence to the point that the overall cost of the algorithm is reduced. For this reason we inserted this option in the algorithm although it is not tested in this paper. It was, however, tested, for the case of unconstrained min-max problems, in a previous work by the authors Vasile (2014). At the end of the minimisation process, the archive A d of the design configurations is updated with the new design solutiond.
If the condition in line 18 of Algorithm 2 holds, the constraint relaxation strategy described in Algorithm 5 is activated. An alternative option, the trade-off strategy in Algorithm 6, can instead be activated if condition in line 21 is true. These alternatives can be visualised in Fig. 1. If the former condition holds, problem A is stopped, the relaxation step in the first block in B is performed finding and finally the constrained min-max problem is updated as in the second block in B. If instead the second condition hold, block A is stopped and then blocks C and D are solved.
The minimisation and restoration steps are alternated until condition in line 9 is satisfied. In particular, the iteration holds until the number of calls to the function has not reached the maximum allowed (n feval < n feval,max ), the number of iterations is below the upper bound (n loop < n loop,max ) and the solutions saved in archives have not converged (max A (σ ) > σ stop ) where max A (σ ) is the maximum standard deviation between all the archives in the last 3 iterations. The violation of at least one of these conditions corresponds to the termination criterion.
Then, the cross-check between all the design vectors archived in A d is performed (line 26). The cross-check procedure is explained in Sect. 3.6 and summarised in Algorithm 7. Finally, the solution is chosen following Algorithm 8 (line 28).
An example of the application of Algorithm 2 to the two functions MWP10 and GFc1 (see Sect. 4.1), without constraint relaxation and trade-off can be visualised in Fig. 3. In sub-figures c and d an initial design guessd and each new designd proposed by the minimisation step are represented as vertical lines. The corresponding worst case scenario for f and c are plotted, with the same colour, as dots and stars respectively. Sub-figure e shows sections of f and c over U for different design configurations. Sub-figure f shows how the space of the maxima of f and c as it appears to the minimisation process. In particular, sub-figure (c) shows that the algorithm is able to find at the first iteration, for this test case, the subspaceD ∈ D in which the design solutions are feasible for any uncertain scenario. The first design guessd (blue line) finds the worst scenario for u = 10. In b indeed, the maximum constraint violation is brought to zero at the second iteration and all the other iterations are then used to minimise the worst case of f working within the design domainD.

Memetic strategy
When problems in Eqs.
(2) to (4) within the min-max Algorithm 2 are global optimisation problems we propose the use of the memetic optimisation algorithm Multi-Population Adaptive Inflationary Differential Evolution Algorithm (MP-AIDEA) (Di . It combines the Darwinian evolution of a number of populations of candidates through a Differential Evolution (DE) strategy with the Lamarckian evolution of the best agent in each converged population through a local search algorithm. In particular, the local refinement is performed if the converged solution is not in the basin of attraction of previous local minima. A number of local restart are then performed before globally restarting.
Algorithm 2 Constrained min-max 1: Initialisation: Algorithm 1 2: ifĀ u f = ∅ ∧Ā uc = ∅ then 3: Initialisation loop: 4: Run if multiple outputs, choose best u a f as in Algorithm 3 7: Update archives as in Algorithm 4 8: end if 9: while n feval < n feval,max ∧ n loop < n loop,max ∧ max A (σ ) > σ stop do 10: Minimisation loop: Run if multiple outputs, choose best u ac : Algorithm 3 17: Update archives: Algorithm 4 18: if relaxation flag ∧ not convergence ∧ satisfy limits on n feval and n loop then 19: Stop An example of a run of MP-AIDEA with n pop = 2 for problem MWP-1&GFc-1 is in Fig. 4 which shows, for the evolving populations, the alternations of DE steps together with the local refinements. For this test case the MATLAB function solver fmincon has been used. Given the constraint on the maximum number of function evaluations, the algorithm is able to perform only one local restart in each population. Colours red, blue and green refer to the first population while colours orange, yellow and brown refer to the second population. As shown in the figure, first the two populations are initialised randomly and evolved independently until convergence is Convergence of MP-AIDEA with two populations in the outer loop of Algorithm 2 for test problem MWP1&GFc1. Coloured areas represent the convergence of the differential evolution steps (different set of colours for different populations) while dots represents optimal solutions of the local search achieved (being the parameter ρ one of the termination criteria). In particular the red (orange) area represent the evolution of the best agent in the population while the blue (yellow) represent the mean value and the green (brown) the worst agent. From the two converged solutions a local search is performed with f mincon until the green and brown points are obtained. From these two local minima then the two populations are locally reinitialised. The overall process is then repeated till convergence.

Constraint relaxation strategy
The min-max problem proposed in this paper imposes quite stringent conditions on the satisfaction of the constraints as constraints need to be satisfied for all possible values of u ∈ U . It is, therefore, possible that no solution d opt is feasible in all U . Since we are interested in the worst case solution for both constraints and objective function, when no feasible d is possible we introduce an automatic relaxation of the constraints.
In order to find a set X ⊂ D that is feasible for all u ∈ U , we first solve the following minmax problem: Problem (5) minimises the maximum violation of the constraints and returns a solution vectors d min,c and u min,c . Vectors d min,c satisfies the constraint: The relaxation strategy is explained in Algorithm 5. If condition on line 18 of Algorithm 2 is satisfied then Algorithm 5 is triggered. Once in Algorithm 5, until condition in line 2 holds, problem (5) is solved (see lines 3 -8 of Algorithm 5), with the iteration between the following minimisation and restoration step max This is an unconstrained min-max formulation where the optimised function is the vector of constraints c. The solution at convergence is the minimum over D of the worst constraint violations in U . In line 10 that value is associated to the relaxation parameter . In this way the algorithm is re-conducted to Eq. (5) by means of the relaxation of the constraint violation and Eq. (1) is translated in:

Scalarisation strategy
This subsection presents a different point of view that is based on the parameter-based scalarisation approach (Kasimbeyli et al. 2019) for the solution of Eq. (1). It is here supposed that the decision maker, who is the final user of this methodology, is able to assign preference weights to the different Quantity of Interests (QoIs) involved in the problem under analysis. The scalarisation approach then has been preferred to the direct multi-objective optimisation because, for a given set of weights, it reduces the Run if multiple outputs, choose best u ac : Algorithm 3 8: Update archives A uc and A c : Algorithm 4 9: end while 10: set = max u∈U max i c i (d, u) 11: Restart Algorithm 2 computational complexity. Problem (11) indeed can be seen as an Epsilon-Constraint Scalarisation (ECS) (Filippi and Vasile 2020;Haimes et al. 1971) formulation of the bi-objective min-max problem: Thus one could be interested in a trade-off between f and and accept larger relaxations of the constraints in favour of a better objective. Note that in the context of uncertainty quantification this implies accepting a higher probability of violating the constraints in favour of a better cost function. We now introduce the assumption that a preference vector ω = [ω f , ω c ] T can be defined a priori and explain the scalarisation approach summarised in Algorithm 6.
In the first part (line 2 -10) the reference points c ideal , f nadir , f ideal and c nadir are calculated. c ideal is the minimum (best) over D of the worst case constraint violations in U : and it is equal to the relaxed constraint . For the corresponding design vector d c-ideal the worst scenario for the objective function is f nadir (lines 3 to 8): In line 9 the unconstrained min-max problem is solved to define the best design configuration d f -ideal that minimises the worst scenarios of the objective function f regardless the constraint violation. d f -ideal is then used in line 10 to calculate the corresponding worst condition for the constraint violation: An example of the reference points for a generic Pareto front applied to the min-max problem is in Fig. 5 where z ideal = [ f ideal , c ideal ] and z nadir = [ f nadir , c nadir ]. In the second part (lines 13 -27) it is instead described the scalarisation procedure. Nadir and ideal points are here used to normalise f and c: Three different methods were considered in this study: ECS, Weighted-Sum Scalarisation (WSS) (Gass and Saaty 1955) and Chebychev/Pascoletti-Serafini Scalarisation (CPSS) where the last one is a smooth combination of Weighted Chebyshev Scalarisation (WCS) (Bowman 1976) and Pascoletti-Serafini Scalarissation (PSS) (Bo et al. 2007). A comparison of the different approaches can be found in Filippi and Vasile (2020). The difference in Algorithm 6 is only within the minimisation loop (lines 13 -22).
The WSS solves the minimisation over the design space D of the weighted sum of the worst case scenarios for f and c in the respective archives (line 18): The CPSS solves two different problems. During the DE step of Inflationary Differential Evolution, the algorithm addresses the following problem: where we have introduced the weight vector ω = [ω f , ω c ] T and the objective vector During the local search step of Inflationary Differential Evolution, instead, the CPSS addresses the following constrained minimisation problem: where z f and z c are the values off andc respectively for the best solution obtained from the DE step in Eq. (20) and from which the local search is initialised. For all three methods, the restoration problem solves the two problems: They are theoretic points that collapse the extreme behaviour of the different solutions in the Pareto front. z ideal is the combination of the best solutions for the different objectives. z nadir represents instead the worst possible combination of points. z utopian is finally defined by means of an from z ideal For the given preference vector, the optimal worst case condition for f and the optimal relaxation of the constraint are given, at convergence, by Eqs. (22) and (23). The three possible scalarisation procedures ECS, WSS and CPSS can be used to reconstruct the Pareto front for the minimisation trade-off between the worst case scenarios max f and max c. Fig. 6 shows a comparison of the three methods, where the combination of the objective function G F f 1 in Table 1 and the constraint function G Fc1 in Table 2 has been chosen as representative test case. We used MP-AIDEA with the following settings: maximum number of function evaluation n feval,max = 3000, number of populations n pop = 2, number of agents in the population N pop = 5, dimension of the bubble for the global restart δ global = 0.1 and DE threshold for each population ρ = 0.1. The whole constrained min-max algorithm has then been run until convergence. For ECS the algorithm has been run for 60 different thresholds ranging between 1 and 6.2 (calculated minimum and maximum constraint violation.) For WSS and CPSS, instead, the trigonometric weights have been used letting θ varying from 0 to π 2 and using a discretisation of 60 interval as well. Fig. 6 shows that WSS is not capable of finding optimal Pareto points in the non convex part of the Pareto front (Das and Dennis 1997) and also, using equally spaced weights, it finds non-equally spaced points in the front. The ECS strategy gives better results. The best performance is, however, obtained with the CPSS. Indeed it allows the user to express a preference through the selection of the weights and, as shown in Fig. 6, with the same number of simulations (60 different descent directions for the Fig. 6 Pareto front corresponding to the trade-off between the two conflicting goals max f and max c. The results refer to test case G F f 1&G Fc1. In particular, the -constraint (EC) approach has been applied with different thresholds while the Chebychev/Pascoleti Serafini (CPSS) and the Weighted-Sum (WSS) scalarisations have been applied with different preference vectors CPSS and epsilon values for the ECS) it is able to find more Pareto optimal solutions than ECS. For these reasons CPSS was chosen and implemented in Algorithm 2

Cross-check
All the optimisation problems in the minmax algorithm require the identification of a global maximum or a global minimum. Since it is proposed to use a memetic algorithm it is possible that some of the maxima or minima in the archive are only locally optimal. Note that the use of a deterministic global optimiser would remove this problem but would introduce a tractability problem due to the potential NP-hard nature of some optimisation problems.
In order to mitigate the occurrence of local minima/maxima in the archives we introduce a cross-check of the solutions following the procedure explained in Algorithm 7. It is performed for each design vectord that can be proposed by the optimiser during the minimisation step and at the end of the whole algorithm (respectively in line 9 and 19 of Algorithm 2). Referring to Algorithm 7, lines 1 -7 regard the objective function f while lines from 8 -14 regard the constraint function c. In both cases, for a givend objectives and constraints are evaluated for all the u vectors in the archivesĀ u f and A uc . We also introduced an option (through local flag) to run a local search from each new pair [ N d, u]. This option slows down Algorithm 2 but improves the quality of the solution if the functions present nested minima/maxima. Finally, line 15, retains the worst values of f and c for the archivesĀ u f andĀ uc for eachd.

Selection of output solution
After the termination criterion in Algorithm 2 is applied and the cross-check over the archives is performed (line 21), the solution for the min-max problem is selected following Algorithm 8. In particular, if a feasible subsetÂ d of the archive A d of the design vectors exists (line 1) the selected solution vector is the one, withinÂ d , min with the cross-check as in Algorithm 7, f andc defined in Eqs. (17) and (18) (17) and (18)  minimising the worst value of f (line 2). If, on the other hand,Â d is an empty set, the design vector that minimises the constraint violation is selected (line 4).

Computational complexity
The computational cost of Algorithm 2 is measured in terms of number of function calls. With reference to the minimisation step, the counter n outer feval takes into account, for both the constrained and the unconstrained min-max problems, the calls to max u∈Ā u f f (d, u a f ) in Eq. (2). The same criterion, then, holds for the constraint Algorithm 7 Cross-Check 1: for all elements u a f ∈Ā u f do 2: if local flag then 3: Compute local maximum f (d, u * a ) s.t. max i∈I c c i (d, u * a ) ≤ from u a f 4: else 5: Compute f (d, u a ) s.t. max i∈I c c i (d, u a f ) ≤ 6: end if 7: end for 8: for all elements u ac ∈Ā uc do 9: if local flag then 10: Compute local maximum max i∈I c c i (d, u * a ) from u ac 11: else 12: Compute max i∈I c c i (d, u ac ) 13: end if 14: end for 15: For eachd Save worst vectors u a f and u ac in the archivesĀ u f andĀ uc .

Algorithm 8 Select Solution -Output
relaxation step in Eq. (9) and for the two trade-off steps in Eqs. (15) and (19). It has to be noted that, as the algorithm proceeds in the search of the global optimum solution, the archivesĀ u f andĀ uc of the uncertainty vectors increase progressively in dimension. Each minimisation step explores a maximum number of possible design configurations which is limited by the input parameter n outer feval,max . However, due to the growth of the archives of the solutions coming from the restoration step, each evaluation of the minimisation loop becomes increasingly more expensive.
With reference to the maximisation step, instead, the cost of the two separate problems in Eqs. (3) and (4) has to be considered. For the former n inner, f feval counts the number of calls to the objective function f (d, u a f ) and it is limited by the input n inner, f feval,max . This holds true also for the two steps in the trade-off strategy in Eqs. (16) and (22) and for the relaxation step in Eq. (10) where the function c is considered instead of f . For the latter, n inner,c feval counts the number of calls to max i c i (d, u) in Eq. (4) where the input n inner,c feval,max is the upper limit. Finally, the parameters n feval,max and n loop,max give an upper limit on the whole cost of Algorithm 2. Note that n feval,max and n loop,max represent an upper limit because, as it is shown in line 9 of Algorithm 2, we use an additional termination criterion that looks at the convergence of the solutions in the archives.
The computational complexity of the different parts of the overall algorithm is as follows: 1. Local Search: the local search uses the Matlab fmincon function. All the alternative algorithms (interior-point, trust-region-reflective, sqp, sqp-legacy, active-set) can be selected. We use here interior-point that works well with both large sparse and small dense problems. The complexity is O(n 3 D ) or O(n 3 U ) depending on which step between minimisation and restoration is considered, where n D is the design and n U the uncertain vector's dimension.

Testing procedure
Algorithm 2 has been tested on the benchmark described and explained in this section. Each test case is a combination of an objective function f and a constraint function c. Depending on the mathematical features of each problem, a local optimiser or a global optimiser have been used for the three problems in Eqs.
(2) to (4). The criteria used to choose the right optimiser is explained in Sect. 4.3. Given the stochastic nature of MP-AIDEA, each optimisation for each problem has been repeated 100 times. Results are then reported in Section 5. For the evaluation of the algorithm's performance, the Success Rate S R is used instead of best value, mean, and variance. The SR was suggested in Vasile et al. (2011) for a generic problem min f and a generic algorithm. It is here generalised to consider also the handling of constraints. The definition of SR is in Sect. 4.2.

Benchmark
The equations of f and c are listed in Tables 1 and 2 respectively. The constraint functions c are more extensively presented in the following Sects. 4.1.1 -4.1.10 and visualised in Figs. 7,8,9,10,11,12 and 13 for the case n D = n U = 1. Table 3 lists lower and upper bounds, dimensions and optimal solutions for the unconstrained problems in Table 1. The same solutions holds also for the constraint min-max problems for which the constraint does not change the global optimum. Table 4  Both f and c are designed to include different structures that can be encounter in practice (Jamil and Yang 2013). In particular, they exhibits the following features: 1. Modality: number of local optima that try to trap the algorithm in the wrong peak. 2. Basin or plateau: a relatively steep decline that surrounds a large area. There is no information to drive the algorithm. 3. Valley: similar to the basin but it is narrow area. 4. Non separability: property related to the coupling between parameters. Nonseparable functions are in general more difficult to optimise. 5. Dimensionality: property related to the number of parameters or dimension of the problem. The search space increases with the dimension, increasing then also its difficulty. 6. Non differentiability: cuspids, corners, tangents and discontinuities are features that make functions non differentiables in some points. Some of the constraint functions present cuspids, corners and discontinuity. In particular discontinuity is an abrupt change in the function values. Discontinuities are classified in jump, infinite, removable, endpoint, or mixed. Some of the constraint functions c present jump discontinuities.

GFf-1
GFf-1 is a modifications of the Rastrigin function where half of the variables are design parameters and the others are uncertain variables.   (25) It is continuous, differentiable, scalable, without valleys and basins, and highly multi-modal with hundreds of local peaks.

GFf-2
GFf-2 is a variation of the saddle-point function MWP-8: where both components d k,R and u k,R are obtained rotating d k and u k respectively by the angle GFf-2 is continuous, differentiable, non-separable, scalable, without valleys and basins, and uni-modal.

GFc-1
GFc-1 is a hyper-plane and it is a linear function in both d and u: where K = − i d i − i u u,i − 0.05 with u u,i the upper bound for the i-th uncertain variable. GFc-1 is continuous, differentiable, separable, scalable, without valleys and basins, and uni-modal.

GFc-2
GFc-2 is a modification of GFc-1. It is a continuous piece-wise linear function where the feasible region is a plateau. It is the intersection of two hyper-planes, the second being at the border between of feasible and infeasible regions.

GFc-3
In GFc-3 there are a jump discontinuities, valleys and plateaus. The feasible area is a narrow multidimensional circle. The function is not differentiable, scalable and uni-modal:

GFc-4
GFc-4 is a modification of the Rastrigin function where a jump discontinuity is introduced: It is highly multi-modal, discontinuous, not differentiable with valleys and plateaus, separable and scalable.

GFc-5
GFc-5 is a modifications of Eq. (31). Here a rotation of the vectors d and u is also introduced.
The rotated components d k,R and u k,R are given by the angle θ k = d i + 2u i .

GFc-6
GFc-6 is a multi-dimensional peak function with high coupling between D and U. It is unfeasible in most of the domain while it is satisfied only in few narrow non linear valleys varying with d.
It presents very narrow non-linear valleys. It is continuous, locally non differentiable, without plateaus, scalable, non separable and multi-modal.

GFc-7
GFc-7 is a multi-dimensional peak functions with high coupling between D and U and narrow unfeasible areas varying with d. (d u (sin(d u,i where d u,i and u u,i are the upper bounds for d i and u i respectively. GFc-7 has very narrow non-linear peaks with large plateaus. It is continuous, locally non differentiable (for the cuspids), without valleys, scalable, separable and multi-modal in D.

GFc-8
GFc-8 is a rotated versions of MWP-8: where d i,R and u i,R are obtained from d i and u i respectively with the rotation angle in Eq. (27). It is non-separable, scalable, continuous, differentiable and uni-modal and without valleys and plateaus.

Success rate
The Success Rate (S R) is adopted here for the performance assessment of Algorithm 2. Its definition is given in Algorithm 9 for a generic algorithm A applied to a generic constrained problem C P on the D × U space. It is defined as the ratio j s n between the index of performance j s and the number of independent experiments n.
First, all the parameters required by Algorithm 2 are fixed (refers to the initialisation in Algorithm 1). The following parameters are then defined:  Tables 3  and 4. In particular, δ k c depends on the uncertain vector u opt,c that is the worst for the constraint function c while δ k f depends on the vector u opt, f that make worst the objective functions f . δ u is necessary to verify the convergence on the maximisation in the inner loop (restoration in Sect. 3) and then to avoid counting as success solution an f k opt close to f ref that is coming from a lucky combination of a wrong maximisation and a wrong minimisation in the outer loop (optimisation in Sect. 3).
Algorithm 9 Success Rate 1: Define the parameters for algorithm A to solve the constrained problem C P; 2: Define the number of repetition n; 3: Define tolerances tol f , tol d and tol u ; 4: Initialise the index of performance j s = 0; 5: for k = [1,2,.., n] do 6: Run A on C P with the defined settings; 7: Compute f k opt ← A(C P(d, u)); 8: Compute c k opt ← A(C P(d, u)); 9: Compute d k opt : optimal solution for the design vector 10: Compute u k opt : optimal solution for the uncertain vector 11:  Table 3 Reference solutions for the test cases in Table 1 ID

Algorithm settings
An important feature of the proposed approach is its modularity in the sense that any optimiser can be plugged in and used for the single optimisation problem in Eqs.
(2) -(4). To enhance efficiencies of Algorithm 2, then, the right combination of optimisation solvers should be selected. An optimal choice would require a prior knowledge of the main features of a given problem. For complex multi-modal functions, we suggest the use of the memetic optimisation solver MP-AIDEA because it has shown to be efficient and effective, on average, on a wide range of problems mixing different characteristics. For continuous uni-modal functions we use instead the Matlab fmincon solver with an interior-point scheme. We give here the parameter settings of MP-AIDEA that have been used for all tests. The number of agents for each population N pop and the maximum number of function evaluations were set to be respectively N pop = max[5, n D ], n outer feval,max = 500n D , n inner,f feval,max = 500n U and n inner,c feval,max = 500n U . The dimension of the bubble for the global restart is δ global = 0.1, the number of populations is n pop = 2 and the convergence threshold of DE is ρ = 0.25.

Results
The results are presented and explained in this section. In particular, four sets of test have been performed. In the first, Algorithm 2 has been combined with the optimiser fmincon, while in the other cases MP-AIDEA has been used. First we consider one unimodal problem. The performance of the algorithm is assessed increasing the dimension of the problem. Then we consider the worst-case complexity analysis on the benchmark presented in Sect. 4.1 with a wide variety of difficulties. A complexity analysis on the algorithm convergence is then presented for a selected test case for different problem dimensions. Finally we apply Algorithm 2 to solve a realistic engineering problem: the design for robustness of a communication satellite.

Uni-modal test problem
For the first set of results, the test case used is given by the combination of the objective function GFf-2 and the constraint function GFc-8. They are both continuous, Table 5 Success rates of GFf-2 and GFc-8 for different problem dimensions (rows) and limits on the maximum number of function evaluations (columns). Optimiser: fmincon. 4e5  6e5  8e5  1e6  2e6  3e6  4e6  5e6  6e6  7e6 10 × 10 0.03 0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 differentiable, unimodal and non-separable. With these features a local optimiser is sufficient to solve Eqs.
(2)-(4) at each iteration. The constraint function c admits only one feasible design vector, which is different from the unconstrained optimum of GFf-2. The local optimiser we used in this test is fmincon. The test functions are devised to be scalable with a predictable value of the exact min-max solution. Results are collected in Table 5 for a number of function evaluations up to 7e6. The table shows up to dimension n D = 40 and n U = 40 the algorithm can achieve S R = 1 within the maximum number of function evaluations. For n D = 50 and n U = 50, 7e6 is not enough and the best result is a success rate of 30%.

Multi-modal test problems
For the second set of experiments, Tables 6, 7, 8, 9, 10, 11 and 12 collect the results for all the test cases given by the combination of objective functions f from Table 1 and constraint functions c from Table 2. The last two columns of each table, n iter,min and n iter,max , collect the minimum and maximum number of loops for which the algorithm achieves S R = 1 (rows with the symbol − correspond to problems for which S R = 1 has not been obtained for any of the 100 runs). For almost all the problems Algorithm 2 converges to the correct solution with an S R = 1 within the maximum number of function evaluations. For some problems (namely GFF1-GFC1, MWP10-GFC4, MWP11-GFC5, MWP11-GFC2) few of the runs did not converge to the correct minimum of f but the S R is still reasonably high.

Convergence complexity
Sections 5.1 and 5.2 show the performance of Algorithm 2 with respect to the worstcase computational complexity where an upper bound on the number of function evaluations (n feval,max ) is fixed for each optimisation step of the approach: n outer feval,max , n inner,f feval,max and n inner,c feval,max . It is here interesting to show an other computational complexity analysis that is related to the order of magnitude of the number of function evaluations needed to converge to the optimal solution. In particular, Table 13 summarises the results for the test case G F f -1&G Fc-1. This test case has been selected as representative because it is scalable in both design D and uncertain U spaces Table 6 GFc-1, MP-AIDEA, 2 populations,  Table 7 GFc-2, MP-AIDEA, 2 populations,  Table 8 GFc-3, MP-AIDEA, 2 populations,  Table 9 GFc-4, MP-AIDEA, 2 populations,   Table 11 GFc-6, MP-AIDEA, 2 populations,  and it is also one of the most difficult within the proposed benchmark. Each row in Table 13  f eval have been determined averaging the sum, for the different algorithm's iterations, of the number of function evaluations at convergence. The optimiser MP-AIDEA has been used. In order to assure convergence in each optimisation step in Algorithm 2, the number of populations has been set equal to the problem dimension n pop = n (the problem is highly multi-modal) and the maximum allowed number of function evaluations of each step has been fixed at n outer feval,max = n inner,f feval,max = n inner,c feval,max = 1e4n while for the whole algorithm it is n feval,max = 2e6n. The remaining input parameters for MP-AIDEA have been fixed as in Sect. 4.3.

Space system design
The min-max approach in Algorithm 2 is finally tested on the design for robustness of a Complex Engineered System (CEdS) under uncertainty. The system under analysis is an observation spacecraft and the goal of the mission is the fire detection within a belt centred at the latitude of 50 deg. The spacecraft is modelled as the network shown in Fig. 14 where the nodes correspond to its subsystems and the links to the coupling between them. The mathematical models that have been used for the nodes are a modification of the ones the authors extensively presented in . The differences are described in the following and are in the explicit definition of a node for the orbital dynamic and in the payload subsystem. Design and uncertain variables are listed in Tables 14 and 15 respectively. Within the orbit node, considering a circular Low Earth Orbit (LEO), the altitude h, inclination i, the minimum elevation angle min at which the ground station is    (Fig. 14). These coupling variables are: the period of each orbit P o , the number of orbits N o the satellite perform due to the shift of the longitude of the ascending node, the time of eclipse for each orbit T ecl , the dynamic pressure p dyn , the mean Earth magnetic field strength K m , the gravitational field K g , the maximum distance to the target D max , and the access time T ac (or total time in view) between the target and the satellite, where the target is the ground station at 22 deg of latitude used for down-link and up-link. With exception of T ac the formulas can be found in . Instead, considering that the satellite ground track is determined by the inclination i and by the longitude of the ascending node L node and that the latter increases by 360 deg in 1346 min (the rotation of the Earth relative to the stars), T ac is calculated considering the total number of orbits i that happens during this period T ac = i T ac,i . Following Wertz et al. (1999), that describes the motion of the satellite as seen from a point on the Earth (the ground station), T ac,i is evaluated as: T ac,i = P 180deg arccos cos λ max,i cos λ min,i with λ min,i and λ max,i the minimum and maximum Earth Central Angle for the i-th orbit.
The Payload System is an infrared camera that is used to detect possible fires and its target is the belt at 50 deg of latitude. Within the payload node the model's parameters are h (shared with the orbit node), P o (coupling parameter), the width for square detector d, the quality factor Q, the operating wavelength λ, the maximum incidence angle of the instrument I A max and the maximum ground sampling distance Y max . The model evaluates the following coupling variables: the data volume DV shared with On Board Data handling (OBDH) and the power requirement P shared with the Power System. The model evaluates also the payload mass and the percentage of coverage area PC of each orbit during which the payload target is seen. In particular, PC is calculated following Wertz et al. (1999) as function of λ max , i and the latitude of the target Lat = 50 deg: where cos φ 1 = − sin λ max + cos i sin Lat sin i cos Lat and cos φ 2 = sin λ max + cos i sin Lat sin i cos Lat The remaining couplings between nodes are the compressed data volume DV c that OBDH send to TTC for down-link to the ground station and the power requirements P of all the nodes (orbit excluded) that the Power sub-system has to make available. Finally, the global outputs of the network are the overall mass M of the satellite, sum of the masses of the components, and the percent coverage PC of payload target land. In the optimisation framework, M is considered to be the performance indicator while PC is the constraint to be satisfied. This mission design problem is translated into the following constrained min-max problem: In order to explore the conflict between f and c, the corresponding Pareto front has been reconstructed. We want to apply here the main min-max method presented in Algorithm 2. The algorithm has then been repeated using 30 different values for the threshold ν through an ECS approach.
The results of Eq. (40) are shown in Fig. 15 for which the optimiser MP-AIDEA has been used with the setting specified in Section 4.3. In particular, Fig. 15a presents the Pareto front reconstructed for the different values of ν, while the shape of the front can be understood looking at Fig. 15b,c,d. There are indeed two different geographical targets on the Earth for the defined mission: the ground station (22 deg of latitude) used for up-link and down-link by TTC and the area that has to be monitored by the payload for possible fire detection (50 deg of latitude).
These two targets are quantified in the respectively node's models by T ac and PC, the latter being the constraint function c and the former having an high impact on the final mass M of the overall spacecraft that is the objective function f . The most Sub-figure d finally shows the percent of land coverage by the payload influential design parameters with regard to PC and T ac in the trade-off within the set of optimal Pareto points are the altitude h = d 1 and the inclination i = d 3 . Fig. 15b shows their optimal values while moving in the front, while Fig. 15c,d show finally the corresponding values of T ac and PC. For low values of ν in the constraint function (left side of the Pareto front) the design solution selects the orbit inclination that maximise the amount of time for the link between the spacecraft' antenna and ground station. This configuration reduces the overall mass M at the expense of the capacity of detect fires (PC). As ν increases, the solutions becomes sub-optimal for Fig. 15c while maximises the area in Fig. 15d.

Conclusion
The paper has presented a new algorithm for the solution of a class of constrained min-max problems. The class of min-max problems emerges naturally from the need to make robust decisions under uncertainty in the case in which constraints need to be always satisfied. The method is based on the alternation of a minimisation and a restoration step. This scheme is fairly optimiser agnostic and we demonstrated its applicability even in the case a simple gradient search is used. For the case in which the min-max solution requires the global exploration of a complex solution space we have proposed the use of a memetic approach based on Inflationary Differential Evolution. Our complexity analysis has revealed that the algorithm is overall of polynomial complexity with maximum exponent equal to 2.
The combination of the proposed solution strategy and memetic global optimiser was extensively tested on a new benchmark of objective and constraint functions with a variety of features that can be encountered in real-life applications.
Results show that the algorithm we propose is successful at identifying the constrained min-max solution with a limited number of calls to objective functions and constraints. Such solution minimises the worst case realisation of the objective function in the uncertain space while guaranteeing its feasibility in all the possible scenarios. The benchmark is complemented by a real case of robust optimisation of a space systems.
In the case in which a feasible solution in all the uncertain domain could not be found, we proposed a constraint relaxation procedure to automatically adapt the admissible region. Finally we proposed a trade-off approach between unconstrained min-max solution and constraint satisfaction based on a combination of Chebyshev and Pascoletti-Serafini scalarisation. This approach is promising for cases in which a user defined relaxation of the constraint is possible as it allows one to explore the optimal trade-off curve between optimality and reliability. The use of this approach will be further developed in future work.