Review on Energy Resilience

Contents lists available at ScienceDirect

International Journal of Production Economics

journal homepage: www.elsevier.com/locate/ijpe

Bi-objective green ride-sharing problem: Model and exact method Yang Yua,b, Yuting Wua, Junwei Wangb,∗ a State Key Laboratory of Synthetic Automation for Process Industries, Institute of Intelligent Systems Engineering, Northeastern University, Shenyang, 110819, PR China b Department of Industrial and Manufacturing Systems Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong

A R T I C L E I N F O

Keywords: Non-linear objective Exact method Pareto-optimal ride Pareto-optimal partition Matrix diagonalization

A B S T R A C T

We investigate the bi-objective green ride-sharing problem (BGRSP) with consideration of the drivers' interests. The first objective is to minimize carbon emissions. The second objective is to maximize average ride profit so that every driver's interest can be satisfied. The average ride profit is the average profit of all used rides and it is non-linear due to the variable number of the used rides. The BGRSP is a nonlinear multi-objective problem. We develop an exact method with three steps to solve the BGRSP. The highlight of the exact method is to cut most of the non-Pareto-optimal solutions and use a decomposition method. First, we define the Pareto-optimal ride and prove that every Pareto-optimal solution of the BGRSP is composed of the Pareto-optimal rides; thus, the solution space is reduced by cutting the non-Pareto-optimal rides. Second, we define the partition (equivalent to the solution of BGRSP) based on the relationship matrix between customers and Pareto-optimal rides which is di- agonalized into several submatrices, and prove that all partitions of the relationship matrix can be obtained by the partitions of the submatrices. Therefore, the larger-scale NP-hard problem is decomposed into several small- scale NP-hard problems, each of which produces partitions of each submatrix. Third, we define the Pareto- optimal partition and prove every Pareto-optimal solution of BGRSP is composed of the Pareto-optimal partitions of each submatrix. Thus, the solution space can be significantly reduced by cutting the non-Pareto-optimal partitions, even by (1-5.5E-42)*100%. The exact method is validated by solving a benchmark instance of pdp_100-lr101 from Li & Lim benchmark with 106 customers and 200 vehicle capacity. The proposed model and method can reduce carbon emissions and make every driver satisfied simultaneously.

1. Introduction

Ride-sharing, such as SuperShuttle in USA, can drastically lower the fares for passengers by improving the utilization of vehicles. Such a transportation mode becomes more popular in many cities in China recently, as it can significantly reduce congestion and greenhouse gas emissions. For example, the “ride-sharing to the airport or railway station” service has emerged in many cities (e.g. Beijing, Shanghai, Shenyang, Chengdu etc.). The potential reason is the developing countries like China suffer the frequent occurrence of fog and haze, and the vehicle exhaust has been regarded as one of the major reasons for the deterioration of air quality in urban areas. Therefore, many Chinese metropolises encourage green ride-sharing and provide information platforms to facilitate the communication between passengers and drivers. Passengers give the ride requirements to the information plat- form and drivers provide the ride-sharing service. One key issue of green ride-sharing is to satisfy the benefit of each driver by maximizing the average ride profit, i.e., the division of total profit by the number of all used rides. However, this issue has never been examined in the

current literature. The goal of this paper is to investigate this green ride-sharing pro-

blem (BGRSP) with consideration of the drivers’ interests by proposing a bi-objective mathematical model. The first objective is to minimize carbon emissions. The second objective is to maximize average ride profit. It is noted that the second objective is non-linear because the ride profit is the division of total profit by the number of all used rides which is variable.

1.1. Related work

The ride-sharing problem is the generalization of the vehicle routing problem (VRP), similar to the Dial-a-Ride Problem (DARP), Pickup and Delivery Problem (PDP), and Vehicle Routing Problem with Time Windows (VRPTW). Thus, this section reviews related studies including the ride-sharing problem, the green VRP and the multi-objective PDP/ VRPTW.

Ride-sharing has received much attention from the transportation optimization community. Furuhata et al. (2013) surveyed the

https://doi.org/10.1016/j.ijpe.2018.12.007 Received 23 October 2018; Received in revised form 4 December 2018; Accepted 9 December 2018

∗ Corresponding author. E-mail address: [email protected] (J. Wang).

International Journal of Production Economics 208 (2019) 472–482

Available online 29 December 2018 0925-5273/ © 2018 Elsevier B.V. All rights reserved.

T

classification of current ridesharing systems. Psaraftis (1983) developed a dynamic programming approach to solve single-vehicle many-to- many DARP. Cordeau and Laporte (2003) presented a tabu search for the multi-vehicle static DARP. Cordeau (2006) proposed a branch-and- cut algorithm to solve the multi-vehicle static DARP. Xiang et al. (2008) presented a fast heuristics for the multi-vehicle dynamic DARP. Parragh et al. (2014) designed a branch-and-price algorithm to solve DARP with split requests and profits. It is noted that the maximization of total profit is the major objective of these studies. Wang et al. (2018) studied the effects of ride-sharing among those calling on taxis in Singapore for similar origin and destination pairs at nearly the same time of day.

However, neither green ride-sharing nor non-linear multi-objective ride-sharing has been considered in literature.

The green vehicle routing problem has been extensively examined. Bektas and Laporte (2011) proposed Pollution-Routing Problem (PRP) with an objective function that minimizes the amount of greenhouse emissions and fuel. Some extensions of PRP have been investigated. Franceschetti et al. (2013) studied a PRP with time-dependent travel times. Demir et al. (2014a) investigated a bi-objective PRP and utilized an adaptive large neighborhood search algorithm to solve it. Koç et al. (2014) investigated the fleet size and mix PRP with consideration of a heterogeneous vehicle fleet. Demir et al. (2014b) reviewed the research on green road freight transportation. Tiwari and Chang (2015) pro- posed a block recombination approach to solve a green VRP. Fukasawa et al. (2016) investigated a PRP that views the speed as a continuous decision variable within an interval. Suzuki (2016) proposed a dual- objective metaheuristic approach to solve a practical PRP from the users’ viewpoint. Zhang et al. (2018) modeled an electric vehicle routing problem with recharging stations for minimizing energy con- sumption and proposed an ant colony algorithm based meta-heuristics to solve the model. Darvish et al. (2018) compared the effect of op- erational decisions not only on costs but also on emissions and re- assessed some well-known logistic optimization problems. In ride- sharing, only Caulfield (2009) estimated the environment benefits of ride-sharing. However, none of these studies has considered exact al- gorithms for the non-linear multi-objective green VRP.

The multi-objective PDP or VRPTW has been investigated using exact algorithms, heuristics or meta-heuristics. In regard to exact al- gorithms for multi-objective linear PDP or VRPTW, Sarpong (2013) used column generation combined with ε-constraint to solve bi-objec- tive vehicle routing problems. Kovacs et al. (2015) obtained the Pareto- optimal solutions for a special case of multi-objective VRPTW. They proposed exact approaches based on the ε-constraint framework to solve the small instances with up to 10 customers. It is noted that for a linear PDP or VRPTW, ε-constraint method and parallel partitioning method (Lemesre et al., 2007) can be used to solve the exact Pareto- optimal solutions by combining the algorithms with the linear pro- gramming. In regard to heuristics or meta-heuristics for the multi-ob- jective PDP or VRPTW, Park (2001) proposed a hybrid genetic algo- rithm incorporating a greedy interchange local optimization algorithm for the vehicle scheduling problem with three conflicting objectives of the minimization of total vehicle travel time, total weighted tardiness and fleet size. Jozefowiez et al. (2009) proposed a bi-objective VRP model where the total length of routes and the balance of routes were minimized, and presented a meta-heuristic method based on an evo- lutionary algorithm to solve the model. Amorim et al. (2012) for- mulated a model of the multi-objective production and distribution planning of the perishable products problem and proposed a heuristic to solve it. Garcia-Najera and Gutierrez-Andrade (2013) solved multi- objective PDP using an evolutionary algorithm. Demir et al. (2014a) presented an adaptive large neighborhood search algorithm combined with a speed optimization procedure to solve the bi-objective PRP with minimization of fuel consumption and driving time. Wang et al. (2015) proposed a multi-objective local search and multi-objective memetic algorithm for solving the multi-objective PDP.

However, there is no exact algorithm for the non-linear multi-

