An Optimization Model of a Sustainable City Logistics Network Design Based on Goal Programming

: This paper investigates the joint optimization problem on the logistics infrastructure investment and CO 2 emission taxes for a sustainable city logistics network design by a goal programming approach where the cost recovery, service level and CO 2 emission reduction goals are involved. The above multi-objective logistics infrastructure capacity investment and CO 2 emission taxes problem is formulated as a bi-level goal programming model. Given the priority structure of the goals, the total deviations from predetermined goals are minimized in the upper level, while the lower level of the model serves as the service route choice equilibrium problem of logistics users. To solve the proposed model, a genetic algorithm is developed, where the method of successive average (MSA) is embedded. The case study focusing on the urban logistics network of Changsha, China demonstrates the effectiveness of the bi-level goal programming model and the genetic algorithm. The ﬁndings reveal that the priority rankings of the goals have a signiﬁcant impact on the joint decisions of CO 2 emission taxes and logistics infrastructure capacity investment. The proposed methodology provides an avenue to balance multiple conﬂicting objectives and obtain an economical and environmental city logistics network.


Introduction
Due to rapid urbanization in recent years, the demand of distribution in city logistics networks has increased dramatically [1,2]. Freight transportation, which plays an essential role along products distribution, is one of main sources of emissions of greenhouse gases (GHGs), e.g., carbon dioxide (CO 2 ), nitrogen oxide, sulphur oxide, and particulate matters [3]. Almost 10% of carbon emissions is attributed to freight transportation [4]. It is evident that freight transportation adversely affects the sustainability of city logistics systems [5]. Therefore, it is significant to implement effective measures (e.g., increasing the coordination between logistics participants [6]) and legislations (e.g., CO 2 emission taxes) to reduce unfavourable environmental impacts caused by increased vehicle pollution emissions so as to improve efficiency and environmental sustainability of city logistics systems.
In terms of constructing a low-carbon urban logistics system, the following effective methods are usually considered, i.e., the design of logistics network and proposals of policies for sustainable logistics management. The former mainly involves optimizing the city logistics network design, while the latter mostly refers to logistics authority effective low-carbon initiative measures, such as CO 2 emission taxes or green subsidies [7].
As we know, it is critical to optimize the logistics network in a city logistics system, as well as the timing of the investment of logistics infrastructure (e.g., logistics park, logistics centre, distribution centre) [8]. When designing a low-carbon city logistics network,

Green Logistics Network Designs
The total costs or the efficiency level of operators are the two major concerns in traditional logistics network design models, while the external environmental costs are usually overlooked [7]. In contrast, the green logistics network design models drew more attention to environment-related factors. Not only was improving the efficiency of logistics service with a decline of logistics costs prioritized, but achieving sustainable and social objectives was as well in the green logistics network design models [28][29][30]. The studies of green network design problems have been investigated intensively in the literature. Harris, Naim [31] demonstrated the necessity to address the economic and environmental objective explicitly by exploring the relationship between logistics cost and environmental effects. Zhang, Zhan [7] explored the regional green logistics network design with low-carbon initiatives under the carbon emission reduction target. Furthermore, some relevant research regarding logistics network design problems subsumed environmental factors such as CO 2 emissions as well as the interactions among different participants, i.e., logistics authorities and logistics users. Bi-level programming is widely used for solving the integrated logistics network optimization models. The upper level for logistics authorities is to determine the optimal set of policies, whereas the selection of logistics service routes of logistics users is processed in the lower level [32][33][34].
An important and effective method management measure to improve the green city logistics system is CO 2 emission taxes charges, which will benefit to switch the transport mode with higher CO 2 emission to the greener transport mode so as to improve social welfare [28,34,35]. Yin, Li [35] addressed the integrated optimization on the toll pricing and capacity of a road network considering congestions by a goal approach. Chen, Yang [36] studied the coastal liner route design problem where goals on CO 2 emission reductions and state subsidy levels are involved.

