A Bi-Level Programming Model for the Railway Express Cargo Service Network Design Problem

Service network design is fundamentally crucial for railway express cargo transportation. The main challenge is to strike a balance between two conflicting objectives: low network setup costs and high expected operational incomes. Different configurations of these objectives will have different impacts on the quality of freight transportation services. In this paper, a bi-level programming model for the railway express cargo service network design problem is proposed. The upper-level model forms the optimal decisions in terms of the service characteristics, and the low-level model selects the service arcs for each commodity. The rail express cargo is strictly subject to the service commitment, the capacity restriction, flow balance constraints, and logical relationship constraints among the decisions variables. Moreover, linearization techniques are used to convert the lower-level model to a linear one so that it can be directly solved by a standard optimization solver. Finally, a real-world case study based on the Beijing–Guangzhou Railway Line is carried out to demonstrate the effectiveness and efficiency of the proposed solution approach.


Introduction
In recent years, the competition between railway and highway transportation has become increasingly intense, especially for the long-distance transportation of high-value freight. The governments expect that more freight flow transported by highways should be diverted to the railways in order to reduce carbon emissions. For example, in Europe, 30% of road freight over 300 km is expected to shift to railways or waterways by 2030. The Ministry of Transport of China also recommends that the railways and waterways should support more freight transportation [1]. However, road freight transport is convenient and efficient currently, and many shippers tend to choose the road transportation mode. Therefore, as one of the largest and busiest railway systems, the China Railway has focused on the development of scheduled train services during past years in an effort to obtain a higher market share in the high-value freight transportation. Among them, the Railway Express Cargo Service Network Design problem (RECSNDP) is the key issue that needs to be addressed in both theory and practice.
The RECSNDP requires the determination of the train operation plans and a set of shipment delivery strategies on the network, which satisfies target customer demand with an acceptable level of service at minimum cost, and without violating capacity restrictions. The prominent features of it are the focus on fast transport, guaranteed pickup and delivery times for customers.
The RECSNDP is normally viewed as a tactical planning problem in which the railway company has to decide which terminals should provide with direct transportation services and at what frequency.  Figure 1 presents a simple railway network. The network contains six stations, named A~F, along the railway line. We generate a potential train operation plan according to the shipment characteristics on the railway network, which includes 12 trains that are denoted by the lines with different styles and colors (see Figure 1). Different line styles represent different levels of speed, i.e., a solid line represents 80 km/h, a dashed line represents 120 km/h, and a dot-and-dash line represents 160 km/h. In this way, a service network is constructed, which includes twenty-one service arcs. Note  characteristics on the railway network, which includes 12 trains that are denoted by the lines with different styles and colors (see Figure 1). Different line styles represent different levels of speed, i.e., a solid line represents 80 km/h, a dashed line represents 120 km/h, and a dot-and-dash line represents 160 km/h. In this way, a service network is constructed, which includes twenty-one service arcs. Note that all the arcs differ from each other even they belong to the same train. For example, though the service arc A01 and service arc A02 belong to the same train, they are regarded as different service arcs. Moreover, even the origin station and destination station of two arcs are the same, they may be different arcs, such as arc A01 and arc A08. In order to express the service network conveniently, we convert the trains into a matrix. In this way, each train can be described by the service network (see Table 1). Meanwhile, all of the service arcs can be identified from the matrix. Table 1. Train operation plan.

