Visualizing Exact and Approximated 3D Empirical Attainment Functions

Most real-world engineering optimization problems are inherently multiobjective, for example, searching for trade-off solutions of high quality and low cost. As no single optimal solution exists for such problems, multiobjective optimizers provide sets of optimal (or near-optimal) trade-off solutions to choose from.The empirical attainment function (EAF) is a powerful tool that can be used to analyze and compare the performance of these optimizers. While the visualization of EAFs is rather straightforward in two objectives, the three-objective case presents a great challenge as we need to visualize a large number of 3D cuboids. This paper addresses the visualization of exact as well as approximated 3D EAF values and differences in these values provided by two competing multiobjective optimizers. We show that the exact EAFs can be visualized using slicing and maximum intensity projection (MIP), while the approximated EAFs can be visualized using slicing, MIP, and direct volume rendering. In addition, the paper demonstrates the use of the proposed visualization techniques on a steel casting optimization problem.


Introduction
When solving engineering optimization problems it is usually not enough to find the best solution with regard to single measure of quality, but other objectives, such as other measures of quality, robustness of the solution or its cost, can be important too.While traditional approaches were used to combine all these heterogeneous objectives into a single one and solve the resulting singleobjective optimization problem, in the last two decades a lot of research has been directed into developing algorithms able to tackle these problems in their true multiobjective form.Most notably, multiobjective evolutionary algorithms (MOEAs) have proven to be effective in solving such problems [1,2].
A MOEA provides a set of trade-off solutions where no solution from the set is better than any other in all objectives.This is called an approximation set and MOEAs typically return a different approximation set in every run.If the performance of a MOEA needs to be analyzed or compared to the performance of another algorithm, several different measures can be adopted [3].Nevertheless, visualization of approximation sets can give an important insight into the properties of the algorithms (or the problem at hand) and is often used in this regard.While simple scatter plots can visualize 2D and 3D approximation sets (approximation sets of higher dimensions that require more sophisticated methods to be visualized will not be discussed in this paper-the interested reader is referred to [4] for a comprehensive review of such methods), scatter plots of multiple approximation sets can become too cluttered and loose their interpretation potential.The empirical attainment function (EAF), which assigns to each vector in the objective space a value signifying how often the vector was attained by the given approximation sets, can be used instead [5].Visualization of 2D EAFs entails plotting of rectangles in different colors (or shades of gray) corresponding to EAF values or differences in EAF values and is rather straightforward to do.
This paper focuses on the visualization of EAFs in the 3D case, which is more demanding than the 2D case since a large number of cuboids or rather unions of cuboids need to be first computed and then visualized.Building on some previous work [6,7], we present how the visualization of exact EAFs cuboids can be done using slicing (showing the 2D images obtained by slicing the cuboids at different angles) and maximum intensity projection (MIP) [8].However, if exactness is not crucial, the 3D objective space can be approximated using a grid of voxels, that is, values in a 3D grid.The paper introduces the idea of visualizing approximated EAFs using slicing, MIP, and direct volume rendering (DVR) [9].A flowchart explaining the connections among these methods is shown in Figure 1.Finally, we demonstrate the use of the proposed methods on steel casting optimization problem.The paper is organized as follows.Section 2 details the background and related work, while Section 3 introduces the benchmark approximation sets that are used in the rest of the paper.Visualization of exact and approximated EAFs is presented in Sections 4 and 5, respectively, and summarized in Section 6.The steel casting use case is depicted in Section 7. The paper concludes with final remarks in Section 8.

Background and Related Work
In addition to recalling some relations and concepts often used in multiobjective optimization, this section presents formal definitions of some new terms related to the EAF, such as attainment anchors and opposites.It also reviews the existing methods for visualization of EAFs and concludes with the presentation of slicing, MIP, and DVR.
As this paper deals with visualization in the objective space, which can be viewed rather independently from the decision space, the following definitions are confined to the objective space.

