A Time-Dependent Fuzzy Programming Approach for the Green Multimodal Routing Problem with Rail Service Capacity Uncertainty and Road Traffic Congestion

School of Management Science and Engineering, Shandong University of Finance and Economics, No. 7366, Second Ring East Road, Jinan, Shandong Province 250014, China Institute for ProductionManagement, WU Vienna University of Economics and Business, Welthandelsplatz 1, 1020 Vienna, Austria Unit of Logistics and Informatics, KTH Royal Institute of Technology, Tekniringen 10, 10044 Stockholm, Sweden School of Traffic and Transportation, Beijing Jiaotong University, Beijing 100044, China


Introduction
With the rapid development of globalization, companies are seeking for specialized partners all over the world to outsource their businesses and also for new markets to make more profit.International trade and global commodity circulation thus get significantly motivated during the last decades [1].Supported by the transportation industry, this prosperity results in a remarkable expansion of the transportation network from a limited region to the entire world.The expanded transportation network extends the distribution channels and enhances the difficulty of transportation organization from the viewpoint of economy, efficiency, reliability, and environmental protection.The traditional unimodal transportation is no longer suitable for this situation.Meanwhile, an advanced transportation mode, namely, the multimodal transportation, emerges in the transportation industry.It combines different transportation modes to generate origin-to-destination routes and integrates their respective advantages in the transportation process.Furthermore, by using standardized containers, multimodal transportation can be highly mechanized and the efficiency of the operations at terminals is greatly improved [2].Many empirical case studies, for example, Bookbinder and Fox [3], Janic [4], and Liao et al. [5], have demonstrated the superiority of the multimodal transportation in the long-haul transportation setting.Nowadays, multimodal transportation is getting more and more popular and plays a crucial role in international trade.Since the transportation chains in this area are usually long combining different actors, effective solutions for their coordination are necessary in order to ensure that the goods are delivered according to the customer's requirements.In this context, the multimodal routing problem has gained interest of both researchers and practitioners [6,7].This operational problem aims at selecting the best routes to move containers from their origins to destinations through the multimodal service network according to customer demands [2].In order to achieve this objective, a number of different factors have to be considered.
First of all, multimodal transportation combines different modes with different operational constraints that have to be represented in the mathematical model.In general, it can be distinguished between time-flexible modes (e.g., road) that can be used whenever they are needed and schedule-based modes (e.g., rail) that operate according to fixed schedules planned in advance.These schedules regulate the route as well as the (un-)loading, departure, and arrival times of the transportation service.Since these schedules cannot be easily changed, they should be considered when planning the routes of the containers and their transshipment between different transportation services; otherwise, the multimodal plans might be infeasible in reality.
Secondly, the actors and stakeholders involved in multimodal transportation systems have usually different interests which are often in conflict.Whereas transportation operators usually tend to minimize only the transportation costs, customers are also interested in the reliability of the routes in order to avoid situations where a plan becomes infeasible due to late arrival of a service or insufficient capacity.Since the multimodal routing is done before the start of transportation, when the information about the transportation services (e.g., traffic situation, travel time, and capacity) is not completely known, possible sources of uncertainty should be considered already during the planning [8].As an example, the scheduled rail services have usually limited capacity due to the limited length of the tracks in the railway stations and limited locomotive power [9].However, the available capacities of rail services are uncertain since they are also dependent on other transportation plans that are difficult to predict [10].If only deterministic capacities are considered in the model, two consequences may happen in practice: the planned routes are either infeasible if the capacities are estimated to be too large, or the planned routes are not the actual best ones if the capacity estimation is too small.Therefore, capacity uncertainty is a factor that should be considered in the decision-making process.In addition to that, road services frequently suffer from various disruptions, such as natural disasters (e.g., bad weather) and manmade disruptions (e.g., traffic accidents and congestion) [11].Among these disruptions, traffic congestion is very common due to the rapidly increasing traffic volume in recent years [12].Traffic congestion results in the travel time uncertainty and delayed arrivals of the containers at the nodes, which disrupts the transshipment if the containers have to take another scheduled transportation service at the same node and also possibly violates the due date constraint at the destination.Thus, congestion is another factor that should be also considered in the decision-making process.
Thirdly, the environmental aspects of transportation have drawn considerable attention from the public and the government [13].Among the various environmental issues, global warming caused by greenhouse gas emissions has been considered as a tremendous threat to human sustainable development.CO 2 accounts for approximately 80% of the entire greenhouse gas emissions [5].During recent decades, large numbers of countries have proposed regulations to ease the global warming by reducing CO 2 emissions.For example, in 2012, the National Development and Reform Commission of China proposed a pilot project on trading carbon emissions [2].Since the transportation sector is one of the biggest contributors to the CO 2 emissions [14,15], the transportation industry should accept this responsibility and consider controlling and optimizing the CO 2 emissions during the transportation planning stage to promote the development of an environmentally friendly transportation.In order to respond to these environmental requirements, in this study, we extend the multimodal routing problem into a green version by considering the CO 2 emissions.
The inclusion of uncertain factors and the combination of multiple objectives (i.e., economic and environmental) in multimodal transport planning is an emerging field which has been studied in recent years [13][14][15].However, the authors usually consider only one source of uncertainty or model the CO 2 emissions in form of charges which limits their influence on the resulting routes.Therefore, we present a multimodal routing planning approach which deals with the following question: How can the sustainability and reliability of multimodal routing be improved by considering environmental factors and uncertainty?
The contributions of the paper are fivefold: (1) A mixed integer nonlinear optimization model for multimodal routing including multiple objectives (i.e., economic and environmental) and multiple uncertainty factors (i.e., road traffic congestion and rail service capacity) is defined.
(2) A fuzzy logic is proposed to model the capacity uncertainty, and a piecewise linear function is used to represent the road traffic congestion.
(3) The influence of CO 2 emissions on the routes is analyzed comparing the CO 2 emissions charging method and a bi-objective approach in which minimizing CO 2 emissions is set as an independent objective.
(4) Since the proposed optimization model is nonlinear, a linearization approach is proposed in order to reduce the complexity of the model.
(5) The planning approach is applied to a real-world case study from China with specific network design and parameters.In this case, the network is limited to road and rail services; however, other modes with fixed schedules (e.g., inland waterway and deep sea 2 Complexity shipping) could be also treated in a similar way as it is the case of rail.
The paper is organized as follows: In Section 2, we review the current literature on multimodal transportation routing problem.In Section 3, we present a detailed modelling methodology.The improvements of our settings compared to the existing multimodal routing literature are indicated by comparing the two sections.In Section 4, based on the proposed modelling framework, we establish a fuzzy chanceconstrained mixed-integer nonlinear programming model to mathematically describe the multimodal routing problem explored in this study.Since the proposed model is nonlinear, in Section 5 we design a linearization-based exact solution strategy to solve the problem.In Section 6, an empirical case study based on the Chinese scenario is given to demonstrate the feasibility of the methods in dealing with the practical problem.Also in the case study section, sensitivity analysis, bi-objective optimization scenario analysis on CO 2 emissions, and fuzzy simulation are adopted to draw conclusions which can further help the decision makers to better organize the multimodal transportation.Finally, conclusions of this study are drawn in Section 7.

