# Visualizing Exact and Approximated 3D Empirical  Attainment Functions

Authors: Tea Tušar; Bogdan Filipič
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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 anapproximation 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 ofcuboids 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 ofvoxels, 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.Figure 1
Flowchart of EAF visualization depending on the preferences.The paper is organized as follows. Section2 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.


## 2. 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 asattainment anchors andopposites. It also reviews the existing methods for visualization of EAFs and concludes with the presentation of slicing, MIP, and DVR.


### 2.1. Multiobjective Optimization Concepts

The multiobjective optimization problem consists of finding the optimum of a function:(1)f:X⟶F,f:(x1,…,xn)⟼(f1(x1,…,xn),…,fm(x1,…,xn)),
where X is an n-dimensional decision space and F is an m-dimensional objective space (m≥2). Each solution x∈X is called a decision vector, while the corresponding element z=f(x)∈F is an objective vector. Without loss of generality we assume that F⊆Rm and all objectives fi:X→R are to be minimized.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).
The objective vectorzA=(z1A,…,zmA)∈Fdominates the objective vector zB=(z1B,…,zmB)∈F; that is, zA≺zB, if
(2)ziA≤ziB∀i∈{1,…,m},thereexistsaj∈{1,…,m}forwhichzjA<zjB.Definition 2 (weak Pareto dominance of vectors).
The objective vectorzA=(z1A,…,zmA)∈Fweakly dominates the objective vector zB=(z1B,…,zmB)∈F; that is, zA⪯zB, if
(3)ziA≤ziB∀i∈{1,…,m}.Definition 3 (strict Pareto dominance of vectors).
The objective vectorzA=(z1A,…,zmA)∈Fstrictly dominates the objective vector zB=(z1B,…,zmB)∈F; that is, zA⪻zB, if
(4)ziA<ziB∀i∈{1,…,m}.Definition 4 (incomparability of vectors).
The objective vectorszA=(z1A,…,zmA)∈F and zB=(z1B,…,zmB)∈F areincomparable; that is, zA∥zB, if
(5)zA⋠zB,zB⋠zA.Definition 5 (approximation set).
A set of objective vectorsZ⊆F is called anapproximation set if zA∥zB for any two objective vectors zA,zB∈Z.Definition 6 (weak Pareto dominance of approximation sets).
The approximation setZA⊆Fweakly dominates the approximation set ZB⊆F; that is, ZA⪯ZB, if every zB∈ZB is weakly dominated by at least one zA∈ZA.Definition 7 (ideal point).
Theideal pointz*=(z1*,…,zm*) represents the minimum possible value in each objective (typically, it cannot be achieved):
(6)zi*=min⁡{fi(x);x∈X}∀i∈{1,…,m}.Definition 8 (minimal element, minimal set).
LetZ⊆F be a set of objective vectors. An objective vector m∈Z is aminimal element of Z if
(7)∀z∈Z,z⪯m⟹z=m.
All minimal elements of Z constitute theminimal set of Z denoted in this paper by Zmin⁡.Maximal elements and the maximal setZmax⁡ can be defined dually.Definition 9 (Pareto front).
The set of Pareto optimal objective vectors known as thePareto front can be formally defined as the minimal set of f(X):
(8)Pf=f(X)min⁡.Definition 10 (Edgeworth-Pareto hull).
LetZ⊆F be an approximation set. The Edgeworth-Pareto hull (EPH) of Z includes all objective vectors weakly dominated by Z:
(9)ZEPH={z∈F;existsz′∈Zsothatz′⪯z}.


### 2.2. Empirical Attainment Function

