Robot Routing Problem of Last-Mile Delivery in Indoor Environments

: With the development of robot technology, trials adopting robots for last-mile delivery are continuing, and the ﬁnal destination of last-mile delivery is further expanding into indoor environments. Unlike existing studies conducted for robot-based last-mile delivery in outdoor environments, two main issues must be solved to enable last-mile delivery in indoor environments using robots. First, it is necessary to reasonably and realistically estimate the robot travel time considering horizontal and vertical movement segments within a given building. Second, optimizing the robot routing problem based on the estimated robot travel time is necessary. In this paper, we proposed a new method to estimate the robot travel time considering robot movement characteristics and an elevator in a building. In addition, we developed a mathematical model of the robot routing problem and problem-speciﬁc heuristic based on a genetic algorithm to quickly solve the proposed mathematical model. It obtained the exact solutions when the problem size was small and near-optimal solutions in the medium-and large-sized problems (average optimality gap: 0.11% and 0.18%, respectively). Through extensive experiments assuming various building structures, it was determined that the proposed model and heuristic can quickly yield realistic solutions for indoor robot-based last-mile delivery


Introduction
E-commerce has increasingly proliferated, and the lockdown measures associated with the COVID-19 pandemic have accelerated this trend [1]. Following this trend, the demand for last-mile delivery, i.e., the final delivery stage of ordered items, has explosively grown. In addition, urbanization expansion has exerted additional pressure on developing practical, innovative, and sustainable city logistic solutions [2]. Last-mile delivery is challenging but crucial because this stage is related to the final contact with customers. Thus, the level of automation is much lower than that of upstream logistic activities, such as warehousing and transportation between terminals.
However, the ongoing COVID-19 pandemic has increased interest in new technology adoption, such as drones and robots, for last-mile delivery purposes since these devices can provide unmanned delivery services. Even though both drones and robots have received attention in regard to contactless delivery, their applicability in terms of delivery coverage differs. For example, drones are the most likely choice in rural areas, while robots are more likely to be chosen in urban areas [1].
This research focuses on the indoor environments inside the multistory buildings in dense urban areas where robot-based delivery seems more suitable. Unlike existing studies conducted for robot-based last-mile delivery in outdoor environments, two main issues must be solved to enable last-mile delivery in indoor environments using robots.
First, estimating the robot's travel time reasonably and realistically is necessary to capture the vertical and horizontal movement within a given building. Extant literature thinks a robot's travel time is the same as a truck's, just distance over speed. However, it is insufficient to simply obtain the robot travel time as in the previous studies when planning indoor robot delivery. Regarding the horizontal movement of robots on any floor, robot

•
A new travel time calculator is provided to estimate robots' reasonable and realistic travel time considering their vertical/horizontal movement in indoor environments.

•
An efficient heuristic was developed to generate solutions quickly, even for large-sized multistory buildings.
The remainder of this paper is organized as follows. Section 2 introduces the related literature. Section 3 provides problem description and mathematical models. Section 4 introduces a heuristic-based solution to manage large-size problems that can frequently occur in practical cases. Section 5 presents experimental results to verify the necessity and performance of the proposed approach. Finally, concluding remarks and future research directions are outlined in Section 6.

