Robert G Ingersoll

en français

Contexte

] ou de la surveillance des infrastructures ferroviaires [Flammini, Pragliola, and Smarra 2016] ont suscité un intérêt croissant des entités commerciales pour ce domaine en pleine expansion. La diversité des applications met en évidence les difficultés théoriques telles que : une méthodologie de conception de contrôle incomplète, des ressources insuffisantes, des limitations physiques dans la communication, etc. L'un des défis de ces applications est la planification de mouvement qui vise à générer une route (en temps réel) dans un environnement complexe, en optimisant un coût donné sous des contraintes kinodynamiques [Hsu, Kindel, et al. 2002] et, simultanément, en maintenant la sécurité et l'intégrité du système en question.

The preparation of this thesis took a long and challenging odyssey, during which the help, friendship, and encouragement of so many people were essential. The following few lines are an attempt to express my gratitude to all of them.

First, I would like to thank my supervising team for guiding me through all aspects of research. I thank (Associate Professor) Ionela Prodan for her trust and kindness, even before the actual beginning of this thesis. I am grateful to (Professor) Sorin Olaru for his enthusiasm for research, insightful advice, and continuous encouragement that made me overcome all the challenging moments. I also thank my co-supervisor, (Professor) Silviu-Iulian Niculescu, for his constant support of my work and ensuring that my thesis is developing optimally. Working under their guidance has been a great opportunity and privilege. No words could ever express how indebted I am for all their support. I look forward to continuing our fruitful collaboration.

As well, I want to thank (Professor) Florin Stoican from the "Politehnica" University of Bucharest for the time he took to discuss new ideas and improvements in my work.

In my eyes (and not only mines), he is an excellent researcher who can always offer constructive and critical suggestions from writing to thinking. Through these suggestions and discussions, I recognized the relevance and realized the potential of my research.

I thank Professor Laurent Lefevre, director of LCIS, for welcoming me into the laboratory in my early days as a PhD student and softening the adaptation to the French culture. I am honored by the interest he has always shown for my work. As well, I am thankful to Professor Didier Dumur, the Head of the Automatic Department (within L2S), and all the professors in the department for their constant support. Furthermore, I would like to thank Dr. Jacques Blanc-Talon, who has regularly attended my work as a representative of DGA, for the fruitful discussions.

I am also grateful to the members that agreed to be part of my thesis committee. I want to thank Koichi Kobayashi, Associate Professor at Hokkaido University, and Mazen Alamir, Directeur de recherche CNRS at GIPSA-lab, for their careful reading of the manuscript and all their valuable and insightful comments. I also thank Vasso Reppa, Associate Professor at TU Delft, and Anthony Tzes, Professor at NYU Abu Dhabi, for accepting to be part of my committee and Antoine Girard, Directeur de recherche CNRS at L2S, for presiding it. Their questions and remarks during the defense are a worthwhile contribution to further the work presented in this thesis.

A special acknowledgment has to go to my former professors and collaborators. My thanks go to my highschool math and physics professors, Claudiu Militaru and Robertino Radu, who, with their classes, constantly increased my interest in the engineering field.

Reaching the end of my doctoral research is a significant achievement for me. To get to this point would not have been possible without the support of my friends. I have to thank Iulian, one of the best friends I had through our time in highschool and beyond, who has always been there to help and encourage me when I needed it. This last period also allowed me to meet many people who made the experience more enjoyable. I want to thank all my friends from both laboratories, LCIS and L2S, who were there to support me during the last years. My warmest thanks go to Thinh for his constant and communicative good mood and positivity, Igysso, for all our discussions on greek history and culture, Silviu, for his essential pieces of advice, Huy, for his calm and serenity, as well as the others Yoann, Thanos, Thomas, Joy, Martin, Jérémy, Vincent, Ugo, Antonello, Andreea, Ion, Gauthier, Dario, Matthieu, Kodjo, Fetra, Baptiste, Jean, Mert, Benjamin, and Merouane, with whom I spent less but valuable time.

One last thank you goes to my family (my parents and brothers), who granted me the support I needed to achieve all my achievements. Without them I would not be the person I am today.

Thank you very much to everyone.

Merci beaucoup à tous.

Mulțumesc mult tuturor.

Qu'il s'agisse de véhicules aériens sans pilote (aéronefs à voilure fixe, drones ou dirigeables), de véhicules terrestres sans pilote ( par exemple, des rovers planétaires), de véhicules de surface sans pilote (navires-drones) ou de véhicules en orbite autour de la terre, [Iagnemma and Dubowsky 2004]) jusqu'aux engins spatiaux en orbite terrestre et aux sondes dans l'espace lointain, les véhicules autonomes sont devenus une ressource exploitable pour une grande variété d'opérations. Que ce soit dans le domaine civil ou militaire, plusieurs véhicules robotisés sont actuellement utilisés (ou en cours de développement) afin de réduire les coûts et d'atténuer le risque humain associé à une mission. Par ailleurs, le niveau d'autonomie peut varier selon les catégories de véhicules téléguidés, certains de ces systèmes nécessitant un certain niveau de supervision humaine à partir d'une base ou d'une station au sol. Alors que la plupart des drones ont été conçus pour intégrer une commande à distance par un opérateur humain, leurs coûts d'exploitation ont encore été réduits en augmentant le niveau d'autonomie. Par exemple, au lieu d'être pilotés à distance, ils peuvent être programmés pour suivre une collection prédéterminée de points de cheminement. Ainsi, leur capacité à manoeuvrer en toute sécurité dans un environnement particulier implique une stratégie de navigation sûre. De plus, lorsque l'environnement particulier contient des zones interdites (par exemple, des obstacles), le véhicule doit être capable de calculer sa trajectoire en temps réel, c'est-à-dire au fur et à mesure du déroulement de la mission, en tenant compte de ses propriétés dynamiques et cinématiques.

Les méthodes de planification de mouvement [START_REF] Lavalle | Planning algorithms[END_REF]] sont bien établies et ont été étudiées sous plusieurs angles . Tous ces sujets couvrent un large spectre qui dépasse le cadre viii Contents de ce manuscrit, mais nous ne retenons que ceux d'entre eux qui font face à une complexité supplémentaire et dans lesquels la présence d'alternatives ou de prises de décisions discrètes est explicitement considérée. Dans le présent travail, l'accent est mis sur les applications spécifiques dans lesquelles des domaines faisables non convexes apparaissent et sont ensuite abordés à l'aide d'un cadre théorique des ensembles approprié.

Comme le suggère le titre du manuscrit, nous nous intéressons au problème de la planification du mouvement dans des environnements contenant une distribution spatiale dense d'obstacles, appelés environnements encombrés.Un exemple classique est celui d'un bateau se déplaçant dans une région de fjords. Afin d'arriver en toute sécurité à sa destination, le bateau doit suivre une trajectoire sans collision, en évitant de se frapper avec ces zones de navigation dangereuses. Par conséquent, l'identification de cette trajectoire implique la prise en compte de contraintes de positionnement non convexes (de sorte que le bateau n'entre pas en collision et évite en toute sécurité les rivages). Plusieurs questions se profilent alors :

• Comment repérer et modéliser ces régions interdites?

• Comment caractériser efficacement les régions faisables au sein desquelles le bateau peut naviguer en toute sécurité?

• Comment sélectionner la trajectoire la plus sécurisée (par rapport à un critère spécifique)?

Essentiellement, dans le manuscrit, la stratégie de navigation repose sur le fait suivant : l'agent (c'est-à-dire le robot ou le véhicule sans pilote) dispose d'une "awareness map", une carte qui contient toutes les informations disponibles sur l'environnement extérieur. Ainsi, la manière dont les obstacles sont modélisés joue un rôle important, impactant fortement les performances de la stratégie de navigation. Dans le présent manuscrit, nous ne nous intéressons pas à la sélection de la "meilleure" représentation interne/effective des régions interdites, mais plutôt à la manière dont ces régions peuvent être décrites efficacement afin de générer avec succès/efficacement des trajectoires (ou au moins de fournir les garanties de la non-existence d'une trajectoire faisable).

Les méthodologies et concepts de planification de mouvement, développés dans ce manuscrit, appartiennent, principalement, à la classe de la planification du mouvement basée sur les MPC(Model Predictive Control). Ils partent de formulations de type MIP (mixed-integer programming) et évoluent progressivement vers des convexifications.

Cette thèse peut être située dans la lignée des recherches de [START_REF] Janeček | Trajectory planning and following for UAVs with nonlinear dynamics[END_REF]], [R. [START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF]] ou [A. [START_REF] Richards | Trajectory planning[END_REF] et, de même, peut être considérée comme une extension de [Ionela Prodan 2012]. En d'autres termes, l'accent est mis sur les stratégies de navigation et de contrôle basées sur l'optimisation, en s'appuyant sur une combinaison de la commande prédictive (MPC) et des méthodes de la théorie des ensembles. Par rapport aux travaux précédents, nous développons des méthodes constructives (basées sur l'interprétation géométrique), offrant une nouvelle perspective dans l'exploration du problème de la planification du mouvement dans des environnements encombrés. La particularité du présent manuscrit vient du fait que nous avons abordé le problème de navigation à partir de différents niveaux de décision.

Contributions

Cette thèse présente plusieurs nouveaux concepts pour la planification de mouvement basée sur l'optimisation d'un véhicule dans un environnement encombré. La plupart des algorithmes, ix Contents techniques et idées présentés ici peuvent être considérés individuellement et certains d'entre eux appartiennent à un domaine plus large des MIP, de la théorie des ensembles ou des systèmes hybrides. Durant cette thèse, nous nous limitons aux problèmes de guidage pour illustrer les concepts théoriques liés à ce cadre. Les contributions de cette thèse peuvent alors être énoncées comme suit :

• fournir une revue de littérature détaillée des résultats de recherche récents et des questions ouvertes dans le domaine de la planification du mouvement dans un cadre mixte en nombres entiers. Ce travail peut être utilisé au profit des communautés de recherche sur le contrôle/commande et l'optimisation, ce qui permet de d'identifier rapidement les sujets de recherche antérieurs, opportuns et pertinents dans le domaine.

• Nous présentons une vue géométrique du problème de l'évitement des collisions en utilisant des surapproximations zonotopiques des obstacles. des obstacles. Nous avons souligné les avantages de choisir une famille particulière famille d'ensembles (zonotopes paramétrés) en ce qui concerne la complexité d'une représentation non convexe du domaine faisable. L'évitement des collisions basé sur le MPC (Model Predictive Control) est comparé à la technique heuristique PRM (Probabilistic Roadmaps). et a confirmé les avantages potentiels. Le problème de la sauvegarde de la topologie initiale en fournissant des conditions de séparation capables d'assurer le maintien des des passages parmi les obstacles est abordé dans le même cadre.

• fournir une solution constructive pour la génération d'un chemin entre deux points dans un environnement encombré, multi-obstacles. L'information globale sur la géométrie des obstacles est considérée comme un point d'entrée pour une procédure d'optimisation convexe qui conduit à un lifting convexe permettant le partitionnement de l'environnement encombré. Ce partitionnement est un élément clé pour décrire un graphe autour des obstacles et, finalement, pour la génération de corridors qui évitent les obstacles. D'un point de vue informatique, la validité de la construction repose sur la faisabilité de l'algorithme de lifting convexe. Cette faisabilité peut être améliorée par une reformulation des obstacles en termes d'un nombre fini de sous-ensembles convexes.. De plus, ce principe permet la généralisation de la construction pour la planification de chemins en présence d'obstacles non convexes. Dans ce dernier cas, un graphe déconnecté dans la génération des couloirs offrira une formelle de l'infaisabilité du problème de planification du chemin pour au moins une paire de points dans l'environnement encombré.

• l'introduction d'une stratégie de navigation fiable dans un environnement encombré, en employant une variation du schéma MPC (Model Predictive Control) générique, fournissant des garanties de faisabilité et de convergence en exploitant l'accessibilité. En outre, nous fournissons une adaptation pour les environnements multi-obstacles variant dans le temps et pour un contexte multi-agents.

x

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d). The corridors for k ∈ {3, 6} . . . . . . . . . . . . . . . . . . . . . . . . . . . (e). The corridor for k = 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10. The solution relying on convex hull. . . . . . . . . . . . . . . . . . . . . . . . (a). The space partitioning considering the convex hull of all possible displacement of the moving obstacle . (a). X-Y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). X-Y-t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c). X-t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d). Y-t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.2. The control scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.3. Simulation results for N p = 10 P = 10 4 I 2d , Q = 10 6 I 2d , R = I d and Q z = 10 5 I 2d .

(a). p z , ϵ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). X-t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiii

List of Tables

Feasibility and convergence aspects for approaches in

Acronyms

MI(L)P

Notation

Throughout this manuscript we use the following standard notations.

The logical operators: ∧ (AND), ∨ (OR) and ¬ (negation).

Let R, Z and N denote the field of real numbers, the set of integers and the set of nonnegative integers, respectively. Notations R n and R m×n denote the vector field and the matrix field of real numbers, respectively. We use the similar notation for the sets of integer and non-negative integers.

For a discrete set I, #I represents its cardinality.

For x ∈ R n we denote x 2 Q = x ⊤ Qx.
Given a compact set S ∈ R n , C X (S) denotes the complement of S over X ∈ R d , and int(S) its interior. Conv(S) is the convex hull of S ∈ R n , Com(S) the space of compact subsets of S, and V(S) the set of its extreme points (possibly infinite collection of points spanning Conv(S)).

We consider the following set operations:

• Minkowski sum of two sets A, B ∈ R n :

A ⊕ B = {a + b | a ∈ A, b ∈ B} • Pontryagin difference of two sets A, B ∈ R n : A ⊖ B = {a ∈ A | a + b ∈ A, ∀b ∈ B}
Any polytope (i.e. a bounded polyhedron) has a dual representation in terms of :

• intersection of half-spaces representation:

P = {x ∈ R n : s ⊤ i x ≤ r i , ∀i ∈ I} • convex hull of extreme points v j ∈ R n , j = 1, . . . , ℓ: P = {x ∈ R n : x = α j v j , α j = 1, α j ≥ 0, ∀j}.
For any polyhedron P ⊂ R d , V(P ) is the (finite) set of its vertices and F k i (P ) is the i-th face of the dimension k < d.

Background and Motivation

General Context

In the last decades, surveillance, patrolling and area assessment activities in various complex environments are becoming ubiquitous. A vast majority of these activities is either dangerous or relatively difficult for human presence due to operating conditions and environmental impact, cost effectiveness, and similar ones. In these situations, the missions can be more efficiently performed either by a robot or, with increased accuracy, by a cooperative team of robots (also interpreted as multi-agent systems) which can partition a given task in order to accomplish it in an effective manner. Applications in precision agriculture [START_REF] Jawad | Energy-efficient wireless sensor networks for precision agriculture: A review[END_REF]], marine support operations [START_REF] Tran | Motion planning for multi-agent dynamical systems in a variable environment[END_REF]], underwater monitoring [Yilmaz et al. 2008] or railway infrastructure monitoring [START_REF] Flammini | Railway infrastructure monitoring by drones[END_REF] lead to a growing interest of commercial entities for this rapidly expanding domain. The diversity of application highlight the theoretical difficulties such as: incomplete control design methodology, insufficient resources, physical limitations in communication, etc. One of the challenges in such application is the motion planning which aims to generate a (real-time) route within a complex environment, optimizing a given cost under kinodynamic constraints [START_REF] Hsu | Randomized kinodynamic motion planning with moving obstacles[END_REF]] and, concurrently, maintaining the safety and integrity of the system at hand.

Roughly speaking, motion planning seems quite elementary, given the fact that humans deal with it instinctively and with no perceptible difficulty in their everyday lives. Literally, this apparent simplicity is a matter of perception, since the basic principles used by human beings for interaction with their environment are utterly challenging for autonomous robots. Indeed, there exists a plethora of trivial methods that can produce satisfactory or, sometimes, quite remarkable results, but their shortcomings can be readily brought out into relief when operating within complex environments. Therefore, in order to develop a safe navigation strategy, the use of nontrivial mathematical tools and tailored algorithmic techniques cannot be avoided.

Ranging from unmanned aerial vehicles (fixed-wing [START_REF] Beard | Small unmanned aircraft: Theory and practice[END_REF], drone [Robin Deits and Russ Tedrake 2015] or blimps), unmanned ground vehicles (e. g., planetary rovers [Iagnemma and Dubowsky 2004]), unmanned surface vehicles (drone ships) [START_REF] Liu | Unmanned surface vehicles: An overview of developments and challenges[END_REF]] to earth-orbiting spacecrafts and deep-space probes, autonomous vehicles have become an exploitable resource for a wide variety of operations. Either in civilian or in military domain, several robotic vehicles are currently employed (or being developed) in order to reduce the costs and to alleviate the human risk associated with a particular mission, i.e., unmanned systems can be made more economical and smaller than their manned counterparts due to the absence of operator-safety constraints. As a side remark, the level of autonomy may differ among the classes of unmanned vehicles, some of these systems require a certain level of human supervision from a base/ground station. For instance, unmanned aerial vehicles are capable of high maneuverability and are extensively exploited within hostile environments. While most UAVs were designed to integrate a remote control by a human operator, their operating costs were further decreased by enhancing the level of autonomy. As an example, instead of being remotely flown, they can be programmed to track predetermined collection of waypoints. Thus, their capability to safely maneuver within a particular environment involves as fundamental part a safe navigation strategy. Moreover, when the particular environment contains forbidden regions (e.g., obstacles), the vehicle should have the capability to compute its path in real-time, i.e., as the mission unfolds, thereby accounting for its dynamic and kinematic properties.

As well, the automotive industry is concentrating its efforts towards the autonomous driving concept. Under this paradigm, the vehicle handles entirely the driving task while the human is playing a passive role. Recently, many companies are developing prototypes able to control the vehicle behavior without human intervention. However, there is a gap between these prototypes and their availability for the general public. There are diverse reasons, but the concept needs to gain in maturity and to ensure a safe and legal interaction with all components of the environment. A first and important step represents the development of the Advanced Driving Assistance Systems (ADAS), which permits the techniques to be perfected and gradually to get closer to Autonomous Driving desideratum [START_REF] Levinson | Towards fully autonomous driving: Systems and algorithms[END_REF]]. Roughly speaking, the term ADAS covers the technological systems which aim to assist the driver, taking over the control of the vehicle in certain situations, e.g., parking lots or highways, necessary to enhance the safety and optimizing the comfort. This type of systems can be viewed as a first generation of semi-autonomous driving and paves the way to automated vehicles [Ballesteros-Tolosana et al. 2017]. The main objective of Autonomous Driving is to compensate the human error in order to achieve an accidentfree traffic [START_REF] Beiker | Legal aspects of autonomous driving[END_REF]]. This means the performances of autonomous vehicle need to be at a superior level compared to that of a human driver in any standard scenario. As a matter of fact, the driving tasks can be divided into basic functional components which can be technically implemented and developed up to a certified level of maturity. By any means, the requirement of a safe motion planning strategy cannot be neglected.

The motion planning methods [START_REF] Lavalle | Planning algorithms[END_REF]] are well established and have been studied in time optimal control, nonlinear control, stabilization, reachability and other related topics. All these topics cover a large spectrum which goes beyond the scope of this manuscript, but we retain only those of them which face additional complexity and in which the presence of alternatives or discrete decisions making is explicitly considered. The emphasis of the present work is placed on the specific applications where non-convex feasible domains appear and are subsequently tackled using an appropriate set-theoretic framework.

Motivation

Over several decades, motion planning has been a paramount research topic in various fields, e.g., control, robotics or, more recent, artificial intelligence. Initially, the design methods dealt with trajectory and/or path planning for holonomic systems performing in obstacle-free environments and were mainly based on optimal control and/or nonlinear programming techniques. The interested reader can find a survey of such numerical methods in, e.g., [START_REF] Latombe | Robot motion planning[END_REF]]. Gradually, the presence of the obstacles could not be neglected anymore, and the research focused on the integration of these notions in the design. As stated, the present thesis focuses on recent techniques which enable such an integrated design (relying, e.g., on model predictive control), but we do not completely disregard earlier methods, aiming to constantly compare our approaches to them. As a side remark, many of these methods were developed in the robotics community and some of them yielded fundamental complexity results.

As the title of the manuscript suggests, we are interested in addressing the motion planning problem within environments containing a dense spatial distribution of obstacles, so-called cluttered environments. A classical example, a boat moving within a fjord region, as illustrated in Figure 1. 1. In order to safely arrive to its destination, the boat should follow a collision-free path, staying away from an undesired crash with those dangerous zones of navigation. Hence, the identification of that paths involves the consideration of non-convex position constraints (such that the boat does not collide and safely avoids the shorelines). Several questions become apparent:

• How to detect and model such forbidden regions?

• How to efficiently characterize the feasible regions within which the boat may safely navigate?

• How to select the most safe (w.r.t. a specific criterion) trajectory/path?

Basically, in the manuscript, the navigation strategy relies on the following fact: the agent (i.e., robot or any unmanned vehicle) has a so-called awareness map which contains all available information about the external medium. Thus, the way in which the obstacles are modeled plays an important role, strongly impacting the performances of the corresponding (design and control) strategy. In what follows, we are not interested in selection of the "best" internal/effective representation of the interdicted regions, but rather in how these regions can be efficiently described in order to successfully/effectively generate trajectories/paths (or at least to provide the guarantees of the non-existence of a feasible trajectory).

One of the inclusive approaches in tackling the enumeration and selection among alternatives in these problems is represented by mixed-integer programming (MIP), i.e., optimization for which part or all of the arguments are required to be integers [START_REF] Jünger | 50 Years of integer programming 1958-2008: From the early years to the state-of-the-art[END_REF]Ionela Prodan, Stoican, et al. 2015]. For many years, MIP problems have been actively bridging the gap between academia, industrial as well as business fields of the construction of the graph is incremental, the algorithm stops when a large enough set of collision-free paths is attained. Thus, a collision-free sample is added as a node in the graph and is linked with the surrounding nodes. The obtained graph is actually a tree.

As PRM, RRT has a variety of versions. Some take into account the equations of motion and generate attainable paths [START_REF] Lavalle | Randomized Kinodynamic Planning[END_REF], others generate only geometric paths which became reference trajectories for a lower level controller. Moreover, some versions are tailored for complex and/or unstable dynamics [START_REF] Leonard | A perception-driven autonomous urban vehicle[END_REF]] or for uncertain dynamics [Weiss et al. 2017].

Potential field approaches rely on the construction of a scalar function (so-called, the potential). This function takes high values when the agent stays within the interdicted zones. In the collision-free workspace the function is decreasing towards the goal configuration (i.e., the potential associated to the destination point is minimal). Thus, the agent may attain the final point moving in the direction of the negative gradient of the potential. [START_REF] Rimon | Exact robot navigation using artificial potential functions[END_REF] provides an historical review on the potential field formulation and how this approach is involved in motion planning. An interesting characteristics is that the potential field formulation is frequently used in decentralized or distributed control strategies [START_REF] Filotheou | Decentralized control of uncertain multi-agent systems with connectivity maintenance and collision avoidance[END_REF]]. An issue of the method is the possibility of traps (local minima), but that can be overcome by using ideas such as a randomized potential field. As a side remark, while MIP-based methods take explicitly into account the constraints and lead to a constrained optimization problem, the potential field-based formulations [Y.-b. Chen et al. 2016] relax the constraints by adding penalty terms in the cost.

Optimization-based control approaches cover, generally speaking, the control design based on optimization techniques, including the classical optimal control, the LMI-based techniques or interpolation-based techniques. A notably useful approach of optimizationbased control for guidance of an agent within multi-obstacle environment is receding horizon technique, also known as Model Predictive Control (MPC). The important assets of the predictive control approach are given by its capability to handle generic state and control constraints, allowing a system to operate closer to its constraints and, thus, often emerge in better performance [START_REF] Murray | Optimization-based control[END_REF]]. Moreover, there exists a variant of this classical MPC approach, MIP-based MPC, which combines the benefits of MIP and MPC.

The main advantage of using MIP for motion planning is the ability to deal with hard constraints induced by obstacle and collision avoidance conditions. Moreover, the algorithms are characterized by completeness (that is, they return a feasible solution whenever it exists) and provide a solution to the general motion planning problem. However, as stated elsewhere, the approach needs to be enhanced and refined for complex and large-scale problems. The class of problems which can be handled in this framework is increasing due to computer performances and the proliferation of software tools (e.g., [START_REF] Cplex | V12. 1: User's Manual for CPLEX[END_REF][START_REF] Gurobi | Robotic swarm control from spatio-temporal specifications[END_REF]), implementing efficiently the state of the art resolution algorithms. Noteworthy, in the literature, beside the software packages, a considerable effort was put in finding a technique which allows the translation/relaxation from non-convex to convex formulations by mitigating the gaps. This is referred as convex relaxation or convexification of the non-convex optimization problem. In the literature there are various works providing methods for convexification under several labels: e.g., convex relaxation [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF], successive convexification [Mao, Szmuk, and Acikmese 2016] or timevarying constraints [Frasch et al. 2013]. The interested reader is referred to Section A.2 for more details on the convex relaxation.

The motion planning methodologies and concepts, developed in this manuscript, belong, mainly, to the class of MPC-based motion planning. They are starting from MIP type of formulations and gradualy move towards convexifications. For the interested reader Section A provides a brief overview of the above mentioned alternatives.

Thesis orientation 1.2.1. Proposed approach

This thesis can be situated in a line with the research of [START_REF] Janeček | Trajectory planning and following for UAVs with nonlinear dynamics[END_REF],[R. [START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF] or [A. [START_REF] Richards | Trajectory planning[END_REF] and, as well, can be viewed as an extension of [Ionela Prodan 2012]. In other words, the focus is placed on optimization-based navigation and control strategies, building on the combination of Model Predictive Control (MPC) and set theoretic methods. Compared to the previously/above-mentioned works, we develop a couple of constructive methods (based on the geometrical interpretation), offering a new perspective in the exploration of motion planning problem within cluttered environments. The particularity of the present manuscript comes from the fact that we tackled the navigation problem from several decision layers.

Firstly, based on a thorough analysis of the recent results in the field, we focus on the representation of the multi-obstacle environment. Thus, having a collection of interdicted/forbidden regions (obstacles) we model them using sets with respect to their intrinsically properties (e.g. convexity and connectivity). The complement of these forbidden region designate the domain of interest, and can be viewed as a combination of (finite) alternatives, which can be readily modeled via MIP techniques. For instance, let us consider the trivial example depicted in Figure 1.2, where an agent needs to pass without collision an obstacle. There are two alternatives, as shown in the figure, and the navigation strategy has to select one of them with respect to an a priori known objective. As can be noted, the selection has to be done from a limited number of alternatives. However, an increase of the number of obstacles and a more dense spatial distribution, so-called cluttered environment, lead to a finite but huge number of alternatives3 , impacting the computational effort. For example, procedures like cell merging [I. Prodan et al. 2012] or logarithmic formulations [J. P. Vielma and Nemhauser 2011] are employed to reduce the complexity but their automatic treatment becomes numerically complex for large numbers of obstacles and/or agents. Therefore, a conceptual update of the approach, which simplifies the eventual complexity, is necessary. To this end, we opted to exploit in this thesis a particular class of convex sets endowed with the symmetry property in order to model the cluttered environment. In other words, considering the obstacles avoidance condition as geometrical problems we revisit the MIP framework, aiming to take advantage of the geometrical insight and properties (as, e.g., symmetry) in order to reduce the complexity and enhance the performance. As a consequence, this permits the reduction the computational cost and an efficient handling of an increased number of constraint. Moreover, in this manner, we will be able to provide a rigorous evaluation of the consequent improvements. However, even though this approach enhances the classical MIP-based motion planning, the physical and computational limitation may affect the performances of a real-time implementations as the non-convexity of the formulation is not leveraged.

Eventually, we orient to the treatment of the navigation problem from a tactical perspective, and splitting the global problem using local information based on graph theory. Hence, the original problem will be addressed aiming to reduce its complexity and, concurrently, to ensure the convergence and the feasibility. To this end, using a set theoretic tool we compute a partition of the cluttered environment taking into account the spatial distribution of the obstacles. Hence, we aim to identify a collection of geometrical paths passing between any pair of obstacles. Next, having this collection and employing graph theory, we focus on generating a feasible reference path which represents a starting point

Contributions

This thesis presents several new concepts for optimization-based motion planning of single vehicle within a cluttered environment. The most of the algorithms, techniques and ideas which are presented here can be considered individually and some of them belong within a broader field of MIP, set-theory or hybrid systems. In this manuscript, though, we restrict ourselves to guidance problems to illustration of the theoretical concepts related to this framework. The contributions of this thesis can then be stated as follows:

providing a detailed literature review of recent research results and open issues in the field of motion planning in a mixed-integer framework. This work can be employed to the benefit of both control and optimization research communities allowing to swiftly identify previous, timely and relevant research topics in the field.

presenting a geometric view for the collision avoidance problem using zonotopic overapproximations of the obstacles. We emphasized the benefits of choosing a particular family of sets (parametrized zonotopes) regarding the complexity of a non-convex feasible domain representation. The MPC(Model Predictive Control)-based collision avoidance is compared to the heuristic PRM (Probabilistic Roadmaps) technique and confirmed the potential advantages. The problem of safeguarding the initial topology by providing that separation conditions able to ensure the maintaining of the passages among the obstacles is addressed within the same framework.

providing a constructive solution for the generation of a path between two points in an environment obstructed by multiple obstacles in a d-dimensional space. The global information on the geometry of the obstacles is considered as an entry point for a convex optimization procedure which leads to a convex lifting allowing the partitioning of the cluttered environment. This partitioning is a key element for describing a graph around the obstacles and ultimately for the generation of corridors which avoid obstacles. From the computational point of view, the validity of the construction relies on the feasibility of the convex lifting procedure. The feasibility can be enhanced by a reformulation of the obstacles in terms of a finite number of convex subsets. Furthermore, this principle allows the generalization of the construction for path planning in the presence of multiple non-convex obstacles. In the latter case, a disconnected graph within the generation of the corridors will offer a formal proof of infeasibility of the path planning problem for at least a pair of points in the cluttered environment.

introducing a safe navigation strategy within a cluttered environment, employing a variation of the generic MPC (Model Predictive Control) scheme, providing feasibility and convergence guarantees by exploiting reachability. Moreover, we provide an adaptation for time-varying multi-obstacle environments and for a multi-agent context.

Organization of the manuscript

The dissertation is organized as follows (see also Figure 1.3 ).

Chapter 2 presents a review of past and present results and approaches in the area of motion planning using MIP (Mixed-integer Programming). Although in the early 2000s MIP was still seen with reluctance as method for solving motion planningrelated problems, nowadays, due to increases in computational power and theoretical advances, its extensive modeling capabilities and versatility are coming to the fore and enjoy increased application and appreciation. This class of control problems involves, essentially, either a selection from a limited number of alternatives or a constrained optimization problem over a non-convex domain. In both situations, MIP has proven to be an efficient modeling technique as it will be shown in the sequel. The material is based on:

• Ioan, D., Prodan, I., Olaru, S., Stoican, F., Niculescu, S. I. (2020). Mixedinteger programming in motion planning. Annual Reviews in Control.

Chapter 3 Maintaining a classical MIP framework for the motion planning problem, this chapter provides: i) various measures for tight zonotopic approximations of the obstacles, guarding the original space topology ii) a rigorous bound for the number of generators such that the complexity is reduced;

iii) propose a navigation strategy with obstacle avoidance guarantees using local zonotopic approximations.

iv) provide a detailed analysis of evolution of the hyperplane arrangements corresponding to a dynamical multi-obstacle environment.

The chapter is based on the following publications:

• Ioan, D. Chapter 4 addresses the path-planning level with a primary objective of global feasibility.

In other words, our method lead either to a feasible path or to a certificate of infeasibility. This certificate ensures that the environment does not allow the construction of a path between the source and the destination. The differential constraints are discarded from the problem formulation as well as the other limitations that may appear in motion planning due to limited steering or energy. Optimality is a secondary objective for the generation of a geometric path and has the potential to explicitly describe a feasible corridor as in [S. Liu et al. 2017]. The main contributions are: i) provide a partitioning of the workspace;

ii) obtain a feasible geometric paths;

The chapter is based on the following publications:

• Ioan, D. Chapter 5 addresses the navigation problem in a multi-obstacle environment by employing convex lifting. The particularity comes from the following fact. We consider an adaptation of the generic MPC (Model Predictive Control) trajectory tracking problem, aiming to guarantee the feasibility and the convergence. Simulation results and proof of concepts illustrations prove the effectiveness of the proposed approach.

The main contributions of this chapter are:

i) providing a detailed analysis of the recursive feasibility of a corridor-constrained MPC;

ii) providing a generic navigation strategy for cluttered environments with feasibility guarantees ("Relay MPC") using the corridors as robustness measure;

iii) introducing the necessary adjustments for time-varying cluttered environments; iv) introducing a reachability-based approach for navigation of multi-agent systems.

The chapter is based on the following publication:

• Ioan, D., Prodan, I., Olaru, S., Stoican, F., Niculescu, S. I. (2020, July).

Navigation in cluttered environments with feasibility guarantees. In IFAC 2020-21st IFAC World Congress.

Chapter 6 completes the thesis with conclusions and discussions of future directions.

Graphically, the previously described chapters can be viewed as shown in Figure 1.3, where the edges suggest dependencies between the various approaches. . . . . . . . . . . . . . . . . . . . . . . . . . . 74 In the first part of this manuscript, the main objective is to identify and summarize the state of the art of MIP-based motion planning. Hence, in the sequel, we place less emphasis on the other control areas employing MIP, even if throughout the chapter we occasionally refer the interested readers to the references covering the broad MIP-based control topics.

■♥tr♦❞✉❝t✐♦♥ ❈❤❛♣t❡r ✶ ▼♦t✐♦♥ ♣•❛♥♥✐♥❣ ❙❡t✲t❤❡♦r❡t✐❝ t♦♦•s ▼✐①❡❞✲✐♥t❡❣❡r t❡❝❤♥✐q✉❡s ▼✉•t✐✲♦❜st❛❝•❡ ❡♥✈✐r♦♥♠❡♥ts ❖♣t✐♠✐③❛t✐♦♥✲❜❛s❡❞ ❝♦♥tr♦• ❈❤❛♣t❡r ✷ ❖❜st❛❝•❡ ❛✈♦✐❞❛♥❝❡ ❈♦♠♣•❡①✐t② ❜♦✉♥❞s )♦♥♦t♦♣✐❝ r❡♣r❡s❡♥t❛t✐♦♥ ▼♦✈✐♥❣ ♦❜st❛❝•❡s ❈♦♥tr♦• ❞❡s✐❣♥ ♣❡rs♣❡❝t✐✈❡ ❈❤❛♣t❡r ✸ P❛t❤ ♣•❛♥♥✐♥❣ ❙♣❛❝❡ ♣❛rt✐t✐♦♥✐♥❣ ❈♦♥✈❡① •✐❢t✐♥❣ ❈♦rr✐❞♦rs •r❛♣❤ t❤❡♦r② ❈❤❛♣t❡r ✹ ◆❛✈✐❣❛t✐♦♥ ❋❡❛s✐❜✐•✐t② ❣✉❛r❛♥t❡❡s ❚r❛•❡❝t♦r② ♣•❛♥♥✐♥❣ ▼✉•t✐✲❛❣❡♥t ♣❡rs♣❡❝t✐✈❡ ❚✐♠❡✲✈❛r②✐♥❣ ❡♥✈✐r♦♥♠❡♥ts ❈❤❛♣t❡r ✺ ❈♦♥❝•✉s✐♦♥s ❋✉t✉r❡ ❉✐r❡❝t✐♦♥s ❈❤❛♣t❡r ✻
During more than 60 years of existence, the field of integer programming was extensively studied in the mathematics community due to its modeling capability and flexibility. In recent years (mainly the last two decades), mostly owing to the growing computational power, the integer programming was brought to the attention of the control and robotic communities. There exists a broad variety of decision making problems that can be dealt through a MIP framework/approach.