Based on the multiobjective concept of goal-attainment [5] we say that an objective vector isattained 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 theattainment surface and its minimal elements are called theattainment anchors.Definition 11 (set boundary).
The boundary of a setZ⊆F is the intersection of the closure of Z, Z¯, with the closure of its complement, (F∖Z)¯:
(10)∂Z=Z¯∩(F∖Z)¯.Definition 12 (attainment surface, attainment anchors).
LetZ⊆F be an approximation set. Theattainment surface of Z is the boundary of the EPH of Z:
(11)SZ=∂ZEPH.
The minimal elements of SZ are calledattainment anchors and are equal to the original approximation set Z:
(12)AZ=SZmin⁡=Z.
Now, imagine that an optimization algorithm is runr times, producing r approximation sets. Then, each objective vector can be attained between 0 and r times.Definition 13 (empirical attainment function).
For algorithm A the empirical attainment function of the objective vectorz∈Z⊆F represents the frequency of attaining z by its approximation sets Z1A,…,ZrA:
(13)αrA(z)=α(Z1A,…,ZrA;z)=1r∑i=1r‍I(ZiA⪯{z}),
where I is the indicator function, defined as
(14)I(b)={1,ifbistrue,0,otherwise,
and ⪯ is the weak Pareto dominance relation between sets.This means that for algorithmA with approximation sets Z1A,…,ZrA an EAF value from the set of frequencies {0,1/r,2/r,…,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 k%-attainment surfaces (alsosummary attainment surfaces) in order to distinguish them from the “ordinary” attainment surfaces of single approximation sets. Vectors in the k%-attainment surface are the tightest objective vectors that have been attained in at least k% of the runs.Definition 14 (summary attainment surfaces).
LetZ1,…,Zr be r approximation sets of algorithm A and αrA(z) its EAF. Let the t/r-superlevel set of αrA be defined as
(15)Lt/r+(αrA)={z∈F;αrA(z)≥t/r}∀t∈{1,…,r}.
Then thesummary (or t/r-)attainment surfaceSt/r is equal to the boundary of Lt/r+(αrA):
(16)St/r=∂Lt/r+(αrA)∀t∈{1,…,r}.
We will use At/r to denote the attainment anchors of St/r:
(17)At/r=St/rmin⁡∀t∈{1,…,r}.
Note how the definition of an attainment surface coincides with the definition of a summary attainment surface ifr=1. In general, the number of all attainment anchors (including duplicates), nA=∑t=1r‍|At/r|, is much larger than the number of vectors contained in the approximation sets, nZ=∑t=1r‍|Zt| (in the context of the EAF computation [10], vectors from the approximation sets and attainment anchors were calledinput andoutput points, resp.). It has been shown in [10] that nA∈O(r2nZ) for three-objective optimization problems.We illustrate these concepts in the two-objective case. Figure2 shows an example of four approximation sets Z1,…,Z4 (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).Figure 2
Four approximation sets in a two-objective optimization problem, their summary attainment surfaces, and all attainment anchors. Colors denote areas with equal EAF values (darker colors correspond to larger values).If two algorithmsA and B are run r times each, they can be compared by computing and visualizing their EAF values αrA and αrB separately. However, there exists a straightforward way to directly compare the algorithms by computingEAF differences, that is, differences in their EAF values [11].Definition 15 (EAF differences).
Assume that algorithmsA and B are run r times each. The EAF differences between algorithms A and B are defined for each objective vector z∈Z as
(18)δrA-B(z)=αrA(z)-αrB(z).
Note that EAF differences need to be computed for all attainment anchors of the combined2r approximation sets. Defined in this way, the differences can adopt values from the set {-1,-(r-1)/r,…,0,1/r,…,1}. Positive EAF differences correspond to areas in the objective space where the algorithm A outperforms the algorithm B, while negative EAF differences denote areas in the objective space where the algorithm B outperforms the algorithm A. Naturally, where the differences are zero, both algorithms attain the area equally well.Let us focus on the example from Figure3, where algorithms A and B 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 A outperformed algorithm B (green hues) and areas where algorithm B outperformed algorithm A (red hues). We can readily notice that algorithm A is more successful in optimizing the first objective, while algorithm B 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.Figure 3
Four approximation sets in a two-objective optimization problem: two produced by algorithmA and two by algorithm B. Colored areas show where algorithm A outperformed algorithm B (green hues) and algorithm B outperformed algorithm A (red hues).We can safely assume that we are interested in a limited portion of the objective space, for example, the one enclosed by reference vectorsr1 and r2, where r1⪻r2. 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 t/r- and (t+1)/r-attainment surfaces, where t∈{1,…,r-1}, can be defined by two vertices: one from the t/r-attainment surface and the other from the (t+1)/r-attainment surface (Figure 4 shows 2D areas between 25%- and 50%-attainment surfaces). While the lower vertices correspond to the attainment anchors of the t/r-attainment surface, the upper vertices are not the attainment anchors of the (t+1)/r-attainment surface but theiropposite. The opposite can be defined for an arbitrary approximation set as follows.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.Definition 16 (approximation set opposite).
Let reference vectorsr1 and r2, where r1⪻r2, define the boundaries of the observed objective space
(19)R={z∈F;zi∈[ri1,ri2]∀i∈{1,…,m}}.
Let Z⊆R be an approximation set enclosed in this space with the attainment surface SZ. TheoppositeOZR of the approximation set Z within the observed objective space R is the maximal set of SZ∩R:
(20)OZR=(SZ∩R)max⁡.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 thehypervolume indicator [12], which is often used to measure the quality of approximation sets.Definition 17 (hypervolume indicator).
Let again reference vectorsr1 and r2, where r1⪻r2, define the boundaries of the observed objective space
(21)R={z∈F;zi∈[ri1,ri2]∀i∈{1,…,m}}.
Thehypervolume indicatorIH of an approximation set Z⊆R can be formulated via the empirical attainment function of Z as
(22)IH(Z)=∫z∈R‍α(Z;z)dz.


### 2.3. Visualization of EAFs

2D EAFs are most often visualized throughbest,median, andworst summary attainment surfaces corresponding to 1/r%, ~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 t/r-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 t, 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 mD case, where m>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 A and B, four plots can be produced, each corresponding to one of the values of αA, αB, δA-B, and δB-A. All attainment anchors are plotted using gray-shaded dots, where the shade corresponds to the value of the chosen function (αA, αB, δA-B, or δB-A). 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 theprosection 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 axis-aligned sampling surfaces of the EPH and use different colors to represent each slice. As the surface of ZEPH is exactly the attainment surface of Z, it should be rather straightforward to apply this method for visualizing attainment surfaces.


### 2.5. 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.


### 2.6. Direct Volume Rendering

Direct volume rendering (DVR) [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 namedirect).First, each voxel value is assigned color and opacity (the RGBA value) by means of a chosentransfer 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, shear-warp 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 at each sample position along the ray. The final RGB values shown in the image are computed from the accumulated RGBA values using the volume rendering integral [23].


## 2.1. Multiobjective Optimization Concepts

The multiobjective optimization problem consists of finding the optimum of a function:(1)f:X⟶F,f:(x1,…,xn)⟼(f1(x1,…,xn),…,fm(x1,…,xn)),
where X is an n-dimensional decision space and F is an m-dimensional objective space (m≥2). Each solution x∈X is called a decision vector, while the corresponding element z=f(x)∈F is an objective vector. Without loss of generality we assume that F⊆Rm and all objectives fi:X→R are to be minimized.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).
The objective vectorzA=(z1A,…,zmA)∈Fdominates the objective vector zB=(z1B,…,zmB)∈F; that is, zA≺zB, if
(2)ziA≤ziB∀i∈{1,…,m},thereexistsaj∈{1,…,m}forwhichzjA<zjB.Definition 2 (weak Pareto dominance of vectors).
The objective vectorzA=(z1A,…,zmA)∈Fweakly dominates the objective vector zB=(z1B,…,zmB)∈F; that is, zA⪯zB, if
(3)ziA≤ziB∀i∈{1,…,m}.Definition 3 (strict Pareto dominance of vectors).
The objective vectorzA=(z1A,…,zmA)∈Fstrictly dominates the objective vector zB=(z1B,…,zmB)∈F; that is, zA⪻zB, if
(4)ziA<ziB∀i∈{1,…,m}.Definition 4 (incomparability of vectors).
The objective vectorszA=(z1A,…,zmA)∈F and zB=(z1B,…,zmB)∈F areincomparable; that is, zA∥zB, if
(5)zA⋠zB,zB⋠zA.Definition 5 (approximation set).
A set of objective vectorsZ⊆F is called anapproximation set if zA∥zB for any two objective vectors zA,zB∈Z.Definition 6 (weak Pareto dominance of approximation sets).
The approximation setZA⊆Fweakly dominates the approximation set ZB⊆F; that is, ZA⪯ZB, if every zB∈ZB is weakly dominated by at least one zA∈ZA.Definition 7 (ideal point).
Theideal pointz*=(z1*,…,zm*) represents the minimum possible value in each objective (typically, it cannot be achieved):
(6)zi*=min⁡{fi(x);x∈X}∀i∈{1,…,m}.Definition 8 (minimal element, minimal set).
LetZ⊆F be a set of objective vectors. An objective vector m∈Z is aminimal element of Z if
(7)∀z∈Z,z⪯m⟹z=m.
All minimal elements of Z constitute theminimal set of Z denoted in this paper by Zmin⁡.Maximal elements and the maximal setZmax⁡ can be defined dually.Definition 9 (Pareto front).
The set of Pareto optimal objective vectors known as thePareto front can be formally defined as the minimal set of f(X):
(8)Pf=f(X)min⁡.Definition 10 (Edgeworth-Pareto hull).
LetZ⊆F be an approximation set. The Edgeworth-Pareto hull (EPH) of Z includes all objective vectors weakly dominated by Z:
(9)ZEPH={z∈F;existsz′∈Zsothatz′⪯z}.


## 2.2. Empirical Attainment Function

Based on the multiobjective concept of goal-attainment [5] we say that an objective vector isattained 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 theattainment surface and its minimal elements are called theattainment anchors.Definition 11 (set boundary).
The boundary of a setZ⊆F is the intersection of the closure of Z, Z¯, with the closure of its complement, (F∖Z)¯:
(10)∂Z=Z¯∩(F∖Z)¯.Definition 12 (attainment surface, attainment anchors).
LetZ⊆F be an approximation set. Theattainment surface of Z is the boundary of the EPH of Z:
(11)SZ=∂ZEPH.
The minimal elements of SZ are calledattainment anchors and are equal to the original approximation set Z:
(12)AZ=SZmin⁡=Z.
Now, imagine that an optimization algorithm is runr times, producing r approximation sets. Then, each objective vector can be attained between 0 and r times.Definition 13 (empirical attainment function).
For algorithm A the empirical attainment function of the objective vectorz∈Z⊆F represents the frequency of attaining z by its approximation sets Z1A,…,ZrA:
(13)αrA(z)=α(Z1A,…,ZrA;z)=1r∑i=1r‍I(ZiA⪯{z}),
where I is the indicator function, defined as
(14)I(b)={1,ifbistrue,0,otherwise,
and ⪯ is the weak Pareto dominance relation between sets.This means that for algorithmA with approximation sets Z1A,…,ZrA an EAF value from the set of frequencies {0,1/r,2/r,…,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 k%-attainment surfaces (alsosummary attainment surfaces) in order to distinguish them from the “ordinary” attainment surfaces of single approximation sets. Vectors in the k%-attainment surface are the tightest objective vectors that have been attained in at least k% of the runs.Definition 14 (summary attainment surfaces).
LetZ1,…,Zr be r approximation sets of algorithm A and αrA(z) its EAF. Let the t/r-superlevel set of αrA be defined as
(15)Lt/r+(αrA)={z∈F;αrA(z)≥t/r}∀t∈{1,…,r}.
Then thesummary (or t/r-)attainment surfaceSt/r is equal to the boundary of Lt/r+(αrA):
(16)St/r=∂Lt/r+(αrA)∀t∈{1,…,r}.
We will use At/r to denote the attainment anchors of St/r:
(17)At/r=St/rmin⁡∀t∈{1,…,r}.
Note how the definition of an attainment surface coincides with the definition of a summary attainment surface ifr=1. In general, the number of all attainment anchors (including duplicates), nA=∑t=1r‍|At/r|, is much larger than the number of vectors contained in the approximation sets, nZ=∑t=1r‍|Zt| (in the context of the EAF computation [10], vectors from the approximation sets and attainment anchors were calledinput andoutput points, resp.). It has been shown in [10] that nA∈O(r2nZ) for three-objective optimization problems.We illustrate these concepts in the two-objective case. Figure2 shows an example of four approximation sets Z1,…,Z4 (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).Figure 2
Four approximation sets in a two-objective optimization problem, their summary attainment surfaces, and all attainment anchors. Colors denote areas with equal EAF values (darker colors correspond to larger values).If two algorithmsA and B are run r times each, they can be compared by computing and visualizing their EAF values αrA and αrB separately. However, there exists a straightforward way to directly compare the algorithms by computingEAF differences, that is, differences in their EAF values [11].Definition 15 (EAF differences).
Assume that algorithmsA and B are run r times each. The EAF differences between algorithms A and B are defined for each objective vector z∈Z as
(18)δrA-B(z)=αrA(z)-αrB(z).
Note that EAF differences need to be computed for all attainment anchors of the combined2r approximation sets. Defined in this way, the differences can adopt values from the set {-1,-(r-1)/r,…,0,1/r,…,1}. Positive EAF differences correspond to areas in the objective space where the algorithm A outperforms the algorithm B, while negative EAF differences denote areas in the objective space where the algorithm B outperforms the algorithm A. Naturally, where the differences are zero, both algorithms attain the area equally well.Let us focus on the example from Figure3, where algorithms A and B 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 A outperformed algorithm B (green hues) and areas where algorithm B outperformed algorithm A (red hues). We can readily notice that algorithm A is more successful in optimizing the first objective, while algorithm B 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.Figure 3
Four approximation sets in a two-objective optimization problem: two produced by algorithmA and two by algorithm B. Colored areas show where algorithm A outperformed algorithm B (green hues) and algorithm B outperformed algorithm A (red hues).We can safely assume that we are interested in a limited portion of the objective space, for example, the one enclosed by reference vectorsr1 and r2, where r1⪻r2. 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 t/r- and (t+1)/r-attainment surfaces, where t∈{1,…,r-1}, can be defined by two vertices: one from the t/r-attainment surface and the other from the (t+1)/r-attainment surface (Figure 4 shows 2D areas between 25%- and 50%-attainment surfaces). While the lower vertices correspond to the attainment anchors of the t/r-attainment surface, the upper vertices are not the attainment anchors of the (t+1)/r-attainment surface but theiropposite. The opposite can be defined for an arbitrary approximation set as follows.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.Definition 16 (approximation set opposite).
Let reference vectorsr1 and r2, where r1⪻r2, define the boundaries of the observed objective space
(19)R={z∈F;zi∈[ri1,ri2]∀i∈{1,…,m}}.
Let Z⊆R be an approximation set enclosed in this space with the attainment surface SZ. TheoppositeOZR of the approximation set Z within the observed objective space R is the maximal set of SZ∩R:
(20)OZR=(SZ∩R)max⁡.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 thehypervolume indicator [12], which is often used to measure the quality of approximation sets.Definition 17 (hypervolume indicator).
Let again reference vectorsr1 and r2, where r1⪻r2, define the boundaries of the observed objective space
(21)R={z∈F;zi∈[ri1,ri2]∀i∈{1,…,m}}.
Thehypervolume indicatorIH of an approximation set Z⊆R can be formulated via the empirical attainment function of Z as
(22)IH(Z)=∫z∈R‍α(Z;z)dz.


## 2.3. Visualization of EAFs

2D EAFs are most often visualized throughbest,median, andworst summary attainment surfaces corresponding to 1/r%, ~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 t/r-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 t, 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 mD case, where m>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 A and B, four plots can be produced, each corresponding to one of the values of αA, αB, δA-B, and δB-A. All attainment anchors are plotted using gray-shaded dots, where the shade corresponds to the value of the chosen function (αA, αB, δA-B, or δB-A). 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 theprosection 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 axis-aligned sampling surfaces of the EPH and use different colors to represent each slice. As the surface of ZEPH is exactly the attainment surface of Z, it should be rather straightforward to apply this method for visualizing attainment surfaces.


## 2.5. 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.


## 2.6. Direct Volume Rendering

Direct volume rendering (DVR) [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 namedirect).First, each voxel value is assigned color and opacity (the RGBA value) by means of a chosentransfer 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, shear-warp 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 at each sample position along the ray. The final RGB values shown in the image are computed from the accumulated RGBA values using the volume rendering integral [23].


## 3. 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.Figure 5
The linear and spherical approximation sets used in all visualization examples in Sections4 and 5.The first set consists of fivelinear approximation sets with a uniform distribution of vectors that satisfy the constraint
(23)∑i=13‍zi=cL.
In order to simulate the behavior of a stochastic algorithm, the values cL follow the normal distribution N(1,0.05). The second set contains fivespherical 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
(24)∑i=13‍zi2=cS2,
where cS~N(0.75,0.05). Each individual approximation set contains 100 objective vectors.The values ofcL and cS 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α5L and α5S to denote EAF values of the linear and spherical approximation sets, respectively. In addition, δ5L-S and δ5S-L will denote differences between those values.


## 4. 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.


### 4.1. Computing the Cuboids

The algorithm for computing EAFs [10] takes as input a series of approximation sets Z1,…,Zr and returns as output a series of sets containing attainment anchors A1,…,Ar that define the corresponding summary attainment surfaces and for which the following holds: Ai⪯Ai+1 for i∈{1,2,…,r-1}. From these sets of attainment anchors and a reference vector r2 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).Algorithm 1:Computing the cuboids.
Input: Sets of attainment anchors A1,…,Ar and a
reference vectorr2, for which A1⪯⋯⪯Ar
andz⪻r2 for all z∈Ai, i=1,…,r
Output: The set of cuboids C
A
r
+
1
←
{
r
2
};
r
1
←
z
*;
foreach adjacent pair of sets Ai and Ai+1, i∈{1,…,r}
do
O← opposite(Ai+1,r1,r2);
foreach  z∈Ai  do
O′←{o∈O;z⪻o};
foreach  o∈O′  do
C←C∪ cuboid(z, o, value(z));
end
end
endFirst, a set containing only the given reference vectorr2 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 r2 can also be computed. Then, the reference vector r1 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 Ai and Ai+1 and computes the opposite of Ai+1 using Algorithm 2. For each attainment anchor z∈Ai all dominated vectors from the opposite are stored in O′. Each pair (z,o), where o∈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.Algorithm 2:Computing the opposite of a 3D approximation set.
Input: Approximation set A and reference vectors r1 and
r2, for which r1⪯z⪻r2 for all z∈Ai,i=1,…,r
Output: The opposite O
O
←
{
r
2
};
foreach  z∈A  do
O′←{o∈O;z⪻o};
O←O-O′;
foreach  o∈O′  do
on1←(z1,o2,o3);
on2←(o1,z2,o3);
on3←(o1,o2,z3);
for  i←1  to  3  do
if ¬ redundant (oni, r1, O) then
O←O∪{oni};
end
end
end
endNote that there might be multiple cuboids stemming from the same attainment anchor (see the analog example of multiple rectangles in Figure4)—we call this aunion 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 Figure6 how to find the opposite of a 3D approximation set. Imagine a single cuboid defined by the reference vectors r1 and r2. The opposite is initialized to {r2}. Every time a vector from Z 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 z1 and z2 (see Figure 6(a)) and vector z3 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 o4 but not o5). Next, for every deleted vector o we create three new vectors in the following way:
(25)on1=(z1,o2,o3),on2=(o1,z2,o3),on3=(o1,o2,z3).Two steps in the construction of the opposite of the approximation setZ.
(a)
Opposite of{z1,z2}
(b)
Opposite of{z1,z2,z3}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 withr1 or collinear with any of the vectors from the existing opposite). In our example, we add to the opposite vectors o6 and o7 but not o8=(2,6.5,5) (not pictured in the figure) as it is collinear with o5 and thus redundant. When these steps have been taken for every vector from Z, the resulting set O represents the opposite of Z. See Algorithm 2 for the algorithmic notation of the described procedure.As we cannot visualize the entire (infinite) objective space, Algorithms1 and 2 use two reference vectors r1 and r2 to limit it. In order to preserve all information, they need to be set in such a way that r1⪯z⪻r2 for all attainment anchors z. While the ideal point is the most sensible choice for r1, we can choose any objective vector dominated by all attainment anchors for r2.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 ofn attainment anchors has the worst-case computational complexity of O(n3), which means that all cuboids between r pairs of attainment surfaces can be computed in O(rn3) time. However, in practice, the loopforeacho∈O′do in Algorithm 2 is executed only a few times for each z∈A, 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 O(nlog⁡n). If this would be possible, the total cost of computing the cuboids would be O(rnlog⁡n). 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.


