# A Bi-Level Probabilistic Path Planning Algorithm for Multiple Robots with Motion Uncertainty

Authors: Jingchuan Wang; Ruochen Tai; Jingwen Xu
Journal: Complexity
 (2020)

---

## Abstract

For improving the system efficiency when there are motion uncertainties among robots in the warehouse environment, this paper proposes a bi-level probabilistic path planning algorithm. In the proposed algorithm, the map is partitioned into multiple interconnected districts and the architecture of proposed algorithm is composed of topology level and route level generating from above map: in the topology level, the order of passing districts is planned combined with the district crowdedness to achieve the district equilibrium and reduce the influence of robots under motion uncertainty. And in the route level, a MDP method combined with probability of motion uncertainty is proposed to plan path for all robots in each district separately. At the same time, the number of steps for each planning is dependent on the probability to decrease the number of planning. The conflict avoidance is proved, and optimization is discussed for the proposed algorithm. Simulation results show that the proposed algorithm achieves improved system efficiency and also has acceptable real-time performance.


---

## Body

## 1. Introduction

In the running of multiple robots in warehouses, motion uncertainties happen unavoidably since there exist some disabled components of robots, the instability of networks, and the interference of people walking. Under this case, robots would not follow the designed paths and the coupled relationship between temporal and space domain for paths is broken. And there is no doubt that other robots are disturbed by the ones where motion uncertainties happen. Finally, the whole multirobot system would be influenced and this brings about chaos or even breakdown of the whole system. It can be seen that this poses a great challenge to the path planning of multiple robots in warehouses. Therefore, taking motion uncertainty into consideration in the path planning of multiple robots is an issue worthy of attention and research.As a kind of method taking motion uncertainty into consideration for the discretized map, Ma et al. [1] designed the path planning algorithm which contains relevant restricted principles that an agent cannot enter the vertex that is occupied by some other agents at the last time step. The conflicts can still occur and this would cause the breakdown of the multirobot system in warehouses. Another kind of method which can deal with the issues caused by motion uncertainties needs to replan the paths of multiple robots when the samplings of robots contain abnormal information. Li et al. [2] proposed a two-stage trajectory planning scheme where a directed graph with cost variable connection is established and cost variable depends on the reachability, the risk of collision, and different state constraints. Gombolay et al. [3] presented a centralized algorithm that handles tightly intercoupled temporal and spatial constraints and scales to environment with large number of robots, which aims at improving the efficiency firstly. However, the above method cannot guarantee the real-time performance. In [4], the authors proposed a hybrid architecture which combines centralized coordination with distributed freedom of action to achieve an appropriate interplay. The architecture takes the new sampling for the state of robots at each time step, and path planning algorithm would be adopted to coordinate paths of multiple robots when sampling information containing motion uncertainties . However, the real-time performance cannot be guaranteed and the efficiency of the whole system is reduced due to too frequent planning. The algorithm introduced in [5, 6] based on generalized probabilistic roadmaps (GPRM) also coordinates paths of multiple robots under motion uncertainty by using the feedback information. In the algorithm, the passive planning strategy is introduced in the multiple traveling salesman problem (MTSP) solution methodology. However, the system efficiency is low and the scalability of the algorithm cannot be guaranteed because far more points need to be sampled. A novel method in [7] based on the sampling considering the motion uncertainty would propagate this delay. Robots’ paths need to be replanned adopting the Pareto-optimal plan repair scheme. However, the problem concerned with the decreased system efficiency is still not solved.As a classical method devoted to solving the uncertainty problem, the Markov decision process (MDP) is introduced to the field of multiple agents. Corresponding multiagent MDPs (MMDPs) have been designed for multiple agent planning. In [8], approximating methods are introduced to be combined with MMDPs for the path planning problem of a team of homogeneous microair vehicles (MAVs) towards a set of goals. This algorithm can guarantee that there are no conflicts. Ma et al. [9] proposed the use of robust plan-execution policies to control how each robot proceeds along its path with a 2-level MDP solver that generates valid plans. To optimize the solutions, a novel optimal solver is designed in [10] for the problem of motion uncertainties in the transition-independent MMDPs. In [11], an equilibrium policy has been proposed in MMDPs. Based on the notion of macroactions, the path planning problem is transferred into the decentralized partially observable Markov decision processes (Dec-POMDPs) in [12, 13]. Multiagent reinforcement learning (MARL) [14–20] solves the learning task in multiagent systems (MASs). Nevertheless, the reasoning agent would deal with learning tasks which grow exponentially according to the number of agents. Therefore, the computational complexity of these methods based on MDPs is so high that the real-time scalability performance cannot be guaranteed.Recently, the planning methods based on homotopy classes have been proposed to solve the problem of motion uncertainties. In [21, 22], the coordination space is established, which is the base for the later path planning. By selecting homotopy classes of the paths which are collision-free, a centralized controller is designed for the path planning problem under motion uncertainties. Furthermore, the encoding based on priorities, which is called the priority graphs, is introduced for the homotopic solutions to solve the coordination problem in [23]. For [21–23], the planning method based on homotrophy classes can be proved that there are no conflicts for the coordinating paths and make sure that all robots will eventually reach their destinations even when some robots are temporarily stopped by a delay disturbance. However, the quality of solutions is decreased to a great extent since this kind of method only considers the next step of motions and all the robots would follow the previous planned paths.Although the above methods ensure that conflict can be avoided among the planning of multiple agents under motion uncertainty, problem still exists. The efficiency of a multirobot system with motion uncertainty cannot be guaranteed. Therefore, the first objective of this paper is to design an algorithm which can plan conflict-free paths for multiple robots even when motion uncertainty happens. And another objective is to improve the system efficiency that includes the makespan of the whole system and the number of planning.To attain these goals, in this paper, the map is partitioned into multiple districts and we propose a bi-level planning architecture that is composed of topology level and route level generating from above map.In the topology level, the passing orders of districts are planned based on the weight of paths in topology map. And in this process, the weight changes dynamically according to the distribution of robots; the aim of planning in the topology level is to achieve the district equilibrium, which can in turn reduce the influence of robots under motion uncertainty and promote the system efficiency. In the route level, a MDP-based probabilistic method is proposed to plan path for all robots in each district separately. The MDP is combined with the probability of motion uncertainty to improve efficiency. And the number of steps for each planning is dependent on the probability to decrease the number of planning. We proved that the designed algorithm is conflict-free. And we also validated the system efficiency improvement of designed algorithm from the optimization discussion and simulation.The remainder of this paper is organized as follows. Section2 provides the problem statement. In Section 3, the path planning algorithm for multiple robots is introduced. Section 4 proves that the proposed algorithm is conflict-free and discusses the optimization. Simulation validation results are reported in Section 5. Section 6 concludes our research.


## 2. Problem Formulation

The problem that this paper attempts to solve is the path planning of multiple AGVs with motion uncertainty in warehouses. In this problem, warehouse environment is discretized and simplified as a graphGV,E whose vertices V correspond to locations and whose edges E correspond to transitions between locations. Figure 1 shows the discretized map.Figure 1
The discretized map of the warehouse. The orange and blue squares represent pickup points of the shelves and storage points, respectively. Other squares represent feasible ways.LetR=r1,r2,…rn be a set of robots. Time is discretized into time steps. To formally describe a plan, the path for a robot ri is a map pi which means that the location pit of ri in each time step tϵN∗ belongs to the vertices V. A robot can either remain stationary or move to an adjacent vertex. Robots are allowed to move simultaneously but restricted to move on the same edge in opposite directions or move at the same vertex simultaneously to ensure that there are no conflicts among multiple robots. Figure 2 gives the illustrations of these restricted conditions.Figure 2
Restricted condition for motions of multiple robots in warehouses.
(a)(b)LetT=t1,…tv be a set of tasks which are generated continually. For each robot that has been assigned a task, it has a start position (pickup point) and a goal position (storage point). And the makespan of all the tasks is M. During the running of robots, motion uncertainty, which means that robots have to stop immediately, could happen at any time. And the probability of the uncertainty is pu.This paper assumes the following about the scheduling of multiple AGVs in warehouses:(i)
All feasible routes are two-way bidirectional streets, which means that two robots can drive side by side and robots can drive in both directions on each feasible way.(ii)
The warehouse is made up of shelves and feasible ways, which is shown in Figure1 (not limited to this scenario).(iii)
The locations of all the robots are known.(iv)
The probability of the uncertainty is known.(v)
The robots in the state of delay will return to normal. And after the recovery, these robots will be put into operation in the warehouse system.(vi)
The task allocation is not taken into consideration in this paper.In this paper, letters of the subscripts which are not in the parentheses are used for identification. Letters or numbers of the subscripts which are in the parentheses are used for index.Therefore, the problem studied in this paper can be formulated in the following equationfd1:(1)minMs.t.pit≠pjt+1,∀i,j∈1,n,i≠j,pjt≠pit+1,∀i,j∈1,n,i≠j,pit≠pjt,∀i,j∈1,n,i≠j.


## 3. The Bi-Level Path Planning Algorithm for Multiple Robots

