A BRANCH-AND-PRICE ALGORITHM FOR A ROUTING AND SCHEDULING PROBLEM FROM ECONOMIC AND ENVIRONMENTAL PERSPECTIVES

. This paper addresses a routing and scheduling problem from two different perspectives: economic and environmental. From economic perspective, we aim to optimize the vehicle routing plan to reduce the operating cost, but from environmental perspective, we aim to optimize the vehicle routing and speed decisions to reduce the carbon emissions. This research can provide two different decision plans under these two different perspectives, and the comparison of the results from the two different perspectives will be very meaningful and helpful to the logistics decision-makers. We formulate the problem using two mixed-integer programming (MIP) models with different objectives. However, this problem is very challenging, with medium-sized instances already difficult for the MIP solver. In order to solve it with larger scale instances, we propose an exact branch-and-price (BAP) algorithm. The BAP algorithm relies on efficiently solving the pricing sub-problem with different objectives. We design two different tailored labeling algorithms to solve it. Extensive computational experiments demonstrate the effectiveness of the proposed BAP algorithm, comparing with the MIP formulation solved by CPLEX with a time limit of 2 h.


Introduction
With the vigorous development of the e-commerce, logistics and freight transportation are particularly important.The main role of logistics in the digital world is to use real-time data to generate vehicle routes, so as to reduce the negative consequences such as congestion, safety and environment [32].This paper studies such a logistics problem: a large number of customers who has demand are usually distributed in a town, a village or a city.Each customer has a different service time horizon (also called time window in the paper) and a different demand quantity.These customers will be serviced by several homogeneous trucks with the same capacity.Therefore, to a certain extent, the problem studied in this paper is similar to a vehicle routing problem with time window (VRPTW) considering carbon emissions [20,30].
Generally, the transportation cost will influence the optimization decisions of logistics company, and the goal is to design a valid route plan to reach the best operations cost [10,28].Except the economic factors of transportation, sustainable transportation is an active and important research topic today, which refers to any green means of transportation with little impact on the environment [22,23].Generally, the transportation sectors are responsible for a large portion of the pollution for Green House Gas (GHG) emissions worldwide [17].According to the work reports by the US Environmental Protection Agency (EPA), the transportation sector contributes 28% of national GHG emissions [8].In other words, the transportation sector of this routing and scheduling problem is to transform the driver to support the customers' required goods by cars or other means of transportation [9].Among these logistics activities, GHG, especially carbon dioxide (CO 2 ) emissions, are the most concerning because CO 2 emissions have direct influences on people's health [1].Based on this motivation, this study develops a routing and scheduling model under economic and environmental perspectives to design the route plan for the logistics company.
In this paper, the aim is to design a reasonable logistics route plan under some constraints, and the objective is to minimize operating cost.Moreover, we also test the objective of minimizing the total carbon emissions, which has a positive linear relationship with fuel consumption [24].Generally, this goal can reduce environmental pollution while optimizing operating costs for the logistics company.In order to verify the rationality of this goal, this paper compares the routing plan of some instances under different goals and compares their operating costs.
The carbon emissions along a route are calculated by a formulation based on the traveling speed.In this paper, we set several speed data for selecting by the vehicles' drivers.Therefore, this paper not only studies the route of visiting customers, but also determines the speed of the vehicles.The consideration of carbon emissions makes the problem be a route and speed joint optimization problem.There is no doubt that our scheduling problem is more complicated than the classical VRPTW.Since VRPTW is NP-hard problem [30], our scheduling problem is of course NP-hard.In order to solve the problem, an exact branch-and-price (BAP) algorithm is proposed in this paper.Experimental results highlight the efficiency of the proposed approach thanks to a comparison of the results with the mixed-integer programming (MIP) formulations solved by CPLEX.
The main contributions of this paper can be denoted as follows: (1) this paper addresses a routing and scheduling problem from economic perspective and environmental protection perspective, and the studied problem has a different objective under different perspectives; (2) an exact branch-and-price approach is proposed to solve the problem, in which two tailored labeling algorithms are used to solve the pricing sub-problem; (3) the experimental results show the effectiveness of the proposed BAP algorithm and can be used as a benchmark for future research of related problem.
The rest of this paper is organized as follows.The relevant literature is investigated in Section 2. Section 3 introduces the problem and mathematical model.Section 4 illustrates the column generation algorithm for its relaxation problem.Section 5 develops a BAP approach in order to solve the problem.The computational experiments are described in Section 6. Section 7 concludes the paper.

