Abstract
Aiming at the local freight train transit system in railway terminal, the optimization of simultaneous delivery and pickup wagon scheme on hybrid siding network is presented. The problem is formulated as a mathematical programming model. The objective function intends to minimize shunting engines’ operating cost, wagons’ travelling cost, and penalty cost for exceeding the retrieval time window. Some constraints indicating the relationship between delivery and pickup wagons, relationship between radial branches and operation batches, relationship between local wagon groups’ retrieved time and outbound train departure time, and relationship between tasks and operation batches are considered in the model. As an NP problem, using traditional method to solve the model is difficult and inefficient. A novel two-stage hybrid optimization procedure is proposed. Firstly, the three-phase approach composed of coding task sequence, dividing task batch, and generating access order sequence population is developed to generate initial solution. Secondly, an asynchronous iteration heuristic is provided. The iteration object is given as the expression form of the solution representing the delivery and pickup wagon sequence and then the initial iteration population is generated. Furthermore, the primary updating iteration population with pending and remaining groups is introduced to optimize iteration object Then, the senior updating iteration population with Lévy flight searching is also introduced to optimize further iteration object. Finally, the experimental scenarios are designed to test the proposed approach. Also, the proposed approach is compared with some other approaches and the performance of the proposed approach is evaluated by some different sized instances.
1. Introduction
Railway freight terminal is a huge and complex system which consists of marshalling yards, handling stations, and railway sidings. Railway freight terminal is the important component of the railway freight transportation chain. Therefore, it is crucial for enhancing the transportation chain effectiveness that the performance of railway freight terminal is optimized.
Delivery and pickup scheme of local wagons is a significant work for the large-scale railway freight terminals. The shunting efficiency of wagons is directly related to the wagons’ turnover rate, the goods’ delivery time, and the accomplishment of the transport production quota. The railway terminal operators must consider how to appropriately organize the delivery and pickup scheme of wagons. Throughout the whole process of wagon shunting, the local wagons will be detached sequentially from different inbound trains, delivered to relevant handling station for loading or unloading, and then retrieved back and attached onto assigned outbound trains. According to the spatial layout of handling stations in a railway terminal, we can divide these handling areas into three types: branch-shaped siding network (BSSN), radial siding network (RSN), and hybrid siding network (HSN). For a branch-shaped siding network, after a shunting locomotive has sent a wagon group to one handling station, the shunting locomotive does not need to go back to the marshalling yard before delivery of another wagon group. However, in a radial siding network, the shunting locomotive must return to the marshalling yard before running on to the next handling station. Furthermore, the hybrid siding network can be divided into some radial sub-branches, and each radial sub-branch is a branch-shaped siding network. In a hybrid siding network, after the shunting locomotive has completed all the operations belonging to the same radial sub-branch, it must return to the marshalling yard before running onto another radial sub-branch. In a radial sub-branch, the working process of shunting engine is the same as that of branch-shaped siding network. In this paper, we discuss optimization of simultaneous delivery and pickup wagon scheme on hybrid siding network and name this type of problem as SDPWS-HSN. Moreover, we develop a two-stage hybrid optimization procedure called as TSHOP.
The remainder of this paper is organized as follows. Section 2 of this paper discusses related earlier research efforts. Section 3 is devoted to the mathematical description of the SDPWS-HSN. In Section 4, we describe the mathematical formulation. Section 5 introduces a two-stage hybrid optimization procedure combining the three-phase approach and asynchronous iteration heuristic for solving SDPWS-HSN. The computational and compared experiments are described in Section 6, and the computational results show the effectiveness of the proposed method. The current work and extensions are summarized in the last section.
2. Literature Review
In this section, we review the related literature on simultaneous delivery and pickup wagons. Literature review indicates that relevant studies can be divided into three groups, that is, vehicle routing problem, vehicle assignment, and delivery and pickup problem. The focus of our literature review will primarily be on the model and approach since it is also what we have taken in this paper.
Vehicle routing problem is the main part of simultaneous delivery and pickup problem and has been widely studied during the last decades. Many variants with different characteristics have been considered and a variety of mathematical and computational techniques have been proposed to solve them. Archetti et al. [1] presented the inventory routing problems with pickups and deliveries. Bettinelli et al. [2] provided the multitrip separate delivery and pickup problem with time windows and solved it by using a branch-price-and-cut algorithm. Hosseinabadi et al. [3] proposed a combinatorial algorithm based on the gravitational emulation local search strategy to solve the open vehicle routing problem. Reil et al. [4] considered the vehicle routing problems with backhauls and time windows. Archetti et al. [5] presented the classical capacitated vehicle routing problem. Zulvia et al. [6] proposed a green vehicle routing problem with time windows and time dependency and developed a multiobjective gradient evolution algorithm. Santos et al. [7] introduced a branch-and-cut-and-price algorithm for the two-echelon capacitated vehicle routing problem. The algorithm relies on a reformulation based on q-routes that combines two important features. Avci and Topalogly [8] developed an adaptive local search solution approach for vehicle routing problem. Cheng et al. [9] established a coordinated optimized model of quantity and route for empty wagon distribution. Guo et al. [10] provided a general model for delivering and fetching wagons with mixed-shaped goods operation sites and introduced an improved simulated annealing algorithm. Liu et al. [11] proposed a two-stage mathematical model to optimize railway freight center stations’ location and wagon flow organization together. Zheng et al. [12] presented an integer programming model for emergency railway wagon scheduling and developed a hybrid biogeography-based optimization algorithm.
The dividing batch is the crucial point of the SDPWS-HSN and essentially belongs to the vehicle assignment. In recent years, many research studies have considered the vehicle assignment problem from different perspectives, using mathematical model method and heuristic algorithm. Choi et al. [13] considered the vehicle assignment problem with demand uncertainty in a hybrid hub-and-spoke network with a single hub. Spliet et al. [14] introduced the time window assignment vehicle routing problems with time-dependent travel times and developed an exact labeling algorithm and a tabu search heuristic. Györgyi et al. [15] proposed a probabilistic approach to solve the dynamic delivery and pickup problems with time window uncertainty. João et al. [16] presented an approach based on integrated flight scheduling and fleet assignment model to minimize generalized total transport network costs. Batsyn et al. [17] introduced the fleet site-dependent vehicle assignment problems with split deliveries and described two types of approaches. Wu et al. [18] presented a stated-choice survey to investigate the attractiveness to customers of two mechanisms for managing fleet volatility: virtual queuing and guaranteed advance reservation, and the results of the study are helpful for the manufacturer to classify the market. Li et al. [19] discussed a particular vehicle assignment problem arising when carrying out loading or unloading operations and proposed a bi-level programming model, and an alternating algorithm by updating predicted control parameters was given. Then, Li et al. [20] focused on the multistage heterogeneous fleet scheduling with fleet sizing decisions, established a mixed integer programming model, and developed a hybrid simulated annealing algorithm. Wang et al. [21] developed a collaborative neuro dynamic optimization approach to solve multivehicle task assignment problems. Qin et al. [22] proposed a dynamic task assignment method for vehicles based on the multiagent reinforcement learning. Zhilyakova et al. [23] presented graph methods for solving unconstrained and constrained locomotive assignment problem on a single-line railway section.
Many studies were addressed from the viewpoint of optimization models and solution methods for the delivery and pickup problem. The delivery and pickup problem has great reference significance for our research. Guo et al. [24] turned the delivery and pickup problem into a single machine scheduling problem and used Hungarian algorithm to calculate the optimal solution of the assignment problem. Adamo et al. [25] presented the problem of cargo pickup path and speed optimization with time window constraints and proposed a branch and bound algorithm. Wang et al. [30] introduced the problem of cooperative green delivery and pickup with the goal of minimizing carbon emissions. An exact solution procedure based on cooperative game theory was provided. Ghilas et al. [27] proposed a branch and bound algorithm to solve the problem of delivery and pickup of cargo with time windows. Haddad et al. [31] discussed the multivehicle one-to-one delivery and pickup problem with split loads and proposed an iterated local search metaheuristic and a branch-and-price algorithm. Danloup et al. [32] presented the problem of delivery and pickup of cargo with transshipment and proposed a large-scale neighborhood search algorithm to solve the problem. Wang et al. [26] proposed a hybrid algorithm based on integrated improved particle swarm optimization algorithm and ant colony algorithm to solve the problem of delivery and pickup optimization based on shared fleet of collaboration enterprise. Capelle et al. [31] built an integer programming model considering the location-routing problem of cargo operations and provided a column-generation solving algorithm. Zhu et al. [32] proposed a synchronous delivery and pickup problem with stochastic demand and developed a new delivery strategy which allows vehicles to cooperate in the delivery of goods. Boysen et al. [33] described the operation process of railway container yard. Li et al. [34] proposed the optimization of cargo train inspection and disassembly and assembly operations at the railway marshalling station and formulated a sorting model based on multimodulation and inspection groups. Shi et al. [35] built a multilayer network flow model to describe the train arrival, disintegration, assembly, grouping, and starting operations of railway marshalling stations and gave a heuristic solution. Jaehn et al. [36] presented the optimization of delivery and pickup wagons for railway sidings and introduced exact algorithms with polynomial runtimes. Gestrelius et al. [37] comprehensively analyzed the arrival, disintegration, assembly, marshalling, and departure process of freight trains in marshalling yard and gave an integer programming model of shunting operation in marshalling yard with the objective function of minimizing the number of wagons towed by the engine and the occupancy of shunting line.
The focus of this paper is the development of some models based on the SDPWS-HSN to aid in making decisions on the delivery and pickup scheme on hybrid siding network. Further, a two-stage hybrid optimization procedure (TSHOP) is proposed for solving the model. A three-phase approach as the first component of TSHOP is developed to generate initial solution. An asynchronous iteration heuristic as the second component of TSHOP is provided to update solution population. The conclusion of the experiment can intuitively provide help for decision makers of railway terminals in formulating development strategies.
3. Problem Formulation
The whole processes of delivery and pickup scheme of local wagons on hybrid siding network are described in this section. We explicate some foundational concepts and analyze the characteristics of the problem. Then, we define the relevant notation and formulate simultaneous delivery and pickup shunting scheme of wagons on hybrid siding network (SDPWS-HSN) as a programming model.
3.1. Understanding Delivery and Pickup Shunting of Wagons on Hybrid Siding Network
Here we discuss the layout of hybrid siding network which consists of marshalling yard and handling stations in railway terminal.
The marshalling yard consists of three kinds of working areas: the arrival yard, the departure yard, and the shunting yard which has some classification tracks with a sophisticated braking system. The arrival yard and shunting yard are often connected by a hump, and an automatic switching system sequentially directs the wagons to their assigned classiication tracks of shunting yard. In handling stations, the local wagons can be loaded or unloaded.
Large quantities of inbound freight trains arriving to the marshalling yard in a nondirect manner are stopped on the arrival yard where local wagons are decoupled and inspected. Considering the relevant constraints, these local wagons need to be rolled into the shunting yard and sorted into their assigned classiication tracks. When the delivery and pickup shunting of wagons is carried out, a shunting locomotive is sent from the engine depot and the local wagons on some given tracks are coupled. After all handling stations are visited, the taken-out wagons are pulled back to shunting yard and sorted into some new outbound trains. Finally, outbound trains must leave from the departure yard. The continuous delivery and pickup shunting begins. The process block diagram illustrating the delivery and pickup shunting of wagons on hybrid siding is shown in Figure 1.