In this section, the bi-level path planning algorithm is introduced. The aim of this structure is that district equilibrium would be taken into consideration in the topology level for reducing the influence of robots under motion uncertainty and distributed path planning in each district can improve system efficiency in the route level even when motion uncertainty occurs. The architecture of the algorithm is shown in Figure3. It is divided into two levels generating from the map partition: topology level and route level. In the topology level, the algorithm plans the order of passing districts. In the route level, the algorithm plans paths for robots in each district separately, considering the motion uncertainty. And the introduction of map partition is as follows.Figure 3
The bi-level architecture of the algorithm.


### 3.1. Map Partition

In this paper, the discretized map in Figure1 is partitioned into multiple districts D=d1,…,dM as shown in Figure 3. For the topology level, the generation of the topology map G′V′,E′ is listed as follows: (1) each vertex v′∈V′ corresponds to district ds; (2) each edge vs′,vs′′∈E′ represents that districts ds and ds′ are connected. Let the horizontal edges and vertical edges be Eh′ and Ev′. For the route level, all feasible ways of each district in the discretized map of Figure 1 are shown in the route map.


### 3.2. Topology Level

In the topology level, the first step of the planning problem introduced in Section3 is to plan the passing orders of districts for the robots which have been assigned tasks. At each time step, there might be multiple tasks assigned to robots. The number of these tasks is N and the robots corresponding to these tasks are Rt=r1,…,rtt≤n. According to the task ti which has been assigned to the robot ri, the districts which contain the start position and goal position of ti are dsiand dgi, respectively (here si and gi represent the index of the corresponding district). To formally describe the plan results, the passing order of districts for a robot ri is a map oi:N∗⟶D.In this paper,TT∈N∗ is defined as the macroscopical time, which is the independent discrete variable in the topology level. With the increase of T, the robot will move from the start district dsi to the goal district dgi. Different from the discrete time step, the macroscopic time represents the temporal space and corresponds to the district in topology level. Therefore, the change of T represents the transfer of the robot in the topology level. And if districts oiT and oiT+1 are different, then the corresponding edge that robot ri will pass in the graph G′V′,E′ is defined as followsfd2:(2)EiT=oiT,oiT+1.In this formulation, the edgeEiT represents the transition from oiT to oiT+1, which is the moving direction for the robot ri. For the robot ri which has been assigned the task, there is the corresponding macroscopical time Tsi and Tgi when ri reaches the start district dsi and goal district dgi, respectively. Then, the planning problem concerned with the districts can be formulated as follows:(3)min∑i=1t∑T=TsiTgiLEiT1+CT,EiTs.t.oiTsi=dsi,1≤i≤t,oiTgi=dgi,1≤i≤t.In this formulation of the optimization objective,L⋅ is the length of the edge EiT, which is the corresponding distance between the center of the two districts. Due to the length difference of the edges in the graph G′V′,E′, the durations of passing theses edges are different. In our problem, the length of the horizontal edge is twice that of the vertical edge. Therefore, the duration of passing the horizontal edge is twice that of the vertical edge. For the other environment map with different layout, the same principle applies. According to this principle, the passing order for the robot ri should satisfy the following:(4)oiT+1=oiT,oiT−1,oiT∈Eh′andoiT−2≠oiT−1,next district,otherwise.Based on this principle, equation (2) can be rewritten as(5)EiT+1=EiT,EiT−1∈Eh′andEiT−2≠EiT−1,next edge,otherwise.In the optimization objective of equation (3), C⋅ represents the crowdedness of the edge EiT at macroscopic time T(here we assume that oiT and oiT+1 are different). And it is defined as follows:(6)CT,EiT=NCoiTT+NCoiT+1T2,where NCoiTT and NCoiT+1T are the normalized crowdedness of the district oiT and oiT+1, respectively. And the crowdedness of each district oiT is denoted as CoiTT, which represents the number of robots in each district at the macroscopic time T. It can be obtained by the passing orders for each robot. Furthermore, these values are changing dynamically since robots are moving. Then, the method of normalization can be adopted to calculate NCoiTT and NCoiT+1T.Furthermore, due to the motion uncertainty of robots themselves, robots would not arrive at the planned district precisely. Therefore, at eachmacroscopic time, there are multiple robots which are assigned new tasks and need new path planning in current warehouse environment plan, the normalized crowdedness in equation (6) should be recalculated according to the distribution and passing orders of all robots in topology level. To be specific, assuming that one robot is located at district ds1 and will move to district ds2 when there are some new tasks and the planning is needed at the macroscopic time T, if the edge corresponding to the transition between ds1 and ds2 belongs to Eh′, then this robot would still stay at this edge at next time T+1 in equation (6). If the edge belongs to Ev′, then this robot would move to next edge at macroscopic time T+1.Therefore, the optimization objective in equation (3) considers the length of the edges which robots would pass and the crowdedness of the related districts, whose aim is to make sure that the paths of robots would not be too long and the districts which have been planned for these robots would not be too crowded. And this design can guarantee the balance among different districts and efficiency of the whole warehousing system.The constraints listed in equation (3) make sure that the algorithm can plan paths for robots from their start position to the goal position in topology level.To solve the optimization problem in equation (3), K shortest path planning algorithm [24] is adopted to obtain the multiple shortest paths for each robot which needs to be planned at the moment. According to the planned paths, the passing order can be generated by using equations (4) and (5). The optimization objective can be calculated based on the rules introduced above, and the optimal solution among all the combinations can be selected, which contains the information concerned with the passing order for the robots assigned with new tasks in the topology level.Figure4 shows an example of the planning results in topology level. And it can be seen that the planning results follow the analysis above; the algorithm in the topology level prefers to choose the paths which are shorter and less crowded since this can minimize the optimization objective mentioned in equation (3).Figure 4
Illustration of the planning result in topology level. Red squares represent the robots which need to be planned. Green squares represent the other robots. The red line shows the planned passing order for each robot assigned with new task. According to the number of robots, the grey shade implies the crowdedness of each district. The fewer the number of robots, the lighter the color.


### 3.3. Route Level

