Developing a Collaborative Planning Framework for Sustainable Transportation

Currently, as being the highest petroleum consuming sector in the world, transportation significantly contributes to the total greenhouse gas emissions in the world. Road transportation not only is responsible for approximately 20% of the total emissions of carbon dioxide in the EU and in the US but also has a steadily increasing trend in contributing to global warming. Initiatives undertaken by authorities, such as Emission cap and trade in the EU, limit the emissions resulted from the actions of the companies and also give economic incentives to companies to reduce their emissions. However, in logistics systems with multiple entities, it is difficult to assess the responsibilities of the companies both in terms of costs and emissions. In this study, we consider a delivery network with multiple customers served by a single carrier, which executes a delivery plan with the minimum transportation cost, and allocate the resulting costs and the emissions among the customers in a fair manner. We develop allocation mechanisms for both costs and emissions. In order to develop a mechanism that provides further reduction of the emissions, we study a setting where the carrier takes the responsibility of the emissions and reflects the resulting inefficiencies while charging the customers.


Introduction and Literature Review
Transportation sector is facing formidable challenges and recently environmental issues lie in the heart of these challenges.Currently, as being the highest petroleum consuming sector in the world, transportation significantly contributes to the total greenhouse gas emissions and hence to the global warming.In the EU, road transportation is responsible for approximately 20% of the total emissions of carbon dioxide.According to European Commission Directorate-General for Energy and Transport [1], in the EU in no other sector has the growth rate of greenhouse gas (GHG) emissions been as high as in transport, which amounted to 23% in 2010 [2] compared to 1990 levels.In the US, transportation is responsible for approximately 28% of the GHG emissions, second to the electricity sector, and compared to 1990 levels it has been increased by 18% [3].This increasing trend in both the EU and the US is mainly due to the increased demand for travel of individuals and increased freight transportation activities.
On the other hand, road transportation can be considered as the most important mode for passenger and freight transportation.In 2010, total goods transport activities in the EU27 (only intra-EU) have amounted to 3,831 billion tonne kilometers.Road freight transport accounts for 45.8% of the total freight volume in the EU [4].In the US, road transportation is responsible for hauling 70% of the total freight in 2012.Even though road transportation plays such a crucial role in passenger and freight transportation, it is also the largest contributor to global warming.Over the next 50 years, the GHG emissions of the transportation sector could grow 80% compared to current levels due to increased transportation activities [5].
Recently, there have been several initiatives to reduce the environmental damage caused by transportation activities.According to European Commission Directorate-General for Energy and Transport [1], the environmental suitability of transport activities is a major concern and immediate action should be taken to mitigate the negative effects of transportation on the environment.According to European 2 Mathematical Problems in Engineering Commission Directorate-General for Energy and Transport [1], the EU has recently adopted a climate and energy package that sets a target of reducing GHG emission in the EU by 20% with respect to 1990.The US transportation sector also considers three potential approaches that could reduce GHG emissions as well as fuel consumption: (i) improvements in vehicle technology, (ii) switching to environmentally friendly fuels, and (iii) transportation/travel demand management (TDM), which corresponds to better utilization of the transportation resources [5].
Based on these facts and figures, it is clear that governments and international organizations should play an active role in achieving sustainable transportation.Several governments have already passed legislations to support sustainable transportation including the EU's Transport Protocol of the Alpine Convention [6].An initiative for reducing GHG emissions is to limit the amount of pollutant emitted by the direct and the indirect actions of a company.This is referred to as "emission cap and trade, " which allows controling pollution using economic means.European Union Emission Trading Scheme [7] is an example to these cap and trade initiatives.In this context, companies have certain limits either allocated or sold by the central authority.Companies must limit the emissions resulted by their actions to their thresholds, which are set by the permits they have.If the thresholds are not sufficient, then companies must acquire additional permits from other companies whose emissions are below their limit.The allowance of acquiring additional permits is referred to as "trade, " which gives economic incentives to the companies to limit their GHG emissions.However, it is quite difficult to assess the companies' emission responsibilities resulted from their indirect actions.
Having this in mind, we develop mechanisms for the fair assessment of the emission responsibilities of the customers in a logistics system.In this context, a single logistics provider (carrier) may serve several customers along a delivery route in order to minimize the transportation costs.Allocating cost among the customers is required to cover the expenses plus the mark-up of the carrier and allocating emissions is required for emission cap and trade purposes discussed above.The GHG emissions are allocated among the customers in units of grams of CO 2 (gCO 2 ).Even though the transportation cost and the emissions are correlated as they both depend on the fuel consumption of vehicle, the relationship is not trivial as there are several other factors contributing to each of these values.For instance, the operating cost of a truck is estimated as 49 cents per mile on average; however, only 24.2 cents of this cost is resulted from the fuel consumption and the rest is due to maintenance, repairs, and depreciation [8].On the other hand, the emissions are largely affected by the payload of the vehicle, as it varies considerably along a delivery route.The cost of emissions depends on the final emission responsibilities of the customers, which determine the total number of permits required by the company as well as the trading cost of a permit on the market.The cost of the permit on the market varies considerably mostly due to the available number of permits.The standard reporting unit is tones of carbon dioxide equivalent (tCO 2e) and there exist several studies to forecast the price per ton of CO 2 [9].In this study, we assume the price as $42 per ton of CO 2 even though this value might fluctuate over time.Note that both per mile cost of transportation and per ton cost emission represent just two parameter values and do not affect the solution methodology at all.On the other hand, the allocation of responsibility is important and affects the cost/emission allocation mechanisms.In the first setting we consider the customers are responsible for both the transportation cost and the emissions directly and the carrier minimizes the transportation cost to serve the customers and allocates both the cost of transportation and the total emissions to the customers individually.In the second setting, the emission responsibility is diverted to the carrier and the carrier decides the minimum total cost plan to serve the customers and allocates a single cost figure, which includes the transportation cost and the emission cost, to the customers.Our motivation for considering the second setting is that it might favor social welfare as the total overall emission might be reduced when the carrier's objective includes this as well.
Unfortunately, neither allocating costs nor allocating emissions are a simple task.The main reason is even though the total cost and the total emission can be easily calculated, it is not trivial to determine which customer is responsible for how much of these totals.We illustrate this on a simple example.Consider a delivery route of a carrier (Figure 1) serving three customer locations.Suppose that the carrier uses a truck with capacity .Suppose that Customers  and  are requesting deliveries with weights equal to 0.1 and Customer  is requesting a delivery with a weight equal to 0.8.Suppose that the cost of traveling between any two locations except from  and  is 1 and negligible between  and .The optimal delivery route is trivial as all demand can be consolidated in one truck (Figure 1).Note that this route also results in the minimum GHG emission, since the largest bulk delivery spends as less time as possible on the vehicle.If we were to allocate both the costs and the emissions proportional to the customers' distance to the depot, then all customers will have equal allocations in both costs and emissions.It is clear that this allocation is not acceptable since Customer  should be allocated a higher emission (due to its load weight) and Customer  should be allocated a higher cost value (due to its location).If we were to allocate both the costs and the emissions proportional to their individual solution values, respectively, then Customers  and  will be allocated the same cost and emission.Again, this is not an acceptable allocation as the cost and the emission contribution of Customer  on this delivery route is less than Customer  and hence should be allocated less in each category (as Customer  has a synergy with Customer ).Finally, if we were to reverse the route, then the contributions of Customer  seem relatively lower as the load delivered to Customer  now spends less time on the vehicle.However, delivery route determination decision is an exogenous decision of the carrier.Hence, the allocation decisions should be independent of the routing decisions.These issues remain even when the carrier takes the responsibility of the emissions and only allocates costs to the customers.Before discussing the related works in the literature, we present a list of desirable properties for an allocation method.First, an allocation method should allocate the entire costs/emissions among the participants, which is referred to as "budget balancedness." Second, each individual should be allocated at most its individual solution value, which is called "individual rationality." Third, if we take this property one step further, the same principle should apply to each subset (coalition) of the members, which is referred to as "stability." Fourth, the carrier's routing decisions should not have an effect on the cost/emission allocations.Otherwise, it would not be proper to charge a customer a higher cost/emission value just because the routing decision did not favor that particular customer.Finally, an allocation mechanism should treat each participant "fairly".In our context, a fair allocation method should satisfy budget balanced and stability conditions (referred to as "core allocation" in the literature) and should be independent of the routing decision of the carrier.Faigle et al. [10], Pál and Tardos [11], Engevall et al. [12], and Geomans and Skutella [13] are some examples from the literature which adopt the same fairness definition.In most cases, a core allocation may not exist; hence, our objective is to identify a good approximation to the core.
Allocating the costs/benefits in a logistics system have been studied in the literature under various settings.Tamir [14], Potters et al. [15], and Kuipers [16] consider a traveling salesperson problem (TSP) and discuss the existence of core allocations.Faigle et al. [10] propose a moat-packing based allocation method for TSP to obtain fair allocations.On the other hand, Göthe-Lundgren et al. [17] and Engevall et al. [12] analyze a vehicle routing problem (VRP) setting from a more practical perspective and propose a computational approach to compute nucleolus.Finally, Özener and Ergun [18] consider a vendor managed inventory setting and propose methods to compute cost-to-serve values under this setting.
Recent papers study the carbon emission and transportation relationship from different perspectives.Bektas and Laporte [19] study the pollution-routing problem, an extension to VRP that takes other factors besides travel distance into account such as greenhouse emissions, fuel, and travel times.Pitera et al. [20] consider a case study on an urban transportation system and propose local search algorithm for the underlying emissions minimization problem.Suzuki [21] proposes a solution approach for a truckrouting problem where the objective is to minimize the fuel consumption and pollutants emission.Xiao et al. [22] study the capacitated VRP with a load-dependent fuel consumption rate and propose a solution approach to this problem using a simulated annealing algorithm with a hybrid exchange rule.Gajanand and Narendran [23] formulate a multipleroute VRP problem to minimize the fuel consumption and show that taking alternative routes may offer reduction in fuel consumption.Xue et al. [24] develop a model of the relationship between traffic flow and vehicle emission under road capacity constraints.Wang et al. [25] focus on the carbon emission rights and propose an auction mechanism which converges to a unique equilibrium.
Finally, Sichwardt [26], Leenders [27], and Naber [28] consider CO 2 allocation among the customers under a VRP setting.Sichwardt [26] analyzes the problem based on the empirical data of a Swedish medium-sized transportation company.The author assumes that the transportation decisions are given based on a combination of the sweep and savings algorithm.After the CO 2 emission calculations on a specific route, the emissions are allocated to the customers based on their distance to depot.Leenders [27] assumes a set of properties as the fairness criteria: (i) individual rationality, (ii) marginality, (iii) efficiency, and (iv) no negative allocation.Based on a comparative analysis, the selected allocation method is the one based on the individual rationality of the shipments.Naber [28] also discusses the vagueness of the term fairness and comes up with a long list of fairness concepts including stability and budget balancedness.The author compares methods such as Shapley Value and nucleolus as well as a proportional-based method.However, the analysis is on a per-route basis and hence the carrier's routing decision may have a huge impact on the final allocation of the emissions.More importantly, as the computations are on an individual route level, the number of customers, which significantly affects the computational complexity of the generic allocation methods, is less than 15.Hence, other than the proportional based method, these approaches are not scalable for real-life sized problems.
Contrary to the studies reviewed above, we propose methods to simultaneously allocate both the transportation costs and the emissions among the customers.We develop an allocation mechanism based on duality and propose an approximation for the Shapley Value.We refer the reader to Owen [29], Kalai and Zemel [30], Pál and Tardos [11], and Geomans and Skutella [13] as examples of the application of the duality-based allocation methods.Our fairness criteria are a more restrictive one compared to the studies discussed above.We computationally show that our methods outperform the proportional methods proposed in those studies.Finally, to the best of our knowledge, this is the first study that considers the setting where the carrier takes over the emission responsibility and allocates a joint cost to the customers.
The remainder of the paper is organized as follows.In Section 2, we provide a formal definition of the problem and list our assumptions.In Section 3, we present solution approaches for the carrier's routing problem under two different settings.In Section 4, we present our proposed solution approaches, dual LP based and Shapley Value, respectively.In Section 5, we computationally demonstrate how our proposed mechanisms perform in comparison with the other allocation schemes.Concluding remarks are provided in Section 6.