3.2. Problem Definition and Notation
To formulate this problem, we present the complete notation for the problem here.
3.2.1. Sets
U is the set of marshalling yard and handling stations in hybrid siding network. Define U = {u∣u = 0, 1, …, }, where is the total number of handling stations. Especially, we use u = 0 to denote marshalling yard.
H is the set of radial branches in hybrid siding network. Define H = {h∣h = 1, 2, …, }, where is the total number of radial branches.
L is the set of inbound freight trains arriving at marshalling yard. Define L = {l∣l = 1, …, }, where is the total number of inbound freight trains throughout the shunting horizon. Let l denote the serial number of inbound freight train. Especially, we use l = 0 to denote wagon groups which are parked on handling station u at the beginning of shunting horizon.
S is the set of operation types. S is written as S = {s∣s = 1, 2}. We denote s = 1 as delivery wagons. Let s = 2 represent pickup wagons.
R is the set of local wagon groups. Define R = {|l ∈ L, h ∈ H, u ∈ U}, which is decoupled from inbound freight train l to handling station u in the radial branch h for loading and unloading. Especially, we define as the detention wagon groups which are parked on handling station u in the branch h at the beginning of shunting horizon.
Ω is the task set of local wagon groups. Define Ω = {∣l ∈ L, h ∈ H, u ∈ U, s ∈ S}, where is the task of delivery wagon group , is the task of pickup wagon group , and is the task of pickup detention wagon group .
W refers to the batch serial number sequences of delivery and pickup wagons. Define W = {∣ = 1, …, }, where is the total number of batches for the whole shunting horizon.
3.2.2. Parameters
is the travel time from handling station u to handling station .
is the number of local wagon group . Especially, we use to denote the number of detention wagon groups at handling station u in radial branch h.
is the loading and unloading time of the wagon group .
is the completion time point of shunting for break-up of the inbound freight train l.
is the departure time window of outbound train coupled by the local wagon group . ET () is the earliest departure time point of outbound train coupled by the local wagon group . LT () is the latest departure time point of outbound train coupled by the local wagon group .
d is the tonnage of shunting engine traction. Here it is expressed as the number of wagons pulled by the shunting engine.
Zu is the capacity of handling station u.
is the unit-minute operating cost of the shunting engine.
is the unit-minute operating cost of the local wagon group.
is the unit-minute penalty cost caused by the failure to retrieve the local wagon group in time.
3.2.3. Activity Variables
is the time point of shunting engine’s departure from marshalling yard in the th batch.
is the time point of shunting engine’s return to marshalling yard in the th batch.
is the time point of local wagon group ’s return to marshalling yard.
is the completion time point of handling operations of the local wagon group .
t (j) is the time when the local wagon group numbered j is taking or sending.
is the number of wagon dispatched to handling station u in the th batch.
is the number of wagon retrieved from handling station u in the th batch.
is the total number of wagons dispatched to all handling stations in the th batch.
is the total number of wagons retrieved from all handling stations in the th batch.
is the total number of wagons pulled by the shunting engine between handling stations u and in the th batch.
3.2.4. Decision Variables
is the batch serial number of task .
is the order number of task .
is the access order sequence of task set R1.
is the binary variable, taking the value 1 if task access order sequence including path in the th batch, and 0 otherwise.
3.3. Objective Function
The objective function aims to minimize the total cost composed of wagons’ operating costs, shunting engines’ operating cost, and penalty cost. It is given by
The first term in the objective function is the shunting engines’ operating cost. It can be obtained by the product of the time consumption of one batch operation and the unit-minute operating cost. The second term in objective function is the wagons’ operating cost. The value can be computed by the product of the wagons’ number of any path, travel time of any path, and the unit-minute operating cost. The third term in objective function is the penalty cost. It can be counted by the difference among the earliest departure time point of departure outbound train and retrieving time point of local wagon groups.
3.4. Constraints
(1)Relationship between pickup wagon and delivery wagon: any local wagon group includes delivery task and pickup task . Obviously, the order number of delivery task must be less than that of pickup task. It can be written as(2)Relationship between radial branch and operation batch: the local wagon groups going to different radial branch must not be divided into the same batch. It is given by(3)Traction number of shunting engine: the total number of wagons pulled by the shunting engine in any path does not exceed the traction number of shunting engine. We have(4)Relationship between task and operation batch: the task of delivery and pickup wagon can only be allocated to one batch number. It can be expressed as(5)Equalization of wagon amount: the sum amount of delivery wagons and detention wagons is equal to the amount of pickup wagons over the entire planning horizon. It is shown as follows:(6)Marshalling yard as the start point and destination point of task: there is only one marshalling yard in the railway freight terminal studied in this paper. So, any task starts from the marshalling yard and eventually returns to the marshalling yard. We have(7)Relationship between local wagon groups’ retrieved time and outbound train departure time: the time point of local wagon group’s return to marshalling yard must be earlier than the latest departure time point of outbound train coupled by the local wagon group. It obviously can be expressed as(8)Handling station capacity: the wagons’ amount for loading and unloading in same station and batch must not exceed the capacity of handling station. We have3.5. Reconstruction of the Model
We simplify the mathematical model by reducing the constraint on the purpose of developing solving method. Constraints (8) and (9) are transformed into the penalty functions and then are put into the objective function as one component. So, the reformulated mathematical model is given bys.t. constraints (2)–(7).
Here the parameters M1 and M2 are defined as penalty coefficients and are two large enough positive numbers. The reformulated objective function (10) is defined as measure function to evaluate the performance of the proposed approach.
4. Solution Algorithm
As an NP problem, it is difficult and inefficient to solve the SDPWS-HSN model by traditional methods. A novel method is proposed to solve this problem. This section proposes a two-stage hybrid optimization procedure, called as TSHOP, to solve the SDPWS-HSN. TSHOP includes two components: three-phase approach for initial solution and asynchronous iteration heuristic for updating solution population. Section 4.1 provides the three-phase approach used to code task sequence, divide task batch, and generate access order sequence population. Section 4.2 develops the asynchronous iteration heuristic used to make a decision on optimal routing of delivery and pickup wagons.
4.1. Initial Solution Generation with TPA
The three-phase approach composed of coding task sequence, dividing task batch, and generating access order sequence population is developed to generate initial solution. This approach is abbreviated as TPA in order to facilitate the expression of the problem. An overview of the framework of TPA is shown in Figure 2. Stage 1. Solution encoding. Step 1.1. Natural number coding sequence. According to the task set Ω, natural number is used to mark all tasks in ascending order. The delivery task is coded as an odd number and the pickup task is encoded as an even number. For the convenience of expression, let Ω be rewritten as Ω = { + 1, …, + Q |l ∈ L, h ∈ H, u ∈ U, s ∈ S}, and Q is the total number of delivery and pickup local wagon groups. Define σ (Ω) as the natural number coding sequence of task set Ω. We have σ (Ω) = {σ ( + 1), …, σ ( + Q)}. Step 1.2. Rearranged coding sequence. Generating random number is executed Q times and Q random numbers correspondingly replace Q natural numbers in σ (Ω). Then, natural number coding sequence σ (Ω) is ranked in ascending order of random numbers. So, the rearranged coding sequence β (Ω) can be given by β (Ω) = {β ( + 1), …, β ( + Q)}. Step 1.3. Updated coding sequence. According to constraint (2) in the SDPWS-HSN model, all rearranged coding sequences are rearranged so that the delivery wagon task is located in front of pickup wagon task with the same inbound freight train l, handling station u, and radial branch h. Any pair of tasks’ failure to meet the above restriction in the rearranged coding sequence shall be interchanged. So, the updated coding sequence γ (Ω) can be obtained and is expressed by γ (Ω) = {γ ( + 1), …, γ ( + Q)}. Stage 2. Batch division. Step 2.1. Dividing batch with radiation branch. According to constraint (3) in the SDPWS-HSN model, two adjacent tasks going to different radial branch cannot be allocated into the same batch. The batch breakpoint is located between two adjacent tasks going to different radiation branches. It means that these tasks before and after breakpoint are assigned to different batches. Step 2.2. Dividing batch with traction number. According to constraints (4) and (5) in the SDPWS-HSN model, the wagon amount allocated to shunting engine in a batch should be less than shunting engine traction d. The batch breakpoint needs to be inserted into the proper location of updated coding sequence. According to the order of updated coding sequence, the wagon amount of tasks is accumulated in turn until it reaches the shunting traction number. Accumulation process is stopped and batch breakpoint is put in the updated coding sequence. Stage 3. Population initialization. Step 3.1. Access order coding sequence. We, respectively, execute the procedure of dividing batch with radiation branch and traction number for updated coding sequence so that the access order coding sequence can be obtained. Step 3.2. Access order coding sequence population. Stage 1 program and stage 2 program are repeatedly executed ℑ times so that ℑ access order sequences are created. The ℑ access order coding sequences are collected to set up access order coding sequence population. The access order coding sequence population is shown in Figure 3.