Logistics Network Design Model with Multi-Objective Optimization
To promote the sustainable urban logistics system, multiple benefit stakeholders and multiple decision objectives are involved. Nevertheless, different decision objectives conflict with each other in most cases. A high-quality solution considering one criterion may be unfavourable for other criteria. This attribute has stimulated research interest in logistics network design based on multi-objective optimization approaches in recent years.
Three methods exist: the weighted sum method, artificial intelligence approaches such as genetic algorithms and neural networks, and multi-objective programming to deal with the multi-objective optimization problems. The weighted sum method, which transforms a multi-objective problem to a single-objective problem, was commonly applied to attain Pareto-efficient solutions [15,37]. Most of the existing studies in low-carbon logistics network design field concentrated on the trade-offs between the total costs and CO 2 emissions. Wang, Lai [38] presented a mixed-integer optimization model for the multi-objective supply chain network problem, which considered the total costs and total CO 2 emissions. Harris, Mumford [39] investigated a capacitated facility location-allocation problem with flexible store allocations for green logistics by an efficient evolutionary multi-objective optimization approach, where financial costs and CO 2 emissions are considered simultaneously. Yuchi, Wang [40] proposed a bi-objective model where the economic and environmental targets are involved for the reverse logistics network design problem arising in the truck tire industry. An optimization model which minimized the total cost and CO 2 emissions was proposed to design a regional timber logistics network [41]. A regional multimodal logistics network design problem considering the uncertainty of logistics demands, the reduction goal of CO 2 emissions, subsidies and multiple stakeholders is investigated in Jiang, Zhang [42]. Considering the carbon cap and carbon cap-and-trade policies, Samuel, Venkatadri [43] introduced robustness into their deterministic mathematical model to investigate the impacts of the quality of returns on the closed-loop supply chain network. A concise but complete review on the combinatorial optimization models proposed for network design problems arising in reverse logistics and waste management supply chain was provided in [44]. Liu, Li [45] presented a review on multi-objective discrete optimization problems (MODOPs) from existing multi-objective metaheuristics for MODOPs, application areas of MODOPs, performance metrics and test instances.
Another important method to solve multi-objective decision problem is the goal programming (GP) approach introduced by Charnes, Cooper [46]. Goal programming serves as an effective and efficient tool which is capable of producing compromise solutions to satisfy the goals specified by decision-makers as much as possible [47]. The GP approach is always capable of providing a good satisfactory solution that conforms to the preference goals defined by users rather than just guaranteeing the theoretical optima, which is consistent with decision-making processes in real-world cases. In addition, unlike the set of Pareto-efficient solutions, a single solution can be offered by the GP, which caters to realistic applications. Bal and Satoglu [48] utilized the GP approach to balance economic, social and environmental targets for the logistics network design of the recovery of waste electric and electronic equipment.
To the best of our knowledge, relevant publications which address both the investment of logistics infrastructures and charges of CO 2 emissions taxes by a multi-objective decision perspective are still scarce. In addition, only a limited number of existing studies employ the GP approach to deal with multi-objective logistics network design models (shown in Table 1). Therefore, this paper applied the GP approach to investigate the joint optimization on CO 2 emission taxes and logistics infrastructure investment in a sustainable urban logistics network.

Network Representation
Several detailed descriptions are made in Figure 1 as preliminaries to model the city logistics network. Figure 1a gives the logistics demand network which records multiple types of logistics demands on OD pairs. Figure 1b represents the logistics service network which embeds multiple logistics service routes to satisfy logistics demands. Figure 1c gives three typical logistics service routes through different logistics nodes and transport modes. m e unit CO2 emission of transport mode mM  (kg/ton-km)

Network Representation
Several detailed descriptions are made in Figure 1 as preliminaries to model the city logistics network. Figure 1a gives the logistics demand network which records multiple types of logistics demands on OD pairs. Figure 1b represents the logistics service network which embeds multiple logistics service routes to satisfy logistics demands. Figure 1c gives three typical logistics service routes through different logistics nodes and transport modes. As shown in Figure 1d, a logistics infrastructure network = ( , ) comprises the logistics nodes set p N as well as the logistics links set p A . More specifically, consists of a number of logistics parks, distribution centers and freight terminals, and p A denotes the physical links between a pair of logistics nodes which includes road segments, rail tracks, and river segments. All logistics transferring actives are conducted at logistics transfer nodes. Moreover, two groups of logistics nodes, namely logistics parks with economies of scale effects and distribution centers with smaller capacity, are considered.
To describe of the model, a virtual arc is introduced, which denotes a logistics transferring operation (e.g., the changes of transport modes). As shown in Figure 2, a logistics service along pair OD is conducted by combined transportation, which embeds a transfer from a waterway to a railway at transfer node H. To simplify the presentation, we denote