Problem Definition
In this section, we provide a formal definition of the problem, explain our emission calculation approach in detail, and finally list our assumptions.
We consider a set of customers  = {1, . . ., } that are served by a common carrier from a single depot location.Including the depot location, the vertex set is  = {0, 1, . . ., }, where 0 represents the depot and edge set is  = {(, ) : ,  ∈ }.The distance along an edge (, ) is denoted by   and these figures are assumed to satisfy the triangle inequality.Note that the proposed solution approaches are valid under an asymmetric distance setting and can be used with minor modifications.The carrier serves customers via a set of delivery routes and each route is assigned to one vehicle from a finite number, , of homogenous vehicles with capacity .Customer  is to be delivered a shipment with a weight denoted by   ,   ≤  for all .Note that a feasible delivery route corresponds to an ordered customer list with total shipment weight less than equal to the truck capacity.
In this context, delivery routes executed by the carrier determine the total transportation costs and the total emissions.It is important to note that the delivery routes determined by the carrier are not necessarily the optimal set of delivery routes to serve the customers.This is the case in many real life delivery systems; nevertheless, the cost to be allocated among the customers corresponds to the total cost of the selected delivery routes.We assume that the carrier is operating a commercial delivery truck with an operating cost of 49 cents per mile or equivalently 30 cents per kilometer [8].Note that these values are just parameters and do not actually affect neither the proposed methodology nor the comparative results.We denote the transportation cost along an edge (, ) by   .Finally, there are two distinct approaches for calculating CO 2 emissions of transportation: energy-based approach and activity-based approach.Similar to the studies discussed above, we employ an activity-based approach.However, the difference of our approach is that we use a piecewise linear function with several break points rather than a linear function in order to better approximate the CO 2 emission values.The general formula is as follows: CO 2 -emission factor per ton-km. (1) Note that the transport volume is the payload of the truck covering that particular distance.The calculation of the average CO 2 -emission factor per ton-km is done assuming a 40-44 ton truck and emission factors presented in Table 1, which is acquired from McKinnon and Piecyk [31].
Based on the data above, the emission value along an edge (, ) with a payload of , which is denoted by   (), is calculated by identifying the average CO 2 -emission factor per ton-km with respect to the payload along the edge and multiplying this factor by the payload volume and the distance.We assume that the truck itself weighs around 14-15 tons and hence the payload volume varies between 0 and 30 excluding the truck's own weight.
In this paper, we consider two different settings based on the emission responsibility (hence the objective function) of the carrier.Our motivation for analyzing two different settings is to assess the performance of these distinct policies from the perspective of fairness in allocating costs and emissions.The settings are as follows.
(i) First setting: the emission responsibility is on the customers and thus the objective of the carrier is to minimize the total transportation costs.
(ii) Second setting: the emission responsibility is on the carrier and thus the objective of the carrier is to minimize the total cost including the transportation costs and the emission costs.
In the first setting, the carrier does not have an incentive to execute low-emission routes and hence treats the emission caused by the delivery routes as a byproduct of the transportation activities.As a result, after the execution of the routes both the total transportation costs and the total emissions are to be allocated among all the customers.
In the second setting, the carrier wants to execute costeffective and low-emission routes.As a result, after the execution of the routes the carrier has to allocate a joint cost among the customers in a fair manner.For this setting, we assume a unit cost of $42 per ton of CO 2 emission and one permit is equivalent to one ton of carbon dioxide.Hence, the emission cost is a step function in terms of tons of emission.However, as the carrier might probably have other activities besides serving these customers and will pay an overall emission cost annually for all its activities, it will be an overestimate to use a step function with intervals of tons.Instead, we assume a linear relationship while estimating the emission cost of the delivery routes of the carrier.

