Design and Analysis of Novel Hybrid Multi-Objective Optimization Approach for Data-Driven Sustainable Delivery Systems

This study presents a novel two-stage solution method designed for sustainable last-mile delivery systems in urban areas. A proposed hybrid solution methodology includes multi-criteria decision-making system to select the most efficient logistics providers by considering different performance indicators, and a mixed-integer linear optimization model for last-mile cargo distributions by drones within metropolitan areas. We present a multi-objective modeling approach by considering time windows for customer services and charging operations of drones and outline important characteristics of the mathematical programming problem to minimize transportation cost (in the meantime carbon dioxide emissions) and total sustainability score of the system by using epsilon constraint method to find out the Pareto frontiers. The main novelty of the proposed solution methodology is the inclusion of many performance indicators of last-mile delivery systems into multi-objective models for design of a sustainable city logistics. Additionally, the proposed model is applied to an illustrative case by using real-life data of one of the metropolitan in Turkey. The approach is shown as comparative analysis of proposed method with other two state-of-art solution methodologies for multi-objective problems, after defining some pre-processing, symmetry breaking steps, valid inequalities, and logic cuts.


I. INTRODUCTION
Companies operating in last-mile deliveries search for new and adaptive systems to survive under harsh competitive conditions. As a nature of last-mile delivery systems, it is the movement of cargoes from transportation hubs to final delivery destinations (mostly individual customers because of increased e-commerce transactions). However, although logistics centers are moved from the city centers to rural areas due to lower real estate costs, more multi-tier distribution systems are developed and expectation to adapt for the last-minute changes (especially on-demand services and urgent delivery requests) or maintenance/break-down cases leads to excess expenditures, longer delivery duration and increased rate of unsatisfied customers. Therefore, proposed systems should ensure redesign and preparation of advanced and dynamic plans for routing problems because customers expect rapid interactions with their suppliers in The associate editor coordinating the review of this manuscript and approving it for publication was Baozhen Yao . delivery or return processes. This situation leads to increase loyalty and satisfaction ratios. Apart from the individual effect of last-mile delivery operations over citizens, there is a significant effect over energy consumption and environmental issues (such as greenhouse gas emissions, carbon footprint, etc.). The efficiency of the last-mile operations not only influences the profitability of retailing but also affects environmental and social performance criteria such as emissions and traffic congestion in regions. For instance, traffic caused by urban cargo delivery accounts for about 10-15% of kilometers travelled in city centers and emits approximately 6% of all transport related greenhouse gas (GHG) emissions [1]. 20-25% of freight vehicle kilometers is related to goods leaving urban areas, and 40-50% is related to incoming goods. The remaining percentage relates to internal exchange (i.e. goods having both their origin and destination within the city) [2].
Developed models or systems should also consider the effect of energy usage and pollution reductions in delivery processes under certain scenarios. One of the innovative solutions for this purpose is the usage of drones or unmanned aerial vehicles in last-mile delivery processes. This concept is getting more attention and popularity after observing some real-life cases (such as Amazon Inc., DHL, etc.). Although, there are many advantages of this new developed system over management of cost, time and carbon emissions, some critical points occur when real-life applications are considered. For example, storage problems of drones at distribution or transshipment centers, bad weather conditions, personal/private data protection, range and capacity limitations of drones, etc. However, when comparing heavy/medium trucks in city centers or at the location of final customers in delivery systems, drones offer significant operational advantages. For example, a drone system can deliver goods directly to the person who ordered it or to pick up the order for cancellation cases in less than 30 minutes and this duration reduces lead time significantly and eliminates route restrictions and many other logistics obstacles. Distribution of cargoes in a combination of both trucks and drones simultaneously is a very complex problem and needs very extensive studies in terms of solution methodologies. Although, vehicle routing problem considering different perspectives (such as; number of vehicles, variety of vehicles, capacity of the vehicles, time windows of customers, type of service provided, number of warehouses, uncertainty in distribution times, etc.) were studied extensively, inclusion of drones in delivery process is a very hot topic and there is still serious gap in literature in terms of modelling and computational side.

