Accepted: June 12, 2024 | Published Online: June 14, 2024

Citation: Niasar MSF, Tayebi AH (2024) Reevaluating Safety and Health Concerns in Transportation for Students Commuting to School, by Exploring Various Metaheuristic Methodologies. Int J Ind Operations Res 7:020.

Copyright: © 2024 Niasar MSF, et al. This This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

The school bus routing problem (SBRP) represents a variant of the well-known vehicle routing problem that combines bus route generation and stop selection. As the organization in charge of student transportation, the municipality seeks to efficiently tackle the bus routing problem while coping with increased complexities and concerns in regard to resource, safety, and health issues. Although significant research has been conducted on SBRP, no study has so far considered these concerns in an integrated approach. Thus, the proposed approach was developed in two ways: First, a combination of mixed load and location-allocation-routing problems was considered. Second, the pre-processing approach was considered to identify the risk factors influencing student safety; correspondingly, the high level of risks was formulated in the model. Given the NP-hard complexity of proposed model, a metaheuristic approach was suggested to solve it for small, medium, and large instances. To this end, different metaheuristic configurations were examined that could explain the mechanism of diversification and structures of neighborhoods selection. The suggested configurations were addressed in m-VND (multi-start combined with VND), m-HA (multi-start combined with heuristic based adaptive layer), p-VND (perturb heuristic combined with VND), and p-HA metaheuristics (perturb heuristic combined with heuristic based adaptive layer).The results illustrated the better performance of the p-HA metaheuristic. Extensive computational experiments were conducted to investigate the effect of risk constraints on the objective function. The finding highlighted that imposing risk constraints contributes to increasing the objective function by 15.07%. Finally, some insights about health and safety risks were proposed that help better control operations and mitigate the encountered risks.

Keywords

SBRP, VND heuristic, p-HA metaheuristic, p-VND metaheurisitic, Safety issue 

Introduction

Transporting students to and from schools in Tehran (with a population of 10 million) involves scheduling daily activities of 4.5 million students. Due to different concentrations and dispersions of students throughout the city, the municipality seeks to implement an efficient bus routing plan. Any step in this regard must be in accordance with safety and health criteria, which include important issues that ought to be addressed in the transportation system. In Tehran, the majority of students go to and from school via the public transportation system, hence the need to make them aware of the risks associated with each mode of public transportation they use. The data from an accident report in Tehran (compiled by author) indicates that bus transportation, due to its strict policies, is safer than using private cars. These policies should also be applied to the school bus routing problem, including when the bus stops to pick up or drop off the students, the time considered for evacuating a bus because of an accident, as well as the time that the students are on the bus. These policies pave the way for establishing rules that facilitate carrying students to their respective schools in a safe and more convenient fashion. Technically, this problem is referred to as the road traffic safety. School buses have special features, including flashing red lights, cross-view mirrors, and stop-sign arms, that make them safer. They are also well equipped with protective seating and high crush standards. Additionally, the municipality is in charge of paying more compensation for students who are injured while using public transport. The system suggested by the French National Council of Transport (CNT) incorporates many criteria for assuring the safety of bus stops: the route leading to the bus stop, pedestrian crossing, visibility of bus stop for driver, size and position of bus stop, quality of place where students wait, and location of bus stop and maneuvers imposed on it. Another interesting work in this area is the model proposed by the Swedish Transport Agency, which calculates the risk of accidents and the level of insecurity for all students who are waiting at bus stops or are on the way to reach bus stops. They formed four kinds of bus stops based on the design and waiting area.

Safety matters are significantly associated with the type of vehicles (e.g., trucks passing students'), waiting/walking area and the local conditions (snow restrictions, etc.).

An important point in the school bus routing problem is to employ capable and physically healthy drivers who have not only perfect driving skills but also the ability to communicate to students safety tips as efficiently as possible. Non-compliance with some safety criteria in school routing problem is partly associated with families. Indeed, overcrowding of students on school buses has become a major concern in recent years, meaning that to reduce the transportation expense, parents allow theses buses to exceed their capacity constraint. The municipalities are also aware of the consequence of this tendency, yet they still give permission to bus drivers to pick up more students than their predefined capacity. Families play a critical role in this regard, as they can help to prevent and control this problem. There are other factors like depreciation and technical defects of the school bus fleet that make the problem even more complicated, and they need to be seriously considered. Parents are always concerned about reports of bus school accidents. This concern will be exacerbated when a given accident is caused by technical problems. A large percentage of Tehran's minibuses are worn out and cannot meet the necessary safety requirements for transferring students. In the period of 2017-2018, the Deputy of Transportation and Traffic of Tehran Municipality announced that all vehicles used for schools are required to undergo a valid technical inspection, and fleet failing this test are not deployed to provide services.

To date, this rule is enforced and all fleet transporting students are required to be accredited through the technical inspection centers. Based on the results of interviews with experts, around 70% of Tehran's minibuses are worn out and must be inspected under specific/technical requirements (compiled by the author), which substantiates that the current minibuses are not appropriate for transporting students. To handle this problem, the public-school transport system is entering in collaboration with the private sector (e.g., the taxi organization). This will help employ new kinds of vehicle in the student transportation system. The above subjects have received special attention in earlier studies, which helps alleviate parents’ worries.

More importantly, to be more efficient in meeting students’ safety and health requirements, it is needed that the families alongside the school committee pay special attention to these issues. This collaboration will lead to a reduction in the number of accidents. Families should ask and teach their children to wait for their bus in a safe place, such as a sidewalk, and open the door after making sure the vehicle stops. Families can remind their children to carefully listen to the advice of the driver and not to stand up from their seats when the vehicle is moving. Also, no student should speak loudly, which will distract the driver and result in an accident. Following these simple, but efficient, rule scan lead to mitigating the risks involved in student transportation. In addition to the steps taken already by the families, if municipality and relevant communities strictly follow their defined rules (e.g., supplying the needed resources), the occurrence of accidents can significantly decline.

There is limited research on student safety transportation. In this regard, Chalki, et al. [1] introduced a new method for transporting students to and from school via a safer and more efficient mechanism. This method was constructed based on devising a safety map that could assign the risk safety features to the considered arcs and nodes. The problem was investigated using different criteria such as traffic flow, speed, and lightening of road. To solve this method, three kinds of algorithms were proposed, of which the genetic algorithm proved superior to the others.

Kotoula, et al. [2] investigated a seven-step method for optimizing the school bus route of private school. They combined the cluster approach with the generic algorithm while taking into account the geographic features of the road network. The seven steps of this method included scenario definition, geocoding, calculating deriving distance, stop clustering, employing the genetic algorithm to develop the route, creating a route per school, and finally route visualization. The paper was aimed at reducing total travel distance and, thus, a significant decrease in NOX and CO ₂ emission.

Transportation network in Tehran can be listed in three different types of traffic areas: 1) Central restricted area, where only public transport (bus, taxi, ambulance, etc.( is allowed; 2) The area of even and odd traffic plan, or the so-called low emission zone (LEZ). (The even-and-odd traffic plan is an effective way to limit access to private vehicles. In order to implement this scheme, the private vehicles, depending on whether the last digit of their license plate is even or odd, can move in this area on alternate days of the week. According to newly approved by law, to decrease air pollution, all private vehicles are allowed to enter this area for a maximum of 10 working days a month). 3) The zone outside these two afore mentioned areas, called the free zone, which can be used by all kinds of vehicles. By enforcing this restriction, it will be possible to control vehicles’ passage from crowded and densely populated places and highways.

