An evaluation of three DoE-guided meta-heuristic-based solution methods for a three-echelon sustainable distribution network

This article evaluates the efficiency of three meta-heuristic optimiser (viz. MOGA-II, MOPSO and NSGA-II)-based solution methods for designing a sustainable three-echelon distribution network. The distribution network employs a bi-objective location-routing model. Due to the mathematically NP-hard nature of the model a multi-disciplinary optimisation commercial platform, modeFRONTIER®, is adopted to utilise the solution methods. The proposed Design of Experiment (DoE)-guided solution methods are of two phased that solve the NP-hard model to attain minimal total costs and total CO2 emission from transportation. Convergence of the optimisers are tested and compared. Ranking of the realistic results are examined using Pareto frontiers and the Technique for Order Preference by Similarity to Ideal Solution approach, followed by determination of the optimal transportation routes. A case of an Irish dairy processing industry’s three-echelon logistics network is considered to validate the solution methods. The results obtained through the proposed methods provide information on open/closed distribution centres (DCs), vehicle routing patterns connecting plants to DCs, open DCs to retailers and retailers to retailers, and number of trucks required in each route to transport the products. It is found that the DoE-guided NSGA-II optimiser based solution is more efficient when compared with the DoE-guided MOGA-II and MOPSO optimiser based solution methods in solving the bi-objective NP-hard three-echelon sustainable model. This efficient solution method enable managers to structure the physical distribution network on the demand side of a logistics network, minimising total cost and total CO2 emission from transportation while satisfying all operational constraints.


Introduction
Location routing models are typical logistics problems which are generally complex NPhard combinatorial optimisation problems (Karp 1972;Nagy and Salhi 2007;Yu et al. 2010). The three-echelon multi-objective location routing problems are NP-hard (Karp 1972;Perl 1983;Perl and Daskin 1985;Daskin et al. 2005;Nagy and Salhi 2007;Hassanzadeh et al. 2009;Yu et al. 2010;Hamidi et al. 2012;Validi et al. 2020). This is because they combine two different conflicting-in-nature problems, viz. facility location problems, and vehicle routing problems (Karp 1972;Daskin et al. 2005). The computational effort required for the solution of such problems grows exponentially with increasing problem size (Erdoğan and Miller-Hooks 2012).
Meta-heuristics are often used to solve these kind of models as a unique solution does not exist for such models. The solution method requires a division of the problem into inter-linked components, viz. facility location, allocation of consumers to facilities, and vehicle routing (Hassanzadeh et al. 2009). This leads to the introduction of multiple phase solutions with specific algorithms (Perl 1983;Wu et al. 2002).
The sustainable three-echelon bi-objective AHP-integrated 0-1 mixed integer locationrouting model reported in Validi et al. (2020) (see "Appendix A") is NP-hard in nature. The model facilitates allocation of distribution centres (DCs) to plants and retailers to DCs. It further facilitates routing of vehicles from plants to DCs, DCs to retailers and retailers to retailers. The model identifies open and closed facilities and the optimised routing pattern throughout the logistics network while minimising CO 2 emissions from transportation and total transportation costs considering disparate operational constraints.
In this article, three solution methods are proposed and compared employing two Design of Experiment (DoE)-guided genetic algorithm (GA)-based optimisers, viz. MOGA-II (Multi-Objective Genetic Algorithm of kind II) (Quagliarella et al. 1998) and NSGA-II (Non-dominated Sorted Genetic Algorithm-II) (Deb et al. 2000), and a particle swarm (PS)-based optimiser, viz. MOPSO (Multiple Objective Particle Swarm Optimisation) (Alvarez-Benitez et al. 2005) in an effort to find the most efficient meta-heuristics based solution method to solve the sustainable three-echelon bi-objective location-routing problem (BO-LRP) introduced in Validi et al. (2020).