Literature Review
The multimodal routing problem is a combinatorial optimization problem which deals with finding the optimal routes for customer demands respecting the constraints of the network [16].Customer demands are represented by transportation orders characterized by their origin, destination, volume of containers, release time at origin, and due date at destination.Whereas the release time is usually considered as fixed, the due date can be either considered as a hard constraint (see, e.g., Verma et al. [17], Sun and Lang [2,18], and Sun et al. [19]) or delayed deliveries can be accepted causing additional penalty costs (see, e.g., Demir et al. [20], Hrušovský et al. [14], and Zhang et al. [21]. However, considering only late arrivals might not be sufficient since early deliveries can also cause problems with the storage of the goods until they are needed.Therefore, various studies consider due date time windows instead of time points [8,22], which is also the case in this study.Moreover, modelling a soft time window is more suitable than a hard one, because customers allow the violation to the time window (i.e., early and delayed deliveries) within a certain degree.
The containers from transportation orders need to be transported in the multimodal service network by various transportation services with their specific characteristics.One of these characteristics is the usually fixed schedule of the high-capacity and environmentally friendly transportation modes (e.g., railway or inland waterway) in contrast to the time-flexible (nonscheduled) road service.However, the majority of the related articles neglects the schedules and considers rail services to be also time-flexible.This facilitates the modelling since only transshipment times between services in terminals need to be considered (Sun et al. [23], Jiang and Lu [24], Cai et al. [25], Zhang et al. [21], and Xiong and Wang [26]).However, this does not correspond to reality where railway services are usually planned well in advance [2,19].Therefore, the schedules are usually fixed in the tactical network design phase (see, e.g., Wieberneit [27], Andersen et al. [28], and SteadieSeifi et al. [29] for this area) and have to be accepted in the operative multimodal routing considered in this paper (Liu et al. [30], Lin [31], Demir et al. [20], and Hrušovský et al. [14]).
The routes can be optimized according to different criteria.Traditionally, the optimization models concentrate on costs or time [32].However, various surveys show that multiple criteria, among others also sustainability and reliability, are important for transportation planners [33,34].In case of sustainability, vast majority of available literature only considers economic and environmental criteria whereas social aspects of transportation are neglected [35,36].As a result, the comparison of transportation costs and CO 2 emissions is decisive for finding the optimal routes.
The CO 2 emissions can be either integrated into the cost objective function as emission charges or modelled as an independent objective in multiobjective optimization.The majority of the articles on optimizing CO 2 emissions in transportation adopts the first approach (e.g., Chang et al. [37], Zhang et al. [21], Wen et al. [38], Sun and Lang [2], Chen et al. [39], Demir et al. [20], and Hrušovský et al. [14]), whereas only a few studies used the second approach (e.g., Sun and Chen [40] and Qu et al. [41]).However, it is difficult to decide which method is better to use since the charging method is easier to implement but might not always be applicable, especially in cases where the emission charges are too low in comparison to transportation costs and therefore do not have any influence on the resulting routes [22].
Another important decision is the choice of the emission calculation method since most of the models require a lot of parameters (e.g., fuel type, vehicle and route characteristics, and traffic condition) that cannot be always estimated in advance [20,[42][43][44][45].As an alternative, the activity-based method that multiplies activity intensity (unit: TEU-km) with the CO 2 emission factor (unit: g/TEU-km) to calculate the total CO 2 emissions [5] yields a better feasibility in transportation practice and has been already applied in various studies [2,5,37,39,41,46].
The reliability of transportation is influenced by various sources of uncertainty (e.g., weather, accidents, congestion, and technical issues) that can cause disruptions in transportation processes.Despite this fact, most of the existing models consider all factors as deterministic and do not take into account uncertainty [47].In order to improve the reliability of transportation plans, the uncertainty should be considered during the planning phase.Despite the variety of possible reasons for disruption, their effects influence either the demand, the travel time, or the supply (capacity) of transportation.Whereas extensive literature about demand uncertainty in multimodal routing is available (see, e.g., Lium et al. [48] and Ho et al. [49]), the consideration of travel time and capacity uncertainty is rather limited [50,51].Moreover, usually the models deal with only one source of uncertainty.However, the consideration of multiple uncertainties at the same time might lead to improved routes and therefore in 3 Complexity our study the travel time uncertainty for road services and capacity uncertainty for rail services should be combined.
Road traffic congestion leads to travel time uncertainty of road services.Only a few studies, for example, Demir et al. [20], Sun and Chen [40], and Uddin and Huynh [11], adopted stochastic programming to formulate the multimodal service network problem under disruptions.However, in practical planning, due to the lack of reliable historical data [52], it is quite difficult to fit the probability distributions to the travel times of all road services in the multimodal service network.Such disadvantage of stochastic programming in addressing the uncertainty is emphasized by Zarandi et al. [52] in a location-routing problem.Furthermore, the fuzzy programming that uses fuzzy variables (usually triangular or trapezoidal numbers) to evaluate the uncertainty fails to represent the variation of the travel time of a road service en route with respect to the time of the day [53].Therefore, it is worthwhile to find an applicable approach (e.g., time-dependent travel time that is widely discussed in the vehicle routing problems [53,54]) to formulate the road traffic congestion and introduce it to the multimodal routing modelling.
The capacity uncertainty has only been considered by a few studies in different transportation planning problems, for example, multiobjective solid minimal cost flow problem by Wang [55], multiobjective multi-item solid transportation problem by Kaur and Kumar [56], and resource allocation problem for containerized cargo transportation networks by Kundu et al. [57].Whereas Wang [55] used stochastic programming to address this issue, Kaur and Kumar [56] and Kundu et al. [57] formulated the capacity uncertainty from a fuzzy programming viewpoint.The fuzzy programming approach has been also used in the context of vehicle routing problem [52,58] and shortest path problem [59], since it can avoid the infeasibility of the stochastic programming in case of lack of historical data [58].However, to the best of our knowledge, there exists no study which would consider capacity uncertainty in multimodal routing problem.
The considered aspects of the multimodal routing approach increase the complexity of the problem.The modelling of constraints for schedule-based services leads to the nonlinearity of the model, which can be seen in the models proposed by Liu et al. [30], Lin [31], and Chang [60].Moreover, the modelling of uncertainty sources might also lead to nonlinearity of the model.As for the problem formulated by a nonlinear model, Xie et al. [61] pointed out that the nonlinearity makes the problem very difficult to solve, while using linearization technique to get the equivalent linear model can enable the problem to be relatively easily solved.The linearization technique should be designed according to the specific nonlinearity of the nonlinear model.Its main idea is to introduce new linear forms to replace the nonlinear components in the model equivalently.
As the literature review shows, the interest of researchers in multimodal routing has been increasing during the last years, considering different aspects of the studied planning problem.However, there still exists potential for improvement especially in the area of sustainability and reliability of transportation.Therefore, our paper contributes to these aspects of the multimodal routing since it introduces and compares two different methods for modelling CO 2 emissions.Moreover, two different sources of uncertainty (travel time and capacity) are considered within the planning process which should reduce the influence of potential disruptions on created transportation plans and improve the reliability.In case of rail service capacity, the fuzzy programming approach is used in order to better represent the possible capacity variations.Although the application of the model to the Chinese case study requires some specific assumptions, the general model can be used for any multimodal service network.
In order to give an overview of the current state of the literature, the most important articles with the covered aspects are shown in Table 1.

Modelling Methodology
As already mentioned, the objective of the multimodal routing problem in this study is to select the best green routes to move containers from multiple transportation orders through the multimodal service network with road traffic congestion and rail service capacity uncertainty.The multimodal service network is composed of nodes and arcs, where nodes represent the terminals serving as origins and destinations of transportation orders as well as transshipment points between rail and road services.The arcs represent the services connecting the nodes.Related information on road and rail services, such as schedules, transportation distance, and capacity, is regarded as the parameters of the nodes and arcs.
In the studied network, rail services are the backbone for multimodal container transportation, while road services play an important role in origin pickup and destination delivery.In some cases, road services also take direct origin-todestination transportation.Figure 1 shows the basic topological structure of the multimodal service network.Overall, the network structure is similar to a hub-and-spoke network with multiple hubs where origins and destinations represent spokes while railway container stations represent hubs.

Modelling Transportation Services.
Service characteristics are dependent on the chosen transportation mode.In case of road service, the trucks are assumed to be fully flexible which means that they are available every time they are needed and it is possible to assign as many trucks as necessary to carry the containers.As a result, the road services are uncapacitated and their departure time can be adjusted according to the customer requirements and the current state of the network.This means that the departure time can be either immediately after the containers are loaded on the vehicle, but can be also postponed in order to avoid traffic congestion.However, it has to be ensured that the truck arrives to its destination before the loading cutoff time of the next planned service in order to avoid infeasibility of the transportation plan.Besides that, too early arrivals should be also avoided as inventory costs can be charged for containers waiting for the next service at a terminal.Therefore, it is important to 4 Complexity determine the planned departure time of a truck from the current node.
In contrast to that, rail services are modelled as block container trains driving between two nodes according to fixed schedules with the following characteristics: (1) Operation time window at the loading/unloading organization station that regulates the operation start time and cutoff time regarding loading/unloading containers on/off the train (2) Scheduled departure/arrival time from/at the loading/unloading organization station (3) Operation period.For the convenience of modelling, same rail services in different operation periods are considered as different services.
Since the loading/unloading operations of both road and rail services are highly mechanized, the time needed for loading/unloading is relatively low.Consequently, we assume that these times can be neglected in the multimodal routing model.Besides that, the capacities of rail services are limited.Taking into account these assumptions, the origin-to-destination transportation process combining road and rail services is depicted in Figure 2 and can be described as follows: (1) The containers will first depart from the origin (node 1) at the planned departure time that is not earlier than the earliest release time defined by the customer, and then arrive at the loading organization station (node 2).
(2) In order to be able to load the containers on the train operated from node 2 to node 3, their arrival time at node 2 should be no later than the operation start time.If the arrival time is earlier than the operation start time, there will be an inventory period from the arrival time to the operation start time at node 2, and the containers should wait until the operation start time and then get loaded.Moreover, there exists a free-of-charge period for inventory costs [64].Only when the inventory period is longer than this period will inventory costs be created.
(3) When loaded on the train, the containers will wait until the scheduled departure time and then leave node 2 along with the train.
(4) The containers will arrive at the unloading organization station (node 3) at the scheduled arrival time of the train.However, they should wait until the operation start time and then get unloaded.Therefore,  (5) After unloading from the train, the transportation organizer can select a planned departure time that is no earlier than the operation start time to move containers to their destination (node 4).In case of an early arrival to the destination, the containers have to be stored again until the planned delivery date in a warehouse provided by the third-party inventory companies that can have different conditions in comparison to warehouses used in intermodal terminals [65].

Modelling Uncertainty.
Considering the uncertainty factors, traffic congestion is an important factor in road transportation that influences the travel times.Congestion can occur either regularly due to insufficient road capacity or unexpectedly due to, for example, an accident.Since the latter one is hard to predict, only regularly occurring congestion is considered in this study.However, this congestion can vary with the time of the day or the section considered, therefore piecewise linear functions can be used to model such situations [53,54].By combining limited historical data and experiences from the experts, such travel time functions can be established.Figure 3 shows an example of the timedependent travel time function.There are two peak hours from 9 am to 10 am and from 6 pm to 8 pm on the arc served by road services.
Besides the road congestion, the available capacity of a rail service can be also uncertain due to various reasons (e.g., order cancellation, vehicle breakdown, and lack of available vehicles).Since the reasons are manifold and the historical data about capacity variation might not be available, using the stochastic programming for capacity uncertainty might not be feasible.As a consequence, capacity uncertainty can be effectively expressed in a fuzzy way by using the triangular fuzzy variable represented by the minimum of the possible capacities, most possible capacity and maximum of the possible capacities.Similar to the traffic congestion issue, fuzzy capacity can be attained by limited historical data and expert experiences.

Modelling and
Optimizing the CO 2 Emissions.As stated in Section 2, there are two approaches to optimize the CO 2 emissions, including the CO 2 emission charging method and the bi-objective optimization method.The optimization model is firstly constructed in Section 4 using the charging method, and the bi-objective method is then introduced in Section 6 where both methods are compared.For calculating the CO 2 emissions, the activity-based method is used.This method calculates the total CO 2 emissions for a transportation service using the formula (EM • Q • L) where EM is the CO 2 emission factor of the service (unit: g/TEU-km), Q is the number of the containers carried by the transportation service (unit: TEU), and L is the corresponding transportation distance (unit: km).

3.4.
Modelling the Customer Demands.Transportation orders are considered to be unsplittable which means that all containers of an order have to use the same route.The release times of orders are fixed whereas the delivery can be done within a due date time window limited by a lower and upper bound.For each transportation order, it is best if the arrival time of the containers falls within the time window.In this case, there will be no extra costs at the destination.However, if the containers arrive at the destination earlier than the lower bound, inventory costs will be created, because the goods need to be stored before being processed at the destination.If the arrival time is later than the upper bound, penalty costs should be charged for the delayed delivery.

Mathematical Model for the Multimodal
Routing Problem

4.2.
A Fuzzy Chance-Constrained MINLP Model.The framework of the mathematical model is derived from the classical arc-node-based formulation that has been widely used in the routing problems of various research fields, for example, telecommunication systems, manufacturing systems, and internet service systems [1].The mathematical model for the multimodal routing problem explored in this study is presented as follows.The proposed model is a combination of the fuzzy chance-constrained programming and the mixed integer nonlinear programming (MINLP).
Objective function: Eqs. ( 1), ( 2), ( 3), (4), and ( 5) are successively the transportation costs en route, loading and unloading costs at the nodes, inventory costs, penalty costs caused by delayed deliveries at destinations, and CO 2 emission costs.Their linear summation represents the generalized costs.The objective function minimizes the generalized costs for multiple transportation orders using multimodal transportation.This setting reflects the objective of satisfying customer demands and minimizing monetary expenditures.Therefore, the multimodal routing in this study is a kind of serviceoriented planning.Subject to Eq. ( 6) is the flow conservation constraint which ensures the continuity of the nodes and arcs on the planned routes, which is a common constraint in the node-arc-based formulation [1].Eq. ( 7) means that no more than one transportation service on one arc can be utilized to serve a transportation order.The combination of the two equations above ensures that containers in a transportation order are unsplittable, which is claimed in Section 3.4.
Eq. ( 8) is the fuzzy chance constraint of the rail service capacity based on the fuzzy credibility measure [66].In this equation, Cr • represents the credibility of the event in • .Eq. ( 8) ensures that the credibility of the event that the entire loaded container volume on a rail service does not exceed its available capacity should not be lower than a predetermined confidence α ∈ 0, 1 .The value of α indicates the decision makers' preference for addressing the corresponding decision-making.
Eq. ( 9) assumes that the arrival time of the containers of transportation order k at the origin equals its earliest release time.
Eqs. ( 10), ( 11), (12), and (13) compute the arrival times of the containers of transportation order k at the nodes on the route by road services.First, we determine the planned departure time of the containers at the current node i by (10).Then, we normalize such departure time into range [0, 24] by using the floor integral function (11) and normalization formula z k i − 24 • n k i , so that we can further generate the corresponding travel time t k ijs of the containers on arc i, j by road service s according to the piecewise linear function (12).Finally, we can gain the arrival time of container flow k at the following node j by (13).
Eq. ( 14) computes the arrival times of the containers of transportation order k at the nodes on the route by rail services.As claimed in Section 3.1, l r i is the effective arrival time of containers at node i.Because we assume that 8 Complexity loading/unloading time can be neglected, y k i is also the earliest departure time of containers at transshipping node.
Eq. ( 15) is the rail service operation cutoff time constraint.It ensures that the arrival time of the containers of transportation order k at node i should not exceed the operation cutoff time of rail service s at the same node, if this rail service has to be used in order to move these containers on arc i, j .Directed arc from node i to node j, and i, j ∈ A.

Sets Representations Γ ij
Set of the rail services on arc i, j .

Ω ij
Set of the road services on arc i, j .

S ij
Set of the transportation services on arc i, j , δ − i Set of the predecessor nodes to node i, and Set of the successor nodes to node i, and δ + i ⊆ N.

Variables Representations w k ijs
Charged inventory time for transportation order k at node i before its containers are moved on arc i, j by service s, unit: h.Unit CO 2 emission costs, unit: ¥/g.em s CO 2 emission factor for service s, unit: g/TEU-km.τ s i Free-of-charge inventory period of service s at node i, unit: h.

M
A significant large positive number.α Confidence in the fuzzy chance constraint.9 Complexity Eqs. ( 16) and ( 17) compute the charged inventory time of the containers of transportation order k at node i before being moved on arc i, j by a rail service and a road service, respectively.In (17), especially if i = o k and s ∈ Ω ij , τ s i = 0.
The above mathematical model for the multimodal routing problem is summarized in Figure 4. Figure 4 classifies all equations (excluding variable domain constraints) into three categories, including "Modelling the network" which realizes the methodology proposed in Section 3, "Setting the objective" which reflects that satisfying customer demands is the core of the routing problem, and "Basic assumptions" which are the basic conditions of the model.This should help the readers to better understand the functions and roles of various equations in the mathematical model.
By using the solution strategy designed in the next section, the proposed model can be used to optimize any green multimodal routing problems that are composed of scheduled rail services with uncertain capacities and time-flexible road services with uncertain travel times.It can plan the best multimodal routes to accomplish multiple transportation orders with soft due date time windows.

An Exact Solution Strategy
5.1.Determining the Solution Strategy.The proposed mathematical model in Section 4.2 is easily understood, because it is established straightforwardly according to the settings of the specific multimodal routing problem.However, the mathematical model contains three categories of nonlinearity, including (1) fuzzy chance constraint as (8), (2) piecewise linear constraint as (12), and (3) general nonlinear components as (3), ( 4), ( 13), ( 14), (16), and (17).Such nonlinearity increases the difficulty of the problem to be effectively solved by the exact solution algorithms (e.g., branch-and-bound algorithm and simplex algorithm).Especially for the large-scale real-world problem, due to the nonlinearity, the solution to the problem will usually get into the local optimum, and massive computational time will be consumed to generate the solution.Therefore, if we can eliminate the nonlinearity of the model by designing the equivalent linear reformulations to the nonlinear equations in the model, then the problem can be solved with the help of the mathematical programming software (e.g., LINGO, CPLEX, and GAMS) where exact solution algorithms can be implemented efficiently.
We devote this study to developing a linearization-based exact solution strategy to address the multimodal routing problem.Although it is widely acknowledged that heuristic algorithms have better performance than the exact solution algorithms in dealing with the multimodal routing, which is known as an NP-hard problem, exact solution algorithms have two significant advantages.The first one is that they can test if the model observes the mathematical logic.The second is that they can provide an exact benchmark to verify the efficiency and accuracy of the heuristic algorithms [67].Fazayeli et al. [68] also highlighted that the development of an exact algorithm and the comparison between the algorithms will be their future work for the location-routing problem on a multimodal transportation network.

Model Linearization
According to (24), we can gain the crisp equivalency to the fuzzy chance constraint (8) by replacing a with ∑ k∈K q k • x k ijs and b with v ijs .However, the crisp equivalency is a piecewise linear function which is still essentially a nonlinear formula.Thus, we conduct following modifications shown as (25) and (26) to transform it into a pure linear function.After such modification, during the simulation process, we can use (25) to substitute (8) if α ∈ 0 5, 1 , (26) to substitute the same equation otherwise.(27).
To linearize the piecewise linear function (27), the following procedures are conducted.

Optimization principle
Minimizing the total generalized costs for multiple orders Eqs.1-5

Customer demand orientation
Optimization object

Mathematical model
Containers in a transportation order are unsplittable.Eqs.6 and 7 The arrival time of the containers of a transportation order at the origin equals its earliest release time.Eq.9  12) is shown as (31).

Linearization of Other Nonlinear Components
(1) Linear Reformulation 1. Eq. ( 3) with nonlinear component ∑ k∈K c k dest • q k • max t L k − y d k , 0 can be linearized by ( 32), (33), and (34) where g k is a newly introduced variable.
If the containers arrive at the destination earlier than (32) and (33), 32) and ( 33).If the containers arrive at the destination later than t According to (32) and (33), g k ≥ 0. Minimization of (32) will then lead to g k = 0.In this case, ∑ k∈K c k dest • q k • max t L k − y d k , 0 is also equivalent to (32) and (33).
Above all, the equivalency above is verified.Similarly, we have the following linear reformulation.
We take z k i = 60 as example to verify the equivalency above.According to (12), n k i = 60/24 = 2 5 = 2.Meanwhile, according to ( 45) and ( 46), there exists Because n k i is a nonnegative integer, finally n k i = 2. Consequently, the equivalency above is tenable.
For detailed proofs of Linear Reformations 4 and 5, readers can refer to our previous study [2].