Assumptions
In this study, the following five basic assumptions are underlined in the city logistics network design model. A1: The city logistics network design model is proposed to support strategic planning and policy evaluation.
A2: In the city logistics system, the logistics authority will decide the joint optimal scheme on the logistics nodes investment and CO2 emission taxes. As shown in Figure 1d, a logistics infrastructure network G p = N p , A p . comprises the logistics nodes set N p as well as the logistics links set A p . More specifically, N p consists of a number of logistics parks, distribution centers and freight terminals, and A p denotes the physical links between a pair of logistics nodes which includes road segments, rail tracks, and river segments. All logistics transferring actives are conducted at logistics transfer nodes. Moreover, two groups of logistics nodes, namely logistics parks with economies of scale effects and distribution centers with smaller capacity, are considered.
To describe of the model, a virtual arc is introduced, which denotes a logistics transferring operation (e.g., the changes of transport modes). As shown in Figure 2, a logistics service along pair OD is conducted by combined transportation, which embeds a transfer from a waterway to a railway at transfer node H. To simplify the presentation, we denote N h s ⊂ N s as the set of logistics transfer nodes and A h s ⊂ A s as the set of virtual logistics transfer arcs.

Network Representation
Several detailed descriptions are made in Figure 1 as preliminaries to model the city logistics network. Figure 1a gives the logistics demand network which records multiple types of logistics demands on OD pairs. Figure 1b represents the logistics service network which embeds multiple logistics service routes to satisfy logistics demands. Figure 1c gives three typical logistics service routes through different logistics nodes and transport modes. As shown in Figure 1d, a logistics infrastructure network = ( , ) comprises the logistics nodes set p N as well as the logistics links set p A . More specifically, consists of a number of logistics parks, distribution centers and freight terminals, and p A denotes the physical links between a pair of logistics nodes which includes road segments, rail tracks, and river segments. All logistics transferring actives are conducted at logistics transfer nodes. Moreover, two groups of logistics nodes, namely logistics parks with economies of scale effects and distribution centers with smaller capacity, are considered.
To describe of the model, a virtual arc is introduced, which denotes a logistics transferring operation (e.g., the changes of transport modes). As shown in Figure 2, a logistics service along pair OD is conducted by combined transportation, which embeds a transfer from a waterway to a railway at transfer node H. To simplify the presentation, we denote

Assumptions
In this study, the following five basic assumptions are underlined in the city logistics network design model. A1: The city logistics network design model is proposed to support strategic planning and policy evaluation.
A2: In the city logistics system, the logistics authority will decide the joint optimal scheme on the logistics nodes investment and CO2 emission taxes.

Assumptions
In this study, the following five basic assumptions are underlined in the city logistics network design model. A1: The city logistics network design model is proposed to support strategic planning and policy evaluation.
A2: In the city logistics system, the logistics authority will decide the joint optimal scheme on the logistics nodes investment and CO 2 emission taxes.
A3: Logistics demand originated from industrial and commercial activities are involved. The customers' points of industry logistics demand come from the peripheral of the city, while those of commercial logistics demand arises from the city center.
A4: Logistics users tend to make their decisions for logistics service routes based on the logistics service disutility which can be evaluated by the service time and monetary cost. The monetary cost refers to the logistics cost and the CO 2 emission cost [35].
A5: The introduction of road-user charging can favor a switch to other greener transport modes, such as railway and waterway. Thus, the CO 2 emissions tax charges are imposed on the transport modes with high unit emissions [29].

Model Formulation
In the city logistics system, the logistics authority as well as logistics users (i.e., carriers) are the two interrelated decision makers. The logistics authority has to decide the optimal logistics infrastructure investments and CO 2 emission taxes charges, so as to realize their objectives as much as possible. Three types of logistics infrastructure investment options are considered in this study. The first two options are to add new physical links or improve existing physical links in current logistics network, and the third one is to locate several logistics terminals in the network. The main objectives specified by the logistics authority are shown as below.
(1) Cost-recovery goal: The total logistics investment cost on the capacity enhancements of logistics infrastructures is covered by the total CO 2 emissions taxes charges in the city logistics network partially or completely.
(2) Service-level goal: The total generalized logistics cost of the city logistics system is affected by the freight flow of the logistics network as well as the decisions made by the logistics authority. After introducing the logistics infrastructure investment and CO 2 emissions taxes charges, a decline of the total generalized logistics cost is expected to exceed a pre-specified level.
(3) Environmental goal: The logistics authority requires that the CO 2 emission decreases by a pre-specified level after the CO 2 emission taxes charging and logistics capacity investment schemes are implemented.
Since above three objectives conflict with each other, favoring one of them usually weakens others. The GP approach is an effective method for us to obtain a good compromise solution among these three conflicting objectives. The goal constraints for corresponding objectives are introduced in Section 4.1.
In addition, the decisions of logistics users and the logistics authority are mutually affected by each other. From the perspective of logistics users, their decisions, which are evaluated by the perceptions of logistics service disutility, are influenced by the decisions made by the logistics authority. Conversely, the logistics authority's decisions are affected by the routes and transport mode choices made by logistics users. These interactions can be described as a bi-level model, which can be summarized in Figure 3. We provide more detailed descriptions about the goal programming decision model based on bi-level programming in following sections.

