Multi-Objective Optimization of Customer-Centered Intermodal Freight Routing Problem Based on the Combination of DRSA and NSGA-III

: The satisfaction of requirements and preferences of shippers is critical to enable the practi-cability of solutions that are derived from intermodal transportation routing problems. This study aims to propose a decision process to help shippers participate better in routing decisions. First, we considered shippers’ requests on transportation cost, timeliness, reliability, and ﬂexibility to construct a multi-objective optimization model. Then, to solve the interactive optimization method that was proposed, NSGA-III was applied to obtain the Pareto front and dominance-based rough set approach to model the preference information. Finally, a case study was conducted and an expert was invited as decision-maker to demonstrate the applicability of the proposed model and the effectiveness of the interactive method for shippers. The results are expected to provide shippers with more rational transportation schemes and insights for the sustainable development of intermodal transportation.


Introduction
The growth of international trade has led to growing interest in intermodal transportation among transportation planners. Intermodal transportation is considered a strategic tool to benefit from the advantages of transport modes, e.g., relieving congestion in road mode [1], environmental benefits [2], better transport reliability, and low cost [3,4]. Different combinations of transportation modes lead to different costs and service levels, which poses a challenge for intermodal operators to transport shippers' freight. Enterprises hope to improve logistics efficiency and effectiveness to maintain their competitiveness under the fierce market competition. Therefore, how to optimize the choice of freight transportation routes in intermodal transportation networks has been highly concerned by enterprise decision-makers and intermodal operators.
Intermodal Freight Routing Problem (IFRP) works to optimize the combination of various transportation modes to form optimal routes for commodities from the origins to the destinations through an intermodal network [5]. The analysis of this problem needs to be based on its several important characteristics, such as multiple objectives, delivery time windows, and uncertainty. These characteristics have been commonly analyzed and discussed in related studies.
From the perspective of freight transportation planning, IFRP is generally a decision that is made by intermodal operators that is based on customer needs. Operators need to consider various factors in transportation operation and management comprehensively in order to achieve effective planning [6]. Previous research on IFRP has focused on helping operators to optimize the overall benefits of network operations around the objectives of transportation cost, time, environmental impact, and transportation risk.
Customer needs are also multidimensional in transportation services. In addition to reducing transport costs and improving timeliness, they also expect more flexible and reliable freight transportation [7]. However, few efforts have been made to integrate these factors into their modeling processes, which may make customers fail to obtain optimal transportation services. Moreover, most previous studies on IFRP applied non-preference, priori, or posteriori approaches to solve their models and did not consider or simply deal with customers' preference information (mainly using assumed weights). Disregarding preferences produces multiple optimal solutions and increases the hesitation of Decision Makers (DMs). The lack of in-depth analysis and modeling of preference information increases the subjectivity of optimal solutions and customers' burden.
In general, existing research facilitates operators' operational decisions but makes it difficult to help operators meet the heterogeneous needs of customers and improve service satisfaction, even when the shippers have the power to make routing decisions.
To fill these research gaps, this study aims to answer the following problems to help customers to obtain the most satisfactory transportation plans in IFRP.
(1) What optimization objectives should be considered to express shippers' requirements for transportation services and how can they be modelled? (2) How to choose the optimization method to truly reflect the DMs' preference and guide the decision-making process?
This study first analyzes and models several important characteristics of IFRP, including scheduled and time-flexible service patterns, time uncertainty, and delivery time windows. Then, based on the analysis of the shippers' interest, this study proposes a multi-objective optimization model of IFRP with transportation cost, timeliness, reliability, and flexibility as objectives. Referring to the relevant literature, the functions of cost and timeliness are determined using generalized cost structure and fuzzy soft time windows, respectively. A reliability modeling method that is based on Monte Carlo and a flexibility modeling method considering service frequency and the latest booking time are proposed for the first time. Third, the solution method should be efficient and reflect the DMs' requests effectively. As such, an interactive optimization approach is proposed to solve IFRP, where the NSGA-III algorithm that is suitable for solving multi-objective optimization problems with constraints is combined with a classical preference modeling method Dominance-based Rough Set Approach (DRSA). This approach can help DMs to obtain satisfactory solutions well while imposing less burden on them and allowing them to learn during the solution process. Finally, we conduct a case study to illustrate the validity of the proposed models.
The rest of this paper is organized as follows. In Section 2, a brief review of the relevant literature is presented. In Section 3, we propose a multi-objective optimization model for IFRP. Then, in Section 4, an interactive multi-objective optimization method combining DRSA with NSGA-III is developed. A case study is conducted in Section 5 to verify the effectiveness of the proposed model and method. The relevant findings are summarized in Section 6. Finally, the conclusions are drawn in Section 7.

Formulation Characteristics of Optimization Models in IFRP
IFRP is typically complex due to its multiple characteristics, such as multiple objectives, multiple transportation service patterns, demand delivery time, uncertainty, and limited network capacity. According to customers' requirements, effective freight transportation planning needs to weigh the cost against multiple service quality indicators to achieve fast, economic, and reliable transportation [6]. Service patterns can be divided into scheduled services and time-flexible services. Usually, rail, sea, and air transport have a fixed departure time, while trucks can travel flexibly in time [5]. If goods are transshipped from a truck to a scheduled service, generally they need to be operated at the latter's specified service time and the service stations [8]. Ziliaskopoulos and Wardell [9] and Moccia, Ropke, and Valentini [10] investigated the optimal path problem through intermodal networks with multiple service patterns. Customers have different demand delivery time requirements and the planned feasible routes need to deliver goods at a specified time [11], the forms of which are usually time points or time windows. Chang [12] described the international intermodal freight routing problem as a multi-objective multimodal multicommodity flow problem with time windows and concave cost. In addition, they summarized transportation mode schedules and the delivery time as two important characteristics in IFRP. Consideration of uncertainty in intermodal transportation planning is one of the key factors in providing adequate support to DMs [13]. Previous studies mainly focused on the uncertainty of time, demand, and network capacity. The latter two factors and network capacity constraints are often correlated with each other. That is, the satisfaction of demand is constrained by the network capacity and when demand increases or network capacity decreases, this constraint is likely to lead to infeasible transportation routes [14,15]. Considering the importance of transportation time uncertainty in IFRP, it has been discussed widely by traffic engineers. Hrušovský et al. [16] and Zhao et al. [17] developed models to measure the travel time uncertainty based on stochastic programming. Zhao, Xue, et al. [18] subsequently introduced stochastic transfer time. Sun and Li [19] suggested the road travel time and the operation time as time uncertainty sources of road-rail intermodal transportation and modelled uncertainty parameters as fuzzy variables.

