A Robust Optimization Model for Multi-Period Railway Network Design Problem Considering Economic Aspects and Environmental Impact ^†

Abstract

The railway network design problem is one of the critical issues in the transportation sector due to its significance and variety of necessary applications. The major issue of this field relates to the decision of whether to increase the railways’ capacity or construct a new route to meet demand. Although the budget is a great concern of the managers for making such a decision, environmental factors should be necessarily included in the decision-making process. Therefore, this research proposes a novel robust bi-objective mixed-integer linear programming (MILP) model to simultaneously minimize the total cost and environmental impact under uncertain conditions and within a given time horizon. The proposed problem addresses strategic and operational decisions through railway project selection and product flow determination. To deal with the bi-objectiveness of the model and tackle the complexity of the problem, a nondominated sorting genetic algorithm (NSGA-II) is employed. The proposed NSGA-II could reach near-optimal Pareto solutions in a reasonable solution time and showed a reliable performance for being employed in large-sized instances. It also indicates that the proposed NSGA-II can be utilized for solving large-sized samples in a very short time.

1. Introduction

Railway projects are considered as one of the major construction projects in the countries and require high executive expenses such that construction costs are now a large part of the total construction credits of various countries. Railway transport is considered to belong among the modes of transport with the lowest negative impact on the environment and society [1]. The European Union (EU) has set a target objective of a modal shift from road freight transport towards other more sustainable transport modes of 30% by 2030 and at least 50% by 2050 for shipments over 300 km [1].

Rail freight transportation plays a critical role in the world’s economy and is an indisputable part of transportation in modern society. Duly, investment in the development or construction of a railway network, which comprises the selection of a series of stations along with the links joining them, is affirmed through the resulting savings in energy, time, and pollution emitted by trains. It is also regarded as a firm step towards sustainable development [2]. Therefore, properly selecting railway plans helps to properly plan the transportation sector and the correct orientation in the country’s development, increases the effectiveness of railway networks, and maximizes social and financial benefits. The railway planning problem can be divided into six subgroups: (i) railway network design problem (RNDP), (ii) crew scheduling problem, (iii) rolling stock assignment, (iv) train routing problem, (v) timetabling problem [3], and (vi) line planning problem [4].

However, there are also some challenges and questions to be addressed. The main challenge is the provision of the resources required to improve the current lines or build new lines. Moreover, pollution reduction is another important challenge which needs to be considered continuously. Therefore, the main question is how we can develop a mathematical model and solution method to consider real-world assumptions and find an optimal solution accordingly. Therefore, in this research, a multi-period railway network design problem is studied considering the economic aspects and robustness of the problem. The main goal is to minimize the total cost of designing the rail transportation network from the government’s point of view. Therefore, we consider various economic aspects of the problem, including direct/indirect costs and operating and construction costs, which must be addressed in many studies. On the other hand, the amount of greenhouse gas emissions associated with rail transport is another factor that has received less attention despite recent environmental laws and concerns [5]. To define the problem space, a set of parameters such as network line capacity [6], which plays a very key role in managing network transport demand, is considered. Moreover, double-track and single-track lines are analyzed independently. Another important issue is the budget available to operate the design project, which is spent on building new lines, improving existing lines, and maintenance. All these significantly impact the capacity parameter, and all of which are considered in the form of a category called capacity improvement projects. On the other hand, to make the model closer to real-world conditions, the cost of lost demand is also assumed, which plays a decisive role in determining the optimal values of the variables and the total cost. Finally, the most critical issue is to consider the degree of uncertainty of the demand parameter in the network, which is formulated based on the interval uncertainty and the proposed linear programming by [7]. Since demand is inherently uncertain, and based on the available information, the uncertainty of the interval in the form of defining the lower and upper bounds in an interval is the best way to deal with it; therefore, the robust optimization approach is employed, which has high efficiency in dealing with uncertainty.

The remainder of this paper is structured as follows. Section 2 reviews the relevant background of the discussed problem. Section 3 presents the problem assumptions and the proposed mathematical modelling of the problem. Section 4 describes the robust approach, and the problem’s robust counterpart is developed. The applied solution methods are discussed in Section 5. In Section 6, the parameter adjustment of the algorithm is investigated. The results and a discussion of the proposed algorithms are presented in Section 7. Finally, the conclusion is given in Section 8.

2. Background