Solution Method for the Carrier's Routing Problem
In this section, we present a solution method that determines the best routing of shipments from the perspective of the carrier.Note that, in the first setting, the carrier's objective is to minimize the cost of transportation.Hence, the routing decision of the carrier does not take into account the emission caused by the resulting delivery routes.The underlying routing problem is a well-studied problem in the literature, VRP.
In order to solve this VRP, we first iteratively generate clusters using a modified K-means clustering approach, and then based on these clusters we generate a number of feasible delivery patterns.Next, we compute the transportation cost and CO 2 emission of the generated patterns, and finally we solve a set partition problem (SPP) to select the delivery patterns to serve the customers with the minimum cost possible.As the number of customers along a route is limited due to the truck capacity restrictions, we are able to solve these TSPs using exact methods and commercial solvers.
Our motivation for choosing this particular solution method among all possible VRP solution methodologies is that we are able to generate a large set of delivery routes, all of which will be used in our allocation methodology.We use a variant of the K-means clustering algorithm, which is discussed in [32,33].In this approach, we partition the customers into a certain number of clusters.In our problem, this number would be equal to the number of available trucks, which is denoted by .After initiating the first  clusters randomly, we implement a variant of the set partition problem (VSPP) to assign customers to these cluster bases.VSPP ensures that the total weight of the shipments assigned to a cluster does not exceed the truck capacity.The details of the overall solution algorithm are as Procedure 1.Note that  represents the iteration number of the outer loop and  represents the iteration number in the inner loop in generating the delivery patterns.Before presenting the formulation of , we introduce the notation used in the formulation.Let   be the binary variable indicating whether a pattern is selected in the final solution or not.Let   be the set of delivery patterns that visit customer .The formulation of  is as follows: The objective is to minimize the total transportation costs.Constraints (3) ensure that each customer is served by exactly one delivery route.Constraint (4) ensures that the number of delivery routes executed should be less than or equal to the total number of vehicles.The last constraints are the binary restrictions of the variables.
Solving the carrier's routing problem would determine the minimum cost solution of delivering to all the customers.Note that, in applying Procedure 1, we calculate not only the transportation costs but also the CO 2 emissions of the delivery routes.Here,  *  corresponds to the optimal objective function value of , which may not correspond the optimal solution to the problem as Procedure 1 does not guarantee optimality.We assume that  *  denotes the total cost to be allocated among the customers.We also need to compute the total emissions to be allocated among the customers.Let  *  denote that value, which is equal to ∑ ∈    *  .Finally, we are able to generate a large set of feasible delivery patterns, not only the ones executed in the final solution.This is very crucial because, for our proposed allocation methods, we require a range of feasible delivery routes, especially to determine a fair allocation method independent of the routing solution preferred and executed by the carrier.
In the second setting, the carrier's objective is to minimize the total cost of the delivery routes, which includes the transportation costs and the emission costs.Under this setting, we modify Procedure 1, specifically the step "solving the TSPs to calculate the cost and the emission of the route." Unfortunately, this modification is a major one.In the original procedure, we solve a mathematical model for TSP with regular assignment constraints as well as subtour elimination constraints.However, when the emission costs are also included in the objective function and since the emission costs depend on the payload of the vehicle on a certain link, we need to keep track of the route traveled and hence the payload on the vehicle at any given time.Therefore, we need an entirely different model to handle this situation.The problem resembles the delivery man problem (DMP) studied by Fischetti et al. [34].In addition to this modification, we calculate the emission values using a piecewise linear approximation with several break points.Hence, we need a further modification of the flowbased model proposed by Fischetti et al. [34].Unfortunately incorporating the piecewise approximation into this flow formulation makes the model too time consuming to be used for large number of delivery patterns.Therefore, we use a step function for assessing the emission cost in this flow-based model's objective function and calculate the emission cost by using our original piecewise approximation after solving the model.The formulation of our flow-based TSP model, FTSP, is as follows: ∈ {0, 1} , ∀, ∈ Ñ, ∀ ∈ {0, 1, . . ., 30} .
In the formulation above Ñ represents the set of customers to be visited by the given TSP pattern  plus the carrier's depot location.The binary variable   denotes whether edge (, ) is selected in the final TSP tour or not.The integer variable   represents the flow volume (weight) along edge (, ).Note that as   values are in kilograms, the flow values along the edges are represented in integer multiples of kilograms and also the weight of the truck is not included in this flow volume.The variable   represents the cost of emission along the edge (, ).Finally, the binary variable   denotes whether on edge (, ) the flow volume falls within an interval between  − 1 and  or not.
The objective function is to minimize the total transportation costs and the emission costs.Constraints (7) and (8) ensure that each customer is visited only once by the delivery pattern.Constraints ( 9), (10), and (11) ensure that the flow conservation within the network is achieved while satisfying the volume delivery requirements of the customers.In order to represent a positive flow on the vehicle along the selected edges of the pattern, we increase the outgoing flow value from the depot by one kilogram, which is negligible in emission calculations.Constraints (10) also guarantee that there will be no subtours in a feasible solution.Suppose that there is a subtour that does not include the depot (otherwise there is only one grand tour, and hence no subtours).For all the nodes in that particular subtour, we sum up the corresponding constraints of type (10) and this would yield a contradiction.Hence, any feasible solution to FTSP does not have subtours.Constraints (12) guarantee that if an edge has a positive flow then that edge has to be covered by the pattern.Constraints ( 13), ( 14), (15), and ( 16) make sure that the flow on an edge falls within a certain interval between 0 and 30 and the emission cost along an edge is computed according to that specific interval.In Constraints (16), the constant, 42, represents unit cost of emission per ton of CO 2 and   ( − 0.5) represents the total amount of emission in tons along edge (, ) if this edge is covered with a payload of  − 0.5 tons.Note that   ( − 0.5) values are parameters that can be computed in a preprocessing step.Also,   ( − 0.5) is an approximation as even though the real flow value, (  − 1)/1000, falls within the interval [ − 1, ], it may not be exactly equal to  − 0.5.The last four constraints are the binary, integer, and nonnegativity restrictions of the variables.Recall that, after solving the model above, we calculate the total emission cost using our original piecewise approximation and update the total cost of the TSP pattern accordingly.The rest follows by replacing   values with   values and solving a modified  to obtain  *  = ∑ ∈    *  .

