STOCHASTIC MEDICAL TOURISM PROBLEM WITH VARIABLE RESIDENCE TIME CONSIDERING GRAVITY FUNCTION

. Medical tourism is a recent term in healthcare logistics referring to travel of patients to receive health services and spending leisure time in a destination country. This transferring of patients leads to access high-quality health services which are cheaper than the original country of patients. During this travel, passengers who are the patients from another country, have this opportunity for complimentary entertainment packages ( e.g. , pleasure tours) in the aftercare period. As far as we know, the term of medical tourism is rarely studied in healthcare logistics and such services are highly important for developing countries. Such facts motivate us to develop a practical optimization model for the Medical Tour Centers (MTCs) for allocation of patients to hospitals in proper time and creation of memorable aftercare time for them. In this regard, the main aim of the proposed model is to maximize the total profit of MTCs through optimal allocation of patients to hospitals while considering an aftercare tour for the passengers. To make the proposed model more realistic, the optimal residence time in attractive places is simulated by a time-dependent gravity function. To address the uncertainty of medical tourism problem, a scenario-based two-stage stochastic optimization approach is extended to encounter different sources of uncertainty existing in surgical success, medical time, restoration restrictions, and the attraction of tourist places. Another novelty of this work is to propose an innovative hybrid meta-heuristic for large-scale instances, which is a combination of Progressive Hedging Algorithm (PHA) and Genetic Algorithm (GA). The model is analyzed by different test problems for small, medium, and large-scale instances where the hybrid meta-heuristic algorithm could solve them with an average gap of 3.4% in comparison with the commercial solver. The results revealed the importance of tourist opinion and public preferences in medical and pleasure tours, respectively, to improve the economic growth in this sector in developing countries.


Introduction
The development of the tourism industry has a prominent role in promoting extensive economic dimensions of countries through earning revenue, creating job opportunities, reducing poverty, and extending social justice and welfare in society [75].The tourism industry counts as the third significant industry after the petroleum and automobile industries due to the significant effect on foreign exchange earnings [18].Over time, the tourism industry has progressed and been divided into different fields, such as cultural tourism, sports tourism, adventure tourism, and medical tourism [42].The global value of its market was approximately $47 billion in 2016.The term "medical tourism" expresses the condition when a person wants to travel abroad to cure diseases and spend leisure time in a host country [13].The people tend to go abroad to receive services due to their affordability with an equal or better quality of care than what is provided in their own country at the same time [51].In fact, the integration of treatment and tourism creates a subset of healthcare and tourism called medical tourism.Hence, tour centers try to facilitate tourism processes to attract more tourists, increase their profits, and improve the economies of their countries by providing low-cost and high-quality services.
Medical tourism is one of the main reasons for globalization and international logistics [17].As defined by the World Trade Organization (WTO), medical tourism is associated with international logistics in health services [8].The elements of the medical tourism supply chain include medical care providers and traveling agencies, accommodating and providing transportation services for allocation of persons to hospitals and other tourism sites [44].To improve customers' satisfaction and increase their profits in a competitive market, the Medical Tour Centers (MTCs) try to provide high-quality services, meeting their needs regarding the therapy of diseases, and creating an opportunity for them to utilize the attractiveness of the host or destination country in a healthcare journey to have a good leisure time [38].A pleasure tour is among the complementary benefits in medical tourism that can lead to competitive advantages, marketing, and profitability [16].Therefore, in addition to providing medical services, tour centers are responsible for transportation, accommodation, and recreation of tourist areas.Various suppliers, each with its specific limitations, provide these services for customers.Therefore, developing a comprehensive model to integrally address the different aspects of medical tourism issues is quite necessary.
COVID-19 pandemic has been huge effect on the tourism industry from the economic and social aspects for healthcare logistics and supply chains [55].One of the important and effective parameters in the revival of the tourism industry in the post Corona era and the reconstruction of the damaged industries is the sense of belonging and dependence on the place.The higher the sense of belonging of to a place, the more it can change decisions about returning to tourist destinations.A sense of place dependence can have a direct impact on tourist loyalty to tourism destinations.Solutions, such as trying to better understand the behavior of tourist, trying to lower costs and increase satisfaction, provide opportunities to attract customers by providing quality amenities and accommodation and reduce incidental costs.Recognizing the conscious and subconscious goals of tourists will improve customers' satisfaction [88].In medical tourism, costs, treatment quality, specialists, level of welfare, and tourism attractions of the host country are essential factors in medical tourism services [72].However, the information gathered from customers' experiences is another essential source in travelers' trip decision-making [86].People tend to select a hospital with high-quality treatment for curing their diseases; so, the customers' feedback creates a competitive environment in the healthcare system, promoting it by attracting more people through high-quality treatments [71].Therefore, as the framework proposed in this study, evaluating patients' interests in hospitals and comparing the feedback of previous clients, MTCs avoid making trips with dissatisfaction.They recommend effective treatment tours to customers at the right time, in the right place, and at the right price.
In an entertainment tour, determining the residence time of the tourists in tourist sites is another significant variable in travel planning, which is often neglected [85].If the residence time in a tourist destination increases, the desirability of visit increases rapidly at first, but it slowly reaches to the peak, while the costs of visit remarkably increase in the visiting period.Therefore, determining the optimal residence time of customers in each touristic place is a crucial measure for MTCs.It is even more critical when they organize trips for patients in a recovery period.It is also important that the treatment time of a patient is not definite and depends on both general conditions of the individual and the quality of treatment, which consequently influences the tour planning complexity [26].Hence, to take the risk of uncertainty into account, different sources like surgical success, overall conditions of the patient, hospital quality, and surgery team qualifications are all considered in the proposed model.This paper presents a novel mathematical model for MTCs to provide low-cost, high-quality medical tours for allocating passengers to hospitals and design aftercare pleasure tours.A two-stage stochastic model is developed to coordinate the treatment and tourism departments and cope with the inevitable uncertainty.The objective function of the proposed model is to maximize the total profit of an MTC to design a medical entertainment tour through distributing patients among various competitive hospitals and touristy cities in a scheduled form.We highlight the main contributions of this paper as follows: -This study presents a novel integrated medical-entertainment mathematical model.-A two-stage stochastic model is investigated in uncertain medical tourism tour plan.
-Public interests and passenger's preference are applied simultaneously.
-The optimal staying time in attractive places is determined using gravity function.
-Hybrid algorithm is proposed to efficiently optimize stochastic mathematical model.
The remainder of the paper is organized as follows.Section 2 reviews the related literature concerning the healthcare logistics considering medical tourism problem to identify the research gap and highlight the importance of this problem in terms of optimization model and solution algorithms.Section 3 describes and formulates the proposed problem.Section 4 explains our solutions for the proposed model.Section 5 provides numerical examples for the suggested model, followed by the computational results of the model to highlight the performance of proposed solution algorithm.Moreover, sensitivity analyses are provided to show the efficiency of proposed model in different cases.Finally, Section 6 presents the conclusions of this work, findings, limitations and recommendations for future studies.