Determination of Optimization Objectives
In previous studies, cost and time were the main optimization objectives for responding to customer needs. Minimizing cost reflects the customer needs for improving the transportation economy and is modelled as an objective in many IFRP studies. The charges for accomplishing transportation orders are usually formulated referring to generalized costs, the necessary components of which are operating costs, handling costs, and storage costs [19].
The time objective involves not only reducing transportation time but more importantly, improving delivery timeliness. Either early or late delivery affects the costeffectiveness of the shippers, especially when enterprises adopt a just-in-time production strategy [20]. Generally, transportation services will be considered timely when goods are delivered at a specified time by the shippers. Time windows have been widely used in modeling the delivery time in IFRP and can be divided into hard time windows and soft time windows. Subsequently, considering that customers can accept the violation of delivery time to a certain degree, Sun and Li [19] further proposed fuzzy soft time windows to model the delivery time and timeliness.
The above studies widely discussed the optimization of the overall operational benefits from the perspective of operators. However, the interests of shippers and operators are different and previous studies failed to estimate the preferences of shippers. The final plans that are given by operators without comprehensive preference analysis for customers are generally not optimal for shippers. The literature on the mode choice problem conducted extensive empirical investigations on shippers' preferences. In those studies, transportation cost, time, service frequency, flexibility, and delay risk, were frequently used to compare transportation services [21][22][23]. The delay risk arises from the transportation time uncertainty and the service frequency is related to flexibility [24]. With the development of intermodal transportation, flexibility and reliability have become important factors for DMs to weigh intermodal services against single-mode services [13]. Therefore, transportation reliability and flexibility should also be considered in customer-centered IFRP.
Transportation reliability, which is generally defined as the on-time delivery ability of transportation services, is considered by modeling the uncertainty in IFRP. Uddin and Huynh [15] proposed that the transportation reliability is significantly affected by the operational uncertainty of transportation networks, i.e., the network capacity uncertainty, and further leads to time uncertainty. Sun and Li [19] defined the transportation reliability as a measure of the success of the planned intermodal routes to deliver containers to customers and investigated the problem of reliable route selection under travel time uncertainty and loading/unloading time uncertainty. Despite the growing number of studies on transportation reliability, it has hardly been modelled as an objective in IFRP.
Few studies have examined the transportation flexibility in IFRP. In mode choice studies, flexibility mainly indicates the ability of the logistics service providers to respond to changes. Khakdaman, Rezaei, and Tavasszy [25] offered a comprehensive definition of flexibility, i.e., the capacity of a logistics system to respond to shipper needs for changes of the service components that are provided (e.g., destination, transit time, or volume) at any point in time during service. A limited number of studies have been done on flexibility, which mainly focused on the response to partial changes, especially transport time. For example, Norojono and Young [26] showed that flexibility is a function of frequency and problem responsiveness. Beuthe and Bouffioux [27] measured the responsiveness in terms of the minimum notice time of transport orders.

Multi-Objective Optimization Approaches
Multi-objective optimization consists of three phases: model construction, optimization (finding the Pareto front), and decision-making (selecting a single optimal solution). According to whether the DMs express preference information and the order of decision and optimization, multi-objective optimization methods can be classified into no-preference, priori, posterior, and interactive methods [28].
For multi-objective optimization of IFRP, most studies have explored methods that focus on the optimization phase rather than the decision phase. Moreover, they commonly did not consider or simply deal with customers' preference information. Xiong and Wang [29] and Peng, Yong, and Luo [30] explored the application of improved genetic algorithms to obtain Pareto sets without considering the decision process. Chen et al. [31] solved the optimization model of IFRP considering container use applying the normalized normal constraint method. In IFRP studies considering the DMs' preferences, priori methods have been more frequently applied, i.e., inputting preference information before optimization. Chang [12] and Lei et al. [32] used the weighting method to integrate multiple objectives, classified the DMs' preferences using weight, and then applied the Lagrangian relaxation technique and the swarm intelligence method to solve the optimization models, respectively. Kaewfak, Ammarapala, and Huynh [33] combined 0-1 goal programming with the analytic hierarchy process to determine weights by a pairwise comparison of the objectives to generate optimal routes. There are also less studied applied posteriori methods in IFRP. Cho, Kim, and Choi [34] evaluated and obtained optimal solutions for IFRP by applying mathematical and multiple attribute decision-making models to the case after obtaining the Pareto front, respectively.
It is evident that non-preference methods are difficult to meet the real needs of customers. The applied priori or posteriori methods, which separate the decision process from the optimization process, tend to impose a greater burden on DMs due to the design limitations. Compared with the above methods, Interactive Multi-objective Optimization (IMO) methods, which have not yet been applied to IFRP, are considered promising methods in multi-objective optimization for solving real-world problems [28] and deserve further research. It is because IMO methods have alternating processes of the decision and the optimization, and often help DMs to obtain satisfactory solutions with less burden on them and allow them to learn during the solution process.
IMO methods aim to incrementally find the optimal solutions for DMs guided by the DMs' preferences. The basic elements of the methods are mainly search engines, preference information, and preference models [35]. The search engine is a decisive factor in solving problems and can be divided into mathematical programming and non-mathematical programming techniques (mainly referring to heuristic methods). The optimal solution to a large-scale or a real-world freight routing planning problem is difficult to obtain by a single exact solution method and heuristic algorithms have proven to be feasible in solving this optimization problem [5], such as genetic algorithms [29], swarm intelligence [32], and memetic algorithms [36]. Preference information can be divided into three categories: expectations, comparison of objective functions, and the comparison of solutions. The frequently used preference models are commonly divided into value functions, dominance relationships, and decision rules [35]. The type choices of preference information and preference model need to be based on the actual situation of the DMs, while reducing the stress of the DMs and controlling the cognitive bias to obtain complete preference information.