objective PDP or VRPTW. In this paper, we develop an efficient exact method to solve the non-linear BGRSP.

1.2. Our contributions

Our main contributions are summarized as follows. The first contribution is to propose the bi-objective green ride-

sharing problem (BGRSP) that minimizes carbon emissions and max- imizes average ride profit simultaneously.

The second contribution is to develop an effective exact method to solve the non-linear BGRSP by cutting most of its non-Pareto-optimal solutions and using a decomposition method. The exact method can solve the ride-sharing benchmark instance of pdp_100-lr101 from Li & Lim benchmark with 106 customers and 200 vehicle capacity.

1.3. Overview of the paper

The paper is structured as follows. Section 2 describes the non-linear BGRSP. Section 3 defines the feasible subset, the ride, and the feasible ride and clarifies the complexities. Section 4 defines the Pareto-optimal rides and demonstrates the relationship matrix between the customers and the Pareto-optimal rides. Section 5 defines the partition based on the relationship matrix which is further diagonalized into several sub- matrices and proves that all partitions of relationship matrix can be obtained based on the partitions of submatrices. Section 6 defines the Pareto-optimal partition and proves that every Pareto-optimal solution of the BGRSP is composed of Pareto-optimal partitions. Section 7 de- scribes the procedure of the developed exact method for the non-linear BGRSP. In the section, we evaluate the developed exact method by solving the benchmark instance with 106 customers and vehicle capa- city of 200. Section 8 presents the conclusions. Appendixes give the PROOF of theorems and the detailed data of the used benchmark in- stance.

2. Non-linear BGRSP

2.1. Problem definition

The BGRSP is described as follows. A complete undirected graph G = {V, E} is given, where V is the set of vertices and E is the set of edges. We have V = {0} W {d}, where vertex 0 represents the depot and W represents the set of customers, each of which having qw units of products (or persons) to be picked up and delivered to d (drop off point) in the preferred time windows. A fleet of identical vehicles is located at the depot. Each vehicle has a capacity Q and runs at a constant speed v. The products (or persons) with each customer is no more than Q. One objective is to minimize the carbon emissions and the other objective is to maximize the average ride profit. The BGRSP is a static many-to-one ride-sharing problem, in which all information is known in advance and there are many pickup points and one drop off point.

2.2. Notations

The following notations are defined.

0: The depot; d: The drop off point; Q: The vehicle capacity; W: The set of the customers; |W|: The number of customers in W; w: Customer w, w∈W; qw: Integer units of product (or persons) in customer w; Wq: The qth subset of W, Wq ⊆ W, q = 1,2, …,2

|W|−1, Wq, q = 1,2, …,2|W|−1 = W; |Wq|: The number of elements in Wq; R: The set of all the feasible rides;

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

473

|R|: The number of feasible rides in R; i: Index of feasible rides, i=1,2, …, |R|; Ti: Feasible ride i; cei: The carbon emissions of ride i; di: The distance of ride i; fi: The consumed fuel of ride i; ci: The cost of ride i, ci = cd * di + cf*fi, where cd is the cost per unit distance and cf is cost per unit fuel; Of: The total fare paid by customers;