Literature review
With the rapid development of environmental awareness in modern logistics industry, more and more scholars studied vehicle routing problem (VRP) with the green issues [2].Bektaş and Laporte [1] first introduced the pollution-routing problem (PRP) in 2011, which is a variant of the VRP.This problem considered the relationship between carbon emissions and speed, was modeled as a mixed-integer linear programming (MILP) model.In 2012, Demir et al. [6] proposed an adaptive large neighborhood search (ALNS) for the PRP to minimize a function combining fuel, emission and driver costs.After that, Kwon et al. [18] considered carbon emission in logistics systems, and proposed a Tabu search (TS) algorithm to minimize the sum of variable operation costs and carbon emission trades.In 2014, Tajik et al. [31] presented a new model for PRP with pickup and delivery, and they introduced a robust counterpart of the MILP model to deal with uncertainty.Based on the previous works, Lin et al. [19] in 2014 systematically reviewed the green VRP (GVRP) from energy saving, emissions reduction and reverse logistics; moreover, they provided a classification of GVRP.Qian and Eglese [25] in 2016 produced routes and schedules for a fleet of delivery vehicles in order to minimize the fuel emissions in a timevarying road network, and they proposed a column generation based TS algorithm to solve this problem.In 2017, Dabia et al. [5] developed an exact branch-and-price (BAP) algorithm for a variant of the PRP, in which the speed of each arc is assumed to be the same.Teoh et al. [32] proposed a data driven multi-objective differential evolution (MODE) algorithm to solve the safe capacitated vehicle routing problems (CVRP) by minimizing the GHG emissions and hazardous risk in 2018.Fathollahi-Fard et al. [10] in 2019 studied the environmental pollution or green emissions into the home health care (HHC) routing and scheduling problems, and proposed a new modified simulated annealing (SA) algorithm to solve this problem.Cai et al. [2] in 2021 developed a hybrid particle swarm optimization (HPSO) algorithm to solve a low-carbon VRP.Luo et al. [22] in 2021 studied a routing and scheduling problem with the consideration of carbon emissions, and they proposed an ant colony optimization (ACO) heuristic for the problem.
Based on the above studies related to VRP with green issues, we can find that it is meaningful to study green logistics.However, when studying green logistics, we should not only focus on energy consumption or exhaust emissions.We also have to pay attention to operational costs.In my opinion, it is very valuable to study green logistics, especially when considering operating costs, it is still possible to conduct research on green logistics.Therefore, this paper considers two objectives, and conducts two independent experiments.One is with the objective of operating cost (in the form of transportation costs), the other is with the objective of carbon emissions.This research can provide two different decision plans under these two different perspectives, which is very interesting and meaningful for a logistics organization.In order to solve the studied problem, we develop an exact BAP approach.There is no doubt that exact algorithm will be very interesting and have theoretical meaning.
Related to the VRP with green issues, some good exact algorithms have been developed in the literature.For example, Qian and Eglese [25] used a column generation based TS algorithm to solve a green VRP in 2016.They designed a column generation method to solve the master problem, and the pricing sub-problem was solved by a TS heuristic.Dabia et al. [5] in 2017 developed an exact BAP algorithm for a variant of the PRP, which is the first time that scholars have not used the algorithm based on meta-heuristics to solve VRP with green issues.In this problem, the master problem is a set-partitioning problem solved by means of column generation, and the pricing sub-problem is a speed-and start-time elementary shortest path problem with resource constraints solved by a tailored labeling algorithm.New dominance rules are developed to discard unpromising labels.In 2019, Yu et al. [34] developed a BAP algorithm for the heterogeneous fleet green VRP with time windows.In this study, in order to solve the pricing sub-problem, a multi-vehicle approximate dynamic programming (MVADP) algorithm based on the labeling algorithm is developed.This research is similar to our study.However, only one speed is considered in this research.multiple speed selections make our study more complicated.In 2021, Wang et al. [33] proposed a BAP algorithm to solve a green location routing problem with multi-type charging infrastructure, where initial feasible columns are given by a hybrid heuristic algorithm, the pricing sub-problems are solved by the label-setting algorithm, and the global lower bound is raised by the Lagrangian lower bound.

Problem description
Throughout this paper, we work on the delivery version of the routing and scheduling problem with homogeneous vehicles from economic perspective and environmental protection perspective.The aim is to optimize the daily routes plan in order to minimize operating cost or carbon emissions.The problem can be defined as follows.Let  = (, ) be a directed graph with a set of nodes  = {0, 1, . . ., ,  + 1} and a set of arcs  = {(, )|,  ∈ ,  ̸ = }.Node 0 and node  + 1 represent the star depot and the end laboratory, respectively.Nodes  = {1, 2, . . ., } represent the customers who has demand from the logistics company.
Each customer  ∈  has a demand   , and each vehicle has the same load capacity .Each customer  ∈  is associated with a service duration   .Each customer  ∈  has a service time window [  ,   ], where   represents the earliest time and   represents the latest time for visiting the customers.Each driver is allowed to arrive before the earliest time   , but the driver must wait until the time is available for the customer and the driver can start working.The driver is prohibited to arrive after the latest time   .The start working time at customer  is denoted by   .The start depot and the end depot have the same time window, meaning the drivers must leave from the start depot and return to the end depot between the earliest time and latest time.There is a fixed planning horizon [0,  ] for the logistics company.
The distance between  and  is denoted as   .The carbon emission between  and  is calculated by a formulation denoted in Section 3.2 based on the distance and speed between  and .The speed   between  and  is an integer variable in the paper.Based on the speed   , it is very easy to calculate the travel time   /  between  and .
The problem is developed to determine a set of routes in order to minimize the objective (operating cost or carbon emissions) under the constraints of time windows, capacity, and the following assumptions: (1) each vehicle leaves from the start depot, deliver the samples to the end depot, and at last returns to the depot, and visits each customer at most once; (2) we set the cost from the start depot to end depot as the fixed cost, which is not considered in the paper; (3) the speed of the vehicle is assumed to be an average speed in each arc.