4.2. Solution Population Updating with AIH
In this section, we propose an asynchronous iteration heuristic, embedding primary updating and senior updating for delivery and pickup routing, called as AIH-PU&SU. An overview of the framework of AIH is shown in Figure 4. The detailed steps are as follows. Stage 1. Setting control parameters of iteration procedure. Let the access order coding sequence population be initial iteration population. The access order coding sequence is made as the iteration object in the iteration procedure. The iteration object will be continually updated using the proposed approaches. The iteration population is composed of ℑ iteration objects. Let ℜ be the max number of iteration throughout iteration procedure. The ξth iteration population is expressed as ωξ ξ = 1, …, ℜ. The iteration object number in each iteration population is kept to ℑ for whole iteration horizon. Stage 2. Primary updating iteration population with pending and remaining group. Step 2.1. Splitting iteration object group in iteration population. Any ξth iteration population throughout whole iteration horizon is divided into two groups, i.e., pending iteration object group and remaining iteration object group . The number of iteration object in the two groups is the same. The pending iteration object group is updated firstly with proposed approach. Then, the remaining iteration object group is updated secondly with different proposed approach. Finally, two implement object groups are combined into a new iteration population. Step 2.2. Updating pending implement object group with random vector group.(1)Pending implement object group. Half of the implement objects in ξth iteration population ωξ are randomly chosen to set up pending implement object group. The pending implement object group is given by = .(2)Random number vector group. ℑ/2 random number vectors with the certain range from 0 to 1 are generated to mark all implement objects in the pending set. The random number vector group is denoted as = . The relationship between pending implement object group and random number vector group is shown in Figure 5.(3)New random number vector group. One vector ρ (γk) is randomly selected from random number vector group . The ith vector ρ (γi) in random number vector group and selected vector ρ (γk) are put into equation (11) to calculate ith new random number vector. The new random number vector can be obtained by Here we put a disturbance parameter into equation (11). The disturbance parameter is a random number between 0 and 1. It intends to achieve a lager searching range. When all vectors in random number vector group perform the above updating procedure, a new random number vector group can be obtained. The updating procedure of new random number vector group is shown in Figure 6.(4)Updating pending implement object group. All random number vectors in new random number vector group are ranked in ascending order of random number value. Each implement object matched with random number vector should be reordered accordingly in the pending implement object group. So, the new pending implement object group can be obtained. Then, each implement object is compared with that of pending implement object group. Then, the implement object with lower function value will be collected into the updating pending implement object group. The updating pending implement object group is written as = {, …, }. The hierarchical structure of this process is outlined in Figure 7. Step 2.3. Updating remaining implement object group with accumulation probability and random vector group.(1)Remaining implement object group. The remaining half of the implement objects in ξth iteration population ωξ are collected to set up remaining implement object group. The remaining implement object group is given by = .(2)Updating remaining implement object group with accumulation probability. For a minimization problem, the solution with lower objective function value has better performance. The cumulative probability of updating pending implement object group is given by The remaining implement object is taken out in turn from remaining group . Meanwhile a random number uniformly distributed between 0 and 1 is generated. If there is , let = . Besides, if there is , let = . The first updating remaining implement object group is written as = {, …, }.(3)Updating remaining implement object group with random vector group. The approach used in stage 2 is adopted to generate the updating remaining implement object group. The updating remaining implement object group is written as = {, …, }. Step 2.4. Constructing primary updating iteration population. The updating pending implement object group and the updating remaining implement object group are merged to set up new iteration population. The primary updating iteration population is written as follows: Stage 3. Senior updating iteration population with Lévy flight searching. Step 3.1. Generating random number vector group. ℑ random number vectors with the certain range from 0 to 1 are generated to mark all implement objects in the primary updating iteration population . The random number vector group is denoted as = . The implement object with the lowest function value is searched from the primary updating iteration population . It is called as the local implement object and denoted as . The random number vector matched with the local optimal implement object can be called as local optimal random number vector and written as . Step 3.2. Updating random number vector group with Lévy flight searching. Lévy flight searching is used to update further the random vector group. The short searching step and long searching step are alternately adopted to avoid falling into the local optimal solution and explore the wider range of solution space in the approach. The ith vector in the random number vector group and local optimal vector are put into Lévy flight searching equation (14) to calculate ith new random number vector. The new random number vector can be obtained by The parameter is called as searching step factor. It is a random number with the range from 0 to 1 and intends to achieve a lager searching range. The parameter Lévy is a random variable with Lévy distribution and can be obtained by The random variable obeys Gaussian distribution with mean 0 and variance 1. The random variable follows Gaussian distribution with mean 0. The variance of can be written as , where is the gamma function. When all vectors in random number vector group perform the above updating procedure, a new random number vector group can be obtained. The updating procedure of new random number vector group is shown in Figure 8. Step 3.3. Achieving senior updating iteration population. All random number vectors in new random number vector group are ranked in ascending order of random number value. Each implement object matched with random number vector should be reordered accordingly in the primary updating iteration population . So, the new iteration population can be obtained. Then, each implement object is compared with that of primary updating iteration population. Then, the implement object with lower function value will be collected into the senior updating iteration population. The senior updating iteration population is written as = {, …, }. The hierarchical structure of this process is outlined in Figure 9. Stage 4. Finding global optimal implement object. Step 4.1. Termination criterion. For the whole iteration horizon, repeat stage 2 to stage 3 until the iteration number arrives at R. Step 4.2. Global optimal implement object. All implement objects are compared by their objective function value. The best of them is identified as the global optimal implement object. So, the optimal handling station accessing sequence can be obtained. Step 4.3. Termination criterion. Iterative optimization is terminated until the maximum number of iterations maxgen is reached. Output the objective function value and the change curve with the number of iterations and record the global optimal implement object.






