Last-mile drone delivery combinatorial double auction model using multi-objective evolutionary algorithms

This study proposes a combinatorial double auction bi-objective winner determination problem for last-mile delivery using drone. Prior studies are limited on solving mixed integer model, which are not efficient for large-scale scenario. However, this is not practical in real cases as the computation time to obtain the solution is longer due to number of combinations of packages and participants anticipated in the last-mile delivery platform. Four multi-objective evolutionary algorithms (MOEAs) with the decomposed winner determination problem model are experimented. This study is able to yield Pareto optimal solutions from multiple runs of mixed linear integer programming (MILP) using different objectives weights in the model. Unmanned aerial vehicle, or drone, has potential to reduce cost and save time for last-mile logistic operations. The result positively shows MOEAs are more efficient than MILP in yielding a set of feasible solutions for undertaking complex winner determination problem models. The percentage of improvement in terms of time spent identifying the best option is almost 100%. This is likely an unprecedented research in drone where combinatorial double auction is applied to complex drone delivery services and MOEAs are used to solve the associated winner determination problem model.


Introduction
Customers would like to buy delivery services are required to give the necessary information such as weights, volume and collection address, to an electronic transportation market (ETM) or an e-commerce platform. Firms can also pass their parcels to a local pickup or drop-off point, such as a SKY-BOX kiosk by aCommerce. The electronic transportation market is open lanes for the participated firms to bids. A firm determines its resources, constraints, or utilities issues before submitting a bid to the firm and both the customer and firm are allowed to submit their bids as packages. This means that the customer is allowed to Q3submit different volumes and addresses of the parcels that are to be delivered as a package bid, if customer is willing to pay only if all the lanes in the packages are assigned. A participated firm submits its proposed optimized routes according to its capacity and utility as a package bid, if the firm is willing to accept the business only if all the lanes in the packages assigned are given to the firm. Hence, this study models a Combinatorial Double Auction (CDA) bi-objective winner determination problem (WDP) for last-mile delivery using drones. CDA is a type of auction that is participated by multiple buyers and multiple sellers and both buyers and sellers are allowed to bid the packages. Allowing a package bid on a set of lanes means that if Lane AtoB ; Lane BtoC f g is submitted, both lanes in this package are adopted if this package is assigned.
In reality, Walmart introduced their last-mile delivery solution, i.e., Walmart To Go includes home delivery and local pickup options for their customers (Sarah 2014). A similar pickup option is also provided by aCommerce with its SKYBOX kiosk, offering their customer options for pickup and drop-off couriers while on the go (aCommerce 2016). In fact, Iwan et al. (2016) carried out an analysis on parcel lockers, where its concept is similar to SKYBOX kiosk and found this local pickup concept is the important directions for last-mile delivery. This shows that large firms have shown interest to use drones in their delivery system. In 2013, Amazon recommended Prime Air service to deliver packages using multi-rotor drones to customers (Amazon.com Inc. 2016). A large logistic firm in German, Deutsche Post DHL also started its Parcelcopter initiative in 2013 by transporting medicine to the island of Juist in the North Sea using their Parcelcopter (DHL International GmbH. 2014). In 2014, Google revealed Project Wing (Stewart 2014) that aimed to deliver larger items than those carried by Prime Air and Parcelcopter. Alibaba, the largest internet retailer in China, also started its drone delivery trials in 2015 (Leo 2015). This study is noted that using drones for last-mile delivery of light-weight parcels is also environment-friend, i.e., less carbon emission as compared with using truck (Goodchild and Toy 2018).
In theoretical perspective, supply chain management and transportation planning, last-mile delivery means delivery from a transportation hub to a final destination. As the transactions of consumer-to-consumer in e-commerce are increasing, last-mile delivery is becoming important in logistic firms. myHermes, one of the logistic firms, which was launched in 2009, already provides consumer-to-consumer services that collect couriers from the door-step of customers (Hermes Parcelnet Ltd. 2017). Last-mile delivery solutions includes home delivery, setting up local pickup and drop-off kiosks for customers, or even collect parcel from the doorstep of customers. Using unmanned aerial vehicles, or commonly known as drones, in last-mile delivery has potential to reduce the cost and time required to deliver packages as the maintenance cost of drones is less than trucks, and the tasks are performed autonomously. Besides that, delivering with drones is unconstrained by traffic conditions, and generally has few obstacles to avoid (Dorling et al. 2016;Marius 2017).
One distinct characteristic of a last-mile delivery system using drone is that the number of participants and the number of lanes will be very large. When the number of participants and number of lanes are very large, the WDP is not an easy problem to be solved using mixed integer linear programming (MILP). Solving a bi-objective WDP using the classical weighted objectives methods need multiple runs with different weights, in order to obtain the Pareto optimal solutions. This is not practical in real cases as the execution cost might not be affordable. This study adopts Multi-Objectives Evolutionary Algorithms (MOEAs), i.e., Non-dominated Sorting Genetic Algorithm II (NSGA-II), reference-point-based Non-dominated Sorting Genetic Algorithm (NSGA-III), Pareto Archived Evolution Strategy (PAES), and Strength Pareto Evolutionary Algorithm (SPEA2), to solve the complex bi-objective WDP model and produce the Pareto front of the solutions in a single run (Zhou et al. 2011). The main reason that the four MOEAs are implemented due to their versatility in solving different types of problems such as knapsack and combinatorial problem (Zhou et al. 2011). These MOEAs are commonly used to tackle real-world applications (Giagkiozis et al. 2015).
To sum up, in today's e-commerce, drones in last-mile deliveries are becoming increasingly important. In instances where ETM must evaluate several objectives, the present paradigm is insufficient. The WDP model must be extended to a bi-objective representation. The scope of the project is the bi-objective WDP model in the context of the Electronic Transportation Market (ETM). However, biobjective models are tough to solve, and large-scale and complicated models are much more challenging. Metaheuristic methods, which include MOEAs, constitute a practical approach to solve large scale instances in complex WDP models. Therefore the objective of this paper is firstly to extend CDA model to a bi-objective representation and applied to last-mile drone delivery services (ETM) and secondly to address the problem using MOEAs. This research is a step forward to make last-mile drone delivery become a common choice in delivery services. The CDAbased WDP model is able to collect a pool of demands and supplies while MOEAs are able to solve this large scale problem in short timeframe.
The contribution of this study is two-fold.
(1) a CDAbased bi-objective WDP model is proposed for last-mile delivery with drones; and (2) MOEAs are applied to find the solutions for the WDP model to handle large and complex scenarios. The results from different MOEAs are compared and analyzed. The convincing solutions indicate that MOEAs are useful for providing solutions for large and complex cases such as drone delivery services; therefore overcoming the inefficiency issue associated with large and complex WDP models in the marketplace. The rest of the study is structured as follows. Section 2 contains a literature review covering drone delivery, CDA-based WDP, and MOEAs. Section 3 describes the methodology used. Section 4 explains the details of our solution using MOEAs. Section 5 contains discussion on data preparation and results. Finally, conclusions and suggestions for further work are presented in Sect. 6.