Carbon emission function
The research addressed in this paper aims to minimize the carbon emissions by optimizing vehicle speeds and routes.Carbon emissions is used to provide an estimate of the GHGs generated by vehicles.In general, the carbon emissions are measured by calculating the fuel consumption, and then multiplying by the carbon emissions conversion factor, namely the carbon emissions have a linear relationship with fuel consumption [24].The carbon emission per unit distance traveled (kilogram per kilometer, kg/km) at speed  is  (), which is developed by the United Kingdom Transport Research Laboratory [15].The emissions function has been used by many researchers, such as Jabali et al. [16], Teoh et al. [32] and so on, which can demonstrate the effectiveness of the emission function.The emissions function  () is provided as follows: where  is the speed of the vehicle in km/h, and the coefficients ,  1 ,  2 ,  3 ,  4 ,  5 and  6 will be different under the vehicles with different types and sizes.The coefficients in this paper are adopted the settings in Hickman et al. [15], and the values of ,  1 ,  2 ,  3 ,  4 ,  5 and  6 are 765, −7.04, 0, 0.006320, 8334, 0, 0, respectively.The vehicle will emit  () kg/km carbon dioxide (CO 2 ) when the vehicle is driven at the speed .It is easy to know that the emissions formula is a function of speed .Based on a survey of the speed limit of French cities, we set three speed selections in this paper, respectively 30 km/h, 45 km/h and 60 km/h.

Mathematical model
As mentioned before, this work focuses on the delivery version of the routing and scheduling problem with homogeneous vehicles.The aim is to optimize the daily routes decisions in order to minimize operating cost or carbon emissions.Based on the problem description, this paper develops a mixed-integer programming (MIP) formulation in this section.In order to model the problem clearly, we introduce two important binary variables.One is the most widely used three-index binary variable which is presented as follows: The other is a variable associated to speed.Based on the above notations, the speed used on arc (, ) is undoubtedly represented by     .However,   is a variable, and   is also a variable, so it is obvious that the function is non-linear.In order to linearize the mathematical model, we generate a new binary decision variable   which is denoted as follows: As mentioned in the previous chapter, this paper studies two different objectives and conducts two independent experiments.One objective is the operating cost, and it is in the form of transportation costs.The equation of this objective is presented as follows: where the objective is the total travel distance, which is the most common objective optimized in routing and scheduling problems.The other objective studied in this paper is carbon emissions, which have a linear relationship with fuel consumption.To a certain extent, The carbon emissions can also reflect operating cost.The formulation of this objective is denoted as follows: It is clear that the difference between two objectives is the carbon emission function.Before introducing the MIP formulation, this paper summarizes the notations used in the MIP model, which are shown in Table 1.
Then, we will introduce the constraints of the mathematical models, and the formulations are presented as follows: It should be noted that the objective of MIP1 is distance, but not about speed.Therefore, we only need to set the fastest speed, then we can meet all the constraints and get the best operating cost (or total distance).In these two MIP models, constraint (3.5a) guarantees that each customer is visited only once.Constraint (3.5b) ensures the flow balance of the vehicles, i.e., the driver visits the customer and then will leave the customer.Constraints (3.5c) and (3.5d) ensure that the vehicles start at the start depot and end at the end depot.Constraint (3.5e) is the flow equation for the demand of customers, and constraint (3.5f) is the capacity constraints.Constraint (3.5g) denotes that the vehicle  cannot arrive at  before   +   +   /  , the reason is that the vehicle  needs the service duration   and travel time from  to .Constraint (3.5h) ensures the time window of the customer .Constraints (3.5i) and (3.5j) ensure that the decision variables   and   are binary.Constraint (3.5k) is relationship between two decision variables   and   .Constraints (3.5l) and (3.5m) ensure the non-negative.

Column generation
The MIP formulation of the studied problem can be directly solved by commercial optimization solvers.However, the computational efficiency of these algorithms degrades significantly as the problem scale or the instance size increases.To overcome this difficulty, we develop a column generation algorithm to solve the relation of routing and scheduling problem in this section.Firstly, we give a set-partitioning-based formulation of the problem; then, we develop the column generation to solve the linear programming (LP) relaxation problem; Finally, we present the labeling algorithm and several effective heuristics to solve the pricing sub-problem.

Set-partitioning formulation for the studied problem
We define Ω as the set of feasible paths.A path is feasible if it satisfies capacity and time window constraints.Let   be a binary variable deciding whether path  is included in the optimal solution or not, define   as the cost (operating cost or carbon emissions) of the path , and let   be a binary variable that denotes the customer  is visited by the path  or not.We formulate the studied problem as a set partitioning model, which is presented as follows: where the objective function  2 (4.1a) minimizes the operating cost or carbon emissions of the chosen paths, constraint (4.1b) guarantees that each customer  ∈  is visited only once, and constraint (4.1c) ensures that the decision variables are binary.We define the LP relaxation of the set-partitioning model (4.1) as the master problem (MP).We use column generation [7] to solve the MP with a small subset Ω ′ ⊆ Ω of feasible paths.The MP with the subset Ω ′ is denoted as the restricted master problem (RMP), and the RMP can be formulated as follows: where the sub-problem that adds feasible routes (also called columns) to the RMP is denoted as the pricing problem.