Research contribution
This research contributes to the body of knowledge by providing two-phase solution methods employing the three prominent optimisers for an integrated bi-objective three-echelon sustainable distribution model and comparing the results obtained from them. Additionally, it investigates the efficacy of the three solution methods. In these solution methods the DoE is coupled to the meta-heuristic optimisers to ensure that optimal feasible results are obtained. The DoE connects the main three-echelon problem to the optimisers by generating an initial population sets for the optimisers, filtering out based on poor design criteria.
To the best of the authors' knowledge, an evaluation of the DoE-guided meta-heuristic based two-phase solution methods for multi-objective three-echelon sustainable distribution model has not been reported in the literature. To bridge this knowledge gap, this article considers the bi-objective three-echelon sustainable distribution model of Validi et al. (2020) and contributes to the current body of the literature by: (a) proposing three DoE-guided two-phase meta-heuristics based solution methods, using the modeFRONTIER ® commercial platform, for solving the above mentioned NP-hard model, and (b) comparing the efficacy of the DoE-guided two-phase disparate meta-heuristic based solution methods to obtain the best stable solution space.
The decision-makers' (DMs') priorities are reflected in comparing and ranking the set of Pareto efficient realistic optimum results by integrating a compensatory multi-criteria decision analysis method, viz. Technique for Order Preference by Similarity to Ideal Solution (TOPSIS).
The remainder of this paper is structured as follows. Section 2 provides a brief survey of literature on sustainable distribution network model and meta-heuristic based solution methods. Section 3 outlines the sustainable three-echelon logistics model, problem statement and relevant data sets. Section 4 explains the DoE-guided two-phase solution methods. Section 5 provides the results and a detailed analysis of these results. Section 6 provides the structure of the optimal routes obtained through the three methods indicating logistical scenarios. Finally, Sect. 7 concludes the paper with an implication to future research.

Literature review
The extant literature on sustainable supply chain (SC) distribution network design is categorised into models and their solution approaches using meta-heuristics. The following sections briefly review articles on these two categories.