The prevalence of corona virus disease calls for more attention to students' health and safety. An alarming report in Iran indicated a case fatality rate of 5% (since the start of the pandemic until April 2021), which is the highest compared to other countries around the world. This finding illustrates that after the pandemic, it is imperative to adopt the right course of action in order to ensure students' safety when they attend the school. It also stresses the need to build an efficient school transport system that respects predefined safety and health concerns. Achieving these objectives requires providing a safety map for school bus routing and drafting appropriate guidelines.

In order to achieve beneficial result, the above issues need to be evaluated in accordance with their impact on safety and health issues, which is compatible with risk assessment. Thus, the risks need to be properly identified and, accordingly, the necessary action ought to be implemented to access them. To understand the real problems and concerns, we attempted to conduct a careful risk analysis of student transportation, which took account of different aspects (including traffic, population density, health, etc.). Risk assessment is one of the forthcoming methods to determine the level of safety in a complex and modern system. In practice, risk is uncertainty about the probability of a loss. In regard to safety, risk is considered as the probability of danger and its consequences.

So far, most studies have dealt with the application of risk assessment in urban transportation, resulting in its present popularity in this industry. Despite the existing knowledge, risk analysis of student transportation has remained in theorical level and has not been operationalized in an integrated fashion with the concept of mixed load and location-allocation in SBRP. In the present study, the risk assessment model was developed and implemented with respect to different features. These features were derived from previous knowledge of the expert team and the available data collected by the researcher. To do so, the potential risks were first identified. Then, risk prioritization was conducted to ascertain the importance of each risk and weight it. The results were used as a basis for the proposed model. Briefly, the contributions of this paper are four-fold:

1)       Providing a risk assessment method to identify and score risks that negatively impact students’ health and safety;

2)       Prioritizing the more threatening risks and formulating them in the model;

3)       Proposing a metaheuristic that investigates the type of diversification in the solution space; and

4)       Analyzing the total transportation cost while relaxing risk constraints in the model.

Literature Review

Limited research exists on student safety in transportation. Chalki, et al. [1] address this gap by introducing a new method for safer and more efficient transportation of students to and from school. The method involves creating a safety map that assigns risk safety features to considered arcs and nodes, considering criteria like traffic flow, speed, and road lighting. To solve this problem, three algorithms are proposed, with the genetic algorithm proving superior to the others. Fernandes, et al. [3] suggest in their review that school-based interventions focusing on the built environment may moderately enhance students' physical activity, health, and active commuting. Ensuring student safety and well-being during school bus transportation poses challenges and concerns for school boards, families, and transportation companies. The transportation of students with disabilities adds further complexities for transportation companies. Ross, et al. [4] conducted a scoping review on transportation services for students with disabilities. The review involves analyzing 20 documents from various perspectives, including the viewpoints of students with disabilities, stakeholders' understanding of essential education, disability rights legislation, challenges and concerns in school transportation, and alternative transportation options. The findings suggest a need for considering alternative practices. Strategies to enhance transportation for students with special needs may include consolidating transportation with other students, providing training en route to school, utilizing inclusive technology, instructing practical skills with medical professionals, and collaborating with non-transportation stakeholders for safe travel planning. In a cross-sectional study, Rothman, et al. [5] explore the correlation between built environment factors and Active School Transportation (AST) safety in seven Canadian communities. They analyze factors like child population density, multi-dwelling housing density, pre-1960 housing density, school design elements (e.g., parking, proximity, crossing guards, cycling infrastructure, sidewalks, car drop-off, traffic calming devices), and road design features (e.g., intersection density, traffic signal density, local road density). Social and environmental factors, such as school population and new immigrants, are also considered. The study suggests that crossing guards, cycling infrastructure, and traffic signs significantly contribute to enhancing AST safety, with their importance varying based on school type, travel culture, and city-specific challenges [5]. Implementing initiatives and new programs in drop-off zones can be a feasible solution to improve school routing safety. Looking at it differently, parents of disabled children may demand a robust response to perceived travel risks, expressing worries about their children's use and access to school transportation services. These concerns stem from issues such as non-compliance with safety regulations and protocols (Falkmer, et al. [6]), insufficient knowledge and skills of drivers, particularly with regard to students with disabilities, and inadequate seating arrangements on buses (Falkmer and Gregersen, [7]). Parents also express concern about the training and regulations provided to bus drivers, particularly their ability to assist disabled children with medical needs during transportation ([6]). Falkmer, et al. ([8]) highlight the primary risk as the potential for severe injury or death resulting from traffic accidents. To ensure passenger safety, factors such as driver expertise, knowledge, and abilities; safe locations for student pick-up and drop-off; efficient traffic management; and prevention of driver fatigue should be carefully considered. A Swedish organization, SAFEWAY2SCHOOL, has developed a technologically advanced program aimed at improving school travel safety. The program utilizes measures such as smart bus stops, GPS tags on buses, and onboard computers to monitor and control student safety comprehensively during bus trips (Falkmer, et al. [8]). Improvements are needed in various aspects of school travel safety, as highlighted by Dubée, et al. [9]. They stress the significance of maintaining appropriate bus schedules and establishing a clear emergency protocol for unforeseen events during school bus trips. Attention should also be directed towards addressing built environment issues, encompassing pick-up and drop-off procedures, bus-to-schoolyard arrangements, and individualized education programs (IEP) for each student. To tackle safety concerns effectively, it is crucial for school staff, boards, service providers, and bus drivers to enhance engagement with students and families, especially those who are vulnerable, and provide support to alleviate their concerns.

Research Method

The method of this research relies on four consecutive steps. In first stage, the main risks and sub-risks are identified and classified. In the second stage, the score of each identified risk is calculated based on the probability and level of impact and, subsequently, risk prioritization is conducted by making pair wise comparisons. In third stage, the risks with higher weights are fed to the model, such that the more threatening risks are assigned the highest priority. Finally, the proposed model is solving by an exact method for small instances and through metaheuristic configurations for all instances (see Figure 1).

The first stage involves constructing a method to identify and classify risks. Specifically, we attempt to identify the hazardous conditions triggered by defects, interaction, and interference of any factors that deviate from the desired objectives.

Risk identification in the risk assessment process is crucial, and failure to perform it will impair performance in the following stages.

Risk identification can be done systematically, experimentally, or inventively. In this study, it is carried out based on a combination of field observation and interview with experts. The data for the survey was gathered from 12 specialists who were systematically dispersed in different regions of the city. For this purpose, Tehran was divided into three zones (north, south, and center zones), and the relevant experts (4 experts per each zone) were interviewed in accordance with each zone.

The researcher went to different districts of city, (e.g., in bus stations), observed the important issues, and then drafted the initial list of hazards. The final list of hazards was prepared and categorized by combining experts’ insights and researcher’s observations. Totally, 13 risks were identified. To facilitate and better understand the identified risks, they were divided in three categories, namely safety, health, and traffic. To select more important risks, it is required to conduct the risk assessment (scoring the risk) and make a trade-off between identified risks and, accordingly, determines their respective weight (risk prioritization). Clearly, the risks with high weights are considered for further analysis. To sum up, in the second stage first we conduct the risk assessment (through the risk analysis) and then prioritize them. Risk prioritization is performed via the Analytic Hierarchy Process (AHP) technique. AHP provides a basis for finding effective solutions in complex decision situations and helps simplify and speed up the decision-making process. AHP is a method of decomposing complex unstructured situations into simpler components, thus producing a hierarchical system problem. This method was developed by Saaty in the mid-1970s and has been broadly considered and refined since then [10]. AHP is run in three stages: 1) Structuring the matrix of judgment between each pair of criteria; and 2) Making pair wise comparisons between the criteria and giving priority to each. 3) After the paired-comparison matrix is generated, it is normalized to match all criteria in the same unit.

