Optimizing the High-Level Maintenance Planning Problem of the Electric Multiple Unit Train Using a Modiﬁed Particle Swarm Optimization Algorithm

: Electric multiple unit (EMU) trains’ high-level maintenance planning is a discrete problem in mathematics. The high-level maintenance process of the EMU trains consumes plenty of time. When the process is undertaken during peak periods of the passenger ﬂow, the transportation demand may not be fully satisﬁed due to the insufﬁcient supply of trains. In contrast, if the process is undergone in advance, extra costs will be incurred. Based on the practical requirements of high-level maintenance, a 0–1 programming model is proposed. To simplify the description of the model, candidate sets of delivery dates, i.e., time windows, are generated according to the historical data and maintenance regulations. The constraints of the model include maintenance regulations, the passenger transportation demand, and capacities of workshop. The objective function is to minimize the mileage losses of all EMU trains. Moreover, a modiﬁed particle swarm algorithm is developed for solving the problem. Finally, a real-world case study of Shanghai Railway is conducted to demonstrate the proposed method. Computational results indicate that the (approximate) optimal solution can be obtained successfully by our method and the proposed method signiﬁcantly reduces the solution time to 500 s.


Introduction
In China, high-speed railway has become a priority option for the long trip due to its convenience and comfortableness, and it account for 60% of rail total passenger traffic. It has long been a difficult problem for the China Railway (CR) that supply enough available Electric Multiple Unit (EMU) trains in tourist rush seasons to fulfill heavy transportation tasks. The EMU trains are the unique vehicles running on the China high-speed railway. In addition, they have a high purchase cost and complicated maintenance regulations. Therefore, the problem becomes much worse.
The gigh-level maintenance planning (HMP) is an important part of the operation and maintenance management for the EMU trains, which covers the third-level maintenance, the fourth-level maintenance and the fifth-level maintenance. Due to the complex regular preventive maintenance, uneven distribution of passenger flow, limited maintenance capacity of workshop and several weeks for maintenance service time, the EMU trains' HMP needs to be scheduled in advance, of which the planning horizon lasts for a natural year or more. The HMP's aim is to provide enough In the CR system, according to the maintenance regulations [1], the EMU trains that have been put into operation in the first period will undergo the high-level maintenance (HM) procedures together in the next few years. When lots of EMU trains are sent to workshop in tourist rush season, it will lead to the travel needs of passengers cannot be met. Not only the operating income but also the traveler satisfaction level on high-speed railway transportation will decrease. Meanwhile, it will increase the maintenance costs when the HM procedures are undergone in advance. Due to the limited capacities of workshops, the maintenance service time will be prolonged when plenty of EMU trains undergo the maintenance procedures together. Therefore, to scientifically formulate the HMP that meets the travel demands and reduces maintenance cost as much as possible is a complicated combinatorial optimization problem that needs to be solved in the field of high-speed railway operation. For the maintenance system, Stuchly et al. [2] and Rezvanizaniani et al. [3] developed a condition based maintenance management system, Shimada [4] introduced the accident prevention maintenance system, and Cheng and Tsao [5] proposed a preventive and corrective maintenances system. These maintenance systems can reduce the maintenance costs and improve the utilization efficiency of EMU train.
Many experts and scholars study the first-level and the second-level maintenances. Maróti and Kroon [6,7] proposed that adjust the operation schedule ahead of time to ensure the maintenance procedures can be undertook in time. Tsuji et al. [8] developed a novel approach based on ant colony optimization and Wang et al. [9] designed an algorithm based column generation to solve the EMU train maintenance plan problem. Giacco et al. [10] integrated the rolling stock circulation problem and short-term maintenance planning, and a mixed-integer linear-programming formulation was proposed.
In other fields, scheduled maintenance planning problems also have been researched. Ziarati et al. [11], Lingaya et al. [12], and Wang et al. [13] studied the locomotive operation and maintenance. Moudania and Félix [14] and Mehmet and Bilge [15] developed the aircraft operation and maintenance. Budai et al. [16] researched the long-term planning of railway maintenance works. In addition, Grigoriev et al. [17] tried to find the length of maintenance plan to minimize the total operation costs.
There are a few relevant literatures available for the EMU trains HMP problem. Lin et al. [18] designed a state function to show the state of EMU train on each day during the planning horizon, In the CR system, according to the maintenance regulations [1], the EMU trains that have been put into operation in the first period will undergo the high-level maintenance (HM) procedures together in the next few years. When lots of EMU trains are sent to workshop in tourist rush season, it will lead to the travel needs of passengers cannot be met. Not only the operating income but also the traveler satisfaction level on high-speed railway transportation will decrease. Meanwhile, it will increase the maintenance costs when the HM procedures are undergone in advance. Due to the limited capacities of workshops, the maintenance service time will be prolonged when plenty of EMU trains undergo the maintenance procedures together. Therefore, to scientifically formulate the HMP that meets the travel demands and reduces maintenance cost as much as possible is a complicated combinatorial optimization problem that needs to be solved in the field of high-speed railway operation. For the maintenance system, Stuchly et al. [2] and Rezvanizaniani et al. [3] developed a condition based maintenance management system, Shimada [4] introduced the accident prevention maintenance system, and Cheng and Tsao [5] proposed a preventive and corrective maintenances system. These maintenance systems can reduce the maintenance costs and improve the utilization efficiency of EMU train.
Many experts and scholars study the first-level and the second-level maintenances. Maróti and Kroon [6,7] proposed that adjust the operation schedule ahead of time to ensure the maintenance procedures can be undertook in time. Tsuji et al. [8] developed a novel approach based on ant colony optimization and Wang et al. [9] designed an algorithm based column generation to solve the EMU train maintenance plan problem. Giacco et al. [10] integrated the rolling stock circulation problem and short-term maintenance planning, and a mixed-integer linear-programming formulation was proposed.
In other fields, scheduled maintenance planning problems also have been researched. Ziarati et al. [11], Lingaya et al. [12], and Wang et al. [13] studied the locomotive operation and maintenance. Moudania and Félix [14] and Mehmet and Bilge [15] developed the aircraft operation and maintenance. Budai et al. [16] researched the long-term planning of railway maintenance works. In addition, Grigoriev et al. [17] tried to find the length of maintenance plan to minimize the total operation costs.
There are a few relevant literatures available for the EMU trains HMP problem. Lin et al. [18] designed a state function to show the state of EMU train on each day during the planning horizon, and a non-linear 0-1 programming model and its solution strategy was proposed. Wu et al. [19] proposed a time-state network to optimize the EMU trains HMP problem. Li et al. [20] presented a forecast method to estimate the maintenance quantity of the EMU trains in future.
Compared to the existing researches, a 0-1 programming model and solution strategy were proposed. In the mathematic model, all necessary regulations and practical constraints were considered.
The remainder of this paper is organized as follows. The problem description of HMP in the CR system is presented in Section 2. In Section 3, a 0-1 programming model is proposed, and then the solution algorithm on the basis of particle swarm algorithm is designed in Section 4. An empirical case is provided to verify the effectiveness of the model and the algorithm in Section 5. The last section gives some conclusions and the possible areas of further research.

