An optimization approach for green tourist trip design

In this paper, the Multi-Objective Multi-Modal Green Tourist Trip Design Problem (MO-MM-GTTDP) as the multi-modal variant of the orienteering problem is investigated. For this problem, a mixed-integer linear model is formulated, containing three objectives: (1) maximizing the total score of the trip, (2) minimizing the total cost of the trip, and, (3) minimizing the total CO2 emission produced by transportations in the trip. Various transportation modes are considered for the tourist to choose to move between points of interest (POIs). The tourist choice may be affected by the transportation time and cost. Moreover, choosing the transportation mode will have an impact on the amount of trip pollutants. The cost of visiting POIs, as well as the cost of transportation between POIs, is considered as the total cost of the tour. The proposed problem is solved for two types of tourists: (1) tourists who give the highest priority to the total score, then the total cost, and at last, the total emission; (2) tourists who would like to compromise between objective functions and select a solution which best suits their situation. For the former type, the problem is solved using lexicographic method and one best solution is proposed. For the latter type, a ε-constraint\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upvarepsilon -\mathrm{constraint}$$\end{document} method is implemented which provides Pareto optimal solutions. In addition, a Multi-Objective Variable Neighborhood Search (MOVNS) algorithm is designed to solve instances of this problem. The performance of this method is evaluated by comparing its results with the ε-constraint\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upvarepsilon -\mathrm{constraint}$$\end{document} and the lexicographic methods. New instances of the problem are generated based on the OP existed benchmark instances. Moreover, the MOVNS is also compared against two other metaheuristic versions of the solution approach using multiple metrics. The conclusion is the high quality of the proposed MOVNS algorithm solutions in practically acceptable computation time (few seconds). Finally, a small case study based on real data on several POIs in the city of Tehran is generated and used to demonstrate the performance of the proposed model and algorithm in practice. For this case study, by using the multi-attribute decision-making method of TOPSIS, the obtained non-dominated solutions are ranked, and the best ones are presented.


Introduction
Nowadays, Tourism has become as one of the significant and even crucial sources of foreign income in a lot of countries, particularly in the context of developing economies (Gautam 2012;Cucculelli and Goffi 2016). This industry has a vital role in economic development, involves multiple sectors, and interlinks with several industries (Agaraj and Murati 2009). Imagine a tourist who wants to visit a touristic city with many attractions, in a couple of hours. According to the tourist's preferences, each POI has a degree of attractiveness and there is a specific visiting cost for each of them. In addition, there are different transportation facility modes with different transferring costs and time that can be selected by tourists to move between POIs. In this case, there are three aspects to consider: the potential utility as well as the cost of visiting each point, the travel cost and the duration of transporting between POIs, and the amount of carbon dioxide released when traversing between POIs (Divsalar et al. 2017).
In general, personalized tourist trip planning is a complicated and time-consuming process which consists in selecting points of interest (POIs) and scheduling of trips (Liao and Zheng 2018), giving a suitable plan to tourists on the sequence of POIs to explore conveniently and costeffectively (Wong , 2011;Leiper 1990). However, tourists often do not have enough time to visit all of the POIs during their day tour (Tsai and Chung 2012). This problem has been known as the ''Tourist Trip Design Problem'' (TTDP). TTDP includes tour route planning for tourists, and maximizing their desirability while considering existing constraints (Vansteenwegen and Oudheusden 2007), such as selecting of POIs they feel are the most appropriate for them ).
Over the past several years, extensive research has been generated on the content of TTDP (Liao and Zheng 2018;Tsai and Chung 2012;Hsu et al. 2012;Lee et al. 2009;Liu et al. 2014;Rodríguez et al. 2012;Zheng et al. 2017). Despite recent achievements, research on the domain of tourist personalized trip planning is still at early stages. By analyzing the existing studies, we find that most of researches have investigated the TTDP without considering the environmental aspects of the tour. However, the unplanned development of this section can have destructive effects on natural resources and local societies. Moreover, tourism is considered as a significant contributor to climate change by making greenhouse gas emission (Becken 2004). In a very recent report by the United Nations World Tourism Organization (UNWTO), mentioned tourism as an influential factor in climate change and global warming's phenomenon providing evidence of the CO 2 emissions from tourism and the implications of the different modes of transport (World Tourism Organization and International Transport Forum 2019). The annual movement of billions of travelers often over large distances made the contribution of the tourism sector to global climate change at approximately 5% of CO 2 emissions. Moreover, tourism's future forecasts show that in all scenarios, it will grow substantially and will account for an increasingly large share of global greenhouse gas emissions, particularly if other sectors manage to achieve absolute emission reductions (Gössling et al. 2013). More specifically, it is shown that tourist arrival has a significant impact on the CO 2 emission of the transportation sector in many countries (Al-Mulali et al. 2015;Ng et al. 2016). These findings indicate that next to governments, and tourism sectors, tourists themselves, pay attention to the environmental quality of the country that they travel (Al-Mulali et al. 2015;Ng et al. 2016). Therefore, it is crucial that any touristic plan should consider the environmental aspects as well. In a TTDP, this subject requires a useful trade-off between maximizing tourist entertainment and minimizing the produced pollution of the tour. Thus, the resulted personalized trip planning problem can be treated as a multi-objective (M-O) problem, which has attracted little attention to date.
Although Vincent et al. (Vincent et al. 2017) made pioneering efforts on this issue from the view of the multimodality of team orienteering problem, their study does not mention some relevant points: First, their research concentrate on attaining the benefits of tourists as uniformly as possible, but it neglects minimizing the cost as well as produced pollutants of the tour. Second, this study considered the choice of POIs, but dismissed sequencing the transportation facilities between POIs, which may depend on the cost and the amount of produced pollutant. Except for this study, the concern of tourism's roles in climate change is nearly neglected in the investigated variants of TTDP. Moreover, to our best understanding, in all the presented M-O variants of the OP there is a trade-off for achieving the maximum utility of the tour by different categories of POIs (e.g., leisure, history, culture, shopping, etc.), and different score for each category. In fact, the cost of the tour has been considered as a constraint.
Knowing the restrictions of previous studies attend to the TTDP, the orienteering problem (OP) and its variants which have been applied to model various versions of the TTDP (Feillet et al. 2005;Gunawan et al. 2016) are investigated in this research. The OP is a single-criterion variant of traveling salesman problem with profits (TSPP) in which a set of vertices with a determined score is given. For all existing pairs of vertices, there is a required travel time. Due to the time limitation, not all vertices can be visited. Besides, the first and final vertices have to be visited, and other vertices can be visited at most once (Divsalar et al. 2013). Next to the tourist trip planning, there are other practical applications mentioned in the literature for the OP and its variants, including, technicians routing, athlete recruitment, military applications, the mobile crowdsourcing problem, the smuggler search problem, the wildfire routing problem, and the integration of vehicle routing, and customer selection problems (Vansteenwegen and Gunawan 2019). In almost any of these applications, minimizing the cost and environmental emissions can be seen as other objectives of the problem.
To conceal the gaps, the present work aims to develop a personalized tourist tour by considering the tour cost as well as the environmental aspects in a multi-modal environment. This problem is complicated because of considering three potentially conflicting objectives and multiple limitations. To dominate the complexity involved, we handle three objectives by using the definition of Pareto front optimality and developing an efficient M-O metaheuristic approach, namely the MOVNS. This problem is solved for two types of tourists: 1) tourists who give the highest priority to the total score, then the total cost, and at last, the total emission; 2) tourists who would like to compromise between objective functions and select a solution which best suits their situation. For the former type, the problem is solved using lexicographic method and one best solution is proposed. For the latter type, a e À constraint method is implemented which provides Pareto optimal solutions. The MOVNS is also able to provide one best lexicographically selected solution. The performance of the presented MOVNS is evaluated according to the M-O evaluating metrics by comparing its results with the results obtained by the implemented epsilon-constraint method, as well as the exact lexicographic method. In addition, a small real example over a number of POIs in the city of Tehran (IRAN) is generated and used to demonstrate the performance of the proposed model and solution algorithms in practice. For this case study, by using the multi-criteria decision-making method of TOPSIS, 1 the obtained non-dominated solutions are ranked, and the best ones are presented to the tourist.
The rest of this study is structured as follows: Sect. 2 manages a review on the literature of TTDP and its variants. Section 3 introduces a mathematical modeling for the MO-MM-GTTDP. Section 4 describes the proposed solution method and presents all its phases in detail. Numerical experiments are performed in Sect. 5, including sensitivity analyses of the model and solution methods. The case study of Tehran is presented and discussed in Sect. 6. Finally, this study is summarized and conclusions and potential developments are provided in Sect. 7.