In the bi-level architecture, when the planning of the topology level has been finished, according to the distribution of the robots, the planning of each district would be conducted for the robots in this district distributedly.In the route level, the route map is made up of all the feasible ways in each district. And the local route planning is to plan the actions of all the robots in the district being planned.In the planning of each district, considering the motion uncertainty of robots and the probabilitypu, the planning step α is introduced that should be intuitively related to pu (α not only refers to the planning step at each planning but also refers to the interval between each two plannings). The aim of this idea is that the planning step can be adjusted dynamically based on pu and some redundant planning in [11, 12] would not be needed anymore. Here, we define the planning step as the follows:(7)α=1pu.It can be seen from this definition that there is an inverse relationship between the planning stepS and the probability of motion uncertainty pu. When pu gets larger, which means that there is higher probability that robots would not follow the previous planned path, the planning step α becomes smaller and the algorithm would plan paths for robots more conservatively. The first aim of computation for probability of motion uncertainty is to guarantee the real-time performance of the planning in route level since the traditional replanning all the paths from the current position to the goal position. However, the planning step is dependent on the probability of the motion uncertainty. Another aim of this is the planning step is adjusted based on pu, which can prevent the situation that frequent replanning might destroy the optimized planned paths.In the route level, the planning of each district can be formulated as the following optimization problem:(8)minπxERS1+λ⋅RS2+⋯+λα−1RSα,s.t.pit≠pjt+1,∀i,j∈1,m,i≠j,pjt≠pit+1,∀i,j∈1,m,i≠j,pit≠pjt,∀i,j∈1,m,i≠j.In the formulation of this optimization objective,πx represents one policy and there are k policies which are π1,π2,…,πk. Each policy contains the information concerned with actions of the robots at each time step in the district. And for each policy, πx=πx,1,πx,2,…,πx,α where πx,1 indicates the policy at the first time step for all the robots in the district being planned. And πx,1=πx,1r1,πx,1r2,…,πx,1rmm≤t where πx,1r1 indicates the policy at the first time step for the robot r1 and m is the number of the robots located in the district being planned. And πxr1=πx,1r1,πx,2r1,…,πx,αr1 represent the α step policy of r1.E⋅ indicates the expectation of the reward which is similar to that of MDP. The aim of designing this objective function is that this can obtain the probabilistic optimal solution since there exists motion uncertainty of robots. λ0<λ<1 is the discount factor. λ has relationship with pu as equation (9). The aim of equation (9) is that the discount factor can be adjusted dynamically based on and when is high; the long-term reward is not reliable anymore and the path planning based on MDP is concerned more on the optimization of total short-term reward.S1,…,Sα represents the state from the start time step when the current district needs to be planned to the time step which is α−1 time steps after the start time step. And for S1, it contains states for all the robots in the district, which means that S1=s1r1,s1r2,…,s1rm. The state s1r1 indicates the state for the robot r1after the action π1,1r1 is taken. R⋅ defines the reward for each state based on the distribution of the robots in the corresponding district.(9)λ=1−pu.Therefore, according to the Markov decision process in [8], the optimization objective in equation (8) implies the expectation of the rewards for the state from the start time step to after α−1 time steps. And it can be rewritten as follows:(10)minπxERS1+λ⋅RS2+⋯+λα−1RSα=minπx∑q=1αλq−1∑Sq′PSq−1πx,q−1Sq′⋅RSq′,where Sq′ indicates the possible state after the action πx,q−1 is taken at the state Sq−1 since there exists the motion uncertainty of robots, RSq′ represents the reward of the robots in this district at the state Sq′, and PSq−1πx,q−1Sq′ is the probability that the state of robots in this district would turn to Sq′ when the action πx,q−1 is taken at the state Sq−1. The probability can be calculated as follows:(11)PSq−1πx,q−1Sq′=∏j=1mPsq−1rj,πx,q−1rjsqrj′,where sq−1rj′is the state of the robot rj after the action πx,q−1,rj is taken at state sq−1rj and Psq−1rj,πx,q−1,rjsqrj′ is the probability that the state of the robot rjwill change into sqrj′ after the action πx,q−1rjis taken at state sq−1rj. And this probability is dependent on the motion uncertainty and can be represented as follows:(12)Psq−1rj,πx,q−1rjsqrj′=pu,πx,q−1rj=move andrjstays still,1−pu,πx,q−1rj=move andrjmoves,1,πx,q−1rj=waiting.For the rewardRSq′ of the policy πx, it is related to the crowdedness of the district and the distances among robots in the district, which can be defined as follows:(13)RSq′=Cdx′⋅ε⋅∑h=1mvaluesqrh′,where Cdx′ is the normalized crowdedness of the district being planned when the policy πx is taken. And the crowdedness Cdx of the district being planned can be calculated as follows:(14)Cdx=m∑i=1m∑j=1mdistanceij,where distanceij is the distance between ri and rj at state Sq′. Then, the crowdedness Cdx needs to be normalized by using the principle of the softmax according to the following formula:(15)Cdx′=eCdx∑h=1keCdh.In equation (13), valuesqrh′indicates the value of the robot rh at state sqrh′. It is defined as the distance between the robot rhand the next district which rhheads to since the passing order has been planned in the topology level. The rule for calculating this distance is based on the Manhattan distance.To avoid the conflicts among multiple robots, the factorε is introduced. And ε in equation (13) is defined as follows:(16)ε=∞,pit=pjt+1,pjt=pit+1orpit=pjt,∃i,j∈1,m,e1/1−pu,pit=pjt,∃i,j∈1,m,1,otherwise.It can be seen from equation (16) that when there would be a conflict between robot ri and rj, ε is set at ∞. When a robot would move to the location at next time step where another robot is staying at for now, then ε is set at e1/1−pu. Otherwise, ε is set at 1.According to equations (13) to (16), the first aim of this design is to prevent that each district would be too crowded since this would cause the coordination space to be too small, which means that once there happens motion uncertainty in some robots, it would be hard for the algorithm to find the solution and the system efficiency would be decreased. Another aim for the set of ε is to prevent the conflicts and make sure that one robot would not choose the action which might cause it to be next to another one. And this is beneficial for the motion coordination of robots since the motion uncertainty of one would cause it to still stay at the previous location, which increases the probability of the conflicts. Therefore, this setting can decrease the probability of the situation like that.To solve the optimization problem in equation (8), the policy set π needs to be generated in the district and the policy with the minimum objective would be selected as the solution. And the principle for the candidate action set of the policies is that robots would not be far away from the next district after taking this action. There are two benefits for this principle. Firstly, the farther the distance from the next district, the larger the objective function in equation (8). Secondly, due to the introduction of this principle, the candidate actions will decrease, which can in turn lower the computational complexity. According to above principle, there are two kinds of candidate action sets when robot is in different locations of each district:(1)
When robot is in any location of district and not coming into the boundary, the candidate action set of robot is denoted asleft,up,still. The action left means that robot drives to the other side of two-way street at next time step, up means that robot follows current side of street, and still means that robot stays still at next time step. For example, in Figure 5, robots r1 and r3 are in the district and not coming into boundary at next time step. Therefore, their candidate action set includes left, up, and still as orange arrows show. And as shown in Figure 5, the boundary between two districts is defined as an area composed of four grids, and one half of the boundary belongs to one district and the other belongs to the neighboring district. Considering that path planning in the route level is conducted separately in each district, robot should only follow one direction (drive to the right) in the process of traversing the boundary to the neighboring district, for avoiding the conflict.(2)
When robot is in the boundary of district or coming into the boundary, the candidate action set of robot is denoted asup,still or left,still according to the rule that robot should only drive to the right in the process of traversing the boundary. For example, robot r2in Figure 5 is coming into the boundary at next time step and its candidate action set isleft,still which means that it needs to drive to the right side of two-way street firstly. r4 and r5 are in the boundary and their candidate action set is up,still.Figure 5
Illustration of the candidate actions. Green squares represent the robots in the district. The blue arrow indicates the next district where these robots are heading to. Therefore, the orange arrows for each robot represent the candidate actions they can take. Besides the actions represented by these orange arrows, staying still is also included in the candidate actions for each time step.According to the above candidate action set, the generation process of robotri and policy πqri is shown in Figure 6. And it can obtain multiple policies for each robot. Then, π can be generated with the combination of multiple policies of each robot in the district.Figure 6
Generation of policyπqri for riin each district. Assuming that the candidate actions in each district include left, up, and still, the generation of policy πqri for robot ri is shown in this figure. In this tree structure, different level of branches corresponds to different steps. Each circle represents the action for ri in the district at the corresponding time step. For example, the action left means that all the robots in the district should go left. A complete policy πqri from the first time step to the αth time step is illustrated in the blue box. Multiple policies for each robot can be generated from this tree structure.


## 3.1. Map Partition

In this paper, the discretized map in Figure1 is partitioned into multiple districts D=d1,…,dM as shown in Figure 3. For the topology level, the generation of the topology map G′V′,E′ is listed as follows: (1) each vertex v′∈V′ corresponds to district ds; (2) each edge vs′,vs′′∈E′ represents that districts ds and ds′ are connected. Let the horizontal edges and vertical edges be Eh′ and Ev′. For the route level, all feasible ways of each district in the discretized map of Figure 1 are shown in the route map.


## 3.2. Topology Level

In the topology level, the first step of the planning problem introduced in Section3 is to plan the passing orders of districts for the robots which have been assigned tasks. At each time step, there might be multiple tasks assigned to robots. The number of these tasks is N and the robots corresponding to these tasks are Rt=r1,…,rtt≤n. According to the task ti which has been assigned to the robot ri, the districts which contain the start position and goal position of ti are dsiand dgi, respectively (here si and gi represent the index of the corresponding district). To formally describe the plan results, the passing order of districts for a robot ri is a map oi:N∗⟶D.In this paper,TT∈N∗ is defined as the macroscopical time, which is the independent discrete variable in the topology level. With the increase of T, the robot will move from the start district dsi to the goal district dgi. Different from the discrete time step, the macroscopic time represents the temporal space and corresponds to the district in topology level. Therefore, the change of T represents the transfer of the robot in the topology level. And if districts oiT and oiT+1 are different, then the corresponding edge that robot ri will pass in the graph G′V′,E′ is defined as followsfd2:(2)EiT=oiT,oiT+1.In this formulation, the edgeEiT represents the transition from oiT to oiT+1, which is the moving direction for the robot ri. For the robot ri which has been assigned the task, there is the corresponding macroscopical time Tsi and Tgi when ri reaches the start district dsi and goal district dgi, respectively. Then, the planning problem concerned with the districts can be formulated as follows:(3)min∑i=1t∑T=TsiTgiLEiT1+CT,EiTs.t.oiTsi=dsi,1≤i≤t,oiTgi=dgi,1≤i≤t.In this formulation of the optimization objective,L⋅ is the length of the edge EiT, which is the corresponding distance between the center of the two districts. Due to the length difference of the edges in the graph G′V′,E′, the durations of passing theses edges are different. In our problem, the length of the horizontal edge is twice that of the vertical edge. Therefore, the duration of passing the horizontal edge is twice that of the vertical edge. For the other environment map with different layout, the same principle applies. According to this principle, the passing order for the robot ri should satisfy the following:(4)oiT+1=oiT,oiT−1,oiT∈Eh′andoiT−2≠oiT−1,next district,otherwise.Based on this principle, equation (2) can be rewritten as(5)EiT+1=EiT,EiT−1∈Eh′andEiT−2≠EiT−1,next edge,otherwise.In the optimization objective of equation (3), C⋅ represents the crowdedness of the edge EiT at macroscopic time T(here we assume that oiT and oiT+1 are different). And it is defined as follows:(6)CT,EiT=NCoiTT+NCoiT+1T2,where NCoiTT and NCoiT+1T are the normalized crowdedness of the district oiT and oiT+1, respectively. And the crowdedness of each district oiT is denoted as CoiTT, which represents the number of robots in each district at the macroscopic time T. It can be obtained by the passing orders for each robot. Furthermore, these values are changing dynamically since robots are moving. Then, the method of normalization can be adopted to calculate NCoiTT and NCoiT+1T.Furthermore, due to the motion uncertainty of robots themselves, robots would not arrive at the planned district precisely. Therefore, at eachmacroscopic time, there are multiple robots which are assigned new tasks and need new path planning in current warehouse environment plan, the normalized crowdedness in equation (6) should be recalculated according to the distribution and passing orders of all robots in topology level. To be specific, assuming that one robot is located at district ds1 and will move to district ds2 when there are some new tasks and the planning is needed at the macroscopic time T, if the edge corresponding to the transition between ds1 and ds2 belongs to Eh′, then this robot would still stay at this edge at next time T+1 in equation (6). If the edge belongs to Ev′, then this robot would move to next edge at macroscopic time T+1.Therefore, the optimization objective in equation (3) considers the length of the edges which robots would pass and the crowdedness of the related districts, whose aim is to make sure that the paths of robots would not be too long and the districts which have been planned for these robots would not be too crowded. And this design can guarantee the balance among different districts and efficiency of the whole warehousing system.The constraints listed in equation (3) make sure that the algorithm can plan paths for robots from their start position to the goal position in topology level.To solve the optimization problem in equation (3), K shortest path planning algorithm [24] is adopted to obtain the multiple shortest paths for each robot which needs to be planned at the moment. According to the planned paths, the passing order can be generated by using equations (4) and (5). The optimization objective can be calculated based on the rules introduced above, and the optimal solution among all the combinations can be selected, which contains the information concerned with the passing order for the robots assigned with new tasks in the topology level.Figure4 shows an example of the planning results in topology level. And it can be seen that the planning results follow the analysis above; the algorithm in the topology level prefers to choose the paths which are shorter and less crowded since this can minimize the optimization objective mentioned in equation (3).Figure 4
Illustration of the planning result in topology level. Red squares represent the robots which need to be planned. Green squares represent the other robots. The red line shows the planned passing order for each robot assigned with new task. According to the number of robots, the grey shade implies the crowdedness of each district. The fewer the number of robots, the lighter the color.