Literature review
Many studies have been conducted on the community's attitudes toward the impacts of tourism on the development of countries from an economic perspective.Medical tourism is also a growing phenomenon in the tourism industry, and its emergence has been taken into consideration both academically and professionally [16].The research on medical tourism has started in the mid-2000s [40].Most of these introductory topics have addressed medical tourism concepts and the factors influencing the development of medical tourism [31,58].We divide the medical tourism literature into three categories: qualitative, statistical, and mathematical models.The medical tourism industry in qualitative and statistical studies has gained considerable attention over recent years.
In qualitative studies, some scholars have investigated criteria, which influence the development of medical tourism.For example, Yu and Ko [83] sought to identify infrastructural factors from the economic perspective to reduce trip expenses and waiting time by integrating organizations that are responsible for the provision of services.Buzinde and Yarnal [11] discussed the privileges of integrating hospitals and airlines to reduce expenses and attract more overseas patients.Han and Hyun [37] investigated the influence of quality, trust, customer satisfaction, and pricing of host countries on the attraction of customers.Heung et al. [39] stated that the obstacles of medical tourism were associated with costs, infrastructures, governmental policies, and governmental support.They reported that the barriers could be eliminated by promoting activities, increasing capacities, and improving service levels.Chen et al. [15] mentioned factors such as insufficient information about medical tourism, inaccurate schedules, and concerns over travel dissatisfaction as the barriers to medical tourists' traveling.Koeshendro et al. [43] identified a few research gaps in medical tourism literature concerning leisure travels in medical tourism, considering the personal opinions of travelers instead of using a common plan and paying attention to medical risks in the travel plan.Momeni et al. [54] introduced factors having significant effects on developing medical tourism, including marketing to attract international patients, global interaction, ethics and social culture, language, state-of-the-art technology and hospitals in a host country, transfer mode of foreign patients, commission, and international cooperation in global markets.Lovelock and Lovelock [47] pointed to the factors involved in medical tourism activities, including medical procedures, the extent of travel for post-treatment recovery, and the complexities of the itinerary.Table 1.Some criteria discussed in qualitative and statistical studies of medical tourism.