HMP Problem at CR
The HMP, a typical discrete system, is a tactical plan that determines when the EMU trains to undergo the high-level maintenance. The length of the planning horizon lasts for about one year. Each train undergoes the high-level maintenance at most once during the planning horizon because of the interval between two adjacent HM processes is longer than the span of planning horizon; the order of the maintenance level for a train is the third-level, the fourth-level, the third-level, the fifth level, the third-level and so on until they are scrapped [1]. Therefore, the level of the high-level maintenance can be deduced beforehand according to the records and regulations of maintenance. According to the maintenance regulations, each train has a time window during which the train can be delivered to workshop on any day. The lower bound of the time window is the earliest date on which a train can be delivered to workshop while the upper bound is the latest date. The detailed problem descriptions can be referred to Lin et al. [18] and Wu et al. [19]. Here, we focus on the generation steps of the time window.
The estimated time of arrival (eta) of the HM can be easily calculated [20]. In this process, the average daily operating mileage is used to describe the daily usage of EMU trains before HM procedures. The notations used in the generation process of the time window are listed in Table 1. The type of the EMU train e; e(g) The maintenance level of the EMU train e; The maximum value of the difference between the actual operating mileage before the HM and the eta e for the train of which the type is m and maintenance level is g; The maximum value of the difference between the eta e and the actual operating mileage before the HM for the train of which the type is m and the maintenance level is g; l e The average daily operating mileage for the EMU train e before HM.
The time window of the start time for the HM can be generated as follows.
Step Step 5. If E = φ, turn to Step 6, otherwise turn to Step 1; Step 6. Set E = E , over. The EMU train's time window is continuous if we set the "day" as the minimum time unit. In this way, the time window can be presented by a time interval.
In addition, the HMP aims to ensure that there are enough well-conditioned trains to meet the passenger transport demand. We set a maximum HM rate to guarantee it. The HM rate is the ratio of the number of trains in HM state to the fleet size.