MIP (Mixed-integer Programming

) is a mathematical optimization problem in which some or all the variables are integers.

As its name indicates, MIP (Mixed-integer Programming) represents a mathematical optimization problem in which the objective is a linear, quadratic function or sometimes a more general criterion to be minimized or maximized, the constraints are linear (or nonlinear) equalities (or inequalities) and there exist some (non-empty) subsets of integer variables playing the role of arguments [START_REF] Williams | Model building in mathematical programming[END_REF][START_REF] Jünger | 50 Years of integer programming 1958-2008: From the early years to the state-of-the-art[END_REF]].

MIP is used to model several design problems and decision processes.

In a larger perspective, MIP is used to model several design problems and decision processes. Consider a typical logistics problem: an airport, which serves on average 50 flights per hour. The airport has only four runways. The task assignment problem that appears is to assign flights to runways, such that the runways are efficiently and uniformly used, while respecting some regulations (e.g., time separation between two consecutive landings/takeoffs or a minimum distance between two runways for simultaneous takeoffs). Another classical situation is described by the well-known traveling salesman problem and its variations, where the salesman wants to visit a number of customers in a minimal time or to cover a minimal distance. This has applications in several domains (e.g, overhauling gas turbine engines or X-ray crystallography [START_REF] Matai | Traveling salesman problem: an overview of applications, formulations, and solution approaches[END_REF]). The above problems can be solved either intuitively, based on experience or by a trial and error method, but for critical situations an accurate mathematical formulation is necessary in view of certification. There are of course many use cases which may employ MIP in applications closer to control and robotics. For example, the optimal power flow in the energy transmission networks [START_REF] Bahiense | A mixed integer disjunctive model for transmission network expansion[END_REF] or the transportation problems in a cluttered environment. Consider again the example of a boat moving within a fjord region. In order to safely arrive to its destination, the boat should follow a given path and avoid collision with the fjords. Thus, the feasible region is non-convex and should be efficiently described in terms of alternatives for the maneuvering.

In the following, a brief classification of the types of problems, which can be modeled through MIP is provided. A first class of problems is designated by those that involve integer quantities (i.e. discrete/quantified inputs or outputs), e.g. the knapsack problem [START_REF] Williams | Model building in mathematical programming[END_REF]]. For this type of problem, MIP may not come as the obvious, natural, first choice, but, usually, it represents a better solution than a classical approach (e.g., use of the classical linear programming and approximate the provided solution to the nearest integer value).

Another MIP-modelisable class of problems, closer to control and robotics, is the one involving logical conditions, extensively treated in: [Alberto Bemporad and Morari 1999;[START_REF] Williams | Model building in mathematical programming[END_REF][START_REF] Scott | Dynamic optimal power flow in microgrids using the alternating direction method of multipliers[END_REF]. For example, in [Alberto Bemporad and Morari 1999], using the notations from [START_REF] Williams | Model building in mathematical programming[END_REF]], boolean algebra tools are aggregated, which allows to transform logical conditions on continuous variables into mixed-integer inequalities (linear inequalities involving continuous and binary variables). These inequalities are mixed with the dynamics to model hybrid systems behaviour.

As well, MIP is a popular modelling tool for sequencing and/or allocation problems (also, named combinatorial problems) [START_REF] Scott | Dynamic optimal power flow in microgrids using the alternating direction method of multipliers[END_REF], including here the typical task assignment problem and its variations (e.g travelling salesman problem [START_REF] Dantzig | Solution of a Large-Scale Traveling-Salesman Problem[END_REF]). This class of problems can be easily extended to networks (and graph theory, as a generalization) problems: resource allocation on a PERT (Project Evaluation and Review Techniques) network [START_REF] Williams | Model building in mathematical programming[END_REF]].

Lastly, but most importantly for this chapter's purpose, MIP turned out to be a captivating method to model non-linearity [Alberto Bemporad and Morari 1999;[START_REF] Vielma | Mixed Integer Linear Programming Formulation Techniques[END_REF] and/or non-convexity [Ionela [START_REF] Prodan | Mixed-integer representations in control design: Mathematical foundations and applications[END_REF][START_REF] Richards | Trajectory planning[END_REF]. A plethora of control engineering problems are naturally and intrinsically characterized by non-linearity and/or non-convexity. For this reason and due to the increasing interest in optimization-based control [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]], the MIP has became an essential technique, which allows to include logical decisions and non-convex constraints in the optimization problem. Therefore, MIP's presence in control can be perceived in: piecewise-affine system identification [START_REF] Roll | Identification of piecewise affine systems via mixed-integer programming[END_REF][START_REF] Bemporad | Identification of hybrid systems via mixed-integer programming[END_REF], assignment problems [START_REF] Alighanbari | Coordination and control of multiple UAVs with timing constraints and loitering[END_REF], persisting exciting control [START_REF] Marafioti | Persistently exciting model predictive control for siso systems[END_REF]], control of hybrid systems [Alberto Bemporad and Morari 1999], fault detection [START_REF] Stoican | Reference governor design for tracking problems with fault detection guarantees[END_REF] or motion planning [Ionela Prodan, Stoican, et al. 2015;A. Richards and J. P How 2002].

MIP Formulations

The academic effort towards formal use of MIP can be trace back to the early '80 (or even earlier), the literature emphasizing at that time the modeling capabilities of MIP. All these formulations share a common characteristic: the encoding of discrete decisions using binary and/or integer variables. These decisions appear in different problems, each using a certain formulation. This section reviews them and provides a brief description of the most used MIP techniques and, concurrently, introduces some basic theoretical notions and tools.

Although generalized disjunctive programming (GDP) is not explicitly used in motion planning, we succinctly present it for the sake of the generality of modeling disjunctions through MIP.

Generalized Disjunctive Programming

Generalized Disjunctive Programming (GDP) has appeared for the first time in [START_REF] Raman | Modelling and computational techniques for logic based integer programming[END_REF] in the research effort to use both quantitative and qualitative information in order to optimally solve chemical engineering problems. To this end, the qualitative information is represented using disjunctions and logic propositions. Compared to MIP, the GDP approach has a relatively more compact formulation, because the logical conditions are not transformed using boolean algebra and inequalities but rather in their natural (logic) form. In other words, GDP represents a combination of algebraic and logical equation as shown by the typical GDP from (2.1):

min x f (x) + k c k (2.1a) s.t. r (x) ≤ 0, x ∈ R n , c k ∈ R, (2.1b
)

j∈J k   Y jk g jk (x) ≤ 0 c k = γ jk   , k ∈ K, (2.1c
)

h(Y ) = true, Y jk ∈ {true, false}, (2.1d)
where r (x) is a generic constraint, which does not depend on logical decisions; c k describes the cost variables, γ jk fixed charges and K is the index set of disjunctions.

g jk (x) ≤ M jk (1 -y jk ), (2.2c
)

j∈J k y jk = 1, k ∈ K (2.2d)
where M jk are "big-M" parameters, represented by sufficiently large constants.

Remark 2.1. This "big-M" formulation For illustration, we consider a common case of modeling the logical implication: "x > 0 -→ α = 1" for a variable x which evolves implicitly within a bounded domain2 of variation X = {x | |x| ≤ x}. The "big-M" formulation simply reverts the logical implication in terms of an inequality x-M α ≤ 0. Intuitively, as long as M is chosen to satisfy M ≥ x, the implication follows naturally. Nonetheless, by considering a too large value of the coefficient M we increase the size of the searching domain. Thus, the idea is to take M as small as possible (but still large enough to fulfill its relaxation role). In our simple example a proper choice is:

M = max x∈X x, i.e., M = x. ♦
Moreover, the cost is reformulated by rewriting each cost variable c k as the product γ jk y jk .

As well, the condition (2.1d) is written in an algebraic form Ay ≤ a.

The apparent advantage of the GDP over MIP is diminished because all the existing GDPdedicated solvers are based on a MIP reformulation. It is enough in general to assume that a direct MIP modeling of the problem at hand could lead to a better and more compact model. In spite of eluding GDP in such a way, there is another possibility to benefit from GDP. The problem can be modeled as a GDP and the resulting MIP reformulation can be parsed using one of the several tools, as detailed in [START_REF] Williams | Model building in mathematical programming[END_REF]].

Geometric viewpoint for MIP

Complex control synthesis or design builds most often on a constrained optimization problems and whenever a non-convex feasible domain needs to be handled, the ability of MIP to mathematically formalize disjunction constraints can be exploited. An active research topic was the MIP efficient description of such non-convex regions. Preliminary results make use of hyperplane arrangements and, hereinafter, we briefly recall these results which involve set-theoretic notions. We refer the readers interested in a deeper technical description to complementary materials referenced in this manuscript(e.g., [J. P. Vielma and Nemhauser 2011;[START_REF] Vielma | Mixed Integer Linear Programming Formulation Techniques[END_REF] or [Ionela Prodan, Stoican, et al. 2015]).

Most sets involved in MIP formulations are polyhedral ones due to their linear constraints representation. These linear constraints lead to polyhedral sets whenever the complements are considered. In what follows, we use the notion of polytope which is a bounded polyhedral set and has a dual representation in terms of intersection of half-spaces or convex hull of extreme points: P = {x :

s ⊤ i x ≤ r i , ∀i ∈ I} = {x : x = α j v j , α j = 1, α j ≥ 0, ∀j ∈ J }.
Consider a finite collection of hyperplanes from R d

H i = x ∈ R d : h i x = k i , i ∈ I (2.3)
where I ≜ {1 . . . N }, and

(h i , k i ) ∈ R 1×d × R.
Each of these hyperplanes divides the space in two disjoint 3 regions:

R + i = x ∈ R d : h i x ≤ k i , (2.4a
)

R - i = x ∈ R d : -h i x ≤ -k i . (2.4b)
A polytope P is a bounded intersection of these half-spaces 4 :

P = i∈I R - i , (2.5)
Its complement (up to its relative interior, see also footnote 1) C(P ) ≜ cl(R d \ P ) over R d represents the union of all regions covering the entire space except P . The fact that the union and intersection operators interchange w.r.t. the complement operator, we write: .6) The region (2.6) is a finite union of convex sets 5 and for a tractable characterization we may employ mixed-integer techniques. Hence, the binary variables (α 1 . . . α N ) ∈ {0, 1} N are introduced to obtain the representation of the polytopic set in the extended space including auxiliary binary variables:

C(P ) = C i∈I R - i = i∈I C R - i = i∈I R + i . ( 2
The relative interiors of these regions are disjoint, but their closures have as a common boundary the hyperplane H i The sign "-" was chosen for the simplicity of notation, any other feasible combination of signs from (3.16a)- (3.16b) could be chosen in order to describe the polytope P . In general (2.6) is non-convex with notable exceptions: the complement of an empty set, the complement of the unconstrained space or the complement of a half-space (possibly in a non-minimal representation).

h i x ≤ k i + M α i , i ∈ I, (2.7a) i∈I α i ≤ N -1, (2.7b)
where M is a big-M coefficient as in Remark 2. 1.

The requirement that a point is outside of the polytopic set P needs to be converted to the condition that the point has to be within the complement of at least one of the half-spaces defining the polytope.

Remark 2.3. The conditions (2.7a)-(2.7b) describe the region (2.6) via a proper combination of binary variables. For instance, the region R + i is expressed by (2.7a) with the following binary variables: .8) When a binary variable takes the value "1" the associated inequality describes, for the limit case (M → ∞, i.e., sufficiently large M ) the entire domain R n . Hence, the condition (2.7b) is necessary in order to ensure that at least one of the binary variables is "0" and, thus, that at least an inequality remains active. ♦ Remark 2. 4. As stated in Remark 2.1, choosing the value of the "big-M" constant may induce redundancy which can lead to increasing complexity in terms of the computational effort. It is straightforward that the constraint (7b) cannot describe the entire R n for a α i = 1, but rather a domain including the region of interest, X. Assuming, for instance, that the union of interdicted regions lies in a bounded cluttered environment X, the value of M is obtained via the following LP:

(α 1 . . . α N ) = (1 . . . 1 0 i 1 . . . , 1). ( 2
M = max i M i , M i = max x∈X (k i -h i x) .
Clearly, we may consider a different value of M for each half-space, that is, M i . In Figure 2.1 an illustrative example is depicted. Specifically, while the constraint is active (i.e., α 1 = 0), the resulting inequality: -h i x ≤ -k i describes the red region. Once α 1 = 1 (inactive/relaxed constraint), the selection of M gives the measure of relaxation. Thus, the idea is to relax sufficiently each half-space such that the remainder of the constraints is not affected.

The above reasoning is not restricted to the complement of one convex polytopic region and can be generalized for the representation of the interdicted region as a finite union of polytopes, see, e.g., [START_REF] Stoican | Hyperplane arrangements in mixed-integer programming techniques. Collision avoidance application with zonotopic sets[END_REF]. Explicitly, in such a case the non-convex region is characterized by:

-h i l x ≤ -k i l + M α i l , ∀i l ∈ I l , (2.9a
) .9b) with I in (2.3) redefined as a union of discrete intervals6 :

i l ∈I l α i l ≤ #I l -1. ( 2
I = l I l .
Consider a collection of obstacles:

P = N j=1 P j
Briefly, gathering the collection of associated support hyperplanes defined as (2.3) we reach the hyperplane arrangement (2.11). Labeling the feasible cells (2.10) into interdicted Σ P = {σ : A(σ) ∩ P = ∅} or allowed Σ X\P = {σ : A(σ) ∩ P = ∅} one can obtain a mixedinteger characterization of the feasible domain:

h i x ≤ k i + M (1 -α i ), (2.12a
)

-h i x ≤ -k i + M α i , (2.12b
)

σ l (i)= ′ + ′ (1 -α i ) + σ l (i)= ′ -′ α i > 0, ∀σ l ∈ Σ P (2.12c)
Even though the construction presented above is a generic one, the binary part of the representation is substantially large as we have a binary variable for each region (3.16b) 7 .

The representation of the interdicted domains is not limited to the polyhedral sets in the wide literature.

O = {(x, y) : (x -x o ) 2 + (y -y o ) 2 ≤ R 2 o } (2.13)
In order to incorporate the corresponding avoidance constraints in the previous polyhedral setting, the obstacles as in (2.13) are approximated with polygons (polyhedra in R 2 , having p vertices). Each interdicted zone can be overapproximated by a set of K p inequalities:

Õ = {(x, y) : (x -x o ) sin 2πm K p + (y -y o ) cos 2πm K p ≤ R o , ∀m = 1, . . . , K p } (2.14)
Using the same reasoning as in (2.7) and the big-M technique, the avoidance constraints are formulated as:

(x -x o ) sin 2πm K p + (y -y o ) cos 2πm K p ≥R o -M β m , ∀m (2.15a) Kp m=1 β m ≤ K p -1. (2.15b)
It is worth to mention that restrictions within the state space of dynamical systems can be directly treated in the modeling stage by including integer auxiliary variables.. In particular is the case of the hybrid systems where switching (a disjunctive type of selection in between modes of functioning) include binary variables or discrete alternatives. The Piece-Wise Affine (PWA) models represent for example an integrated approach, where non-convex domains and local dynamics find a unitary representation. This unitary modeling needs to be differentiated from the classical motion planning framework where the

Y i1 ∨ Y i2 ↔ Y i2 = ¬Y i1
where g i (x) = x 1 ± x 2 ± 6 and we associate: Y i1 = true with R - i (the below half-space) and Y i2 = true with R + i (the above half-space). The interdicted and the allowed region (both illustrated in Figure 2.4) can be described using the following boolean function:

-convex region: h(Y ) = ¬Y 11 ∧ Y 21 ∧ ¬Y 31 ∧ Y 41 -non-convex region: h(Y ) = Y 11 ∨ ¬Y 21 ∨ Y 31 ∨ ¬Y 41
The MIP formulations were briefly sketched using basic set-theoretic notions. For a more extensive presentation the interested reader is referred to the well-known works about polyhedral and hyperplane arrangement notions [START_REF] Ziegler | [END_REF][START_REF] Kuhn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF]] or to the more recent monographies, as [Ionela [START_REF] Prodan | Mixed-integer representations in control design: Mathematical foundations and applications[END_REF]], which offer a more detailed mathematical descriptions.

Optimization-based control framework.

Prerequisites

Ahead of proceeding with typical MIP implementations in motion planning, some clarifications are in order. For motion planning there exists numerous applications that will be treated in detail in Section C. An important aspect is that all those applications have a dynamical behavior of a point of mass in 2D or 3D spaces governed by dynamical mathematical models (either linear or nonlinear).

Up front, let us consider a generic model of the controlled dynamical system designated by an ordinary differential equation (ODE): .16) where x(t) ∈ R nx denotes the state vector, u(t) ∈ R nu the control input vector, and

ẋ(t) = f (x(t), u(t), w(t)) ( 2
w(t) ∈ R nw the disturbance 8 . The mapping f (•, •, •) : R nx × R nu × R nw → R nx
is a continuous function admitting an equilibrium point (i.e., f (x, ū, 0) = 0; without any loss of generality we may assume f (0, 0, 0) = 0).

Remark 2.5. The model (2.16) describes a dynamic in continuous time. As will be detailed below, a considerable part of the literature is based on discrete time dynamics. Hence, we need to adopt, concurrently, keeping similar notations, the discrete time counterpart: .17) The correlation between (2.16) and (2.17) is done using one of the various discretization techniques [START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF].

x(k + 1) = f (x(k), u(k), w(k)) ( 2
♦
In what follows, a relevant feature of the previously-mentioned models is the way of selecting the input variables and consequently the decision domain. In the literature we can distinguish these choices with respect to the linearity/ non-linearity of the models. Thus, for linear models the references concentrate towards the models in which the state is composed by the position and the speed of the agent and the input is represented by the acceleration of the agent, see, for example: [Papen et al. 2017;[START_REF] Yu | Energy-Efficient Trajectory Planning for a Mobile Agent by Using a Two-Stage Decomposition Approach[END_REF]Z. Liu et al. 2017;[START_REF] Bellingham | [END_REF]A. Richards and J. P How 2002]. This type of model implies actuation on each direction and is usually based on the kinematics equations describing the behavior of the controlled systems. Regarding the nonlinear models, the favored variant consists of models, in which the components of the state are the position and heading angle, and the input is represented by speed and steering angle,e.g., [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF]. A usual approach is to work on trajectory tracking which hides the nonlinearities induced by obstacle avoidance but moves the difficulty in the trajectory generation step.

From the control perspective, the earliest methodologies have been based either on optimal control or on non-linear programming, but the controlled systems have operated in obstacle-free environment. Recently, this type of restricted environment is no longer an appropriate assumption from practical reasons (especially, in a multi-agent context).

Hence, the control community has had to propose novel techniques or to adapt the old ones. First adjustments consisted in restricting the controlled system to track predetermined trajectories around the obstacles, but the the resulting control strategies are relevant only for particular dynamics and environments [LaValle 2006].

As a result of the collision avoidance constraints, MIP framework is frequently encountered within optimization-based control

While the focus of this chapter is MIP and its implications in motion planning, ultimately, the obtained formulation comes as the result of optimization-based control strategy. Thus, in order to provide a general overview we briefly recall in basic terms the finite-time optimal control formulation and a series of concepts related to MPC (Model Predictive Control).

Essentially, the optimal control problem [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF][START_REF] Diehl | Lecture notes on optimal control and estimation[END_REF]] consists in finding a control input u ⋆ ∈ R nu such that the system (2.16) converge towards a trajectory x ⋆ (t) ∈ R nx , all by minimizing a specific criterion: min

x(t),u(t) J = min x(t),u(t) H (x(t f )) + t f t 0 L (x(t), u(t)) dt (2.18a) s.t. ẋ(t) -f (x(t), u(t), w(t)) = 0 (2.18b) x(t 0 ) -x 0 = 0 (2.18c) g (x(t), u(t)) ≤ 0 (2.18d)
where t ∈ [t 0 , t f ], (t ft 0 ) denotes the time horizon length, x 0 the initial state, and

g(•,
•) is a vector-function incorporating the physical constraints (e.g., terminal or stage constraints of the controlled system). The OCP (optimal control problem) (2.18) can take various forms according to: i) the choice of functions H(•) and L(•, •), e.g. minimum time problem or terminal control problem [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF]];

ii) the selection of the horizon ([t 0 , t f ]), e.g., [START_REF] Garg | Direct trajectory optimization and costate estimation of finite-horizon and infinite-horizon optimal control problems using a Radau pseudospectral method[END_REF]].

Model predictive control

The main drawback of the generic optimal control problem (2.18) comes from the fact that it is a functional optimization problem, i.e., it involves calculus of variations instead of ordinary calculus for selecting the optimal input function. In order to mitigate this shortcoming, MPC may represent an interesting practical solution. It restrains the set of possible inputs such that the optimization is implemented over a finite set of "decision variables" (by optimising over a finite horizon). Moreover, MPC avoids solving a feedback problem but rather a sequence of finite-time open-loop problems.The feedback is used indirectly (the most recent measurements providing the initial condition for the next openloop optimization) [START_REF] Maciejowski | Predictive control: with constraints[END_REF]] through a receding horizon optimization policy.

The main elements that need to be designed whenever a MPC technique is employed are:

the cost function/ the optimization criterion;

the (internal) model of the system (see Section 2.2.1);

the state and input constraints (and their representation).

Each of these concepts have their importance, but we will concentrate next on the representation of the constraint sets.

By considering that a discretization and linearization part can be handled efficiently at a lower level and that motion planning task can concentrate on the navigation and constraint satisfaction, we consider, in the sequel, agents described by linear dynamics:

x k+1 = Ax k + Bu k , (2.19a
)

y k = Cx k (2.19b)
In the classical MPC studies, the state and input constraints (2.21b) are assumed to be convex containing the origin.

Often, throughout the manuscript, a standard formulation with a quadratic cost is considered for the MPC problem:

J (x k , u k ) =     x k+Np|k 2 P V f (x k ,u k ) + Np-1 l=1 x k+l|k 2 Q + Np-1 l=1 ∆u k+l|k 2 R V(x k ,u k )       (2.20) s.t. x k+l|k = Ax k+l-1|k + Bu k+l-1|k , (2.21a
)

x k+l|k ∈ X , u k+l|k ∈ U , (2.21b)
x k+l|k ∈ C X (P).

( .21c) where N p is the prediction horizon. In addition to these basic ingredients, we considered here additional constraints (2.21c), which need to be efficiently formulated using above presented MIP-techniques. The agent dynamics are described by (2.19) with x ∈ R d the state vector, u ∈ R du the input vector and the matrices A, B of appropriate dimension. Furthermore, in (2.20) -( 2.21) N p is the prediction horizon, the matrices P (terminal cost penalty), Q (output error penalty) and R (control move penalty) are positive semi-definite and of appropriate dimensions. The sets X and U are compact sets from R d and R du , respectively. The notation x k+l|k stands for the value of x at time instant k + l, predicted upon the information available at time k ∈ N.

Remark 2.6. For further implications, the value of P from the terminal cost V f (x k , u k ) is selected such that the function V (x) = x ⊤ P x is a Lyapunov function for a prestabilizing state-feedback law u = Kx with K ∈ R d×du in a vicinity of xi [Mayne et al. 2000]. As well, the values of Q and R from the cost per stage V(x k , u k ) are chosen in order to enhance the tracking performances.

Standard MIP-based problems in motion planning

The motion planning problems [LaValle 2006] are encountered in several enginnering fields, from robotics to automotive and aerospace and have been studied in time optimal control, nonlinear control, stabilization, reachability and other related frameworks. To review all these topics and their particularities is out of the scope of this manuscript, but we aim to retain those of them in which the presence of alternatives or discrete decisions making is explicitly considered. The emphasis is placed on the specific applications where non-convex feasible domains appear and are subsequently encoded by mixed-integer techniques.

Task assignment

Obstacle In the schematic view of Figure 2.5, the motion planning sub-areas, which can be efficiently formulated through a MIP approach, can be delineated: i) Task assignment (TA) is the strategic decision of allocating one objective to a particular subsystem (who-goes-where?), these objectives can be interchanged and their suitability is measured by a specific criterion. TA is a discrete decision-making process where the number of alternatives within the pairs (resource objective) association is countable.

ii) (Collision-free) Path planning is the construction of a route in the position space without an explicit parametrization in time and without explicitly considering the dynamics. The collisions in this framework bring the alternatives into the design and open the way to a mixed integer formulation.

iii) (Collision-free) Trajectory planning represents the problem of constructing a function which associates a time interval to a path. It takes into account the particularities of the agent's dynamics and the resulting function is generated in an open-loop manner. The time-parametrization and the dynamical constraints remain a continuum problem but the preservation of the collision-free attributes of the solution inherits the mixed integer structure.

iv) (Online) Obstacle and collision avoidance is the problem of finding input control signals which minimize a performance criterion as in (18a) while simultaneously avoiding collisions within closed-loop strategy. This real-time task is intended to deal with the uncertainties based on the feedback mechanism all by integrating the collision avoidance constraints.

We have graphically illustrated in Figure 2.5 the functional relationship among the items described above. A noteworthy aspect is that for the task assignment and path planning sub-areas the particularities of the dynamics have, generally, almost no influence in problem statement, in contrast with the other sub-areas. Besides, TA and path planning are not tackled in most of the motion planning literature as stand-alone topics, but together with or as part of the others sub areas. As depicted in Figure 2.5 (the dotted lines), path planning is usually implicitly included in trajectory planning and TA is either a part of obstacle avoidance or included in path planning.

As well, collision-free trajectory planning and collision avoidance represent similar approaches, the difference between them being the manner of interaction with navigation and control. While the obstacle avoidance intervenes into the navigation-control loop, the trajectory planning only generates the reference (to be tracked) for this loop.

A considerable part of the control literature includes also the formation control in motion planning tasks, although it represents, in fact, a control level which gathers all the above sub-areas of motion planning, and (as, e.g., TA) makes sense only in a multi-agent context.

For the sake of completeness, we treat this topic in Section 2. 3.5, where connectivity maintenance problems in multi-agent systems are presented.

In Table 2.1 we depicted the classification of the references with respect to the motion planning problem predominantly treated therein. As well, we included those reference addressing formation control, connectivity maintenance and corner cutting avoidance problems (topics covered in Section 2. 3.5).

Remark 2.7. Before detailing the motion planning problems, let us highlight the differences between the path and trajectory planning with respect to their time-parametrization (see, e.g., [START_REF] Beard | Small unmanned aircraft: Theory and practice[END_REF] or [L. Yang et al. 2014] for the detailed studies). The choice between these two methodologies (within the low level feedback control scheme) depends on particularities of the application, e.g. for a fixed wing UAV, a trajectory planning methodology can have some undesired consequences [START_REF] Beard | Small unmanned aircraft: Theory and practice[END_REF]. These particularities are recalled and discussed in Section C. ♦ Remark 2.8. Within the motion planning and control, a pre-condition of the use ofMIP approaches is the existence of an awareness map which entails the global knowledge about the environment and the strategic move within a predefined goal and a local/global information on the uncertainties and conflicts. Thus, representation of the environment has a significant impact on the performances of the motion planning strategies (except the task assignment problem where the impact is often associated to economic and performance considerations). ♦

Task assignment

Let us present a standard MIP formulation for the task assignment problem, emphasizing its importance in the context of motion planning for a multi-agent system.

Task assignment refers to a strategic decision of allocating one objective to a particular subsystem (it responds to the question "who-goes-where?").

Consider the undesired occurrence of an earthquake, which affects 10 buildings. An emergency committee must assign for each building a rescue team, in order to minimize the intervention time, taking into account the number of rescue teams and their necessary time to arrive at respective buildings. This example can be easily reformulated for the case of target tracking in a multi-agent system. At large, a task assignment algorithm must provide an assignment which optimizes a specific criterion.

As a generalization, let us consider N agents and M targets. In order to minimize the amount of time to reach all the targets, we should optimally assign agents to targets. To this end, we define a binary variable x ij , which equals one if the agent i visits the target j and zero otherwise. Also, we consider a cost associated to each combination agent-target:

c ij .
We take into account that each agent must attain at least one target, and, at the same time, each target must be reached by an agent (constraints (2.22b) and ( 2.22c), respectively). Thus, we obtain a standard MIP formulation of task assignment problem, described in (2.22).

min

x ij N i=1 M j=1 c ij x ij (2.22a) s.t. N i=1 x ij ≥ 1, ∀j (2.22b) M j=1 x ij ≥ 1, ∀i (2.22c
)

x ij ∈ {0, 1}, ∀i, j (2.22d)
Since we have an integer requirement on the variables x ij , the problem (2.22) belongs to the class of ILP (integer linear programming) (a subclass of MIP, where all the variables are integer.) A noteworthy aspect is given by a special characteristic of (2.22): solving its LP relaxation9 allows obtaining the solution of ILP original problem.

In specific situations, LP relaxations may lead to an integer solution, as in Figure 2.6 methods, e.g., in [Alighanbari, Kuwata, and J. P. How 2003] a comparison between the Tabu search and the MIP formulation is done, and is concluded that for higher dimensions the Tabu search exhibits better performance. However, the tangible theoretical advantage of the MIP formulation (2.22) is that it guarantees the global optimum. This has a substantial importance in some critical situations.

Remark 2.9. In a similar manner as the task assignment, the notorious travelling-salesman problem (TSP) appears first time in [START_REF] Dantzig | Solution of a Large-Scale Traveling-Salesman Problem[END_REF] posed as MILP. For the classical TSP, the formulation ( 19) is used by converting the constraints (19b) and (19c) into equalities and imposing that N = 1 and M ≥ 2. ♦

Besides task assignment, a plethora of resource allocation problems can be readily formulated as MIP. Many of these problems are used in several applications within related motion planning domains such as aircraft maintenance (or repair) [START_REF] Bajestani | Scheduling a Dynamic Aircraft Repair Shop with Limited Repair Resources[END_REF], crew scheduling(and flight retiming) [START_REF] Mercier | An integrated aircraft routing, crew scheduling and flight retiming model[END_REF] etc.

There exist approaches building on the idea of state space (workspace) organization in terms of available resources. More precisely, these approaches partition the environment in equally-shaped cells and the cells are viewed as a shared resource. In this manner, the collision and obstacle avoidance problems are transformed in resource allocation with mutual exclusion.

The collision avoidance problems can be formulated as a resource allocation problem with mutual exclusion.

In the line of developments, in [X. Wang et al. 2015] the collision avoidance constraints are given by the fact that each cell cannot be visited by two distinct agents at the same time, i.e. the agents do not use the same resource simultaneously. The formulation is an extension of the task assignment problem (2.22), the novelty comes from the evaluation and resolution of the MIP in an online manner.

Another reference using this "resource allocation" approach is [Haghighi, Sadraddini, and Belta 2016] where a robotic swarm has the objective of respecting specific and complex patterns in a bidimensional workspace. These patterns are described using the spatiotemporal logic (SpaTeL). SpaTeL consists in propositions indicating the number of agents which can access a cell at a certain time. These SpaTeL formulas are converted in mixedinteger linear constraints, resulting a formulation comparable with the task assignment (2.22) but more complex. A similar approach is employed in [Z. Liu et al. 2017], but the difference is given by the inclusion of the communication constraints in the motion planning problem under the form of signal temporal logic (STL).

Path planning

As mentioned in Section 2.3, a path planning strategy involves the generation of a route without an explicit parametrization 10 . Further, this route is called the planned path.

Path planning refers to the construction of a route in the position space without an explicit parametrization in time and without explicitly considering the system dynamics.

Definition 2.2. A path between two points x 0 and x f within the navigation space is defined by a map γ :

[0, 1] → X with γ(0) = x 0 , γ(1) = x f (2.24) ♦
There are in the literature two major classes of methods to describe the function γ. The first one, the explicit methodology, consists in providing the value of γ for the entire interval [0, 1] [Janeček, Klaučo, Kalúz, et al. 2017]. The second class of methods is the most encountered alternative, the function takes its values at a finite number of elements in the domain γ ∈ [0, 1] from a given set of, so-called, waypoints [R. J. M. Afonso, Galvão, and K. H. Kienitz 2013].

Remark 2. [START_REF]An iterative solution for the MINLP[END_REF]. Despite the lack of explicit time parametrization, it should be forseen that the waypoints represent either the nodes of graph in a dynamic-free setting, or the boundary conditions for intervals of time. ♦

For further use, we define the set of all waypoints describing the planned path: W = {x w 0 . . . x w Nw } and consider an additional condition to Def. 2.2:

∀x w k , ∃! θ k ∈ [0, 1] such that γ(θ k ) = x w k .
The methodologies involved in the selection and determination of the waypoints are outside the scope of this paper, although they are evoked in Section A. In the path planning problem, MIP is not employed in the generation of the path, but rather to address the problem of distribution of the waypoints (i.e., their ordering along the path), satisfying some criteria. In other words, MIP plays an instrumental role for optimizing a path obtained with one of the well-established techniques from the robotic field (e.g., sampled based methods -see Section A ). The criterion of this optimization procedure has to include/cover costs given by the physical and "economical11 " limitation.

One of the most used approaches treating this type of problem is based on a particular case of the task assignment (2.22), the travelling-salesman problem (TSP). Under this form we can identify the MIP usage in path planning and a classical example in this sense is to be found in [A. [START_REF] Richards | Aircraft trajectory planning with collision avoidance using mixed integer linear programming[END_REF][START_REF] Schouwenaars | Safe trajectory planning of autonomous vehicles[END_REF]]. In this prototype problem, they deal with a robot which should visit N w waypoints while minimizing the cost of the maneuver. Hence, the MI constraints which enforce the visiting of the considered waypoints are the following:

∀i ∈ {1, . . . , T }, ∀k ∈ {1, . . . , N w }, (2.25a
)

x i -x w k ≤ M (1 -b ik ), (2.25b)
T i=1 b ik = 1, ∀k and b ik ∈ {0, 1}, (2.25c)
where T is the number of time steps, and b ik is a binary variable that indicates whether or not the waypoint k was visited at the time step i. The constraint (2.25c) ensures that each of the waypoints is visited once by the agent. The ordering of the waypoints along the interval depends on the selected criterion (e.g., the minimum time to visit all the waypoints). In practice, the attainability of the waypoints is affected by various factor and, thus, in most of the situations, reaching a vicinity of the waypoint is a reasonable objective. In such cases, (2.25) needs to be adapted,e.g., the constraint (2.25b) is replaced by

H k (x i -x w k ) ≤ h k + M (1 -b ik )
, where H k and h k are given by the supporting hyperplanes describing the considered vicinity region.