Review Summaries
As we can see from Sections 2.1-2.4, the existing literature on IFRP has made sustainable achievements, both in terms of modeling and optimization methods. However, the literature review also shows the existence of the research gaps as follows.

1.
Studying IFRP from the shippers' perspective and systematically considering the optimization objectives to improve the rationality of routing selection for customers 2.
Optimizing the customer demand on transportation reliability and flexibility by modeling rationally.

3.
Applying an interactive approach to solve the multi-objective optimization of IFRP, combining the preference models with heuristic algorithms to help the DMs to obtain the most preferred solution more efficiently and effectively.

Problem Description
A single-commodity intermodal freight routing problem aims to select an optimal complete route from the origin to the destination through an intermodal transportation network [5]. To better demonstrate the problem, a simplified network diagram containing 8 nodes and 17 directed line segments is proposed in Figure 1. Nodes O and D denote the origin and the destination of the commodity flow, respectively, and T1-T6 denote the transshipment nodes in complete routes from O to D. The directed line segments represent the transport channels containing the available transportation services between two related nodes. Services at nodes and line segments have different attribute values such as cost and time, which determine different cost and service benefits of feasible paths.

Modeling Important Characteristics and Process
As mentioned earlier, we consider and model the important characteristics of IFRP as follows.

Multiple Transportation Service Patterns
Transportation services can be divided into scheduled services and time-flexible services.
Referring to the analysis of Moccia et al. [10], a scheduled transportation service has a fixed route, the expected departure time at the loading terminal with a service time window, i.e., a closed time interval from the operation start time to the operation cutoff time, and the expected arrival time at the unloading terminal. The service time window [Dl − , Dl + ] of a scheduled transportation service has the following constraint to forwarding goods: (1) if the batch of goods arrives at the loading terminal of s before the operation start time Dl − , they must be stored first until the operation starts. (2) if the arrival time is later than the operation cut-off Time Dl + , these goods cannot be shipped by the service s.
Due to the consideration of travel time uncertainty, the expected arrival time for the scheduled services will not exist in this study. Time-flexible transportation services are not subject to planning restrictions, the cargo can be loaded onto vehicles immediately after arriving at the vehicles' loading terminals. If goods need to be stored in the service terminals, generally these terminals provide for free storage time and chargeable storage time.

Time Uncertainty
This study considers travel time uncertainty and operation time uncertainty at nodes, which are mainly modelled as random variables or fuzzy variables in related studies [19]. Referring to Zhao, Xue, et al. [18], we modelled these two uncertain parameters as random variables.

Demand Delivery Time
To better express the customer needs, this study adopts the fuzzy soft time windows that were proposed by Sun and Li [19] to model the delivery time requirements of shippers and explain the timeliness. The method uses trapezoidal fuzzy numbers to model the delivery time windows, allowing customers to express a satisfactory time range and a larger allowable time range. The fuzzy membership degree indicates the customer satisfaction with the delivery time, i.e., the timeliness level. First, let the shippers choose a satisfactory time interval and when the moment completing a freight order falls into this interval, the satisfaction will reach the maximum. Then, choosing the earliest and latest allowable time points, four parameters of a fuzzy soft time window, T = (T 1 , T 2 , T 3 , T 4 ), are formed. The fuzzy membership degree function is then used to measure the customer satisfaction, denoted by µ(T) (see Figure 2).

Nodal Time in Intermodal Transportation
The nodal process of intermodal transportation in this study consists of three parts: the shipment process at the origin, the delivery process at the destination, and the transshipment process at the other nodes. Combined with the previous problem characteristics, we analyzed that the nodal time can be divided into variable time and fixed time. The variable time can be further divided into the variable operation time and the storage time.
Taking the time in the transshipment nodes as an example, there are two types of transshipment processes in Figure 3, including from the transportation service r to the time-flexible service s and r to the scheduled service s − . The nodal time for r to s is not only a variable time, but also a variable operation time, consisting of the unloading time for r and the loading time for s. It is because the fixed time and storage time are only relevant for the scheduled services. The fixed time is the time from the operation start time Dl − to the departure time for service s − . The storage time happens when the cargo is stored in the service terminals of the scheduled services due to waiting for the operation start. It is variable and determined by the unloading time of service r, drayage time and the operation start time of service s − as Figure 3, but at the origin by the release time that is specified by the shippers and the operation start time of the scheduled services. The rest of the nodal time for r to s − is defined as the variable operation time. Similarly, the time division at the origin or destination can be obtained. In this study, the variable operation time at each node was modelled as a random variable.

Assumption
The modeling is required to follow several assumptions to make it rigorous:

1.
All containers are ISO standardized and are measured in 40 ft equivalent units (FEU) 2.
Containers of batches are prepared for loading at their release time and delivered when completing the unloading operation at the destination. 3.
The scheduled transportation services always have a fixed departure time and service time windows. Storage cost occurs when storage time exceeds the free time and is proportional to storage days.

4.
Travel time for all modes and variable operation time at each node are random variables. The values of the latter are taken only concerning transportation modes involved.