=a w i{1, customer is in ride 0, otherwisew i,

;

=x i{1, ride is selected 0, otherwisei

;

=

=

O c x

x f i

R i i

i R

i

1 | |

1 | | : The average profit of ride, where =O c xf i

R i i1

| | is total

profit and = xi R

i1 | | is the number of all used rides.

2.3. Model

The non-linear BGRSP is modeled based on set-partition formulation as follows.

= ce xmin ,

i

R

i i 1

| |

(1)

=

=

O c x

x max ,f i

R i i

i R

i

1 | |

1 | |

(2)

= = =

a x w W1, 1,2, ..., | |. i

R

w i i 1

| |

, (3)

Objective (1) minimizes the carbon emissions. Objective (2) max- imizes the average ride profit. Equation (3) guarantees each customer is served exactly once.

Objective (2) is non-linear and thus the model cannot be solved by the ε-constraint method combined with the algorithms for linear pro- gramming. Therefore, we develop an effective exact method to obtain the Pareto-optimal solutions.

3. Feasible subset, ride, and feasible ride and the complexities

To rapidly solve the non-linear BGRSP, we first define the feasible subsets of W, the ride and the feasible ride.

DEFINITION 3.1. (Feasible subset, Wq). Wq of W is to satisfy:

(1) total product (or persons) in Wq is not more than Q; and (2) the customers in Wq can be delivered to the drop off point together

by one vehicle.

For the set of W, there are 2|W|-1 non-empty subsets. Using DEFINITION 3.1, we cut the infeasible subsets.

LEMMA 3.1. The number of the feasible subsets of W is at mostC W| | 1 +C W| |

2 + …+C W

Q W | | min{ , | |}.

PROOF. Case 1: |Wq| = 1, at most C W| | 1 feasible subsets

Case 2: |Wq| = 2, at most C W| | 2 feasible subsets.

And so on. Case Q Wmin{ , | |}: |Wq| = Q Wmin{ , | |}, since qw is a positive integer,

an arbitrary feasible subset contains at most Q Wmin{ , | |}customers. This means there are at most C W

Q W | | min{ , | |} feasible subsets.

W Q W| | min{ , | |}q , so these cases cover all the possibilities. Thus, we can conclude that the number of feasible subsets of W is at most C W| |

1 +C W| | 2 +…+C W

Q W | | min{ , | |}. □

So LEMMA 3.1 means when <Q W| | at least + + ++ +C C C...W

Q W Q

W W

| | 1

| | 2

| | | | non-feasible subsets are cut. Thus the com-

putation time of producing all feasible rides can be greatly saved by

cutting the infeasible subsets. Subsequently, all the rides are produced based on the feasible sub-

sets.

DEFINITION 3.2. (Ride of Wq). A ride of Wq is formed by a vehicle departing from the depot (0), picking up and delivering all customers in Wq to the drop off point (d).

LEMMA 3.2. The number of the rides of Wq is |Wq|!. PROOF. |Wq|! sequences are in Wq, i.e., the permutation of custo-

mers in Wq. According to DEFINITION 3.2, adding 0 and d into an ar- bitrary sequence forms a ride of Wq. Therefore, there are |Wq|! rides for Wq. □

For example, the subset Wq with 2 customers labeled as 1 and 2, there are 2! rides, i.e., 0- > 1- > 2- > d and 0- > 2- > 1- > d.

DEFINITION 3.3. (Feasible ride of Wq). A feasible ride of Wq satisfies all the constraints of the customers in Wq. The feasible ride is shared by all the customers in the ride.

Definition 3.3 means the more constraints, the fewer feasible rides.

LEMMA 3.3. The number of the feasible rides of Wq is at most |Wq|!. PROOF. When all rides are feasible, the number is maximal.

LEMMA 3.4. The number of all the feasible rides of W is at most + + +C C C Q W1! 2! .. min{ , | |}!W W W

Q W | | 1

| | 2

| | min{ ,| |} .

PROOF. Combining LEMMA 3.1 and LEMMA 3.3. □

Table 1 demonstrates the feasible subsets, 15 rides and 8 feasible rides in a simple case with three customers labeled as 1, 2 and 3.

DEFINITION 3.4. (Feasible solution of ride-sharing). A feasible solution of ride-sharing is a subset of feasible rides (ST), if and only if ST satisfies:

(1) =CT CT i i T T ST, , ,i i i i , where CTi is the customer set of feasible ride i (Ti); and

(2) == CT W T ST,iST i i1| | .

Definition 3.4 means that each feasible solution of ride-sharing is composed of one or more feasible rides; thus, non-feasible rides can be cut. Based on the feasible rides in Table 1, Table 2 shows 6 feasible solutions of the simple case.

Using 6 feasible solutions we can obtain the exact Pareto-optimal solutions for non-linear BGRSP. However, the number of produced feasible solutions using all the feasible rides may be enormous. We propose the definition of Pareto-optimal ride and prove that every Pareto-optimal solution of non-linear BGRSP is composed of Pareto- optimal rides. Thus, non-Pareto-optimal rides are cut so that the solu- tion space can be effectively reduced.

Table 1 Feasible subsets and rides.

Wq Rides of Wq Feasible

{1} 0- > 1- > d ✓ {2} 0- > 2- > d ✓ {3} 0- > 3- > d ✓ {1,2} 0- > 1- > 2- > d ✓

0- > 2- > 1- > d – {1,3} 0- > 1- > 3- > d ✓

0- > 3- > 1- > d – {2,3} 0- > 2- > 3- > d ✓

0- > 3- > 2- > d ✓ {1,2,3} 0- > 1- > 2- > 3- > d ✓

0- > 1- > 3- > 2- > d – 0- > 2- > 1- > 3- > d – 0- > 2- > 3- > 1- > d – 0- > 3- > 1- > 2- > d – 0- > 3- > 2- > 1- > d –

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

474

4. Pareto-optimal rides and relationship matrix between customers and Pareto-optimal rides

4.1. Pareto-optimal rides

We define the dominance relation.

DEFINITION 4.1. (Dominance relation of ride, t). Let Rq be the set of all the feasible rides of Wq and t t R, q. t t t tdominates ( )t if and only if t t t tand1 1 2 2, where at least one inequality is strict; t tand1 2 are carbon emissions and cost of t respectively; and t tand1 2 are carbon emissions and cost of t respectively.

DEFINITION 4.2. (Pareto-optimal ride of Wq). A ride x Rq is Pareto- optimal ride of Wq if and only if

¬ x Rq such that x xt .

LEMMA 4.1. The number of Pareto-optimal rides of Wq is no more than |Wq|!.

PROOF. When all feasible rides of Wq are Pareto-optimal, the number is maximal.

In practice, the number of Pareto-optimal rides of Wq is so far less than |Wq|!, especially when |Wq|! is large.

DEFINITION 4.3. (Dominance relation of solution, s). Let S be the set of all the solutions and s s S, . s s s sdominates ( )s if and only if s s s sand1 1 2 2 , where at least one inequality is strict; s sand1 2 are carbon emissions and average profit of rides of srespectively; and s sand1 2are carbon emissions and average profit of rides of s respectively.

DEFINITION 4.4. (Pareto-optimal solution of BGRSP). A solutions Sis Pareto-optimal solution of BGRSP if and only if ¬ s Ssuch that s ss .

THEOREM 4.1. In every Pareto-optimal solution, all rides are Pareto-optimal rides.

PROOF. See Appendix A.

According to THEOREM 4.1, the Pareto-optimal solutions can be exactly obtained using Pareto-optimal rides. Thus, the solution space can be greatly reduced by cutting the non-Pareto-optimal rides. The Pareto-optimal rides of a simple case are shown in Table 3.

POTi: Pareto-optimal ride i; cei: carbon emissions of POTi; ci: cost of POTi; Customers: customers in POTi.

In Table 3, according to the customers in POTi, i.e., the column of Customers, the relationship matrix between the customers and the Pareto-optimal rides can be obtained, as shown in Table 4.

4.2. Relationship matrix between customers and Pareto-optimal rides

The notations used in relationship matrix are defined as follows.

A: Relationship matrix, = ×A

a a a a a a

a a a (

...

... ... ... ...

... )

N

N

W W W N

W N

1,1 1,2 1,

2,1 2,2 2,

| |,1 | |,2 | |,

| | ;

aw,n: Relationship value, aw,n = 1 if customer w is contained in Pareto-optimal ride n, 0 otherwise; |W|: Number of rows, i.e., the number of customers; N: Number of columns, i.e., the number of Pareto-optimal rides;

Rw: Row w, Rw = [aw,1 aw,2 ... aw,N]. Rw expresses which rides con- tain customer w; Cn: Column n, Cn = [a1,n a2,n … a|W|,n]

T. Cn expresses which custo- mers are contained in column n.

DEFINITION 4.5. (Customer set of a column). Customer set of Cn is expressed as CSn = {w| aw,n = 1}.

CSn means that if aw,n = 1 then Cn contains customer w. Based on the relationship matrix, a partition of customer set can be

expressed as a set (SC) of columns.

5. Partition expression based on matrix and relationship matrix diagonalization

5.1. Partition expression based on matrix

DEFINITION 5.1. (Partition of W, expressed as a set (SC) of columns). A set (SC) of columns is a partition of W if and only if SC satisfies:

(1) =CS CS n n C C SC, , ,n n n n ; and (2) == CS W C SC,nSC n n1| | .

Comparing DEFINITION 5.1 and DEFINITION 3.4, a partition based on relationship matrix is equivalent to a solution of BGRSP. Thus, using the relationship matrix, all solutions (i.e., partitions) can be produced.

THEOREM 5.1. The time complexity of producing all partitions based on relationship matrix is a non-polynomial function regarding the number of customers.

PROOF. See Appendix B.

We develop a recursive algorithm to produce all partitions based on the relationship matrix. However, according to THEOREM 5.1, the al- gorithm cannot solve larger-scale instances. It takes 4.25 h to obtain all partitions of the submatrix with 68 customers. Thus, we propose an innovation method based on generalized matrix block diagonalization to decompose the larger-scale problem into several small-scale pro- blems.

5.2. Diagonalization of relationship matrix

A is a generalized block diagonal structure with a total of D diagonal submatrices, each of which is denoted as Ad and has the individual size |W|d ×Nd, d = 1,2, …,D, whereas the other nondiagonal matrix entries are all zero (Wang et al., 2012). A is expressed as Fig. 1.

= ×A

A A

A

( . . . )

D

W N

1

2 | |

Obviously, == W W| | | |d D

d1 , where |W|d is the number of rows in Ad. In addition, == N Nd

D d1 , where Nd is the number of columns in Ad.

The following notations are defined after matrix diagonalization:

Table 2 Feasible solutions produced by all the feasible rides.

Feasible solution

Expression based on set-partition

Expression based on feasible rides

1 {{1},{2},{3}} {0- > 1- > d, 0- > 2- > d, 0- > 3- > d} 2 {{1},{2,3}} {0- > 1- > d, 0- > 2- > 3- > d} 3 {0- > 1- > d, 0- > 3- > 2- > d} 4 {{2},{1,3}} {0- > 2- > d, 0- > 1- > 3- > d} 5 {{3},{1,2}} {0- > 3- > d, 0- > 1- > 2- > d} 6 {{1,2,3}} {0- > 1- > 2- > 3- > d}

Table 3 Pareto-optimal rides of a simple case.

POTi cei ci Customers POTi cei ci Customers

1 2.8 42 1 11 3.2 51 5 2 2.9 45 1,2 12 4.1 60 5,6 3 4.1 58 1,2,3 13 3 45 6 4 4 59 2,3,1 14 3.2 47 7 5 3.2 48 1,3 15 4 58 7,8 6 3.1 49 3,1 16 3.6 53 7,9 7 2.4 37 2 17 2.3 35 8 8 2.7 41 2,3 18 2.8 53 8,9 9 2.1 31 3 19 2 30 9 10 2 31 4

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

475

= = + = + = + = + = + = +

= + = + = + = + = + = +

= + = + = + = + = + = + ×

A d A

a a a a a a

a a a : Submatrix ,

...

... ... ... ... ...

... ;d d

i d W i i

d Ni i d W i i

d Ni i d W i i

d Ni Nd

i d W i i

d Ni i d W i i

d Ni i d W i i

d Ni Nd

i d W i W d i

d Ni i d W i W d i

d Ni i d W i W d i

d Ni Nd W d Nd

1 1| | 1, 1

1 1 1 1| | 1, 1

1 2 1 1| | 1, 1

1

1 1| | 2, 1

1 1 1 1| | 2, 1

1 2 1 1| | 2, 1

1

1 1 | | | | , 1

1 1 1 1 | | | | , 1

1 2 1 1 | | | | , 1

1 |

= W| |i d

i1 1 : Total number of rows contained in the first d-1 sub-

matices. += W| | 1i d

i1 1 means the order of the first row in Ad;

= Ni d

i1 1 : Total number of columns contained in the first d-1 sub-

matices. += N 1i d

i1 1 means the order of the first column in Ad;

Rwd: The part of Rw in Ad, Rwd = + + += = =a a a[ ... ]w N w N w N N, 1 , 2 ,id i id i id i d11 11 11 , where

= += a 1n

N w N n1 ,

d i d

i1 1 . The entries of Rw not in Ad are all 0;

BRwd: The part of row w before Ad, BRwd = =

[ 0 0 ... 0 ] Ni

d i1

1 ;

ARwd: The part of row w after Ad, ARwd = = +

[ 0 0 ... 0 ] Ni d

D i1

;

Rw: Row w, Rw = [PBdw Rwd PAdw] =

= = + = + = +

= +

a a a[ 0 0 ... 0 ... 0 0 ... 0 ]

i d Ni

w i d Ni w i

d Ni w i d Ni Nd

i d D Ni1

1 , 1

1 1 , 1 1 2 , 1

1

1

;

Cnd: The part of Cn in Ad, Cnd = + + += = =a a a[ ... ]W n W n W W n

T | | 1, | | 2, | | | | ,i

d i i

d i i

d i d1

1 1 1

1 1 , where

= += a 1t

W W t n1

| | | | ,

d i d

i1 1 . The entries of Cn not in Ad are all 0;

BCnd: The part of column n before Ad, BCnd =

=

[ 0 0 ... 0 ] W

T

| |i d

i1 1

;

ACnd: The part of column n after Ad, ACnd = = +

[ 0 0 ... 0 ] W

T

| |i d D

i1

;

Cn: Column n, Cn = [BCnd Cnd ACnd] =

= = + = + = +

= +

a a a[ 0 0 ... 0 ... 0 0 ... 0 ]

i d W i

i d W i n i

d W i n i d W i W d n

i d D W i

T

1 1| |

1 1| | 1, 1

1| | 2, 1 1| | | | ,

1| |

;

W_Ad: Customer set of Ad. W_Ad = w w A{ | }d . DEFINITION 5.2. (w Ad). Customer w Ad if and only if

= += a 1n

N w N n1 ,

d i d

i1 1 .

DEFINITION 5.3. (n Ad). Column n Ad if and only if

= += a 1t

W W t n1

| | | | ,

d i d

i1 1 .

After matrix diagonalization, the following properties can be con- cluded.

PROPERTY 5.1. w A d d w A, if thend d , i.e., customer w is in one and only one submatrix.

PROPERTY 5.2. n A d d n A, if thend d , i.e., column n is in one and only one submatrix.

PROPERTY 5.3. =W A W A d d_ _ ,d d , where W_Ad is the customer set of Ad.

PROPERTY 5.4. == W A W_dD d1 . Theorem 5.2and THEOREM 5.3 are proposed to prove that all

partitions of relationship matrix can be produced using the partitions of submatrices.

THEOREM 5.2. Every partition (p_A) of A can be produced by the partitions of D submatrices, i.e., =

= p A p A_ _

d

D d

1 , where p A_ dis an

existing partition of submatrix d. PROOF. See Appendix C. THEOREM 5.3. == p A p A_ _dD d1 , where p A_ dis an arbitrary

partition of Ad and p A_ is an existing partition of A. PROOF. See Appendix D.

Combining THEOREM 5.2 with THEOREM 5.3, it is concluded that all the partitions of matrix A can be produced using the partitions of submatrices. Thus, we can decompose the larger-scale problem into several small-scale problems. We develop a recursive algorithm to produce all the partitions of matrix A using the partitions of sub- matrices.

The result of block diagonalization of relationship matrix of Table 4 is shown in Table 5, where 4 submatrices are obtained.

According to THEOREM 5.3, we know the number of partitions of A equals = N P_d

D d1 , where N_Pd is the number of partitions of Ad. If

= N P_d D

d1 is large, then it is impossible to produce all partitions of A. Therefore, we propose the definition of Pareto-optimal partitions and prove that every Pareto-optimal solution of BGRSP is composed of Pareto-optimal partitions. Thus, non-Pareto-optimal partitions are cut so that the solution space can be greatly reduced.

6. Pareto-optimal partition of submatrix

DEFINITION 6.1. (Dominance relation of partition, p). Let P be the set of all the partitions of submatrix and p p P, . p p p pdominates ( )p if and only if p p p p p p, , and1 1 2 2 3 3 where at least one inequality is strict; p p p, ,1 2 3are carbon emissions, cost, and rides number of p respectively; and p p p, ,1 2 3are carbon emissions, cost, and rides number of p respectively.

Table 4 Relationship matrix between customers and Pareto-optimal rides of the simple case.

Customer Rides

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 3 0 0 1 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1

1

2

| |

. .

.

D W N

A A

A

A

Fig. 1. Generalized block diagonalization of relationship matrix.

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

476

DEFINITION 6.2. (Pareto-optimal partition of a submatrix). A partition p P is Pareto-optimal partition of a submatrix if and only if ¬ p Psuch that p pp .

THEOREM 6.1. In every Pareto-optimal solution of BGRSP (i.e., Pareto-optimal partition of matrix A), any partition of submatrix is Pareto-optimal.

PROOF. See Appendix E. According to THEOREM 6.1 and THEOREM 5.3, the number of the

solutions based on Pareto-optimal partitions is = N POP_d D

d1 , where N_POPd is the number of Pareto-optimal partitions of Ad.

Table 6 shows partitions and Pareto-optimal partitions of each submatrix for a simple case. Using partitions, 56 = 7*1*2*4 solutions are produced. However, using Pareto-optimal partitions, 2 = 2*1*1*1 solutions are produced.

After obtaining all the solutions based on Pareto-optimal partitions, we can solve the Pareto-optimal solutions based on the ε-constraint method.

7. The exact method with three steps for non-linear BGRSP and computational experiments

7.1. The exact method with three steps for non-linear BGRSP

We develop an exact method with three steps to solve the non-linear BGRSP. The first step produces all the Pareto-optimal rides. The second step produces the relationship matrix between customers and Pareto-

optimal rides, performs the matrix diagonalization, and produces the Pareto-optimal partitions of each submatrix. The third step obtains the Pareto-optimal solutions of BGRS. The brief description of the exact method is given in Table 7.

In substep 3.2, the improved ε-constraint method is utilized to ob- tain the Pareto-optimal solutions. The ε-constraint method involves a constraint on one objective and optimizes the second objective (Demir et al., 2014a). The time complexity of ε-constraint method for bi-ob- jective is O(MN), where M is the number of Pareto-optimal solutions and N is the number of all solutions. To reduce the time complexity, we propose an improved exact algorithm by performing ε-constraint in the solution space of cutting some non-Pareto-optimal solutions. To cut non-Pareto-optimal solutions in advance, we define the Pareto-optimal solution with the minimal carbon emissions (POS1) and the Pareto- optimal solution with the maximal average profit (POS2).

DEFINITION 7.1. (POS1). POS1 is the Pareto-optimal solution with the minimal carbon emissions.

DEFINITION 7.2. (POS2). POS2 is the Pareto-optimal solution with the maximal average profit.

Both POS1 and POS2 are used to cut non-Pareto-optimal solutions before performing ε-constraint method, as shown in Fig. 2.

The time complexity to obtain POS1 and POS2 is O(N). The time complexity to cut the solutions dominated by POS1 and POS2 is O(N). Assume the number of solutions dominated by POS1 or POS2 is C, the time complexity of obtaining left M-2 Pareto-optimal solutions is O((M- 2)*(N-C)). To sum, the total time complexity of the improved ε-con- straint method based on cutting is O(M(N-C)-2(N-C)+2N) = O(M(N- C)). Therefore, when M or C is large, the improved ε-constraint method based on cutting can greatly save the running time.

7.2. Method to calculate the carbon emissions of a ride

The amount of carbon emissions of a vehicle on the arc (i,j) is ap- proximately expressed as:

=e rfP ,ij ij (4)

where f is the fuel consumption index parameter, r is the carbon emissions index parameter, and Pij is the energy consumed on the arc (i,j) approximately expressed as (see Bektaş and Laporte, 2011):

+ +P a w f d v d( ) ,ij ij ij ij ij ij ij ij 2

(5)

where wij is the empty vehicle weight, fij is the load carried by the vehicle on this arc, and aij and βij are expressed as the two following equations respectively.

= + +a a g gCrsin cos ,ij ij i (6)

where a is the acceleration (m/s2), g is the gravitational constant (9.8 m/s2), θij is the arc angle, Cr is the coefficients of rolling resistance. aij is an arc specific constant, and aij(wij+ fij)dij represents the load-

Table 5 Block diagonalization of relationship matrix of a simple case.

Customers Rides

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 3 0 0 1 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1

Table 6 Partitions and Pareto-optimal partition of submatrix.

POT cei ci Customers d Partitions_Ad N_Pd POPd N_POPd

1 2.8 42 1 1 {{1},{2},{3}} 7 {{1,2,3}} 2 2 2.9 45 1,2 {{1},{2,3}} {{2,3,1}} 3 4.1 58 1,2,3 {{1,2},{3}} 4 4 59 2,3,1 {{1,2,3}} 5 3.2 48 1,3 {{2,3,1}} 6 3.1 49 3,1 {{1,3},{2}} 7 2.4 37 2 {{3,1},{2}} 8 2.7 41 2,3 9 2.1 31 3 10 2 31 4 2 {{4}} 1 {{4}} 1 11 3.2 51 5 3 {{5},{6}} 2 {{5,6}} 1 12 4.1 60 5,6 {{5,6}} 13 3 45 6 14 3.2 47 7 4 {{7},{8},{9}} 4 {{7,9},{8}} 1 15 4 58 7,8 {{7},{8,9}} 16 3.6 53 7,9 {{7,8},{9}} 17 2.3 35 8 {{7,9},{8}} 18 2.8 53 8,9 19 2 30 9

d: submatrix d; Partitions_Ad: partitions in Ad; N_Pd: number of partitions in Ad; POPd: Pareto-optimal partitions in Ad; N_POPd: number of Pareto-optimal par- titions in Ad.

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

477

induced energy requirements.

= C A0.5 ,ij d (7)

where Cd is the drag coefficient, A is the frontal surface area of the vehicle (m2), ρ is the air density (kg/m3). βij is a vehicle specific con- stant, and βij = 0.5CdAρ represents the speed-induced energy require- ments.

Given the type of vehicle and environment, the aij, wij, βij, r and f are all constants.

7.3. Software and hardware specifications

The algorithms to implement the three steps in the exact method for BGRS were coded in C#. All tests were executed on the computer with Intel Core(TM) i7-4790 processor at 3.60 GHz under the Microsoft Windows 7 operation system using 8.00 GB RAM.

7.4. Experiment parameter

We test the developed exact method using the real instances of ride- sharing to airport in China. There are at least three companies pro- viding PDCA in Shenyang, such as Zhongshan, Shuntian, and Jiantong companies (Yu et al., 2016, 2017). Ride-sharing to airport has the fol- lowing characteristics: 1) the vehicle capacity is small; and 2) the ar- rival airport time window is tight. The developed exact method can obtain the Pareto-optimal solutions for the instance of ride-sharing to airport with 330 customer points and 528 customers. The instance can cover all instances in the current ride-sharing to airport of Zhongshan company in Shenyang (Wang et al., 2018b). Also, we use the general benchmark instance to evaluate the performance of the exact method for the BGRSP.