The pricing sub-problem
The pricing sub-problem constructs a feasible route with a minimum reduced cost, using the dual values obtained from the LP solution of the RMP.If the constructed route has negative reduced cost, its corresponding column is added to the RMP.Otherwise, the LP procedure will be terminated with an optimal solution to the continuous relaxation of the MP.The pricing problem searches for the routes with a negative reduced cost, and its objective function is defined as follows: where c is the reduced cost of path , and   is the dual variable associated with the formulation (4.1b).

The pricing problem with the objective of operating cost
As mentioned in Section 3.3, the objective is operating cost, and we use the total distance to present the operating cost, which is not about speed.Therefore, we only need to set the fastest speed, then we can meet all the constraints and get the best operating cost (total distance).In other words, we only need to set one speed level (the fastest speed), which is enough to meet all the constraints and have no influences on the total distance.In order to solve this problem, we design a labeling algorithm, which will be introduced in the next section.

The pricing problem with the objective of carbon emissions
For the pricing problem with the objective of carbon emissions, it involves a speed optimization problem (SOP).Then, we will introduce the SOP in detail.For a feasible path  = ( 0 ,  1 , . . .,  ℎ ), the problem that calculates the optimal carbon emissions   is defined as a SOP.Let   := {  | = 0, 1, . . ., ℎ} be the vertex sets of path , let Arc  := {(  ,  +1 )| = 0, 1, . . ., ℎ − 1} be the arc sets of the path , and let   be a three-index binary variable that denotes the arc (, ) ∈ Arc  is traveled with a speed level  ∈ .As mentioned above, we use three speeds, namely 30, 45, and 60 km/h, which corresponding to the speed level 1, 2 and 3, namely  = {1, 2, 3}.  is a decision variable that presents the start working time at vertex  ∈   .The SOP is formulated as a MIP formulation, we name the MIP formulation as MIP3, which is presented as follows: where the objective function  5 (4.4a) minimizes the carbon emissions of the path , constraints (4.4b) and (4.4c) guarantee no violation for the time windows, and constraint (4.4d) ensures that the decision variables are binary.Since the SOP is optimized for the feasible solutions, capacity constraint need not be considered.This problem is highly related to time, it is very hard to find a best speed for each edge.Therefore, the SOP is difficult to be solved only by using a dynamic programming algorithm.In this paper, we solve the SOP by solving the MIP3 by using a commercial solver.The MIP solver is enough to solve this problem.

The labeling algorithm
The labeling algorithm has been successfully applied into many VRP-related researches [5,13,21,27].In general, the labeling algorithm is performed where the labels are extended from the start depot (i.e., node 0) to its successors.Recently, many scholars have been using bidirectional search to speed up the labeling algorithm, namely labels are extended both forward from the start depot (i.e., node 0) to its successors, and backward from the end depot (i.e., node  + 1) to its predecessors [4].It should noted that the efficiency of labeling algorithm largely depends on the dominance rules.
In this paper, in order to solve these two different pricing problems, we design two different labeling algorithms.For the pricing problem with the objective of operating cost, we design a bidirectional labeling algorithm to solve it.For the pricing problem with the objective of carbon emissions, we propose a forward labeling algorithm to solve it.The detailed introductions are below.

A tailored bidirectional labeling algorithm for the pricing problem with the objective of operating cost
As mentioned before, the cost of path is distance.Thus, we can set only one speed level (the fastest speed, namely 60 km/h in this paper) to meet all constraints, which will not influent the distance.In this section, we design a bidirectional labeling algorithm to solve the pricing problem.In the bidirectional labeling algorithm, we set a fixed time   =  /2 as the critical resource, and the forward and backward labels are not allowed to extend beyond   .Then we will introduce the bidirectional labeling algorithm in detail.
Firstly, we define the notations of forward label   .A forward label consists of   = (︀ , , , , c, , V )︀ .The related attributes are presented as follows:

𝑖:
The last node visited by   ; : Sum of demand when leaving node  in label   ; : Departure time using the fastest speed at node  in label   ; : The distance along label   ; c: Reduced cost of label   ;  : Set of nodes visited by vehicle in order along label   ; V : Set of nodes that are unreachable from node  in label   .
The forward labeling algorithm stores all of the feasible routes extended over the start depot ( = 0) to its successors.The search is restricted to elementary paths by discarding extensions to any vertex  ∈  .When the label )︀ is computed by using the following rules: The initial forward label is   = (0, 0, 0, 0, 0, {0},  ∖ {0}).In the process of forward extension, all labels should not be dominated labels.The forward dominance rule is denoted as follows.