Literature Review
The interest in freight distribution with robots has surged over the past several years, producing numerous field tests and studies [4][5][6][7]. Since this research is related to last-mile delivery in indoor environments, we aimed to focus on service robots that can move and deliver indoors.
Continental illustrated a pack of robotic dogs exiting an autonomous vehicle and entering a building for final delivery [8]. A robotic dog is an autonomous, four-legged robot that can manage obstacles in indoor environments. The robotic dog can traverse challenging terrains such as curbs and stairs and can deliver a package in front of the door [9]. Additionally, there exists a bipedal robot, approximately the size and shape of a small adult human, that can navigate indoor environments with the help of light detection and ranging (LIDAR) technology and other sensors [10]. Even though both concepts have yet to be pilot tested with customers, they have indicated the possibility of employing robots to deliver ordered items in buildings.
Along with the development of relevant robot technology, academia has shown interest in last-mile delivery using robots for the past couple of years [11]. Most related research utilized the existing vehicle routing problem (VRP) to adopt robots in outdoor environments [12][13][14][15][16][17]. Last-mile delivery problems have been solved with various VRP formulations via mixed integer linear programming (MILP), and efficient algorithms have been developed to address larger problems due to the inherent complexity of the VRP. Of course, it might occur that the outdoor environment remains a more general solution for delivery services because technical difficulties are encountered indoors where the global positioning system (GPS) is inoperable [18].
However, owing to advances in simultaneous localization and mapping (SLAM) and the control of robots [19][20][21], research on VRP with robots in indoor environments has been initiated. The studies [22,23] attempted to optimize a task allocation problem (e.g., patrol and delivery services) and plan routes for indoor robots. They suggested a general framework and algorithm, referred to as the fixed destination multi-depot multiple traveling salesman problem. Their study demonstrated that the approach should differ between indoor and outdoor routing. Additionally, [24] introduced an operation design of a robot logistics system for the hotel industry. They assumed that the robot could serve dual jobs, such as delivering and retrieving dishes within a hotel. They calculated the minimum number of required robots, considering a hotel with a long hallway on each floor, and then built a mathematical model to minimize the number of unassigned tasks to robots within the unit execution cycle, e.g., 15 min. Table 1 summarizes the recent research related to robot-based delivery problems. Although the adoption of robots for last-mile delivery is an emerging field in both the industry and academia, almost all attention in academia remains focused on the outdoor routing problem. In other words, there is a lack of research focusing on the indoor routing problem, as indicated in Table 1. Although [24] considered vertical movement using an elevator in their travel time calculation procedure, they did not capture the possibility of intermediate stops or any characteristics of horizontal robot movement (e.g., acceleration, deceleration, and cornering).
Unlike outdoor truck routing problems, the indoor robot routing problem (IRRP) requires a more plausible and sophisticated travel time calculation. Applying the existing Appl. Sci. 2022, 12, 9111 4 of 20 method to IRRP will reduce the reliability of the derived results. To the best of our knowledge, this paper is the first attempt to tackle the IRRP while considering the realistic travel time of robots moving both horizontally and vertically during last-mile delivery. A general but powerful GA-based heuristic also increases the practicality of the research. Our research will help companies that want to try indoor robot-based delivery services.

Problem Description and Mathematical Models
This section describes the problems and introduces the mathematical models of the IRRP. The network representation and travel time calculation processes are explained first, and a mathematical model of the IRRP is then formulated.

Travel Time Calculator Considering Robot Movement
As mentioned in Section 1, one of the two essential issues to consider robots for delivery in indoor environments is the realistic estimation of the travel time. The travel time between two points within buildings requires a few seconds up to a few minutes; thus, if errors are accumulated, then the gap between the actual travel time and the estimated travel time can increase. Therefore, a different approach is needed when a vehicle moves indoors, in contrast to the simple estimation approach of the travel time based on the average speed of the robot. In fact, in the case of truck routing problems, assuming the average speed seems reasonable. Since trucks usually require several hours to move between outdoor locations, any delay occurring due to unexpected road conditions can be offset via acceleration before reaching the destination. However, in the case of indoor robot-based delivery, robots slowly move between indoor locations similar to a human operator walking. Additionally, robots might move much slower if corners are encountered. Moreover, using an elevator inside a building results in an extra delay in the journey. All travel times of all robots within buildings are counted in seconds rather than in minutes or hours. Thus, there is insufficient opportunity to compensate for any delay resulting from corners or elevators within buildings.
In this subsection, we divide the indoor movement process of the robot into segments along the horizontal and vertical directions, calculate the travel time for the horizontal (t h ij ) and vertical (t v ij ) movement segments separately, and then sum these travel times to determine the final travel time (t ij ) of the robot.
First, which factor imposes the most notable influence on the robot travel time required for indoor horizontal movement?
The presence of corners is the greatest obstacle to horizontal robot movement across the floor. In a typical building, corridors and corners are continuously connected, facilitating movement. Therefore, corners mainly affect the travel time because they limit the speed of robot movement. After initiation, the robot increases its speed during movement, and when a corner is encountered, the robot slows down to round the corner. Additionally, after rounding the corner, the robot again accelerates and moves through the hallway.
As expected, the more corners between any two points, the more the robot must accelerate and/or decelerate to round these corners, which leads to a delay in travel time. Of course, the lag time to round a corner can comprise only a few seconds, but the lag time can be extended if there are many corners between two points.
Before estimating the travel time required for the horizontal movement segment of the robot indoors, the following assumptions apply:

•
The indoor pathway of the robot consists of corridors and corners.

•
The robot remains stationary at the starting point before movement occurs.

•
The robot immediately accelerates after starting, and when it reaches an appropriate speed, the robot moves while maintaining this speed.