Literature review
According to Gavalas et al. (Gavalas et al. 2014), variants of tourist trip design problem can be categorized as singletour and multiple-tour TTDP. The majority of published researches in the literature of TTDP model the problem based on OP and its variants ). The OP is categorized as a variant of TSPP which is a biobjective problem. In the TSPP, the time budget constraint of the OP is also considered as an objective function. Therefore, the goals in the TSPP are maximizing the collected profit while the travel cost is minimized (Gavalas et al. 2014). Bérubé et al. (Bérubé et al. 2009) present the first exact Pareto solutions for TSPP by using the e À constraint method in 2009. At the same year, the M-O Orienteering Problem (MOOP) is introduced by Schilde et al. (Schilde et al. 2009) as a variant of the OP, in which each POI can be in a different category (e.g., history, culture, shopping, leisure) and a different score has been considered for each category. They presented an adaptation of the Pareto Ant Colony Optimization (P-ACO) (Doerner et al. 2004) metaheuristic and a M-O extension of Variable Neighborhood Search (P-VNS) introduced by Chao et al.  for the first time. TTDP with clustered points of interest is a recent extension of this problem, where clusters representing different types of attraction sites. A fuzzy greedy randomized adaptive search procedure is used to select the candidates in the construction phase. Computational results show the effectiveness and efficiency of their proposed approach and its suitability to be part of a personalized electronic tourist guide devices (Expósito et al. 2019). In 2015, Yu-Han Chen et al. (Chen et al. 2015) proposed an Ant Colony algorithm for MOOP with Time Windows (MOOPTW) to find a set of nondominated solutions. In the MOOPTW, there is a different specified time interval for visiting each POI. An artificial bee colony algorithm is also presented for the bi-objective MOOP by Martin-Moreno & Vega-Rodrigues (2018) (Martin-Moreno and Vega-Rodriguez 2018). They compared their method with the P-ACO and P-VNS of Schilde et al. (Schilde et al. 2009) and proved a higher-quality performance over the benchmark and real-word instances. Recently, Zhenga et al. (Zheng and Liao 2019) have introduced a modified version of MOOP, which considers a compromise between the fairness of individual members and the total obtained score of the group. In this problem, a personalized trip is proposed for a heterogeneous group of tourists modeled as a common M-O problem. In this study, an ant colony optimization heuristic is combined with the differential evolution algorithm as the solution approach. The Time Windows version of the OP (OPTW) is introduced by Marisa et al. in 1992 (Kantor andRosenwein 1992), and for the first time is used to model the TTDP by Vansteenwegen et al. (Vansteenwegen et al. 2011). They argued that time windows can have a significant impact on the regular OP. Another variant of the OP is the Time Dependent OP (TDOP) proposed by Fomin and Lingas (Fomin and Lingas 2002). In the TDOP, the cost of moving between POIs depends on the time that this travel is started. Therefore, the TDOP has been suggested to model the Multi-Modal cases. In 2016, two metaheuristic algorithms are presented for solving the M-O TDOP by Mei et al. (Mei et al. 2016): A M-O Memetic Algorithm (MOMA) and a M-O Ant Colony System (MACS). Liao et al. in 2018 (Liao andZheng 2018) introduce a time-dependent stochastic TDOP to design a personalized route where travel times or wait times are not static. This study designs a hybrid heuristic algorithm-based on random simulation (RSH2A). They have evaluated their solution method by a case study in China. Another bi-objective variant of the OP is recently introduced by Dutta et al. (Dutta et al. 2020), where they maximize the visited clusters' profits as well as the customers' satisfaction in the open Set Orienteering Problem (SOP). The SOP is a clustered version of the OP, which can be applied to model the TTDP (Archetti et al. 2018).
An interesting multi-tour variant of the OP is the Orienteering Problem with Hotel Selection (OPHS), which is introduced by Divsalar et al. (Divsalar et al. 2013;Divsalar et al. 2014). In the OPHS, they considered a set of attractions with an associated profit and a set of non-rated hotels, with the goal of determining a number of connected trips to visit a some of available attractions and maximizing the total collected profits. They proposed two methods to solve the OPHS, a VNS, and a memetic algorithm. Another multi-tour extension of the OP is the Team Orienteering Problem (TOP), in which the goal is to find K tours (or paths) in a limited time to maximize the total score given that each nonterminal node can be visited at most once in all K tours. The TOP is introduced by I-Ming et al.  and is used to model the TTDP variant by ). Lately, Tsakirakis et al. (Tsakirakis et al. 2019) introduced the Similarity Hybrid Harmony Search (SHHS) algorithm as a new variant of Harmony Search (HS) for TOP. They proposed the static and dynamic versions of the method. In the static version, the values for the parameters are predefined, and in the dynamic one, parameters are adjusted during the solve procedure. The results of applying the proposed method for benchmark instances of TOP show the superiority of the dynamic version.
In a research by Gunawan et al. (Gunawan et al. 2017), an Iterated Local Search as well as a hybrid Simulated Annealing and Iterated Local Search (SAILS) is presented for the TOP with Time Windows (TOPTW). They tested their proposed methods on benchmark TOPTW instances and concluded from the results that the proposed ILS and SAILS improves the solutions in terms of quality. In 2018, Hu et al. (Hu et al. 2018) introduced the M-O TOP with Time Windows (MO-TOPTW), where multiple profits for each point are considered. They developed a M-O evolutionary algorithm on the basis of a decomposition technique and constraint programming to solve this problem. The results show an improvement by comparing the results of applying benchmark instances with best-known solutions. Another version of MO-TOPTW is introduced by Hapsari et al. (Hapsari et al. 2019), which is different in objective functions. In this research besides the maximizing of the scores, the available time of the tour is considered as a second objective instead of a constraint. The authors used an Adjusted ILS (AILS) to solve the proposed model on a real case and compared it with some other metaheuristics. The result demonstrates the higher performance of their method in terms of computational time and providing itineraries with higher total score. In 2020, Aghdasi et al. (2019) proposed a new variant of MO-TOPTW for rescue applications with five objectives. They present an efficient solution by combining M-O evolutionary algorithms (MOEAs) with learning algorithms. The results show that their algorithm decreases the gap with best-known solutions of TOPTW by up to 42% In 2017, the Multi-Modal Team Orienteering Problem (MM-TOPTW) as new functional variant of the TOPTW was introduced by Vincent et al. (Vincent et al. 2017). In this research, several transportation modes are considered for tourists in the TTDP. The goal is to achieve the maximum value of score without violating the travel cost and time limitations. In addition, time windows are considered for the opening and closing time of each POI. They proposed a two-level particle swarm optimization with multiple social learning terms for the MM-TOPTW. The authors implemented their proposed algorithm for solving 56 benchmark VRPTW instances. The results showed promising solutions within an acceptable computation time.
Due to the wideness of the existing OP variants and applications, for an all-inclusive literature review, we refer the readers to the fairy recent surveys on the OP variants by Gunawan et al. (Gunawan et al. 2016) and Vansteenwegen et al. (Vansteenwegen et al. 2011).
With an overall look to the presented literature, very few number of research among the OP variants and its extensions discuss the M-O and the multi-modal variant of the orienteering problem or, more specifically, the TTDP. A detailed comparison of the mentioned M-O versions of OP variants are presented in Table 1. This table summarizes the difference in the various objective functions of each M-O OP as well as their proposed solution methods. These objectives can be divided in two main types: The ones which focus on the node selection, and the ones which depend on the arc selection. The current study is among the very few ones which consider both node-and arc-based objective functions simultaneously. In addition to the tourist's preferences and the total cost of the tour, that contains the visiting and the transportation cost, giving a choice to the tourist to choose his/her transportation facility to move between POIs can be considered as an effective factor in the problem. Another factor, which has not been noticed in the literature, is the environmental aspect of the tour. As emphasized in Sect. 1, due to the major effect of tourism on air pollution, and climate change of the touristic cities, it is certainly worth to consider this factor in a TTDP. Accordingly, we consider tourist preferences, the cost, and the amount of CO 2 emission simultaneously, and introduce a new variant of TTDP, which we call the M-O Multi-Modal Green TTDP (MO-MM-GTTDP). As a result, because of the compromise between three objectives, the MO-MM-GTTDP is a complex problem to solve. To handle these conflicting objectives, a Pareto-based local search method is proposed to achieve diversified and realistic routes for this problem based on an effective M-O optimization approach. This is among the first works in which a MOVNS is developed for an OP variant. The advantage of using such a structure is its flexibility in using the power of local moves in tackling various objectives simultaneously. Furthermore, we consider two types of tourists as the potential users of this method. The first type has a fixed rank order of the objective functions, and the second type is interested in receiving a list of indifferent solutions. In this research, for the first type of tourists the lexicographic method is used to provide one best solution. Motivated by the fact that the user will have a fuller knowledge and better understanding of the solution space when he/she knows the set of all non-dominated solutions (Haimes and Li 1988), for the second type of tourists, the MOVNS is used to find a list of Pareto solutions. However, it is always beneficial to help the user to rank and clarify the best route. To this end, in our case study, the TOPSIS method is applied to the obtained non-dominated solutions from a real example, made over real data from the city of Tehran.