Cost/Emission Allocation Methods
In this section, we describe our proposed allocation methods both for cost and emission allocation.Before describing each method, note that the cooperative game under a VRP setting may have an empty core, which follows directly from the proof on the existence of the core allocation under a TSP setting, and hence both the cost and the emission allocation problems are not guaranteed to have a budget balanced and stable cost allocation.Having this in mind, our objective is to find the most fair allocation method, which we define as being exactly budget balanced and as stable as possible (also note that stability implies individual rationality).Hence, for all the allocation methods, we calculate the percent instability values for all possible groups of customers and decide which allocation method yields better solutions based on these instability values.Finally, ease of computation is another criterion in choosing the best allocation method as the applicability of the proposed methods directly depends on their computational times and the related works in the literature mostly focus on the theoretical contribution and ignore this aspect of the problem.

An Allocation Method Based on the Distance to Depot.
First, we describe an allocation method, which allocates the costs and the emissions proportional to the customers' distances to depot.It is natural to assume that the contribution of a customer to both transportation cost and CO 2 emission of a deliver route is correlated to the customer's distance to the depot, even though this correlation might be stronger in cost contribution than that of emission contribution as the weight of the shipment is ignored.The obvious advantages of such a cost/emission allocation method are as follows: (i) it is a simple/fast allocation method that requires no complex calculations.(ii) it provides individually rational allocations, and (iii) obtained allocations do not depend on the routing solution of the carrier.The disadvantage of this method is that it completely ignores the synergies among the customers and it is very unlikely to obtain approximate core, hence fair, allocations using this proportional method.A particular disadvantage of this method in allocating CO 2 emission is that this method completely ignores the shipment weights of the customer, which is a very important factor in determining the CO 2 emissions.We compute   , the cost allocated to customer , and   , the emission allocated to customer , as Procedure 2.
Note that for our second setting, we replace  *  value with  *  value in Procedure 2 to obtain cost allocations.