Definition 1 (Pareto dominance of vectors
All minimal elements of  constitute the minimal set of  denoted in this paper by  min . Maximal elements and the maximal set  max can be defined dually.Definition 9 (Pareto front).The set of Pareto optimal objective vectors known as the Pareto front can be formally defined as the minimal set of f(): Definition 10 (Edgeworth-Pareto hull).Let  ⊆  be an approximation set.The Edgeworth-Pareto hull (EPH) of  includes all objective vectors weakly dominated by :

Empirical Attainment Function.
Based on the multiobjective concept of goal-attainment [5] we say that an objective vector is attained when it is weakly dominated by an approximation set, that is, when it is contained in the set's EPH.The surface delimiting the attained vectors from the ones that have not been attained by an approximation set is called the attainment surface and its minimal elements are called the attainment anchors.
Definition 11 (set boundary).The boundary of a set  ⊆  is the intersection of the closure of , , with the closure of its complement, ( \ ): Definition 12 (attainment surface, attainment anchors).Let  ⊆  be an approximation set.The attainment surface of  is the boundary of the EPH of : The minimal elements of   are called attainment anchors and are equal to the original approximation set : Now, imagine that an optimization algorithm is run  times, producing  approximation sets.Then, each objective vector can be attained between 0 and  times.
Definition 13 (empirical attainment function).For algorithm A the empirical attainment function of the objective vector z ∈  ⊆  represents the frequency of attaining z by its approximation sets   1 , . . .,    : where I is the indicator function, defined as and ⪯ is the weak Pareto dominance relation between sets.
This means that for algorithm  with approximation sets   1 , . . .,    an EAF value from the set of frequencies {0, 1/, 2/, . . ., 1} can be assigned to every vector in the objective space.Of course, in practice the EAF cannot be computed for every objective vector and only attainment surfaces with a constant EAF value (and their corresponding attainment anchors) are of interest.They are called %attainment surfaces (also summary attainment surfaces) in order to distinguish them from the "ordinary" attainment surfaces of single approximation sets.Vectors in the %attainment surface are the tightest objective vectors that have been attained in at least % of the runs.
We will use  / to denote the attainment anchors of  / : Note how the definition of an attainment surface coincides with the definition of a summary attainment surface if  = 1.In general, the number of all attainment anchors (including duplicates),   = ∑  =1 | / |, is much larger than the number of vectors contained in the approximation sets,   = ∑  =1 |  | (in the context of the EAF computation [10], vectors from the approximation sets and attainment anchors were called input and output points, resp.).It has been shown in [10] that   ∈ ( 2   ) for three-objective optimization problems.
We illustrate these concepts in the two-objective case.Figure 2 shows an example of four approximation sets  1 , . . .,  4 (each consisted of four objective vectors indicated with crosses) and their summary attainment surfaces.Dots denote all attainment anchors and colors are used to emphasize areas with equal EAF values (darker colors correspond to higher values).If two algorithms  and  are run  times each, they can be compared by computing and visualizing their EAF values    and    separately.However, there exists a straightforward way to directly compare the algorithms by computing EAF differences, that is, differences in their EAF values [11].
Definition 15 (EAF differences).Assume that algorithms  and  are run  times each.The EAF differences between algorithms  and  are defined for each objective vector z ∈  as Note that EAF differences need to be computed for all attainment anchors of the combined 2 approximation sets.Defined in this way, the differences can adopt values from the set {−1, −( − 1)/, . . ., 0, 1/, . . ., 1}.Positive EAF differences correspond to areas in the objective space where the algorithm  outperforms the algorithm , while negative EAF differences denote areas in the objective space where the algorithm  outperforms the algorithm .Naturally, where the differences are zero, both algorithms attain the area equally well.
Let us focus on the example from Figure 3, where algorithms  and  solving a two-objective optimization problem produced two approximation sets each.The lines show the overall summary attainment surfaces.The colored areas emphasize areas of the objective space where algorithm  outperformed algorithm  (green hues) and areas where algorithm  outperformed algorithm  (red hues).We can readily notice that algorithm  is more successful in optimizing the first objective, while algorithm  attains better lower values of the second objective.This small example nevertheless shows that visualization of 2D EAF differences can be very helpful when analyzing and comparing the results of two algorithms.We can safely assume that we are interested in a limited portion of the objective space, for example, the one enclosed by reference vectors r 1 and r 2 , where r 1 7r 2 .Then, areas with equal EAF values/differences correspond to unions of rectangles in 2D and unions of cuboids in 3D [7].If we allow overlapping, then each rectangle or cuboid between the /and ( + 1)/-attainment surfaces, where  ∈ {1, . . ., −1}, can be defined by two vertices: one from the /attainment surface and the other from the (+1)/-attainment surface (Figure 4 shows 2D areas between 25%-and 50%attainment surfaces).While the lower vertices correspond to the attainment anchors of the /-attainment surface, the upper vertices are not the attainment anchors of the ( + 1)/attainment surface but their opposite.The opposite can be defined for an arbitrary approximation set as follows.
Definition 16 (approximation set opposite).Let reference vectors r 1 and r 2 , where r 1 7r 2 , define the boundaries of the observed objective space Let  ⊆  be an approximation set enclosed in this space with the attainment surface   .The opposite    of the approximation set  within the observed objective space  is the maximal set of   ∩ : Finding the opposites of attainment anchors in 2D is rather trivial [7].The computation of opposites (and the corresponding cuboids) for the 3D case is more demanding and is presented in detail in Section 4.1.
Finally, let us recall the formal definition of the hypervolume indicator [12], which is often used to measure the quality of approximation sets.
Definition 17 (hypervolume indicator).Let again reference vectors r 1 and r 2 , where r 1 7r 2 , define the boundaries of the observed objective space The hypervolume indicator   of an approximation set  ⊆  can be formulated via the empirical attainment function of  as 2.3.Visualization of EAFs.2D EAFs are most often visualized through best, median, and worst summary attainment surfaces corresponding to 1/%, ∼50%, and 100%-attainment surfaces, respectively.The first attempt to visualize 3D summary attainment surfaces is documented in [13].At the time, an efficient algorithm for computing 3D attainment anchors was not yet available, so separate summary attainment surfaces were approximated by discretizing the objective space using a grid.Assume the /-attainment surface needs to be visualized.The algorithm examines the objective space for each objective separately.First, it finds the intersections between all approximation sets and the lines stemming from the 2D grid of the remaining two objectives.Then, the intersections on each of these lines are counted and if they amount to , a marker is drawn at the corresponding height.Combining these markers by considering all objectives yields informative visualization of the chosen 3D summary attainment surface.Years later, an efficient algorithm for computing attainment anchors and their EAF values for the 3D case was finally provided [10,14].The attainment function tools from [14] were used to compute all EAF values in this paper.While [10] offers an estimation of the number of attainment anchors also for the general D case, where  > 3, no algorithm for computing 4D EAFs exists yet.As a consequence, we do not discuss visualization of EAFs beyond 3D.
Recently, the idea of visualizing EAF differences in addition to EAF values was introduced in [11].When visually comparing algorithms  and , four plots can be produced, each corresponding to one of the values of   ,   ,  − , and  − .All attainment anchors are plotted using grayshaded dots, where the shade corresponds to the value of the chosen function (  ,   ,  − , or  − ).On top of this, the best, median, and worst overall summary attainment surfaces are drawn.Examples from [11,15] show that this kind of visualization can be very useful in exploratory analysis.
Building on this work we have performed some initial research on visualizing 3D EAFs using slicing and maximum intensity projection in [6,7].These two methods, followed by direct volume rendering, are introduced next.

2.4.
Slicing.This is a simple method where spatial data is sliced using a plane.Then, the intersection of the plane with the data is visualized in 2D.In its most general form, slicing can be applied to any kind of spatial data; that is, the data does not have to be represented by voxels.
While the plane used in slicing is most often axis-aligned, the nature of multiobjective optimization problems yields the need for a plane that always contains one axis and the origin of the objective space (which usually corresponds to the ideal point) and intersects all summary attainment surfaces.Therefore (similarly to the prosection plane used in [4,16]), the slicing plane cuts the objective space at an angle  ∈ (0 ∘ , 90 ∘ ).
While not used for visualizing attainment surfaces, the interactive decision maps (IDMs) are related to slicing using axis-aligned planes [17].They visualize a number of axisaligned sampling surfaces of the EPH and use different colors to represent each slice.As the surface of  EPH is exactly the attainment surface of , it should be rather straightforward to apply this method for visualizing attainment surfaces.

Maximum Intensity Projection.
Maximum intensity projection (MIP) is a volume rendering method for spatial data represented by voxels.The method inspects voxels in direction of parallel rays traced from the viewpoint to the projection plane and takes the maximum value encountered in the voxels along a ray as the projection value for the ray.
The method, originally called maximum activity projection (MAP) [8], was proposed for 3D image rendering in nuclear medicine and tested in tomographic studies.It was later accepted not only in medical imaging, but also in scientific data visualization in general.
The advantages of MIP are its simplicity and efficiency and the ability of achieving high contrast, which arises from the fact that maximum voxel values are projected.On the other hand, as a limitation, the resulting projections lack the sense of depth of the original data (see [18,19] for attempts to remedy this issue).Moreover, the viewer cannot distinguish between left and right and front and back.As an improvement, animations are usually provided, consisting of a sequence of MIP renderings at slightly different viewpoints, which results in the illusion of rotation.[20] is another volume rendering method based on the voxel representation of data that is often used for visualization in medicine and science.It can generate very appealing results by employing advanced shading techniques and can achieve real-time interactive rendering.DVR is able to produce images of volumetric data without the need to explicitly extract geometry or surfaces (hence the name direct).

Direct Volume Rendering. Direct volume rendering (DVR)
First, each voxel value is assigned color and opacity (the RGBA value) by means of a chosen transfer function.The transfer function can be a simple ramp, a piecewise linear function, or an arbitrary table (the problem of specifying a good transfer function is a research area on its own [21,22]).Then, one of the following techniques is used to project the RGBA values to the screen: ray casting (for only ray casting is used in this paper, any mention of DVR can be understood to refer to ray-casting DVR from now on), splatting, shearwarp factorization, or texture-based volume rendering.In ray casting [9], viewing rays are traced from the viewpoint through the volume, accumulating color and opacity values

Benchmark Approximation Sets
In order to show the outcome of different visualization methods, we need some approximation sets to visualize.The two sets of approximation sets used in this paper are not products of optimization algorithms but rather benchmark approximation sets with known properties used for visualization purposes [7].They are shown in Figure 5.
The first set consists of five linear approximation sets with a uniform distribution of vectors that satisfy the constraint In order to simulate the behavior of a stochastic algorithm, the values   follow the normal distribution (1, 0.05).The second set contains five spherical approximation sets with a nonuniform distribution of vectors where only few vectors are located in the middle of the approximation set, while most of them are near its corners.The vectors from the spherical approximation sets satisfy the constraint where   ∼ (0.75, 0.05).Each individual approximation set contains 100 objective vectors.The values of   and   were chosen so that the sets are intertwined; in one region, the linear approximation sets dominate the spherical ones, while in others, the spherical sets dominate the linear ones.This assures that there will always be visible EAF differences between the two sets of approximation sets.We will use   5 and   5 to denote EAF values of the linear and spherical approximation sets, respectively.In addition,  − 5 and  − 5 will denote differences between those values.

Visualizing Exact EAFs
Exact EAFs can be represented as unions of cuboids with values corresponding to either EAF values of one algorithm or EAF differences between two algorithms.This section will first show the procedure for computing the cuboids from the given attainment surfaces and then present visualization of these cuboids using slicing and MIP.

Computing the Cuboids. The algorithm for computing
EAFs [10] takes as input a series of approximation sets  1 , . . .,   and returns as output a series of sets containing attainment anchors  1 , . . .,   that define the corresponding summary attainment surfaces and for which the following holds:   ⪯  +1 for  ∈ {1, 2, . . .,  − 1}.From these sets of attainment anchors and a reference vector r 2 that limits the observed objective space, the cuboids can be computed.We propose a very simple algorithm for this purpose that uses the first set of attainment anchors and the opposite of the second set of attainment anchors to compute overlapping cuboids (see Algorithm 1).
First, a set containing only the given reference vector r 2 is added at the end of the input set of attainment anchors, so that any cuboids ranging from the last set of attainment anchors to r 2 can also be computed.Then, the reference vector r 1 needed for computation of the opposites is set to be equal to the ideal point.Next, the main loop iterates through all adjacent pairs of sets   and  +1 and computes the opposite of  +1 using Algorithm 2. For each attainment anchor z ∈   all dominated vectors from the opposite are stored in   .Each pair (z, o), where o ∈   , represents the two vertices required to define the cuboid.In addition, the cuboid is given a value: either the EAF value or the EAF difference in z.
Input: Approximation set  and reference vectors r 1  Note that there might be multiple cuboids stemming from the same attainment anchor (see the analog example of multiple rectangles in Figure 4)-we call this a union of cuboids.This is not a problem since all such overlapping cuboids have the same value (if they did not have the same value, another attainment anchor would have already split the cuboid).The result of Algorithm 1 is a set of cuboids (or unions of cuboids).
Next, we show with the help of Figure 6 how to find the opposite of a 3D approximation set.Imagine a single cuboid defined by the reference vectors r 1 and r 2 .The opposite is initialized to {r 2 }.Every time a vector from  is "cut into" the existing cuboid, the vectors strictly dominated by it are removed from the opposite and new vectors are added to the opposite.Assume we have already performed this step for vectors z 1 and z 2 (see Figure 6(a)) and vector z 3 is next in line (see Figure 6(b)).First, we delete from the opposite all vectors that are strictly dominated by the current vector z (this means that in our example we delete vector o 4 but not o 5 ).Next, for every deleted vector o we create three new vectors in the following way: Finally, each of these vectors is added to the opposite if it is not redundant.This means that it has to contribute to the opposite (it cannot be coplanar with r 1 or collinear with any of the vectors from the existing opposite).In our example, we add to the opposite vectors o 6 and o 7 but not o 8 = (2, 6.5, 5) (not pictured in the figure) as it is collinear with o 5 and thus redundant.When these steps have been taken for every vector from , the resulting set  represents the opposite of .See Algorithm 2 for the algorithmic notation of the described procedure.As we cannot visualize the entire (infinite) objective space, Algorithms 1 and 2 use two reference vectors r 1 and r 2 to limit it.In order to preserve all information, they need to be set in such a way that r 1 ⪯ z7r 2 for all attainment anchors z.While the ideal point is the most sensible choice for r 1 , we can choose any objective vector dominated by all attainment anchors for r 2 .
The presented approach for computing the cuboids is very simple and easy to implement but comes at high computational cost.Computing the cuboids between two sets of  attainment anchors has the worst-case computational complexity of ( 3 ), which means that all cuboids between  pairs of attainment surfaces can be computed in ( 3 ) time.However, in practice, the loop foreach o ∈   do in Algorithm 2 is executed only a few times for each z ∈ , meaning that in practice the computational complexity is quadratic in the attainment anchors set size (the check for redundancy still takes linear time).
Another possible approach would be to deal only with nonoverlapping cuboids, which is somewhat similar to the problem of computing the hypervolume indicator in 3D.
In future work we wish to investigate if a more efficient algorithm for our problem could be found by, for example, modifying the space-sweep algorithm from [24], which computes the hypervolume indicator in the 3D case with the computational complexity of only ( log ).If this would be possible, the total cost of computing the cuboids would be ( log ).However, time complexity is not the only important aspect in the computation of cuboids.We would also like to study which method would yield a lower number of cuboids as this considerably affects the tractability of their visualization.

Visualizing Cuboids Using Slicing.
In this paper, the slicing plane is aligned with one of the three axes (assume for now that this is  3 ), includes the origin r 1 , and cuts through the underlying plane ( 1  2 ) at some angle  ∈ (0 ∘ , 90 ∘ ).In this way, each slice always captures all attainment surfaces.Slicing through the objective space containing a large number of cuboids means that only those cuboids that intersect the slicing plane are visualized.When a cuboid is sliced, the result is a rectangle in 3D (see the example in Figure 7).Two steps are needed to compute this rectangle in 2D.First, we need to calculate the projected cuboid vertices z  and o  yielding the rectangle in 3D, and second, we need to rotate the vertices by angle − to get z  and o  .Depending on the angles this is done in the following way: When slicing through a union of cuboids stemming from the same vertex z, some of the projected vertices o  can be redundant.In order to decrease the number of total rectangles to plot, it is therefore sensible to keep only nonredundant vertices o  (those that are nondominated with regard to the inverted weak dominance relation ⪰).
The cuboids can be visualized by slicing them at different angles, thus showing different parts of the objective space.We present visualization using slicing of our benchmark approximation sets at angles  = 5 ∘ and  = 45 ∘ .The EAF values   are presented on the same plot (see Figure 8).From the plots of EAF values it is easy to distinguish linear sets (green hues) from the spherical ones (red hues) as the shape of the sets is readily visible.We can also see that the two best spherical approximation sets are better than the remaining three by a solid margin.While these plots provide a lot of information by themselves, the comparison between the sets is best visualized using EAF differences.They nicely show regions in the objective space where linear sets outperform the spherical ones and vice versa.
Note that in order to fully explore the entire objective space, slicing planes at several different angles should be performed.Also, one might want to consider slicing planes aligned to the axis  1 or  2 , too.However, as our benchmark approximation sets are rather symmetrical, there is no need to do this here.

Visualizing Cuboids Using MIP.
As explained in Section 2.5, MIP shows only the maximum value encountered on rays from the viewpoint to the projection plane.This means that it is not sensible to use this method for visualizing EAF values as most of the objective space has the maximum value and the resulting visualization would not be very informative.However, MIP seems a good choice for visualizing EAF differences where the highest values are of most interest as they represent the largest differences between the algorithms.
MIP could be used to visualize cuboids of different values by sorting the cuboids so that those with higher values would be put on top of those with lower values.However in general, the 3D plotting tools capable of visualizing cuboids render them in a sequence that tries to maintain some notion of depth (cuboids near the viewpoint are shown in front of the cuboids further away) and do not allow for custom sorting of the cuboids according to their values.Therefore, this sorting can only be done after the viewpoint has been set and the cuboids are already projected to 2D.At that time we can choose to visualize them in the order of ascending values, which (although a bit cumbersome) effectively achieves the MIP visualization of the cuboids.
Another difficulty in using MIP for visualizing exact EAF differences is that we need to visualize a large number of cuboids, which can be challenging to do.While many of them are completely covered by cuboids with larger values and could therefore be spared without altering the final image, removal of such cuboids is not trivial as it depends also on the position of the viewpoint.
Nevertheless, we present visualization of exact EAF differences using MIP in Figure 9.The differences  − 5 and  − 5 need to be visualized separately in order to avoid overlapping of cuboids.The MIP visualizations are very informative; we can easily see which parts of the objective space are better attained by the linear approximation sets and which by the spherical approximation sets.As is typical with MIP, while we are able to "see through the cuboids, " we loose the sense of depth.Although it is inevitable to lose some information when projecting 3D structures onto 2D, this can be amended by combining two visualization techniques: MIP and slicing.Together they give a good idea of the 3D "cloud" of cuboids.

Visualizing Approximated EAFs
If we wish to visually compare the outcome of two algorithms and are not interested in the details of such a visualization, we can use approximated instead of exact EAFs.Approximation means that the objective space is discretized into a grid of voxels (similarly to what was proposed in [13]).This section first presents how such discretization is performed and then illustrates visualization of approximated EAFs using slicing, MIP, and DVR.

Discretization into Voxels.
A voxel is a single point on a regularly spaced 3D grid of V 1 × V 2 × V 3 voxels.Recall that, based on the reference vectors r 1 and r 2 , the 3D observed objective space  is defined as If  is a cube, V 1 = V 2 = V 3 ; otherwise, some care must be taken to ensure that the grid of voxels is truly regular.Either  must be normalized prior to approximation or the number of voxels in each dimension V  must be set to be proportional to  2  −  1  for all  ∈ {1, 2, 3}.Note that a voxel represents only a single point not a volume.How this missing information is reconstructed depends on the visualization method.
Computing voxel values when the cuboids have already been computed is rather straightforward.Iterating over all cuboids, the voxels that are "covered" by the cuboid adopt its value.If we are not interested in visualizing the exact EAFs (and therefore have not computed the cuboids), we can compute the voxel values directly from the set of approximation sets.The value of each voxel is set to either EAF value or EAF difference by counting how many approximation sets weakly dominate it.

Visualizing Voxels Using Slicing.
For visualization using slicing we assume that the whole volume of the voxel has the same value as its center.Therefore, slicing of voxels is done in the same way as slicing of cuboids (see Section 4.2).In order to avoid showing plots that are very similar to those presented in Figure 8, Figure 10 shows only the slices of  − at angle  = 5 ∘ produced using different discretizations.
We can see that by refining the discretization we get results that are increasingly similar to the exact EAFs.As we believe that for the purpose of this study the discretization into 128 3 voxels suffices, we will use this discretization in the remainder of the paper.

Visualizing Voxels Using MIP.
Visualizing voxels using MIP is trivial with a tool supporting such visualization, as there are no additional parameters to set.We use a volume rendering engine called Voreen [25,26] to produce the MIP images from Figure 11 and all DVR images.and we can see that although these are approximations of the images presented in Figure 9, the results are quite similar.

Visualizing Voxels Using DVR.
Visualization using DVR is a bit trickier as we need to define the transfer function that assigns color and opacity to each voxel value.In our case, voxel values are discrete as they equal either the EAF values or the EAF differences, which makes this task easier.
Figures 12 and 13 show visualization using DVR of the EAF differences between the two benchmark approximation sets  − 5 and  − 5 , respectively.The first five plots in both figures show the voxels with a single value from {1/5, . . ., 5/5}.The transfer functions used to obtain these plots are simple piecewise linear functions that set the opacity of voxels of the desired value to 1 and the opacity of all other voxels to 0. In this way, we are able to visualize each  of the values separately, which is reflected in the nice and informative plots in Figures 12 and 13.The biggest advantage of DVR is that, by setting the transfer function "the right way, " it is possible to see inside the visualized volume.This was attempted on the plots shown in Figures 12(f) and 13(f), which show completely opaque voxels with value 5/5 and almost transparent voxels with value 1/5 (both plots are rotated to give a nicer view of inner voxels).
In contrast to MIP, DVR can be used for visualizing EAF values, too. Figure 14 shows an example of such a visualization for the spherical approximation sets, where only the values   5 ∈ {1/5, 4/5} are shown (opacity is set to 0.3 for these two values and to 0 for all other values).

Discussion
Let us summarize the properties of all presented visualization methods (see also Table 1).Slicing is a rather simple method that enables visualization of exact EAF values and differences.Its biggest advantage is its accuracy; it enables a detailed analysis of the approximation sets in the same way that can be done for 2D EAFs.Its biggest drawback is that it is able to visualize only one slice at a time.Generally, the approximation sets are not as symmetrical as in our benchmark case; therefore, multiple slices are needed to sufficiently inspect the EAFs.Also, note that the meaning of the angle  changes if the observed objective space has different ranges  (if  is not a cube, it is not cut exactly in half by the plane at angle  = 45 ∘ ).While slicing can be used also for visualizing approximated EAFs, there is no need to do so as it works well for exact EAFs.
It is not reasonable to use MIP for visualizing EAF values since usually a large portion of the objective space has the maximum EAF value and such plots are not very informative.While MIP can be used for visualizing exact EAF differences, it works much better on the approximated case, for which it was conceived.Visualizing 3D cuboids is rather impractical (especially for a very large number of cuboids) and requires sorting of cuboids with regard to their value in order to produce MIP image.However, MIP can easily be used to visualize approximated EAF differences as it has no parameters to set.Its biggest advantage is that it combines all values in a single image giving precedence to largest differences (which are the most important when visualizing EAF differences), while its the biggest disadvantage is the lost sense of depth.
Finally, DVR can be used to visualize approximated EAFs.It produces nice and informative visualizations of both the EAF values and differences.With an insightful setting of the transfer function it is even possible to "look through" the cloud of cuboids.The need to define the transfer function is at the same time the biggest disadvantage of DVR as it might be demanding for a user not familiar with the method.Nevertheless, it is much easier to find the right transfer function for EAFs than, for example, for some medical data, as EAFs take only discrete values.

Steel Casting Use Case
This section shows how the described visualization methods can be used on a real-world optimization problem with three objectives from a previous study [27].

The Steel Casting Problem.
Continuous casting of steel is a complex metallurgical process where molten steel is cooled and shaped into semimanufactures of desired dimensions.Cooling is done using water in the primary and secondary cooling subsystems.Primary cooling is performed in the mold, while secondary cooling comprehends wreath cooling at the exist from the mold, and spray cooling of the strand.
The goal is to set the parameters of this process (casting speed, mold outlet coolant temperature, and wreath and spray coolant flows (see Table 2)) in such a way that the quality of the cast steel is as high as possible.The quality of steel is defined by the following three objectives: distance from desired metallurgical length (the length of the liquid core in the strand), distance from desired shell thickness (thickness of the solid shell at mold exit), and distance from desired surface temperature at a predefined point in the strand.
where   is the value of the observed variable and  *  is its desired value.
As real-world experimentation with parameter settings is expensive and time-consuming and could also be dangerous, we have simulated it using a numerical model of steel casting based on a meshless technique for diffusive heat transport [28].One simulation of the steel casting process takes approximately 2 minutes on a standard desktop computer.
Two instances of this optimization problem were studied: discrete and continuous.The discrete problem instance was set by domain experts as presented in Table 2 in such a way that all possible solutions (9639 in total) could be explored.While the discrete problem instance enables the rough exploration of the objective space, the continuous problem instance where any value within the variable bounds could be chosen is the one we wish to solve.

Results of Optimization Algorithms.
The discrete problem instance was solved using the exhaustive search (ES) algorithm.These solutions form the Pareto front of the discrete problem instance and we wanted to see whether an algorithm tackling the continuous problem instance could come close to this Pareto front.
We chose the algorithm DEMO (differential evolution for multiobjective optimization) [29], which is a MOEA algorithm that uses differential evolution [30], to explore the decision space for the continuous problem instance.DEMO was run five times, each time exploring 3200 solutions.More detailed information on the employed experimental setup can be found in [27].The majority of the explored solutions were infeasible.For example, out of 9639 solutions found by ES only 1242 were feasible and of those only 72 were mutually nondominated.DEMO was solving an extended problem and therefore found more nondominated solutions: on average, 645 per run (although their number varied significantly over different runs).
First, we visualize all final (feasible and nondominated) solutions found by both algorithms in Figure 15.We can see that ES found two distinct subsets of solutions.A more detailed analysis of results revealed that they correspond to two of the three mold outlet coolant temperatures (the remaining one always produces infeasible solutions).Since DEMO was not bound by this discretization, it was able to find solutions also around these subsets.Both algorithms found a few solutions with low distances from desired shell thickness that are rather "detached" from the rest.They actually lie on a larger disconnected region of the Pareto front of which just this minor part is feasible.
This visualization is able to show that DEMO is able to reach the Pareto front of the discrete problem and is therefore a good choice for solving this problem.However, there are two aspects we are interested in that are hard to infer from this visualization alone (they could, of course, be computed from the solutions).First, we wish to see whether different runs of DEMO produce similar results.Visualizing all five approximation sets together with different markers or visualizing one approximation set at a time does not provide a good means for comparison.Second, although the approximation sets found by DEMO look linear, they are in fact slightly convex making it very hard to see (even by rotating the objective space) whether results by ES are in fact dominated by those by DEMO.We will try to find the answers to these questions using the proposed visualization methods.
In order to use the proposed visualization methods to visualize the results of ES and DEMO, we need to address two issues.The first issue is the uneven ranges of the attainment surfaces, which are important when computing voxels.As we can infer from Table 3, the observed objective space is equal to [0, 1] × [0, 2] × [0, 7.5].Moreover, the feasible results found by the two algorithms actually lie in an even smaller cuboid contained in [0, 0.6] × [0.5, 2] × [0, 7.5].If the number of voxels was proportional to the objective ranges, the first two objectives would be discretized too roughly.Therefore, we choose to normalize the objective vectors with regard to [0, 0.6] × [0.5, 2] × [0, 7.5].The second issue is the uneven number of runs of both algorithms (ES was run once).To leave the meaning of the EAFs differences intact, we copy the results by ES to get five runs in total.This seems a reasonable approach as the deterministic ES would actually produce five equal results if it was run five times.Now, visualization methods as described in Sections 4 and 5 can be used on these results.Figure 16 presents the DVR of approximated EAF differences between the two algorithms,  DEMO−ES 5 and  ES−DEMO

5
, which gives us a general idea of their performance.However, this does not answer our two questions.To check the repeatability of DEMO, we need to inspect the difference between the best and worst  17) suggests that although DEMO's approximation sets were of considerably different cardinality, they attain a similar portion of the objective space.This is further confirmed by a more detailed inspection using slicing at angles  = 25 ∘ and 45 ∘ (see Figure 18), which shows that only a small border of the attained objective space is attained less than five times (denoted by light green hues).
Finally, MIP can be used to check whether DEMO ever finds better solutions than those found by ES.If this was not the case, each solution found by ES would have exactly one union of cuboids (albeit small) for which  ES−DEMO 5 = 5/5.As accuracy is important in this case (small cuboids can be omitted if the discretization into voxels is not fine enough), we use MIP on exact EAF differences.Figure 19 clearly shows  that although  ES−DEMO 5 = 5/5 for some solutions, this does not hold for all of them, meaning that DEMO was actually able to find solutions that dominate those by ES.This was of course possible only because DEMO was solving an extended instance of the problem.