To do so, 4 expert analysts were interviewed in order to score and weight the risk factors, and the average weight were given to the AHP for analysis.

All analysts were asked to answered two questions: What is the score of each risk (according to its impact and level of probability)? Which criterion is more important and stronger (pair wise comparison)? The aggregated results of risk [assessment] matrix and AHP are given in the following tables. Table 1 reveals the score for each risk independently, and Table 2 illustrates the weight of each risk based on pair wise comparison.

In terms of safety, the results showed that the size of the bus stop has the highest importance. Regarding health, the results indicated that the highest weight pertains to the density of population and prevalence of corona virus. Lastly, concerning traffic, the highest weight was related to traffic volume.

It was observed (from Table 1) that the density of population is higher in central and south zones; traffic volume is greater in the center and south due to density of businesses and government departments and agencies; and the location of bus stops is a common problem in most parts of Tehran.

Problem Description and Mathematical Model

The SBRP investigated in this paper is a simplification of the well-known vehicle routing problem (VRP), where a set of schools, one type of students, and a set of garages and identical buses-each of them with a fixed capacity-are incorporated. Each student is assigned to a bus stop, according to his/her maximum allowable distance. Afterward, the bus starts from the garage to pick up the students already assigned to it, and then goes to its assigned schools. Each student has a relevant school to which he/she is assigned. Since the problem considers both student allocation and routing problem while respecting the mixed load effect, it is possible to carry the students of different schools in the same bus. In order to make the model more practical, the problem is restricted to three different zones (north, center, and south), and risk analysis is performed accordingly. Subsequently, the zone with a high-risk score is customized in the model. The purpose is to optimize the total travel time. The most important assumptions of our problem are as follows:

1) Each bus starts from the garage to pick up the students from the stop, and drops them off at their respective schools;

2) Students must be allocated to their respective stops while maximum walking distance is satisfied;

3) Safety (size of bus stop), healthy factors (prevalence of corona virus and population density), and traffic concerns (traffic volume) are formulated in the model. In other words, the proposed problem features two concerns on the arc (traffic volume and population density) and health and safety issues on the nodes (prevalence of corona virus and size of bus stop), corresponding to each zone;

4) An upper time window is considered for each school, such that a bus should arrive at a given school and drop-off the student before the school’s latest time window;

5) A lower time window is considered for each stop, such that a bus should start its service to pick up the student after the stop’s earliest time; and

6) The maximum number of students allocated to each stop must not exceed the limit (Table 3).

The mathematical programming formulation of the school bus routing problem is as follows:

$$Min{\displaystyle \sum _{i\in N}^{}{\displaystyle \sum _{j\in N}^{}{\displaystyle \sum _{k\in K}^{}{t}_{ij}{X}_{ijk}^{}}}}\text{ (1)}$$

S.t.

${\displaystyle \sum _{j\in N}^{}{X}_{jik}}={\displaystyle \sum _{j\in (P^{+}\cup P^{-})}^{}{X}_{ijk}}=y_{ik}$                                    $\forall i\in P^{+},k\in K\text{ (2)}$

$${\displaystyle \sum _{j\in (P^{+}\cup P^{-})}^{}{X}_{jik}}={\displaystyle \sum _{j\in N}^{}{X}_{ijk}}=h_{ik}\forall i\in P^{-},k\in K\text{ (3)}$$

$${\displaystyle \sum _{i\in G}{\displaystyle \sum _{j\in P^{+}}^{}{X}_{ijk}}}\le 1\forall k\in K\text{ (4)}$$

$${\displaystyle \sum _{j\in P^{-}}{\displaystyle \sum _{i\in G}^{}{X}_{jik}}}\le 1\forall k\in K\text{ (5)}$$

$${\displaystyle \sum _{i\in G}^{}{\displaystyle \sum _{j\in G}^{}{X}_{ijk}}}=0\forall k\in K\text{ (6)}$$

$${\displaystyle \sum _{k\in K}^{}{y}_{ik}}\le 1\forall i\in P^{+}\text{ (7)}$$

${\displaystyle \sum _{k\in K}^{}{Z}_{il}^k}\le s_{il}\forall l\in S,j\in P^{-}\text{ (8)}$

$Z_{il}^k\le y_{ik}\forall l\in O_i,i\in P^{+},k\in K\text{ (9)}$

$y_{ik}\le {\displaystyle \sum _{l\in S}{Z}_{il}^k}\forall i\in P^{+},k\in K\text{ (10)}$

${\displaystyle \sum _k{D}_{jl}^k}\le q_{jl}\forall l\in S,j\in P^{-}\text{ (11)}$

$D_{jl}^k\le h_j^k\forall l\in W_j,j\in P^{-},k\in K\text{ (12)}$

$$h_{jk}\le {\displaystyle \sum _{l\in S}{D}_{jl}^k}\forall j\in P^{-},k\in K\text{ (13)}$$

$${\displaystyle \sum _{i\in P^{+}}^{}{Z}_{il}^k}={\displaystyle \sum _{j\in P^{-}}^{}{D}_{jl}^k}\forall l\in S,k\in K\text{ (14)}$$

$${\displaystyle \sum _{i\in P^{+}}^{}{\displaystyle \sum _{k\in K}^{}{Z}_{il}^k}}=1\forall l\in S\text{ (15)}$$

${\displaystyle \sum _{l\in S}^{}{Z}_{il}^k}\le ms\forall i\in P^{+},k\in K\text{ (16)}$

$$L_{ik}^{}=0\forall i\in G,k\in k\text{ (17-a)}$$

$$L_{ik}^{}+{\displaystyle \sum _{l\in S}^{}{Z}_{jl}^k}{H}_i\le L_{jk}^{}+bigM(1-X_{ijk})\forall i\in P,j\in P^{+},k\in k\text{ (17-b)}$$

$$L_{ik}^{}-{\displaystyle \sum _{l\in S}^{}{D}_{jl}^k}\le L_{jk}^{}+bigM(1-X_{ijk})\forall i\in P,j\in P^{-},k\in k\text{ (17-c)}$$

$${\displaystyle \sum _{l\in S}^{}{Z}_{il}^k{H}_i\le }{L}_{ik}^{}\le C\forall i\in P^{+},k\in k\text{ (17-d)}$$

$$T_{ik}^{}+ap.{\displaystyle \sum _{l\in S}^{}{Z}_{il}^k}+ad.{\displaystyle \sum _{l\in S}^{}{D}_{il}^k}+t_{ij}^{}\le T_{jk}^{}+bigM(1-X_{ijk})\forall i\in P,j\in P,k\in K,i\ne j\text{ (18)}$$

$$T_{ik}+t_{ij}^{}\le T_{jk}^{}+bigM(1-X_{ijk})\forall i\in G,j\in P^{+},k\in K\text{ (19)}$$

$$T_{ik}^{}\le T_{jk}^{}+bigM(1-Z_{il}^k)\forall i\in P^{+},j\in P^{-},l\in S,k\in K\text{ (20)}$$

