PAINT–SiCon: constructing consistent parametric representations of Pareto sets in nonconvex multiobjective optimization

We introduce a novel approximation method for multiobjective optimization problems called PAINT–SiCon. The method can construct consistent parametric representations of Pareto sets, especially for nonconvex problems, by interpolating between nondominated solutions of a given sampling both in the decision and objective space. The proposed method is especially advantageous in computationally expensive cases, since the parametric representation of the Pareto set can be used as an inexpensive surrogate for the original problem during the decision making process.


Introduction
Multiobjective optimization concerns simultaneously optimizing multiple conflicting objectives (see, e.g., [1,2]). Multiobjective optimization has its roots in economics (see [3]). However, many effective applications can be nowadays found, among others, in mechanics, aerospace and automotive engineering, environmental science, medicine and even biology [4][5][6][7][8]. For a typical multiobjective optimization problem there exist multiple Pareto optimal solutions, i.e., solutions where none of the objectives can be improved without impairing at least one of the remaining objectives. However, in most cases, only one solution has to be selected for implementation. The person, or the group, appointed for choosing such a solution is called the decision maker. It is often assumed that the decision maker chooses the solution through an enlightened screening of the range of the available values of the objective functions. Naturally, a better decision can be taken when the decision maker is made aware of an as large as possible set of alternative Pareto optimal solutions.
To this purpose, Pareto front and Pareto set approximation methods become advantageous. Following [9], we mean by a Pareto front approximation a set in the objective space that can be used as a surrogate for the actual Pareto front in decision making. A Pareto set approximation is a similar approximation in the decision space. Several approximation methods can be found in the literature (see, among many others, [9][10][11][12]), but only few of them are effective for the application to nonconvex problems and even fewer involve the surrogate set in the actual decision making process (see, e.g., [13,14]).
A subclass of the Pareto set approximation methods consists of the so-called continuation methods, which translate mathematically the intuition of continuity and regularity that is apparent to human eye even when visualizing also a relatively small set of Pareto optimal solutions. These methods, rather than generate an as large as possible set of Pareto optimal solutions, attempt to trace explicitly the Pareto curves and surfaces. Continuation can be performed either locally or globally. Local continuation consists in considering a given Pareto optimal solution x 0 and then attempt to trace a mesh of new Pareto optimal solutions in a small neighborhood of x 0 . Iterating this process starting from one of the newly found solutions, one can generate the local branch of the Pareto set passing through x 0 (see [10,[15][16][17], and the references therein). Global continuation, on the other hand, is intended to draw consistently connections among a given set of optimal solutions, in order to build a higher-dimensional structure, which is usually a simplicial complex (see [18]), globally approximating the Pareto set. These ideas are reminiscent of "connecting the dots" children's games. Alternatively, a global continuation could start from a geometrical decomposition of the feasible set (e.g., a Delaunay tessellation, a cubic tessellation) and next produce a linear approximation of the portion of the Pareto set passing through each tessellation cell [12,19,20]. The strengths and weaknesses of these methods emerge when dealing with highly nonlinear and nonconvex problems, where the Pareto front can be disconnected, or composed by different local branches intersecting one another, as illustrated in Fig. 1.
By keeping these problems in mind, we have developed a new PAINT-SiCon global continuation method. The input of the PAINT-SiCon algorithm is a sampling of the decision space. The algorithm interpolates between the nondominated points of the sample both in the decision and in the objective space, respecting the domination structure of the objective space. Here, by respecting the domination structure, we mean that the constructed set is inherently nondominated as defined in [13]. Furthermore, the method detects connected components of the set of Pareto optimal solutions. The advantages of the method include that it is parameter free and especially developed to deal with nonconvex problems. Finally, the parametric representation of the approximation is computationally inexpensive to handle and can be viewed as a surrogate problem, like e.g., in [21]. Thus, the method facilitates decision making by embedding the use of the approximation in any interactive multiobjective optimization method.
The algorithm proposed in this paper draws ideas from two Pareto set approximation algorithms, i.e., PAINT [14,21] and SiCon [20]. The PAINT method differs from the new PAINT-SiCon in the fact that the former works only in the objective space, producing an interpolation of the Pareto  On the other hand, SiCon produces an approximation of all the local Pareto optimal solutions, which could possibly include some dominated points. In addition, SiCon requires the computation of the whole Delaunay tessellation of the sample points and derivative information on the objective functions and this, for high-dimensional problems, can be computationally exhaustive. The PAINT-SiCon algorithm proposed in this paper attempts to combine the strengths of both algorithms while avoiding their respective shortcomings. Especially, the new PAINT-SiCon method will be able to approximate both in the objective and decision spaces, detect disconnectedness in the Pareto set and is computationally less expensive than the SiCon method. In addition, the new PAINT-SiCon method includes a new way of including the knowledge about the constraints of the problem into the construction of the approximation, which is something that neither forebear method can do. This paper is organized as follows: First, in Sect. 2 we recall relevant mathematical background, fix the notation and give an overview of the two forebear methods PAINT and SiCon (focussing on details about their strengths and weaknesses). In addition, the paper contains an electronic supplementary material, where we review terminology related to computational and differential geometry that is necessary for this paper, but is not often discussed in multiobjective optimization literature. In Sect. 3, the mathematical basis for the new PAINT-SiCon method is established, while the technical details of the algorithm are given in Sect. 4. Section 5 describes the workings of the method on two critical nonconvex problems. Conclusions and perspectives are presented in Sect. 6.

Preliminaries
In this section, we review the relevant mathematical background from the perspective of this paper and give an overview of the two forebear approximation methods PAINT and SiCon. Further background material, specifically related to computational and differential geometry that are not often discussed in multiobjective optimization literature, is reported in Sect. 1 of the Supplementary Material (SM). In the Supplementary Material, we review the definitions that are needed in this research and set the notation and terminology. We start with the basic definitions of multiobjective optimization.

Multiobjective optimization
Let X ⊆ R n and consider a vector function -A general multiobjective optimization problem (MOP) consists in optimizing (maximizing or minimizing) jointly the scalar functions f 1 , . . . , f k : X −→ R: -The domain X of the functions is called the feasible set, while any point x ∈ X is called a (feasible) solution. The functions f 1 , . . . , f k : X −→ R are called the objective functions or simply the objectives. The space R n is called the decision space, while R k is called the objective space. The vectors z = f (x) ∈ R k , for any x ∈ X , are called outcomes.
Apart from pathological cases, it is not possible to find a solution which is optimal for all objectives taken separately. On the other hand, adopting a natural partial ordering on the objective space R k , it is usually possible to find infinite maximal elements, which in this context are named Pareto optimal solutions.
f i (y) for all i = 1, . . . , k, and f j (x) > f j (y) for at least one j ∈ {1, . . . , k}. The relation of domination induces a partial ordering on X .
-A solution x ∈ X is said Pareto optimal if there does not exist any solution which dominates x . The image f (x ) of a Pareto optimal solution x is called a Pareto optimal outcome. -The set of Pareto optimal solutions is called the Pareto set (PS), while the set of Pareto optimal outcomes is called the Pareto front (PF). Therefore PS is a subset of the decision space R n and of the feasible set X , while PF is a subset of the objective space R k .

PAINT: PAreto front INTerpolation
PAINT is a continuation method operating in the objective space. The PAINT method was proposed in [21], and it is based on the concept of an inherently nondominated PF approximation introduced in [13] and on the mathematical theory presented in [14]. PAINT defines a polyhedral complex using as nodes a given set of Pareto optimal outcomes, in a way that, by construction, there do not exist two vectors in the whole complex one dominating the other. As a first step, the method produces the full Delaunay tessellation based on the given Pareto optimal outcomes, and defines a basic approximation of the PF as the set of the given solutions. Next the method selects sequentially simplexes from the Delaunay tessellation that can be added to the approximation of the PF without impairing its inherent approximation. The resulting approximation of the PF has been proven accurate for a number of test problems. For instance, Fig. 2 shows an example of applying the PAINT method to approximate the PF of the standard 3-objective DTLZ2 problem [22]. The dots in the figure are the given input to the PAINT method, and the triangles and the quadrilaterals are the approximation constructed by PAINT. The actual PF of the DTLZ2 problem is the part of the ball surface that is in the positive orthand. In this case, the PAINT method correctly traces the shape of the PF. The main advantage of PAINT is that the approximation constructed can effectively support the decision making process. The surrogate problem is a multiobjective mixed-integer linear problem, therefore is almost computationally inexpensive and can be solved with any standard method. 1 Solving the approximated problem gives a number of candidate Pareto optimal solutions that can be appraised and compared by the decision maker. The selected approximate solutions can be used as reference points for finding the final preferred solution for the original problem, applying for instance the achievement scalarizing problem [23]. The analysis of the approximation can guide the experimenter for a deeper investigation of the problem, involving the computation of brand new Pareto optimal solutions in particularly interesting regions of the PF. Successively the PAINT method can be applied to the augmented set of solutions. 2 Unfortunately, as already noted in [21], the functional relation between decision variables and the objective values is not involved in the construction of the approximation. Therefore, the behavior of the approximation structure in the decision space can be unpredictable in the presence of nonconvexity. Figure 3a illustrates how this can lead the PAINT method to fail to recognize disconnectedness in the PS. The right-side panel in the Fig. 3a shows the outcomes given a sample of points in the objective space (the dots), where the nondominated ones have been given a red color. The red and gray line segments represent the PF approximation constructed by PAINT. However, visualizing the corresponding connecting segments in the decision space, it is clear that some of them (the grey lines) are unreasonable, because their extrema are too far apart. Indeed, such connections have to cross cells which nodes are dominated. Based only on dominance relations in the objectives space, PAINT method has produced a unique connected component for the PS, while clearly the true PS of the problem has three separate components.
Furthermore, although the problem in Fig. 3 is an unconstrained problem, there still appears disconnectedness in the PS that PAINT cannot detect. Even more severe issues may emerge if the multiobjective optimization problem is constrained as the PAINT method does not deal with constraints when defining the approximation. We will show in the next section how PAINT-SiCon can cope with constrained problems with disconnected PSs.

SiCon: Singular continuation
Singular Continuation [19,20], referred to as SiCon, is a method setting forth from the characterization of Pareto optimality based on first and second order derivatives of the objective functions given by Smale in [24] (see also [25,26]) and further using piece-wise linear approximations of implicitly defined manifolds (see, e.g., [27]). According to the SiCon method, linearized algebraic equations of the PS are determined for all the simplexes of a Delaunay tessellation of the feasible set based on a given sample of points. The solution of these equations is a polygon or a polyhedron representing a linear approximation of the part of the PS contained in the tessellation simplex. Taking the union of the simplexes gives a global approximation of the set of local Pareto optimal solutions. The outcome is a quadratically precise approximation of the PS, and the functional relation between PS and PF guarantees that the approximation is topologically consistent.
The method can approximate both the Pareto critical set 3 or the set of local Pareto optima. This can be a drawback, because there is the possibility that some globally dominated solutions to be included in the PS. Furthermore, the method is quite computationally intensive,

The mathematical basis of PAINT-SiCon method
The new method we propose here is an extension of the PAINT method aimed at reproducing the results of SiCon, while not being so computationally expensive, not requiring derivative information and also handling the constraints of the problem. The main issue about PAINT, which hinders consistent representations of the PS in the nonconvex case, is the missing connection between the decision space and the objective space.
As a first step, we consider the structure (i.e., the simplicial complex) emerging from the PAINT method in the objective space and then backwards propagate it on the decision space, i.e., considering a simplex having two outcomes f (x 1 ) and f (x 2 ) of nondominated sample points x 1 , x 2 ∈ F (where F ⊂ X is a finite sample set) as vertices one may define a new simplex in the decision space having the sample points x 1 and x 2 as vertices. Similar backpropagation can be performed with simplexes of any dimension. Back-propagating a simplex is, thus, defined for the rest of the paper as forming the simplex, the vertices of which are the inverse images of the vertices of the original simplex. In this way, we can guarantee that the approximation in the decision space complies with the domination structure. However, as observed above, some of these connections can be unreasonable, therefore we need a compelling criterion for removing the invalid simplexes. We propose a criterion based on a Delaunay tessellation of the decision space. We will show that with a dense enough sampling of the feasible set, a subset of the Delaunay tessellation restricted to the feasible set X can be used to extract the connected components of the PS. Then, only the simplexes that are defined by sample points in a single connected component are validated. The criterion can be generalized for constrained problems, by using restricted Delaunay tessellations to the feasible set X ⊆ R n [28,29]). The simplexes of a restricted tessellation have always the circumcenter inside the restricting subset. Figure 5 illustrates the difference of the Delaunay tessellation constrained to a set and the normal Delaunay tessellation.
We start with some definitions and, then, we build such tessellations characterizing the connected components of PS.
Definition 1 (Simplicial covering) Consider a simplicial complex C. Given any X ⊂ X , consider the union of the simplexes in C intersecting X : We refer to C ∩ X as the simplicial covering of X (based on the simplicial complex C).
Definition 2 (Simplicial patch) Consider a simplicial complex C and a set X ⊆ X . The connected components of C ∩ X are referred to as the simplicial patches of the parts of X (based on the simplicial complex C).
By examining the simplicial covering of the Pareto critical set in Fig. 3b, i.e., the green triangles, one clearly notices that the simplicial patches can separate the different connected components of the PS. This is clearly a valuable feature that should be investigated more deeply. Can we say that this is true in general, or at least that this is a generic 4 property? More precisely we state: Definition 3 (Faithful simplicial covering) For a simplicial complex C, a simplicial covering C∩X of a set X ⊆ X is called faithful if every simplicial patch contains exactly one connected component of X .
Therefore, a faithful simplicial covering of a set X is able to separate its connected components. At this point we want to prove that under suitable hypothesis on the objective functions, a simplicial covering based on the set of nondominated points is a faithful simplicial covering for PS. In order to do that we need to assume some regularity on the functions in study.