Improved Linear Reformulation.
For the sake of readability, the improved linear reformulation derived from the mathematical model in Section 4.2 is summarized as follows.

Empirical Case Study Based on the Chinese Scenario
6.1.Case Description.In this section, we present a real-world case to demonstrate the feasibility of the proposed model and exact solution strategy in optimizing the practical problem.This empirical case is oriented on a Chinese scenario.In this scenario, several consignments of containers need to be transported from the western inland cities (e.g., Lanzhou and Hohhot) to the eastern sea ports (e.g., Qingdao) through the road-rail multimodal service network.
In the multimodal service network, "road services" refer to the container trucks and "rail services" refer to the block container trains.The topological structure of the multimodal service network is shown in Figure 5.In the multimodal service network, there are 10 periodically operated rail services (block container trains) and 15 road services (container trucks).The schedules of the rail services and the timedependent travel times of the road services are separately presented in Appendix A (available here) and Appendix B (available here) in the supplementary material.Note that all the values in Appendix A (available here) and Appendix B (available here) are expressed as real numbers; for example, 10 : 30 is rewritten as 10.5 and 5 : 00 at the second day of the planning horizon is 29 (24+5).Other times in different days can be expressed in the same way.
There are 10 transportation orders within a short-term planning horizon.The containers in these orders need to be transported from Lanzhou and Hohhot (inland cities in Northwest China) to Jiaozhou and Huangdaogang (sea ports in Qingdao, Shandong Province) and Lianyungang (a sea port in Jiangsu Province) through the multimodal service network shown in Figure 5.The detailed information on the multiple transportation orders is given in Appendix C (available here) in the supplementary material.Note that in Appendix C (available here), if the destination inventory costs equal zero, it means that the receivers have their own inventory facilities and there is no need to rent them from other warehouse companies.
6.2.Parameter Setting.The unit transportation costs en route (unit: ¥/TEU) regarding the two transportation services are determined by (48).
As can be seen in ( 48), the unit transportation costs for rail services en route include two parameters: c 1 rail which is the costs related to the container volume (unit: ¥/TEU) and c 2 rail which is the costs related to the turnover of the containers (the multiplication of the container volume and corresponding transportation distance, unit: ¥/TEU-km).For road services, such costs only relate to one parameter c 2 road which is similar to c 2 rail (unit: ¥/TEU-km).The values of all these parameters above are regulated by the Ministry of Transport of China and China Railway Corporation [70].20 ft ISO containers are utilized to contain goods in the empirical case.According to 2015 China Railway Statistical Bulletin [71], the national railways consumed 15.69 million tons of the standard coal and accomplished 33,503.67hundred million ton-kilometers of freight transportation.The unit CO 2 emission rate of the standard coal is 0.69, and approximately 92% carbon is conversed into CO 2 .Assuming that the 20 ISO containers are all fully loaded (each container weights 24 tons), we can then calculate the CO 2 emission factor of rail service as the following formula [72] where 44 is the molecular weight of one CO 2 molecule and 12 is the atomic weight of one C atom.Similarly, by using the data from 2015 China Transportation Industry Development Statistical Bulletin [73], we can also gain the CO 2 emission factor of road services which is 1064 g/TEU-km.The values of all the parameters in the empirical case are summarized as shown in Table 3.In addition, the values of the unit inventory costs of transportation services at the nodes (c s i ) and the free-of-charge inventory periods of services at the nodes (τ s i ) are all presented in Figure 5.When the nodes are specifically the railway stations, all the corresponding c s i = 3 125 ¥/TEU and τ s i = 48 h.When the nodes are the origins of the transportation orders, c s i varies based on the locations of the nodes and 6.3.Simulation Environment.Described by the linearized programming model shown in Section 5.3, the empirical case can be solved by any exact solution algorithms, for example, branch-and-bound algorithm and simplex algorithm.In this study, we adopt the branch-and-bound algorithm to solve the problem.The algorithm is conducted by the mathematical programming software LINGO 12 designed by LINDO Systems Inc., Chicago, IL, USA [75], on a ThinkPad Laptop with Intel Core i5-5200U 2.20 GHz CPU 8 GB RAM.The scale of the empirical case is indicated by Table 4.