Dominance 1 (Forward Dominance rule). A Label
) ) and at least one of the above inequalities is strictly satisfied.
The pseudo-code of the forward labeling is described in Algorithm 1.Before describing the algorithm, we need to define two sets ℒ   and ℒ   .Let ℒ   be the set of forward labels waiting for extension, and define ℒ   as the set of forward labels waiting for merging.The termination condition of the forward labeling algorithm is to determine whether ℒ   is an empty set.In order to speed up the algorithm, we set the departure time as a strict judgment resources.If the departure time is larger than  /2, we will not extend this label (Lines 6, 7).In the labeling algorithm, we ensure that the newly extended labels and all labels in the ℒ   and ℒ   are not dominated.If any label is dominated, then the label will be discarded (Lines [15][16][17][18].The output of the labeling algorithm is set ℒ   .The difference of backward labeling and forward labeling is that the backward labeling starts at the end depot, namely  + 1 to its predecessors.Similarly, we define   = (︀ , , , , c, , V )︀ as the backward label.The related attributes are presented as follows: Algorithm 1: The forward labeling algorithm.

𝑖:
The first node visited by   ; : Sum of demand when departing from node  in label   ; : Departure time using the fastest speed at node  in label   ; : The distance along label   ; c: Reduced cost of label   ;  : Set of nodes visited by vehicle in order along label   ; V : Set of nodes that are unreachable from node  in label   .
The initial backward label is   = ( + 1, 0, , 0, 0, { + 1},  ∖ { + 1}).When the label is computed by using the following rules: Similarly, in the process of backward extension, all labels should not be dominated labels.The backward dominance rule is denoted as follows.

Dominance 2 (Backward Dominance rule). A Label
) and at least one of the above inequalities is strictly satisfied.
The extension process of the backward labeling algorithm is similar to Algorithm 1, there are two sets ℒ   and ℒ   , respectively are the set of backward labels waiting for extension, and the set of backward labels waiting for merging.The line 6 in Algorithm 1 will be modified as  ≤  /2.After get the set of ℒ   , then we can merge the forward labels and backward labels.
Algorithm 2: The merging process.
Merge   and   , and new label is ; We define the fixed label as , and let ℒ  be the set of fixed paths with negative reduced cost.The process of the merging is presented in Algorithm 2. There are some constraints about merging as follows: -The last visited node of   and the first visited node of   is same, and the same node between   and   can only be this node.-The new path  will not violate the capacity and time window constraints.
-The reduced cost of new path  will be smaller than 0, in which c = c   + c   +     .

A tailored forward labeling algorithm for the pricing problem with the objective of carbon emissions
In this section, a tailored forward labeling algorithm is proposed for the pricing problem with the objective of carbon emissions.Because the pricing problem involves a SOP, which is closely related to the departure time.If we get a backward label and the departure time of the first node has changed, the whole path will be calculated again.Therefore, in this section, we only use forward labeling algorithm to solve the pricing problem.
We first define some notations related to a label.Given a label   = (︀ , ,   ,   , , c, , V )︀ , where the tuple represents the state of the path associated with the label   .The related attributes are presented as follows:

𝑖:
The Set of nodes visited by vehicle in order along label   ; V : Set of nodes that are unreachable from node  in label   .
Before presenting the extension of the labels, we review the relationship between speed and carbon emissions.As mentioned in Section 3.2, three speeds can be selected by the caregiver, namely 30 km/h, 45 km/h and 60 km/h.According to the parameters setting in Section 3.2, we can calculate that  (30) = 8.0432,  (45) = 26.546,and  (60) = 62.974.Therefore, this paper concludes that  (30) <  (45) <  (60).In this paper, we use three notations to denote this three speeds, and shown as follows:  min = 30 km/h,  mid = 45 km/h and  max = 60 km/h.Thus, it concludes that  ( min ) <  ( mid ) <  ( max ) in this paper.
The forward labeling algorithm stores all of the feasible routes extended over the start depot to its successors.The search is restricted to elementary paths by discarding extensions to any vertex  ∈  .When the label is computed by using the following rules: Before updating the other notations tuple )︀ , this paper first proposes a proposition, which is presented as follows.Proof.The departure time  ≥   +   , and   =   +   , so the departure time   is the earliest departure time.  is the departure time of best carbon emissions at vertex  in label   , which means it is impossible to reduce carbon emissions in this path up to vertex .Therefore, no matter what vertex the label extends after vertex , the carbon emissions of the path up to vertex  will not decrease nor increase, and the earliest departure time will also not be changed.Thus, in label   , the vertex  seems a new "depot" for the vertex expanded after .So the vertex  is thus as optimal vertex.)︀ should be re-optimized for the route from the last optimal vertex to vertex  using MIP3 ((4.4a)−(4.4d)).

Based on the
Proof. is the best carbon emissions of the label   , and the related best departure is   .And as mentioned before,  ( min ) <  ( mid ) <  ( max ).So the minimal carbon emissions of arc (, ) is    ( min ).Thus if   +   / min ≤   holds, it is easy to get that the minimal carbon emissions of the path up to vertex  is  ′ =  +    ( min ), and other notations  ′  and c′ will be updated as presented in Proposition 4.2.As for the other situation, if   +   / min >   , then it will not be allowed to use the speed  min on the arc (, ).Therefore, we cannot decide to adjust the speed of the path up to vertex  or arc (, ) in order to achieve the lowest overall carbon emissions up to vertex .
If   +   / min >   , the speed of the path up to vertex  should be re-optimized.In this paper, the problem is solved by a MIP formulation ((4.4a)-(4.4d))using CPLEX.We can get the best carbon emissions  ′ and the corresponding  ′  .The reduced cost is updated by the formulation c′ =  ′ − ∑︀ ∈     , where  is the current path up to vertex .
After generating the new labels, the dominance rules are used to reduce the unpromising labels that cannot lead to the optimal solution in order to save the storage space and calculation time of the computer.The dominance rule is shown as follows: )︁ is dominated by the other label ) ) and at least one of the above inequalities is strictly satisfied.
The pseudo code of the forward labeling algorithm is similar to Algorithm 1.The difference is that Lines 6, 7 in Algorithm 1 are deleted.Finally, we will extend all the labels in ℒ   with negative reduced cost to the end depot, and add the paths to the column generation process.