Literature review
This section presents a literature review on drone delivery, CDA, and MOEAs techniques. Ducret (2014) presented parcel deliveries and urban logistics and found last-mile delivery is still understudied, which is generally more expensive, less efficient, and the most polluting section in a logistic chain. The main problems of last-mile delivery include uncertainty in terms of traffic that influences the delivery process and the substantial pollution emission (Tavana et al. 2017). Tavana et al. (2017) suggested that using drone can reduce the influence of the aforementioned issues and, therefore, the service provided is faster and more reliable. In fact, technological advancement has made drone delivery a reality, which is able to change the common choice of delivery options as drone is operated autonomously and unconstrained by road traffic conditions (Dorling et al. 2016;Marius 2017). Besides that, Goodchild and Toy (2018) estimated the carbon dioxide emission and vehicle miles travelled levels and presented that using drones, rather than trucks, can reduce carbon dioxide emission if the service zones are close to depot or the carried load is light in weight. Ghelichi et al. (2021) advocated delivering medical supplies to rural and suburban regions using drones. He looked at the crucial aspects of restricted cargo capacity, limited operational coverage range of each drone, energy consumption, and timeliness for delivering medicinal goods within a certain time frame. Meanwhile, Moshref-Javadi and Winkenbach's (2021) comprehensive review of use cases of drones for delivery revealed that the majority of the reviewed models were designed for applications in e-commerce and healthcare/emergency services, while other applications, such as food and mail deliveries, remain underrepresented in the academic discussion.