5.
The intermodal transportation network is assumed to be uncapacitated, which means all nodes, links, and terminals have sufficient capacity to complete a given transport or transshipment task.

Defining Symbols and Parameters
The symbols including sets, indices, and variables in the intermodal freight routing problem in this study are defined in Table 1, in which they are classified into four types, all the cost parameters are in dollars and the time parameters are in hours if not specified. Using these symbols, a multi-objective optimization model is developed in the next section.

B
Set of freight batches that need to be transported The origin of all batches of freight d The destination of all batches of freight t b 0 Release time of freight batch b at the origin claimed by the shipper µ The timeliness level of transportation routes The soft delivery time window of freight batch b specified by shippers An allowable maximum value for changes of expected transportation time claimed by shippers for freight batch b ε r i Free storage time set by the service terminal of service r at node i AT l

Symbols for the Transportation Network
The time completing freight transportation in route l E(AT l ) Expectation value of AT l , involved with multiple random variables.

Constructing the Objectives of the Optimization Model
The goal of the optimization model is to minimize the total cost and maximize the objectives of timeliness, reliability, and flexibility.
Cost objective: Minimize the total cost, including the cost for the departure operation at the origin and the arrival operation at the destination, travel cost, transshipment operation cost at other nodes, and the chargeable storage cost at nodes due to waiting for shipment onto scheduled transportation services. The cost model is as shown in Equation (1).
Timeliness objective: Maximize the timeliness level. First, AT needs to be obtained, which consists of the release time of a freight batch, the used time at the origin, and the destination, travel time, and transshipment time at other nodes. Then, combining that with the delivery time windows of freight batches, the previously mentioned µ(x) function was used to calculate the satisfaction level of the delivery time, i.e., timeliness level. Taking stochastic time uncertainty into account, the time completing the freight transportation is also random. In this study, we used the random expectation function E(x) to handle the time variable. The objective model is as shown in Equation (2).
Reliability objective: Maximize the possibility of on-time delivery of feasible routes, i.e., the reliability of delivery time. As reviewed before, this study considers the transportation reliability with stochastic travel time and the stochastic operation time at nodes. It is difficult to quantify the reliability of intermodal transportation under multiple parameter uncertainties, and related studies generally use simulation methods [37] or estimate the overall distribution [38] to quantity.
In this study, A Monte Carlo method was used to evaluate the on-time delivery probability. The applied steps are: (I) for route l, the calculation of AT l is from a random vector t T that is composed of multiple time variables, which obey the specified distributions. (II) based on all the variables' distributions, generate M random arrays t k T , k = 1, 2, . . . , M. (III) for route l, calculate AT l for M random arrays. If there are AT l values of M l arrays falling in the interval E(AT l ) − δ b , E(AT l ) + δ b , then according to the law of large numbers, the following equation is used to estimate the value of θ for route l: Thus, the reliability objective model is shown in Equation (3).
Flexibility objective: Based on the previous review, this study defines the flexibility objective as maximizing the ability to adapt to customer demand changes for transportation time before transportation and proposes a modeling approach that is determined by two factors, the service frequency f and the latest booking time To of the transportation services.
Specifically, the importance of these two parameters is treated as the same and the extreme value normalization method is used to make the two parameters dimensionless. The model for the flexibility objective is shown in Equation (4).

Presenting the Constraint Sets of Optimization Model
Equation (5) is the flow conservation constraint that ensures the continuity of the feasible routes. Equations (6) (9) and (10) represent the time at the origin, destination, or other nodes. Equation (11) indicates the moment for freight arriving at node i. For equation (12), it represents the permissible delivery time window constraint that AT of freight transportation must satisfy. Equation (13) means that the freight transportation must satisfy the constraint of the service time window of a scheduled transportation service if it needs to be shipped by it. Finally, Equations (14) and (15) are the chargeable storage time for the freight batches at the origin or transshipment nodes.

An Interactive Optimization Approach
In multi-objective optimization, many Evolutionary Multi-objective Optimization (EMO) algorithms such as NSGA-II, have difficulties in handling many objectives due to too many non-dominated solutions, not keeping diversity well, the difficult representation of trade-off surface, and so on. NSGA-III was proposed to release those difficulties. It maintains the diversity of non-dominated solutions and has been shown to work well from testing and other problems with 3-15 objectives [39]. Preference models should retain complete preference information, and the processing should be easy to understand for DMs. DRSA requires only a simple classification of a sample of solutions by DMs, then deduces transparent and easy to interpret decision rules for DMs. This "glass box" rather than "black box" model improves the quality of the interaction and allows the DMs to receive the resulting recommendations well. Moreover, DRSA has been successfully applied to IMO [40]. More importantly, DRSA complements any multi-objective optimization method nicely, especially EMO methods for obtaining the optimal set [28]. Greco, Matarazzo, and Słowiński [41] proposed the basic idea of DRSA that was applied in interactive evolutionary multi-objective optimization, i.e., DRSA-EMO, which successfully integrates these two methods. Therefore, this study explores an IMO method that is based on the combination of NSGA-III and DRSA to solve the constructed model.
According to the relevant steps of DRSA-IMO [40] and DRSA-EMO [41], the complete interaction process can be divided into two phases: the computation phase and the dialogue phase. In the computational phase, a Pareto set is obtained by solving the NSGA-III algorithm after iterations. Then we enter the external interaction phase. To reduce the cognitive burden of DMs, no more than 15 solutions need to be shown to DMs in this study. We designed it that if the number of solutions in the set is less than six, the DMs are asked to indicate the most preferred solution as the final solution directly. Otherwise, they still need to classify the solutions of the representative set. Next the DRSA method is utilized to derive the decision rules and then the only one rule that is recommended by DMs needs to be transformed into constraints for the objective values, which are added to the algorithm for the next round of solving.
The following interaction steps are designed for the study.
Step 1. Initialize the algorithm parameters and set the interaction counter it = 0.
Step 2. Generate an initial population that is based on the coding rule. Then run NSGA-III for 500 generations and output the final Pareto front X and objectives values for every x ∈ X.
Step 3. If the number of objects in X is less than 15, then go on; Otherwise, randomly select 15 solutions as the input for the next step.
Step 4. Present these solutions combined with a parallel coordinate plot to DM.
Step 5. If the number of solutions is not more than five, DM will be required to specify the most satisfied one; otherwise, go on.
Step 6. The DM selects some solutions as "better".
Step 7. Apply the DRSA method in the "better" solutions to derive a set of decision rules "if C(x) ≤ a 1 , T(x) ≥ a 2 , R(x) ≥ a 3 , and F(x) ≥ a 4 , then, the solution x is better".
Step 8. Present these rules to the DM and require him or her to select the most important rule in those.
Step 9. Transform the conditional part of the final rule into a new constraint for the proposed problem.
Step 10. Update it = it + 1 and go to Step 2.
In the interactive process, the DMs can express preferences for the solution set output by the optimization algorithm, which can then be fed back to the algorithm to guide the convergence, thereby helping obtain the ideal solution for the DMs. Moreover, the process is reversible; DMs can always go back to the Pareto front that is considered in the previous iteration and continue from that. This is a learning-oriented optimization concept, i.e., the interactive process allows the DMs to understand their preferences and the "shape" of the Pareto set [40]. In this way, the advantages of the two methods are combined.