Train
Speed In addition, the service arc can be rewritten as (the origin station, the destination station, travel time, train ID). For example, if the distance between A and B is 500 km, then the arc A→B provided by train 1 can be expressed as (A, B, 500/80, 1). In practice, there may exist a series of train operation plans that match the shipment characteristics and we can generate these plans using above mentioned method. Clearly, one plan corresponds to a service network.
The second step is to assign the shipments to the service network based on the train operation plan. For example, shipment 1: A→E is considered on the complete service network. The possible transportation strategies that can be adopted by the carrier include (but are not limited to) the following: In Strategies (1) and (2), the shipment need to wait with the train. Whereas in Strategies (3)~(8), the shipment need to be reclassified once. In contrast, a waiting (at B) and once reclassification operations (at C) is needed in Strategy (9). Waiting on the same train will take one to two hours, but the reclassification will take almost eight hours. The total costs of a service network can be expressed as the sum of the service network setup costs and cargo reclassification delay costs. The first term in the total costs comes from the train operation and the second term comes from the shipments moving over the service network. For the train operation, there exists a fixed cost at the origin station and variable costs involving the operating mileages and intermediate stops. The express cargos might be transferred from one train to another train during its itinerary. The reclassification will increase the workloads of intermediate stations.
We use the cargo reclassification delay costs to calculate the workloads. The detailed calculation process is shown in the latter sections.
The express train service is usually of high cost, fast speed and low frequency, while the general train service has a lower speed at a lower cost and is always provided more frequently. Providing more express train service will increase fixed costs, but the general train service cannot guarantee the prescribed transit time; although non-stop train service save the time for shipments, it requires more train services. Furthermore, selecting the optimal transportation strategies for all the shipments while respecting the prescribed transit time constraints and capacity constraints is an optimization problem [24]. In conclusion, the RECSNDP is a complicated combinatorial optimization problem, which is generally not possible to get the optimality for real-life instances.

Mathematical Model
This section aims to provide a mathematical description for the rail express cargo service network design problem. To facilitate the model formulation, we make the following assumptions throughout this paper: Assumption 1. In this paper, we research the express cargo service network between the railway hubs and ignore the detailed operations within the hub. Each hub including marshalling stations, logistics centers and other stations can be viewed as a station.

Assumption 2.
Each shipment should not be split during the transportation process. This means that each shipment can only choose a single route (a train service chain). However, railway operators can decide whether to transport the whole volume of a shipment or just a portion.

Notations
The indices, sets, input parameters and decision variables used in the optimization model are listed in Table 2.  The fixed costs of the train (i,j,l,P), where i denotes the origin station, j denotes the destination station, l denotes the speed level, and P denotes the stopping patterns; Binary parameters, if the service arc m is provided by the train (i,j,l,P), α lP ijm = 1, otherwise, α lP ijm = 1; The costs associated with stopping at hub k for a train; R g Unit transportation income obtained by railway operators from carrying goods g; ξ g Continuous variable, the proportion of shipment g that can be transported.

Mathematical Formulation
The objective function and constrained conditions can be expressed with a bi-level programming formulation as follows: where f lP ij and Z(Y) are solved by the lower-level program; ∀m ∈ E, l and P come from the UP (5) f lP ij ∈ R ∀i, j ∈ V, l and P come from the UP (13) where y lP ij come from the upper-level program. The objective of upper-level Formulation (1) is to minimize the sum of the costs of train operation and maximize the sum of the income obtained by the railway operator for carrying the goods. The first term in the objective function is the total operation cost of all express trains that depend on the value of y lP ij . The second term is the income solved by the lower program. Moreover, decision variable y lP ij domains are specified by Constraints (2). Finally, Constraints (3) state the logic relationship between decision variables y lP ij and f lP ij . The objective of lower-level Formulation (4) aims to maximize the transportation income of the railway operator, which is also the second term of Formula (1). The income is equal to the total transportation fees paid by the shippers minus the total handling (i.e., transfer) costs. Constraints (5) guarantee that the total volumes on service arcs are less than the capacity provided by the train operation plan from the upper-level formulation. Constraints (6) ensure that the time consumption of transporting the cargo is less than prescribed transit time requested by the shippers. Constraints (7) describe the logical relationship between the two groups of decision variables. Constraints (8)~(10) are the well-known flow balance constraints. Specifically, Constraints (8) are the flow balance constraints at the departure nodes, Constraints (10) are the flow balance constraints at the arrival nodes, and Constraints (9) are the flow balance constraints at the intermediate (passing) nodes. Finally, decision variable domains are specified by Constraints (11), (12) and (13).
The lower-level program of above mathematical model is inherently nonlinear and is difficult to solve (to optimum). For the solution of nonlinear models, heuristic algorithms, e.g., simulated annealing algorithm or genetic algorithm, are generally used to find the (approximate) optimal solution. However, in this bi-level model, the lower-level program should serve as a solid aid for the upper-level program, which requires the lower-level program to obtain an optimal flow distribution plan quickly and accurately. In this case, the heuristic algorithm is unable to meet the demands. Therefore, the linearization techniques are adopted to convert the nonlinear model into a Linear Program (LP) so that exact solution approaches can be used to solve it to its global optimality quickly.