Problem description and mathematical formulation
In this section, the Multi-Objective Multi-Modal Green Tourist Trip Design Problem (MO-MM-GTTDP) is formally formulated as a mixed-integer linear problem. Consider G = (V, E), as an undirected graph. In this graph, V is the set of N vertices, and E is the set of its edges. We consider the POIs as the vertices of the graph G from 1…N. Each POI is associated a score S i based on the tourist's preferences as well as a visit cost showed by VC i . K is the set of types of transportation facilities to move between POIs. Each transportation mode produces a certain amount of emission per distance measure. P ijk represents the amount of carbon dioxide emission of traveling from i to j using the transportation mode k, and C ijk is the cost of this transport. t ijk is the travel time of passing edge ði; j; kÞ 2 E. The start and the end node of the visit tour are fixed to nodes 1 and N, respectively. Each node is visited at most once. X ijk is the binary variable set to 1 if edge ði; jÞ is included the solution, and traversed by transportation mode type k 2 K; 0 otherwise. Then, the mathematical model is presented using Eqs.
(1) to (11). d þ S ð Þ ¼ fði:jÞ 2 E : i 2 S:j 6 2 S} shows all the outgoing edges from set S and for the simplicity d þ fig ð Þ is presented as d þ i ð Þ. In a similar way, d À S ð Þ ¼ fði; jÞ 2 E : i 6 2 S; j 2 S shows all the incoming edges to set S. Then, using the graph notation, the mixed-integer model is presented as follows: Equation (1) maximizes the total profit of visited POIs. Equation (2) minimizes the total cost of the tour, including the cost of visiting POIs and the cost of moving between POIs. Equation (3) minimizes the carbon dioxide emissions produced by transportation facilities. Constraint (4) limits the total travel time within the time budget T max . Constraints (5) ensure that in each path between i and j, only one transportation mode is selected. Constraints (6) & (7) ensure that node 1 is the start point, and node N is the end point of the tour. In a similar way, constraint (8) is the flow constraint and ensures that each node is visited at most once. Finally, constraints (9) and (10) ensure that there are no sub-tours. The values of decision variables are determined by constraints (11).

Solution strategies
We assume that there are two types of tourists regarding the priorities on the objective functions. The first type gives the highest priority to maximize its happiness, then, next priority to minimize its cost, and the least priority to minimize the emission. Taking these priorities into consideration, we have used the lexicographical method to solve instances of the MOGTTDP. The second type of tourists would like to benefit from having a list of Pareto solutions.
In the following sections first, the lexicographical method is described, and then, the epsilon-constraint method as one of the most reliable M-O methods is implemented to provide Pareto solutions using the presented mathematical model. Moreover, the obtained results from this method are used to evaluate the proposed metaheuristic. Next to the epsilon-constraint method, an efficient metaheuristic is developed to solve the MO-MM-GTTDP instances, and to obtain Pareto solutions. This method is the M-O version of the VNS method (MOVNS) and is based on combining several local moves within a systematic change procedure, avoiding to trap in a local optimum. Because of considering three objective functions simultaneously, designing special moves to focus on improving each objective and the transportation facility selection, the MO-MM-GTTDP is more challenging than the OP. Thus, significant computational efforts are assigned to evaluate combinations of vertices and the selected transportation facility to move between them. We have also used an optional module as the last step of the presented MOVNS so that the best solution according to the lexicographical priorities of objective functions can be provided to the tourist of the first type.

Lexicographic optimization method
The lexicographic method is a sequential elimination method in which the decision maker is first asked to rank order the objectives in terms of their importance (Chankong and Haimes 2008). In the case of TTDP, we assume that the first type of tourists gives the higher rank to maximize total score, the second rank to minimize total cost, and the third rank to minimize total emission. Therefore, in our implementation, the proposed model is solved in three steps. In the first step, the model is optimized, while the second and third objective functions (Eqs. (2) and (3)) are ignored. Thereby, the optimal solution is achieved where the maximum possible score (maxZ1) is obtained. In the second step, the first objective is imposed as the following constraint to the model, and the second objective is optimized.
The optimal solution obtained at this stage (minZ2) is then used as the right-hand side of the following constraint. This constraint is also added to the model of the previous step. Then, Eq. (3) is used as the objective function and the resulted model is solved to optimality.
This way, we end up with one best solution for the tourist of the first type. Moreover, an extra module has been added to the MOVND approach to lexicographically sort the found Pareto list and give one best solution as an output. This module is explained in Sect. 4.3.6.