The basic parameters are listed in Table 8.

7.5. Computational experiments for non-linear BGRSP using benchmark instance

The instance is adopted from Li & Lim PDPTW benchmark pdp_100- lr101 with 106 customers and 1496 product. The detailed data is given in Appendix F. The vehicle capacity (Q) is 200, weight of empty vehicle is 2000 kg, and frontal surface area (A) is 2 m2. Other parameters are the same as those in Table 8. Depot (0) coordinate is (40, 50), the same as that of pdp_100-lr101.

In the first step, 4339 Pareto-optimal rides are produced. The run- ning time is 62 s.

In the second step, 43,200 solutions based on Pareto-optimal par- titions are produced. The running time is 77 s. The computation results of the second step are given in Table 9.

N_POTd is the number of Pareto-optimal rides in Ad; N_Pd is the number of partitions in Ad; N_POPd is the number of Pareto-optimal partitions in Ad; =Ratio of cut 100%

N P N POP N P

_ _ _

d d d

. From Table 9, we can see that approximately 99% non-Pareto-op-

timal partitions are cut. If using all partitions of submatrix, then the number of solutions is = += N P_ 7.82 E 45d d1

15 , which states it is im- possible to produce all feasible solutions. However, if using Pareto- optimal partitions, the number of solutions is == N POP_ 43,200d d1

