Sustainable distribution system design: a two-phase DoE-guided meta-heuristic solution approach for a three-echelon bi-objective AHP-integrated location-routing model

This article introduces a sustainable integrated bi-objective location-routing model, its two-phase solution approach and an analysis procedure for the distribution side of three-echelon logistics networks. The mixed-integer programming model captures several real-world factors by introducing an additional objective function and a set of new constraints in themodel that outbound logistics channels find difficult to reconcile. The sustainablemodel minimises CO2 emissions from transportation and total costs incurred in facilities and the transportation channels. Design of Experiment (DoE) is integrated to themeta-heuristic based optimiser to solve the model in two phases. The DoE-guided solution approach enables the optimiser to offer the best stable solution space by taking out solutions with poor design features from the space and refining the feasible solutions using a convergence algorithm thereby selecting the realistic results. Several alternative solution scenarios are obtained by prioritising and ranking the realistic solution sets through a multi-attribute decision analysis tool, Technique for Order Preference by Similarity to Ideal Solution (TOPSIS). The robust model provides the decision maker the ability to take decisions on sustainable open alternative optimal routes. The outcomes of this research provide theoretical and methodological contributions, in terms of integrated bi-objective location-routing model and its two-phase DoE-guidedmeta-heuristic solution approach, for the distribution side of three-echelon logistics networks. B Arijit Bhattacharya A.Bhattacharya@uea.ac.uk; arijit.bhattacharya2005@gmail.com Sahar Validi S.Validi@hud.ac.uk P. J. Byrne pj.byrne@dcu.ie 1 Logistics Transport and Supply Chain Subject Group, School of Business, University of Huddersfield, Queensgate, Huddersfield HD1 3DH, UK 2 Norwich Business School, University of East Anglia, Norwich Research Park, Norwich NR4 7TJ, UK 3 Dublin City University Business School, Dublin City University, Glasnevin, Dublin 9, Ireland


Introduction
Sustainable logistics operations and their effect on carbon emissions and costs are a current strategic challenge in modern supply chains (SCs) (Du et al. 2017;McKinnon et al. 2015;Brandenburg and Rebs 2015;Kumar et al. 2012;Srivastava 2007). Logistics activities have been identified as one of the most significant sources of air pollution and greenhouse gas emissions (Sheu and Li 2014;Wang et al. 2011). Efficient and effective design of the outbound logistic system is one of the critical success factors for sustainable movement of products to multiple retailers and consumers (Lopes et al. 2008). Therefore, sustainable location-routing problem (LRP) is of major concern in logistics networks design (Srivastava 2007;Seuring and Müller 2008).
A sustainable transportation strategy within a logistics network promotes an approach that seeks to achieve mutually reinforcing benefits for the economy, the environment and society (Ilbery and Maye 2005). The Paris Agreement (United Nations Treaty Collection 2017) sets out a global action plan to reduce member states' carbon emission through a legally binding global climate deal and encourages businesses to significantly reduce carbon emissions from their operations. Transportation activities on the demand side of the SC via roadways leave harmful effects on human health and the environment. Therefore, the transportation routes should be planned so as to minimise CO 2 emissions and costs resulted from transportation activities.
Driven by the current laws and regulations and competitive business opportunities on the sustainable transportation mechanism, this article addresses three inter-linked aspects of sustainable location-routing. First a sustainable bi-objective three-echelon integrated optimisation model is proposed, next a Design of Experiment (DoE)-guided meta-heuristic-based robust solution approach is provided to solve the computationally NP-hard integrated model, and the Decision-Makers' (DMs') prioritisation and ranking of the realistic solution sets, and various routing scenarios of the three-echelon sustainable location-routing are illustrated.
The model addresses all three echelons on the distribution side of the logistics network responsible for carbon emission. The weighting procedures of Analytic Hierarchy Process (AHP) (Saaty 1994) is infused into the bi-objective 0-1 integer programming in order to provide flexibility in the vehicle selection decision-making process. The computational complexity of the NP-hard model necessitates the solution approach to divide the model into two inter-connected phases. While Phase-I deals with the sustainable transportation mechanism from the processing plants to the multiple distribution centres (DCs) and DCs to multiple retailers, Phase-II uses the outcome of Phase-I to locate the non-dominated Pareto realistic optimal routes among retailers. The DoE-guided Multiple-Objective Genetic Algorithm of kind II (MOGA-II) optimiser enables the model to obtain the best set of realistic solutions. Disparate scenarios of the realistic vehicle routes are captured with alternative possible outcomes from both the phases. The final set of optimal and realistic solutions is then obtained from the synergistic effect resulted from both phases. The best set of realistic sustainable optimal solutions is then obtained.
The remainder of this paper is organised as follows. Section 3 formulates the sustainable three-echelon location-routing model. The two-phase solution approach to the NP-hard model is explained in Sect. 4. Next section implements the model and its two-phase solution approach on a three-echelon dairy processing logistics network based in Ireland. The results and analysis of the implemented model are illustrated in Sect. 6. A critical discussion on the results is also included in this section. Finally, Sect. 7 concludes the paper with an implication of the proposed approach on sustainable location-routing on the demand side of logistics network.