Optimization Result and Multimodal Route Illustration.
The optimization result of the empirical case and the performance of the exact solution strategy are shown in Table 5.
The structure of the minimal generalized costs for the multiple transportation orders is illustrated by Figure 6 where C1 to C5 are successively the transportation costs en route, loading/unloading operation costs, inventory costs, penalty costs, and CO 2 emission costs, which also correspond to (1) to (5) in the objective function.
The best routes for the containers in the multiple transportation orders are given in Table 6.According to the planned routes, the accomplishments of 4 transportation orders satisfy their due date time windows exactly.6.5.Case Discussion 6.5.1.Sensitivity of the Multimodal Routing with respect to the CO 2 Emissions.First of all, we analyze the effect of the CO 2 emissions on the multimodal routing.Let the unit 13 Complexity emission costs vary from 50 to 100 ¥/ton, the sensitivity of the multimodal routing result with respect to this variation can be seen in Figure 7.
Figure 7 shows that the CO 2 emission costs increase linearly due to the fact that the emission volume is a constant (87.1 tons), while the rest of the generalized costs stay constant with the unit emission costs increasing.Even though we continue the simulation until the unit emission costs reach up to 1000 ¥/ton, which is obviously infeasible in practice (according to the formulation of the Chinese Academy of Environmental Planning [73], the feasible unit emission costs should reach up to 100 ¥/ton at 2030), the sensitivity still keeps the same tendency as shown in Figure 6.It means that the best routes in the empirical case do not modify when the emission costs vary and are hence insensitive to this variation.Therefore, charging for CO 2 emissions within a reasonable range is not particularly helpful to promote the green multimodal transportation, at least in this empirical case.6.5.2.Bi-Objective Optimization regarding CO 2 Emissions.Furthermore, we modify the proposed model into a biobjective optimization model with the objectives of minimizing the generalized costs (Z 1 = Eq 1 + Eq 2 + Eq 3 + Eq 4) and of minimizing the CO 2 emissions (Z 2 = Eq 51). minimize Using the widely utilized weighted sum method [43] defined in ( 51) and ( 52), we can solve the bi-objective 14 Complexity optimization problem and generate its Pareto solutions shown in Figure 7.
Especially, in Figure 8, compared with the last Pareto solution, the second last one can lower the CO 2 emissions from 87.1 to 83.7 by 3.90% and cause only a slight increase of the generalized costs from 1571.2 to 1586.0 by 0.94%.By comparing Figure 8 with Figure 7, we can conclude that a bi-objective optimization approach is much more effective than charging for CO 2 emissions in helping decision makers to control and optimize the CO 2 emissions in the multimodal service network.6.5.3.Sensitivity of the Multimodal Routing with respect to the Confidence.Confidence α in the fuzzy chance constraint ( 8) is set by the decision makers and reflects their preference for the reliability of the multimodal routing such that all the loaded volume of a rail service does not exceed its carrying capacity.Its value has a remarkable effect on the optimization result of the multimodal routing, which can be investigated by using sensitivity analysis.Let α vary from 0.1 to 1.0 with a step size of 0.1, we can obtain the sensitivity of the multimodal routing result (indexed by the minimal generalized costs) with respect to such variation, which can be seen in Figure 9.
We can also draw some helpful insights from Figure 9 summarized as follows.
(1) Overall, larger values of confidence α will lead to larger generalized costs for the best routes.
(2) The variation of the generalized costs with respect to the confidence is stepwise.In this case study, the multimodal routing is sensitive to the confidence when its value exceeds 0.3.
(3) The economy and reliability of the multimodal routing cannot reach their respective optimum simultaneously.The economy will be sacrificed if the decision makers decide to improve the reliability of the planned best routes, and vice versa.
6.5.4.Fuzzy Simulation.During the decision-making process regarding the multimodal routing, it is challenging for decision makers to determine an objective confidence value.In order to help them select reasonable confidence values, fuzzy simulation is adopted in this study.The fuzzy simulation simulates the actual deterministic transportation scenario by randomly generating the carrying capacities of rail services according to the membership degree of the triangular fuzzy numbers.The simulation process is shown in Figure 10 [8].
After the simulation, we can first gain deterministic carrying capacities v ijs for i, j ∈ A, s ∈ Γ ij .Eq. ( 8) can be consequently transformed into a deterministic carrying capacity constraint as (53).
Then, we can check if the planned given by the fuzzy programming model satisfy (15) or not.If the planned routes satisfy the constraint, we define them as a successful plan; otherwise, a failed plan which means that the planned routes are infeasible due to the violation of the capacity constraint when adopted to move containers in the practical transportation.If the transportation starts and a plan is found to be infeasible due to the insufficient actual capacities of some rail services, an alternative plan should be generated according to the following principle in order to reduce the additional costs and CO 2 emissions caused by the alternative: the rail service with insufficient capacity should be preferentially used by the transportation orders with larger container volumes initially assigned to it as many as possible on the condition that (53) must be satisfied, that is, the transportation orders with larger container volumes should  15 Complexity be accomplished by the initial plan as far as possible, and the rest transportation orders usually with smaller container volumes assigned to the rail service should be accomplished by the origin-to-destination road services.The fuzzy simulation should be conducted several times.In this study, we ran the simulation 50 times, and the fuzzy simulation result is shown in Figure 11.
As we can see from Figure 11, the successful times of the planned routes increase remarkably when the confidence value changes from 0.7 to 0.8.Although the generalized costs increase approximately by 1.29%, the success ratio of the planned routes considerably increases from 26% to 100% by more than 3.8 times when the confidence value is set to 0.8 instead of 0.7.
Note that when α = 0 5, according to ( 25) and ( 26), the fuzzy chance constraint can be rewritten as ijs which is the deterministic capacity constraint widely employed by the current literature (e.g., Sun and Lang [2] and Sun et al. [19]).So we know that the successful ratio of the multimodal routing with deterministic capacity consideration is only 12% in the 50 simulations for the specific case.
As we can see from Figure 11, the successful ratio of the multimodal routing can be significantly improved by more than 8.3 times by considering the capacity uncertainty and setting the confidence to 0.8, 0.9, or 1.0.Therefore, it is of great value to formulate the capacity uncertainty when planning the multimodal routes, so that higher reliability of the multimodal routing decision-making can be guaranteed.We can further achieve the optimization results of the 50 deterministic problems formulated by the linearized model with (25) and (26) replaced by (53).The optimization results are presented in Appendix D (available here) in the    17 Complexity supplementary material.Such results can be considered as the actual best routes.By comparing the actual best routes with the planned ones, we can find the confidence value that best matches the actual transportation.
In the 50 fuzzy simulations, the best actual routes whose generalized costs are equal to ¥ 1,575,559 emerge 38 times that account for 76% of the entire simulation times (see in Figure 12).That is to say, such routes are the most likely best ones in the practical transportation.The planned best routes given by the fuzzy programming model with confidence values of 0.8, 0.9, and 1.0 match the most likely best actual routes better than others.Above all, the best values of the confidence in this empirical case are recommended to be 0.8, 0.9, and 1.0.6.6.Managerial Implications.It should be noted that the optimization results and conclusions in this section are drawn based on a specific case whose setting is shown in Sections 6.1 and 6.2.The results and conclusions are not universal and are sensitive to the case itself and might change if the case changes.Therefore, the above analysis and discussion in this section only provide a procedure and a demonstration on how to adopt the proposed methods to deal with the practical problems.Although the cases might be different, we can always use such procedure to address them to draw corresponding optimization results and conclusions.
Despite this fact, the results show important implications for the practical implementation of the model: (1) Although the multimodal service network is complex, the linearized optimization model can deliver the optimal solution in a very short time so that new plans for incoming orders can be found very fast.
(2) The objective of the optimization is to minimize the total costs consisting of different cost categories Step 1: Step 2: Step 3: Step 4: Step 5: Step 6: Step 7: Step 8: Generate a random number Calculate its membership degree: , if v ijs ≤ v ijs ≤ v ijs   18 Complexity (i.e., transportation, transshipment, inventory, penalty, and emission costs).However, there is a clear dominance of transportation costs over all other cost categories which makes it difficult to influence the optimal solution by changing other cost factors than transportation costs.
(3) Based on the previous implication, it is clear that the charging method for CO 2 emissions does not improve the environmental sustainability of the multimodal routes if realistic values are used for emission costs.Therefore, rather a bi-objective approach should be used in order to be able to compare the costoptimal and emission-optimal solution.In this case, often a solution can be found where one objective is significantly improved whereas there is only a small increase for the other one.
(4) The minimal values of total costs can be usually achieved under the deterministic setting since no buffer times or capacities need to be included.However, uncertainty plays an important role in practice and therefore should be considered also in planning.
Although the total costs are increasing with the increasing confidence level, this also reduces the risk of disruptions during transportation which can be seen in a lower ratio of failed plans.As a result, more transportations can be performed as originally planned and the need for re-planning for failed plans is minimized.Since such re-planning usually leads to additional costs, reducing the share of failed plans is of high interest.