Criteria References
The integration of supply chain elements Yu and Ko [83]; Buzinde and Yarnal [11]; Some statistical studies have been conducted on tours and medical tourism considering various factors, including the integration of supply chain elements in terms of economic and social factors [75,87], the effect of capacity and price [67], the quality of treatment [48], behavioral theory, and customer satisfaction [60].Table 1 shows a few criteria mentioned in qualitative and statistical studies, some of which are considered in this study.
As mentioned earlier, the research on the healthcare logistics is very rich and thematically, we can divide the literature of these studies into two categories.One type of the studies covers medical and healthcare resources and services in hospitals to provide service at the right time [3,50,59,81].In the second type of studies, the patients stay at home to receive health services known as home care services, which may be used for elderly and old people to keep and help them at their homes [5,28] instead of hospitals, or public health services.Recently, many optimization models and solutions were devoted to address this challenging decision-making problem for routing and scheduling of caregivers and patients [29].In this study, we develop a mathematical model to handle health services in hospitals for curing of patients who do not live in the host country.The aim of the travelers not only is to cure their disease, but also have the leisure time in the attractive host country.
Despite the rapid growth and increasing interest of researchers and experts in this field, most articles related to this area fall in the qualitative and statistical categories, while providing mathematical models in medical tourism problems is quite scarce.For instance, Ahmadimanesh et al. [2] proposed a mathematical model in supply chain network design for the medical tourism problem.Their proposed model determined the optimal number of treatment units and final capacity for accommodation.Rezaeiahari and Khasawneh [65] developed a mathematical model for medical tourists' scheduling in destination medical centers to minimize deviations from patients' preferred time from the beginning of treatment.They assumed that service providers were not available during certain periods.Rezaeiahari and Khasawneh [66] proposed another mathematical model to minimize the difference between patients' scheduling and patients' preferred times and minimize patients' flow times in medical centers.They proposed a multi-period medical tourism problem and assumed that the duration of treatment is uncertain.To solve the problem, they applied a simulation approach with two optimization algorithms, Tabu Search (TS) and Simulated Annealing (SA).Motevalli-Taher and Paydar [56] developed a multi-objective mathematical model to reduce the outbreak of Corona Virus during Covid-19 pandemic.They modeled a mathematical model to reduce the number of patient from outside of city to prevent from the spread of virus to other cities.
As reported in Table 2, the development of a medical and/or entertainment tour design and trip planning problem in an uncertain environment, especially at operational and tactical levels is rarely studied in the literature.The base formulation of the proposed problem can be considered as the integration of orienteering problem (OP) and traveling salesman problem (TSP) as two well-studied combinatorial optimization problems.
The OP is a challenging problem considering a certain number of nodes with a specified score.OP is generally aimed at finding the path with the maximum specified length to maximize the scores by visiting several nodes.Therefore, the OP is the combination of the Knapsack problem (KP) and TSP [79].This problem is mentioned with other titles, such as the selective TSP, bank robber problem, and maximum collection problem [80], in the literature.Team orienteering problem (TOP) is an extension of OP whose goal is determining routes to maximize the total collected scores in a limited time.A team gathers ratings by visiting nodes in the same period [76].Many scholars have studied OP, TOP, and the development of these models [73,80].
Some scholars have reviewed and classified various studies and applications concerning OP and TOP, such as the tourist trip design problem [36].Yu et al. [84] developed TOP considering time-window and multiple transportation modes in the tourist trip design application.Lin and Vincent [45] investigated a mathematical model with some mandatory nodes that must be visited as crucial nodes.Abbaspour and Samadzadegan [1] presented a mathematical model for the OP of a passenger in an individual journey.The time dependency and multi-modal travel vehicles were considered in the tourist touring.In addition, in order to solve the tourist trip planning problem and find the shortest trajectory, genetic algorithm (GA) was employed with variable chromosome lengths.Salazar-Aguilar et al. [69] addressed a model, which included the team navigation problem for scheduling different tasks.In the provided model, there were two categories of nodes, namely obligatory and arbitrary ones.The purpose of the problem was to maximize the total collected score of non-compulsory points in the team.Zheng et al. [91] modeled an individual daily tourist travel planning and hotel allocation intending to maximize the passenger utility function, assuming that the utility rate decreases over time.They applied a heuristic algorithm to solve the problem.Zheng and Liao [90] proposed a mathematical model with distinct nodes in a tourist personal travel path problem.The problem was developed as a bi-objective model, which was aimed at maximizing the total group utilities as the first objective function and maximizing the minimum utilities gained from other group members as the second objective function.
Feillet et al. [30] introduced TSP with profits.They surveyed this problem, which does not require visiting all nodes, but its goal is to ensure the amount of collected profit.Another consideration in this regard is to maximize the objective function, which restricts the total cost of the journey.Another procedure in a practical tour problem is to consider the difference between travel cost and collected profit [34].Tang and Miller-Hooks [77] investigated the probability of the prior tour problem considering the uncertainty in service times, travel times, and travel costs to improve each customer service.The problem was modeled by chance-constraint stochastic programming and solved using exact and heuristic approaches.Evers et al. [26] addressed an OP in which the weight of each arc was stochastic with hard capacity constraints.The importance of stochastic parameters affected the expected value of the total benefit collected in the nonlinear model.Angelelli et al. [4] integrated stochastic routing problem and probabilistic TSP considering the cost of each arc and the benefit of each node.A mathematical linear integer model in a two-stage stochastic problem was solved using the branchand-cut algorithm and some heuristic approaches.Vincent et al. [82] proposed a TOP with time-dependent scores considering time windows and time budget in the tour trip design problem.
In the literature regarding the attractive OP, the tendency to select each location is determined by comparing the attraction probability of one location rather than the others [25].Freeman et al. [32] developed an attractive function to investigate the impacts of vicinity and time on concert tour planning.In an attractive TSP, the attraction is gained based on gravity function among different places to maximize the attraction of places [25].The model proposed in this article considers the optimal assignment of patients to hospitals and the determination of residence to suggest tour packages to customers among various hospitals and attractive places based on their interests.Competitive attractiveness among hospitals is presented because of customer behavior by the Huff gravity function [32].The attractiveness of each hospital depends on the quality of the hospital's treatment services compared to other hospitals.Besides, the attractiveness of places is considered as a gravity function that varies with the amount of time each person resides there.
Table 2 summarizes some of the previous studies related to medical, tourism, and OPs.Accordingly, the literature on medical tourism, OP, and gravity function are reviewed simultaneously.The previous studies have rarely utilized mathematical models in the medical tourism problem.Very few studies have focused on the scheduling of medical tourism problems in medical centers [65,66], and none of these studies has been applied to consider gravity function based on the quality of medical services in medical tour planning.Despite the fact that the OP in tourist trip design has recently attracted the attention of researchers, few studies have considered the TOP in tourist trip designs [69,82,84], among which Vincent et al. [82] has merely considered time dependency and time budget in their model.The attractive tourist tour problems have been studied to some extent [32,90,91].However, the issue of assigning and orienteering medical tourists, especially at operational and tactical levels, has not been explored using mathematical models so far.Generally speaking, considering the tourism industry problems integrated with the health system in the mathematical model has received scant attention in the academic world.To bridge the existing gaps, this paper has proposed a mathematical model to assign and schedule patients in medical and pleasure tour planning based on the quality of medical services and time-dependent gravity functions.Individual and public opinions are also considered in the model to integrate the behavior of customers in the trip design.Despite the inevitable importance of uncertainty, it is repeatedly neglected in the medical and tourism planning problems [4,26,66].Exploring the related literature shows that the duration of treatment is strongly stochastic [66], due to therapeutic risks.Customers' interests in visiting tourist spots is another source of uncertainty, which directly influences the quality of a tour.Hence, a two-stage stochastic optimization framework is employed to investigate the medical-pleasure tour design in an uncertain environment in which the therapeutic and entertainment tour planning is integrated with the gravity function.The gravity function enables the model to increase the attraction of the journey with the help of a higher level of trip customization (e.g., residence).
Based on the above explanations, the integration of the medical-pleasure tour into TOP to assign and schedule travelers, consideration of gravity function both in medical and tourism centers, the assumption of uncertainty in the duration of treatment and weights of entertainment places, and proposing a hybrid solution approach, collectively constitute a significant contribution in this context.