An Allocation Method Based on the Individual Solutions.
Second, we describe an allocation method, which allocates the costs and the emissions proportional to the solutions where the customers are served individually by the carrier.It is clear that the cost allocations will be exactly the same as the proportional allocation based on the distance to the depot.However, the emission allocation will be different as the weight of the shipments becomes a factor in this allocation method.Hence, we expect this allocation method to perform similar to the one discussed above except that the emission allocations might be better in terms of fairness criteria.The details of the procedure are as Procedure 3. Again, for our second setting, we replace  *  value with  *  value in Procedure 3 to obtain cost allocations.

An Allocation Method
Based on Duality.Third, we introduce our first proposed method for allocating the costs and the emissions, which is based on duality.The relationship between core allocation and LP duality is a well-established one since Owen [29].In this allocation method, we use the dual of the LP relaxation of  to obtain cost allocations.Note that, from the optimal values of the corresponding dual variables, we obtain weights and by multiplying with  *  and  *  we obtain the actual allocations.
There are several advantages of the duality-based allocation method: (i) it is a simple/fast allocation method that only requires an LP solution, (ii) obtained allocations do not depend on the routing solution of the carrier, (iii) and it assesses the synergies among the customers better than other allocation methods discussed and hence is more likely to yield fair allocations.
The LP relaxation of  is as follows: Let   be the dual variables associated with constraints (22), let V be the dual variable associated with constraint (23), let   be the dual variables associated with constraints (24), and finally let   be the set of customers visited by pattern , and then the dual of the  LP is as follows: Solving the dual LP above in polynomial time yields the optimal values for the dual variables   , which are the variables to be used in determining the cost allocations.Note that we need to solve another, modified dual LP in order to obtain dual variables for allocating the emissions.In that case, we just replace   values with   values to obtain ỹ *  .The details of the procedure are as Procedure 4.
To obtain duality-based cost allocations for our second setting, we replace  *  value with  *  value and   values with   values.