Conclusions
Nowadays, the multimodal transportation has been widely promoted all over the world, and its routing optimization has attached great importance by both researchers and practitioners.Meanwhile, the public concern on environmental protection raises, and environmentally friendly development gets tremendous highlights in almost every industry.Therefore, it is valuable to combine multimodal routing with the environmental views.In addition to that, different uncertainty factors can influence the performance of multimodal transportation and therefore should be considered already in the planning phase.However, all these factors increase the complexity of multimodal routing modelling and therefore make the transportation planning very challenging in practice.
In order to respond to these challenges, we developed a nonlinear mixed-integer optimization model covering multiple objectives (i.e., economic and environmental) and multiple uncertainty factors (i.e., traffic congestion and rail service capacity).Moreover, in order to facilitate the solution of this problem, a linearization approach was proposed.The application of the proposed model to a real-world case study showed that it can deliver optimal results in a relatively short time which is very helpful for the practical transportation planning.Besides that, it was shown that the inclusion of environmental criteria is meaningful only when the bi-objective method is used, since the weight of the emission costs is very low in comparison to transportation costs if the charging method is used.Regarding the reliability of transportation, including the uncertainty factors into planning and increasing confidentiality levels significantly reduce the risk of disruptions during the execution phase but also lead to increased total costs.However, the reduced risk of disruption also reduces the need for re-planning and therefore helps to avoid additional costs connected with this process.
Compared with the models proposed by other studies listed in the literature review section, our model can comprehensively formulate the issues of road traffic congestion and rail service capacity uncertainty that actually exist in the routing decision-making process and thus have better feasibility in dealing with the practical problem.Another improvement is that by modifying the objective function, our model can easily use two different approaches, including charging for CO 2 emissions and bi-objective optimization, to separately optimize the CO 2 emissions and identify which one is more effective in planning green routes.Although the model is nonlinear, it can be linearized, which enables the problem to be effectively solved by any exact solution algorithm on any mathematical programming software.
This study mainly focuses on the modelling and the utilization of a classical exact solution algorithm based on a linearized model and does not involve any intelligent algorithms.The exact solution strategy can efficiently optimize the presented case, but its performance is doubtful when the scale of the case is increasingly expanded.As a result, for further improvement to this study, we will test the feasibility of the exact solution strategy in dealing with a larger-scale empirical case.If the test result is not satisfactory, we will try to design some advanced intelligent solution algorithms (e.g., heuristic algorithms) with higher solution efficiency, so that the problem can be effectively optimized and routing suggestions can be provided to the decision makers in time when the problem scale gets larger.