Epsilon-constraint method
To have a benchmark for evaluating the performance of the proposed metaheuristic, first an epsilon-constraint approach is implemented to solve the instances of the MO-MM-GTTDP. Therefore, in this section this method is briefly introduced.
Assume the following M-O problem: .., f m x ð Þ are the m objective functions, and D is the feasible region. A M-O problem usually consists of a set of solutions called the Pareto optimal front, where each Pareto optimal solution demonstrates a trade-off between different objectives. In fact, it is not possible to improve all objectives functions at the same time. In M-O optimization rather than the single-objective, comparing solutions is more complex. Two subjects have to be investigated in M-O optimizations: Pareto dominance and Pareto optimality (Abounacer et al. 2014).
In the e-constraint method, a part or the entire efficient set for a M-O problem is produced. This method can provide a representative subset of the Pareto set, which in most cases is efficient. In this method, the decision maker chooses only one objective as its main one to be optimized. The other objective functions are added as constraints to be less than or equal to given epsilon values (Khalili-Damghani et al. 2013). Then, by variation in the epsilon values of the constrained objective functions, the efficient solutions of the problem are obtained. Moreover, to ensure that only non-dominated solutions will remain in the obtained Pareto front, the augmented version of the e-constraint method is implemented (Mavrotas 2009;Mavrotas and Florios 2013) which avoids the production of weakly Pareto optimal solutions. For a detail description of the econstraint method, one can refer to Mavrotas (2009) (Mavrotas 2009).
For the implementation of the e-constraint method, the first objective (1 in the model) is used as the main objective function, and two other objectives (2 and 3 in the model) are placed in the set of constraints. Since the main objective is maximized and two others are minimization functions, in our implementation these constraints are first set to a large value and then decreased step by step in each iteration. This method is implemented in CPLEX 12.6 to solve the created instances of the MO-MM-GTTDP in 20 iterations as the stopping condition. Computational examination demonstrates that the exact method sometimes requires large computing times. Therefore, a metaheuristic algorithm is proposed and implemented to solve real size instances of the MO-MM-GTTDP.

The multi-objective variable neighborhood search method
As mentioned in the literature section, the VNS has been successfully applied to many variants of the OP, such as TOP (Archetti et al. 2007), TOP with time windows (Labadie et al. 2012), and OP with Hotel Selection (Divsalar et al. 2013). Therefore, in this paper, the M-O version of the VNS algorithm (MOVNS) is developed to tackle the MO-MM-GTTDP, which is in fact modeled as a M-O OP.
Considering that the presented problem is M-O with three conflicting objectives; this is the first use of the MOVNS algorithm for solving an OP variant. In this algorithm, the optimizing concept is considering two main views: First, how to define the set of pre-established neighborhoods, and second, how to explore using these neighborhood structures (Duarte et al. 2015). A general overview of the applied MOVNS is presented in Fig. 1. According to Algorithm 1, after some preprocessing, at the first step of the algorithm, three initial solutions are generated. This initialization phase is explained in Section b. Then the improvement phase is started with the ''shaking'' phase followed by applying the local search part which is called the ''MOVND''. 2 Both shaking and the MOVND are presented in detail in Sections c and d, respectively. In the MOVND part of this algorithm, three main sets of neighborhoods are used, which correspond to three objective functions. In Section e, the ''update'' method as the third part of the improvement phase is described. In this method, the available Pareto solutions are updated so that finally dominated solutions are deleted, and the only effective (non-dominated) Pareto solutions are remained and moved to the next iteration. Obviously, these steps are repeated as long as the maximum iteration number is not exceeded.

Solution representation
To implement the proposed MOVNS, each solution of the MO-MM-GTTDP is represented as a two-dimensional array containing the list of selected nodes as well as the transportation mode that has been selected for moving between the selected nodes. The sequence of visiting nodes starts at the origin node (1) and ends at the destination node (n). An example of the solution representation of an instance of the problem with 20 nodes and 3 modes of transportation facilities is presented in Table 2.

Initialization
In the initialization part of the algorithm, due to the existence of three objective functions of the problem, three initial solutions are constructed; each of them corresponds to one objective function when creating a solution. In the preprocessing step, four matrices, including the distance between every pairs of nodes (matrix d ), the travel time (matrix t ) and the cost (matrix c ) of moving between every pairs of nodes with any transportation mode, as well as the amount of co 2 emission that is produced by moving between every pairs of nodes using each transportation mode (matrix p ), are generated. It should be mentioned that there are three types of constant values taken as the input and used in creating these matrices: constant-time, constant-cost, and constant-pollutant. Obviously, these constants depend on the transportation mode. Moreover, the Score and the visiting cost values for each node are used in the description of the algorithm as vectors of S(i) and VC(i).
All three initial solutions are made based on the famous nearest neighbor strategy. However, in each of them, a different selection measure is used. In the first initial solution, among the feasible non-included nodes, the selection priority for visiting node(i) is the maximum ratio of S(i) over the distance of the last visited node to node(i). After selection, node(i) is added to the end of the visited nodes list, and the process continues until no more node is feasible to visit. It should be noted that in this solution, it is assumed that the fastest transportation mode is selected for moving between POIs. In this case, the fastest facility is the one with the minimum constant-time. Considering the second objective of the model, for the second initial solution, the denominator of the mentioned ratio for node selection is cost (i). The cost of node(i) is equal to VC(i) plus the transportation cost of moving between the last visited node to node(i). For this solution, the cheapest transportation mode, the facility with the minimum constant-cost, is selected for moving between POIs. Creating the third initial solution is similar to the first except that the cleanest transportation mode is selected. An overview of the process of generating the first initial solution is presented in Fig. 2.

Shaking phase
According to the general structure of the algorithm (Algorithm 1), the Multi À ObjectiveShake is the first step in the improvement phase. The main goal of the shaking phase is to diversify the solution search space by skipping from local optimum. Considering the M-O nature of the problem, since both the sequence of visited nodes and transportation modes may affect the solution quality in terms of any of the objective functions, in the shaking phase of the proposed algorithm, two shaking moves are implemented: Node-Remove and Change-TrpMode. The first move is to remove some of the visited nodes randomly, and then to apply the local search to improve the solution. The number of nodes to remove for each solution is calculated by multiplying a constant parameter (a) by the number of visited nodes in the solution. At the same time when removing node (i) from the visited node list, the transportation mode for moving between node (i) to node (i?1) is removed, and the same transportation mode used to move between node ði À 1Þ and node ðiÞ is used to move between node ði À 1Þ and node ði þ 1Þ. The second move in the shaking phase, Change-TrpMode, randomly changes the transportation modes for some of the arcs between visited nodes while keeping the solution feasible. The number of arcs for which the transportation modes are changed is again obtained by multiplying the constant parameter (a) by the number of visited nodes in the solution. It should be mentioned that the first shake move is applied on each solution in the current Pareto list. However, the second move is only applied on half of the solutions which are selected randomly. Algorithm 3 (Fig. 3) represents the pseudocode of the shaking phase.

