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

Authors: Said Ali Hassan; Yousra Mohamed Ayman; Khalid Alnowibet; Prachi Agrawal; Ali Wagdy Mohamed
Journal: Complexity
 (2020)

---

## Abstract

This article proposes a new problem which is called the Stochastic Travelling Advisor Problem (STAP) in network optimization, and it is defined for an advisory group who wants to choose a subset of candidate workplaces comprising the most profitable route within the time limit of day working hours. A nonlinear binary mathematical model is formulated and a real application case study in the occupational health and safety field is presented. The problem has a stochastic nature in travelling and advising times since the deterministic models are not appropriate for such real-life problems. The STAP is handled by proposing suitable probability distributions for the time parameters and simulating the problem under such conditions. Many application problems like this one are formulated as nonlinear binary programming models which are hard to be solved using exact algorithms especially in large dimensions. A novel binary version of the recently developed gaining-sharing knowledge-based optimization algorithm (GSK) to solve binary optimization problems is given. GSK algorithm is based on the concept of how humans acquire and share knowledge during their life span. The binary version of GSK (BGSK) depends mainly on two stages that enable BGSK for exploring and exploitation of the search space efficiently and effectively to solve problems in binary space. The generated simulation runs of the example are solved using the BGSK, and the output histograms and the best-fitted distributions for the total profit and for the route length are obtained.


---

## Body

## 1. 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. The target of the TAP is to select the route with maximum profitability for the advisory company.The 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. The 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. Therefore, 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 routeIn 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.The 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.The rest of this paper is organised as follows: Section2 describes in detail the mathematical model of the Travelling Advisor Problem (TAP) including the definition of problem variables, constraints, and the objective function. The proposed model is a nonlinear binary model with a dimension depending on the number of candidate workplaces.The Stochastic Travelling Advisor Problem (STAP) simulation procedure is explained in detail inSection 3. The 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 inSection 4; the case study is implemented in a very important service sector which is the occupational health and safety. This section reveals also the role and advices offered by the Occupational and Health Advisor to organizations and workplaces.InSection 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. These 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.


## 2. Mathematical Model for the Travelling Advisor Problem (TAP)

The new Travelling Advisor Problem (TAP) is defined on a graphG 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]). The 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.The mathematical model for such a deterministic TAP with known and fixed parameters representing the travelling and advising times is formulated as follows.


### 2.1. Decision Variables

Let(1)xim=1,if workplaceiis visited by the advisor on positionmof his route,iandm=1,2,...,n=total number of workplaces,0,otherwise.


### 2.2. Constraints

#### 2.2.1. Position Constraints

Each positionm in the advisor route has at most one workplace:(2)∑i=1nxim≤1,m=1,2,...,n.


#### 2.2.2. Workplace Constraints

Each workplacei can be in one position of the advisor route or not visited:(3)∑m=1nxim≤1,i=1,2,…,n.


#### 2.2.3. 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:(4)∑i=1nxim+1≤∑i=1nxim,m=1,2,…,n−1.(i)
If∑i=1nxim+1=1, then ∑i=1nxim=1,m=1,2,...,n−1(ii)
If∑i=1nxim+1=0, then there is no restriction on the value of ∑i=1nxim=1,m=1,2,…,n−1


#### 2.2.4. Working Hours/Day Constraints

The total time spent by the advisor in travelling and site inspection should be within the maximum working hours/day =T = 8 hours:(5)∑i=1nt0,ixi1+∑m=1n∑i=1ntixim+∑i=1n∑j=1j≠inti,j.∑m=1n−1xim.xjm+1≤T,where(i)
t0,i  = transportation time between the headquarter and workplace i, where i = 1, 2, …, n(ii)
ti,j  = transportation time between the two adjacent workplaces i and j, where i, j = 1, 2, …, n(iii)
ti  = inspection time for workplace i, where i = 1, 2, …, nThis 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.


#### 2.2.5. Binary Constraints

All the decision variables are 0-1:(6)xim=0or1,i,m=1,2,…,n.


#### 2.2.6. The 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:(7)Maxz=∑m=1n∑i=1npixim,in which pi  = profit of visiting workplace i, where i = 1, 2, …, n.The optimum solution will produce two distinct situations:(1)
If∑m=1nxim=1i = 1, 2, …, n, then all the n-workplaces are visited by the advisor in one working day and the problem is completed.(2)
If∑m=1nxim=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.


## 2.1. Decision Variables

Let(1)xim=1,if workplaceiis visited by the advisor on positionmof his route,iandm=1,2,...,n=total number of workplaces,0,otherwise.


## 2.2. Constraints

### 2.2.1. Position Constraints

Each positionm in the advisor route has at most one workplace:(2)∑i=1nxim≤1,m=1,2,...,n.


### 2.2.2. Workplace Constraints

Each workplacei can be in one position of the advisor route or not visited:(3)∑m=1nxim≤1,i=1,2,…,n.