Model Formulation
In the city logistics system, the logistics authority as well as logistics users (i.e., carriers) are the two interrelated decision makers. The logistics authority has to decide the optimal logistics infrastructure investments and CO2 emission taxes charges, so as to realize their objectives as much as possible. Three types of logistics infrastructure investment options are considered in this study. The first two options are to add new physical links or improve existing physical links in current logistics network, and the third one is to locate several logistics terminals in the network. The main objectives specified by the logistics authority are shown as below.
(1) Cost-recovery goal: The total logistics investment cost on the capacity enhancements of logistics infrastructures is covered by the total CO2 emissions taxes charges in the city logistics network partially or completely.
(2) Service-level goal: The total generalized logistics cost of the city logistics system is affected by the freight flow of the logistics network as well as the decisions made by the logistics authority. After introducing the logistics infrastructure investment and CO2 emissions taxes charges, a decline of the total generalized logistics cost is expected to exceed a pre-specified level.
(3) Environmental goal: The logistics authority requires that the CO2 emission decreases by a pre-specified level after the CO2 emission taxes charging and logistics capacity investment schemes are implemented.
Since above three objectives conflict with each other, favoring one of them usually weakens others. The GP approach is an effective method for us to obtain a good compromise solution among these three conflicting objectives. The goal constraints for corresponding objectives are introduced in Section 4.1.
In addition, the decisions of logistics users and the logistics authority are mutually affected by each other. From the perspective of logistics users, their decisions, which are evaluated by the perceptions of logistics service disutility, are influenced by the decisions made by the logistics authority. Conversely, the logistics authority's decisions are affected by the routes and transport mode choices made by logistics users. These interactions can be described as a bi-level model, which can be summarized in Figure 3. We provide more detailed descriptions about the goal programming decision model based on bi-level programming in following sections.

Specification of the Goal Constraints
• Cost-recovery constraint The cost-recovery constraint implies that the total CO 2 emissions taxes charges covers a pre-specified percentage of the total investment cost on expanding the capacities of the logistics infrastructures, which is defined as, The left-hand of inequality (1) is the investment cost, while the right-hand represents the total CO 2 emissions taxes charges. Parameter φ denotes the percentage of the total investment cost which is covered by the total CO 2 emissions charges. If φ is smaller than 1.0, inequality (1) is a partial cost-recovery constraint. If φ equals 1.0, the total investment cost becomes fully covered by CO 2 emissions charges. The residual part of the investment cost, that is, is subsidized by other revenue. The level of capacity improvement and the CO 2 emissions charges determine that whether inequality (1) is satisfied. The logistics authority hopes that inequality (1) will be satisfied as closely as possible. In order to balance the left-and right-hand of inequality (1), two variables are defined to denote positive and negative deviations. Subsequently, we can obtain the equivalent constrain as Equation (2): where the overachievement and the underachievement of the cost-recovery goal are represented by d + 1 and d − 1 , respectively. • Service-level constraint The service-level constraint states that the changes of total generalized logistics cost of city logistics system, before and after introducing the joint schemes on logistics infrastructure investment and CO 2 emissions taxes charges, should reach an expected level, which is shown as follows.
The generalized logistics cost can be calculated according to the lower-level decision model of logistics users in Equations (11)- (14). Equation (3) requires that the total generalized logistics cost after carrying out joint optimization schemes of investment and CO 2 emissions taxes is at most α 1 times before the joint schemes.
By introducing positive and negative deviational variables d + 2 and d − 2 of total generalized logistics cost regarding the corresponding target value, the goal constraint on the improvement of service level of city logistics system can be expressed as: The environmental constraint states that the unit of CO 2 emission of the city logistics network before and after the joint schemes of logistics infrastructures investment and CO 2 emission taxes charges is no more than a predefined threshold, namely, the total amount of the CO 2 emissions after the joint scheme never exceeds a desired level without the joint scheme. This can be expressed as: By introducing positive and negative deviational variables d + 3 and d − 3 of the CO 2 emissions per unit shipment for the associated target value, the constraint (5) can be represented as Equation (8): The parameter α 2 is a positive constant that specifies the decrease threshold of the CO 2 emission of unit shipment.

