Patent Publication Number: US-2023134026-A1

Title: Method and device for coordinating vehicle routes in confined areas

Description:
TECHNICAL FIELD 
     The present disclosure relates to the field of automatic vehicle control. In particular, it proposes techniques for coordinating the movement of multiple vehicles in an environment with conflict zones. 
     BACKGROUND 
     It is believed that fully automated vehicles (AV) have the potential to drastically change the transport industry, both in terms of increased safety and efficiency [1]. The most drastic improvements are expected when a substantial part of the vehicles on public roads are fully automated, as in e.g., “robo-taxis” and hub-to-hub transports on highways. Unfortunately, managing the unpredictable conditions on public roads in a reliably safe manner has proven to be harder than initially expected, and the current state-of-the-art exhibits a lack of production-level maturity. 
     However, confined areas, such as mines, ports and logistic centers, lack many of the difficult aspects of public road driving, and present use cases for near-future, large-scale deployment of automated vehicles. Within this context, the AVs form a component in transport solutions for commercial operations, where the material-flow is largely handled without human involvement. 
     One of the challenges in such systems is the efficient coordination of multiple AVs use of mutually exclusive resources, such as intersections, narrow roads, work-stations (e.g. crushers, loading/unloading spots etc.) and, in the case of electrified AVs, charging-stations. Poor coordination can lead to substantial decreases in energy-efficiency and productivity, reducing the benefits of automation. In detail, failure to coordinate the traffic situation efficiently has the following implications:
         Energy waste: With poor planning, the vehicles expend more energy than required to perform the transport missions. An efficient planner, on the other hand, prevents vehicles from coming to unnecessary stops at, e.g., intersections, by pre-emptive speed adjustments. Moreover, an efficient planner harvests potential energy better by, e.g., prioritizing vehicles going downhill in the MUTEX management, avoiding unnecessary deceleration. Energy efficiency is particularly relevant for Battery-Electric Vehicles (BEV), where poor performance leads to more frequent charging stops where no value is added to the transport solution.   Underutilization: With poor planning, vehicles are often stuck in states where no value is added. Examples are vehicles waiting to enter charging stations or loading/unloading areas. An efficient planner avoids this by performing MUTEX management and speed planning in such a way that waiting times are reduced, leading to more value increasing time spent in the transport solution and a smaller vehicle fleet.   Deadlocks: Some road network traffic layouts are prone to deadlocks, a situation where the vehicles cannot continue their missions without manual intervention in the form of reset. With a deadlock prone traffic layout, such situations will inevitably occur unless a traffic planner restricts access to certain areas at key instants.   Safety: With poor planning, it is possible for the system to be near dangerous states (e.g., risk of collision or similar) such that lower-level emergency systems are activated by small disturbances. An efficient planner avoids this by including sufficient margins.   Physical failure states: With poor planning, vehicles may enter states from which they cannot recover. An example is the complete stop of a heavily loaded vehicle in a steep incline. If stopped, the engine/motor and transmission in use can be such that it cannot be re-started. An efficient planner avoids this by pre-emptive speed planning.   Hardware wear: An efficient planner avoids unnecessary wear of the machinery, which for instance could be induced when heavy vehicles are decelerated with the service breaks repeatedly.       

     The problem of handling mutual exclusive resources has been addressed for industrial robots [16], [17], where different scheduling algorithms have been explored. The coordination of multiple AVs, however, adds a different aspect to the challenge, where for example, the dynamics of the vehicles and the road topology play a significant part in the optimization problem. The coordination of automated vehicles at intersections has been widely discussed in the literature recently, see [2] for a comprehensive survey. In general, the problem is difficult to solve, and was formally shown to be NP-hard in [15]. Often relying on simplifying assumptions and heuristics, a number of methods has been presented that solve the problem using, e.g., hybrid system theory [3], reinforcement learning [4], scheduling [5], model predictive control (MPC) [6], [7] or direct optimal control (DOC) [8], [9]. 
     Coordination of AVs in confined areas has some distinct differences compared to the intersection scenarios often found in the literature. For instance, the full site-layout of confined areas is typically known at the planning stage, and it can often be expected that no non-controlled actors will disturb the execution of a plan once its formed. The motion of each vehicle can therefore be planned from the start of a transport mission to its end. Planning for confined sites is therefore benefited by methods that can handle long planning horizons. This in contrast to the intersection coordination context found in the literature, where a cutout around the intersection proper is most often considered, with vehicles arriving at speed at the intersection cutout [6], [7]. Moreover, in confined areas, the motion plan must consider coordination of the vehicles&#39; use of a number of resources besides intersections, e.g., different variety of collision zones, and the quality of the plan has a direct and significant impact on the profit margins of the transport solution provider. 
     Direct optimal control is contemplated as a favorable type of method for solving the optimization problem. In the present disclosure, a two-stage heuristic for obtaining the optimal crossing order for each collision zone is utilized. The idea proposed in [9] is adapted to suit the goals in this application. The problem is formed in the space frame as it allows for optimizing the motion over the full vehicle path, as well as, allowing a more convenient handling of the collision constraints that occur in the collision zones. Furthermore, the collision constraints for the merge-split collision zone are defined such that the zone is not excluded for only one vehicle at a time, and the algorithm is modified to be able to handle these zones along with being able to handle multiple collision zones between vehicles. Although the approach focuses on confined sites, the method of handling the mutual exclusion zones is extendable to other scenarios as well (e.g., public road applications). 
     From similar scenarios that have been considered, the authors in [11] and propose an optimization approach for handling merge scenarios, which is a subset of the merge-split collision zone and [12] uses a game theoretical strategy for optimizing traffic flow through multiple intersection collision zones. The handling of multiple intersections and zones of different types were also identified in the survey [2] as topics for further work in this field. 
     SUMMARY 
     One objective of the present disclosure is to make available more widely applicable vehicle coordination techniques. A further objective is to enable planning over the entire site (a full plan covering the entire navigable environment of the vehicles) with at most a feasible effort increase. A further objective is to enable planning at all moments with at most a feasible effort increase. In particular, it is desired to remove the time horizon of the OCP and allow it to extend to the end of each vehicle&#39;s assigned transport mission. A further object is to enable planning for environments with multiple conflict zones (CZs). A further object is to enable planning for environments with at least one dwelling zone (e.g., charging, tipping, servicing), which is such that special vehicle dynamics apply to any vehicle therein, and/or at least one merge-split or narrow-road CZ. 
     At least some of these objectives are achieved by the invention as defined by the independent claims. The dependent claims relate to embodiments of the invention. 
     In a first aspect of the invention, there is provided a method of planning routes for a plurality of vehicles operating in a common environment which includes at least one CZ, wherein movements of each vehicle are controllable by a control signal. The method comprises: obtaining a predefined objective function J; solving a first optimization problem for a first objective function V i (T i ) derived from the predefined objective function, to obtain a vehicle crossing order at each CZ, wherein the first optimization problem is subject to safety constraints; and solving an optimal-control problem (OCP) for the predefined objective function subject to the vehicle crossing order at the CZs and subject to the safety constraints, to obtain a control signal u i  for each of the vehicles. The OCP includes a dynamic vehicle model representing evolution with respect to path length p i ,in which time t i (p i ) and path speed v i (p i ) are state variables. The OCP includes such a model, for example, if the objective function for which it is being solved relates to a dynamic vehicle model representing evolution with respect to path length (p i ), in which time (t i (p i )) and path speed (v i (p i )) are state variables. 
     The invention according to the first aspect enables the following advantages compared with the prior art. Because the OCP is formulated with path length (distance along the path) as the independent variable, rather than time, it becomes possible to plan over the entire site at all moments without significant effort increase. The use of path length also avoids certain numerical difficulties which affected the prior art, relating to requirements for computational resource and the complexity in formulating the OCP. Further, it becomes possible to handle CZs of other types than intersections, namely, the merge-split combination, charging stations, loading/unloading stops, tunnels and other single-lane passages. It becomes possible to account for all CZs present at a site. The use of path length as independent variable further allows the site state and topography (e.g., the road slope, the friction conditions on different parts, etc.) to be accounted for. 
     In some embodiments (First Approach), the this is implemented by letting the CZ entry and CZ exit times be decision variables of the first optimization problem. Alternatively, the first objective function is a parametric local optimum of the predefined objective function, wherein the parametric local optimum is parametrized by tentative CZ entry and CZ exit times and is independent of the safety constraints. 
     In other embodiments (Second Approach), the route planning is implemented by letting the pairwise relative vehicle crossing orders (b) as well as a state trajectory and control signal for each of the vehicles (W) be decision variables of the first optimization problem. 
     In a second aspect of the present invention, there is provided a method of controlling a plurality of vehicles operating in a common environment which includes at least one CZ, wherein movements of each vehicle is controllable by a control signal, the method comprising: sensing current positions P of all vehicles; obtaining a predefined objective function J; solving a first optimization problem for a first objective function V i (T i ) derived from the predefined objective function, to obtain a vehicle crossing order at each CZ, wherein the first optimization problem is subject to safety constraints and subject to the sensed vehicle positions; solving an optimal-control problem (OCP) for the predefined objective function subject to the vehicle crossing order at the CZs and subject to the safety constraints, to obtain a control signal u i  for each of the vehicles; and feeding the control signal to the vehicles. 
     In a third aspect of the present invention, there is provided a device comprising memory and processing circuitry which is configured to perform the above method of planning routes for a plurality of vehicles operating in a common environment which includes at least one CZ, or to perform the above method of controlling these vehicles. The invention further relates to a computer program containing instructions for causing a computer, or the device in particular, to carry out the above methods. The computer program may be stored or distributed on a data carrier. As used herein, a “data carrier” may be a transitory data carrier, such as modulated electromagnetic or optical waves, or a non-transitory data carrier. Non-transitory data carriers include volatile and non-volatile memories, such as permanent and non-permanent storage media of magnetic, optical or solid-state type. Still within the scope of “data carrier”, such memories may be fixedly mounted or portable. 
     The second and third aspects of the invention generally share the effects and advantages of the first aspect, and they can be implemented with an equal degree of technical variation. 
     The term “conflict zone” may include a multiple-exclusion (MUTEX) zone. 
     The term “path length” refers to the vehicle&#39;s position constrained to a predefined path. As such, even though the vehicle is designed for two-dimensional motion, the path length is a one-dimensional quantity, which can be represented, e.g., as the distance p i  from a reference point on the path. Similarly, path speed v i  can be understood as the time derivative of the path length. 
     The term “route” may refer to a representation indicating a vehicle&#39;s position (e.g., a point in space) at different times. A “path” may refer to the set of positions that a vehicle has traversed or will traverse, without indicating when this occurred or will occur. As such, a route may be considered to be a function from time to space, whereas a path is a static quantity. 
     Generally, all terms used in the claims are to be interpreted according to their ordinary meaning in the technical field, unless explicitly defined otherwise herein. All references to “a/an/the element, apparatus, component, means, step, etc.” are to be interpreted openly as referring to at least one instance of the element, apparatus, component, means, step, etc., unless explicitly stated otherwise. The steps of any method disclosed herein do not have to be performed in the exact order described, unless explicitly stated. 
     Embodiments of the disclosed invention may enable one or more of the following:
         State explicit performance objectives and ability to determine how these should be traded off between each other.   State explicit system constraints (e.g., individual transport mission deadlines, speed and acceleration limits etc.).   Include explicit safety constraints with tunable margins.   Incorporate physical constraints and effects (e.g., the vehicle dynamics and propulsion system efficiency, including the effect increased energy consumption induced by uphill driving, or utilize the potential energy acquired by a vehicle at a hill).   Incorporate the energy state of a vehicle to automatically determine when a BEV should be charged.   Introduce explicit priorities among the vehicles.   Introduce feedback: the traffic plan can be recomputed at a high rate to handle deviations from the current plan.   Have flexibility in tuning, all the above can be changed on-the fly, leading to situations where the overall system behavior can be changed at-will (e.g., changing prioritization between throughput and energy efficiency depending on customer needs).   Perform site-wide optimization considering complete vehicle trajectories and the site topology and all conflict zones, simultaneously.       

    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Aspects and embodiments are now described, by way of example, with reference to the accompanying drawings. 
         FIG.  1 A  illustrates an intersection CZ. 
         FIG.  1 B  illustrates a merge-split CZ. 
         FIG.  1 C  illustrates a narrow road section. 
         FIG.  2 A  is an illustration of the heuristic approach for when an original problem is a non-convex problem with a disconnected solution space. 
         FIG.  2 B  shows the MIQP approximation of the original problem, which approximation allows for finding a region where the solution of the original problem could be. 
         FIG.  2 C  shows through the fixed-order NLP), the approximate optimal solution is acquired. 
         FIG.  3    shows a mock-up confined site area together with a legend of the line styles to be used in  FIGS.  4  to  12   . 
         FIG.  4    shows uncoordinated crossing for the merge-split CZ for the red vehicle (solid curve) and blue vehicle (dashed curve). 
         FIG.  5    shows coordinated crossing for the merge-split CZ for the red vehicle (solid curve) and blue vehicle (dashed curve). 
         FIG.  6    shows uncoordinated crossing of an intersection CZ for the red vehicle (solid curve) and green vehicle (dash-dotted curve). 
         FIG.  7    shows coordinated crossing of an intersection CZ for the red vehicle (solid curve) and green vehicle (dash-dotted curve). 
         FIG.  8    shows speed profiles for the coordinated vehicles (from top to bottom: red, blue, green, black). 
         FIG.  9    shows the time gap between the red and blue vehicles in the merge-split CZ for the coordinate and uncoordinated cases. 
         FIG.  10    shows uncoordinated and coordinated crossing of the intersection zone. An intersection of the depicted trajectory with the intersection zone is equivalent to both vehicles being in the CZ at the same time. 
         FIG.  11    shows uncoordinated crossing for the merge-split CZ for the red and blue vehicle. The grey box is the occupancy time in the zone for the red vehicle. While it is allowed for the blue vehicle to enter the CZ whilst the red vehicle is in it, an intersection between the vehicles translates to a collision between the vehicles 
         FIG.  12    shows coordinated crossing for the merge-split CZ for the red and blue vehicle. The grey box is the occupancy time in the zone for the red vehicle. The blue vehicle enters the CZ whilst the red vehicle is in and keeps the desired time gap. 
         FIG.  13    is a flowchart of a method for planning routs of a plurality of vehicles, according to an embodiment of the invention. 
         FIG.  14    shows an environment with multiple intersection CZs, merge-split CZs and dwelling-type CZs, in which a plurality of vehicles operates, as well as a device suitable for controlling said vehicles is deployed. 
         FIG.  15    is a sequence diagram showing an exchange of signals between a vehicle and a device configured to control a group of vehicles to which the vehicle belongs. 
         FIG.  16    shows a mock-up confined site area together with a legend of the line styles to be used in  FIGS.  17  to  19   . 
         FIG.  17    shows speed profiles for the coordinated vehicles (thicker line) and uncoordinated vehicles (thinner line, occasionally overlapping). 
         FIG.  18    shows, for the  1   st  and  3   rd  vehicles, uncoordinated and coordinated crossing of the intersection zone. An intersection of the depicted trajectory with the intersection zone is equivalent to both vehicles being in the CZ at the same time. 
         FIG.  19    shows, for the  9   th  and  10   th  vehicles, uncoordinated and coordinated crossing of a narrow road collision zone. An intersection of the depicted trajectory with the CZ is equivalent to both vehicles being in the CZ at the same time. 
     
    
    
     DETAILED DESCRIPTION 
     The aspects of the present disclosure will now be described more fully hereinafter with reference to the accompanying drawings, in which certain embodiments of the invention are shown. These aspects may, however, be embodied in many different forms and should not be construed as limiting; rather, these embodiments are provided by way of example so that this disclosure will be thorough and complete, and to fully convey the scope of all aspects of the invention to those skilled in the art. Like numbers refer to like elements throughout the description. 
     Problem Formulation 
     For the confined site coordination problem, N a  vehicles that need to be coordinated efficiently through different types of exclusion zones without collisions are considered. Assuming a confined site means that no unexpected traffic participants are present, i.e., no other vehicles than the planned ones, no pedestrians, etc. Furthermore, all coordinated vehicles are assumed to be automated, capable to communicate, their full paths and common exclusion zones are known a priori and overtakes do not occur. 
       FIG.  14    shows an environment  1400  in which a plurality of vehicles  1460  operate. The environment  1400  includes roads  1430  as well as multiple intersection CZs  1410 , merge-split CZs  1410  and dwelling-type CZs (dwelling zones)  1420 . The environment  1400  may be a confined area with no other traffic participants than the vehicles  1460 . The vehicles  1460  may be conventional vehicles, partially autonomous vehicles or fully autonomous vehicles; they may be cars, trucks, buses, construction equipment; they may be single vehicles or multi-unit vehicle combinations. Also visible in  FIG.  14    is a device  1450  suitable for controlling said vehicles  1460 . The device  1450  includes processing circuitry  1452  and is configured to communicate over a wireless interface  1451  with corresponding wireless interfaces  1461  on the vehicles  1460 . As indicated for one example vehicle  1460 , the vehicles  1460  circulating in the environment  1400  may further include vehicle-carried processing circuitry  1462 . 
     Vehicle Model 
     The motion of the vehicles along their path is described by 
         {dot over (p)}   i ( t )= v   i ( t ) 
         {dot over (x)}   i ( t )= f   i ( p   i ( t ),  x   i ( t ),  u   i ( t )) 
       0≤ h   i ( p   i ( t ),  x   i ( t ),  u   i ( t )),
 
     where p i (t) ∈   is the position, x i (t) ∈    n  the vehicle state, u i (t) ∈    m  the control input, and i is a vehicle index. The state is subdivided as x i (t)=(v i (t), z i (t)), with the speed along the path v i (t) ∈   and z i (t) ∈    n−1  collecting possible other states. The functions f i  and h i  describe respectively the dynamics and constraints that capture, e.g., actuator and speed limits. Both functions are assumed to be smooth. 
     Vehicle Model in Space 
     For confined site optimization a main interest is to optimize the trajectories of the vehicles over the whole site, i.e., over their full paths. However, the time it takes a vehicle to traverse a path is unknown. Therefore, the necessary time horizon is undefined, making the coordination problem difficult to formulate with time as the independent variable. Due to this, the problem is reformulated “in space”, using that 
     
       
         
           
             
               
                 dp 
                 i 
               
               dt 
             
             = 
             
               
                 
                   
                     v 
                     i 
                   
                   ( 
                   t 
                   ) 
                 
                 ⁢ 
                     
                 gives 
                 ⁢ 
                     
                 dt 
               
               = 
               
                 
                   dp 
                   i 
                 
                 / 
                 
                   
                     v 
                     i 
                   
                   ( 
                   t 
                   ) 
                 
               
             
           
         
       
     
     The spatial dynamics are 
     
       
         
           
             
               
                 
                   
                     
                       dt 
                       i 
                     
                     
                       dp 
                       i 
                     
                   
                   = 
                   
                     1 
                     
                       
                         v 
                         i 
                       
                       ( 
                       
                         p 
                         i 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     4 
                     
                       , 
                       TagBox[&#34;,&#34;, &#34;NumberComma&#34;, Rule[SyntaxForm, &#34;0&#34;]] 
                     
                     5 
                     
                       , 
                       TagBox[&#34;,&#34;, &#34;NumberComma&#34;, Rule[SyntaxForm, &#34;0&#34;]] 
                     
                     6 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 dx 
                 i 
               
               
                 dp 
                 i 
               
             
             = 
             
               
                 1 
                 
                   
                     v 
                     i 
                   
                   ( 
                   
                     p 
                     i 
                   
                   ) 
                 
               
               ⁢ 
               
                 
                   f 
                   i 
                 
                 ( 
                 
                   
                     p 
                     i 
                   
                   , 
                   
                     
                       x 
                       i 
                     
                     ( 
                     
                       p 
                       i 
                     
                     ) 
                   
                   , 
                   
                     
                       u 
                       i 
                     
                     ( 
                     
                       p 
                       i 
                     
                     ) 
                   
                 
                 ) 
               
             
           
         
       
       
         
           
             
               0 
               ≤ 
               
                 h 
                 ⁡ 
                 ( 
                 
                   
                     p 
                     i 
                   
                   , 
                   
                     x 
                     i 
                   
                   , 
                   
                     u 
                     i 
                   
                 
                 ) 
               
             
             , 
           
         
       
     
     where the position p i  is the independent variable. 
     Conflict Zone Modelling 
     A conflict zone (CZ) is described by the entry and exit position [p i   in , p i   out ] on the path of each vehicle. From the known positions, the entry and exit times through the CZ (t i   in , t i   out ) can be defined. In the present disclosure, two types of conflict zones are considered as depicted in  FIG.  1   . Let  ={I 1 , I 2 , . . . , I r     0   } denote a set of all intersections (total number: r 0 ) in the confined site and    r ={q r,1 , q r,2 , . . . , q r,l } denote the set of vehicles that cross intersection I r . In the intersection-like CZ, it is desired to only have one vehicle inside the CZ, i.e., not allowing the vehicle j to enter the CZ before vehicle i≠j exits the CZ, or vice-versa. The order in which the vehicles cross the intersection I r  is denoted  =(s r,1 ,s r,2 , . . . ,  ), where S r,1 , s r,2 , . . . are vehicle indices and  ={ , . . . ,  }. A sufficient condition for collision avoidance for the r-th intersection CZ can be formulated as 
     
       
         
           
             
               
                 
                   
                     
                       
                         t 
                         
                           s 
                           
                             r 
                             , 
                             i 
                           
                         
                       
                       ( 
                       
                         p 
                         
                           s 
                           
                             r 
                             , 
                             i 
                           
                         
                         out 
                       
                       ) 
                     
                     ≤ 
                     
                       
                         t 
                         
                           s 
                           
                             r 
                             , 
                             
                               i 
                               + 
                               1 
                             
                           
                         
                       
                       ( 
                       
                         p 
                         
                           s 
                           
                             r 
                             , 
                             
                               i 
                               + 
                               1 
                             
                           
                         
                         in 
                       
                       ) 
                     
                   
                   , 
                   
                     i 
                     ∈ 
                     
                       
                         [ 
                         
                           1 
                           , 
                           
                             
                               
                                 ❘ 
                                 &#34;\[LeftBracketingBar]&#34; 
                               
                               
                                 𝒬 
                                 𝓇 
                               
                               
                                 ❘ 
                                 &#34;\[RightBracketingBar]&#34; 
                               
                             
                             - 
                             1 
                           
                         
                         ] 
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     where t is determined from (4). 
     In the merge-split CZ case, let  ={M 1 , M 2 , . . . , M w     0   } denote a set of all merge-split zones (total number: w 0 ) in the site and Z w ={z w,1 , z w,2 , . . . , z w,h } denote the set of vehicles that cross the merge-split CZ M w . It is desired to have the vehicles in the zone at the same time, instead of blocking the whole zone, and thus increasing throughput efficiency. This requires having rear-end collision constraints once the vehicles have entered the CZ safely. In this case, the order in which the vehicles enter the zone is denoted as  =(s w,1 , s w,2 , . . . , s w,|Z     w     | ), and  ={ , . . . ,  }. The collision avoidance requirement for this w-th CZ leads in having to enforce the following constraints: 
     
       
         
           
             
               
                 
                   
                     
                       
                         t 
                         
                           s 
                           
                             w 
                             , 
                             i 
                           
                         
                       
                       ( 
                       
                         p 
                         
                           s 
                           
                             w 
                             , 
                             i 
                           
                         
                         in 
                       
                       ) 
                     
                     + 
                     
                       Δ 
                       ⁢ 
                       t 
                     
                   
                   ≤ 
                   
                     
                       t 
                       
                         s 
                         
                           w 
                           , 
                           
                             i 
                             + 
                             1 
                           
                         
                       
                     
                     ( 
                     
                       
                         p 
                         
                           s 
                           
                             w 
                             , 
                             
                               i 
                               + 
                               1 
                             
                           
                         
                         in 
                       
                       + 
                       c 
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   t 
                   
                     
                       s 
                       
                         w 
                         , 
                         i 
                       
                     
                     , 
                     
                       k 
                       i 
                     
                   
                 
                 + 
                 
                   Δ 
                   ⁢ 
                   t 
                 
               
               ≤ 
               
                 
                   t 
                   
                     s 
                     
                       w 
                       , 
                       
                         i 
                         + 
                         1 
                       
                     
                   
                 
                 ( 
                 
                   
                     p 
                     
                       
                         s 
                         
                           w 
                           , 
                           i 
                         
                       
                       , 
                       
                         k 
                         i 
                       
                     
                   
                   - 
                   
                     p 
                     
                       s 
                       
                         w 
                         , 
                         i 
                       
                     
                     in 
                   
                   + 
                   
                     p 
                     
                       s 
                       
                         w 
                         , 
                         
                           i 
                           + 
                           1 
                         
                       
                     
                     in 
                   
                   + 
                   c 
                 
                 ) 
               
             
             , 
           
         
       
       
         
           
             
               k 
               
                 s 
                 
                   w 
                   , 
                   i 
                 
               
               in 
             
             ≤ 
             
               k 
               i 
             
             ≤ 
             
               k 
               
                 s 
                 
                   w 
                   , 
                   i 
                 
               
               out 
             
           
         
       
       
         
           
             
               
                 
                   
                     t 
                     
                       s 
                       
                         w 
                         , 
                         i 
                       
                     
                   
                   ( 
                   
                     p 
                     
                       s 
                       
                         w 
                         , 
                         i 
                       
                     
                     out 
                   
                   ) 
                 
                 + 
                 
                   Δ 
                   ⁢ 
                   t 
                 
               
               ≤ 
               
                 
                   t 
                   
                     s 
                     
                       w 
                       , 
                       
                         i 
                         + 
                         1 
                       
                     
                   
                 
                 ( 
                 
                   
                     p 
                     
                       s 
                       
                         w 
                         , 
                         
                           i 
                           + 
                           1 
                         
                       
                     
                     out 
                   
                   + 
                   c 
                 
                 ) 
               
             
             , 
           
         
       
       
         
           
             i 
             ∈ 
             
               
                 
                   [ 
                   
                     1 
                     , 
                     
                       
                         
                           ❘ 
                           &#34;\[LeftBracketingBar]&#34; 
                         
                         
                           Z 
                           w 
                         
                         
                           ❘ 
                           &#34;\[RightBracketingBar]&#34; 
                         
                       
                       - 
                       1 
                     
                   
                   ] 
                 
               
               . 
             
           
         
       
     
     The time Δt above is added for an additional constant time gap. However, only including a time gap would lead to the problem that at low speeds the vehicles would be very close to each other, unless Δt is not very large. Large Δt, on the other hand, would lead to too large distance between the vehicles at high speeds. Therefore, a constant c is included as well. Assuming that the vehicles do not reverse, the velocity of a moving vehicle is larger than zero, meaning c&gt;0. 
     Method—First Approach 
     Optimal Control Problem 
     A direct optimal control approach is used to solve the problem. The independent variable is discretized as p i =(p i,1 , . . . p i,N     i   ), the input approximated using zero order hold such that u(p)=u i,k , p ∈[p i,k ,p i,k+1 [, and equations (4, 5) are (numerically) integrated on this grid, giving the “discretized” state transition relation 
     
       
         
           
             
               
                 [ 
                 
                   
                     
                       
                         t 
                         
                           i 
                           , 
                           
                             k 
                             + 
                             1 
                           
                         
                       
                     
                   
                   
                     
                       
                         x 
                         
                           i 
                           , 
                           
                             k 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
                 ] 
               
               = 
               
                 F 
                 ⁡ 
                 ( 
                 
                   
                     x 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     u 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     p 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     p 
                     
                       i 
                       , 
                       
                         k 
                         + 
                         1 
                       
                     
                   
                 
                 ) 
               
             
             , 
           
         
       
     
     where F denotes the integration of (4, 5) from p i,k  to p i,k+1 . Using this, the problem of finding the optimal vehicle trajectories that avoid collisions can be formalized as: 
     (Optimal Coordination Problem) Obtain the optimal state and control trajectories X*={x 1   + , . . . x N     a     + },  ={u 1   + , . . . , u N     a     + }, given the initial state X 0 ={x 1,0 , . . . x N     a     ,0 }, by solving the optimization problem 
     
       
         
           
             
               
                 
                   
                     min 
                     
                       
                         x 
                         
                           i 
                           , 
                           k 
                         
                       
                       , 
                       
                         u 
                         
                           i 
                           , 
                           k 
                         
                       
                       , 
                       
                         𝒪 
                         𝒥 
                       
                       , 
                       
                         𝒪 
                         ℳ 
                       
                     
                   
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         N 
                         a 
                       
                     
                     
                       
                         J 
                         i 
                       
                       ( 
                       
                         
                           x 
                           
                             i 
                             , 
                             k 
                           
                         
                         , 
                         
                           u 
                           
                             i 
                             , 
                             k 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     subject 
                     ⁢ 
                         
                     to 
                     ⁢ 
                         
                     initial 
                     ⁢ 
                         
                     guess 
                     ⁢ 
                         
                     
                       x 
                       
                         i 
                         , 
                         0 
                       
                     
                   
                   = 
                   
                     
                       x 
                       ^ 
                     
                     
                       i 
                       , 
                       0 
                     
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   system 
                   ⁢ 
                       
                   dynamics 
                   ⁢ 
                       
                   
                     ( 
                     9 
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   state 
                   ⁢ 
                       
                   and 
                   ⁢ 
                       
                   input 
                   ⁢ 
                       
                   constraints 
                   ⁢ 
                       
                   
                     ( 
                     6 
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     d 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     safety 
                     ⁢ 
                         
                     constraints 
                     ⁢ 
                         
                     
                       ( 
                       7 
                       ) 
                     
                   
                   , 
                   
                     
                       ( 
                       8 
                       ) 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     e 
                   
                   ) 
                 
               
             
           
         
       
     
     Note in particular that this involves finding the crossing orders  ,  , which makes the problem combinatorial and difficult to solve. 
     Problem (10) can be solved as Mixed Integer Nonlinear Program (MINLP). However, finding a solution to such problems is known to be difficult, especially when the constraints or the objective function are non-convex. Therefore, a common procedure is to apply an approach where the integer part of the solution is obtained first using a heuristic, and the continuous part of the solution thereafter obtained by solving the nonlinear program (NLP) that results from fixing the integers to the values found with the heuristic. This approach is followed in the present disclosure. The heuristic that is used is similar to that of [9] and approximates the integer solution of (10) by solving a Mixed Integer Quadratic Problem (MIQP). With the integer solution given, the state and control trajectories are obtained by solving the “fixed-order coordination” NLP, i.e., Problem (10) with fixed crossing orders  ,  . 
     Crossing Order Heuristic 
     The crossing order heuristic follows the procedure proposed in [9] where the coordination problem is approximated as time-slot scheduling MINLP in T, where T=(T i , . . . , T N     a   ) and T i  the set of entry and exit times for vehicle i for all the CZ it encounters. With this formulation, the integer part of the solution, i.e., the crossing orders, are treated in T space. The problem reads: 
     
       
         
           
             
               
                 
                   
                     min 
                     
                       T 
                       , 
                       
                         𝒪 
                         𝒥 
                       
                       , 
                       
                         𝒪 
                         ℳ 
                       
                     
                   
                       
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         N 
                         a 
                       
                     
                       
                     
                       
                         V 
                         i 
                       
                       ( 
                       
                         T 
                         i 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     11 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       
                         subject 
                         ⁢ 
                             
                         to 
                         ⁢ 
                         
                             
                              
                         
                         ⁢ 
                         
                           T 
                           i 
                         
                       
                       ∈ 
                       
                         dom 
                         ⁢ 
                            
                         
                           ( 
                           
                             V 
                             i 
                           
                           ) 
                         
                         ∀ 
                       
                     
                     = 
                     1 
                   
                   , 
                   … 
                      
                   , 
                   
                     
                       N 
                       a 
                     
                     ⁢ 
                         
                     and 
                     ⁢ 
                         
                     
                       ( 
                       7 
                       ) 
                     
                   
                   , 
                   
                     ( 
                     8 
                     ) 
                   
                   , 
                 
               
               
                 
                   ( 
                   
                     11 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
     
     where dom(V i ) is the domain of the function V i . The definition of V i (T i ) follows next. 
     The Vehicle problem: The function V i (T i ) is the optimal value function of the following parametric optimal control problem, denoted the “Vehicle Problem”: 
     
       
         
           
             
               
                 
                   
                     
                       V 
                       i 
                     
                     ( 
                     
                       T 
                       i 
                     
                     ) 
                   
                   = 
                   
                     
                       min 
                       z 
                     
                         
                     
                       
                         J 
                         i 
                       
                       ( 
                       
                         
                           x 
                           i 
                         
                         , 
                         
                           u 
                           i 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     12 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     subject 
                     ⁢ 
                         
                     to 
                     ⁢ 
                         
                     
                       ( 
                       6 
                       ) 
                     
                   
                   , 
                   
                     ( 
                     9 
                     ) 
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       b 
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     12 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     ξ 
                     i 
                     in 
                   
                   = 
                   
                     t 
                     i 
                     in 
                   
                 
               
               
                 
                   ( 
                   
                     12 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     ξ 
                     i 
                     out 
                   
                   = 
                   
                     t 
                     i 
                     out 
                   
                 
               
               
                 
                   ( 
                   
                     12 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
     with ξ i  being the parameter over which optimization is performed. Note that dom(V i ) thereby is the set of parameters T i  for which (12) has a solution. 
     An MIQP-based heuristic: The time-slot scheduling MINLP (11) is still a difficult problem, with non-convex objective function and constraints. Motivated by the availability of efficient solvers, the following MIQP approximation of (11) is considered: 
     
       
         
           
             
               
                 
                   
                     
                       
                         min 
                            
                       
                       T 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           1 
                         
                         
                           N 
                           a 
                         
                       
                         
                       
                         
                           1 
                           2 
                         
                         ⁢ 
                         
                           T 
                           i 
                           T 
                         
                         ⁢ 
                         
                           
                             ∇ 
                             2 
                           
                           
                             V 
                             i 
                           
                         
                         ⁢ 
                         
                           T 
                           i 
                         
                       
                     
                   
                   + 
                   
                     
                       
                         ( 
                         
                           
                             ∇ 
                               
                             
                               V 
                               i 
                             
                           
                           - 
                           
                             
                               
                                 ∇ 
                                 2 
                               
                               
                                 V 
                                 i 
                               
                             
                             ⁢ 
                             
                               T 
                               i 
                               
                                 [ 
                                 0 
                                 ] 
                               
                             
                           
                         
                         ) 
                       
                       T 
                     
                     ⁢ 
                     
                       T 
                       i 
                     
                   
                 
               
               
                 
                   ( 
                   
                     13 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     subject 
                     ⁢ 
                         
                     to 
                     ⁢ 
                         
                     
                       ( 
                       7 
                       ) 
                     
                   
                   , 
                   
                     ( 
                     8 
                     ) 
                   
                   , 
                 
               
               
                 
                   ( 
                   
                     13 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
     
     where the derivatives are evaluated at T i   [0] , which is obtained from the first step of Algorithm 1 below. Note that besides using the first and second order Taylor expansion of the objective function and constraints, respectively, the constraint dom(V i ) is removed. This is motivated by the conclusions in [9]. 
     The solution to the MIQP problem provides an approximately optimal time-slot schedule T and therefore an approximate crossing orders  ,  . 
     A comment on differentiation of V i (T i ): From the solution of the parametric NLP (12), the first and second order derivatives w.r.t. the parameter can be calculated using the parametric sensitivity analysis [14]. In particular, this can be done at a low additional cost when the optimal problem is solved. Note that for a solution to problem (12), feasible entry and exit times, i.e., feasible values for the parameters, for each CZ that vehicle i encounters are necessary. These times can, for example, be computed by solving an optimization problem without the safety constraints. 
     Fixed-Order NLP 
     With the found crossing order, the integer part of the solution is obtained, making (10) a NLP. Obtaining the optimal state and control trajectories is thus found through solving the fixed-order coordination problem 
     
       
         
           
             
               
                 
                   
                     min 
                     
                       
                         x 
                         
                           i 
                           , 
                           k 
                         
                       
                       , 
                       
                         u 
                         
                           i 
                           , 
                           k 
                         
                       
                     
                   
                      
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         N 
                         a 
                       
                     
                     
                       
                         J 
                         i 
                       
                       ( 
                       
                         
                           x 
                           
                             i 
                             , 
                             k 
                           
                         
                         , 
                         
                           u 
                           
                             i 
                             , 
                             k 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     14 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     subject 
                     ⁢ 
                         
                     to 
                     ⁢ 
                         
                     
                       ( 
                       
                         10 
                         ⁢ 
                         b 
                       
                       ) 
                     
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       c 
                     
                     ) 
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       d 
                     
                     ) 
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       e 
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     14 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       𝒪 
                       𝒥 
                     
                     = 
                     
                       
                         𝒪 
                         ^ 
                       
                       𝒥 
                     
                   
                   , 
                   
                     
                       𝒪 
                       ℳ 
                     
                     = 
                     
                       
                         𝒪 
                         ^ 
                       
                       ℳ 
                     
                   
                 
               
               
                 
                   ( 
                   
                     14 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
     The two stage approximation approach is summarized in Algorithm 1 and depicted in  FIG.  2   . 
     Algorithm 1 
     
         
         Input: N a ,  ,  ,   Z w , vehicle paths 
         Output: X*,  *
       1. For all i, solve NLP (10) without safety constraints (7), (8) to obtain feasible entry and exit times T i   [0] .   2. For all i, solve NLP (12) and calculate the first- and second-order sensitivities (∇V i , ∇ 2 V i ).   3. Solve MIQP (13) for the approximate crossing orders  ,      4. Solve the fixed-order NLP (14) using  ,   and obtain X*,  .   
     
       
    
     Method—Second Approach 
     A Practical Reformulation of the Collision Constraints 
     A common way to handle constraints such as (7) and (8) is to introduce auxiliary binary variables and use the “big-M” technique. For example, an equivalent representation to the constraint (7), with b s     r     ,i,i+1  ∈{0,1}, i ∈  and a sufficiently large M is 
         t   s     r,i   ( p   s     r,i     out )− t   s     r,i+1   ( p   s     r,i+1     in )≤ b   s     r     ,i,i+1   M,  
 
         t   s     r,i+1   ( p   s     r,i+1     out )− t   s     r,i   ( p   s     r,i     in )≤(1− b   s     r     ,i,i+1 ) M.  
 
     In the case where b s     r     ,i,i+1 =0, the vehicle i+1 is constrained to cross the CZ after the vehicle i, with the opposite being true if b s     r     ,i,i+1 =1. All integer variables for all CZs are connected in b ∈   2   r     0     +w     0   . This is to say, b corresponds to pairwise relative crossing orders. 
     Optimal Control Problem 
     In the Second Approach too, a direct optimal control approach is used to solve the problem. The independent variable is discretized as p i =(p i,1 , . . . , p i,N     i   ), the input approximated using zero order hold such that u(p)=u i,k , p ∈[p i,k , p i,k+1 [, and equations (4, 5) are (numerically) integrated on this grid, giving the “discretized” state transition relation 
     
       
         
           
             
               
                 [ 
                 
                   
                     
                       
                         t 
                         
                           i 
                           , 
                           
                             k 
                             + 
                             1 
                           
                         
                       
                     
                   
                   
                     
                       
                         x 
                         
                           i 
                           , 
                           
                             k 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
                 ] 
               
               = 
               
                 F 
                 ⁡ 
                 ( 
                 
                   
                     x 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     u 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     p 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     p 
                     
                       i 
                       , 
                       
                         k 
                         + 
                         1 
                       
                     
                   
                 
                 ) 
               
             
             , 
           
         
       
     
     where F denotes the integration of (4, 5) from p i,k  to p i,k+1 . As seen, the problem of finding the optimal vehicle trajectories that avoid collisions can then be formalized as the above-stated Optimal Coordination Problem. 
     The Optimal Coordination Problem (10) can be expressed as a Mixed Integer Nonlinear Program (MINLP), where the crossing order correspond to the combinatorial (“integer part”) and the state and control trajectories corresponds to the “NLP part”. In essence, the Optimal Coordination Problem can be stated as 
     
       
         
           
             
               
                 
                   
                     
                       
                         min 
                         
                           𝒲 
                           , 
                           b 
                         
                       
                     
                     
                       
                         J 
                         ⁡ 
                         ( 
                         𝒲 
                         ) 
                       
                     
                   
                   
                     
                       
                                 
                         
                           s 
                           . 
                           t 
                           . 
                             
                         
                       
                     
                     
                       
                         
                           g 
                           ⁢ 
                           
                             ( 
                             𝒲 
                             ) 
                           
                         
                         = 
                         0 
                       
                     
                   
                   
                     
                         
                     
                     
                       
                         
                           h 
                           ⁢ 
                           
                             ( 
                             𝒲 
                             ) 
                           
                         
                         ≤ 
                         0 
                       
                     
                   
                   
                     
                         
                     
                     
                       
                         
                           c 
                           ⁢ 
                           
                             ( 
                             
                               𝒲 
                               , 
                               b 
                             
                             ) 
                           
                         
                         ≤ 
                         0 
                       
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     where  ={X,  }, J( )=Σ i=1   N     a   J i (w i ), g( ), h( ) gather all equality and inequality constraints, and c( , b)=c w ( )+Cb are the integer constraints for the combinatorial part of the problem with C a matrix that captures the influence of the integer variables. 
     Since finding a solution to MINLP problems is known to be difficult, a common way of obtaining a solution is through a two-stage approximation procedure. In this disclosure, under the Second Approach, this way is followed and an alternative heuristic is proposed. 
     Crossing Order Heuristic 
     The crossing order heuristic is based on solving an MIQP that is assembled from a quadratic approximation of (15). The quadratic approximation is formed in a similar fashion as QP sub-problems are formed in SQP methods. In essence, (15) can be reformulated as: 
     
       
         
           
             
               
                 
                   
                     
                       min 
                       
                         𝒲 
                         , 
                         b 
                       
                     
                         
                     
                       
                         
                           1 
                           2 
                         
                         [ 
                         
                           
                             
                               Δ𝒲 
                             
                           
                           
                             
                               b 
                             
                           
                         
                         ] 
                       
                       T 
                     
                     ⁢ 
                     
                       
                         H 
                         ⁡ 
                         ( 
                         
                           𝒲 
                           , 
                           λ 
                           , 
                           μ 
                         
                         ) 
                       
                       [ 
                       
                         
                           
                             Δ𝒲 
                           
                         
                         
                           
                             b 
                           
                         
                       
                       ] 
                     
                   
                   + 
                   
                     
                       ∇ 
                       𝒲 
                     
                     
                       
                         
                           J 
                           ⁡ 
                           ( 
                           𝒲 
                           ) 
                         
                         T 
                       
                       [ 
                       
                         
                           
                             Δ𝒲 
                           
                         
                         
                           
                             b 
                           
                         
                       
                       ] 
                     
                   
                   + 
                   
                     J 
                     ⁡ 
                     ( 
                     
                       𝒲 
                       ** 
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     16 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   	 
                   
                     
                       
                         
                           s 
                           . 
                           t 
                           . 
                               
                         
                       
                       
                         
                           
                             
                               g 
                               ⁡ 
                               ( 
                               
                                 𝒲 
                                 ** 
                               
                               ) 
                             
                             + 
                             
                               
                                 ∇ 
                                 𝒲 
                               
                               
                                 
                                   
                                     g 
                                     ⁡ 
                                     ( 
                                     
                                       𝒲 
                                       ** 
                                     
                                     ) 
                                   
                                   T 
                                 
                                 [ 
                                 
                                   
                                     
                                       Δ𝒲 
                                     
                                   
                                   
                                     
                                       b 
                                     
                                   
                                 
                                 ] 
                               
                             
                           
                           = 
                           0 
                         
                       
                     
                     
                       
                           
                       
                       
                         
                           
                             
                               h 
                               ⁡ 
                               ( 
                               
                                 𝒲 
                                 ** 
                               
                               ) 
                             
                             + 
                             
                               
                                 ∇ 
                                 𝒲 
                               
                               
                                 
                                   
                                     h 
                                     ⁡ 
                                     ( 
                                     
                                       𝒲 
                                       ** 
                                     
                                     ) 
                                   
                                   T 
                                 
                                 [ 
                                 
                                   
                                     
                                       Δ𝒲 
                                     
                                   
                                   
                                     
                                       b 
                                     
                                   
                                 
                                 ] 
                               
                             
                           
                           ≤ 
                           0 
                         
                       
                     
                     
                       
                           
                       
                       
                         
                           
                             
                               
                                 c 
                                 𝒲 
                               
                               ( 
                               
                                 𝒲 
                                 ** 
                               
                               ) 
                             
                             + 
                             
                               
                                 ∇ 
                                 𝒲 
                               
                               
                                 
                                   
                                     
                                       c 
                                       𝒲 
                                     
                                     ( 
                                     
                                       𝒲 
                                       ** 
                                     
                                     ) 
                                   
                                   T 
                                 
                                 [ 
                                 
                                   
                                     
                                       Δ𝒲 
                                     
                                   
                                   
                                     
                                       b 
                                     
                                   
                                 
                                 ] 
                               
                             
                             + 
                             Cb 
                           
                           ≤ 
                           0 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     16 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
     
     where 
       H( , λ, μ)=blkdiag({H i } i=1   N     a   ,  0   r     0     +w     0     ,r     0     +w     0   )
 
     is a positive-definite block diagonal matrix with H i (w i , λ i , μ i )=∇ w     i     2   (w i , λ i , μ i )=∇ w     i     2 J i (w i )−∇ w     i     2 λ i   T g(w i )−∇ w     i     2 μ i   T h(w i ), where λ i , μ i  are the dual variables and zeros of appropriate size for the integer variables, and Δ = − **, with a solution guess (solution conjecture)  **. The constant term J( **) can be discarded. For the heuristic used in the Second Approach, a simplification is made by which the dual variables (λ i , μ i ) are equal to zero. This results in that the hessian only includes the expansion of the cost function, i.e., H i (w i )=∇ w     i     2 J i (w i ). The solution guess  ** can be obtained, for example, by solving the optimization problem (10) without safety constraints (10e), or by a forward simulation of the vehicles with, for example, an LQR controller. 
     The MIQP problem (16) can be compactly rewritten as 
     
       
         
           
             
               
                 
                   
                     
                       
                         min 
                         
                           𝒲 
                           , 
                           b 
                         
                       
                     
                     
                       
                         
                           
                             
                               
                                 1 
                                 2 
                               
                               [ 
                               
                                 
                                   
                                     𝒲 
                                   
                                 
                                 
                                   
                                     b 
                                   
                                 
                               
                               ] 
                             
                             T 
                           
                           ⁢ 
                           
                             H 
                             [ 
                             
                               
                                 
                                   𝒲 
                                 
                               
                               
                                 
                                   b 
                                 
                               
                             
                             ] 
                           
                         
                         + 
                         
                           
                             J 
                             T 
                           
                           [ 
                           
                             
                               
                                 𝒲 
                               
                             
                             
                               
                                 b 
                               
                             
                           
                           ] 
                         
                         + 
                         α 
                       
                     
                   
                   
                     
                       
                               
                         
                           s 
                           . 
                           t 
                           . 
                              
                         
                       
                     
                     
                       
                         
                           
                             A 
                             eq 
                           
                           [ 
                           
                             
                               
                                 𝒲 
                               
                             
                             
                               
                                 b 
                               
                             
                           
                           ] 
                         
                         = 
                         
                           b 
                           eq 
                         
                       
                     
                   
                   
                     
                         
                     
                     
                       
                         
                           
                             A 
                             ineq 
                           
                           [ 
                           
                             
                               
                                 𝒲 
                               
                             
                             
                               
                                 b 
                               
                             
                           
                           ] 
                         
                         ≤ 
                         
                           b 
                           ineq 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     where J now contains all the first order terms, α contains the linear terms and where the constraints in (16) are grouped as A eq , A ineq , b eq , b ineq  in (17). The solution to the MIQP problem provides an approximately optimal crossing orders  ,   that is obtained from the values of the integer variables b. The   component of the solution, which corresponds to the state trajectory and control signal for each of the vehicles, can be discarded from this step. 
     Fixed-Order NLP 
     With the found crossing order, the integer part of the solution is obtained, making (10) a NLP. Obtaining the optimal state and control trajectories is thus found through solving the fixed-order coordination problem 
     
       
         
           
             
               
                 
                   
                     min 
                     
                       
                         x 
                         
                           i 
                           , 
                           k 
                         
                       
                       , 
                       
                         u 
                         
                           i 
                           , 
                           k 
                         
                       
                     
                   
                      
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       
                         N 
                         a 
                       
                     
                     
                       
                         J 
                         i 
                       
                       ( 
                       
                         
                           x 
                           
                             i 
                             , 
                             k 
                           
                         
                         , 
                         
                           u 
                           
                             i 
                             , 
                             k 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     18 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     subject 
                     ⁢ 
                         
                     to 
                     ⁢ 
                         
                     
                       ( 
                       
                         10 
                         ⁢ 
                         b 
                       
                       ) 
                     
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       c 
                     
                     ) 
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       d 
                     
                     ) 
                   
                   , 
                   
                     ( 
                     
                       10 
                       ⁢ 
                       e 
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     18 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       𝒪 
                       𝒥 
                     
                     = 
                     
                       
                         𝒪 
                         ^ 
                       
                       𝒥 
                     
                   
                   , 
                   
                     
                       𝒪 
                       ℳ 
                     
                     = 
                     
                       
                         𝒪 
                         ^ 
                       
                       ℳ 
                     
                   
                 
               
               
                 
                   ( 
                   
                     18 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
     The two-stage approximation approach is summarized in Algorithm 2 and illustrated in  FIG.  2   . 
     Algorithm 2 
     
         
         Input: N a ,      Z w , vehicle paths 
         Output: X*,  
       1. For all i, obtain a solution guess w i ** by, e.g., solving NLP (10) without safety constraints (7), (8).   2. Calculate and form the approximation terms H, J, α.   3. Solve MIQP (17) for the approximate crossing orders  ,      4. Solve the fixed-order NLP (18) using  ,   and obtain X*,  .   
     
       
    
     Simulation Results 
     In this section, a simulation example showing the operation and performance of the coordination algorithm will be presented. 
     Simulation Setup 
     The vehicles are modelled as a triple integrator  =u, whereby the spatial model is: 
     
       
         
           
             
               
                 [ 
                 
                   
                     
                       
                         
                           dt 
                           i 
                         
                         
                           dp 
                           i 
                         
                       
                     
                   
                   
                     
                       
                         
                           dv 
                           i 
                         
                         
                           dp 
                           i 
                         
                       
                     
                   
                   
                     
                       
                         
                           da 
                           i 
                         
                         
                           dp 
                           i 
                         
                       
                     
                   
                 
                 ] 
               
               = 
               
                 
                   
                     [ 
                     
                       
                         
                           
                             1 
                             / 
                             
                               v 
                               i 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           
                             1 
                             / 
                             
                               v 
                               i 
                             
                           
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                     
                     ] 
                   
                   ⁢ 
                   
                     x 
                     i 
                   
                 
                 + 
                 
                   
                     [ 
                     
                       
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                       
                       
                         
                           
                             1 
                             / 
                             
                               v 
                               i 
                             
                           
                         
                       
                     
                     ] 
                   
                   ⁢ 
                   
                     u 
                     i 
                   
                 
               
             
             , 
           
         
       
     
     where a i  is the acceleration and u i  is the jerk. The dynamics are discretized using multiple shooting and an Explicit Runge-Kutta-4 (ERK4). 
     The constraints are chosen as bounds on the speed and longitudinal acceleration ( v   i ≤v i,k ≤ v   i , a i,k ≤ā i,lon ) in order to obey speed limits and physical constraints. The vehicles are expected to manoeuvre on curved roads, thus, it is needed to consider the lateral forces they experience. As the one-dimensional model that is used in this disclosure does not account for lateral motion, the following constraint is enforced, as similarly proposed in [8]: 
     
       
         
           
             
               
                 
                   
                     ( 
                     
                       
                         a 
                         
                           i 
                           , 
                           k 
                         
                       
                       
                         
                           a 
                           _ 
                         
                         
                           i 
                           , 
                           lon 
                         
                       
                     
                     ) 
                   
                   2 
                 
                 + 
                 
                   
                     ( 
                     
                       
                         
                           
                             κ 
                             i 
                           
                           ( 
                           
                             s 
                             
                               i 
                               , 
                               k 
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           v 
                           
                             i 
                             , 
                             k 
                           
                         
                       
                       
                         
                           a 
                           _ 
                         
                         
                           i 
                           , 
                           lat 
                         
                       
                     
                     ) 
                   
                   2 
                 
               
               ≤ 
               1 
             
             , 
           
         
       
     
     where ā i,lat  is the lateral acceleration limit and κ i (p i,k ) is the road curvature, that is assumed to be available at every point along the path. 
     In this disclosure, under the First Approach, it is of interest to minimize the change of acceleration and jerk in order for the vehicles to operate in a comfortable way. Minimizing the acceleration and jerk would also lead to a reduction in fuel/energy consumption. This requirement leads to the following objective function for each vehicle: 
     
       
         
           
             
               
                 
                   J 
                   i 
                 
                 ( 
                 
                   
                     x 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     u 
                     
                       i 
                       , 
                       k 
                     
                   
                 
                 ) 
               
               = 
               
                 
                   
                     ∑ 
                     
                       k 
                       = 
                       1 
                     
                     N 
                   
                     
                   
                     ( 
                     
                       
                         ( 
                         
                           
                             
                               P 
                               i 
                             
                             · 
                             
                               a 
                               
                                 i 
                                 , 
                                 k 
                               
                               2 
                             
                           
                           + 
                           
                             
                               Q 
                               i 
                             
                             · 
                             
                               u 
                               
                                 i 
                                 , 
                                 k 
                               
                               2 
                             
                           
                         
                         ) 
                       
                       · 
                       
                         
                           Δ 
                           ⁢ 
                           
                             p 
                             
                               i 
                               , 
                               k 
                             
                           
                         
                         
                           v 
                           
                             i 
                             , 
                             k 
                           
                         
                       
                     
                     ) 
                   
                 
                 + 
                 
                   
                     R 
                     i 
                   
                   · 
                   
                     t 
                     
                       f 
                       , 
                       i 
                     
                   
                 
               
             
             , 
           
         
       
     
     with Δp i,k =p i,k+1 −p i,k , P i , Q i , R i  the appropriate weights and t f,i  being the final time. 
     For evaluation purposes, a mock-up confined area with four vehicles as depicted in  FIG.  3    will be investigated. There are in total two merge-split CZs and five intersection CZs. Every vehicle starts from an initial velocity of 50 [km/h] and the chosen bounds are v i =3.6 [km/h],  v   i =90 [km/h], ā i,lon =4 [m/s 2 ], ā i,lat 2 [m/s 2 ]. The weight coefficients for the optimal problem are: P i =1, Q i =1, R i =10. The total shooting points for this scenario is N=100 and is same for all vehicles. In the intersection-like scenario the CZ is created with 5 meter margin ahead of and behind the collision point, whereas in the merge-split scenario the margin is 15 meters for both the entry and exit collision point. For the merge-split CZ it is also desired to keep at least 0.5 s margin between the vehicles. 
     The Casadi [18] toolkit with the IPOPT [19] solver will be utilized for solving the optimization problems (12, 14) and Gurobi for solving the MIQP problem (13). 
     Under the Second Approach, the following trade-off between minimization of the total travel time and the squares of the longitudinal acceleration and longitudinal jerk is considered: 
     
       
         
           
             
               
                 J 
                 ⁡ 
                 ( 
                 
                   
                     
                       x 
                       i 
                     
                     ( 
                     
                       p 
                       i 
                     
                     ) 
                   
                   , 
                   
                     
                       u 
                       i 
                     
                     ( 
                     
                       p 
                       i 
                     
                     ) 
                   
                 
                 ) 
               
               = 
               
                 
                   
                     ∫ 
                     
                       p 
                       
                         i 
                         , 
                         0 
                       
                     
                     
                       p 
                       
                         i 
                         , 
                         
                           N 
                           i 
                         
                       
                     
                   
                   
                     
                       ( 
                       
                         
                           
                             p 
                             i 
                           
                           ⁢ 
                           
                             
                               
                                 a 
                                 i 
                               
                               ( 
                               
                                 p 
                                 i 
                               
                               ) 
                             
                             2 
                           
                         
                         + 
                         
                           
                             Q 
                             i 
                           
                           ⁢ 
                           
                             
                               
                                 j 
                                 i 
                               
                               ( 
                               
                                 p 
                                 i 
                               
                               ) 
                             
                             2 
                           
                         
                       
                       ) 
                     
                     ⁢ 
                     
                       1 
                       
                         
                           v 
                           i 
                         
                         ( 
                         
                           p 
                           i 
                         
                         ) 
                       
                     
                     ⁢ 
                     
                       dp 
                       i 
                     
                   
                 
                 + 
                 
                   
                     R 
                     i 
                   
                   ⁢ 
                   
                     
                       t 
                       i 
                     
                     ( 
                     
                       p 
                       
                         i 
                         , 
                         
                           N 
                           i 
                         
                           
                       
                     
                     ) 
                   
                 
               
             
             , 
           
         
       
     
     where P i , Q i , R i  are the appropriate weights. This objective is integrated with Forward Euler, leading to the “discretized” objective expression 
     
       
         
           
             
               
                 
                   J 
                   i 
                 
                 ( 
                 
                   
                     x 
                     
                       i 
                       , 
                       k 
                     
                   
                   , 
                   
                     u 
                     
                       i 
                       , 
                       k 
                     
                   
                 
                 ) 
               
               = 
               
                 
                   
                     ∑ 
                     
                       k 
                       = 
                       1 
                     
                     N 
                   
                     
                   
                     ( 
                     
                       
                         ( 
                         
                           
                             
                               P 
                               i 
                             
                             ⁢ 
                             
                               a 
                               
                                 i 
                                 , 
                                 k 
                               
                               2 
                             
                           
                           + 
                           
                             
                               Q 
                               i 
                             
                             ⁢ 
                             
                               j 
                               
                                 i 
                                 , 
                                 k 
                               
                               2 
                             
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           Δ 
                           ⁢ 
                           
                             p 
                             
                               i 
                               , 
                               k 
                             
                           
                         
                         
                           v 
                           
                             i 
                             , 
                             k 
                           
                         
                       
                     
                     ) 
                   
                 
                 + 
                 
                   
                     R 
                     i 
                   
                   ⁢ 
                   
                     t 
                     
                       i 
                       , 
                       N 
                     
                   
                 
               
             
             , 
           
         
       
     
     with Δp i,k =p i,k+1 −p i,k .For evaluation purposes, the mock-up confined area depicted in  FIG.  16    will be investigated. There are in total ten vehicles and two merge-split CZs, two narrow road CZs and sixteen intersection CZs. Every vehicle starts from an initial velocity of 50 km/h, and furthermore vehicle 5 and vehicle 10 start from a nonzero initial time for the purpose of the interacting vehicles to occupy the CZ at the same time. The respective initial times for those vehicles are t 0,5 =56 and t 0,10 =11.3 seconds. The chosen velocity and acceleration bounds are  v   i =3.6 [km/h],  v   i =90 [km/h], ā i,lon =4 [m/s 2 ], ā i,lat =2 [m/s 2 ]. The weight coefficients for the optimal problem are: P i =1, Q i =1, R i =10. The total shooting points for this scenario is N=100 and is same for all vehicles. In the intersection-like scenario the CZ is created with 5 meter margin ahead of and behind the collision point, where as in the merge-split and narrow road scenarios the margin is 15 meters for both the entry and exit collision point. For the merge-split CZ it is also desired to keep at least 0.5 s margin between the vehicles. The CasADi toolkit [18] with the IPOPT solver [19] is used for solving the optimization problem (14) and Gurobi for solving the MIQP problem (12). 
     Discussion of Results 
     In the following, the results for this simulation scenario are presented and, in particular, the uncoordinated and coordinated results for one merge-split zone and one intersection zone, both circled over in  FIG.  3   . The uncoordinated results are obtained as the individual optimum of each vehicle, i.e., the speed profiles are obtained by solving the optimization problem without any safety constraints. 
       FIG.  4   , in which path coordinate is plotted as a function of time, shows uncoordinated crossing for the merge-split CZ for the red vehicle (solid curve) and blue vehicle (dashed curve) introduced in  FIG.  3   . Similarly,  FIG.  5    shows coordinated crossing for the merge-split CZ for the red vehicle (solid curve) and blue vehicle (dashed curve).  FIG.  6    shows uncoordinated crossing of an intersection CZ for the red vehicle (solid curve) and green vehicle (dash-dotted curve).  FIG.  7    shows coordinated crossing of an intersection CZ for the red vehicle (solid curve) and green vehicle (dash-dotted curve). 
       FIGS.  11  and  12    illustrate the uncoordinated and coordinated crossing of the merge-split CZs, respectively. In the figures, the path trajectories of the involved vehicles are depicted. The path trajectories have an offset of their entry position to the merge-split CZ, meaning that, for both vehicles, position of zero indicates the entry position to the CZ. The grey rectangle depicts the CZ, with entry and exit times for the red vehicle, since that is the first vehicle that enters the zone in both cases. As mentioned, in this type of CZ it is desired to have multiple vehicles inside the zone whilst keeping a minimum gap between them. This means that while it is allowed for both trajectories to be inside the zone, it is not allowed for the trajectories to intersect whilst in the zone. In the uncoordinated case, it can be noticed that the blue vehicle enters the CZ later than the red vehicle, however, intersects the red vehicle and exits the zone first. The interpretation of this behavior is that inside the CZ the blue vehicle “overruns” the red vehicle, meaning that a collision occurs. In the coordinated case, the vehicles are capable of being inside the CZ at the same time and keeping at least a minimum designated gap between each other. Note that an intersection of the trajectories outside the CZ is not relevant as they are no longer on a shared road. By allowing the vehicles to be both in the CZ, and not blocking the whole zone for one vehicle, the throughput is increased meaning that the vehicles could get to their end destinations sooner. This could translate to increased efficiency of the site. To supplement the figures,  FIG.  9    depicts the time gap between the vehicles while they are in a merge-split CZ. In the figure, the moment where the blue vehicle “overruns” the red vehicle in the uncoordinated case can be noticed, as well as the fact, that the desired minimum gap is kept for the coordinated case. 
     In the scenarios where two (or more) vehicles intersect, it is desired to only have one vehicle in the zone. The reasoning for this decision is that the CZ in this case occupies a relatively small patch of road and blocking off the whole zone for one vehicle will not result in major loss of throughput efficiency.  FIG.  10    (First Approach), just like  FIG.  18    (Second Approach), depicts the uncoordinated and coordinated vehicle trajectories along with the intersection zone. For this dependency, a collision is defined if the trajectory intersects the intersection zone, which is equivalent to the two vehicles being in the zone at the same time. It is noticeable that in the uncoordinated case the intersection zone is intersected, whilst the coordinated vehicles satisfy the collision constraints as they successfully avoid being in the CZ at the same time. 
     For the narrow road collision zones, the collision avoidance is defined in the same way as in the intersection zone since the road conditions allow for only one vehicle to occupy the zone at a time. In  FIG.  19    (Second Approach), the uncoordinated vehicle trajectory intersects the collision zone which, as in the intersection zone scenario, results in both vehicles being inside the zone at the same time and is not desired. The coordinated case manages to avoid concurrent occupancy of the zone. 
       FIG.  8    (First Approach) depicts the vehicles&#39; speed profiles throughout the site. As the objective function is to minimize the change of acceleration and jerk, it is noticeable that the speed profiles are smooth throughout the whole site and would not be challenging to follow in an application case. In the speed plots, the solid black lines represent the curvature constraints. It is worth noting that in the uncoordinated case for the whole site, two more collisions occur: one in the other merge-split CZ between the blue and green vehicle and another in the intersection CZ of the blue and black vehicle. When the vehicles are coordinated, these collisions are avoided. For sake of brevity, these results will not be depicted here, as the critical behavior is similar to the above-mentioned collision zones. 
       FIG.  17    (Second Approach) depicts the vehicles&#39; speed profiles throughout the site for the uncoordinated and coordinated case. As the algorithm is aware of all CZs the vehicles encounter, it is able with small changes of the speed profile to avoid collisions that can happen further away. With the objective function aiming to minimize the acceleration and jerk, it is noticeable that the speed profiles are smooth throughout the whole site and that they would not be challenging to follow in an application case. In the speed plots, the solid black lines touching the speed curves from above represent the curvature constraints. Under the Second Approach too, more collisions occur in the uncoordinated case for the complete site, but when the vehicles are coordinated these collisions are avoided. 
     The simulation scenario is implemented in MATLAB on a 2.90 GHz Intel Xeon computer with 32 GB of RAM. Under the First Approach, the total computational time for solving the example is 40.6313 seconds. The main computational load comes from the Casadi toolkit and setup. To be precise, the total time all solvers take is 1.815 seconds, of which the MIQP requires 0.08 seconds and the main Optimal Control Problem (OCP) 0.709 seconds and the remainder is the required time to solve the individual vehicle problem. The computational effort besides the solver time is the time Casadi takes to set up the environment. The total solver time is, however, still fairly high especially for real-time vehicle application and needs to be further improved. Under the Second Approach, the total time all solvers take is 1.973 seconds, of which the MIQP requires 0.140 seconds and the main Optimal Control Problem (OCP) 1.833 seconds and the remainder is the required time to solve the individual vehicle problem. 
     Conclusions and Future Theoretical Work 
     The present disclosure has presented two optimal-control based approaches for coordination of automated vehicles for sites where multiple collision zones of different types could occur. The approaches optimize the trajectories of the vehicles for their whole path length, taking all collision zones into account. To compute the crossing order in the collision zones, both approaches rely on an MIQP-based heuristic. Simulation results demonstrate the capability to coordinate the vehicles without collisions. 
     Implementations and Performance 
       FIG.  13    is a flowchart of a method  1300  for controlling a plurality of vehicles  1460  (see  FIG.  14   ) operating in a common environment  1400  which includes at least one CZ, wherein movements of an i th  vehicle is controllable by a control signal u i . The execution of the method  1300  may be coordinated by a central processing resource being either fixedly installed or vehicle-carried. In the following description, by way of example, it will be assumed that the execution of the method  1300  is coordinated by processing circuitry  1452  fixedly installed in the device  1450 . The device  1450  can delegate the execution of some steps of the method  1300 , or parts thereof, to vehicle-carried processing resources (e.g., processing circuitry  1462 ) by transmitting a corresponding delegation signal D to the vehicles  1460  over the wireless interface  1451  and receiving a processing result message E in return. 
     In a first step  1310  of the method  1300 , current positions of all vehicles are sensed. The vehicle positions may be self-reported positions represented by a signal P (see  FIG.  15   ) transmitted from each vehicle  1460  to the device  1450 . The transmission of the signal P may be event-triggered, periodic or triggered by a request (not shown) transmitted from the device  1450  to the vehicles  1460 . 
     In a second step  1320 , a predefined objective function J is obtained. The objective function may be related to a dynamic vehicle model representing evolution with respect to path length (p i ), in which time (t i (p i )) and path speed (v i (p i ) are state variables. 
     In a third step  1330 , a first optimization problem for a first objective function V i (T i ) derived from the predefined objective function is solved, to obtain a vehicle crossing order at each CZ, wherein the first optimization problem is subject to safety constraints. The safety constraints for an intersection zone  1410  or a dwelling zone  1420  may include a mutual exclusion requirement. The safety constraints for a merge-split zone  1410  may include a minimum longitudinal spacing requirement, Δt or Δt+c. 
     In a fourth step  1340  of the method  1300 , an OCP for the predefined objective function is solved subject to the vehicle crossing order at the CZs (obtained from the third step  1330 ) and subject to the safety constraints, to obtain a control signal u i  for each of the vehicles. It is understood that OCP includes a dynamic vehicle model representing evolution with respect to path length p i , in which time t i (p i ) and path speed v i (p i ) are state variables. 
     In a fifth step  1350 , the control signals u i  are fed to the vehicles (signal U in  FIG.  15   ) for execution. 
     Under the First Approach, the CZ entry and CZ exit times may be decision variables of the first optimization problem. Further, the first objective function V i (T i ) can be a parametric local optimum of the predefined objective function J, wherein the parametric local optimum is parametrized by tentative CZ entry and CZ exit times and is independent of the safety constraints. In this case, the tentative CZ entry and CZ exit times are decision variables in the first optimization problem. 
     Optionally, the first objective function is approximated  1334  by a polynomial, such as a quadratic polynomial. In particular, the polynomial may be a Taylor expansion around an approximate set of CZ entry and CZ exit times; see equation 13a. The approximate set of CZ entry and CZ exit times may be obtained  1332  by a preceding step of optimizing the predefined objective function independently of the safety constraints. The optimization problem expressed in terms of the polynomial may be solved as a mixed-integer quadratic program (MIQP). 
     Under the Second Approach, the CZ entry and CZ exit times may be decision variables of the first optimization problem that is solved in the third step  1330 . Further, the first objective function can be a quadratic approximation of the predefined objective function J. 
     Optionally, the first objective function is transformed  1334  by a variable substitution (     − **) for those decision variables that represent state trajectories and control signals for the vehicle. The variable substitution corresponds to subtracting a solution conjecture ( **). The solution conjecture ( **) can be obtained  1332  by optimizing the predefined objective function independently of the safety constraints. The first optimization problem with the quadratic approximation may be solved as a mixed-integer quadratic program (MIQP). 
     In one embodiment, independently of whether the First or Second Approach is taken, a subset comprising only the obtaining  1320 , solving  1330 , and OCP solving  1340  steps of the method  1300  can be implemented in isolation, for the purpose of planning routes for said plurality of vehicles  1460 . 
     The invention can be embodied as a software component that returns motion plans (e.g. speed profiles) for all vehicles from their current position to the end of their mission. The core of this software component is the generation of the motion plans, done such that no collisions occur and such that the physics of reality is considered. As seen above, the motion plans are generated by (approximately) solving an optimization problem. This optimization problem is formulated and solved based on an optimal control formulation of the site-coordination problem, using procedures from what is known as Numerical Optimal Control (NOC). NOC is a field where mathematical optimization is used to compute control commands for a dynamical system (e.g. a mechanical system such as trucks/machines) such that its evolution optimizes a performance criterion while satisfying certain constraints. 
     The solution to the problem yields the motion profiles for each involved vehicle, from its current position on the site to the terminal point of its transport mission, which
         are consistent with what the vehicles can do given the available knowledge of the environment and the limitations imposed by the designer (satisfy the vehicle dynamics and the vehicle constraints),   are free of collisions and satisfies MUTEX requirements,   maximizes performance in terms of the metric specified under the two conditions above.       

     As regards computational tractability, the following can be noted. In contrast to OCPs found in most other fields, the site coordination problem is combinatorial, as the order in which the vehicles on the site utilize the CZs (or MUTEX zones) must be determined to produce a solution. This means that the size of the space of possible solutions grows exponentially and that its solution therefore is computationally intractable for anything but small problems, in the general case. For illustration, consider a k-way crossing with n j  vehicles approaching along a j-th lane. The total number of crossing combinations are 
     
       
         
           
             
               
                 ( 
                 
                   
                     n 
                     1 
                   
                   + 
                   
                     n 
                     2 
                   
                   + 
                   … 
                   + 
                   
                     n 
                     k 
                   
                 
                 ) 
               
               ! 
             
             
               
                 
                   n 
                   1 
                 
                 ! 
               
               ⁢ 
               
                 
                   n 
                   2 
                 
                 ! 
               
               ⁢ 
                  
               … 
               ⁢ 
                  
               
                 
                   n 
                   k 
                 
                 ! 
               
             
           
         
       
     
     For a total of 8 vehicles, all approaching a four-way intersection with two from each direction, there would be 2520 distinct solutions, differing in which order the vehicles use each CZ (the crossing order). As the number of CZs and vehicles increase, the number of solutions grows rapidly. 
     Algorithms to handle this type of problems typically operate by systematically ruling out large parts of the possible crossing orders. This process involves solving so-called “relaxations” of the original problem. These relaxations are not combinatorial in themselves but otherwise inherits their complexity-class from the non-combinatorial part of the original problem. In the case of site-coordination problems, these are non-convex, nonlinear programs (NCNLPs), which typically are non-trivial to solve, even though so called local solutions often can be obtained. This means that, for a typical problem, many relaxations must be obtained to find a solution, and this number grows with the number of possible crossing orders in the original problem&#39;s solution-space. Hence, with the site coordination problems considered, most real problems become practically unsolvable in the context within they must be used, as the plan will be irrelevant once it is obtained. 
     This is one of the reasons for which the invention includes an approximating separation, performed in two steps:
         (1) The construction of the “combinatorial part” of the solution through a “site-aware” heuristic, i.e., for each CZ the determination of which vehicle goes first, second, third etc., within the set of contending vehicles, taking into account the vehicle limitations, the site topology etc. approximately; and   (2) the solution of a site-wide optimal control problem where the MUTEX order is fixed. This is a standard OCP and computationally tractable optimization problem, solvable with commercial tools within the time limits admitted by the usage context, whose solution gives the optimal motion profiles for all vehicles given the current vehicle state and the heuristically obtained MUTEX order.
 
Here, the site-aware heuristics may include two stages:
   (automatically) simplified representations of the performance objectives, physical constraints and vehicle constraints; and   a formulation of a simplified representation of the problem, using these simplifications, as a Mixed-Integer Quadratic Program (MIQPs).
 
Indeed, while the size of the combinatorial solution space remains the same (e.g., the number of distinct ways the 8 vehicles could cross the intersection in the example above), there are several well developed commercial tools to handle MIQPs efficiently. This allows the invention to handle significantly larger (and therefore relevant) problems within the timespans admitted by the intended usage context. The excellent scalability is one of the prominent features of the invention, allowing it to handle coordination of a desired number of vehicles in real time or near real time conditions.
       

     The principal steps of the invention can be summarized as follows:
         1. The central site planner obtains the current position and velocity of each automated vehicle on the site via wireless communication.   2. The central site planner finds a crossing order by a site-aware heuristic based on the solution of a Mixed-Integer Quadratic Program (MIQP).
           a. If the MIQP solution fails, the planner sets the crossing order of the previous cycle to be the current.   
           3. The central site planner solves an OCP for all vehicles where the vehicles are positioned according to the measurements obtained in step  1  and the crossing order is that produced in step 2.
           a. If the OCP solution fails, the planner sets a system fault and produce plans for all vehicles to stop immediately.   
           4. The central planner distributes the plan produced in step 3 to all vehicles.       

     Already on this general level, the invention differs from [6] notably in that the non-approximated problem, and therefore the OCP solved in step 3, is formulated with path length (distance along the path) as the independent variable, rather than time. This enables planning over the entire site (the full plan) at all moments without significant effort increase and avoids certain numerical difficulties, both issues preset in [6]. If one were to follow the lines of [6], severe penalties arise in one or both of computational resource requirements and problem formulation complexity. Notionally ‘using space rather than time’ as the independent variable is therefore an important step towards being able to treat site planning problems with optimization, which is not trivially realized. 
     Under step 2, the site-aware heuristic may operate as follows:
         For each automated vehicle, compute a quadratic approximation of that vehicle&#39;s optimal performance metric as a function of the previous CZ timeslot schedule (initialized by e.g. first-come first-served), this amounts to so solving one small OCP (‘greedy’ problem) for each vehicle and applying tools from parametric sensitivity analysis.   Attempt to solve, using the central site planner, an MIQP constructed with the quadratic approximations from a). This gives a tentative CZ schedule.   From the CZ schedule, extract the crossing order.
 
The MIQP is thus constructed on a high-level abstraction of the original OCP (the problem of finding the optimal time-slot schedule) and does not deal with the positions or velocities of the vehicles directly. The resulting MIQP has small relaxations but necessitates the solution of one small OCP per vehicle.
       

     With regard to step 4, it is noted that the plan produced by the central site planner need not be followed precisely by the vehicles on the site. If this is the case, replanning is necessary. Embodiments of the invention may therefore include feedback by periodically solving the site wide OCP based on fresh measurements of the system state. Also, the crossing order could be updated using a local search heuristic, or re-computed from scratch with either of the two heuristics (approaches). 
     The above-described solution can be parallelized in one or more of the following ways:
         Parallel computation of each vehicle&#39;s cost and quadratic approximation. Central solution of MIQP and subsequent solution of site wide OCP.   On-board computation (i.e., using processing resources carried by the vehicle) of each vehicle&#39;s OCP, data is sent to a server which solves the MIQP site wide OCP.   On-board computation of each vehicle&#39;s OCP, data is sent to server which computes the MUTEX utilization order, server and vehicles solve the problem cooperatively, distributed computations.   On-board computation of each vehicle&#39;s OCP, data is sent vehicle to vehicle, MUTEX utilization order is negotiated vehicle to vehicle, decentralized computation.   Site-wide OCP: The solution can be parallelized by use of tailored linear-algebra which leverages the numerical structure of the problem. For instance, parts of the computations needed to solve a common problem can be delegated by the central site planner to processing resources carried by the vehicles. This may relate to computationally offloading the central site planner without necessarily ensuring that the delegation causes each vehicle to carry out ‘its own’ computations.       

     The aspects of the present disclosure have mainly been described above with reference to a few embodiments and examples thereof. However, as is readily appreciated by a person skilled in the art, other embodiments than the ones disclosed above are equally possible within the scope of the invention, as defined by the appended patent claims.