$$T_{ik}\ge a_i-(1-y_{ik})bigM\forall i\in P^{+},k\in K\text{ (21)}$$

$$T_{ik}\le b_i+(1-h_{ik})bigM\forall i\in P^{-},k\in K\text{ (22)}$$

$${\displaystyle \sum _{i\in P^{-}}{\displaystyle \sum _{k\in K}^{}{X}_{ijk}}}\le Pg\forall j\in G\text{ (23)}$$

$$R_{ik}=0\forall i\in G,k\in K\text{ (24)}$$

$$R_{ik}^{}+A{r}_{ij}{C}_{ij}{X}_{ijk}\le R_{jk}^{}+bigM(1-X_{ijk})\forall i\in G,j\in P^{+},k\in K,i\ne j\text{ (25)}$$

$$R_{ik}^{}+A{r}_{ij}{C}_{ij}{X}_{ijk}\le R_{jk}^{}+bigM(1-X_{ijk})\forall i\in P,j\in P,k\in K,i\ne j\text{ (26)}$$

$$R_{ik}\le T_r\forall i\in P,k\in K\text{ (27)}$$

$y_{ik}\in \{0,1\}\forall i\in P^{+},k\in K\text{ (28)}$

$X_{ijk}\in \{0,1\}\forall i,j\in N,i\ne j,k\in K\text{ (29)}$

$Z_{il}^k\in \{0,1\}\forall i\in P^{+},l\in S,k\in K\text{ (30)}$

$D_{jl}^k\in \{0,1\}\forall j\in P^{-},l\in S,k\in K\text{ (31)}$

$$h_{ik}\in \{0,1\}\forall i\in P^{-},k\in K\text{ (32)}$$

$$r_{lk}\in \{0,1\}\forall l\in S,k\in K\text{ (33)}$$

The objective function (1) minimizes the total travel time traversed by all buses. Constraints (2) require that a bus entering the stop node should leave it as well. In constraints (3), the same constraints for school node are exposed (3). Constraints (4) ensure that a bus cannot start more than once from its home location (here, garage).Similarly, constraints (5) specify that it is not allowed for a bus to arrive at its final location (garage) more than once, which implies that some buses could remain unused. Constraints (6) specify that direct transfer from garage to garage is not possible. Constraints (7) enforce that each stop is visited no more than once. Constraints (8) enforce that each student is taken from the stop to which he/she walks. Constraints (9) represent that picking up a student from a non-visited stop by bus k is not possible. Constraints (10) guarantee that stops are not visited unnecessarily. Constraints (11) ensure that each student is delivered to his/her respective school. Constraints (12) guarantee that whenever a student is assigned to a bus, the school associated with this student is also visited by the same bus. Constraints (13) guarantee that schools are not visited unnecessarily. Constraints (14) impose that the number of pickup and delivery students in each route is equal. Constraints (15) state that each student should be picked up exactly once. Constraints (16) ensure that number of students allocated to each allowable stop must not be more than m _s . The next four sets of constraints (17-a), (17-b), (17-c) and (17-d) are load constraints. Thus, constraints (17-b) state that when a node i is followed by a pickup node , the number of students after visiting node j is greater than or equal to the summation of the number of students after servicing node i and the number of students picked up in node j . Similar to constraints (17-b), inequality (17-c) proves that when node i is followed by delivery node , the number of students after visiting node j is greater than or equal to the number of students after visiting node i minus the number of students delivered to node j . In practice, constraints (17-b) and (17-c) determine the load on a bus only after it leaves each node on its route. Constraints (17-d) specify the capacity of buses. Constraints (18)-(23) are time-related constraints. The arrival time of each bus to a node in p is calculated according to constraints (18). Constraints (19) are similar to constraints (18), but they are designed for the routes from garage to stop. Constraints (20) ensure that students are picked up by a bus before they are delivered. Constraints (21) and (22) indicate the time window for stops and schools, respectively. Constraints (23) restrict the number of available parking places in each garage. Constraints (24)-(27) ensure that the global risk is at most equal to the risk threshold Tr with big M (a sufficiently large number). Finally, variables and their types are presented in constraints (28)-(33).

Solution Methodology

The SBRP has already been proven to be an NP-hard problem; hence, for large instances, it is required to adopt a metaheuristic approach. It is necessary to explain how to design a metaheuristic that is capable of considering both diversification and intensification strategies and how to investigate this behavior through the solution space. In practice, to make an algorithm perform more efficiently, an appropriate trade-off must be achieved between intensification and diversification mechanisms. In this regard, there are some concerns pertaining to neighborhood structure and exploration mechanisms that need to be properly addressed. Undoubtedly, choosing the accurate order of operators will contribute positively to reducing computing time on the way to reach the optimal solution. Thus, it is crucial to choose the right order of neighborhoods in the search space. More importantly, we need to make a reasonable compromise between quality of solution and algorithm computing-time. Concentrating on neighbors which can ensure high quality may be desirable, but it results in more computing time. Conversely, focusing on smaller neighbors can hinder the perturbation effect in the solution space-hence diminishing the effect of exploration. To date, the adaptive mechanism of neighborhood selection is the more well-known approach in a broad range of routing problems. However, the recent work of Sörensen, et al. [11] revealed the extent to which the ALNS adaptive layer could help improve heuristic performance. Unsurprisingly, they showed that the benefit of using the adaptive selection mechanism is not considerable. They recommended that when researchers intend to use the ALNS, a special course of actions needs to be undertaken instead of blindly copying the algorithm from other studies. This shows that while designing the algorithm, it is imperative to find elements that boost its performance.

In this regard, in order to achieve a good performance, we need to have a close look at the following issues and implement them, if possible.

1) Rewarding the metaheuristic based on both its performance and executing time;

2) Paying careful attention to the difference made in the objective function while executing the operator;

3) Allowing the selection of operators that exhibit a poor performance, even with a low probability;

4) Weighing the neighborhood-based degree of intensification and diversification in the solution space;

5) Eliminating some useless and worse neighborhoods; and

6) Considering the time interval between the execution of each operator in the solution space.

In many studies, operators are weighted based on their performance in the past iterations. According to Sörensen’s findings, the adaptive mechanism alone cannot guarantee to find good solutions, and additional features need to be taken into account. Therefore, the most important task is to provide a well-defined structure for selecting both the order of local search operators and different approaches of the exploration.

It is important to make a metaheuristic analysis of the diversification mechanism in the solution space, type of neighborhoods (well-known or more specified), and role of each neighborhood in the performance (whether time or cost) of the algorithm. In this regard, we present different metaheuristic configurations as follows. For diversification, we introduce two kinds of structures: A multi-start structure and a perturb-and-improve structure, abbreviated as m- and p-, respectively. Concerning the selection of operator, we consider two kinds of neighborhood structure: Traditional order and systematic order (selection based on adaptive mechanism). Totally, we suggest four kinds of metaheuristics: Multi-start heuristic with a fixed order of neighborhoods (m-VND), multi-start heuristic with a systematic order of neighborhoods (m-HA), perturb heuristic with a fixed order of neighborhoods (p-VND), and finally perturb heuristic with a systematic order of neighborhoods (p-HA).