### 2.2.3. 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:(4)∑i=1nxim+1≤∑i=1nxim,m=1,2,…,n−1.(i)
If∑i=1nxim+1=1, then ∑i=1nxim=1,m=1,2,...,n−1(ii)
If∑i=1nxim+1=0, then there is no restriction on the value of ∑i=1nxim=1,m=1,2,…,n−1


### 2.2.4. Working Hours/Day Constraints

The total time spent by the advisor in travelling and site inspection should be within the maximum working hours/day =T = 8 hours:(5)∑i=1nt0,ixi1+∑m=1n∑i=1ntixim+∑i=1n∑j=1j≠inti,j.∑m=1n−1xim.xjm+1≤T,where(i)
t0,i  = transportation time between the headquarter and workplace i, where i = 1, 2, …, n(ii)
ti,j  = transportation time between the two adjacent workplaces i and j, where i, j = 1, 2, …, n(iii)
ti  = inspection time for workplace i, where i = 1, 2, …, nThis 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.


### 2.2.5. Binary Constraints

All the decision variables are 0-1:(6)xim=0or1,i,m=1,2,…,n.


### 2.2.6. The 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:(7)Maxz=∑m=1n∑i=1npixim,in which pi  = profit of visiting workplace i, where i = 1, 2, …, n.The optimum solution will produce two distinct situations:(1)
If∑m=1nxim=1i = 1, 2, …, n, then all the n-workplaces are visited by the advisor in one working day and the problem is completed.(2)
If∑m=1nxim=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.


## 2.2.1. Position Constraints

Each positionm in the advisor route has at most one workplace:(2)∑i=1nxim≤1,m=1,2,...,n.


## 2.2.2. Workplace Constraints

Each workplacei can be in one position of the advisor route or not visited:(3)∑m=1nxim≤1,i=1,2,…,n.


## 2.2.3. 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:(4)∑i=1nxim+1≤∑i=1nxim,m=1,2,…,n−1.(i)
If∑i=1nxim+1=1, then ∑i=1nxim=1,m=1,2,...,n−1(ii)
If∑i=1nxim+1=0, then there is no restriction on the value of ∑i=1nxim=1,m=1,2,…,n−1


## 2.2.4. Working Hours/Day Constraints

The total time spent by the advisor in travelling and site inspection should be within the maximum working hours/day =T = 8 hours:(5)∑i=1nt0,ixi1+∑m=1n∑i=1ntixim+∑i=1n∑j=1j≠inti,j.∑m=1n−1xim.xjm+1≤T,where(i)
t0,i  = transportation time between the headquarter and workplace i, where i = 1, 2, …, n(ii)
ti,j  = transportation time between the two adjacent workplaces i and j, where i, j = 1, 2, …, n(iii)
ti  = inspection time for workplace i, where i = 1, 2, …, nThis 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.


## 2.2.5. Binary Constraints

All the decision variables are 0-1:(6)xim=0or1,i,m=1,2,…,n.


## 2.2.6. The 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:(7)Maxz=∑m=1n∑i=1npixim,in which pi  = profit of visiting workplace i, where i = 1, 2, …, n.The optimum solution will produce two distinct situations:(1)
If∑m=1nxim=1i = 1, 2, …, n, then all the n-workplaces are visited by the advisor in one working day and the problem is completed.(2)
If∑m=1nxim=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.


## 3. Stochastic Travelling Advisor Problem Simulation Procedure

The 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.The 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]. The 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]. The 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. The most famous application is the distribution of activity times in project networks like that of the advising time in STAP [8].The probability density function of beta distribution is given by(8)probx|α,β=xα−1∗1−xβ−1Bα,β,where the beta function is given by(9)Betaα,β=∫01tα−11−tβ−1dt.Alpha (α) and beta (β) can be any positive value greater than zero. They are called the shape parameters of the beta distribution [6]. These 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]. There 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. This method is known as program evaluation and review technique (PERT) analysis or PERT method [11]. The PERT distribution 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]. Therefore, PERT has become a classic tool for estimating uncertain activities durations [13].The 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]):(10)μ=16a+4m+b,σ2=136b−a2.The PDF of PERT distribution can be symmetric, right-skewed, or left-skewed according to the values of the parametersa, m, and b [16]. This is better clarified with visual illustrations of the PERT Distribution using special computer programs (for example, Rosenmai [10] and EasyFit [17]).


### 3.1. PERT-Beta Distributions

The 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]:(11)α=μ−ab−aμ−ab−μσ2−1,β=b−μb−aμ−ab−μσ2−1.The 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. The 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(RAND,α,β,a,b) function [19]. In Excel, the RAND() function is used to obtain uniformly distributed probability values. The 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. The 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. The steps of the solution procedure are shown in Figure 1.Figure 1
Steps of the simulation procedure of the stochastic TAP.


## 3.1. PERT-Beta Distributions

The 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]:(11)α=μ−ab−aμ−ab−μσ2−1,β=b−μb−aμ−ab−μσ2−1.The 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. The 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(RAND,α,β,a,b) function [19]. In Excel, the RAND() function is used to obtain uniformly distributed probability values. The 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. The 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. The steps of the solution procedure are shown in Figure 1.Figure 1
Steps of the simulation procedure of the stochastic TAP.