Theoretical background
Literature on conventional LRPs is rich in terms of both methodologies and applications. Prominent reviews can be found in Madsen (1983), Nagy and Salhi (2007), Dekker et al. (2012), Demir et al. (2014), Prodhon and Prins (2014), Eskandarpour et al. (2015), Drexl and Schneider (2015), Vega-Mejía et al. (2017). The purpose of LRP methodologies is the simultaneous determination of the location of facilities and the routes of vehicles for product transportation (Laporte et al. 1988;Nagy and Salhi 2007;Yu et al. 2010). A broad range of LRP applications is found in the literature. A synopsis of these applications is illustrated in Table 1.
The nature of the LRP problems stimulates practitioners to use a range of optimisation techniques. Stochastic (Laporte et al. 1989;Albareda-Sambola et al. 2007) and deterministic (Albareda-Sambola et al. 2005) optimisation models are reported in the literature. For example, a mixed-integer programming technique is reported for minimising the sum of the fixed costs and distribution costs (Diabat and Simchi-Levi 2009) in LRP. Belenguer et al. (2011) propose an exact approach using a 0-1 linear model based on a branch-and-cut algorithm for solving the LRP with capacity constraints on DCs and vehicles. Rath and Gutjahr (2014) report a three-objective warehouse location-routing problem for disaster relief. Karaoglan et al. (2011) propose an exact algorithm based on a branch-and-cut technique. Berger et al. (2007) propose a branch and price algorithm. Integer-linear programming is very common in solving LRP. Applications of integer-linear programming are found in multi-level locationrouting-inventory (Ambrosino and Scutellà 2005), "road-train" routing (Semet 1995) and Eulerian location (Ghiani and Laporte 1999) etc. LRP has been viewed also as a Lagrangian relaxation model (Cappanera et al. 2004) and mixed-integer programming (Alumur and Kara 2007).
The optimisation formulations of efficient LRPs are NP-hard by nature (Nagy and Salhi 2007;Yu et al. 2010). One of the characteristics of these NPhard problems is that the solution grows exponentially with increasing problem size (Erdogan and Miller-Hooks 2012). Therefore, in order to tackle this situation, the application of metaheuristics in LRP is necessary (Golden and Skiscim 1986;Tuzun and Burke 1999;Prins et al. 2007;Bräysy et al. 2009;Prins et al. 2009). Examples on the use of meta-heuristics in LRP are abundant. This includes the use of particle swarm optimisation (Liu et al. 2012), tabu search (Russell et al. 2008), simulated annealing (Stenger et al. 2012), greedy randomised adaptive search procedure (Nguyen et al. 2012), variable neighbourhood search algorithms (Derbel et al. 2011), ant colony optimisation (Ting and Chen 2013) and honey bees mating optimisation . The use of genetic algorithms in LRPs are reported in Zhou and Liu (2007), , Jin et al. (2010), Karaoglan et al. (2012). Some variants of genetic algorithms are also reported in LRP literature (Hwang 2002;Zhou and Liu 2007;Jin et al. 2010). Table 2 lists some of the prominent literature on the use of meta-heuristics in LRPs.
It is reported that outbound logistics is more complicated than inbound logistics in a logistics network due to the product values, customer delivery requirements and stringent regulations (Wu and Dunn 1995). The environmental impact of product transportation to DCs and retailers can be reduced substantially with effective logistics management (Wu and Dunn 1995). Conventional logistics system does not serve this purpose effectively as it does not consider the environmental impact of the outbound logistics mechanism. Consideration of the environmental aspect of LRP in outbound logistics requires substantial changes in LRP formulation procedures.
Although literature on general sustainable operations is growing, literature specifically on sustainable location-routing is scant. Reverse logistics and green-vehicle routing methodologies are found in (Zhu et al. 2008;Erdogan and Miller-Hooks 2012). Wang et al. (2005) reports a trade-offs model between the cost factors and the environmental impact of a supply chain. A two-layer multi-objective sustainable location-routing model with its solution approach is reported in Validi et al. (2014aValidi et al. ( , b, 2015. Literature on conventional LRPs is rich in terms of both methodologies and applications (Laporte et al. 1988;Demir et al. 2014;Prodhon and Prins 2014;Drexl and Schneider 2015). The purpose of LRP methodologies is the simultaneous determination of the location of facilities and the routes of vehicles for product transportation (Laporte et al. 1988;Yu et al. 2010). A broad range of LRP applications is found in the literature (Perl and Daskin 1985;Alumur and Kara 2007;Govindan et al. 2014). A wide range of optimisation techniques, including applications of integer-linear programming are found in literature (Semet 1995;Ambrosino and Scutellà 2005).   Ant colony optimisation Bell and McMullen (2004), Bin et al. (2009) Bi-objective mixed-integer non-linear programming model with a modified archived multi-objective simulated annealing meta-heuristic algorithm Rayat et al. (2017) Two-stage stochastic programming model Rezaee et al. (2017)

Contribution and objectives
Driven by growing social awareness, current laws and regulations and competitive business opportunities, this article contributes to the current body of knowledge by addressing a number of inter-linked decisions in logistics network design, location decisions, routing decisions and sustainability. This article contributes to the state-of-the-art of location-routing literature in the following aspects evolved from the literature review: (a) How a bi-objective optimisation model for three-echelon distribution networks with sustainable low-cost outcome can be developed that captures several real-world factors difficult to reconcile by organisations?
(b) How are decision makers' preferences accommodated with regard to costs and CO 2 emission in optimising the logistics network design based on varying organisational strategies? (c) What is the effective approach to solve the proposed bi-objective three-echelon locationrouting optimisation model?
To address these research aspects a sustainable three-echelon bi-objective integrated location-routing optimisation model is proposed that addresses the demand side of three echelon SCs. The weighting procedures of Analytic Hierarchy Process (AHP) (Saaty 1994) is infused into the bi-objective 0-1 integer programming optimisation model in order to reflect the DMs' priorities in the vehicle selection decision-making process. The computational complexity of the NP-hard model necessitates the solution method to divide it into two inter-connected phases. Phase-I deals with the transportation mechanism from the processing plants to the multiple Distribution Centres (DCs) and from DCs to multiple retailers. Phase-II uses the outcome of Phase-I to locate the non-dominated Pareto realistic optimal routes among retailers. A Design of Experiment (DoE)-guided meta-heuristic-based robust solution approach, integrated with the Multi-Objective Genetic Algorithm of kind II (MOGA-II), is adopted to solve the computationally NP-hard integrated model. The final set of optimal and realistic solutions is obtained from the combination of both phases. The realistic solutions are prioritised through Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) (Hwang and Yoon 1981) by reflecting DMs' priorities. To the authors' knowledge, literature is scant in combining three inter-linked aspects of sustainable location-routing.

The proposed sustainable integrated model
The purpose of this three-echelon bi-objective sustainable model is to minimise the level of CO 2 emission caused from transportation and minimise a combination of costs on the demand side of distribution networks. The proposed model is formulated based on a set of realistic assumptions. The model considers three key players on the distribution side of a SC, viz., plants, DCs and retailers. Two fleets of vehicles/trucks are considered for transporting the products throughout the SC network. A fleet of trucks transport products from plants to DCs, and another fleet of trucks transport products from DCs to retailers and then from retailers to other retailers. Each route may be a combination of different types of roads. In every country different speed limits apply to different types of roads. Speeds in different types of routes are captured in the model by the use of an appropriate variable.
The proposed model allocates DCs to the plants and retailers to DCs, and routes the vehicles from plants to DCs, DCs to retailers and retailers to retailers thereby serving the demand-side of the SC. The problem is to identify open and closed facilities and the optimised routing pattern throughout the network while minimising: (a) the CO 2 emissions from transportation, and (b) the total costs of operating facilities, satisfying demand at each facility and total transportation costs. The DMs' subjective opinions are considered to identify the optimal location-routing pattern. The model is generic as it is applicable to any three-echelon logistics network.
The routes on the demand side of the SC start from plants to DCs, DCs to retailers and then from retailer to retailer. The sustainable element of the location-routing design is addressed in the model by introducing two sustainable components, viz. (i) a sustainable objective function, and (ii) an AHP-integrated sustainable constraint.
The three-echelon bi-objective AHP-integrated 0-1 mixed integer location-routing model is associated with 10 operational constraints. A set of 0-1 integer decision variables defines open/close plants and DCs (F s and E j ) and all real feasible routes from plants to DCs, DCs to retailers and routes connecting retailers to retailers (V s jt n k , L jit n k and O ii t n k ). A set of parameters ( p t n k s j , p t n k ji and p t n k ii ) represents the amount of CO 2 emission in each real feasible route from plants to DCs, DCs to retailers, and among retailers. Two sets of variables represent the sum of fixed costs ( f s and f j ) and sum of the variable costs (v s and v j ). Another set of parameters represents the cost of serving DCs from plants, retailers from DCs and retailers from other retailers (c t n k s j , c t n k ji and c t n k ii ) ( Table 3).

Objective functions
The objective function (1) of the model minimises total CO 2 emission from transportation between facilities i.e. plants to DCs, DCs to retailers and retailers to retailers: The objective function (2) minimises costs of operating a facility, cost of serving DCs and retailers and vehicle-routing costs:

Constraints and decision variables
Each DC is connected to only one plant, each retailer is connected to only one DC and each retailer is connected to only one other retailer (if on a multi-route): Constraint (4) limits the length of each multi-stop route. One route comprises a start point from a DC to a retailer and then serving another retailer. Multi-stops are considered from the DCs to the retailers and from one retailer to another retailer:

t n k ji
Cost of serving retailer i ∈ I using vehicle k ∈ K from DC j ∈ J Q s j Fraction of demand at plant s ∈ S that is served by DC j ∈ J c t n k ii Cost of serving retailer i ∈ I using vehicle k ∈ K and t n ∈ T from retailer i ∈ I Q ji Fraction of demand at DC j ∈ J that is served at retailer i ∈ I p t n k s j CO 2 emission from transportation t n ∈ T n , from plant s ∈ S to DC j ∈ J using vehicle k ∈ K and t n ∈ T Q ii Fraction of demand at retailer i ∈ I that is served by retailer i ∈ I p t n k ji CO 2 emission from transportation from DC j ∈ J to retailer i ∈ I using vehicle k ∈ K and t n ∈ T Each route on the demand side of the logistics network is connected to a facility: s∈S j∈J V s jt n k ≥ 1 ∀k ∈ K , ∀t n ∈ T (routes from plant to DCs) , j∈J i∈I L jit n k ≥ 1 ∀k ∈ K , ∀t n ∈ T (routes from DCs to retailers) , and i∈I i ∈I O ii t n k ≥ 1, ∀k ∈ K , ∀t n ∈ T (routes from retailers to retailers) .
A route entering a node (i.e., DC and retailer) must exit the same node: j∈J i∈I L jit n k − i∈I j∈J L i jt n k 0, ∀k ∈ K , ∀t n ∈ T (from DCs to retailers) and O i it n k 0,∀k ∈ K , ∀t n ∈ T (from retailers to retailers). (6) The sustainable three-echelon model considers that a route can operate out of only one facility: s∈S j∈J V s jt n k ≤ 1, ∀k ∈ K , ∀t n ∈ T (routes from plants to DCs) , j∈J i∈I L jit n k ≤ 1, ∀k ∈ K (routes from DCs to retailers) i∈I i ∈I O ii t n k ≤ 1, ∀k ∈ K , ∀t n ∈ T (routes connecting retailers).
The flow of the products from the supply nodes (plants and DCs) into the facilities is ensured: There is a restriction of throughput at each facility up to the allowable maximum limit at each node that links the flow variables to the facility location variables: A retailer must be assigned to a facility if the route leaves the facility: Constraint (10) is applied to multi-stop routes connecting an open DC to a retailer and then through the first retailer serving other retailer(s). Assuming that the retailers have no supply of products, open DCs and routes connecting DCs to served retailers are included in this constraint.
Capacity for each vehicle is restricted: The sustainable constraint is formulated by infusing the pair-wise comparison matrices of the AHP technique into the mixed-integer programming model to introduce flexibility in the form of prioritising the vehicles to be used for the transportation: where w mn is associated with the priorities of the DMs regarding the type of vehicles (T n ) and B m contributes to the parameters of the objective functions. The decision variables are: T n 1, if vehicle/truck t n ∈ T n is selected to transport the products 0, if not (18)

Case of a dairy processing firm
A case of a dairy processing industry's three-echelon SC, operating in the east of Ireland, is considered for implementing the model and its solution approach. The SC has 2 processing plants, 6 DCs and supplies 22 retailers. The model considers variable costs as the total costs of serving each DC at each plant and each retailer at each DC. The variable costs of the processing plants (v s ) is defined as: where a s variable cost (e)/unit of product shipped from a plant to a DC; r j demand (units of product) at each DC. One 'unit' refers to a two-litre carton of milk. Variable costs at the DCs (v j ) is defined as: where a j variable cost (e)/unit of product shipped from a supply point (here DC) to a consumption point (i.e., retailer); r i demand at each retailer. Average demand at each retailer is assumed as 2/3 of the total population at the location of the retailers. Heavy duty Diesel fuelled fully loaded refrigerated vehicles are considered for the transportation activities. An average speed on the road is assumed.
The volume of burnt diesel is calculated using Eq. (21). 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 [Eq. (22)]: Litres of diesel burnt in each path Fuel efficiency (L/km) × Distance (km) CO 2 emission from a diesel vehicle (kg) Litres of diesel burnt × 2.64 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. Equation (23) considers e1.53/l and e11.50/h as the cost of Diesel in Ireland and average wage of a heavy-duty vehicle driver: Cost of serving a route e 1.53 × Litres of diesel bunt per km + e (11.5 × Disance (km) z .
The variable 'z' represents the average speed in each path. The CO 2 emission during transportation between the DCs and retailers and among the retailers, and corresponding costs of serving each route ( p t n k ji and c t n k ji ) are computed. The costs of serving each route among retailers ( p t n k ii and c t n k ii ) are determined. For illustrative purpose three different types of vehicle with different levels of CO 2 emission and costs, are considered for transportation.
The weight matrix (w mn ) is determined using the pair-wise comparison matrix of AHP. The right hand side matrix of constraint (12), i.e.,B m , is found considering the average round off values of the CO 2 emission. The boundary values of CO 2 emission and costs of serving the routes are found. Tables 16,17,18,19,20,21,22,23 and 24 of appendix provide the collated data utilised for the dairy processing industry's three-echelon SC.   Table 4 illustrates a snapshot of the technical details for the two phases.
TOPSIS is utilised to analyse the selected feasible optimal solution sets (designs) from Phase-I and Phase-II. It prioritises the results with refined result sets facilitating identification of open routes among retailers along with the type of the vehicles for transportation and its costs and CO 2 emission. Phase-I of the sustainable bi-objective three-echelon locationrouting model considers a set of realistic assumptions, as illustrated in Table 5.
Vehicle routes among retailers are not considered in Phase-I. This is because of the computational difficulty for the NP-hard nature of the model. Therefore, the demands of the unserved retailers are met through other retailers in Phase-II. In order to formulate Phase-II of the model the following set of realistic assumptions are considered (Table 6).  i ∈ I and retailer node i ∈ I , then retailer i ∈ I must be assigned to retailer node i ∈ I and consequently to facility (serving retailer i ∈ I ) Constraint 9: AHP-infused constraint considering the stakeholders' priorities Non-negativity constraints Integer constraints Outcomes from Phase-II Routes connecting served retailers to un-served retailers from Phase-I As served retailers don't supply products, the routes connecting them to open DCs from Phase-I are considered in Phase-II

Results and analysis
MOGA-II optimiser is selected to implement the sustainable bi-objective location-routing model as it has smart multi-search elitism and rapid convergence. Table 7 elucidates the optimiser parameter set up.

Phase-I: Results and analysis
Phase-I is implemented with an initial population of 50 different designs in the DoE table consisting of 10 DoE sequences. Phase-I is executed with 250 generations that generates 12,500 real feasible results. A statistical solution summary (Table 8) on the computed maximum and minimum levels of CO 2 emission and costs based on the DoE tables is obtained.  If a route opens up to cover a retailer, the demand transporting through the route is known The number of vehicles required to transport products in each route are considered as an objective function coefficient   MOGA-II optimiser's performance is analysed using convergence plots (Fig. 2a, b). MOGA-II starts converging with mild fluctuations after 65 generations. The solution for CO 2 emissions converges in a steady manner as compared with costs.
The optimiser generates a feasible space of solutions guided by the DoE tables. Figure 3a, b illustrates the characteristics plot of CO 2 versus costs on the realistic result tables. In these plots colours and diameters of the bubbles are related to the objective functions of the model. In Fig. 3a F2 values refer to costs, whereas in Fig. 3b F1 values refer to CO 2 emission. The colour schemes of the plots indicate the feasible solutions, infeasible solutions and solutions with errors. These colour schemes are illustrated in the legends of the plots. Further, the plots contain unique identity (ID) numbers in green colours. The red coloured bubbles in the plots are not realistic in nature as these solutions involve high CO 2 emission and high costs. Therefore, selection of alternative solutions is confined within the blue coloured optimum realistic solutions of Fig. 3.
In order to evaluate the realistic results, 20 out of the 5540 realistic solutions are selected guided by box-whiskers. This selection refers to the results table represented by the 4-D bubble plots in Fig. 3a, b. The red coloured solutions, in the 4-D bubble plots, are not realistic in nature as those involve high CO 2 emission and high costs. Therefore, selected realistic solutions are taken from the blue coloured bubbles for further evaluation using TOPSIS.
The sustainable bi-objective location-routing model for a three-echelon logistics network is a strategic decision-making procedure. In such strategic decision-making processes, it is desirable to elucidate the ranking of the set of selected feasible realistic optimal designs according to the stakeholders' preferences. Nine different weight matrices are considered to compare the selected results. It is found that extreme decision matrices are those that represent  a situation where the stakeholder believes minimising CO 2 emissions is more important and minimisation of costs is less important to the stakeholder, or vice versa. Therefore, using a moderate weight matrix W 5 0.5 0.5 in TOPSIS the selected results are prioritised. 20 results out of 412 realistic results generated by MOGA-II are selected using TOPSIS and the top three selected designs are shown in Table 9.
A routing scheme for the first ranked result is obtained (Table 10). Further, in order to perform a test on the robustness of the model a scenario analysis (Table 11) on the closed routes is conducted. This scenario analysis for Phase-I of the model shows the effect on the CO 2 emission and total costs if the closed routes are open considering the moderate TOPSIS weight matrix, viz., W 5 0.5 0.5 . Performance evaluation of MOGA-II shows that the 20 selected results lie on the Pareto frontier ( Fig. 4) that follows the Pareto optimality, and it is very efficient.  A feasible space of solutions guided by the DoE tables is obtained (Fig. 5a, b) illustrating linear relationships between CO 2 emission and costs. The objective functions of the model solved in Phase-II do not consider plants and DCs as they were already covered in Phase-I. Therefore, the fixed and variable costs of operating plant and DCs and the CO 2 emissions from transportation (plants to DCs) are not considered in these objective functions. The feasible and optimal results are chosen from the blue coloured bubbles of Fig. 5 guided by a set of history diagrams and box-whisker plots.

Phase-II: Results and analysis
A total of 20 selected results are prioritised by TOPSIS using weight matrix W 5 0.5 0.5 , and the top three are picked for further analysis. These three results cover the rest of the retailers. Phase-II of the solution approach is designed to find the optimal realistic designs for serving the following un-served retailers from Phase-I, viz. retailers #08, 13, 14 and 15. Considering these three results, open routes among retailers are identified along with the type of the vehicles used for transportation and their corresponding costs and CO 2 emission (Table 13).

Final results
The final result is a combination of the Phase-I and Phase-II results. Table 14 presents the sustainable location-routing plan while Table 15 shows the number of vehicles and the    It is noted that the selected 20 different optimal and feasible design IDs for Phase-II of the integrated model (as shown in Table 12) lie on the Pareto front (Fig. 6).
The selected 20 different optimal and feasible results (designs) for Phase-II of the integrated model lie on the Pareto front (Fig. 6). The realistic optimum routes from the processing plants, through DCs, to the retailers are mapped in Fig. 7. Figure 7 is an illustrative example of the final results for the proposed sustainable three-echelon bi-objective location-routing model for the demand side of the supply chain.

Managerial implications
The supply chain distribution system consists of processing plants, DCs and retailers. The considered case is representative of many other real-world logistics networks with intermediate DCs between production facilities and retailers. With the help of this three-echelon bi-objective AHP-integrated location-routing model logistics managers will be able to allocate DCs to the plants and retailers to DCs, and routes the vehicles from plants to DCs, DCs to retailers and retailers to retailers. By using the model and its solution approach logistics managers can identify open and closed facilities and the optimised routing pattern throughout the network. The optimised routing pattern provides two advantages to managers while satisfying the demand at each facility and total transportation costs, viz. minimisation of (a) CO 2 emissions from transportation, and (b) total costs of operating facilities. The logistics managers can utilise the benefits from this model as it captures several real-world factors that organisations find difficult to reconcile. The proposed model bridges the gap between the existing distribution solutions and the environmental impact of distribution. The model considers capacities of the processing plants, DCs, retailers and the routes connecting them and allocates plants to DCs, DCs to retailers and retailer to retailers suitably. The routing of the vehicles serves the demand-side of the SC and simultaneously minimises the CO 2 emissions from the vehicles during transportation. Both the fixed and variable costs have been optimised. Further, the model eliminates organisational estimation from contrasting objectives (i.e. cost versus environmental impact). The sustainable model has been designed to be robust with an ability for logistics managers to take immediate decisions by opening alternative closed feasible routes during SC disruptions.

Conclusions
From a practical viewpoint, sustainable business development has become the norm for general organisational development. Such organisational activity is generally driven by necessity rather than for altruistic reasons. Such developments are generally driven by a combination of dictate (e.g. legislative rules and polluter pay principles), competitive pressures (e.g. consumers and competitors) and for identified efficiencies (e.g. reduced cost of operations). This article contributes to this growing sustainable business development domain through the formulation of a sustainable location-routing solution approach for a three-echelon logistics network. More specifically, an innovative sustainable integrated bi-objective 0-1 mixed- integer programming method has been presented based for the three-echelon dairy processing logistics network. The implemented solution approach consists of a DoE guided efficient genetic algorithm based meta-heuristic with TOPSIS being used for the purpose of ranking sets of selected results. The solution approach has been developed to allow decision makers to develop routing plans based on varying organisational strategies. For example, decision makers can input into the model their own organisational preferences in terms of the influence of costs with respect to environmental performance. Organisations that face strong environ-  mental influences (e.g. consumer pressure, legislation, etc.) can place a heavier weight on this factor through AHP in the model. Likewise, an organisation that favours cost or even a more balanced approach can influence these in similar fashions. Each setup of the model will provide a set of results that are attuned to the organisations strategic preferences.
The model presented in this paper has been specifically tested in the case of the threeechelon dairy logistics network. Although this model has been illustrated in this way, the underlying DoE-guided unique solution approach and the analysis procedure do not have any geographical nor specific structural limitations in terms of the number of entities in each of the three echelons (e.g. processing plants, DCs, retailers). In short, the sustainable bi-objective three-echelon location-routing model, its solution approach and analysis procedure can be implemented to any three-echelon logistics network.

Future scope of research
The presented model clearly presents a novel procedure for evaluating the bi-objectives, viz. cost and environmental impact associated with CO 2 emissions for logistics distribution system. Although the model is innovative in nature and extends the current LRP domain, there are a number of extensions which would add additional value. In the first instance a fourth echelon to the model, incorporating grocery shops/super markets would be a useful extension to the current presented model. Furthermore, an extension to the current representation of environmental impact as CO 2 emissions would provide for a more comprehensive depiction. Such an extended representation could include the introduction of carbon taxes, carbon caps and trading and carbon offsets in the proposed model using appropriate objective functions. In addition, collaborative partnership among the multiple processing plants, distribution centres and retailers also provides opportunity for future research. The proposed model does not currently consider variability of demand under dynamic conditions. The variability in demand could possibly have a significant effect on the selection of the type of vehicle based on capacity. Concerns over the NP-hardness of the current model limited the inclusion of demand variability in the proposed methodology, but would be a useful area for further research. Inclusion of uncertainty issues can also be another area for future research directions. Uncertain issues like catastrophes can be used in the dynamic model to detect the vulnerability of individual nodes (i.e., open and closed routes).  National primary and secondary routes (dual carriageways included) 100 80 Regional and local roads 80 50 Built up areas (town and city) 50 30

Table 21
Total distance and the road types for the routes between plants and DCs (Validi 2014) DC ( j   (Validi et al. 2014a(Validi et al. , b, 2015Validi 2014) Vehicle type CO 2 emission Costs Vehicle-1 (T 1 ) Medium Medium Vehicle-2 (T 2 ) L o w H i g h