NSGA-III
Deb and Jain [39] proposed NSGA-III on the basis of NSGA-II and subsequently extended it to the problem with constraints [42]. Unlike NSGA-II which uses crowding distance to select individuals of the critical layer in the elite retention process, NSGA-III maintains diversity by a set of uniformly distributed reference points. The application framework in this study is shown in Figure 4. Step 1. Input algorithm parameters, including the size of population N, maximum number of iterations maxGen, variation rate, crossover rate, etc.
Step 2. Generate uniformly distributed reference points on the inner and outer hyperplanes.
Step 3. Determine the encoding and decoding rules, initialize the population P t , and set the initial value of interaction counter it = 0.
Step 4. Generate the offspring population Q t using the selection operator, the improved crossover operator, and the variation operator.
Step 6. Combine the parent population P t and the offspring population Q t , the size of the new population is 2N.
Step 7. Fast non-dominated sorting of the combined population to obtain nondominated layers F 1 , F 2 , · · · , F L .
Step 8. Put individuals of layers starting with F 1 into a new population S t layer by layer, until the size of S t will be no less than N for the first time if all the individuals in layer l are put there. If the size of S t = N, then finish step 10; Otherwise, go to Step 9.
Step 9. Normalize the objective values of all the individuals in S t and reference points so that both sets are in a range. Then, associate each population member of S t with a reference point and use the niche-preservation operation to obtain the remaining individuals for S t from the layer l.
Step 10. S t become the new parent population. If Gen ≤ maxGen, then go back to Step 4; otherwise, stop.
The adaption parts for the proposed problem in this study are now described, which include two aspects.

Coding Rule and the Choice of Crossover and Variational Operators
In this study, we use the priority-random coding method that was proposed by Hong et al. [43], which guarantees that any path that is obtained by decoding corresponds to a complete path from the origin to the destination in a given intermodal network. The method first assigns a number to each of all the nodes, followed by a priority (a non-repeating integer from 1 to n) to each and uses the priority sequence as decision variables. Subsequently, based on the backward and forward relationships among the nodes in the given network, the backward nodes are selected one by one from the starting node. The inter-node transportation services are randomly chosen and a complete route is obtained finally.
A group of decision variables that is constituted by all the nodes' priority sequences in one route must take different values from each other as a way to avoid the inability of comparison between multiple backward nodes of a selected node. Therefore, referring to the above study, the partial matching crossover operators [44] and the reversal variation operators are chosen to perform the crossover and variation process.

Adaptative Improvement for the Interaction Process
NSGA-III is utilized to solve the proposed multi-objective problem model. In the interaction process, DRSA is used to obtain the DMs' threshold preference for objectives, then the preference information needs to be added into the problem model in the form of constraints to guide a new round of iterations. Individual vectors not satisfying any value constraints of objectives are treated as infeasible solutions.

DRSA
DRSA is an inference method for partially inconsistent and preference-ranked data that has been successfully applied to IMO. It aims to obtain a representation of the DM's preferences based on some exemplary decisions that are given by the DM for easily understandable "if . . . , then . . . " decision rules. According to Greco, Matarazzo, and Slowinski [40], exemplary decisions in DRSA can be preferences for dominance relationships among solutions, or classifications of solutions, such as "better" or "otherwise". In this study, we used DRSA for ranking problems, i.e., use the classification of solutions to give exemplary decisions. This section briefly introduces DRSA, more details can be available in Greco, Matarazzo, and Slowinski [45], and Slowinski, Greco, and Matarazzo [46].
DRSA is a method for deriving decision rules based on a set of basic concepts and mathematical relationships. Therefore, important concepts need to be understood first.
(1) Decision table   The DRSA method starts with an information table that is called the decision table. The form of the data table is a four-tuple information system DS = U, R, V, f . U is a nonempty infinite set of objects, U = {x 1 , x 2 , . . . , x n } and R is the set of attributes, R = C∪D. The subsets C and D represent the set of conditional attributes and the set of decision attributes, respectively, D = ∅. V r denotes the nonempty set of attribute values for attribute r ∈ R, i.e., the value domain of attribute r, V = ∪ r∈R V r . f : U × R → V , is an information function that specifies the property value of each object x in U.
(2) Upward union and downward union The derivation of the decision rules requires determining the classification level of the decision attributes D and then selecting the approximate sets of the upward union and downward union of classes as the input of the derivation process according to the trends of the objective functions. The set of decision attributes D can divide U into a finite set of classes Cl, Cl = {Cl t , t ∈ T}, T = {1, . . . , n}. For ∀x ∈ U, x can only belong to one class Cl t . The upward union and downward union of classes can be defined as Cl ≥ t = ∪ s≥t Cl s and Cl ≤ t = ∪ s≤t Cl s . As for ∀r ∈ T, if r > s, the objects in Cl r are better than objects in Cl s .