## 4. 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].The 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. They 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 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].The 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. The five factories shown in Table1 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.Table 1
Workplaces for the case study.
#FactoryAddressProfit ($) × (102)1El Nasr Co. For Intermediate ChemicalsKm 28, Cairo Alex Desert Rd., Industrial Zone, Abou-Rawash, Giza. Beside Smart Village3.12Eva PharmaEl-Sadat, Zawya Abou Muslim, Giza Governorate4.73Eastern Company450 Al Haram, Al Omraneyah Al-Sharqeyah, El-Omraniya, Giza Governorate5.54El Helal and Golden StarIndustrial Zone A, 10th of Ramadan, Cairo Governorate4.25Juhayna GroupSecond 6th of October Industrial Zone, Giza Governorate3.5Figure 2
Location of the candidate workplaces.


### 4.1. Estimating the Probabilistic Times

The 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. Figure3 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). The 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.Figure 3
Parameters of the stochastic travelling time.Table 2
PERT-beta distribution parameters of the travelling times from the starting location.
t0,iambμσαβt0,12235453.1250.2083.853.08t0,22240452.10.154.321.41t0,3610124.16670.254.232.33t0,470901055.1830.31673.883.03t0,54560651.5750.0924.331.67Similarly, 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]. The 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.Table 3
PERT-beta distribution parameters of the advising times.
tiambμσαβt12.75343.130.041.274.29t21.8522.752.10.021.014.22t33.7545.254.170.061.014.22t44.656.55.180.11.354.31t51.451.521.580.010.443.89Figure 4
Examples of beta BDF and CDF for travelling and advising times. (a) BDF fort0,1; (b) BDF for t0,2; (c) BDF for t0,3; (d) BDF for t0,4; (e) BDF for t0,5; (f) CDF for t0,1 and t0,2; (g) BDF for t1; (h) CDF for t1.
(a)(b)(c)(d)(e)(f)(g)(h)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 Table4 and 10 examples of the generated travelling times between workplaces are shown in Table 5.Table 4
An example of the generated random values for travelling and advising times.
Workplaceti,j (hours)12345ti (hours)00.590.620.21.611.053.04110.780.681.581.251.97520.741.590.993.8731.361.394.77142.21.5295Table 5
Examples of generated travelling times between workplaces.
#t1,2t1,3t1,4t1,5t2,3t2,4t2,5t3,4t3,5t4,510.780.681.581.250.741.590.991.361.392.220.670.641.531.150.811.230.851.561.272.6030.620.641.541.240.781.410.811.251.492.2840.810.651.661.120.81.590.821.351.442.0450.700.721.541.280.791.220.991.541.442.1260.800.711.441.210.831.390.661.481.491.8670.800.581.531.300.751.350.911.421.482.3380.720.671.611.300.751.460.891.631.352.1590.680.751.491.250.821.320.971.351.382.49100.760.741.561.190.821.600.971.491.432.19The simulation process is repeated 30 times, the generated random values for the used parameters according to the BERT-beta distribution for each simulation run are substituted in the mathematical model described inSection 2, and the optimum solution is obtained for each case for further evaluation.


## 4.1. Estimating the Probabilistic Times

The 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. Figure3 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). The 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.Figure 3
Parameters of the stochastic travelling time.Table 2
PERT-beta distribution parameters of the travelling times from the starting location.
t0,iambμσαβt0,12235453.1250.2083.853.08t0,22240452.10.154.321.41t0,3610124.16670.254.232.33t0,470901055.1830.31673.883.03t0,54560651.5750.0924.331.67Similarly, 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]. The 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.Table 3
PERT-beta distribution parameters of the advising times.
tiambμσαβt12.75343.130.041.274.29t21.8522.752.10.021.014.22t33.7545.254.170.061.014.22t44.656.55.180.11.354.31t51.451.521.580.010.443.89Figure 4
Examples of beta BDF and CDF for travelling and advising times. (a) BDF fort0,1; (b) BDF for t0,2; (c) BDF for t0,3; (d) BDF for t0,4; (e) BDF for t0,5; (f) CDF for t0,1 and t0,2; (g) BDF for t1; (h) CDF for t1.
(a)(b)(c)(d)(e)(f)(g)(h)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 Table4 and 10 examples of the generated travelling times between workplaces are shown in Table 5.Table 4
An example of the generated random values for travelling and advising times.
Workplaceti,j (hours)12345ti (hours)00.590.620.21.611.053.04110.780.681.581.251.97520.741.590.993.8731.361.394.77142.21.5295Table 5
Examples of generated travelling times between workplaces.
#t1,2t1,3t1,4t1,5t2,3t2,4t2,5t3,4t3,5t4,510.780.681.581.250.741.590.991.361.392.220.670.641.531.150.811.230.851.561.272.6030.620.641.541.240.781.410.811.251.492.2840.810.651.661.120.81.590.821.351.442.0450.700.721.541.280.791.220.991.541.442.1260.800.711.441.210.831.390.661.481.491.8670.800.581.531.300.751.350.911.421.482.3380.720.671.611.300.751.460.891.631.352.1590.680.751.491.250.821.320.971.351.382.49100.760.741.561.190.821.600.971.491.432.19The simulation process is repeated 30 times, the generated random values for the used parameters according to the BERT-beta distribution for each simulation run are substituted in the mathematical model described inSection 2, and the optimum solution is obtained for each case for further evaluation.