The framework of BAP algorithm
To optimally solve the routing and scheduling problem, we develop a branch-and-price (BAP) algorithm, which is the leading exact algorithm for the routing and scheduling problems [3,12,14,26].The BAP algorithm is composed by embedding the column generation procedure into the branch-and-bound framework.That is because that the objective value of the RMP may not be integers, but it can provide the lower bound (LB) of the studied problem at the node of the search tree.The flowchart of the BAP algorithm is shown as follows.
As shown in Figure 1, it is clear that there are some initial columns for the RMP.In this paper, we use the route Depot-Customer-Depot as the initial columns.Therefore, the number of the initial columns is the number of customers of the instance.At next, we will introduce the branching and search strategies of the BAP algorithm.

Branching and search strategies
After the column generation, we first use the integer branching strategy to save computation time.The best integer solution at the root node is obtained by solving the RMP as a 0-1 integer programming problem.The objective value of the integer solution will be the initial upper bound of the studied problem.If the upper bound is equal to the lower bound.The problem is solved and accordingly the proposed BAP algorithm terminates.
In the proposed BAP algorithm, the branch-and-bound tree is explored according to a best bound first strategy.As analysed in Feillet [11], branching on a fractional   variable poses some difficulties and is impractical, and branching on   = 0 would need a more complicated modification of the pricing problem in order to check column .In contrast, including or forbidding arcs in the solution of the pricing problem will be easily achieved.According to Reihaneh and Ghoniem [26], branching on arcs is given higher priority, this empirically results in the best improvement in the lower bound.Therefore, a branching strategy on arcs is adopted in this paper [21].
Let   be the set of all columns that contain arc (, ) ∈ , ,  ∈  .The sum of the flows on arc (, ) is equal to ∑︀ ∈   .If there exists at least an arc (, ) with fractional ∑︀ ∈   , then we branch on the value ∑︀ ∈   which is the closest to the midpoint 0.5.Two new child nodes are generated accordingly by forcing arc (, ) in one node and forbidding arc (, ) in the other node.In the former case, all columns containing arcs (,  ′ ) and ( ′ , ) with  ′ ̸ =  and  ′ ̸ =  are deleted when solving the pricing problem.In the latter case, columns using arc (, ) have to be removed when solving the pricing problem.

Computational experiments
We conduct extensive computational experiments of the proposed BAP algorithm, and compare its performance with that of a MIP formulation (provided in the Sect.3.3) solved by a state-of-the-art optimization solver CPLEX.All of the algorithms in this paper are coded in C++ programming language.The RMPs and MIPs are both solved by CPLEX 12.10.0.Computational experiments are conducted on a PC with an Inter(R) Core(TM) i7-7700 CPU @3.60 GHz and a 16 GB RAM, under a Linux operating system.Computation times are reported in seconds on this machine.Each instance is optimal or not solved within a time limit of 7200 Cpu s.

Problem instances and experimental setup
In this paper, we generate the test instances based on the classical Solomon VRPTW benchmark instances [30].The Solomon VRPTW instances are very famous and widely used by a large number of scholars such as Dabia et al. [5], Shi et al. [29], Yu et al. [34], etc.According to the geographical distribution characteristics of the instances, the Solomon VRPTW instances can be divided into three categories namely C-type instances (clustered customers), R-type instances (uniformly distributed customers) and RC-type instances (a mix of R and C types).In this paper, we set up three scale instances, with 10, 25, and 50 customers respectively, and we name them small-scale, medium-scale and large-scale instances.
In order to guarantee the correctness of the new instances, we modified the Solomon's benchmark instances with reference to the works of Dabia et al. [5] and Yu et al. [34].In the basis of the Solomon VRPTW benchmark instances, the rules of generating the test instances of the proposed problems are as follows: -we set the coordinate of the end depot as (30, 40); -the planning horizon was set as 12 hour (h), and all customers' time windows were scaled accordingly using the coefficient 12/ 0 , therefore, the time window [  ,   ] of customer  in the Solomon instances was modified as [  * (12/ 0 ),   * (12/ 0 )]; -in the tests, the distance   was not changed, but we set the unit of distance as kilometers (km); -the service time was set to 0.75 h for all customers; -the other parameters were not changed.
As for the speed settings of the test instances, we set three kinds of speed namely 30 km/h, 45 km/h and 60 km/h to test the influence of different speed to the carbon emissions.
In this paper, we use the proposed BAP algorithm and MIP formulation to do two different experiments.Experiment 1 is conducted from the economic perspective, and with the objective of operating cost (distance).Experiment 2 is conducted from the environmental protection perspective, and with the objective of carbon emissions.We also compare some results with these different objectives.