## 3.3. Route Level

In the bi-level architecture, when the planning of the topology level has been finished, according to the distribution of the robots, the planning of each district would be conducted for the robots in this district distributedly.In the route level, the route map is made up of all the feasible ways in each district. And the local route planning is to plan the actions of all the robots in the district being planned.In the planning of each district, considering the motion uncertainty of robots and the probabilitypu, the planning step α is introduced that should be intuitively related to pu (α not only refers to the planning step at each planning but also refers to the interval between each two plannings). The aim of this idea is that the planning step can be adjusted dynamically based on pu and some redundant planning in [11, 12] would not be needed anymore. Here, we define the planning step as the follows:(7)α=1pu.It can be seen from this definition that there is an inverse relationship between the planning stepS and the probability of motion uncertainty pu. When pu gets larger, which means that there is higher probability that robots would not follow the previous planned path, the planning step α becomes smaller and the algorithm would plan paths for robots more conservatively. The first aim of computation for probability of motion uncertainty is to guarantee the real-time performance of the planning in route level since the traditional replanning all the paths from the current position to the goal position. However, the planning step is dependent on the probability of the motion uncertainty. Another aim of this is the planning step is adjusted based on pu, which can prevent the situation that frequent replanning might destroy the optimized planned paths.In the route level, the planning of each district can be formulated as the following optimization problem:(8)minπxERS1+λ⋅RS2+⋯+λα−1RSα,s.t.pit≠pjt+1,∀i,j∈1,m,i≠j,pjt≠pit+1,∀i,j∈1,m,i≠j,pit≠pjt,∀i,j∈1,m,i≠j.In the formulation of this optimization objective,πx represents one policy and there are k policies which are π1,π2,…,πk. Each policy contains the information concerned with actions of the robots at each time step in the district. And for each policy, πx=πx,1,πx,2,…,πx,α where πx,1 indicates the policy at the first time step for all the robots in the district being planned. And πx,1=πx,1r1,πx,1r2,…,πx,1rmm≤t where πx,1r1 indicates the policy at the first time step for the robot r1 and m is the number of the robots located in the district being planned. And πxr1=πx,1r1,πx,2r1,…,πx,αr1 represent the α step policy of r1.E⋅ indicates the expectation of the reward which is similar to that of MDP. The aim of designing this objective function is that this can obtain the probabilistic optimal solution since there exists motion uncertainty of robots. λ0<λ<1 is the discount factor. λ has relationship with pu as equation (9). The aim of equation (9) is that the discount factor can be adjusted dynamically based on and when is high; the long-term reward is not reliable anymore and the path planning based on MDP is concerned more on the optimization of total short-term reward.S1,…,Sα represents the state from the start time step when the current district needs to be planned to the time step which is α−1 time steps after the start time step. And for S1, it contains states for all the robots in the district, which means that S1=s1r1,s1r2,…,s1rm. The state s1r1 indicates the state for the robot r1after the action π1,1r1 is taken. R⋅ defines the reward for each state based on the distribution of the robots in the corresponding district.(9)λ=1−pu.Therefore, according to the Markov decision process in [8], the optimization objective in equation (8) implies the expectation of the rewards for the state from the start time step to after α−1 time steps. And it can be rewritten as follows:(10)minπxERS1+λ⋅RS2+⋯+λα−1RSα=minπx∑q=1αλq−1∑Sq′PSq−1πx,q−1Sq′⋅RSq′,where Sq′ indicates the possible state after the action πx,q−1 is taken at the state Sq−1 since there exists the motion uncertainty of robots, RSq′ represents the reward of the robots in this district at the state Sq′, and PSq−1πx,q−1Sq′ is the probability that the state of robots in this district would turn to Sq′ when the action πx,q−1 is taken at the state Sq−1. The probability can be calculated as follows:(11)PSq−1πx,q−1Sq′=∏j=1mPsq−1rj,πx,q−1rjsqrj′,where sq−1rj′is the state of the robot rj after the action πx,q−1,rj is taken at state sq−1rj and Psq−1rj,πx,q−1,rjsqrj′ is the probability that the state of the robot rjwill change into sqrj′ after the action πx,q−1rjis taken at state sq−1rj. And this probability is dependent on the motion uncertainty and can be represented as follows:(12)Psq−1rj,πx,q−1rjsqrj′=pu,πx,q−1rj=move andrjstays still,1−pu,πx,q−1rj=move andrjmoves,1,πx,q−1rj=waiting.For the rewardRSq′ of the policy πx, it is related to the crowdedness of the district and the distances among robots in the district, which can be defined as follows:(13)RSq′=Cdx′⋅ε⋅∑h=1mvaluesqrh′,where Cdx′ is the normalized crowdedness of the district being planned when the policy πx is taken. And the crowdedness Cdx of the district being planned can be calculated as follows:(14)Cdx=m∑i=1m∑j=1mdistanceij,where distanceij is the distance between ri and rj at state Sq′. Then, the crowdedness Cdx needs to be normalized by using the principle of the softmax according to the following formula:(15)Cdx′=eCdx∑h=1keCdh.In equation (13), valuesqrh′indicates the value of the robot rh at state sqrh′. It is defined as the distance between the robot rhand the next district which rhheads to since the passing order has been planned in the topology level. The rule for calculating this distance is based on the Manhattan distance.To avoid the conflicts among multiple robots, the factorε is introduced. And ε in equation (13) is defined as follows:(16)ε=∞,pit=pjt+1,pjt=pit+1orpit=pjt,∃i,j∈1,m,e1/1−pu,pit=pjt,∃i,j∈1,m,1,otherwise.It can be seen from equation (16) that when there would be a conflict between robot ri and rj, ε is set at ∞. When a robot would move to the location at next time step where another robot is staying at for now, then ε is set at e1/1−pu. Otherwise, ε is set at 1.According to equations (13) to (16), the first aim of this design is to prevent that each district would be too crowded since this would cause the coordination space to be too small, which means that once there happens motion uncertainty in some robots, it would be hard for the algorithm to find the solution and the system efficiency would be decreased. Another aim for the set of ε is to prevent the conflicts and make sure that one robot would not choose the action which might cause it to be next to another one. And this is beneficial for the motion coordination of robots since the motion uncertainty of one would cause it to still stay at the previous location, which increases the probability of the conflicts. Therefore, this setting can decrease the probability of the situation like that.To solve the optimization problem in equation (8), the policy set π needs to be generated in the district and the policy with the minimum objective would be selected as the solution. And the principle for the candidate action set of the policies is that robots would not be far away from the next district after taking this action. There are two benefits for this principle. Firstly, the farther the distance from the next district, the larger the objective function in equation (8). Secondly, due to the introduction of this principle, the candidate actions will decrease, which can in turn lower the computational complexity. According to above principle, there are two kinds of candidate action sets when robot is in different locations of each district:(1)
When robot is in any location of district and not coming into the boundary, the candidate action set of robot is denoted asleft,up,still. The action left means that robot drives to the other side of two-way street at next time step, up means that robot follows current side of street, and still means that robot stays still at next time step. For example, in Figure 5, robots r1 and r3 are in the district and not coming into boundary at next time step. Therefore, their candidate action set includes left, up, and still as orange arrows show. And as shown in Figure 5, the boundary between two districts is defined as an area composed of four grids, and one half of the boundary belongs to one district and the other belongs to the neighboring district. Considering that path planning in the route level is conducted separately in each district, robot should only follow one direction (drive to the right) in the process of traversing the boundary to the neighboring district, for avoiding the conflict.(2)
When robot is in the boundary of district or coming into the boundary, the candidate action set of robot is denoted asup,still or left,still according to the rule that robot should only drive to the right in the process of traversing the boundary. For example, robot r2in Figure 5 is coming into the boundary at next time step and its candidate action set isleft,still which means that it needs to drive to the right side of two-way street firstly. r4 and r5 are in the boundary and their candidate action set is up,still.Figure 5
Illustration of the candidate actions. Green squares represent the robots in the district. The blue arrow indicates the next district where these robots are heading to. Therefore, the orange arrows for each robot represent the candidate actions they can take. Besides the actions represented by these orange arrows, staying still is also included in the candidate actions for each time step.According to the above candidate action set, the generation process of robotri and policy πqri is shown in Figure 6. And it can obtain multiple policies for each robot. Then, π can be generated with the combination of multiple policies of each robot in the district.Figure 6
Generation of policyπqri for riin each district. Assuming that the candidate actions in each district include left, up, and still, the generation of policy πqri for robot ri is shown in this figure. In this tree structure, different level of branches corresponds to different steps. Each circle represents the action for ri in the district at the corresponding time step. For example, the action left means that all the robots in the district should go left. A complete policy πqri from the first time step to the αth time step is illustrated in the blue box. Multiple policies for each robot can be generated from this tree structure.