Remark 2.11. In the majority of sampled-based methods (which create a meshed graph of points by sampling of the positioning space) [LaValle 2006], the problem of finding the shortest path is tackled using one of the various algorithms (Dijkstra, A ⋆ are common choices). However, in [START_REF] Taccari | Integer programming formulations for the elementary shortest path problem[END_REF]] several MIP formulations for the elementary shortest path problem are introduced as extensions of TSP. ♦

Besides the approach relying on (2.25), there exist a few others works, which treat MIPbased path planning in different manners. For instance, in [START_REF] Vitus | Tunnel-milp: Path planning with sequential convex polytopes[END_REF]], a Tunnel-MILP approach is employed. This algorithm divides the global motion planning problem in three main tasks. The first task consists in finding a path as in Def. 2.2, while ignoring the vehicle's dynamics constraints. Next, the path is used along with a convex decomposition of the space to generate a sequence of N R convex polytopes from the start to the goal:

F i = {x ∈ R 2 : A i x ≤ b i }, ∀i ∈ {1, . . . , N R }.
In the third task, this sequence (a tunnel of polytopes) is used to constrain the position of the agent p(t) to an optimal and dynamically feasible route from the initial point to the goal inside the tunnel. More precisely, a MILP problem is employed, forcing the agent to remain in one of the regions of the tunnel at all times. This leads to some OR-constraints 12 , similar with the GDP formulation (2.1), which can be readily rephrased using MI techniques: .26) where p(t) is the position of the considered agent.

i∈{1,...,N R } A i p(t) ≤ b i , ∀t (2
The same idea is used in [R. [START_REF] Afonso | Waypoint trajectory planning in the presence of obstacles with a tunnel-MILP approach[END_REF], where the procedure to find the path relies on generalized Voronoi diagrams. Here, a number of points is considered: on the facets of the obstacles (modeled as convex bounded sets) and on the boundaries of the navigation environment in order to obtain a Voronoi graph. Subsequently, the nodes of the graph that are inside the prohibited region are removed along with the edges that connected these nodes to the other ones in the graph. After that, using a Delaunay triangulation the space is divided into triangles whose intersection is empty, except for their sides. The triangles that are crossed by the path are computed and merged to form convex polytopes, attempting to have the minimal number of convex regions along the tunnel. These polytopes are used further to impose the obstacle avoidance constraints as in the tunnel-MILP formulation (2.26).

In a similar approach, [R. [START_REF] Deits | Efficient mixed-integer planning for UAVs in cluttered environments[END_REF] treats the problem of path planning in a multi-obstacle environment for a rotary-wing UAV. The underlying idea is based on a mixed-integer optimization which assigns polynomial trajectories to convex regions (known to be obstacle free). The paths are defined as piece-wise polynomial functions in time with vector-valued coefficients. Having an a priori (offline) chosen degree of polynomials and number of pieces, the optimization problem returns the coefficients of each polynomial which ensure a collision-free trajectory. This polynomial parametrization is possible because the model of the considered system 13 is characterized by differential flatness. Regarding the obstacle-free convex regions, they are obtained offline using IRIS (Iterative Regional Inflation by Semidefinite programming [Robin Deits and Russ Tedrake 2015]), a technique for greedy convex segmentation of the free space. The proposed approach is tested for various simulations scenarios with satisfactory results, although sometimes the convex segmentation does not cover the entire obstacle-free space. Moreover, a comparison with the classical 14 MIP method is done, this alternative being characterized through a more complex integer program. An interesting aspect of this work is the Interpreted as XOR.

The nano-quadcopter Crazyflie 2.0 From the authors point of view: use the obstacle faces ("hyperplanes") directly in the generation of safe region dependence of the generated paths on the considered convex safe region segmenting the space.

Besides the above works,we can cite references which treat, but only marginally, the problem of path planning. For example, in [Janeček, Klaučo, Kalúz, et al. 2017] or in [Janeček, Klaučo, and Kvasnica 2017] a toolbox is described which includes various helpers to generate a reference path. Thus, a circular trajectory of a known radius or a polygonic reference passing through the given set of waypoints are provided, but the use of MIP is kept to a minimum.

Trajectory planning

The trajectory planning pertains to the problem of determining both the path and how to move along it. Thus, a trajectory planning strategy returns a path which is explicitly parametrized in time. In this section we consider first the works dealing with general problems and, subsequently, we focus on classical applications in the field of autonomous vehicles (i.e., autonomous overtaking, merging junctions or lane changing).

Trajectory planning refers to the construction of a function generated in openloop, which associates a time to each point on the path and takes into account the particularities of the system dynamics.

Definition 2.3.

A trajectory between two points x 0 , x f from the navigation space is given by the continuous function γ

: [t 0 , t f ] → X with γ(t 0 ) = x 0 , γ(t f ) = x f . (2.27)

♦

Recalling the path planning task, the specification of the continuous function describing the trajectory can be done in two different ways: explicitly, for the entire interval 0, 1}, and by using a discrete set of values (waypoints). For the second one, the set of all waypoints is modified as follows: W = {(x w 0 , t w 0 ) . . . (x w Nw , t w Nw )} and an additional condition to the Def. 2.3 is considered: γ(t w k ) = x w k , ∀x w k . Nevertheless, the waypoint method is rarely used in MIP-based trajectory planning (and in the related literature). Hence, whenever it is not specified otherwise, the trajectory is given explicitly as the continuous function γ(t).

Trajectory planning in various tasks

One of the popular references treating the problem of trajectory planning is [M. [START_REF] Earl | [END_REF] which introduces an iterative MILP algorithm. More precisely, considering circular obstacles (see Section 2. 1.2, eq. (2.13)) and a traditional nonholonomic (car-like) vehicle, the authors propose an algorithm which guarantees obstacle avoidance over the entire trajectory and distributes avoidance times efficiently, resulting in smaller MILP formulations.

The proposed algorithm needs the discretization in time of the continuous nonholonomic model. The non-uniform discretizations are exploited in order to bypass an extensive computational effort corresponding to the resolution of large MILPs. Support for nonuniform discretizations in time allows the use of intelligent time step selection algorithms for the generation of more efficient MILP formulations. Thus, the idea behind of the iterative MILP algorithm is to find a convenient distribution of the avoidance times, and circumventing the growing of the specific MIP problem complexity (induced by the enumeration of the constraints at each sampling time). Therein, two different problems are considered: trajectory-generation with obstacle-avoidance requirements and minimumtime trajectory-generation problems.

Iterative MILP algorithms address the issues of MILP coping with large-scale models.

The work of [M. G. Earl and R. D'Andrea 2002] studies the cooperative control of multivehicle systems. Building on the requirements of a robotic competition, the authors model these simplified competition rules using a hybrid system. Further, the control problem is stated in an optimization framework, using MILP. The considered simplified competition involves two teams of robots, the attackers and the defenders, on a playing field with a region at its center called the defense zone. The attackers are drones directed toward the defense zone. The objective for the defenders is to block the attackers from entering the defense zone by intercepting each attacker before it enters the zone. Once an attacker enters the defense zone or is intercepted by a defender it remains stationary for the rest of the competition. While pursuing its objective defenders must avoid collisions with other defenders and obstacles as well as avoid entering the defense zone which is off limits to defending robots. The control strategy is implemented with a centralized controller with perfect knowledge of the system, perfect access to all states, and with the ability to transmit control signals to the defenders instantaneously. The controller needs to figure out the inputs to provide each defending robot so that the objective is achieved. Using MPC, they obtained a set of control inputs that minimize the number of attackers that enter the defense zone over the duration of the drill and, in addition, is consistent with the system dynamics (robot dynamics) and the constraints (no collisions, etc.). The obtained model is a system composed of continuous and discrete states (a hybrid system) with linear dynamics subject to inequality constraints and logical rules. In other words they obtained a MLD (mixed-logical dynamical) system [Alberto Bemporad and Morari 1999].

The resulting optimization problem is a MILP which can be readily solved using one of the state-of-the-art solvers (e.g., ILOG CPLEX or GUROBI).

Multi-robot systems can be modeled as hybrid (MLD) systems using MIP techniques. The strategy for one team is fixed and modeled by state machines, while for the other the behavior is controlled using the iterative MILP method. On a side note, the approach therein was developed independently from a similar one introduced in [A. [START_REF] Richards | Aircraft trajectory planning with collision avoidance using mixed integer linear programming[END_REF].

The trajectory planning topic is treated in [START_REF] Cetin | Hybrid mixed-logical linear programming algorithm for collision-free optimal path planning[END_REF] through the generation of collision-free trajectories for the reconfiguration of spacecraft formations, aiming also for an optimal fuel consumption. In order to model the spacecraft and their corresponding safety regions, unrotated cubes are used. Moreover, the trajectories to be followed are discretized in time using cubic splines and, thus, the generic problem is translated to an optimization problem. The resolution led to a trajectory parameterized by the spacecraft positions and velocities at a set of waypoints. To this end, the "big-M" technique is used to write the parametrized optimization as a MILP whose solution can be obtained either using standard MILP solvers (see Section C) or using the concept of a sequential linear programming. These two alternatives are compared on two standard validation tests whose aim is to swap the position within a spacecraft with minimum fuel consumption. The comparison leads to an ample discussion on the feasibility of MILP and on the methods necessary to shorten the resolution time.

In [A. [START_REF] Richards | Performance evaluation of rendezvous using model predictive control[END_REF]] MILP has been used for open-loop vehicle trajectory design, enabling the inclusion of non-convex constraints such as plume impingement avoidance. Recent advances in the field of autonomous vehicles bring to light the issues like vehicle lane change and overtaking on highways in the assisted driving framework. In order to perform such maneuvers, it is fundamental to compute suitable and comfortable trajectories that take into account the vehicle limitations as well as safety restrictions.

Trajectory planning in autonomous vehicles tasks

Safety restriction and vehicle limitations require efficient MIP formulations.

In [Molinari, Anh, and Re 2017] an efficient MIP formulation for the autonomous overtaking problem is introduced. The considered vehicles have a simple dynamic model, the non-convex feasible region being represented using hyperplanes arrangement similar with [Ionela [START_REF] Prodan | Mixed-integer representations in control design: Mathematical foundations and applications[END_REF]](as in Section 2.1). A complete formulation for trajectory generation with collision avoidance guarantees is presented for the case of an agent surrounded by a number of agents. The MPC criterion contains a desired reference state, which ensures that the overtaking takes place. They also considered two methods for reducing the numbers of binary variables: logarithmic formulation and cell merging [START_REF] Stoican | Hyperplane arrangements in mixed-integer programming techniques. Collision avoidance application with zonotopic sets[END_REF]. The illustrative example therein and a comparative analysis between the two reduction methods from complexity point of view, are validated in IPG CarMaker (an accurate vehicle simulator).

In [Ballesteros-Tolosana et al. 2017], the lane change and overtaking maneuvers aim to generate trajectories able to ensure the comfort of passengers. Mathematically, the problem is formulated as an Optimal Control Problem (OCP) due to its handling kinematic and collision avoidance constraints (the last ones in terms of hyperplane arrangements as in Section 2.1). In order to mitigate the drawback of the substantial computational effort associated with that kind of representation, a pre-analysis step is presented therein, a step consisting in the enumeration of the all possible overtaking configurations and of their resulting compact MIP formulation. The so-obtained non-linear constrained optimal control problem is solved using a multiple-shooting approach which leads to improvements in the computational burden when compared to the pre-analysis step.

The presence of MIP in merging junctions problem can be found in a plethora of works.

For instance, [START_REF] Fayazi | Optimal scheduling of autonomous vehicle arrivals at intelligent intersections via MILP[END_REF] deal with the optimal scheduling of autonomous vehicle arrivals at intersection, eliminating the need of traffic signals. The idea is to design an intersection controller able to coordinate the flow of vehicles through intersection, scheduling the intersection access in safe conditions and receiving information from all subscribing vehicles. Moreover, each vehicle informs the controller about its movement and desired schedule. Therefore, the optimization-based controller should find the optimal sequence of vehicles crossing the intersection and their corresponding time of intersection access, minimizing the difference between the current time and the expected arrival time of the last vehicle passing the intersection. The resulting optimization problem has to take into account physical constraints (speed limit and maximum acceleration) and a safety "window" between two consecutive vehicles accessing the intersection. These constraints lead to disjunctions and, respectively if-then statements, which are modeled through MIP (big-M formulation).

Likewise, [X. [START_REF] Huang | Speed trajectory planning at signalized intersections using sequential convex optimization[END_REF] tackled the problem of speed trajectory planning at signalized intersection. The idea is to optimize the vehicle speed trajectory over multiple intersections in order to minimize the fuel consumption and the travel time. In order to obtain a mathematical formulation of the problem , the authors assume that the signal traffic state is known and the influence of a lead vehicle is ignored. The main advantage of this formulation (with respect to, e.g., [START_REF] Fayazi | Optimal scheduling of autonomous vehicle arrivals at intelligent intersections via MILP[END_REF]) is the additional turnings at intersections (including a turning speed constraint). After a comprehensive description of the fuel consumption model, the effect of turning and acceleration model, a MIP formulation is stated. The presence of integer variables models the crossing of the intersection without violating red light, indicating the active green phase window.

Also, [Y. [START_REF] Zhang | Modelling and traffic signal control of a heterogeneous traffic network with signalized and non-signalized intersections[END_REF]] presents a control strategy for a heterogeneous traffic network. Heterogeneity comes from the presence of both signalized and non-signalized15 intersections. After the validation of the proposed model of the heterogeneous traffic network, a control strategy is stated using a MIP formulation. Moreover, [Y. [START_REF] Zhang | Traffic light scheduling for pedestrians and vehicles[END_REF]] proposes a traffic signal scheduling strategy which takes into account both vehicles and pedestrians presence. After developing a mathematical model for the pedestrians flow and for the vehicle traffic network (more than one intersection/junction), the problem of scheduling is stated as a MIQP (mixed-integer quadratic program).

Obstacle and collision avoidance

Obstacle and collision avoidance is the problem of finding input control signals which minimize a performance criterion while simultaneously avoiding collisions through a closed-loop strategy.

An adequate collision and obstacle avoidance strategy is critical in order to ensure the safety and the integrality of both the system and the environment. It is generally accepted that the collision avoidance problem is a challenging task due to the presence of the non-convex feasible domain. Describing this non-convex domain has computational and structural implications, generally leading to a trade-off between computational efficiency and control performance. Since MIP is able to explicitly model the non-convexity as in Section 2.1, it has become over time a proper approach for expressing collision avoidance problem.

The MIP approaches in feedback control with obstacle and collision avoidance problems can be distinguished in two main categories. The first class of MIP approaches is widespread in non-MIP (robot) motion planning literature and consists in choosing a partitioning method which does not depend on the particularities of the environment (e.g the obstacles' form). A standard partitioning is the one formed by equally-sized square cells [Haghighi, Sadraddini, and Belta 2016;X. Wang et al. 2015]. Using this approach, the collision-avoidance problem becomes, in most of the cases, a resource allocation problem, as presented in the Section 2. 3.1. The second MIP approach or the geometric approach relies on set theory and on the capability of mixed-integer techniques to efficiently encode the description of non-convex sets (see Section 2.1). Even though different, both approaches share a common characteristic: the workspace is decomposed in cells and subsequently the MIP is employed in representing disjunctive constraints.

Remark 2.12. In the literature there are works which make a distinction between obstacle and collision avoidance. More precisely, the collision avoidance refers only to the mutual collision avoidance within multi-agent systems, and obstacle avoidance may refer to the collision avoidance either with stationary obstacles or with moving obstacles which have an a priori known trajectory, in most of the cases a periodic one. As well, some others references consider that obstacle avoidance refers to the static environment, while collision avoidance covers any kind of moving obstacles. ♦

Static obstacles

The work of [START_REF] Culligan | Online trajectory planning for UAVs using mixed integer linear programming[END_REF]] presents a path planner using MILP to solve a receding horizon optimization problem for unmanned aerial vehicles (UAVs). The MILP formulation has two important components: the hard constraints for obstacle and multi-vehicle avoidance and an approximation of vehicle dynamics with extensions to 3D. Moreover, some enhancements of the MILP framework are discussed in order to provide decreased resolution time and also increase the capability of the path planner. These improvements consist of various techniques like a variable time step size, linear interpolation points, and horizon minimization. A noteworthy aspect, the concept of variable time steps is extended to the receding horizon, non-iterative MILP formulation. Variable time step sizing allows the simulation horizon time to be lengthened without increasing resolution time dramatically. Horizon minimization decreases resolution time by removing unnecessary obstacle constraints from the the problem (similar with [Janeček, Klaučo, Kalúz, et al. 2017]).

In [START_REF] Schouwenaars | Safe trajectory planning of autonomous vehicles[END_REF]] the author presents a framework for safe online trajectory planning of unmanned vehicles through partially unknown environments. A MPC framework is employed, MILP playing an instrumental role in order to incorporate the collision avoidance constraints (similar with the techniques in Section 2.1). An interesting aspect is that the agents can be controlled either through a standard (velocity) control system, or by using a maneuver scheduler that allows the implementation of a maneuver from a discrete set of possible ones. This hybrid control architecture is applied and enhanced for a particular type of dynamics, corresponding to a small-scale helicopter. For feasibility gurantees a beneficial concept was considered: terminal feasible invariant set, a set in which an agent can remain for an indefinite period of time with anti-collision guarantees. Effectively, these sets are computed online, being represented as affine constraints on the last stage of the prediction horizon. Via these sets, it is provided an a priori known backup plan that is dynamically feasible and obstacle-free and, thus, the feasibility and safety can be guaranteed. The proposed strategy was tested on an unmanned Boeing aircraft using scalable loiter circles as feasible invariant sets. From the multi-agent perspective, the control strategy is a distributed one, each agent only computing its own trajectory while accounting for the latest planned behaviour of the agents from its proximity. The potential conflicts are solved in real-time such as to preserve the feasibility guarantees. In order to illustrate the benefits of the considered strategy, the algorithm was run over a scenario involving a fleet of small-scale helicopters which aimed at maintaining wireless connectivity in a cluttered environment.

Obstacle and collision avoidance constraints are usually imposed at the sampling time without regards to the intra-sample behavior of the agent. Hence, it is possible for an agent to "cut the corner" of an obstacle in continuous-time while apparently respecting the constraints in discrete-time. The idea employed in the literature is to consider additional constraints which ensure that the segment between two consecutive positions does not cut the obstacle [START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF]A. G. Richards and Oliver Turnbull 2013]. [Stoican, Ionela Prodan, and Grøtli 2018] provides a treatment of the multi-obstacle case within a hyperplane arrangement setting with exact and over-approximated representations.

Mobile obstacles

The mobile obstacles in the MIP framework are not so common as the static environments due to the associated computational burden. In most of the references treating the timevarying environments, the interdicted mobile regions are modeled as rectangular/cubical exclusion regions [A. [START_REF] Richards | Aircraft trajectory planning with collision avoidance using mixed integer linear programming[END_REF][START_REF] Richards | Trajectory planning[END_REF]. In fact, this representation is a particular case of the polytopic representation (see Section 2.1), but with a limited number of corresponding mixed-integer linear constraints. Thus, the exclusion region constraints in a MIP formulation for two moving obstacles within a 2D environment (x-y coordinates) are given as:

x 1 -x 2 ≥ d -M α d1 (2.28a) x 2 -x 1 ≥ d -M α d2 (2.28b
)

y 1 -y 2 ≥ d -M α d3 (2.28c) y 2 -y 1 ≥ d -M α d4 (2.28d) 4 i=1 α di ≤ 3 (2.28e)
where d is the safety distance (length of the edge of square exclusion region), α di are binary variables, while (2.28e) ensures that at least one of the above constraints is active.

It is worth to mention that the rectangular exclusion zone can be used in order to model other non-convex constraints, e.g. [A. [START_REF] Richards | Trajectory planning[END_REF][START_REF] Culligan | Online trajectory planning for UAVs using mixed integer linear programming[END_REF]] Similar to trajectory planning, the collision avoidance problem can be found/applied in tasks regarding the autonomous vehicles. For instance, in [Mukai, Natori, and Fujita 2017] the problem of merging for vehicles on a motor way is tackled using MIP along with a receding horizon strategy (MPC). Thus, interdicted regions are described using logical statements (AND/OR), which are further modeled as mixed -integer constraints, using the "big-M" method.

Furthermore, [Molinari, Anh, and Re 2017] treat the autonomous overtaking problem using MPC in an efficient MIP formulation. The considered vehicles have a simple dynamic model, the non-convex feasible region being represented as in Section 2.1. The MPC criterion contains a desired reference state, which ensures that the overtaking take place.

As well, [START_REF] Bali | Merging vehicles at junctions using mixed-integer model predictive control[END_REF] propose a method for vehicle merging scenarios in junctions with relative cost prioritization. The method is based on MPC, employing MIQP optimization. The scheme provides optimal control properties while maintaining safety and recursive feasibility. The latter properties are ensured through positive control invariance of simple time headway constraints. For examples with two vehicles, tunable prioritization and gap acceptance are verified and presented on a decision graph. Priorities are then demonstrated to be respected in an example with four vehicles.

Equivalent MLD formulation for obstacle avoidance

The mixed logical dynamical (MLD) systems were introduced for the first time in [Alberto Bemporad and Morari 1999] and represent a relevant framework for modeling and controlling systems which incorporates linear dynamic equations, logic rules, and operating constraints. These are described by linear dynamic equations subject to linear inequalities involving real and integer variables as in (2.29).

x k+1 = Ax k + Bu k + B 2 δ k + B 3 z k (2.29a
)

y k = Cx k + Du k + D 2 δ k + D 3 z k (2.29b) E 2 z k +E 3 δ k ≤ E 1 u k + E 4 x k + E 5 (2.29c)
where x k , y k , u k denote state,output and, respectively, input vectors of the system. The δ k ∈ {0, 1} n δ and z k ∈ R nz are auxiliary logic and continuous variables, respectively.

Let us consider a system described by a single agent which has to move in a non-convex domain. The agent is described by the dynamics16 :

x k+1 = Ax k + Bu k (2.30a
)

y k = Cx k + Du k (2.30b)
We examine a basic example, the agent has to stay outside a polytope P = {x : Sx < r}:

x k+1 / ∈ P.

We aim to describe the considered system via a MLD formalism, more precisely, to rewrite the above condition in the form of (2.29c).

Firstly, it is straightforward that (2.30a) and (2.29a) are equivalent if we consider the matrices B 2 and B 3 null. Next, writing the imposed condition (x k+1 / ∈ P ) in a similar manner as in (2.7a) and replacing x k+1 according to the dynamic, we compactly write: .31) where M is a sufficiently large constant (according to the "big-M" formulation), N is the number of linear constraints describing the polytope

-M I N 1 ⊤ α k ≤ SA 0 x k + SB 0 u k + -r N -1 , ( 2
P = {x ∈ R d | Sx ≤ R, S ∈ R N ×d , R ∈ R N }
, and I N ∈ R N is the identity matrix. By adding constraints on input (S u u k ≤ r u )and considering the null matrix O N ∈ R N , the MLD formulation of the considered system is:

x k+1 = Ax k + Bu k (2.32a)   -M I N 1 ⊤ O N   α k ≤   SA 0 0   x k +   SB 0 -S u   u k +   -r N -1 r u   (2.32b)
An important notion in relationship with a MLD system is the well posedness according to Definition 1 from [Alberto Bemporad and Morari 1999]. Namely, this ensures the existence and unicity of a trajectory in forward time. The formulation (2.32) fulfills these requirements and can be used as a alternative MIP formulation. ♦

In [START_REF] Ritter | Adaptive Observation Strategy for Dispersion Process Estimation Using Cooperating Mobile Sensors[END_REF]] the equivalent MLD formulation is further extended. A multi-vehicle cooperating system is modelled and controlled such that the system reaches some target locations in a reasonable time avoiding any possible collision among the agents. Thus, MLD-system represents a valid compact formulation for modeling multi-vehicle system, because all existing constraints (dynamics, collision avoidance, measurement etc.) are embedded in a complete configuration, which allows to use existing methods to solve the control problem.

MIP reformulation for drift counteraction optimal control

Another form in which MIP is employed in motion planning is within DCOC (drift counteraction optimal control) or optimal exit-time control [START_REF] Zidek | Optimal and receding horizon drift counteraction control: Linear programming approaches[END_REF]. The main objective in such problem is to satisfy prescribed constraints for as long as possible [START_REF] Zidek | Receding horizon drift counteraction and its application to spacecraft attitude control[END_REF]]. In other words, DCOC is a particular optimal control problem, aiming to determine a sequence of control inputs that maximizes the first exit-time from a given set. There exist various applications in motion planning where DCOC represents an useful tool, handling well the systems with finite resources (fuel or energy). In [START_REF] Zidek | Receding horizon drift counteraction and its application to spacecraft attitude control[END_REF]] the MILP formulation of DCOC with application in spacecraft attitude control is fully addressed and for its resolution a LP-based iterative procedure is introduced.

As well, in a similar approach, [Maia and R. [START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF] present an implementation of the shifting prediction horizon for an MPC controller, tackling the obstacle avoidance problem, using binary variables and implicitly MIP. Moreover, [A. [START_REF] Richards | Performance evaluation of rendezvous using model predictive control[END_REF] MILP optimization is used to effect a variable horizon length, leading to guaranteed finite-time completion.

Connectivity maintenance and formation control

In today's complex and various environments, the vast majority of the activities is too difficult and time-consuming to be handled by only one agent/robot/entity. Thus, in order to perform these activities with increased accuracy, redundancy and in a reduced time, cooperative teams of robots/agents may be employed. In this manner, the key elements of risk for the safety and integrity of systems are mitigated at the expense of an extensive augmentation of the systems to be supervised and controlled. Factors like the large scale, geographical distribution, high failure rates and heterogeneity of network systems are becoming decisive in consideration of the suitable control architecture.

An important issue in the control of multi-agent systems is the coordination of a cooperative team of agents (robots) in order to accomplish a given "mission" in an efficient manner. In several cases, these teams of agents are required to converge (and maintain) a specific spatial configuration. Therefore, formation control is often a prerequisite for other applications in motion planning.

Connectivity maintenance refers to providing constraints which guarantee a connected communication graph among the agents at all times.

More precisely, connectivity maintenance for a formation of agents refers, generally, to the ensemble of feedback laws, path/trajectory planning, collision avoidance and task assignment which guarantee that agents can communicate between themselves (to ensure data flows for, e.g., control decisions, data gathering, distributed estimation, etc.).

Remark 2. [START_REF] Relay | [END_REF]. Even though connectivity maintenance and formation control are not different labels for the same phenomenon, the underlying idea is similar, the difference being the global objective. A key difference between the two notions is the relevance of relative distance between agents: for connectivity maintenance the unobstructed communication is the key aspect, while the inter-distance or relative velocity are relevant only in their effect on the "line of sight" between agents (e.g., distance less than a communication range threshold). ♦

Formation control aims to provide constraints guaranteeing desired relative distances and velocities among a group of agents.

There are two main approaches in the literature which tackle the formation control problems [START_REF] Qu | Cooperative Control of Dynamical Systems: Applications to Autonomous Vehicles[END_REF]]: the leader-follower design and leaderless approach. The former consists in designating one of the agents as a leader moving in a certain way with the remaining agents tracking the leader in order to maintain formation around it. The later involves the coordination of the agents through a global consensus for attaining the global objective.

Choosing between these two approaches depends on the particularities of the team/group of agents, i.e., their capability of communication and sensing. Moreover, these aspects affect further the control architecture, as will be detailed in Section B.

Corner cutting avoidance conditions

The corner cutting avoidance problem is an important but often overlooked part in motion planning. Obstacle and collision avoidance constraints are usually imposed at the sampling time without regards to the intra-sample behavior of the agent(s). Hence, it is possible for an agent to "cut the corner" of an obstacle while apparently respecting the constraints.

Corner cutting aims to provide constraints guaranteeing intra-sample collision avoidance.

The idea is conceptually simple: an agent avoids corner cutting iff its next position lies in its visibility region (the union of all rays spanned from the current position not intersecting any of the obstacles). Note that the dual notion (the under-shadow region -the region hidden from the viewpoint of an agent) may be used to enforce full coverage of the feasible space, e.g., the gallery problem.

Since describing exactly such regions is impractical, mixed-integer formulations (in either exact or over-approximated form) are used. Within the hyperplane arrangement framework presented earlier, these constraints employ (either implicitly or explicitly) three sign tuples: the current and future positions of the agent, σ, σ + ∈ Σ X\P , and the coordinates of the obstacle(s), σ • ∈ Σ P .

We are aware of results from [Maia and R. K. Specifically, in [Maia and R. K. Galvao 2009], to avoid cutting a single obstacle, there has to exist at least one half-space containing the obstacle which does not contain the current and successor positions of the agent: .33) Since such conditions often appear in MPC problems where both current and successor sign tuples are decision variables, (2.33) is rewritten such as to avoid nonlinearities:

∃i s.t. σ(i) = σ + (i) = 0. ( 2
i [1 -σ • (i)] σ + (i) < N + i [1 -2σ • (i)] σ(i), (2.34)
for all possible values σ • ∈ {0, 1} N . The idea is that from all constraints (2.34) at least one reduces to 17 (2.33), i.e., the one for which σ • = σ + .

[A. [START_REF] Richards | Inter-sample avoidance in trajectory optimizers using mixed-integer linear programming[END_REF] improves on [START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF]] by reducing the number of constraints (2.34) from 2 N to a more manageable N . This is done by forcing two consecutive positions x, x + to respect the same constraint, with our notation:

-h ⊤ i x ≤ -k i + M σ + (i) and -h ⊤ i x + ≤ -k i + M σ + (i) for all i = 1 . . . N
. This implies that there exists at least one index i s.t. both x and x + lie on the same side of the hyperplane (and thus on the opposite side from the obstacle), similar with (2.33). [R. [START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF] proposes a logarithmic scheme to reduce the number of binary variables involved in the selection of the active hyperplanes.

The common shortcoming of these approaches is that they do not easily handle the multiagent multi-obstacle case. [START_REF] Stoican | Exact and overapproximated guarantees for corner cutting avoidance in a multiobstacle environment[END_REF] generalizes (2.33) with: .35) As stated earlier, when both σ, σ + are variables, (2.35) becomes nonlinear due the bi-linear term appearing in it. Hence, a relaxation such as in (2.34) may be used: .36) for all σ • ∈ Σ X\P , i.e., (2.36) reduces to (2.35) when σ • = σ + .

i |σ •,j (i) -σ(i)| • |σ •,j (i) -σ + (i)| > 0, ∀σ •,j ∈ Σ P . ( 2
i |σ •,j (i) -σ • (i)| • |σ •,j (i) -σ + (i)| > -|σ • -σ|, ∀σ •,j ∈ Σ P , ( 2
17 Both (2.33) and ( 2.34) assume a single obstacle (Σ P = {σ •,1 }) where, by convention, σ

•,1 (i) = 1, ∀i. This implies that i σ(i) ≤ N -1, i σ + (i) ≤ N -1.

Connectivity maintenance conditions

The optimal dynamic formation control is a popular topic in the multi-agent systems field, and MIP-based methods provide a suitable tool for the modeling of the corresponding constraints. In what follows, we briefly present the topics where the connectivity is induced by the internal factors (as, e.g., communication).

The mixed-integer constraints enforce that each agent is within the communication range of at least one other agent. Moreover, pairs isolation need to be avoided, i.e., if an agent is enforced to stay in a communication range of another, the latter cannot be constrained to stay in the range of the former, but in the range of a third one. In order to formulate the corresponding mixed-integer constraints we proceed as in [R. J. Afonso, Maximo, and Galvão 2020]. Hence, we consider a binary variable ξ i,j,k ∈ {0, 1} which is 1, if the i-th agent is inside of the communication range of j-th agent at the timestep k, and 0 otherwise. Next, we can write the following constraint, ensuring the communication connectivity among the N agents:

N -1 i=1 N j=1 ξ i,j,k = N -1, 0 ≤ k ≤ T, (2.37a
)

v i ,v j ξ i,j,k ≤ #S -1, 0 ≤ k ≤ T, ∀S ⊂ V, (2.37b
)

ξ i,j,k = 0, 1 ≤ i ≤ N, 1 ≤ j ≤ i, 0 ≤ k ≤ T. (2.37c)
The constraint (2.37b) prevents the possible pairs isolation, imposing for each sub-graph S of the original communication graph V that at least one node has to be able to communicate with a node outside the graph. Moreover, alongside with (2.37a) it guarantees that only connected communication graph can be obtained (for a detailed discussion, see Theorem 2 from [R. J. Afonso, Maximo, and Galvão 2020]), i.e. the connectivity is maintained.

Furthermore, [X. [START_REF] Sun | Optimal dynamic formation control of multi-agent systems in environments with obstacles[END_REF] is focusing on the coverage control problem.

The agent team works in a 2D mission space, maximizing coverage in this space. The formation control is approached from a leader-follower setting perspective, the leader moves on a given trajectory while the remaining agents must maintain the formation. Thus, the formation becomes dynamic as soon as the leader starts moving along the given trajectory, and it must adapt to any environmental change (a new mission, a new composition of the team or some obstacles detected). Firstly, it is presented a brief formulation for a general optimal formation problem as an optimization problem, whose objective function depends on the agent spatial positions. The constraints of the optimization problem are given by a feasible space for positions (included in the 2D mission space) and by the condition that an undirected graph must be connected. This graph is modeling the desirable links among agents. The feasible space can be either convex (without obstacles) or non-convex (with obstacles). This aspect complicates the way in which the graph constraint should be treated. Hence, these two cases must be tackled separately. For a convex feasible space, the optimal dynamic formation problem is rewritten as a MINLP, by introducing a set of flow variables over the undirected graph (associating to every link between two agents an integer flow amount). The problematic constraint is transformed into a set of mixed-integer inequalities. Also, the computationally demanding re-solving of the MILP for all time instants is avoided, providing a sufficient condition for maintaining optimal formation in a certain time interval. For a non-convex feasible space, the problem is solved in a different approach. At each moment, when the connectivity is lost, a new undirected connected graph is constructed such that the effort to maintain the initial formation is minimal. An algorithm is developed to construct the new graph, and this is used as an input for CPA (Connectivity Preservation Algorithm). An illustrative example of this second situation is considered therein, where the initial formation is computed using the MILP from the convex feasible space case, obstacles being considered afterwards.

Connectivity induced by external constraints

This subsection introduces the types of objectives/ goals encountered in the formation control, goals enforced by the external18 factors.

One of the classical examples where connectivity maintenance is externally induced is the problem of efficient coverage of a specific area by mobile sensors. For instance, [START_REF] Ritter | Adaptive Observation Strategy for Dispersion Process Estimation Using Cooperating Mobile Sensors[END_REF]] tackled this problem for a mobile sensors platform composed of a number of sensor-equipped autonomous vehicles. An adaptive observation strategy for on-line estimation of the state of a dispersion process is proposed, based on the model of the process and the multi-agent collaborative systems. The multi-vehicle cooperating system is modelled using MLD systems [Alberto Bemporad and Morari 1999]; and controlled such that the target location is reached in a reasonable time while simultaneously avoiding any possible collision among the agents. The target positions of the agents are received from the PDE19 -based state estimation, more precisely the optimal measurement locations. Thus, the MLD-systems represent a good alternative for modeling multi-vehicle systems, because all existing constraints (dynamics, collision avoidance, measurement etc.) are embedded in a complete configuration, which allows to use standard methods to solve the control problem.

As well, a standard problem is given by the rendezvous maneuver, a group of agents (usually, a fleet of spacecraft) has to achieve a specific formation in a minimum time. [Papen et al. 2017] treat this problem for a fixed-wing-UAV fleet in an environment containing both static and dynamic obstacles. The control strategy is implemented through a distributed architecture (specifically, MPC applied for each agent), an agent knowing only its own trajectory (defined by its current position and velocity). By introducing binary variables to model the obstacles and to eliminate some nonlinear constraints (like velocity and acceleration constraints), the MPC optimization problem is transformed into a MILP. Noteworthy, the wake vortex of turbulent air formed behind every UAV is modeled as a dynamic obstacle because within these regions it is arduous to maintain control over the UAVs. To solve the resulting MILP is necessary a considerable time, but with some relaxation of the constraints and fine-tuning the complexity is limited.

Conclusions and open problems

In the preceding material we provided our evaluation on the state-of-the-art for MIP-based motion planning and also we aim to identify active topics and open problems in this field.

It is important to mention that the valuable insights in the description of a non-convex feasible region represent an useful modeling tool not only for the motion planning but also for broader control fields.

As mentioned before, although the history of MIP starts almost 60 years ago, the interest of the control and robotic community on this topic is relatively recent and the research in this area is quite active. It is obvious that there exists a substantial progress in all of the topics mentioned in Section 2. 3. However, there are many points that can be further enhanced since the developing of new and performing methods to provide exact solution for MIP is exponentially growing. For the interested reader, Section C provides an overview of software tools that have been recently employed in resolution of MIP.

As open and active problems regarding mainly obstacle and collision avoidance topic we can mention the trade off between conservativeness and complexity in non-convex region modeling and representation. Even though there exist several valuable improvements on classical MIP formulations, the complexity remains an arduous issue imposing restrictions: only a small size problem can be solved in a real-time manner. The improvements may be accomplished by exploiting the underlying combinatorial structure of the MIP formulation. Nevertheless, whenever the problems are inherently non-convex and/or they involve alternative choices, then mixed-integer representations provide an useful and powerful tool, but we have to proceed with caution evaluating the structural properties that may lead to compact formulations. As well, important progress may be expected in the direction of the MIP resolution algorithms, which can be adapted and improved in order to attenuate the inherent complexity of the mixed-integer approaches by capitalizing on the different particularities of the problems to be solved. Additionally, another aspect that can be viewed as an open issue and, in some sense, is generated by the previous ones is that there are not so many works in the field of MIP-based motion planning which validate the theoretical results over, at least, experimental platforms, even though there are many applications able to benefit from their use, e.g., search and rescue, environmental measurements tasks, area coverage and the like.

This chapter presented a review of recent research and developments in MIP-based motion planning. In addition to the results reviewed above, there are many other publications treating MIP formulations in feedback decision making that could not find a direct link with the scope of the motion planning in spite of our best effort. Covering all the variety of application and often inhomogenous problem formulations involving MIP is a task that can be challenging even for a review.

Mixed-integer description of non-convex feasible regions

. . 107

The present chapter addresses the problem of collision avoidance in a multi-obstacle environment and focuses on its representation in optimization-based control problems. As mentioned before (Section 2.2), the design problem is commonly stated in the literature in terms of a constrained optimization problem over a non-convex domain which can be efficiently characterized using hyperplane arrangements (Section 2.1). In what follows, we consider additional structural assumptions by the use of zonotopic over-approximation and highlight their benefits when introduced in the obstacle avoidance problem.

As stated in Section 2.1, the obstacles are often modeled as convex polytopic regions in which case the resulting non-convex and non-connected feasible domain can be efficiently described through its associated hyperplane arrangement [Ionela [START_REF] Prodan | Mixed-integer representations in control design: Mathematical foundations and applications[END_REF]] which allows to characterize the domain as a union of convex cells. Next, a mixedinteger formulation is used to characterize the union of cells and is introduced in a standard constrained optimization problem. Techniques like cell-merging [I. Prodan et al. 2012] or logarithmic formulations [J. P. Vielma and Nemhauser 2011] may be employed to reduce the complexity of the formulation but nonetheless their numerical complex observe a combinatorial increase with the number of obstacles and/or agents.

Therefore, the following question appears:

"Can the number of cells for the feasible region be decreased by over-approximating the obstacles safeguarding the path-planning features of the initial problem? "

In the first part of the present chapter, we propose to address this question by considering zonotopes and their properties for such obstacle over-approximation and cell-counting while providing an explicit measure of the problem complexity in terms of the total number of cells.

Due to their intrinsic properties, zonotopes are widely used in control, e.g, for collision detection [Guibas, A. Nguyen, and L. Zhang 2003], reachability analysis [Althoff, Stursberg, and Buss 2010], fault diagnosis [START_REF] Puig | Fault diagnosis and fault tolerant control using setmembership approaches: Application to real case studies[END_REF]] or guaranteed state estimation [Alamo, Bravo, and Camacho n.d.]. Moreover, there are toolboxes like CORA [START_REF] Althoff | CORA 2015 Manual[END_REF]] able to manage zonotopic sets representation.

We start from exact formulations (analytic formulation of the zonotope's volume, exhaustive separation conditions, etc) and relax them to a linear constrained form (which depends on a collection of centers and scaling factors of the zonotopes representation). Particular attention is given to the control-related issues. Foremost, we handle the approximation such that separation between any two obstacles is guaranteed (thus avoiding changes in the domain topology and subsequently, infeasibility in the motion planning procedure). Lastly, the collision avoidance problem is formulated as a predictive control optimization problem whose constraints are activated only when the obstacles are contained in the reachable set of the agent.

Specifically, some of the noteworthy aspects of the approach can be resumed as follows: i) we provide tight zonotopic approximations of the obstacles and various measures for these over-approximations;

ii) we provide an exact bound for the number of generators such that the complexity is reduced;

iii) we ensure separation among the zonotopic approximations, i.e., preserving the original space topology; iv) we formulate an optimization problem correlated with reachability properties for constraint activation.