Results of Experiment 1 and Experiment 2
As mentioned in Section 6.1, this paper conducts two different experiments from the economic perspective and environmental protection perspective, respectively.This section will present the experimental results of these two experiments.
Tables 2-4 report our computational results of Experiment 1 for solving the base MIP formulation in Section 3.3 (using CPLEX vs. solving the set partitioning formulation using the proposed BAP algorithm).Tables 5-7 report our computational results of Experiment 2 for solving the base MIP formulation in Section 3.3 (using CPLEX vs. solving the set partitioning formulation using the proposed BAP algorithm).In these 6 tables, Columns 1 and 2 report the information of instances, and respectively are the name and the size of instance.For CPLEX, columns 3-5, respectively, report the best upper bound (BUB1) obtained by the CPLEX solver, be optimal or just an upper bound, the Cpu time (CpuT1), and the optimality gap (Gap1) when the calculation is terminated.For the proposed BAP algorithm, columns 6-9 focus on the performance of the proposed BAP algorithm at root node stage and report the following: (i) the lower bound (LB) solved by the column generation; (ii) upper bound (UB) corresponding to the best integer solution at the root node by solving the RMP as a 0-1 integer programming problem using CPLEX; (iii) the optimality gap (Gap2 = (UB − LB)/LB * 100%) at the root node based on the UB and LB; and (iv) the Cpu time (CpuT2) used at the root node stage.Columns 10-14 report the overall performance of the proposed BAP algorithm, which is presented as the following: (i) the best upper bound (BUB2) of the proposed BAP algorithm; (ii) the optimality gap (Gap3 = (BUB2 − LB)/LB * 100%) compared with the LB at termination or within a time limit of 7200 Cpu s; (iii) the optimality gap (Gap4 = (BUB2 − BUB1)/BUB1 * 100%) compared with the BUB1 at termination or within a time limit of 7200 Cpu s; (iv) the number of explored branch-and-bound tree nodes (Nodes); and (v) the total Cpu time (CpuT3).
An entry "-" in the table means that the algorithm is not able to obtain an associated bound within the time limit of 7200 s.As for the proposed BAP algorithm, when the time limit 7200 s is reached, the algorithm will not be terminated until it finishes processing the current branch-and-bound node.In these 6 tables, in column "BUB2", there are many "-", which means that the best upper bound not be found in branch-and-bound explored stage, and we set the UB at the root node by solving the RMP as a 0-1 integer programming problem using CPLEX as the BUB.average gap (Gap1); and (iv) the average Cpu calculation time.Columns 7-14 summarize the performance of the proposed BAP algorithm and report in the following: (i) the number of solved instance and total instance (N/T); (ii) the ratio (Ratio2) between solved instances and total instances; (iii) the average gap (Gap2) between UB and LB; (iv) the average Cpu time (CpuT2) used at the root node stage; (v) the average gap (Gap3) between BUB2 and LB; (vi) the average gap (Gap4) between BUB2 and BUB1; (vii) the average number of explored branch-and-bound tree nodes (Nodes); and (viii) the average total Cpu time (CpuT3).
As for Experiment 1 with the instance scale of 10 customers, it is clear that the MIP solver can solve all the instances optimally, and the gap between the proposed BAP algorithm and MIP model is 0, which illustrates that the proposed BAP algorithm can also solve all the instances with 10 customers optimally.Besides, the average Cpu time of MIP solver is 7.99 s, but the Cpu time of the BAP algorithm is only 0.03 s, which fully demonstrates the efficiency of the BAP algorithm.The experimental results of the small-scale instances are very illustrative, because the small-scale instances can accurately test the accuracy and effectiveness of the proposed algorithm.With the increases of the scale of increases, the solving ability of the MIP solver becomes weaker.The MIP solver can only solve 16 out of 56 instances (28.57%) with the scale of 25.However, the proposed MIP model can solve all the instances with this scale, and the gap between the proposed BAP algorithm and MIP model is −1.61%, which highlights the solving ability and efficiency of the BAP algorithm.This conclusion is more obvious in 50-scale instances.We can see that the MIP solver can only solve 5 out of 56 instances (8.93%), and the average Cpu time is 6614.85s, which denotes that the MIP solver has difficulty in handling large-scale instances.But the proposed BAP algorithm solves 41 out of 56 instances (73.12%), the gap between the proposed BAP algorithm and MIP model is −21.08%, and the Cpu time of BAP algorithm is 2661.59s, which obtains a great advantage in comparison with MIP solver.To sum up, the proposed BAP algorithm is effective and efficient for solve the studied problem from economic perspective, and has the ability to solve the instances with relatively large scale up to 50.
Similarly, three scale instances are used to test the proposed BAP algorithm in Experiment 2 .As shown in Table 8, for the small-scale instances with 10 customers, the MIP solver solves 52 instances out of 56 instances (92.86%), and the average Gap1 is 1.09%, so the results of MIP solver are very accurate and comparable.From Table 8, the proposed BAP algorithm solves all the small-scale instances with 10 customers, and the average Gap4 between the BUB2 and BUB1 is 0, which illustrates the efficiency and effectiveness of the proposed BAP algorithm.As for the medium-scale instances with 25 customers, the overall BAP algorithm itself produced global optimal solutions in Cpu times that are significantly shorter than the MIP solver, and the MIP solver can only solve 10 out of 56 instances (17.86%) while the proposed BAP can solve all medium-scale instances, which further demonstrates the efficiency and effectiveness of the proposed BAP algorithm.As for the large-scale instances, the MIP solver seems very weak and can only solve 5 out of 56 instances (8.93%), even the MIP solver cannot give a UB for many large-scale instances (32 out of 56 instances, which is shown in Tab. 8).While the proposed BAP algorithm can solve 35 out of 56 instances (62.5%), which demonstrates the solving ability of the proposed BAP algorithm for large-scale cases.And the average gap of root node (Gap2) and overall gap (Gap3) are 3.35% and 1.48%, respectively.For the large-scale instances, the gaps are relatively small and acceptable, which further illustrates the effectiveness of the proposed BAP algorithm.
Besides, via the comparison between Experiment 1 and Experiment 2 , we can see that the proposed BAP algorithm can both solve many instances, but there are large differences in the Cpu time.The solution time in Experiment 1 is much faster than that in Experiment 2 .There are two main reasons.First but also most important is that the problem in Experiment 2 is complicated than the problem in Experiment 1 .The BAP algorithm in Experiment 2 involves a SOP solved by MIP solver.Thus the BAP algorithm will call CPLEX to solve MIP3 model for many times, which will take a lot of solution time.Second is that the BAP algorithm in Experiment 1 uses bidirectional search to speed up the labeling algorithm, but in Experiment 1 only uses the forward labeling algorithm.
For a logistics company or organization, reducing operating cost is the most important, though environmental issues deserve everyone's attention.The research in this paper gives the scheduling plan from two perspectives, which can give decision makers more choices when making the scheduling plan.