Linearization Technique
Clearly, Constraints (5), (6), (8), (9) and (10) are nonlinear constraints because they involve the production of two decision variables. For the purpose of linearization, we need to introduce an auxiliary continuous decision variable r m g and an auxiliary binary decision variable s mn g . The definition of them as follows: s mn g = x m g x n g ∀g ∈ G, m ∈ E, n ∈ E, t m = s n In this way, Constraints (5), (6), (8), (9) and (10) can be converted to: ∀m ∈ E, l and P come from the UP (16) In addition, we add the following constraints into the original lower-level model: Accordingly, the objective of lower-level model (4) can be rewritten as follows: After above process, there still exists a nonlinear term in Formula (26), i.e., s mn g ξ g . Similarly, we need to introduce another auxiliary continuous decision variable t mn g to linearize this term. The definition of it is as follows: As a result, the Formula (27) can be converted into: Furthermore, the following constraints should also be satisfied: By replacing ξ g x m g with r m g in Constraints (5), (8), (9) and (10), replacing x m g x n g with s mn g in Constraints (6), replacing s mn g ξ g with t mn g in Constraints (26), and adding Constraints (21)~(25) and (29)~(31) to the lower-level model, the original non-linear model can be transformed to a linear one, which means we can directly use a standard optimization solver (e.g., CPLEX or Gurobi) to solve the lower-level model.
The solution strategy of the bi-level model is as follows. Firstly, because of the upper-level model is easier to get the reasonable options owe to the few constraints it contains, we can generate a number of train operation plans based on the practical experience. Then, we substitute each plan into the lower-level program to get the objective function value of the whole bi-level programming. Finally, by comparing the objective function values among all these plans, we select the best solution.

Numerical Example
Here, we make a numerical example analysis based on the railway corridor called "Beijing-Guangzhou Railway Line" (see Figure 2), which is one of the most important north-south rail lines in China. The route, from Beijing to Guangzhou, has a total length of 2324 km. It connects six provincial capitals or municipalities and many large-and medium-sized cities, and it is one of the busiest major railways in China. The goods transported to the south on the Beijing-Guangzhou railway mainly include coal, steel, petroleum, timber and export materials. Relatively speaking, the proportion of high value-added goods to the north is higher than to the south. Therefore, this paper uses the northward direction of the Beijing-Guangzhou railway as an example for analysis. This route passes through many important railway hubs, including Guangzhou railway hub (H1), Changsha railway hub (H2), Wuhan railway hub (H3), Zhengzhou railway hub (H4), Shijiazhuang railway hub (H5), and Beijing railway hub (H6). (29)~(31) to the lower-level model, the original non-linear model can be transformed to a linear one, which means we can directly use a standard optimization solver (e.g., CPLEX or Gurobi) to solve the lower-level model. The solution strategy of the bi-level model is as follows. Firstly, because of the upper-level model is easier to get the reasonable options owe to the few constraints it contains, we can generate a number of train operation plans based on the practical experience. Then, we substitute each plan into the lower-level program to get the objective function value of the whole bi-level programming. Finally, by comparing the objective function values among all these plans, we select the best solution.