Conclusions
When comparing the results of multiobjective optimization algorithms it is important to be able to visualize them as this can provide new information regarding the algorithms or the given problem.If the algorithms are stochastic, the EAF can be used to describe how well the algorithms attain the objective space with their multiple approximation sets.Visualization of EAFs is rather straightforward in 2D but presents a challenge in 3D as multiple cuboids need to be visualized.This paper presented how these cuboids can be computed and visualized using slicing and MIP.If accuracy of visualization is not crucial, the EAF values and differences can be approximated by discretizing the objective space into a grid of voxels.In this way, visualization can be performed using slicing, MIP, and DVR.We have shown how all these methods perform on two sets of benchmark approximation sets and discussed their advantages and disadvantages.
In addition, we demonstrated the use of slicing, MIP, and DVR on a real-world optimization problem solved by exhaustive search and MOEA algorithm.We have shown that these powerful visualization methods are able to give a new insight regarding the performance of the algorithms that cannot be otherwise seen using solely "standard" visualization of approximation sets.
In the future, we wish to find a more efficient way to compute the cuboids either by improving the provided algorithm or by adjusting existing efficient algorithms for hypervolume calculation to suit our needs.

Figure 4 :
Figure 4:  The area between the 2D 25%-and 50%-attainment surfaces from the example in Figure2.Vectors from the opposite of the attainment anchors of the 50%-attainment surface are denoted by squares.