An Approximation for the Shapley
Value.Finally, we introduce an allocation method based on an approximation for the Shapley Value.The Shapley Value is a generic allocation method introduced by Shapley [35].In a nutshell, Shapley Value is the weighted average of the marginal contribution of each member to each subset of the collaboration.Let   () be the marginal cost of adding customer  to the subset  and let   () be the marginal emission of adding customer  to the subset , and then the Shapley Value calculations are as follows: The advantages of the Shapley Value allocation method are as follows: (i) obtained allocations do not depend on the For  = 1, . . .,  Let   be the set of the customers including the nearest 5 customers and customer

End For
Procedure 5: Approximation for the Shapley Value.
routing solution of the carrier, (ii) it assesses the synergies among the customers, and (iii) it satisfies additional properties such as equal treatment of equals, additivity, and dummy properties.A major disadvantage of the Shapley Value is that it requires an exponential effort to compute the allocation unless an implicit method is developed (as the number of subsets grows exponentially).Hence, in its current form Shapley Value is not a suitable method for the large instances of the problem.Therefore, we propose an approximation for the Shapley Value, in which we calculate the weighted average of the marginal contribution of each customer to the subsets that include only the nearest 5 customers.By doing that, we limit the number of subsets and hence limit the computational effort that is required to calculate the Shapley Value.The details of the procedure are as Procedure 5.
To obtain the Shapley Value cost allocations for our second setting, we update   () values and repeat Procedure 5.