In both p-VND and p-HA metaheuristics, the initial solution is first built (construction phase). Afterward, the initial solution undergoes the improvement phase (second phase) until local optima is found. Finally, the diversification heuristic is applied to escape the current solution (the stuck solution) from the local optima. The diversification stage comes to play when the solution reaches out to local optima. To put it simply, this structure uses the construction phase once and repeats the improvement phase from the perturbed solution for a number of iterations. In contrast, in both m-VND and m-HA metaheuristics, the initial solution is built in the first phase and is transferred to the second phase for further improvement. In case local optima are discovered, the solution restarts again from the initial solution to reach more diversification.

The solution repeats both the construction phase and the improvement phase until the maximum number of iterations is determined.

The main difference between these VND and HA configuration is in the way each selects neighborhoods: The HA metaheuristics arrange the neighborhoods based on their performance, while the VND metaheuristics are constructed based on a fixed order, from small to large neighborhoods.

According to VND metaheuristics, all operators are examined based on the first-improvement strategy, which accepts all feasible movements that contribute to improving the current solution. The neighborhoods are ordered and explored from small to larger (more complex) ones. In case improvement occurs, the solution starts the intensification from the first operator. This procedure is terminated when the local search cannot further be improved by any of the operators. In the HA metaheuristics (adaptive mechanism), the operators are selected based on their performance. Thus, after executing the move, whether the solution is improved or worsened, the next operator needs to be selected via the roulette wheel mechanism.

To better analyze the exploration process, we introduced two kinds of diversification configurations, called multi-restart and perturbation mechanisms.

The structures of multi- and perturbation heuristics combined with the adaptive layer are given in Algorithm 1 and Algorithm 2, respectively. The allocation of students to bus stops is carried out in the construction phase, and the improvement phase (through the student allocation sub-problem), which is performed to check the feasibility of the solution. The mechanism for assigning students to allowable bus stops is comprehensively detailed in the study by Fallah, et al. [12].

Construction Phase

The construction phase is executed sequentially in three main steps: 1) Allocating students to potential bus stops (student allocation step); 2) Assigning the allowable bus stop to the available garages, preferably the closest garage (stop allocation step); and 3) Determining the route between the potential bus stops. In the first stage, each student is assigned to the nearest bus stop such that maximum walking distance and safety are observed. The safety issue requires that the number of allocated students to each allowable bus stop should not exceed the predefined value. The mechanism of the assigning students to allowable bus stops is fully described in the study by Fallah, et al. [12]. In the second step, the stops are assigned to the closest garages as long as the total number of allocated stops reaches the threshold (For each instance, the value of threshold is obtained through the total number of available stops divided by the number of available garages). For each garage, if the number of assigned stops exceeds this threshold, the stop is accordingly assigned to the next closest garage. This process continues until all eligible stops (i.e., those to which students are already assigned to) find the irrespective garage. In the third step, routing is conducted through a greedy randomized adaptive search procedure (GRASP), because it helps overcome the myopic behavior of the greedy heuristic as it creates a balance between intensification and diversification. In the routing procedure, two issues are emphasized for inserting non-visiting nodes: Generating new routes from garages and, correspondingly, adding non-visited stops to the current route. For each new route from the garage, all non-visited-stops are sorted in increasing order of their criterion function (explained further) and added in the list U. Then, the restricted candidate list (so called RCL), which contains the α first nodes from the list U, is created. After that, one node is selected randomly from the RCL and inserted as the first node of the current route. The selected node is removed from the list U and saved as the current node. New non-visited stops in U are selected and added to the current route on the basis of the procedure explained below.

1) An eligible list is created to sort all non-visited stops in increasing order of the criterion function value. It is worth mentioning that the eligible list (Le) is applied to nodes that respect predefined constraints.

2) In this step, the list RCL-e is constructed that comprises the α first nodes. If the considered RCL-e is empty, new routes are generated. Otherwise, one stop is selected from the RCL-e and inserted in the current route; it is then saved as the current node and removed from the list U. In case a new route needs to be generated, the bus (i.e., route) returns to its respective school(s) (dropping off the students) and goes back to the closest garage. That means it is not mandatory for the bus to go back to the garage from which its operation started.

In this study, two different kinds of criterion functions are considered. The first one relies on the mechanism of selecting the next non-visited stop on the basis of the closest distance. The second criterion function is introduced as cf = (n + m)/d, where n demonstrates the number of eligible student(s) and also m represents the number of students in the candidate stop that have a common school with students who have already been picked up by the bus; finally, d shows the distance between the last visited node (garage or stop) and the candidate stop. Each criterion function is calculated for each non-visited node in the list U. This procedure is completed and updated when a new request is received. The construction phase serves two objectives: 1) Concentrating on the selection of the next non-visited node based on the problem characteristics; and 2) Generating two different kinds of initial solutions. The important point is that in visiting each node, the health risk needs to be investigated more carefully. In order to make the problem closer to reality, the health risk is identified and formulated based on the recent statistical data pertaining to the prevalence of the corona virus in any zone. If the visited node is positioned in an area with a high prevalence rate, it belongs to the highly risky or high-danger zone. In case the incidence rate is slightly lower in an area, it is considered as part of the orange or danger zone. Finally, an area where the situation is relatively normal and rather safe is considered within the yellow zone. Region classification is conducted based on the results of risk analysis. Accordingly, the north of Tehran is marked as the yellow zone, the center corresponds to the orange zone, and the south is associated with the red zone.

It is worth mentioning that when the bus visits the nodes, depending the location of the stop, the penalty demand is proportionately applied. More precisely, if the bus stop is in the yellow zone, the number of each student is considered as one demand; when it passes the orange zone, each student is considered as 1.1 demands; and if the bus crosses the red zone, each student is considered as 1.2 demands. This simple strategy requires that the bus is obliged to pick up a smaller number of students in areas with a high prevalence rate, thus paying special attention to social distancing. The important item is the size of α: If it is small, the construction heuristic concentrates on the greedy mechanism, but if it is large (equal to N), the solution adopts a random behavior.

Improvement Phase

Once an initial solution is constructed, it is delivered to the next stage for further improvement. As mentioned already, two improvement strategies are suggested: 1) Selecting the neighborhoods using a fixed traditional order; 2) Selecting the neighborhoods based on their performance. The performance mechanism of each neighborhood relies on the execution time of the operators and their role in changing the objective function. In a nutshell, the first strategy justifies the VND heuristic as a simple but effective algorithm, and the second strategy helps make an appropriate trade-off between intensification and diversification mechanisms.

Variable Neighborhood Descent (VND)

The improvement phase is based on the variable neighborhood descent (VND), a species of the variable neighborhood search (VNS) metaheuristic method [13]. The VND heuristic needs a hierarchical order of neighborhoods. This ordering is usually achieved through putting simpler, less complex, and less perturbative neighborhoods at the beginning of the list. To escape from the local optima, one must choose larger neighborhoods when no further improvement can be achieved from the small heuristics.

Small neighborhoods intensify fewer solutions and, therefore, can be searched in less time than large neighborhoods. Larger neighborhoods are only recalled when all smaller neighborhoods become useless, i.e., when the current solution is a local optimum in regard to all smaller neighborhoods. As soon as any local search improves the current solution, the move is executed and VND intensifies the search from the first neighborhood in the list. The VND stops its operation when the solution reaches the local optimum for all considered neighborhoods.

In order to find the right order of neighborhoods, different combinations of neighborhoods are tested, and the most promising order is incorporated in the VND algorithm.