Drone delivery literature
Still, different methods have been examined to resolve the routing problems in drone deliveries. Sundar and Rathinam (2013) optimized the routing problem of single drone with multiple depots and multiple targets where the depots allow drone to refuel. Ferrandez et al. (2016) showed another way of trucks and drones collaboration. It aims to optimize the delivery network. Dorling et al. (2016) considered the issues of battery weight, payload weight, and drone reuse in a vehicle routing problem. Wang et al. (2019) even combined truck delivery, truck-based drone delivery and independent drone delivery platforms. Not to mention Sawadsitang et al. (2019) proposed a multi-objective and three-stage stochastic optimization model for the drone delivery scheduling. There are also multi-drone networks modeled using auction-based concept (Shin et al. 2019). Aforementioned, these studies conclude that with the current technologies, drone delivery still has limitations, but has some advantages over truck delivery in certain scenarios. However, the marketplace of e-commerce is a significant area for drone delivery to work well. There are a lot of transactions in e-commerce that the goods are light in weights (Wang 2016).
Previous research on drone delivery has mostly focused on overcoming routing issues. Unlike drone routing problems presented, this study proposes to use the WDP method to assign a pool of customers and logistic firms. A large enough pool of customers and logistic firms, despite whether drone delivery is used in tandem with trucks and the routing problem optimized by the logistic firms is more practical in the economic sense. Prior studies focused on routing from the technical aspect and the economic aspect of drone delivery is yet to be fully explored. As the use of drone in last-mile delivery is increasingly gaining importance in today's e-commerce era, the need to look at drone transport service procurement is timely and vital. The purchase of transportation services is a WDP concern. Lafkihi et al. (2019) have emphasized the necessity to create new transportation models in order to increase the efficiency and efficacy of the procurement process. As a result, it is critical to investigate novel methods, such as combinatorial auction drone delivery, to minimize logistical costs as well as other negative environmental externalities. The next subsection presents a discussion on a variety of methods to solve WDP in transportation service procurement.

Combinatorial double auction-based winner determination problem
There are many methods to solve WDP in transportation service procurement. In this aspect, transportation services need optimization-based procurement owing to the enormous economy prospect (Caplice and Sheffi 2003). The commonly used mechanism in solving WDP is auction. CDA is a type of auction mechanism that has both double auction (DA) and combinatorial auction (CA) properties. DA is an auction which considers a multi-seller and multibuyer environment. In this environment, a third party acts as auctioneer, collecting the bids from sellers and buyers. The auctioneer could also be an independent platform or cloud-based ETM and regulates the price that allows clearance or matching between bids. The DA mechanism has also been extended for a multi-attribute setting by Cheng et al. (2016) in the context of perishable supply chain trading. In addition, CA allows carriers to bid on multiple lanes in a package. Song and Regan (2003) applied simulation in CA and found that both shippers and carriers could gain significant benefits. This is because CA solves transportation procurement problems that account for carriers' economy of scope (Sheffi 2004). Ignatius et al. (2010) leveraged the imprecision of bid prices as a means to provide a better flexibility to WDP formed using CA. Furthermore, Ignatius et al. (2011) formalized a Fuzzy Combinatorial Auction Winner Determination Problem that allows the auctioneer to estimate its ''true'' revenue despite price uncertainties. Ignatius et al. (2014) were to consider multi-objective properties for winner determination in CA and provide solution comparisons across weighted objective methods, pre-emptive goal programming, and compromise programming. Accounting for the uncertainty of demand in transportation procurement, Ma et al. (2010) proposed a two stage stochastic integer programming for CA-based WDP. Later on, Remli and Rekik (2013) proposed a constraint generation algorithm to similar WDP with consideration of uncertainty in demand.
CDA has just been applied to transportation service procurement recently. Based on both DA and CA properties, CDA considers a multi-shippers and multi-carriers environment, where they are allowed to bid on multiple lanes in a package. Similar to DA, the auctioneer is a third party agent or an ETM. There are only few studies on CDA in transportation service procurement. For instance, Motlagh et al. (2010) derived a linear CDA model in transportation service procurement that allows carriers to constraint the range of volume of loads to be carried. The derived CDA model is compared with CA using a test procedure. The results show that CDA is superior to CA in generating revenue, and it remains relatively stable for reduced market clearing flexibility. The CDA is also used in business to business trading in a centralized marketplace, or in a multi-agent coordination system with artificial intelligence.
In the vein of techniques to solve CDA model, Xia et al. (2005) showed how a general CDA problem could be reduced to a combinatorial single-sided auction, which is a multi-dimensional knapsack problem. He contrasted the branch-and-bound method with the intelligent search method and discovered that theoretically, the LP relaxation bounds outperform the specific search strategies. In addition, Triki (2021) proposed a combinatorial auction model for the procurement of occasional drivers in the freight transportation. The model was solved by using two heuristic, one based on decomposition and the other on a cost-comparison greedy approach Hammami et al. (2021) proposed an exact non-enumerative and heuristics methods to solve the combinatorial bid construction problem with stochastic pricing in truckload transportation services procurement auctions. They run simulation on instances up to 50 contracts. Furthermore, Tafsiri and Yousefi (2018) used CDA to allocate resources in the cloud computing market, formulating the CDA model as an integer LP model and solving the model with a heuristic resource allocation algorithm with a quasi linear time complexity.
Other applications included native vegetation offsets in Nemes et al. (2008), a trust incentive problem in Wang et al. (2010), and cloud computing in Samimi et al. (2016) and Tafsiri and Yousefi (2018). In sum, accounting for cases that ETM needs to consider more than one objective, this study extends the CDA model to a bi-objective representation, and tackle the problem with MOEAs. A review on MOEAs is presented in the next subsection.