II. LITERATURE REVIEW
The survey of Lin et al. [3] gives and analyses an exhaustive literature review over the vehicle routing problems and classifies them into different application domains. According to this survey, green vehicle routing problems should concern energy usage of the systems to reduce fuel consumption that directly affects the greenhouse gas emissions and increases transportation efficiency. In this perspective, Bektas and Laporte [4] presents extended version of classical vehicle routing problem (VRP) by considering the amount of greenhouse emissions, travel times and total transportation costs. In their study, they developed a mathematical model with comprehensive objective functions to indicate the vehicle loads and speeds. In addition, Conrad and Figliozzi [5] presents a routing problem for electric vehicles to decide on the locations of charging stations as the Recharging Vehicle Routing Problem (RVRP). They try to decide on recharging operations of the electric vehicles during the transportation operations in which vehicles can be recharged in the locations of customers instead of dedicated stations. Therefore, there can be some savings in the total transportation times because there can be both loading activities and recharging of the vehicles instead of spending time at the dedicated stations for recharging. Erdogan and Miller-Hooks [6] introduces the Green Vehicle Routing Problem (G-VRP) by considering the recharging of the vehicles with limited fuel capacities. In G-VRP problem, vehicles are eliminated from the risk of running out of fuel and service time of each customer and maximum service duration restriction is posed on each route. After this study, Schneider et al. [7] proposes an updated version of G-VRP by including customer time windows for the delivery operations in which electric vehicles are used. Murray and Chu [8] introduces the vehicle routing problem with drones (VRPD) in which a drone collaborates with a truck to distribute customer parcels within minimum delivery time. A mixed integer linear programming (MILP) formulation and a heuristic adopting ''Truck First, Drone Second'' idea are proposed and tested on the instances with 10 customers. Ponza [9] presents a solution methodology for VRPD problems by using both MILP model and simulated annealing algorithm. Ferrandez et al. [10] proposes an optimization model of a truck-drone system in bicycle delivery networks by using the K-means algorithms to find the most efficient launch locations as well as using a genetic algorithm to assign the truck route between starting nodes. Wang et al. [11] introduces a more general problem considering multiple trucks and drones with the objective of minimizing the total duration of the delivery mission. Carlsson and Song [12] generalizes the VRPD in that two vehicles can meet anywhere, at not only depot and customers. Dorling et al. [13] proposes Drone Delivery Problem (DDP) by considering two VRP based models that are Minimum-Time-DDP (MT-DDP) and Minimum-Cost-DDP (MC-DDP). Agatz et al. [14] studies the VRPD with the objective of minimizing total logistics cost by using route-first, cluster-second heuristics based on local search and dynamic programming. Wang and Sheu [15] proposes a MILP model and a branch-and-price algorithm for VRPD. Different from the classical VRP, there are two types of vehicles in proposed model: a drone may have multiple times of flying and landing, each of which may be associated with a different truck; and a truck may launch and collect multiple drones at different times and locations. Karak and Abdelghany [16] develops a methodology that extends the classic Clarke and Wright algorithm to solve the hybrid vehicle-drone routing problem. Kitjacharoenchai et al. [17] presents a MILP formulation with the objective of minimizing the arrival time of both trucks and drones at the depot after completing the deliveries. A new heuristics algorithm is also developed to solve large sized problems containing up to a hundred locations. Compared to the existing studies, the main novelty of this study is the inclusion of large set of sustainability indicators into model. Most of the existing studies ( [9]- [17]) deal with only a single objective function (mostly total transportation cost or time) for vehicle routing problems by ignoring multi-criteria assessment in a perspective of multi-objective cases. In addition, multiple truck and drone combinations are presented in most of the studies, but none of them may consider re-charging operations of drones at nodes and they only allow operations of drones under fullcapacity battery conditions that leads to misrepresentation of drone capabilities in terms of flight ranges. This study also considers recharging of the drones and traces the battery levels and positions of individual drones for different logistics VOLUME 8, 2020 providers. Most of the existing studies, there is dependency to fly back to the main hub (depot or truck) after drop operation is completed. However, this study enables to either fly back to depot to deliver the next parcels or fly directly to another customer for pickup based on their charge-levels. In addition, most of the existing models are solved by using small-scale illustrative data sets, there is no testing and validation of proposed models in real-life cases.
Considering a single objective function in getting a feasible solution set by using exact solution methodologies for such complicated problems may not be so realistic in real-life cases. Therefore, usage of the meta-heuristic methodologies can be alternative solution. Dokeroglu et al. [47] shares an extensive survey about new generation meta-heuristics algorithms. Apart from classical methodologies, (e.g. the genetic algorithms, ant colony optimization, particle swarm optimization, simulated annealing techniques), nature-inspired population-based optimization algorithms are discussed such as Cuckoo Search [48]; Firefly Algorithm [49]; Grey Wolf Algorithm [50]; Whale Optimization [51]. If harsh competition conditions are considered, objective of minimizing total operational cost gets much more importance, but inclusion of sustainability factors (such as; time, GHG emissions, customer satisfaction rates, etc.) will lead to more realistic models and solution sets for decision makers. Therefore, inclusion of other objective functions will convert the problem into bi-or multi-objective problems. There are several methodologies to get non-dominated solution sets or Pareto frontiers for multi-objective problems, which are weighted sum and epsilon constraint models in general. In this study, the epsilon-constraint method is used because the weighting method may sometimes miss some of the efficient solutions in some feasible regions [18]. In the epsilon-constraint method, one of the objectives is taken as a single objective function and others are included into a model as constraint. Therefore, a set of solutions can be found by improving one objective and worsening others. Several approaches are proposed in the literature for the epsilon-constraint method. The augmented epsilon-constraint (AUGMECON) method by Mavrotas [18], its improved version (AUGMECON2) by Mavrotas and Florios [19] proposes solution methods to find a set of exact solutions in multi-objective problems by adding slack variables into the objective functions that are taken as constraints and penalty term into the main objective function. Fattahi and Turkay [20] presents a novel one-direction search (ODS) method to find an exact non-dominated frontier of bi-objective MILPs. Moreover, multi-objective evolutionary algorithms (MOEAs) are widely used in getting solutions for large-scale real-life problems. In these algorithms, several conflicting objective functions are included in order to increase the performance of decision-making problems. Therefore, as given in study of Zhang et al. [52], MOEAs are classified into three main segments that are; • algorithms based on Pareto-dominance (e.g. nondominated sorting genetic algorithm II (NSGA-II) [45]; strength Pareto evolutionary algorithm 2 (SPEA2) [53]; Pareto envelope-based selection algorithm II (PESA-II) [54] and cellular multi-objective genetic algorithm (cellular MOGA) [55]); • algorithms based on specific indicators (Zitzler and Kunzli [56]); • algorithms based on decomposition methods (e.g. the evolutionary algorithm based on decomposition (MOEA/D) [57]; MOEA/D-DE [58]) If sustainability factors are considered, multi-criteria decision-making models should be used for pairwise comparisons that can allow decision makers to weight coefficients and compare alternatives with relative ease [21]. Noorizadeh [22] proposes data envelopment analysis (DEA) to select green suppliers to reduce carbon emissions. According to Kannan et al. [23], goal programming (GP) is a useful method to solve multi-criteria decision-making problems since it can easily be applied to solve complex problems. Jolai et al. [24] obtains the efficiency ranks of suppliers and chooses the higher ranked ones with using a fuzzy multi-criteria method. Then, they calculate order quantities of suppliers using a multi-objective mathematical model. Ku et al. [25] presents an approach combined fuzzy analytic hierarchy process (FAHP) with fuzzy goal programming (FGP) methods for solving supplier selection problem. Boran et al. [26] uses a hybrid approach combined intuitionistic fuzzy sets with Technique for Order Preference by Similarity to Ideal Solution (TOPSIS) to apply a decisionmaking problem in supplier selection. Liu and Zhang [27] proposes a new method that is a combination of entropy weight and an improved Elimination and Choice Translating Reality (ELECTRE) III method. Weights of each indicator are determined based on entropy and all alternative suppliers are ranked according to their strong abilities. Koksalan and Ozpeynirci [28] proposes an interactive approach combined the utilities additives discriminants (UTADIS) method with approach of Koksalan and Ulu [29] to sort non-reference and reference alternatives. UTADIS is used for estimating the additive utility function that uses alternatives assigned to categories. To prioritize supply chain risks, Prasanna and Kumanan [30] proposes an approach combining AHP with Preference Ranking Organization Method for Enrichment Evaluation (PROMETHEE) for plastic industry. Resat and Unsal [43] presents a novel two-stage hybrid solution method designed for sustainable supply chain systems and their application in packaging industry. This model includes both analytic hierarchy process (AHP) to manage sustainability indicators and multi-objective aggregate planning model to design a sustainable supply chain network. However, this model highly depends on expert opinions and lacks of data verification and assessment scale in multi-criteria decision systems.
In addition, routing problems under real-time conditions includes many uncertainties due to their dynamic environments. Therefore, there are many studies to deal with such uncertainties to design more reliable and risk-free systems. For example, Meng et al. [37] proposes a novel methodology for linear models to deal with reliability conditions robustly and efficiently. Also, Meng and Keshtegar [38] shares some approaches to get feasible solution sets for reliability based design optimization models under some probabilistic cases. Lee et al. [39]; Gomez et al. [40] and Ehmke et al. [41] present vehicle routing problems with some uncertain time windows. Their computational experiments shares some solution methodologies to deal with probabilistic travel time and demand distributions for vehicle routing problems. Resat [42] presents a simulation based optimization model for the risk management of freight transportation by considering sustainability aspects.
The main novelties of this study are to: • Design a multi-objective optimization framework in order to provide alternative solution sets with different possible routes and schedules; • Develop a MILP model for drone based routing problem with time windows for last-mile delivery systems; • Consider sustainability concept by integrating as much as possible performance indicators of logistics providers into proposed model by using TOPSIS method; • Implement a hybrid methodology into the proposed model and validate it by using real-life data sets.
• Share a comprehensive analysis of proposed model by comparing with two state-of-art solution methodologies for multi-objective problems in terms of different performance metrics.
The rest of the paper is structured as: After this comprehensive literature review and problem definition section, details of proposed methodology for multi-objective hybrid model will be given in Section III. Computational details and data used in the model are given in Section IV. Details of multicriteria decision analysis model and developed a MILP model are shared in Section V. After obtaining outputs of study for decision makers, some future works and assessment of the findings will be given in Section VI.

