# Optimization of Simultaneous Delivery and Pickup Wagon Scheme on Hybrid Siding Network of Railway Terminal

Authors: Bing Li; Ce Yun; Guanglan Xu; Xiaoqian Zhou; Hua Xuan
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## 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.


---

## Body

## 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. Section2 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 Figure1.Figure 1
Illustration of the process for delivery and pickup wagons on hybrid siding network.


### 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, …, u¯}, where u¯ 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, …, l¯}, where l¯ 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 = {lhu|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 0hu 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 Ω = {lhus∣l ∈ L, h ∈ H, u ∈ U, s ∈ S}, where lhu1 is the task of delivery wagon group lhu, lhu2 is the task of pickup wagon group lhu, and 0hu2 is the task of pickup detention wagon group 0hu.W refers to the batch serial number sequences of delivery and pickup wagons. Define W = {w∣w = 1, …, w¯}, where w¯ is the total number of batches for the whole shunting horizon.


#### 3.2.2. Parameters

πu,u′ is the travel time from handling station u to handling station u′.Zlhu is the number of local wagon group lhu. Especially, we use Z0hu to denote the number of detention wagon groups at handling station u in radial branch h.π˜lhu is the loading and unloading time of the wagon group lhu.πl0 is the completion time point of shunting for break-up of the inbound freight train l.ETlhu,LTlhu is the departure time window of outbound train coupled by the local wagon group lhu. ET (lhu) is the earliest departure time point of outbound train coupled by the local wagon group lhu. LT (lhu) is the latest departure time point of outbound train coupled by the local wagon group lhu.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.η1 is the unit-minute operating cost of the shunting engine.η2 is the unit-minute operating cost of the local wagon group.η3 is the unit-minute penalty cost caused by the failure to retrieve the local wagon group lhu in time.


#### 3.2.3. Activity Variables

πw is the time point of shunting engine’s departure from marshalling yard in the wth batch.π˜w is the time point of shunting engine’s return to marshalling yard in the wth batch.π¯lhu is the time point of local wagon group lhu’s return to marshalling yard.π⌢lhu is the completion time point of handling operations of the local wagon group lhu.t (j) is the time when the local wagon group numbered j is taking or sending.χwu is the number of wagon dispatched to handling station u in the wth batch.χ˜wu is the number of wagon retrieved from handling station u in the wth batch.χw is the total number of wagons dispatched to all handling stations in the wth batch.χ˜w is the total number of wagons retrieved from all handling stations in the wth batch.χwu,u′ is the total number of wagons pulled by the shunting engine between handling stations u and u′ in the wth batch.


#### 3.2.4. Decision Variables