Definition 4
We say that the vector function f : X −→ R k is sufficiently regular if 1. the function f is C ∞ , 2. the PS is a stratified set of dimension k − 1 5 and 3. there exists > 0 such that two distinct connected components of the PS are at least apart each one from the other.

Remark 1
The stratification of the global Pareto set PS, stated in point 2. of Definition 4, can be false in degenerate cases but we conjecture that, in generic cases, it is true. It has been proven that generically the local PS is stratified [26,31,32] but a proof for the global PS is still missing. In addition, a viable analytic criterion for checking this property is not available, yet. This issue has been discussed in [19,[24][25][26]33]. In addition, we are going to assume that the input points in the decision and objective spaces are in general position. Thus, all the polytopes are in fact simplexes. This can be done, because if this is not the case, then the points can be perturbed infinitesimally to guarantee general position, as is shown in [34].

Proposition 1 If f : X −→ R k is sufficiently regular there exists a faithful simplicial covering of PS.
Proof If > 0 is the minimum distance between two connected components of PS, then a simplicial covering where all the simplexes are strictly smaller than 2 separates the connected components of PS. Indeed, consider such a simplicial covering and two points p, q ∈ PS belonging to adjacent simplexes, i.e., to simplexes having at least a common node. Then, if the diameters of the simplexes are all smaller than 2 , we have that | p − q| < 2 2 = . Therefore p and q cannot belong to different connected components of PS. As a result, any two points belonging to distinct connected components cannot belong to adjacent simplexes, and the covering of the PS is faithful.
Since the PS is not known, we need now some recipe for defining such coverings for an a priori not precisely known set X . In particular we want to know if the simplicial covering of an approximation of X covers faithfully the unknown approximated set X . For that we define the following: Definition 5 Given a set X ⊆ X and a finite set of points Λ = x 1 , . . . , x a , we say that Λ is a δ-approximation of X if max max i.e., d H (X , Λ) is the Hausdorff distance between the sets X and Λ.
With the above definition we can give sufficient conditions for a simplicial covering to be faithful. This is given by the following theorem.
Theorem 1 Let X ⊆ X and let δ > 0 be the minimum Hausdorff distance between two connected components of X . Let Λ be a sample of points where the maximum diameter of an empty ball is δ > 0, and let D be a Delaunay tessellation based on Λ. Let F ⊆ Λ be the points at a distance smaller than δ from a point of X . Then F is a δ -approximation of X and the subset of simplexes of D having all vertices in F is a simplicial covering of X . If δ < δ 2 , then the simplicial covering is also faithful for X .
Proof The proof can be divided in the following steps: 1. The maximum diameter of the Delaunay tessellation cells is smaller or equal than δ . 2. F is actually a δ -approximation of X . Indeed it is clear that for each point s ∈ X there exists a point v in Λ, and therefore in F, at a distance smaller that δ from s. Assume by absurdum that for a point s ∈ X there is no point of Λ in a ball of radius δ . The existence of a ball with radius δ not intersecting Λ is in contradiction with the hypothesis on Λ. 3. If a simplex τ in the tessellation intersects X , then its vertices are in X . Indeed every point in τ is within a δ distance from each of its vertices. Therefore, C, i.e., the set of simplexes with all vertices in X , is a simplicial covering of X . 4. Consider two points from two distinct connected components in X , s and q. Their respective distance is larger than δ. If there exists two adjacent simplexes, one containing s and one containing q, there should exist at least a common vertex v for two simplexes, with v ∈ X . But for the triangle inequality it should be which cannot hold if δ > 2δ . Therefore C is a faithful simplicial covering of X .
On the basis of Theorem 1, we claim that the connected components of the PS of a generic objective vector function can be distinguished consistently by a simplicial covering based on the nondominated set of a sufficiently fine approximation.