MOVND
The local search part in the improving phase of the presented MOVNS algorithm is called the M-O variable neighborhood decent (MOVND). Again, according to the M-O nature of the problem, the MOVND contains three sub-parts that improve each solution according to each objective function. In fact, in this algorithm, three ''VND'' 3 structures are used to improve the solution in terms of three objective functions, separately: VND-1, VND-2, and VND-3. As can be seen in Algorithm 4, in the overall structure of the MOVND, each VND includes local moves specifically designed to improve the solution in terms the corresponding objective function locally. More specifically, VND-1 (similarly VND-2 and VND-3) tries to improve the value of z 1 (symmetrically z 2 and z 3 ) and find the local optimum with respect to z 1 (uniformly z 2 and z 3 ) by using the implemented neighborhood structures. One main difference of the M-O over the single-objective optimization technique is that it is not enough only to check whether the final solution obtained by the VND is improved compared to the incumbent solution, but it needs to consider every explored solution during the search and give them a chance to be part of the Pareto front (Duarte et al. 2015).
As shown in Fig. 4, all different solutions explored during a local move search are stored in the list S. After each VND, an Update function is applied to compare the solutions in S with the ones that already exist in the Pareto list. This results in an updated Pareto list after each local move. The detail of the update process is presented in the next section. In the remainder of this section, each VND is explained separately.
As mentioned earlier, the VND-1 is designed to improve the first objective function of the model. So, concerning the time budget of the tour, two local search moves are used to raise the total score directly or to reduction the total travel time with a hope to be able to visit more POIs and increase the total score: insert-Node and Two-Opt. A brief description of each move is given as follows: Insert-Node: this move is increasing the total score of the tour by inserting nodes to the visited node list. Among non-visited nodes for which the insertion does not make the tour infeasible, the best insertion position of node (i) is found with the maximum ratio of the score over the length (i) where length (i) for the insertion of node (i) between node (j) and node (j þ 1) is efficiently calculated using length (i) ¼ t j;i þ t i;jþ1 À t j;jþ1 . It should be mention that the transportation mode for moving between node (j) and node (i) as well as between node (i) and node (j þ 1) are assumed to stay the same as the former mode which is used between node (j) and node (j þ 1) in the current solution.
Two-Opt: Due to its efficiency, this move is one of the most commonly used local neighborhood structures in routing problems. It reduces the travel time in the tour by replacing two existing edges with two new ones. In the implementation, for every possible pair of nodes, among the visited nodes it is checked whether the reversal in the order of nodes between them produces a travel time saving, and then the pair of nodes with the highest saving is chosen. After selection, the order of the nodes, as well as the transportation mode for moving between them, is reversed. To check out the feasibility of this move, the time matrix ðmatrix t Þ is being used.
VND-2 is designed with respect to the second objective. The only local move used in this is called the Change-Node.
Change-Node: It tries to minimize the total cost of the tour by replacing a visited node by a none-visited node. As the first condition, node (j) is replaced with node (i) if this replacement decreases the total cost of the tour. This cost includes the visit cost, and transportation cost. If this change is feasible, then the difference between the scores is checked out, and with the increase in the total score, node (j) is replaced with node (i). If such a node is not found (i.e., the score of node (i) is smaller than the score of node (j)), then the replacement happens for the node with the maximum of the diff C =diff S ratio. diff C represents the difference of node (j) and node (i) cost (including the visit and the transport), and diff S is equal to the difference of those nodes' score. Figure 5 represents the pseudocode of Change-Node.
The third part of the MOVND is the VND-3. This part concentrates on the third objective function and includes one local move which tries to decrease the produced co 2 emissions over all the tour, according to z 3 . This local move is called Change-TrpMode. It is designed to change the current transportation mode between POIs with a less pollutant transportation mode while respecting the time budget constraint. It starts from the first arc of the tour and the transportation mode used for moving between the first and second nodes of the visited node list. It greedily changes the transportation mode to reduce co 2 emissions. To be more specific, in each application of this local move, the transportation mode of an arc with the maximum ''diff p '' value is changed. diff p represents the difference between produced CO 2 emissions by using the new transportation mode instead of the current one for an arc. The pseudocode of the Change-TrpMode is presented in Algorithm 6 (Fig. 6).

Update Pareto list
Two decision sets will be non-dominated or indifferent if under the existing state of knowledge no logical basis can be found for ranking one as 'better' than another when all commensurable and non-commensurable objectives are being considered (Haimes and Hall 1974). The task of keeping track of such solutions in the presented method is done with the Update Pareto List module. According to the presented M-O model, this module has an important role in our MOVNS algorithm. As shown in Fig. 1, in an iteration of the algorithm, the Update function is called in four places of the MOVNS. This function is designed to compare all the obtained solutions after application of each VND with the existing Pareto solutions. Thereby, every time new found solutions are added to the current list of solutions we ensure that weak optimal solutions are not chosen, and non-dominated solutions are not left out. More specifically, these comparisons are performed using the corresponding objective values of each solution. As a result of this phase, all none-dominated solutions are remained in the Pareto list and the duplicated solutions as well as the dominated ones are removed from the list.

Pareto solution selection
As explained, two types of tourists are considered regarding the priorities on the objective functions. For the first type, who has fixed and predetermined priorities on the objective functions, and is not willing to sacrifice the value of an objective with a higher priority to gain any value of the objective with lower priority, a lexicographical method is used. To provide such a solution to the tourist using our MOVNS approach, an optional module is implemented as a post-processing step. In this step, the final Pareto list obtained by the MOVNS is lexicographically sorted first based on their score objective value, then the total cost value, and last their total emission value. Thereby, the tourist is provided by one single solution for which the lexicographical order of objective functions is satisfied.
The second type of tourists who might give different priorities to the objective functions would like to benefit from having a list of Pareto (non-dominated) solutions and selecting a solution from this list. This type of tourists may use an MCDM approach helping them in selecting a route among non-dominated solutions. In our case study, the TOPSIS method is being used which gives the tourist the opportunity to have paired comparisons on objective functions.

Experiments and results
The proposed MOVNS algorithm is coded in Visual C??. The mathematical model and the lexicographic method as well as the e À constraint method are implemented in CPLEX 12.6 to evaluate the performance of the presented metaheuristic algorithm. First, new benchmark instances are generated to evaluate the performance of solution methods. These instances are described in Sect. 5.1. Generally, the performance of solution methods in solving the M-O models is evaluated by a set of standard criteria. The selected evaluating criteria in this research are described in Sect. 5.2. Before applying the MOVNS method on the benchmark instances, parameters of this algorithm should be tuned. In Sect. 5.3, it is described that how the Taguchi experimental design (Taguchi et al. 2005) is used for parameter tuning in this research. Section 5.4 is divided in two parts. First, the results of the proposed MOVNS are  compared with the results of the epsilon constraint and the lexicographic method over all generated benchmark instances and analyzed based on the evaluation criteria. Then, the MOVNS is compared against two other metaheuristics including a basic VNS as well as an iterated local search.

Test instances
Because there exist no test instances for the MO-MM-GTTDP in the literature, a set of standard instances are generated to evaluate the work of the proposed algorithm. These instances are generated using the existing instances for the OP (Tsiligirides 1984) by changing the number of nodes and available time for the tour. In addition, some new data are added to these instances, including the visit cost of POIs, travel time and cost for moving between POIs by any of available transportation modes, and also the amount of CO 2 emissions produced by any of available transportation modes. The visit cost of each POI is generated randomly with a uniform function in the interval of [10,150]. In all generated instances, three transportation modes are considered (e.g., Bicycle, Bus, and Taxi). The average travel time per distance unite for moving between POIs by using each transportation mode is assumed to 6, 4, and 2 units, respectively. Moreover, the cost of moving between POIs using each transportation mode is set as 0.3, 1, and 5 units per distance unite, respectively. The pollution constant per distance unit using each transportation mode is also assumed to be 0, 0.069, and 0.17, respectively.
The set of benchmark instances in this research is generated based on set1 of OP instances of Chao, et al.,  which is available on [www.mech.kuleuven.be/ cib/op]. In total, 30 instances are generated. Table 3 demonstrates the main characteristics of these instances. In general, Instances with 5 to 32 available POIs with a time budget of 10 to 80 are generated and used for evaluation of the algorithm. It should be mentioned that larger instances are unlikely to be solved using the implemented epsilonconstraint method in CPLEX.

Evaluation criteria
Considering that the outcome of the M-O algorithms is usually an approximation of the Pareto optimal set, an important issue is how to compare the performance of these algorithms with each other (Zitzler and Thiele 1999). Traditionally, there are two goals for a M-O procedure: 1) a good convergence to the Pareto optimal front, and 2) a good diversity in obtained solutions. Three types of metrics are usually used for performance measurement, which evaluates 1) spread of solutions on the known Pareto optimal front, 2) convergence to the known Pareto optimal front, and 3) certain combinations of convergence and spread of solutions (Wang et al. 2011). In this section, next to the Cpu-time, six other criteria are selected for the performance evaluation of the proposed MOVNS and briefly described. SID, SI, and NPS belong to the first set of above explained metrics. MID, and IGD are the metrics from the second set which find the convergence quality, and the Hyper Volume (HPV) is computed to measure the combination of convergence and diversity of the results obtained from both MOVNS, and the e À constraint method. Since the MO-MM-GTTDP has three objective functions, the set of effective Pareto solutions are assumed to be as: f 1 1 ; f 1 2 ; f 1 3 À Á ; . . .; f n 1 ; f n 2 ; f n 3 À Á .