III. METHODOLOGY
A two-step hybrid approach is proposed in this study to be able to add the largest number of indicators that can be used in the mathematical model of sustainable last-mile delivery processes. The main steps of the methodology are shared in Figure 1.
The proposed methodology includes mainly three steps: First one is the collection of relevant data ( kαi ) for logistics provider selection in last-mile delivery systems from final customers. The main critical point in data collection part is to assess criteria for logistics provider selection of final customers in cargo deliveries. Key performance indicators of the delivery operations affecting customer decisions and used in the proposed model are indicated in Table 1.
The defined indicators given in Table 1 are obtained by using semi-structured interview method. The surveys are conducted in both verbal and written forms in a group of 250 customers and collected scores of different cargo delivery com- gives the score table of each logistics provider per performance indicator and this table will be the main input of the TOPSIS system. Then by following Steps 1-2 in Algorithm 1, normalized decision matrix is constructed. Decision makers firstly define associated weights (χ α ) based on their relevant variances by following Steps 3-5. Then, by multiplying normalized matrix with associated weights, weighted normalized decision matrix is created under Step 6. Final performance scores of the logistics providers (π k ) are searched at the point where it is closest to the positive ideal point ( + α ) and as well as furthest position to negative ideal point ( − α ). Then, positive and negative ideal points are found (Step 7) and best feasible solution set is shared with decision makers by using their relative proximity.
In the second step of proposed model, a complex decision problem is structured at hierarchical levels and decision alternatives are generated to decrease complexity in multi-criteria problems by using Algorithm 1. Decision makers can reach alternative solution sets by performing different associated weights. Algorithm 1 given below gets final scores of different last-mile delivery companies under defined performance indicators [31].
After collecting final performance scores of different companies (π k ), final scores of companies are given into mathematical model as parameter. Pareto frontier for sustainable distribution channels in city centers are tried to be obtained. As indicated in study of Resat [32], green zones should be created within the city centres in order to achieve more To satisfy this condition, drone and conventional truck combinations are allowed in the proposed model and there are many assumptions used to simplify and reduce the complexity of the operations that are listed as below.