Computational Study
We carried out a computational study to demonstrate the performance of our proposed allocation methods.We perform our experiments on randomly generated instances with varying characteristics.These instances are all generated on a region of 1,000 × 1,000 unit square and the Euclidian distances between two locations are divided by 10 to obtain the distances in kilometers.Note that, for calculating the transportation cost, we multiply the distances in kilometers by the operating cost per kilometer multiplier, which is taken as $0.30.
We generate a total of 20 instances using different parameter configurations.Half of these instances have 25 customers and the rest has 50 customers.We generate a certain portion of the customer locations in clustered regions, representing metropolitan areas.There are 3 clusters in 25customer instances and 4 clusters in 50-customer instances.Among all the customers, 40% of them are generated within the cluster regions and the remaining ones are randomly placed over the entire map.The load-carrying capacity of the trucks is assumed to be equal to 30 tons.In half of the instances, low volume instances, the shipment weight is uniformly distributed between 5% and 25% of the truck capacity.In the rest of the instances, high volume instances, the shipment weight is uniformly distributed between 10% and 30% of the truck capacity.Low volume and high volume instances are equally distributed among 25-and 50-customer instances.
A major decision in our computational analysis is the number of iterations in our iterative VRP solution procedure.Recall that  represents the iteration number of the outer loop and  represents the iteration number in the inner loop in generating the delivery patterns.Our objective is to generate as many patterns as possible; however, the computational time of the procedure is also a concern.In our computations, we perform two different sets of runs with  = 50 and  = 100, whereas  is always set to 20.In solving , we set the time limit to 3600 seconds (1 hour) and the optimality gap to 1% for the termination of the run.Finally, all the computational experiments are carried out on a 64bit Windows Server with two 2.4 Ghz Intel Xeon CPU's and 24 GB RAM.The algorithms are implemented using C++ and CPLEX Concert Technology.We test the performance of the proposed allocation methods with respect to both solution quality and computational time.Recall that the quality of an allocation method is assessed using our fairness criteria, budget balanced, and stability.As we cannot allow any budget deficiency or surplus, all the allocation methods have to satisfy the budget balanced property.Hence, we need to compute the maximum percent deviation from stability of each allocation method over all the subsets of the customers.Then, we call this maximum percentage deviation as "percent instability" of allocation method and assume that the higher this value is, the less fair the allocation is.Unfortunately, the problem comes with the fact that there are 2  subsets and this number is huge when  = 25 or  = 50.Therefore, we use the "approximate stability assessment" procedure developed by Özener et al. [18] to assess the stability of the allocation methods.The purpose of this method is to generate only a fraction of all the possible subsets and, however, at the same time make sure that these subsets are the important ones.The definition of the "importance of a subset" is having the highest possibility of being instable as well as having a high percentage instability value.In that procedure, the inherent assumption is that the customers that are closer to each other have the highest synergy potential and hence constitute an important subset.In their procedure, they pick a random point on the map and assign probabilities to the customers based on their distance to that base point.Then, they generate a subset based on these probabilities and check the percent instability of that subset.We modify this procedure for our application.In our modified procedure, we discard the generated subsets if the customers cannot be served by a single delivery route (total shipment volume exceeds the truck capacity) or their total shipment volumes are less than 75% of the truck capacity.In that sense, we assume that the highest synergy customers subsets are the ones that correspond to the efficient delivery routes.We generate 25,000 subsets for each of the 25customer instances and to 50,000 for each of the 50-customer instances.Note that the procedure does not allow duplicate subsets.
Table 2 summarizes the performance of the allocation methods when  = 50.The first four columns correspond to 25-customer instances and the next four columns correspond to 50-customer instances.Within each four-column group, each column represents an allocation method: proportional allocation method PR, individual solutions based allocation method IS, duality based method DU, and Shapley Value SV, respectively.The first three rows summarize the performance of the methods on allocating costs and the next three rows summarize the performance of the methods on allocating CO 2 emissions.The first row, "Subset, " presents the average number of instable subsets (over all the generated subsets) with each cost allocation method.The row "Ave" presents the average of the average percent instability of each cost allocation method over all the generated subsets whereas the next row "Max" presents the average of the maximum percent instability of cost allocation methods over all the generated subsets.The next three rows present the same performance measures for the CO 2 emission allocations.
The computational experiments reveal that the dualitybased method provides significantly better allocations especially compared to the proportional and individual solution based allocation methods.The average instability of dualitybased cost allocation method over all instances is equal to 3.59% on 25-customer instances and 3.48% on 50-customer instances.The same values for proportional, individual solution based and Shapley Value allocation methods are 7.26%, 7.26%, and 5.64%, respectively, for 25-customer instances and 10.14%, 10.14%, and 6.36%, respectively, for 50-customer instances.Hence, we conclude that the superior performance of the duality-based method is even more significant in larger instances.When we compare the number of instable subsets, again duality-based method outperforms the other methods.When the CO 2 emission allocations are considered, we observe the same pattern, the superiority of the duality-based method, except that proportional and individual solution based allocation methods do not yield the same values.
Table 3 presents the same results except that the classification of the 20 instances is now based on the load volumes.Hence, the first four columns present the results for the low volume instances and the next four columns present the results for the high volume instances.Even though the ranking of the allocation methods does not change with respect to the instance characteristics based on the volume of the loads, we observe that the overall performances of all the allocation methods are slightly worse in low volume instances.The reason for such a difference in performance is due to the fact that, in the low volume instances, the synergies among customers are higher, which lead to more challenging instances in terms of stability.
Tables 4 and 5 present the results when  = 100, whereas in the former table the classification is based on the number of customers; in the latter one, the classification is based on load volumes.From both tables, we observe that the performance of the duality-based method improves whereas the performances of the other allocation methods do not change significantly.Here, one might expect the exact same performance from all the allocation methods except the duality-based method as their procedures are not affected.The differences are due to selecting a different seed in random variable generation in each of these computations.Table 6 summarizes the computational times of the proposed allocation methods.The proportional method and individual solution based method clearly do not take even a second to compute and hence their computational times are not presented.The first row presents the computational times of the methods on 25-customer instances when  = 50, the second row presents the same values for 50-customer instances when  = 50, the third row for 25-customer instances when  = 100, and finally the last row for 50-customer instances when  = 100.We observe that duality-based method does not take longer than a second to compute even on the 50-customer instances when  = 100.Even though the Shapley Value method has the highest computational time of all the allocation methods, it is still within the acceptable limits.Finally, the total computational time values also include the stability assessment procedure run times, which are significant as expected.However, we should stress that the assessment computations should not be included in the running time of the allocation methods.We conclude that our proposed methods are not only effective but also computationally efficient.
Under our second setting, we provide additional incentives to decrease the emissions by switching the emission responsibility to the carrier.From Table 7, we observe that the second setting in fact decreases the total emission values.The average percentage improvement is 5.39% and the maximum percentage improvement is 14.60%.
Tables 8 and 9 present the results for the joint cost allocation methods under both classification rules, whereas in the former table  = 50 and in the latter  = 100.From both tables, we observe that the relative performances of the methods are similar with respect to each other; however, the instability values of all the methods are worse compared to their performance in allocating cost and emissions individually.This is inline with our expectations as in the joint cost setting opportunities for cost saving for the subsets might be higher, which leads to more restrict stability conditions.Finally, Table 10 summarizes the computational times of the joint cost allocation methods.Compared with the previous performance of the allocation methods, we observe that there is an increase in terms of CPU seconds.This due to the additional complexity of the TSP procedure.Fortunately, the increase is not significant except the stability assessment procedure run time, which is not included in the actual running time of the algorithms.
Based on our detailed computational study, we conclude that the Shapley Value approximation performs better than proportional and individual solution based allocation methods.The duality-based method significantly outperforms all other methods both in allocating costs and CO 2 emissions with respect to the average and the maximum percent instability as well as the number of instable subsets.Finally, the average instability value of the duality-based allocation method is always below 6% regardless of difference instance characteristics.