15 . Thus the solution space is reduced by (1-

+ 43,200

7.82 E 45 )*100%= (1–5.525E-

42) *100% by cutting non-Pareto-optimal partitions. As a result, the Pareto-optimal solutions can be obtained rapidly.

In the third step, we use the improved ε-constraint method to obtain the exact Pareto-optimal solutions. There are 185 Pareto-optimal so- lutions. The running time is 0.08 s. The 185 Pareto-optimal solutions mean that the solutions with the minimal carbon emissions and the maximal average ride profit are not identical. To operate bi-objective green ride-sharing, one of Pareto-optimal solutions needs to be selected

Table 7 Brief description of the exact method for BGRS.

Step Substep Description

1 1.1 Produce all feasible subsets of W. 1.2 Produce all sequences of each subset by permutation of the customers in the subset. 1.3 Produce all rides of each subset by adding 0 and d into each sequence of the subset. 1.4 Obtain the feasible rides by judging if a ride satisfies the constraints of the customers in the ride. 1.5 Obtain Pareto-optimal rides of each subset.

2 2.1 Obtain the relationship matrix (A) of customers and Pareto-optimal rides 2.2 Perform matrix (A) diagonalization. 2.3 Obtain all partitions of each submatrix. 2.4 Obtain Pareto-optimal partitions of each submatrix.

3 3.1 Produce the solutions of BGRS using Pareto-optimal partitions of each submatrix. 3.2 Solve Pareto-optimal solutions of BGRS by the improved ε-constraint method based on cutting.

Carbon Emission

Average Profit

POS1

POS2

Solutions dominated by POS1 or POS2 Solutions left Cutting operation

Fig. 2. POS1, POS2, and solutions dominated by POS1 or POS2.

Table 8 The parameters of the experiments.

Parameter name Value

(v) speed of the vehicles 40 km/h weight of a customer or product 100 kg (β) maximum ride time proportional to the direct ride time 1.5 (Cd) drag coefficient 0.45 (Cr) coefficients of rolling resistance 0.015 (ρ) air density 1.225 kg/m3

(f) fuel consumption index parameter 1/44000 J/g (r) carbon emissions index parameter 3.125 g/g (0) depot coordinate (20, 30) (d) drop off point coordinate (80, 50) (Q) capacity 4 weight of empty vehicle 1400 kg (A) frontal surface area of the vehicle 1 m2

(cd) cost per unit distance 1 RMB (cf) cost per unit fuel 5 RMB/kg (Of) total fare paid by customers 5000 RMB

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

478

according to the preference of manager.

8. Conclusions

In this paper, we propose a bi-objective mathematical model for the BGRSP; in particular, two objectives of the model are minimization of the carbon emissions and maximization of the average ride profit. Based on the analysis of its characteristics, we develop an exact method to

solve the model by cutting most of the non-Pareto-optimal solutions and using a decomposition method. The exact method has been vali- dated by solving the ride-sharing benchmark instance of pdp_100-lr101 from Li & Lim benchmark with 106 customers and vehicle capacity of 200. Our exact method can be used to solve the static many-to-one non- linear ride-sharing problem/dial-a-ride problem/PDP. Future work will consider the static many-to-many green ride-sharing problem with many delivery points. The method in this paper to divide the subsets using the capacity Q does not work for the many-to-many problem; therefore, the difficulty will be how to produce subsets and obtain the Pareto-optimal rides of each subsets. Furthermore, we will employ the proposed algorithm to solve other important problems with similar multiple objectives (Han et al., 2015; Wang et al., 2016, 2018a).

Acknowledgments

This research is supported by the National Natural Science Foundation of China (71571156, 71831006, 71571037, 71601089 and 71620107003), by the Fundamental Research Funds for the Central Universities (N170405005), by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. [T32-101/15-R]), by the University of Hong Kong through the Seed Fund for Basic Research (Grant No.201611159213), and by the open project funded by State Key Laboratory of Synthetical Automation for Process Industries (PAL-N201802).

Appendix A. Proof of THEOREM 4.1

THEOREM 4.1. In every Pareto-optimal solution, all rides are Pareto-optimal rides.

PROOF. Let a set of feasible rides (ST) be a Pareto-optimal solution of BGRSP covering customer set of W. Assume at least one ride is non-Pareto- optimal in ST.