This section reviews the studies related to the railway network design problem (RNDP). In an early study. In an early study [8], the authors indicated that railway projects could be completed in different phases, either gradually or all at once. To solve the problem, they employed a heuristic time sequence method. They also reduced the model size significantly by ignoring the small demands from the origin–destination (OD) matrix and preloading them onto their shortest path. This size reduction had little effect on the optimality level of the problem. Lin et al. [9] studied a railway network design problem, considering the project budget and routing costs. In their model, decisions such as establishing a new railway or extending the old railways were made. Yan et al. [10] extended the railway network design model in multiple stages. They also investigated the effects of investment on the railway project choice.

Arnold et al. [11] developed a mixed-integer linear programming model which aimed to determine the optimal location of rail terminals. They employed a heuristic approach to solve their model. Pazour et al. [12] designed a network model for high-speed railways and highways to reduce road congestion. The model objective was the minimization of travel time, and the model was uncapacitated. Sun et al. [13] developed a bi-level programming model for a transport network by considering multi-periods and demand uncertainty in the model. The model was able to optimize the form of the transport network and the construction time at the same time.

Laporte et al. [14] designed a railway network considering the network robustness. The study focused on designing railway lines so that the utility function is optimized. Marín and García-Ródenas [15] developed a model for rapid transit networks which aimed to maximize the transit demand and minimize the travel time. They also considered budget constraints and limitations for the physical networks. In another study, they studied an integrated model of network design and line planning without decreasing the system efficiency in the case of zero failure [16]. They solved the problem by considering the concept of robustness from users’ and operators’ viewpoints.

In the context of rapid transit network design, a Simulated Annealing (SA) algorithm was utilized to maximize trip coverage in the problem [17]. Garcia et al. [18] designed an infrastructure railway network and developed a robust model of the problem as well. They employed a greedy randomized adaptive search method for solving their problem. Canca et al. [19] presented a railway rapid transit network design in line with line planning. They considered and minimized the total cost, including construction costs, fleet purchase costs, train operational costs, and personnel costs. Syedvakili et al. [4] developed a mixed-integer model for designing a railway network from the government’s viewpoint. The model made decisions such as constructing a new line or improving the existing line, and the objective function was to minimize the total costs. A joint multi-objective, decision-making model for the capacity optimization and allocation of an urban railway network was offered by Wang et al. [20] using multi-source data. They applied their proposed methodology to the Beijing metro system, considering the operational costs of trains and the costs related to the waiting time of passengers. Seyedvakili et al. [21] proposed a multi-period railway network design model for the long-term planning of the railway network. Using multi-period railway network design problems can help decisionmakers consider various development projects, even when constructing new lines and improving existing lines.

Many countries are currently seeking effective methods to reduce carbon emissions. The transportation industry plays a significant role among the sources of air pollution. Lin et al. [22] analyzed the carbon emissions of railway and highway transport, estimated the environmental benefit of building a railway, and recommended that the railway network design be modeled as a bi-level programming problem, considering emissions.

In terms of a solution method, many studies have employed multi-objective metaheuristics to solve the problem by generating Pareto frontiers [23,24,25].

Considering the literature, we define our main contributions as follows:

• The possibility of expanding the existing capacity of the rail transport network through construction projects;
• Considering demand shortage and time horizon in the problem;
• Considering greenhouse gas emissions as a result of rail transportation operations and minimizing the total emissions in line with the total costs;
• Studying the uncertainty of demand using the robust optimization approach;
• Solving the problem bi-objectively using exact and meta-heuristic approaches.

3. Problem Definition

3.1. Problem Assumptions

The main assumptions of the problem are: (i) a set of new or existing directed arcs is considered in the network, (ii) the budget is limited, (iii) a planning period (time horizon) is defined, (iv) each arc or rail line has a limited capacity, (v) the capacity of any existing arc or rail line can be increased by defining a new project at the beginning of the time horizon, (vi) each project has a cost that is calculated based on the available budget, (vii) lost demand is considered in the problem, and (viii) the rail freight demand between each pair of origin–destination nodes (OD) is uncertain.

3.2. Mathematical Modelling

The proposed model investigates the effects of two important factors, including planning period and pollution reduction, which have not been addressed previously. The mathematical modelling for the proposed railway network design problem is presented as follows:

Sets

$A$	Set of directed arcs (rail lines),
$A^{\prime }$	Set of new directed arcs (rail lines),
$K$	Set of origin–destination (OD) pairs,
$P$	Set of capacity improvement projects,
$T$	Set of time periods.

Parameters