Figure 5 :
Figure 5: The linear and spherical approximation sets used in all visualization examples in Sections 4 and 5.

3 (b) Opposite of {z 1 , z 2 , z 3 }Figure 6 :
Figure 6: Two steps in the construction of the opposite of the approximation set .

Figure 7 :
Figure 7: A cuboid sliced using the plane aligned with  3 that slices the  1  2 plane under angle.The cuboid vertices z and o are projected to 3D rectangle vertexes z  and o  , respectively.

Figure 8 :
Figure 8: Slices of the exact 3D EAF values and differences for the benchmark approximation sets under two angles.Darker colors represent higher EAF values/differences.

5
and   5 are shown separately, while the EAF differences  − 5 and  − 5

5 Figure 9 :
Figure 9: MIP of exact 3D EAF differences between the benchmark approximation sets.

Figure 11
Figure 11 shows the MIP for  − 5 and  − 5

5 Figure 11 :
Figure 11: MIP of approximated 3D EAF differences between the benchmark approximation sets.

Figure 12 :
Figure 12: DVR of approximated 3D EAF differences between the linear and spherical benchmark approximation sets  − 5 .

Figure 15 :
Figure 15: The best solutions found by ES on the discrete problem instance and by DEMO on the continuous problem instance (all five runs shown).

Figure 18 :
Figure 18: Slices of exact 3D EAF values of DEMO at two angles.

Table 1 :
Summarized properties of the presented visualization methods.

Table 2 :
Steel casting process parameters.

Table 3 :
Variables defining the optimization objectives.

Table 3
details the variables defining these optimization objectives.Their bounds and desired values were determined by domain experts.If a parameter setting yields values outside the boundary constraints presented in Table3, it is deemed infeasible.The objectives of the feasible solutions are computed by