wlhus is the batch serial number of task lhus.λlhus is the order number of task lhus.ψlhus is the access order sequence of task set R1.yuu′w is the binary variable, taking the value 1 if task access order sequence including path u,u′ in the wth 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(1)MinF=∑w=1w¯π˜w−πwη1+∑u=1u¯∑u′=1u¯∑w=1w¯πuu′yuu′wχwuu′η2+∑l=0l¯∑h=1h¯∑u=1u¯maxπ¯lhu−ETlhu,0η3.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 grouplhu includes delivery task lhu1 and pickup task lhu2. Obviously, the order number of delivery task must be less than that of pickup task. It can be written as(2)λlhu1<λlhu2,∀lhu∈R.(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)wlhus≠wl′h′u′s′,∀lhus,l′h′u′s′∈R,h≠h′.(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)χwu.u′≤d,∀w∈W,u.u′∈U.(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)∑w=1w¯wlhus=1,∀lhus∈R.(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)∑w=1w¯χw+∑h=1h¯∑u=1u¯ZOhu=∑w=1w¯χ¯w.(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)youw=yuow=1,∀u∈U,w∈W.(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)π¯lhu≤LTlhu,∀lhu∈R.(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 have(9)Z0hu+∑w=1w¯χwu−χ˜wu≤Zu,∀u∈U.


### 3.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 by(10)MinFnew=F+M1∑l=1l¯∑h=1h¯∑u=1u¯maxπ¯lhu−LTlhu,0+M2∑h=1h¯∑u=1u¯maxZ0hu+∑w=1w¯χwu−χ˜wu−Zu,0,s.t. constraints (2)–(7).Here the parametersM1 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.


## 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 Figure1.Figure 1
Illustration of the process for delivery and pickup wagons on hybrid siding network.


## 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, …, u¯}, where u¯ 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, …, l¯}, where l¯ 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 = {lhu|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 0hu 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 Ω = {lhus∣l ∈ L, h ∈ H, u ∈ U, s ∈ S}, where lhu1 is the task of delivery wagon group lhu, lhu2 is the task of pickup wagon group lhu, and 0hu2 is the task of pickup detention wagon group 0hu.W refers to the batch serial number sequences of delivery and pickup wagons. Define W = {w∣w = 1, …, w¯}, where w¯ is the total number of batches for the whole shunting horizon.


### 3.2.2. Parameters

πu,u′ is the travel time from handling station u to handling station u′.Zlhu is the number of local wagon group lhu. Especially, we use Z0hu to denote the number of detention wagon groups at handling station u in radial branch h.π˜lhu is the loading and unloading time of the wagon group lhu.πl0 is the completion time point of shunting for break-up of the inbound freight train l.ETlhu,LTlhu is the departure time window of outbound train coupled by the local wagon group lhu. ET (lhu) is the earliest departure time point of outbound train coupled by the local wagon group lhu. LT (lhu) is the latest departure time point of outbound train coupled by the local wagon group lhu.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.η1 is the unit-minute operating cost of the shunting engine.η2 is the unit-minute operating cost of the local wagon group.η3 is the unit-minute penalty cost caused by the failure to retrieve the local wagon group lhu in time.


### 3.2.3. Activity Variables

πw is the time point of shunting engine’s departure from marshalling yard in the wth batch.π˜w is the time point of shunting engine’s return to marshalling yard in the wth batch.π¯lhu is the time point of local wagon group lhu’s return to marshalling yard.π⌢lhu is the completion time point of handling operations of the local wagon group lhu.t (j) is the time when the local wagon group numbered j is taking or sending.χwu is the number of wagon dispatched to handling station u in the wth batch.χ˜wu is the number of wagon retrieved from handling station u in the wth batch.χw is the total number of wagons dispatched to all handling stations in the wth batch.χ˜w is the total number of wagons retrieved from all handling stations in the wth batch.χwu,u′ is the total number of wagons pulled by the shunting engine between handling stations u and u′ in the wth batch.


### 3.2.4. Decision Variables

wlhus is the batch serial number of task lhus.λlhus is the order number of task lhus.ψlhus is the access order sequence of task set R1.yuu′w is the binary variable, taking the value 1 if task access order sequence including path u,u′ in the wth batch, and 0 otherwise.


## 3.2.1. Sets

U is the set of marshalling yard and handling stations in hybrid siding network. Define U = {u∣u = 0, 1, …, u¯}, where u¯ 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, …, l¯}, where l¯ 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 = {lhu|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 0hu 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 Ω = {lhus∣l ∈ L, h ∈ H, u ∈ U, s ∈ S}, where lhu1 is the task of delivery wagon group lhu, lhu2 is the task of pickup wagon group lhu, and 0hu2 is the task of pickup detention wagon group 0hu.W refers to the batch serial number sequences of delivery and pickup wagons. Define W = {w∣w = 1, …, w¯}, where w¯ is the total number of batches for the whole shunting horizon.


## 3.2.2. Parameters

πu,u′ is the travel time from handling station u to handling station u′.Zlhu is the number of local wagon group lhu. Especially, we use Z0hu to denote the number of detention wagon groups at handling station u in radial branch h.π˜lhu is the loading and unloading time of the wagon group lhu.πl0 is the completion time point of shunting for break-up of the inbound freight train l.ETlhu,LTlhu is the departure time window of outbound train coupled by the local wagon group lhu. ET (lhu) is the earliest departure time point of outbound train coupled by the local wagon group lhu. LT (lhu) is the latest departure time point of outbound train coupled by the local wagon group lhu.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.η1 is the unit-minute operating cost of the shunting engine.η2 is the unit-minute operating cost of the local wagon group.η3 is the unit-minute penalty cost caused by the failure to retrieve the local wagon group lhu in time.


## 3.2.3. Activity Variables

πw is the time point of shunting engine’s departure from marshalling yard in the wth batch.π˜w is the time point of shunting engine’s return to marshalling yard in the wth batch.π¯lhu is the time point of local wagon group lhu’s return to marshalling yard.π⌢lhu is the completion time point of handling operations of the local wagon group lhu.t (j) is the time when the local wagon group numbered j is taking or sending.χwu is the number of wagon dispatched to handling station u in the wth batch.χ˜wu is the number of wagon retrieved from handling station u in the wth batch.χw is the total number of wagons dispatched to all handling stations in the wth batch.χ˜w is the total number of wagons retrieved from all handling stations in the wth batch.χwu,u′ is the total number of wagons pulled by the shunting engine between handling stations u and u′ in the wth batch.


## 3.2.4. Decision Variables

wlhus is the batch serial number of task lhus.λlhus is the order number of task lhus.ψlhus is the access order sequence of task set R1.yuu′w is the binary variable, taking the value 1 if task access order sequence including path u,u′ in the wth 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(1)MinF=∑w=1w¯π˜w−πwη1+∑u=1u¯∑u′=1u¯∑w=1w¯πuu′yuu′wχwuu′η2+∑l=0l¯∑h=1h¯∑u=1u¯maxπ¯lhu−ETlhu,0η3.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 grouplhu includes delivery task lhu1 and pickup task lhu2. Obviously, the order number of delivery task must be less than that of pickup task. It can be written as(2)λlhu1<λlhu2,∀lhu∈R.(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)wlhus≠wl′h′u′s′,∀lhus,l′h′u′s′∈R,h≠h′.(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)χwu.u′≤d,∀w∈W,u.u′∈U.(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)∑w=1w¯wlhus=1,∀lhus∈R.(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)∑w=1w¯χw+∑h=1h¯∑u=1u¯ZOhu=∑w=1w¯χ¯w.(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)youw=yuow=1,∀u∈U,w∈W.(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)π¯lhu≤LTlhu,∀lhu∈R.(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 have(9)Z0hu+∑w=1w¯χwu−χ˜wu≤Zu,∀u∈U.


## 3.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 by(10)MinFnew=F+M1∑l=1l¯∑h=1h¯∑u=1u¯maxπ¯lhu−LTlhu,0+M2∑h=1h¯∑u=1u¯maxZ0hu+∑w=1w¯χwu−χ˜wu−Zu,0,s.t. constraints (2)–(7).Here the parametersM1 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. Section4.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 Figure2.Stage 1. Solution encoding.Step 1.1. Natural number coding sequence.According to the task set Ω, natural number is used to mark all taskslhus in ascending order. The delivery task lhu1 is coded as an odd number and the pickup task lhu2 is encoded as an even number. For the convenience of expression, let Ω be rewritten as Ω = {lhus + 1, …, lhus + 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 σ (Ω) = {σ (lhus + 1), …, σ (lhus + Q)}.Step 1.2. Rearranged coding sequence.Generating random number is executedQ 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 β (Ω) = {β (lhus + 1), …, β (lhus + 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 lhu1 is located in front of pickup wagon task lhu2 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 γ (Ω) = {γ (lhus + 1), …, γ (lhus + 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 Figure3.Figure 2
Illustration of the three-phase approach for initial solution.Figure 3
A diagram of the access order coding sequence population.


### 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 Figure4. 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 ωξremaining and remaining iteration object group ωξremaining. The number of iteration object in the two groups is the same. The pending iteration object group ωξremaining is updated firstly with proposed approach. Then, the remaining iteration object group ωξremaining 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 ωξremaining = γ1Ω,…γℑ/2Ω.(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ρξpending = ργ1,…ργℑ/2. 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 ρξpending. 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(11)ρijnew=ρij+ϕρij−ρkj.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ρξpending perform the above updating procedure, a new random number vector group ρξnew 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ρξnew 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 ω¯ξpending = {γ1updatingΩ, …, γℑ/2updatingΩ}. 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 ωξremaining = γℑ/2+1Ω…γℑΩ.(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(12)pγjupdating=∑i=1j1/γiupdating∑i=1ℑ/21/γiupdating,j=1,…,ℑ2.The remaining implement objectγℑ/2+iΩ is taken out in turn from remaining group ωξremaining. Meanwhile a random number ϕ uniformly distributed between 0 and 1 is generated. If there is ϕ≤qγ1updating, let γℑ/2+inewΩ = γ1updatingΩ. Besides, if there is qγjupdating<ϕ≤qγj+1updating, let γℑ/2+inewΩ = γj+1updatingΩ. The first updating remaining implement object group is written as ω¯ξremaining = {γℑ/2+1newΩ, …, γℑnewΩ}.(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ω¯ξremaining = {γℑ/2+1updatingΩ, …, γℑupdatingΩ}.Step 2.4. Constructing primary updating iteration population.The updating pending implement object groupω¯ξpending and the updating remaining implement object group ω¯ξremaining are merged to set up new iteration population. The primary updating iteration population is written as follows:(13)ω¯ξprimary=γ1primaryΩ…γℑprimaryΩ.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ω¯ξprimary. The random number vector group is denoted as ρξprimary = ργ1primary,…,ργℑprimary. The implement object with the lowest function value is searched from the primary updating iteration population ω¯ξprimary. It is called as the local implement object and denoted as γoptimalprimary. The random number vector matched with the local optimal implement object can be called as local optimal random number vector and written as ργoptimalprimary.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. Theith vector ργiprimary in the random number vector group ρξprimary and local optimal vector ργoptimalprimary 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(14)ρijnew=ρijprimary+ϕ⋅Lévy⋅ρijprimary−ρoptimal,jprimary.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(15)Levy=ℵ℘1/ε.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 Varℵ=Γ1+ε×sinπε/2/Γ1+ε/2ε×2ε+1/21/ε, where Γ is the gamma function.When all vectors in random number vector groupρξprimary perform the above updating procedure, a new random number vector group ρξnew 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ρξnew 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 ω¯ξupdating. 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 ω¯ξsenior = {γ1seniorΩ, …, γℑseniorΩ}. 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 atR.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.Figure 4
Illustration of asynchronous iteration heuristic.Figure 5
Pending implement object group and random number vector group.Figure 6
New random number vector group.Figure 7
Updating pending implement object group.Figure 8
New random number vector group.Figure 9
Senior updating iteration population.


## 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 Figure2.Stage 1. Solution encoding.Step 1.1. Natural number coding sequence.According to the task set Ω, natural number is used to mark all taskslhus in ascending order. The delivery task lhu1 is coded as an odd number and the pickup task lhu2 is encoded as an even number. For the convenience of expression, let Ω be rewritten as Ω = {lhus + 1, …, lhus + 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 σ (Ω) = {σ (lhus + 1), …, σ (lhus + Q)}.Step 1.2. Rearranged coding sequence.Generating random number is executedQ 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 β (Ω) = {β (lhus + 1), …, β (lhus + 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 lhu1 is located in front of pickup wagon task lhu2 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 γ (Ω) = {γ (lhus + 1), …, γ (lhus + 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 Figure3.Figure 2
Illustration of the three-phase approach for initial solution.Figure 3
A diagram of the access order coding sequence population.


## 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 Figure4. 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 ωξremaining and remaining iteration object group ωξremaining. The number of iteration object in the two groups is the same. The pending iteration object group ωξremaining is updated firstly with proposed approach. Then, the remaining iteration object group ωξremaining 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 ωξremaining = γ1Ω,…γℑ/2Ω.(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ρξpending = ργ1,…ργℑ/2. 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 ρξpending. 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(11)ρijnew=ρij+ϕρij−ρkj.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ρξpending perform the above updating procedure, a new random number vector group ρξnew 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ρξnew 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 ω¯ξpending = {γ1updatingΩ, …, γℑ/2updatingΩ}. 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 ωξremaining = γℑ/2+1Ω…γℑΩ.(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(12)pγjupdating=∑i=1j1/γiupdating∑i=1ℑ/21/γiupdating,j=1,…,ℑ2.The remaining implement objectγℑ/2+iΩ is taken out in turn from remaining group ωξremaining. Meanwhile a random number ϕ uniformly distributed between 0 and 1 is generated. If there is ϕ≤qγ1updating, let γℑ/2+inewΩ = γ1updatingΩ. Besides, if there is qγjupdating<ϕ≤qγj+1updating, let γℑ/2+inewΩ = γj+1updatingΩ. The first updating remaining implement object group is written as ω¯ξremaining = {γℑ/2+1newΩ, …, γℑnewΩ}.(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ω¯ξremaining = {γℑ/2+1updatingΩ, …, γℑupdatingΩ}.Step 2.4. Constructing primary updating iteration population.The updating pending implement object groupω¯ξpending and the updating remaining implement object group ω¯ξremaining are merged to set up new iteration population. The primary updating iteration population is written as follows:(13)ω¯ξprimary=γ1primaryΩ…γℑprimaryΩ.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ω¯ξprimary. The random number vector group is denoted as ρξprimary = ργ1primary,…,ργℑprimary. The implement object with the lowest function value is searched from the primary updating iteration population ω¯ξprimary. It is called as the local implement object and denoted as γoptimalprimary. The random number vector matched with the local optimal implement object can be called as local optimal random number vector and written as ργoptimalprimary.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. Theith vector ργiprimary in the random number vector group ρξprimary and local optimal vector ργoptimalprimary 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(14)ρijnew=ρijprimary+ϕ⋅Lévy⋅ρijprimary−ρoptimal,jprimary.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(15)Levy=ℵ℘1/ε.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 Varℵ=Γ1+ε×sinπε/2/Γ1+ε/2ε×2ε+1/21/ε, where Γ is the gamma function.When all vectors in random number vector groupρξprimary perform the above updating procedure, a new random number vector group ρξnew 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ρξnew 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 ω¯ξupdating. 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 ω¯ξsenior = {γ1seniorΩ, …, γℑseniorΩ}. 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 atR.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.Figure 4
Illustration of asynchronous iteration heuristic.Figure 5
Pending implement object group and random number vector group.Figure 6
New random number vector group.Figure 7
Updating pending implement object group.Figure 8
New random number vector group.Figure 9
Senior updating iteration population.


## 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 Section5.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 Figure10. 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.Figure 10
Hybrid siding network in railway terminal.Table 1
Travel time of shunting engine between a pair of handling stations in railway terminal (unit: min).
u0123456789101100181920191920262527292211802027272538444345474021920024282739454446484132027240262240464547494241927282602739454446484151925272227039454446484162038394039390282826494272644454645452802624554882543444544442826023544792745464746462624230564910294748494848495554560511122404142414142484749510Table 2
Handling station capacity (unit: wagons per day).
u1234567891011Zu1836222425282622212132Table 3
Data of arrival train and wagon flow in railway terminal.
Natural number coding of taskSerial number of trainCompletion time point of train break-upHandling stationsNumber of local wagon groupOperation typesHandling timeCompletion time point of handling operationsDeparture time window of outbound train coupled by the local wagon group1106:00111155—[9:20, 10:00]2106:00111255—[9:20, 10:00]3106:00312160——4106:00312260——5106:00615175—[12:00, 13:00]6106:00615275—[12:00, 13:00]7106:00811155——8106:00811255——9206:2628140——10206:2628240——11206:26410150——12206:26410250——13206:26712160——14206:26712260——15206:26915175—[11:40, 12:40]16206:26915275—[11:40, 12:40]17307:3059145——18307:3059245——19307:3067135——20307:3067235——21307:301013165—[11:00, 12:10]22307:301013265—[11:00, 12:10]23408:1035130——24408:1035230——25408:10612160——26408:10612260——27408:101015175——28408:101015275——29408:10119145—[11:30, 12:00]30408:10119245—[11:30, 12:00]310—3723506:30—320—61427007:00—330—91527507:00[10:20, 11:00]340—11824006:20—


### 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.Table 4
Testing scope of three parameters.
ParameterRangeIteration population size ℑ{50, 80, 100}Iteration maximum number ℜ{100, 120, 160}Disturbance parameterφ{0.3, 0.6, 0.9}Table 5
Testing results of parameter combination.
No.ℑℜφOFVCPU time1501000.33909.7360.752501200.33897.4963.183501600.33822.6161.334501000.63847.3155.735501200.63870.0758.626501600.63655.8660.467501000.93927.0762.538501200.93909.7360.739501600.93897.4961.4110801000.33954.3784.6511801200.33943.3481.2312801600.33898.6584.3513801000.63899.1783.6614801200.63903.5683.5715801600.63843.7780.6216801000.93878.0482.5417801200.93922.4181.7118801600.93854.8785.43191001000.33785.26108.69201001200.33764.51108.85211001600.33671.27109.12221001000.63783.75110.31231001200.63742.61108.09241001600.63753.71107.11251001000.93821.01109.66261001200.93782.56109.73271001600.93756.90110.36In Table5, 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 Figure10. 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 byRAT = [OFV (TSHOP) − OFV (CPLEX)]/(OFV (CPLEX)) × 100%.The results of performance comparison are displayed in Table6.Table 6
Performance for CPLEX and TSHOP applied to small-size problem.
Tasks sizeCPLEXTSHOPRAT (%)OFVCPU timeOFVCPU time83655.86569.733655.8660.460.00106676.331493.726917.81109.373.621210127.373315.4510796.42153.286.6114#3600.013617.36197.33#Average6819.852515.7810393.14146.545.24For 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 byRAT1 = [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 byRAT2 = [OFV (TPA-ABCA) − OFV (TSHOP)]/OFV (TSHOP) × 100%.Here we chose control parameter combination from Table4 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.Table 7
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 16 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.614246.9416077.1815268.66240.07106.13199.920.130.07501600.614047.1616380.3215366.64350.05167.62311.830.170.091001200.616978.2419411.9218526.36212.2696.50165.300.140.091001600.616601.6819522.218267.66331.52150.28272.930.180.10Average15468.5017847.916857.33283.47130.13237.500.150.09Table 8
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 24 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.623829.4027264.2027191.00347.47153.76308.450.140.14501600.622288.1827085.6627010.16517.65227.14478.960.220.211001200.626591.7431805.4831543.98299.11135.31275.840.200.191001600.625677.2432238.2230776.40477.93180.57470.890.260.20Average24596.6429598.3929130.39410.54174.20383.530.200.18Table 9
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 34 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.639470.5656531.4449835.36658.54212.67449.650.300.21501600.639011.4455818.1648902.52851.15345.59727.260.300.201001200.645658.9664629.4058852.12475.71192.64389.790.290.221001600.643847.1464231.2258071.56745.12297.26635.570.320.24Average41997.0360302.5653915.39682.63262.04550.570.300.22From 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 Table10. 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.Table 10
Comparison for three algorithms on average deviation rate.
Instance with 16 item tasksInstance with 24 item tasksInstance with 34 item tasksRAT10.150.200.30RAT20.090.180.22Now 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.Table 11
The optimal solution for TPA-GA, TPA-ABC, and TSHOP applied to the instance with 34 item tasks.
Natural number coding of taskAccess order coding sequencesHandling station accessing sequencesCompletion time of taskBatch serial numberTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOP152911611406:0008:1006:261112631363707:1508:5207:11122323313333908:1008:5207:35222411251546908:3609:3207:352325313014311709:0210:1408:112426242119310609:0211:0508:572427152220910609:4912:1009:3234283213261610:1512:5709:323539252561610:1513:5209:32353102632666611:1514:3010:473631129157119811:5715:1711:3246312211617109512:0816:3212:174641322342910111113:1317:2112:5847414191330671114:0218:0913:435841520191666914:3718:5514:325851630203116315:1919:3015:1968617335496316:0819:3015:39786181662396316:0820:4515:397861927724108317:0421:1316:098962028831108318:1922:0816:0989621111114119:0622:5216:36910622912224119:2623:4217:3191062313431731120:1100:0818:11101072422321231020:5500:0819:021110725122422431021:2200:3820:071110726101727251021:5001:0020:071110727143328791022:3501:4621:22121172817261856523:2002:1222:1013118291827956200:0503:0122:37131183032810310200:2704:1623:17131283149832801:2705:0400:011313932341012112402:0705:4400:45141310337182585602:5706:1101:24151311348142687603:5206:5602:24151411We 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.Figure 11
Comparison for solution performance with iteration number of algorithms under different size problems. (a) Instance with 16 item tasks. (b) Instance size with 24 item tasks. (c) Instance with 34 item tasks.
(a)(b)(c)


## 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 Figure10. 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.Figure 10
Hybrid siding network in railway terminal.Table 1
Travel time of shunting engine between a pair of handling stations in railway terminal (unit: min).
u0123456789101100181920191920262527292211802027272538444345474021920024282739454446484132027240262240464547494241927282602739454446484151925272227039454446484162038394039390282826494272644454645452802624554882543444544442826023544792745464746462624230564910294748494848495554560511122404142414142484749510Table 2
Handling station capacity (unit: wagons per day).
u1234567891011Zu1836222425282622212132Table 3
Data of arrival train and wagon flow in railway terminal.
Natural number coding of taskSerial number of trainCompletion time point of train break-upHandling stationsNumber of local wagon groupOperation typesHandling timeCompletion time point of handling operationsDeparture time window of outbound train coupled by the local wagon group1106:00111155—[9:20, 10:00]2106:00111255—[9:20, 10:00]3106:00312160——4106:00312260——5106:00615175—[12:00, 13:00]6106:00615275—[12:00, 13:00]7106:00811155——8106:00811255——9206:2628140——10206:2628240——11206:26410150——12206:26410250——13206:26712160——14206:26712260——15206:26915175—[11:40, 12:40]16206:26915275—[11:40, 12:40]17307:3059145——18307:3059245——19307:3067135——20307:3067235——21307:301013165—[11:00, 12:10]22307:301013265—[11:00, 12:10]23408:1035130——24408:1035230——25408:10612160——26408:10612260——27408:101015175——28408:101015275——29408:10119145—[11:30, 12:00]30408:10119245—[11:30, 12:00]310—3723506:30—320—61427007:00—330—91527507:00[10:20, 11:00]340—11824006:20—


## 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.Table 4
Testing scope of three parameters.
ParameterRangeIteration population size ℑ{50, 80, 100}Iteration maximum number ℜ{100, 120, 160}Disturbance parameterφ{0.3, 0.6, 0.9}Table 5
Testing results of parameter combination.
No.ℑℜφOFVCPU time1501000.33909.7360.752501200.33897.4963.183501600.33822.6161.334501000.63847.3155.735501200.63870.0758.626501600.63655.8660.467501000.93927.0762.538501200.93909.7360.739501600.93897.4961.4110801000.33954.3784.6511801200.33943.3481.2312801600.33898.6584.3513801000.63899.1783.6614801200.63903.5683.5715801600.63843.7780.6216801000.93878.0482.5417801200.93922.4181.7118801600.93854.8785.43191001000.33785.26108.69201001200.33764.51108.85211001600.33671.27109.12221001000.63783.75110.31231001200.63742.61108.09241001600.63753.71107.11251001000.93821.01109.66261001200.93782.56109.73271001600.93756.90110.36In Table5, 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 Figure10. 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 byRAT = [OFV (TSHOP) − OFV (CPLEX)]/(OFV (CPLEX)) × 100%.The results of performance comparison are displayed in Table6.Table 6
Performance for CPLEX and TSHOP applied to small-size problem.
Tasks sizeCPLEXTSHOPRAT (%)OFVCPU timeOFVCPU time83655.86569.733655.8660.460.00106676.331493.726917.81109.373.621210127.373315.4510796.42153.286.6114#3600.013617.36197.33#Average6819.852515.7810393.14146.545.24For 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 byRAT1 = [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 byRAT2 = [OFV (TPA-ABCA) − OFV (TSHOP)]/OFV (TSHOP) × 100%.Here we chose control parameter combination from Table4 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.Table 7
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 16 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.614246.9416077.1815268.66240.07106.13199.920.130.07501600.614047.1616380.3215366.64350.05167.62311.830.170.091001200.616978.2419411.9218526.36212.2696.50165.300.140.091001600.616601.6819522.218267.66331.52150.28272.930.180.10Average15468.5017847.916857.33283.47130.13237.500.150.09Table 8
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 24 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.623829.4027264.2027191.00347.47153.76308.450.140.14501600.622288.1827085.6627010.16517.65227.14478.960.220.211001200.626591.7431805.4831543.98299.11135.31275.840.200.191001600.625677.2432238.2230776.40477.93180.57470.890.260.20Average24596.6429598.3929130.39410.54174.20383.530.200.18Table 9
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 34 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.639470.5656531.4449835.36658.54212.67449.650.300.21501600.639011.4455818.1648902.52851.15345.59727.260.300.201001200.645658.9664629.4058852.12475.71192.64389.790.290.221001600.643847.1464231.2258071.56745.12297.26635.570.320.24Average41997.0360302.5653915.39682.63262.04550.570.300.22From 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 Table10. 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.Table 10
Comparison for three algorithms on average deviation rate.
Instance with 16 item tasksInstance with 24 item tasksInstance with 34 item tasksRAT10.150.200.30RAT20.090.180.22Now 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.Table 11
The optimal solution for TPA-GA, TPA-ABC, and TSHOP applied to the instance with 34 item tasks.
Natural number coding of taskAccess order coding sequencesHandling station accessing sequencesCompletion time of taskBatch serial numberTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOP152911611406:0008:1006:261112631363707:1508:5207:11122323313333908:1008:5207:35222411251546908:3609:3207:352325313014311709:0210:1408:112426242119310609:0211:0508:572427152220910609:4912:1009:3234283213261610:1512:5709:323539252561610:1513:5209:32353102632666611:1514:3010:473631129157119811:5715:1711:3246312211617109512:0816:3212:174641322342910111113:1317:2112:5847414191330671114:0218:0913:435841520191666914:3718:5514:325851630203116315:1919:3015:1968617335496316:0819:3015:39786181662396316:0820:4515:397861927724108317:0421:1316:098962028831108318:1922:0816:0989621111114119:0622:5216:36910622912224119:2623:4217:3191062313431731120:1100:0818:11101072422321231020:5500:0819:021110725122422431021:2200:3820:071110726101727251021:5001:0020:071110727143328791022:3501:4621:22121172817261856523:2002:1222:1013118291827956200:0503:0122:37131183032810310200:2704:1623:17131283149832801:2705:0400:011313932341012112402:0705:4400:45141310337182585602:5706:1101:24151311348142687603:5206:5602:24151411We 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.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 Figure10. 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 byRAT = [OFV (TSHOP) − OFV (CPLEX)]/(OFV (CPLEX)) × 100%.The results of performance comparison are displayed in Table6.Table 6
Performance for CPLEX and TSHOP applied to small-size problem.
Tasks sizeCPLEXTSHOPRAT (%)OFVCPU timeOFVCPU time83655.86569.733655.8660.460.00106676.331493.726917.81109.373.621210127.373315.4510796.42153.286.6114#3600.013617.36197.33#Average6819.852515.7810393.14146.545.24For 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 byRAT1 = [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 byRAT2 = [OFV (TPA-ABCA) − OFV (TSHOP)]/OFV (TSHOP) × 100%.Here we chose control parameter combination from Table4 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.Table 7
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 16 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.614246.9416077.1815268.66240.07106.13199.920.130.07501600.614047.1616380.3215366.64350.05167.62311.830.170.091001200.616978.2419411.9218526.36212.2696.50165.300.140.091001600.616601.6819522.218267.66331.52150.28272.930.180.10Average15468.5017847.916857.33283.47130.13237.500.150.09Table 8
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 24 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.623829.4027264.2027191.00347.47153.76308.450.140.14501600.622288.1827085.6627010.16517.65227.14478.960.220.211001200.626591.7431805.4831543.98299.11135.31275.840.200.191001600.625677.2432238.2230776.40477.93180.57470.890.260.20Average24596.6429598.3929130.39410.54174.20383.530.200.18Table 9
Performance for TSHOP, TPA-GA, and TPA-ABCA applied to the instance with 34 item tasks.
Parameter setOFVCPU time (s)RAT1 (%)RAT2 (%)ℑℜφTSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCA501200.639470.5656531.4449835.36658.54212.67449.650.300.21501600.639011.4455818.1648902.52851.15345.59727.260.300.201001200.645658.9664629.4058852.12475.71192.64389.790.290.221001600.643847.1464231.2258071.56745.12297.26635.570.320.24Average41997.0360302.5653915.39682.63262.04550.570.300.22From 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 Table10. 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.Table 10
Comparison for three algorithms on average deviation rate.
Instance with 16 item tasksInstance with 24 item tasksInstance with 34 item tasksRAT10.150.200.30RAT20.090.180.22Now 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.Table 11
The optimal solution for TPA-GA, TPA-ABC, and TSHOP applied to the instance with 34 item tasks.
Natural number coding of taskAccess order coding sequencesHandling station accessing sequencesCompletion time of taskBatch serial numberTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOPTPA-GATPA-ABCATSHOP152911611406:0008:1006:261112631363707:1508:5207:11122323313333908:1008:5207:35222411251546908:3609:3207:352325313014311709:0210:1408:112426242119310609:0211:0508:572427152220910609:4912:1009:3234283213261610:1512:5709:323539252561610:1513:5209:32353102632666611:1514:3010:473631129157119811:5715:1711:3246312211617109512:0816:3212:174641322342910111113:1317:2112:5847414191330671114:0218:0913:435841520191666914:3718:5514:325851630203116315:1919:3015:1968617335496316:0819:3015:39786181662396316:0820:4515:397861927724108317:0421:1316:098962028831108318:1922:0816:0989621111114119:0622:5216:36910622912224119:2623:4217:3191062313431731120:1100:0818:11101072422321231020:5500:0819:021110725122422431021:2200:3820:071110726101727251021:5001:0020:071110727143328791022:3501:4621:22121172817261856523:2002:1222:1013118291827956200:0503:0122:37131183032810310200:2704:1623:17131283149832801:2705:0400:011313932341012112402:0705:4400:45141310337182585602:5706:1101:24151311348142687603:5206:5602:24151411We 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.Figure 11
Comparison for solution performance with iteration number of algorithms under different size problems. (a) Instance with 16 item tasks. (b) Instance size with 24 item tasks. (c) Instance with 34 item tasks.
(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.



---

*Source: 6713617-2022-05-11.xml*