### 4.2. Visualizing Cuboids Using Slicing

In this paper, the slicing plane is aligned with one of the three axes (assume for now that this isf3), includes the origin r1, and cuts through the underlying plane (f1f2) 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 Figure7). 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
(26)φz=arctan(z2-r21z1-r11),φo=arctan(o2-r21o1-r11),
this is done in the following way:
(27)z′={(r11+z2-r21tanφ,z2,z3),ifφz≥φ,(z1,r21+(z1-r11)tanφ,z3),ifφz<φ,o′={(o1,r21+(o1-r11)tanφ,o3),ifφo≥φ,(r11+o2-r21tanφ,o2,o3),ifφo<φ,z′′={(r11+z2-r21sinφ,z3),ifφz≥φ,(r11+z1-r11cos⁡φ,z3),ifφz<φ,o′′={(r11+o1-r11cos⁡φ,o3),ifφo≥φ,(r11+o2-r21sinφ,o3),ifφo<φ.
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 ⪰).Figure 7
A cuboid sliced using the plane aligned withf3 that slices the f1f2 plane under angle. The cuboid vertices z and o are projected to 3D rectangle vertexes z′and o′, respectively.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 α5L and α5S are shown separately, while the EAF differences δ5L-S and δ5S-L are presented on the same plot (see 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.
(a)
Slice ofα5L at φ = 5°
(b)
Slice ofα5S at φ = 5°
(c)
Slice ofδ5L-S and δ5S-L at φ = 5°
(d)
Slice ofα5L at φ = 45°
(e)
Slice ofα5S at φ = 45°
(f)
Slice ofδ5L-S and δ5S-L at φ = 45°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 axisf1 or f2, too. However, as our benchmark approximation sets are rather symmetrical, there is no need to do this here.


### 4.3. Visualizing Cuboids Using MIP

As explained in Section2.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 Figure9. The differences δ5L-S and δ5S-L 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.MIP of exact 3D EAF differences between the benchmark approximation sets.
(a)
MIP ofδ5L-S
(b)
MIP ofδ5S-L


## 4.1. Computing the Cuboids

The algorithm for computing EAFs [10] takes as input a series of approximation sets Z1,…,Zr and returns as output a series of sets containing attainment anchors A1,…,Ar that define the corresponding summary attainment surfaces and for which the following holds: Ai⪯Ai+1 for i∈{1,2,…,r-1}. From these sets of attainment anchors and a reference vector r2 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).Algorithm 1:Computing the cuboids.
Input: Sets of attainment anchors A1,…,Ar and a
reference vectorr2, for which A1⪯⋯⪯Ar
andz⪻r2 for all z∈Ai, i=1,…,r
Output: The set of cuboids C
A
r
+
1
←
{
r
2
};
r
1
←
z
*;
foreach adjacent pair of sets Ai and Ai+1, i∈{1,…,r}
do
O← opposite(Ai+1,r1,r2);
foreach  z∈Ai  do
O′←{o∈O;z⪻o};
foreach  o∈O′  do
C←C∪ cuboid(z, o, value(z));
end
end
endFirst, a set containing only the given reference vectorr2 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 r2 can also be computed. Then, the reference vector r1 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 Ai and Ai+1 and computes the opposite of Ai+1 using Algorithm 2. For each attainment anchor z∈Ai all dominated vectors from the opposite are stored in O′. Each pair (z,o), where o∈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.Algorithm 2:Computing the opposite of a 3D approximation set.
Input: Approximation set A and reference vectors r1 and
r2, for which r1⪯z⪻r2 for all z∈Ai,i=1,…,r
Output: The opposite O
O
←
{
r
2
};
foreach  z∈A  do
O′←{o∈O;z⪻o};
O←O-O′;
foreach  o∈O′  do
on1←(z1,o2,o3);
on2←(o1,z2,o3);
on3←(o1,o2,z3);
for  i←1  to  3  do
if ¬ redundant (oni, r1, O) then
O←O∪{oni};
end
end
end
endNote that there might be multiple cuboids stemming from the same attainment anchor (see the analog example of multiple rectangles in Figure4)—we call this aunion 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 Figure6 how to find the opposite of a 3D approximation set. Imagine a single cuboid defined by the reference vectors r1 and r2. The opposite is initialized to {r2}. Every time a vector from Z 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 z1 and z2 (see Figure 6(a)) and vector z3 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 o4 but not o5). Next, for every deleted vector o we create three new vectors in the following way:
(25)on1=(z1,o2,o3),on2=(o1,z2,o3),on3=(o1,o2,z3).Two steps in the construction of the opposite of the approximation setZ.
(a)
Opposite of{z1,z2}
(b)
Opposite of{z1,z2,z3}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 withr1 or collinear with any of the vectors from the existing opposite). In our example, we add to the opposite vectors o6 and o7 but not o8=(2,6.5,5) (not pictured in the figure) as it is collinear with o5 and thus redundant. When these steps have been taken for every vector from Z, the resulting set O represents the opposite of Z. See Algorithm 2 for the algorithmic notation of the described procedure.As we cannot visualize the entire (infinite) objective space, Algorithms1 and 2 use two reference vectors r1 and r2 to limit it. In order to preserve all information, they need to be set in such a way that r1⪯z⪻r2 for all attainment anchors z. While the ideal point is the most sensible choice for r1, we can choose any objective vector dominated by all attainment anchors for r2.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 ofn attainment anchors has the worst-case computational complexity of O(n3), which means that all cuboids between r pairs of attainment surfaces can be computed in O(rn3) time. However, in practice, the loopforeacho∈O′do in Algorithm 2 is executed only a few times for each z∈A, 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 O(nlog⁡n). If this would be possible, the total cost of computing the cuboids would be O(rnlog⁡n). 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.


## 4.2. Visualizing Cuboids Using Slicing

In this paper, the slicing plane is aligned with one of the three axes (assume for now that this isf3), includes the origin r1, and cuts through the underlying plane (f1f2) 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 Figure7). 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
(26)φz=arctan(z2-r21z1-r11),φo=arctan(o2-r21o1-r11),
this is done in the following way:
(27)z′={(r11+z2-r21tanφ,z2,z3),ifφz≥φ,(z1,r21+(z1-r11)tanφ,z3),ifφz<φ,o′={(o1,r21+(o1-r11)tanφ,o3),ifφo≥φ,(r11+o2-r21tanφ,o2,o3),ifφo<φ,z′′={(r11+z2-r21sinφ,z3),ifφz≥φ,(r11+z1-r11cos⁡φ,z3),ifφz<φ,o′′={(r11+o1-r11cos⁡φ,o3),ifφo≥φ,(r11+o2-r21sinφ,o3),ifφo<φ.
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 ⪰).Figure 7
A cuboid sliced using the plane aligned withf3 that slices the f1f2 plane under angle. The cuboid vertices z and o are projected to 3D rectangle vertexes z′and o′, respectively.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 α5L and α5S are shown separately, while the EAF differences δ5L-S and δ5S-L are presented on the same plot (see 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.
(a)
Slice ofα5L at φ = 5°
(b)
Slice ofα5S at φ = 5°
(c)
Slice ofδ5L-S and δ5S-L at φ = 5°
(d)
Slice ofα5L at φ = 45°
(e)
Slice ofα5S at φ = 45°
(f)
Slice ofδ5L-S and δ5S-L at φ = 45°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 axisf1 or f2, too. However, as our benchmark approximation sets are rather symmetrical, there is no need to do this here.


## 4.3. Visualizing Cuboids Using MIP

As explained in Section2.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 Figure9. The differences δ5L-S and δ5S-L 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.MIP of exact 3D EAF differences between the benchmark approximation sets.
(a)
MIP ofδ5L-S
(b)
MIP ofδ5S-L


## 5. 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.


### 5.1. Discretization into Voxels

A voxel is a single point on a regularly spaced 3D grid ofv1×v2×v3 voxels. Recall that, based on the reference vectors r1 and r2, the 3D observed objective space R is defined as
(28)R={z∈F;zi∈[ri1,ri2]∀i∈{1,2,3}}.
If R is a cube, v1=v2=v3; otherwise, some care must be taken to ensure that the grid of voxels is truly regular. Either R must be normalized prior to approximation or the number of voxels in each dimension vi must be set to be proportional to ri2-ri1 for all i∈{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.From the observed objective spaceR the grid of v1×v2×v3 voxels is constructed so that the voxel at grid position (k1,k2,k3) has the following coordinates:
(29)((r12-r11)(2k1-1)2v1,(r22-r21)(2k2-1)2v2,(r32-r31)(2k3-1)2v3),
where ki∈{1,…,vi} for i∈{1,2,3}.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.


### 5.2. 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 Section4.2). In order to avoid showing plots that are very similar to those presented in Figure 8, Figure 10 shows only the slices of δ5L-S and δ5S-L at angle φ=5∘ produced using different discretizations.Slices of exact and approximatedδ5L-S and δ5S-L at angle φ = 5° showing the approximation error.
(a)
Using cuboids
(b)
Using 643 voxels
(c)
Using 1283 voxels
(d)
Using 2563 voxelsWe 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 1283 voxels suffices, we will use this discretization in the remainder of the paper.


### 5.3. 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.MIP of approximated 3D EAF differences between the benchmark approximation sets.
(a)
MIP ofδ5L-S
(b)
MIP ofδ5S-LFigure11 shows the MIP for δ5L-S and δ5S-L and we can see that although these are approximations of the images presented in Figure 9, the results are quite similar.


### 5.4. 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.Figures12 and 13 show visualization using DVR of the EAF differences between the two benchmark approximation sets δ5L-S and δ5S-L, 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.DVR of approximated 3D EAF differences between the linear and spherical benchmark approximation setsδ5L-S.
(a)
DVR ofδ5L-S = 1/5
(b)
DVR ofδ5L-S = 2/5
(c)
DVR ofδ5L-S = 3/5
(d)
DVR ofδ5L-S = 4/5
(e)
DVR ofδ5L-S = 5/5
(f)
DVR ofδ5L-S∈{1/5,5/5}DVR of approximated 3D EAF differences between the spherical and linear benchmark approximation setsδ5S-L.
(a)
DVR ofδ5S-L = 1/5
(b)
DVR ofδ5S-L = 2/5
(c)
DVR ofδ5S-L = 3/5
(d)
DVR ofδ5S-L = 4/5
(e)
DVR ofδ5S-L = 5/5
(f)
DVR ofδ5S-L∈{1/5,5/5}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 Figures12(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. Figure14 shows an example of such a visualization for the spherical approximation sets, where only the values δ5S∈{1/5,4/5} are shown (opacity is set to 0.3 for these two values and to 0 for all other values).Figure 14
DVR of approximated 3D EAF values of the spherical benchmark approximation sets forα5S∈{1/5,4/5}.


## 5.1. Discretization into Voxels

A voxel is a single point on a regularly spaced 3D grid ofv1×v2×v3 voxels. Recall that, based on the reference vectors r1 and r2, the 3D observed objective space R is defined as
(28)R={z∈F;zi∈[ri1,ri2]∀i∈{1,2,3}}.
If R is a cube, v1=v2=v3; otherwise, some care must be taken to ensure that the grid of voxels is truly regular. Either R must be normalized prior to approximation or the number of voxels in each dimension vi must be set to be proportional to ri2-ri1 for all i∈{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.From the observed objective spaceR the grid of v1×v2×v3 voxels is constructed so that the voxel at grid position (k1,k2,k3) has the following coordinates:
(29)((r12-r11)(2k1-1)2v1,(r22-r21)(2k2-1)2v2,(r32-r31)(2k3-1)2v3),
where ki∈{1,…,vi} for i∈{1,2,3}.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.


## 5.2. 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 Section4.2). In order to avoid showing plots that are very similar to those presented in Figure 8, Figure 10 shows only the slices of δ5L-S and δ5S-L at angle φ=5∘ produced using different discretizations.Slices of exact and approximatedδ5L-S and δ5S-L at angle φ = 5° showing the approximation error.
(a)
Using cuboids
(b)
Using 643 voxels
(c)
Using 1283 voxels
(d)
Using 2563 voxelsWe 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 1283 voxels suffices, we will use this discretization in the remainder of the paper.


## 5.3. 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.MIP of approximated 3D EAF differences between the benchmark approximation sets.
(a)
MIP ofδ5L-S
(b)
MIP ofδ5S-LFigure11 shows the MIP for δ5L-S and δ5S-L and we can see that although these are approximations of the images presented in Figure 9, the results are quite similar.


## 5.4. 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.Figures12 and 13 show visualization using DVR of the EAF differences between the two benchmark approximation sets δ5L-S and δ5S-L, 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.DVR of approximated 3D EAF differences between the linear and spherical benchmark approximation setsδ5L-S.
(a)
DVR ofδ5L-S = 1/5
(b)
DVR ofδ5L-S = 2/5
(c)
DVR ofδ5L-S = 3/5
(d)
DVR ofδ5L-S = 4/5
(e)
DVR ofδ5L-S = 5/5
(f)
DVR ofδ5L-S∈{1/5,5/5}DVR of approximated 3D EAF differences between the spherical and linear benchmark approximation setsδ5S-L.
(a)
DVR ofδ5S-L = 1/5
(b)
DVR ofδ5S-L = 2/5
(c)
DVR ofδ5S-L = 3/5
(d)
DVR ofδ5S-L = 4/5
(e)
DVR ofδ5S-L = 5/5
(f)
DVR ofδ5S-L∈{1/5,5/5}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 Figures12(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. Figure14 shows an example of such a visualization for the spherical approximation sets, where only the values δ5S∈{1/5,4/5} are shown (opacity is set to 0.3 for these two values and to 0 for all other values).Figure 14
DVR of approximated 3D EAF values of the spherical benchmark approximation sets forα5S∈{1/5,4/5}.


## 6. Discussion

Let us summarize the properties of all presented visualization methods (see also Table1). 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 R 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.Table 1
Summarized properties of the presented visualization methods.
Slicing
MIP
DVR
Exact EAF values
+ Enables detailed analysis of the approximation sets− Visualizes only one slice at a time
− Not sensible to use since usually a large portion of the objective space has the maximum EAF value
− Not possible to use
Exact EAF differences
+ Enables detailed analysis of the approximation sets+ Capable of simultaneously visualizing positive and negative differences without overlapping − Visualizes only one slice at a time
+ No parameters to set+ Shows all values in a single image − No sense of depth − Feasible only for a limited number of cuboids (impractical)
− Not possible to use
Approximated EAF values
− No need to use it on the approximation as it works well for exact EAFs
− Not sensible to use since usually a large portion of the objective space has the maximum EAF value
+ Nice visualizations+ Enables looking through the volume and preserves the sense of depth − Requires definition of the transfer function
Approximated EAF differences
− No need to use it on the approximation as it works well for exact EAFs
+ No parameters to set+ Shows all values in a single image − No sense of depth
+ Nice visualizations+ Enables looking through the volume and preserves the sense of depth − Requires definition of the transfer functionIt 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.


## 7. 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].


### 7.1. 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 Table2)) 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 desiredmetallurgical length (the length of the liquid core in the strand), distance from desiredshell thickness (thickness of the solid shell at mold exit), and distance from desiredsurface temperature at a predefined point in the strand. Table 3 details the variables defining these optimization objectives. Their bounds and desired values were determined by domain experts.Table 2
Steel casting process parameters.
Parameter
Lower bound
Upper  bound
Disc. step
Casting speed (m/min)
1.5
2.0
0.01
Mold outlet coolant temp. (°C)
33
35
1
Wreath coolant flow (m3/h)
10
40
5
Spray coolant flow (m3/h)
25
65
5Table 3
Variables defining the optimization objectives.
Variable
Lower  bound
Upper  bound
Desired  value
Metallurgical length (m)
10
11
10
Shell thickness (mm)
11
15
13
Surface temperature (°C)
1115
1130
1122.5If 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
(30)fk=|yk-yk*|∀i∈1,2,3,
where yk is the value of the observed variable and yk* 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 andcontinuous. 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.


### 7.2. 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 Figure15. 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.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).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 Table3, 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 Sections4 and 5 can be used on these results. Figure 16 presents the DVR of approximated EAF differences between the two algorithms, δ5DEMO-ES and δ5ES-DEMO, 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 summary attainment surfaces. The narrow layer between α5DEMO=1/5 and α5DEMO=4/5 as well as α5DEMO=1/5 and α5DEMO=5/5 (see Figure 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).DVR of approximated 3D EAF differences between the two algorithms for all values in{1/5,…,5/5}.
(a)
DVR ofδ5DEMO-ES
(b)
DVR ofδ5ES-DEMODVR of approximated 3D EAF values of DEM.
(a)
DVR ofα5DEMO∈{1/5,4/5}
(b)
DVR ofα5DEMO∈{1/5,5/5}Slices of exact 3D EAF values of DEMO at two angles.
(a)
Slice ofα5DEMO at φ = 25°
(b)
Slice ofα5DEMO at φ = 45°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δ5ES-DEMO=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 δ5ES-DEMO=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.Figure 19
MIP of exact 3D EAF differences between ES and DEMO,δ5ES-DEMO.


## 7.1. 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 Table2)) 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 desiredmetallurgical length (the length of the liquid core in the strand), distance from desiredshell thickness (thickness of the solid shell at mold exit), and distance from desiredsurface temperature at a predefined point in the strand. Table 3 details the variables defining these optimization objectives. Their bounds and desired values were determined by domain experts.Table 2
Steel casting process parameters.
Parameter
Lower bound
Upper  bound
Disc. step
Casting speed (m/min)
1.5
2.0
0.01
Mold outlet coolant temp. (°C)
33
35
1
Wreath coolant flow (m3/h)
10
40
5
Spray coolant flow (m3/h)
25
65
5Table 3
Variables defining the optimization objectives.
Variable
Lower  bound
Upper  bound
Desired  value
Metallurgical length (m)
10
11
10
Shell thickness (mm)
11
15
13
Surface temperature (°C)
1115
1130
1122.5If 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
(30)fk=|yk-yk*|∀i∈1,2,3,
where yk is the value of the observed variable and yk* 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 andcontinuous. 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.


## 7.2. 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 Figure15. 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.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).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 Table3, 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 Sections4 and 5 can be used on these results. Figure 16 presents the DVR of approximated EAF differences between the two algorithms, δ5DEMO-ES and δ5ES-DEMO, 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 summary attainment surfaces. The narrow layer between α5DEMO=1/5 and α5DEMO=4/5 as well as α5DEMO=1/5 and α5DEMO=5/5 (see Figure 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).DVR of approximated 3D EAF differences between the two algorithms for all values in{1/5,…,5/5}.
(a)
DVR ofδ5DEMO-ES
(b)
DVR ofδ5ES-DEMODVR of approximated 3D EAF values of DEM.
(a)
DVR ofα5DEMO∈{1/5,4/5}
(b)
DVR ofα5DEMO∈{1/5,5/5}Slices of exact 3D EAF values of DEMO at two angles.
(a)
Slice ofα5DEMO at φ = 25°
(b)
Slice ofα5DEMO at φ = 45°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δ5ES-DEMO=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 δ5ES-DEMO=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.Figure 19
MIP of exact 3D EAF differences between ES and DEMO,δ5ES-DEMO.


## 8. 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.



---

*Source: 569346-2014-09-17.xml*