Specification of the Lower-Level Decision Model of Logistics Users
As aforementioned, the logistics users' disutility is estimated by the logistics service time, transportation cost and CO 2 emission taxes, which can be defined as Equation (9): x k a , y m a ) + e m l a y m a , ∀a ∈ A, m ∈ M, k ∈ K (9) Considering the different transportation modes, the logistics service time on link a for different modes is measured in Equation (10). A US Bureau of Public Roads-type function is applied to compute service time for HGVs and LGVs [33,49,50], while the estimation of service time spent by railways and waterways consider both the transport service time and the shift interval.
Note that t m0 a is the free-flow transport service time, and t md is the average shift interval for transport mode m.
For each given logistics infrastructure investment options and CO 2 emission charges (namely given the vector x k a , y m a , ∀a ∈ A, m ∈ M, k ∈ {1, 2} , the following model is solved to obtain the equilibrium shipments for logistics users in city logistics network [51]:

Specification of the Upper-Level Multi-Objectives Decision Model of Logistics Authority
The objective function of the GP model consists of undesired absolute deviations from the specified goal levels, where both underachievement and overachievement values are considered. The GP model is aimed to minimize these deviations under a predefined priority structures of the different goals [49].
In this paper, the priority factors of the three goals follow a decreasing order, namely, the ith goal takes precedence over the (i + 1)th goal, which is represented as follows: where P i is the ith priority factor, and the operator of "<<" represents the relative importance among various objectives of logistics authority based on the priority of decision structure. Hence, the upper-level multi-objectives decision problem can be formulated as follows: and is subject to: x k a = {0, 1} , ∀a ∈ A, k ∈ K (21) where the arc flow vectors v = (v m a , ∀a ∈ A, m ∈ M) are obtained by solving the lower-level decision model of logistics users, which is the above UE model shown as in Equations (11)-(14).
The objective function Equation (16) aims to minimize the total positive deviations from the three goals, which states that the overachievement of each goal constraint should be as low as possible. Equations (17)- (19) are the goal constraints previously defined. Equation (20) means that the CO 2 emissions tax charge is non-negative and does not exceed the upper bound of charges. Equation (21) implies that logistics investment option is a 0-1 variable. Equation (22) demonstrates that at least one variable of a pair of positive and negative deviational variables is zero. Equation (23) means that the deviational variables are non-negative.

Solving Algorithm
This section provides an analysis on the corresponding algorithm, which solves the above proposed bi-level decision model based on goal programming. A hybrid genetic algorithm (HGA) embedded with the method of successive average (MSA) is given based on the characteristic of proposed model [50,[52][53][54][55][56].
(1) Initialize the important parameters of GA, which consist of the population size (pop_size), the probability of crossover (Pc) and the probability of mutation (Pm).
(2) The stopping criterion of the hybrid algorithm is defined as a certain number of generations (Max_gen) is reached.
(3) Initialize the iteration number n as zero and generate the initial population Pop n = G n k | k = 1, 2, · · · , pop_size} . Let the individual G n k be the k-th decision scheme in the n-th generation; that is, G n k = X n k , Y n k ) is the decision vector where X n k , Y n k denote the chosen investment options and CO 2 emissions tax charges, respectively.
Step 2. Calculate arc flow vectors. Given a joint scheme G n k = X n k , Y n k ) , the flow v m a of arc a by transport mode m can be obtained by solving the lower-level UE problem in Equations (9)-(14) based on the method of successive average (MSA).
Calculate the value of the objective function of upper-level decision model (Z U ) n k on the basis of Equations (16)- (23).
The fitness value f it n k of each G n k of Pop n is calculated based on the objective function value (Z U ) n k .
Sort the value fitness f it n k of individual by descending order among population pop n . Step 6. GA Operations. Three operators including selection, crossover and mutation based on the value of f it n k are conducted to generate a new population pop n . Afterwards, n is increased by 1.
Check the termination condition. If the iteration number reaches max_gen, then the algorithm terminates and outputs the final solution (X * k , Y * K ), v m a * } . Otherwise, go to step 2.

