Reference point based evolutionary multiobjective optimization with dynamic resampling for production systems improvement
Abstract
Purpose
Stochastic simulation is a popular tool among practitioners and researchers alike for quantitative analysis of systems. Recent advancement in research on formulating production systems improvement problems into multiobjective optimizations has provided the possibility to predict the optimal tradeoffs between improvement costs and system performance, before making the final decision for implementation. However, the fact that stochastic simulations rely on running a large number of replications to cope with the randomness and obtain some accurate statistical estimates of the system outputs, has posed a serious issue for using this kind of multiobjective optimization in practice, especially with complex models. Therefore, the purpose of this study is to investigate the performance enhancements of a reference point based evolutionary multiobjective optimization algorithm in practical production systems improvement problems, when combined with various dynamic resampling mechanisms.
Design/methodology/approach
Many algorithms consider the preferences of decision makers to converge to optimal tradeoff solutions faster. There also exist advanced dynamic resampling procedures to avoid wasting a multitude of simulation replications to nonoptimal solutions. However, very few attempts have been made to study the advantages of combining these two approaches to further enhance the performance of computationally expensive optimizations for complex production systems. Therefore, this paper proposes some combinations of preferencebased guided search with dynamic resampling mechanisms into an evolutionary multiobjective optimization algorithm to lower both the computational cost in resampling and the total number of simulation evaluations.
Findings
This paper shows the performance enhancements of the referencepoint based algorithm, RNSGAII, when augmented with three different dynamic resampling mechanisms with increasing degrees of statistical sophistication, namely, timebased, distancerank and optimal computing buffer allocation, when applied to two realworld production system improvement studies. The results have shown that the more stochasticity that the simulation models exert, the more the statistically advanced dynamic resampling mechanisms could significantly enhance the performance of the optimization process.
Originality/value
Contributions of this paper include combining decision makers’ preferences and dynamic resampling procedures; performance evaluations on two realworld production system improvement studies and illustrating statistically advanced dynamic resampling mechanism is needed for noisy models.
Keywords
Citation
Ng, A., Siegmund, F. and Deb, K. (2018), "Reference point based evolutionary multiobjective optimization with dynamic resampling for production systems improvement", Journal of Systems and Information Technology, Vol. 20 No. 4, pp. 489512. https://doi.org/10.1108/JSIT1020170084
Download as .RISPublisher
:Emerald Publishing Limited
Copyright © 2018, Emerald Publishing Limited
1. Introduction
Realworld optimization problems very often involve multiple objectives that have to be considered simultaneously. In terms of production systems improvement, there are almost always at least two objectives to be considered – the targeted condition (Rother, 2009) and the cost of the improvement. As in many realworld multiobjective problems (MOP), it is obvious that these two objectives are in conflict with each other, so that an improvement in one objective can only be obtained at the expense of degradation of the other objective. MOPs can be readily solved by a priori approaches which transform the problems into some singleobjective ones by formulating the objective functions into single, weightedsum functions. However, this is usually not applicable if the decision maker (DM) does not explicitly know how to weigh the various objectives before any optimal alternatives are known. At the same time, it is not easy to understand the effects of the weights, in terms of correlation and nonlinear outcomes, meaning that a small change in weights can alter the solution dramatically. In contrast to any a priori approaches in which the DM has to explicitly determine their preference regarding the objectives before the optimization process, a posteriori approaches aim to find the entire set of best “tradeoff”, or socalled Paretooptimal solutions, so that the DM can decide which solution to implement after the optimization has been completed. The goal of a posteriori optimization approach is therefore to find a converged set of solutions that also feature wide diversity, to spread as much as possible in the objective space and form an optimal tradeoff curve/surface, or socalled efficient frontier (EF). This approach is desirable, as it allows the DM to obtain a complete picture about the problem under study, e.g. the relationship between the decision variables and the objectives (Bandaru et al., 2017) and provides the DM with many alternatives to choose from. However, it becomes a critical issue if the optimization process involves computationally expensive function evaluations, e.g. stochastic simulation runs on largescale, complex models.
Regarding production systems analysis, stochastic simulation is not only a popular tool for the evaluations of long, complex and realworld production systems, but probably the only feasible choice, especially when the processing times and downtimes follow nonexponential or nonnormal distributions (Negahban and Smith, 2014). Stochastic simulation is the only available choice for researchers and practitioners in industry alike if more complex flows and other types of variability (e.g. setups) are included in the study of unbalanced production lines. As claimed by Tempelmeier (2003), “If quantitative performance evaluation is carried out at all (in the industry), then in almost any case simulation is the only tool used.” Wang and Chatwin (2005) summarize the three major weaknesses of analytical/mathematical methodologies, compared to computer simulation: analytical evaluation is impractical when it encounters stochastic elements, such as many random and nonlinear operations that exist in virtually any manufacturing system; due to randomness in a dynamic system which changes with time, the mathematical modeling of a complex dynamic system requires many simplifications which may cause the models to become invalid; analytical methods are not sufficient for optimization because mathematical models can only be built with simplifying assumptions that may affect the accuracy of the optimization results. Additionally, they also point out that in many cases, one must resort to simulation even though the system under study is analytically tractable because some performance measures of the system have values that can only be computed by running simulations. Having listed all these advantages offered by simulation, it is still an important issue that stochastic simulations rely on running a large number of replications to cope with the randomness and obtain some accurate statistical estimates at a specified level of confidence (Chen and Lee, 2010). While the recent advance in computing technology has greatly ameliorated this issue, Kleijnen (2015), however, correctly points out that modelers also demand greater modeling complexity when the technology advances. Therefore, it is not difficult to understand that the computational cost grows very quickly, with regard to handling a stochastic MOP in which the aim is to both generate a widespread Paretooptimal solution set and, at the same time, ensuring the statistical soundness of the optimal solutions. For instance, a simulation model that takes a CPU 3 min to execute a replication run will cost 250 h, if 10 replications are required and the optimization takes 5,000 evaluations, which for any DM would simply be a prohibitively long waiting period. As a matter of fact, it is not surprising that in a real industrial decisionmaking situation that involves a short lead time, either the simulation accuracy or the optimality of the solutions, or even both, have to be compromised for the DM if such a stochastic MOP approach has to be used. To address this issue, many techniques that aim to obtain more accurate results with a reduced total number of replications and/or evaluations have been proposed in the literature. In terms of reducing the number of replications, most of these techniques involve the idea that replications should be dynamically allocated to solutions that are believed to be nearoptimal. This kind of dynamic allocation of replications is also commonly called dynamic resampling (DR), as the mechanism attempts to identify which solutions should be given more samples during the optimization process. Optimal computing buffer allocation (OCBA) (Chen et al., 2008) is the best example of these advanced statistical techniques that can estimate the optimal number of replications to the solution. These can be combined with various evolutionary algorithms to optimize the accuracy of the optimal solutions. For example, a confidencebased DR integrated into the wellknown NSGAII algorithm to improve the confidence of the Pareto ranking was proposed in Syberfeldt et al. (2010).
On the other hand, in terms of reducing the total number of evaluation for faster convergence toward the EF, a guided search that considers the preference(s) of the DM during the optimization seems to be a promising approach proposed by EMO researchers (Branke, 2008). While it may be impractical for a DM to completely specify the tradeoffs between the different objectives as suggested in a priori approaches, it is also reasonable to assume that a DM has at least some rough idea about what he/she would prefer. In the literature related to multicriteria decisionmaking (MCDM), there is an abundance of methods that a DM can use to interactively specify his/her preference(s) regarding the objectives. With such interactive methods, procedures of some iterative algorithms can be repeated so that a DM can progressively provide the preference information required to find the most preferred solution. Miettinen et al. (2008) identify three ways in which a DM can provide preference information to the interactive methods: tradeoff information; classificationbased and reference point (RP) approach. In general, tradeoff information concerns determining, to a certain quantity, how much improving a certain objective will deteriorate another one. Such a tradeoff relationship can be measured objectively when the move from one feasible solution to another one is considered, socalled objective tradeoff. It can also be acquired subjectively when the DM determines how much of the value of one objective can be sacrificed to desirably improve another one at a specified quantity. Therefore, two different schemes, whether to determine at each iteration the objective tradeoffs and present them to the DM, or ask the DM to provide the subjective tradeoff information, to find a Paretooptimal solution, are developed in different tradeoffbased interactive methods. The latter, subjective methods usually assume the DM has an implicit value function, v, so that the aim of the iterative algorithms is to use the progressively acquired tradeoff information to approximate and maximize v.
With classification methods, the DM indicates his/her preferences by classifying which objective function should be improved and which one could be impaired or relaxed from its current values (Miettinen et al., 2008). For example, the satisficing tradeoff method STOM, proposed by Nakayama (1995), asks the DM to classify the objective functions into three classes: I^{≤}, values should be improved; I^{≥} values can be relaxed and I^{=}, those values are accepted. The DM has to specify the desirable aspiration levels for the functions in I^{≤} and the upper bounds for functions in I^{≥}.
In comparison, while classification methods assume that some objective function must be allowed to become worse, a reference point to be selected in the RP approach can be chosen more freely. The RP approach requires the DM to specify RP(s) in terms of aspiration and reservation levels for all objective functions. Paretooptimal solutions can then be generated through maximizing some achievement scalarizing function (ASF) which is a nonlinear aggregation of the objective functions. As an interactive approach, the DM is supposed to learn about the decision situation and to explore different interesting parts of the Paretooptimal set through freely modifying the RPs. This learning aspect through interactively varying the preferences is more flexible, compared to the two other methods. On a 2D/3D plane, the RP approach would also allow the DM to easily verify whether the obtained results correspond to his/her preferences visually (Mkaouer et al., 2013). At the same time, the RP approach makes it possible to find Paretooptimal solutions that correspond to multiple preference conditions simultaneously. Additionally, RPbased methods are applicable to MOPs with a large number of objectives (>3), a large number of variables, with linear or nonlinear constraints (Chikumbo et al., 2012), partly because the preference relation gives a finer order of vectors of objective space than only the Pareto dominance relation (LópezJaimes and Coello Coello, 2014). With such motivations, several algorithms with the capability to focus the search toward the preference areas in the objective space, specified by RPs, have been proposed.
As reviewed in this introduction, many algorithms that consider the preferences of the DM are targeted at converging faster toward Paretooptimal solutions. There are also significant efforts in proposing various DR procedures to avoid wasting an abundance of simulation replications to nonoptimal solutions. However, to our best knowledge, very few attempts have been made to illustrate the advantages of combining both of these two approaches to further enhance the performance of the algorithm, especially if the problems under study involve computationally expensive evaluations, such as production systems simulation. Therefore, this paper aims to propose a combination of a preferencebased guided search with DR into an evolutionary multiobjective optimization (EMO) algorithm to lower both the computational cost in resampling and simulation evaluations. Specifically, based on the above discussions about the advantages and flexibility offered by the interactive RP approach, we carry out an experimental, comparative study of performance enhancements of RNSGAII augmented with three DR mechanisms when applied to realworld production systems improvement. These DR mechanisms are: timebased, distancerank and OCBA (the order follows the increasing degree of statistical sophistication). The results will also be compared with the optimizations generated with the original RNSGAII and NSGAII.
In the remainder of this paper, we first introduce how a production systems improvement problem can be formulated into a MOP in Section 2, including an outline of how a RPbased algorithm can enhance the accuracy and efficiency when a lowbudget computational cost (i.e. a limited number of resamples) is sought or allowed. Then, in Section 3, we describe the three DR mechanisms selected for examination in this paper. Section 4 introduces RNSGAII and how it can be augmented with the selected DR mechanisms. Experimental results when applied to the DRaugmented RNSGAII in two realworld production system simulation models are presented in Section 5, followed by conclusions and future work in Section 6.
2. Production systems improvement as multiobjective problems
Consider a general MOP that consists of m objective functions, f_{i}(x), i = 1,…, m, which can be minimized and maximized with x as a decision (design) vector, comprising n decision variables, x_{i}, within their respective lower bounds (x_{i}^{L}) and upper bounds (x_{i}^{U}). Mathematically:
For a continuous MOP, we call
In many MOPs, where the objectives f_{i}(x) are in conflict with each other, finding a single best optimal design is impossible because improving one objective would deteriorate the other (Deb, 2004). This gives rise to the concept of EF that denotes the best tradeoff in ℤ^{m} with respect to f_{i}(x). The definition of dominance is essential for an optimization to find and compute x that constituent EF: considering only Min{f_{1}(x),…, f_{m}(x)}, a design vector x_{1} is said to dominate another solution vector x_{2}, denoted as x_{1}≼x_{2}, iff f_{1}(x_{1}) ≤f_{i}(x_{2}) ∀i ∈ {1,…, m}, but ∃j∈{1,…, m} s.t. f_{1}(x_{1}) < f_{i}(x_{2}). A decision vector x* ∈Φ is Paretooptimal if ∄x that dominates x*. In other words, the Paretooptimal set, PS, consists of x* that are nondominated to each other in ℤ^{m}. Equivalently, an objective vector z*∈ℤ^{m} is Paretooptimal if the decision vector correspondent to it is Paretooptimal (Miettinen, 1999). With the above formulation, it is not difficult to understand that an optimizationbased design approach will involve a process for finding the PS that forms EF so that a DM can make a final selection on which z* to choose. In the literature, this is referred to as the postoptimality decisionmaking process that involves some higherlevel information, very commonly including some preferences of the DM which were not formulated into the objectives in the MOP.
Production systems improvement problems can be formulated readily as MOPs, as shown earlier in Pehrsson et al. (2011, 2013) and more recently in Pehrsson et al. (2016). The basic idea proposed was based on an observation that many decisionmaking situations in production system improvement projects can be effectively formulated into MOPs. The primary objective of the MOP is usually related to a key performance index (KPI), such as system throughput or manufacturing leadtime; the second objective of the MOP can be formulated as a function of investments, or steps of changes, or socalled improvement actions, needed to improve the KPI of the system (Pehrsson et al., 2016). For example, if the production system throughput (τ) is the primary KPI as the improvement target, so that f_{1}(x) = maximize(τ), then we can define the total number of changes, i.e. improvement actions, as minimize(f_{2}(x)) in the MOP. If we consider three types of discrete, twolevel decision variables: {α_{i}, β_{i}, γ_{i}}, i = 1,…, n, where n = number of workstations so that each variable can either be set to its original value or to some improved value, to represent the possible improvements for availability, processing times and repair times, respectively. For a workstation i, such that, e.g. α_{i} = {0.9, 0.95}, β_{i} = {65 s., 60 s} and γ_{i} = {20 min., 10 min.}, then mathematically, the second objective in the MOP can be formulated as a summation function:
α^_{i} = {0, 1}, β^_{i} = {0, 1}, γ^_{i} = {0, 1} as the binary variables which represent whether correspondent decision variables, α_{i}, β_{i} and γ_{i} retain their original values or are changed to some improved values. For instance, with the above example, if α_{i} is not changed (i.e. 0.9), then α^_{i} = 0. Otherwise α^_{i} = 1 if α_{i} is improved to 0.95. Similarly, β^_{i} = 1 if β_{i} is cut to 60 s and γ^_{i} = 1 if γ_{i} is reduced to 10 min. It is very important to note that improved values for β_{i} and γ_{i} mean their values are reduced.
With this formulation, the second objective function is a discrete function that varies in the range [0, 3n] because the maximum number of changes can only be 3n. Figure 1 illustrates this biobjective problem graphically, showing how the Paretooptimal solutions in the objective space can be used to support the decisionmaking in regarding the selection of the optimal (minimal) changes required to improve the system to achieve the target throughput.
In Figure 1, the original throughput, τ_{0}, represents the current condition when no improvement has been made, i.e. f_{2}(x) = 0. If the target of the KPI is τ_{t}, then in the objective space, z*_{1} denotes the optimal number of changes (c_{t}) required to improve the production system from τ_{0} to τ_{t}. As illustrated, z*_{t} dominates z_{1} as they have the same f_{2}(x) value. However, a higher τ can be achieved with z*_{t}. On the other hand, in this MOP, τ_{max} can also be a very interesting value for the DM because it denotes the maximum τ that can be obtained for the production line under study when all possible changes that can be made in this MOP are considered. Note that τ_{max} can be further increased if the processing times of the workstations, β_{i}, can be further reduced to some values not considered in the MOP. However, as reducing physical processing times can only be done to a certain extent and α_{i} is always <1 (max. availability 100 %) and γ_{i} ≥ 0, there is always a maximum practical τ in a real production system. Any solutions beyond z*_{max}, e.g. z_{2}, which give a τ ≤ τ_{max} with f_{2}(x) value > c_{max} are dominated by z*_{max} because they require more, unnecessary changes and would certainly be considered as undesirable by the DM.
In a realworld decisionmaking process, finding z*_{t} is a crucial process especially if the DM already has a clear target KPI value in mind, despite the fact that how this preference information in the objective space can be provided for the optimization process is an interesting issue that has induced research into different methods. On the other hand, although the DM may have the target τ_{t} clearly in mind, finding a single z*_{t}, which could be done with a singleobjective optimization process can be undesirable because the DM may want to explore more possible solutions around (c_{t}, τ_{t}), for several reasons. First, it may be preferable to select a slightly higher value than τ_{t} to add a safety factor to τ_{t} to compensate some possible evaluation error of the model (simulation) evaluation (see also the discussion below on the estimation of f_{i}(x)) or to address what is called preference correlation. This is a situation where a DM would prefer to have one performance measure higher or lower than intended in order to compensate for another performance measure that is being higher or lower than normal (Rosen et al., 2007). Such a preference correlation situation can occur easily if the DM has explored the existence of some optimal tradeoff solutions around τ_{t} that exert “knee” points in MOP literature (Branke et al., 2004; Deb and Gupta, 2011), shown as z*_{k} in Figure 2. As illustrated in Figure 2a, despite that τ_{t} is the target, but as the cost (required changes, c_{t}–c_{k}) of advancing from z*_{k} to z*_{t} is high, the small gain with τ_{t} – τ_{k} may not justify so that the DM might relax τ_{t} and select z*_{k} as the chosen solution for implementation. On the other hand, for the “knee” case in Figure 2b, the cost of increasing τ_{t} to τ_{k}, represented as c_{k}  c_{t}, is so small so that the DM would be very comfortable selecting z*_{k}. In both cases, the DM has gained better knowledge about the behavior of the system under study and the decision can be made with greater confidence which could not be easily obtained without finding more tradeoff solutions. In this way, the DM adjusts τ_{t} according to some broader information gained in a multiobjective optimization process, which cannot be easily acquired from singleobjective optimization, and is referred to as “knowledgedriven optimization” (Bandaru et al., 2017).
With a simulationbased optimization approach, all function evaluations are done through the simulation model, giving only an estimation of f_{i}(x). In other words, F(x) can only be estimated by F^(x) through the performance values obtained from simulation replication runs. Stochastic simulation will result in different F^(x) in each replication, so the estimated output objective function values used for the optimization are actually F^(x) = {(
3. Guided search using reference points
Using a reference point (RP) to indicate a DM’s aspiration levels in different objectives is perhaps the most important way of providing preference information in an interactive MCDM procedure (Branke, 2008). As suggested by Wierzbicki (1980, 1982, 1986), the RP approach can be used to find a weakly, єproperly or Paretooptimal solution closest to a supplied RP of aspiration level, reflecting desirable values for the objective function, based on optimizing an achievement scalarizing function (ASF). Given a RP, z¯, for MOP of minimizing (f_{1}(x), f_{2}(x), …, f_{m}(x)) as in equation (1), the following singleobjective optimization is solved:
Here, w_{i} is the ith component of a chosen weight vector used for scalarizing the objectives, in which each component of the reference point represents the aspiration levels that the decision maker requires for each objective. For a chosen RP, the closest Paretooptimal solution, in the sense of the weightedsum of the objectives, is the target solution found by a RP procedure. If the DM is interested in biasing some objectives more than others, a suitable weight vector can be used with the RP selected.
When using a RP search procedure, the RP can be given interactively, preferably with the approximated nadir and ideal objective vectors presented to the DM to illustrate the ranges of ℤ^{m}. The interactive procedure starts with the DM specifying a RP, z¯. By solving the ASF with a given z¯ and w_{i}, a z* will be located on EF as shown in Figure 3. If the DM is not satisfied with the Paretooptimal solutions generated, a new reference point can be suggested and the above procedure is repeated. In this way, reference points can be created, adapted or deleted interactively during the optimization run. By repeating the procedure with different reference points, the DM tries to evaluate the region of Paretooptimality, instead of one particular Paretooptimal point.
To make a RP procedure more interactive and useful in practice, Wierzbicki (1980) also suggested the generation of m more new RPs by using the following equation:
where d = ‖z¯ − z*‖ and e_{j} is the jth unit vector for j = 1,…, m. By using this calculation, z¯_{j} represents the newly generated RP along the jth coordinate direction vector. The new m RPs can be used to find m more z* as shown in Figure 3. In this way, the new m RPs will form new ASFs so that m additional z*, instead of only one single Paretooptimal solution, can be found. As d is proportional to the distance between z¯ and z*, an interesting effect with this RP approach is that a RP which is further away from the EF (e.g. z¯_{1} in Figure 3) results in larger spaced auxiliary reference points and thus can yield solutions with wider spread (s_{1}), compared to z¯_{2} which is closer to EF and finds lessspaced auxiliary reference points (hence solutions with less spread, s_{2}).
Depending on the ASF selected, the solutions found by a RPbased method can be weakly or єproperly Paretooptimal. For example, to avoid the generation of weakly Paretooptimal solutions, a sufficiently small positive scalar ρ, can be incorporated into equation (5):
On the other hand, one common approach to avoid using the weights, w_{i}, as additional preference information, is to calculate w_{i} by normalizing each objective function (Miettinen, 1999; LópezJaimes and Coello Coello, 2014):
Here z** ∈ ℤ^{m} is called a utopian, infeasible objective vector whose components are formed by
It is interesting to note that without using the absolute value of (f_{i}(x)  z¯_{i}) in the ASF function [equation (3)], the RP approach can be used for both feasible RP and infeasible RP (infeasible in every component i). A feasible RP means that ∃ z∈ℤ^{m} s.t. z ≼ z¯, i.e. there is at least one solution that dominates z¯. In other words, an infeasible point is a z¯ that cannot be achieved from any solution in the feasible objective space. If a reference point is feasible and is not a Paretooptimal solution, the DM may then be interested in knowing about solutions which are Paretooptimal and close to the reference point. On the other hand, if z¯ is an infeasible RP, the DM would be interested in finding Paretooptimal solutions which are close to it. As suggested in the MCDM literature (Miettinen, 1999), a DM should never be pessimistic when choosing a RP to indicate his/her preference. Therefore, starting with some infeasible RP is recommended. Figure 3 also illustrates that an infeasible RP which is further away from the EF, such as z¯_{1}, is desirable in the early interactive process, to explore a wider region of the EF. In general, as mentioned previously, how a DM can interactively vary the preferences during the optimization process to learn the behavior of the system under study through exploring the objective space is an important motivation for the use of a RP approach in the guided search, compared to other MCDM methods.
As described in Goulart and Campelo (2016), any Paretobased algorithm can be adapted to incorporate the RP approach. The earliest use of RP to guide an EMO algorithm was proposed by Fonseca and Fleming (1993, 1998). In the current paper, we decided to investigate the combined effects of the DR and RP approaches on an existing, efficient and wellknown EMO algorithm, namely, NSGAII. Therefore, the RP extension of NSGAII, i.e. RNSGAII (Deb et al., 2006), is further extended with three different DR mechanisms for running the experiments and then comparing the results with the original RNSGAII (no DR) and NSGAII (with or without DR but no RP). In the next section, we introduce the three DR mechanisms before presenting the details of the DRaugmented RNSGAII in Section 5.
4. Dynamic resampling
The simplest resampling mechanism is to distribute the total amount of the replication budget equally among the evaluated solutions. This is a popular method among simulationbased optimization practitioners because it requires minimal implementation effort. If r denotes replication runs for a single solution, x, and B is the total replication budget, then for such a static resampling mechanism, N = B/r solutions can be evaluated. While increasing r can improve the accuracy of evaluating F^(x), this will reduce N. At the same time, the replications that are assigned to evaluate nonoptimal solutions will be wasted. In comparison to such a static resampling mechanism, the basic idea of DR is to dynamically allocate r_{i} for a solution i, in a dynamic basis. An intuitive DR procedure would be to reduce the standard error until it drops below a certain userdefined threshold. The required sampling budget for obtaining the reduction in threshold standard error (e_{t}) can be calculated as:
However, as the sample mean changes as new samples are added, this oneshot sampling allocation might not be optimal. The number of fitness samples drawn might be too small to reach the error threshold, in case the sample mean has been shown to be larger than the initial estimate. On the other hand, a oneshot strategy might add too many samples, if the initial estimate of the sample mean is too big. Therefore, DR is often done sequentially. Through this sequential approach, the number of required samples can be determined more accurately. Therefore, DR techniques can be classified into either oneshot sampling strategies or sequential sampling strategies. Some other sequential samplingbased algorithms proposed can be found in Branke and Schmidt (2004) and BartzBeielstein et al. (2007).
In the following subsections, we introduce the three different DR mechanisms based on various combinations of optimization time, distance to the RP and rank in the nondomination sorting of the solutions. In the following, we denote:
B = maximum overall number of simulation runs, i.e. total budget, and B_{t} = current overall number of simulation runs and B_{F} = final samples for replicating the final set of solutions.
The number of sample(s) that should be allocated for a solution x, calculated in a DR iteration, is b_{s}. Minimum and maximum number of samples for an individual solution: b_{min} and b_{max}. In other words, for static resampling, each solution will have a constant sampling budget, b_{s} = b_{min} = b_{max}.
Acceleration parameter for the increase of the sampling budget: a > 0 – increasing a decreases the acceleration of the sampling budget, decreasing a increases the acceleration.
The calculated normalized sampling need, s, is discretized in the following way which guarantees that b_{max} is already assigned for s < 1:
4.1 Timebased resampling
Timebased dynamic resampling or TDR (Siegmund et al., 2013) is a generally applicable DR mechanism which can be used with any optimization algorithm. It assumes that the accuracy of the objective values should be increased gradually when the optimization converges so that the maximum number of replications should be used only toward the end of the optimization run. This hypothesis was proven in Aizawa and Wah (1993, 1994) for singleobjective EAs. Timebased dynamic resampling allocates a small sampling budget at the beginning and a high sampling budget towards the end of the optimization. The strategy of this resampling algorithm is to support the algorithm, when the solutions in the population are close to the Paretofront, and to save sampling budget at the beginning of the optimization when the solutions are still far away from the Paretofront. The normalized timebased resampling need, s^{T}, is calculated as:
Timebased resampling is a DR technique that does not consider the variance of the solutions. It makes only one sampling decision for each solution according to the elapsed time, i.e. a oneshot allocation strategy. However, if a solution survives into the next generation of an evolutionary algorithm, the decision is made again. The total number of samples available for optimization is limited by the final samples, B_{F}, added to the last population. For timebased resampling, this means the elapsed time needs to be adjusted so that 100 per cent of the time has already been reached when B − B_{F} samples have been used. The modified allocation function is shown in equation (10). This function is used in all the other resampling algorithms that apply the elapsed optimization time as one of their resampling criteria.
4.2 Distance rankbased dynamic resampling (DR2)
An earlier attempt to combine multiple different resampling criteria in a resampling algorithm to combine the rankbased dynamic resampling and the distancebased dynamic resampling (DDR) strategies, so called distancerank dynamic resampling (DR2), was reported in Siegmund et al. (2015). It uses four different factors to determine the resampling allocation for individual solutions: Paretorank, time, reference point distance and progress. Although the elapsed optimization time can be combined with any other resampling criterion with little effort, as seen with ranktimebased dynamic resampling, the Paretorank and the reference point distance are two truly different resampling criteria. Therefore, the term hybrid is emphasized for the DR2 algorithm.
The allocation function for rankbased DR is given in equation (11). R_{s} stands for the Paretorank of solution s and R_{max} is the maximum Paretorank of the population. Equation (12) shows the distancebased allocation of the DDR sampling algorithm. The reference point distance of solution s is denoted as d_{s} and δ¯ is the current minimum reference point distance achieved by any solution in the population. The mechanisms regarding how we considered the elapsed optimization time and the population progress toward the reference point in combination with the reference point distance in equation (12), are described in Siegmund et al. (2013).
Equation (13) describes the sampling allocation of DR2. Like ranktimebased resampling, the minimum of both the normalized sampling need of distancebased resampling and rankbased resampling is used to create a combined sampling allocation. However, the distancebased sampling need is not calculated individually for each solution. Instead, DR2 identifies the solution s_{m} closest to R, d_{m} = min_{k}{d_{k}}, and the normalized sampling need for
4.3 Distancebased OCBAm dynamic resampling (DOCBAm)
With the distance information available toward the preferred area, a resampling algorithm can allocate samples to enable an optimal comparison between solutions regarding distance. A comparison of solutions within a stochastic simulation scenario requires considering information about their objective variance. With this information, a statistically significant calculation can be made regarding which one of two solutions is closer to the reference point. For this purpose, we proposed the use of a Ranking and Selection algorithm called OCBAm (Chen et al., 2008) to determine the solution subset of size m which contains the closest solutions to the RP with a desired accuracy. The Optimal Computing Budget Allocation is an R&S algorithm, which allocates samples to the considered designs according to a scalar objective value. While the original OCBA procedure by Chen et al. (2000) allocates samples to identify the best design according to the single scalar objective value, the OCBAm algorithm described in (Chen et al., 2008) extended the original algorithm to identify the best m designs. As in our case we treat the scalar objective value for OCBAm as the distance to the RP and its variance according to all calculated objective value samples so far, the variation is therefore called distancebased OCBAm dynamic resampling or DOCBAm in below.
The sampling allocation problem that is to be solved by DOCBAm is stated in equation (14) (Siegmund et al., 2016). The corresponding notation is given in Table I. It is solved asymptotically, which means DOCBAm recalculates the allocation iteratively, after having received new objective value evaluation data.
where:
Consequently, in the list of solutions sorted by reference point distance, there are m − 1 solutions that are closer to z¯ than for solution m. The index k_{m+1} belongs to the next best solution in the ordered list. This solution and all other worse solutions are not selected into the next population. Note also that n_{0} is for bounding N_{k}; it is related to b_{min} but if b_{min} = 1 then n_{0} should be bigger to estimate the standard deviation.
Based on the Karush–Kuhn–Tucker condition, the allocation is optimal if the following holds for all solutions i and j, i ≠ j (Chen et al., 2008):
This allocation can be achieved by selecting one solution arbitrarily (let it be solution s) and putting all other in relation to it, as in the equation below:
The allocation for solution i, with respect to the problem constraint,
Equation (17) can be used because N_{i} should be assigned its correct share of N_{total}. For more details of the iterative DOCBAm procedure, readers are referred to Siegmund et al. (2016).
5. DRAugmented RNSGAII
RNSGAII is based on the nondominated sorting genetic algorithm II, NSGAII (Deb et al., 2002), which is a widely used and representative EMO algorithm. NSGA II sorts the solutions in population and offspring into different nondominated fronts. Those selected comprise all solutions in those best fronts that completely fit in size into the next population – if population size is P and the number of solutions in those selected best fronts is Q, where P ≥ Q. However, in many cases Q > P, so that not all fronts can be accommodated into the next population. Then a subset of solutions which have the largest distances to their neighbors will be selected into the next population. This selection mechanism is called crowding distance and it guarantees that the result population will be diverse. The parental selection for offspring generation follows this fitness hierarchy. After selection has been completed, offspring solutions are generated by tournament parental selection, crossover and mutation. The offspring are evaluated and the selection step is performed again. The RNSGAII algorithm replaces the crowding distance operator by the distance to userspecified RP(s). Solutions that are closer to a RP will receive a higher selection priority. As shown previously, the DM defines the RPs in areas of interested to find more alternative solutions in some limited regions of the objective space. However, in order not to lose all diversity, a clustering mechanism with a minimum objective vector distance ε is used. Solutions that are too close to each other (d < ε) are considered with a lower priority as illustrated in Figure 4. The same fitness hierarchy, dominance first, then reference point distance and εclustering, is used for parental selection.
The RPs can be created, adapted or deleted interactively during the optimization run, taking effect each time a new generation is started. In case all RPs have been removed, RNSGAII will continue the optimization as an instance of the regular NSGAII algorithm. As RNSGAII uses nondomination sorting, it tends to prioritize population diversity before convergence to the RPs. Therefore, extensions have been proposed which limit the influence of the Paretodominance and allow the algorithm to focus faster on the RPs (Siegmund et al., 2012).
The RNSGAII algorithm is described in pseudocode in Algorithm 1; the steps that are augmented with a DR mechanism are highlighted. As RNSGAII is a populationbased EA, the DR algorithm should aim at supporting the correct selection of the best solutions for the next generation.
Compared to its predecessor, NSGAII, less reporting has been done on the application of RNSGAII in solving realworld problems. One application that gained some publicity was reported in Chikumbo et al. (2012) and concerned solving a land use management problem in New Zealand, which comprised 14 objectives and 3,150 integer variables. Some other applications using RPbased EMO, in particular RNSGAII that can be found are logistics network design (Rajabalipour et al., 2014), reversible logic circuit design (Wang and Wang, 2014) and workflow scheduling in computing grid (Garg, 2011; Garg and Singh, 2012; Navaz and Ansari, 2012) or cloud computing (Verma and Kaushal, 2015). In manufacturing, there are very few reported applications of RNSGAII, except, for instance, tool sequence and parameter optimization of rough milling (Churchill et al., 2013a, 2013b). To our best knowledge, the application studies presented in the next section are among the first ones testing RNSGAII on complex production system models.
6. Industrialscale applications and results
This section addresses the experimental setup and results of using the DRaugmented RNSGAII on two realworld industrial application studies, both found in the automotive industry. The optimizations have the same aim – to effectively identify the exact areas of optimal improvements required to reach the desired target condition within a limited simulation budget, B = 5000.
The first model represents a complex assembly line comprising 71 workstations (i.e. n = 71), 72 buffers, pallets for carrying two variants for assembly and disassembly processes. The reduction for each α_{i} (processing time) is 20%; improvement in availability for β_{i} from its current value to 0.98 and repair times (γ_{i}) reduced by 50%, i = 1, , n. The second complex production system is an automotive engine machining line that includes multiple parallel sections, gantry robots, machining centers and assembly stations which conduct multiple operations. Apart from some maintenance issues that affect the availabilities of the machines, multiple variants have to be processed on the line and variations in the weekly volume contribute to the higher variability of the system compared to the previously mentioned assembly line. The crucial issue for the company was the poor production capacity of this line, far below the capacity required to maintain production; a lean improvement project was in order. However, given the size, complexity and variability of the line, locating where and what to improve, besides the effect of performing the improvements if only traditional Lean tools, e.g. value stream mapping, were to be used, was believed to be extremely difficult. For these two cases, it was obvious that not a single, but multiple improvements had to be made to achieve the targeted throughput levels. The first model has been described with more details in Bernedixen et al. (2015) where discussions on the importance of avoiding any manually tedious and errorprone optimization setup by using some software support are presented. In addition, the engine machining line under study is the same one that is presented in Ng et al. (2014). While that article is focused on the relationship between bottleneck improvements using NSGAII, the current paper mainly focuses on the results of applying the DRaugmented RNSGAII. Various versions of the simulation models have been built, but the latest complete models for the optimization studies were developed using FACTS Analyzer (Ng et al., 2007, 2008). Figure 5 illustrates the FACTS Analyzer model developed for the assembly line using standard modeling objects provided by the software, including multiple Operation, Buffer, Assembly, Disassembly and CONstant WorkInProcess (CONWIP) Control objects. As a typical discreteevent simulation model built for studying improvement, the complexity of the model is not in the products flow, but rather in the scale of the model that many attributes can be subjected to change. Using trialanderror simulation experiments to evaluate all these possible changes is not only too timeconsuming but also a waste of the potential offered by the valid model.
The five algorithm configurations in the experiments were as follows:
NSGAII with static resampling and 10 samples (Static10).
RNSGAII  Static10.
RNSGAII  Timebased Dynamic Resampling 110.
RNSGAII  Distance Rankbased Dynamic Resampling (DR2) 110.
RNSGAII  Distancebased OCBAm Dynamic Resampling (DOCBAm) 110.
Two different experiment setups were run:
An optimistic reference point closes to the ideal point, z¯_{1}.
z¯_{2}, a reference point above the maximum achievable throughput, τ_{max}.
The reference points for the assembly line were z¯_{1} = (0, 100) and z¯_{2} = (100, 110) and for the engine machining line model were z¯_{1} = (0, 100) and z¯_{2} = (400, 100).
As stated, both optimization studies were allowed to run with B = 5000. The results from five replicated optimizations are illustrated using median attainment surfaces (Knowles, 2005). The results from the two experiment setups for the two industrial application cases with the five algorithm configurations are shown in Figures 610, where the minImpmaxTP plots represent the two objectives used in all the optimizations, namely, Minimize(Improvement Actions) as defined in equation (2) and Maximize(Throughput).
For the assembly model with two different RPs (Figures 6 and 7), distinguishing the performance among the three DR mechanisms seems to be difficult. Unexpectedly, the difference between RNSGAII Static10 and NSGAII Static10 is not significant. However, what can be concluded here is that, in general, DRbased algorithms perform significantly better than Static Resampling (Static10). For the case with RP z¯_{2}, DOCBAm performs best, and there is no clear advantage with using Timebased DR or DR2 (Figure 7). RNSGAII Static10, however, is clearly better than NSGAII for the solutions close to z¯_{2}.
On the other hand, the result variances from using different DR mechanisms are more pronounced for the engine machining line. In Figure 8, the attainment surface results for the optimistic RP (z¯_{1}, closer to the ideal point) have shown that the best performance can be achieved by the RNSGAII with DR2 as the DR mechanism, followed by DOCBAm and timebased DR. The higher throughput performance obtained with RNSGAII DR2 can be statistically verified using the confidencebased dominance postoptimality analysis, first introduced in Ng et al. (2008) and later used for a systematic comparison of production control mechanisms in the objective space using NSGAII (Ng et al., 2012). The basic idea of confidencebased dominance is to consider the statistical significance when determining whether an objective function of a solution, i.e. a design vector x_{1}, is better than another solution vector x_{2}, before a dominance relationship (e.g. x_{1} ≼ x_{2}, iff f_{1}(x_{1}) ≤f_{i}(x_{2}) ∀i ∈ {1,…, m}, but ∃j ∈ {1,…,m} s.t. f_{1}(x_{1}) < f_{i}(x_{2})) can be concluded as described earlier in Section 2. In comparison to the minCTmaxTP, i.e. Minimize(Cycle Time) against Maximize(Throughput) analysis considered in (Ng et al., 2012), the analysis for the minImpmaxTP plots in the current paper is simpler because only the objective function value of system throughput is subjected to uncertainty from stochastic simulation. In other words, the Welch difference test used in (Ng et al., 2012) and other simulation literature is applied only to f_{2}(x) as below:
where f_{2}(x_{1}) and f_{2}(x_{2}) represent the mean values of the objective function, f_{2}(x), (i.e. maxTP) of the two solutions, x_{1} and x_{2}, respectively, under the statistical difference test; σ_{1} and σ_{2} are their standard deviations, calculated using n_{1} and n_{2} as the numbers of replications, respectively, which set by the algorithms for generating the solutions. The studentt distribution value to be set for a (1 − α)% confidence has to be estimated with a degree of freedom f^, computed and rounddown using the following equation:
Provided that the simulation outputs conform to normality, the use of the Welch difference test is extremely useful for the comparison of two solutions, x_{1} and x_{2}, generated with different numbers of replications (i.e. n_{1}≠ n_{2}) because of the applications of DR during the optimization processes. By making a complete pairwise comparison with all the Paretooptimal solutions, z*_{t}, generated with different algorithms in different, multiple optimization runs, the outstanding performance of DR2 in finding solutions that have higher throughput in the region closer to z¯_{1} can be statistically verified by the visualization of the significantly nondominated solutions as shown in Figure 9. In this figure, only the solutions that are statistically better in their minImpmaxTP values, among all the generated z*_{t}, after the Welch difference test, are visualized. Error bars showing the ± 95% standard errors of the significantly nondominated solutions obtained from the confidencebased dominance postoptimality analysis are also plotted with these solutions for representing the variability of the simulation outputs.
For the second scenario with z¯_{2} as the RP, similar results have been achieved, as shown in Figure 10. However, in this case, the performance differences between the DR mechanisms on RNSGAII are insignificant and therefore we can assume they perform equally well. The significantly better results of DR algorithms, especially DR2 on the engine machining line are due to the higher stochasticity of the model (and the real line) caused by the higher variantdependent processing times at the workstations and lower machine availabilities. In general, RNSGAII outperforms NSGAII, while the DR algorithms on RNSGAII outperform static resampling.
7. Conclusions and future work
This paper has discussed and illustrated the advantages of combining the RPbased approach and DR mechanisms to enhance the performance of EMO algorithms for the optimizations of production system improvement that involve stochastic, computationally expensive simulation evaluations. Specifically, experimental, comparative evaluations of the performance enhancements on RNSGAII augmented with three different DR mechanisms, with increasing degrees of statistical sophistication, namely, timebased, distancerank and OCBA, which were applied to an assembly line and an engine machining line have been conducted. The conflicting optimization objectives were the throughput of the lines and the necessary improvements required to enhance them. The best results illustrated in the attainment surface plots were achieved by using RNSGAII with DR2 as the DR mechanism, followed by DOCBAm and timebased DR. As a general conclusion, DRbased algorithms perform significantly better than algorithms with static resampling. In addition, the more stochasticity that the simulation models exert, the more the statistically advanced DR could significantly enhance the performance of the optimization process. As a result, decision makers can find better tradeoff solutions with lower computational budgets. Future work involves extensive experiments on both hypothetical and realworld models that exert various levels of stochasticity using the DRbased algorithms described in this paper.
Figures
Notation for DOCBAm

Mean of distance of the objective vector of sol. s_{k} to z¯, d_{k}=d(s_{k}, z¯) 
σ_{k}  Standard deviation of the stochastic variable d_{k} 
N_{k}  No. of samples for solution k 
σ¯_{km+1}  Standard deviation of the mean d¯_{k} s.t.

n_{0}  Minimum number of samples for all solutions 
N_{total}  Total budget samples allocated for this iteration 
S  The set of all considered solutions (e.g. parents and offspring in an EMO) 
S_{m}  A subset of S which contains the best m solutions (closest to z¯) 
References
Aizawa, A.N. and Wah, B.W. (1993), “Dynamic control of genetic algorithms in a noisy environment”, in Proceedings of the Fifth International Conference on Genetic Algorithms, pp. 4855.
Aizawa, A.N. and Wah, B.W. (1994), “Scheduling of genetic algorithms in a noisy environment”, Evolutionary Computation, Vol. 2 No. 2, pp. 97122.
Bandaru, S., Ng, A.H.C. and Deb, K. (2017), “Data mining methods for knowledge discovery in MultiObjective optimization: Part a – survey”, Expert Systems with Applications, Vol. 70, pp. 139159.
BartzBeielstein, T., Blum, D. and Branke, J. (2007), “Particle swarm optimization and sequential sampling in noisy environments”, Metaheuristics – Progress in Complex Systems Optimization, Operations Research/Computer Science Interfaces Series, Springer, Boston, MA. Vol. 39, pp. 261273.
Bernedixen, J., Pehrsson, L., Ng, A.H.C. and Antonsson, T. (2015), “Simulationbased multiobjective bottleneck improvement: towards an automated toolset for industry”, in Winter Simulation Conference (WSC), IEEE, pp. 21832194.
Branke, J. and Schmidt, C. (2004), “Sequential sampling in noisy environments”, in Proceedings of the 8th International Conference on Parallel Problem Solving from Nature, 2004, Birmingham, Lecture Notes in Computer Science, Vol. 3242, pp. 202211.
Branke, J., Deb, K., Dierolf, H. and Osswald, M. (2004), “Finding knees in multiobjective optimization”, in International Conference on Parallel Problem Solving from Nature, Springer Berlin Heidelberg.
Branke, J. (2008), “Consideration of partial user preferences in evolutionary multiobjective optimization”, in Multiobjective Optimization: Interactive and Evolutionary Approaches, Branke, J., Deb, K., Miettinen, K., Slowinski, R. (Eds.), SpringerVerlag, pp. 157178.
Chen, C.H., He, D., Fu, M. and Lee, L.H. (2000), “Efficient simulation budget allocation for selecting an optimal subset”, Informs Journal on Computing, Vol. 20 No. 4, pp. 579595.
Chen, C.H., Lin, J., Yucesan, E. and Chick, S.E. (2008), “Simulation budget allocation for further enhancing the efficiency of ordinal optimization”, Discrete Event Dynamic Systems, Vol. 10 No. 3, pp. 251270.
Chen, C.H. and Lee, L.H. (2010), Stochastic Simulation Optimization – An Optimal Computing Budget Allocation, World Scientific Publishing Company.
Chikumbo, O., Goodman, E. and Deb, K. (2012), “Approximating a multidimensional pareto front for a land use management problem: a modified MOEA with an epigenetic silencing metaphor”, IEEE Congress on Evolutionary Computation, CEC 2012, pp. 1015.
Churchill, A.W., Husbands, P. and Philippides, A. (2013a), “Multiobjective tool sequence and parameter optimization for rough milling applications”, IEEE Congress on Evolutionary Computation, CEC 2013, pp. 14751482.
Churchill, A.W., Husbands, P. and Philippides, A. (2013b), “Tool sequence optimisation using preferential multiobjective search”, in Proceeding of the Fifteenth Annual Conference Companion on Genetic and Evolutionary Computation Conference Companion – GECCO, Amsterdam, The Netherlands, 2013, pp. 181182.
Deb, K., Agarwal, S., Pratap, A. and Meyarivan, T. (2002), “A fast and elitist multiobjective genetic algorithm: NSGAII”, IEEE Transactions on Evolutionary Computation, Vol. 6 No. 2, pp. 182197.
Deb, K. (2004), MultiObjective Optimization Using Evolutionary Algorithms, John Wiley and Sons Ltd.
Deb, K., Sundar, J., Reddy, U.B. and Chaudhuri, S. (2006), “Reference point based multiobjective optimization using evolutionary algorithms”, International Journal of Computational Intelligence Research, Vol. 2 No. 3, pp. 273286.
Deb, K. and Gupta, S. (2011), “Understanding knee points in bicriteria problems and their implications as preferred solution principles”, Engineering Optimization, Vol. 43 No. 11, pp. 11751204.
Fonseca, C.M. and Fleming, P.J. (1993), “Genetic algorithms for multiobjective optimization: formulation, discussion, and generalization”, International Conference on Genetic Algorithms, pp. 416423.
Fonseca, C.M. and Fleming, P.J. (1998), “Multiobjective optimization and multiple constraint handling with evolutionary algorithms  Part I: a unified formulation”, IEEE Transactions on Systems, Man, and Cybernetics  Part A: Systems and Humans, Vol. 28 No. 1, pp. 2637.
Garg, R. (2011), “MultiObjective optimization to workflow grid scheduling using reference point based evolutionary algorithm”, International Journal of Computer Applications, Vol. 22 No. 6, pp. 4449.
Garg, R. and Singh, A.K. (2012), “Reference point based multiobjective optimization to workflow grid scheduling”, International Journal of Applied Evolutionary Computation ( Computation), Vol. 3 No. 1, pp. 8099.
Goulart, F. and Campelo, F. (2016), “Preferenceguided evolutionary algorithms for manyobjective optimization”, Information Sciences, Vol. 329, pp. 236255.
Kleijnen, J.P.C. (2015), “Design and analysis of simulation experiments”, 2nd edition. International Series in Operation Research and Management Science, Springer AG: Switzerland.
Knowles, J. (2005), “A summaryattainmentsurface plotting method for visualizing the performance of stochastic multiobjective optimizers”, in Proceedings of the 5th International Conference on Intelligent Systems Design and Applications (ISDA), 810 September 2005, Wroclaw, Poland, IEEE Computer Society, pp. 552557.
LópezJaimes, A. and Coello Coello, C.A. (2014), “Including preferences into a multiobjective evolutionary algorithm to deal with manyobjective engineering optimization problems”, Information Sciences, Vol. 277, pp. 120.
Miettinen, K. (1999), Nonlinear Multiobjective Optimization, Kluwer Academic Publishers.
Miettinen, K., Ruiz, F. and Wierzbicki, A.P. (2008), “Introduction to multiobjective optimization: Interactive approaches”, in Multiobjective Optimization: Interactive and Evolutionary Approaches, Branke, J., Deb, K., Miettinen, K., Slowinski, R. (Eds.), SpringerVerlag, pp. 2758.
Mkaouer, M.W., Kessentini, M., Bechikh, S. and Tauritz, D.R. (2013), “Preferencebased multiobjective software modelling”, in Proceedings of 1st International Workshop on Combining Modelling and SearchBased Software Engineering, CMSBSE 2013, pp. 6166.
Nakayama, H. (1995), “Aspiration level approach to interactive multiobjective programming and its applications”, in Advances in Multicriteria Analysis, Pardalos, P.M., Siskos, Y., Zopounidis, C. (eds.), Kluwer Academic Publishers, Dordrecht, pp. 147174.
Navaz, S. and Ansari, U. (2012), “An evolutionary algorithm in grid scheduling by multiobjective optimization using variants of NSGA”, International Journal of Scientific and Research Publications, Vol. 2 No. 9, pp. 15.
Negahban, A. and Smith, J.S. (2014), ” “Simulation for manufacturing system design and operation: Literature review and analysis”, Journal of Manufacturing Systems, Vol. 33 No. 2, pp. 241261.
Ng, A.H.C., Bernedixen, J. and Pehrsson, L. (2014), “What does multiobjective optimization have to do with bottleneck improvement of production systems?”, in Proceedings of the Swedish Production Symposium, SPS’14, Swedish Production Academy, Gothenburg.
Ng, A.H.C., Bernedixen, J. and Syberfeldt, A. (2012), “A comparative study of production control mechanisms using simulationbased multiobjective optimization”, International Journal of Production Research, Vol. 50 No. 2, pp. 359377.
Ng, A.H.C., Svensson, J. and Urenda, M. (2008), “MultiObjective simulation optimization for production systems design using FACTS analyzer”, in Proceedings of the 2’nd Swedish Production Symposium (SPS’08), Stockholm, Sweden, 23 November 2008.
Ng A.H.C., Moris, M.U., Svensson, J., Skoogh, A. and Johansson, B. (2007), “FACTS analyzer: an innovative tool for factory conceptual design using simulation”, in Proceedings of the 1st Swedish Production Symposium (SPS07), Gothenburg, 2830 august 2007, pp. 18.
Pehrsson, L., Ng, A.H.C. and Bernedixen, J. (2011), “Multiobjective production system optimisation including investment and running costs”, in Evolutionary MultiObjective Optimization in Product Design and Manufacturing, L. Wang, A. Ng, K. Deb (Eds), Springer, pp. 431454.
Pehrsson, L., Ng, A.H.C. and Stockton, D.J. (2013), “Industrial cost modelling for multiobjective optimisation of production systems”, International Journal of Computer and Industrial Engineering, Vol. 66 No. 4, pp. 10361048.
Pehrsson, L., Ng, A.H.C. and Bernedixen, J. (2016), “Automatic identification of constraints and improvement actions in production systems using multiobjective optimization and postoptimality analysis”, International Journal of Manufacturing Systems, Vol. 39, pp. 2437.
Rosen, S.L., Harmonosky, C.M. and Traband, C.M. (2007), “A simulation optimization method that considers uncertainty and multiple performance measures”, European Journal of Operational Research, Vol. 181 No. 1, pp. 315330.
Rother, M. (2009), Toyota Kata: Managing People for Improvement, Adaptiveness and Superior Results, McGrawHill.
Siegmund, F., Ng, A.H.C. and Deb, K. (2013), “A comparative study of dynamic resampling strategies for guided evolutionary multiobjective optimization”, in Proceedings of the Congress on Evolutionary Computation, Cancún, México, pp. 18261835.
Siegmund, F., Ng, A.H.C. and Deb, K. (2015), “Hybrid dynamic resampling for guided evolutionary multiobjective optimization”, in Proceedings of the 8th International Conference on Evolutionary MultiCriterion Optimization, Guimarães, Portugal. Lecture Notes in Computer Science, pp. 366380, Vol. 9018.
Siegmund, F., Ng, A.H.C. and Deb, K. (2016), “A ranking and selection strategy for preferencebased evolutionary multiobjective optimization of variablenoise problems”, IEEE Congress on Evolutionary Computation, CEC 2016, Vancouver, Canada, July 2429.
Siegmund, F., Bernedixen, J., Pehrsson, L., Ng, A.H.C. and Deb, K. (2012), “Reference pointbased evolutionary multiobjective optimization for industrial systems simulation”, in Proceedings of the Winter Simulation Conference 2012, Berlin.
Simon, D. (2013), Evolutionary Optimization Algorithms: Biologically Inpsired and PopulationBased Approaches to Computer Intelligence, John Wiley and Sons, NJ pp. 557558.
Syberfeldt, A., Ng, A.H.C., John, R.I. and Moore, P. (2010), “Evolutionary optimisation of noisy multiobjective problems using confidencebased dynamic resampling”, European Journal of Operation Research, Vol. 204 No. 3, pp. 533544.
Rajabalipour, C.H., Islam, M.N. and Desa, M.I. (2014), “A polarbased guided multiobjective evolutionary algorithm to search for optimal solutions interested by decisionmakers in a logistics network design problem”, Journal of Intelligent Manufacturing, Vol. 25 No. 4, pp. 699726.
Tempelmeier, H. (2003), “Practical considerations in the optimization of flow production systems”, International Journal of Production Research, Vol. 41 No. 1, pp. 149170.
Wang, Q. and Chatwin, C.R. (2005), “Key issues and developments in modelling and simulationbased methodologies for manufacturing systems analysis, design and performance evaluation”, The International Journal of Advanced Manufacturing Technology, Vol. 25 Nos 11/12, pp. 12541265.
Wang, X. and Wang, X. (2014), “Reference pointbased evolutionary multiobjective optimization for reversible logic circuit synthesis”, in The 7th International Conference on Biomedical Engineering and Informatics (BMEI), Dalian, China, 2014, pp. 955959.
Verma, A. and Kaushal, S. (2015), “Costtime efficient scheduling plan for executing workflows in the cloud”, Journal of Grid Computing, Vol. 13 No. 4, pp. 495506.
Wierzbicki, A.P. (1980), “The use of reference objectives in multiobjective optimization”, in: G. Fandel and T. Gal, (eds), Multiple Criteria Decision Making Theory and Applications, SpringerVerlag, Berlin, pp. 468486.
Wierzbicki, A.P. (1982), “A mathematical basis for satisficing decision making”, Mathematical Modelling, Vol. 3 No. 5, pp. 391405.
Wierzbicki, A.P. (1986), “On the completeness and constructiveness of parametric characterizations to vector optimization problems”, OperationsResearchSpektrum, Vol. 8 No. 2, pp. 7387.
Acknowledgements
This work was partially financed by the Knowledge Foundation (KKS), Sweden, through the BlixtSim project (20112014). The authors gratefully acknowledge their provision of research funding and thank the industrial partners, Volvo Car Corporation and AB Volvo, for their collaborative supports during the project.