Zonotopic framework

In order to efficiently describe the non-convex feasible region for a dynamical (mobile) agent, we recall the notions and techniques in Section 2.1 and the references therein, which involve polyhedral sets and hyperplane arrangements. .1) and it has some interesting properties. Let us consider two sets X and Y , then:

h Q (η) = sup q∈Q η T q. ( 3
i) h αX (η) = αh X (η), ∀α ≥ 0, ii) h X⊕Y (η) = h X (η) + h Y (η), iii) h AX (η) = h X (A ⊤ η), where A is a linear transformation, iv) the inclusion X ⊆ Y holds iff h X (η) ≤ h Y (η), ∀η ∈ R d . ♦
Furthermore, if Y is a polytope in the half-space form:

Y = {x ∈ R n : s ⊤ i x ≤ r i , i = 1 . . . n h }, (3.2)
then condition iv) has to be checked only for a finite number of vertices, i.e. the condition becomes:

h X (s i ) ≤ r i , ∀i. (3.3)
Additionally, if X is also a polytope in the extreme vertices form:

X = {x ∈ R n : x = α j v j , α j = 1, α j ≥ 0, j = 1 . . . n v } (3.4)
the inclusion condition is rewritten as:

s ⊤ i v j ≤ r i , ∀i, j. (3.5)
In what follows, we will use extensively zonotopes as a particular class of polytopes, endowed with a third representation due to their symmetry. Owing to all their properties, zonotopes provide a good balance between numerical complexity and fidelity of representation.

Definition 3.2 (Zonotopes-[Kuhn 1998]).

A zonotope is a centrally symmetric polytope, which can be described as a Minkowski sum of line segments. In its generator representation a zonotope Z(G, c) is described by center c ∈ R d and generator matrix G = g 1 . . . g m ∈ R d×m :

Z(G, c) = {c + m k=1 ξ k g k : ξ ∞ ≤ 1}. (3.6)
Zonotopes own several properties of practical interest [START_REF] Fukuda | From the zonotope construction to the Minkowski addition of convex polytopes[END_REF]]:

i) are closed under linear transformation:

RZ(G 1 , c 1 ) = Z(RG 1 , Rc 1 ); (3.7)
ii) are closed under Minkowski sum:

Z(G 1 , c 1 ) ⊕ Z(G 2 , c 2 ) = Z( G 1 G 2 , c 1 + c 2 ); (3.8)
iii) are symmetric, up to their center:

-Z(G, c) = Z(G, -c); (3.9)
iv) their volume has an explicit formulation [START_REF] Gover | Determinants and the volumes of parallelotopes and zonotopes[END_REF]:

Vol(Z(G, c)) = 1≤j 1 <j 2 ...j d ≤m det(G j 1 ...j d ) , (3.10)
where G j 1 ...j d denotes the matrix composed from columns of indices j 1 . . . j d from G.

v) their corresponding half-space representation [Althoff, Stursberg, and Buss 2010]: .11) .12) where a pair

Z(G, c) = 1≤j 1 <...j d-1 ≤m {x ∈ R d : h ⊤ i (x -c) ≤ k i }. ( 3
h i ⊥g j l , ∀j l ∈ {j 1 . . . j d-1 }, k i = j l / ∈{j 1 ...j d-1 } h ⊤ i g j l . ( 3
(h i , k i ) ∈ R d × R corresponds to each sequence of d -1 generators 1 ≤ j 1 < j 2 . . . j d-1 ≤ m.
Adding the center from generator representation, a zonotope is formulated in terms of its halfspaces as in (3.12).

Remark 3.1. For general convex sets (as for the polytopes) the operations i) and ii) may demand a considerable computational effort and they are numerical unstable for greater dimension. Nevertheless, a minor drawback for (3. 7) and (3.8) is that their result can be higher order than the initial zonotopes. This inconvenient can be avoided using methods for enclosing a given zonotope within one of a lower order (at the price of overapproximations), see, e.g., [START_REF] Althoff | Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization[END_REF]. ♦

Using the support functions properties (see Definition 3.1) in combination with the definition of a zonotope means that: a) the inclusion Z(G, c) ⊆ Y , with Y defined as in (3.2), is valid iff: .13) b) the inclusion of a polytopic set X, defined as in (3.4), into a zonotope X ⊆ Z(G, c) holds iff:

s ⊤ i c + j s ⊤ i g j ≤ r i ∀i, j. ( 3
h ⊤ i (v j -c) ≤ k i ∀i, j. (3.14)

Hyperplane Arrangements

Consider a finite collection of hyperplanes from R d :

H = {H i } i∈I l (3.15) with H i = {x ∈ R d : s ⊤ i x = r i }.
Each of these hyperplanes divides the space in two disjoint regions:

R + i = x ∈ R d : s ⊤ i x ≤ r i , (3.16a)
R - i = x ∈ R d : -s ⊤ i x ≤ -r i . (3.16b)
Based on the Definition 2.1 of hyperplane arrangements, as stated in Section 2.1, several aspects regarding complexity features can be discussed.

Definition 3.3. An arrangement 1 A is considered to be in general position iff any pair of hyperplanes is intersecting, i.e., they are not parallel an do not share the same normal.

♦

As a side remark, if a hyperplane arrangement is in general position, then the number of cells remain constant whenever there is a perturbation in the position of the hyperplanes.

Remark 3.2. For a given hyperplane arrangement A(H), the number of feasible cells is bounded by the Buck's formula [START_REF] Buck | Partition of space[END_REF]], w.r.t. the space dimension d and the number of hyperplanes N = #I, as follows:

γ(N ) ≤ d k=0 N k (3.17)
with equality satisfied if the arrangement is in general position. ♦ Example 3.1. Consider the hyperplane arrangement depicted in Figure 3.1a, described by N=3 hyperplanes H i , i = 1, 2, 3. These hyperplanes divide the entire space (R 2 , d = 2) in 7 feasible cells, each characterized by a unique sign tuple σ. Hence, the feasible sign tuples are: (+, +, +), (-, +, +), (-, +, -), (-, -, -), (+, -, -), (+, -, +) and (-, -, +).

The last tuple characterizes the shadowed region

A = 3 i=1 R σ(i) i = R - 1 R - 2 R + 3 .
We note that the remaining sign tuple (+, +, -) is unfeasible and describes

A = 3 i=1 R σ(i) i = R + 1 R + 2 R - 3 = ∅.
In line with Definition 3.3 and since the number of cells checks the equality from (3.17), the arrangement in Figure 3.1a is in general position.

Each of the regions of the hyperplane arrangement, as in (2.11), can be either bounded or unbounded. Therefore, we denote: the total number of regions as r(A) and the number of bounded regions as b(A). For instance, for Example 3.1 we have r(A) = 7, and b(A) = 1, as it can be seen in Figure 3.1a.

H 3 H 1 H 2 A(-, -, +) (a) Hyperplane arrangement • ∅ • H 1 • H 1 ∩ H 2 • H 2 • H 2 ∩ H 3 • H 3 • H 3 ∩ H 1 rank = 0 rank = 1 rank = 2
(b) The central (sub-) arrangements Before proceeding further, some auxiliary notions are required. Then the characteristic polynomial of A is defined as

Definition 3.4. A (sub-)arrangement B ⊆ H is called central if

H i ∈B H i = ∅.
χ A (t) = B⊆A,B central (-1) #B t d-rank(B) (3.18)
Furthermore, the total number and the number of bounded regions characterizing the arrangement are: .19b) Remark 3.5. The equations (3.19) are a particularization of the Zaslavsky's Theorem [START_REF] Zaslavsky | Facing up to arrangements: Face-count formulas for partitions of space by hyperplanes: Face-count formulas for partitions of space by hyperplanes[END_REF]]. ♦ Remark 3.6. For a hyperplane arrangement with N hyperplanes in general position, the bounds are .20) Let us revisit Example 3.1 and note that the central sub-arrangements are:

r(A) = (-1) d χ A (-1), (3.19a) b(A) = (-1) rank(A) χ A (1). ( 3
r(A) = d k=0 N k , b(A) = N d -1 . ( 3
{R 2 , H 1 , H 2 , H 3 , H 1 ∩ H 2 , H 2 ∩ H 3 , H 1 ∩ H 3 }.
Each of these elements leads to a term in (3.18). In Figure 3.1b we depict the corresponding rank of each element. For instance, if B = H 1 ∩ H 3 then #B = 2 and rank(B) = 2 (i.e., the two hyperplane span R 2 ). That leads to the following term in Whitney's formula:

(-1) 2 t d-2 = 1.
Similarly, we obtain the other components of the formula (3.18) as follows:

χ A (t) = (-1) 0 t 2-0 + (-1) 1 t 2-1 + (-1) 1 t 2-1 + (-1) 1 t 2-1 + (-1) 2 t 2-2 + (-1) 2 t 2-2 + (-1) 2 t 2-2 = t 2 -3t + 3 (3.21)
Further we use (3.19) and obtain the total number of regions r(A) = (-1) 2 χ A (-1) = 7 and the number of bounded regions b(A) = (-1) 2 χ A (1) = 1.

Example 3.2. Consider the hyperplane arrangement in Figure 3.2a, which builds up the arrangement in Figure 3.1a by adding a fourth hyperplane H 4 . Hence, the number of feasible sign tuples doubles (containing elements as:(+, +, +, -), (+, +, +, -)). The highlighted region is described by the tuple (-, -, +, -), i.e.,

A = 3 i=1 R σ(i) i = R - 1 R - 2 R + 3 R - 4 .
We note that the tuple (+, +, -, +) The notions presented above have an instrumental role in counting the regions without explicitly computing them. Hence, by employing combinatorial notions and avoiding a burdensome geometrical analysis, one is able to provide a measure of the complexity of the representation. As seen in Remark 3.6, for arrangements in general position the computational effort is negligible, having explicit formulas. However, for more complex cases, counting regions need a different, more practical, approach. For example, efficient algorithms are proposed in [START_REF] Avis | Reverse search for enumeration[END_REF]] and [START_REF] Geyer | Efficient mode enumeration of compositional hybrid systems[END_REF].

describes A = 3 i=1 R σ(i) i = R + 1 R + 2 R - 3 = ∅, and
H 3 H 1 H 2 H 4 H 4 A(-, -, +, -) (a) Hyperplane arrangement H 3 H 1 H 2 H 4 H 4 (b) Perturbed hyperplane arrangement

Illustrative examples

In the following subsection, we consider two more complex multi-obstacle environments. The methods in YALMIP [START_REF] Lofberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF]] have been used in order to count the feasible cells.

Consider the collection of obstacles (blue regions in Figure 3. Before the main developments we caution that a zonotope is centrally symmetric. Hence, not any convex shape can be rigorously approximated by an enclosing zonotope and this will be the main source of conservatism hereafter (and a potential direction of improvement [J. [START_REF] Scott | Constrained zonotopes: A new tool for set-based estimation and fault detection[END_REF]). In particular, notice that a symmetric approximation (zonotope) with a pre-defined topology (pre-imposed fixed direction of the generators) can be arbitrarily conservative. Nonetheless, in the literature a variety of algorithms provides over-approximations [Guibas, A. Nguyen, and L. Zhang 2003], [Althoff, Stursberg, and Buss 2010] and many of them consider pre-defined families of zonotopes (e.g., by fixing the directions of the generators). In corner cases, the approximation may be arbitrarily bad, i.e., the ratio between approximated and actual volume tends towards an large value. For example in Figure 3.4, this ratio tends to infinity for small values of the parameter ϵ, more precisely we have:

lim ϵ→0 V □ V △ = ∞,
where V □ stands for the volume of the enclosing zonotopic over-approximation, and V △ is the actual volume of the convex shape.

Zonotopic framework

Considering Definition 3.2 we refer to a family of zonotopes parametrized after their centers c ℓ ∈ R d and scaling factors ∆ ℓ ∈ R m×m applied to a common generator "seed" (an a priori given matrix G ∈ R d×m ):

Z(G∆ j , c j ), j = 1 . . . N o . (3.23)
∆ ℓ is a diagonal matrix whose diagonal elements may be either equal or distinct 4 . The k-th diagonal element is noted as δ j k .

Having a common generator seed:

i) provides linear inclusion constraints in term of the scaling and center parameters, ii) allows measures for the obstacle over-approximation in terms of ℓ 1 , ℓ ∞ norms and the volume, iii) expresses explicitly the complexity of the representation (i.e., the number of cells of the associated arrangement).

Set inclusions for parametrized zonotopes

Replacing g k with g k • δ j k in (3.12), the half-space representation of the j-th zonotope from (3.23) is given by:

h i s.t. h i ⊥ g k , ∀k ∈ {k 1 . . . k d-1 }, (3.24a
)

k i (∆ j ) = k / ∈{k 1 ...k d-1 } h ⊤ i g k δ j k , (3.24b)
where i enumerates the

p(d, m) = m d -1 (3.25)
combinations of d -1 distinct generators selected from the list of m available ones (i.e.,

1 ≤ k 1 < • • • < k d-1 ≤ m).
Remark 3.7. Note that h i remains unchanged with respect to (3.12) as long as the normal on {g k } k∈{k 1 ,.

Z(G∆ j , c j ) ⊆ Y : s ⊤ i c j + m k=1 s ⊤ i g k • δ j k ≤ r i , ∀i, (3.27a
)

X ⊆ Z(G∆ j , c j ) : h ⊤ i (v k -c j ) ≤ k i (∆ j ), ∀k. (3.27b) Remark 3.8.
Using X in its vertex representation (3.4) means that (3.27b) is translated into finding the tight enclosure of a set of points in R d by a zonotope. [Guibas, A. Nguyen, and L. Zhang 2003] enforces a priori given points (v k in (3.27b)) to stay on the zonotope boundary. Hence, the "inclusion" condition is written as:

v k = c j + m i=1 ξ ijk g i , -δ j i ≤ ξ ijk ≤ δ j i , ∀i, k. (3.28)
The advantage of condition (3.28) is that it uses directly the generator representation but requires a large number of inequalities. ♦

Measures for zonotope approximations

The overall goal is to provide adequate over-approximations (3.23) for the collection (3.22).

That is, seek a (intrinsically symmetric) zonotope Z(G∆ j , c j ) enclosing the (usually nonsymmetric) polytope P j such that a specific measure parametrized after c j , ∆ j is minimized:

(∆ j , c j ) * =arg min ∆ j ,c j C(∆ j , c j ) (3.29a) s.t. P j ⊆ Z(G∆ j , c j ). (3.29b)
In the following proposition we consider several measures.

Proposition 3.1. For the cost defined as in (3.29), the following measures are available: i) zonotope volume Vol(Z(G∆ j , c j )): .32)

C(∆ j , c j ) = 1≤k 1 <•••<k d ≤m   det(G k 1 ...k d ) • k∈{k 1 ,...,k d } δ j k   ; (3.30) ii) generator sum m k=1 g k δ j k : C(∆ j , c j ) = ||∆ j || 1 = m k=1 δ j k ; (3.31) iii) largest generator max k=1...m g k δ j k : C(∆ j , c j ) = ||∆ j || ∞ = max k=1...m δ j k . ( 3
Proof. i) Replacing g k with g k • δ j k in (3.10) leads to Vol(Z(G∆ j , c j )) = 1≤k 1 <•••<k d ≤m det(G k 1 ...k d • ∆ k 1 ...k d j ) .
Noting that det(M • N ) = det(M ) • det(N ) and that ∆ j is a diagonal matrix with positive elements gives (3.30).

ii) Assuming without loss of generality that ||g k || 1 = 1, the cost reduces immediately to (3.31).

iii) Similar reasoning with case ii), assuming ||g k || ∞ = 1.

Remark 3. 9. Volume (3.30) is a sum of polynomial terms δ j k , thus, non-linear. Imposing equality among the scaling factors (δ j 1 = . . . = δ jm = δj ), leads to a simplified volume formulation: Vol(Z(G∆ j , c j )) =

1≤k 1 <•••<k d ≤m det(G k 1 ...k d ) • δd j , (3.33)
which can be used instead of (3.30)

Complexity bounds within a zonotopic framework

In this section, we present a detailed analysis of the problem complexity in terms of the total number of cells. Specifically, having a collection as in (3.26), which is generated by the zonotopic approximations (3.23) and induces the hyperplane arrangement A( H), we aim to:

i) provide a cell count with explicit dependence on m, the number of generators and N o , the number of obstacles;

ii) compare it against the original arrangement A(H).

Hyperplane arrangement induced by the zonotopes

As shown in (3.26), the 2 • p(d, m) • N o hyperplanes are partitioned in 2N o families: in each family there is a hyperplane with normal h i but with a different offset ±k i (∆ j ). This particularity allows the following proposition.

Proposition 3.2. A( H) has the characteristic polynomial:

χ(A) = d k=0 (-1) k t d-k • p(d, m) k • (2N o ) k , (3.34)
Proof. The left-side of (3.38) comes from (3.35) and the right side from (3.20) with n → n * o .

While the previous results hold for R d , it is worthwhile to particularize them for d ∈ {2, 3}.

Corollary 3.2. For the case d = 2, we have that:

1. the total and the bounded number of cells for A( H) is: ii) Introducing p(2, m) in (3.38) reduces it to a second-order inequality:

r(A) = 1 + m • 2N o + m(m -1) • 2N 2 o , (3.39a) b(A) = 1 -m • 2N o + m(m -1) • 2N 2 o ; (3.39b)

the largest m s.t. A( H) has fewer cells than A(H) is:

m = (N o -1) + (N o -1) 2 + n * o (n * o + 1) 2N o . ( 3
m 2 • 2N 2 o -m • 2N o (N o -1) - n * o (n * o + 1) 2 ≤ 0. (3.41)
The inequality (3.41) has two real solutions (one negative and one positive) which means that the largest m verifying the inequality is the closest (at its left) from the positive solution, hence, (3.40).

Corollary 3.3. For the case d = 3, we have that:

i) the total and the bounded number of cells for A( H) is:

r(A) = 1 + η • 2N o + η(η -1) • 2N 2 o + 1 6 η(η -1)(η -2)8N 3 o , b(A) = 1 -η • 2N o + η(η -1) • 2N 2 o - η(η -1)(η -2) 6 8N 3 o ;
where η = 1 2 m(m -1); ii) the largest m s. t .43) where η * represents the real solution of a solvable third order equation (3.45).

. A( H) has fewer cells than A(H) is :

m = 1 + √ 1 + 8η * 2 . ( 3
Proof. For d = 3, the number of support hyperplanes for a zonotope becomes p(3, m) = m ii) Introducing p(3, m) in (3.38) reduces it to a third-order inequality in η as it follows:

2(8N 3 0 -12N 2 0 + 12N 0 )η -3(8N 3 0 -4N 2 0 )η 2 + 8N 3 0 η 3 ≤ n * o + 1 6 n * o (n * o -1)(n * o + 1). (3.44)
Noting λ for the right term of (3.44) we rewrite it as: .45) Therefore, we can compute the discriminant: .46) We note that ∆ is negative for the so-defined positive values of N 0 and λ. This means that the third order equation has one real solution and two complex ones. The computing of them involves calculating:

8N 3 0 η 3 -12N 2 0 (2N 0 -1)η 2 + 8N 0 (2N 2 0 -3N 0 + 3) -λ ≤ 0. ( 3
∆ = N 6 0 [-359424 -34560λ -1728λ 2 + (829440 + 82944λ)N 0 -(1022976 + 27648λ)N 2 0 + 663552N 3 0 -221184N 4 0 + 16384N 6 0 ]. ( 3
∆ 0 = N 4 0 (-432 + 192N 2 0 ), (3.47a
)

∆ 1 = N 6 0 (-17280 -1728λ + 41472N 0 -13824N 2 0 ), (3.47b
) .47c) Thus, the real solution is given by:

C = 3 ∆ 1 ± 1728N 2 0 ∆ 2 . ( 3
η * = 1 24N 3 0 12N 2 0 (2N 0 -1) + C + ∆ 0 C . (3.48)
Introducing (3.48) in the η definition, we have a second-order equation in m with two real solutions (one negative and one positive). Thus, the largest m verifying the inequality is the closest (at its left) from the positive solution, (3.43) .

Illustrative examples

Let us revisit the previous examples and proceed to find the zonotopic over-approximations of the obstacles in Figure 3. .50) Maintaining the structure of Table 3.1, we delineate in Table 3.2 some parameters of interest. Thus, we indicate the total computing time corresponding to each considered method (t sol ), the relative modification of the number of cells ( ∆γ(N ) γ(N ) ) and a couple volume specifications: the volume of the over-approximation (V ) and the relative error with respect to the volume of the polytopic obstacles ( ∆V V ).

G 1,3 ∈ 1 0 1 0 1 1 , -1 0 -1 0 1 1 , (3.49) G 4,5 ∈      1 0 0 0 1 0 0 0 1   ,   1 0 0 1 0 1 0 1 0 0 1 1      . ( 3
Chapter 3. Mixed-integer description of non-convex feasible regions Measure Constraints G by it. For instance, in [S. [START_REF] Olaru | Positive invariant sets for fault tolerant multisensor control schemes[END_REF]] in order to find the equation describing this hyperplane, the definition of the hyperplane is rewritten as:

t sol #H γ * (N ) ∆γ(N ) γ(N ) (%) t γ * (N ) #Σ P V ∆V V (%) d = 2 δ 1 ( 3
h(x) = γ ⊤ 1 x = 0, (3.52)
where γ ⊤ = γ 0 γ 1 . . . γ d . Next, the value of γ is obtained by solving the following optimization problem:

min γ 2 (3.53a) s.t. 1 v ⊤ k γ ≥ 1, ∀v k ∈ V(P i ), (3.53b
)

1 ν ⊤ ℓ γ ≤ 1, ∀ν ℓ ∈ V(P j ); (3.53c)
where P i and P j are the polytopes to be separated. ♦ Remark 3. [START_REF] Relay | [END_REF]. Alternatively to Remark 3.12, a separating hyperplane can be obtained as follows. We compute an obstacle-based space partitioning, {X i } i=1:No , as in (4.20), H ij is the common facet between any two neighboring cells. This space-partitioning method is exhaustively presented in Chapter 4.

♦

The first approach in overcoming the potential overlapping of the zonotopic over-approximation consists in considering an additional constraint in (3.29). More precisely, this constraint enforces that the approximation to stay above or below a separation hyperplane. The decision between above or below is made based on the position of the center of the polytope relative to this given separation hyperplane. In other words, we impose:

Z(G j ∆ j , c j ) ⊂ R + ij or Z(G j ∆ j , c j ) ⊂ R - ij , ∀i. (3.54)
The reformulation of (3.54) as an LP is straightforward: all vertices of the zonotopic approximation have to stay on the same side of the separation hyperplane as its center. Thus, the optimization problem (3.29) takes the form (3.55) with (3.55c) formulated recalling the inclusion condition of a parametrized zonotope into a polytope (3.27).

min Proof: The linearity of separation constraints (3.55c) directly leads to the convexity of the feasible space regarding zonotopic over-approximation.

∆ j ,c j C(∆ j , c j ) (3.55a) s.t. P j ⊆ Z(G∆ j , c j ) (3.55b) Z(G j ∆ j , c j ) ⊂ R ± ij , ∀i. (3.
Remark 3. 14. Certainly,(3.55c) enforces the corridors but it represents an additional structural constraint and impacts the feasiblity of the zonotopic overapproximations. An example in this sense is depicted in the Figure 3.9, where we take a triangle in R 2 and proceed to over-approximate with a parametrized zonotope with G given by the Euclidean unit vectors. The problem (3.55) is infeasible. ♦

Despite the limitation highlighted in Remark 3.14, there exists a strong constructive result for the existence of a zonotopic over-approximation which satisfies the corridor constraints. The idea is to generate a face of the zonotopic over-approximation h ⊤ 1i x ≤ w 0i with i = 1, . . . , d+1 be the "d-1"-dimensional facet of this simplex. Let P 2 , ..., P d+2 be zonotopic sets, each one having within the half-space representation one and only one of the constraints h ⊤ 1i x ≤ w 0i +ϵ. Obviously for ϵ → 0, P 2 , . . . , P d+1 remain zonotopic while the feasible set allowing a zonotopic approximation of P 1 reduces to P 1 itself. But P 1 is not a zonotopic set and thus no over-approximation exists 5 . Proof: We can construct the matrix G such that the resulting over-approximation to have the facets parallel with the "d + 1" separating hyperplanes. The obtained zonotopic over-approximation can be adjusted via the parameters in (3.23).

The algorithm combines the advantages given by Theorem 3.2 and Corollary 3.3. This algorithm performs well in simulations, but does not offer theoretical guarantees on the separation between the zonotopic over-approximations.

to provide feasibility guarantees. To this end, we rewrite (2.20) in a simpler form:

u * M P C d+1 = arg min u N 1 i=1 x k+i|k -x1|k 2 Q 1 + N 2 i=N 1 +1 x k+i|k -x2|k 2 Q 2 + . . . + N d+1 i=N d +1 x k+i|k -x3|k 2 Q 3 (3.56) s.t. x k+l|k = Ax k+l-1|k + Bu k+l-1|k , (3.57a)
x k+l|k ∈ X , u k+l|k ∈ U , (3.57b)

x k+l|k / ∈ P(k) (3.57c)
where

P(k) = {P i ∈ P : i s.t. Path(x k , x f ) ∩ Z i = ∅} is

♦

Remark 3.17. The horizons N 1 , N 2 and N 3 are such that:

-x1|k is reachable in N 1 steps from x k ; -x2|k is reachable in N 2 -N 1 steps from x 1|k ; -.... -xd|k is reachable in N d+1 -N d steps from x d-1|k .
Also the number of subsums in the cost function is related with the piece-wise structure of the path which can be generated by ensuring the existence of a corridor and the feasibility of the obstacle avoidance. ♦ Form (3.57) provides feasibility guarantees but does not certify the convergence. In order to impose terminal costs and terminal constraints in such an MPC strategy, one needs to have a global feasibility guarantee for the path planning task. Obviously such a requirement cannot be reached with local information but can receive an efficient solution at the tactical decision stage as discussed in Chapter 4.

Remark 3.18. Constraint x / ∈ P(k), is reformulated with a restricted number (of local) obstacles with the same approach as in Section 2.1. Thus, the complexity of enumeration is substantially reduced, i.e. the number of cells generated by the reduced set of hyperplanes is significantly reduced compared to the MIP one generated in the general case (3.15). ♦

Hereinafter, the number of anti-collision constraints is reduced in term of number of obstacles and the only computational advantage can be obtained by decreasing the complexity of the obstacles' representation. Thus, the zonotopic can provide a practical benefit.

Parametrized zonotopic representation within time-varying multi-obstacle environments.

In many practical applications, as autonomous overtaking [Ballesteros-Tolosana et al. 2017] or coordination of vehicles at intersections [Hult et al. n.d.], we have to deal with dynamical obstacles. Therefore, the resulting hyperplane arrangement is time-varying. Since a re-computation of the arrangement at each time instant may be impractical in terms of complexity, the question to be addressed is:

How can one use hyperplane arrangements for the resolution of the collision avoidance problem in dynamical environments?

A possible solution may be based on parametrized polyhedra [D. Olaru S. D. 2004] and/or by choosing a particular class of polyhedra, zonotopes [START_REF] Kuhn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF]]. The fundamental idea of such an approach is to concentrate the generators of the shape in a compact form (notable by exploiting the property of symmetry -in the zonotopic case) and to separate them from the parameters affected by the time-variation.

As was stated above, the collection (3.15) can evolve due to the dynamical behavior of the environment, or in other words whenever the obstacles change their position in time.

For this reason and in order to avoid a burdensome re-computation, a pre-analysis is necessary.

Let us consider a parameter p ∈ R np and a linear dependence of the polyhedral set description:

H i (p) = {x ∈ R d : s ⊤ i x = r i -(s p i ) ⊤ p}, (3.58)
Remark 3.19. Specifically,in (3.58) we separate the time-varying component of the hyperplanes in (3.15):

r i (k) = r i -(s p i ) ⊤ p(k).
For compactness, k was omitted in (3.58). ♦ Hence, we have a collection of hyperplanes in the lifted space H * = {H * i } ∈ R d+np with:

H * i = {x ∈ R d+np : s ⊤ i (s p i ) ⊤ x p = r i }, (3.59)
The arrangement A(H * ) is composed of cells A(σ * ). The parametrization (3.58) leads to a dependence of the domain of existence of each cell A(σ * ) w.r.t. the parameter p. Therefore, we need to project on the parameter space in order to describe the range of variation of the respective configuration. Conversely, by cutting the polyhedron ( 3.59) at a particular value p leads to a parameterized form of the arrangement A(σ, p).

Recalling the halfspace representation of a zonotope (3.11), we note that the hyperplanes defining a zonotope are strongly dependent on the value of the zonotope center. Therefore, these hyperplanes can be written under the following form: .60) Consequently, the center c is a suitable choice as the parameter p in (3.58) or as one of its component. Noteworthy, in the context of the translation motion, c is the only time-varying characteristic of a zonotope.