•
The robot decelerates ahead of its destination and stops upon reaching the destination.

•
The robot decelerates to the predetermined speed ahead of a corner and accelerates after rounding the corner.

•
The acceleration and deceleration ratios (the acceleration and deceleration rates per unit time) of the robot remain constant.
To explain how the travel time is estimated in this research, certain situations were assumed. Figure 1 shows the example of indoor robot movement from origin to destination in a multi-story building. First, let us assume the robot only moves on a single floor (n th floor) to explain the calculation procedure for the horizontal movement. There are two corners between the starting node i and the destination node, the elevator. Then, we can divide the robot's path into L 1 , L 2 , and L 3 based on the two corners. Figure 2 shows a graph of the robot's moving speed and travel time at that moment. The three sections we divided are further divided into three parts again (acceleration (A), constant (U), and deceleration (D), respectively) to follow the above assumptions applied for the horizontal movement segment.


The robot remains stationary at the starting point before movement occurs.  The robot immediately accelerates after starting, and when it reaches an appropriate speed, the robot moves while maintaining this speed.  The robot decelerates ahead of its destination and stops upon reaching the destination.  The robot decelerates to the predetermined speed ahead of a corner and accelerates after rounding the corner.  The acceleration and deceleration ratios (the acceleration and deceleration rates per unit time) of the robot remain constant.
To explain how the travel time is estimated in this research, certain situations were assumed. Figure 1 shows the example of indoor robot movement from origin to destination in a multi-story building. First, let us assume the robot only moves on a single floor ( floor) to explain the calculation procedure for the horizontal movement. There are two corners between the starting node and the destination node, the elevator. Then, we can divide the robot's path into , , and based on the two corners. Figure 2 shows a graph of the robot's moving speed and travel time at that moment. The three sections we divided are further divided into three parts again (acceleration ( ), constant ( ), and deceleration ( ), respectively) to follow the above assumptions applied for the horizontal movement segment.   The mathematical notations in Figures 1 and 2 are defined as follows: The rate of acceleration; The constant speed; The rate of deceleration; The portion of the moving distance between two nodes for = 1,2,3; The moving distance at the different rates for = 1,2,3 and = , , ; The travel time corresponding to the moving distance for = 1,2,3 and = , , ; The safe robot speed obtained via deceleration to safely round a corner; The maximum speed that a robot can maintain to travel through a corridor.
The following basic ideas were applied to calculate the travel time of the robot between two nodes.
1. First, the travel time is calculated based on the acceleration and deceleration sections. The portion of the moving distance between two nodes for i = 1, 2, 3; L i j The moving distance at the different rates for i = 1, 2, 3 and j = A, U, D; T i j The travel time corresponding to the moving distance L i j for i = 1, 2, 3 and j = A, U, D; V sa f e The safe robot speed obtained via deceleration to safely round a corner; V max The maximum speed that a robot can maintain to travel through a corridor. The following basic ideas were applied to calculate the travel time of the robot between two nodes.

1.
First, the travel time is calculated based on the acceleration and deceleration sections. This can be easily calculated from zero to a specific speed (in the acceleration section) or from a specific speed to zero (in the deceleration section). 2.
Next, the distance between the acceleration and deceleration sections is calculated based on the travel time in the acceleration and deceleration sections. As shown in Figure 2, we can easily calculate the moving distance by multiplying the travel time and speed. 3.
The distance encompassing the constant-speed section can be determined by subtracting the distance between the acceleration and deceleration sections from the total distance. 4.
Finally, the travel time in the constant-speed section can be calculated by dividing the length of the constant-speed section by the constant speed.
In any section composed of a corridor and a corner, the above procedure can be applied repeatedly. With the use of the notations and procedure defined above, the moving distances and travel times for the given example can be calculated, as summarized in Table 2. Table 2. Moving distance and travel time of the robot for the given example (horizontal).
Moving distance for the acceleration and deceleration sections: 3. Moving distance for the constant-speed section:

Travel time for the constant-speed section:
Travel time for the acceleration and deceleration sections: Moving distance for the acceleration and deceleration sections: 3. Moving distance for the constant-speed section:

Travel time for the constant-speed section:
Travel time for the acceleration and deceleration sections: Moving distance for the acceleration and deceleration sections: 3. Moving distance for the constant-speed section:

Travel time for the constant-speed section:
With the use of the equations in Table 2, the robot travel time required for horizontal movement (t h ij ), considering the given number of corners (c ij ) between two nodes i and j, can be generalized as follows: where: Once the travel time required for horizontal movement is calculated with Equation (1), the vertical movement segment (L 4 ) of the robot should be considered because most buildings in cities are multistory buildings. Usually, robots can move between different floors by taking the elevator, as shown in Figure 1.
As expected, it is difficult to calculate the travel time related to the elevator due to its variability. More specifically, it is difficult to capture the elevator waiting time and the number of intermediate stops. Although we can simply calculate the travel time required for vertical movement by dividing the height by the average elevator speed, similar to other studies [22][23][24], there remains a gap between the mathematical model and the actual situation.
Therefore, we adopted three scenarios (i.e., peak, normal, and off-peak scenarios) and binomial distribution assumptions to address these scenarios.
Before explaining how to calculate the moving distances and travel times required for vertical movement under the various scenarios, the following notations are defined: h The distance between floors; L 4 The moving distance of vertical movement between the n th floor and m th floors; w ω The average elevator waiting time under scenario ω, for ω ∈ {peak, normal, off-peak}; s The time spent at a stop;  Figure 3.
Here, since it is a matter of choosing whether to stop or not to stop during vertical movement on an elevator from the n th floor to the m th floor, which can be regarded as a Bernoulli trial, the number of intermediate stops can be expressed with the following binomial distribution. where: the maximum possible number of intermediate stops is r = |m − n| − 1, the stop probability under scenario ω is determined as  The following basic ideas are applied to calculate the travel time of the robot between the two floors in the above example: The following basic ideas are applied to calculate the travel time of the robot between the two floors in the above example: 1.
First, the density of the building is assessed, and the number of floors to be traversed by the robot is determined, after which the moving distance associated with vertical movement between the n th and m th floors is calculated. 2.
Next, the journey time (T 4 ) is calculated to traverse the distance between floors on an elevator considering the probability of intermediate stop occurrence.

3.
Finally, the total travel time associated with vertical movement between the n th and m th floors can be determined by summing the initial elevator waiting time (w ω ), the time required for two fixed stops (s), and the journey time (T 4 ) calculated in step 2.
With the use of the notations and procedure defined above, the moving distance (L 4 ) and travel time between the two floors (the n th and m th floors) can be calculated as follows: In Equation (4), the mean value of the intermediate stop occurrence is calculated by multiplying r and p ω according to the properties of the binomial distribution.
As a result, t ij , i.e., the travel time between any two nodes i and j within the building, can be calculated by summing the results obtained with Equations (1) and (4). Note that Equation (4) is only needed when two nodes are located on different floors.

Indoor Robot Routing Problem (IRRP)
An indoor network can be defined as a directed graph with a reasonable travel time. Let G = (V, A) be a graph where V = {0, 1, . . . , n} comprises a vertex set, and A = {(i, j)|i, j ∈ V, i = j} is an arc set. The travel time, t ij , is associated with the arc set, A. Additionally, vertex 0 represents a depot at which K homogeneous robots are housed, while the remaining vertices correspond to n customers. Each customer exhibits a nonnegative delivery quantity, and each robot exhibits a fixed capacity.
To formulate the IRRP as an MTVRP model, the mathematical model reported in [25] was modified.
The operating conditions are assumed as follows: • One or more robots are allowed to be launched together.

•
The robot starts and ends at the depot.

•
Multiple trips are allowed, as all demands must be satisfied.
• Since robot storage can be retrofitted to accommodate two or more compartments, a single robot can visit more than one customer up to the maximum capacity. • Robots are considered fully autonomous and equipped, which suggests that loading and unloading operations do not require much time, which can be neglected.

•
The robot possesses a sufficient battery capacity to cover one building.
Before we formulated the model, the following parameters and variables are defined: Objective function (5) minimizes the total travel time of the robots. Constraints (6)-(8) are flow conservation constraints, ensuring that each robot trip starts and ends at the depot. Additionally, these constraints guarantee that each customer is served exactly once. Constraint (9) defines the carrying load quantity flows, and Constraint (10) enforces capacity limitations on the robots. When a robot visits customer i, it should at least satisfy the demand of customer i, but the robot cannot carry in excess of the capacity. Note that these two constraints also eliminate subtours [26]. Finally, Constraints (11) and (12) specify the decision variables, where x k ij is defined as a binary variable and can only equal 0 or 1, and the carrying load quantity u k i is ensured to be nonnegative.