5. Case Study
In this section, we try to evaluate the performance of the TSHOP provided in previous section to solve the SDPWS-HSN by some traditional measures such as objective function and execution time. The scenarios for SDPWS-HSN are described in Section 5.1. The key control parameters of TSHOP are set in Section 5.2. In Section 5.3, many small-size and large-size contrast experiments are developed to test the performance of the TSHOP compared with CPLEX solver and some proposed swarm intelligence algorithms. Section 5.4 reports the relationship between the iteration number and the quality of the solution in three sets of experiments.
5.1. Scenario Settings of the Experiments
In order to design the experiments scheme, the scenarios settings must be firstly set up. We choose Zhengzhou railway terminal in China as the tested target object and generate the tested dataset. However, the collected data still have some inconsistencies caused by the total number of local wagons, the time interval of the inbound trains, wagon number of each wagon group being dispatched to various handling stations, loading or unloading time of wagon group being dispatched to various handling stations, and so on. Therefore, the data have to be cleaned.
A hybrid siding network consisting of 11 handling stations is shown in Figure 10. The travel time f between a pair of handling stations is listed in Table 1 where 0 represents the marshalling station. The handling station capacity is affected by the effective length of the handling siding. The handling station capacity is listed in Table 2. During the whole shunting horizon, 4 freight trains arrive at the marshalling station one by one. Data of arrival train and wagon flow in railway terminal are displayed in Table 3.