## 4. Discussion and Analysis

### 4.1. Conflict Avoidance

In this paper, we propose a bi-level path planning algorithm. The architecture of the proposed method is divided into topology level and route level. And in the topology level, the passing order of districts is planned for robots. There is no conflict in topology level because the districts are allowed to traverse for multiple robots at the same time. Therefore, the conflict avoidance of the proposed algorithm in the paper mainly depends on whether or not it could plan conflict-free path in route level. And we try to prove it in the following proposition.Proposition 1.
The path planning in route level could plan conflict-free paths for all robots in the district.Proof.
In route level of the proposed method, the path for robots in the district is planned for many times andα steps with MDP at each time. Therefore, we try to use mathematical induction to prove the conflict-free path can always be planned by the proposed method in the district at each time. And we firstly describe some symbols of some variables to facilitate the expression of proof.
For any districtds∈D with m robots, the route map can be denoted as GpVp,Ep. And robots in the district can be denoted as Rp=r1,…,rm. nrerepresents the nreth path planning in each district. In the nreth planning: the beginning of time step is denoted as tnre, and the planning step is α. The solution with MDP is πnreand planned path of πnre for any robot ri is pitnre,pitnre+1,…,pitnre+α. The proof detail is as follows:(a)
Firstly, whennre=1, prove that the path planning is conflict-free:At the first time of path planning, robots are leaving the boundary in each district. These robots of each district are heading to different directions and located in the different sides of two-way street in the map. Therefore, there always exists one solutionπ1 that makes sure that all robots go on with their current side of two-way street in the α time steps and no robots from the opposite direction run in the same side of two-way street. And the planned paths for robots meet equation (17). It can be concluded that reward of the above solution π1 according to equation (16) will not be ∞ and π1could be feasible solution of the optimization objective, which could prove that the proposed method could plan conflict-free path for all robots when nre=1.(17)pit≠pjt+1andpjt≠pit+1,∀i,j∈Rp,tnre≤t≤tnre+1,pit≠pjt,∀i,j∈Rp,tnre≤t≤tnre+1.(b)
Secondly, whennre=k1, suppose that the planned path in the k1th planning is conflict-free and prove the planned path in the k1+1th planning is also conflict-free:Considering that the path planning in thek1th planning is conflict-free, the ending location pitk1+α of planned path for any robot ri∈Rp satisfies equation (17), which means that there does not exist the situation of direction conflict and location. On the other hand, the ending location pitk1+α of the k1th path planning is equivalent to the beginning location pitk1+1 of the k1+1th path planning. If there exists v1i,v2i∈Rp where v1i represents the next location from with the action of going up and v2i represents the next location from pitk1+1 with the action of going left, equation (18) can be concluded. It means that each robot rican go left or up without conflict with other robots. Therefore, it can construct a candidate solution πk1 that makes robots heading to the same direction run in the same side of two-way streets in the α steps of path. And when robots heading to the different direction run in the same side of two-way street, the planned path will not have conflict with other robots. Therefore, it can be proved that the planned path in the k1+1th planning is conflict-free.(18)phtnre≠v2i,∀h∈Rpandh≠iifpjtnre=v1i,∃j∈Rpandj≠i,phtnre≠v1i,∀h∈Rpandh≠iifpjtnre=v2i,∃j∈Rpandj≠i,and when robots go to the boundary in the last time planning, considering that the robots always follow one direction when they are coming into boundary or traversing on the boundary, the district-distributed way of the proposed method in route level will not have conflict between robots from different districts.
To summarize, there always would be conflict-free solution for path planning in the route level with the proposed method.


### 4.2. Optimization

For a multirobot system with each robot subject to uncertainty, directly computing the optimal solution is very difficult and it must recompute many times during the running of robots. Therefore, it is hard to directly prove the optimization or near optimization. In this section, we will discuss the optimization of the proposed method. Considering that the architecture of proposed method in this paper is divided into two levels and paths are planned separately, we discuss the optimization of proposed method from following two parts:(1)
For path planning in topology level: in the topology level, the optimal orders of districts for all robots can be obtained and the reason is as follows. Firstly, K shortest path planning algorithm [24] is adopted to obtain the multiple orders of passing districts for each robot. And combinations of orders of passing districts among robots can be generated. When K is large enough, all possible combinations of orders of passing districts can be obtained. Secondly, equation (3) is put up with in the topology level, which aims to minimize the corresponding distance between the center of the two districts and total crowdedness between districts among all the robots assigned with new tasks. And above optimization objective can be calculated based on all possible combinations of orders of passing of districts for robots.(2)
For path planning in route level: the proposed method takes path planning for all robots in each district for many times and plans paths ofα steps for the robots at each time until all robots leave the district. So, it is difficult to directly prove the optimization. In this part, we compare above method with existing method in [8] to discuss the optimization of the proposed method. Similar to the proposed method, the existing method also takes path planning in each district for many times. The difference is that at each time planning, the existing method plans the path of static αglobal steps with the MDP method and executes the first α steps of them. The optimization objective is shown in equation (19). αglobal can make sure that the destination of planned path is the location where robots are leaving from the current district. In the existing method in [8], α is also intuitively related to pu following equation (7). The existing method in [8] has proved the optimization.(19)minπx∑q=1αλq−1∑Sq′PSq−1πx,q−1Sq′⋅RSq′+∑q=1αglobal−αλα+q−1∑Sq′PSq−1πx,q−1Sq′⋅RSq′.For discussing the optimization of proposed method, we try to discuss the similarity of firstα steps for solution at each time planning compared with existing method in [8] in the following two situations:(a)
When the motion certainty probabilitypu is low pu≤0.2:In this situation,α is close to αglobal according to (7). And it is obvious that the first α steps of solution from the proposed method and existing method in [6] are close to the same. For example, Figure 7 shows a situation where three robots are located at different states, and Table 1 shows the different solutions of path planning for the situation in Figure 7 using above two methods when pu is different (αglobal is set to 10 for the existing method in [8]). As shown in Table 1, when pu is 0.1 or 0.2, first α steps of the proposed method are the same as the existing method in [6].(b)
When the motion certainty probabilitypuis high (pu>2):In this situation, consider that both the proposed method and existing method in [8] choose the MDP method where the discount factor λ exists in optimization objective and has relationship with pu according to equation (9). Therefore, when pu is high, λα in (19) will be lower. This means the last αglobal−α steps of solution have little influence on the computation for equation (19). And in this situation, the first α steps of the proposed method are close to the existing method in [6]. For example, in Table 1, when α is close to αglobal according to (19), it is obvious that the solution from the proposed method is close to first α steps of solution from the existing method in [8]. Figure 7 shows a situation where three robots are located at different states, and Table 1 shows the different solutions of above two methods when pu is different. As shown in Table 1, in the example of Figure 7, when pu is 0.4 or 0.5, although the solution of above two methods is not exactly the same, the solution from the proposed method is close to the first α steps of solution from the existing method in [8].Figure 7
Illustration of the example of trying to plan path for three robots in the district. The original location of robots is shown in the figure. And the blue arrow indicates the next district where these three robots are heading to.v1∼v11 represent different locations in the district and v6 is the destination where robots are leaving from the current district.Table 1
The solution of path planning for three robots which are in different locations as shown in Figure7.
The solution of the proposed methodThe firstα−step solution of the existing method in [8]pu=0.1α=10λ=0.9r1:v0,up,v2,up,v5,left,v4,upv7,left,v6,leftr2:v12,up,v9,up,v7,left,v6,leftr3:v13,up,v11,up,v10,left,v8,stillv8,left,v7,left,v6,leftr1:v0,up,v2,up,v5,left,v4,upv7,left,v6,leftr2:v12,up,v9,up,v7,left,v6,leftr3:v13,up,v11,up,v10,left,v8,stillv8,left,v7,left,v6,leftpu=0.2α=5λ=0.8r1:v0,up,v2,up,v5,left,v4,upv7,leftr2:v12,up,v9,up,v7,left,v6,upr3:v13,up,v11,up,v10,left,v8,stillv8,leftr1:v0,up,v2,up,v5,left,v4,upv7,leftr2:v12,up,v9,up,v7,left,v6,upr3:v13,up,v11,up,v10,left,v8,stillv8,leftpu=0.3α=3λ=0.7r1:v0,up,v2,up,v5,upr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftr1:v0,up,v2,up,v5,upr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftpu=0.4α=3λ=0.6r1:v0,up,v2,up,v5,upr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftr1:v0,up,v2,up,v5,leftr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftpu=0.5α=2λ=0.5r1:v0,up,v2,upr2:v12,up,v9,upr3:v13,up,v11,upr1:v0,up,v2,upr2:v12,up,v9,upr3:v13,up,v11,up