Mathematical Model of the HMP Problem
In this section, we propose a 0-1 programming model for the HMP problem. The constraints of the model include the maintenance interval, the passenger transportation demand, and the capacity of workshop. The objective function is to minimize the mileage losses of all EMU trains.

Notations
The all notations that used in the model are listed in Table 2. Table 2. Notations used in the model.

Notations Definition Indices e
The index of EMU trains; m The index of types for the EMU trains; t The index of dates during the planning horizon; g The index of the maintenance level; The set of all types; m ∈ M T The set of all dates during the planning horizon; t ∈ T G The set of all maintenance levels; g ∈ G Input parameters The maintenance level of the EMU train e; e(m) The type of the EMU train e; c The unit penalty fee for the unused mileage before the HM; The length of the planning horizon; |E| The number of the EMU trains which need to be maintained during the planning horizon;

Notations Definition
Decision Variables x t e Binary variable, indicates whether the EMU train e selects the t-th day to start the HM procedures during the planning horizon, x t e = 1 if yes, x t e = 0 otherwise; y t e Binary variable, indicates whether the EMU train e in the e(g)-th level maintenance state on the t-th day during the planning horizon, y t e = 1 if yes, y t e = 0 otherwise; z t e Binary variable, indicate whether the EMU train e in the e(g)-th level delivery interval on the t-th day during the planning horizon, z t e = 1 if yes, z t e = 0 otherwise;

Optimization Objective
The objective function of the mathematical model is to maximize the service efficiency of the EMU trains, i.e., to minimize the unutilized mileage. In this way, the objective function can be presented as follows.

Constraints Analysis
According to Section 2, each train e can choose one and only one delivery date during the time window. This is the uniqueness constraint.
Any time t out of the time window for the train e cannot be selected.
The θ t is a variable according to travel demand which should be guaranteed. Theoretically, the value of θ t is different for each day during the planning horizon. To describe this requirement, a set of constraints established as follows.
Because of the various itineraries, the number of the trains with the specific type in the HM state must be less than the given threshold value on the t-th day. Constraints in this respect can be expressed as follows.
The number of trains in each level of the HM state should not be exceeded the capacity of workshops. A set of constraints can be listed as follows.
Meanwhile, restricted by the limited resources, only a few trains are permitted to enter the workshop over several days. This situation can be described in the form of mathematical inequalities as follows.
Symmetry 2018, 10, 0349 6 of 14 In addition, the logical relationships between those three sets of decision variables are presented as follows.
Finally, all of the decision variables are binary variables.
x t e , y t e , z t e ∈ {0, 1} ∀ e ∈ E, t ∈ T

