Stochastic Travelling Advisor Problem Simulation with a Case Study: A Novel Binary Gaining-Sharing Knowledge-Based Optimization Algorithm

,


Introduction
A new problem which we are going to call the Travelling Advisor Problem (TAP) in network optimization is defined for an advisor who wants to settle on the foremost profitable route for visiting some or all candidate workplaces each associated with a corresponding profit. He begins from a predetermined starting location and next wants to visit each chosen workplace exactly once within the day working hours. e target of the TAP is to select the route with maximum profitability for the advisory company.
e Travelling Advisor Problem is a new and different version of the well-known Travelling Salesman Problem (TSP) studied in operations research and theoretical computer science. e general idea of TSP is defined for a list of cities (represented as nodes) and their pairwise distances, and the task is to find the shortest possible tour that starts at the hometown, visits each city exactly once, and returns back home. Droste [1] stated that the number of different tours is very large, so one might not think to solve the problem by simply calculating the length of each possible tour. erefore, suitable algorithms are needed to solve these situations.
Instead, in the Travelling Advisor Problem (TAP), nodes represent a set of workplaces and an advisor's job is to visit a subset of the workplaces and spend some time in each for giving his advice. In quite a nutshell, an advisor starts his working day from his company's headquarter, and then he proceeds to travel to and advise a chosen workplace and then transfers to a new one until his working hours are over. During a workday of maximum eight hours, each travelling time between one workplace and the next in addition to the advising time spent at each visited workplace is considered in addition to the travelling time spent from the headquarter to the first visited workplace. Each workplace will yield a different profit than the other, meaning that the goal is to maximize the advisor's total profit.
TAP differs from TSP in the following main points: (i) In the TSP, the time is open till completing visiting of all customers, while in the TAP, the available time is limited by the day working hours (ii) In the TSP, the salesman will visit all the customers, while in the TAP, the advisor will determine a route containing some or all the workplaces which optimizes the problem objective function within the available limited day time (iii) In the TSP, no time is consumed in customer places, or it is immaterial, while in the TAP, the advising time at a workplace is a basic factor in the day time limit constraint, and hence, it directly affects the choice of the optimum solution of the problem (iv) In the TSP, the objective is to complete the route while minimizing the total travelling time, while in the TAP the objective function is to maximize the total profit in the chosen route In simple practical application problems where variations in travelling times and advising times are not so effective or the decision maker wants only to have a rough picture of the considered problem, these parameters can be considered as deterministic quantities. In general practice, the travelling and advising times are random and very difficult to be exactly measured or evaluated, considering that both the travelling and advising times should be independent and continuously distributed random variables. e stochastic version of the Travelling Advisor Problem has a wide range of real-life applications in various service and advisory fields. Such fields show up in many consulting domains such as health and safety, industry, agriculture, business, education, telecommunications, investing, quality assurance, social and community services, pollution, medical, tourism, marketing, sales, advertising, sports, arts, and cooking. e rest of this paper is organised as follows: Section 2 describes in detail the mathematical model of the Travelling Advisor Problem (TAP) including the definition of problem variables, constraints, and the objective function. e proposed model is a nonlinear binary model with a dimension depending on the number of candidate workplaces.
e Stochastic Travelling Advisor Problem (STAP) simulation procedure is explained in detail in Section 3. e new problem proposes PERT-beta probability distributions for the travelling and advising times and resolving the problem under such conditions, and the steps of the simulation solution procedure are also explained.
A real practical application case study is explained in Section 4; the case study is implemented in a very important service sector which is the occupational health and safety.
is section reveals also the role and advices offered by the Occupational and Health Advisor to organizations and workplaces.
In Section 5, a novel binary version of a recently developed gaining-sharing knowledge-based optimization algorithm (GSK) is introduced to solve the TAP. GSK cannot solve the problem with binary space; therefore, binary gaining-sharing knowledge-based optimization algorithm (BGSK) is proposed with two new binary junior and senior stages. ese stages allow BGSK to explore and exploit the search space of the problem efficiently; an example of the experimental results of one of the simulation runs of STAP is also presented. Section 6 gives the conclusions and the suggested points for future research studies, respectively.

