## 1 Introduction  
Multi-robot systems play a vital role in next-generation factories, urban search and rescue, and package delivery, and they are anticipated to be applied in space exploration [1, 2]. One of the key ingredients to a multi-robot system is the motion planning module. The motion planning module should find a set sequence of valid configurations that move multiple robots from starting points to goal points safely [3]. However, in unstructured environments, the multi-robot motion planning algorithm being both optimal, complete, real-time, and flexible is still a valuable and challenging problem to be solved, especially in the presence of disturbances. As mentioned in previous publications, the multi-robot motion planning (MMP) problem is generally defined as a multi-objective optimization problem with different types of constraints [4, 5]. The methods for solving MMP problems are mainly classified into centralized and decentralized methods [3]. Centralized methods plan the motion trajectories of robots jointly [6]. In [7], [PERSON] et al. proposed a two-level algorithm called conflict-based search (CBS) to find optimal paths for multiple robots. The core of the algorithm is to maintain a constraint tree to resolve conflicts and plan the optimal path for each robot with constraints. In [8], [PERSON] et al. improved the CBS algorithm based on explicit estimation search. Their method utilizes online learning to acquire the unacceptable estimated cost of the solution and adopts explicit estimation search to resolve conflicts. Although these methods guarantee optimality and real-time performance, the paths are discrete nodes and the application is limited to holonomic robots. Since the above methods do not consider the robot kinematic and dynamical characteristics, the paths generally need to be further processed before they can be applied to non-holonomic robots. To deal with the MMP problem of non-holonomic robots, a two-stage pipeline, i.e., reference path planning and trajectory generation, is applied in a series of works [9, 10, 11]. In these works, collision-free paths are planned for all robots firstly. Then, the reference paths are treated as an initial guess of the MMP problem. The nonlinear optimization problem [11] or quadratic problem [9, 10] is formulated to obtain collision-free dynamically feasible trajectories. Compared with solving the problem without an initial guess, the two-stage pipeline can achieve better efficiency and quality. Decentralized methods typically decouple large-scale multi-robot motion planning problems into multiple smaller-scale multi-robot motion planning problems [6]. The most famous decentralized method is the velocity obstacle method [12, 13, 14, 15]. The main idea of the velocity obstacle method is to choose a suitable velocity in the feasible velocity set. As this method only considers local obstacle avoidance, the global performance is poor and may lead the robot into deadlock. Another well-known decentralized method is distributed model predictive control (DMPC) [4, 16, 17]. In [17], [PERSON] et al. proposed a distributed model predictive control (DMPC) method for multi-robot motion planning to generate trajectories in real-time. The DMPC method improves computational efficiency but reduces the quality of the solution. More importantly, multi-robot systems operating in the real world are often subject to various disturbances. Disturbances lead to inaccurate estimates of the states of multi-robot systems. To ensure robustness in multi-robot motion planning, disturbances need to be considered in the design of the algorithm. The disturbances are often modeled as chance constraints [18, 19, 20]. In [21], [PERSON] et al. defined a chance constraint called probabilistic safety barrier certificates (PrSBC) and solved an optimization problem with the PrSBC. In [22], [PERSON] et al. studied a probabilistic variant of the reciprocal velocity obstacle method called probabilistic reciprocal velocity obstacle (PRVO).