Heuristic-Based Solution
As mentioned in Section 1, commercial optimization solvers cannot always provide the optimal solution for large-sized MILP problems within reasonable computation times. The same problem should be addressed in our case because the problem size can easily become large when considering more floors and customers within buildings. Therefore, an efficient heuristic is needed to address the computational tractability in actual applications, which can quickly generate optimal or near-optimal solutions in terms of the objective value of the IRRP. In this regard, we proposed a heuristic comprising two basic steps: (1) initial solution generation using nearest-neighbor insertion and (2) solution improvement with the GA. The GA is a popular metaheuristic method that can find solutions for NP-hard problems within a reasonable time [27]. This study employed the floor-splitting heuristic to determine whether the initial solution obtained in Step 1 can be updated with the GA. Additionally, certain well-known and validated approaches, such as roulette wheel selection, partially mapped crossover operator for reproduction, and improving rate-based swap operator for mQtation, were adopted.
Each step of the proposed heuristic approach can be explained as follows:

•
Step 1: Initial solution generation with the nearest-neighbor insertion heuristic.

Step 1.1 (preparation):
The travel time is calculated from the depot to every node with the travel time matrix obtained with the travel time calculator presented in Section 3.1 (as shown in Figure 4a).
Step 1.2 (starting node selection): From the current set of nodes (an initial set containing all nodes), the node with the longest travel time from the depot is chosen. However, if there are no more nodes in the set of nodes, Step 1 is terminated.

(loading capacity check and trip generation):
If the loading capacity of the robot is not violated by considering the demand of the chosen node, the chosen node is added to the current trip (initially, there are no nodes constituting the trip), and the chosen node is eliminated from the set of nodes for update purposes. The process then proceeds to Step 1.4. Conversely, this trip is terminated by inserting the depot (i.e., node 0) at both the beginning and end of this trip, and the chosen node is eliminated from the set of nodes for update purposes. The process proceeds to Step 1.2.
Step 1.4 (nearest-neighbor insertion): The nearest (closest) neighbor node is selected among the nodes chosen during Step 1.2 in terms of the travel time with the travel time matrix, as shown in Figure 4b. The process proceeds to Step 1.3 with this chosen neighbor node. Figure 4c shows an initial solution containing five trips generated by Step 1. Robots are assigned to these trips once a final solution is obtained in Step 2.

•
Step 2: Solution improvement using the GA.
Step 2.1 (initial population creation): In this research, the chromosome was represented similarly to the initial solution, as shown in Figure 4c, and the population consisted of one hundred chromosomes. The initial solution obtained in Step 1 was added to the population first. Then, the other ninety-nine chromosomes were randomly generated.
Step 2.2 (roulette wheel selection for genetic operations): Two parent chromosomes must be selected to generate offspring that could potentially represent better solutions. The chromosome with the highest fitness value (i.e., the lowest total travel time, which is the objective function value of the IRRP) was chosen as one of the two parent chromosomes. The fitness value of any chromosome can be easily calculated by adding the travel times of all trips and determining the inverse value. The other parent chromosome can be chosen with the well-known roulette wheel selection method considering the following probability (p i ) for each chromosome i [28]: where f i →the fitness value of chromosome i.
taining all nodes), the node with the longest travel time from the depot is chosen. However, if there are no more nodes in the set of nodes, Step 1 is terminated.
Step 1.3 (loading capacity check and trip generation): If the loading capacity of the robot is not violated by considering the demand of the chosen node, the chosen node is added to the current trip (initially, there are no nodes constituting the trip), and the chosen node is eliminated from the set of nodes for update purposes. The process then proceeds to Step 1.4. Conversely, this trip is terminated by inserting the depot (i.e., node 0) at both the beginning and end of this trip, and the chosen node is eliminated from the set of nodes for update purposes. The process proceeds to Step 1.2. Step

(nearest-neighbor insertion):
The nearest (closest) neighbor node is selected among the nodes chosen during Step 1.2 in terms of the travel time with the travel time matrix, as shown in Figure 4b. The process proceeds to Step 1.3 with this chosen neighbor node. Figure 4c shows an initial solution containing five trips generated by Step 1. Robots are assigned to these trips once a final solution is obtained in Step 2.