Mathematical Model for the Travelling Advisor Problem (TAP)
e new Travelling Advisor Problem (TAP) is defined on a graph G with a set of n nodes V representing the workplaces (customers) and an additional node denoting the Health and Safety Agency headquarter (HQ) where the advisor starts the job and a set of arcs representing the travelling times between two distinct workplaces (Pinter [2]). e time of inspection in a workplace and the travelling time between two workplaces are specified.
While the advisor should start his route at the headquarter (HQ), he ends his route at the last visited workplace. e mathematical model for such a deterministic TAP with known and fixed parameters representing the travelling and advising times is formulated as follows.

Consecutive Position Constraints.
A position (m+1) cannot exist in the advisor tour unless the preceding position m exists, and this is achieved by the following set of constraints: where (i) t 0,i � transportation time between the headquarter and workplace i, where i � 1, 2, . . ., n (ii) t i,j � transportation time between the two adjacent workplaces i and j, where i, j � 1, 2, . . ., n (iii) t i � inspection time for workplace i, where i � 1, 2, . . ., n is is a quadratic inequality in two variables, the first part is for travelling from the headquarter to the first position in the route, and the second part is the inspection times at workplaces, and the third part is the travelling time between different positions in the route.

Binary
Constraints. All the decision variables are 0-1: 2.2.6. e Objective Function. It is formulated for maximizing the total profits of the Occupational Health and Safety Agency gained by visiting the workplaces during the working day time limit: in which p i � profit of visiting workplace i, where i � 1, 2, . . ., n. e optimum solution will produce two distinct situations: (1) If n m�1 x m i � 1 i � 1, 2, . . ., n, then all the n-workplaces are visited by the advisor in one working day and the problem is completed.
(2) If n m�1 x m i � 0 for any i, then the corresponding workplace i is not visited by the advisor in the first working day. In this case, it is needed to eliminate the visited workplaces, adding one more day and repeat the procedure for another working day.

Stochastic Travelling Advisor Problem
Simulation Procedure e problem so far has been concerned with deterministic travelling and advising times. However, solutions of the deterministic models deteriorate once applied in real-life problems where these times are stochastic [3]. In this section, the Stochastic Travelling Advisor Problem (STAP) will be handled, by proposing probability distributions for the time constraints and resolving the problem under such conditions. e probability distribution function (PDF) used for time simulations should be continuous and limited between two-time intercepts and have a unique mode in its defined range therefore making it ideal to choose here, as the beta PDF satisfies all the conditions stated [4]. Probability distribution of travelling times using beta distribution has been proposed and validated by many researchers [5]. e beta distribution is very flexible and commonly used to represent where the uncertain variable is a random value between zero and a positive value [6]. e beta distribution is flexibly described over the interval [0, 1] and the beta density function is a very flexible form to represent outcomes like probabilities [7]. Applications include modelling random variables that have a finite range a to b. e most famous application is the distribution of activity times in project networks like that of the advising time in STAP [8].
e probability density function of beta distribution is given by where the beta function is given by Complexity 3 Alpha (α) and beta (β) can be any positive value greater than zero. ey are called the shape parameters of the beta distribution [6]. ese two parameters work together to determine whether the distribution is symmetrical, positively skewed, or negatively skewed. Depending on the values of α and β, the beta distribution can take a variety of shapes [9]. ere are many graph packages that allow to experiment with different values of α and β and visualize how the shape changes (see, for example, Rosenmai [10]).
In PERT method, three-point estimate for activity duration is an estimate that includes optimistic, most likely, and pessimistic estimate. is method is known as program evaluation and review technique (PERT) analysis or PERT method [11]. e PERTdistribution is a family of continuous probability distributions defined by the minimum or best (a), most likely (m), and maximum or worst (b) values that a variable can take. It is a transformation of the four-parameter beta distribution with an additional assumption on its expected value.
PERT introduces uncertainty into the estimates for activity duration and is well suited for those situations where there is insufficient background information to specify accurately time. Its three-time estimates become the framework on which the probability distribution curve for the activity is erected [12]. erefore, PERT has become a classic tool for estimating uncertain activities durations [13]. e three estimates of the activity duration enable the expected mean time μ, as well as the standard deviation σ and variance σ 2 to be derived mathematically (Shankar et al. [14] and Golpîra [15]): e PDF of PERT distribution can be symmetric, rightskewed, or left-skewed according to the values of the parameters a, m, and b [16]. is is better clarified with visual illustrations of the PERT Distribution using special computer programs (for example, Rosenmai [10] and EasyFit [17]).