Model Construction
On the basis of above analysis, a 0-1 programming model for the EMU train HMP problem is proposed as follows.  x t e = 1 ∀ e ∈ E x t e = 0 ∀ e ∈ E, t / ∈ [WS e , WE e ] and t ∈ T x t e , y t e , z t e ∈ {0, 1} ∀ e ∈ E, t ∈ T We can see from the HMP model that the number of all decision variables equals to 3 × |T| × |E| that is the product of two factors: the time span of the planning horizon and the number of trains. But the search space will reach |T|ˆ|E| according to the model. It is too complicated to be solved by using CPLEX or Gurobi within a reasonable time. Thus, a meta-heuristic solution strategy based on the particle swarm optimization (PSO) algorithm is designed to address this problem.

Modified Particle Swarm Optimization Algorithm
PSO algorithm is a population based stochastic optimization technique motivated by social behavior of organisms. PSO algorithm has the advantage of fast convergence speed and high accuracy solution and it is easy to be applied in most areas [21], which makes it attract great attention from researchers. The algorithm can also be used in solving the combinatorial optimization problem [22,23]. In this section, we present a modified particle swarm optimization (MPSO) algorithm to solve the HMP model based on analysis and preprocess.

Processing of Model Constraints
The constraint conditions of the HMP model need to be processed before applying the MPSO algorithm. Constraints (2), (3), (8)- (12) are the logical relationships, and they can be observed by the specific encoding rules (see the latter section), while the others can be removed by the penalty function method. For the value of the penalty factor, it is necessary to combine the actual application of the EMU trains and the strength of the constraint. Among the Inequations (4)-(7), the Inequation (7) has the strongest constraint and the Inequation (4) has the weakest constraint. The relationship of the penalty coefficient is λ 4 > λ 3 > λ 2 > λ 1 . The optimization model can be presented as follows.

General Particle Swarm Optimization Algorithm
In general PSO algorithm, the basic update equations of the velocity and position of the particles are as follows: where V i (r) and P i (r) denote the velocity and the position, respectively, for particle i in the r − th iteration. HB i (r) denotes the best position in the history for particle i by the end of the r − th iteration; GB(r) denotes the best position in the history for all of the particles by the end of the r − th iteration.
ω denotes the inertia weight; c 1 denotes the self-learning factor; c 2 denotes the social learning factor; ξ and η are the random numbers in [0, 1].

Inertia Weight
In order to make the particles have a better search ability in the early stage of evolution and have a better development ability in the later stage of the evolution, the linear time-varying inertia weight is adopted in this paper. The inertia weight can be calculated as follows.
where ω(r) denotes the inertia weight at the r-th iteration. ω max and ω min denotes the maximum inertia weight and the minimum inertia weight, respectively. And MAXR denotes the maximum number of the evolution iterations.

Learning Factor
In the same way, in order to make the particles strengthen the global search ability in the early stage and converge to the global optimum in the later period, we decrease the self-learning factor and increase the social learning factor continuously during the process of optimization. The calculation formulae are as follows.
where c 1 (r) and c 2 (r) denote the self-learning factor and the social learning factor in the r − th iteration. c 1 and c 1 denote the initial value and the final value for the self-learning factor; c 2 and c 2 denote the initial value and the final value for the social learning factor; they are the constants.

Update Equations
The particle continuously updates its position in the search space at an unfixed speed. The velocity represents the variation of position in magnitude and direction like the definition in classical physics, and it has the same dimension as the position. Let hbest i (r) denote the corresponding fitness value of HB i (r). Let gbest(r) denote the corresponding fitness value of GB(r). The fitness value can be calculated by the formula (13). The update equations of the position and velocity for the particle i are as follows.
V i (r + 1) = ω(r)·V i (r) + c 1 (r)·ξ(r)·(HB i (r) − P i (r)) + c 2 (r)·η(r)·(GB(r) − P i (r)) (19) where ξ(r) and η(r) are the random numbers in [0,1]. Each dimension of the velocity of a particle is limited to an interval [−V max , V max ], and if it is out of the interval, we set the boundary value of the interval as the actual velocity component. Similarly, the components of the position vector are limited to the time window for each train. Let pbest i (r) denote the fitness value for P i (r).