Theorem 2 Let G be a grid of points in X with maximum diameter of an empty ball equal to δ, and let F be the subset of nondominated points. If for every Pareto optimal solution there is a nondominated point in F at a distance smaller than δ > 0, and the minimum separation distance among two connected components of PS is larger than 5δ, then the set of Delaunay simplexes having at least a node in F is a faithful simplicial covering of PS.
Proof Although, we assume that every point in PS can be approximated by a nondominated point in the sample G, i.e., by a point in F, we cannot say that all points of G in a δneighborhood of PS are nondominated. In practice this is what one observes, e.g., Fig. 7. So F is not a δ-approximation of PS, however, in our hypothesis, we have that every point in PS is contained in a simplex with at least a vertex in F. Therefore, if the connected components of PS are distant at least 4δ, we have that the simplicial covering C composed by the simplexes having at least one vertex in F is faithful for PS.
Indeed, simplicial covering C must be contained in a tubular neighborhood of PS of radius 2δ 6 Indeed, C is composed by simplexes intersecting PS and some of the empty simplexes sharing at least one vertex with a nonempty simplex. If the sizes of those simplexes is smaller than δ, the farthest point from PS is at a distance smaller than 2δ. Therefore, the tubular neighborhoods of distinct connected components cannot intersect, and the covering is faithful. Remark 2 Thus, by building a Delaunay tessellation of the decision space and selecting the simplexes having at least a nondominated vertex, we have a correct representation of the connected components of the PS. Therefore, we could consider the simplexes defined by the PAINT method and validate only those contained in the simplicial patches of the covering. However, this strategy requires the full computation of the Delaunay tessellation, which is likely to be exhaustive in the interesting cases. Thus, we propose to detect the connected components by defining a graph which nodes are the nondominated points and which line segments (also called links and one-dimensional polytopes, in this context) are Delaunay polytopes. This structure is considerably less expensive than the construction of the full tessellation. At this point, if one validates only the Delaunay polytopes among the PAINT simplexes, practical experience suggests that a great number of valid simplexes would be rejected. We propose then to accept all PAINT simplexes whose extrema belong to the same connected components. This is illustrated more in detail in the next section.