## 5. Proposed Methodology

Metaheuristic algorithms have been developed to solve the complex optimization problem with continuous variables. Mohamed et al. [32] recently proposed a novel gaining-sharing knowledge-based optimization algorithm (GSK), which is based on the ideology of acquiring knowledge and sharing it with others throughout their lifetime. The 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]). The proposed methodology is described as follows.


### 5.1. Gaining-Sharing Knowledge-Based Optimization Algorithm (GSK)

A constrained optimization problem can be formulated mathematically as follows:(12)MinfX;X=x1,x2,…,xDim,s.to.giX≤0;i=1,2,…,m,X∈αp,βp;p=1,2,…,Dim,where f denotes the objective function; X=x1,x2,…,xDim are the decision variables; giX are the inequality constraints, and αpandβ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.The 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. The 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. These 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. These people have the experience to judge other people and can categorize them (good or bad). The 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. The initial population is created randomly within the boundary constraints as follows:(13)xtp0=αp+randpβp−αp,wheret is for the number of populations and randp 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:(14)DimJ=Dim×Genmax−GGenmaxk,DimS=Dim−DimJ,wherek>0 denotes the knowledge rate that controls the experience rate, DimJ and DimS represent the dimension for the junior and senior stages, respectively, Genmax 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. Thus, individuals are updated as follows:According to the objective function values, the individuals are arranged in ascending order. For everyxtt=1,2,…,NP, select the nearest best xt−1 and worst xt−1 to gain knowledge and also choose randomly xr to share knowledge. Therefore, 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. The updated individual can be determined as follows:The individuals are classified into three categories (best, middle, and worst) after sorting individuals into ascending order (based on the objective function values). Best individual =100p% (xbest), middle individual = Dim−2∗100p% (xmiddle), and worst individual = 100p% (xworst). For every individual xt, choose two random vectors of the top and bottom 100p% individual for gaining part and the third one (middle individual) is chosen for the sharing part. Therefore, the new individual is updated through the pseudocode presented in Figure 6.Figure 5
Pseudocode for junior gaining-sharing knowledge stage. Note:kf>0 is the knowledge factor.Figure 6
Pseudocode for senior gaining-sharing knowledge stage. Note:p∈0,1 is the percentage of best and worst classes.


### 5.2. 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. The working mechanism of BGSK is presented in the following sections.


#### 5.2.1. Binary Initialization

The initial population is obtained in GSK using equation (18) and it must be updated using the following equation for binary population:(15)xtp0=roundrand0,1,where the round operator is used to convert the decimal number into the nearest binary number.


#### 5.2.2. Binary Junior Gaining and Sharing Stage

The binary junior gaining and sharing stage is based on the original GSK withkf=1. The individuals are updated in original GSK using the pseudocode (Figure 2) which contains two cases. These two cases are defined for binary stage as follows:Case 1. When fxr<fxt: there are three different vectors (xt−1,xt+1,xr), which can take only two values (0 and 1). Therefore, a total of 23 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. The results of each possible combination are presented in Table 6.Subcase (a): ifxt−1 is equal to xt+1, the result is equal to xrSubcase (b): whenxt−1 is not equal to xt+1, then the result is the same as xt−1 by taking −1 as 0 and 2 as 1Table 6
Results of the binary junior gaining and sharing stage of Case 1 withkf=1.
xt−1xt+1xrResultsModified resultsSubcase (a)00000001111100011111Subcase (b)1001110121010−1001100Table 7
Results of the binary junior gaining and sharing stage of Case 2 withkf=1.
xt−1xtxt+1xrResultsModified resultsSubcase (c)1100311000110111000011−20Subcase (d)0000000100210010−100001−101010001001000110110101111110211011−10110121111111The mathematical formulation of Case 1 is as follows:(16)xtpnew=xr;ifxt−1=xt+1,xt−1;ifxt−1≠xt+1.Case 2. When fxr≥fxt: there are four different vectors (xt−1,xt,xt+1,xr), that consider only two values (0 and 1). Thus, there are 24 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.Subcase (c): ifxt−1 is not equal to xt+1, but xt+1 is equal to xr, the result is equal to xt−1Subcase (d): if any of the conditions arisext−1=xt+1≠xrorxt−1≠xt+1≠xrorxt−1=xt+1=xr, the result is equal to xt by considering −1 and −2 as 0, and 2 and 3 as 1The mathematical formulation of Case 2 is as follows:(17)xtpnew=xt−1;ifxt−1≠xt+1=xr,xt;otherwise,