## 4.1. Conflict Avoidance

In this paper, we propose a bi-level path planning algorithm. The architecture of the proposed method is divided into topology level and route level. And in the topology level, the passing order of districts is planned for robots. There is no conflict in topology level because the districts are allowed to traverse for multiple robots at the same time. Therefore, the conflict avoidance of the proposed algorithm in the paper mainly depends on whether or not it could plan conflict-free path in route level. And we try to prove it in the following proposition.Proposition 1.
The path planning in route level could plan conflict-free paths for all robots in the district.Proof.
In route level of the proposed method, the path for robots in the district is planned for many times andα steps with MDP at each time. Therefore, we try to use mathematical induction to prove the conflict-free path can always be planned by the proposed method in the district at each time. And we firstly describe some symbols of some variables to facilitate the expression of proof.
For any districtds∈D with m robots, the route map can be denoted as GpVp,Ep. And robots in the district can be denoted as Rp=r1,…,rm. nrerepresents the nreth path planning in each district. In the nreth planning: the beginning of time step is denoted as tnre, and the planning step is α. The solution with MDP is πnreand planned path of πnre for any robot ri is pitnre,pitnre+1,…,pitnre+α. The proof detail is as follows:(a)
Firstly, whennre=1, prove that the path planning is conflict-free:At the first time of path planning, robots are leaving the boundary in each district. These robots of each district are heading to different directions and located in the different sides of two-way street in the map. Therefore, there always exists one solutionπ1 that makes sure that all robots go on with their current side of two-way street in the α time steps and no robots from the opposite direction run in the same side of two-way street. And the planned paths for robots meet equation (17). It can be concluded that reward of the above solution π1 according to equation (16) will not be ∞ and π1could be feasible solution of the optimization objective, which could prove that the proposed method could plan conflict-free path for all robots when nre=1.(17)pit≠pjt+1andpjt≠pit+1,∀i,j∈Rp,tnre≤t≤tnre+1,pit≠pjt,∀i,j∈Rp,tnre≤t≤tnre+1.(b)
Secondly, whennre=k1, suppose that the planned path in the k1th planning is conflict-free and prove the planned path in the k1+1th planning is also conflict-free:Considering that the path planning in thek1th planning is conflict-free, the ending location pitk1+α of planned path for any robot ri∈Rp satisfies equation (17), which means that there does not exist the situation of direction conflict and location. On the other hand, the ending location pitk1+α of the k1th path planning is equivalent to the beginning location pitk1+1 of the k1+1th path planning. If there exists v1i,v2i∈Rp where v1i represents the next location from with the action of going up and v2i represents the next location from pitk1+1 with the action of going left, equation (18) can be concluded. It means that each robot rican go left or up without conflict with other robots. Therefore, it can construct a candidate solution πk1 that makes robots heading to the same direction run in the same side of two-way streets in the α steps of path. And when robots heading to the different direction run in the same side of two-way street, the planned path will not have conflict with other robots. Therefore, it can be proved that the planned path in the k1+1th planning is conflict-free.(18)phtnre≠v2i,∀h∈Rpandh≠iifpjtnre=v1i,∃j∈Rpandj≠i,phtnre≠v1i,∀h∈Rpandh≠iifpjtnre=v2i,∃j∈Rpandj≠i,and when robots go to the boundary in the last time planning, considering that the robots always follow one direction when they are coming into boundary or traversing on the boundary, the district-distributed way of the proposed method in route level will not have conflict between robots from different districts.
To summarize, there always would be conflict-free solution for path planning in the route level with the proposed method.


## 4.2. Optimization

For a multirobot system with each robot subject to uncertainty, directly computing the optimal solution is very difficult and it must recompute many times during the running of robots. Therefore, it is hard to directly prove the optimization or near optimization. In this section, we will discuss the optimization of the proposed method. Considering that the architecture of proposed method in this paper is divided into two levels and paths are planned separately, we discuss the optimization of proposed method from following two parts:(1)
For path planning in topology level: in the topology level, the optimal orders of districts for all robots can be obtained and the reason is as follows. Firstly, K shortest path planning algorithm [24] is adopted to obtain the multiple orders of passing districts for each robot. And combinations of orders of passing districts among robots can be generated. When K is large enough, all possible combinations of orders of passing districts can be obtained. Secondly, equation (3) is put up with in the topology level, which aims to minimize the corresponding distance between the center of the two districts and total crowdedness between districts among all the robots assigned with new tasks. And above optimization objective can be calculated based on all possible combinations of orders of passing of districts for robots.(2)
For path planning in route level: the proposed method takes path planning for all robots in each district for many times and plans paths ofα steps for the robots at each time until all robots leave the district. So, it is difficult to directly prove the optimization. In this part, we compare above method with existing method in [8] to discuss the optimization of the proposed method. Similar to the proposed method, the existing method also takes path planning in each district for many times. The difference is that at each time planning, the existing method plans the path of static αglobal steps with the MDP method and executes the first α steps of them. The optimization objective is shown in equation (19). αglobal can make sure that the destination of planned path is the location where robots are leaving from the current district. In the existing method in [8], α is also intuitively related to pu following equation (7). The existing method in [8] has proved the optimization.(19)minπx∑q=1αλq−1∑Sq′PSq−1πx,q−1Sq′⋅RSq′+∑q=1αglobal−αλα+q−1∑Sq′PSq−1πx,q−1Sq′⋅RSq′.For discussing the optimization of proposed method, we try to discuss the similarity of firstα steps for solution at each time planning compared with existing method in [8] in the following two situations:(a)
When the motion certainty probabilitypu is low pu≤0.2:In this situation,α is close to αglobal according to (7). And it is obvious that the first α steps of solution from the proposed method and existing method in [6] are close to the same. For example, Figure 7 shows a situation where three robots are located at different states, and Table 1 shows the different solutions of path planning for the situation in Figure 7 using above two methods when pu is different (αglobal is set to 10 for the existing method in [8]). As shown in Table 1, when pu is 0.1 or 0.2, first α steps of the proposed method are the same as the existing method in [6].(b)
When the motion certainty probabilitypuis high (pu>2):In this situation, consider that both the proposed method and existing method in [8] choose the MDP method where the discount factor λ exists in optimization objective and has relationship with pu according to equation (9). Therefore, when pu is high, λα in (19) will be lower. This means the last αglobal−α steps of solution have little influence on the computation for equation (19). And in this situation, the first α steps of the proposed method are close to the existing method in [6]. For example, in Table 1, when α is close to αglobal according to (19), it is obvious that the solution from the proposed method is close to first α steps of solution from the existing method in [8]. Figure 7 shows a situation where three robots are located at different states, and Table 1 shows the different solutions of above two methods when pu is different. As shown in Table 1, in the example of Figure 7, when pu is 0.4 or 0.5, although the solution of above two methods is not exactly the same, the solution from the proposed method is close to the first α steps of solution from the existing method in [8].Figure 7
Illustration of the example of trying to plan path for three robots in the district. The original location of robots is shown in the figure. And the blue arrow indicates the next district where these three robots are heading to.v1∼v11 represent different locations in the district and v6 is the destination where robots are leaving from the current district.Table 1
The solution of path planning for three robots which are in different locations as shown in Figure7.
The solution of the proposed methodThe firstα−step solution of the existing method in [8]pu=0.1α=10λ=0.9r1:v0,up,v2,up,v5,left,v4,upv7,left,v6,leftr2:v12,up,v9,up,v7,left,v6,leftr3:v13,up,v11,up,v10,left,v8,stillv8,left,v7,left,v6,leftr1:v0,up,v2,up,v5,left,v4,upv7,left,v6,leftr2:v12,up,v9,up,v7,left,v6,leftr3:v13,up,v11,up,v10,left,v8,stillv8,left,v7,left,v6,leftpu=0.2α=5λ=0.8r1:v0,up,v2,up,v5,left,v4,upv7,leftr2:v12,up,v9,up,v7,left,v6,upr3:v13,up,v11,up,v10,left,v8,stillv8,leftr1:v0,up,v2,up,v5,left,v4,upv7,leftr2:v12,up,v9,up,v7,left,v6,upr3:v13,up,v11,up,v10,left,v8,stillv8,leftpu=0.3α=3λ=0.7r1:v0,up,v2,up,v5,upr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftr1:v0,up,v2,up,v5,upr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftpu=0.4α=3λ=0.6r1:v0,up,v2,up,v5,upr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftr1:v0,up,v2,up,v5,leftr2:v12,up,v9,up,v7,leftr3:v13,up,v11,up,v10,leftpu=0.5α=2λ=0.5r1:v0,up,v2,upr2:v12,up,v9,upr3:v13,up,v11,upr1:v0,up,v2,upr2:v12,up,v9,upr3:v13,up,v11,up


## 5. Simulation Result and Analysis