H i (c) = {x ∈ R d : s ⊤ i x = r i + s ⊤ i c}, ( 3
Next, we rewrite (3.59) for the case p ← c ∈ R d :

H * i = x c ∈ R 2d : s i -s i x c = r i (3.61) The arrangement A(H * ) with H * (c) = {H * i (c)} is composed of feasible cells A(σ * ), σ * ∈ Σ * .
An interesting feature of this approach is that it allows to compute the domain of existence of each cell A(σ * ) w.r.t. the center position by projecting on the parameter (center) space:

Dom{A(σ * )} = c : ∃x s.t. x c ∈ A(σ * ) (3.62)
All these validity domains can be computed offline and represent polyhedral domains in

R d .
The above reasoning can lay an instrumental role in the next multi-agent scenario. Let us consider two zonotopes:

Z 1 = Z(G 1 , c 1 ), Z 2 = Z(G 2 , c 2 )
with their compact halfspace representations:

Z 1 = {x : H 1 (x -c 1 ) ≤ w 1 }, Z 2 = {x : H 2 (x -c 2 ) ≤ w 2 } (3.63)
A cell A(σ ′ ) from the hyperplane arrangement induced by these two zonotopes is given by:

A(σ ′ ) = {x ∈ R n : H σ i 1 (x -c 1 ) ≤ w σ i 1 , H σ i 2 (x -c 2 ) ≤ w σ i 2 (3.64)
where

σ ′ = σ 1 σ 2 with σ i sign tuple w.r.t. Z i , i ∈ {1, 2}.
Note that the important aspect at the above described multi-agent scenario is not the absolute position of the centers, but their relative position. Thus, we can assume, without loss of generality, that the obstacle Z 1 is fixed and Z 2 is a moving obstacle. Hence we can construct the correspondent collection of hyperplanes Ĥ(c) = { Ĥ}:

Ĥ = x c 2 ∈ R 2n : H σ i 1 0 H σ j 2 -H σ j 2 x c 2 ≤ w σ i 1 + H σ i 1 c 1 w σ j 2 (3.65)
Similarly to (3.62), we compute the domain of existence of each cell A(σ ′ ) in the parameter c 2 space (validity domain) [D. Olaru S. D. 2004]:

Dom(A(σ ′ )) = proj c 2 A(σ ′ * )

Illustrative example

Consider the hyperplane arrangement induced by the following two zonotopes. Using this method for a simple case (d = 2, "rectangular" zonotopes Figure 3.13a) we obtain that only 81 tuples are feasible (∃c 2 such that the tuple is feasible) from 2 8 = 256 possible ones. Moreover, the symmetry of the zonotopes brings a benefit with respect to the number of valid domains Figure 3.13. Figure 3.13a presents one particular position of the obstacles for a given value of the parameters, while Figure 3.13b depicts the number of regions describing the complement of the obstacles as a function of the inter-distance between their centers. . . . . . . . . . . . . . . . . . . . . . . . . . . 128 The previous chapters discussed the interest in both control and robotics communities for the navigation through multi-obstacle environments and its many applications including, among others, monitoring or surveillance, autonomous overtaking or precision agriculture. We recall that, from a mathematical point of view, the main difficulty arises from the non-convexity of the feasible regions in the motion space and consequently in the lack of connectivity in the solution space.

Path planning based on convex lifting

To tackle the navigation problem, most of the approaches divide the problem in two main sub-tasks: path/trajectory generation and trajectory tracking. Moreover, these (sub-)tasks are usually viewed as independent or decoupled problems. For instance, the classical sample-based1 approaches are prone to focus on the first task, neglecting the second one and are, thus, simplifying the problem. On the other hand, the optimizationbased strategies (as those developed in Chapter 3), e.g. mixed-integer formulations or the convexification techniques, potential field methods [Y.-b. Chen et al. 2016] and settheoretic approaches, merge the planning and tracking tasks at the expense of a higher computational complexity. We have shown in Chapter 3 that this last category may have drawbacks especially in the case of cluttered/congested multi-obstacle environments, the computational effort being strongly correlated with the number and structural distribution of the obstacles. Even if symmetry and other topological properties are exploited, the entire optimization-based framework builds on the hypothesis that a solution from an initial point to a target point in the state space exists. This hypothesis is fundamental for the feasibility of the receding horizon optimization and the subsequent developments. Nevertheless, its fulfilment and certification, particularly in a time-varying environment is not guaranteed. The present chapter aims to provide a complete solution to this "initial to final point" path planning problem based on a convex optimization construction (despite the non-convexity of the original problem) by means of a lossless convexification which produces the roadmap of corridors. The properties of the solution will be exploited in order to obtain a certified path from any initial state to a target state, if it exists by a simple graph search.

Needless to say, and regardless of the method, an important aspect is the way in which the environment (obstacles) is represented. A popular practice is the use of convex sets, either polytopes or ellipsoids. Naturally, we maintain the workspace representation employed in Section 3. 3.2, namely polytopic/zonotopic sets.

Philosophically, this chapter is closely related to our results regarding obstacle and collision avoidance presented in Chapter 3 (especially, Section 3.3). The essence of the approach is motivated by the solution introduced in Section 3. 3

which concluded with:

There are two ways to elude the overlapping among the over-approximations. Both share a common feature: introduce separating hyperplane constraints.

Evidently, the following question arises:

Can we pre-determine those separating hyperplanes in order to improve the navigation strategy?

The answer is positive but needs to consider planning and tracking as distinct although interacting tasks. Once the path-planning with collision avoidance is guaranteed, both graph-based and optimization-based approaches can be merged to develop a navigation strategy with feasibility guarantees (this will be discussed in Chapter 5).

The content of the chapter addresses the path-planning level with a primary objective of global feasibility. In other words, our method leads either to a feasible path or to a certificate of infeasibility. This certificate ensures that the environment does not allow the construction of a path between the source and the destination. The dynamical constraints are discarded from the problem formulation as well as the other limitations that may appear in motion planning due to limited steering or energy 2 . The problem becomes a geometrical one and optimality is a secondary objective for the generation of a geometric path. Is worth noting that obtaining a geometric path has the potential to explicitly describe a feasible corridor as in [S. Liu et al. 2017]. We will exploit this property in the practical navigation task.

Before entering in the main development, few positioning remarks can be useful from a mathematical point of view. The solution will exploit the convex lifting notion, which has been previously employed in constrained control and PWA (piecewise affine) control implementations [N. [START_REF] Nguyen | Convex lifting: Theory and control applications[END_REF] but has not been used (to the best of our knowledge) in motion planning or obstacle avoidance. Our work is the first to establish a link between the convex lifting, the polyhedral partitions and the path selection in the navigation space. This versatile optimization-based approach for the construction of a partition starting from the obstacles can be understood as a convexification procedure for the characterization of the non-convex feasible regions in the motion space.

The main contributions in what follows are threefold: i) provide a partitioning of the navigation space based on convex lifting; ii) construct feasible corridors based on a graph of interconnections in the multi-obstacle environment;

iii) propose a path with obstacle avoidance guarantees.

Particularly important from a structural point of view, our solution is not restricted to R 2 , nor R 3 and provides a generic path-generation technique in any finite dimensional space with obstacles.

From obstacle-based partitioning to path generation

This section will first introduce some prerequisites, necessary in handling space partitioning and convex lifting. Next, we will detail the chosen approach, presenting the reasoning behind our choice.

Prerequisities

Consider a finite dimensional output space R d and a finite number of non-overlapping regions P j ∈ Com(R d ), j∈ I = {1, . . . , N o } describing obstacles:

P = No j=1 P j ; P i ∩ P j = ∅, ∀i =j. (4.1)
Through the boundedness assumptions and due to the navigation objectives w.r.t. P j , j ∈ I, the union of these obstacles ( 4.1) will be considered to lie in a bounded cluttered environment X:

P ⊂ int(X) ⊂ R d (4.2)
Consequently, the obstacle-free/feasible domain is C X (P) ≜ X \ P. In order to fix the ideas, Figure 4.1 depicts such an environment (for graphical easiness in R 2 ). Given the collection of obstacles P, a corridor between two points x 0 , x f ∈ int(C X (P)) is enabled by the existence of two continuous functions:

γ : [0, 1] → C X (P) (4.3) ρ : [0, 1] → R >0 (4.4) satisfying γ(0) = x 0 , γ(1) = x f (4.5) γ(θ) ⊕ B 0,ρ(θ) ⊂ C X (P), ∀θ ∈ [0, 1] (4.6)
Based on (4.3)-( 4.6) the corridor is defined as: The above definitions provide a mathematical description for the corridors. The objective of this chapter is their effective construction in a cluttered environment.

Π = {x ∈ R d : ∃θ ∈ [0, 1] s.t. x ∈ γ(θ) ⊕ B 0,ρ(θ) }. ( 4
Remark 4. 1. Corridor construction encompasses classical path planning. Indeed, the corridor can be understood as a compact family of feasible paths between the initial and final point. The selection of a path within a corridor is conceptually simpler in comparison with the (MIP based) direct search of a path avoiding the obstacles. Additionally, given a path within a corridor, the distance to the boundary of the corridor represents a robustness margin for the trajectory planning. ♦

The notion of distance in the context of a multi-obstacle environment is essential for the characterization of the "nearest" obstacles and indirectly leads to the partitioning of the cluttered environment according to the distribution of obstacles.

Definition 4.2. A family of sets {X

i } i∈I verifying: i) X = No i=1 X i , ii) int(X i ) int(X j ) = ∅, ∀i = j ∈ I, iii) P i ⊂ int(X i ), ∀i ∈ {1 . . . , N o }
is called a partition of X induced by the collection of obstacles P. ♦

Proposed Approach

As mentioned above, there exists in the literature a broad variety of sampled-based methods which aim to solve the navigation problem. The main drawback of this kind of methods is the probabilistic completeness [START_REF] Barraquand | A random sampling scheme for path planning[END_REF]]. More specifically, the probability that the algorithm returns a feasible solution tends to 1 if the number of sample points is sufficiently large (approaching ∞), as it was empirically shown in [Hsu, Latombe, and Kurniawati 2007]. However, these probabilistic completeness proofs do not provide any guarantee on the time in which the algorithm finds the optimal path (if there exists one). We place our work in opposition to this line of develoopments and we target as main objective the finite-time guarantee of resolution. The procedure we proposed answers the global feasibility of the path planning and develops a method characterized by completeness, eliminating the risks of an heuristic/unpredictable behavior.

Therefore, in the sequel, we divide the generic navigation problem in three sub-problems characterized in mathematical terms as follows:

P1) Given the collection of convex obstacles P, describe a partition of the cluttered environment around them. Provide a constructive algorithm for the case of the polytopic obstacles. (Section 4.2) P2) Given any two points in the cluttered environment x 0 , x f ∈ int(C X (P)), construct a corridor or provide a certificate of infeasibility. (Section 4. 3.1) P3) Given a non-empty corridor, select according to some performance criteria a continuous path π : [0, 1] → Π guaranteeing collision avoidance, i.e. π(θ)

∩ P = ∅, ∀θ ∈ [0, 1]. (Section 4.2)
Since there are in the literature works which treat in detail problem P3) using à priori defined corridors, see, e.g., [START_REF] Faulwasser | Model predictive path-following for constrained nonlinear systems[END_REF], we focus hereinafter on providing a detailed solution of the problems P1) and P2), respectively.

Typically, the motion planning methods employ a path planner and rely on a graph as a result of the workspace partitioning, e.g., by using a grid of square/cubic cells [X. Wang et al. 2015] or Voronoi diagrams [START_REF] Sugihara | Approximation of generalized Voronoi diagrams by ordinary Voronoi diagrams[END_REF]]. Therefore, before the main development, it is worth showing that the problems P1)-P2) are not trivial. An intuitive, but not completely adequate idea for partitioning of the cluttered environment is related to a Voronoi-like construction [START_REF] Sugihara | Approximation of generalized Voronoi diagrams by ordinary Voronoi diagrams[END_REF]Aurenhammer 1991].

We recall that a Voronoi region V i corresponding to a collection of points c i ∈ R d is described by:

V i = x ∈ R d | d(x, c i ) ≤ d(x, c j ), ∀j = i (4.
c 2 c 3 c 4 c 5 V 6 V 5 V 2 V 3 V 1 V 4 Figure 4.4.: An example of classical Voronoi diagram in R 2 .
As can be seen in Figure 4.4, the Voronoi graph is a polyhedral partition of the workspace which satisfies i)-ii) (Definition 4. 2), but relative to a given set of points and not w.r.t. regions as will be the case of the obstacles. Indeed, the points in the classical Voronoi partitioning can be replaced by the collection of obstacles. Hence, the construction of such diagrams relies on the selection of certain points in the workspace: either mass (or Chebyshev) centers of the obstacles or a set of points approximating the boundary of the obstacles [Rubens Junqueira Magalhaes Afonso, R. K. H. Galvao, and Karl Heinz Kienitz 2013], as depicted in Figure 4. [START_REF]Construction of a partitioning for the cluttered environment X ⊃ P[END_REF]. Nevertheless, the resulting diagrams either do not satisfy the characteristics of a partition (as in Definition 4.2) or the resulting graph needs an additional processing step in order to remove redundant nodes and edges and whose overall complexity is not negligible. Moreover, these solutions are generally defined only for the planar case.

Convex lifting as a tool for space partitioning

As mentioned in the introduction of this chapter, the problem to be solved is essentially non-convex and retrieving a convex or convexity-based counterpart will represent an important step forward from the constructive point of view (practically eluding the MIP formulations). The convexification is a well understood concept and it can be seen as a transformation of the original problem in a different space of design parameters, where the design enjoys convexity properties. Finding this novel space of parameters is non-trivial task, often concerns conceptually different (larger) dimensions and can involve approximations. When the convexification is done without loss of precision in terms of the set of feasible solutions or optimality, the lossless convexification term is used. In the present framework, the design will involve a convex lifting, which represents a search for a convex construction in a higher dimensional space, retrieving by projection desirable properties. This section presents in-depth this essential technique for an obstacle-based space partitioning, addressing, at the same time, the particularities that may hinder the applicability of the proposed navigation strategy.

Convex Lifting. Background.

From a geometrical perspective, the lifting can be viewed as the inverse operation with respect to the orthogonal projection. Besides the control theory, this technique was exploited in various fields, such as: signal processing, mechanics or geometry. As will be emphasized by the mathematical definitions, the operation consists in lifting a (finite) collection of polyhedra onto a higher dimensional space. More precisely, it means the construction of a convex surface in R d+1 whose orthogonal projection on R d is methodically related in terms of containment or polyhedral partition to the original collection of polyhedra.

Up to now, the control literature resort to the convex lifting in the PWA/constrained control implementations, exploiting the capability of the lifting techniques in resolution of the inverse parametric linear/quadratic programming problem [START_REF] Gulan | Efficient embedded model predictive vibration control via convex lifting[END_REF]N. A. Nguyen, Sorin Olaru, et al. 2017]. In such problems, the collection of polyhedra is representing a polyhedral partition (or cell complex5 ). However, the convex lifting has a conceptually larger applicability.

Despite the fact that the problem of the existence was intensively studied by the mathematics community with a plethora of structural results [Aurenhammer 1991], there exist only a few works, which provide constructive solutions for the convex lifting. Hence, before introducing our contribution, we briefly recall the existing results in order to highlight the differences between the approach therein and our development. For further details regarding the concept (of convex lifting) and its implications (in quadratic programming) and/or for completeness, the interested reader is referred to [N. A. Nguyen 2015]. Definition 4.4. Given a cell complex {X i } i∈I of a polytope X with I a finite subset of N, the function z : X → R is called a convex PWA lifting with:

z(x) = a ⊤ i x + b i , x ∈ X i , and
a i ∈ R d , b i ∈ R, ∀i ∈ I (4.9)
if z(x) is continuous and convex over X.

Remark 4. 4. The lifting z(x) satisfies the continuity condition if For each pair of neighboring regions (X i , X j ),

a ⊤ i x + b i = a ⊤ j x + b j , ∀x ∈ X i ∩ X j as long as X i , ∀i ∈ I is a compact (polyhedral) set. Remark 4.5. The function z(x) being convex over {X i } i∈I means that z(x) > a ⊤ j x+b j , x ∈ X i \ X j , ∀i = j ∈ I ♦ Using Definition 4.
add continuity conditions ∀v ∈ V(X i X j ):

a ⊤ i v + b i = a ⊤ j x + b j (4.10)
add convexity conditions ∀v ∈ V(X i ) \ V(X j ): The feasibility of the optimization problem (4.12) is a necessary and sufficient condition for the existence of the convex lifting corresponding to the given partition.

a ⊤ i v + b i ≤ a ⊤ j x + b j + ξ (4.
Remark 4. 6. The existence of a convex lifting is guaranteed by the feasibility of (4.12) and this property is not affected by the value of the constant ξ > 0. This constant has exclusively a numerical conditioning influence (by avoiding the singularities through its positiveness and which should not be chosen large in order not to affect the numerical representations). ♦

Algorithm 4 is based on the vertex representation of the partition components. Note that the vertex enumeration may imply an important computational effort (e.g., the algorithm proposed in [START_REF] Avis | Reverse search for enumeration[END_REF]) is strongly dependent on the dimension d, number of vertices and number of facets of the given polytope. This drawback can be alleviated by using an algorithm based on the halfspace representation [N. A. Nguyen, Sorin Olaru, et al. 2017]. We do not dwell here on this computational part of the construction and preserve for the clarity of presentation a dual (generators) representation. Moreover, for the problem at hand (navigation in cluttered environments) this issue is not crucial, in most of the cases the enumeration is not necessary, the obstacles being already described explicitly as a collection of vertices.

z(x) > a ⊤ j x + b j , x ∈ X \ X j , holds. (4.15)
Proof. i) Since the sets X i ⊂ X are compact and cover X without gaps ( X i = X), there are pairs (i, j) with X i ∩ X j = ∅. For each x ∈ X i ∩ X j from the uniqueness of z(x) w.r.t parameters a i and b i it follows a ⊤ i x + b i = a ⊤ j x + b j and thus the lifting (4.14) is continuous.

ii) A necessary and sufficient condition of convexity is that inequality

λz(x 1 ) + (1 -λ)z(x 2 ) ≥ z(λx 1 + (1 -λ)x 2 )
holds for any pair

x 1 , x 2 ∈ X and λ ∈ [0, 1].
Without loss of generality let us assume that x 1 ∈ X i , x 2 ∈ X j and λx 1 +(1-λ)x 2 ∈ X k . Thus, the inequality becomes:

λ(a ⊤ i x 1 + b i ) + (1 -λ)(a ⊤ j x 2 + b j ) ≥ a ⊤ k (λx 1 + (1 -λ)x 2 ) + b k . Rewriting the right-hand side as λ(a ⊤ k x 1 + b k ) + (1 -λ)(a ⊤ k x 2 + b k )
and applying (4.15) shows that the inequality holds, thus (4.14) is convex.

Theorem 4.1. A piecewise affine lifting for a collection of obstacles P = No j=1 P j with int(P i ∩ P j ) = ∅, ∀i = j is continuous and convex iff (a i , b i ) satisfy:

a ⊤ i v + b i ≥ a ⊤ j v + b j + ϵ, ∀v ∈ V(P i ) ∩ int(X i ), ∀i = j, (4.16a)
a ⊤ i v + b i ≥ a ⊤ j v + b j , ∀v ∈ V(P i ), ∀i = j, (4.16b)
a ⊤ i v + b i ≤ M, ∀v ∈ V(P i ), ∀i. (4.16c)
In (4.16) ϵ, M > 0 are suitably chosen and V(P i ) denotes the collection of extreme points of P i .

Proof. Let us recall that P i ⊂ X i are convex. Thus, with a suitably chosen ϵ, and checking (4.16c), (4.16b) guarantees the inequality:

a ⊤ i x + b i ≥ a ⊤ j x + b j , x ∈ X \ X j , (4.17)
Whenever the obstacles are disjoint, i.e. P i ∩ P j = ∅, the vertices of the obstacle are strictly inside the interior of the set X i , the inequality holds strictly as in (4.16a).

All these conditions imply the convexity and by adding the continuity properties (through Lemma 4.1) for the lifting (4.14), the proof is complete.

Taking P i as polyhedral sets (i.e, having a finite number of extreme points), allows to obtain the lifting as the result of the following convex optimization problem: .18) Based on the solutions of ( 4.18), we define the following "d+1"-dimensional polyhedron:

min a i ,b i No i=1 a i b i ⊤ 2 2 s.t. (4.16a) -(4.16c) hold. ( 4
L = x z ∈ R d+1 : a ⊤ i -1 x z ≤ -b i , i ∈ I . (4.19)
Projecting the facets of L on X a polyhedral partition can be obtained. {X i } i=1:No .

Corollary 4.1. The polyhedral partition {X i } i=1:No obtained by convex lifting for disjoint obstacles has the following properties:

i) P i ⊂ int(X i ), ∀i, ii) X i ∩ P j = ∅, ∀j = i.
Proof. Projecting the facets 6 of L on X, we obtain the polyhedral partition {X i } i=1:No :

X i = proj(F d-1 i (L), X ), ∀i, (4.20)
which allows proving the properties: i) The feasibility of ( 4.18) and the vertex-representation of P i , ∀i, implies that for all x ∈ P i we have:

x z(x) = x a ⊤ i x + b i ∈ L.
Thus, the lifting

corresponding to P i is included in F d-1 i (L).
By projecting over X, it directly leads to inclusion P i ⊂ X i . ii) Let us suppose that there is a y ∈ X i ∩ P j . According to i) we have P j ⊂ X j . That leads to y ∈ X j and, from feasibility of (4.16a):

a ⊤ j y + b j ≥ a ⊤ i y + b i + ϵ. (4.21)
Moreover, y ∈ X i ∩ P j means that y ∈ X i ∩ X j . Thus, via (Lemma 4.1), we have:

a ⊤ i y + b i = a ⊤ j y + b j . (4.22)
Combining ( 4.21) and ( 4.22), gives ϵ ≤ 0, which contradicts that ϵ > 0.

Algorithm 5 Construction of a partitioning for the cluttered environment X ⊃ P, Input: P = No j=1 P j with P i ∩ P j = ∅, ∀i = j, X and a given constants M > 0 . Output: {X i } i∈I a polyhedral partition of X 1: For each pair of obstacles (P i , P j ) impose the constraints:

a ⊤ i v + b i ≥ a ⊤ j v + b j + ϵ, ∀v ∈ V(P i ), ∀i = j, (4.23a)
a ⊤ i v + b i ≤ M, ∀v ∈ V(P i ), ∀i. (4.23b)
2: Solve the convex optimization problem: Compute the projections X i of the facets of L on X Remark 4.9. In contrast with Algorithm 4 (i.e., the constructive solution in [N. A. Nguyen 2015]), the method summarized in Algorithm 5 is characterized by the following aspects:

min a i ,b i No i=1 a i b i ⊤ 2 2 s.t
6 Facet F d-1 i (L) is the (d -1
)-order face of polyhedron L, i.e., the inequalities describing L remain the same except the i-th which is converted to an equality.

Non-liftable partition

For the convex-non-liftable cluttered environments, the solution can be analysed in the perspective of the structural result presented in [N. [START_REF] Nguyen | Fully inverse parametric linear/quadratic programming problems via convex liftings[END_REF]], which can be summarized as follows:

For any polyhedral partition there always exists one subdivision such that the internal boundaries of this partition are preserved and the new partition is convexly liftable.

Transposed in the present framework of obstacles included in the partition's regions, it follows that obstacles can be subdivided and represented as collections of convex subsets, thus enabling convex lifting.

The next result is formally stating the existence of a finite refinement of the cluttered environment leading to a convex lifting and the associated polyhedral partition.

Theorem 4.2. Given a collection of obstacles P = No j=1 P j with int(P i )∩int(P j ) = ∅, ∀i = j ∈ Z [1,No] and P ⊂ X, there exists at least one representation of these obstacles as finite unions of convex polyhedral sets P = Ño j=1 Pj = P with int( Pi ) ∩ int( Pj ) = ∅, ∀i = j ∈ Z [1, Ño] which guarantees the existence of a polyhedral partition of X = No k=1 Xk such that for any j ∈ {1, . . . , Ño } there exists one and only one k ∈ N [1,No] satisfying Pj ⊂ Xk and int(

Pj ∩ Xl ) = ∅, ∀l ∈ N [1, No] \ {k}.
Sketch of Proof : Since the obstacles are polytopic sets, hyperplanes intersecting an obstacle allow the cutting of the respective obstacle and and its representation in terms of a union of two convex sets. Thus, we obtain a decomposition of the initial setting. According to the theorem in [N. [START_REF] Nguyen | Fully inverse parametric linear/quadratic programming problems via convex liftings[END_REF]], the hyperplane arrangement generating the decomposition can be chosen such that the corresponding partition, induced by the new obstacle setting, is convexly liftable. The finite number of cutting is related to the notion of k-stress and can be proved using the hyperplane arrangements which are known to lead to convex liftable partitions. Corollary 4.2. Given a collection of obstacles P = No j=1 P j with int(P ) i ∩ int(P ) j = ∅, ∀i = j with P ⊂ X, and their associated sets of supporting hyperplanes H(P i ), ∀j ∈ {1, . . . , N o }, the decomposition of P by H(P) induces a convexly liftable polyhedral partition.

Proof. The reasoning relies on Theorem 4.4.10 from [N. A. Nguyen 2015] and on Theorem 3 in [START_REF] Aurenhammer | A criterion for the affine equivalence of cell complexes inR d and convex polyhedra inR d+ 1[END_REF]], which guarantee the feasibility of the convex lifting (4.18) for any partition resulted from a hyperplane arrangement technique. Definition 4.3. Theorem 4.2 guarantees the existence of a liftable partition starting from a finite decomposition of the obstacles. The particularity of the resulting partition is that there exist at least one index i ∈ N [1,No] such that int(P i ∩ Xj ) = ∅ and int{P i ∩ Xl } = ∅ for some j, l ∈ N [1, Ño] and Xj = Xl .

Corollary 4.3. Given a collection of obstacles for which no convexly liftable partition exists according to

Proof. Let us consider P i ∈ P as a union of two convex sets from P, i.e., P i = Pj ∪ Pl . Computing the lifting associated with P we obtain the polyhedral partition X for which we have the bijective correspondence Pk ⇐⇒ Xk , ∀k. When we report the resulting partition X to the original set of obstacles, the partition cell corresponding to P i is X i = Xj ∪ Xl , and P i intersect both Xj and Xl . This reasoning can be generalized for P i as a finite union of more than two sets.

As a consequence of Corollary 4.3, some boundaries of Xi , i ∈ {1, . . . , Ño } are intersecting the original obstacles. In order to resort this property one simply needs to merge the original obstacles and their corresponding sets in the partition. Explicitly, for each P i ∈ {P 1 , . . . , P No } it can be obtained a set of indices I i = {1 ≤ j ≤ Ño | int(P i ∩ Xj ) = ∅} and finally a corresponding region X i = j∈I ⟩ Xj . Therefore, by relaxing the bijective correspondence between obstacles and partition cells and the convexity of these cells Corollary 4.3 guarantees the existence of a partition which can be employed in a path planning strategy. The difference comes from the fact that instead of using points on the boundaries of the obstacles as in Example 4.2, we split them in finite unions of convex sets.

Remark 4.14. Practically, the technique based on Corollary 4.3 represents only one of the possible solution for the identification of the non-convexly liftable cell complex/polyhedral partition. However, this adjustment will increase the complexity of the implementation. Therefore, the complete refinement (inherited from the hyperplane arrangement) may not be necessary in practical applications and only represents an upper bound for the feasibility of present convex-lifting based methodology.

In order to avoid the unnecessary increasing of the complexity (Remark 4.14) we propose Algorithm 6. This procedure bypasses the splitting of all obstacles and concentrates only on those obstacles which lead to non-liftability. In other words, we can iteratively7 run the Algorithm 5 until the optimization problem (4.14) becomes infeasible. To construct a partition, we apply the splitting procedure only for the obstacles which are not in the set generating the last convexly liftable partition.

Algorithm 6 Treatement of convexly non-liftable obstacles via Corollary 4. 31: Find the sets of supporting hyperplanes H(P i ), ∀i ∈ {1, . . . , N o } 2: Construct iteratively the "maximal" set S ⊂ P such that a convex liftable partition { Xi } i=1:No and P = S exists. Find the hyperplanes from H(S) which intersect P i : { Hk } 1≤k

5:

Split P i using { Hk } 1≤k :

P i = P i \ { Hk ⊕ B 0,ϵ } 1≤k , 0 < ϵ ≪ 1 (4.25) 6:
P = P ∪ {P i } 7: end for 8: Compute the partition for P 9: Extract the cells Xi , ∀i ∈ {1, . . . , Ño } 10: Compute the partition cells X i , ∀i ∈ {1, . . . , Ño } such that P i ⊂ X i , P i ∈ P Example 4. 6. Let us consider the cluttered environment depicted in Figure 4.13a and containing six obstacles. If we apply Algorithm 5 for this environment, the lifting (4. 14) is infeasible, and, thus, the partition is convexly non-liftable. The highlighted obstacle (yellow) in Figure 4.13d is the obstacle that can be seen as the source of the non-liftablity. Employing Algorithm 6 we split this particular obstacle in order to obtain a liftable partition as the one in Figure 4.14a.

Remark 4.15. As can be seen in Figure 4.14, Algorithm 6 represents a relaxation of the theoretical upper bound given by Theorem 4.2. However, the number of local cuttings can be high, depending on the neighbours and the geometry of the partition and can subsequently increase the complexity of last steps in Algorithm 6. That is, the resulting partition cell corresponding to the divided obstacle impacts the post-processing step for the path planning procedure and can be further improved.

♦

Let us focus next on a variant of Algorithm 6 which selects only the supporting hyperplanes associated to the neighbors of the problematic obstacle (from the liftability perpective). As depicted in Figure 4.14b, such an approach leads to a considerably reduced number of cuttings and, implicitly, a diminution of the number of convex region composing the cell associated to the divided obstacle. We do not propose a new specific algorithm implementing this strategy but remark that the approach leads to a simple modification in Step 4 of Algorithm 6 where only the hyperplanes of the neighbour obstacles are used for cutting.

This idea can be further refined and we can replace the pairs of hyperplanes by the separating hyperplane obtained from the iterative lifting procedure. Since we have the separation between each pair of obstacle from the "maximal" liftable subset S, we can use the facets of the partition in order to proceed with the splitting. Following this approach, one can operate the cutting as in (4.26) and update the collection of obstacles before computing a new partition.

Algorithms 6 and 7 lead to convex liftable solutions but it results into the cells that need to be merged in order to cover exactly the original set of obstacles, i.e.,

X i = j∈I X ij .
For such a region X i the convexity is not guaranteed. However, the resulting partition satisfies the conditions in Definition 4.2 and, as we will see in the next subsection, solves the path planning problem without needlessly increasing the complexity.

Algorithm 8 Treatement of convexly non-liftable obstacles via hybrid method II 1: Find the maximal subset of obstacles S ⊂ P generating a convex liftable partition { Xi } i=1:No and P = S 2: for P i / ∈ S do 3:

Find the cells { Xj } j̸ =i such that P i ∩ Xj = ∅ 4:

Find the facets of { Xj } j̸ =i which intersect P i : { Fk } 1≤k

5:

Split P i using { Fk } 1≤k :

P i = P i \ { Fk ⊕ B 0,ϵ } 1≤k , 0 < ϵ ≪ 1 (4.27) 6:
P = P ∪ {P i } 7: end for 8: for X k such that ∃ I k ⊂ {1... Ño } with Pi ⊂ X k , ∀i ∈ I k do 9:
Run Algorithm 5 for X = X k and P = The above solutions (Algorithms 6-8) for the convexly non-liftable cases relax the convexity condition on the cells of partition. Thus, a natural question is whether or not there exists a method able to identify a polyhedral partition w.r.t. to the considered polyhedral obstacle setting. In order to tackle this ultimate challenge, we propose an approach which relies on selecting those pairs of neighboring obstacles which are problematic from the lifting perspective, and merge each pair into an altered obstacle. For instance, the newly formed component replacing the pair (P i , P j ) is a convex hull P ij = Conv{P i ∪ P j } which verifies the condition P ij ∩ P k = ∅, ∀k / ∈ {i, j}. Thus, we have a modified collection of obstacles, P, which contains P ij instead of P i , P j . We employ the lifting procedure Algorithm 5 for P, and, further, the cell X ij corresponding to P ij can be partitioned w.r.t. the obstacles P i , P j . But P i ∪ P j is nothing else than a particular sub-set of obstacles that can be tackled by convex lifting in itself, independently. The entire procedure is summed up in Algorithm 9.

Remark 4.17. Algorithms 6-8 proceed by sub-cuttings in order to enhance the convex liftability while Algorithm 9 proceeds by merging and outputs convex partitions of a cluttered environment. ♦

Remark 4.18. The replacement of non-convex obstacles by a union of convex obstacles enables the construction of a partition as in Section 4.2. However, the existence of nonconvex obstacles lead can lead to infeasibility in path planning procedure (as will be discussed in Section 4.3) Remark 4.19. As can be seen in the case of the cluttered environment in Figure 4.18a, in the construction of a partition for a non-convex obstacle setting, the importance of the choice of convex decomposition and the "trimming" procedure is increasing when compared to the convexly non-liftable case, from the perspective of path planning.

From paths to corridors

Geometric path generation based on polyhedral partitions

The partitioning introduced in the previous section induces a graph structure which allows to compute an obstacle-avoiding path. .) to associate weights to the nodes (f : N → R), not to the edges. Consequently, the classical shortest-path algorithms are slightly modified but can still be applied in order to select a path in between an initial and a final node. ♦

Starting from the partition {X i } i=1:No of the workspace X, our goal is to construct an associated graph in order to generate feasible paths through X. For the construction of the graph, we have to select the nodes, the edges and the associated weights from the constructive parameters of the compact sets X i (vertices and faces). The principle is simple and relies in the use of the faces of the polyhedra in the partition. It is known that these faces form a lattice (combinatorial structure) for a polyhedron. In Table 4.2 we delineate the existing possibilities for selecting graph components. While for d = 2 the differences among the alternatives are negligible, for higher dimensions the choice of the graph representation may have a significant effect on the characteristics of the corridor (4.7).

Domain of f

N E 1 E F 0 (X i ) F 1 (X i ) 2 F d-2 (X i ) F d-1 (X i ) 3 N F 1 (X i ) F 0 (X i ) 4 F d-1 (X i ) F d-2 (X i ) Table 4.1.: Alternative selection of graph triple (N , E, f ). Proposition 4.2. Γ {F d-2 (X i )} X i ∈X , {F d-1 (X i )} X i ∈X , f is a connected planar graph.
Proof. As {X i } i=1:No is a polyhedral partition of X, it satisfies the properties from Definition 4.2, i.e. the polytopic regions X i share common boundaries and cover X without gaps. Since the regions X i are polytopic, the facets F d-1 (X i ) are disjoint except on their where I f i = 1, . . . , N f i with N f i the number of facets of X i . Similarly, we obtain (j ⊥ , x ⊥ )(x f ). These two auxiliary nodes are added to the graph by linking them with the incident nodes of the containing edge and removing this edge. A new graph Γ1 (x 0 , x f ) which preserves the properties of Γ 1 is thus obtained. A graph search algorithm (e.g. Dijkstra's Algorithm [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF]) is employed and the shortest path between the nodes induced by x 0 and x f is obtained.

Problem (4.30) is non-convex: it is actually a MINLP problem, but its complexity is polynomial in the number of vertices of the graph and as such particularly attractive for path planing in comparison with the NP complexity of the MIP problems described in Chapter 2. Moreover, for d = 2 case the solution can be iteratively computed. Thus, we find the projection of x 0 ∈ X i on each of the facets (edges) of X i and we verify whether the segment linking x 0 and its projection intersects the obstacle P i . Afterwards, we select the projections corresponding to the segments which do not collide with the obstacle and, finally, we find the minimum among them. We summarized the procedure in the following algorithm (which can be readily extended for higher dimensions):

Algorithm 10 An iterative solution for the MINLP (4.30)

1: Find X i s.t. x 0 ∈ X i ; 2: for each F j (X i ) do 3:
Find x ⊥ j = proj(x 0 , F j (X i ));

4:

P = Conv({x 0 , x ⊥ j }); 5: if P ∩ P i = ∅ then 6: d j = dist(x 0 , x ⊥ j ); 7: else 8:
d j = Inf; 9:
end if 10: end for 11: [d min , j min ] = min j d j ;

12: x ⊥ = x ⊥ j min ;

Remark 4.21. For further use, we denote the shortest path through the graph between x i and x f as Path(x 0 , x f ) = (x 0 = x i , x1 , . . . , xn , xn+1 = x f ). This represents an ordered set of points where no segment defined by a pair of consecutive points cuts any of the obstacles. It is not a path in the sense stated in problem P3), but is a sufficient condition for the existence of a corridor (4.7). ♦ Proposition 4.3. Any polyhedral partition {X i } i=1:No provides a corridor [START_REF]Treatement of convexly non-liftable obstacles via a hybrid method[END_REF] for a given pair x 0 , x f ∈ int(C X (P)).