Problem definition
A host country in the medical-entertainment tourism industry should have many potential attractive places for tourists, as well as low expenses of healthcare, eminent doctors at the international level, prosperous progress in surgical procedures (such as heart disease treatment, cosmetic surgery, and infertility treatment), and short waiting time for allocation of patients to hospitals.Besides, for an attractive tourism trip design, companies try to simultaneously diminish overall expenses and satisfy customers' expectations in both medical and entertainment dimensions.Therefore, a crucial issue in an MTC is a recommendation system to design a memorable trip plan for passengers.Furthermore, for efficient management, strengthening the infrastructure of medical tourism, as well as assigning and transferring patients, are essential measures in medical tourism.
The integrated medical-entertainment tour design is a rather novel business in the medical tourism industry, suggesting suitable travel plans to patients during their travels.The objective of a travel agency is to maximize the total profit of the company according to customers' preferences in a journey.The attraction in travel acquires both through satisfaction in therapeutic services and attractive places.
The present study proposes a two-stage stochastic programming framework.Under this framework, the decisions are categorized into two classes, with and without scenario indices.In the first stage, decisions are made before stochastic events, and no scenario is taken into account.The second stage is the reaction to the stochastic events to mitigate the risk of the decisions made at the first stage [52].Therefore, decisions are Table 2.The contribution of this study versus other studies.decisive at the first stage.At the second stage, the outcomes depend on different events, which have taken place based on the likelihood of scenarios [9].In this study, the first stage involves the assignment decision variables to allocate patients to hospitals with limited capacities.At the second stage, all decision variables, including the logistics, transportation, tour planning, and other operational decisions, depend on the scenario that happens.Another issue which is considered is that customers have a limited time.Patients are allowed to stay in the host country according to the type of visa issued.Hence, depending on the percentage of recovery gained in the hospital, patients wish to attend various entertainment tourist programs regarding the required time and tourism areas' attraction.The timing of a journey, therefore, depends on the probability of improvement over which the scenarios are defined.
Given the uncertain parameters through a pre-estimated set of discrete scenarios, the model aims to determine an optimal detailed schedule of a medical-entertainment tour for each patient to meet different objectives of an MTC.The main objective of MTC is to maximize the total profit obtained from the total revenue of medical and entertainment tours minus the total costs of treatment, transportation, visiting, and accommodation.In order to embed the patients' interests and their expectations of the whole journey, which are a combination of satisfaction with the therapeutic and tourism sectors, two new constraints are introduced.The first one depends on the performance of a hospital among other suggested hospitals based on the calculation of average public opinions extracted from public rating systems.The second one relies on each patient's residence in each attractive place and is identified by a time-dependent gravity function in the proposed model.In the following sections, gravity functions will be explained more.

Medical tour:
− Each patient is introduced as a tuple (patient number and home country) with a unique ID. − Each patient travels from the home country to the host country to treat an illness.− At the beginning of a journey, the medical tour commences, and the treatment process begins to cure the disease.− Each patient is allocated to only one hospital for the treatment.− The capacity of each hospital is limited.− The attractiveness of each hospital is determined by its performance compared to other hospitals.

Transportation:
− The planning of travel and transportation initiates from the home country.− At the end of a journey, the patient must return to their home country.− The travel costs, from the home to the host country and vice versa, from the hospital to touristy cities, and between touristy cities are included in the transportation costs.

Entertainment tour:
− An entertainment tour begins after finishing the treatment process at the hospital.− In the trip planning, each patient should visit at least one node.− The resident time (days) a person stays at each node is a decision variable, and the model determines it according to the proposed gravity function.− An MTC should undertake transportation costs, treatment costs, and accommodation costs throughout the time.The first two equations are the components of the first stage, consisting of the total revenue of treatment (RT) and total transportation cost from home countries to hospitals.The next ones are the components of the second stage, including the total revenue of visiting sites (RS), the total cost of treatment (CT), the total cost of transportation (TC) except for the total cost of travel from the home country to the host country, and the total cost of accommodation (CA), under different scenarios.The total revenue of treatment and total revenue of visiting sites are formulated in equations (3.1) and (3.2), respectively.Equation (3.3) gives the total cost of treatment under scenario , which depends on the medical time.Equation (3.4) expresses the sum of transportation costs under scenario  from the origin node to the hospital, from the hospital or site node to another site node (for sightseeing), and the return of the patient from the site to the origin node.Equation (3.5) calculates the total accommodation cost under scenario , which depends on the patient's residence in each touristy node.The difference between the revenues and costs under all possible scenarios is formulated as the expected profit function.Equation (3.6) maximizes the profit function in the medical tour problem.

Gravity function
In medical-entertainment tour design, a fixed number of hospitals is considered, and patients should visit one of them.The allocation of patients to hospitals depends on the relative utility function of each hospital compared to others, computed by average feedback of patients cured in them [71].Hence, the gravity function is applied to determine the general relative attraction of each hospital.Skellern [71] suggested an additive utility function of hospitals from each patient's perspective, which considers the hospital-specific and patient-specific options.Therefore, the gravity function of hospitals and the patients' preferences are considered such that the customers' expectations are fulfilled more than those of people in previous experiences.In equation (3.21),  ℎ and  ℎ represent the attractiveness of hospital ℎ versus other hospitals and the general gravity function of hospital ℎ, respectively.
In this paper, the gravity function of tourist attractiveness is formulated as a nonlinear function.The attraction of tourist site visits increases with a descending slope over time.Therefore, the negative exponential distribution function in Figure 1, which is commonly applied to estimate the duration of a particular event, is used in this study according to equation (3.22).The attraction of each tourism node for each passenger depends on how long it takes a passenger to stay in that node.It is calculated using the exponential cumulative distribution function in equation (3.23), as shown in Figure 2. In this figure,  = 1 means that the tourist will remain at least a unit.Equation (3.24) shows the time-dependent gravity function of attraction in a tourist site at a selected node.

Linearization of proposed model
It is advisable to convert the nonlinear model to a linear one since it improves the problem-solving speed and efficiency remarkably.The attraction of tourist sites is formulated as a nonlinear function.The attraction of a journey depends on the residence time the tourists stay at the nodes.In literature, there is a well-known linearization technique named the piecewise linear method, which transforms a nonlinear function into a linear one [53].The rise in the number of segmentations makes the estimation of the model more accurate while increasing the computational complexity.Figure 3 demonstrates the nonlinear gravity function of visiting places against its linear approximation.
For the linearization of the nonlinear function, equation (3.25) is used to calculate the slope of each line (  ) and estimate the value of gravity function.
According to Figure 3,   <  −1 , therefore, the gravity function increases with a descending slope.One of the significant contributions of this model is the variable residence time in each node, which is realistically