(3) Dominance relation and approximations
The approximations of the upward or downward unions are associated with the dominance relationships between objects. If an instance is described as x P-dominates y with respect to P, then it means x q y for all q ∈ P, which is denoted by xD P y. As for P ⊆ C, ∀x, y ∈ U, The P-dominating set and P-dominated set of x are D + P (x) = {y ∈ U|yD P x} and D − P (x) = {y ∈ U|xD P y}. The P-lower approximation of the subset Cl ≤ t and Cl ≥ t can be defined as (16) and (17): The P-upper approximations are (18) and (19): The P-boundary of Cl ≥ t and Cl ≤ t can be defined as (20) and (21): (4) Decision rule Based on the above concepts, decision rules are further introduced. The form of decision rule r is "If φ then ψ", where φ is the conditional part and ψ for the decision part. The φ is the ensemble of conditional value pairs, where the values assumed by conditional criteria are specified. The decision part specifies the assignment of one or more decision classes. When only one attribute is considered in the decision set, two types of decision rules can be used as follows.

•
Certain D ≥− decision rules, providing lower profiles of objects belonging to P Cl ≥ t : • Certain D ≤− decision rules, providing upper profiles of objects belonging to P Cl ≤ t : i f f (x, q 1 ) ≤ r q 1 and f (x, q 2 ) ≤ r q 2 and . . . f (x, q n ) ≤ r q n , then x ∈ Cl ≤ t , where {q 1 , q 2 , . . . , q n } = P c ⊆ C.

(5) Quality of approximation
With given data sets, we can derive a set of decision rules using the DRSA method. Then the credibility levels of these decision rules can be measured using cr. The formulation is shown in Equation (22), where the operator "{ . . . }" denotes the set of objects that satisfy the internal conditions, and "| . . . |" means the number of objects in the internal set.
Based on the above knowledge, these needed steps for using the DRSA method to solve the problem in this study are as follows and presented in Figure 5. Step 1. Obtain the information for a decision table based on a representative sample of Pareto front output by the chosen algorithm. The table includes solutions and all the objective values of each solution.
Step 2. Use all the objectives as conditional attributes and invite experts as DMs to analyze the table data. They are required to classify each object in the sample solutions as "better" or "otherwise", corresponding to Cl 1 and Cl 2 .
Step 3. Obtain the P-lower approximate of Cl 1 according to the above decision table, P ⊆ C. Then, derive a set of decision rules for P Cl ≤ t and calculate the values of cr. The obtained rules are then shown for DMs in the interactive process.