#### 5.2.3. Binary Senior Gaining and Sharing Stage

The working mechanism of binary senior gaining and sharing stage is the same as the binary junior gaining and sharing stage with value ofkf=1. The individuals are updated in the original senior gaining-sharing stage using pseudocode (Figure 3) that contains two cases. The two cases were further modified for binary senior gaining-sharing stage in the following manner:Case 1. Whenfxmiddle<fxt: it contains three different vectors (xbest,xmiddle,xworst), and they can assume only binary values (0 and 1), and thus a total of eight combinations are possible to update the individuals. These total eight combinations can be classified into two subcases ((a) and (b)) and each subcase contains only four different combinations. The obtained results of this case are presented in Table 8.Subcase (a): ifxbest is equal to xworst, then the obtained results are equal to xmiddleSubcase (b): on the other hand, ifxbest is not equal to xworst, then the results are equal to xbest with assuming −1 or 2 equivalent to their nearest binary value (0 and 1, respectively)Table 8
Results of binary senior gaining and sharing stage of Case 1 withkf=1.
xbestxworstxmiddleResultsModified resultsSubcase (a)00000001111100011111Subcase (b)1001110121010−1001100Table 9
Results of binary senior gaining and sharing stage of Case 2 withkf=1.
xbestxtxworstxmiddleResultsModified resultsSubcase (c)1100311000110111000011−20Subcase (d)0000000100210010−100001−101010001001000110110101111110211011−10110121111111Case 1 can be mathematically formulated in the following way:(18)xtpnew=xmiddle;ifxbest=xworst,xbest;ifxbest≠xworst.Case 2. Whenfxmiddle<fxt: it consists of four different binary vectors (xbest,xmiddle,xworst,xt), and with the values of each vector, a total of sixteen combinations are presented. The sixteen combinations are also divided into two subcases ((c) and (d)). Subcases (c) and (d) further contain four and twelve combinations, respectively. The subcases are explained in detail in Table 9.Subcase (c): whenxbest is not equal to xworst and xworst is equal to xmiddle, then the obtained results are equal to xbestSubcase (d): if any case arises other than (c), then the obtained results are equal toxt by taking −2 and −1 as 0 and 2 and 3 as 1The mathematical formulation of Case 2 is given as follows:(19)xtpnew=xbest;ifxbest≠xworst=xmiddle,xt;otherwise.The pseudocode of BGSK is presented in Figure7.Figure 7
Pseudocode for BGSK.


### 5.3. Example of Results

The STAP is solved using the proposed novel BGSK algorithm and the values of parameters are presented in Table10. 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. Moreover, Figure 8 shows the convergence graph of the solutions of STAP using BGSK. From the figure, it can be observed that after the 25th iteration, it converges to the global optimal solution, which shows the robustness of the BGSK.Table 10
Numerical values of parameters.
Parameters of BGSKConsidered valuesNP200k10kr0.9p0.1kf1Maximum number of iterations100Figure 8
Convergence graph of BGSK.The obtained optimum solutions for each simulation run including total profit, route length (hours), and visited workplaces are presented in Table11. Figure 9 illustrates the histogram for the count of runs of the visited workplaces in the optimum solutions.Table 11
Optimum solutions for the 30 simulation runs.
#Total profit ($)Route length (hours)Visited workplaces110.27.042, 3210.27.745, 3310.26.703, 249.007.211, 2510.26.823, 269.007.013, 579.007.685, 389.006.933, 5910.26.833, 21010.26.713, 2118.205.012, 51210.27.333, 21310.28.003, 2148.205.295, 21510.26.743, 2169.007.233, 51710.26.563, 2189.007.113, 51910.26.912, 3208.607.713, 12110.27.113, 22210.27.172, 3239.007.575, 3249.007.523, 52510.26.693, 22610.26.813, 2279.007.413, 5289.007.855, 3299.007.785, 33010.27.433, 2Figure 9
Histogram for the count of runs of the visited workplaces in the optimum solutions.From Figure9, 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. The 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. The histogram for the total profit and the route length is obtained. Table12 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.Table 12
Goodness-of-fit summary for the total profit and the route length.
DistributionKolmogorov–SmirnovAnderson–DarlingChi-squaredStatisticRankStatisticRankStatisticRankTotal profitBeta0.47633311.88431.02221Normal0.3446223.746512.12892Pert0.533334171.24410.193Uniform0.2881319.80212N/ATriangularNo fitRoute lengthBeta0.1428810.5266110.603542Normal0.1844931.084230.190061Pert0.1721820.7523921.05593Triangular0.2660451.674445.06074Uniform0.2355848.92515N/AThe histograms and the best-fitted distributions for total profit and route length for the results are shown in Figure10. 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. The descriptive statistics for both distributions are presented in Table 13.Figure 10
Histogram and distribution fitting. (a) Total profit ($) and (b) route length (hours).
(a)(b)Table 13
Descriptive statistics for the total profit and route length.
StatisticTotal profitRoute lengthDistributionNormalBetaSample size3030Range22.99Mean9.57337.0633Variance0.504090.42964Std. deviation0.709990.65547Coef. of variation0.074160.0928Std. error0.129630.11967Skewness−0.43748−1.5419