MOEAs
The WDP model needs to consider more than one objectives, which often are in contradiction in the real cases and the model is not possible to yield on solution that could minimize or maximize both objectives. A set of trade-off solutions is known as the non-dominated Pareto-optimal solutions (Deb 2001). The classical weighted objectives method is considered to solve multi-objectives problems (Geoffrion 1968). However, the solution set is obtained through multiple runs with different objectives weights, which is a time-consuming and laborious process. Population-based multi-objectives algorithms are able to produce a set of Pareto optimal solutions in single run and have been widely accepted as a useful method for solving multi-objectives problems (Nayyar et al. 2018a, b;Zhou et al. 2011). Metaheuristic methods, which include MOEAs, constitute a practical approach to solve large scale instances in complex WDP models.
MOEAs are developed from single objective evolutionary algorithm. Examples include the Multi-Objective Genetic Algorithm (MOGA), Non-dominated Sorting Genetic Algorithm (NSGA) and Niched Pareto Genetic Algorithm (NPGA). These methods convert a single-objective genetic algorithm to a MOEA by adding the necessary operators. In general, the two common features of MOEA operators are assigning a fitness to each population member based on a non-dominated sorting procedure and preserving diversity among the solutions of the same non-dominated front. To produce better solutions, the elitism strategy is adopted. Some elitist MOEA examples include Strength Pareto Evolutionary Algorithm (SPEA) (Zitzler and Thiele 1998), Pareto-archived Evolution Strategy (PAES) (Knowles and Corne 1999) and elitist objective genetic algorithm (Rudolph 2001).
There are also researchers that study on MOEA based on decomposition, MOEA/D where the MOP is decomposed to subproblems to be optimized simultaneously (Zhang and Li 2007). For dynamic cases, Jiang and Yang (2016) developed a steady-state and generational evolutionary algorithm. Some examples of applications of MOEAs include data mining, spectrum assignment problem and power flow problem. Coello (1999) introduced the strength and weakness of those MOEAs that were popular while Zhou et al. (2011) carried out a survey on MOEAs and presented the state of the art of MOEAs.
This study employs four widely used MOEAs, namely NSGA-II, NSGA-III, PAES and SPEA2, to tackle the CDA-based WDP model for last-mile drone delivery. A review on these MOEAs and their associated characteristic features are described as follows.
NSGA-II uses a fast non-dominated sorting procedure for ranking the solutions in its selection and crowding distance assignment (Deb et al. 2002). In NSGA-II, for each solution, two entities is calculated: • The number of solutions which dominate the solution (Domination count) • A set of solutions that a particular solution dominates The solution that has zero domination count occupies the first front. The domination count is reduced accordingly for the subsequent fronts. Deb and Jain (2013) suggested a reference point-based many objectives NSGA-II framework evolutionary algorithm, namely NSGA-III. Jain and Deb (2013) extended NSGA-III for solving generic constrained problems. The basic framework of NSGA-III is similar to that of NSGA-II, but with significant changes in its selection operator Jain and Deb 2013). The crowding distance operator is replaced with methods that maintain diversity among the population members by supplying and adaptively updating a number of well-spread reference points.
PAES uses its Pareto archive for evolution with a local search and starts from a population of one to archive from the previous solutions and a copy of that chromosome is mutated to form new candidate solutions after the evaluation of one chromosome (current solution) (Knowles and Corne 1999). The current and new candidate solutions are compared, if one solution dominates the other, the nondominated solution is absorbed into the archive. For neither solution dominates one another, the candidate solution is compared with a reference population of previously achieved non-dominated solutions. If the comparison fails to favor one candidate solution over the existing ones, the solution which resides in the least crowded region of space is selected.
In addition, SPEA2 uses a fine-grained fitness assignment strategy, a density estimation technique, and an enhanced archive truncation method (Zitzler et al. 2001) in its algorithm. In the fitness assignment state of SPEA2, a strength value is assigned to each member in the population and the archive (an external set) according to the number of solutions it dominates (Zitzler et al. 2001). The raw fitness value is calculated according to the strength value. The density estimation technique used is an adaption of the k-th nearest neighbor method (Dehnad 1987). The fitness value is calculated by adding both the raw fitness value and the density estimation value. SPEA2 has a fixed number of individuals in its archive form SPEA. The dominated individual may be copied to the archive if the number of non-dominated individuals is smaller than the fixed number of individuals in the archive. The aforementioned MOEAs are widely applied in different fields. Examples include an inventory control system (Chołodowicz and Orłowski 2017) and an operational planning problem (Alexandre et al., 2017). In both applications, NSGA-II and SPEA2 were used and compared. NSGA-II and PAES were used in a partial flexible job shop scheduling problem (Rabiee et al. 2012). In general, these MOEAs outperform each other in different fields.
The bibliometric network diagram for the three domains of drone deliveries, CDAs, and MOEAs is shown in Fig. 1. Vos Viewer was used to construct the diagram. To begin, we conduct a query in the Database Scopus using key terms ((''Transportation procurement'') AND (''drone delivery'' OR ''combinatorial double auction'' OR MOEA OR ''Multiple Objective Evolutionary Algorithm'')). However, only one item was found for this query, making the bibliometric network diagram insufficient. As a result, we modified the query to include the new keyword: (transport*) AND (''drone delivery'' OR ''combinatorial double auction'' OR MOEA OR ''Multiple Objective Evolutionary Algorithm''. The system searches for articles that are linked to ''transport*'' and meet any of the other criteria, such as ''drone delivery'' OR ''combinatorial double auction'' OR MOEA OR ''multiple objective evolutionary algorithm.'' There were 149 articles identified in the search results. Then we used the appropriate filter to narrow down the publications that were relevant. Only English papers were chosen (148 papers met this criterion), and then only articles were chosen (exclude conference papers, notes, reviews etc.). At this point, 88 papers fulfilled our journal article requirements. The third filter was used to hunt for relevant articles by looking at the article title and abstract. Finally, 57 papers were found to be eligible for bibliometric analysis. Appendix A contains the details of 57 papers.
As shown in the picture, CDA research is the most traditional and the first to occur (year 2016), followed by drone delivery research in 2019 and MOEAs regaining prominence in recent years, and finally freight transportation (yellow in color) is the focus in recent years. The three domains are interconnected and overlap as creative development, implying that there is tremendous opportunity for research projects combining the three in the context of transportation procurement, which is an economic subfield of freight transportation. Overall technical limitations in existing studies include that the CDA model is mostly single objective, and few scholars have investigated the biobjective model. Furthermore, previous studies had largely considered the CDA model as linear and solved the problem on a small scale. Due to these inadequacies, this research is conducted to examine a bi-objective model and upscale the model in a large scale setting. Large and complicated circumstances would make it difficult to solve a linear model. Therefore, in this research we also propose a decomposition model with two levels of variables, namely the assignment of packages and the proportion of loads, and then use MOEAs to discover solutions to the model.
Considering last-mile delivery with drone, CDA-based bi-objectives WDP model is unprecedented. MOEAs are applied to find the solutions for the WDP model to handle large and complex scenarios. The main reason of employing MOEAs is to tackle realistic, large scale scenarios in the transportation service procurement field. This is vital as last-mile delivery services has gained attention from logistic companies to cope with the increasing number of transactions in e-commerce, especially in C2C market. Some companies like Amazon and Alibaba have conducted run trials on using drone to deliver goods. This research could be a step forward to make last-mile drone delivery become a common choice in delivery services. The CDA-based WDP model is able to collect a pool of demands and supplies while MOEAs are able to solve this large scale problem in short timeframe. The convincing solutions indicated by the metaheuristic methods would be able to provide solutions for large and complex scenarios, therefore overcoming the inefficiency issue associated with large and complex WDP models in the marketplace.

Combinatorial double auction
Referring to the multi-objective combinatorial auction for transportation service procurement in Ignatius et al. (2014), this study extends the single-objective MILP model to a weighted bi-objective representation. The first objective is to maximize the profit of the ETM. The profit is calculated by summing all payments from the assigned customers' (buyers) packages minus the total cost of the assigned firms' (sellers) packages, as follows: The second objective is to maximize the marketplace fairness. This objective aims to maximize the numbers of assigned packages as a way to maximize the satisfaction of the participants in the marketplace. The second objective is formulated as follows: The weighted bi-objective model that maximizes both the profit and market fairness is a s ij x s ij À Ry s 0 8s; i; j ð Þ ð8Þ a s ij x s ij À a s ij y s ¼ 0 8s; i; j ð Þ ð9Þ x s ij ¼ R þ 8s; i; j ð Þ ð11Þ y c ¼ 0; 1 f g 8c ð12Þ where R in Eqs. (5) and (8) is an arbitrary large value. Note that obj 1 and obj 2 are used to normalize the importance of both objectives to avoid the condition whereby either one of the objectives is insignificant to affect the solution, e.g., a situation where the profit is 10000 while marketplace fairness is only 10. Inequality (4) ensures the supply of assigned packages is greater than the demand for all lanes. Inequality (5) to Inequality (7) impose constraints upon the volume of load pertaining to the range submitted by the firms while Eqs. (8) and (9) ensure that all the loads in the assigned customers' packages are accepted.

Decomposed model
There are two levels of variables in this problem. Letx be the vector of x c ij and x s ij whileỹ as the vector of y c and y s such that i 2 I and j 2 J for 8c 2 C and 8s 2 S. Letx andỹ be the vectors of decision variables determining the proportion of loads assigned and the 'winning' packages, respectively. An example ofỹ that contains a 6-package bid from carriers and a 6-package bid from a shipper is shown in Fig. 2. Figure 2 shows that c 2; c 5; c 6 f gpackage bid from customers (buyers) and s 1; s 2; s 6 f gpackage bid from firms (sellers) are assigned. Notice that variables in vector y are the first level binary independent variables, while variables in vectorx are the dependent variable corresponding to those inỹ. This study decomposes the MILP model in two stages: i. Assign the packages (values ofỹ) ii. Find the proportion of loads for the assigned packages (values ofx) The decomposed model becomes: Inequality (16) ensures that vectorỹ is a feasible solution (i.e., supply is greater than demand).
Note that v 2 S given that y v ¼ 1, and k 2 C given that is the revenue from the assigned firms' package bid, y v and X k2 C\y k ¼1 ð Þ is the cost from the assigned customers' package bid, y k . Inequality (20) ensures the supply is greater than demand while Eq. (21) imposes constraints upon the ranges ofx so that it is within the ranges submitted by firms. Notice that this study assumes all customers want all their submitted volumes of loads to be delivered if they are assigned (i.e., the proportion of loads for customers are 1 in all lanes).

MOEAs
Four MOEAs are engaged in stage 2 of the proposed decomposed model. A general framework that depicts the fitness selection strategy of each MOEA is shown in Fig. 3. In this paper, MOEAs evolveỹ only with respective selection features and methods. Objective 2 scoring as shown in Eq. (15) can be computed directly onỹ. However, to get objective 1 scoring,x need to be optimized (stage 2) during fitness evaluation as shown in Fig. 4. Generally,x is optimized by the following steps: 1. If stage 1 is feasible (supply is greater or equal demand in all lanes),ỹ is forwarded to stage 2; else continue to evolveỹ 2. In stage 2, setx to the lower bound in a submitted bid 3. In each lane, the prices submitted in bids are sorted in an ascending order 4. In each lane, if the sum of the minimum supply is lower than the demand,x is adjusted within the upper bound according to the sorted prices. 5. Afterx is optimized, the profit is calculated The pseudo-code for optimizingx is:

Data generation and results analysis
This section will describe the generation of instances, the evaluation indicators of Pareto solutions, the results of instances and the analysis of results.   i.e., the demand of shipper m for all m 2 M. Each a m ij is generated with a discrete uniform distribution in the range [1000, 10000], subject to a uniformly generated a 50 and i 6 ¼ j. Notice that m at the top right of the matrix denotes shipper m. Then, A m is separated to several packages, 7 5 s such that s 2 S randomly uses an exponential distribution for each m 2 M. Note that a non-zero a m ij is selected randomly and store in A s . In every selection of a m ij , the selected a m ij is stored to the next package subject to probability of an exponential distribution. By setting k ¼ 4, the probability is calculated by 1 À e Àkt , where t is the number of lanes stored in the current package.

Generation of instances
Then, according to A m , the supply of carrier n, . . b ll 2 6 6 6 4 3 7 7 7 5 n is generated for all n 2 N. As this is a multi-shipper and multi-carrier problem, as long as lane i; j ð Þ is opened, each b n ij is generated with a discrete uniform distribution in the range [1000, 10000], subject to a uniformly generated a 50 and i 6 ¼ j. Again, n at the top of the matrix denotes the index of carrier n. B n is separated into several packages, 7 5 c such that c 2 C randomly uses an exponential distribution for each n 2 N.
Assuming that the distance is longer when i À j j j is larger,p s ij and p c ij are then generated using a normal distribution with mean i À j j j15 and i À j j j10 with a standard deviation of 2 for the corresponding lanes and packages, respectively. The lower and upper bounds that carriers could accept for the corresponding lanes and packages, i.e., L c ij and U c ij , are set to 40 and 100% of the quantity of bids submitted, respectively. A summary of the CDA parameters are shown in Table 1.
This study generates the transportation data for different numbers of lanes with a computer of 1.80 GHz processor and 8 GB RAM, and with MATLAB and Java software packages. Each simulated data set is fed into the weighted bi-objective CDA model, which is solved using the weighted MILP method. The simulated data set is fed into the decomposed model and solved using NSGA-II, NSGA-III, PAES, and SPEA2. The results of different weighted objectives solved by MILP are shown in Table 2. The settings of the MOEAs are presented in Table 3 (based on settings used in the corresponding published studies).

Results of instances and analysis
This study presents the results of different MOEAs in Fig. 5a-c, with respect to the number of lanes. Notice that for each MOEA with the same population size, the results of 10 runs and with different numbers of fitness evaluations (NFE) (i.e., 1000, 5000 and 10,000) are shown. The approximated time required to run 11 different weight combinations for the bi-objective model using MILP are shown in Table 4. The durations needed to run each MOEA are recorded in Table 4.

Analysis of results
This study categorizes the scenarios into three categories, i.e., a small scale case of 12 lanes with a total of 105 package bids; a medium scale case of 42 lanes with a total of 593 package bids; and a large scale case of 90 lanes with a total of 1795 package bids. Notice that the MILP results are treated as the Pareto optimal solutions. For the small  scale case, the time needed for 1 simulation to run 10,000 NFE for each MOEA is approximately the same as that required for 11 points from MILP with different weight combinations. Except PAES, the three MOEAs are able to produce more than 11 points on the Pareto optimal front. As such, NSGA-II, NSGA-III and SPEA2 are more efficient that MILP in solving the small scale case. For the medium scale case, the time and solution set from NSGA-II and SPEA2 are similar. While NSGA-III takes a longer time (few more seconds), its solution set is more diversified, and there are some solutions that fall at the Pareto optimal front. Comparing with 11 MILP runs MILP (more than 5 h), the time needed for 1 run of NSGA-II, NSGA-III, and SPEA2 is extremely fast (not more than 30 s). Furthermore, the solution sets from MOEAs except PAES are close to the Pareto optimal front.
For the large scale case, the outcome is similar to those of the medium scale case, except that the distance of the solution sets produced by NSGA-II, NSGA-III and SPEA2 is further from the Pareto optimal front. However, considering the time needed to run 11 MILP runs (more than 10 h), these MOEAs are more efficient in producing the feasible results.

Evaluation indicators of Pareto
This study further investigates the performance of these MOEAs by randomly run 30 seeds for each MOEAs with respective scenarios and calculate some performance metrics such as generational distance (GD), hyper volume (HV) and spacing. GD measures the closeness of the solution set to the PF (Van Veldhuizen 1999;Van Veldhuizen and Lamont 2000;Yen and He 2014). It reflects the convergence of the solutions. The nearer the GD value to 0, the closer the solutions are to the PF. Let say our solution set, S contain P number of solutions such as S ¼ z 1 ; z 2 ; . . .; z P f g , this study calculates the GD as   On the other hand, HV consider both closeness of the solution set with PF and diversity of the solution set (Zitzler and Thiele 1998;Van Veldhuizen 1999;Yen and He 2014). It calculates the area of the enclosed space of the solution set with a reference point. For example, referring to the bi-objective problem, the equation to calculate the hyperarea of the solution set and origin (as reference point) is where a z p À Á is the rectangle area bounded by origin and f z p À Á for bi-objective problem. Considering the difficulty to choose reference point to calculate hyperarea, Auger et al. (2009) suggested an algorithm that helps us to decide reference point in a more proper way. Using true PF as a reference, David and Veldhuizen (1999) proposed a hyperarea ratio metric defined as In this metric, a high value means the solution set is close to the reference point and have good diversity.
Lastly, spacing measure the distribution of the solution set (Schott 1995;Van Veldhuizen and Lamont 2000;Yen and He 2014). The equation is given as below where d p is the distance in objective space of the solution p and the nearest point in PF, r such as r 2 PF and P is the number of individuals in the solution set. A value of 0 for this metric means that all members in solution set are equidistantly spaced. Figure 6a-c shows the graphs of GD, HV and spacing metrics for each scenarios.
From the graph of GD metrics, the distance between the solution set and reference points from NSGA-II and SPEA2 are similar throughout the evolution for all cases. Meanwhile, NSGA-III has similar ''closeness'' to PF with NSGA-II and SPEA2 after 10,000 NFE for small scale case but is slightly closer to reference points for medium and large scale cases. Besides that, it needs fewer NFE to get solution set that has similar GD value compared to NSGA-II and SPEA2. From the HV metrics, this study found that NSGA-II and SPEA2 are again similar throughout the evolution while NSGA-III generally has higher HV value than both NSGA-II and SPEA2. NSGA-II and SPEA2 can obtain similar ''diversity'' in small scale case after 10,000 NFE. This shows that the diversity of NSGA-III is better than NSGA-II and SPEA2. From approximation set in Fig. 5a-c, it is clear that except small scale case, NSGA-III is more diversified than NSGA-II and SPEA2. From the spacing metrics, NSGA-III tend to be more evenly spread than NSGA-II and SPEA2 at early evolution stage (NFE around 1000) but after 10,000 NFE, NSGA-II and SPEA2 are more evenly spread than NSGA-III.
Generally, a summary of the analysis is as follows: • NSGA-II and SPEA2 produce similar solutions sets with similar computational time. • NSGA-III is more diversified, and its solution set is generally closer to the Pareto front than those of NSGA-II and SPEA2, but it requires a slightly longer computational time. • PAES produces poor results.

Conclusions
The CDA-based WDP model proposed by Motlagh et al. (2010) has been extended to a bi-objective representation, and applied to last-mile drone delivery services. The resulting MILP formulation has been decomposed and solved using MOEAs. The main reason of employing MOEAs is to tackle realistic, large-scale scenarios in the transportation service procurement field. This is vital as last-mile delivery services has gained attention from logistic firms to cope with the increasing number of transactions in e-commerce, especially in the consumer-toconsumer market. Some firms like Amazon and Alibaba have conducted run trials on using drone to deliver goods. Our study is a step forward to make last-mile drone delivery become a common choice in delivery services. In practices, CDA-based WDP framework is able to collect a pool of demands and supplies while MOEAs are able to solve this large-scale problem in a short timeframe. As last-mile delivery business is competing speed of delivery to the end user, failure to meet customers' expectation might cause business loss.
The contributions of this study are (1) propose a new framework, that is to use decomposed bi-objectives CDAbased WDP model in last-mile drone delivery system; (2) show the need to implement metaheuristic methods to the proposed framework for better efficiency; (3) compare the results for four common MOEAs and analyze the results with common metrics. In general, NSGA-II, NSGA-III, and SPEA2 are able to produce feasible solutions very rapidly. The accuracy rates deteriorate slightly in large scale cases based on the default MOEA settings. The performance could be improved by tuning the parameters and evolution operators such as crossover and mutation operators in the respective algorithm according to a specific scenario for further research.
Although this study has reached its aims, there is an unavoidable limitation regarding the data simulation procedures. As there are no open sources data that matches with the scenarios of CDA-based WDP in the proposed framework, data is simulated randomly using common distributions such as uniform distribution, normal distribution and exponential distribution. It would be better representing real situation if the data is the real data or it is generated according to the real distribution.
Future research might concentrate on increasing the 'intelligence' of knowledge-based operators in order to achieve a more accurate solution set in less time. This may be done via operators like the initial population generating function, crossover function, and mutation function, among others. Not to mention performing a sensitivity study on the parameters for the algorithms that are used for tuning. Besides, the number of objectives can be further increased, e.g., marketplace confidence and reputation of the participants. This is because decision-makers generally wish to factor in the possibility of revocation. Furthermore, in the representation model accounting for the robustness feature, a fuzzy environment may be applied to the pricing or the load volume.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Human or Animal Rights This article does not contain any studies with human participants or animals performed by any of the authors.
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://creativecommons. org/licenses/by/4.0/.