Figure 1 :
Figure 1: Basic topological structure of the multimodal service network.

Figure 2 :
Figure 2: Diagram of a simple road-rail multimodal route.

Figure 3 :
Figure 3: Time-dependent travel time formulated by piecewise linear function.

i
variable.x k ijs = 1 if service s on arc i, j is used for transportation order k; x k ijs = 0 otherwise.y k i Arrival time of the containers of transportation order k at node i. z k i Planned departure time of the containers of transportation order k from node i. n k i A nonnegative integer variable that indicates how many periods (24 hours) z k i exceed.t k ijs Travel time of road service s on arc i, j when used for transportation order k, unit: h.t k ijs = f ijs z k i where f ijs is the piecewise linear function of the travel time with respect to the corresponding planned departure time.Note that in this function, the input departure time z k i should fall within range [0, 24].Otherwise, it should be first normalized into such range before being input into the function.Parameters Representations l s i Loading/unloading operation start time of rail service s at node i. u s i Loading/unloading operation cutoff time of rail service s at node i. v ~ijs Fuzzy available capacity of rail service s on arc i, j and v ~ijs = v min ijs , v M ijs , v max ijs where v min ijs < v M ijs < v max ijs , unit: TEU.d ijs Travel distance of service s on arc i, j , unit: km.c ijs Unit transportation cost when service s is used to move containers on arc i, j , unit: ¥/TEU.c s Unit loading/unloading operation costs of service s, unit: ¥/TEU.c s Unit inventory costs of transportation service s at node i, unit: ¥/TEU-h.ck dest Unit destination inventory costs for transportation order k, unit: ¥/TEU-h.c k pen Unit penalty costs for transportation order k caused by delayed delivery at its destination, unit: ¥/h.c co 2 Technique 5.2.1.Crisp Equivalent and Linearization of the Fuzzy Chance Constraint.From the definition of the fuzzy creditability, (24) can be obtained where a is a deterministic number while b is a triangular fuzzy number described by b 1 , b 2 , b 3 and b 1 < b 2 < b 3 .[69]:

26 5. 2 . 2 .
Linearization of the Piecewise Linear Function.Similarly to(24), the piecewise linear function f ijs • is also a nonlinear expression where • represents z k i − 24 • n k i for the sake of conciseness.Let P ijs represent the set of the linear segments of f ijs • , p the index of the segments (p = 1, … , P ijs ), time interval corresponding to the linear segment f p ijs • .Then, f ijs • has a universal formulation shown as

Step 1 . 28 Step 2 .
Define ∀k ∈ K, ∀ i, j ∈ A, ∀s ∈ Ω ij Define two nonnegative variables λ kp− ijs and λ kp+ ijs , and distribute them to the lower bound tt p− ijs and upper bound tt p+ ijs of the corresponding departure time interval tt p− ijs , tt p+ ijs as their weights.Hence, we have the following equations.

Figure 4 :
Figure 4: Framework of the mathematical model.

Figure 5 :
Figure 5: Multimodal service network for the empirical case.

Figure 6 :
Figure 6: Structure of the minimal generalized costs.

Figure 8 :Figure 9 :
Figure 8: Pareto solutions to the bi-objective optimization problem.

Figure 12 :
Figure 12: Ratios of the optimization results in the 50 fuzzy simulations.

Table 1 :
Literature review of the articles in the multimodal routing problem.
U k are separately the lower bound and upper bound of the due date time window, respectively.The rest of the symbols in the mathematical model and their respective representations are summarized as follows in Table2.
4.1.Notations.In this study, G = N, A, S represents a roadrail multimodal service network in form of a directed graph, where N, A, and S separately denote the node set, arc set, and transportation service set in the network.For each transportation order k (∀k ∈ K where K is the transportation order set), its five attributes are represented as origin o k , destination d k , volume q k (unit: TEU), earliest release time t k release at origin, and due date t k due = t L k , t U k where t L k and t

Table 2 :
Notations.Indices of the nodes, i, j, and h ∈ N. s, rIndices of the transportation services, s and r ∈ S. h, iDirected arc from node h to node i, and h, i ∈ A. i, j

Table 3 :
[74]es of the parameters in the empirical case.Road service; * * rail service.*** Suggested by the research group at the Chinese Academy of Environmental Planning[74].

Table 4 :
Scale of the empirical case problem.

Table 5 :
Optimization result and strategy performance.

Table 6 :
Best routes in the empirical case.Figure 7: Sensitivity of the multimodal routing with respect to the CO 2 emissions.