## 5.1. Gaining-Sharing Knowledge-Based Optimization Algorithm (GSK)

A constrained optimization problem can be formulated mathematically as follows:(12)MinfX;X=x1,x2,…,xDim,s.to.giX≤0;i=1,2,…,m,X∈αp,βp;p=1,2,…,Dim,where f denotes the objective function; X=x1,x2,…,xDim are the decision variables; giX are the inequality constraints, and αpandβ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.The 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. The 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. These 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. These people have the experience to judge other people and can categorize them (good or bad). The 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. The initial population is created randomly within the boundary constraints as follows:(13)xtp0=αp+randpβp−αp,wheret is for the number of populations and randp 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:(14)DimJ=Dim×Genmax−GGenmaxk,DimS=Dim−DimJ,wherek>0 denotes the knowledge rate that controls the experience rate, DimJ and DimS represent the dimension for the junior and senior stages, respectively, Genmax 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. Thus, individuals are updated as follows:According to the objective function values, the individuals are arranged in ascending order. For everyxtt=1,2,…,NP, select the nearest best xt−1 and worst xt−1 to gain knowledge and also choose randomly xr to share knowledge. Therefore, 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. The updated individual can be determined as follows:The individuals are classified into three categories (best, middle, and worst) after sorting individuals into ascending order (based on the objective function values). Best individual =100p% (xbest), middle individual = Dim−2∗100p% (xmiddle), and worst individual = 100p% (xworst). For every individual xt, choose two random vectors of the top and bottom 100p% individual for gaining part and the third one (middle individual) is chosen for the sharing part. Therefore, the new individual is updated through the pseudocode presented in Figure 6.Figure 5
Pseudocode for junior gaining-sharing knowledge stage. Note:kf>0 is the knowledge factor.Figure 6
Pseudocode for senior gaining-sharing knowledge stage. Note:p∈0,1 is the percentage of best and worst classes.


## 5.2. 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. The working mechanism of BGSK is presented in the following sections.


### 5.2.1. Binary Initialization

The initial population is obtained in GSK using equation (18) and it must be updated using the following equation for binary population:(15)xtp0=roundrand0,1,where the round operator is used to convert the decimal number into the nearest binary number.


### 5.2.2. Binary Junior Gaining and Sharing Stage

The binary junior gaining and sharing stage is based on the original GSK withkf=1. The individuals are updated in original GSK using the pseudocode (Figure 2) which contains two cases. These two cases are defined for binary stage as follows:Case 1. When fxr<fxt: there are three different vectors (xt−1,xt+1,xr), which can take only two values (0 and 1). Therefore, a total of 23 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. The results of each possible combination are presented in Table 6.Subcase (a): ifxt−1 is equal to xt+1, the result is equal to xrSubcase (b): whenxt−1 is not equal to xt+1, then the result is the same as xt−1 by taking −1 as 0 and 2 as 1Table 6
Results of the binary junior gaining and sharing stage of Case 1 withkf=1.
xt−1xt+1xrResultsModified resultsSubcase (a)00000001111100011111Subcase (b)1001110121010−1001100Table 7
Results of the binary junior gaining and sharing stage of Case 2 withkf=1.
xt−1xtxt+1xrResultsModified resultsSubcase (c)1100311000110111000011−20Subcase (d)0000000100210010−100001−101010001001000110110101111110211011−10110121111111The mathematical formulation of Case 1 is as follows:(16)xtpnew=xr;ifxt−1=xt+1,xt−1;ifxt−1≠xt+1.Case 2. When fxr≥fxt: there are four different vectors (xt−1,xt,xt+1,xr), that consider only two values (0 and 1). Thus, there are 24 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.Subcase (c): ifxt−1 is not equal to xt+1, but xt+1 is equal to xr, the result is equal to xt−1Subcase (d): if any of the conditions arisext−1=xt+1≠xrorxt−1≠xt+1≠xrorxt−1=xt+1=xr, the result is equal to xt by considering −1 and −2 as 0, and 2 and 3 as 1The mathematical formulation of Case 2 is as follows:(17)xtpnew=xt−1;ifxt−1≠xt+1=xr,xt;otherwise,


### 5.2.3. Binary Senior Gaining and Sharing Stage