Case 1: |ST| = 1, ST contains only one ride, i.e., ST = {T1}. According to the assumption, T1 is non-Pareto-optimal ride. Let POT1 be one of Pareto-optimal rides for T1. Let NST be the set containing POT1, i.e., NST = {POT1}. Customer set of POT1 is W, therefore, NST is a solution of BGRSP. Since |ST| = 1, the carbon emissions of ST (ST1) is T11 and average profit of rides of ST (ST2) is Of-T12. Also, the carbon emissions of NST (NST1) is POT11 and average profit of rides of NST (NST2) is Of -POT12. Therefore, NST1- ST1 = POT11-T11and NST2- ST2 = POT T( )1 12 2 . Since POT1 tT1, POT T POT T0 and 01 1 1 11 1 2 2 where at least one inequality is strict. Thus, NST1- ST1 0 and NST2- ST2 0 where at least one inequality is strict. Therefore, NST sST. But this contradicts the assumption that ST is a Pareto-optimal solution.

Case 2: |ST| 2, ST contains at least two rides. Assume Ti is non-Pareto-optimal ride. The carbon emissions of ST (ST1) is + += = +T T Tt i

t i t i ST

t1 1

1 | |

1 1 1

and average profit of rides of ST (ST2) is + += = +O T T T ST

( ) | |

f t i

t i t i ST

t1 1

2 2 1 | |

2 , where Ti1and Ti2are carbon emissions and cost of Ti respectively. Let POTi be one of Pareto-optimal rides for Ti. Let NST be the new set of rides, i.e., NST = {T1, …, Ti-1, POTi, Ti+1, …, T|ST|}. Customer set of POTi is same as that of Ti, therefore, NST satisfies the two conditions of feasible solution of BGRSP. The carbon emissions of NST (NST1) is + += = +T POT Tt

i t i t i

ST t1

1 1

| | 1 1 1

and the

average profit of rides of NST (NST2) is + += = +O T POT T

ST ( )

| | f t

i t i t i

ST t1

1 2 2 1

| | 2 . Therefore, NST1- ST1 = POTi1-Ti1and NST2- ST2 =

POT T ST

( ) | |

i i2 2 . Since POTi tTi, POT T POT T0 and 0i i i i1 1 2 2 where at least one inequality is strict. Thus, NST1- ST1 0 and NST2- ST2 0 where at least one inequality is strict. Therefore, NST sST. But this contradicts the assumption that ST is a Pareto-optimal solution.

These cases cover all the possibilities. Since the assumption that a Pareto-optimal solution of BGRSP is composed of non-Pareto-optimal rides led to a contradiction, in every Pareto-optimal solution, all rides are Pareto-optimal. □

Appendix B. Proof of THEOREM 5.1

THEOREM 5.1. The time complexity of producing all partitions based on relationship matrix is a non-polynomial function regarding the number of customers.

PROOF. Set-partition is a well-known NP-hard. The number of partitions of set W is expressed as = =F W P W S(| |) (| |, )S W

1 | | , where P(|W|,S) is the

count of partitioning |W| elements into S subsets. P(|W|,S) is the Stirling numbers of the second kind (Knopfmacher and Mays, 2006). F(1) to F(10) are 1, 2, 5, 15, 52, 203, 877, 4140, 21147 and 115975. When |W| > 4, < < F W2 (| |)W| | . Since there are F(|W|) partitions, the time complexity of producing all partitions is at least > >O F W O( (| |)) (2 )W| | . Thus, producing all partitions is non-polynomial function regarding the number of customers. □

Appendix C. Proof of THEOREM 5.2

THEOREM 5.2. Every partition (p_A) of A can be produced by the partitions of D submatrices, i.e., = =

p A p A_ _ d

D d

1 , where p A_ dis an existing

Table 9 Computation results of the second step for the instance.

Ad Customers in Ad N_POTd N_Pd N_POPd Ratio of cut %

A1 1,2 3 2 1 50 A2 3,4,5,6,7,8,9,10,11,12 1600 247,746 3 99.999 A3 13,14,15,16,17,18,19,20,21,22 569 86,601 2 99.998 A4 23,24,25,26,27,28,29,30,31,32 403 53,011 2 99.996 A5 33,34,35,36,37,38,39,40 318 5640 3 99.947 A6 41,42,43,44,45,46,47,48,49,50,51 632 223,638 4 99.998 A7 52,53,54,55,56,57 65 210 1 99.524 A8 58,59,60,61,62,63 75 232 3 98.707 A9 64,65,66,67,68,69,70 51 211 1 99.526 A10 71,72,73,74,75,76,77,78 139 1908 5 99.738 A11 79,80,81,82,83,84 73 263 2 99.240 A12 85,86,87,88,89,90,91 189 1182 5 99.577 A13 92,93,94,95,96,97,98 151 1152 1 99.913 A14 99,100, 101,102,103,104,105 70 388 2 99.485 A15 106 1 1 1 0

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

479

partition of submatrix d.

PROOF. Let p_A be an arbitrary partition of A and SC_pA be the set of columns of p_A. For every Ad, the set (SC_Ad) of C n A,n d can be obtained from SC_pA. Thus =

= SC pA SC A_ _

d

D d

1 . Since p_A is a partition of A, SC_Ad satisfies

=CS CS n n C C SC A, , , _n n n n d and =CS w w A C SC A{ | }, _n d n d. For arbitraryn Ad, = + + +

= = = =

= +

C a a a[ 0 0 ... 0 ... 0 0 ... 0 ]n W

W n W n W W n W

T

| | | | 1, | | 2, | | | | ,

| |i d

i i d

i i d

i i d

i d

i d D

i1 1

1 1

1 1

1 1

1

. Let |SC_Ad| be the number of columns in

SC_Ad. Then SC_Ad can be expressed as:

1 1 1 | | 1, | | 1, | | 1,

1 1 1

1 1 1 | | 2 , | | 2 , | | 2 ,

1 1 1

1 1 | | | | , | | | | , | |

1 1

1

1

0 0 ... 0 ... ... ... | |

0 0 ... 0 ...

... _

... ... ... ...

d d d W j W k W li i i

i i i

d d d W j W k W li i i

i i i

d d W W j W W k Wi d i d i

i i

d

i i

d

W

a a a

a a a SC A

a a a 1 | | ,

1

1

| | | _ |

0 0 ... 0 ... ... ... | |

0 0 ... 0

d W ld

i

d

D

i i d

W SC A

W

, where j k l n C SC A, , { | _ }n d .

Let =

+ + +

+ + +

+ + +

×

= = =

= = =

= = =

SC pA

a a a a a a

a a a _ (

...

... ... ... ...

... ) _d

W j W k W l

W j W k W l

W W j W W k W W l

W SC A

| | 1, | | 1, | | 1,

| | 2, | | 2, | | 2,

| | | | , | | | | , | | | | ,

| | | |

i d

i i d

i i d

i

i d

i i d

i i d

i

i d

i d i d

i d i d

i d

d d

1 1

1 1

1 1

1 1

1 1

1 1

1 1

1 1

1 1

. Then SC_pAd is a set of columns in Ad.

According to DEFINITION 4.5, in SC_Ad, = = = + += =CS w a w W W W{ | 1, | | 1, ..., | | | | }n w n i d

i i d

i d, 1 1

1 1 . Obviously in SC_pAd, =CS CSnd n.

Since CSnd = CSn, SC_pAd satisfies: (1) =CS CS n n C C SC pA, , , _nd n d nd n d d; and (2) =CS w w A C SC pA{ | }, _nd d nd d. Thus SC_pAd is an existing partition of Ad.

Therefore, = =

p A p A_ _ d

D d

1 . □

Appendix D. Proof of THEOREM 5.3

THEOREM 5.3. == p A p A_ _dD d1 , where p A_ dis an arbitrary partition of Ad and p A_ is an existing partition of A.

PROOF. == p A p A p A p A p A_ _ _ ... _ ... _dD d d D1 1 2 . Let SC_Ad is the set of columns ofp A_ d. Let = =

SC SC A_ d

D d

1 .

(1) Since p A_ dis a partition of Ad, =CS CS n n C C SC A, , , _nd n d nd n d d. According to PROPERTY 5.3, =W A W A d d_ _ ,d d , so =CS CS n A d d n A, , ,nd n d d d. n Ad, =