Cpu-time
The first performance evaluation criterion is the computational time of the algorithm. For NP-hard problems, it is an essential criterion for evaluating the solution methodology, specifically for the desired application of the MO-MM-GTTDP in tourist trip planning.

Maximum spread index (MSI)
Zitzler and Thiele (Zitzler and Thiele 1999) present the MSI for the first time in 1999. It shows the diversity of effective Pareto solutions. Obviously, the larger value of MSI indicates that solutions are scattered in a wider space, which shows the higher performance of the algorithm. This criterion is calculated using Eq. (14).
In this equation, K is the number of objective functions and f max j and f min j are the highest and lowest values obtained by the algorithm for the j th objective function.

Spacing index (SI)
The second evaluation criterion measures the uniformity of the non-dominated solutions in the solution space. To calculate this measure, first, the minimum value of the Euclidean distance between every pair adjoining Pareto solutions in the solution space is calculated by Eq. (15).
In Eq. (15), n is the number of Pareto solutions, K is the number of objective functions, and f i k is the amount of k th objective function in the i th solution. Then, the SI value is obtained by Eq. (16) as follows: where d is the average of d i s. The lower this criterion, the better the uniform distribution of the solutions.

Number of Pareto solution (NPS)
The third criterion shows the number of approximate nondominated Pareto solutions that evaluate the performance of the algorithm. Obviously, having a higher value for the NPS means the higher quality of the M-O algorithm.

Mean ideal distance (MID)
The MID metric finds the mean deviation of Pareto solutions from the ideal solution calculated by Eq. (17) (Maghsoudlou et al. 2016). In this equation, n shows the number of Pareto solutions, j is the index of the objective function, and f j Max and f j Min represent the highest and the lowest values of j th objective, respectively.f j Best shows the ideal value of the corresponding objective. We consider the best value of each objective function among the solutions in the Pareto list as the corresponding ideal value. The lower the value of MID, the higher the quality of the Pareto list.

Inverted generational distance (IGD)
This performance indicator evaluates the quality of an obtained solution set based on the distance between the solution and a reference set (Ishibuchi H, Masuda H, Tanigaki Y, Nojima Y. Modified distance calculation in generational distance and inverted generational distance. InInternational conference on evolutionary multi-criterion optimization 2015). Equation (18) formulates the IGD indicator, where p* demonstrates the set of reference point, PF denotes the non-dominated solutions generated by the algorithm, dist(p, PF) shows the nearest distance from p to solutions in PF, and the distance from p to the solution y in PF is calculated by dðp; yÞ ¼ P m j¼1 ðp j À y j Þ 2 (Sun et al. 2018). Here, m shows the dimension of each Pareto solution. In our case, PF is the list of non-dominated solutions obtained by the MOVNS, and p* is the Pareto front obtained by epsilon-constraint. Therefore, the IGD is only calculated for the results obtained by our MOVNS.

Hyper volume (HPV)
The hyper volume indicator is a measure of the quality of a set PF = {p 1 ,p 2 ,...,p n } of n non-dominated objective vectors produced in a run of a M-O optimizer. Assuming a minimization problem involving m objectives, this indicator consists of the measure of the region which is simultaneously dominated by PF and bounded above by a reference point r 2 R m such that r ! (max p p 1 ,...; max p p m ), where p = (p 1 ,...; p m ) 2 P & R m , and the relation C applies component-wise. As illustrated in Fig. 7, this region consists of an orthogonal polytope, and may be seen as the union of n axis-aligned hyper-rectangles with one common vertex (the reference point,r) (Fonseca et al. 2006). In this research, the approximation method introduced by Deng, and Zhang (Deng and Zhang 2020) is used to calculate the HPV for the results obtained by both the MOVNS, and the epsilon method. The reference point is the worst objective values obtained by any of these methods. In our calculations, all values are normalized in [0,1] and (1.1, 1.1, 1.1) is used as the reference point. Fig. 7 Region show the HPV of a given Pareto front, and a reference point (r) (Fonseca et al. 2006)

Taguchi-based Parameters tuning
One of the important factors that may have an impact on the efficiency of the metaheuristic algorithm is the predetermined value of parameters of the algorithm. Taguchi experimental design is a well-known method for achieving a good value for parameters of metaheuristics without performing the full factorial possible combinations of experiments. This method is on the basis of two crucial factors: the orthogonal array (OA) and the signal to noise (S/N) ratio. The OA represents a matrix that is carrying the experimental scheme depend on different factors' levels (Azadeh et al. 2017). In the S/N ratio, S and N represent the term ''signal'' as a useful value and ''noise'' as an unpleasant value, respectively. The Taguchi method tends to minimize the variation (S/N ratio) by minimizing the effect of noise and find the best level of the influential controllable parameter based on the robustness (Tsai et al. 2007). The first step in using this method is to determine the controllable parameters of the proposed algorithm. A positive side of the proposed MOVNS is that there are only two parameters to be set before running the algorithm: the iteration number and the alpha coefficient. As described in Sect. 3.d, the Alpha coefficient is a constant value that has a direct impact on the number of removable nodes as well as changeable transportation modes during the shaking phase of the algorithm. Five levels are considered for each of these two parameters, which are presented in Table 4. These values are selected based on some preliminary experiments.
To tune the parameters of the proposed method, the above-described evaluating criteria (MID, MSI, SI, IGD, HPV, NPS and CPU) are used as responses of each experiment design to determine the most appropriate level for each parameter. To do these analyses, the Minitab software is used. It should be noted that all these experiments are implemented on three randomly selected instances with different problem sizes. Before this selection, all generated instances are classified in to three categories small, medium, and large size according to the number of POIs and the available time budget. Table 5 shows the classified instances and the ones selected for the parameter tuning. Results obtained from selected instances are normalized and averaged before applying the Taguchi method. Figures 8, 9 and 10 illustrate the results of the parameter tuning obtained from the Minitab. According to the S/N ratio plots (Figs. 8,9 and 10), it is inferred that the best value of the iteration number and the alpha coefficient for small instances are 50 and 0.5, for medium instances are 50 and 0.8, and for large instances are 100 and 0.4 respectively.

Results and discussion
In this section, the detailed results of the proposed MOVNS are presented and analyzed in two main subsections. First, all instances are solved by the MOVNS and its performance is compared with both thee À constraint, and the lexicographic methods. Second, two other versions of the