Concluding Remarks
In this paper, we study the cost and emission allocation problem of a logistics system.Road transportation is considered as the largest contributor to the world's GHG emissions and consequently to global warming.In order to mitigate the negative effects of the transportation activities to the environment, certain initiatives are undertaken, such as emission cap and trade in the EU.These initiatives try to give incentives to companies to reduce their emissions resulted from their direct and indirect actions.However, it is quite difficult to assess the companies emission responsibilities resulted from their indirect actions.
Having this in mind, we develop mechanisms for the fair assessment of the emission responsibilities of customers (recipients) in a logistics system.We consider a delivery network with multiple customers served by a single carrier and develop allocation mechanisms for both costs and emissions.Our assessment criterion is the fairness of the allocations.We also analyze the setting where the carrier takes over the emission responsibility and allocates a joint cost to the customers, which provides additional incentives to decrease the emissions.Based on a detailed computational study, we observe that our proposed methods, a duality-based allocation method and an approximation for the Shapley Value, perform significantly better than the allocation methods discussed in the literature.Our computational analysis also establishes that the proposed methods are computationally efficient and can be implemented for real-life sized problems.

Table 1 :
Carbon emission factors (gCO 2 /ton-km) for 40-44 ton trucks with varying payloads and levels of empty running.
Initiate a set of delivery patterns  = {} For ℎ = 1, . . .,  Initiate the cluster centers by randomly choosing  locations out of all customer locations For  = 1, . . .,  Assign customers to cluster centers by solving VSPP For  = 1, . . .,  Generate a delivery pattern  for cluster  and add this pattern to set  Calculate the centroid for all locations assigned to this cluster center Relocate the cluster center to the location closest to center of gravity of the cluster Remove the duplicated patterns in  For  = 1, . . ., || Solve a TSP to calculate the cost and the emission of the route   and   , respectively End For Solve  to select the final  delivery patterns Procedure 1: Solving the carrier's routing problem.

Table 2 :
Instability values of the proposed methods on 25-and 50-customer instances when  = 50.

Table 3 :
Instability values of the proposed methods on low and high volume instances when  = 50.

Table 4 :
Instability values of the proposed methods on 25-and 50-customer instances when  = 100.

Table 5 :
Instability values of the proposed methods on low and high volume instances when  = 100.

Table 6 :
Computational times of the proposed methods in CPU seconds.

Table 7 :
Percentage improvement values in the emission values under the second setting.

Table 8 :
Instability values of the proposed methods for joint cost allocation when  = 50.

Table 9 :
Instability values of the proposed methods for joint cost allocation when  = 100.

Table 10 :
Computational times of the proposed methods for joint cost allocation in CPU seconds.