A. ASSUMPTIONS
• Logistics operations start from the distribution centres (DC) located outside of the city centres, a single type of medium truck (with capacity of a ton) from distribution centres to transfer points carries the cargo. When trucks come to transfer points (TP), there will be transfer from conventional vehicles to drones; • Standard type of octocopter drone with lithium-based batteries is used in predefined green zones (highly crowded and carbon-dense areas) to reduce environmental effects due to heavy road transportation activities in order to carry a single parcel from transfer points to final destination of customers; • Drones can either fly back to transfer points to deliver the next parcels or fly directly to another customer for pickup depending on their charge-levels; • Single type of product (a standard parcel) is considered for e-commerce transactions under last-mile deliveries with a dimension of cargo box 42*42*38 cm (Width-Length-Height) and it should be up to 5 kg due to drone capacities; • The capacities of drones are also assumed as constant. Each drone can carry 5 kg cargo with a speed of 30 km/h in delivery processes. Range of each drone with full battery is assumed as around 5 miles (app. 8 km); • The degree of extra warehousing requirement for drone delivery dominates the comparison of drone and truckbased scenarios, because non-operating drones should also be stored in some specific warehouses. For servicing an urban area with on-demand delivery, two main approaches have been proposed. The first is to locate distribution centers such that all of the service area is within delivery range of a distribution center. The second is to establish way-stations such that drones can fly from one to another and exchange batteries in a series of hops from distribution center to customer destination [33]. In this study, way-stations (transfer points) are considered and drones are stored at these places and delivery operations start from these stations; • Due to lower battery lives of drones and flight ranges, there will be some options to recharge these drones during the routing. In our problem, it is assumed that each drone can be recharged at only recharging stations. In addition, there is no limitation in the recharging ratios that means that drones can be recharged less than the maximum level. Drone batteries are recharged with constant ratios and each battery has a constant capacity; • Decision makers define associated weights based on their relevant variances; • All of the customers have to be visited once during a day; • All the routes are started in the transfer points and ends again in the transfer points (not exactly the same one); • Different logistics providers can use and store their drones at different transfer points; • The total demands of the customers in the green zones should not exceed the truck capacities (back order option is eliminated); • The deliveries should be managed within time windows of the customers; • Transportation times are directly depending on Euclidean distances between nodes; • The distances between the routes are taken as symmetric (d ij = d ji ); • Demands of the customers are known beforehand.