Comparing with epsilon and lexicographic methods
After parameter tuning, all mentioned instances are solved using both implemented methods. All tests were conducted on a laptop with Intel core i7, 3.1 GHz CPU and 8.00 GB RAM. The results of solving instances by the presented MOVNS algorithm as well as the e-constraint method are summarized in Table 6. The first column in this table represents the instance's number. The second column gives the instance's name, which includes the number of POIs and the available time budget. Then the value of each evaluation metric is presented for both MOVNS and econstraint methods, respectively. As explained in Sect. 5-2.6 the IGD is only calculated for the MOVNS. Because of the stochasticity of the MOVNS, each instance is solved for 30 times and the best and the average results over these 30 runs are reported.
Comparisons of the results of both algorithms based on each evaluation metric are illustrated in Figs. 11,12,13,14,15 and 16. From Fig. 11, it can be seen that in terms the MID metric, the best results of the MOVNS are always better than or equal to the results obtained by the epsilon method, and the average results are very much close, and competing, and in many cases better than the epsilon method.
The MSI criterion is presented in Fig. 12. According to this figure, the value of the MSI in the results of the MOVNS is approximately the same as the epsilon-constraint method for most of the instances. More precisely, as expected, in few instances with less time budget, the econstraint method provides Pareto solutions with more extensive solution space, but in instances with 25 POIs or higher, the presented MOVNS algorithm represents a higher value of MSI. Generally, according to the obtained MSI criterion over the results of all instances, it can be concluded that our MOVNS algorithm has acceptable performance to provide Pareto solutions with wide solution space. Figure 13 illustrates the uniformity of the set of Pareto solutions obtained by each of the two methods. As can be seen, in most cases, the MOVNS is superior in terms of spacing index criterion. Figure 14 shows the number of Pareto solutions obtained by each of the implemented algorithms. It is worth mentioning that this measure is relatively unfair because we fixed the maximum number of Pareto solutions that can be obtained by the e À constraint method on 20. However, the comparing procedure in some instances that have less than 20 Pareto solutions in the e À constraint method is considerable.      HPV is one of the most trustworthy metrics in comparing M-O methods which consider both convergence and diversity at the same time. From Fig. 15, it is evident that the MOVNS performs well from the HPV metric point of view. In all instances, the HPV obtained from the MOVNS results is either the same or higher than the results of the econstraint method. Moreover, it is visible that in all cases the average and the best results of the MOVNS are very close to each other, which shows the robustness of the MOVNS in providing high-quality solutions.
The computational time of both algorithms is compared over the generated instances in Fig. 16. As the problem size increases, the Cpu-time of the e À constraint method becomes more than 4 h. Consequently, the exact algorithm is not applicable to real-world problems.
To have a better view on the effectiveness of the presented MOVNS in practice, an example of the obtained Pareto solutions for the largest generated instance with 32 POIs and the time budget of 80 using both MOVNS and e À constraint method is plotted in Fig. 17. Although, it is a 3-dimensional plot, it can be seen that the MOVNS is able to find higher number of non-dominated solutions which are almost uniformly distributed along the Pareto front. Moreover, when we compared both Pareto fronts in terms of non-dominancy, only 1 solution from the MOVNS results was removed (movns-rem), while 6 solutions from the epsilon results were dominated by the MOVNS (epsrem). Therefore, it can be concluded that the MOVNS performs effectively for this instance.
So far it is illustrated that the MOVNS is effective in finding Pareto solutions. However, as explained in Sect. 4.3.6, for the tourists of type 1, an optional module is used to post-process the found Pareto list and provide one best solution for the tourist. To evaluate the efficiency of the MOVNS plus this post-processing in providing lexicographic best solution, its results are compared with the exact lexicographic method which is explained in Sect. 4.3 and implemented in CPLEX. Each instance was solved 30 times by the MOVNS and once by the epsilon method. The best lexicographic optimization solution obtained by the first run is selected as the MOVNS result. As a result, for all instances, the MOVNS was able to find exactly the same solution as achieved by the exact method, which shows that the presented method is indeed very effective in finding the lexicographic best solution. These results are presented in Table 7.

Comparing with other metaheuristics
To evaluate the proposed MOVNS and show its performance over other existing methods in the literature, two other versions of the solution method are compared with the proposed MOVNS.
The first metaheuristic to compare with our method is a basic variable neighborhood search, called the ''BASIC_VNS,'' in which the main properties of the proposed VNS become inactive. Namely, the following changes are applied: 1) Only the first initial solution method is used and the other two initial solution strategies are ignored (See Sect. 4.3.2; 2) in the shaking phase, the Change-TrpMode is ignored and only the Node-Remove shake is applied (See Sect. 4.3.3; 3), in the MOVND part of the algorithm, only the first VND is kept active (including the first two local moves), and VND-2 and VND-3 are skipped (See Sect. 4.3.4). Fig. 17 A 3D plot of the Pareto fronts obtained by MOVNS and e-constraint for instance (32,80) As described in the literature section (Table 1), there are other multi-objective variants of the OP in the literature with various objective functions. One of the more recent solution methods is the Iterated Local Search method of Hapsari et al. (2019). The ILS method is one of the few methods from the literature (according to Table 1) in which arc-based objective functions are considered (similar to the current study). Moreover, it is a local move-based method and, in this sense, is closer to the proposed MOVNS. Therefore, as the second metaheuristic to compare with the proposed MOVNS, we redesigned our method using the ILS structure. More specifically, the ''Initialization'', and ''Shaking'' phases are kept the same, and the same local moves as presented in Sect. 4.3.4, are placed in an ILS structure, where in one iteration, all moves are applied consecutively and not in a VND structure. Then, the ''Update'' procedure is only applied once at the end of the iteration. This adapted algorithm is called the ''ILS.'' All instances are solved 30 times by each of the three algorithms, the ''BASIC-VNS'', the ''ILS'', and the ''MOVNS''. Then, for each instance and each performance metric, the Wilcoxon's rank sum test (Wilcoxon 1945) is conducted two times: first, between the results of the MOVNS and the BASIC-VNS, and second, between the results of the MOVNS, and the ILS under the significance level of 5%. Table 8 shows the number of instances that for each performance metric the MOVNS outperforms each of the other two algorithms.
The best value of each performance metric for each instance and each method is presented in Table 9. In this table, for each instance the best value of the three algorithms is bolded out.
From these two tables, it can be observed that, except for the MID and the CPU, for the other metrics the MOVNS outperforms the ILS. If we look at the definition of the MID, as the deviation from the ideal solution which is in our case the best solution obtained by each method, it can be concluded that although on average the MID values of the ILS results are lower, it might be because of the lower quality of the considered ideal solution rather than the higher quality of the obtained Pareto set.
The comparison of the MOVNS and the BASIC-VNS illustrates the effectiveness of the methodological contributions of the proposed solution method. Moreover, each of the main parts of the MOVNS method is analyzed further in Sect. 5.5.2.

Sensitivity analyses
In this section, two types of sensitivity analyses are performed. The first type is on the problem definition by using the mathematical model, where we investigate the consideration of the economic budget as a hard constraint rather than considering the total cost as an objective to be minimized. Second type is on the solution algorithm, where  An optimization approach for green tourist trip design 4323 the importance of some of the main modules of the proposed MOVNS is checked and discussed.

Analyses on the problem definition
We use an example to show that the user would benefit from considering the total cost as an objective function, thereby having a list of Pareto solutions rather than imposing the budget as a hard constraint to the model. We make our reasoning based on the medium size instance with 20 POIs and 40 as its time budget. The Pareto list obtained by applying the e-constraint is presented in Table 10. First, we assume that the user has a hard budget limit of 300 for his/her trip. According to the solutions, the best available solution would be number 4, where the total score of 15 is gained by spending about 259 out of the 300 budget. However, if the Pareto solutions were provided with the user, he/she might select the 5th solution, where spending almost 20 more than the predetermined budget would lead to a higher score. As another example, if the available budget of the user would be 400, again the best obtained solution of the model becomes number 5 with the total score of 20. However, one might realize that spending about less than 5 cost unit more than the available budget would provide a much higher score by selecting solution number 8. With these analyses, it is evident that considering the total cost as the objective function of the model may have benefits for the user compared to imposing the budget as a hard constraint to the model.