Numerical Example
Here, we make a numerical example analysis based on the railway corridor called "Beijing-Guangzhou Railway Line" (see Figure 2), which is one of the most important north-south rail lines in China. The route, from Beijing to Guangzhou, has a total length of 2324 km. It connects six provincial capitals or municipalities and many large-and medium-sized cities, and it is one of the busiest major railways in China. The goods transported to the south on the Beijing-Guangzhou railway mainly include coal, steel, petroleum, timber and export materials. Relatively speaking, the proportion of high value-added goods to the north is higher than to the south. Therefore, this paper uses the northward direction of the Beijing-Guangzhou railway as an example for analysis. This route passes through many important railway hubs, including Guangzhou railway hub (H1), Changsha railway hub (H2), Wuhan railway hub (H3), Zhengzhou railway hub (H4), Shijiazhuang railway hub (H5), and Beijing railway hub (H6).

Parameters of the Model
The freight demand, volume, prescribed transit time and costs for each railway hub are shown in Table 3. As can be seen in the table, there is a total of 58 flows, with the flow ID from F01 to F58 corresponding to g = 1 to 58. The freight cost g R is based on the Chinese Railway Tariff (see the data published on the website of 12306) and the transportation distance of shipment g . The car volume and freight prescribed transit time are analyzed and speculated based on historical data.

Parameters of the Model
The freight demand, volume, prescribed transit time and costs for each railway hub are shown in Table 3. As can be seen in the table, there is a total of 58 flows, with the flow ID from F01 to F58 corresponding to g = 1 to 58. The freight cost R g is based on the Chinese Railway Tariff (see the data published on the website of 12306) and the transportation distance of shipment g. The car volume and freight prescribed transit time are analyzed and speculated based on historical data. Column o g in Table 3 shows the origin hub of the flow and column d g indicates the destination hub of the flow. The unit of car volume N g is car, the unit of freight transportation prescribed time T g is hour, and the unit of freight cost R g is CNY/car. The parameters of all potential train service arcs in the railway network are listed in Table 4. The columns s m and t m in Table 4 represent the origin and destination of the potential train service arc m, respectively. The column τ l m represents the running time on arc m of the train at the different speed levels, respectively. Specifically, l = 1 denotes the speed of the train at 80 km/h, l = 2 denotes the speed of the train at 120 km/h and l = 3 denotes the speed of the train at 160 km/h. The last column is the mileages for each potential train service arc. According to the railway operation engineer suggestions, ten train operation plans are generated for the upper model. We present them by assigning values to parameter α lP ijm y lP ij (see Appendix A).
In this way, there exists lots of train service arcs in each plan, e.g., the service arc A01 from hub H1 to hub H2 will be provided by train t1, t2, t4, t6 and t12. In order to facilitate solving the model, let (the origin of train service arc, the destination of train service arc, time consumption and train) denote the train service arc, e.g., (H1, H2, 10.8, t1) represents that the service arc from H1 to H2 takes 10.8 h and it is provided by train t1 in train operation plan I. The cost of express cargo train includes two parts, the first one is the fixed costs, and the second is the variable costs. Specifically, the fixed cost includes the expenditure for the train at the departure station and the running expenditure, and the variable cost is the expenditure of stops along the itinerary. At the departure railway hub, the expenditure for the train in the three speed levels is 4000 CNY, 4500 CNY and 5000 CNY, respectively. The running fee rate for the three speed levels is 1, 1.2 and 2 per kilometer. The expenditure of the train at each stop is 10% of the costs at the departure station.
When a car is transferred from one train to another train, it usually takes approximately ten hours according to the statistics of the main railway hubs in China's railway system. The reclassification delay τ k mn , m(t) = n(t), k ∈ V is influenced by capacity of the railway hub k. However, the transfer time τ k mn , m(t) = n(t), k ∈ V between two arcs belonging to the same train is significantly less than the transfer time τ k mn , m(t) = n(t), k ∈ V between two arcs belonging to different trains. Table 5 shows the values of the time delay τ k mn , k ∈ V in hub k. The unit of parameter τ k mn , k ∈ V is hour. For example, a shipment from hub H1 to hub H4 takes train service arc (H1, H2, 10.8, t1) from hub H1 to hub H2, takes train service arc (H2, H3, 7.6, t2) from hub H2 to hub H3, and takes train service arc (H3, H4, 8.3, t2) from hub H3 to hub H4 under the train operation plan I. Because of the goods flow needs to be reclassified in hub H2, it will take 11.1 h. However, the goods flow is forced to stop in hub H3 on the train t2, another 1.3 h is needed. Therefore, it will take 10.8 + 11.1 + 7.6 + 1.3 + 8.3 = 39.1 h from hub H1 to H4 for the shipment. In addition, the maximum capacity for one express cargo train is set as 50 cars, i.e., C lP m = 50.