Solution approach
Due to complexity theory for NP-hard optimization problems [80]; nowadays, there is a great deal of interest in development of heuristics and meta-heuristics for decision-making problems like online learning [89], scheduling [6,27,49], multi-objective optimization [46], healthcare [20,23], transportation [22,62], etc.In a similar way, this study develops a hybrid meta-heuristic algorithm for solving the stochastic medical tourism tour problem.In this regard, a Progressive Hedging Algorithm (PHA) is applied to tackle the uncertainty by a decomposition approach to analyze each stochastic scenario.In the proposed hybrid meta-heuristic, the GA is an approach to solve sub-problems while PHA decomposes the main model into different sub-problems.
The main reasons to justify the selection of GA and PHA are first to address the uncertainty in an efficient way while solving the main model effectively in large-scale networks.In the literature review, many studies [35,41,64] addressed the uncertain decision-making problems of PHA.As such, the evolutionary algorithms and different variations of GA were employed successfully in the literature on routing and OP optimization [12,19,61,78].These facts motivate our efforts in the development of a hybrid meta-heuristic based on the GA and PHA for the proposed stochastic tourism routing problem.

PHA for stochastic mixed-integer program
Two-stage stochastic programming is a branch of linear stochastic programming.In this approach, the model is run for stochastic variables in each scenario, considering the occurrence probability of the scenarios [63].The two-stage stochastic mixed-integer programming is commonly applied as follows: in which the probability of each scenario is specified by vector .Therefore, variable  is determined before the uncertainty is realized.Then, the model decides based on the occurrence probability of events to maximize the expected value of profits in the second stage.Equations (4.4)-(4.6)present the recourse action to react to the stochastic events in the second stage.
Step 6. Stop algorithm if the values of first stage variables in all scenarios are equal; otherwise, repeat steps 2-6.

GA
GA is a well-known random search algorithm where main idea is inspired by nature.In nature, better next generations are created by the composition of parents with better chromosomes, which is simulated in GA by the crossover operator.Mutations in chromosomes may also happen that probably improve the next generation [21].
In order to find the best possible solution in large-sized problems at acceptable times, it is essential to solve the model with heuristic, meta-heuristic, or hybrid algorithms.In this paper, an efficient hybrid algorithm is developed, which combines PHA and GA to cope with the uncertainty and complexity of the proposed model simultaneously.

Initial population
To generate an initial population according to the number of patients, closed routes are created randomly, which include some cells, a home country, a hospital, and a variable number of tourism sites.Moreover, given the variable length of nodes, the residence period is allocated corresponding to each tourist site.The set of solutions to a problem, which are obtained either randomly or algorithmically, is called the initial population, in which the number of preliminary solutions equals the population size.Each member of the initial population is named a chromosome.Each chromosome includes a set of genes.Figure 4 shows the representation of a chromosome for our GAs.In this figure, the trip starts from the home or origin country, where the origin country is country number 2, and terminates to the same country.The hospital number 3 is randomly selected for the treatment of patient.One of the important parts of this chromosome is the tourism places that the patient should be visited.The number of the places the tourism that can be visited is variable from one to total number of tourism places.Therefore, the tourism places are coded as cells so that the number of places can be different.As shown in Figure 4, tourism places 3 and 4 are selected for visiting among 5 different places.In the end, the patient comes back to the original country.For each tourism site, which is selected, the staying time is determined randomly.Hence, the tourism stays 1-and 3-time unit to visit tourism places of 3 and 4, respectively.

Crossover
The crossover operation is performed in each iteration in the following way.Parents are selected based on the probabilities of crossover.In crossover operation, two offspring inherit part of the genes from two selected parents.The number of tourism sites in each chromosome is variable.Therefore, a sub-route is selected randomly from each parent of crossover.The sub-route includes at least one tourist site and at most the maximum number of tourist sites of the selected parent.Figure 5 represents an example of the crossover operator, in which the sub-routes of 1-5, and 3 are selected randomly from the parents 1 and 2, respectively.Then, each selected subroute is inserted into the other parent, and the duplicated places are removed.For instance, in the sub-route of parent 1, node 5 is the duplicated node in parent 2, therefore; it is deleted, and only node 1 is added to offspring 2. To find the best possible place for a sub-route, a greedy heuristic approach, the best-insertion method, is used.This method determines both the route in which the sub-route is inserted and the two consecutive tourist sites between which the sub-route is inserted. 1 and   denote the first and last tourist sites in the subroute, while ℎ  and ℎ +1 represent tourist sites in a route in the offspring.The following explains how to calculate the payoff of inserting a sub-route between sites [10]: where (ℎ  • ℎ +1 ) denotes the profit of traveling between two consecutive sites.The algorithm searches through the whole offspring and inserts the sub-route in the place with the largest payoff value.Therefore, a new offspring is generated, and this operation is used again for the other parent to generate the second offspring.As shown in Figure 5, the best place for replacing the sub-route of parent1 is between nodes 2 and 5 of offspring2.Also, the best place for putting the sub-route of parent2 is after node 2 in offspring1.

Mutation
The mutation operator provides the property of being random and escaping from local optimum points.In this problem, the candidates for mutation are selected randomly based on the occurrence probability of mutation.In the proposed algorithm, mutation is performed on the assignment of a patient to a hospital, which is swapped with the random selection among hospitals.Besides, mutation changes the tourism route.Hence, one point of visiting places is randomly picked and deleted; if all the tourist points in the set are selected, only we remove the chosen node.Otherwise, on the unselected nodes, some are randomly selected and replaced with the taken nodes.In Figure 6, we represent an example of a mutation in the proposed algorithm.For the selected chromosome, the hospital number is randomly changed from hospital 3 to hospital 1.In addition, for the visiting route, which is 3-4, the visiting node of 3 was selected and removed.Also, among 3 unselected nodes which are 1, 2, and 5; first, two nodes are randomly chosen from 3 nodes.In Figure 6, as seen by the orange color, nodes 5, and 1 were randomly selected and added.We mention that the number of nodes, which we add was randomly selected.
Finally, the initial population, crossover, and mutation are combined and sorted according to the fitness function (objective function).Arranged members are selected as long as the original population and the rest of the population are removed.The steps of the proposed GA are described below: Step 1. Set parameters for the model.
Step 2. Generate a primary population with random chromosomes.
Step 3. Evaluate the initial population with fitness function.
Step 4. Sort population based on the value of fitness function.