Management decision analysis
In this paper, we design a BAP algorithm to solve the routing and scheduling problem from two different perspectives.A large number of computational and experimental results show the effectiveness of the proposed algorithm.In order to make readers better understand the differences in transportation costs and carbon emissions between the two different perspectives, this section compares the results from different perspectives, and the comparison results are presented in Table 9.
As can be seen from Table 9, from the perspective of environment, carbon emissions can be reduced by 45.04%, 43.95% and 42.27% when the transportation cost is only increased by 1.27%, 1.57% and 6.33%.This result shows that the mathematical model from the perspective of environment can achieve more balanced results than the mathematical model from the perspective of economic cost.Moreover, this conclusion shows that the mathematical models from different perspectives proposed in this paper can better compare the differences of transportation cost and carbon emission, and enable decision makers to make better decisions suitable for the operation of the company.

Conclusions
In this paper, we studied a routing and scheduling problem from economic and environmental perspective.From economic perspective, we aim to optimize the vehicle routing plan to reduce the operating cost, but from environmental perspective, we aim to optimize the vehicle routing and speed decisions to reduce the carbon emissions.This research can provide both a vehicle scheduling plan with the minimal operating cost and an environment-friendly scheduling plan.We formulate the problem as two MIP models with different objectives, and try to use CPLEX solver to solve the MIP model.However, this problem is very challenging, with medium-sized instances already difficult for the MIP solver.In order to solve the studied problem with larger scale instances, we propose an effective BAP algorithm to precisely solve this problem, where the master problem and the pricing sub-problem are solved by a column generation algorithm and a labeling algorithm, respectively.The BAP algorithm relies on efficiently solving the pricing sub-problem.As for the pricing problem with operating cost objective, we design a tailored bidirectional labeling algorithm to solve it.As for the pricing problem with environmental objective, we design a tailored forward labeling algorithm to solve it.Extensive computational results show that the proposed BAP algorithm outperforms a state-of-the-art MIP optimization solver, which highlights the effectiveness and efficiency of the proposed BAP algorithm.
This research can provide two different decision plans under these two different perspectives, and the comparison of the results from the two different perspectives shows that the mathematical model from the perspective of environment can achieve more balanced results than the mathematical model from the perspective of economic cost.The comparison of the results from the two perspectives will be very meaningful and helpful to the decision-makers.
At last, we discuss some future research directions.One direction is to consider traffic information into our model.Indeed, traffic has a huge impact on vehicle speed.The modeling and algorithmic questions in this direction are both challenging.Another direction is consider some uncertainty factors to increase the interest of this research.

Proposition 4 . 1 .
For a label   = (︀ , ,   ,   , , c, , V )︀ , we define   =   +   as the earliest departure time.The earliest departure time corresponds to the minimal carbon emissions of the path up to vertex .The vertex  in label   with the earliest departure time   and minimal carbon emissions is thus as optimal vertex.

Table 2 .
Computational results on 10-customers instances of Experiment 1 .

Table 3 .
Computational results on 25-customers instances of Experiment 1 .

Table 4 .
Computational results on 50-customers instances of Experiment 1 .

Table 5 .
Computational results on 10-customers instances of Experiment 2 .

Table 6 .
Computational results on 25-customers instances of Experiment 2 .

Table 9 .
Management decision analysis.