Encoding Rules and Initial Solution
It is a crucial step to make the particle of the MPSO and the solution of a certain problem correspond with each other. We use a particle to represent an overall HMP for all of the EMU trains. According to Equations (2) and (3), we set the dimensionality of a particle to |E|. The value of each dimension represents the start time of the HM for the corresponding EMU train. The detailed description is shown in Figure 2 with the help of schematic diagram.
It is a crucial step to make the particle of the MPSO and the solution of a certain problem correspond with each other. We use a particle to represent an overall HMP for all of the EMU trains. According to Equations (2) and (3), we set the dimensionality of a particle to E . The value of each dimension represents the start time of the HM for the corresponding EMU train. The detailed description is shown in Figure 2 with the help of schematic diagram.

Algorithm Steps
Step 1. Generate the time window for all of the EMU trains that need to be maintained during the planning horizon, turn to Step 2.
Step 2. Initialization. Assign values to related parameters including the size of particle swarm I , max  , min  , 1 c , 1 c , 2 c , 2 c , MAXR and max V . Generate the initial solution (0) i P according to Section 4.3.1, and generate the initial velocity (0) i V randomly. Set 0 r  , turn to Step 3.
Step 3. Calculate the fitness value () i pbest r of each particle, turn to Step 4.
Step 4. Compare the fitness values of ()  In Figure 2, t e denotes the start time of the HM for the e − th EMU train. Therefore, x t e = 1 when the value of t e is determined, and the t in

Algorithm Steps
Step 1. Generate the time window for all of the EMU trains that need to be maintained during the planning horizon, turn to Step 2.
Step 2. Initialization. Assign values to related parameters including the size of particle swarm I, ω max , ω min , c 1 , c 1 , c 2 , c 2 , MAXR and V max . Generate the initial solution P i (0) according to Section 4.3.1, and generate the initial velocity V i (0) randomly. Set r = 0, turn to Step 3.
Step 3. Calculate the fitness value pbest i (r) of each particle, turn to Step 4.
Step 4. Compare the fitness values of P i (r) with HB i (r). If pbest i (r) < hbest i (r), then HB i (r) = P i (r), hbest i (r) = pbest i (r). Turn to Step 5.
Step 7. Update V i (r) and P i (r) according to Formulas (19) and (20). If any dimension in V i (r) out of [−V max , V max ], we set the boundary value of the interval as an actual value. If any dimension in P i (r) out of [WS e , WE e ], we set the boundary value of the time window as an actual value, turn to Step 8.
Step 9. Make GB(r) feasible according to the HMP model, and output GB(r) and gbest(r). Over.