The PAINT-SiCon method
The input of the PAINT-SiCon algorithm is a sampling Λ of the feasible set X and, if the sampling is fine enough and the functions are sufficiently regular, the algorithm outputs a PS approximation which is a simplicial complex and separates the connected components. We start by evaluating the objective functions on a given set of points, then we extract the nondominated set and build the approximation of the PS according to the PAINT method. We then back-propagate the approximation in to the decision space and validate the backpropagated simplexes by checking whether the vertices of the link are in the set of vertices of the same simplicial patch based on the Delaunay complex restricted to the feasible set S. The method described above is computationally less expensive than the SiCon method, but the computational burden is still not entirely negligible, if we require to compute the full Delaunay tessellation in the decision space. This difficulty can be overcome by checking the Delaunay condition only for the links we are interested. This choice considerably simplifies the procedure and makes the method efficient and appealing. Algorithm 1 describes the functionality of the PAINT method in more detail. Algorithm 2 partitions the Delaunay and similarly for higher order simplexes. 6: Divide the nondominated points into connected components C 1 , . . . , C m w.r.t. the Delaunay tessellation by calling Algorithm 2 with the nondominated points as input. Algorithm 2 Separating the Connected components of the PS 1: Consider as input the nondominated points given by Algorithm 1. 2: Set D 1,Λ X as such 1-dimensional Delaunay simplexes (based on the sampling Λ) restricted to X that both of the end points x 1 , x 2 ∈ Λ are nondominated in Λ i.e., and set componentnumber = 0 3: while D 1,Λ X = ∅ do 4: Set componentnumber = componentnumber + 1. 5: Select 1-dimensional simplex conv(x 1 , x 2 ) ∈ D 1,Λ X and set C componentnumber = {x 1 , x 2 } and set for all conv(x, y) ∈ D 1,Λ X do 8: if p ∈ V then 9: Set tessellation defined on the nondominated set F into connected components. A straightforward way to do step 1. of Algorithm 2 is to construct the Delaunay tessellation of the complete sample set Λ restricted to the feasible set S and then considering only the one-dimensional simplexes. After this, one can check whether both of the end points are nondominated points and, if so, include the simplex in D 1,Λ X . However, this approach is computationally expensive, as only the complexity of computing the restricted Delaunay tessellation in the decision space is exponential in the number of dimensions.
We have, however, developed a new optimization-based approach to finding the onedimensional Delaunay simplexes, which are needed for Algorithm 2. By definition, a simplex conv(x 1 , x 2 ) belongs to the restricted Delaunay tessellation [28] if and only if there exists an open ball B(x, r ) such that the center x ∈ X , B(x, r ) ∩ Λ = {x 1 , x 2 } and B(x, r ) ∩ Λ = ∅. Since the feasible set is often described in multiobjective optimization problems by mathematical constraints, a natural way of checking whether the simplex is in this Delaunay tesselation is by solving the optimization problem max t s.t.
If the optimal value of the problem (5) is greater than zero, then the simplex defined by these points is Delaunay; otherwise, it is not Delaunay. The problem can be written as x ∈ X and t ∈ R, (6) which is linear if the constraints setting X are linear. Thus, one needs to solve a(a −1) (where a is the number of nondominated points in the sample F) times problem (6) for finding out the restricted Delaunay one-dimensional simplexes of nondominated points in the sample. If the constraints are linear, the optimization problems can be solved with interior point methods in polynomial time [35]. Even if the case of nonlinear constraints, this approach still saves the computational cost of checking the Delaunay property.
After the above, the algorithm starts with a random simplex in the collection D X and first sets that the vertices of this one-dimensional simplex are in this component. Then, the algorithm starts inspecting which sample points are connectable (using the one-dimensional simplexes in D X ) from each of the sample points already in the sample. The whole component has been found, when there are no other sample points that are connected to sample points in the component with one-dimensional simplexes in D . Finally, the algorithm sets the nondominated sample points that are not connected to any other sample points as components that have merely one sample point. Figure 6 shows the Delaunay one-dimensional simplexes of nondominated points (the labeled points) for an arbitrarily defined sample set. In this example, the collection D X would thus contain one-dimensional simplexes conv(x 5 , x 14  Once the PAINT-SiCon approximation A has been constructed, it can be represented with almost the same surrogate problem as the PAINT surrogate. Let each simplex S ∈ A be represented as a row A i in a matrix A so that S = conv(x A i,1 , . . . , x A i,b ), where x i are solutions in the sample X and b ∈ N is an integer. With this matrix, the PAINT-SiCon surrogate problem can be written as In the above problem also the Pareto optimal solutions x are approximated and not only the Pareto optimal outcomes z. Thus, the decision maker does not have to find the corresponding solution to the actual problem to get information about the values of the decision variables. In addition, because the PAINT-SiCon method validates the polytopes in the decision space, the feasible values of the variable x in Problem (7) does not contain values that are far away from the connected components of the PS, and, because of continuity, the values of the objective functions are accurately approximated by the linearly approximated values given by the PAINT-SiCon surrogate problem.

Examples
In this section, examples of applying the PAINT-SiCon method are shown. Our test problems are the two-objective L&H 2×2 with two-dimensional decision space and the three-objective L&H 3×3 with three-dimensional decision space that have both been developed by A. Lovison. In both problems there is a PS composed by two distinct connected components, which give rise to two local PFs intersecting each other. For testing the algorithm, we have implemented it using the Octave [36] scripting language. For the sake of clarity, we give the explicit construction of these problems, which could be useful in order to construct problems with similar features. Basically we combine smooth bump functions, defined as Gaussians, and operate projections and rotations for orienting the fronts in the desired directions.L&H 2×2 is defined by summing two Gaussians, one centered in (0, 0), the other in (0, −1.5), the first bump with a larger amplitude, the second appreciably higher and sharper. The graph of the sum is projected on the (x, z) plane, such that the axes of the bumps are projected one above the other. Finally the plane is rotated by 45 • such that the axes of the Gaussians become oriented as the vector (1, 1). More formally: It can be seen that the results are qualitatively similar, but we notice that PAINT-SiCon did not made use of derivative information and did not require the computation of the full Delaunay tessellation of the decision space. Figure 8 illustrates the application of the PAINT-SiCon method to the L&H 3×3 problem. The top panels in Fig. 8 show the back-propagated PAINT simplexes (on the left) and the PAINT simplexes in the objective space (on the right). The bottom panels in the figure show the PAINT simplexes that have been validated both in the decision space (on the left) and in the objective space (on the right). For this problem also, the PAINT-SiCon method has correctly approximated the shape of the PF i.e., the approximation also consists of two connected components.
Surrogate problem (7) for this multiobjective optimization problem have parameter a = 13 describing the simplexes and b = 2 describing the maximal number of vertices in a simplex. 9 The matrix A characterizes the simplexes in a way described in Sect. 2.2 e.g., when the first simplex has as its vertices solutions x 17 and x 18 then the first line of the matrix A is [17,18]. Problem (7) is a mixed-integer linear multiobjective optimization problem and it can now be used as a surrogate of the original problem in e.g., interactive decision making or any a posteriori multiobjective optimization method.

Conclusions and topics of further research
In this paper, we have introduced the PAINT-SiCon PS approximation method that is an extension of the PAINT method using the ideas of the SiCon method. The method produces consistent parametric representations of the PS and the PF based on a sampling in the feasible set. These parametric representations can be used as a computationally inexpensive surrogate for the original multiobjective optimization problem. The method has been especially developed for nonconvex problems and can detect separate connected components of the PF unlike the PAINTmethod. In addition, the PAINT-SiCon method is computationally less expensive than the SiCon method.
As a topic of further research, great interest is in the definition of refinement strategies. The question remains if one can utilize the PS approximation constructed by the PAINT-SiCon method to produce new Pareto optimal solutions without building explicitly the Delaunay tessellation of the feasible set. Clearly, the points near the PS approximation are natural candidates for improving the approximation. Producing well distributed candidate points in the neighborhoods of the PAINT-SiCon structure, leading to more and more accurate approximations of the PS is still an open problem. Another open question is how to prevent getting stagnated in the neighborhoods of the first approximation detected, possibly missing a global branch of Pareto optimal points located in a still unexplored region of the feasible set. Clearly this question is a multiobjective version of the exploration-exploitation dilemma of global optimization.
Finally, when the method is used on real world applications, checking the smoothness conditions on such functions is usually impracticable. However, as noticed above, these conditions are not restrictive and we can expect that they are generally met, at least after an arbitrarily small perturbation.