Prior to implementing any intra or inter local search operator, both the cost of the solution and its feasibility in regard to the predefined constraints must be checked. If the local search operator finds better solutions and the considered constraints are satisfied, the corresponding move is executed; otherwise, it is discarded. The local search is terminated when the solution is stuck in a local optimum and cannot be further improved by executing any of the local search operators. As mentioned already, both m-VND and p-VND metaheuristics utilize the VND in the improvement phase.

Modified Adaptive Neighborhood Selection

In the second improvement strategy, following the well-known method of ALNS (adaptive large neighborhood search, proposed by Ropke and Pisinger [14]) though with a slight modification, we select the neighborhoods based on their performance. In this strategy, the invocation of the neighborhoods happens in the order of their effectiveness with respect to the problem at hand.

More practice, the search process is divided into a number of segments, each consisting of a number of consecutive iterations n . The roulette-wheel mechanism is used to select the next neighborhood h _i for the operation in each iteration n as follows:

$\frac{w(h_i)}{{\displaystyle \sum _{i\in I}^{}w(h_i)}}$                                                                                                 (34)

Where w ( h _i ) demonstrates the weight of each employed neighborhood. This probability depends on the weight of each operator. At the beginning, each neighborhood has the same weight, but at the end of each segment (after a number of consecutive iterations), all heuristics’ weights are updated according to the following formula:

$w_{i,j+1}=(1-\gamma )w_{i,j}+\gamma \frac{\mu (h_i)}{{\varphi }_{ij}}$                                                      (35)

It means when a segment ends, the new weights are calculated based on the accumulated score of each heuristic in that segment.

In the above formula, w _i,j signifies the weight of neighborhood i in segment j, $\mu (h_i)$ shows the accumulated score of neighborhood i in the last segment, and ${\varphi }_{ij}$ refers to the number of times neighborhood i is repeated in the last segment. The value of 0 < g < 1 is an important factor that helps control the behaviors of the adaptive mechanism in the proposed algorithm. When g is set close to 0, the previous weight of the heuristic is considered; however, if g is set to 1, the weight of the heuristic is considered only based on new accumulated scores.

At each iteration, the score of neighborhoods used is updated as follows: σ ₁ when the applied heuristic results in a global best solution thus far; σ ₂ when the selected heuristic provides a smaller solution than the current solution; and σ ₃ in case the applied heuristic results in a worse solution than the current solution.

The modification introduced in this study is related to the mechanism of neighborhood scoring. More precisely, instead of assigning a fixed score to the neighborhoods that create better or worse solutions, the score of each neighborhood is calculated based its performance and role in the intensification and diversification mechanism. This method dynamically controls the selection of operators based on a combination of three partial functions, which are combined to determine the score of heuristic i in each iteration n as follows:

$$\mu (h_i)=0.5*{\varphi }_n(f_{{}_1}(h_i))+0.5*{\varphi }_n(f_{{}_2}(h_i,h_k)+{\delta }_n{f}_{{}_3}(h_i)\text{ (36)}$$

The first criterion is expressed by the following formula:

$f_{{}_1}(h_i)=\frac{I_{{}_n}(h_i)}{T_n(h_i)}$                                                                                                         (37)

Where $I_{{}_n}(h_i)$ represents the change in the objective function when the heuristic is applied, and $T_n(h_i)$ is the time it takes the heuristic to explore.

The second criterion investigates the dependency among the operators with respect to their performance, and it is expressed by the following formula:

$f_{{}_2}(h_i)=\frac{I_{{}_n}(h_i,h_k)}{T_n(h_i,h_k)}$                                                    (38)

Where $I_{{}_n}(h_i,h_k)$ represents the change in the objective function, and $T_n(h_i,h_k)$ represents the time it takes the heuristic it to be recalled after the heuristic k.

The third criterion, given below, expresses the time elapsed since the heuristic was last executed. It shows how long the heuristic has been inactive and examines the possibility of using it.

$f_{{}_3}(h_i)={\tau }_{{}_n}(h_i)$                                                                          (39)

The measures f₁ and f₂ are used to control intensification, and the measure f₃ is used to control diversification. ${\varphi }_{{}_n}$ and δ _n parameters in this model are used to weigh f₁, f₂ and f₃. Specifically, ${\varphi }_{{}_n}$ is the intensification parameter that controls and weighs f₁, and f₂ and δ _n is used to weigh f₃. In each iteration, if there is an improvement in the objective function, the value of ${\varphi }_{{}_n}$ is increased and the value of δ _n is reduced. Conversely, the value of ${\varphi }_{{}_n}$ decreases and the value of δ _n increases when the objective function shows a worse performance. This concept is expressed as follows:

${\varphi }_{{}_n}(h_j)=\left\{\begin{array}{l}0.99ifimprove\\ \max (({\varphi }_{{}_{n-1}}-0.01),0.01)otherwise\end{array}\right\}$                                                 (40)

${\gamma }_{{}_n}(h_j)=1-{\varphi }_{{}_n}(h_j)$                                                                                    (41)

Compared to the traditional VND, our adaptive model is novel in many ways: 1) The weight assigned to each neighborhood is consistent with its performance, time of operation, as well as dependency among operators; 2) Our method is more compatible with the size of the problem; 3) It can create a balance between diversification and intensification in the search space; and4) The bad move, too, has a chance to be selected, though with a lower probability.

Type of Neighborhoods

Two types of neighborhoods are incorporated in this study which are addressed in the special and common to the SBRP. The common neighborhood resembles the well-known VRP, whereas special cases are compatible with problem characteristics such as single load, mixed load, and health concepts. Since two different strategies laid out to select the neighborhoods, the mechanism of exploration is different in each case. In the adaptive layer, neighborhoods are selected based on their performance, while the VND approach demands that the neighborhoods be recalled in order of complexity, from small to large (more complex) ones.

To obtain the right combination of neighborhoods through the VND approach, we performed a pilot study (see Section calibration).

Remove-Insert Intra Routes Operator

This operator removes one stop from the candidate route and inserts it in another location within the route without checking the feasibility constraint.

2-Opt Intra Route Operator

This operator attempts to remove and reinsert two edges between stops along the same route.

Remove-Insert Operator (Considering Health Feature)

This operator attempts to remove one stop and reinsert it in the same zone (whether on the same route or not). It performs this move in each zone (red, orange, and yellow) independently, thus (1) helping to avoid further spread of corona virus between the zones and (2) increasing intensification in the solution space.

Similar School-Remove-Insert

The operator separates bus stops and adds them on the routes where the schools related to the students of those stops are already assigned. This move will prevent additional insertion of schools on the route, and this local search helps local intensification.

Convert Single to Mixed Load

This operator serves to reduce the number of single load routes as well as remove the stops from single-load routes and insert them on mixed-load routes. This procedure is in line with the concept of mixed load effect.

Convert Mixed Load to Single Load

This operator removes the stop from the mixed-load route and insert it in another route in order to remove the stop whose students travel along-distance relative to existing school(s).

Replace

In order to minimize total traveling time, this operator chooses a non-visited stop from the list and then replaces it with a stop that is already included on the route. The list of non-visited stops is sorted in a descending order of the number of students who can reach that stop. Before applying any move, the cost check as well as a student allocation sub problem must be executed.

Remove

This operator is executed by removing a stop from the current route in order to decrease the cost of the route while maintaining solution feasibility. Like the above move, before applying any move, a student allocation sub problem must be solved. If it contributes to achieving a feasible solution, the move is performed. Due to the triangular inequality characteristic, it is not essential to perform cost check since a stop removal will always produce a reduction in the travel cost.