. Analyses on the solution method
To show the effectiveness of the main modules in the MOVNS, especially on the importance of considering cost, and the CO 2 emission in the decision-making process, two sensitivity analyses are conducted. First analyses are performed on the initialization phase, where we analyze the sensitivity of the solution quality to ignoring the second, and the third initial solutions where the cheapest and the cleanest transportation modes are selected, respectively (See Sect. 4.3.2). Then, we perform a similar sensitivity analysis on the effectiveness of the VND-2 and VND-3 (explained in Sect. 4.3.4) on the solution quality. For this purpose, a medium size instance is randomly selected with 25 POIs and 50 as its time budget. Each time, one of the mentioned algorithm parts is skipped in the solution method, and the instance is solved for 30 times.
The best values of each of the 3 performance metrics (MID, MSI, HPV) are calculated over these runs. The obtained values are then compared with the corresponding ones achieved using the full version of the algorithm. Table 11 compares these values. As can be seen, for all metrics best values are obtained when the full algorithm is used. This is specifically recognizable for the HPV in which both convergence and diversity of the solutions are taken into account.

Case study
In this section, the proposed solution approach has been tested over a real dataset obtained from the city of Tehran. Tehran is the capital of Iran and one of the main important touristic places in the country. It has numerous famous touristic attractions. According to available data from 2018, annually, 1,200,000 foreign tourists visit Tehran. As a small test set, 10 POIs are selected from all the existing POIs in Tehran. These POIs are chosen from four different categories: park and green spaces, handicraft centers, museums, and shopping centers. Table 12 presents the name of each of these POIs as well as their visiting cost.  The corresponding score for each POI and are generated randomly in the range of [0, 10]. The available time budget for the tour is considered 4 h. Due to the available transportation facilities in Tehran, we proposed three transportation modes to the tourist for moving between each pair of POIs: the Subway, BRT, 4 and Taxi. All the data related to the distance and available transportation modes between all pairs of nodes, the cost of traveling as well as the amount of CO 2 emission which is produced by each mode are shown in Table 13. Time values are the average value over different times of the day that obtained using Google Map to relatively considering the traffic. Moreover, the amount of produced CO 2 emission is based on the available information on the [www.co2nnect.org] that is 0.069 kg/km for BRT, 0.17 kg/km for Taxi, and 0 kg/km for Subway per person. The presented MOVNS approach is applied to solve the real data instance. The results contain 46 non-dominated solutions. All of these non-dominated solutions are illustrated in Fig. 18. For some of these non-dominated solutions, more details, including the proposed route, are presented in Table 14.
As can be seen in Fig. 18, all Pareto solutions are clustered in three groups displaying the impact of each objective function. This means that one group members have higher scores than the other groups. The second and third groups have lower travel costs as well as less pollution produced during the tour. Therefore, it can be concluded the proposed solution method (MOVND) has acceptable performance to provide Pareto solutions with considering all the aspects of the problem.

Selecting routes among the non-dominated solution
For the type 1 tourist who would not like to reduce his/her happiness (score) to reduce the cost or to reduce the CO 2 emission, the Lexicographic method is used. This gives one best solution out of the 46 non-dominated ones with objective values of: Z 1 ¼ 52,Z 2 ¼ 56500, andZ 3 ¼ 4:20.
However, for the second type of tourists, considering the high number of obtained non-dominated solutions, choosing one best tour for the tourist may still be a difficult task. Selecting a solution among the non-dominated ones is itself a challenging problem in the area of M-O optimization (MOP), which is not in the scope of the current research. However, to make the procedure of tour recommendation more practical, in the rest of this section, a simple but practical multi-attribute decision-making (MADM) technique is used to indicate which solutions might be more interesting for the tourist among the 46 obtained solutions. As a result, at the end of this section, three routes are proposed to the tourist with considering his/her preferences.
In the first step, the existing criteria are weighed with the paired comparison method. This paired comparison method is introduced by Saaty (Saaty 2008) as one of the high performance multi-criteria decision-making techniques in 2008. Table 15 shows the paired comparison matrix between every two criteria, which is filled with the 9 scale standard range of the pairwise comparison. It should be noted that in this case study, the paired comparison matrix is filled artificially. In the second step, we rank all the obtained solutions by using the TOPSIS method. TOPSIS is one of the MADM methods developed originally by Ching-Lai Hwang and Yoon (Hwang and Yoon 1981) in 1981. This method is one of the most common multi-criteria decision-making techniques, which compares a set of alternatives by identifying weights for each criterion, normalizing scores for each criterion, and calculating the geometric distance between each alternative and the ideal alternative. The ideal alternative in this method is the one with the highest score in every criterion (Rezki and Aghezzaf 2017). As can be imagined, this ideal alternative does not necessarily exist among the list of alternatives. For more details on the TOPSIS method, one may refer to the research by Triantaphyllou et al. (Triantaphyllou 2000).
In the current case study, all 46 non-dominated solutions are considered as alternatives in the TOPSIS method, numbered from A 1 to A 46 . Each alternative A i has three criteria value corresponding to each objective value of Z 1 , Z 2 , and Z 3 . The final rank of each non-dominated solution depends on the paired comparison matrix. In practice, this matrix can be generated according to the tourist preferences over objective functions.
In the presented matrix, the priority is given to the tourist desirability by maximizing the total score. Then, the environmental aspect of the tour is considered with a lower weight. The minimum weight is given to the travel cost. Indeed, these decisions may have an impact on the selection of transportation modes. Applying the TOPSIS method, the normalized vector of weights for criteria is found as 0:222518; 0:126834; 0:650648 ½ . The weighted normalized decision matrix (V) is then calculated by multiplying the vector W by the decision matrix. Equations 19 and 20 show the positive and the negative ideal alternatives obtained from the V matrix, respectively.   An optimization approach for green tourist trip design 4327 solutions, respectively. Equation 21 calculates the relative closeness from the negative ideal for each alternative.
Then, all 46 alternatives are ranked based on the corresponding CL i value. Table 16 represents the three best non-dominated solutions after this ranking. These three solutions are graphically presented in Figs. 19,20 and 21. According to the TOPSIS method, Fig. 19 shows the route with the first rank. As we see in Table 16, this solution has the lowest value of traveling cost. Considering the Paired Comparison matrix of three objective functions, the Taxi is not proposed for moving between POIs in this route. Figures 20 and 21 represent the second and third routes, respectively. The second route has more score by suggesting a POI with more visiting cost as well as score. It should be mentioned that the total produced pollutant is the same in these two routes, Because of the facility mode, which produced more pollutants has a higher moving cost per unit of distance. In the third route, the value of total score (Z 1 ) is very different from two other routes. Obviously, as the value of the first objective function increases, the total cost and the produced pollutant are also increased. We can conclude from the obtained solutions that all the three aspects of the problem have been considered. Then, based on the tourist's priorities, the three best routes are presented.
Tourists as users of this system are considered as two types. For the first type, a lexicographic optimization method provides one best solution with the first priority to total score, second priority to total cost, and third priority to total emission of the proposed tour. For the second type, a Pareto list of solutions is provided. We propose a MOVNS algorithm to cope with this problem. On the basis of the trade-off between improving three objective functions, to generate combinations of the selected POIs and the transportation facility modes, high-quality non-dominated solutions are obtained in short computation times. We solve the problem by implementing it in CPLEX using both the lexicographic and the e-constraint method as exact solution approaches. We use our 30 generated instances based on benchmark OP instances and compare two methods by using four M-O evaluating criteria. The results represent the effectiveness of our proposed metaheuristic algorithm in finding Pareto solutions as well as lexicographic best solution.   Furthermore, we implemented the proposed MOVNS algorithm on a small real case in Tehran to prove the applicability of the method. In addition, we used the TOPSIS method as an effective multi-criteria decisionmaking approach to summarize and suggest the three best tours to the tourist.
This model can be more realistic by considering time windows or the time dependent travel times using various transportation facilities in future works. A considerable challenge would be dealing with problems with multiple tours. Moreover, presenting a comprehensive algorithm which combines MOVNS approach with a multi-criteria decision-making method to simplify the selection for tourist can be considered in a future research.