PERT-Beta Distributions.
e built-in beta distributions provided by most software systems are parameterized by two shape parameters (α and β) and two location parameters (a and b). In the PERT context, the desired beta distributions are specified by two statistics (mean and variance) and the same two location parameters (a and b). Hence, in order to carry out simulation with the software, one needs to first convert the PERT mean and variance and the extreme values a and b into the associated shape parameters α and β [18]. When these PERT mean and variance formulas are substituted into beta formulas, one gets the unique beta distribution parameters α and β needed for each time activity as a beta distribution which varies between (a and b) and fits for modelling the activity duration [13].
Beta distributions defined on the interval [a, b] in this fashion are referred as PERT-beta distributions because they exhibit means and variances as specified by the PERT mean and variance formulas. PERT-beta distributions have the following formulas for location parameters (a and b) and shape parameters (α and β) [19]: e best, most likely, and worst travelling times are estimated based on experience and with the help of Google Maps to get traffic data from the starting position to each candidate workplace and from each workplace to another in any given day and time. e minimum, most likely, and maximum advising times are also estimated based on the experience of the advising group.
PERT-beta distribution is used to generate random values for the travelling and advising times using the BETAINV(RAN D, α, β, a, b) function [19]. In Excel, the RAND() function is used to obtain uniformly distributed probability values. e generated random values for the travelling and advising times are introduced into the mathematical model, the model is solved using an appropriate method, and the results are gathered including the optimum solution, the route duration, and the expected total profit. e process of generating other random variables is repeated many times while the corresponding run results are gathered in each time. After enough number of runs, the final concluded results are worked out. e steps of the solution procedure are shown in Figure 1

Practical Application Case Study
One of the most important areas of the service sector is chosen in this case study, which is the occupational health and safety field.
Workers in all countries around the world are exposed to many hazards while working in their workplaces [20]. Among these hazards are the chemical hazards [21] and the biological hazards [22]. e purpose of health and safety advisor is to facilitate and promote safety, health, and wellness across all function areas of organizations, to collaborate with both staff and management to ensure regulatory compliance and to undertake reviewing, audits, and developing strategies for best practices [23].
Occupational Health and Safety (OH&S) Acts and Regulations are laws that govern workplaces. ey outline the rights and responsibilities of the employer, the worker, and the supervisor to ensure working environments are healthy and safe [24]. Egyptian law regulations in this field are enrolled under the International Labour Organization and the World Health Organization (WHO) [25,26].
It is the responsibility of the advising agency to check up the risks at the workplace and those that can happen, and 4 Complexity these inspections help distinguishing and preventing risk [27]. Health and safety advisors make sure that people comply with all rules of safety and health and that work environments are not the source of illness, wounding, or dying [28]. Advisors accomplish their work by visiting workplaces to ensure that work is accurately performed in compliance with regulations, Marriott et al. [29]. Other responsibilities are to predict possible hazards, providing advices on health and safety, recording violations of regulations, gathering illegal evidence, follow-up evidence in court, writing remarks and reports, and providing training and educational support to employers [30].
e Faculty of Engineering, Cairo University, located in Giza, Giza Governorate, is chosen as the starting position where the advising group will start the route as one of its research centers is responsible for this advising project. e five factories shown in Table 1 are the candidate workplaces to be visited, and the table shows also the profit for the visit of each workplace. Using a Google Maps tool "My Map," a location map is created to pinpoint the factory locations for better illustration as shown in Figure 2.