Proof. As Proposition 4.2 states, there exists a connected graph Γ induced by the partition {X i } i=1:No . By connecting any two points to the graph, the connectivity is preserved. Thus, by using a search algorithm for the extended graph Γ we identify a continuous piecewise affine function γ generated by the edges composing Path(x 0 , x f ). Moreover, we consider ρ as a piecewise constant function defined for each edge of Path(x 0 , x f ) as the minimum Hausdorff distance between the edge and the obstacles around it. The set Π, described by functions γ and ρ is a corridor in the sense of (4.7).

Algorithm 11 Continuous path generation

Input: Γ(x 0 , x f ) Output: a path π : [0, 1] → X with π(θ) ∩ P = ∅, ∀θ 1: Find Path(x 0 , x f ) using Proposition 4.2.

2: Determine γ(θ) from Definition 4.1 for Path(x 0 , x f ).

3: Find ρ(θ) such that:

ρ(θ) ≤ min P i ∈P d H (P i , γ(θ)), ∀θ ∈ [0, 1]; (4.31)
4: Find a path π(θ) inside the corridor characterized by γ(θ), ρ(θ) as:

π(θ) = arg min ν C(ν) s.t ν(θ) ∈ γ(θ) ⊕ B 0,ρ(θ) , ∀θ ∈ [0, 1].
Each step of Algorithm 11 represents a further refinement along the trajectory generation procedure.

Step 1 provides the graph structure of the workspace, from which, in Step 2, γ(θ), a PWA (continuous) function is given (under the construction proposed in this section: a collection of edges). To these, in Step 3 is attached, ρ(θ), a continuous width function which provides a measure of the acceptable deviation from the nominal γ(θ).

Step 4 replaces the feasible solution γ(θ) by an optimization-based selection π(θ). The cost C(ν) is often taken as the path length but may be replaced/augmented by a cost describing energy/performance/smoothness. For illustration purposes we revisit the obstacle collection shown in Figure 4.9 to which we apply Algorithm 11.

First, we construct the associated graph (depicted in Figure 4.19a) and find a path γ, as shown in Figure 4.23. Next, we provide an approximation of the corridor width ρ (gray area in Figure 4.23 is the corridor, as defined in (4.7)). To compute the corridor width we sampled the continuous parameter θ and introduced it in (4.31). We sidestep Step 4 of the algorithm by choosing π = γ. This path is provided as reference to a standard path tracking mechanism which (green line with diamond markers) is shown to respect the constraints (there is no intersection with the obstacles and the destination is successfully reached). For illustration of the ultimate path tracking task we considered a standard double integrator dynamic and applied an MPC (Model Predictive Control) strategy.

We depicted in the following table the corresponding runtimes (measured in seconds) for the d = 2 case. We focus on the d = 2 case due to its relevance in the motion planning context. Thus, we have considered 10 randomly generated10 multi-obstacle environments for each obstacle count and we indicate the computational time for each procedure involved in the motion planning. As can be seen in the table and in Fig. 4.21, the computational time corresponding to LP (4.18) and to the construction of the graph Γ are strongly dependent on the number of the obstacles. Regarding the resolution of the MINLP (4.30), the influence of #P is a bit smaller, whereas for the Dijkstra's Algorithm is not relevant 4.30) t Dijkstra Let us recall the example illustrated in Figure 4.18. As stated, by replacing the nonconvex obstacles with their corresponding convex subsets we obtain the partition depicted in Figure 4.22a. Further, we are able to construct the graph as in Figure 4.22b, which need to be adjusted. That is, we eliminate from the graph those nodes and edges intersecting the initial non-convex obstacle setting and, as well the remaining isolated component.

#P t LP t Γ t (
In this manner, we obtain the graph depicted in Figure 4.22c. One can remark that the procedure identifies all the feasible passageways for the non-convex multi-obstacle environment, and, thus, the graph can be further employed in a trajectory generation strategy.

Corridors as patchy polyhedra

Remark 4.23. The corridor Π can be written as

Π = Nc i=1 Π i with Π i = {x ∈ R d : ∃ θ ∈ [0, 1] s.t. x ∈ γ i ( θ) ⊕ B 0,ρ i ( θ) } (4.32)
where γ i (0) = xi-1 and γ i (1) = xi . ♦

For illustrating the construction of the corridor, we revisit the obstacle collection shown in Figure 4.9. Hence, we provide an approximation of the corridor width ρ (red area in Figure 4.23 is the corridor defined in (4. 7)). To compute the corridor width we sampled the continuous parameter θ and used the corresponding values in (4.31).

Remark 4.24. Since we deal with a polytopic representation of the environment, for further use we consider the polytopic underapproximation of the corridors (blue area in Figure 4.23). The construction of the corridor is based only on the topological characteristics, but the method proposed here is not unique. For example, instead of using the Hausdorff distance, the method in [S. Liu et al. 2017] relies on finding the maximal ellipsoid including a segment from a given path and not intersecting the obstacles.

Conclusions

This chapter presented a constructive solution for the generation of a path between two points in an environment obstructed by multiple obstacles in a d-dimensional space. The The present chapter provides a control-based perspective of the navigation strategy. More precisely, relying on the ideas and concepts proposed in Chapter 4, we proceed further in providing global feasibility guarantees and developing an effective control strategy.

Navigation in cluttered environments

Basically, in a first stage (covered by Chapter 4) we neglected the dynamical constraints and the physical limitations that may appear in the motion planning in order to generate a feasible geometric path. As stated, the resulting path ensures the avoidance of obstacles and has the potential to explicitly describe a feasible corridor similar to concepts handled in the literature (e.g. [S. Liu et al. 2017] or in [START_REF] Faulwasser | Nonlinear model predictive control for constrained output path following[END_REF]). At a second stage, using the geometric path and the corridor as starting points, we find some appropriate trajectory respecting the agent's dynamics and constraints using a MPC strategy.

Additionally, this chapter is concerned with path planning in congested time-varying environments. Our approach exploits the partitioning of the feasible space according to the obstacles' distribution via convex lifting, see Section 4.2. While the procedure is well established for static environments, it needs adaptations in a dynamical context. The second part of this chapter covers the adaptation of the generic methodology for timevarying multi-obstacle environments providing the techniques necessary for the case of a a priori known dynamic behavior. Not in the least, we illustrate the results with geometrical examples.

Moreover, the third part of the present chapter pertains to navigation of a multi-agent system within a cluttered environment. Thus, we introduce a methodology tailored to the presence of a team of agents, based on the techniques presented in the previous chapters (Chapters 2 -4).

The main contribution is threefold: i) a detailed analysis of the recursive feasibility of a corridor-constrained MPC;

ii) a generic navigation strategy for cluttered environments with feasibility guarantees;

iii) an investigation on the evolution of the space partitioning (and, implicitly, of the convex lifting) corresponding to a time-varying multi-obstacle environment.

Motivation

The main motivation for this chapter's work comes from the approaches presented in the previous chapters, each one with their own advantages and inconveniences. Note that, for obtaining a collision free trajectory for an agent through the cluttered environment, we employed different strategies, all within an MPC framework:

1. Classical MIP-based trajectory planning: Specifically,the MPC problem (2.20) with the obstacle avoidance constraints was employed using the MIP formalism (2.7).

As stated in Section 2.1, the shortcomings of this method are: i) it becomes numerically complex for large numbers of obstacles;

ii) it needs a large prediction horizon in order to attain its destination/target.

However, the main advantage is given by the following theoretical guarantee: if there exist a feasible trajectory, then there exists a value of the prediction horizon such this feasible solution is reached. ii) we have a manageable complexity of the obstacles due to the zonotopic framework as in Section 3. 3.2 iii) the convergence to the target point is guaranteed, only if a piece-wise linear path is available and the MPC is successively tracking a sequence of d nodes on this path.

The inconveniences are coming from the pre-processing of the obstacles approximation and the dependence on the path-planning mechanism before the effective trajectory tracking.

Geometrical path computation via space partitioning:

Based on space partitioning (via convex lifting, see Chapter 4), we compute a geometrical path further used as a reference for a MPC trajectory tracking strategy. The asset of this approach lies in that the the MPC strategy is unburdensome (convex formulation), but the optimality and convergence guarantees cannot be readily certified. In short, the strategy relies on the path planning for the obstacle avoidance while the navigation concentrate on the performance of the tracking.

Corridor-constrained trajectory planning:

The approach is based on Algorithm 11 in Section 4. 3.1, where the fourth step is replaced with an MPC algorithm enforcing corridor-constrained trajectory 1 . The benefits of this approach are: i) the inherent non-convex constraints are replaced with convex ones;

ii) the complexity of the convex MPC strategy is reasonable as long as a fixed structure is used.

Nevertheless, the absence of feasibility and convergence guarantees (due to the length of the path in comparison with the prediction horizon which prevents from employing terminal constraints) for this approach is troublesome and it has to be thoroughly addressed.

An illustrative example can provide a qualitative comparison among the above strategies and paves the way for the forthcoming improvements of this chapter.

Illustrative example:

Consider the cluttered and complex environment depicted in Figure 5.1a, containing N o = 33 obstacles, and an agent described by the dynamics (5.2) in R 4 with:

A = O 2 I 2 O 2 -µ M I 2 , B = O 2 M -1 I 2 , ( 5.1)
where µ = 3 and M = 60. The agent's state is composed of position and velocity components, grouped in the state vector x = p x p y v x v y ⊤ , whereas the input is given by the accelerations u = a x a y

⊤

. Both state and input are constrained to : In Table 5.1 we delineate some noteworthy computational characteristics: N goal -the number of discrete time instants to attain a neighborhood/vicinity of the final point, t goal -the total time to compute the trajectory, and ℓ t (m) -the length of the resulting trajectory. Since the comparison is a qualitative one within an MPC framework, note that all the MPC parameters, employed in all four above-detailed approaches, are identical, excepting length, but it functions significantly better from the computational perspective due to the convex formulation.

X = {x ∈
We remark that t goal corresponding to the last three approaches is quite similar, but there is a substantial variation regarding the length of the trajectory. Within the four examples, the last approach performs better in the sense of the parameters delineated in Table 5. 1.

Nevertheless, the choice of the tuning parameters needs an application-dependent analysis.

To sum up, in the absence of strong results regarding the feasibility and convergence, the classification of the approaches is not relevant as long as they come with different certificates. Hence, in the following section, we aim to develop a method enhancing the corridor-based (the fourth) approach and providing the missing theoretical guarantees which can recommend it from performance, guarantees and computational perspectives..

Relay MPC. Feasibility guarantees

This section introduces the formulation of an MPC controller governing the motion of a mobile agent within polyhedral corridors and its proof of recursive feasibility and convergence. The strategy proposed in the sequel exploits the existence of a feasible path and its associated corridors resulted from the partitioning of the environment. First we recall the necessary concepts and, next, we present the formulation associated to a single compact corridor (a segment from (4.7)), see Figure 5.4. Then the extension to the entire corridor (4.7) is tackled.

Prototype MPC with recursive feasibility guarantees

In what follows, the aim is to obtain a collision-free trajectory inside a corridor Π i for an agent described by the following LTI dynamics:

x k+1 = Ax k + Bu k , ∀k (5.2)
with x k ∈ R d the state vector, u k ∈ R du the input vector and the matrices A, B of appropriate dimensions. Also, the physical constraints lead to the compact sets X ⊂ R d and U ⊂ R du , respectively.

Therefore, the MPC problem to be solved at each time step throughout a corridor Πi can be formulated as:

P( Πi , N p ,X f , xi , X ) : min u J (N p , xi , x k , u) (5.5a) s.t. x k+l+1|k = Ax k+l|k + Bu k+l|k , (5.5b
)

u k+l|k ∈ U , ∀ℓ = 1 : N p -1 (5.5c)
x k+l|k ∈ Πi (5.5d)

x k+Np|k ∈ X f (x i ), (5.5e)
where the constraint (5.5e) is based on the existence of an invariant set containing the origin X 0 f which can further be parameterized w.r.t the shift of the fixed point, leading to X f (x i ).

Having the formulation (5.5), the question to be addressed is how to select the parameters N p , X f , X in order to ensure the recursive feasibility for a given pair (x i , Π i ). The recursive feasibility is one of the fundamental challenges in MPC literature. Basically, it represents the property that for all initial feasible states and for all optimal sequences of control inputs the MPC optimization problem remains feasible [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]]. [START_REF] Löfberg | Oops! I cannot do it again: Testing for recursive feasibility in MPC[END_REF]] offers a broad overview on this topics, especially from the computational perspective.

The selection of the parameters in (5.5) can be viewed as an additional design step, which builds on the, so-called, backward reachable set (BRS).

The N-step (BRS) is the set of all states that can reach a final position xi in N stepsassociated to the system described by LTI dynamics (5.2): .6) Since we have "hard" constraints on state (the ones given by the corridor), we have to compute the N -step BRSs taking into account its limitation: Ri N ⊆ Πi . In fact, the computation of a N -step BRSs remains valid as long as the controller operated within the same corridor constraints. From this observation, the computation of N -step BRSs can be done iteratively as in Algorithm 12.

R i N (x i ) = A -N   xi ⊕ N -1 j=0 -A j BU   . ( 5

Algorithm 12

Computing N-step BRS for (5.2) taking into account the state constraints Πi Input: xi , A, B, U , Πi Output: Ri As its name suggests, the idea behind the "Relay MPC" strategy is to ensure the transitions from the current segment of corridor ( Πi ) to the next ( Πi+1 ). For that reason, we choose the terminal sets X f i as the intersections between two consecutive segments. Apply MPC strategy solving P( Πi , N p , X f (x i ), xi , X )

6:

until x k+1|k ∈ X f (x i )

7:

x 0 = x k+1|k 8:

Update the parameters of P 9: end for Proposition 5.2. If Path(x 0 , x f ) exists and the control law based on P(•) (5.5) is recursively feasible then the convergence x0 → x f is guaranteed.

Proof. The existence of Path(x 0 , x f ) leads to the construction of the corridor Π = Πi . The recursive feasibility and the selection of terminal sets as in ( 1) ensure the transition from Πi to Πi+1 . Recursively, for any point in X f (x 1 ) the MPC problem with horizon N i p is feasible with respect to final constraints X f (x 2 ). On the same argument any point in X f (x i-1 ) is (recursively) feasible for a terminal constraint in X f (x i ) and thus the terminal set at the end of the path X f (x Nc ) can be reached in finite time. Moreover, for the last segment we have X f (x Nc ) = RNc 0 which directly leads to x f = xNc being reachable, and, by consequence, it proves the asymptotic convergence of the scheme.

Using dynamics (5.1) we revisit the obstacle collection depicted in Figure 4.9 to which we apply Algorithm 13. As well, in Figure 5.7 the values of the acceleration and velocity are plotted along the simulation horizon.

As a side remark, for the same initial and final position an MIP based MPC strategy can be employed, considering the entire set of obstacles and imposing a large prediction horizon (e.g., N p = 40). Despite the large prediction horizon, the resulting trajectory does not converge to the final position, the agent remaining on the boundary of one of the obstacles due to lack of vision on a long-horizon and the benefits of a short-time increase of the cost (the pink trajectory in Figure 5.6).

Regarding the offline part, as stated above, the computational effort related to the reachability is substantial, e.g., for the trajectory in Figure 5.6 is around 3 minutes, but this can be improved by replacing the polyhedral representation of the sets with a zonotopic one [START_REF] Althoff | CORA 2015 Manual[END_REF]] or considering upper bounds. Intuitively, there is a connection between the length of a corridor segment (∆γ i = γ i (1)γ i (0) ) and the length of the minimal prediction horizon(N p min ). Therefore, we depict in Table 5.3 the values corresponding to (5.9). Remark 5. 1. In practice, whenever the length of a segment of the corridor is greater than a user-defined value, that segment can be split Table 5.3.: The values of ∆γ i and N p min for the corridors depicted in Figure 5. 6. fact, the dimension of the corridor. That is, the function ρ(•), as in (4.7), becomes ρ(•) = ρ(•) + ∆ρ. In Table 5.4 we delineate: N goal -the number of steps to attain a neighborhood of the final point, t goal -the total time to compute the trajectory and ℓ t -the trajectory length. As it can be seen in Table 5.4, the computing time t goal has similar values (slight modifications), while N goal and ℓ t increase with the values of ∆ρ. This behavior is counter-intuitive, but it can be explained by the fact that the decreasing in terms of steps (or distance) gained along of the corridor are wasted with the maneuvers associated to changing the segments of the corridor (activation of inputstate constraints). However, by moving all complex operations from on-line to off-line, we note that the computing time is kept reasonable, allowing a comparison with the classical sampled-based methods. Table 5.4.: The impact of corridor dimension (ρ(•)) on the computational features.

Adjustments/Adaptation for time-varying cluttered environments

Before the main developments on-time varying cluttered environments we will establish a series of restriction on these displacements which have to be counteracted by the control of each agent. Obviously, a unbounded degree of freedom in the motion of obstacles cannot be counteracted by constrained control, thus, the following apparently restrictive assumptions are borrowed from the various mechanical robots applications, so-called capsule robots [START_REF] Chernousko | On the optimal motion of a body with an internal mass in a resistive medium[END_REF]] or vibro-robots [START_REF] Vartholomeos | Dynamics, design and simulation of a novel microrobotic platform employing vibration microactuators[END_REF]. This class of mechanical robotic system is characterized by a one-dimensional motion, and their dynamics is generated by the periodic movement of an internal mass.

Remark 5.2. In the sequel, in order to concentrate the scope of the work we have examined the case in which the dynamical behavior of the environment consists only in the translation motion of the obstacles, neglecting any rotational motion behavior. Thus, the following assumptions are necessary.

Assumption 5.1. A movement of an obstacle P i ⊂ P refers to a displacement of its center on a given direction

δ i ∈ R d .
scratch of the entire corridor, and we compute only those segments of the corridor affected by the movement of the obstacle. Algorithm 14 sum up the entire procedure.

Algorithm 14 Continuous path generation within time-varying environment

Input: Γ(x 0 , x f ) (for the initial setting), K = {1, . . . , k max } Output: a path π : [0, 1] × K → X with π(θ, k) ∩ P(k) = ∅, ∀θ, k ∈ K 1: Find Path(x 0 , x f ) using Proposition 4.2 for Γ(x 0 , x f ) . Find ρ(θ, k) such that:

ρ(θ, k) ≤ min P i (k)∈P(k) d H (P i (k), γ(θ)), ∀θ ∈ [0, 1];
(5.12)

5:

Find a path π(θ, k) inside the corridor characterized by γ(θ), ρ(θ, k) as:

π(θ, k) = arg min ν C(ν) s.t ν(θ) ∈ γ(θ) ⊕ B 0,ρ(θ,k) , ∀θ ∈ [0, 1].
6: end for Remark 5. 3. The set K = {1, . . . , k max } is obtained from (5.10), i.e., k max is corresponding to the displacement of the obstacle center given by τ max * δ j taking into account a sampling time. This sampling time may differ from the one corresponding to the agent. In this case a detection mechanism is needed to identify the position of the moving obstacle and generating the path, accordingly. ♦

As mentioned before, each step of Algorithm 14 can be viewed as a further refinement along the trajectory generation procedure. While the first steps of the algorithm provide the graph structure and γ(θ) (a PWA (continuous) function) associated to the initial static workspace, Step 4 and Step 5 are taking into account the time-varying behavior, producing ρ(θ, k), a piece-wise width function that provides a measure of the permissible deviation from the nominal γ(θ).

Remark 5. 4. Similarly to the static case (Algorithm 11), Steps 4 and 5 can be improved by choosing a standard approach in motion planning field, by reformulating the problem in terms of the weights associated to a collection of functions. ♦

For the sake of illustration, we recall the obstacle collection shown in Figure 5.8 to which we apply Algorithm 14 and we depict in Figure 5.9 the resulted corridors and how they differ w.r.t. the one generated by the static environment. Note that we consider a sampling time for the moving obstacle such that K = {1, . . . , 10}.

As stated above, the moving obstacle has an influence only on a part of the corridor. More precisely, up to k = 5 only one segment of the corridor is altered, after this instant (k ≥ 6) we have two segments of the corridor which need to be updated.. The remaining segments are identical with the ones in Figure 5.9a. Since the corridor is constructed as in Section 4. 3.3 it is clear that the maximal number of segments affected by a moving obstacle at a given time instant is two. Therefore, Algorithm 14 can be improved in order to skip those segments of corridor whose width remain the same.

A potential drawback of the method summarized by Algorithm 14 is that it may spawn narrow segments of corridor as the one in Figure 5.9e. This may generate issues in i) the strategy employed for path/trajectory planning ii) the number of moving obstacles However, the current study remains preliminary and does not offer clear guidelines for the selection of the more appropriate method.

Extensions to multi-agent systems.

In most of the developments of the previous chapters, we have considered only one agent within the multi-obstacle environment. As stated elsewhere, there exist however applications which can be arguable better performed by a team of agents. By consequence, a second strong requirement become apparent, i.e., the avoidance of collisions among those agents. Thus, besides the (usually) time-invariant constraints induced by the obstacles avoidance, in the navigation problem formulation there exists a set of time-varying constraints generated by the anti-collision requirement.

As a side remark, in what follows, we assume that each agent has access to the so-called awareness map which contains all available information about the multi-obstacle environment. Moreover, we presume that there exists between any two agents a communication link, which allows a bidirectional transfer of information, e.g., vehicle-to-vehicle (V2V) communication [START_REF] Darbha | Benefits of V2V communication for autonomous and connected vehicles[END_REF]. Under these assumption, we are able to consider the following adaptation of the navigation strategy.

As stated in the introduction of Chapter 5, the method herein has two phases: planning and tracking. However, for sake clarity, we consider a particular breakdown of the navigation strategy into decision levels, as depicted in Figure 5.11.

Hence, the strategical layer consists of all procedures involved in the characterization of the environment according to the method presented in Chapter 4. In other words, this level builds a processed awareness map, in the sense that the procedure leads to a map containing all the possible corridors (similar with a roadmap) within the cluttered environment. Practically, starting from the obstacles distribution Algorithm 5 is employed, the graph Γ is constructed as in Section 4. 3.1 and, by contrast to the single-agent strategy (Algorithm 11), the corridors are constructed for the entire graph Γ using (4.31).

The next decision layer, the tactical one, covers the remaining steps of the planning phase. Firstly, for each agent the awareness map is adjusted by adding the initial and final positions to the graph Γ as in (4.30) and finding the shortest path through Γ between those two positions. Thus, to each agent a sequence of corridors is attached/allocated.

Finally, at the operational level, the effective control strategy is employed. Specifically for the method herein, for each agent a variant of the Relay MPC (Algorithm 13) is considered, as will be further detailed.

Considering the decision breakdown in Figure 5.11, there are two possible modes to intervene for satisfying the collision avoidance requirement. Clearly, we can eliminate the risk of collision from the planning phase (more precisely, tactical level) by selecting paths with disjoint sets of nodes (and, implicitly, of edges). Nevertheless, this may lead to trajectories too far (in terms of performance or consumed energy) from the optimal ones. In order to elude this shortcoming and, at the same time, to streamline the tracking phase, we can impose in the planning phase that the paths of the two agents share a minimal number of common nodes, as resumed by the next remark.

Conclusions

The chapter presented a constructive solution for the generation of collision-free trajectories between two points in an environment containing multiple obstacles in a d-dimensional space. This builds on the geometry of the obstacles and the convex lifting procedure describing a graph around the obstacles. This graph represents a key element in order to generate collision-free trajectories employing MPC controllers with recursive feasibility guarantees and convergence in between an initial and a final position.

Not in the least, the present chapter concerned the path planning in time-varying congested environments. Basically, we extended the results for time-varying multi-obstacle environments, considering, as well, re-configuration and adjustments of the graph. While the procedure is well established for static environments, a dynamical evolution of the context can support a certain degree of robustness by adjusting the width of the navigation corridors. Whenever the reconfiguration of the path is needed, the computation will be impacted accordingly. The second part of the present chapter covered the adaptation of the generic methodology providing an explicit resolution in the case of a a priori known dynamical behavior.

Conclusions

Conclusions

The present manuscript had as main objective to establish and to shed light on the safe navigation strategy in the presence of constraints. Concepts from control theory, optimization and computer science have been gathered and used to create valuable techniques that were then applied for navigation problems involving multi-obstacle environments. Thus, we used and built on a combination of optimization-based control and set-theoretic tools to obtain safe navigation strategies.

In the first part of the manuscript, an in-depth theoretical study of the concepts of MIP (Mixed-integer Programming) is presented. This study provided a detailed review of breakthrough research results and open issues in the field of multi-agent motion planning in a mixed-integer framework. As mentioned in the respective context, although the history of MIP starts almost 60 years ago, the interest of the control and robotic community on this topic is relatively recent and the research in this area is quite active. Clearly, there exists a substantial progress in all MIP-related topics. However, there are many points that can be further enhanced since the developing of new and performing methods to provide exact solution for MIP is exponentially growing. We identified the open and active problems regarding mainly the obstacle and collision avoidance topic, e.g., trade off between conservativeness and complexity in non-convex region modeling and representation.

Even though there exist several valuable improvements on classical MIP formulations, the complexity remains an arduous issue imposing restrictions, foremost among them, only reduced-dimension and complexity can be solved in a real-time manner. Additionally, another aspect that can be viewed as an open issue and, in some sense, is generated by the previous ones is that there are not so many works in the field of MIP-based motion planning validating their theoretical results. This issue is still happening even though there are many applications able to benefit from the theoretical results, e.g., search and rescue, environmental measurements tasks, area coverage and the like. Moreover, our aim was that this work to be employed to the benefit of both control and optimization research communities allowing to swiftly identify previous, timely and relevant research topics in the field.

We have to emphasize that the navigation within cluttered environment represents a broad and stimulating domain [START_REF] Latombe | Robot motion planning[END_REF]]. Thus, we concentrated our effort on a specific thematic related to obstacle avoidance, a complex and challenging issue. This problem is commonly stated in the literature in terms of a constrained optimization problem over a non-convex domain. Describing this non-convex domain has computational and structural implications, generally leading to a trade-off between computational efficiency and control performance. Preliminary results make use of hyperplane arrangements to characterize these regions Ionela [START_REF] Prodan | Mixed-integer representations in control design: Mathematical foundations and applications[END_REF]. The improvements may be accomplished by exploiting the underlying combinatorial structure of the MIP formulation. Nevertheless, whenever the problems are inherently non-convex and/or they involve alternative choices, mixed-integer representations provide an useful and powerful tool, but we have to proceed with caution evaluating the structural properties that may lead to compact formulations We opted to address this issue by using zonotopic over-approximations. With their particular representation, zonotopes balance between complexity of representation and numerical difficulty. Firstly we provide the necessary tools to overestimate the shapes of interest with predefined complexity and we focus on the benefits of the zonotopes regarding the decrease of the computational burden of general collision avoidance problem. Afterwards, the zonotopic approximations have been exploited from the control perspective, treating the eventual issues generated by their use in the control problem like the lost or the narrowing of the passageways among the obstacles. These last topics were not addressed in the open literature to the best of our knowledge. Therefore, we provide a set of conditions able to safeguard the initial domain topology, maintaining the passageways and, simultaneously, the advantages given by the zonotopic representation of the multi-obstacle environment. As well, we extend the results and we treat the problem of collision avoidance in a time-varying multi-obstacle environment. Although, the previous formalism is adequate for a static multi-obstacle environment, it may be impractical in a dynamical context. Nevertheless, to overcome this issue we make use of a technique which introduces an additional analysis step, by evaluating the problem in a lifted space (using a suitable parametrization).

Not in the least, we developed a constructive solution for the generation of a path between two points in an environment obstructed by multiple obstacles in a d-dimensional space.

The global information on the geometry of the obstacles is considered as an entry point for a convex optimization procedure which leads to a convex lifting allowing the partitioning of the cluttered environment. This partitioning is a key element for describing a graph around the obstacles and ultimately for the generation of corridors which avoid obstacles. From the computational point of view, the validity of the construction relies on the feasibility of the convex lifting. It was shown that feasibility can be improved by a reformulation of the obstacles in terms of a finite number of convex subsets. Furthermore, this principle allows the generalization of the construction for path planning in the presence of multiple non-convex obstacles. In the latter case, a disconnected graph within the generation of the corridors will certify the infeasibility of the path planning problem for at least a pair of points in the cluttered environment. Finally, we have addressed the navigation problem in a multi-obstacle environment, extending the convex lifting approach and providing feasibility guarantees. Moreover, we consider an adaptation of the generic MPC (Model Predictive Control) for corridor-constrained trajectory tracking problem, aiming to guarantee the feasibility and convergence. Also, we have extended the results for timevarying multi-obstacle environments, considering re-configurations and adjustments of the graph, which generates the geometric path.

Future directions

First of all, we consider that the most relevant topic to be prospectively treated is the development of distributed control strategies for task partitioning and guaranteed constraint validation for the multi-agent systems (MAS). While in the present manuscript the focus was on modeling and on numerical performance issues, future work should to integrate all of these into the actual multi-level control block. Indeed, eventually the multi-agent system has to have a task to solve. This will require that the task partitioning phase has been optimally treated, all by handling communication constraints and the explicit introduction of the constraint and cost constructions from before. The approaches should exploit the existing consensus strategies such that stability and performance are guaranteed a priori at the level of the MAS.

Secondly, a further development is the deployment of the methodological and theoretical results.The ideal test of the theoretical results should include an experimental platform (e.g., aerial units like nano-quadcopters with complex scenarios in terms of obstacles, in order to cover the possible safe evolution of civilian UAV in a crowded urban area). These challenging benchmarks preserve a certain complexity of the low level control and can operate in an indoor environment or in the open field. There are many applications, both civilian and military, which can benefit from their use, as for example search and rescue, communications, environmental measurements tasks, area coverage and the like.

As well, a topic that need to be considered for improvement is mixed-integer representation for a time-varying cluttered environment. As stated, each agent (i.e., robot or any unmanned vehicle) has a so-called awareness map which contains all available information about the external medium. Obstacles and other agents (denoted by a priori given safety region) need to be observed; their shape to be estimated and, ultimately, introduced into the awareness map. This will be done through parametrized mixed-integer descriptions. E.g., if the shape is fixed and the future trajectory of an agent / obstacle is known (or predicted), then the map will be updated instead of being re-computed from scratch at each environmental modification. Furthermore, if the agent is part of the formation, the networked estimation and reconstruction of the scene propose challenges from the distributed decision making perspective. Finally, we consider that set-theoretic notions like positive and controlled invariance can be further improved to minimize online computations. These tools and notions are used to characterize safety regions around agents, to validate stability and performance criteria (both in the sense of tracking error and estimation error).

A. Alternatives approaches to MIP-based motion planning

This section briefly presents the existing alternatives to MIP which are extensively used in motion planning problems. We delineate in Table A.1 the state-of-the-art references for each of these alternatives.

A.1. Alternative formulations

Alternative formulations for MIP range from graph-based approaches, potential field formulations to the corresponding optimization problems relaxations.

Graph-based approaches

In contrast to the MIP approaches where the discrete decisions are encoded in a mathematical formalism and are solved as such, the graph-based approaches reduce these discrete decision to the search of the shortest path between nodes in a graph. Although these techniques can be applied for any MIP-based motion planning problems (see Section 2.3), in this subsection we concentrate on collision-free path/trajectory planning.

Remark A. The graph-based approaches in path planning are classified by graph construction. A first category builds the graph based on the cell decomposition methods, considering an explicit representation of the (multi-obstacle) environment. Thus, the similarity with the MIP approach is undeniable, polytopic sets being a popular decomposition primitive.

For instance, polytopic sets were already used in [Lozano-Pérez and Wesley 1979]. More precisely, the graph nodes are defined by the vertices of the polytopic obstacles, and by the initial and the final point of the path. This graph is named the visibility graph because the link between the nodes is given by a straight line which does not intersect any obstacles.

In other words two linked nodes can "see" each other. Thus, the collision-free path is the shortest path through this graph between the initial and final position of the agent. This approach performs well from a theoretical perspective if the considered agent is a point. For the general case when the dimensions are not negligible, the above method is extended by considering an artificial "growing" of the obstacles (commonly, via Minkowski sum). The main drawback of the above presented method is the computational burden, specifically in a complex environment (large number of obstacles with complicated forms) where the number of vertices is exhaustive.

The second category tries to mitigate this issue avoiding an explicit representation of the environment. This type of approach is also referred in the literature as sampled-based algorithms (see [START_REF] Lavalle | Planning algorithms[END_REF]) for a complete literature review on this topic). Hence, the explicit representation is eliminated and the focus is put on a collision checking module. This module checks the feasibility of the trajectories connecting random sampled points from the obstacle-free space. The graph is constructed such as the link between two nodes represents a feasible path in the (multi-obstacle) environment.

Remark A.2. An interesting notion for sampled-based algorithms is the probabilistic completness [START_REF] Barraquand | A random sampling scheme for path planning[END_REF][START_REF] Ladd | Measure theoretic analysis of probabilistic path planning[END_REF]. That is, the probability that the algorithm returns a feasible solution tends to one if the number of sample points is sufficiently large (→ ∞). This is empirically demonstrated in [Hsu, Latombe, and Kurniawati 2007]. ♦

In the literature, there are two important sampled-based algorithms: PRM (Probabilistic RoadMaps) [Hsu, Latombe, and Kurniawati 2007] and RRT (Rapidly-exploring Random Tree) [Weiss et al. 2017]. The difference between them is given by the method of constructing the graph. The former (PRM) is a multiple-query method in the sense that after the construction of the roadmap (a rich set of feasible paths) it answers queries by computing an optimal path through the graph. Henceforth, the PRM is an useful method if an awareness map of the environment is available [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF]. There are a plethora of variants for the PRM, each representing a valuable improvement [START_REF] Ladd | Measure theoretic analysis of probabilistic path planning[END_REF][START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF].