B. MODEL FORMULATION
The main purpose of this study is to create a multi-objective MILP model to find all possible tours covering all customers exactly once while minimizing total transportation distance (cost), as well as minimize total carbon emissions and maximize total supplier score in terms of sustainability.  Inventory cost keeping a drone at transfer points κ The pi constant ζ Overall power efficiency of the drone η Number of rotors ρ The density of air g The (1) m∈M n∈N x mn ≥ maxCap conv y n , ∀n ∈ N (12) i∈I n∈N k∈K v∈V k∈K v∈V i∈I x mn ≥ 0, ∀m ∈ M , ∀n ∈ N The objective function (1) tries to minimize total cost of the last mile delivery operation and includes five parts. First and second parts show total transportation cost occurred inside and outside of the green zones respectively. These expressions use total distance travelled in/outside of green zones in determined time duration. Third part ensures total transfer cost for transfer of cargo from conventional vehicles to drones, forth part is about total cost of opening transfer points on the borders of green zones. Last part gives the total inventory holding cost of unused drones kept in transfer points. Objective function (2) tries to maximize total sustainability score of the system by using individual sustainability scores of the logistics providers. Objective function (3) shows the third objective function by minimizing carbon emissions of the drone delivery system in green zones. Constraints (4) ensure that each vertex between customers will be visited exactly once by each drone, so this constraint enforces the drones to complete their tours. Constraints (5) ensure that incoming and outgoing flows has to be equal to guarantee the elimination of sub tours for each drone. Constraints (6) indicate that there are at most ''max'' drones going out from transfer points and Constraints (7) show at most ''max'' drones can come into the transfer points. Constraints (8) show the time balance of arrival of the drones to the nodes. The difference between arrival times of the nodes has to be more than the summation of total travel time and service time within the node for each drone. Constraints (9) guarantee that leaving time of nodes should be within the customer time-windows (the opening and closing time of nodes). Constraints (10) ensure the balance of total energy consumption of the drones while travelling between the nodes and total energy load of the drones should be until the maximum level of the battery. Constraints (11) ensure that drones of each logistics provider start with full loaded batteries initially at transfer points. Constraints (12) show the distribution of the cargo from distribution centres to the transfer points and this constraint forces the model to give a decision about opening transfer points. Constraints (13) ensure that total carried cargo has to fulfill the total demand of the customers in the green zone. Constraints (14) indicate that total number of visits of supplier k from transfer point n can not exceed the total number of customers. Constraints (15)

IV. COMPUTATIONAL EXPERIMENTS
An illustrative example is carried out by using real-life data obtained from final customers. Aim of this case study is to represent two-stage hybrid optimization approach and demonstrate detailed calculations for sustainable delivery systems. While designing the last-mile delivery system, not only cost minimization but also environmental and sustainability effects are considered.

A. DATA
In this section, all necessary data for comprehensive analysis of design and development of proposed sustainable delivery system in the illustrative example are shared. Firstly, locations of the customers are obtained from webbased mapping service of one of the biggest search engine providers (Google Maps Platform) and given in Table 2. The real-coordinates of the customers are converted into some integer values and x-and y-axis of the graph are located at the nodes where their coordination levels are set at lowest (yaxis) and on right-hand side (x-axis) of the sample network. Euclidean distances between distribution centres, transfer points and customers are calculated by using Constraint (21).
Secondly, the energy usage and environmental impacts of drone delivery system depend substantially on the range of drones and manner of implementation. Some specific parameters used for drone services in this study are listed as: • Standard type of octocopter drone with lithium-based batteries is used in case studies to carry a single parcel from transfer points to final destination of customers, and then returning empty to the transfer points; • Cargo capacities carried by standard drone systems can range between 0,3 to 5 kg and their speed range is between 20 to 120 km/h [34]. Therefore, it is assumed that each drone can carry 5 kg cargo with a speed of 30 km/h in delivery processes. Range of each drone with full battery is assumed as around 5 miles (app. 8 km); • Drones are expected to fly up to max. 120 minutes with full battery capacity and not carrying any payload [33], therefore it is assumed that their flight capacity can be max. 100 minutes in the air; • In the objective of minimization of greenhouse gas emissions and energy use in the delivery systems by using, drone energy calculations are included as a one of the objective functions. The masses of the drone body, battery, and package, the gravitational constant, and the total drag force are used in the calculations and parameters are taken from the study of Stolaroff et al. [33].
Thirdly, some operational limitations and scalars are included into model. Such as; number of available drones per each logistics provider, warehouse capacities for drones, time windows for each customer, charge consumption and recharge rates, some cost parameters for direct delivery and transshipment operations. The proposed MILP model for this problem is written in GAMS modelling environment and solved with IBM ILOG CPLEX 12,1 [35]. Both models are executed on a computer  with Intel Core I5 2520 M CPU with 2,50 GHz dual core processor, and with 4,00 GB of RAM. An optimality gap of 1% is set for the solutions.