Diversification Stage

The diversification strategies are used to generate new initial solutions and, consequently, escape from local optima and explore different parts of the search space. This contributes to creating a new promising starting point as an input for the local search block. In fact, the diversification mechanism will navigate the algorithm to search for unexplored parts of the solution. In order to discover the effect of diversification, we must have the mechanism investigate the degree of perturbation. This shows that, the size of perturbation is crucial since it can result in supercharge improvement or worsen the current solution. Previous evidence suggests that lower deterioration increases the risk of getting trapped in local optima. On the other hand, large degrees of deterioration can mislead the path leading toward the global optima. Therefore, the design of diversification heuristics must be compatible with the size of perturbation and status of exploration. In this regard, two different diversification structures are proposed: 1) Perturbing part of the solution space (randomly or constructively); 2) Generating new initial solutions compatible with the multi-start procedure.

The multi-start metaheuristic attempts to escape from local optima via repeating a number of iterations both in the construction and the improvement phases. In contrast , the perturb configuration implements the construction phase only once and repeats the local search and perturbed heuristic for consecutive iterations.

During the perturb configuration, the destroy and repair heuristics are employed. First the best solution found is ruined partially by the destroy operator. Afterward, the repair heuristic helps reconstruct the new current solution, hence preparing it for the local search block. In the destruction phase, we introduce two kinds of heuristics: 1) Random-destroy heuristic, which starts to destroy random part of the solution space; 2) Constructed-destroy operator, which attempts to destroy the part of the solution with the highest removal gain. When a random part of the solution is selected, there is no prior knowledge about which part of the solution needs to be destroyed. In contrast, when the constructed-destroy heuristic comes to play, that part of the solution space is destroyed that brings the lowest gain (depending to previous performance). In the following, we only rely on the explanation of perturb configuration, which encompasses the destroy and repair operators.

Destroy Phase

Random destroy heuristic

In this method, the random part of solution space (i.e., random bus stops) are selected, removed and added to the list U. The size of perturbation is controlled by a value of p * q , where p denotes the number of routes constructed in the improvement phase and q is the size of perturbation, which is controlled between q _min and q _max . If the perturbation mechanism leads the search toward reaching the global best solution in the next improvement phase, the value of q is set in q min so that the intensification strategy is preserved. However, if no global best solution appears, the value of q is increased by 10% in order to maintain diversification. This process continues until the value of q reaches q _max . In this case, the value of q is set to q _min again.

Constructed destroy heuristic

In this stage, the heuristic attempts initially to select the stops based on the ratio of , where the value of d _i represents the amount of demand in each stop, while the value of s _i is represented by the following formula. This ratio is sorted for all bus stops in descending order. Afterward, k stops from the top of the list are selected and added to list U in order to select the bus stop with the lowest demand and the highest removal gain. Meanwhile, the priority is to remove the stops and add them in a location that will reduce the cost. The value of k is controlled between  and , where n is the number of stops and  and  are minimum and maximum percentage of stops to be removed in the diversification stage.

Each time diversification is applied, one of the aforementioned destroy operators is randomly nominated.

Repair phase

The repair operation is simple and effective. Accordingly, all removed stops in list U are added to the solution via the GRASP procedure, which helps complete the current solution. This procedure is followed without violating the considered constraints. In this stage, all students of a given stop in list U remain fixed, and there is no necessity to perform student allocation. The newly generated solution will be the input of the local search heuristic to be considered for further improvement.

Instance generation

The mechanism of generating some features is similar to what was described in Fallah, et al. [12] (generating stops, school, students, and garage); however, here we only discuss the issues related to risk consideration. In this regard, our generated instance is distributed in three zones: North, center and south, which are represented by a, b, and c areas, respectively. In order to incorporate the risk concept into the current problem, we take account of characteristics such as density of population and traffic volume (the risk with a high score). To put it more precisely, the high-density population risk and high traffic volume risk are limited to the center and south (zones B and C), which highlights the importance of considering the effect of risk analysis in the problem under study. For these areas, the traffic volume and population density risk factors are randomly generated in an interval of (1.1 to 1.40), respectively. For other areas (with low and medium levels of risks), both traffic volume and population density risks are set to 0in order to mitigate the complexity of the problem. In order to visualize the prevalence distribution of COVID-19, three health coefficients are considered for the north (= 1), center (= 1.1), and south (= 1.2) zones. This coefficient is multiplied by the number of students to be picked up or dropped off. Accordingly, when the bus picks up a student from the north zone, he/she is equivalent to 1, but when a student is picked up in the south zone, he/she is equivalent to 1.5. Therefore, for a bus with a capacity of 25 passengers, we can only pick up a maximum of (25/1.2 = 20) students in the southern zone, hence adhering to social distancing rule in the bus. In other words, for the zone with a high prevalence of COVID-19, the bus picks up a limited number of students lower than its default capacity.

Computational Results

Our computational experiments were carried out in three different stages outlined below:

In the first stage, the heuristic parameters were tuned (see Table 4) by conducting a full factorial experiment on a subset of instances.

In the second stage, the performance of four proposed algorithms on all instances (small, medium, and large) was compared.

In the third stage, some risk analyses were undertaken in order to show the effect of risk constraints on the objective function.

Calibration of Metaheuristic

The metaheuristics presented earlier have parameters that need to be adjusted and calibrated so as to make a reasonable compromise between the quality of the solution and computing time. Setting the controllable parameters of metaheuristic algorithms is of utmost importance while designing these algorithms. It helps better understand the behavior of each metaheuristic. In order to tune the controllable parameters of multi-start and perturb configurations, the parameters were evaluated under different values and the best values were chosen. The main parameters are described in Table 4. A full factorial experiment (combining all the parameter settings shown in Table 4) was performed to solve 10 instances (4 instances from set S, 4 instances from set M, and 2 instances from set L) by using the four suggested metaheuristics. For each instance and each specific parameter, 5 runs were executed. In some of the runs of the analysis, all neighborhoods in the improvement phase were deactivated for all metaheuristics. In that case, the algorithm behaves like a random restart. Two performance measures were investigated: the average solution cost and the average computational time. It is important to note that the maximum number of iterations is not included in the list of analyzed parameters, since a large number of iterations will logically provide better solutions in longer computational time. Thus, in this stage (calibration stage), the number of iterations remains fixed in 350 for all calculations, and only the maximum number of non-improvement iterations is considered in our analysis. Table 5 presents the P values of the F-tests. Asterisk values indicate the significance of each parameter on both the objective function and the solution time of the algorithm. According to Table 5, important parameters of the algorithm that significantly affect both the solution quality and computational time include: The majority of local search operators, minimum and maximum percentage of routes and stops to be destroyed, maximum number of non-improvement iterations. Meanwhile, similar school removal and mixed-to-single-load- operators are the parameters that only have a significant effect on computing time.

The best values of the parameters are given in Table 6. Analyzing the considered neighborhoods reveals that, on average, the effects of similar school removal (for all metaheuristic configurations) and mixed load-to-single-load operators (for m-VND and p-HA metaheuristics) on the quality of solutions are less than other operators, although they generate a slight improvement in the results when activated. Consequently, these two operators are not included for further analysis.

It worth mentioning that we independently investigated the effect of different numbers of iterations on the quality of solution. Other parameters were fixed at their best value and were taken from Table 6.