For the case when the environment is not a priori known, the RRT method is more suitable. In this approach the construction of the graph is incremental, the algorithm stops when a large enough set of collision-free paths is attained. Thus, a collision-free sample is added as a node in the graph and is linked with the surrounding nodes. The obtained graph is actually a tree. As PRM, RRT has a variety of versions. Some take into account the equations of motion and generate attainable paths [LaValle and Kuffner 2001], others generate only geometric paths which became reference trajectories for a lower level controller. Moreover, some versions are tailored for complex and/or unstable dynamics [START_REF] Leonard | A perception-driven autonomous urban vehicle[END_REF]] or for uncertain dynamics [Weiss et al. 2017].

Besides the classical graph-based approach, there exists in the literature a multitude of methods which combine the standard graph algorithms with advanced control strategies. For instance, [START_REF] Berntorp | Positive invariant sets for safe integrated vehicle motion planning and control[END_REF]] presents a method for real-time integrated motion planning which uses feedback control, positive invariant sets, and equilibrium trajectories of the closed-loop system. In order to generate the collision-free trajectories the method employs, in an offline manner, a graph search over reference paths, each being associated with a constraint admissible positive invariant set. Next, they use pre-designed unconstrained linear quadratic controllers to track the reference paths.

As well, [START_REF] Altché | Partitioning of the free spacetime for on-road navigation of autonomous ground vehicles[END_REF] addresses the problem of trajectory planning using an approach which consists in partitioning the feasible ("collision-free") region, while allowing to decompose the NP-hard problem as a path-finding problem in a well-designed graph followed by a (simple) optimization phase ("in MPC fashion") for a quadratic convex cost function. Also, [START_REF] Franzè | The obstacle avoidance motion planning problem for autonomous vehicles: A low-demanding receding horizon control scheme[END_REF] deals with the problem of obstacle avoidance in an unknown environment (considered agents are UGVs -autonomous ground vehicles). The proposed approach consists of two parts: an offline part which computes the ellipsoidal approximation of the one-step controllable sets for all possible scenarios (these approximations guarantee the existence of a feasible path through multi-obstacles environment) and a online part which involves the development of a MPC-based strategy in order to keep the agent in that sequence of ellipsoidal sets.

Potential field formulation

While MIP-based methods take explicitly into account the constraints and lead to a constrained optimization problem, the potential field-based formulations [Y.-b. Chen et al. 2016] relax the constraints by adding penalty terms in the cost. Essentially, the potential field approach relies on construction of a scalar function (so-called, the potential). This function takes high values when the agent stays within the interdicted zones. In the collision-free workspace the function is decreasing towards the goal configuration (i.e., the potential associated to the destination point is minimal). Thus, the agent may attain the final point moving in the direction of the negative gradient of the potential. [START_REF] Rimon | Exact robot navigation using artificial potential functions[END_REF] provides an historical (and more detailed) review on the potential field formulation and how this approach is involved in motion planning. An interesting characteristics is that the potential field formulation is frequently used in decentralized or distributed control strategies [START_REF] Filotheou | Decentralized control of uncertain multi-agent systems with connectivity maintenance and collision avoidance[END_REF].

A.2. Optimization problem relaxation

As was stated above, one of the most significant capabilities of MIP is to handle non-convex constraints in non-convex optimal control problems. A natural way to tackle this type of problems is by extending the methods and techniques used for the convex optimal control. Often, the MI formulations are solved heuristically (e.g., [START_REF] Quaritsch | Networked UAVs as aerial sensor network for disaster management applications[END_REF]] applies genetic algorithms) or relaxed by iterative solving (the optimization problem is broken into "reasonable" sub-problems which are solved iteratively). As an example, [START_REF] Xu | Energy-efficient UAV communication with multiple GTs based on trajectory optimization[END_REF]] employs binary variables to model the link between an UAV and a ground sensor but relaxes the formulation through time allocation tactics and channel communication pre-scheduling. Thus, in the literature, a considerable effort was put in finding a technique which allows the translation/relaxation from non-convex to convex formulations without any major gap. This is referred as convex relaxation or convexification of the non-convex optimization problem. In the literature there are various works providing methods for convexification under several labels: e.g., convex relaxation, succesive convexification or time-varying constraints.

As its name suggests, the basic idea of successive convexification is to solve a non-convex optimal control problem via a sequence of convex sub-problems. The non-convexity comes either from having non-linear dynamics [Mao, Szmuk, and Acikmese 2016] and/or from non-convex state (and/or control) constraints [Mao, Daniel Dueri, et al. 2017]. In both cases the same technique is applied: linearization, commonly, using a first order Taylor approximation (in a successive manner). Therefore, a preliminary condition on the functions which generate the non-convexity is necessary: they have to be differentiable. In an iterative manner, the linearization is done about the solution obtained at the previous step.

Although the linearization procedure leads to convex formulations, it also introduces two new issues, namely artificial infeasibility1 and approximation error. These two drawbacks are addressed in the literature, a variety of algorithms was developed and, recently, a convergence analysis was elaborated [X. [START_REF] Liu | Solving Nonconvex Optimal Control Problems by Convex Optimization[END_REF].

For instance, [D. [START_REF] Dueri | Trajectory optimization with inter-sample obstacle avoidance via successive convexification[END_REF]] deals with the problem of the trajectory optimization for autonomous vehicles in an environment containing cylindrical and ellipsoidal obstacles.

The approach employs the Successive Convexification technique which is used to solve the non-convex optimal control problem via a convergent sequence of convex optimization problems. It considers a discrete-time, finite horizon constrained optimization problem with a number of non-convex state constraints. Several assumptions are necessary in order to employ the technique at hand. The first assumption can be readily satisfied and involves that the obstacles boundaries should not be in contact with the ones of any other constraint. A second problematical assumption builds on a dynamic with a finite number of stationary points. Having the general formulation for the non-convex problem, the successive convexification technique is applied in order to obtain the sequence of convex sub-problems, each being linearized in an iterative procedure. This linearization results in a convex problem but, at the same time, introduces two drawbacks: approximation error and artificial infeasibility. To mitigate these two issues of the convexification, the authors introduce trust regions and penalty functions, respectively. The drawbacks of the convexification procedure and the ways to alleviate them are treated in-depth in the references such as [START_REF] Harris | Lossless convexification of non-convex optimal control problems for state constrained linear systems[END_REF]].

An interesting alternative for handling the non-convexity of the collision avoidance problem is based on time-varying constraints. The idea is mentioned in [Frasch et al. 2013] and is used in, e.g., [Janeček, Klaučo, Kalúz, et al. 2017] and [START_REF] Yu | Energy-Efficient Trajectory Planning for a Mobile Agent by Using a Two-Stage Decomposition Approach[END_REF]]. In few words, the non-convex domain is decomposed in a sequence of convex regions and switching instants are introduced. At each moment the agent should stay in one of the convex regions. For instance, in [Janeček, Klaučo, Kalúz, et al. 2017] the switching instants are steps of the prediction horizon in the MPC controller. Moreover, this approach is coupled with a heuristic black box which establishes the get-around direction and the sequence of convex sub-domains. This method can be seen as a particularization of the successive convexification. The difference comes from the way of obtaining the sequence of convex sub-problems/sub-domains. Time-varying constraints approaches consider an a priori known number of sub-domains, since the successive convexification is an iterative method, the sequence is growing until a feasible solution is obtained.

In [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF] a different approach for convex relaxation is presented. For instance, therein the authors treat the aircraft conflict problem 2 . Basically, the provided formulation is based on complex numbers representation and it results in a tight convex relaxation for the inherently non-convex optimization problem. It is worth mentioning that the above reference includes a comprehensive literature review on the formulation of the air conflict problem as an optimal control problem using mixed-integer techniques (as MILP or as MINLP) 3 . Coming to the application, the aircraft separation condition is stated: the relative position of two aircraft should be greater than a certain threshold. As it was expected, this condition leads to a non-convex feasible domain which is modeled using a binary variable (actually, the feasible region is composed of two convex ones and the solution should be in one of these two convex subsets). The control actions (speed variation rate and heading deviation angle) admit a natural representation in the complex number form. Even so, the non-convexity is not eliminated (the disjunction constraint is maintained) but the considered formulation is useful towards the convex relaxation approach. This approach is extensively treated in [START_REF] Coffrin | The QC relaxation: A theoretical and computational study on optimal power flow[END_REF],

where the non-convex constraint is tackled by deriving the corresponding convex hull, the problem being transformed in a MIQCP (mixed-integer quadratically constrained program). A further relaxation is possible by entirely omitting the non-convexity. An algorithm which included the relaxation is presented and tested (with excellent results) on two classical benchmark problems for conflict resolution.

As well, [X. [START_REF] Huang | Speed trajectory planning at signalized intersections using sequential convex optimization[END_REF] solve the MI optimization problem using sequential convex optimization method (i.e., search the local optimum by forming convex sub-problems), avoiding the eventual curse of dimensionality. Likewise, [Papen et al. 2017] solve a MILP by using a relaxation of the constraints and fine-tuning the complexity in order to limit the computation time w.r.t. a classical MILP resolution.

B. Control architectures for MIP-based motion planning

There are three well-established classes of control architectures, and they have been extensively studied in various application domains: centralized, distributed and decentralized.

The last two methods require the local controllers to optimize over only their local inputs having similar computational burden. The difference between these two is given by the impact of communication, decentralized control requires no communication among the agents.

To identify/discuss the control architectures weaknesses and strengths with respect to the motion planning field is not the scope of this manuscript, thus, we only focus on how they are employed along with MIP. In Table B.1, we delineate a classification of MIP-based motion planning references w.r.t. the control architectures and number of agents involved in the formulations.

It is worth mentioning that the control strategies employed in other architectures than the centralized way are optimization-based, with predilection MPC. Thus, in this section an emphasis is put on the specific MPC implementations. Nevertheless, we have not neglected the references where MIP is used specifically in distributed/decentralized non-MPC strategies due to its capability to formulate task allocation problems. These are extensively treated in Section 2.3.1.

The control architectures involving MIP naturally evolve from centralized to decentralized and distributed strategies.

B.1. Centralized

Due to its theoretical simplicity, the centralized approach is the most used way of controlling a multi-agent system. In this architecture, the multi-agent system is treated as a whole, equating an extended single-agent system. The physical restrictions (e.g. communication limitations) are completely ignored, each agent having complete knowledge on the behavior/actions of the others, all information being available in the single global controller [START_REF] Rawlings | Model predictive control: Theory and design[END_REF]]. However, this methodology is limited, not only because of the undeniable physical constraints, but also due to numerical difficulties, induced by the substantial complexity of the extended system.

As stated in Section 2.3, in the literature, MIP may appear at different levels of control.

For instance, in [M. [START_REF] Chen | Multi-vehicle collision avoidance via hamilton-jacobi reachability and mixed integer programming[END_REF] the problem of collision avoidance is treated for a multi-agent system (minimum N = 3 agents). They consider N agents, each having a similar dynamic and N -1, so-called, danger zones: Z ij . Every agent makes use of two controllers: a "liveness controller" helps to complete the agent own objective (e.g. reach a target) and a "safety controller" has to keep away the agent from the compared to centralized architecture: e.g., decrease of complexity and scalability. However, a loss of performance and of global stability may become apparent, being hard to be ensured compared to the centralized approach. Decoupling into independent nodes, design of robust control strategies, consensus seeking, all of these try to solve the problem, but still, with limited success [START_REF] Cao | An overview of recent progress in the study of distributed multi-agent coordination[END_REF]].

Within the MIP framework, most of the distributed control approaches using MPC strategy have been proposed. The features of MPC allow to handle explicitly the interactions between the different subsystems/agents. For example, in [START_REF] Schouwenaars | Safe trajectory planning of autonomous vehicles[END_REF]], a distributed MPC strategy is employed for navigation of a fleet of vehicles through a partially unknown cluttered environment.

Due to the inherent problem of complexity, and implicitly the lack of scalability, MIP was not a popular method for distributed architectures. However, the control community has given a particular attention to the distributed MIP resolution algorithms. For instance, [START_REF] Testa | A finite-time cutting plane algorithm for distributed mixed integer linear programming[END_REF] propose an algorithm for resolution of a MILP where the constraints are distributed among the agents. Likewise, [START_REF] Vujanic | A decomposition method for large scale MILPs, with performance guarantees and a power system application[END_REF]] provides a decomposition method particularly useful for large scale MILPs, based on Lagrangian duality. Moreover, there are some works, which formulate the problem using MIP techniques, and for its resolution an heuristic method is employed, e.g., [START_REF] Van Parys | Online distributed motion planning for multi-vehicle systems[END_REF] use Alternating Direction Method of Multipliers (ADMM).

B.3. Decentralized

During the last decades, decentralized control of multi-agent systems has gained a significant amount of attention due to the great variety of its applications, including multi-robot systems, transportation, multi-point surveillance and biological systems. An important topic of research is multi-agent navigation in both the robotics and the control communities, due to the need for autonomous control of multiple robotic agents in the same workspace. Important applications of multi-agent navigation arise also in the fields of air-traffic management and in autonomous driving by guaranteeing collision avoidance with other vehicles and obstacles.

As seen in the previous subsection, a main motivation for the modification of the control architecture is the computational burden of the centralized problem. A first alteration was to decompose ("distribute") the problem and solve the resulting sub-problems in expectation of achieving a consensus. Proceeding further and eliminating the aim of reaching consensus leads to a decentralized strategy ("everybody for itself"). Under this strategy, each agent owns a controller, which acts without taking into account any information about the behaviour of the other agents. Moreover, the exchange of information is limited (in most of the cases, reduced at a minimum).

A basic illustration for the decentralized approach is the merging junctions1 problem [START_REF] Bali | Merging vehicles at junctions using mixed-integer model predictive control[END_REF]. In a centralized (or in a well constructed distributed) approach, whenever two (or more) agents arrive at an intersection, there is a prioritization w.r.t a well defined criterion. For the decentralized approach, with limited communication and with no guarantee of optimality, in the most probable scenario the system reaches a deadlock. In other words, the decentralized approach cannot provide, usually, the theoretical guarantees characterizing the centralized and distributed approaches. However, in practice, the decentralized methods may lead to efficient strategies, avoiding the computational demand of the others methods; but an understanding of the sensitive aspects is absolutely necessary.

For instance, recalling the previously-mentioned (Subsection B.1) example from [X. Wang et al. 2015] -the problem of collision avoidance in a multi-robot system. Using a decentralized control strategy, they assume that each agent chooses its trajectory independently (from an a priori known set), without informing the rest of the systems. In this case the collision avoidance strategy should consider all possible trajectories combinations and has to return the initial time delays such that the movement to be finished in the shortest time without collisions. The resulting optimization problem is MIP, as in the centralized case. Note: (✓) treated, (×) not treated

C. Control applications using MIP

This section covers some critical details for both the computer simulations and the hardware implementation of the solutions and MIP approaches presented throughout this manuscript. It is worth mentioning that there does not exist some clear and uncontested guidelines capable of generating the most efficient MIP formulation for a given problem. The performance1 of a formulation is customarily strongly dependent on the specific software tools or hardware platform. In the sequel, these aspects are briefly documented, the emphasis being placed on differences between the existing alternative implementations and their influence over the practical performance in motion planning.

C.1. Software for MIP

As was stated above, MIP is a powerful tool for planning and control problems due to its modeling capability and, additionally, the availability of specialized solvers. In the last decade, a consistent effort was put on developing MIP-specialized solvers in order to mitigate the numerical issues generated by the presence of integer/binary variables.

The optimization modeling languages are toolboxes which convert the mathematical formulation in a solvable form for the solvers. [Janeček, Klaučo, Kalúz, et al. 2017;[START_REF] Bemporad | MIQP. M: A Matlab function for solving mixed integer quadratic programs[END_REF][START_REF] Culligan | Online trajectory planning for UAVs using mixed integer linear programming[END_REF] Python [START_REF] Welder | Spatio-temporal optimization of a future energy system for power-to-hydrogen applications in Germany[END_REF] Julia [START_REF] Lubin | Polyhedral approximation in mixed-integer convex optimization[END_REF] Modelling Language Yalmip [Y. [START_REF] Zhang | Modelling and traffic signal control of a heterogeneous traffic network with signalized and non-signalized intersections[END_REF]Mukai, Natori, and Fujita 2017] GAMS [START_REF] Lee | New algorithms for nonlinear generalized disjunctive programming[END_REF] CVX/CVXPY [D. [START_REF] Dueri | Trajectory optimization with inter-sample obstacle avoidance via successive convexification[END_REF] Pyomo [START_REF] Legg | A stochastic programming approach for gas detector placement using CFD-based dispersion simulations[END_REF] JuMP [START_REF] Welder | Spatio-temporal optimization of a future energy system for power-to-hydrogen applications in Germany[END_REF]] Solvers CPLEX [Janeček, Klaučo, Kalúz, et al. 2017;[START_REF] Stoican | Hyperplane arrangements in mixed-integer programming techniques. Collision avoidance application with zonotopic sets[END_REF] GUROBI [Janeček, Klaučo, Kalúz, et al. 2017] SCIP [START_REF] Berthold | Solving mixed integer linear and nonlinear problems using the SCIP Optimization Suite[END_REF] Prior to proceeding further, we should emphasize the difference between optimization modeling languages and solvers. The former designates a toolbox/package/library which interacts with the latter: the mathematical model of the constrained optimization problem is put into an internal form which is then solved, and whose subsequent result is retrieved and displayed. Note that most modern tools do pre-processing steps (which may reduce the number of binary variables) or may even, automatically, put the problem in a MI form (e.g., in YALMIP when specifying a complementarity condition or a bilevel program).

AMPL [A.
Diverse programming languages and online resources facilitate the specification and resolution of MIP problems. In addition to classical Matlab, recently the attention of the optimization-based control community is moving towards other advanced programming languages as Python or Julia, which have become more accessible to the broader scientific and engineering community. In Table C.2 we summarize these programming languages, the modeling languages and the solvers which may be joined in order to solve MIP.

Without being exhaustive, there are some popular optimization modeling tools: YALMIP [START_REF] Lofberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF]], MPT [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF], AMPL [START_REF] Fourer | AMPL: A mathematical programming language[END_REF], CVXPY [START_REF] Diamond | CVXPY: A Python-Embedded Modeling Language for Convex Optimization[END_REF], PYOMO [START_REF] Hart | Pyomo-optimization modeling in python[END_REF] or JuMP [START_REF] Dunning | JuMP: A Modeling Language for Mathematical Optimization[END_REF]. All2 of them are open-source modeling languages which allows users to express a wide range of optimization problems (not exclusively, MIP) in a high-level (almost algebraic or pseudo-code) syntax. As depicted in Table C.2, each modeling language is developed taking into account the specific features of a programming language. Note that an hierarchy among these modeling tools is strongly dependent on the experience of the user/researcher, the preference for one of them has a negligible impact regarding the resolution performances (e.g., computational burden). These performances are influenced rather by the choice of the solver w.r.t. the formulation.

Regarding the solvers, there exists a broad variety of options, we mention here that the most used in the field of MIP-based motion planning are [START_REF] Cplex | V12. 1: User's Manual for CPLEX[END_REF]], [START_REF] Gurobi | Robotic swarm control from spatio-temporal specifications[END_REF]] or [START_REF] Mosek | The MOSEK optimization toolbox for MATLAB manual[END_REF]]. For instance, [X. [START_REF] Huang | Speed trajectory planning at signalized intersections using sequential convex optimization[END_REF]Y. Zhang, Su, C. Sun, et al. 2017] or [Mukai, Natori, and Fujita 2017], use GUROBI, while [A. [START_REF] Richards | Aircraft trajectory planning with collision avoidance using mixed integer linear programming[END_REF]Schouwenaars, A. Richards, et al. 2001] or [M. [START_REF] Earl | [END_REF] utilize CPLEX.

Remark C.1. The majority of the solvers (mentioned in Table C.2) is able to deal with quadratic objectives and/or constraints, elements which are influential in many control strategies and/or applications, see, e.g., (2.18). ♦

While a review of how the resolution techniques employed by the solvers is well beyond the scope of this paper, we mention only the core techniques: branch-and-bound and cutting-plane algorithms. There exists a multitude of variations, each with its strengths and weaknesses. For instance, branch-and-cut method [START_REF] Stubbs | A branch-and-cut method for 0-1 mixed convex programming[END_REF] gathers the advantages of branch-and-bound and cutting-plane algorithms, decreasing the number of nodes to be explored in the search tree by iteratively introducing constraints to cut the feasible region.

In general, the solvers may be classified using different criteria, e.g. convex/non-convex, heuristic/deterministic. There are in the literature more detailed surveys, e.g. [START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF]], treating this topic, but what is relevant from the current paper's perspective is the following aspect. Some of the currently available and reliable solvers may employ heuristics in order to accelerate the standard algorithms. This is a necessary requirement especially for complex (large) problem formulations and real-time resolution. As a word of caution, the performances (particularly, computational times) may vary considerably from a solver to another due to the use of heuristics. Hence, the concept of "the best MIP solver" is pointless, in our opinion. On a more positive note, we have observed that for any given problem there can be found a solver, capable of handling it. Rather, the user should test with his/her own solver, observe the behavior and choose accordingly.

Besides these powerful commercial solvers, there exists a variety of non-commercial/opensource solvers able to provide reasonable performances, in some cases better than commercial ones. An important characteristic of this kind of solvers is that their adjustment and adaptation to the challenges encountered in the real-world applications and to the real-time control requirements can be done in a faster and more straightforward manner than commercial ones which, in addition, have to take into account commercial considerations, balance between free/payed functionalities, etc.. For instance, SCIP (https://scip.zib.de) [START_REF] Achterberg | SCIP: solving constraint integer programs[END_REF]] was at the beginning a MILP solver which implemented the branch-and-bound algorithm with various heuristics, while the later versions were able to solve MINLPs with quadratic objective, non-convex MINLPs or MISDPs (mixed integer semidefinite programs).

Moreover, a part of the control community has concentrated its attention on techniques to adapt the standard MIP resolution algorithms. For instance, [START_REF] Bemporad | An efficient branch and bound algorithm for state estimation and control of hybrid systems[END_REF] propose an efficient branch-and-bound algorithm, enhancing the tree exploring strategy. The application therein involves the control and state-estimation for a MLD-system [Alberto Bemporad andMorari 1999]. Similarly, [Alberto Bemporad 2015] provides an algorithm which combines the classical branch-andbound with nonnegative-least-squares (NNLS) methods, in order to solve MIQP problems generated by the hybrid MPC applications. The idea is further developed in [START_REF] Naik | Embedded mixed-integer quadratic optimization using accelerated dual gradient projection[END_REF] where the NNLS is replaced by accelerated-dual-gradient projection algorithm.

There are also some other works which exploit the particularities of the problem structure within branch-and-bound algorithm. For example, [START_REF] Feng | Branch-and-lift algorithm for obstacle avoidance control[END_REF]] propose a variation of the branch-and-bound algorithm, branch-and-lift algorithm, which has better performances in resolution of a classical obstacle avoidance problem. As well, [START_REF] Hespanhol | A Structure Exploiting Branch-and-Bound Algorithm for Mixed-Integer Model Predictive Control[END_REF] provide an iterative variation branch-and-bound algorithm which exploits the block-sparse optimal control structure of the problem and also the information at previous time steps.

As a side remark, besides the standard MIP resolution algorithms, there exist in the literature some heuristic techniques. Without being exhaustive, we mention only two3 of these:

i) ADMM (alternating direction method for multipliers)-based methods [START_REF] Kanno | Alternating direction method of multipliers as a simple effective heuristic for mixed-integer nonlinear optimization[END_REF][START_REF] Takapoui | A simple effective heuristic for embedded mixedinteger quadratic programming[END_REF]]. Although an algorithm for solving convex optimization problems, ADMM turns out to be an effective method to approximately solve some nonconvex problems as well. The idea behind the heuristic method is to use multiple restarts of ADMM with random initial points, in most of cases this provides an acceptable solution with small computational cost. This technique is frequently used in optimal (power) flow problems, e.g. traffic signal control [START_REF] Timotheou | Distributed traffic signal control using the cell transmission model via the alter-nating direction method of multipliers[END_REF] or network [P. [START_REF] Scott | Dynamic optimal power flow in microgrids using the alternating direction method of multipliers[END_REF][START_REF] Abboud | Distributed caching in 5G networks: An alternating direction method of multipliers approach[END_REF]]. However, [Van Parys and Pipeleers 2016] employ ADMM in motion planning for a multi-vehicle systems. As well, for this technique there exist a multitude of variations and adaptations. For instance, [START_REF] Stellato | Embedded mixed-integer quadratic optimization using the OSQP solver[END_REF]] develop a novel branch-and-bound algorithm, tailored for a ADMM-based solver.

ii) FP (Feasibility Pump) [START_REF] Fischetti | The feasibility pump[END_REF]. A heuristic method for finding a feasible solution of a given MIP, FP aims to minimize the difference between the solution of the LP-relaxation and the one of the original MIP. For instance, [START_REF] Miertoiu | Feasibility Pump Algorithm for Sparse Representation under Laplacian Noise[END_REF] use and adapt the algorithm for sparse representation.

Table C. 3.: Alternative and heuristic techniques for MIP resolution Alternative References Branch-and-bound variation [START_REF] Bemporad | An efficient branch and bound algorithm for state estimation and control of hybrid systems[END_REF][START_REF] Bemporad | MIQP. M: A Matlab function for solving mixed integer quadratic programs[END_REF]Alberto Bemporad 2015;[START_REF] Naik | Embedded mixed-integer quadratic optimization using accelerated dual gradient projection[END_REF][START_REF] Feng | Branch-and-lift algorithm for obstacle avoidance control[END_REF][START_REF] Hespanhol | A Structure Exploiting Branch-and-Bound Algorithm for Mixed-Integer Model Predictive Control[END_REF] ADMM [START_REF] Kanno | Alternating direction method of multipliers as a simple effective heuristic for mixed-integer nonlinear optimization[END_REF][START_REF] Takapoui | A simple effective heuristic for embedded mixedinteger quadratic programming[END_REF][START_REF] Timotheou | Distributed traffic signal control using the cell transmission model via the alter-nating direction method of multipliers[END_REF]P. Scott and Thiébaux 2014;[START_REF] Abboud | Distributed caching in 5G networks: An alternating direction method of multipliers approach[END_REF]] FP [START_REF] Fischetti | The feasibility pump[END_REF][START_REF] Fischetti | Feasibility pump 2.0[END_REF][START_REF] Miertoiu | Feasibility Pump Algorithm for Sparse Representation under Laplacian Noise[END_REF] Remark C.2. Apart of solvers and modeling languages, there exist in the literature works as [Janeček, Klaučo, Kalúz, et al. 2017] which provide a toolbox for MPC-based control for obstacle avoidance problem. The toolbox is developed in an object-oriented manner, allowing to readily set up the associated MPC problem and, afterwards, providing an efficient formulation of the underlying optimization problem for any inexperienced user.

In other words, the OPTIPLAN toolbox encapsulates all aspects related to formulation and resolution of the MIP problem. ♦

C.2. Hardware platforms

Before detailing, it is worth mentioning that some of the existing works takes into account the particularities of the hardware platforms developing specific methods, whereas most of the remaining presents generic methods able (or not) to adapt to the constructive constraints.

There exists a broad variety of robotic platforms which are used in academic and/or commercial applications: aerial, surface or underwater vehicles. Frequently, these robots are involved in activities which are unsafe or troublesome for humans. Even though different levels of autonomy are possible, we can easily remark that the control community has considered with predilection the unmanned vehicles. The main justification comes from the elimination/mitigation of human risk. This aspect has beneficial consequences on cost effectiveness and, in most of the cases, on accuracy of the operations.

In Table C.4, we depicted the classes of unmanned vehicles which have been used as hardware platforms for MIP-based motion planning problems: UAV (unmanned aerial vehicle), USV (unmanned surface vehicle) or UGV (unmanned ground vehicle), UUV (unmanned underwater vehicle). The specific characteristics of each class (or, sometimes, sub-class) of vehicles lead to various challenges in the design of the motion planning strategy. For instance, the USVs/UGVs are moving in a bidimensional workspace, whereas the UAVs and UUV in a 3D workspace thus leading to a higher complexity of the control problem. Another classification is given by the ability to stop and/or go backwards, for example, fixed-wing UAV need to maintain a minimum velocity (to avoid stalling), but the quadcopters/helicopters (rotary wing UAVs), having more degrees of freedom, can maintain an arbitrary velocity (up to being stationary in midair, i.e. hovering).

For instance, in [START_REF] Culligan | Online trajectory planning for UAVs using mixed integer linear programming[END_REF]] an experimental validation of the MILP framework is done: test flights on an indoor quadrotor testbed demonstrate the reliability of the approach for the optimal path planner. For example, using the MILP path planner to create a plan ten seconds into the future, the quadrotor can navigate through an obstacle-rich field with MILP solve times under one second. Simple plans in obstacle-spare environments are solved in less than 50ms. A multi-vehicle test is also used to show non-communicating deconfliction trajectory planning using MILP.

Many applications in precision agriculture, disaster management and target tracking assume a collaboration between an UAV and ground-based sensors. The UAV serves as a mobile sink: it prolongs the sensors' lifetime (by canceling their need to communicate with a base station [START_REF] Xu | Energy-efficient UAV communication with multiple GTs based on trajectory optimization[END_REF]) and reduces operational costs (by canceling the need of direct human supervision [START_REF] Jawad | Energy-efficient wireless sensor networks for precision agriculture: A review[END_REF]). Such applications impose energy-based limitations in the motion planning procedure, either induced by path length [START_REF] Khan | Mobile collector aided energy reduced (MCER) data collection in agricultural wireless sensor networks[END_REF] or by communication requirements [START_REF] Xu | Energy-efficient UAV communication with multiple GTs based on trajectory optimization[END_REF]]. Furthermore, many works simplify the motion planning by assuming predefined path primitives (e.g., the UAV is constrained to move in straight, parallel lines [C. Wang et al. 2015] or in spirals [START_REF] Yue | Path Planning for UAV to Collect Sensors Data Based on Spiral Decomposition[END_REF]). Not in the least, when the environment is cluttered or uneven, the communication links may be weakened or lost as a result of signal attenuation [START_REF] Jawad | Energy-efficient wireless sensor networks for precision agriculture: A review[END_REF]]. Thus, bounds on the communication time at a waypoint have to be considered, which are difficult to handle by fixed-wing UAVs. As stated elsewhere, the result is a nonlinear (in cost and constraints) constrained optimization problem which is often

D. Illustrative example for parameterized zonotopic representation

In what follows, we introduce an example in order to highlight the applicability of the concepts discussed in Section 3.3.3. This example can be interpreted as a surveillance mission with respect to two intruders. More precisely, an agent has to monitor these two intruders and but also to remain outside their region of visibility. The agent is described by the LTI dynamics (2.19) in R d , d = 3, with the following parameters: .1) where µ = 3 and M a = 60. The agent's state is composed from position and velocity components

A = O d I d O d -µ M I d , B = O d M -1 a I d , ( D
x = p x p y p z v x v y v z ⊤ ,
whereas the input is given by the acceleration u = a x a y a z ⊤ .

Both state and input are constrained: X = {x : -15 ≤ x i ≤ 15, ∀i = 1 . . . 2d, x 2d-1 = p z ≥ 0} and U = {u : -3 ≤ u i ≤ 3, ∀i = 1 . . . 2d}. For the monitoring ability, the agent is equipped with downwards facing visual sensors [START_REF] Papatheodorou | Collaborative visual area coverage[END_REF], able to cover a conical field of view. Hence, the agent is able to survey (in a xyplane) an ellipsoidal region centered in p x p y , whose dimension depends on the value of p z and on the characteristics of the visual sensor (hereinafter, we consider the field of view angle ϕ).

The two intruders (and their visibility zones) are modeled using two convex bodies with uncorrelated motions, whose trajectories are depicted in

P k = Z(G∆ 1 , c 1 (k)) ∪ Z(G∆ 2 , c 2 (k)) (D.
2)

The scenario depicted in Figure D.1 will be considered for the test cases such that the motion of the intruders (denoted by red and blue agents) is periodic and restricted to the

3 . 3 .

 33 (a). Hyperplane arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). Perturbed hyperplane arrangement . . . . . . . . . . . . . . . . . . . . . . . Hyperplane arrangement for polytopic obstacles. . . . . . . . . . . . . . . . . .

3 . 4 .

 34 (a). d = 2 and #H = 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). d = 3 and #H = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Arbitrarily conservative approximation. . . . . . . . . . . . . . . . . . . . . . . 3.5. Zonotopic approximation for the multi-obstacle environment from Figure 3.3. . .

3 . 6 . 3 . 9 .

 3639 (a). Volume with G = G 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). Volume with G = G 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Illustration of obstacle avoidance for d = 2. . . . . . . . . . . . . . . . . . . . . (a). Polytopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). Zonotopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7. Illustration of obstacle avoidance for d = 2. . . . . . . . . . . . . . . . . . . . . 3.8. Zonotopic approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a). initial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). G = G 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The feasibility of (3.55). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  (a). infeasible for G = G 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). feasible for G 2 Q 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10. An example in R 2 for the proof of Theorem 3.3 . . . . . . . . . . . . . . . . . 3.11. Maintaining corridors procedure over Figure 3.8 . . . . . . . . . . . . . . . . .

12 .

 12 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . Z * ( P(k)) .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13. Illustrative example of the parametrization of the hyperplane arrangements.. . .

  (a). Zonotopes(3.63) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). Number of validity domains .. . . . . . . . . . . . . . . . . . . . . . . . . . 3.14. Illustration of obstacle avoidance . . . . . . . . . . . . . . . . . . . . . . . . .(a). polytopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi (b). zonotopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. A 2D cluttered environment . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. A corridor Π as in Definition 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Partition of a non-convex set . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. An example of classical Voronoi diagram in R 2 . . . . . . . . . . . . . . . . . . 4.5. Space-partitioning via Voronoi techniques. . . . . . . . . . . . . . . . . . . . .(a). The Voronoi graph for the starred points . . . . . . . . . . . . . . . . . . . (b). Partition via Voronoi construction algorithm . . . . . . . . . . . . . . . . . 4.6. Space-partitioning using separating hyperplanes. . . . . . . . . . . . . . . . ..4.7. For d = 2: a cell complex and one of its affinely equivalent polyhedron. . . . . .