The working mechanism of binary senior gaining and sharing stage is the same as the binary junior gaining and sharing stage with value ofkf=1. The individuals are updated in the original senior gaining-sharing stage using pseudocode (Figure 3) that contains two cases. The two cases were further modified for binary senior gaining-sharing stage in the following manner:Case 1. Whenfxmiddle<fxt: it contains three different vectors (xbest,xmiddle,xworst), and they can assume only binary values (0 and 1), and thus a total of eight combinations are possible to update the individuals. These total eight combinations can be classified into two subcases ((a) and (b)) and each subcase contains only four different combinations. The obtained results of this case are presented in Table 8.Subcase (a): ifxbest is equal to xworst, then the obtained results are equal to xmiddleSubcase (b): on the other hand, ifxbest is not equal to xworst, then the results are equal to xbest with assuming −1 or 2 equivalent to their nearest binary value (0 and 1, respectively)Table 8
Results of binary senior gaining and sharing stage of Case 1 withkf=1.
xbestxworstxmiddleResultsModified resultsSubcase (a)00000001111100011111Subcase (b)1001110121010−1001100Table 9
Results of binary senior gaining and sharing stage of Case 2 withkf=1.
xbestxtxworstxmiddleResultsModified resultsSubcase (c)1100311000110111000011−20Subcase (d)0000000100210010−100001−101010001001000110110101111110211011−10110121111111Case 1 can be mathematically formulated in the following way:(18)xtpnew=xmiddle;ifxbest=xworst,xbest;ifxbest≠xworst.Case 2. Whenfxmiddle<fxt: it consists of four different binary vectors (xbest,xmiddle,xworst,xt), and with the values of each vector, a total of sixteen combinations are presented. The sixteen combinations are also divided into two subcases ((c) and (d)). Subcases (c) and (d) further contain four and twelve combinations, respectively. The subcases are explained in detail in Table 9.Subcase (c): whenxbest is not equal to xworst and xworst is equal to xmiddle, then the obtained results are equal to xbestSubcase (d): if any case arises other than (c), then the obtained results are equal toxt by taking −2 and −1 as 0 and 2 and 3 as 1The mathematical formulation of Case 2 is given as follows:(19)xtpnew=xbest;ifxbest≠xworst=xmiddle,xt;otherwise.The pseudocode of BGSK is presented in Figure7.Figure 7
Pseudocode for BGSK.


## 5.2.1. Binary Initialization

The initial population is obtained in GSK using equation (18) and it must be updated using the following equation for binary population:(15)xtp0=roundrand0,1,where the round operator is used to convert the decimal number into the nearest binary number.


## 5.2.2. Binary Junior Gaining and Sharing Stage

The binary junior gaining and sharing stage is based on the original GSK withkf=1. The individuals are updated in original GSK using the pseudocode (Figure 2) which contains two cases. These two cases are defined for binary stage as follows:Case 1. When fxr<fxt: there are three different vectors (xt−1,xt+1,xr), which can take only two values (0 and 1). Therefore, a total of 23 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. The results of each possible combination are presented in Table 6.Subcase (a): ifxt−1 is equal to xt+1, the result is equal to xrSubcase (b): whenxt−1 is not equal to xt+1, then the result is the same as xt−1 by taking −1 as 0 and 2 as 1Table 6
Results of the binary junior gaining and sharing stage of Case 1 withkf=1.
xt−1xt+1xrResultsModified resultsSubcase (a)00000001111100011111Subcase (b)1001110121010−1001100Table 7
Results of the binary junior gaining and sharing stage of Case 2 withkf=1.
xt−1xtxt+1xrResultsModified resultsSubcase (c)1100311000110111000011−20Subcase (d)0000000100210010−100001−101010001001000110110101111110211011−10110121111111The mathematical formulation of Case 1 is as follows:(16)xtpnew=xr;ifxt−1=xt+1,xt−1;ifxt−1≠xt+1.Case 2. When fxr≥fxt: there are four different vectors (xt−1,xt,xt+1,xr), that consider only two values (0 and 1). Thus, there are 24 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.Subcase (c): ifxt−1 is not equal to xt+1, but xt+1 is equal to xr, the result is equal to xt−1Subcase (d): if any of the conditions arisext−1=xt+1≠xrorxt−1≠xt+1≠xrorxt−1=xt+1=xr, the result is equal to xt by considering −1 and −2 as 0, and 2 and 3 as 1The mathematical formulation of Case 2 is as follows:(17)xtpnew=xt−1;ifxt−1≠xt+1=xr,xt;otherwise,


## 5.2.3. Binary Senior Gaining and Sharing Stage