5.2. Parameter Setting
The test work has been executed using a computer with Intel(R) Core(TM) i7-10510U CPU @ 1.80 GHz 2.30 GHz. The hybrid heuristic combining asynchronous iteration procedure TSHOP has been coded in Microsoft Visual C++.
Here we give unit-minute operating cost of the shunting engine η1 = 16, unit-minute operating cost of the local wagon group η2 = 1.4, unit-minute penalty cost caused by the failure to retrieve the local wagon group in time η3 = 9, and the tonnage of shunting engine traction d = 40.
Three key parameters are closely related to the performance of TSHOP. The three parameters are, respectively, the iteration population size ℑ, the maximum number of iterations ℜ, and disturbance parameter φ. Here they are tested to achieve a better combination configuration. The testing scope of three parameters is listed in Table 4. We collected 27 sets of parameters from Table 4 and adopted the experiment scenario mentioned in Section 5.1 to test them. Each set is tested 10 times. Also, we use two measures of performance. The first one is the average value of the objective functions obtained by executing TSHOP approach 10 times. In order to facilitate representing the measure, it is abbreviated to OFV. The second measure of performance is the average CPU time to execute the TSHOP approach 10 times. The testing results of parameter combination are illustrated in Table 5.
In Table 5, all instances are executed with no CPU time limit imposed. The performance of TSHOP approach improves with ℑ and ℜ increasing. However, ℜ has more influence on the average value of the objective functions and ℑ has an obvious impact on the CPU time. Moreover, the parameter φ needs to arrive at a balance value. The too high or too low value of φ affects the performance of solution. From the testing results, we find the best parameter set as ℑ = 50, R = 160, and φ = 0.6.
5.3. Performance Evaluation
5.3.1. Performance for Small-Size Instance
To evaluate the performance of the proposed heuristic, we firstly give four sets of contrast experiments with 8 item tasks, 10 item tasks, 12 item tasks, and 14 item tasks using the same hybrid siding network composed of 11 handling stations shown in Figure 10. They all are small-size instances. Here the CPLEX solver can be used to solve them. The four sets of contrast experiments are oriented to compare the performances of the TSHOP approach with that of CPLEX solver.
The computer execution time of TSHOP approach and CPLEX solver is strictly limited within 3000 seconds. Each contrast experiment is executed 10 times. We use two measures to evaluate the performance. They are, respectively, the average objective function value and average CPU time. Then, in order to facilitate the analysis of the problem, we provide the deviation rate as one analytical indicator. The deviation rate is abbreviated as RAT. It represents the deviation rate from TSHOP to CPLEX on optimal solution and can be calculated by RAT = [OFV (TSHOP) − OFV (CPLEX)]/(OFV (CPLEX)) × 100%.
The results of performance comparison are displayed in Table 6.
For the first three sets of contrast experiments with 8 item tasks, 10 item tasks, and 12 item tasks, the objective function value of CPLEX solver is slightly better than that of TSHOP. But the CPU time of CPLEX solver is obviously higher than that of TSHOP. When the contrast experiments arrive at 14 item tasks, the CPU time of CPLEX solver goes further away from the reasonable interval. The feasible solution cannot be found by CPLEX solver within 3000 seconds and is denoted as the symbol “#.” Thus, one obvious conclusion is that the performance of TSHOP is significantly prominent when compared to CPLEX solver.
5.3.2. Performance for Large-Size Instance
When the increase of task number goes a step further, the CPU time of CPLEX solver will go further away from the reasonable interval. The swarm intelligence algorithm has outstanding performance advantages in solving large-size instance. In many swarm intelligence algorithms, the genetic algorithm (GA) and artificial bee colony algorithm (ABCA) have the characteristics of few parameters, strong search ability, and easy implementation. GA and ABCA can only be adopted to implement the iterative procedure in order to find optimal solution. Also, TPA procedure in TSHOP needs to be introduced to generate initial solution. The TPA procedure is drawn into GA and ABCA to generate two new approaches written as TPA-GA and TPA-ABCA. The two approaches will be compared with TSHOP approach proposed in this paper.
In this section, the three sets of large-size contrast experiments are oriented to compare the performances of the TSHOP approach with the TPA-GA approach and TPA-ABCA approach. In the first set of contrast experiment, there are 16 item tasks for delivery and pickup local wagon group. In the second set of contrast experiment, there are 24 tasks. In the third set of contrast experiment, the task amount increases to 34.
We use similarly two measures to evaluate the performance. The first measure of performance is the average objective function value obtained by the TSHOP approach, the TPA-GA approach, and TPA-ABCA approach. The second is the CPU time to run TSHOP program, TPA-GA program, and TPA-ABC program.
In order to facilitate the analysis of the problem, we present two analytical indicators. The first indicator denoted as RAT1 represents the deviation rate from the TSHOP to TPA-GA on the optimal solution. Let OFV (TPA-GA) be the average value of objective functions obtained by running TPA-GA approach 10 times. Let OFV (TSHOP) be the average value of objective functions obtained by running TSHOP approach 10 times. RAT1 is calculated by RAT1 = [OFV (TPA-GA) − OFV (TSHOP)]/OFV (TSHOP) × 100%.
The second indicator denoted as RAT2 represents the deviation rate from the TSHOP to TPA-ABC on the optimal solution. Let OFV (TPA-ABC) be the average value of objective functions obtained by running TSHOP program 10 times. RAT2 is calculated by RAT2 = [OFV (TPA-ABCA) − OFV (TSHOP)]/OFV (TSHOP) × 100%.
Here we chose control parameter combination from Table 4 for testing three compared algorithms. Each instance is solved 10 times using the parameter listed in Table 4, and the average results are recorded in Tables 7–9.
From above three sets of testing instances, we can conclude as follows:(1)For the three sets of different size instance, the performance of TSHOP is significantly prominent when compared to TPA-GA and TPA-ABCA. Although the average CPU time of TSHOP is more than that of TPA-ABCA and TPA-GA, the performance of TSHOP is better than others when the major criterion OFV is adopted. The average deviation rate from the TSHOP to TPA-ABCA and TPA-GA on the optimal solution, respectively, reaches 16% and 22%. The performance comparisons are shown in Table 10. In sum, we can conclude that the TSHOP is effective to obtain the good delivery and pickup wagon scheme for SDPWS-HSN problem.(2)At the instance with 16 item tasks stage, it can be observed that the TSHOP approach exhibits slight superiority on the major criterion OFV and analytical indicators RAT. When the increase of task amount goes a step further into 24 item tasks stage, it is an interesting discovery that the superiority shows a steady upward trend. After entering the stage of 34 item tasks, an obvious observation is that the superiority begins to appear a accelerated upward trend. Meanwhile, the CPU time for obtaining solutions by TSHOP approach is still maintained in a reasonable range.
Now for the largest-size instance with 34 item tasks, we find the optimal implement objects, i.e., the optimal handling station accessing sequences. The TSHOP, TPA-ABCA, and TPA-GA are all executed by the best parameter set as ℑ = 50, R = 160, and φ = 0.6. The result is listed in Table 11.
We focus on the performance of TPA-GA, TPA-ABCA, and TSHOP program under 160 iterations. Here the objective function value of global optimal implement objects obtained by the TPA-GA program executed with 160 iterations is 64074.6. The objective function value of global optimal implement objects obtained by the TPA-ABCA program executed with 160 iterations is 58081.8. The objective function value of global optimal implement objects obtained by the TSHOP program executed with 160 iterations is 46528.6. Then, the objective function values obtained by three approaches for solving the SDPWS-HSN problem are compared. The objective function values obtained by TSHOP program are obviously better than those solved by TPA-GA program and TPA-ABCA program.
5.4. Correlation between Solution Performance and Iteration Number
Compared with iteration population size ℑ and disturbance parameter φ, the iteration number ℜ has more influence on the solution performance. To further study the relationship between solution performance and iteration number, TSHOP program, TPA-GA program, and TPA-ABCA program are all set up with the same default parameter setting (ℑ = 50, φ = 0.6) and tested by the three sets of different size instances. The performance is shown in Figure 11. It can be seen that the superiority of TSHOP appears an accelerated upward trend along with the increasing of iteration times.