In order to evaluate the proposed method, simulations were implemented under the Matlab environment with Intel 3.60 GHz Core i7-4790U CPU and 8G RAM. The map of warehouse is shown in Figure1.The motion uncertaintypu ranges from 0 to 1/2 to simulate the situations under different motion modes. The number of robots is set at 10, 20, 30, 40, and 50. The number of tasks is 1000.The makespan, number of planning steps, and computation time are contrasted with those of the existing method in [21, 22] under a different number of robots. And for the existing method in [21, 22], the probability corresponding to the disturbance intensity whether during the following second the robot will be prevented from moving is set to pu and the lower bound which represents the best possible travel time is set to the shortest distance from start to end for each robot. The average of 20 repeated simulations is shown in Tables 2–4. The tasks have been generated before simulations.Table 2
Simulation results in a warehouse environment of 50 × 34 4-neighbor grids: makespan (standard deviation).
puRobot numberExisting method in [21, 22]Proposed methodDecrease percentage (%)0.11019794.5272.419022.4187.13.92013582.6285.612210.8193.510.13011394.5294.29297.9203.818.4409019.3314.87044.1209.421.9507768.4332.95818.5216.725.10.21022956.7280.422061.4193.53.92015636.6290.114041.7199.410.23013103.2301.510666.1204.518.6409019.3317.68702.4210.422.2509587.2334.97123.3220.725.70.31025895.6287.524859.8197.24.02016879.3298.415140.7200.110.43014374.5309.711657.7208.718.94012287.3338.69510.4218.522.65010625.4344.87820.4229.526.40.41028190.4294.127034.6196.54.12018724.5302.616777.2200.410.43016394.6315.313246.8206.819.24014274.8334.710991.6218.523.05012463.4354.79085.8230.527.10.51031852.4300.630514.6199.24.22019952.6310.917857.6201.510.53017365.4325.613979.1208.719.54015284.7344.811692.8221.823.55013561.2356.39777.6231.527.9Table 3
Simulation results in a warehouse environment of 50 × 34 4-neighbor grids: number of planning (standard deviation).
puRobot numberExisting method in [21, 22]Planning percentage (%)Proposed methodPlanning percentage (%)0.11019794.3272.41002587.0102.413.62013582.6285.61001721.7104.114.13011394.5294.21001348.2107.314.5409019.3314.81001056.6109.915.0507768.4332.9100907.7112.615.60.21022956.7280.41005405.0104.324.52015636.6290.11003510.4106.425.03013103.2301.51002815.9109.626.44011185.6317.61002427.9112.127.9509587.2334.91002108.5115.929.60.31025895.6287.51008825.2106.535.52016879.3298.41005708.1108.537.73014374.5309.71004558.2110.839.14012287.3328.61003937.3113.441.45010625.6344.81003456.6116.544.20.41028190.4294.110012544.1108.646.42018724.5302.61008170.5110.748.73016394.6315.31006861.8113.651.84014274.8334.71006034.4116.854.95012463.4354.71005324.3119.958.60.51031852.4300.610017728.9110.358.12019952.6310.910010750.3112.960.23017365.4325.61008764.9115.662.74015284.7344.81007647.1118.465.45013561.2356.31006736.8122.468.9Table 4
Simulation results in a warehouse environment of 50 × 34 4-neighbor grids: computation time (standard deviation).
puRobot numberExisting method in [21, 22]Proposed method0.1100.0015s0.0001540.0083s0.000356200.0021s0.0002020.1287s0.000391300.0026s0.0002830.2176s0.000514400.0033s0.0003620.3025s0.000698500.0041s0.0004510.3786s0.0008240.2100.0016s0.0001560.0041s0.000322200.0022s0.0002090.095s0.000371300.0027s0.0002910.1756s0.000432400.0034s0.0003690.2476s0.000599500.0042s0.0004680.3012s0.0007740.3100.0017s0.0001620.0029s0.000316200.0023s0.0002180.0523s0.000354300.0028s0.0004010.1026s0.000407400.0035s0.0003790.1424s0.000558500.0043s0.0004810.1672s0.0007320.4100.0019s0.0001710.0026s0.000301200.0025s0.0002260.0278s0.000334300.0030s0.0003110.0478s0.000386400.0037s0.0003970.0623s0.000514500.0045s0.0004920.0745s0.0007020.5100.0021s0.0001780.0021s0.000289200.0027s0.0002350.0312s0.000318300.0032s0.0003230.0558s0.000371400.0039s0.0003980.0797s0.000501500.0047s0.0005060.0976s0.000688


### 5.1. Makespan

Table2 shows the simulation results of average makespan of the existing method and the proposed method. Standard deviations are also shown in Table 2. Decrease percentage of makespan (the last column in Table 2) is calculated according to the following formula:(20)decrease percentage=makespan1−makespan2makespan1×100%,where makespan1 and makespan2 represent the makespan of the existing method in [21, 22] and the proposed method, respectively. Graphical summaries of simulation results of makespan are shown in Figure 8.Figure 8
Decrease percentage of average makespan in a warehouse environment of 50 × 34 4-neighbor grids.


#### 5.1.1. Average Time

Compared to the existing method in [21, 22], it can be seen that the proposed method in this paper effectively decreases the makespan. When the number of robots is relatively small, the decrease percentage is low because the robots are located sparsely in the warehouse environment and the motion uncertainties of robots would not influence too many other normal robots, which means that most of the normal robots could still follow the previous planned paths. As the number of robots increases, the decrease percentage of makespan is greatly improved because the number of influenced robots would increase further, which means that the proposed method would replan to find more optimal paths based on the probability pu for more robots to avoid robots which are in the motion uncertainty. Under this circumstance, by combining MDP, redundant waiting of the method in [21, 22] is not necessary. When the value of pu increases, the proposed method could still improve the system efficiency and has a smaller variance compared to the existing method in [21, 22]. The percentage of decrease can be up to 27.9% when the number of robots is 50 and pu is 0.5.


#### 5.1.2. Variance

Compared to the existing method in [21, 22], the variance of the proposed method is smaller. Furthermore, with the increase of the number of robots, the variance in the results of the proposed method increases less than that of the existing method in [21, 22]. The comparison between the results of the proposed method and the existing method in [21, 22] means that the proposed method has less randomness and is more stable. The reason is that the method in this paper would optimize paths of all the robots depending on the probability pu.In the proposed method, unlike the existing method in [21, 22] that would solve the problem of motion uncertainty through a waiting strategy, the proposed method to coordinate paths of robots under motion uncertainties adopts the strategy of replanning new paths to choose the probabilistic optimal solution. Therefore, when motion uncertainty happens, all the normal robots would follow the paths that could reach the destination as soon as possible from the global view instead of adopting the strategy of going forward and waiting. Therefore, compared to the existing methods in [21, 22], the decrease of the efficiency induced by the motion uncertainties could be lowered. In this way, the efficiency of the whole system could be improved when motion uncertainties happen by adopting the proposed method. Compared to the existing methods, the efficiency could be increased, which is beneficial for the overall running of the system.


### 5.2. Number of Planning Steps

Percentage of planning (listed in Table2) is calculated according to the following formula:(21)planning percentage=Nremakespan×100%,where Nre and makespan represent the number of planning and the makespan of corresponding method, respectively. The makespan of each method is given in Table 2. Graphical summaries of simulation results of planning percentage are shown in Figure 9. It can be seen from Table 3 that the percentage of planning of the proposed method is much more smaller than that of the method in [21, 22]. It has to plan at each time step since the method of [21, 22] is a kind of reactive method. Therefore, the percentage is 100%. On the contrary, the proposed method in the route level would only plan when robots have finished the path of planning step α or when the motion uncertainty happens. Therefore, the percentage of planning would decrease, compared with the existing methods in [21, 22].Figure 9
Planning percentage in a warehouse environment of 50 × 34 4-neighbor grids. (a)pu = 0.1. (b) pu = 0.2. (c) pu = 0.3. (d) pu = 0.4. (e) pu = 0.5.
(a)(b)(c)(d)(e)Another conclusion we can get from Table3 is that when the probability pu is small, the percentage reduction compared with the method in [21, 22] is larger and when the probability pu increases, the percentage reduction becomes smaller. The reason for this is that when pu is small, there would not be too many time steps when motion uncertainties happen, which means that there would not be too many time steps when planning is needed. Therefore, the percentage of planning is relatively small, which causes the percentage reduction to be larger. When pu becomes larger, the time steps when motion uncertainties happen become more, which causes more planning. Therefore, the percentage of planning increases, which causes the percentage reduction to be smaller. However, it is noteworthy that even if the percentage of planning becomes larger when puincreases, the percentage of the proposed method is still fewer than that of the method in [21, 22].


### 5.3. Computation Time

Table4 shows the average of the computation time under different pu. The statistics in the table demonstrate that the proposed method could satisfy the requirement of the real-time performance. The reason that the proposed method is higher than that of the existing method in [21, 22] is that the existing method only plans for the next time step while the proposed method needs to plan the whole path to the destination. Therefore, the computation of the proposed method would be larger than that of the existing method. However, as shown in Figure 10, the trend of computation time is more close to linear relationship for the proposed method and increasing speed will reduce slowly with robot number increasing. On the other hand, even if the computation time is somewhat larger, the proposed method could still plan conflict-free paths for all the robots in real time, which can be applied for the lifelong scheduling in warehouses.Figure 10
Computation time in a warehouse environment of 50 × 34 4-neighbor grids. (a)pu = 0.1. (b) pu = 0.2. (c) pu = 0.3. (d) pu = 0.4. (e) pu = 0.5.
(a)(b)(c)(d)(e)The real-time performance of the bi-level probabilistic path planning algorithm for multiple robots with motion uncertainty is guaranteed by the bi-level planning architecture and the distributed planning in each district. The topology level would plan the passing orders of districts for all robots assigned with tasks and the path planning is conducted in each district for route level. The planning in each district only needs to take the actions of the current district being planned into consideration. Therefore, the number of planning depends onpu, which means that the computation is adjusted by pu.