Case Study
In this section, we implement the proposed method to solve a practice problem. The detailed description of the case can be found in the literature [20]. The proposed model is solved by the commercial optimization solver, e.g., Gurobi, as well as the MPSO algorithm. The exact method is coded in Python 2.7 and implemented within Spyder 3.1.4 and the MPSO algorithm is implemented in C++. All the computational experiments are conducted on the computer with Intel Core i5-6200U CPU and 8 GB RAM.
According to the literature [20], some parameters are valued as follows. Inv = 115, |E| = 60, |T| = 533. To protect data confidentiality, we can only generate the time window in advance for each train and use an ID number replace the train. The initial conditions of trains when |T| = 0, which include the train ID, type, the average daily operating mileage, the time window, the maintenance level and the maintenance service time, are listed in Table 3. From Table 3, it can be seen that there are 60 EMU trains will undergo HM procedures during the planning period. Among them, there are 101 candidate dates at most, and the minimum values are only 36 candidate dates.
During the planning horizon, the HM rate is valued as follows: All of these trains undergo the HM procedures in the factory. The maximal capacity is ten, i.e., b 1 + b 2 + b 3 = 10. Meanwhile, the receiving capability per day for a factory is limited, e.g., According to the demand, d t m1 = d t m2 = 4, d t m3 = 10 (t ∈ T). Therefore, the Formula (13) can be converted to the following form.
In addition, the values of other parameters are as follows. c = 0.001, λ 1 = 100, λ 2 = 500, λ 3 = 800, λ 4 = 1000; I = 1000, ω max = 1.2, ω min = 0.8, c 1 = 2.5, c 1 = 0.5, c 2 = 0.5, c 2 = 2.5, V max = (WE e -WS e +1)/2 MAXR = 1000. Based on the data given above, the HMP model is solved by the proposed algorithm. The program runs for 500 s. The returned optimal fitness value is 3,213,121. The curve of the optimal fitness value in the iterative process is depicted in Figure 3.   Based on the data given above, the HMP model is solved by the proposed algorithm. The program runs for 500 s. The returned optimal fitness value is 3,213,121. The curve of the optimal fitness value in the iterative process is depicted in Figure 3. As can be seen in Figure 3, in the first 300 iterations, the algorithm has a strong search capability, which can effectively avoid the occurrence of the premature phenomena; and from about the 300th to the 690th iteration, the development ability of the algorithm is strengthened, which is helpful to search the optimal solution; the fitness value remains the same in the last 200 iterations, indicating that the (approximate) optimal solution for the HMP problem has been generated. The (approximate) optimal solution is listed in Table 4, and the first column is the train's ID; the second column is the start date of HM procedures denoted by e t .   As can be seen in Figure 3, in the first 300 iterations, the algorithm has a strong search capability, which can effectively avoid the occurrence of the premature phenomena; and from about the 300th to the 690th iteration, the development ability of the algorithm is strengthened, which is helpful to search the optimal solution; the fitness value remains the same in the last 200 iterations, indicating that the (approximate) optimal solution for the HMP problem has been generated. The (approximate) optimal solution is listed in Table 4, and the first column is the train's ID; the second column is the start date of HM procedures denoted by t e . We present the daily HM rates from the (approximate) optimal solution, and compare those with the predefined maintenance rate thresholds (see Figure 4). We present the daily HM rates from the (approximate) optimal solution, and compare those with the predefined maintenance rate thresholds (see Figure 4). From Figure 4, we can see that the HM rate of the optimal solution remains below the threshold in each period. Therefore, the proposed HMP model and the algorithm can meet the travel demands. However, the HM rate from the unscheduled solution fluctuates sharply.
In addition, to compare the performance of the proposed algorithm, we solve the HMP model by the Gurobi solver because the HMP model is linear. The detailed numerical comparison results of two solution methods, i.e., the Gurobi and the MPSO algorithm, are shown in Table 5. We can conclude that our proposed approach is efficient and effective from Table 5. This result demonstrates that the proposed algorithm is more efficient than the commercial optimization solver with reducing 84.31% in the solution time consumption. From Figure 4, we can see that the HM rate of the optimal solution remains below the threshold in each period. Therefore, the proposed HMP model and the algorithm can meet the travel demands. However, the HM rate from the unscheduled solution fluctuates sharply.

Conclusions
In addition, to compare the performance of the proposed algorithm, we solve the HMP model by the Gurobi solver because the HMP model is linear. The detailed numerical comparison results of two solution methods, i.e., the Gurobi and the MPSO algorithm, are shown in Table 5. We can conclude that our proposed approach is efficient and effective from Table 5. This result demonstrates that the proposed algorithm is more efficient than the commercial optimization solver with reducing 84.31% in the solution time consumption.

Conclusions
This paper researches the electric multiple unit train's HMP problem. A 0-1 integer programming model and a modified particle swarm optimization algorithm are proposed. The objective function of the model is to minimize the unutilized mileage for all trains, and the model considers the necessary regulations and practical constraints, including passenger transport demand, workshop maintenance capacity, and maintenance regulations. A real-world instance demonstrates that the proposed method can efficaciously obtain the (approximate) optimal solution (see Table 4), and the solution strategy significantly reduces the solution time to 500 s (see Table 5). This result also demonstrates that the proposed algorithm is more efficient than the commercial optimization solver with reducing 84.31% in the solution time consumption. Optimize the workshop's overhaul process to shorten the maintenance service time is needed in future research.