The working mechanism of binary senior gaining and sharing stage is the same as the binary junior gaining and sharing stage with value ofkf=1. The individuals are updated in the original senior gaining-sharing stage using pseudocode (Figure 3) that contains two cases. The two cases were further modified for binary senior gaining-sharing stage in the following manner:Case 1. Whenfxmiddle<fxt: it contains three different vectors (xbest,xmiddle,xworst), and they can assume only binary values (0 and 1), and thus a total of eight combinations are possible to update the individuals. These total eight combinations can be classified into two subcases ((a) and (b)) and each subcase contains only four different combinations. The obtained results of this case are presented in Table 8.Subcase (a): ifxbest is equal to xworst, then the obtained results are equal to xmiddleSubcase (b): on the other hand, ifxbest is not equal to xworst, then the results are equal to xbest with assuming −1 or 2 equivalent to their nearest binary value (0 and 1, respectively)Table 8
Results of binary senior gaining and sharing stage of Case 1 withkf=1.
xbestxworstxmiddleResultsModified resultsSubcase (a)00000001111100011111Subcase (b)1001110121010−1001100Table 9
Results of binary senior gaining and sharing stage of Case 2 withkf=1.
xbestxtxworstxmiddleResultsModified resultsSubcase (c)1100311000110111000011−20Subcase (d)0000000100210010−100001−101010001001000110110101111110211011−10110121111111Case 1 can be mathematically formulated in the following way:(18)xtpnew=xmiddle;ifxbest=xworst,xbest;ifxbest≠xworst.Case 2. Whenfxmiddle<fxt: it consists of four different binary vectors (xbest,xmiddle,xworst,xt), and with the values of each vector, a total of sixteen combinations are presented. The sixteen combinations are also divided into two subcases ((c) and (d)). Subcases (c) and (d) further contain four and twelve combinations, respectively. The subcases are explained in detail in Table 9.Subcase (c): whenxbest is not equal to xworst and xworst is equal to xmiddle, then the obtained results are equal to xbestSubcase (d): if any case arises other than (c), then the obtained results are equal toxt by taking −2 and −1 as 0 and 2 and 3 as 1The mathematical formulation of Case 2 is given as follows:(19)xtpnew=xbest;ifxbest≠xworst=xmiddle,xt;otherwise.The pseudocode of BGSK is presented in Figure7.Figure 7
Pseudocode for BGSK.


## 5.3. Example of Results

The STAP is solved using the proposed novel BGSK algorithm and the values of parameters are presented in Table10. 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. Moreover, Figure 8 shows the convergence graph of the solutions of STAP using BGSK. From the figure, it can be observed that after the 25th iteration, it converges to the global optimal solution, which shows the robustness of the BGSK.Table 10
Numerical values of parameters.
Parameters of BGSKConsidered valuesNP200k10kr0.9p0.1kf1Maximum number of iterations100Figure 8
Convergence graph of BGSK.The obtained optimum solutions for each simulation run including total profit, route length (hours), and visited workplaces are presented in Table11. Figure 9 illustrates the histogram for the count of runs of the visited workplaces in the optimum solutions.Table 11
Optimum solutions for the 30 simulation runs.
#Total profit ($)Route length (hours)Visited workplaces110.27.042, 3210.27.745, 3310.26.703, 249.007.211, 2510.26.823, 269.007.013, 579.007.685, 389.006.933, 5910.26.833, 21010.26.713, 2118.205.012, 51210.27.333, 21310.28.003, 2148.205.295, 21510.26.743, 2169.007.233, 51710.26.563, 2189.007.113, 51910.26.912, 3208.607.713, 12110.27.113, 22210.27.172, 3239.007.575, 3249.007.523, 52510.26.693, 22610.26.813, 2279.007.413, 5289.007.855, 3299.007.785, 33010.27.433, 2Figure 9
Histogram for the count of runs of the visited workplaces in the optimum solutions.From Figure9, 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. The 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. The histogram for the total profit and the route length is obtained. Table12 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.Table 12
Goodness-of-fit summary for the total profit and the route length.
DistributionKolmogorov–SmirnovAnderson–DarlingChi-squaredStatisticRankStatisticRankStatisticRankTotal profitBeta0.47633311.88431.02221Normal0.3446223.746512.12892Pert0.533334171.24410.193Uniform0.2881319.80212N/ATriangularNo fitRoute lengthBeta0.1428810.5266110.603542Normal0.1844931.084230.190061Pert0.1721820.7523921.05593Triangular0.2660451.674445.06074Uniform0.2355848.92515N/AThe histograms and the best-fitted distributions for total profit and route length for the results are shown in Figure10. 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. The descriptive statistics for both distributions are presented in Table 13.Figure 10
Histogram and distribution fitting. (a) Total profit ($) and (b) route length (hours).
(a)(b)Table 13
Descriptive statistics for the total profit and route length.
StatisticTotal profitRoute lengthDistributionNormalBetaSample size3030Range22.99Mean9.57337.0633Variance0.504090.42964Std. deviation0.709990.65547Coef. of variation0.074160.0928Std. error0.129630.11967Skewness−0.43748−1.5419


## 6. Conclusions and Points for Future Research Studies

The 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)
The stochastic nature appears in the presentation of the travelling and the advising times. The simulation procedure is used to generate stochastic parameters following PERT-beta distribution and then gather results obtained from the repeated simulation runs.(3)
The 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. The 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. The 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. The 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 factorkf=1. The 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)
The nonlinear binary mathematical model and the solution method are used to mimic and solve a real case study example. The 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)
The 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.The 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, education, telecommunications, investing, quality assurance, social and community services, pollution, medical, tourism, marketing, sales, advertising, sports, arts, and cooking.(4)
To check the performance of the BGSK algorithm in solving different types of complex optimization problems, and further works can be investigated by the extension of BGSK.



---

*Source: 6692978-2020-12-01.xml*