## 5.1. Makespan

Table2 shows the simulation results of average makespan of the existing method and the proposed method. Standard deviations are also shown in Table 2. Decrease percentage of makespan (the last column in Table 2) is calculated according to the following formula:(20)decrease percentage=makespan1−makespan2makespan1×100%,where makespan1 and makespan2 represent the makespan of the existing method in [21, 22] and the proposed method, respectively. Graphical summaries of simulation results of makespan are shown in Figure 8.Figure 8
Decrease percentage of average makespan in a warehouse environment of 50 × 34 4-neighbor grids.


### 5.1.1. Average Time

Compared to the existing method in [21, 22], it can be seen that the proposed method in this paper effectively decreases the makespan. When the number of robots is relatively small, the decrease percentage is low because the robots are located sparsely in the warehouse environment and the motion uncertainties of robots would not influence too many other normal robots, which means that most of the normal robots could still follow the previous planned paths. As the number of robots increases, the decrease percentage of makespan is greatly improved because the number of influenced robots would increase further, which means that the proposed method would replan to find more optimal paths based on the probability pu for more robots to avoid robots which are in the motion uncertainty. Under this circumstance, by combining MDP, redundant waiting of the method in [21, 22] is not necessary. When the value of pu increases, the proposed method could still improve the system efficiency and has a smaller variance compared to the existing method in [21, 22]. The percentage of decrease can be up to 27.9% when the number of robots is 50 and pu is 0.5.


### 5.1.2. Variance

Compared to the existing method in [21, 22], the variance of the proposed method is smaller. Furthermore, with the increase of the number of robots, the variance in the results of the proposed method increases less than that of the existing method in [21, 22]. The comparison between the results of the proposed method and the existing method in [21, 22] means that the proposed method has less randomness and is more stable. The reason is that the method in this paper would optimize paths of all the robots depending on the probability pu.In the proposed method, unlike the existing method in [21, 22] that would solve the problem of motion uncertainty through a waiting strategy, the proposed method to coordinate paths of robots under motion uncertainties adopts the strategy of replanning new paths to choose the probabilistic optimal solution. Therefore, when motion uncertainty happens, all the normal robots would follow the paths that could reach the destination as soon as possible from the global view instead of adopting the strategy of going forward and waiting. Therefore, compared to the existing methods in [21, 22], the decrease of the efficiency induced by the motion uncertainties could be lowered. In this way, the efficiency of the whole system could be improved when motion uncertainties happen by adopting the proposed method. Compared to the existing methods, the efficiency could be increased, which is beneficial for the overall running of the system.


## 5.1.1. Average Time

Compared to the existing method in [21, 22], it can be seen that the proposed method in this paper effectively decreases the makespan. When the number of robots is relatively small, the decrease percentage is low because the robots are located sparsely in the warehouse environment and the motion uncertainties of robots would not influence too many other normal robots, which means that most of the normal robots could still follow the previous planned paths. As the number of robots increases, the decrease percentage of makespan is greatly improved because the number of influenced robots would increase further, which means that the proposed method would replan to find more optimal paths based on the probability pu for more robots to avoid robots which are in the motion uncertainty. Under this circumstance, by combining MDP, redundant waiting of the method in [21, 22] is not necessary. When the value of pu increases, the proposed method could still improve the system efficiency and has a smaller variance compared to the existing method in [21, 22]. The percentage of decrease can be up to 27.9% when the number of robots is 50 and pu is 0.5.


## 5.1.2. Variance

Compared to the existing method in [21, 22], the variance of the proposed method is smaller. Furthermore, with the increase of the number of robots, the variance in the results of the proposed method increases less than that of the existing method in [21, 22]. The comparison between the results of the proposed method and the existing method in [21, 22] means that the proposed method has less randomness and is more stable. The reason is that the method in this paper would optimize paths of all the robots depending on the probability pu.In the proposed method, unlike the existing method in [21, 22] that would solve the problem of motion uncertainty through a waiting strategy, the proposed method to coordinate paths of robots under motion uncertainties adopts the strategy of replanning new paths to choose the probabilistic optimal solution. Therefore, when motion uncertainty happens, all the normal robots would follow the paths that could reach the destination as soon as possible from the global view instead of adopting the strategy of going forward and waiting. Therefore, compared to the existing methods in [21, 22], the decrease of the efficiency induced by the motion uncertainties could be lowered. In this way, the efficiency of the whole system could be improved when motion uncertainties happen by adopting the proposed method. Compared to the existing methods, the efficiency could be increased, which is beneficial for the overall running of the system.


## 5.2. Number of Planning Steps

Percentage of planning (listed in Table2) is calculated according to the following formula:(21)planning percentage=Nremakespan×100%,where Nre and makespan represent the number of planning and the makespan of corresponding method, respectively. The makespan of each method is given in Table 2. Graphical summaries of simulation results of planning percentage are shown in Figure 9. It can be seen from Table 3 that the percentage of planning of the proposed method is much more smaller than that of the method in [21, 22]. It has to plan at each time step since the method of [21, 22] is a kind of reactive method. Therefore, the percentage is 100%. On the contrary, the proposed method in the route level would only plan when robots have finished the path of planning step α or when the motion uncertainty happens. Therefore, the percentage of planning would decrease, compared with the existing methods in [21, 22].Figure 9
Planning percentage in a warehouse environment of 50 × 34 4-neighbor grids. (a)pu = 0.1. (b) pu = 0.2. (c) pu = 0.3. (d) pu = 0.4. (e) pu = 0.5.
(a)(b)(c)(d)(e)Another conclusion we can get from Table3 is that when the probability pu is small, the percentage reduction compared with the method in [21, 22] is larger and when the probability pu increases, the percentage reduction becomes smaller. The reason for this is that when pu is small, there would not be too many time steps when motion uncertainties happen, which means that there would not be too many time steps when planning is needed. Therefore, the percentage of planning is relatively small, which causes the percentage reduction to be larger. When pu becomes larger, the time steps when motion uncertainties happen become more, which causes more planning. Therefore, the percentage of planning increases, which causes the percentage reduction to be smaller. However, it is noteworthy that even if the percentage of planning becomes larger when puincreases, the percentage of the proposed method is still fewer than that of the method in [21, 22].


## 5.3. Computation Time

Table4 shows the average of the computation time under different pu. The statistics in the table demonstrate that the proposed method could satisfy the requirement of the real-time performance. The reason that the proposed method is higher than that of the existing method in [21, 22] is that the existing method only plans for the next time step while the proposed method needs to plan the whole path to the destination. Therefore, the computation of the proposed method would be larger than that of the existing method. However, as shown in Figure 10, the trend of computation time is more close to linear relationship for the proposed method and increasing speed will reduce slowly with robot number increasing. On the other hand, even if the computation time is somewhat larger, the proposed method could still plan conflict-free paths for all the robots in real time, which can be applied for the lifelong scheduling in warehouses.Figure 10
Computation time in a warehouse environment of 50 × 34 4-neighbor grids. (a)pu = 0.1. (b) pu = 0.2. (c) pu = 0.3. (d) pu = 0.4. (e) pu = 0.5.
(a)(b)(c)(d)(e)The real-time performance of the bi-level probabilistic path planning algorithm for multiple robots with motion uncertainty is guaranteed by the bi-level planning architecture and the distributed planning in each district. The topology level would plan the passing orders of districts for all robots assigned with tasks and the path planning is conducted in each district for route level. The planning in each district only needs to take the actions of the current district being planned into consideration. Therefore, the number of planning depends onpu, which means that the computation is adjusted by pu.


## 6. Conclusion

For improving the system efficiency when there are motion uncertainties among robots in the warehouse environment, this paper proposes a bi-level probabilistic path planning algorithm. In the proposed algorithm, the map is partitioned into multiple interconnected districts and the architecture of proposed algorithm is composed of topology level and route level generating from above map: in the topology level, the district crowdedness is taken into consideration to optimize the passing orders of districts and achieve the district equilibrium; in the route level, path planning is conducted in each district with a MDP-based probabilistic method, where the planning step is related to the probability of the motion uncertainty to improve the efficiency and guaranteeing the real-time performance. Furthermore, the conflict avoidance of designed algorithm is proved and optimization is discussed. The simulation results validate that the designed algorithm can improve the system efficiency and the percentage can be up to 27.9% in the designed warehouses environment of this paper, which also accords with optimization discussion. It can be also seen that the number of planning has been reduced to a large extent.Our future work will involve the further optimization of paths by taking the interaction of flow among different districts into consideration since the actions are only obtained by the distributed planning in each district. Another aspect is that the real experiments are needed to verify the effectiveness of the proposed method.



---

*Source: 9207324-2020-06-05.xml*