Step 2: Solution improvement using the GA.
Step 2.1 (initial population creation): In this research, the chromosome was represented similarly to the initial solution, as shown in Figure 4c, and the population consisted of one hundred chromosomes. The initial solution obtained in Step 1 was added to the population first. Then, the other ninety-nine chromosomes were randomly generated.
Step 2.2 (roulette wheel selection for genetic operations): Two parent chromosomes must be selected to generate offspring that could potentially represent better solutions. The chromosome with the highest fitness value (i.e., the lowest total travel time, which is The proportion of the fitness value of each chromosome to the sum of the fitness values of all chromosomes in the whole population can be computed except for the best option. This represents the probability of the selection of an individual chromosome. A roulette can be partitioned according to the cumulative sum of p i , starting from the solution with the highest value. Then, a fixed point is randomly generated between zero and one. Since the area of each chromosome was proportional to each probability to be selected, a chromosome with a corresponding generated fixed point was finally chosen as the second parent chromosome. Step 2.3 (genetic operations): Once the two parent chromosomes were selected via the roulette wheel selection in Step 2.2, genetic operations such as crossover and mutation operations were applied. Genetic operations play an important role in searching for a broader solution space. The partially mapped crossover (PMX) operator and improving rate-based swap mutation (ISM) operator were adopted.
Step 2.3.1: PMX operator for the crossover operation. The crossover operator reproduces child chromosomes. However, a conventional n-point crossover is not suitable to generate child chromosomes in our case because every node must be visited once in the route plan. Therefore, the PMX operator was chosen as the crossover operator [29]. After choosing two random cut points on the parent chromosomes, the segment between the cut points was exchanged. One arbitrary point was generated for each parent and employed to establish a segment. This is called the mapping section. Then, the remaining parts with no conflicts (i.e., exhibiting no shared nodes) were maintained at their position. However, if a conflict occurs (i.e., the same node applies), the gene from the mapping section fills the position. The crossover operator created two new child chromosomes from the two parent chromosomes. Details on the PMX operator can be found in [29].
Step 2.3.2: ISM operator for the mutation operation.
The mutation operator widely expands the exploration space by introducing attributes not present in the parent chromosomes. In this research, we employed swap mutation, which generates a certain probability for each gene, thus replacing it with another gene if the probability is lower than the mutation rate. What is important when selecting the other gene is that the demand must remain the same to preserve the length of the trip and the feasibility of the solution. For example, if the selected gene exhibits a demand of 2, all the genes with a demand of 2 from the chromosome are chosen to compile a set. From the established set, one gene is randomly chosen to implement the swap mutation operation. When initiating the GA, the initial mutation rate is set to 0.2, which is increased to 0.3 from the 50th generation on and to 0.5 from the 100th generation on to avoid falling into a local optimum [30].
Through the above genetic operations, two child chromosomes are obtained. Among the four chromosomes, i.e., the two parent chromosomes and two child chromosomes, the best chromosome with the lowest objective function value is chosen and sent to the new population of the next generation. Until the new population is fully filled, both Steps 2.2 and 2.3 are repeated. Once the new population is created, the next generation starts.
In this research, the GA was subject to stopping conditions such as 200 total iterations or 30 iterations with an unchanged solution. Finally, once the final solution was obtained with the proposed heuristic approach, we converted the final solution into the optimal robot route plan. Figure 5 shows an example of how the final solution can be converted into the final robot route.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 robot route plan. Figure 5 shows an example of how the final solution can be conve into the final robot route. As shown in Figure 5, trips are assigned to both robots considering their availab For example, Trip 3 can be assigned to Robot 2 since Robot 1 is still delivering pa according to Trip 1. The overall flow of the proposed heuristic is shown in Figure 6. As shown in Figure 5, trips are assigned to both robots considering their availability. For example, Trip 3 can be assigned to Robot 2 since Robot 1 is still delivering parcels according to Trip 1. The overall flow of the proposed heuristic is shown in Figure 6. As shown in Figure 5, trips are assigned to both robots considering their availability. For example, Trip 3 can be assigned to Robot 2 since Robot 1 is still delivering parcels according to Trip 1. The overall flow of the proposed heuristic is shown in Figure 6.

Experimental Results
In this section, numerical experiments are described that validated the proposed approach. First, an illustrative example was prepared to explain indoor delivery environments by focusing on multiple floors, their layout, and robot properties. Then, the performance of the proposed heuristic approach was compared to that of the commercial optimization tool Gurobi Optimizer 9.0 (http://www.gurobi.com, accessed on 10 August 2022). All experiments were performed in a Python 3.7.1 (https://www.python.org/, accessed on 10 August 2022) environment with a 3.10-GHz Intel Core i5 processor and 8 GB of RAM (HP Elite Desktop, Intel Core i5 3.1 GHz, 8GB RAM, 500 GB HDD, Keyboard/Mouse, WiFi, 17in LCD Monitor, DVD, Windows 10 (Renewed), Korea, Incheon)

Illustrative Example
In this research, the illustrative example was based on an actual university building containing six floors, as shown in Figure 7a. Among the many rooms in the building, twelve rooms were selected as the demand nodes (i.e., customers). Figure 7b shows the indoor layout of the first floor of the building, and the remaining floors exhibit almost the same layout. We assumed that the robots have an exact indoor map and know how many corners there are. To generate demand information for this example, each demand node exhibited a nonnegative delivery demand (i.e., the number of parcels) that was randomly generated to attain a value of one or two. Additionally, the distance between the nodes was calculated based on the actual positions of the twelve selected rooms. All the detailed information is summarized in Table 3. Table 3. Customer data in the building.

Floor
Customer ID Delivery Demand To generate demand information for this example, each demand node exhibited a nonnegative delivery demand (i.e., the number of parcels) that was randomly generated to attain a value of one or two. Additionally, the distance between the nodes was calculated based on the actual positions of the twelve selected rooms. All the detailed information is summarized in Table 3.
Based on the number of robots and their capacities, it was assumed that two robots were employed to provide delivery services with three storage compartments. Although the number of robots and their capacities can vary depending on the robot type and given indoor environments, we set these values following previous studies [12,31,32]. Additionally, the average speed of the robot was 1 m/s, considering the average walking speed of people indoors [33,34]. The applied maximum acceleration, maximum deceleration, and safe speed were set to 0.3 m/s 2 , 0.6 m/s 2 , and 0.5 m/s, respectively [35].
Moreover, the indoor delivery robots are acknowledged to possess sufficient battery power for longer than several hours when fully charged [36]. In regard to the time related to the elevator, both the time spent (s) at a stop and the time spent (s ) when an intermediate stop occurs were set to 14 s in accordance with past related research [37]. Detailed informa-tion on all system parameters, including the above parameters, is provided in Table 4. Note that all the parameters can be flexibly adjusted when the indoor environments change, and these parameters will become more reliable if historical data are available. Table 3. Customer data in the building.

Floor
Customer ID Delivery Demand  Based on the parameters above, the travel time matrix among the nodes must be generated by applying the proposed travel time calculator. This matrix was derived by considering the characteristics of robot movement (e.g., acceleration, deceleration, and cornering) and the additional time required for traveling on the elevator (e.g., elevator waiting time, stops, and occurrence probability of intermediate stops) with the equations presented in Section 3.1. The derived travel time matrix assuming the normal density scenario (ω = normal) is presented in Table 5. Next, the IRRP was solved to obtain the delivery route for the two robots, and the obtained routes are listed in Table 6. The total travel time until all the demands were fulfilled by the two robots reached 4561.4 s.  However, what if the travel time matrix is generated with a simple conventional travel time calculator? Will a difference occur in terms of the total travel time from the IRRP solution? To answer these questions, we generated the travel time matrix simply by dividing the moving distance by the average speed of the robot without considering any acceleration, deceleration, or cornering aspects of the robots.
With this simple travel time matrix, the total travel time decreased to 3108 s, revealing a 1381.4 s difference (approximately 32% shorter than the proposed travel time calculator). Although the total travel time decreased with the simple travel matrix, the obtained routes to solve the IRRP can be viewed as naive and unrealistic. The difference in the total travel time should not be ignored when considering the battery power of the robot.
Moreover, the naive routes resulting from the simple travel time matrix can cause delivery failure, which can lead to customer dissatisfaction because a much longer time is required for the robots to finish the delivery. Additionally, the accumulation of the differences in the total travel time can cause a serious problem if delivery services in multiple buildings are required within a given number of hours. For instance, delivery delays are more likely to occur in buildings assigned later in the delivery schedule.

Performance Validation of the Heuristic Approach
In the previous example of two robots and a six-story building with twelve customers, Gurobi, a commercial solver, quickly obtained the optimal solution. However, if the problem size is increased by adding more customers and robots, we must assess whether Gurobi can obtain the optimal solution rapidly.
Three cases with different numbers of floors and robots were designed to evaluate and compare the performance between the proposed heuristic approach and Gurobi. The small-sized problem case consisted of six floors and two robots, and the medium-sized problem case was composed of nine floors and three robots. Finally, the large-sized problem case involved twelve floors and four robots.
It was further assumed that there can occur at least zero to at most two corners between any pair of nodes, and the normal density scenario was applied. Ten test instances were randomly generated in each of the three cases. The proposed heuristic was coded in Python 3.7.1 on the same computer. All system parameters and data generation methods followed the illustrative example above to streamline our analysis. Additionally, the computational time was limited to sixty minutes (i.e., 3600 s) because there are many buildings where delivery services should be completed on the next day. All robot route plans should be generated during the night. Thus, it seems unacceptable if more than an hour is needed to generate a plan for a single building. The objective values and computation time of both the proposed heuristic and Gurobi for the test instances are summarized in Table 7.
As indicated in Table 7, the performance of the proposed heuristic was investigated by dividing the whole process into two main steps: (1) initial solution generation via nearest-neighbor node insertion and (2) solution improvement via the GA. This approach determined the contribution of the GA-based solution improvement step to the performance of the final solution.
As demonstrated in Table 7, in the small-sized problem case, the average computation time of Step 1 of the proposed heuristic was 0.002 s, while Gurobi obtained the optimal solution within 12.61 s on average. Thus, Step 1 of the proposed heuristic seems to be very swift, but the quality of the solution still requires improvement because the average percentage gap with the Gurobi solution in terms of the objective value was 5.78%. In other words, the quality of the initial solution resulting from Step 1 is not yet satisfactory and requires improvement. That is why Step 2 of the proposed heuristic was developed to improve the solution.
Step 2 of the proposed heuristic certainly took longer, namely, 52.89 s on average, but it obtained the exact solutions to the Gurobi solutions in every instance. Although it is more effective to employ Gurobi for small-sized problems in terms of the computation time, the need for heuristics increased in regard to medium-and large-sized problems considering both the computation time and solution quality.
As indicated in Table 7, the proposed heuristic required 143.16 s (for the medium-sized problem) and 200.09 s (for the large-sized problem) on average to derive the solutions. In contrast, Gurobi required 813.94 s (for the medium-sized problem) and 2196.32 s (for the large-sized problem) on average and sometimes failed to obtain a solution within the running time (e.g., instances 14 and 17 for the medium-sized problem and instances 22,24,26,27,29, and 30 for the large-sized problem).
Additionally, the proposed heuristic could obtain near-optimal solutions in the mediumand large-sized problem cases. The average optimality gap between the proposed heuristic and Gurobi in the medium-and large-sized problem cases reached 0.11% and 0.18%, respectively.
In summary, the proposed heuristic seems to be a viable tool to generate a suitable and feasible route more rapidly than commercial optimization software with increasing problem size.

Conclusions
During the COVID-19 pandemic, the growth of non-face-to-face services has been steep, and many new innovations have occurred in logistics. A typical example is last-mile delivery automation using robots. Currently, this process occurs in the test phase, and there are still a few steps left to commercialize, but this situation suggests a new normal, not a one-time happening attributed to the pandemic [38,39]. We expected that robot delivery could further constitute a probable option for last-mile delivery in urban settings with many multistory buildings, and relevant research was performed. Therefore, this research also aimed to enhance the practicality of last-mile delivery services involving indoor robots carrying multiple items. The proposed heuristic procedures might allow practitioners to generate suitable robot route plans within a reasonable time, which is helpful for the implementation of robot-based indoor delivery services in the real world.
In this research, we first investigated the importance of the travel time, which plays an essential role in solving the IRRP. If there are corners and elevators between the origin and destination points within buildings, the impact of corners and elevators should be considered in route planning. Next, three cases were assessed by changing the number of floors, robots, and customers to evaluate the validity of the proposed heuristic: small-(6 floors and 2 robots), medium-(9 floors and 3 robots), and large-sized problems (12 floors and 4 robots). The experimental results indicated that the proposed heuristic could generate satisfactory solutions (i.e., 0%, 0.11%, and 0.18% optimality gaps on average, respectively) within a short time in each case. Moreover, solutions were all obtained much faster than with Gurobi, a commercial optimization software, even in the large-sized problem case, which could not obtain optimal solutions in 6 out of 10 problem instances within an hour.
Regarding future research directions, the solution obtained with the proposed approach can be validated via test simulations using robotics middleware frameworks (e.g., ROS 2 with Gazebo and Navigation 2). In addition, the mathematical model and heuristic can be revised to manage dynamic indoor situations such as avoiding obstacles and sudden reliability problems (e.g., delivery failure due to breakdown or battery discharge). They are important issues to resolve when implementing this study in real-world environments. Moreover, crowdedness, another uncertain factor of buildings, might be an important topic for further studies. Considering the crowdedness of hallways could upgrade our travel time calculator more realistically.