Loop:
Step 5. Select parents for crossover.
Step 6. Generate offspring by applying crossover operation.
Step 7. Evaluate offspring via the fitness function.
Step 8. Select parents for mutation.
Step 9. Create new solutions by applying mutation on selected parents.
Step 10.Calculate the value of fitness function for mutant population.
Step 11.Merge the last population, crossover, and mutation.
Step 12. Sort merged population and truncate it to the number of population size.
Step 13.Go to Step 5 until the ending conditions are satisfied.

Hybrid algorithm
Given the complexity of the large-sized problem, a meta-heuristic approach is applied.In order to simplify the scenario-based problem, a hybrid PHA-GA algorithm is proposed.The problem is first reformulated and broken down into  (number of scenarios) sub-problems under the PHA algorithm, in which a GA is the core solver.In this algorithm, two GAs are coded to solve the problem.One (GA1) is based on the fitness function of the first stage of PHA, while the other (GA2) is formed based on the second stage.Hence, the first solution is processed in the GA1, and the convergence occurs on the GA2 to reach the best solution.All of the constraints and equations in the model except constrains (3.8), (3.17) and (3.18) are considered in the chromosomes to provide the whole of search space.However, there are some constraints to be added to the model to confirm the feasibility of the model.Hence, in the proposed algorithm, the constraints (3.8), (3.17) and (3.18) are added as a penalty costs to the fitness function to satisfy the feasibility of the solutions.If the constraints do not fulfill, the value of fitness function increases with penalties and does not provide feasible solution.The steps of the hybrid PHA-GA algorithm are explained below.
Step 1. Set the first fitness function in the GA based on the first stage of PHA.
Step 2. Run the GA to find the best solution for each scenario, which is applied as the initial population of the next GA.
Step 3. Set the fitness function of GA based on the penalty function given in the second stage of the PHA.
Step 4. Update the penalty function and repeat step 3 and 4 until the terminate condition is fulfilled.

Taguchi
The effectiveness of meta-heuristic calculations incredibly relies upon the right selection of parameters.For tuning the parameters, several strategies are accessible to effectively structure the experiments.One of the techniques, which is broadly utilized in most researches, is full factorial design.This technique tests all mix of conceivable variables.In each case, when the quantity of components altogether expands, this strategy is not cost-effective since the number of required trials is growing considerably [7].Taguchi experimental design technique is recommended to lessen the number of investigations.For a large number of decision variables,   Taguchi method determines orthogonal arrays to decrease the number of experiments.To demonstrate the variant of the response variable in finding the best level of each factor in experiments, the proportion of signal to noise (S/N) is calculated as follows:
In this paper,   and  point to the objective function in experiment  and the number of replication, respectively.The main parameters of the proposed PHA-GA algorithm that need to be tuned are as follows: the population size of GA1 and GA2, maximum iteration in the GA1 and GA2, and the coefficient () in the second GA.Table 3 describes the factors and levels for tuning parameters in the GA-PHA algorithm.In this study, we consider five factors in three levels; therefore, this method reduces 3 5 experiments to orthogonal array L 27 with 27 experiments.Table 4 shows the 27 experiments, which indicates the levels of each factor in each experiment and control factor as a result.
To gather reliable data, each experiment was performed three times and the average of the outputs is considered as a result.The main effects of means and S/N ratio are depicted in Figure 8.The results of the selected criteria for calibration are reported in Table 5.

Computational results and sensitivity analysis
To survey the computational experiment on medical-entertainment tour design, the proposed model is examined based on a real example of central medical services in the Middle East, and the problems with various sizes are solved.We code the proposed model under the OPL and scripts accessed via IBM ILOG CPLEX version 12.8, and code the hybrid approach, which is the consolidation of PHA and GA via MATLAB R2019a.The codes are run in a computer with an Intel R ○Core TM i7-4500U CPU at 2.39 GHz and an 8.00 GB RAM.

Case description
In this paper, several patients are considered who are willing to travel abroad to receive a range of medical services.In addition, some tourist attractions are available, with which the patients can arrange a trip to have leisure time after treatment.An integrated system of medical tourism centers, including transportation companies, hospitals, and tour centers, wants to assign patients based on their preferences to serve the appropriate services.Therefore, patients are first assigned from their own country to the hospital in the host country to receive medical services.Then, after the treatment procedures are completed, one or more recreational tours are formed through which the passengers (patients) visit the touristy cities.This strategy is obtained through the coordination and integration of various elements in providing services.Therefore, the impacts of healthcare quality services, residence time in tourist sites, and customers' specifications in medical-pleasure tour planning are investigated as an integrated system, and the effects of integration on the companies' profitability are examined.Figure 9 illustrates the position of each patient, hospital, and tourist site in an example.In addition, the navigation of each patient in the medical-entertainment tour and the departure time of each traveler from the hospital and each touristy city are demonstrated.The capacity of hospitals, traveling cost, cure cost, and accommodation cost are known.The proposed two-stage stochastic model is run, whose detailed results are presented in Table 6.Table 7 describes the scheduling and sequence of visiting places on the trips.The arrival time of each patient under each scenario is also presented in Table 7.The results show that in medical tour planning, patients are assigned to hospitals for treatment based on the quality of medical centers and the preferences in the desired hospital, which has higher priority in planning than gaining profit from the medical tour.For example, the first preference of patient 3 is to be cured in hospital 3 in Mashhad.Moreover, patient 2 prefers to go to hospitals 2, 3, and 1, respectively.The tour plan shows that the treatment of patient 2 should be arranged in the first intended medical point of the customer.Therefore, respect to the customers' opinions and public ratings are important considerations in this study to provide high-quality medical tour services.On the other hand, in pleasure tour planning, if a patient is assigned to a touristy city, which is more attractive for him/her, the residence of the patient increases compared to the residence in a less interesting city from his/her perspective.For example, in Table 7, patient 1 should visit two tourist sites, Tehran and Esfahan, under both scenarios 1 and 2, but the residences under the two scenarios are different.In scenario 1, the residence in each site is very short, but in scenario 2, the residence is longer because the individual attraction of the two scenarios is disparate, and the interest of patient 1 to visit these cities in scenario 1 is lower than those in scenario 2. Consequently, the residence of a visitor depends on the rate of individual interest; if the interest in a node is higher, the residence   8, as well as the gap between the two procedures.The values of the recourse problem and PHA method are approximately equal with a difference of only 0.01.Table 9 lists the objective value of each problem without considering scenario reports separately.It also presents the results of the objective function in the wait and see (WS) problem, which is considered a problem without scenarios.Moreover, the average value of the objective function is computed by running each scenario, and the expected value of perfect information (EVPI) and the value of the stochastic solution (VSS) are calculated to show the necessity of using the stochastic method in the problem.Equations (5.1) and (5.