Key Operators Design of Genetic Algorithm
Since the design of GA operators and the settings of parameters have a great impact on the performance of the algorithm, key elements including the coding method, selection method and crossover and mutation methods are elaborated as follows. Step 6. GA Operations. Three operators including selection, crossover and mutation based on the value of n k fit are conducted to generate a new population n pop . Afterwards, n is increased by 1.
Check the termination condition. If the iteration number reaches max_gen, then the algorithm terminates and outputs the final solution Otherwise, go to step 2.

Key Operators Design of Genetic Algorithm
Since the design of GA operators and the settings of parameters have a great impact on the performance of the algorithm, key elements including the coding method, selection method and crossover and mutation methods are elaborated as follows.

Coding Method
As shown in Figure 4, there exists two parts in the chromosome of the individual.  Part 1 means the selection of the candidate investment projects on infrastructures, in which the value of 1 represents the corresponding logistics infrastructure will be invested, otherwise, the value is 0. The values of y1 and y2 means the unit CO 2 emissions taxes charges for HGV and LGV, respectively.

Selection Method
In this study, we employed the method of roulette wheel during the selection process. The rank-based evaluation function Eval(V i ) = a(1 − a) i−1 , ∀i = 1, 2, 3, · · · , pop_size is also adopted.
It should be pointed out that populations are sorted in the descending order based on their fitness, i.e., the first individual (i.e., i = 1) has the best fitness, whereas the last (i.e., i = pop_size) has the worst one. Additionally, the sum of the value of all individuals is equal to 1, i.e., ∑ pop_size i=1 The selection is conducted based on spinning the roulette wheel, which runs pop_size times. Chromosomes are selected to generate a new population in four steps as below.
Step 3. The i-th chromosome Step 4. Steps 2 and 3 are iteratively applied pop_size times to obtain chromosomes as parents for next generation.

Crossover and Mutation Operators
(1) Crossover operator After selecting the operator, a new generation will be obtained by crossover and mutation operations. Considering the features of the chromosome in this study, we applied two different crossover methods. As shown in Figure 5, the two-point crossover operator is implemented for Part 1, while a convex combination method is applied for Part 2. For a given pair of parents Y 1 and Y 2 , we can obtain two children C 1 and C 2 by the crossover oper- i , where λ 1 , λ 2 > 0, λ 1 + λ 2 = 1. A representation crossover operation is shown in Figure 5.

Crossover and Mutation Operators
(1) Crossover operator After selecting the operator, a new generation will be obtained by crossover and mutation operations. Considering the features of the chromosome in this study, we applied two different crossover methods. As shown in Figure 5, the two-point crossover operator is implemented for Part 1, while a convex combination method is applied for Part 2. For a given pair of parents Y1 and Y2, we can obtain two children C1 and C2 by the crossover operator, which is shown in Figure 5.
(1) + 2 (2) and 2 = 2 (1) + 1 (2) , where 1 , 2 > 0; 1 + 2 = 1. (2) Mutation operator Let r be a random positive real number no more than 1. If r < Pm, a mutation is implemented. Given a parent, = ( , ) is chosen for mutation. Figure 6a depicts the twopoint mutation operation for vector X where the genes from position 1 to position 2 in the (2) Mutation operator Let r be a random positive real number no more than 1. If r < Pm, a mutation is implemented. Given a parent, V = (X, Y) is chosen for mutation. Figure 6a  (2) Mutation operator Let r be a random positive real number no more than 1. If r < Pm, a mutation is implemented. Given a parent, = ( , ) is chosen for mutation. Figure 6a depicts the twopoint mutation operation for vector X where the genes from position 1 to position 2 in the chromosome are reversed. With regard to the second part (i.e., vector Y), we try to update as + ⋅ . If + ⋅ is infeasible, then we replace as ′ , which is randomly selected within the interval [0, ] until the new individual is feasible. If the above operations fail to achieve a feasible solution over certain number of iterations, is set as zero. Consequently, the new child is ′ = + ⋅ as shown in Figure 6b. The mutation operation for vector Y.

Data Input
In this section, a real-world case based on the logistics network design of Changsha City, China is employed to evaluate the proposed model and solution algorithm. There are 129 links, 70 nodes as well as 32 OD pairs of logistics demand in the city logistics network. The corresponding data of the logistics network is provided in Appendix A.
Two types of logistics demand originated from industry and commercial operations are included in this logistics network. Since the logistics demands of the corresponding 32 OD pairs are unavailable, we generate this missing information by an aggregation method according to the data related to the major indicators of the National Economy and Social Development published by the Changsha Statistical Bureau. The logistics demand between each of the 32 OD pairs is provided in Table 2. Table 2. Logistics demand OD pairs (10,000 tons/week). Destination  35  36  37  38  39  40  41  42  52  19  13  18  11  42  41 28 24