Estimating the Probabilistic Times.
e minimum (best), most likely, and maximum (worst) travelling times are calculated using Google Maps and past practical experience to get traffic data from the starting location to each workplace and from each workplace to another. Figure 3 shows an illustrative example for estimating the parameters (a and b) for the travelling time from the starting location to Workplace number 3 (Eastern Company). e most likely time m is estimated using experience, and it is in the interval [a and b] and depends mainly on the chosen day of the week and time of the day. In the same manner, other time parameters are estimated for all the travelling times of the case study. Table 2 shows the PERT-beta distribution parameters of the travelling times from the starting location to the workplaces.
Similarly, the parameters (a, m, and b) for PERT distribution of the advising time at each workplace are estimated based on the experience of the advisory group. Table 3 represents the PERT-beta distribution parameters of the advising times. Figure 4 expresses some examples of the beta BDF and CDF for travelling and advising times [31]. e distributions are right-skewed for all the five stochastic advising times, as m values are closer to the a values than they are to the b values.
Using such PERT-beta distributions, 30 simulation random values are generated using Excel software for all the travelling and advising times. An example of the generated values for the starting travelling and for the advising times is shown in Table 4 and 10 examples of the generated travelling times between workplaces are shown in Table 5. e simulation process is repeated 30 times, the generated random values for the used parameters according to  Complexity 5 the BERT-beta distribution for each simulation run are substituted in the mathematical model described in Section 2, and the optimum solution is obtained for each case for further evaluation.