Performance of hybrid PHA-GA
To analyze the complexity of the proposed model, 15 test problems with different sizes were generated.To generate such random test problems, Table 10 shows the range of parameters for our test problems.
Figure 10 shows the convergence behavior of algorithm to reach the optimal solution.The model was run for a test problem by considering 4 origin countries, 8 hospitals, 10 tourism sides, and 40 patients in 100 iterations, after 20 iterations the algorithm achieved to the steady state condition.
To compare the effectiveness of the proposed algorithm with that of the exact method, 15 various instances with small, medium, and large sizes are investigated.Table 11 shows the 15 different sizes of the problem and represents each instance's objective function and computational time, which are yielded by solving the model using PHA in CPLEX and hybrid PHA-GA.In Table 11, the second to the fifth columns represent the dimensions of the test problems.Moreover, the sixth and seventh columns provide the results of the objective function and computational time obtained from PHA, while the eighth and ninth columns list the corresponding values yielded by hybrid PHA-GA.The last column of Table 11 shows the gap between the two approaches.The results show that in the medium-and large-sized test problems, the CPLEX is not able to gain the feasible solution in a reasonable time.Hence, for medium and large-scale problems, employing hybrid PHA-GA results in a better objective function and lower computational time compared to the PHA.
The computational times of 15 instances solved by CPLEX, GA, and hybrid PHA-GA are illustrated in Figure 11.In small-and medium-sized problems, the CPLEX solved problems in more reasonable times compared to others methods.In medium-sized problems, the computational time of GA and hybrid PHA-GA are approximately equal.Applying hybrid PHA-GA in large-scale problems results in more reasonable computational time than the other methods.Therefore, in large-scale problems, the hybrid PHA-GA is a more efficient approach than GA by attaining solutions in a shorter computational time.In small-and medium-sized test problems, the difference between GA and hybrid PHA-GA is not noteworthy, and in the small-sized problems, the proposed approach finds the optimal solutions with negligible gaps.
To represent the reliability of the algorithm, in Table 12, we show the best, worst, and average solutions, as well as the standard deviation of the results after 10 times run.To check the results, we compared the results of CPLEX with the suggested algorithm.It was found that the gap values did not exceed 3.4% for the considered algorithm, which addresses the accuracy.Table 13 illustrates the details about the performance of hybrid PHA-GA versus GA.A comparison of the objective functions and computational times shows that the value of the objective function of WS solved by GA is more than that solved by hybrid PHA-GA in each test problem.Furthermore, the computational time of GA is significantly longer compared to the hybrid PHA-GA.Hence, the hybrid PHA-GA finds the solution in a shorter computational time.The last row of Table 13 shows the average gap of objective functions and computational times, which are 42.5% and 20.8%, respectively.Therefore, applying hybrid PHA-GA considering  scenarios is highly reasonable in this study because it provides low-cost, high-quality solutions compared with calculating the expected value of solutions by the iterative running of GA for each scenario.Figures 12 and 13, respectively, illustrate the results of the objective function in each iteration for the test problems 4 and 5, which are run under three scenarios by hybrid PHA-GA with the occurrence probability of 0.15, 0.8, and 0.05, respectively.The results indicate the variability of the objective function in each iteration to converge the average value of the scenarios to the RP problem.In Figure 12, the stop condition is fulfilled after ten iterations, and the average value of all scenarios is very close to the RP.Furthermore, Figure 13 shows the variability of the objective function of each scenario in each iteration to converge the expected value of scenarios to the value of RP.

Sensitivity analysis
Figure 14 shows the percentages of costs, revenues, and profits in medical and tourist travel in various test problems with different dimensions.The results indicate that the percentage of medical and tourism costs increases with the test problems' size, while the percentage of transportation costs decreases.Increasing the problems' dimension results in higher percentages of tourism revenues and lower percentages of medical revenues.Besides, the ratio of tourism profits to medical profits increases by more than 9%, which indicates that the augmentation of the number of tourists, in terms of economies of scale, has a significant impact on the rise in the percentage of tourism profits.
The sensitivity analysis is also performed on medical and tourism preferences.In Figure 15, the horizontal axis indicates the percentage of treatment preferences of previous visitors and the percentage of individual preferences in tourism trip to decide on patients' travel plan.By considering the higher percentage of customers'    16.The effect of scenario on the objective function.
both in therapeutic attractiveness and in tourism, the profit of MTC is reduced.As shown in this figure, the medical preferences are not significant on profit but the reduction in profits by increasing the effect of tourism attractiveness is more significant than medical Therefore, increasing attention to travelers' preferences in MTC to provide medical services is justifiable.Although tourism services reduce total profit, it can be compensated by providing satisfaction to improve the market of medical tourism.Therefore, it is better the MTCs pay attention to the customers' preferences on medical services versus they should focus on the public attraction for recommending the visiting tourism sides.
The sensitivity analysis is also performed on the scenarios provided in Figure 16.As shown in the figure, due to the problem's stochastic nature, as the number of scenarios increases, the objective function grows with a decreasing slope.However, the number of scenarios does not significantly affect the objective function after considering ten scenarios, while it leads to a rise in the complexity and computational time for solving the problem.This figure also indicates that by increasing the number of scenarios, the objective function's value gradually becomes stable.