V. RESULTS
In this section, details of the illustrative example are given and main concepts and analysis of solution on an illustrative case are discussed. Case study is started to find final sustainability scores of four different logistics providers (L1-L2-L3-L4) operating in district of Izmir, Turkey. The data collection methodology is applied more than 250 final customers who take cargo delivery services from these four suppliers. Therefore, final score table of each provider per each performance indicator ( i∈I kαi |I | ) are given in Table 3. The scores and associated weights for defined indicators given in Table 3 are obtained by using semi-structured interview method for different logistics providers. The most important indicator of logistics provider selection process in last-mile deliveries is the delivery duration of the providers after orders are set by the customers. Safety of cargoes and private data follow this criteria. Least important ones are related with environmental conditions such as different types of pollution.
After final score table is created for logistics providers, data set is given to Algorithm 1 and final scores of the suppliers (given in Table 4) are obtained for the mathematical model.
Graphical representation of outputs of the proposed mathematical model are given in Figure 2. In the illustrative case, VOLUME 8, 2020  four distribution centers (DC1-DC4) are defined at the outside of the city centers. Medium trucks are used in delivery of total cargo between distribution centers and transfer points (Node 0-Node 7-Node 46-Node 64). After cargoes are arrived to transfer points, individual delivery of the parcels is managed by using drones.
Details of routes of individual drones of different logistics providers and their corresponding travelled distances under illustrative case are shared in Table 5. The distances given in bold indicate that drones are charged at those routes by using charging stations. For example, logistics provider #3 follows path of Node 7 -Node 16 -Node 15 -Node 46 in Green Zone 1. However, total distance of this route is around 14 km that exceeds the range of drones (8 km), therefore there should be recharging process (at Node 15) to be successfully delivered to final customer (Node 46).
In Figure 3, details of Pareto frontiers of three objective functions are given. Hundred (100) different solution sets are obtained from proposed model for decision makers and different strategies can be followed to satisfy requirements of the customers as given in Table 6. Table 6 includes four different scenarios. Firstly, feasible solution sets are shared if the problem considered as single objective and corresponding values are obtained for different objective functions separately. Then, two different Pareto solution sets are chosen randomly and their comparison are given. If any stakeholder would like to spend much more money (from 0,15 to 3,68 mil. Euro), they can reduce their carbon emissions 93,1% and increase total sustainability score of the system 502,7%.
As given in Table 6, the proposed model shares better results compared to current status of delivery operations in Izmir. In terms of the economic assessment, the worst solution obtained from the proposed model gives around 8% better solution (from 4,59 to 4,24 mil. Euro) than the current status and 19% worse status is obtained even if the worst scenario in terms of the carbon emissions (from 68,9 to 82,1 kg-CO 2 e) is considered. The main reason behind these values is the  frequent usage of diesel-engine trucks in last-mile deliveries currently. The carbon emissions rate of the current status (case of truck usage) are calculated by using the equation given in Resat and Turkay [36]. The nadir points are constructed from the worst objective values over the efficient set of a multi-objective optimization problem. In the proposed multi-objective model, there are three objective functions. Two of them (total cost and carbon emission rates) are related with minimization and one of them (total supplier score) is about maximization, so maximum values of minimization functions and minimum value of maximization function are taken as nadir points and are shared in Table 6.
Moreover, there is an experiment to indicate computational complexity of the proposed model as given in Table 7. In this experiment, different number of nodes (including transfer points) are considered and their relevant CPU times are collected. As indicated in Table 7, an order for the growth rate of computational effort function in terms of number of variables is around 2 (shows O(n 2 )). Also, total CPU times in getting the Pareto set of proposed problem with 75 nodes by using NSGA-II and GoNDEF are around 3.145,87 s (39% less than proposed model) and 7.235,67 s (41% higher than proposed model), respectively.
Additionally, more sensitivity analysis is made to understand the effect of different associated weights of sustainability indicators that are given in Table 3. Therefore, three important indicators with highest weights over process according to answers of more than 250 final customers are chosen which are lead times, traceability and safety of cargoes. As given in Figure 4, their associated weights are increased 50%, 100% and 200%, and their impact over cost and total GHG  emissions are shared. While increasing their relative weights, weights of the other indicators are decreased in a ratio of their initial weights. For example, if the weight of safety of cargoes increased from 0,10 to 0,20, this additional 0,10 in the calculations will be reduced from initial weights of others by considering their initial weights, such as initial weight of lead time is 0,20 and after this change occurs in safety of cargoes, weight of lead time is reduced to 0,175 to satisfy the summation of all weights to one.
As given in Table 8, when associated weights of the performance indicators are increased, final sustainability scores of the logistics providers change in a ratio of their relevant scores that are given in Table 3. By using these final sustainability scores in proposed hybrid model, different Pareto frontiers are obtained for total cost vs total GHG emissions as given in Figure 4. The main outcomes of this analysis are that when companies save some money over their operations, impact of associated weights of lead time and safety of cargoes indicators will be higher over GHG emissions. Also, if decision makers focus on same GHG emission rates (such as range of 150-200 kg-CO2e) for different weight ratios, effect of indicator of safety of cargoes will be higher compared to other two indicators. If companies increased the weight of their safety scores, this leads to much more spending to satisfy the same level of GHG emissions in their system.
In order to justify the main advantageous of proposed system, a comparison analysis is made by using other two well-known solution methods for multi-objective optimization problems in terms of different multi-objective performance metrics. The main idea behind selection of these two solution methodologies is that although proposed solution methodology shares exact solution approach, how Pareto solution sets will be changed if another exact solution algorithm is used or apart from the exact solution algorithms, VOLUME 8, 2020 if one of the meta heuristics approaches, in which nondominated sorting technique is used to provide the solution as close to the Pareto-optimal solution as possible, is used. One of these chosen methodologies (GoNDEF) proposes an exact method to generate all non-dominated points of multiobjective mixed-integer linear programs that is provided by Rasmi and Turkay [44]. This method handles multi-objective problems in a combination of some sub-problems and creates potential sub-Pareto frontiers by testing some predefined integer solution sets and eliminating dominated solutions from non-dominated ones. Final non-dominated Pareto solution set is obtained by adding these non-dominated facets of the subproblems together. Other methodology (NSGA-II) is one of the state-of-art techniques that shares an evolutionary algorithm for multi-objective problems shared by Deb et al. [45]. However, design of NSGA-II algorithm directly depends on various parameters that may have different impacts on its computational effectiveness. Therefore, there are some assumptions over the numeric values of these control parameters in order to minimize perturbation in the design of problem.
• Population size is set as 100 in order to create more balanced sample set with other two methodologies • Number of runs in order to evaluate the problem is chosen as 1.000 due to the computational efforts of large-scale problems. In addition, as given in Figure 5, convergence level of the problem reaches its maturity level around value of 1.000 and after this level marginal increase in the convergence ratio may be negligible.
• The mutation probability is taken as 0,000023 (0,023%) that is obtained from 1/total number of variables in multi-objective problem.
• Probability of crossover is taken as one of the highest values which is 0,98 (98%). Riquelma et al. [46] presents a comprehensive survey about performance metrics of multi-objective problems and shares ranking of the most commonly used performance indi-  cators in comparison studies. Therefore, hypervolume (HV) and generational distance (GD) in an aspect of accuracy; spread delta (SD) measure and inverted generational distance (IGD) from diversity perspective are used in comparison analysis. Therefore, these four different performance indicators are investigated and compared by using three methodologies.
Moreover, generational distances (GD) and inverted generational distances (IGD) between optimal Pareto frontiers of proposed methodology and other solution methodologies are obtained separately to calculate the level of convergence. In the measurement of the GD and IGD, decision makers should measure the distance between the obtained Pareto set and a reference one. Due to the complexities in getting true Pareto set for real-life cases, Pareto set obtained by proposed model is assumed as true Pareto frontier and it is used in calculation of GD and IGD. As given in Figure 6, GoNDEF shows poor quality in approximation with higher GD value in all combinations of bi-objective functions. Lastly, spread delta metrics are compared with different solution methodologies to examine how evenly the solutions are distributed among the approximation fronts in objective space. NSGA-II shows generally wide and uniform spread out of the solutions across the Pareto front due to lower spread value, however GoNDEF gets better position in a combination of total GHG emissions and supplier score.
As given in Figure 7, volumes between generated nondominated points of the Pareto frontier and nadir points are compared with each other. GoNDEF shows a better coverage of the region compared to proposed model (blue region), however proposed model (purple region) slightly dominates objective space obtained by Pareto approximation set of NSGA-II (gray region). For example; although proposed model gets 93,8% of the objective space for the case of minimization of total cost and total GHG emissions, GoNDEF and NSGA-II present 0,95% higher (94,7%) and 2,86% lower (91,1%) coverage compared to proposed model, respectively.
As given in Table 9, several alternative scenarios are considered in terms of the true Pareto sets. If the Pareto set of proposed model is considered as true Pareto set, NSGA-II shows poor quality in approximation with higher GD value, however, when other two scenarios are considered, quality level of proposed model in terms of approximation is better than other two methods with lower GD ratios. However, although proposed model shows more fidelity if solution set of GoNDEF is considered as reference frontier in a case of IGD parameter measurement, GoNDEF gets better position with lower rates in other two options. For the case of spread delta, proposed model show better role compared to GoNDEF, but worse role than NSGA-II. Lastly, when volumes between generated nondominated points of the Pareto frontier for three objective functions and relevant nadir points are compared with each other, GoNDEF shows a slightly better coverage than then proposed model and NSGA-II is the worst one.
Non-dominated solution sets including 100 different solutions are obtained from these two known methodologies and are compared with the solution set obtained by using proposed methodology as given in Figure 8. Although, difference in the obtained Pareto sets from three methodologies are changing on average maximum 3,1% less in NSGA-II (−1% in cost; −3,1% in supplier score and −2,1% in GHG emissions) and 4.6% higher in GoNDEF (2% in cost; 3% in supplier score and 4,6% in GHG emissions) compared to proposed methodology, changes in Pareto sets can be negligible based on these deviations. Then, results obtained from three different methodologies (under Point 1 in Figure 8) are compared with each other. For this illustrative example, proposed model has less total cost and GHG emissions and higher supplier score than other two methodologies. Also, although NSGA-II shows better performance in cost and supplier score functions, GoNDEF has better GHG emissions rate. When details of obtained variables in solution sets are investigated, such conclusions are observed.
• In the proposed model, solution set includes higher usage rate of logistics provider 2 (L2) (which has the highest sustainability score) in overall operations in a unit of individual transactions (from node i to j), this ratio is equal to around 31% of all number of transactions (474) used in the system. This ratio is around 29% in NSGA-II and 28% in GoNDEF. This situation leads to higher sustainability score in the system. Usage rates of L1 in the solution sets are around 21%, 22%, and 21% in the proposed, NSGA-II and GoNDEF, respectively.
• In terms of the cost calculations, three critical cases are observed. One of them is about the distribution of cargoes from distribution centers (m) to transfer points (n). In the GoNDEF, total distance travelled by medium tracks (around 174,8 km) is 6,4% higher than the proposed model. One of the results of this case is supply of cargoes from DC1 to Node 0 (positions are given in Figure 2). Second important case is about the distances travelled by drones within the green zones. Although total number of transactions is the highest in NSGA-II (513 unit), total travelled distance of all drones are highest in GoNDEF (715 km) (this value is 678 km in proposed model and 699 km in NSGA-II). Thirdly, another factor affecting cost is number of drones waiting in warehouses. Surprisingly, although GoNDEF has highest total travelled distance of all drones, ratio of waiting drones at transfer points for all logistics providers is also highest with ratio of (4,2%). Also, majority of this non-operational drones in solution set is from logistics provider 2.
• Another critical approach is about recharging operations. In the proposed model, there are 74 charging station visits occurred. This solution leads to major benefit for proposed model that is the elimination of number of redundant return routes of drones (benefits over distance travelled and number of transactions) because of range limitations and masses of cargoes carried by drones from GHG emissions perspective.