Proposed Methodology
Metaheuristic algorithms have been developed to solve the complex optimization problem with continuous variables. Mohamed et al. [32] recently proposed a novel gainingsharing knowledge-based optimization algorithm (GSK), which is based on the ideology of acquiring knowledge and sharing it with others throughout their lifetime. e original GSK solves optimization problems over continuous space, but it cannot solve the problem with binary space. So, a new variant of GSK is introduced to solve the proposed TAP. A novel binary gaining-sharing knowledge-based optimization algorithm (BGSK) is proposed over discrete binary space with new binary junior and senior gaining and sharing stages. On the other hand, there are many constraint handling techniques in the literature (Deb [33]; Cello [34]; Muangkote et al. [35]. In this paper, the augmented Lagrangian method is used to handle the constraints, in which a constrained optimization problem is converted into an unconstrained optimization problem (Long et al. [36]; Bahreininejad [37]). e proposed methodology is described as follows.

Gaining-Sharing Knowledge-Based Optimization Algorithm (GSK).
A constrained optimization problem can be formulated mathematically as follows:   6 Complexity where f denotes the objective function; X � [x 1 , x 2 , . . . , x Dim ] are the decision variables; g i (X) are the inequality constraints, and α p and β p are the lower and upper bounds of decision variables, respectively, and Dim represents the dimension of individuals. If the problem is in maximization form, then consider minimization � −maximization.
e human-based algorithm GSK is of two stages: junior and senior gaining and sharing stages. All persons acquire knowledge and share their views with others. e people from early-stage gain knowledge from their small networks such as family members, relatives, and neighbours want to share their opinions with others who might not be from their networks, due to curiosity of exploring others. ese people may not have the experience to categorize the people. In the same way, the people from the middle or later age enhance their knowledge by interacting with friends, colleagues, social media friends, etc., and share their views with the most suitable person, so that they can improve their knowledge. ese people have the experience to judge other people and can categorize them (good or bad). e process mentioned above can be formulated mathematically in the following steps: Step 1: to obtain the starting solution of the optimization problem, the initial population must be obtained. e initial population is created randomly within the boundary constraints as follows: where t is for the number of populations and rand p denotes uniformly distributed random number between 0 and 1.
Step 2: at the beginning, the dimensions of the junior and senior stages should be computed through the following formula: where k( >0) denotes the knowledge rate that controls the experience rate, Dim J and Dim S represent the dimension for the junior and senior stages, respectively, Gen max is the maximum number of generations, and G denotes the generation number.
Step 3: junior gaining-sharing knowledge stage: in this stage, the early aged people gain knowledge from their small networks and share their views with the other people who may or may not belong to their group. us, individuals are updated as follows: According to the objective function values, the individuals are arranged in ascending order. For every x t (t � 1, 2, . . . , NP), select the nearest best (x t−1 ) and worst (x t−1 ) to gain knowledge and also choose randomly (x r ) to share knowledge. erefore, to update the individuals, the pseudocode is presented in Figure 5.
Step 4: senior gaining-sharing knowledge stage: this stage comprises the impact and effect of other people (good or bad) on the individual. e updated individual can be determined as follows: e individuals are classified into three categories (best, middle, and worst) after sorting individuals into ascending order (based on the objective function values). Best individual � 100 p% (x best ), middle individual � Dim − 2 * 100p% (x middle ), and worst individual � 100 p% (x worst ). For every individual x t , choose two random vectors of the top and bottom 100 p% individual for gaining part and the third one (middle individual) is chosen for the sharing part. erefore, the new individual is updated through the pseudocode presented in Figure 6.

Binary Gaining-Sharing Knowledge-Based Optimization Algorithm (BGSK).
To solve problems in discrete binary space, a novel binary gaining-sharing knowledge-based optimization algorithm (BGSK) was proposed. In BGSK, the new initialization and the working mechanism of both stages (junior and senior gaining-sharing stages) are introduced over binary space, and the remaining algorithms remain the same as the previous one. e working mechanism of BGSK is presented in the following sections.   8 Complexity

Binary Initialization.
e initial population is obtained in GSK using equation (18) and it must be updated using the following equation for binary population: x 0 tp � round(rand(0, 1)), (15) where the round operator is used to convert the decimal number into the nearest binary number.

Binary Junior Gaining and Sharing
Stage. e binary junior gaining and sharing stage is based on the original GSK with k f � 1. e individuals are updated in original GSK using the pseudocode (Figure 2) which contains two cases. ese two cases are defined for binary stage as follows: there are three different vectors (x t−1 , x t+1 , x r ), which can take only two values (0 and 1). erefore, a total of 2 3 combinations are possible, which are listed in Table 3. Furthermore, these eight combinations can be categorized into two different subcases ((a) and (b)) and each subcase has four combinations. e results of each possible combination are presented in Table 6.
Subcase (a): if x t−1 is equal to x t+1 , the result is equal to x r Subcase (b): when x t−1 is not equal to x t+1 , then the result is the same as x t−1 by taking −1 as 0 and 2 as 1 e mathematical formulation of Case 1 is as follows: Case 2. When f(x r ) ≥ f(x t ): there are four different vectors (x t−1 , x t , x t+1 , x r ), that consider only two values (0 and 1). us, there are 2 4 possible combinations that are presented in Table 7. Moreover, the 16 combinations can be divided into two subcases ((c) and (d)). Subcases (c) and (d) have four and twelve combinations, respectively.
the result is equal to x t by considering −1 and −2 as 0, and 2 and 3 as 1 e mathematical formulation of Case 2 is as follows:

Binary Senior Gaining and Sharing
Stage. e working mechanism of binary senior gaining and sharing stage is the same as the binary junior gaining and sharing stage with value of k f � 1. e individuals are updated in the original senior gaining-sharing stage using pseudocode (Figure 3) that contains two cases. e two cases were further modified      for binary senior gaining-sharing stage in the following manner: it contains three different vectors (x best , x middle , x worst ), and they can assume only binary values (0 and 1), and thus a total of eight combinations are possible to update the individuals. ese total eight combinations can be classified into two subcases ((a) and (b)) and each subcase contains only four different combinations. e obtained results of this case are presented in Table 8.
Subcase (a): if x best is equal to x worst , then the obtained results are equal to x middle Subcase (b): on the other hand, if x best is not equal to x worst , then the results are equal to x best with assuming −1 or 2 equivalent to their nearest binary value (0 and 1, respectively) Case 1 can be mathematically formulated in the following way: Case 2. Whenf(x middle ) < f(x t ): it consists of four different binary vectors (x best , x middle , x worst , x t ), and with the values of each vector, a total of sixteen combinations are presented. e sixteen combinations are also divided into two subcases ((c) and (d)). Subcases (c) and (d) further contain four and twelve combinations, respectively. e subcases are explained in detail in Table 9.
Subcase (c): when x best is not equal to x worst and x worst is equal to x middle , then the obtained results are equal to x best Subcase (d): if any case arises other than (c), then the obtained results are equal to x t by taking −2 and −1 as 0 and 2 and 3 as 1 e mathematical formulation of Case 2 is given as follows: e pseudocode of BGSK is presented in Figure 7.

Example of Results.
e STAP is solved using the proposed novel BGSK algorithm and the values of parameters are presented in Table 10. BGSK runs over personal computer Intel ® CoreTM i5-7200U CPU @ 2.50 GHz and 4 GB RAM and coded on MATLAB R2015a. To get the optimal solutions, 30 independent runs are performed.   x best x worst x middle Results Modified results Subcase (a) x best x t x worst x middle Results Modified results Subcase (c) 10 Complexity Moreover, Figure 8 shows the convergence graph of the solutions of STAP using BGSK. From the figure, it can be observed that after the 25 th iteration, it converges to the global optimal solution, which shows the robustness of the BGSK. e obtained optimum solutions for each simulation run including total profit, route length (hours), and visited workplaces are presented in Table 11. Figure 9 illustrates the histogram for the count of runs of the visited workplaces in the optimum solutions.
From Figure 9, it is resulted that visiting workplace numbers 3 and 2, respectively (Eastern Company and Eva Pharma), is the most occurring solution in 40% of the cases. e second occurring solution is to visit workplaces 3 and 5 with a chance of 20% only.
Histogram, goodness of fit, graph fitting, and descriptive statistics of the simulation results for the obtained are instructed. e histogram for the total profit and the route length is obtained. Table 12 shows the goodness-of-fit summary of the simulation results for the total profit and the route length. From the table and the shape of the PDF of the distributions, it is concluded that the goodness-of-fit test is in favor of the normal distribution for the total profit and the beta distribution for the route length. e histograms and the best-fitted distributions for total profit and route length for the results are shown in Figure 10.     A normal distribution is the most suitable to fit the total profit with a mean of $957.33, a minimum of $820, and a maximum of $1020. A beta distribution is the most suitable to fit the route length with a mean of 7.0633, a minimum of 5.01, and a maximum of 8 in hours. e descriptive statistics for both distributions are presented in Table 13.

Conclusions and Points for Future
Research Studies e main conclusions for this paper can be summarized as follows: (1) A new application problem called Stochastic Travelling Advisor Problem (STAP), where an advisor is likely to select the most profitable route for visiting a subset of candidate workplaces. Continuous probability distributions are used to model the stochastic variables for the problem in order to be a closer representation of real-world applications. (2) e stochastic nature appears in the presentation of the travelling and the advising times. e simulation procedure is used to generate stochastic parameters following PERT-beta distribution and then gather results obtained from the repeated simulation runs. (3) e problem arises extensively in the advising domain for the service sector. It differs from the deterministic version in the fact that the travelling and advising times are stochastic in nature which is the case in real-life problems. (4) Travelling Advisor Problem (TDP) looks like the famous Travelling Salesman Problem (TSP) but with basic distinct differences, mainly, the limited available route time, visiting only a subset of the places, considering the time consumed in customer places as a basic component, and the objective of maximizing the total profit. (5) A nonlinear binary mathematical model is formulated for the given stochastic problem. e binary decision variables represent the allocation of chosen workplaces into the positions of the proposed route, and the nonlinearity appears in the limited time constraint. e final objective is to maximize the total profit of the chosen route. (6) A real application case study in the field of occupational health and safety as one of the very important fields of the service sectors is presented. e stochastic model is formulated and the simulation run problems are solved, and the output histograms and the best-fitted distributions for the total profit and for the route length are obtained. (7) To obtain the solution of the proposed STAP nonlinear binary programming model, a novel binary gaining-sharing knowledge-based optimization algorithm (BGSK) is introduced which includes the two main binary junior and senior gaining-sharing stages with the knowledge factor k f � 1. e BGSK is the first variant of the developed GSK, and the proposed algorithm is applied with an augmented Lagrangian method to handle the constraints. (8) e nonlinear binary mathematical model and the solution method are used to mimic and solve a real case study example. e optimal solution is obtained which is more profitable by 10.8% than the solution given by the management based on selecting workplaces with the highest profit under the time limit constraint. (9) e obtained results by BGSK show its robustness and convergence and prove that it can find the global optimal solution of STAP. BGSK gives better results when compared to the health and safety agency management.
e points for future research studies can be stated as follows: (1) To propose other stochastic mathematical models' formulations for the same problem starting with the design of the decision variables and compare the effectiveness of computations for each model (2) To augment the proposed Stochastic Travelling Advisor Problem (STAP) with its variations as STAP with time window (STAPTW), STAP with multiple advisors (mSTAP), Multiobjective Stochastic Multiple Travelling Advisor Problem (MOSmTAP), and other variations. (3) To apply the same problem formulation to other similar advisory fields that can show up in many other consulting domains such as industry, agriculture, business,

Data Availability
Actually, data are available on request. e interested researchers can contact the corresponding author.

Conflicts of Interest
e authors declare that they have no conflicts of interest.