(a)

(b)

(c)
6. Conclusions
This study aims to develop a model and approach to aid in making decisions on delivery and pickup scheme of wagons on the hybrid siding network of the railway terminal. We name the problem as simultaneous delivery and pickup wagon scheme on hybrid siding network in railway terminal (SDPWS-HSN). In this paper, a mathematical programming model has been discussed for SDPWS-HSN. As an NP problem, it is difficult and inefficient to solve SDPWS-HSN model by using traditional methods. A novel heuristic algorithm is proposed to address this problem.
In this paper, we propose a two-stage hybrid optimization procedure, abbreviated as TSHOP, to solve the SDPWS-HSN. Firstly, a three-phase approach composed of coding task sequence, dividing task batch, and generating access order sequence population is developed to generate initial solution. Secondly, the asynchronous iteration heuristic is provided to update the solution population. Finally, we evaluate the performance of the TSHOP in terms of two traditional measures: objective function and execution time. The results show that the proposed heuristic algorithm is able to obtain better solutions with reasonable CPU time.
Future research will focus on some new heuristics to solve the SDPWS-HSN. Another prospective research will be the extension of the SDPWS-HSN model by considering some further realistic aspects and relaxing some hard constraints.
Data Availability
All data, models, and codes generated or used during the study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by the Philosophy and Social Science Applied Research Major Project of Henan Province Education Department of China (grant no. 2022-YYZD-24), Henan Province Philosophy and Social Science Planning Project (grant no. 2021BJJ087), Henan Province Science and Technology Research Plan Project of China (grant no. 202102310310), and National Natural Science Foundation of China (grant nos. U1604150 and U1804151).