Data Input
In this section, a real-world case based on the logistics network design of Changsha City, China is employed to evaluate the proposed model and solution algorithm. There are 129 links, 70 nodes as well as 32 OD pairs of logistics demand in the city logistics network. The corresponding data of the logistics network is provided in Appendix A.
Two types of logistics demand originated from industry and commercial operations are included in this logistics network. Since the logistics demands of the corresponding 32 OD pairs are unavailable, we generate this missing information by an aggregation method according to the data related to the major indicators of the National Economy and Social Development published by the Changsha Statistical Bureau. The logistics demand between each of the 32 OD pairs is provided in Table 2. Table 2. Logistics demand OD pairs (10,000 tons/week). Destination  35  36  37  38  39  40  41  42   52  19  13  18  11  42  41  28  24  53  29  19  28  16  64  61  42  35  54  14  9  14  8  32  31  21  18  55  34  22  32  19  74  71  49  41 (Source: Data adapted from references [28,49]).
The VOT is USD $10 per hour [57,58]. The transfer cost and time at the logistics transfer nodes are set according to [3,59]. The average unit emissions of HGVs, LGVs, railways and waterways are 0.132, 0.283, 0.022 and 0.016 (kg/ton-km), respectively [3,60]. The parameters of the GA are given as below. The crossover rate (Pc) and the mutation rate (Pm) are set to 0.5 and 0.1, respectively. The maximal number of generations (i.e., stopping criterion) is set to 150. These input data are considered as the base case in the following analysis unless otherwise specified.

Convergence of the Hybrid Genetic Algorithm Based on Frank-Wolf
The evolution progress of the proposed algorithm is depicted as Figure 7. The proposed algorithm converges after 150 generations. The average CPU searching time for each scenario is about 1015 s. Following observations are made from Figure 7 where the convergence processes of the GA under those three scenarios are depicted. First, a dramatic decline of the logarithmic values of the objective function Z U (x, y, d) under three different scenarios within about 40 iterations is shown in Figure 7. After 50 iterations, the logarithmic values of Z U (x, y, d) under all scenarios become stable. In addition, the logarithmic value of Z U (x, y, d) with the highest setting of priority factor p 2 converges faster than other scenarios. Furthermore, the logarithmic values of Z U (x, y, d) with three different priority structures can converge to different levels. Overall, these results reveal that the settings of priority factor P i have a large impact on the converged objective value while its influence of solution time can be neglected.

Effects of Solution Structure of Goals on Model
In this section, we will further assess the impacts of the priority structures of goals on optimal decisions of green city network design. The following three scenarios are discussed in this study. Scenario 1 (T1): Cost-recovery goal >> Service-level goal >> Environmental goal; Scenario 2 (T2) Service-level goal >> Environmental goal >> Cost-recovery goal; and Scenario 3 (T3) Environmental goal >> Service-level goal >> Cost-recovery goal (scenarios are shown in Table 4). The cost-recovery and service-level goal are considered as the primary goal in Scenarios T1 and T2 respectively, while Scenario 3 considers the environmental goal as the primary goal.
From Table 4, we can see that the solutions with these three different priorities structures differ each other considerably. This observation suggests that the optimal solution is largely dependent on the priority structures of goals. Given the candidate logistics infrastructure investment options and the thresholds of the cost-recovery goal, service-level goal and environmental goal, it can be seen that the investment projects numbers and total equivalent uniform per week cost are different. That is, 11 projects with 14,880 USD/week under T1, 12 projects with 16,500 USD/week under T2 and 13 projects with 17,700 USD/week under T3. The corresponding unit shipment CO 2 reduction rates are 26.6%, 35.5% and 40.7% under T1, T2 and T3, respectively. The unit CO 2 emissions tax charges for HGV and LGV under T3 are 0.275 and 0.252, more than that of T2 (0.261 and 0.239) and T1 (0.272 and 0.241). This means that a higher priority rank in an environmental goal would result in much more improvement in unit shipment CO 2 emissions, which also requires higher more investment costs and higher CO 2 emissions tax charges. A comparison on the deviational variables of scenarios T1, T2 and T3 to analyse the influence of ranking structures of the goals on the achievement levels is shown in Table 5. In Table 5, a zero value of a deviational variable implies that the associated goal is fully achieved, whereas a positive value indicates the overachievement level of the corresponding goal. Furthermore, it is obvious that all second priority goals are achieved exactly for all three different scenarios. Moreover, the environmental goal has a deviation of 8.4% from its goal level under scenario T1. When the service-level goal is regarded as the primary (i.e., T2) and the service-level goal and environmental goal can be exactly achieved, while the overachievement level of the cost-recovery goal is 17.3% of its target value. In regard to the priority structure in scenario T3, the last two goals are both completely achieved, while the overachievement level of the cost-recovery goal is 22.7% from its target value. The above observations implies that goal priority structures have a large impact on the achievements level of each goal.