$B$	The available budget for projects,
${CAP}_{ij}$	The capacity of arc $(i,j)$,
${CAP}_{ijt}^p$	The added capacity to the arc $(i,j)$ as a result of project p at period $t$,
${CC}_p$	Fixed cost of conducting project $p$,
${CO}_{ijt}$	The unit operational cost for arc $(i,j)$ at period $t$,
${GHG}_{ij}$	The amount of greenhouse emissions on the arc $(i,j)$,
${{CU}_t}^{od}$	The cost of lost demand for each pair of ODs at period $t$,
${D_t}^{od}$	The demand between each pair of ODs at period $t$,
${\delta }_{ij}^p$	1 if arc $(i,j)$ belongs to the project p, otherwise 0.

Decision variable

${{ud}_t}^{od}$	Lost demand for each pair of ODs at period $t$,
$x_{ijt}^{od}$	The flow of each OD in arc $(i,j)$ at period $t$,
$y_p$	1 if project $p$ is conducted, otherwise 0.

Regarding the sets, parameters, and variables, the proposed mathematical modelling is defined.

$$\mathrm{minimize} TC=\sum _{t\in T}\sum _{od\in K}\sum _{(i,j)\in A\cup A^{\prime }}{CO}_{ijt}{x}_{ijt}^{od}+\sum _{t\in T}\sum _{od\in K}{{CU}_t}^{od}{{ud}_t}^{od}$$

(1)

$$\mathrm{minimize} GHG=\sum _{t\in T}\sum _{od\in K}\sum _{(i,j)\in A\cup A^{\prime }}{GHG}_{ij} x_{ijt}^{od}$$

(2)

Subject to

$$\sum _{p\in P}{CC}_p{y}_p\le B,$$

(3)

$$\sum _j{x}_{ijt}^{od}-\sum _j{x}_{jit}^{od}={D_t}^{od}-{{ud}_t}^{od}\hspace{1em}\hspace{1em}\forall i=o;od\in K;t\in T,$$

(4)

$$\sum _j{x}_{jit}^{od}-\sum _j{x}_{ijt}^{od}={D_t}^{od}-{{ud}_t}^{od}\hspace{1em}\hspace{1em}\forall i=d;od\in K;t\in T,$$

(5)

$$\sum _j{x}_{ijt}^{od}-\sum _j{x}_{jit}^{od}=0\hspace{1em}\hspace{1em}\forall i\ne od;od\in K;t\in T,$$

(6)

$$\sum _{od\in K}{x}_{ijt}^{od}\le {CAP}_{ij}+\sum _{p\in P}{CAP}_{ijt}^p{\delta }_{ij}^p{y}_p\hspace{1em}\hspace{1em}\forall i,j\in A\cup A^{\prime };t\in T,$$

(7)

$$x_{ijt}^{od},{{ud}_t}^{od}\ge 0\hspace{1em}\hspace{1em}\forall od\in K;t\in T,$$

(8)

$$y_p\in \left\{\mathrm{0,1}\right\}\hspace{1em}\hspace{1em}\forall p\in P.$$

(9)

Objective Formula (1) includes minimizing the total cost of the network (transportation costs and lost demand). Objective Formula (2) seeks to minimize the total volume of greenhouse gas emissions in the rail transport network. Constraint (3) expresses the budget limitation for the total cost of the projects. This budget covers all types of budgets needed for capacity improvement projects. Constraints (4) to (6) define the network equilibrium relationships that transfer demand from origins to destinations, considering the possibility of losing the demand in each period. Constraint (7) indicates the capacity constraint that prevents the high flow in the arcs (railway blocks) and increases capacity by carrying out projects on the arcs. The effect of increased capacity is incorporated in a given time period. Constraints (8) and (9) show the domain of the variables. The presented model is deterministic, while traffic demand is an uncertain parameter. Hence, an approach based on robust optimization is presented in the next section to address the inherent uncertainties.

4. Robust Counterpart of the Mathematical Model

Robust optimization searches for near-optimal solutions which are feasible with high probability [26]. Moreover, it proposes a linear formulation that contributes to the suggested mixed-integer linear programming (MILP) of the problem, which is usually neglected by other techniques such as stochastic programming. Hence, in this research, the model proposed by [7] is further explained for the linear optimization problem, in which the objective function is minimized and the uncertainty coefficients exist in both the objective function and the constraints.

The optimization problem is assumed as follows:

$$\begin{gathered} \mathrm{minimize} c^{T}x \\ \mathrm{subject} \mathrm{to}\dots \\ Ax\le b \\ l\le x\le u \end{gathered}$$

(10)

Moreover, the uncertainty level is defined as:

The constraint coefficient $a_{ij},j\in N=\left\{\mathrm{1,2},\dots ,n\right\}$ is modelled as an independent random variable that is symmetrical; however, ${\widehat{a}}_{ij},j\in N$ are of unknown distribution. They take value in the interval [$a_{ij}-{\widehat{a}}_{ij},a_{ij}+{\widehat{a}}_{ij}$], where ${\widehat{a}}_{ij}$ indicates the deviation from the nominal coefficient $a_{ij}$. Similarly, the objective coefficients $c_j,j\in N$ assume the value in the interval of [$c_j-d_j,c_j+d_j$], where $d_j$ is the deviation from the nominal coefficient $c_j$. Since the objective function is minimization and the purpose of robust models is to obtain the maximum regret, only one side of the proposed interval is considered. Accordingly, it is supposed that $c_j$ takes value in [$c_j,c_j+d_j$].

Now, ${\Gamma }_i$ is defined as below to formulate the robust counterpart:

Let constraint i be as $a_i^{T}x\le b_i$. $J_i$ is denoted as a set of uncertain coefficients in row i. For each row of i, a parameter of ${\Gamma }_i$ (it is not necessarily an integer number) is defined such that $[0,\left|J_i\right|]$. In fact, the role of ${\Gamma }_i$ in constraints is to justify the suggested approach’s robustness and the solution’s conservatism level. It has been proven that all coefficients are unlikely to become uncertain simultaneously [7]. Therefore, it can change to the maximum value of [${\Gamma }_i$−$⌊{\Gamma }_i⌋$]$a_{it}$. In other words, only a subset of coefficients is allowed to affect the solutions adversely. This assumption ensures that, if the same condition occurs, the optimally robust solution will be feasible. Furthermore, due to the symmetric distribution of the variables, even if the number of changing coefficients violates $⌊{\Gamma }_i⌋$, the optimal solution will remain feasible with a very high probability. Hence, ${\Gamma }_i$ is considered as the protection level for constraint i. ${\Gamma }_0$ controls the robustness level in the objective function. Therefore, it is intended to determine the optimal solutions when ${\Gamma }_0$ number of coefficients are changed in the objective function and have the greatest impact on the solution. Generally, higher values of ${\Gamma }_0$ raise the conservatism level against the higher cost it should be paid in the objective function.${\Gamma }_0$ is necessarily an integer number; however, the rest of ${\Gamma }_i$s can be an integer or non-integer. On this basis, the nominal linear robust counterpart of the problem is obtained as follows [7]:

$$\begin{gathered} {\mathrm{minimize}\hspace{1em}c}^{T}x+\underset{\left\{s_0|s_0\subseteq J_0,\left|s_0\right|\le {\mathrm{\Gamma }}_0\right\}}{\max }\left\{\sum _{j\in s_0}{d}_j\left|x_j\right|\right\} \\ \mathrm{subject} \mathrm{to} \\ \sum _j{a_{ij}x}_j+\underset{\left\{s_i\bigcup \left\{t_i\right\}|s_i\subseteq J_i,\left|s_i\right|\le ⌊{\mathrm{\Gamma }}_{\mathrm{i}}⌋,t_i\in J_i\backslash s_i\right\}}{\max }\left\{\sum _{j\in s_i}{\widehat{a}}_{ij}\left|x_j\right|+({\Gamma }_i-⌊{\Gamma }_i⌋){{\widehat{a}}_{it}}_{\mathrm{i}}\left|{x_t}_{\mathrm{i}}\right|\right\}\le b_i, \\ l\le x\le u, \end{gathered}$$

(11)

To transform the model into a linear optimization model, Lemma 1 is needed.

Lemma 1.

For the vector $x*,$ the protection level of constraint $i$ is calculated through Equation (12):

$${\beta }_i(x*,{\mathrm{\Gamma }}_i)=\underset{\left\{s_i\bigcup \left\{t_i\right\}|s_i\subseteq J_i,\left|s_i\right|\le ⌊{\mathrm{\Gamma }}_i⌋,t_i\in J_i\backslash s_i\right\}}{\max }\left\{\sum _{j\in s_i}{\widehat{a}}_{ij}\left|{x*}_j\right|+({\mathrm{\Gamma }}_i-⌊{\mathrm{\Gamma }}_i⌋){{\widehat{a}}_{it}}_{\mathrm{i}}\left|{{x*}_t}_{\mathrm{i}}\right|\right\},$$

(12)

which is equal to the optimum value of the objective function of the linear model in Equation (13).

$$\begin{gathered} {\beta }_i\left(x*,{\mathrm{\Gamma }}_i\right)=\mathrm{m}\mathrm{a}\mathrm{x}\sum _{j\in J_i}{\widehat{a}}_{ij}\left|{x*}_j\right|z_{ij} \\ \mathrm{subject} \mathrm{to} \\ \sum _{j\in J_i}{z}_{ij}\le {\mathrm{\Gamma }}_i, \\ 0 \le z_{ij} \le 1\hspace{1em}\hspace{1em}\hspace{1em} \forall i,j\in J_i. \end{gathered}$$