4 . 9 .

 49 (a). A cell complex {X i } i=1:.... . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). The lifting X .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.8. Space-partitioning via convex lifting for a 2D environment. . . . . . . . . . . .(a). P = i=1:No P i for N o = 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). {X i } i=1:No . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The lifting L as in Corollary 4.1 and embeddings of the obstacles. . . . . . . . . 4.10. Space-partitioning {X i } i=1:No from (4.18). . . . . . . . . . . . . . . . . . . . .

  (a). N o = 31, d = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). N o = 5, d =3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.11. Beyond the convex liftable case.. . . . . . . . . . . . . . . . . . . . . . . . .

  (a). convex non-liftable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). non-convex obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12. A counter-example for liftability of a cell-complex. . . . . . . . . . . . . . . . .

  (a). Convex obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). A non-liftable partition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13. A convex non-liftable obstacle setting, its corresponding maximal subset inducing the convex liftable partition . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  (a). The convex non-liftable setting . . . . . . . . . . . . . . . . . . . . . . . . . (b). The "maximal" set S ⊂ P and its corresponding partition . . . . . . . . . . (c). The obstacle P i / ∈ S divided using entire H(S) . . . . . . . . . . . . . . . . (d). Output of Algorithm 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14. Output of a variant of Algorithm 6 over the environment in Figure 4.13 . . . .

  (a). The obstacle P i / ∈ S divided using a subset of H(S) . . . . . . . . . . . . . (b). The resulting partition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15. Partition of a convexly non-liftable cluttered environment via Algorithm 7 . . .

  (a). The obstacle P i / ∈ S divided using the facets of the partition cells . . . . . . (b). The corresponding partition . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16. Algorithm 8 run over a convexly non-liftable cluttered environment . . . . . . .

  (a). The obstacle P i / ∈ S divided using the facets of the partition cells . . . . . . (b). The partition via Algorithm 8. The highlighted facets are provided by Steps 8-13 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17. Algorithm 9 run over a convexly non-liftable cluttered environment. . . . . . .

  (a). The "altered"(after merging) collection of obstacles and the corresponding partition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). The partition corresponding to the "altered" collection of obstacles . . . . . (c). The final polyhedral partition . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18. Convex lifting of the originally non-convex obstacles and their associated partitions.100 (a). The initial non-convex obstacle setting . . . . . . . . . . . . . . . . . . . . . (b). The partition corresponding to the replacing union of convex regions . . . (c). The partition of the the originally non-convex obstacles . . . . . . . . . . . 4.19. Graph Γ(N , E, f ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  (a). d = 2 . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). d = 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20. The shortest path Path(x i , x f ) and a feasible trajectory within the corridor. . . . 4.21. Computing times w.r.t. the number of obstacles. . . . . . . . . . . . . . . . . . 4.22. The graph construction corresponding to the environment in Figure 4.18 . . . .(a). The partition via convex lifting . . . . . . . . . . . . . . . . . . . . . . . . . (b). The unprocessed graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c). The final graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23. Path(x i , x f ) and Π = Nc i=1 Π i .. . . . . . . . . . . . . . . . . . . . . . . . . . .

5. 1 .

 1 Simulation results for the environment Figure 4.19 . . . . . . . . . . . . . . .

  (a). A cluttered environment. Start position and destination. "Active" obstacles. (b). Classical MIP-based trajectory planning . . . . . . . . . . . . . . . . . . . . (c). Zonotopic obstacle over-approximation for trajectory planning via Algorithm 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d). Geometrical path computation via space partitioning . . . . . . . . . . . . (e). Corridor-constrained trajectory planning . . . . . . . . . . . . . . . . . . . 5.2. Graph obtained via convex lifting. . . . . . . . . . . . . . . . . . . . . . . . .

Π

  i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5.4. Illustration for a segment from (4.7). . . . . . . . . . . . . . . . . . . . . . . 5.5. Algorithm 12 illustration for dynamics (5.2) with d = 2 . . . . . . . . . . . . .

  (a). The terminal set xi ⊕ λX 0 f . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). The N-step BRSs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6. Algorithm 13 over the example in Figure 4.9 . . . . . . . . . . . . . . . . . . . 5.7. Input and velocity during Algorithm 13 . . . . . . . . . . . . . . . . . . . . . 5.8. A time-varying multi-obstacle environment. . . . . . . . . . . . . . . . . . . .

2 5. 9 .

 29 (a). The initial setting (and its space partitioning).The trajectory of the moving obstacle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). The superposed settings and the initial space partitioning as in Section 4.Variation of corridor width as in Algorithm 14 . . . . . . . . . . . . . . . . . .

  (a). The initial corridor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). The corridor for k = 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c). The corridor for k = 5 . .

  . . . . . . . . . . . . . . . . . . . . . . . . . (b). The corresponding corridor . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11. The decision layers for a multi-agent system. . . . . . . . . . . . . . . . . . . . D.1. Trajectories of the two intruders. . . . . . . . . . . . . . . . . . . . . . . . . .

2. 1 .

 1 Classification of the research in motion planning using MIP. . . . . . . . . . .

3. 1 .

 1 Constructive parameters and computing time of the hyperplanes arrangements given in Figure 3.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Performance of the zonotopic over-approximation techniques. . . . . . . . . . . 3.3. Simulation results for (2.20) with N p = 10, P = 10I 2d , Q = I 2d , R = I d . . . . . . 3.4. Results for the agent trajectories in Figure 3.7. . . . . . . . . . . . . . . . . . . 3.5. Simulation results for the collision avoidance . . . . . . . . . . . . . . . . . . .

4. 1 .

 1 Alternative selection of graph triple (N , E, f ). . . . . . . . . . . . . . . . . . . 4.2. A breakdown of the runtime associated to the path planning procedure. . . . . .

5. 1 .

 1 Results for the agent trajectories in Figure 5.1. Computational aspects. . . . . .

3 .

 3 The values of ∆γ i and N p min for the corridors depicted in Figure5.6. . . . . .5.4. The impact of corridor dimension (ρ(•)) on the computational features.. . . . .

Figure 1 . 3 .:

 13 Figure 1.3.: Dependencies among the thesis chapters

Figure 2 . 5 .:

 25 Figure 2.5.: Motion planning sub-areas and their interdependence

1 (

 1 Support function -[START_REF] Schneider | Convex bodies: the Brunn-Minkowski theory[END_REF]). For a given set Q ⊂ R d , the support function, evaluated at η ∈ R d , is:

Figure 3 . 1 .

 31 Figure 3.1.: Hyperplane arrangement

Remark 3 . 4 .Theorem 3 . 1 (

 3431 Notation #B denotes the number of hyperplanes and rank(B) the rank of the intersection.♦ Whitney's theorem). Let A be an arrangement in an d-dimensional space.

Figure 3 .

 3 Figure 3.2.: Hyperplane arrangement

Figure 3 . 4 .:

 34 Figure 3.4.: Arbitrarily conservative approximation.

  .40) Proof. For d = 2, the number of support hyperplanes for a zonotope becomes p(2, m) = m 1 = m. i) Introducing p(2, m) in (3.35) directly leads to (3.39).

3 .

 3 The simulations were done for different variants of the generator seed G in (3.23) for both d = 2 and d = 3 dimensions:

Figure 3 . 3 Corollary 3 . 4 .

 3334 Figure3.10.: An example in R 2 for the proof of Theorem3.3

  Figure 4.1.: A 2D cluttered environment Definition 4.1.Given the collection of obstacles P, a corridor between two points x 0 , x f ∈ int(C X (P)) is enabled by the existence of two continuous functions:

. 7 )Figure 4 . 2 .

 742 Figure 4.2.: A corridor Π as in Definition 4.1

8 )Example 4 . 1 .

 841 where d(x, y) stands for the distance between x ∈ R d and y ∈ R d . Consider the following collection of coordinates in R 2 : and d(x, y) = (xy) ⊤ (xy), the Euclidean distance. Employing a standard Voronoi construction algorithm[Aurenhammer 1991], we depicted in Figure4.4 the corresponding Voronoi regions.

Figure 4 . 5 .

 45 Figure 4.5.: Space-partitioning via Voronoi techniques.

Example 4 . 2 .

 42 We delineate in Figure4.5 an example of Voronoi graph construction based on a set of points approximating the boundary of the obstacles. The stars indicate the points taken at the faces of the obstacles and boundaries of the environment for the construction of the Voronoi graph. Running the construction algorithm for this set of points leads to a separation within the workspace, which partially verifies the conditions of a partition in Definition 4.2. While the components of the partition cover the entire workspace and fulfill the facet-to-facet property, we do not have a bijective association

11 ) 3 :

 113 Solve the convex optimization problem (

Remark 4 . 7 .

 47 Note that the construction of convex liftings based on Algorithm 4 is limited to polytopic partitions. The presence of unbounded polyhedra can be handled by adequate introduction of rays and lines in the generators representation. Despite their geometrical Lemma 4.1. A PWA lifting is continuous. Moreover, it is convex if:

Figure 4 .

 4 Figure 4.12.: A counter-example for liftability of a cell-complex.

Remark 4 . 13 .

 413 The technique in Corollary 4.3 is a formal construction which encapsulate the intuition of Example 4.2. In both cases we compute a partition (here via convex lifting, while Example 4.2 uses a Voronoi technique) w.r.t. the obstacles components and afterwards we merge the cells of the resulting partition for the the obstacle setting.

1 .

 1 Run Algorithm 8 for X = X k and P = Algorithm 8 returns a valid partition regardless of the initial convexly non-liftable obstacle setting.Proof. Basically, Algorithm 8 uses the divide-and-conquer paradigm 8 starting with Step 8 and exploiting the recursion from Step 10. Since Algorithm 5 is feasible for any two disjoint sets, Step 9 ensures the finiteness and feasibility of the algorithm.

Remark 4 .

 4 16. Similarly to Algorithm 7, the partition identified by Algorithm 8 contains a couple of non-convex cells corresponding to the divided obstacles. However, the partitions resulted from these algorithms are different.♦

Definition 4 . 6 .

 46 A weighted graph Γ (N , E, f ) is defined by the triple (N , E, f ) with N the set of nodes, E the set of edges, and f : E → R a function which associates to each edge a real positive value (its weight). Remark 4.20. A variation of the Definition 4.6 can employ the weight function f (

Remark 4 .

 4 22. Solving a parameterized optimization problem, as is required in Steps 3 and 4, is often difficult if not downright infeasible. Standard approaches are to sample the continuous variables (θ in our case) or to project along a basis function and reformulate the problem in terms of the associated weights.♦

Figure 4 .

 4 Figure 4.21.: Computing times w.r.t. the number of obstacles.

Contents 6 . 1 . 6 . 2 .

 6162 Conclusions . . . . . . . . . . . . . . . . . . . . . . . .. . 129 Future directions. . . . . . . . . . . . . . . . . . . . . . . 131

R 4 :

 4 |x| ≤ 15 15 0.35 0.35 ⊤ } U = {u ∈ R 2 : |u| ≤ 10 10 ⊤ } In Figure 5.1 we depicted the agent trajectory, generated by each of the four methods enumerated in the beginning of Section 5.1. The agent has the same initial and final position, both shown in Figure 5.1a: x 0 = -13.23 -12.

Figure 5

 5 Figure 5.3.: Π = Nc i=1

N 1 :f 3 :

 13 Solve the following LP:maxλ λ s.t. xi ⊕ λX 0 f ⊂ Πi , λ < 1. Ri 0 = xi ⊕ λX 0 for j=1:N do Ri j = A -1 Ri j-1 ⊕ (-BU ) ∩ Πi .4: end for 2. (ON-LINE) Apply Algorithm 13

  (x 0 , x f ) = {x 0 , x1 , . . . , xNc }, {N i p min } i=1:Nc , {X f (x i )} i=1:Nc , 1: x 0 = x0 2: for each Πi ⊂ Π do

2 :

 2 Determine γ(θ) from Definition 4.1 for Path(x 0 , x f ). 3: for k ∈ K do 4:

Figure D. 1 .

 1 Since the motion is within the xy-plane, the resulting interdicted regions can be defined as zonotopes by taking a common generator matrix 1 scaling factors ∆ 1 = diag( 2.5 3 ) and ∆ 2 = diag( 1 0.5 ), while the centers, c 1 (k) and c 2 (k) respectively, follow the trajectories depicted in Figure D.1. For further use we denote:

  1.1. Illustration of a boat moving within a fjord region.. . . . . . . . . . . . . . . . 1.2. Illustration of a trivial obstacle avoidance alternatives. . . . . . . . . . . . . . . 1.3. Dependencies among the thesis chapters . . . . . . . . . . . . . . . . . . . . . 2.1. Illustration "big-M", M computed as in Remark 4. . . . . . . . . . . . . . . . . 2.2. Polytopic obstacles (in gray) and the corresponding hyperplane arrangements. . . 2.3. A circular obstacle and its polytopic approximation. . . . . . . . . . . . . . . . 2.4. The illustrative example for GDP formulation. {L i } i=1:4 hyperplanes . . . . . . 2.5. Motion planning sub-areas and their interdependence . . . . . . . . . . . . . . 2.6. The solution of the LP relaxation is exactly the solution of ILP. . . . . . . . . .

3.1. Hyperplane

arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

(a). Hyperplane arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b). The central (sub-) arrangements . . . . . . . . . . . . . . . . . . . . . . . .

3.2. Hyperplane arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 5 .

 5 1.

1. Introduction Contents 2.1. MIP Formulations

  . . . . . . . . . . . . . . . . . . . . . . 14 2.1.1. Generalized Disjunctive Programming . . . . . . . . . . 14 2.1.2. Geometric viewpoint for MIP . . . . . . . . . . . . . . 15 2.1.3. GDP formulation for region description . . . . . . . . . 20 2.

2. Optimization-based control framework

  . . . . . . . . . . . 21 2.2.1. Prerequisites . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.2. Model predictive control . . . . . . . . . . . . . . . . . 22 2

.3. Standard MIP-based problems in motion planning

. . . 24 2.3.1. Task assignment . . . . . . . . . . . . . . . . . . . . . 27 2.3.2. Path planning . . . . . . . . . . . . . . . . . . . . . . 29 2.3.3. Trajectory planning . . . . . . . . . . . . . . . . . . . 32 2.3.4. Obstacle and collision avoidance . . . . . . . . . . . . . 35 2.3.5. Connectivity maintenance and formation control . . . . 39 2.4. Conclusions and open problems

. . . . . . . . . . . . . . 43

,

  Prodan, I., Stoican, F., Olaru, S., Niculescu, S. I. (2019, July).

	Complexity bounds for obstacle avoidance within a zonotopic framework. In 2019 American Control Conference (ACC) (pp. 335-340). IEEE.
	• Ioan, D., Olaru, S., Prodan, I., Stoican, F., Niculescu, S. I. (2019, July). Parametrized hyperplane arrangements for control design with collision avoid-ance constraints. In 2019 IEEE 15th International Conference on Control and Automation (ICCA) (pp. 1591-1596). IEEE.

,

  Olaru, S., Niculescu, S. I., Prodan, I., Stoican, F. (2019, June).

	• Ioan, D., Olaru, S., Prodan, I., Stoican, F., Niculescu, S. I. (2019). From Obstacle-Based Space Partitioning to Corridors and Path Planning. A Convex Lifting Approach. IEEE Control Systems Letters, 4(1), 79-84.
	Navigation in a multi-obstacle environment. from partition of the space to a zonotopic-based MPC. In 2019 18th European Control Conference (ECC) (pp. 1772-1777).

2. Mixed-integer programming in motion planning Contents 3.1. Zonotopic framework

  . . . . . . . . . . . . . . . . . . . . 46 3.1.1. Prerequisites . . . . . . . . . . . . . . . . . . . . . . . 46 3.1.2. Hyperplane Arrangements . . . . . . . . . . . . . . . . 48 3.1.3. Illustrative examples . . . . . . . . . . . . . . . . . . . 51 3.1.4. Zonotopic framework . . . . . . . . . . . . . . . . . . 53 3.1.5. Set inclusions for parametrized zonotopes . . . . . . . . 54 3.1.6. Measures for zonotope approximations . . . . . . . . .

55 3.2. Complexity bounds within a zonotopic framework

  . . . 56 3.2.1. Hyperplane arrangement induced by the zonotopes . . . 56 3.2.2. Illustrative examples . . . . . . . . . . . . . . . . . . . 59 3.2.3. Numerical simulations. Comparisons. . . . . . . . . . .

61 3.3. Control design perspective

  . . . . . . . . . . . . . . . . . 64 3.3.1. Corridor constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.3.4. Neighborhood-aware MPC implementation . . . . . . .

	64
	3.3.2. From polyhedral to zonotopic MPC for collision avoidance 68
	3.3.3. Parametrized zonotopic representation within time-varying multi-obstacle environments.

73 3.4. Conclusions

  The logical function h(Y ) corresponds to logical decisions in terms of boolean algebra and it is expressed in conjuctive normal form (CNF-"product of sums"). The idea of GDP can be summarized as follows: If Y jk = true, then the constraints g jk (x) ≤ 0 and c k = γ jk are imposed. Otherwise, they are ignored.It is worth mentioning that the disjunction j∈J k represents, in fact, an exclusive relationship. That is, in a disjunction indexed by k only a boolean variable Y jk can be true. Some formulations state this requirement explicitly (as a separate constraint), others include it implicitly in the boolean function Ω(Y ).

	The problem (2.1) can be written as MIP, employing the binary variables y jk ∈ {0, 1} instead of the boolean ones Y jk and replacing the constraint (2.1c) with:
	min x	f (x) +	k	y jk γ jk	(2.2a)
	s.t. r (x) ≤ 0,		(2.2b)

  Remark 2.2. It is worth discussing the value of the coefficient M .

	1 ([Williams 2013]) consists in choosing a very large positive number M which plays the role of a relaxation constant ([J. P. Vielma and Nemhauser 2011; A. Richards and J. How 2005]). The binary variables play the role of a "switch", activating/deactivating the corresponding constraints. This is a powerful artifact able to encode logical conditions but one should proceed with caution in selection of the value of M. A too large value may hinder the resolution of the MIP. For a detailed analysis regarding the "big-M" technique, we refer the interested reader to, e.g. [J. P. Vielma and Nemhauser 2011; Hooker 2011]. ♦

  For instance in [M. G. Earl and R. D'Andrea 2005], circular obstacles are considered as the one depicted in Figure 2.3, each being determined by its radius R o and by the coordinates of its center (x o , y o ):

  The considered problem is trajectory-generation with obstacle-avoidance requirements and minimum-time trajectory-generation problems. These problems involve vehicles that are described by mixed logical dynamical [AlbertoBemporad and Morari 1999] equations (a hybrid representation). The algorithms use fewer binary variables than standard MILP methods, and require less computational effort.Similarly, in [Matthew G. Earl and Raffaello D'Andrea 2005] the iterative MILP methods are motivated on problems derived from an adversarial game between two teams of robots.

Also, in [M.

[START_REF] Earl | [END_REF]

iterative MILP algorithms are presented, algorithms which address the issue of the MILP coping with large-scale models.

  Galvao 2009; A. Richards and O. Turnbull 2015;[START_REF] Galvao | motion planning using MIP Recalling that R k is centered in 0: CA k x 0 ∈ P j ⊕ R k , (3.68) where k =[END_REF] which discuss over-approximated corner cutting constraints and provide constructive details: [Maia and R. K.Galvao 2009] provides the initial construction; [A.[START_REF] Richards | Inter-sample avoidance in trajectory optimizers using mixed-integer linear programming[END_REF]] and [R. J.[START_REF] Afonso | Reduction in the number of binary variables for inter-sample avoidance in trajectory optimizers using mixed-integer linear programming[END_REF] improve it by reducing the number of necessary constraints and by reducing the number of necessary binary variables, respectively.[START_REF] Stoican | Exact and overapproximated guarantees for corner cutting avoidance in a multiobstacle environment[END_REF] goes further by providing exact and over-approximated descriptions of the under-shadow (and of its complement, the visible) region generated by multiple agents within a multi-obstacle environment.

Contents 4.1. From obstacle-based partitioning to path generation

  . . 78 4.1.1. Prerequisities . . . . . . . . . . . . . . . . . . . . . . 78 4.1.2. Proposed Approach . . . . . . . . . . . . . . . . . . . 81 4.

2. Convex lifting as a tool for space partitioning

  . . . . . . 84 4.2.1. Convex Lifting. Background. . . . . . . . . . . . . . . 84 4.2.2. Space partitioning via Convex Lifting . . . . . . . . . . 86 4.2.3. Beyond the convex liftable case . . . . . . . . . . . . . 91 4.

3. From paths to corridors

  . . . . . . . . . . . . . . . . . . . 101 4.3.1. Geometric path generation based on polyhedral partitions 101 4.3.2. Path generation based on non-convex polytopic partitions 105 4.3.3. Corridors as patchy polyhedra . . . . . . . . . . . . . . 107 4.

4. Conclusions

  Remark3.10. Taking the volume as the best measure for zonotope approximation,(3.30) is nonlinear but exact and (3.31),(3.32) are linear but coarse. Within this hierarchy is natural to consider the Hausdorff distance[START_REF] Ziegler | [END_REF]] as an intermediate measure (in terms of complexity and precision). This direction is not followed in this work. As it is challenging to obtain an explicit expression of the distance in terms of the zonotope parameters the problem remains open, but Proposition 3.1 offers a list of relevant candidates in this respect.♦

	.	♦
	Solving (3.29) with (3.30) is burdensome in higher dimensions. This was the reason for considering the simpler (but less accurate) measures (3.31) -the Manhattan norm and (3.32) -the infinity norm.

Table 3 .

 3

	d = 3	δ ∞ Vol Vol δ 1 δ ∞ δ 1 δ ∞ Vol Vol	.27) (3.27) (3.27) (3.27) δ j = δ (3.28) (3.28) (3.27) (3.27) (3.27) (3.27) δ j = δ	G 1 G 2 G 3 G 1 G 2 G 3 G 1 G 2 G 3 G 1 G 2 G 3 G 1 8,44 8,13 8,02 8,27 8,19 8,01 8,19 9,40 9,19 9,30 9,08 8,88 9,05 G 2 8,24 G 3 8,44 G 1 8,68 G 2 8,50 8,61 G 3 G 4 9,82 G 5 10,50 120 62480 42 505 28 225 42 534 42 441 28 225 42 441 42 510 28 225 40 530 42 479 28 225 42 531 42 504 28 225 38 472 42 431 28 210 42 532 60 8400 G 4 9,71 60 8000 G 5 10,59 120 51396 G 4 11,02 60 8400 G 5 11,79 84 24528 G 4 10,49 60 7220 G 5 11,07 120 49200	20,53 -46,30 27,45 5,25 -46,30 5,25 21,72 -46,30 26,49 14,32 -46,30 26,73 9,41 3,77 10,05 8,72 3,64 10,90 -85,31 9,26 -86,01 -10,12 -85,31 -57,11 -87,37 -13,96	9,53 3,81 10,09 8,19 3,91 8,19 9,66 3,75 10,27 9,94 3,84 11,28 196 101 159 402 123 266 105,98 1145,09 4952 197 101 167 374 175 374 199 101 169 481 193 272 287 143 235 523 161 370 934 100,22 1127 932,12 9432 105,07 934 413,07 2218 90,54 1124 867,08 9436	376,98 410,07 368,93 897,92 583,33 897,92 368,41 410,07 374,99 897,92 583,33 628,72 373,78 410,07 383,41 902,85 583,33 661,92 1857,46 2019,34 2623,3 5852,55 1857,46 1908,2 2623,3 5852,55	71,7 86,78 68,04 308,98 165,69 308,98 67,8 86,78 70,8 308,98 165,69 186,37 70,25 86,78 74,63 311,23 165,69 201,49 323,56 360,47 498,19 1234,56 323,56 335,13 498,19 1234,56

2.:

Performance of the zonotopic over-approximation techniques.

  a significantly reduced set of obstacles, where Path(x k , x f ) denote a set of points from the feasible space connecting x k and x f .

Remark

3.16

. In order to identify the set Path(x k , x f ) various method can be employed, ranging from the earlier method of visibility graph

[START_REF] Lozano-Pérez | An Algorithm for Planning Collision-free Paths Among Polyhedral Obstacles[END_REF]

to Voronoi-like techniques, as in [R.

[START_REF] Afonso | Waypoint trajectory planning in the presence of obstacles with a tunnel-MILP approach[END_REF]

.

111 5.2. Relay MPC. Feasibility guarantees

  . . . . . . . . . . . . . . 115 5.2.1. Prototype MPC with recursive feasibility guarantees . . 115 5.2.2. Relay MPC . . . . . . . . . . . . . . . . . . . . . . . 118

5.3. Adjustments/Adaptation for time-varying cluttered en-

vironments . . . . . . . . . . . . . . . . . . . . . . . . . .

121 5.4. Extensions to multi-agent systems

  .. . . . . . . . . . . .

126 5.5. Conclusions

  4 and the continuity/convexity conditions, the construction of z(x) is summarized by Algorithm 4. Construction of a Convex Lifting for a given polyhedral partition {X i } i∈I of a polytope X Input: {X i } i∈I and a given constant ξ > 0. Output: (a i , b i ), ∀i ∈ I 1: Extract all pairs of neighboring regions in {X i } i∈I 2:

	Algorithm 4

Table 5 .

 5 1.: Results for the agent trajectories in Figure5.1. Computational aspects.

	Approach Figure N goal t goal (sec) ℓ t (m) 1 5.1b 340 10542 20.41 2 5.1c 457 6.1457 24.52 3 5.1d 415 5.15 21.77 4 5.1e 390 4.35 21.14
	Approach Fig. 1 5.1b non-convex Convexity 2 5.1c non-convex 3 5.1d convex 4 5.1e convex	Convergence Obstacle avoidance ✓ ✓(path and navigation) ✓ ✓(if a feasible path is available) ✓ ✓(path), ×(navigation) × ✓(path and navigation)

Table 5 .

 5 2.: Feasibility and convergence aspects for approaches in Table5.1.

  . In this way, one can manage the trajectory

	i ∆γ i N p min	1 0.77 4.02 0.09 1.27 4.50 2.39 0.53 2 3 5 8 10 11 11 24 13 14 29 19 13

Table A .

 A 1.: Classification of the alternatives of MIP in motion planning Method/Technique

Table B .

 B 1.: Classification of MIP approaches in motion planning Agent(s) Centralized Decentralized Distributed References(e.g.) Matthew G. Earl and Raffaello D'Andrea 2005; Janeček, Klaučo, Kalúz, et al. 2017; Haghighi, Sadraddini, and Belta 2016], [X. Sun and Cassandras 2015; Ritter et al.

	1 . . . N 1 . . . N N = 1	✓ × ✓	✓ × ×	× ✓ ×	2014] [X. Wang et al. 2015] [Papen et al. 2017] [Molinari, Anh, and Re 2017; Ragi and Mittelmann 2017]

N = 3 / N ≥ 3,partially ✓ × × [M.

Chen, Shih, and Tomlin 2016] 1 . . . N ✓ × × [

Table C .

 C 1.: Existing software for MIP implementation Programming Language Modeling Language Solvers(e.g.)

	MATLAB	Yalmip	CPLEX, GUROBI, SCIP, CBC, GLPK, LPSOLVE, INTLINPROG, MOSEK, XPRESS
		AMPL	CPLEX, GUROBI
	PYTHON	CVXPY	GUROBI, CVXOPT, Ele-mental, XPRESS, GLPK, CBC, MOSEK
		Pyomo	GUROBI, GLPK, CPLEX, Cbc, IPOPT
	Julia	JuMP	CPLEX, GUROBI, SCIP

Table C .

 C 2.: Software for MIP implementation-examples Software References(e.g.) Programming Language MATLAB

  [START_REF] Richards | Trajectory planning[END_REF] A. Richards and J. P How 2002]

See, e.g., Figure1.1, where the complexity is fairly amplified.

It is also referred to as "big-D" or "big-R".

Due to the physical limitations.

Sets of indices.

There exists a variety of technical procedures like cell merging[Ionela Prodan, Stoican, et al. 2015] or logarithmic formulations [J.P. Vielma and Nemhauser 2011] which may be employed to reduce the complexity of the formulation.

The disturbance will be considered bounded w(t) ∈ W, and W ⊂ R nw to be a convex and compact set.

We replace(2.22d) with the relaxation x ij ∈ [0, 1]

In the sense that the time is not explicitly associated to the arch of curve in the output or space.

For example, in the sense of the economical MPC.

A non-signalized intersection follows the first-come-first-serve (FCFC) principle.

For the clarity of the presentation, in what follows we will neglect the equations(2.29b) and (2.30b) but the reasoning can be readily adapted.

W.r.t the multi-agent system.

Partial Differential Equation

For compactness, whenever clear from the context we use notation A.

If not explicitly stated otherwise, we consider the later case since the former is a simplification of the later.

m(m -1). i) Introducing p(3, m) in(3.35) directly leads to(3.42).

An exemplification in R 2 is provided in Figure3.10

Approaches related to the construction of a graph structure.

They can be addressed using the techniques developed in Chapter 5.

For simplicity the cell complex can be understood as a union of non-overlaping polyhedra enjoying a facet-to-facet property.

By adding obstacles one by one in the lifting procedure.

An often used pattern, "divide et impera" relies on the idea to decompose a given problem into two or more similar and simpler sub-problems, to solve them in turn, and to compose their solutions to solve the given problem.

More precisely, we generate randomly the vertices describing the obstacles

The trajectory point have to stay inside the given corridor.

A solution of the original non-convex problem can become infeasible for the sequence of the convex sub-problems.

According to the air traffic rules, the aircraft have to be separated by at least 5 nautical miles horizontally and 1000 ft vertically, otherwise they are in conflict.

 3 A more detailed review on the air conflict detection and resolution can be found in[START_REF] Kuchar | A review of conflict detection and resolution modeling methods[END_REF]].

Unsignalized intersections.

Computational time, feasibility etc.

Except AMPL.

The most encountered.

For the sake of clarity, we considered two "boxes", but the reasoning and simulation are generic and apply also for more complex shapes of the obstacle, in as much as they are (approximated by) zonotopic sets.

Acknowledgements

The work presented in this thesis was carried out as a collaboration between two laboratories, L2S from CentraleSupélec and LCIS from Grenoble INP, with funding from the DGA, Ministry of Defense.

to which correspond the following bounds: .35b) Proof. To apply Whitney's theorem and obtain the bounds (3.20), the key is to count how many central arrangements of rank k there are. Hence, we need to avoid selecting parallel hyperplanes from the available p(d, m) • 2N o . We consider each of the 2N o families of p(d, m) hyperplanes and select as follows: from the 1st family we select j - 1 ; in the 2nd family there remain mj - 1 available (the others are parallel to those already selected and thus unsuitable) from which we select j + 1 ; the procedure continues up to the 2N 0 -th family from which we select j +

No . Enumerating all sequences j - 1 , . . . , j + No which add to k and taking all the possible combinations for a fixed sequence provides the total number of central arrangements of rank k: .36) Writing explicitly the combinatorial terms for a sequence j -

No ! as the coefficient from the multinomial theorem [START_REF] Hazewinkel | Multinomial coefficient, Encyclopedia of Mathematics[END_REF]] allows to rewrite (3.36) as .37) Introducing (3.37) in (3.18) leads directly to (3.34) which, applied as in (3.19), leads to (3.35), thus concluding the proof.

These results allow to further derive a bound for the number of generators in the zonotopic representation and to provide a rigorous bound for d = 3.

Corollary 3.1. Assuming n * o support hyperplanes in (3.22), for any m ∈ N + which verifies: .38) the arrangement A( H) has fewer cells than A(H).

Algorithm 3 Active obstacles selection

Input: P, N pred , A, B, C, U Output:

The following simulation scenario highlights the link between Algorithm 3 and the Theorem 3.3. Let us consider a new configuration of the obstacles (depicted in Figure 3.14), for which one cannot find an over-approximation maintaining all corridors open. By selecting at most"d + 1" obstacles and using Theorem 3.3 we find a proper approximation with all corridors unobstructed. We proceed by applying a MPC optimization problem over the non-convex constraints described by the earlier defined topologies. Further, using the Algorithm 3 we select only the obstacles intersecting the trajectory and for the zonotopic case we use Algorithm 1 in order to obtain separation. We delineate in the Table 3.5 the number of steps in which the agent attain a vicinity of the origin (denoted as N goal ), the corresponding computing time (t goal ) and the MPC parameters.

Conclusions

In this chapter, the obstacle avoidance problem was treated with the aim to reduce the complexity of the original problem. We first recalled the set theory notions directly related to the complexity evaluation from a geometrical perspective. Basically, a particular family of sets (parametrized zonotopes) were employed in order to efficiently overapproximate the initially polytopic region of interest. To this end, we consider as an explicit measure the problem complexity (number of cells) and we concentrate exclusively on over-approximation aspect and on its benefits on the MPC(Model Predictive Control)based collision avoidance problem. Additionally, we compared this optimization-based

Appendices

danger zones of the other ones. The architecture of the "safety controller" is centralized and it guarantees, for N = 3, that every agent avoids entering into any other danger zones. Therein, MIP is used to provide the higher level control logic to synthesize the cooperative "safety controller". The objective is to group the agents in a pair such that the corresponding avoidance maneuvers may not lead to a dangerous configuration for the other agents.

For the generic obstacle and collision avoidance problem (in MI-MPC framework Section 2. 3.4), in the corresponding OCP for the global system the dynamic behavior of the individual agent is coupled by means of the cost function and the constraints. Also, perfect knowledge of each agent dynamics (described by equations) is available to all the other agents. Consequently, the global model will be used in a predictive control context which permits the use of non-convex constraints for collision avoidance behavior.

A non-MPC example, [X. Wang et al. 2015] the problem studied in this paper is the collision avoidance in a multi-robot system. The approach is quasi-similar with one of previously-mentioned papers [Haghighi, Sadraddini, and Belta 2016]. The workspace is decomposed in equally-shaped cells and each cell cannot be visited by two robots at the same time. Each robot has to complete its own task by choosing from a set of possible trajectories. These trajectories are described by a sequence of adjacent cells and by the crossing time (in which an agent passes through respective cell). In this formulation the collision-avoidance problem becomes a resource allocation problem (the cells can be considered a shared resource). Assuming that the control of the agents is independent (i.e. they cannot pause their movement in order to give or to receive priority), the idea presented in this paper is to compute an initial delay time for each agent such that no collision to appear. In the centralized control strategy the trajectory is chosen by a central unit. Hence, the control algorithm has to return, besides the initial time delay, an "optimal" trajectory, the objective remaining unchanged. The associated optimization problem is MIP, binary variables modeling the disjunction, resulted from the corresponding resource allocation problem. Noteworthy, the objective is a min max one because it is assumed that the robots work in parallel and the shortest (min) time for finishing the movement is the time of the slowest (max) robot. The objective can be reformulated as a standard minimization through the tools from, e.g., [START_REF] Scott | Dynamic optimal power flow in microgrids using the alternating direction method of multipliers[END_REF].

While the complexity of MIP formulations increase in worst case situations exponentially w.r.t the number of binary variables, the reliability of MIP for real-time implementation decreases in the case of higher dimensions of the global systems. The same drawback is valid for the increase of the number of agents regardless of the formulation's efficiency.

B.2. Distributed

The underlying idea of distributed control approaches [START_REF] Maestre | Distributed model predictive control made easy[END_REF]] is to divide the global control problem into a certain number of sub-problems, each involving a specific collection of local controllers or agents. Hence, each agent does not have access at the global information, but can be partially informed on the behavior of the other components of the local subsystem.

In large-scale multi-agent systems, where the agents are dispersed within a given workspace, it is more convenient to handle a set of smaller and/or simpler problems than treating the complex global system. The overall control strategy is given by the behavior of the local controllers, which may have a cooperative interaction. There are several advantages ). This holds likewise for input and for sets X xy , U xy , which are projections of the initial sets on the xy-subspace. Furthermore, M represents a sufficiently large constant, ensuring the penalty of the constraints avoidance to be further translated as reference for the tracking on the z-coordinates.

The existence of the common "seed" G in (D.2) allows us to rewrite the control problem (D.4) obtaining a parametrization with respect to the positions of the centers. Further using a similar procedure like in the case of (2.21c), we rewrite (D.5c) as mixed-integer constraints like in Section 2.1 and the resulting MI problem can be solved using, e.g., CPLEX [START_REF] Cplex | V12. 1: User's Manual for CPLEX[END_REF]].

Regarding the block MPC z , the control strategy relies on resolution of (D.6). where A z , B z are selected from the extended system (D.1). The reference zref|k has the following components: p ref z|k which is the sum between the reference given by the visual area