Conclusions and Future Research
In this paper, an optimization model based on GP approach is proposed to address the joint optimization on logistics infrastructure investment and CO 2 emissions taxes setting problem in the city logistics network, in which three decision objectives (i.e., cost-recovery objective, service-level objective and environmental one) are considered. We formulate this multi-objective optimization problem as a bi-level GP model. The upper level of the model is a GP formulation, whereas the lower level denotes a logistics service route choice user equilibrium problem. To solve the bi-level GP model, a hybrid genetic algorithm embedded with the method of successive average (MSA) is designed. The optimization model is empirically verified by a practical case study of the sustainable city logistics network in Changsha City, Hunan Province, China.
The computational experiments demonstrate the effectiveness and efficiency of the genetic algorithm for solving the bi-level GP model. The logarithmic values of the objective function under all scenarios can converge within 50 iterations. In addition, results of the case study show that the priority structures of the goals have a significant impact on decisions, as well as the achievement level of each specified goal. Higher priority ranking of the environmental goal or service-level goal can lead to a larger improvement at logistics system performance in terms of total generalized cost and unit shipment CO 2 emissions. Meanwhile, the corresponding schemes also require more investment and higher CO 2 emissions taxes charges. The achievement level of each goal will be conducive to improve the decision-making process by adjusting the priority structures for different goals. The proposed model and algorithm in this study provides a powerful decision support for designing a sustainable city logistics network.
Some potential directions for further research are as follows. First, some other case studies from real-world logistics networks can further verify the results of this paper and the performance of the proposed model. Second, various random factors exist in real-life cases, such as the logistics demand fluctuation, which plays a key role in the design of logistics network and management policies making. Consequently, it is essential to develop a stochastic model or robust model to cope with uncertainties in the city logistics system.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
Sets M set of all transport modes, including the heavy goods vehicles (HGVs), light goods vehicles (LGVs), railways and waterways represented by "1", "2", "3" and "4", respectively W set of origin-destination (OD) pairs in the logistics network set of logistics service routes between OD pair K set of logistics infrastructure investment options. "K = 1" represents improving capacities of existing infrastructures and "K = 2" represents adding new infrastructures A set of arcs, A = A 0 ∪ A 1 A 0 set of existing arcs without logistics infrastructure investment A 1 set of arcs after logistics infrastructure investment General Variables q w the logistics demand between OD pair w ∈ W (tons/week) t m a transport time over arc a ∈ A M s by transport mode m (hour) µ m a disutility on arc a ∈ A M s by transport mode m δ m ar indicator variable which implies whether link a is on route r by transport mode m h m wr freight flow on route r ∈ R w between OD pair w ∈ W by transport mode m (tons/week) v m a freight flow on logistics service arc a ∈ A by transport mode m (tons/week) Decision Variables x k a binary variable which equals to 1 if the logistics infrastructure investment option k over arc a is chosen and to 0 otherwise X vector of binary variable x k a , X = (x k a , ∀a ∈ A, k ∈ K) y m a emission taxes per unit CO 2 emission of transport mode m over link a ($/kg) Y vector of variable y m a , Y = (y m a , ∀a ∈ A, m ∈ M) Constants q m w potential demand between OD pair w ∈ W by transport mode m (tons/week) t m0 a free freight flow transport service time on arc a ∈ A by transport mode m (hour) t md average shift interval of transport mode m (hour) l m a length of arc a ∈ A by transport mode m (km) c m a unit transport cost on arc a ∈ A served by mode m ($/ton-km) Cap m a the maximum capacity on arc a ∈ A by mode m (tons/km) I k a the fixed cost of potential logistics infrastructures project k over link a ($/ton) τ vot value of time (VOT) ($/ton-hour) ch max the maximum of unit CO 2 emissions taxes charge ($/kg) e m unit CO 2 emission of transport mode m ∈ M (kg/ton-km)