Solution architecture and selected algorithms
4.2 Selected algorithms
For the purpose of this project, two algorithms were selected for the suite. One for the search mechanism and the other for rigorous statistical R&S analysis.
The selected search algorithm is the cross-entropy method for multi-objective op-timisation (MOO CEM) (Bekker, 2012) while the statistical, ranking and selection algorithm is the MMY procedure (Yoon,2018). MOO CEM is believed to be a good search mechanism for large-scale MOSO problems which are often computationally demanding. MOO CEM is, in effect, known to be a relatively fast converging meta-heuristic (Bekker, 2012). MMY, on the other hand, is believed to be an effective and an efficient modern R&S procedure, ideal for MOSO problems.
4.2.1 The MOO CEM metaheuristic
The MOO CEM is a multi-objective adaptation of the CEM metaheuristic described in Section2.5.1.3. The algorithm was developed by Bekker(2012). The author dedicates this section to the full description of the MOO CEM as it is the principal algorithm of the MOOSolver library. The description is based on the works ofBekker & Aldrich (2011).
The algorithm adapts the CEM for MOO problems by using a number of key con-cepts. The first one is that of a Working matrix. The MOO CEM metaheuristic uses a Working matrix that consists of N rows and n+m+1 columns, where N is an arbitrary number of solutions (i.e. the population size), n is the number of decision variables (DVs) and m is the number of objectives. Sample values of the first DV are stored in column 1, the second DV in column 2, and so on up to column n. The objective function values for objective 1 are stored in column n + 1, for objective 2 in column
n + 2, and for objective m in column n + m. The last column is used to store the Pareto ranking value of each solution. The structure is shown in Table4.1.
Table 4.1: Structure of the working matrix.
Decision variables Objectives Rank
x11 x12 · · · x1n f11 f12 · · · f1m ρ1
... ... ... ... ... ... ... xN 1 xN 2 · · · xN n fN 1 fN 2 · · · fN m ρN
The metaheuristic generates new populations by creating a sequence of probability density functions (pdfs) for each DV in every iteration (this is the “specified” mechanism alluded to in phase one of the CEM in Section2.5.1.3). To form a sample vector xifrom a pdf hi(·; ˆVt−1), a truncated normal distribution is used for each decision variable. For the n DVs defined over ranges [li, ui], li is the lower limit and ui the upper limit of DV
Using truncated distributions makes it easy to contain the search. As required by the CEM, an arbitrarily large value for σi (i.e. the standard deviation) is initially assigned, using σi = 10 · (ui− li). The first n columns of the working matrix are filled with sample values from each applicable truncated normal distribution.
Next, each of the objective functions is evaluated using the row vectors X1i, ..., XN i
of Table4.1. This yields two or more performance vectors fj(x) with 1 < j ≤ m.
The best combinations of objective functions are found by doing a Pareto-ranking using Algorithm1 in Section2.1.
The values of the decision variables in the non-dominated set (i.e. the elite matrix ) provided by Algorithm 1 are used to construct a histogram for each decision variable.
The histogram concept is the second key concept used to adapt the CEM into a MOO algorithm. The histograms provide guiding information for the MOO CEM algorithm and are maintained while the algorithm is searching for non-dominated solutions.
The concept works as follows: for a decision variable that is defined in the range [li, ui], the lower boundary of the first class is set equal to li, and the upper boundary of the last class is set equal to ui. Next, the upper boundary of the first class is set equal to the minimum value of the decision variable xi in the elite matrix, i.e. min(Elite(·, i)).
The lower boundary of the last class is equal to the maximum value of the decision variable in the elite matrix, namely max(Elite(·, i)). A number of equal-sized classes are formed between these two boundaries using (max(Elite(·, i)) − min(Elite(·, i)))/r if r of these classes are formed, resulting in a total number of r + 2 classes (see Figure 4.2).
Figure 4.2: Example of a histogram for decision variable xi and r = 3.
The class boundaries for the histogram of decision variable xi are recorded in a vector Ci = {ci1, ci2, ..., ci(r+2), ci((r+2)+1)}, with ci1 = li and ci((r+2)+1) = ui. Note that Ci contains r + 3 elements because the histogram has r + 2 classes, and that the class widths of the first class ([ci1, ci2]) and the last class ([ci(r+2), ci((r+2)+1)]) can be different from each other and from the widths of the r classes.
The elite matrix has the same columns as the working matrix shown in Table 4.1, and the values in column i, 1 ≤ i ≤ n are used to determine frequency values for decision variable xi. The decision variable values are classified according to the following rule:
xij belongs to the class [cik, ci(k+1)) if cik ≤ xij < ci(k+1), 1 ≤ k ≤ r + 2. The histogram frequency values are recorded in a vector Ri = {τi1, τi2, ..., τi(r+1), τi(r+2)}, where τi1 is
equal to the frequency count of decision variable xi in the range [ci1, ci2), τi2 represents the count in range [ci2, ci3), and so on.
In preparation for the next iteration of the algorithm, the new population of possible solutions is formed proportionally according to the class frequencies for each decision variable: suppose the elite matrix contains Er rows and there are τik occurrences in class [cij, ci(j+1)) for a given decision variable xi, then bN τik/Erc values will be created from this class range for this variable (the population size is N and 1 ≤ k ≤ r + 2).
When the proportional numbers do not add up to N due to the rounding down of the proportion calculation, the small difference is arbitrarily added to the last class.
When generating observations from a class range [cij, ci(j+1)], temporary µ0ik and σik0 values are used. These values are associated with the specific histogram class ranges, so for the class [cik, ci(k+1)) for decision variable xi, the parameter estimators are µ0ik = cik+ U (ci(k+1)− cik), whereas σik0 = (ci(k+1)− cik), 1 ≤ k ≤ r + 2 and U is a uniformly distributed random number.
To prevent premature convergence, the histogram frequencies are adjusted during each iteration t with a preset probability of typically 0.1 − 0.3. To do so, the maximum frequency over all classes is determined for a given decision variable. The frequency in each class is then subtracted from this value, resulting in an inverted histogram as shown in Figure 4.3. This ensures that search ranges that were given small proportions of population candidate allocations receive higher proportions of allocations, while search ranges with high proportions of population allocations receive fewer allocations after frequency inversion.
The algorithm will readjust the frequencies according to the rankings returned by the candidates so that a class that does not contribute to the elite matrix effectively becomes eliminated as the search progresses. The histogram concept also allows for the accommodation of discontinuous search spaces.
To ensure exploitation in the MOO context, the process described above is repeated ol times as an outer loop of the algorithm. This is the third concept used to adapt the CEM for MOO. After each loop, the elite matrix is ranked again and the number of classes of the histograms is incremented. Increasing the number of classes as the search progresses makes it possible to maintain good combinations of decision variable values as the resolution of the decision variable spaces becomes finer.
0 10 20 30 40 50
ci1 ci2 ci3 ci4 ci5 ci6
Classes
Frequency
Figure 4.3: The inverted histogram of Figure4.2.
The parameter vectors (µi, σi) are updated using (4.1) and (4.2) respectively, and the values in the DV columns of the elite matrix.
ˆ
µi,t = αmµ˜i,t+ (1 − αm)ˆµi,t−1 (4.1) ˆ
σi,t = αmσ˜i,t+ (1 − αm)ˆσi,t−1 (4.2) where t is the iteration count index and αm = 0.7 in all cases. This process is continued until the σi-value of each decision variable has decreased below a common threshold epsilon (). On algorithm termination, the elite matrix should contain the solutions members of the Pareto set, as well as their associated objective functions values.
To support exploration and exploitation of the search, the initial ranking threshold th is relaxed and a value of two is selected. When a new loop starts and a new popula-tion is formed, the elite matrix is trimmed and the threshold is set to one. When the algorithm terminates, the existing elite matrix is refined a final time, and the threshold then used is zero, which means all solutions selected are non-dominated. This con-stitutes the fourth and final concept. The algorithm is presented in pseudo-code as Algorithm9.
Algorithm 9 MOO CEM metaheuristic
1: Set Elite = ∅, t = 1, ol= 1.
2: Initialise decision variable vectors xi = ∅, 1 ≤ i ≤ n.
3: For each decision variable xi, 1 ≤ i ≤ n, initialise a histogram class vec-tor Ci = {ci1, ci2, ..., ci(r+2), ci((r+2)+1)} and histogram frequency vector Ri =
15: Rank the performance values using the Pareto ranking of Algorithm1with a relaxed th = 2 to obtain an updated elite matrix Elite.
16: Form new histogram class vectors Ci and histogram frequency vectors Ri based on Elite, 1 ≤ i ≤ n.
17: Use the values in Elite and compute ˜vit for all 1 ≤ i ≤ n.
18: Smooth the vectors ˜vit for all 1 ≤ i ≤ n, using (4.1) and (4.2).
19: If all σi,t > or less than the allowable number of evaluations has been done, increment t and reiterate from Step 4.
20: Rank the elite matrix using the Pareto ranking of Algorithm1 with th = 1.
21: Increment ol.
22: If ol is less than the allowable number of loops, return to Step 2.
23: Rank the elite matrix using the Pareto ranking of Algorithm1with th= 0 to obtain the final elite matrix.
4.2.2 The MMY procedure
ProcedureMMY is the multi-objective variant of the MY procedure presented in Section 2.4.1.1. It is the first MOO R&S procedure with the indifference-zone approach in the literature, that guarantees correct selection following the Bayesian probabilistic theory (Section 2.4.1.1). The algorithm was selected to equip MOOSolver with a modern, efficient multi-objective R&S procedure.
Table 4.2 provides the notation used in procedure MMY, which is presented in Algorithm10.
Table 4.2: Notation for procedureMMY.
M the number of systems in the problem;
S the feasible solution set, i.e., S = {1, ..., M };
I the set of systems that are still in competition;
H the number of objectives;
K the objective set, i.e., K = {1, ..., H};
Ni the total number of simulation replications assigned to system i;
X¯ik(Ni) the sample mean of system i for objective k based on Ni observations;
Sp the observed Pareto set based on ¯Xik (i ∈ S and k ∈ K);
Spc the observed non-Pareto set based on ¯Xik (i ∈ S and k ∈ K);
n0 the number of simulation replications at the first stage;
δk∗ the indifference-zone value for objective k;
P∗ the minimum required value for P(CS).
Following are some definitions used in the procedure. Let
δijk = max{δ∗k, ¯Xjk(Nj) − ¯Xik(Ni)} (4.3) and dxe denotes the smallest integer greater than x. Consider a pair of systems (i, j) where system i is observed as non-dominated and system j can be any other system in S. This pair (i, j) (i ∈ Sp and j ∈ S, j 6= i) is relevant to Steps 4 and 5 in Algorithm 10.
For each pair (i, j) (i ∈ Sp and j ∈ S, j 6= i), let
K1 = {k | | ¯Xjk− ¯Xik| ≤ δ∗k, k ∈ K} (4.7)
Algorithm 10 ProcedureMMY
1: Select the probability of correction requirement P∗ = 1 − α, the indifference-zone value δ∗k for each objective k ∈ K, and the first-stage sample size n0 ≥ 10. Set I = {1, 2, ..., M } and β = Mα.
2: Simulate n0 replications for all M systems, and calculate sample means ¯Xik(n0) and sample variances Sik2(n0) (i ∈ S and k ∈ K). Let Ni= n0.
3: Observe the Pareto set Sp and the non-Pareto set Spc based on the sample means X¯ik(Ni) (i ∈ S and k ∈ K) using Algorithm 1 without the indifference-zone where h3 is the solution to (4.12).
8: Delete system j from I if conditions in (4.6) are satisfied.
9: If |I| = 0, then stop and present the current Pareto set Sp as the final solution set.
Otherwise, for each system i ∈ Sp∩ I, that is, systems in Sp that were not deleted from I in Step 6, add system j ∈ S (j 6= i) to I if it does not satisfy conditions (4.4) or (4.5). Similarly, for each system j ∈ Spc∩ I, that is, systems in Spcthat were not deleted from I in Step 8, add the corresponding system i ∈ Sp to I if it does not satisfy (4.6). Go to Step 10.
10: Take one additional observation Xi,k,Ni+1 from each system i ∈ I, and set Ni ← Ni+ 1 (∀i ∈ I). Set I = {1, 2, ..., M } and update ¯Xik(Ni) and Sik2ik(N i) for all
and
where Φ denotes the cumulative distribution function (cdf) of the standard normal dis-tribution. Note that K1 and k0 should be defined for every pair of (i, j) (i ∈ Sp and j ∈ S, j 6= i). Step 4 in Algorithm10 deals with (i, j) pairs when K1 = K, that is, system i and j are observed to be indifferent to each other, while Step 5 considers the case when K1 6= K.
The constants h1 (in Step 4) and h2 (in Step 5) are the solution to the following equations, respectively: Nj− 1 degrees of freedom, respectively. The system that dominates system j with the maximum probability can be found with
i = arg max
Note that such i should be defined for every j ∈ Spc. This pair of systems (i, j) (i ∈ Sp, j ∈ Scp) is considered in Step 7 in Algorithm10. The constant h3 in the same step of the algorithm is the solution to
Nj− 1 degrees of freedom, respectively.
4.2.2.1 The relaxed Pareto set approach
An important particularity of theMMY procedure is that it outputs an approximate relaxed Pareto set with respect to the indifference-zone values selected by the user. In effect, one would expect an IZ-based MOO R&S procedure to include the IZ concept during its Pareto ranking step; however, this is not the case in procedure MMY as pointed out in Step 3 of Algorithm 10. The relaxed Pareto set approach is, in effect, one that is not strict about the IZ concept with regards to solutions that should be considered as non-dominated. Yoon (2018) found in her work that using a strict IZ regime can cause the R&S algorithm to run indefinitely when comparing a pair of systems whose true means are very close. Hence, the researcher (i.e. Yoon (2018)) proposed going instead for a relaxed Pareto set in order to avoid this. This approach therefore makes the procedure, one that is safer and more flexible in practice.
Moreover, the researcher (i.e. Yoon (2018)) also argues in her work that allowing the final Pareto set to possibly consider as Pareto-optimal solutions, non-dominated solutions (under no IZ regime) that would otherwise be seen as dominated under a strict IZ regime, gives the decision-maker a more comprehensive final set.
To illustrate the concept, consider Figure 4.4 where the solutions in red are non-dominated solutions and those in black are non-dominated ones. The IZ value for both objectives in this figure is selected as 0.5. A relaxed Pareto set (Figure4.4(c)) contains all non-dominated systems that do not have indifferent systems (solutions in red), at least one system from a group of indifferent systems (solutions in green); and regarding systems classified as non-dominated without the IZ concept, but dominated under the IZ regime (solutions in blue), a relaxed Pareto set may or may not contain them.
Obj. 2
Figure 4.4: Pareto set examples (Yoon,2018).
4.2.2.2 MMY implementation challenge
The major challenge in implementing the MMY procedure is the hi values calcula-tion (i = 1, 2, 3). hi values are solutions to double integral equations that are time-consuming to solve as they are solved numerically. It is therefore impractical to dy-namically calculate these values when the procedure is being executed. Instead, these values are calculated beforehand and saved in tables (i.e. an h1table contains h1values etc.) which the procedure can then look-up during execution. Generating these tables, however, is another challenge. It can be observed in (4.9), (4.10) and (4.12), that they all account for the total number of systems M being considered by the procedure. This means that every hi table is unique to a specific number of total systems being consid-ered by the procedure. In other words, for a problem where e.g. M = 8, the procedure will look-up a different set of h tables than one where e.g. M = 9.
Now, consider using the results of a metaheuristic that uses the Pareto approach directly as inputs to the R&S procedure. This would mean that the R&S procedure must have a set of h tables for every possible number of solutions in the approximated Pareto set obtained by the metaheuristic. This is in effect a difficult thing to do, in addition to being impractical as the set of solutions could be relatively very large.
This is what led the author to consider the interactive approach and actually limit the number of systems that the user can select from the approximate Pareto set. In this study, this number was limited to up to 10, meaning that nine set of h tables had to be calculated beforehand (i.e. 1 < M ≤ 10). Note that the author kept P∗ constant with a value of 90.
4.3 Chapter summary
In this chapter, the solution approach proposed for the purpose of this study was described. Additionally, the selected algorithms for the MOO optimisation suite were also presented and described in detail.
In the next chapter, the development and implementation of the MOO optimisation suite are presented.