Conclusion and future studies
This is the first study that considers a medical tourism travel plan in the attractive TOP.A novel mathematical model was presented that integrated the medical, transportation, and tourism industries to maximize the profit of MTCs.In this problem, a two-stage scenario-based model was applied.In the first stage, based on the customers' expectations and feedback of previous patients, passengers were assigned to the hospitals.In the second stage, the decisions were made according to the probability of recovery in hospitals, which depended on the stochastic treatment time.Another novelty of the problem was to decide on the period of residence by a proposed time-dependent gravity function in each tourist site, which depended on the interests of travelers.A real example was investigated to design an integrated medical-entertainment tour for patients from different countries with various characteristics and interests.In small-and medium-sized problems, the model was solved by CPLEX, a meta-heuristic approach, which integrated PHA and was developed to solve the medium-and large-sized scenario-based problem in a reasonable time.The various instances were tested to validate the approach, the results showed that the average gap between exact and hybrid algorithm was 3% highlighting the accuracy.By comparing GA and hybrid PHA-GA to show the impact of scenarios, the results revealed that the average values of objective function and computation time were improved by 42.5% and 20.8%, respectively.Therefore, by using hybrid PHA-GA, the performance of the solution significantly enhanced because this method applied evolutionary features and decomposed the problem.
From the managerial insight, integrating the medical-entertainment tour design is justifiable and profitable for MTCs because the profit of pleasure tours in comparison with medical tours is higher.Therefore, the development of regions for establishing infrastructure to serve the tourists from the economic and social is reasonable.In addition, we advise to MTCs to pay attention to the advantage of the opinions of tourists in therapeutic tour to increase customer satisfaction and increase the market share.On the other hand, for pleasure tours, the public attractions provide profitability for MTCs.Furthermore, the profit of MTCs was compensated by increasing the number of customers in terms of economy of scale and aggregating the pleasure tour with medical services.Therefore, from a managerial perspective, in the long-term, the integration of pleasure tours with medical services promotes profit.
Although this study was among the first studies for the development of a stochastic medical tourism problem, there are some limitations that can be considered as the path of numerous studies in the future, -The model can be developed as a multi-objective model with more suppliers of different objectives, e.g., maximizing total score and minimizing the waiting time to plan the travel.-The model can be extended as a medical-pleasure tour design with other types of attractive utility functions, and the effects of each utility function can be investigated.-Other considerations, such as time window, time dependency, and multiple periods can be considered.

Figure 1 .
Figure 1.The distribution of gravity function.

Figure 2 .
Figure 2. The cumulative distribution of gravity function.

Figure
Figure 3. Linearization of gravity function.

Figure 6 .
Figure 6.A representation of mutation operator.

Figure 7
Figure 7 represents a flowchart of the hybrid algorithm, showing the sequence of steps in it.

Figure 7 .
Figure 7.The flowchart of the proposed hybrid algorithm.

Figure 9 .
Figure 9.The position of patients and the schematic navigation of patients in medicalentertainment tour in scenario 1.

Figure
FigureThe convergence behavior of algorithm.

Figure 11 .
Figure 11.The computational time of small-and large size instances.

Figure 12 .
Figure 12.The value of the objective function in different scenarios in hybrid PHA-GA in instance 4.

Figure 13 .
Figure 13.The value of the objective function in different scenarios in hybrid PHA-GA in instance 5.

Figure 14 .
Figure 14.The percentage of costs, revenues, and profits in medical and tourism industry in different instances.

Figure 15 .
Figure 15.The effect of medical and tourism preferences on profit.
should be after the traveling time of the patient to the host country.Constraint(3.14)expresses the scheduling of treatment when the patient is assigned to the hospital.Constraint (3.15) yields the arrival time, leaving time, and the residence time of each patient visits sites based on the gravity function.Constraint (3.16) specifies the maximum duration of the journey for each patient.Constraints (3.17) and (3.18) guarantee that the passengers' preferences are met in the hospital and tourist sites, respectively.Constraints(3.19)and (3.20) define the type of variables.
∑︁ ∈    =   ∀ ∈ , ∀ℎ ∈ , ∀ ∈  (3.9)  )︀ ∀ ∈ , ∀ ∈ , ∀ ∈  ∪ (),  ̸ = , ∀ ∈  (3.15) guarantees the allocation of each patient to a tourist site after being cured at the hospital.Constraint (3.10) represents the connectivity of each path and the balance of nodes.Constraint(3.11)indicates the allocation of each patient to each tourist site.Constraint (3.12) states that each patient at the end returns to the origin node.Constraint(3.13)shows that the beginning of the service time for treatment if the patient is assigned to the hospital 3. Linearization of gravity function.considereda decision variable.Therefore, equation (3.26) determines the time-dependent attraction of a visiting place once the node is visited.Equations (3.27)-(3.30)areadded for the linearization of the piecewise function.(ST  ) =   ST  +  (ST −1 )LST  < ST  ≤ UST .6) Equation (4.7) gives the expected value of the objective function in the second stage, considering the probabilities of scenarios.The common two-stage stochastic programming is defined as equations (4.7)-(4.11).
.11)Equations (4.12)-(4.17)are the extensive form of a two-stage stochastic program indicating the variables in the first stage depending on the scenarios.Moreover, equation(4.15)isadded to the model and introduced as non-anticipative constraints.These constraints express the independence of first-stage decision variables from the scenario.Owing to the complexity of the two-stage stochastic problem, the primary model is reformulated by denoting the first-stage decisions with scenario indices.The reason behind this reformulation is that the non-anticipativity constraints limit the problem to split scenarios from the first-stage variables.This way, the new formulation decomposes into  independent sub-problems, which are more tractable in terms of solving complexity.

Table 3 .
The factors and levels of proposed hybrid algorithm.

Table 4 .
The orthogonal array L 27 for experiment.

Table 5 .
The best value of parameters of hybrid algorithm gained from Taguchi.there increases.Therefore, MTCs tend to create memorable leisure tours for customers to visit different touristy cities based on their interests.The results of the objective function in the two-stage scenario-based stochastic model and PHA are presented in Table

Table 6 .
Assignment and navigation of patients in medical tour plan under different scenarios.

Table 7 .
The timing of each customer in journey under different scenarios.

Table 8 .
The results of objective function.

Table 10 .
Sources of the range of parameters.

Table 11 .
The performance of solution procedures with small-and large size instances.
Notes.: The number of origin counties for service; : The number of hospitals; : The number of tourist cities;  : The number of patients; OM: Out of memory.

Table 12 .
Results of different test problems.

Table 13 .
The performance of hybrid PHA-GA by considering scenario versus GA without scenario.