Case Study
Based on the demand of S shipyard in Shanghai, China, for importing large amounts and multiple categories of shipbuilding materials from Italy, this study selected Genoa and Shanghai as the origin and the destination of the routing problem, respectively. To know the impact of different needs on the choice of transportation services and to verify the validity of the proposed approach better, five freight batches were considered, and every batch contained a 40-foot equivalent unit or equivalent cargo that was converted to cubic meters for air transport. The range of operation plan was one month, and the different freight information of batches that were considered included the release time at the origin, delivery time windows, and the potential requirements for cost and service quality. To exclude the impact of force majeure of the epidemic on the transportation industry, the paper focused on operation data of the constructed transportation network in 2018 and took August as an example. Table 2 shows the time information of the freight batches, where the figures are obtained by converting the time unit hour, 0:00 on 1 August in 2018 represents 0 and 1 h is 1 unit.
After analyzing the transportation channels between the two cities, international rail, ocean shipping, and air transport were taken as the main international transportation modes and several representative nodes for every mode were selected. Considering road and rail networks are well-developed in Europe and the Chinese mainland, we selected those as modes connecting the international transport services to the origin and the destination. Inland river transport on the Yangtze River was also considered between some cities. Therefore, a transportation network consisting of 19 nodes and 312 inter-node transport services was constructed as Figure 6. The numbers on the directed arcs represent transportation services, such as "20-21" on the arc between Genoa and Rome means that there are two transportation services between the two nodes, numbered 20 and 21, respectively. Due to the difficulty of obtaining all the services' data, only international transportation services were considered as scheduled modes in this study. Considering that efficiency bottlenecks in intermodal transportation plans usually exist in the international segment [34], two parameters for flexibility objective, i.e., the service frequency and the latest booking time were also obtained from the service data of international transport. As for the time uncertainty, we set the variable operation time at each node to obey uniform distribution and inter-node travel time to obey triangular distribution. The allowable maximum value for the fluctuation of the expected transportation time for reliability objective was 12 in all batches.
The collected data contained information of services at nodes and inter-node transportation services. Specifically, the operating schedules of international transportation modes were obtained from several international freight companies including China Railway Container, COSCO Shipping Lines, Air China, and China Eastern. The transit time data for the road and rail services were estimated from the ratios of miles to speed (km/h). The road and rail networks were downloaded from the Open Street Map (OSM: https://www.openhistoricalmap.org, accessed on 26 August 2021) website. The information that was related to the freight rates and other transit time were provided by the Jincheng Logistics website (shipping.jctrans.com, accessed on 20 August 2021). The operating time information at the nodes was obtained by consulting with the relevant freight companies. For convenience, we numbered the nodes and recorded those in Table 3. Besides, considering the important impact of international transport sectors on the transportation schemes and the apparent differences in costs and service levels between the modes, this study divided all the route types into the following three categories: air-mode, sea-mode, and international rail-mode services.  The model that was constructed was solved by the proposed interactive approach. The running of NSGA-III was implemented by Python and an experienced manager that was in charge of procurement at S shipyard was invited to be the DM. With the classification of the solutions from DM, jMAF software (http://ww.cs.put.poznan.pl/jblaszczynski/Site/ faq.html, 18 November in 2021) was applied to derive the decision rules that were based on DRSA.
After the first running, the results showed that Batch 3 had more feasible solutions and multiple service types, so this batch was used as an example to describe a complete interactive solving process. In this study, we set the initial population size N = 100, the maximum number of iterations maxGen = 500.
This problem is a discrete problem with constraints and a limited data size, the number of Pareto solutions that were generated by the iterations was small, so it is not suitable for selecting the representative solution sets from Pareto fronts by methods that are suitable for large-scale problems, such as cluster analysis. Thus, this study used a simple randomized approach to determine the schemes that were shown to the DM. First, we ran NSGA-III to generate the Pareto front without preference for Batch 3, which contained 18 solutions. Then, 15 schemes were randomly selected and presented to the DM. We also provided the domain of objectives values of solutions and used a parallel coordinate plot to visualize the sample distribution to help the DM make decisions, see Figure 7. The DM classified these schemes into two categories: "better" and "others". The evaluation of the given set from the DM is shown in Table 4, where the figures inside the "[]" represent the transportation service numbers. After obtaining the DM's preferences, two decision rules that matched the options that were evaluated as "better" were output by the DRSA method, as shown below: Rule1.1: if C 3 (l) ≤ 6346.12, R 3 (l) ≥ 0.74, then the route l is "better", l ∈ L. Rule1.2: if T 3 (l) ≥ 1, R 3 (l) ≥ 0.74, then the route l is "better", l ∈ L.
The credibility cr of the above two decision rules were 1 and the DM preferred rule 1.1. Therefore, the following constraints were added to the original optimization problem. C ≤ 6346.12, R ≥ 0.74 (23) The above discussion showed that the DM preferred transportation services with lower transport costs and higher reliability. Based on this rule, the routes satisfying the constraints mainly involved international rail or sea transport. Next, we utilized the rule to generate the second sample of non-dominated solutions and presented it to the DM again. The set contained seven feasible solutions and, according to the design, these still needed to be classified by the DM. The evaluation is shown in Table 5. The above preference information was applied to the DRSA method, and we can get the following decision rule. The cr of that was 1.
Rule 2.1: if C 3 (l) ≤ 5803.06, R 3 (l) ≥ 0.84, then the route l is "better". The rule was set as the default option for the DM. In this interaction, it can be seen that the DM further made requirements for reliability and cost. The rule was then transformed into a constraint of the optimization problem, and the third Pareto set was obtained in Table 6. For this time, the number of solutions was less than five, so the DM could select an optimal solution directly. The route belongs to rail-international rail-track combination scheme. The delivery time of Batch 3 was the longest one. Within the set time window, the services with air, sea, or international rail are all competitive. However, the results of the two interactive processes indicate that the lower cost and higher reliability are necessary. The international rail-mode plans have a cost advantage compared to air-mode services and a higher reliability than sea-mode services, which made the solving enter into the third interactive round. In this process, the DM eventually chose the scheme with a low cost but less reliability, which indicates that the DM values the cost objective more when reliability reaches a certain level. We can also see that the alternative routes with international rail have high flexibility, although it was not mentioned in the decision rules.

Results
Similarly, the proposed interactive method was applied to solve the routing problems for the other four batches of shipments. After the first round of algorithm running, the Pareto sets of all batches and the distribution can be shown in Figure 8.
First, it can be seen that reliability and flexibility have a positive relationship, and the cost objective has certain conflicting relationships with the reliability and flexibility. Transportation services with higher reliability also tend to have higher flexibility values. But when the cost objective is optimal, the flexibility and reliability are generally less. Further analysis shows that this situation is significantly related to the types of transportation services. In Batches 2-5 with multiple service types, the air-mode services have obvious advantages in reliability and flexibility objectives but are costly. The sea-mode plans have the opposite characteristics, while the international rail-mode services fall in between. Second, the timeliness is not obviously related to the other objectives and all the service types are likely to have high timeliness. It is because the ranges of shippers' satisfaction time windows of batches are different. Under the shorter delivery time requirements such as Batch 1 and 2, the air-mode services have the highest timeliness. Other service types achieve high timeliness when the specified delivery time is extended or the range of that is expanded.
Then, the subsequent steps were performed for the remaining batches according to the designed process. The information, including the number of non-dominated solutions that were obtained by algorithm iterations, representative decision rules (rdr), and the preferred solutions, are summarized in Tables 7 and 8. A and B in Table 8 denote the optimal vector and the worst vector of all the objectives values of the Pareto set for each batch, respectively.
It can be seen that for Batch 1, the alternative air-mode services have high timeliness and reliability. The DM emphasized more flexibility when expressing preferences. This may be because the DM's requirement for this batch was more urgent, and the specific demand time was uncertain, so on-time delivery and highly flexible services were needed. The airmode services are the best choice for this type of demand. In Batch 3, there are many types of transportation services that are available, the shipper emphasized both cost and reliability objectives and gradually increased the constraints in two interactions. An international rail-mode service plan was chosen lastly. For Batch 5 with the largest permissible delivery time range, the DM first proposed the requirement for cost and timeliness, then reinforced the cost reduction constraints. Finally, the DM obtained the least costly "sea" scheme. In the above three batches of freight, the air-mode, international rail-mode, and sea-mode services showed their respective advantages. For Batch 2, the DM imposed high timeliness and flexibility requirements, so the air-mode services were preferred again. As the constraint was further strengthened, the DM also hoped to reduce the cost. The above processes show that the DM firstly pursued timely delivery and more flexibility in service time, secondly lower costs.   For the longest delivery time request of Batch 4, the alternatives generally have high timeliness. After two interactions, the DM's requirement for reliability reached one and finally chose the rail-international rail-rail plan. It may be because the DM proposes transportation plans in advance and the importance of the reliability of delivery time. The satisfaction of requests of these two batches in the interactive process shows the rich preference expression process of the DM. The DM guided the convergence of feasible non-dominated solution sets by expressing preference information during the interaction process and after several interactions, the most satisfactory transportation routes for the different demands were finally determined.
Let us compare the solutions that were obtained by the proposed interactive method with the solutions that were given by applying a linear-value-based priori method and a posterior method, which was based on the technique for order preference by similarity to an ideal solution. The weights of the objectives were set based on the preference information of the DM for five batches at the first interaction, as shown in Table 9. Table 9. The weights setting of the objectives in the adopted prior and posterior methods.

Batch No.
Weights of Objectives C T R F A comparison of the final solutions that were obtained by the three methods is shown in Figure 9. First, besides high timeliness and flexibility commonly for Batch 2, the plan that was given by the interactive method also has a lower cost. This is because the applied prior and posterior methods can only obtain the DM's preference information once and cannot have any effect on the selection process if the DM hopes to supplement information. The final solution of the interactive method for Batch 3 is more costly but more reliable. Further analysis reveals that the DM made adjustments to the requirements of the cost and reliability values during the two interactions. The restriction and adjustment on the objective values from the DM are hard to reflect in the weight.
For Batch 1, the final solution of the posterior approach has less flexibility but a cost advantage than the result of the priori approach. This illustrates that the ability to extract weight information of different methods differs. This situation will further affect the complete expression of the DM's preference information.
Through solving the multi-objective optimization of IFRP for the five batches of freights using the proposed interactive method, it can be seen that the DM can avoid the pressure of comparing the objective weights and freely emphasize single or multiple objectives based on their needs. Moreover, it allowed the DM to strengthen, supplement, or even modify their preferences through multiple interactions.

Conclusions
In the rapidly evolving global business environment, understanding customer preferences and requirements becomes one of the key factors to the success of any transportation service [47]. In this study, we focused on the multi-objective modeling and optimization of a customer-centered intermodal transportation routing problem. We modelled the cost, timeliness, flexibility, and reliability as optimization objectives, and propose an interactive multi-objective optimization approach that is based on the combination of NSGA-III and DRSA to help shippers effectively find the most preferred solutions.
In the case study, we used the solution process and the method comparison for five freight batches to illustrate the impact of different requests on transportation service selection, and to verify the validity and rationality of the proposed method. The following findings can be deduced.
(1) The reliability objective and flexibility objective have a positive relationship and the cost objective has negative relationships with reliability and flexibility. This result greatly relates to the quantified methods of objectives in this study and the types of transportation services. The values of reliability and flexibility are determined by the latest booking time, service frequency, and the fluctuation range of the total transportation time. Compared with other service types, air-mode plans have the shortest service time and are booked flexibly. As such, these services perform better in flexibility and reliability but at a high cost. In contrast, sea-mode services have low transportation costs but the lowest reliability and flexibility, while international rail-mode services fall somewhere in between.
(2) The timeliness objective is not sensitive to the variations of the other objectives. The result differs from the finding that was observed by Sun and Li [19], which considers a conflicting relationship between the timeliness and economic objectives. A possible explanation may be related to the comparing basis of the two objectives. The objective values except for timeliness of Pareto solutions that were compared in this study are greatly influenced by service types, and there may be multiple types of schemes with optimal timeliness within one satisfaction time interval. The latter compared solutions under the only road-rail transportation and the relationships between the objectives were given only relating to the applied data in the case.
(3) Compared to other approaches that were applied in IFRP, the proposed interactive approach reduces the decision burden on the shippers, allowing them to be more engaged and better able to express their true expectations so that they can obtain satisfactory results with higher effectiveness.
Specifically, this interactive approach can handle the expected trade-offs between the multiple requirements of real decision-making scenarios. Shippers can avoid the pressure of comparing goal weights and achieve learning during the interaction process. More importantly, this approach provides shippers a chance to freely emphasize certain goals and reinforce, complement, or even modify their goal preferences, thus helping shippers to fully express their preferences.
Theoretically, this study provides a reference for IFRP research on how to make freight decisions around customer needs. The proposed optimization model and interactive approach enrich the multi-objective planning research of IFRP. Practically, it is expected to help stakeholders to achieve a balance between the cost and the service level more systematically. It is beneficial for intermodal operators to improve services to customers and to better understand the market demand so that they can optimize and promote sustainable development. Additionally, it can also help shippers to enhance the flexibility in choosing transportation services and procure services that better meet their interests, thereby reducing the adverse effects of irrational plans, such as higher inventory costs or production delays.
There are some limitations to this study. The modeling of reliability considers only the stochastic uncertainty of travel time and variable operation time at nodes. The modeling of flexibility focuses on the responsiveness of demand for transportation time changes before transportation. Both the reliability and flexibility of transportation include other aspects, such as network disruptions for reliability and changes of demand volume during transportation for flexibility. Due to data limitations, this study only covers a portion of these factors. Therefore, future efforts could be made to further refine these two objectives. Second, this study pays less attention to the selection of the DMs and the performance evaluation of the interactive approach. It is worthwhile to combine the method with group decisions and conduct user experience validation as an alternative evaluation approach, to further extend the application of this approach and demonstrate the user participation benefits in it.
Author Contributions: Conceptualization, methodology, software, data curation, writing-original draft preparation, C.S.; writing-review and editing, funding acquisition, H.W. and M.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement: No applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.