Computational Performance of the Proposed Model
This section aims to test the computational efficiency and effectiveness of the proposed lower-level models, including the nonlinear model (NLP) and the linear model (LP). We test six different instances with different number of flows (from 10 to 58) under train operation plan I. All of the instances data come from Table 3. The parameters for each instance are listed in Table 6. The other parameters can be found in Section 5.1. Table 6. Problem instances.

Instance ID
The termination rules for all models are setup with a solver relative gap of 0 in Spyder 3.1.4.

Computational Performance of the Proposed Model
This section aims to test the computational efficiency and effectiveness of the proposed lowerlevel models, including the nonlinear model (NLP) and the linear model (LP). We test six different instances with different number of flows (from 10 to 58) under train operation plan I. All of the instances data come from Table 3. The parameters for each instance are listed in Table 6. The other parameters can be found in Section 5.1. Table 6. Problem instances.

Computational Results and Analysis
All the above instances based on the Beijing-Guangzhou railway line can be solved to optimality after running for about 10 min. The optimal objective function value and transportation demand satisfaction rate of each train operation plan are listed in Table 7. By comparing the optimal objective function values for each train operation plan, we find the plan I is the most profitable. A total of 764.03 cars were involved in the calculation of the 58 flows, and the results show that all cars could be delivered to their destinations within the prescribed transit  Table A1.
The detailed service arcs selection for each shipment is shown in Table 8. There is only one flow, F51, be reclassified that transfer from train t2 to train t6 in hub H4. The total reclassification delay is 76 car-hours (7.6 × 10). The total transfer time without reclassification is 176.963 car-hours, and most of them, approximately 98%, occur in hub H2. The operating frequency ( f lP ij ) of each train in Plan I is shown in Table 9. In this experiment, we set the value of C lP m to be 50. The actual number of cars in each train can be calculated by C lP m f lP ij . Take train t1 as an example, the actual number of cars is 44 cars (50 × 0.88) in t1. Meanwhile, the operating frequency becomes 1 in this scenario accordingly.

Conclusions
In this paper, a bi-level programming model for express cargo service network design is proposed. The upper-level model forms the optimal decisions in terms of the service characteristics, and the lower-level model carries out the flow distribution according to the train operation plan and service commitments. According to the characteristics of high value-added goods, the prescribed transit time constraint is added to the model, considering both the transportation time on the arcs and the transfer time between arcs. In addition, the arc capacity constraint and flow balance constraint are considered. Moreover, linearization techniques are used to convert the lower-level model to a linear one so that it can be directly solved by a standard optimization solver. Finally, a real-world case study is carried out based on the Beijing-Guangzhou Railway Line, which is one of the busiest railway lines in China. The example consists of six railway hubs and five segments, and we consider the northward direction.
The results indicate that the bi-level model has significance in real-life application to the rail express cargos service network design problem. Although the train operation plans come from the upper-level model were developed manually, the results of flow distribution from the lower-level model can serve as a solid aid in improving the quality of train operation planning. In the future, researchers can focus on the joint optimization of the train operation plan and the flow distribution scheme with the three-dimensional bin-packing method. It is also an interesting future work of applying multiobjective optimization methods (e.g., Pareto-optimization approaches) to solve the railway express cargo service network design problem.
Author Contributions: The authors contributed equally to this work.