= = +

C C[ 0 0 ... 0 0 0 ... 0 ]n W

nd W

T

| | | |i d

i i d D

i1 1

1

, so CSn=CSnd. Thus,

=CS CS n n C C SC, , ,n n n n . (2) Since p A_ dis a partition of Ad, =CS W A C SC A_ , _nd d nd d. n Ad, =

= = +

C C[ 0 0 ... 0 0 0 ... 0 ]n W

nd W

T

| | | |i d

i i d D

i1 1

1

, so CSn=CSnd. Thus,

=CS W A n A_ ,n d . Then, = = =

CS n A C SC CS n A W A, (i.e., ) , _n n d

D n d

d

D d

1 1 . == W A W_dD d1 according to PROPERTY 5.4.

Thus, =CS W C SC,n n .

Since SC satisfies the above two conditions, SC is an existing partition of A. Therefore, we can conclude == p A p A_ _dD d1 . □

Appendix E. Proof of THEOREM 6.1

THEOREM 6.1. In every Pareto-optimal solution of BGRSP (i.e., Pareto-optimal partition of matrix A), any partition of submatrix is Pareto-optimal.

PROOF. Let p_A be an arbitrary Pareto-optimal partition of matrix A. Assume that in p_A at least one submatrix's partition is not Pareto-optimal. Case 1: D = 1, A contains only one submatrix, where p_A = p A_ 1be the non-Pareto-optimal partition of A1. Letpop A_ 1be one of Pareto-optimal

partitions of A1 and np_A = pop A_ 1. The carbon emissions and average profit of rides of p A_ is p A_ 1 and O p A

p A _

_ f 12

13 , where p A_ 12and p A_ 13 is the cost

and rides number of p A_ 1 respectively. The carbon emissions and average profit of rides of np A_ is pop A_ 11 and O pop A

pop A _

_ f 12

13 . Sincepop A_ 1 p p A_ 1, then

pop A p A pop A p A pop A p A_ _ , _ _ , and _ _1 1 1 1 1 11 1 2 2 3 3 where at least one inequality is strict. Thus, if =pop A p A_ _1 11 1, then

pop A p A pop A p A_ _ , and _ _1 1 1 12 2 3 3 where at least one inequality is strict, and so > O pop A

pop A O p A

p A _

_ _

_ f f12

13

12 13

, therefore, pop A_ 1 s p A_ 1; else

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

480

<pop A p A_ _1 11 1, then pop A p A pop A p A_ _ , and _ _1 1 1 12 2 3 3, and so O pop A

pop A O p A

p A _

_ _

_ f f12

13

12 13

, therefore, np A_ s p A_ . But this contradicts the as- sumption that p_A is Pareto-optimal partition of matrix A.

Case 2: D 2, A contains at least two submatrices, i.e., = =

p A p A_ _ d

D d

1 , where p_Ad is a partition of Ad. Assume a non-Pareto-optimal partition

exists in Ad, i.e., p_Ad. The carbon emissions of p_A (p_A1) is + += = +p A p A p A_ _ _s d

s d s d D

s1 1

11 1 1 and the average profit of rides of p_A (p_A2) is + +

+ + = = +

= = +

O p A p A p A

p A p A p A

( _ _ _ )

_ _ _ f s

d s d s d

D s

s d

s d s d D

s

1 1

2 2 1 2

1 1

3 3 1 3 , where p A p A p A_ , _ , and _d d d1 2 3 p A_ d2 are carbon emissions, cost and rides number of p_Ad respectively.

Letpop A_ dbe one of Pareto-optimal partitions of Ad. Let = = = +

np A p A pop A p A_ _ _ _ s

d s d

s d

D s

1

1

1 , according to THEOREM 5.3 np_A is a partition of A.

The carbon emissions of np_A (np_A1) is + += = +p A pop A p A_ _ _s d

s d s d D

s1 1

11 1 1 and the average profit of np_A (np_A2) is + +

+ + = = +

= = +

O p A pop A p A

p A pop A p A

( _ _ _ )

_ _ _ f s

d s d s d

D s

s d

s d s d D

s

1 1

2 2 1 2

1 1

3 3 1 3 . Sincepop A_ d pp_Ad, then pop A p A pop A p A_ _ , _ _d d d d1 1 2 2 pop A p A, and _ _d d3 3where at least one in-

equality is strict. Thus, if =pop A p A_ _d d1 1, then pop A p A pop A p A_ _ and _ _d d d d2 2 3 3 where at least one inequality is strict, and so

> + +

+ +

+ +

+ + = = +

= = +

= = +

= = +

O p A pop A p A

p A pop A p A

O p A p A p A

p A p A p A

( _ _ _ )

_ _ _

( _ _ _ )

_ _ _ f s

d s d s d

D s

s d

s d s d D

s

f s d

s d s d D

s

s d

s d s d D

s

1 1

2 2 1 2

1 1

3 3 1 3

1 1

2 2 1 2

1 1

3 3 1 3 , therefore, np A p A_ _s . If <pop A p A_ _d d1 1, then

pop A p A pop A p A_ _ and _ _d d d d2 2 3 3, and so > + +

+ +

+ +

+ + = = +

= = +

= = +

= = +

O p A pop A p A

p A pop A p A

O p A p A p A

p A p A p A

( _ _ _ )

_ _ _

( _ _ _ )

_ _ _ f s

d s d s d

D s

s d

s d s d D

s

f s d

s d s d D

s

s d

s d s d D

s

1 1

2 2 1 2

1 1

3 3 1 3

1 1

2 2 1 2

1 1

3 3 1 3 , therefore, np A p A_ _s . But

this contradicts the assumption that p_A is Pareto-optimal partition of matrix A. These cases cover all the possibilities. Since the assumption that a Pareto-optimal solution of BGRSP can be composed of non-Pareto-optimal

partitions led to a contradiction, in every Pareto-optimal solution, partitions are all Pareto-optimal. □

Appendix F. Data of instance with 106 customers used in test

NO. X Y Earliest time Latest time Demand NO. X Y Earliest time Latest time Demand

1 21 24 18 28 28 54 28 18 143 153 5 2 22 22 18 28 2 55 55 54 144 154 20 3 24 12 31 41 5 56 4 18 144 154 10 4 15 10 32 42 20 57 10 43 145 155 9 5 6 38 32 42 16 58 37 31 155 165 14 6 15 30 34 44 26 59 31 67 155 165 3 7 27 69 34 44 10 60 37 31 155 165 14 8 47 16 35 45 25 61 61 52 156 166 9 9 35 40 37 47 16 62 55 60 157 167 21 10 53 52 37 47 11 63 45 10 157 167 18 11 41 37 39 49 16 64 42 7 167 177 11 12 25 24 39 49 20 65 25 30 169 179 9 13 2 60 51 61 5 66 26 27 170 180 9 14 60 12 54 64 31 67 32 12 171 181 2 15 14 37 54 64 11 68 11 31 171 181 16 16 35 17 60 70 7 69 6 68 178 188 10 17 31 52 60 70 27 70 55 45 186 196 6 18 37 47 60 70 6 71 65 55 197 207 11 19 8 56 61 71 27 72 2 48 197 207 27 20 62 77 61 71 20 73 30 60 204 214 27 21 27 43 62 72 9 74 47 47 204 214 15 22 15 47 65 75 16 75 45 65 206 216 12 23 24 58 78 88 19 76 49 73 207 217 25 24 19 21 78 88 10 77 53 12 210 220 6 25 30 5 81 91 8 78 45 30 212 222 9 26 45 20 82 92 11 79 23 3 222 232 7 27 50 35 83 93 16 80 25 21 223 233 8 28 64 42 83 93 9 81 20 20 224 234 11 29 20 65 87 97 12 82 22 27 225 235 28 30 55 5 88 98 29 83 63 23 226 236 2 31 55 5 89 99 29 84 57 29 230 240 13 32 11 14 89 99 18 85 35 69 241 251 3 33 49 11 99 109 18 86 56 39 242 252 36 34 11 14 99 109 18 87 56 39 242 252 36 35 40 60 101 111 21 88 63 65 243 253 8 36 15 77 103 113 9 89 55 20 249 259 6 37 49 42 103 113 13 90 46 13 249 259 8 38 26 52 104 114 9 91 46 13 249 259 8 39 10 20 105 115 19 92 65 35 263 273 3 40 15 60 106 116 5 93 5 30 267 277 9 41 12 24 116 126 13 94 30 25 269 279 7 42 57 68 117 127 15 95 15 19 270 280 19 43 44 17 118 128 7 96 41 49 271 281 25 44 20 50 121 131 9 97 17 34 272 282 16 45 5 5 123 133 20 98 13 52 275 285 9 46 67 5 123 133 31 99 65 20 292 302 2 47 20 26 123 133 9 100 26 35 296 306 12 48 40 25 125 135 9 101 53 43 299 309 8

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