(13)

Lemma 1 was proved in Bertsimas and Sim (2004). By inserting the dual of Equation (13) in the robust counterpart model, it is formulated as follows:

$$\begin{array}{ll}{\mathrm{minimize}c}^{T}x+z_0{\mathrm{\Gamma }}_0+{\displaystyle {\displaystyle \sum }_{j\in J_0}}{p}_{0j}& \\ \mathrm{subject} \mathrm{to}& \\ {\displaystyle {\displaystyle \sum }_j}{a_{ij}x}_j+z_i{\mathrm{\Gamma }}_i+{\displaystyle {\displaystyle \sum }_{j\in J_i}}{p}_{ij}\le b_i& \forall i\\ z_0+p_{0j}\ge d_j{y}_j& \forall i\in J_0,\\ z_i+p_{ij}\ge {\widehat{a}}_{ij}{y}_j& \forall i\ne 0,j\in J_i,\\ p_{ij}\ge 0& \forall i,j\in J_i,\\ y_j\ge 0& \forall j,\\ z_i\ge 0& \forall i,\\ {-y}_j\le x_j\le y_j& \forall j,\\ l_j\le x_j\le u_j& \forall j.\end{array}$$

(14)

Assume that we have $x-b=A$ equation, instead of $Ax\le b$. Then, the robust model could be gained using the following inequality (Hatefi and Julai, 2014) [27]:

$$\stackrel{-}{A}-{{\mathrm{\Gamma }}_{}\widehat{A}}_{}\le x-b\le \stackrel{-}{A}+{\mathrm{\Gamma }}_{}\widehat{A}$$

(15)

Therefore, using the demand uncertainty interval $\widetilde{{D_t}^{od}}\in [\stackrel{-}{{D_t}^{od}}-\widehat{{D_t}^{od}},\stackrel{-}{{D_t}^{od}}+\widehat{{D_t}^{od}}]$, the optimal robust optimization model is presented as follows:

$$\mathrm{minimize} TC={\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{od\in K}{\displaystyle \sum }_{(i,j)\in A\cup A^{\prime }}{CO}_{ijt}{x}_{ijt}^{od}+{\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{od\in K}{{CU}_t}^{od}{{ud}_t}^{od}$$

(16)

$$\mathrm{minimize} GHG={\displaystyle \sum }_{t\in T}{\displaystyle \sum }_{od\in K}{\displaystyle \sum }_{(i,j)\in A\cup A^{\prime }}{GHG}_{ij}{x}_{ijt}^{od}$$

(17)

subject to

$$\begin{array}{l}\left(\stackrel{-}{{D_t}^{od}}-{{\mathrm{\Gamma }}_t}^{od}\widehat{{D_t}^{od}}\right)-{{ud}_t}^{od}\le {\displaystyle {\displaystyle \sum }_j}{x}_{ijt}^{od}-{\displaystyle {\displaystyle \sum }_j}{x}_{jit}^{od}\le \left(\stackrel{-}{{D_t}^{od}}+{{\mathrm{\Gamma }}_t}^{od}\widehat{{D_t}^{od}}\right)-{{ud}_t}^{od}\\ \forall i=o;od\in K;t\in T,\end{array}$$

(18)

$$\begin{array}{l}\left(\stackrel{-}{{D_t}^{od}}-{{\mathrm{\Gamma }}_t}^{od}\widehat{{D_t}^{od}}\right)\le {\displaystyle {\displaystyle \sum }_j}{x}_{jit}^{od}-{\displaystyle {\displaystyle \sum }_j}{x}_{ijt}^{od}\le \left(\stackrel{-}{{D_t}^{od}}+{{\mathrm{\Gamma }}_t}^{od}\widehat{{D_t}^{od}}\right)-{{ud}_t}^{od}\\ \forall i=d;od\in K;t\in T,\end{array}$$

(19)

Equations (3) and (6)–(9).

5. Solution Method

In this section, the methods for solving the robust mathematical model are presented. According to the bi-objectiveness of the model, an Epsilon constraint method is supposed to solve the model optimally. Moreover, a nondominated sorting genetic algorithm (NSGA-II) is employed to solve large-sized problems.

5.1. Exact Method: ε-Constraint (EC)

The ε-constraint (EC) is a popular method to resolve multi-objective problems by generating Pareto frontiers. It considers one of the objectives as the main objective of the problem and transfers the other objective functions into the constraints in each step [28]. The superiority of ε-constraint among other exact methods is its ability to control the number of generated Pareto solutions. In other words, the problem is solved by considering different values of ${\epsilon }_i$, and a Pareto solution is reported each time. For the proposed problem, the EC method is employed through Model (20):

$$\begin{array}{l}\mathrm{minimize} f1\left(x\right)\\ x\in X,\\ f2\left(x\right)\le {\epsilon }_2,\\ \dots \\ fn\left(x\right)\le {\epsilon }_n.\end{array}$$

(20)

The main steps of the EC method are: (i) Selecting an objective as the main objective function of the problem; (ii) The problem is solved in a single-objective mode for each objective function, and the optimal value of each objective function is obtained; (iii) The difference between the maximum and minimum values of the second objective function is divided into several predetermined parts. A table of values ${\epsilon }_2,\dots ,{\epsilon }_{n }$ is then generated; (iv) The problem is solved using the main objective function and ${\epsilon }_2,\dots ,{\epsilon }_{n }$; (v) Pareto solutions are reported.

5.2. Metaheuristic Approach: NSGA-II

Metaheuristic algorithms provide computational intelligence paradigms that are particularly used for complex optimization problems [29]. In the literature, NSGA-II is known as one of the most popular and widely employed algorithms used to deal with multi-objective optimization models [30]. Besides the functionality of NSGA-II, it can be considered as a standard for many other multi-objective optimization algorithms. The unique approach of NSGA-II in resolving multi-objective optimization problems has been used repeatedly to create a novel algorithm. Undoubtedly, this algorithm is one of the most basic multi-objective evolutionary optimization algorithms; therefore, it is employed in this research. To represent a viable solution of the proposed NSGA-II, a binary string is defined to show the solution representation. Figure 1 presents the specified string for the NSGA-II algorithm.

The chromosome in Figure 1 contains $p$ cells. A randomly generated number between 0 and 1, which is assigned to each cell: cells with a value below 0.5 change to 0, and those above 0.5 are rounded up to 1. The cells with a value of 1 show that the project related to that cell is conducted. By specifying the conducted projects, the arcs for each project are known. Then, each of these arcs is allocated a flow value randomly between their half of the capacity and their full capacity. Note that the arcs with the lower costs and rate of emission are of high priority when being allocated. Given the conducted projects and arc flows, the lost demand is calculated. All these values together form one initial solution for the NSGA-II algorithm. After generating the initial solutions, they may violate one or more problem constraints. In this case, the algorithm finds a neighborhood solution to satisfy the constraints. Then, a single-point crossover and mutation are implemented on the generated initial solutions to form a new population for the next generation. The pseudo-code of the proposed NSGA-II is presented as below Algorithm 1 [31].

Algorithm 1. Pseudo-code of proposed NSGA-II.

1. Initialize population 2. Generate N random solutions and insert into population 3. For (i=1 to MaxGen) do 4. Generate Childpop of size N 5. Select Parents from Population 6. Create Childern from Parents 7. Mutate Childern 8. Combine Population and ChildPop into CurrentPop whth size 2N 9. FOR each individual in CurrentPop DO 10. Assign rank based on Pareto non-dominated sort 11. End FOR 12. Generate sets of non-dominated vectors along Pareto Front 13. LOOP by adding solutions to next generation of Population starting from the best front 14. UNTIL N solutions found and determine crowding distance among the points of each front 15. End FOR 16. Print results

6. Parameter Adjustment: Taguchi Method

The parameters of the NSGA-II algorithm require an adjustment to attain higher-quality results. For this purpose, the Taguchi method is applied to obtain the optimal value of the NSGA-II parameters. Accordingly, three levels have been considered for each parameter of the NSGA-II given in

Table 1. Levels of NSGA-II parameters.

Parameter	Value in Each Level
Crossover rate	0.6, 0.7, 0.8
Mutation rate	0.1, 0.15, 0.2
Initial population	200, 300, 400

. All these different levels solved the problem, and the results were analyzed using MINITAB software. The results of the Taguchi parameter adjustment are shown in Figure 2 and Table 1 and

Table 2. Adjusted optimal values for NSGA-II parameters.

Crossover Rate	Mutation Rate	Initial Population	# of Generation
0.7	0.15	400	50

# is nomination for number.

.

7. Computational Results

In this section, the computational outputs are presented and analyzed. To do so, 15 problem instances are designed in small, medium, and large sizes, which are then solved with the suggested exact method and metaheuristic algorithm.

Table 3. Different size instances.

Problem Size	Problem No.	# of Directed Arcs	# of New Directed Arcs	# of Capacity Improvement Projects	# of Time Periods
Small	1	5	2	2	1
	2	6	3	2	1
	3	7	4	3	1
	4	8	5	4	2
	5	10	5	4	2
Medium	6	12	6	5	3
	7	13	7	5	3
	8	15	8	6	4
	9	17	10	7	4
	10	20	12	8	4
Large	11	25	15	9	5
	12	30	17	9	5
	13	50	20	10	6
	14	75	25	10	6
	15	100	30	10	6

# is nomination for number.

and

Table 4. Value of parameters.

Parameter	Distribution Function	Parameter	Distribution Function
${CAP}_{ij}$	U (1000, 10,000)	${{CU}_t}^{od}$	U (100, 300)
${CAP}_{ijt}^p$	U (100, 1000)	${D_t}^{od}$	U (10, 150)
${CC}_p$	U (1000, 5000)	${\delta }_{ij}^p$	U (0, 1)
${CO}_{ijt}$	U (200, 500)	$B$	U (1,000,000, 10,000,000)
${GHG}_{ij}$	U (10, 100)	-	-

illustrate the size of the different problems and value of the parameters, respectively.

7.1. Validation of the NSGA-II

To validate the performance of NSGA-II, the results of the exact method are compared to the algorithm. All these 15 problem instances were solved using ε-constraint (EC), and each problem’s Pareto frontier was extracted. The EC method was implemented in GAMS software using a CPLEX solver with a 3600 s time limitation for implementation. Then, the results of the proposed NSGA-II were compared to the EC, and the performance of the algorithm was evaluated. In larger problems, the EC was not capable of solving the problems optimally in the proposed time limitation. The results of Problems 1 and 6 are displayed as a problem of small- and medium-sized Pareto frontiers. The metaheuristic was coded in MATLAB software and implemented on an operating system Intel^® Core™ i7-2600 CPU @ 3.40 GHz, 4.00 GB and 64-bit. Figure 3 and Figure 4 represent the Pareto frontier of Problems 1 and 6, respectively. The figures show the conflict between the model objectives and the trade-off between them. In other words, they have an inverse relationship. The larger the first objective, the lower the second one.

According to Figure 3 and Figure 4, the NSGA-II algorithm generates Pareto frontiers with a negligible gap to the EC. In other words, it is successful in achieving near-optimal solutions and showed a reliable performance.

7.2. Measuring the Performance of NSGA-II

In this section, a number of metrics are used to evaluate the algorithm’s performance with higher accuracy. Mean of ideal distance (MID), space metric (SM), diversification metric (DM), and simple additive weighting (SAW) are defined and calculated for the Pareto frontiers obtained using EC and NSGA-II. Equations (21) to (24) present the formula for MID, SM, DM, and SAW, respectively.

$$MID=\sum _{i=1}^n\frac{\sqrt{{\left(\frac{f_{1i}-f_1^{best}}{f_{1.total}^{max}-f_{1.total}^{min}}\right)}^2+{\left(\frac{f_{2i}-f_2^{best}}{f_{2.total}^{max}-f_{2.total}^{min}}\right)}^2}}{n}$$

(21)

$$SM=\frac{\sum _{i=1}^{n-1}\left|\stackrel{-}{d}-d_i\right|}{\left(n-1\right)\stackrel{-}{d}}$$

(22)

$$DM=\sqrt{{\left(\frac{\max f_{1i}-\min f_{1i}}{f_{1.total}^{max}-f_{1.total}^{min}}\right)}^2+{\left(\frac{\max f_{2i}-\min f_{2i}}{f_{2.total}^{max}-f_{2.total}^{min}}\right)}^2}$$

(23)

$$SAW=\frac{DM+SM+\frac{1}{MID}}{3}$$

(24)

where n is the number of Pareto solutions (i = 1,2, …, n), $f_{i.total}^{max}$ and $f_{i.total}^{min}$ are the maximum and minimum value of objective i in all iterations, ($f_1^{best}$,$f_2^{best}$) is the coordination of the ideal point, $d_i$ is the Euclidean distance, and $\stackrel{-}{d}$ is the average of $d_i$. Using Equations (21) to (24), the value of the metrics was calculated in each of the 15 problems and is presented in

Table 5. The values of metrics for EC and NSGA-II.

	DM		MID		SM
Problem	EC	NSGA-II	EC	NSGA-II	EC	NSGA-II
1	1.84	1.75	0.84	0.9	0.33	0.39
2	1.52	1.4	1.17	1.19	0.69	0.8
3	1.83	1.75	1.08	1.17	0.49	0.64
4	1.32	1.16	1.18	1.25	0.32	0.46
5	0.98	0.77	0.95	1.03	1.51	1.69
6	1.79	1.64	0.99	1.01	0.84	0.9
7	1.23	0.94	1.25	1.28	1.79	1.85
8	0.99	0.86	1.04	1.15	0.91	0.99
9	1.74	1.63	0.63	0.66	0.62	0.81
10	1.52	1.02	1.52	1.63	0.9	0.94
11	1.44	1.72	1.45	1.4	1.4	1.39
12	0.69	0.73	1.56	1.48	0.88	0.81
13	0.85	1.06	1.44	1.11	1.66	1.42
14	0.98	1.57	1.98	1.57	1.49	1.4
15	1.39	1.48	0.88	0.79	1.83	1.81
Average	1.34	1.30	1.20	1.17	1.04	1.09

. In Figure 5, the values of the metrics for small- and medium-sized problems have been compared for both approaches. Accordingly, NSGA-II is rigorous in finding Pareto frontiers close to the optimal frontiers using the EC method. Moreover, for a general analysis of the NSGA-II performance, a comparison of the SAW measure has been presented in Figure 6. This figure also proves the closeness of the Pareto solutions using NSGA-II to the optimal Pareto solutions. As a result, the proposed NSGA-II can be reliable in solving large-sized problems that the EC is not able to solve optimally in the specific time limit.

Notably, EC is better than NSGA-II in some cases because it is not capable of solving the problem optimally within the proposed time limit. For instance, in terms of the MID metric, the suggested NSGA-II algorithm shows a better performance than EC when the size of the problem increases. On such an occasion, the proposed NSGA-II can be highly helpful and reliable. Moreover,

Table 6. Computational time (s).

Problem	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
EC	0.2	4.3	13.1	17.2	66.9	112.4	549.7	1233	2416	3491	-	-	-	-	-
NSGA-II	13.5	29.4	44.1	60	94.5	112.9	163.2	169.2	181.7	221.5	245.7	266.3	297.6	310.7	350.1

and Figure 7 illustrate the computational time of the methods for each problem. Evidently, by increasing the problem size, the computational time of the EC method increases exponentially, until it cannot solve the large-sized problem in 3600 s. However, the proposed NSGA-II is capable of solving large-sized problems in a reasonable time.

8. Conclusions, Managerial Insights, and Future Works

This paper studied a multi-period railway network design by developing a novel robust bi-objective mixed-integer linear programming (MILP) model. Various economic aspects related to this problem, including direct/indirect costs and operating and construction costs were considered in the model. Moreover, the model was able to lessen the environmental impact. To overcome the demand uncertainty, a robust optimization approach based on interval-based uncertainty was employed. Then, the robust counterpart model of the problem was developed and discussed. The model was solved optimally using the ε-constraint (EC) method. Due to the complexity of the robust model, an efficient metaheuristic, namely, the nondominated sorting genetic algorithm II (NSGA-II), was utilized to deal with large-sized problems. Computational results proved the efficiency of NSGA-II in small- and medium-sized instances as it could find near-optimal Pareto solutions. As a result, the proposed NSGA-II can be utilized for solving large-sized samples in a very short time.

This study could also assist managers and decisionmakers in this area. First, by generating the Pareto frontiers, we illustrated the conflict between cost and emission. This means that the higher amount of cost leads to a lower amount of emission. Here, depending on the viewpoint of the managers, each of the Pareto solutions can be selected as the solution to this problem. If the manager insists on lower expenses, a Pareto solution from the left-hand side of the chart is selected. In the case of being an environmentally friendly manager, a Pareto solution from right-hand side makes more sense. Therefore, it has facilitated the decision making for the managers according to their point of view and policy. Second, the discussion of whether to strengthen the existing facilities or establishing a new one has always been controversial. Based on real data, this study can make such a decision. For instance, a capacity increase is suggested by the model when there are numerous lanes available but cannot fully meet the demand. On the other hand, when there is an insufficient number of facilities compared to the demand, establishing a new lane is recommended.

According to the main limitations of the study, the following suggestions can be considered for future research directions:

• Developing other uncertainty handling techniques, such as fuzzy programming [32], to be compared to the proposed robust optimization;
• Applying other well-known algorithms, such as the Multi-Objective Grey Wolf Optimization (MOGWO) algorithm [30] and Multi-objective Particle Swarm Optimization (MOPSO) to be compared to the proposed NSGA-II;
• Implementing neural networks, machine learning, and deep learning tools to deal with effective demand forecasts [33,34,35];
• A real case study can be conducted to show the applicability of this model in the real world and provide deeper insights to this problem.