Sustainable distribution network models
A distribution network on the demand side of a logistics network facilitates manufacturers in delivering products to their consumers through a number of logistics intermediaries. Distribution networks tend to be dominated by transportation means that primarily use fossil fuel. As a result, a significant focus for sustainable distribution network studies has been on minimising fuel consumption and reducing carbon footprint. The design of the transportation network plays a pivotal role in achieving this goal. The extant literature introduces different models to address this need. For example, two and three echelon sustainable distribution network models are introduced in Validi et al. (2014aValidi et al. ( , b, 2015Validi et al. ( , 2020. Wang et al. (2011) report a multi-objective sustainable SC network model that minimises total cost and CO 2 emission. Daryanto et al. (2019) propose a three-echelon SC model that reduces carbon emission through optimising the number of deliveries and delivery size. Chaabane et al. (2012) propose a sustainable SC using a mixed-integer linear programming that provides some underpinning of optimal SC strategies. Zhang et al. (2017) report a multi-criteria optimisation model and propose a bio-inspired algorithm to design a sustainable supply chain (SC) network. Varsei and Polyakovskiy (2017) propose a multi-objective mixed-integer programing model to design a sustainable SC. Some comprehensive review articles are reported in this fast growing subject area. For example, a review on sustainable SC network design models, solution approaches and applications is found in Eskandarpour et al. (2015). Brandenburg and Rebs (2015) present a literature review and provide seven modelling guidelines to facilitate model-based SC research.

Meta-heuristic based solution methods
A review of the available literature reveals that the three echelon LRPs are computationally NP-hard. LRPs combine facility location and vehicle routing problems. These are two conflicting problems that contribute to the hardness of NP-hard problems. Other characteristics, such as the number of objective functions, the integer variables, problem size and inclusion of substantially factors add significantly to the complexity and difficulty of finding a solution space. Due to the lack of unique solutions to these LPRs, heuristics, metaheurists and hyper-heuristics, in multi-phase solution methods, are used to find a feasible solution space. Multi-phase solution methods usually break down the problem into individual components forming two or three connected problems to be solved (Karp 1972;Golden and Skiscim 1986;Tuzun and Burke 1999;Prins et al. 2007;Bräysy et al. 2009;Prins et al. 2009;Caballero et al. 2007;Hassanzadeh et al. 2009;Daskin et al. 2005;Yu et al. 2010;Validi et al. 2020;Dai et al. 2019).
One of the major challenges in designing an evolutionary-based meta-heuristic solution method is the details of its design, in particular parameter setting (Hoos 2011;Barbosa and Senne 2017). A range of methods from more conventional approaches with low-mutation rates to uniform crossover choices and experimental designs have been introduced and adopted in practice (Eiben and Smit 2011). Existing parameter setting methods start with generating different settings of parameters and then evaluate them by testing and comparing their performances (Castelli et al. 2012). Huang et al. (2020) categorise the existing parameter setting methods into three main groups: (i) simple generate-evaluate methods: consisting of a generation step for configuring parameters followed by an evaluation phase (ii) iterative generate-evaluate methods: which create the parameters and improve them gradually in a repeating process, and (iii) high-level generate-evaluate methods: consisting of a sophisticated group of methods that quickly generate an 'elite' group of parameters and then carefully select the best of them. The main advantage of the latter methods is in adopting a search algorithm instead of applying a random method in the generating phase and then instead of evaluating, selecting the best set of parameters amongst the generated set. Design of Experiments (DoE) belongs to this third group of parameter setting methods.

Research gaps
A review of the literature reveals that a growing number of attempts have been made to formulate two or three echelon logistics/distribution network models. Consequently, a variety 1 3 of meta-heuristics based solution methods have been developed and adopted for solving the existing and new NP-hard models. Yet, to date the extant literature has not reported: (i) a DoE-guided two-phase meta-heuristic solution methods capable of solving an integrated bi-objective AHP-integrated sustainable three-echelon distribution model, and (ii) a comparative study on the efficacy of DoE-guided two-phase meta-heuristic solution methods.

The model and preliminaries
The sustainable three-echelon bi-objective AHP-integrated 0-1 mixed integer sustainable distribution network model introduced in Validi et al. (2020), which serves the demandside of the logistics network, is considered. The summarised model is presented in "Appendix A". A summary of the main details of the model are presented in Fig. 1.
The model allocates DCs to the plants and retailers to DCs. The model seeks to identify open and closed facilities and the optimised routing pattern throughout the transportation network while minimising CO 2 emissions from transportation and total transportation costs. Two sustainable components are present in the model: (i) a sustainable objective function and (ii) an AHP-integrated sustainable constraint. The model considers 10 operational constraints in total. Integer decision variables define open/close plants and DCs. The 1. Each demand node on one route 2. Limits the length of each route 3. Assigns each route to one facility 4. Any route entering a node must exit the same node 5. A route can operate out of only one facility 6. Defines the flow into a facility from the supply points (in terms of demand) 7. Restricts throughput at each facility to the maximum allowed at that site and links the flow variables and facility location variables 8. If path leaves customer node and facility , then customer i must be assigned to facility j 9. AHP-integrated constraint (sustainability constraint), considering the DMs' priorities Integer constraints Non-negativity constraints Fig. 1 The details of the model 1 3 decision variables also define all real feasible routes from plants to DCs, DCs to retailers and routes connecting retailers to retailers. The parameters in the model represent the amount of CO 2 emission in each real feasible route from plants to DCs, DCs to retailers, and among retailers. The costs of serving DCs from plants, retailers from DCs and retailers from other retailers are defined by parameters in the model.

Problem statement and data set
A dairy processing industry's three-echelon logistics network, operating in the east of Ireland, is considered to test the model using the proposed solution methods. This logistics network has two milk processing plants, six DCs and supplies twenty-two retailers ( Table 20). The geographical locations of the facilities and retailers enable the model to measure the traversed distances, costs of transportation and the associated CO 2 emission. The fixed costs associated with the processing plants and DCs are illustrated in Tables 21  and 22. These costs refer to the total fixed costs of operating each facility in a cycle time of 2-3 days. The variable cost at each DC is provided in Table 22. The variable cost is assumed to be the cost that each plant or DC should bear in order to meet the demand of the consumer. Variable cost at each plant is illustrated in Table 21. Table 23 depicts the capacity of each DC. Average demand at each retailer is estimated as 2/3 of the total population at the geographical location of the retailer. Heavy duty vehicles are considered for deliveries throughout the network. These vehicles are refrigerated. In Ireland, like other countries, different speed limits are set depending on the nature of the roads. In order to avoid cumbersome computational process an average speed on the road is assumed. Table 24 represents the stipulated speed limits according to the "Road Traffic Act 2004" and the corresponding average speeds are assumed for the purpose of the model.
Diesel fuelled fully loaded vehicles are considered. The fuel efficiency is considered as 0.35 based on the report of the UK's Department of Energy and Climate Change (2010) and Nylund and Erkkilä (2005). Guidelines to DEFRA's (2008) greenhouse gas conversion factors aid in calculating the CO 2 emission from the diesel vehicles.
The cost of serving each of the routes from the plants to the DCs and from the DCs to the retailers is the sum of fuel costs and driver's wage. This cost considers the fuel cost and driver's wage. The average cost of diesel in Ireland as €1.53/l and the average wage of a heavy-duty vehicle driver in Ireland as €11.50/h. The speed of the vehicle affects the cost of serving the route. Each route includes a combination of different types of roads, e.g., motorways, national primary and secondary routes, regional routes and local roads. Therefore, average speed in each path is considered.
The CO 2 emission from a diesel vehicle (in kg) is determined using the information presented in Table 25. Table 1 shows the CO 2 emission from the burnt fuel for travelling between the plants and the DCs using the designated routes and corresponding cost of serving each route.
In a similar manner, the CO 2 emission during transportation of the products between the DCs and retailers and corresponding costs of serving each route between the DCs and retailers are computed. The CO 2 emission during transportation of the products among the retailers is also calculated. In addition, the cost of serving each route among retailers are determined.
The AHP-integrated constraint of the model considers the DMs' priorities in terms of vehicle types. Three different types of vehicle are considered for transportation of the 1 3 products for illustrative purpose. The vehicles have different levels of CO 2 emission and costs. Table 26 illustrates the truck types based on these two attributes. The DMs are asked to compare the vehicles using Saaty's nine-point scale (Saaty 1977(Saaty , 1994. The weight matrix (Table 27) is determined using the pair-wise comparison matrix of AHP. The boundary values (i.e. lowest and highest) of CO 2 emission are found considering the average round off values of the CO 2 emission during transportation and costs of serving the routes (Table 28).

Solution approach
The three echelon bi-objective combinatorial optimisation location routing model is NPhard. A two-phased solution method is designed to find the feasible solution area for the model. This solution method divides the main model into two interlinked problems considering different elements of the model. modeFRONTIER ® have been selected as the main solution platform for executing the solution method. modeFRONTIER ® is a multi-disciplinary and multi-objective software developed by ESTECO SpA and is chosen to solve the LRP model. modeFRONTIER ® is a multi-criteria decision-analysis optimisation and design environment. This platform consists of a variety of evolutionary based algorithms and a Design of Experiment (DoE)-guided solution approach (ESTECO 2020). DoE connects the evolutionary-based meta-heuristics to the optimisation model by generating an initial population. A variety of distributions and designs are considered in the main DoE search algorithm to generate sets of parameters and then the best set of parameters are selected as the initial set of parameters. The details of the DoE in each of the proposed solution methods are explained in Sect. 4.1.

Proposed solution methods
The sustainable three-echelon bi-objective distribution network model is solved in two phases, viz. Phase-I and Phase-II. Based on the outcome of Phase-I, Phase-II is modified and then solved. The final results are obtained when both phases are complete. The process of solving the model and analysing results in the two-phase method is illustrated in Fig. 2.  Fig. 2 The scope and the relationship between the two phases of the solution methods by connecting the optimisation model to the selected optimisers. DoE considers a combination of designs and distributions to search for the best set of parameters. Design of experiment sequence, Random design, Sobol sequence, Uniform Latin hypercube, Incremental space filler, and Constraint satisfaction design are the search algorithms used. A population of 50 have been selected by the DoE process to enhance the performance of the meta-heuristics and consequently improve the optimisation process. The details of the designs are presented in Table 2.
Amongst many different combinations of designs tested, this combination of designs proved to be the most efficient for the three-echelon sustainable distribution model. The same initial population has been used for all three optimisers.
Two GA-based optimisers, viz. MOGA-II and NSGA-II, and one PS-based optimiser, viz. MOPSO, are considered in solving the sustainable three-echelon distribution network problem. The chosen optimisers have disparate requirements and distinctive specifications. Therefore, they have different set up details (Table 3).
Based on experimentation and analysis of convergence characteristics and the number of realistic results for each of the three optimisers, the final number of generations has been set to 250 for the MOGA-II and MOSPSO optimisers and to 50 for the NSGA-II optimiser. 250 generations for the MOGA-II and MOSPSO optimisers, generates a larger feasible solution area with clear convergence within this limit. Convergence for the NSGA-II optimiser has been seen to occur much faster and through experimentation NSGA-II was found to converge in 50 generations. After the 50th generation the optimiser doesn't improve the results any further. For this reason, the number of generations in NSGA-II is set to 50 generations, as to run to 250 generations was surplus to requirements.

Results and analyses
The three optimisers have been tested with the same data set using modeFRONTIER ® , a multi-disciplinary optimisation platform. Obtained results from Phase I and Phase II are further analysed and final outcomes of the model are presented. The analysis process is illustrated in Fig. 3. modeFRONTIER ® is the main platform for implementing the solution methods and performing Analysis of Variance (ANOVA). The main advantages of the embedded ANOVA in modeFRONTIER ® are: (i) its direct access to the feasible solution space and the results obtained from the three optimisers simultaneously, and (ii) the test results are part of the DoE that guide the solution methods. Using the optimisation platform, the initial results obtained through the solution methods are refined, and a set of realistic results are identified. In the process of refining results, unrealistic and identical results are eliminated. Unrealistic results are those that don't show any improvement when compared to the previous generations and those that are operationally unrealistic. This is followed by a performance study on each optimiser. Next, ANOVA test is performed to compare means of multiple groups for the optimised data. A set of results are selected for further analysis. Pareto

Refining results
Eliminating unrealistic and identical results to find the set of realistic results

ANOVA test
Comparing the means of two or more groups of realistic results

Selected results
Selecting top 20 result, within each set of realistic results, for further analysis

Pareto efficiency
Checking the Pareto efficiency of selected results

TOPSIS
Considering the priorities of DMs to rank the selected results

Outcomes Geographical maps
Geographically mapping the ranked selected results

Scenario Analysis
Analysing the effect of opening closed routes on the value of objective functions Fig. 3 Analysis process efficiency is examined on the selected results. Finally, the selected results are ranked using TOPSIS (Hwang and Yoon 1981). TOPSIS, a multi-criteria decision-making method, has been considered for ranking data due to its compensatory nature. This compensatory method defines weight and normalises scores for each criterion by selecting an alternative ('results' here) with the shortest distance from the ideal and the farthest distance the non-ideal result (Zavadskas et al. 2006;Yoon and Hwang 1995). When compared with non-compensatory methods, the compensatory capability of TOPSIS provides realistically well-balanced trade-offs between criteria (Asgharpour 1998), which work very well in ranking the selected results of this bi-objective three-echelon LRP model considering the DMs' priorities.

Phase-I results
A statistical summary on the feasible results offered by the three optimisers, are presented in Table 4. These statistical summaries assist in analysing: (a) the total results table which consists of all feasible results obtained by each optimiser, (b) the refined realistic results table which consists of realistic and non-identical results, and (c) the selected results table which consists of selected results from the three lowest sets of feasible results from the Bubble plots.
The minimum value for both CO 2 emissions and costs is obtained by NSGA-II in 50 generations (Table 4). This does not necessarily yield a minimum value for the total costs. Figures 4, 5 and 6 illustrate the feasible real solution space for MOGA-II, MOPSO, and NSGA-II respectively. Each bubble in these figures represents a realistic result. In these plots, the colour and diameter of bubbles are representatives for values of the objective functions where the bubble size represents costs and the colour represents CO 2 emission. Colours range from dark blue to red, and lowest to highest value of the objective function respectively. The diameters of the bubbles range from small to large and lowest to highest

Performance study of the optimisers in Phase-I
The performance of the optimisers regarding their convergence is investigated through convergence plots. Figures 7, 8 and 9 present the convergence plots for MOGA-II, MOPSO and NSGA-II respectively with regard to the objective functions. All optimisers converge in a steady manner. The convergence of all of these evolutionary optimisers is studied and compared based on their final results. The initial population set is generated by DoE for these evolutionary optimisers. This DoE-guided

ANOVA for Phase-I
One-way ANOVA is performed for both the total CO 2 emissions and total costs of transportation to compare the means of multiple groups for the optimised data. Tables 5 and 6   illustrate the ANOVA results for the three optimisers with respect to CO 2 emissions and costs respectively. The most important assumption requested by ANOVA is that the standard deviations within each group are the same. It is found that Hartley and Bartlett's tests are both true. These two statistical tests verify that the standard deviations within each group is the same, therefore the most important assumption requested by ANOVA is valid. It is noted that as the F-ratio increases, the p value decreases.

Selection of Phase-I results
Within the feasible solution space, unrealistic and identical results are eliminated to identify the refined realistic results. A set of 20 results from the refined realistic results are selected for further analysis. The first three shades of blue in Figs. 4, 5 and 6 represent the lowest three sets of values for the objective functions. Therefore, the selected results, with respect to the objective function values, are chosen from these three sets of results. A low value for one objective function may not necessarily yield a low value for the other objective function and vice versa. The realistic results are shown within the feasible solution space in Figs. 10, 11 and 12 for MOGA-II, MOSPO, and NSGA-II respectively.
The values for objective functions offered by optimisers are presented with reference to the result ID of results. The selected results are highlighted in bright green. In the selection of the three-echelon BO-LRP results, DMs' priorities are considered. Therefore, in each set of the selected results at least one result represents extreme decision-making events.
Selected results are ranked in order to prioritise the optimised results. The ranking process is explained in the following section.

Ranking selected results from Phase-I
The selected sets of 20 results from MOGA-II, MOPSO and NSGA-II are ranked using TOPSIS. Ranking is performed to prioritise the best results for different DMs. Nine weight matrices are defined for TOPSIS. Each matrix represents a type of DM. The DMs select relevant weights from Saaty's nine-point scale (1977,1994). DMs have been asked to prioritise options of transportation based on two criteria viz., CO 2 emission and cost, thereby introducing flexibility in the decision-making process. The decision-makers' priorities for CO 2 emission and cost of these matrices are as follows: w 1 = (0.1, 0.9), w 2 = (0.2, 0.8), w 3 = (0.3, 0.7), w 4 = (0.4, 0.6), w 5 = (0.5, 0.5), w 6 = (0.6, 0.4), w 7 = (0.7, 0.3), w 8 = (0.8, 0.2) and w 9 = (0.9, 0.1). We consider the weight matrix of w 5 as it has equal importance on CO 2 emission and cost. The first three ranked results by TOPSIS using the weight matrix w 5 for the optimisers are presented in Table 7.
As shown in Table 7, the best result with regard to objective function values is offered by NSGA-II. MOPSO is in the second best and MOGA-II the third optimiser.

Pareto efficiency for Phase-I
Pareto efficiency is examined to evaluate the performance of the optimisers. The 20 selected results obtained from MOGA-II, MOPSO and NSGA-II optimisers are separately examined with regard to their Pareto efficiency and presented in Figs. 13, 14 and 15.
As is evident from Figs. 13, 14 and 15 the selected results from MOGA-II and NSGA-II optimisers follow the Pareto optimality and are strongly efficient. MOSPO shows a reasonably strong Pareto efficiency on the selected results. In the selection process of results, extreme decision-making events are considered as well. Therefore, a small number of results representing these events exist in each selected results table. These results do not affect the Pareto efficiency of the selected results as they do not

Phase-II results
The first three ranked results obtained in Phase-I using weight matrix w 5 (Table 7) is used as a basis for implementing Phase-II. The three optimisers MOGA-II, MOPSO, and NSGA-II used in Phase-I are used to implement the model in Phase-II. In order to analyse the results obtained from Phase-II: (i) results are refined and a set of realistic results are identified, (ii) a performance study is conducted on each optimiser, (iii) Analysis of Variance (ANOVA) test is performed to compare means of multiple groups for the optimised

Refinement of Phase-II results in feasible real solution space
All results from the MOGA-II, NSGA-II, and MOPSO optimisers in Phase-II are feasible and real. A statistical summary of these results is illustrated in Table 8. These results assist in analysing: (a) the total results table consisting of all feasible real results obtained from the optimisers, (b) the refined realistic results table consisting of realistic and non-identical results, and (c) the selected results table consisting of selected results from the three lowest sets of feasible results. As presented in Table 8, the minimum value for both CO 2 emissions and costs is obtained by NSGA-II in 50 generations. This does not necessarily yield a minimum value for the total costs.
MOGA-II, NSGA-II and MOPSO generate real feasible solutions spaces in Phase-II.   green. Bubble size represents costs while colour represents the CO 2 emission. Colours range from dark blue to red, lowest to highest value of the objective function respectively. The diameter of bubbles ranges from small to large, lowest to highest value of the objective function respectively.

Performance study on optimisers in Phase-II
The performance of the optimisers regarding their convergence is studied comparatively through plots presented in Figs. 19, 20 and 21 illustrate the convergence plots for MOGA-II, MOPSO and NSGA-II, with respect to the objective functions. It is evident from plots 'a ' and 'b' of Figs. 19,20 and 21 that all the optimisers are converging in a steady manner similar to Phase-I. MOGA-II and MOPSO optimisers converge in 250 generations while NSGA-II in 50 generations.

ANOVA for Phase-II
The means of multiple groups of the optimised data are compared through one-way ANOVA with respect to the objective functions. These comparisons are presented in Tables 9 and 10.

3
The most important assumption in ANOVA is that the standard deviations within each group are the same. Hartley and Bartlett's tests are true in this case. These two statistical tests verify that the standard deviations within each group are the same, which means that the assumption of ANOVA is valid. It is also noted that as the F-ratio increases, the p-value decreases. The p-values of the ANOVA table for CO 2 emission and costs (Tables 9 and 10) in MOGA-II optimiser are zero. This suggests that there are significant differences between the groups. At least one sample mean is significantly different from the other sample means.

Selection of results in Phase-II
After eliminating the un-realistic and identical results from the feasible solutions space, the refined realistic results are considered for further analysis. The analysis process performed in Phase-I to select results is adapted here. Figures 22, 23 and 24 exhibit the bubble plots on selected results for MOGA-II, MOSPO and NSGA-II respectively.

Ranking selected results from Phase-II
Similar to the process of ranking results in Phase-I, selected results from Phase-II are ranked. The first three ranked results by TOPSIS using weight matrix w 5 for all three optimisers are presented in Table 11.

Pareto efficiency in Phase-II
The 20 selected results obtained from MOGA-II, MOPSO, and NSGA-II optimisers are separately examined with regard to their Pareto efficiency. Pareto frontiers are presented in Figs. 25, 26 and 27.

Optimal three-echelon distribution network
The final result from the three-echelon BO-LRP is the combined outcome from Phase-I and Phase-II. These results concern facility and vehicle-routing decisions on the demand side of the three-echelon logistics network in the east of Ireland.

Final results
The sustainable three-echelon bi-objective distribution network results offer an optimal facility location and vehicle-routing decision. With regard to facility location decisions, optimum open and closed DCs are offered in each result. With respect to vehicle-routing decisions: (i) optimum routes connecting plants to DCs, DCs to retailers and connection in between retailers, (ii) type of vehicle and (iii) the number of heavy good vehicles (HGVs) required for transporting products in each route are offered. Table 12 illustrates the 'facility location', 'type of truck' and the 'routing pattern' for the disparate optimisers. NSGA-II offers the best results for the three-echelon BO-LRP in two inter-linked phases (Table 12). NSGA-II converges and offers the best results in 50 generations while MOGA-II and MOPSO converge in 250 generations.
With regard to vehicle-routing decisions, the quantity transported in each open route and the numbers of HGVs required for transporting the defined load of product are presented in Tables 13, 14 and 15 with respect to the MOGA-II, MOPSO and NAGA-II first ranked  2,3,4,5,6,7,11,13,18,21 DC 5 serving retailers: 1,9,10,12,16,17,19,20,22  results. The performance of the optimisers in Phase-I and Phase-II shows that NSGA-II is very efficient in solving the bi-objective NP-hard three-echelon sustainable LRP. A schematic presentation of the final results obtained from MOGA-II, MOPSO and NSGA-II are illustrated in Fig. 28a-c. 1 3

Scenario analysis with final results
Disparate scenario analysis on the first ranked results obtained from the three optimisers is presented in Tables 16, 17 and 18.  Tables 16, 17 and 18 are examples of the guidance available to DMs for locating the feasible and realistic optimal distribution routes, considering the trade-offs with respect to the objective functions, if a closed route is forcibly opened.

Conclusions
Three two-phase solution methods are proposed and employed as the sustainable three-echelon bi-objective distribution network model is constrained to a considerable degree and impossible to solve in single phase. Three-echelon LRPs deal with two types of decisions in designing the logistics networks, viz. facility location and vehicle routing decisions.  The open DCs and the routes connecting them to the served retailers from Phase-I are included in Phase-II. The objective function and constraint elements relating to open DCs and routes connecting them to retailers from Phase-I are considered in Phase-II. This is the link for connecting the two phases in order to reach to final optimal result(s) for the three-echelon BO-LRP. Two GA-based and one PS-based meta-heuristic optimisers have performed efficiently in reaching an optimum result for the BO-LRP model.
The solution methods to the three-echelon sustainable distribution network model are DoE-guided. The optimisers have been set with the same parameter values as possible to compare their performances in solving the three-echelon BO-LRP. An experimentally generated initial population for the optimisers by DoE is proven to be efficient in offering the optimum set of results for the optimisation model. The convergence of the optimisers is examined with regard to the objective functions for each optimiser. One way ANOVA is performed for all results obtained from each optimiser for the objective functions. The Hartley and Bartlett's statistics tests verify that the standard deviations within each group is the same. Based on the statistical selection criteria a set of results are identified and consequently ranked by TOPSIS. Pareto efficiency of the selected results is studied. Selected results from all the three optimisers are proved to be strongly Pareto efficient.

3
The final results are obtained through the execution of two interconnected phases. The final results of the model consist of (i) information on the open/close DCs, (ii) the vehicle routing patterns connecting the plants to DCs, (iii) the vehicle routing patterns connecting the open DCs to the retailers, (iv) the vehicle routing patterns connecting the retailers to retailers, and (v) the number of trucks required in each route to transport the products. By obtaining the above optimal setting, the physical distribution network on the demand side of the SC network can be structured with the main aim of minimising the total cost and minimising the total CO 2 emitted from transportation while satisfying the operational constraints.
Scenario analysis is performed on the final results. This scenario analysis provides guidance to DMs when a closed route is forcibly opened. Various scenarios depict the amount of CO 2 emitted and the total costs from a closed route if forced to be open. This provides support for SC network resilience.
Future research may focus on the use of meta-heuristics and assessment of their efficiencies in solving an NP-hard model with uncertainties using the DoE-guided two-phase solution method. Further, solution of the NP-hard sustainable distribution network model using hyper heuristic is another area of research to be focused in future.  --

Acknowledgements
The authors sincerely convey thanks to the Guest Editors of the special issue entitled "New Trends in Multiobjective Optimization Techniques and Applications" and anonymous reviewers for their constructive comments.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.