481

49 20 40 127 137 12 102 37 56 302 312 19 50 49 58 128 138 10 103 56 37 302 312 3 51 49 58 128 138 10 104 18 18 305 315 20 52 16 22 141 151 26 105 18 24 308 318 13 53 57 48 142 152 16 106 36 26 320 330 25

The data in the Appendix is obtained by sorting Earliest time of Li & Lim PDPTW benchmark pdp_100-lr101.txt (http://www.sintef.no/projectweb/top/pdptw/li-lim- benchmark). Some Earliest time are increased, but the time window width is not changed. Xcoord and Ycoord of each customer is not changed. Demand = ABS(Demand of pdp_100-lr101), to avoid Demand is negative. The demands of each customer are delivered to drop off point with coordinate of (80, 50) in the time window. The customers whose delivery time windows overlap are possibly delivered by a vehicle. Depot (0) coordinate is (40, 50) same as that of pdp_100-lr101.

References

Amorim, P., Günther, H.O., Almada-Lobo, B., 2012. Multi-objective integrated production and distribution planning of perishable products. Int. J. Prod. Econ. 138 (1), 89–101.

Bektaş, T., Laporte, G., 2011. The pollution-routing problem. Transport. Res. Part B 45 (8), 1232–1250.

Caulfield, B., 2009. Estimating the environmental benefits of ride-sharing: a case study of Dublin. Transport. Res. Transport Environ. 14 (7), 527–531.

Cordeau, J., Laporte, G., 2003. A Tabu search heuristic for the static multi-vehicle dial-a- ride problem. Transport. Res. Part B 37 (6), 579–594.

Cordeau, J., 2006. A branch-and-cut algorithm for the dial-a-ride problem. Oper. Res. 54 (3), 573–586.

Darvish, M., Archetti, C., Coelho, L.C., 2018. Trade-offs between environmental and economic performance in production and inventory-routing problems. Int. J. Prod. Econ.

Demir, E., Bektaş, T., Laporte, G., 2014a. The bi-objective pollution-routing problem. Eur. J. Oper. Res. 232 (3), 464–478.

Demir, E., Bektaş, T., Laporte, G., 2014b. A review of recent research on green road freight transportation. Eur. J. Oper. Res. 237 (3), 775–793.

Franceschetti, A., Honhon, D., Woensel, T.V., Bektaş, T., Laporte, G., 2013. The time- dependent pollution-routing problem. Transp. Res. Part B Methodol. 56, 265–293.

Fukasawa, R., He, Q., Song, Y., 2016. A disjunctive convex programming approach to the pollution-routing problem. Transp. Res. Part B Methodol. 94, 61–79.

Furuhata, M., Dessouky, M., Ordonez, F., Brunet, M., Wang, X., Koenig, S., 2013. Ridesharing: the state-of-the-art and future directions. Transp. Res. Part B Methodol. 57, 28–46.

Garcia-Najera, A., Gutierrez-Andrade, M.A., 2013. An evolutionary approach to the multi- objective pickup and delivery problem with time windows. In: IEEE Congress on Evolutionary Computation, pp. 997–1004.

Han, B., Zhang, W., Lu, X., Lin, Y., 2015. On-line supply chain scheduling for single- machine and parallel-machine configurations with a single customer: minimizing the makespan and delivery cost. Eur. J. Oper. Res. 244 (3), 704–714.

Jozefowiez, N., Semet, F., Talbi, E.G., 2009. An evolutionary algorithm for the vehicle routing problem with route balancing. Eur. J. Oper. Res. 195 (3), 761–769.

Knopfmacher, A., Mays, M., 2006. Ordered and unordered factorizations of integers. Math. J. 10 (1), 72–89.

Koç, Ç., Bektaş, T., Jabali, O., Laporte, G., 2014. The fleet size and mix pollution-routing problem. Transp. Res. Part B Methodol. 70, 239–254.

Kovacs, A.A., Parragh, S.N., Hartl, R.F., 2015. The multi-objective generalized consistent vehicle routing problem. Eur. J. Oper. Res. 247 (2), 441–458.

Lemesre, J., Dhaenens, C., Talbi, E.G., 2007. Parallel partitioning method (ppm): a new

exact method to solve bi-objective problems. Comput. Oper. Res. 34 (8), 2450–2462. Park, Y.B., 2001. A hybrid genetic algorithm for the vehicle scheduling problem with due

times and time deadlines. Int. J. Prod. Econ. 73 (2), 175–188. Parragh, S.N., Jorge, P. de S., Bernardo, A.L., 2014. The dial-a-ride problem with split

requests and profits. Transport. Sci. 49 (2), 311–334. Psaraftis, H.N., 1983. An exact algorithm for the single vehicle many-to-many dial-a-ride

problem with time windows. Transport. Sci. 17 (3), 351–357. Sarpong, B.M., 2013. Column Generation for Bi-objective Integer Linear Programs:

Application to Bi-objective Vehicle Routing Problems. Doctoral dissertation. Institut d’Optique Graduate School.

Suzuki, Y., 2016. A dual-objective metaheuristic approach to solve practical pollution routing problem. Int. J. Prod. Econ. 176, 143–153.

Tiwari, A., Chang, P.C., 2015. A block recombination approach to solve green vehicle routing problem. Int. J. Prod. Econ. 164, 379–387.

Wang, J.W., Dou, R., Muddada, R.R., Zhang, W., 2018a. Management of a holistic supply chain network for proactive resilience: theory and case study. Comput. Ind. Eng. 125, 668–677.

Wang, J.W., Muddada, R.R., Wang, H., Ding, J., Lin, Y., Liu, C., Zhang, W., 2016. Toward a resilient holistic supply chain network system: concept, review and future direction. IEEE Systems Journal 10 (2), 410–421.

Wang, J.W., Yu, Y., Tang, J., 2018b. Compensation and profit distribution for cooperative green pickup and delivery problem. Transp. Res. Part B Methodol. 113, 54–69.

Wang, J., Zhou, Y., Wang, Y., Zhang, J., Chen, C.L., Zheng, Z., 2015. Multiobjective ve- hicle routing problems with simultaneous delivery and pickup and time windows: formulation, instances, and algorithms. IEEE Trans. Cybern. 46 (3), 582–594.

Wang, Y., Tian, Z., Feng, C., Feng, S., Zhang, P., 2012. Performance analysis of gen- eralized block diagonal structured random matrices in compressive sensing. In: International Symposium on Communications and Information Technologies, pp. 79–797.

Wang, Y., Zheng, B., Lim, E.P., 2018. Understanding the effects of taxi ride-sharing—a case study of Singapore. Comput. Environ. Urban Syst. 69, 124–132.

Xiang, Z., Chu, C., Chen, H., 2008. The study of a dynamic dial-a-ride problem under time-dependent and stochastic environments. Eur. J. Oper. Res. 185 (2), 534–551.

Yu, Y., Tang, J., Li, J., Sun, W., Wang, J., 2016. Reducing carbon emission of pickup and delivery using integrated scheduling. Transport. Res. Part D 47 (2016), 237–250.

Yu, Y., Lou, Q., Tang, J., Wang, J., Yue, X.H., 2017. An exact decomposition method to save trips in cooperative pickup and delivery based on scheduled trips and profit distribution. Comput. Oper. Res. 87 (2017), 245–257.

Zhang, S., Gajpal, Y., Appadoo, S.S., Abdulkader, M.M.S., 2018. Electric vehicle routing problem with recharging stations for minimizing energy consumption. Int. J. Prod. Econ. 203, 404–413.

Y. Yu et al. International Journal of Production Economics 208 (2019) 472–482

482