As in the above, 12 instances were used for analysis and the solution was tested for different numbers (200, 250, 300, 350, 400, 450, 500, 550, and 600) of iterations. The analyses demonstrated that the best number of iterations which provides an ideal compromise between computing time and quality of solution is 450.

Generally, establishing the promising order of neighborhoods employed in a VND heuristic may have a substantial role in the quality of the solution. We tested different orders of neighborhoods on a subset of instances to understand which combination would work better. On average, the best order of neighborhoods incorporated in the VND heuristic is given as follows.

Metaheuristic Comparison

Having determined the optimal parameter settings for each solution approach, now we compare the proposed metaheuristics by testing them on small, medium, and large instances. To test and make a fair comparison between the four proposed metaheuristics in term of the quality of solution and computing time, we ran each metaheuristic 5 times over a fixed number of 450 iterations for each instance.

The experimental analysis was executed on 100 instances consisting of three subsets named Set S, Set M, and Set L. Set S covered 25 instances (5-10 stops), set M consisted of 30 instances (20-25 stops), and finally set L comprised 30 instances (40-50 stops). Additionally, we considered a maximum walking distance of 5, 10, 15, and 20 in our calculations. (For more details on the results of each experiment, see Appendix 1).

The aggregated results are shown in Table 7 and Table 8. Table 7a highlights the percentage gap between the best-found solution of each metaheuristic after 5 runs and the optimal solution (GAMS/CPLEX solver). Table 7b presents the percentage difference between the four metaheuristic configurations. Specifically, after all runs were executed, the metaheuristic with a better quality of solution (on average) was considered as the basis and the performance of other metaheuristics was compared against this one. In both Table 7a and Table 7b, the first column shows the type of metaheuristic, while the next three columns present the set of instances (from small to large).

Analyzing the results in Table 7a reveals that p-HA features the lowest percentage gap with the optimum solutions (exact method) in small instances. Also, m-VND ranked second in terms of this gap.

Comparing the four proposed metaheuristics (Table 7b), we realize that on average p-HA outperforms the other configurations. It seems that, on average, the combination of perturbation and adaptive mechanisms could lead to a better performance. Accordingly, p-HA exploits two mechanisms that are implemented simultaneously: Perturbation heuristic and selecting local search operator based on its performance. Hence, p-HA is considered as the basis algorithm against which comparisons are made with the other proposed metaheuristics. To better understand the performance of the four metaheuristics, we compared them in small, medium, and large instances. It was noted that m-VND, followed by m-HA, has a lower percentage gap than the basis metaheuristic; p-VND showed a higher deviation from the basis metaheuristic, suggesting this combination is not consistent with the considered instances. Therefore, this combination is of little interest. On the other hand, in terms of computing time, Table 8 shows that m-HA is the slowest metaheuristic and relying on its combination requires more computational time.

Risk Factor Effect

This supplementary analysis is undertaken to investigate the effect of each risk constraint on the objective function. Due to its better performance compared to the other three proposed metaheuristics, the p-HA configuration was chosen for this purpose. The experiment was performed on a subset of 12 instances. The main purpose of this experiment was to determine the extent to which each risk constraint affects the suggested objective function. More precisely, we sought to understand how much the value of the objective function would increase if all the risk constraints were active in the model, and how much it could diminish if the risk constraints were relaxed. The results revealed that in case of omitting the safety risk constraints (constraints 1), the value of the objective function decreases by up to 6.53%, and if all risk constraints (safety, healthy and traffic) are relaxed, the objective function decreases by 15.07%. This significant difference must be considered by the municipality or any organization in the private sector. The results also showed that the size of the problem has a positive effect on the risk constraints. Moreover, risk constraint relaxation significantly affects the objective function in all instances, particularly large ones (Figure 2).

Conclusions and Future Research

This paper dealt with the school bus routing problem by combining both bus stop selection and route generation while considering a mixed-load context and some risk constraints. In particular, depending on the walking distance from their home, students are assigned to potential bus stops. At the same time, these stops are inserted in the bus route to take the students to the school such that the travel time is minimized. This study also considered the safety, traffic and health risk constraints, including maximum number of students to be assigned to a bus stop, as well as traffic and health risk constraints (imposed on the arcs and nodes). To solve the generated small, medium, and large instances of the problem (totally 100 instances) in a reasonable computing time, we devised and suggested four different metaheuristics that are capable of generating different diversification strategies in the solution space. Overall, on average the p-HA metaheuristic outperformed other metaheuristics in regard to the quality of solution. As for the computing time, p-VND metaheuristic, followed by p-HA (though with a slight difference), showed a faster performance than others. It seems, therefore, that p-HA can be a more reliable metaheuristic. On the other hand, m-HA configuration cannot provide the solution in a reasonable time and, hence, requires more computing time. These metaheuristic configurations move in two directions, placing as much emphasis on the mechanism of neighborhood selection as on the status of diversification. It can be concluded that while designing the metaheuristic, we need to take account of the mechanism of diversification, kind of neighborhoods, and neighborhoods’ execution time. We conducted some risk analyses and identified three major risks. Clearly, these risks restriction will be more salient after student transportation situation returns to its normal situation. With the reopening of schools, the safety of school buses will be a major concern for many transport companies and students’ parents. In this regard, taking safety instructions seriously can help efficiently mitigate and control such concerns, thus reducing the probability of accidents. Regarding the safety risk, many accidents occur for students when they are waiting for the bus or are getting on or off the vehicle. However, following a few simple but effective measures will protect students against many injuries:

1.       Students should be asked to wait for their school service in a safe place, such as on the sidewalk and away from traffic.

2.       Ensure that students are not allowed to walk in the bus until it stops and its doors are opened.

3.       Remind students that they should never run behind the bus as they may not be seen while the driver is entering the stop.

4.       Teach students to choose a suitable chair, sit on it, and do not stand while the vehicle is moving, as a sudden brake might cause them to fall.

5.       Students must be asked not to speak loudly. Noisy conditions caused by students increase the threat of driver distraction and disturb his concentration on outside circumstances such as red lights and police warnings.

6.       Students must never talk to the driver as it may lead to an accident.

7.       A risk committee should be set up to analyze and evaluate the potential risks (in terms of probability and severity).

8.       Appropriate actions should be introduced to control and mitigate defined risks.

9.       A cost analysis should be performed for the proposed actions.

10.    The results of risk analysis and the measures taken should be shared with other related stakeholders. 

Concerning population density risk in some areas, it seems necessary that the municipality assign more buses to communicate students in a safer and more convenient fashion. In this regard, it can enter a joint collaboration with the private sector to handle this issue much better.

Regarding the healthy risk as a short-term action, the focus should be on the social distancing approach and explaining the health requirements for students while boarding the bus.

Future research can be aimed at including stochastic parameters in the model. In this paper, all input parameters of the case study (e.g., arrival time of bus, speed of bus) were considered to be deterministic. However, in real-world applications, there are different factors that make these parameters stochastic. For example, owing to an unpredicted situation, a bus may have a delay on the route. Moreover, the expected speed of the bus is not usually a fixed parameter due to traffic conditions. As a future work, these parameters can be considered within uncertain conditions. Another interesting idea is to consider an objective function to minimize the total student journey, traveling time from home to stop, waiting time in the stop, traveling time on the bus, as well as pick-up time at the stop and drop-off time at the school. Finally, it is more rewarding to investigate both hard and soft time window constraints at the same time.

---