VI. CONCLUSION
As sustainability concept is getting much more importance in global competition conditions, companies propose innovative and challenging solutions for daily life operations. Last mile cargo deliveries occupy one of the highest importance in logistics activities. Therefore, any improvement or solution methodology may lead to significant improvements in financial, environmental as well as social aspects. Drone delivery systems or drone routing problems are most popular and highly appreciated concepts in current conditions. This study focuses on novel solution methodology for hybrid drone vehicle routing problem in an aspect of sustainability. The novelty of the proposed model is the inclusion of large set of sustainability indicators in the management of city logistics activities by considering delivery systems via drones. The proposed model firstly integrates opinions of many stakeholders (final customers) in the assessment of the logistics providers operating in last-mile delivery systems. After collecting enough data for the proposed model, hybrid solution methodology including TOPSIS and MILP optimization tools takes and processes real-life data sets. Apart from the existing studies, proposed model can satisfy recharging operations of drones throughout the operations and not specifically dependent to the home depots or warehouses of the drones. At the end of system, different Pareto thresholds are obtained for decision makers and sensitivity analysis of the solutions are shared. One of the obtained solutions are compared with the current status of the delivery system in city of Izmir and improvements are shared in terms of total cost and GHG emissions, as well as traceability of the individual drones of the different logistics providers can be satisfied. In addition, proposed model is compared with two well-known solution methodologies for multi-objective problems and obtained Pareto solution sets are compared in terms of hypervolumes, generational distance, spread delta and inverted generational distances. Although proposed model gets better coverage of objective region than NSGA-II, GoNDEF covers solution regions more dominantly. Moreover, proposed model is superior to the other methods for all performance metrics except IGD. The decision makers can use this model in a ranking of the logistics providers in last-mile deliveries by integrating feedback of the customers under different performance indicators. In addition, obtained Pareto results will enhance the decision making processes and stakeholders may get some initiatives or priorities to support environmental conditions in city centers.
One of the expected studies in the near future is inclusion of all characteristics of the drones to create more realistic case and enlargement of the constraints to cover more operational restrictions (such as; topographical conditions; climatic conditions (wind, etc.); different type of drone usage; etc.). Additionally, some more sensitivity analysis over the variety of associated weights that are used in getting the final sustainability scores of logistics providers should be done to analyze how these parameters are important for the cases.
H. GIRAY RESAT received the bachelor's degree in chemical engineering from Bogazici University and the Ph.D. degree from the Industrial Engineering Department, Koc University, Istanbul, Turkey.
Since 2016, he works as an Assistant Professor at the Izmir University of Economics, Izmir, Turkey, and specializes in sustainable manufacturing and transportation solutions for industrial clusters. He works on the development of systematic approaches to systems engineering problems and addresses development of approaches to integrate data (both qualitative and quantitative) with physical models and development of mathematical programming algorithms to optimize system performance. The core research area is related to intermodal / multimodal logistic area and its applications in real life. In addition, it relates with integration methods for the intermodal logistics incorporating economics, service level, and environmental concerns. VOLUME 8, 2020