2.1 Definitions and context

The description of the concept of free flight [Reference Schultz, Shaner and Zhao15] already refers to aircraft conflict management as a “stochastic control problem” due to the several sources of uncertainty. Airspace is broken down into three regions with increasing levels of restriction as aircraft fly from low-density areas to the vicinity of terminals. Together with the description of the concept a principal analysis of functions, communication and data to be shared is provided. A “mixed focus” configuration is outlined defining the roles of aircraft crews, air traffic controllers and airline operators. Several conflict determination and resolution methods are presented categorised according to aiding logic type, control variables and whether they are single or continuous solutions. Three variables are identified as possible ways to avoid conflicts: speed, altitude and heading. Although the concept of free flight is clearly outlined, no safety or performance analysis is presented.

The feasibility of the free flight concept of operations for urban air mobility has already been investigated [Reference Yang and Wei32]. The aim of the framework elaborated is to enable safe and efficient flights and a high density of operations, building on the already established benefits of free flight. The problem of conflict resolution between aircraft while guiding them to their destination was formulated as a multi-agent Markov decision process and solved by the Monte Carlo tree search algorithm. The effect of creating airspace sectors and flights with higher priority was also investigated using numerical simulations. The scenario (2D level-flying aircraft in a high-density environment) and the concept (attempting to limit conflict probability by introducing coordination) are similar to those of this paper. An important difference however is that in the study by Yang and Wei [Reference Yang and Wei32] conflicts are resolved by heading changes, potentially leading to unsolvable situations. Also, the effect of the method applied can only be seen by a posteriori analysing the simulation data, thus, no assurances can be made before deploying the strategy.

With the transition to free flight, the traditionally centralised air traffic control becomes a distributed one [Reference Blom, Bakker, Obbink and Klompstra33]. To analyse the emergent behaviour of an integration of systems with the novel behaviour Petri net modelling and Monte-Carlo simulations are used. The overall goal is to understand how to design the free flight framework for optimal safety and capacity by understanding the process of multiple conflict resolution and providing an unbiased estimator for aircraft collision risk considered to be rare events. The conflict resolution itself is based on the well-known Modified Potential Field method. Even though the density of aircraft simulated is higher than usually encountered in general aviation, it still does not compare to projected PAV (Personal Aerial Vehicle) density.

2.2 Sense and avoid methods

The robustness of a collaborative sense and avoid algorithm in the terminal control area is investigated by Allignol et al. [Reference Allignol, Barnier, Durand, Manfredi and Blond34], using an adaptation of the ORCA (Optimal Reciprocal Collision Avoidance) algorithm [Reference van den Berg, Guy, Lin and Manocha35], changing the heading of UASs (Unmanned Aerial System) to avoid conflicts with aircraft on a fixed trajectory. Conflict and separation are defined based on the distance to other aircraft for a given lookahead time. Probability theory and random variables are used in this context to assess the effect of navigation accuracy on the performance of the algorithm. This analysis shows how important it is to take the effect of noise and uncertainty into account. However, this analysis only considers normal distribution for position and velocity, since the distribution is a result of inaccuracies, not intentional stochastic design.

A probabilistic conflict detection algorithm is described by Hao et al. [Reference Hao, Zhang, Cheng, Liu and Xing36]. The underlying assumption is that waypoints and arrival times are constrained, therefore the reachable set of a given aircraft and the probability distribution of future positions changes continuously with time. Using a truncated spatio-temporal Brownian bridge motion model the conflict probability of multiple aircraft can be calculated in space-time. This method is only applicable to trajectory-based operations, not in a free flight scenario.

A probabilistic collision avoidance algorithm for unmanned aircraft can be formulated based on Markov decision processes (MDPs) and partially observable Markov decision processes (POMDPs) [Reference Temizer, Kochenderfer, Kaelbling, Lozano-Perez and Kuchar37]. The key advantage of this method is that it can accommodate a wide range of sensor types. After a state-space reduction, POMDP solvers can be used to produce optimal policies which can be implemented in real-time. Even though the solution is optimal for different sensor configurations, this algorithm does not consider the higher level, ensemble related performance of the system.

For traffic management and coordination, a novel air traffic flow model is presented by Alam et al. [Reference Alam, Chaimatanan, Delahaye and Feron38] based on functional airspace blocks. Distributed 4D trajectory planning is achieved based on input from airlines, airports and Air Navigation Service Providers (ANSPs) and centralised planning by the central flow management unit (CFMU). The overall goal of the algorithm is to minimise the interaction between trajectories by applying route/departure-time allocation techniques. Departure-time techniques adjust planned departure time within bounds, while route modification consists of waypoint repositioning. This enables a MIP (Mixed-Integer Optimisation Problem) formulation. Interactions are identified using a four-dimensional grid structure and iteratively minimised using a hybrid method with local search (LS) in the inner-loop and simulated annealing (SA) in the outer. The presented distributed method, however, is only applicable for pre-planning purposes.

Another distributed, 4D path planning algorithm is introduced by Egorov et al. [Reference Egorov, Kuroda and Sachs27] for an unstructured airspace, emphasizing the importance of encounter aware planning. The necessary characteristics of a traffic management system are described as robust, reliable and scalable. During the off-line, pre-planning phase, trajectories are optimised considering the probability distribution of the types of conflicts according to the number of vehicles affected. As a metric encounter shift is introduced. The efficiency of the method is proven by simulations for a multi-modal traffic region. The method, however, is only applicable off-line for strategic deconfliction, and thus, it cannot adjust for actual operational deviations.

2.3 Airspace and system-level analysis

The stochastic simulator presented by Lecchini et al. [Reference Lecchini Visitini, Glover, Lygeros and Maciejowski39] is able to take flight related uncertainties into account. In this stochastic environment conflict resolution is formulated as a constrained optimisation problem, which is solved by the Markov Chain Monte Carlo method. A drawback of the applied algorithm is that it can only be applied off-line.

Having a prescribed structure of the airspace can provide an a priori separation procedure for aircraft, and thus, limit the probability of conflicts. The relationship between capacity and structure is analysed by Sunil et al. [Reference Sunil, Hoekstra, Ellerbroek, Bussink, Vidosavljevic, Delhaye and Nieuwenhuisen40, Reference Sunil, Ellerbroek, Hoekstra, Vidosavljevic, Amtzen, Bussink and Nieuwenhuisen41] for four structures: full mix, layers, zones and tubes, each further constraining the traffic flow movement together with possible resolution methods. For the first three, separation is defined based on distance, for the 4D tubes time-based separation is used. Simulation results show that airspace structure has an influence on conflict probability, the best safety, efficiency and stability metric can be achieved by segmenting traffic into altitude bands according to heading range.

The introduction of structure to our framework, and studying its effect on safety and efficiency is part of our future plans. The current paper, however, focuses on how to apply constraints and structure just when they are necessary. These constraints are applied by changing the speed distribution. Also, it is unlikely that PAVs could occupy several altitude bands, therefore, introduction of structure into the UAM airspace may be difficult.

Since all analysis in this paper is performed in 2D, the methodology may be compared to land-based systems. Self-organisation and emergent formations of a swarm of autonomous, car-like vehicles was investigated by Vanualailai et al. [Reference Vanualailai42, Reference Kumar, Vanualailai, Sharma, Chaudhary and Kapadia43]. The papers investigate the emergent dynamics of a Lagrangian swarm. It is shown using computer simulations that the resulting formations are cohesive, well-spaced and bounded around the centroid of the formation. The focus of the papers, however, is on the patterns arising and formation keeping. Therefore, the results are not directly comparable.

Macroscopic modelling, MFD (Macroscopic Fundamental Diagram) and aggregate level description are well-established methods in the field of urban ground mobility. MFDs describe the relations between traffic flow, density and velocity [Reference Geroliminis and Daganzo44]. MFDs can be used to apply high-level control to regulate ensemble dynamics. The control goal is to maximise throughput and stabilise traffic flow. A feedback regulator-based gating method is described by Keyvan-Ekbatania et al. [Reference Keyvan-Ekbatania, Kouvelasb, Papamichaila and Papageorgiou45]. For control design, a basic dynamics model is used, developed based on the NFD (Network Fundamental Diagram) of the entire network. For a heterogeneous urban network, an optimal perimeter traffic flow controller can be derived analytically, as discussed by Aalipour et al. [Reference Aalipour, Kebriaei and Ramezani46]. These papers show how important modelling aggregate dynamics and traffic flows are for proper control design. However, these models are generally applicable for ground traffic networks.

The papers referenced provide an overview of the research related to general air traffic management, conflict probability calculation and conflict resolution. However, no methods have yet been advised to bound the conflict probability without having hard restrictions on airspace structure. This paper presents a novel coordination method in which constraints on the trajectory, namely, changed in speed are applied only when they are necessary. This is in line with the general concept of free flight. In this paper, the regulation of conflict probability is investigated by modifying ensemble speed distribution using the distributed stochastic task allocation algorithm. This enables a real-time, distributed coordination method. Because of the appropriately defined parameters the current state of the agent can be efficiently disseminated, therefore the communication burden of the algorithm is low.

3.1 Homogeneous Markov Chains (HMC)

The probabilistic guidance law based on Markov chains was introduced an later elaborated by Açikmese and Bayard [Reference Acikmese and Bayard28, Reference Acikmese and Bayard47]. A predefined desired formation (desired density distribution) of agents was achieved using an open-loop (no feedback from the actual state), stochastic, decentralised algorithm, formulated the following way. Let $\boldsymbol{{x}}\left( t\right) $ be the vector of probabilities corresponding to the tasks/states/speeds at time t such that $\textbf{1}^T \boldsymbol{{x}}\left( t\right) = 1$ , with a given desired distribution $\boldsymbol{\Theta}$ . The probability vector evolution is given by Equation (1).

(1) \begin{equation} \boldsymbol{{x}} (t+1) = \boldsymbol{{M}}(t),\ \boldsymbol{{x}}(t) \end{equation}

The matrix $\boldsymbol{{M}}$ is the Markov matrix containing the state transition probabilities. The next definite state is given as $\boldsymbol{{y}}\left( t+1\right) $ according to Equations (2) and (3), where $\boldsymbol{{y}}$ is a vector with a single 1 as the jth element, and all other elements 0.

(2) \begin{gather} \sum_{l = 1}^{j-1} \boldsymbol{{x}} \left[l \right] \leq z \leq \sum_{l = 1}^{j} \boldsymbol{{x}} \left[l \right] \end{gather}

(3) \begin{gather} \boldsymbol{{y}}\left( t+1\right) \left[j \right] = 1 \end{gather}

In Equation (2) z is a random number with uniform probability distribution such that $z \in \left[0, 1 \right] $ holds.

The Markov matrix $\left(\boldsymbol{{M}} \right) $ is synthesised using the Metropolis-Hastings algorithm, as given by Equation (4), where $\boldsymbol{{I}}$ denotes the identity matrix.

(4) \begin{equation} \boldsymbol{{M}}[i, j] = \begin{cases} \boldsymbol{{K}}[i, j]\ \boldsymbol{{F}}[i, j] & \text{if } \ i \neq j\\ \\[-7pt] \boldsymbol{{K}} [j, j] + \sum_{l \neq j} \left(\boldsymbol{{I}} - \boldsymbol{{F}} \left[l, j \right] \right) \boldsymbol{{K}} \left[l, j \right] & \text{if } \ i = j \end{cases} \end{equation}

Since not all state transitions are physically possible (due to acceleration constraints), an adjacency matrix $\left( \boldsymbol{{A}}_a\right) $ defined in Equation (5) indicates permitted transitions (Equation (6), $\odot$ representing the Hadamard product).

(5) \begin{equation} \boldsymbol{{A}}_a[i, j] = \begin{cases} 1 & \text{if transition from } i \text{ to } j \text{ is possible} \\ 0 & \text{if transition from } i \text{ to } j \text{ is not possible} \end{cases} \end{equation}

(6) \begin{equation} \left(\textbf{1} \textbf{1}^T - \boldsymbol{{A}}_a^T \right) \odot \boldsymbol{{M}} = \textbf{0} \end{equation}

The elements of the proposal matrix follow Gaussian distribution (Equation (7)).

(7) \begin{equation} \boldsymbol{{K}}\left[j, i \right] = \begin{cases} g\left( \|c_j - c_i \| ; 0, \delta_i \right) & \text{if } \ \boldsymbol{{A}}_a [i, j]=1 , i \neq j\\ \\[-7pt] 0 & \text{if } \ \boldsymbol{{A}}_a [i, j] = 0\\ \\[-7pt] 1 - \sum_{l \neq j} \boldsymbol{{K}} \left[l, j \right] & \text{if } \ j = i \end{cases} \end{equation}

In Equation (6) $c_i$ is the vector describing the states and the function $g\left(z; \eta, \delta \right) $ (Gaussian distribution) is defined as in Equation (8).

(8) \begin{equation} g\left(z; \eta, \delta \right) = \frac{1}{\delta \sqrt{2 \pi}} \mathrm{e}^{-\frac{\left(z-\eta \right)^2}{2 \delta^2}} \end{equation}

The elements of the acceptance matrix are subject to the constraint in Equation (9).

(9) \begin{equation} 0 \leq \boldsymbol{{F}} [i, j] \leq \min \left(1, \boldsymbol{{R}}[i, j] \right) \end{equation}

Acceptance matrix elements are determined as in Equation (10), with $\alpha \in \left(0, 1 \right] $ being a design parameter.

(10) \begin{equation} \boldsymbol{{F}} [i, j] = \alpha \min \left(1, \boldsymbol{{R}}[i, j] \right) \end{equation}

In Equations (9) and (10) the matrix $\boldsymbol{{R}}$ is defined according to Equation (10), where $\boldsymbol{\theta}$ is the vector of desired (steady-state) distribution.

(11) \begin{equation} \boldsymbol{{R}}[i, j] =\frac{\boldsymbol{\theta}[i] \boldsymbol{{K}}\left[j, i \right] }{\boldsymbol{\theta}\left[j \right] \boldsymbol{{K}}[i, j] } \end{equation}

Using Equations (4)–(11) the Markov matrix can be constructed, which can be used to determine states at the next time step (Equation (1)).

Speed at the next time step is determined according to Equation (12), where $\boldsymbol{{V}}$ is a prescribed speed vector.

(12) \begin{equation} v(t+1) = \boldsymbol{{V}}^T \boldsymbol{{y}}(t+1) \end{equation}

The main advantage of this method is that only the desired state needs to be transmitted to agents, who realise the Markov chain independently, no communication about the state is necessary. As no feedback is incorporated, the Markov matrix is constant in time, which leads to a homogeneous Markov chain (HMC). Asymptotic convergence of swarm distribution is guaranteed to desired distribution if the spectral radius condition is satisfied. The algorithm was extended to handle collaborative sub-swarms (composite swarms) using consensus techniques. A Linear Matrix Inequality (LMI) formulation of the problem has also been elaborated ensuring convergence with a prescribed exponential rate, also minimising a cost function. Scalability was demonstrated by applying the method with the desired distribution given as a greyscale image [Reference Acikmese and Bayard47].

3.2 Inhomogeneous Markov Chains (IMC)

The swarm guidance method based on Markov chains was later extended by Bandyopadhyay et al. [Reference Bandyopadhyay, Chung and Hadaegh29] by incorporating feedback from the current distribution state into the transition matrix synthesis. This lead to a closed-loop strategy by generating inhomogeneous Markov chains. Markov matrix changes in time due to the feedback, thus, the Markov chain is inhomogeneous (IMC). The algorithm was proven to converge and minimise transition costs. Apart from the Markov chain construction algorithm (Equation (4)), the same speed/task/bin allocation method can be used as in the case of homogeneous Markov chains, the optimal Markov chain is constructed using the following algorithm. The feedback gain is the Hellinger distance (for the kth agent at time instant t given by Equation (13)) between the current and desired distribution, if it is sufficiently low, the transition matrix is changed to identity to avoid unnecessary transitions.

(13) \begin{equation} \xi_t^k = D_H \big(\boldsymbol{\Theta}, \mu_t^k \big) = \frac{1}{\sqrt{2}} \sqrt{\sum_{i } \left( \sqrt{\Theta [i]} - \sqrt{\mu_t^k [i]}\right) ^2} \end{equation}

(14) \begin{equation} \mu_t^* [i] = \frac{1}{m_t} \sum_{k = 1}^{m_k} r_t^k [i] \end{equation}

In Equation (13) $\mu_t^k$ denotes the estimation of agent k at time t of the total agent distribution $\left(\mu_t^* \right) $ defined as in Equation (14), $m_t$ denoting the total number of agents at time t, and $r_t^k$ is the current task of agent k at time t.

As an estimation method, Bayesian Consensus Filtering [Reference Bandyopadhyay and Chung48] may be used. Given the feedback $\left( \xi\right) $ and the transition cost matrix $\left(\boldsymbol{{C}} \right) $ , the problem can be formulated as a linear program subject to constraints giving the optimal Markov matrix as in Equation (15).

(15) \begin{equation} \boldsymbol{{M}}_t^k [i, j] = \begin{cases} \varepsilon_M \ \xi_t^k \ \boldsymbol{\Theta} \left[j \right] \left(1 - \frac{\boldsymbol{{C}}_t [i, j] }{C_{t, \max} + \varepsilon_C} \right) & \text{if } \ i \neq j\\ \\[-7pt] 1 - \sum_{l \neq i} \boldsymbol{{M}}_t^k \left[i, l \right] & \text{if } \ i = j \end{cases} \end{equation}

Parameters $\varepsilon_M \in \left(0, 1 \right] $ and $\varepsilon_C \in {\rm I\!R^+}$ in Equation (15) are constant design parameters. The technique is also extended with a time-dependent method to facilitate initial convergence and reduce undesirable transitions after sufficient convergence and a deterministic procedure to escape transient bins.

3.3 Local information consistency – Inhomogeneous Markov chains

The feedback-based swarm guidance algorithm was further developed Jang et al. [Reference Jang, Shin and Tsourdos30] to function without having to rely on global information or an estimation of the distribution of the entire swarm. This is called Local Information Consistency Assumption (LICA) as opposed to the general Global Information Consistency Assumption (GICA). By using local information only the algorithm remains scalable with the number of agents involved, an important aspect when utilised for traffic management on a large scale. Since the method is fundamentally applied to probabilistic swarm guidance in a bin partitioned world, in the paper local is defined in terms of neighbouring bins (transition is possible within a single time step) as opposed to a definition based on communication range used in this paper. It is assumed that aircraft cooperate, and share navigation data with all others within their communication range, which are considered to be the local quantities. Local current and desired distribution are defined as given by Equations (16) and (17) respectively. The convergence of local quantities is equivalent to that of global quantities if the topology graph of the bins is at least bidirectional and strongly connected.

(16) \begin{equation} \overline{\mu}_t^*[i] = \frac{1}{m_{k, neighbour}} \sum_{k = 1}^{m_{k, neighbour}} r_t^k [i] \end{equation}

(17) \begin{equation} \overline{\Theta} [i] = \frac{\Theta [i] }{\sum_{\text{neighbour bins} \ n} \Theta \left[n \right]} \end{equation}

The feedback is defined according to Equation (18), saturated to $\left[\varepsilon_\xi, 1 \right] $ , $\alpha$ denotes the sensitivity parameter. Markov matrix synthesis is similar to Equation (15), and is given by Equation (19).

(18) \begin{equation} \overline{\xi}_t^k [i] = \left(\frac{\left\vert \overline{\Theta} [i] - \overline{\mu}_t^* [i] \right\vert }{\overline{\Theta} [i]} \right)^\alpha \end{equation}

(19) \begin{equation} \boldsymbol{{M}}[i, j] = \begin{cases} \varepsilon_M \ \max \left( \overline{\xi}_t^k [i], \overline{\xi}_t^k \left[j \right]\right) \boldsymbol{\Theta} \left[j \right] \left(1 - \frac{\boldsymbol{{C}}_t [i, j] }{C_{t, \max} + \varepsilon_C} \right) & \text{if } \ i \neq j\\ \\[-7pt] 1 - \sum_{l \neq i} \boldsymbol{{M}}_t^k \left[i, l \right] & \text{if } \ i = j \end{cases} \end{equation}

An important feature of this LICA-based method is that using local information only, and enables an asynchronous implementation. Algorithms for Markov chain synthesis under flux upper limits, and, for local information based quorum model are also reported by Jang et al. [Reference Jang, Shin and Tsourdos30].

The above-described algorithms were implemented with a behaviour/task allocation perspective to achieve the traffic coordination goals introduced. Within the context of the paper, $\boldsymbol{{x}}$ is the probability vector corresponding to potential speed states, $\boldsymbol{\Theta}$ is the desired speed state distribution of vehicles, $\boldsymbol{{A}}_a$ is the matrix representing the acceleration constraints, and $\boldsymbol{{V}}$ is the vector of potential speeds within limits a and b.

4.1 Simulator platform

The primary aim of the simulations is to determine whether distributed task allocation algorithms are applicable for stochastic traffic flow management to enable a high-level optimisation layer to control ensemble dynamics and characteristics using a single parameter.

As it is not common in air traffic management research, to carry out the analysis and study the emergent dynamics, a custom simulation environment was created in Matlab.

This way the fidelity and complexity of the simulation modules could be appropriately decided with full control over all the parameters and the possibility to easily extend the framework. Since there are no external data flows to be considered, all the vehicles are synthetic and the point of the simulation is to explore the system rather than individual vehicles, the simulation time is different than real-time.

The architecture of the simulation framework is based on an agent-based modelling approach. The system is solved in a time-based manner, with fixed time steps. Each vehicle is modelled as an agent which moves and makes decisions according to the implemented logic and the location of surrounding vehicles and the information transmitted by them. The speed for the next simulation step is decided based on the task allocation algorithm run individually by each agent, respecting the acceleration constraints. Then, a motion step is taken by the vehicle. A simulation step is completed once all agents made a decision and a movement step.

In case of centralised task allocation algorithms the Markov matrix is generated centrally and passed on to each agent. For decentralised algorithms, it is generated based on the state information of the surrounding vehicles. Since for the communication model a disk radio is used, the decision of whether the information needs to be taken into account from another vehicle becomes a distance calculation. When task allocation is asynchronous, agents are put in a random order, then synthesise the Markov matrix and update the state sequentially.

During the simulations conflict events are detected and stored centrally, with the ID of all vehicles involved. The time from the sensing range to the conflict range is measured. Vehicle states and speeds are also recorded, together with state changes. Based on these, the conflict rate and throughput metrics, the number of transitions and the distribution distance from the desired can be calculated.

Since the results for iterations and different parameters are independent in steady-state, the simulations can be run in parallel. Therefore, to run the simulations a high-performance computing facility was used. The general run time was of the order of a few hours for the set of distribution parameters.

In this paper, a free flight scenario is analysed with en-route traffic, maintaining level flight, therefore it is sufficient to handle the problem in 2D. Conflicts are defined as stated in Section 1 and elaborated below. As a behaviour (task), speed changes are considered only, conflict resolution via heading changes are preserved for the deterministic low level.

4.2 Simulation scenario

The simulations are performed in a $100\,\mathrm{nmi}^2$ region $\left(10\,\mathrm{nmi} \times 10\,\mathrm{nmi} \right) $ as an “infinite world”, in order to maintain the total number of agents.

An “infinite world” means that once an agent crosses the border of the simulation area it reappears on the opposite side, with all other parameters (state, speed, heading, etc.) unchanged. Aircraft density was chosen in line with urbanisation trends and market penetration predictions based on the analysis by Sunil et al. [Reference Sunil, Hoekstra, Ellerbroek, Bussink, Nieuwenhuisen, Vidosavljevic and Kern49]. Considering the area simulated, the expected density and that the number of vehicles need to the integers, simulations were performed with the following agent numbers: $165, \ 211, \ 258, \ 305$ . The resulting densities of vehicles match the scenarios analysed by Sunil et al.

Within the simulation, aircraft models need to be representative of urban air mobility vehicles. Considering the diversity of such vehicles, instead of modelling a specific type of vehicle, a surrogate model is developed based on kinematics and certain dynamic characteristics, which makes the analysis more general.

For aircraft kinematics a simple integrator model is considered as given by Equations (20)–(23), where $\psi$ denotes the heading, $\Delta t$ is the discrete time step and $\left( x, y\right) $ is the position of the aircraft. The speed of the aircraft is denoted by v, $v_x$ and $v_y$ are the velocity components in the x and y direction.

(20) \begin{align} v_x (t) &= v (t) \cos (\psi) \end{align}

(21) \begin{align} v_y (t) &= v (t) \sin (\psi) \end{align}

(22) \begin{align} x (t) &= x \left(t-1 \right) + v_x (t)\ \Delta t \end{align}

(23) \begin{align} y (t) &= y \left(t-1 \right) + v_y (t)\ \Delta t \end{align}

The discrete time step was chosen based on rough approximations and estimations of the dynamics of personal aerial vehicles [Reference Gerboni, Venrooij, Nieuwenhuizen, Joos, Fichter and Bülthoff50, Reference Geluardi, Venrooij, Olivari, Bülthoff and Pollini51]. Since the aim of the middle layer is to determine the prescribed behaviour of the bottom level and provide a reference to the guidance algorithm, therefore the time step size was selected considering the settling time of the step response dynamics, which is $\Delta t = 5\,\mathrm{s}$ for both longitudinal velocity transfer functions in the cited references. The modelled dynamics are generally representative of a wide range of vehicles, independently of their architecture and design.

4.3 Speed distributions

For the speed distributions, the truncated exponential distribution (Equation (24)) and the truncated normal distribution (Equation (25)) were selected. Truncated exponential distribution is given by Equation (24). Lower limit is a, upper limit is b, rate is $\lambda$ .

(24) \begin{equation} f\!\left(z;\ \lambda, a, b \right) =\frac{\lambda \mathrm{e}^{- \lambda z}}{\mathrm{e}^{- \lambda a} - \mathrm{e}^{- \lambda b}} \end{equation}

Truncated normal distribution is given by Equation (25). Lower limit is a, upper limit is b, mean is $\eta$ , standard deviation is $\sigma$ .

(25) \begin{equation} f\!\left(z;\ \eta, \sigma, a, b \right) =\frac{\phi\!\left(\frac{z - \eta}{\sigma} \right) }{\sigma \!\left(\Phi\left(\frac{b-\eta}{\sigma} \right) - \Phi\!\left(\frac{a-\eta}{\sigma} \right) \right) } \end{equation}

The probability density function of the standard normal distribution is given by Equation (26), the cumulative distribution of the standard normal distribution is in Equation (26) (where erf is the error function).

(26) \begin{equation} \phi \left(\xi \right) = \frac{1}{\sqrt{2 \pi}} \mathrm{e}^{\left(\!-\frac{1}{2} \xi^2 \right) } \end{equation}

(27) \begin{equation} \Phi \left(\xi \right) = \frac{1}{2} \left(1 + \textrm{erf} \left(\frac{\xi}{\sqrt{2}} \right) \right) \end{equation}

As a distribution parameter the rate $\left(\lambda \right) $ is considered in case of the exponential distribution and the mean $\left(\eta \right) $ in case of the normal distribution.

Speed limits were chosen taking into consideration the operation in an urban environment. Therefore the minimum speed is set as $a = 15\,\mathrm{Kn}$ (below this speed, manoeuvres are considered low-speed [Reference Gerboni, Venrooij, Nieuwenhuizen, Joos, Fichter and Bülthoff50, Reference Perfect, Jump and White52]), the maximum speed is set as $b = 180\,\mathrm{Kn}$ . Thus, the average roughly matches the speed forecast by Perfect et al. [Reference Perfect, Jump and White52] for the decision variable limits, and the range is in harmony with the characteristic speeds [Reference Elevate2, Reference Sunil, Ellerbroek, Hoekstra, Vidosavljevic, Amtzen, Bussink and Nieuwenhuisen41].

Decision variable limits are selected so that the mean speed values change between 20 and $90\,\mathrm{Kn}$ in order to retain the necessary freedom and flexibility. For the exponential distribution the range was chosen as $\lambda \in \left[0.0025, 0.2 \right] $ . For the normal distribution the standard deviation is kept constant at $\sigma = 10$ , and the mean is adjusted as $\eta \in \left[5, 90 \right] $ .

The dynamics model has a relationship with the task allocation adjacency matrices since it governs the available subsequent states (speeds), which is based on the level flight acceleration limits of the aircraft. Considering passenger comfort, the acceleration limit is specified as $1\mathrm{Kn/s} \left(\approx 0.5\mathrm{m/s^2} \right) $ [Reference Hoberock53], which seems to match the average acceleration capabilities of personal aerial vehicles [Reference Gerboni, Venrooij, Nieuwenhuizen, Joos, Fichter and Bülthoff50]. Variation of acceleration limits with speed is not considered here. Given the time step of $\Delta t = 5\mathrm{s}$ the proposal and adjacency matrices are defined to enable transitions to states with a speed deviation of $\pm 5\mathrm{Kn}$ . Based on dynamic limits, the minimum speed was chosen as $a = 15\mathrm{Kn}$ , the maximum speed as $b = 180\mathrm{Kn}$ and the acceleration limit for the $\Delta t = 5\mathrm{s}$ timestep as $5\mathrm{Kn}$ . Thus, the available states were created by dividing the range between the minimum and the maximum speed into $5\mathrm{Kn}$ wide segments $\left(15, 20, \ldots, 180 \right)$ .

4.4 Description of the environment

Since the scenario considered is personal aerial vehicle traffic in an urban environment reduced separation margins are used. For levelled flight only horizontal separation is meaningful, which was defined as $R_C = 0.135\mathrm{nmi}$ [Reference Sunil, Ellerbroek, Hoekstra, Vidosavljevic, Amtzen, Bussink and Nieuwenhuisen41]. Since cooperation between aircraft is assumed and separation is achieved by a distributed algorithm, communication range has to be defined. The nominal range of ADS-B is greater than $90\text{nmi}$ (air-to-air) for low interference (FRUIT – False Replies Unsynchronised to Interrogator Transmission) conditions [54]. For a high-traffic-density scenario, such a high range is neither necessary nor probable (because of the high interference due to high communication density). It is assumed that the communication range scales similarly to separation limits (standard ATM separation for ADS-B equipped aircraft [55] of $5\mathrm{nmi}$ reduced to $R_C = 0.135\mathrm{nmi}$ ), thus, the data broadcast range used in this paper is $2.5\mathrm{nmi}$ . For the LICA based algorithm, local states and quantities are the ones within this range. To verify that data dissemination is possible, the density capacity of the ADS-B protocol was estimated. Media access is based on TDMA: each agent transmits at most $6.2$ messages of $120\mathrm{\mu s}$ on average every second [56], therefore ideally 1,344 aircraft can be handled simultaneously. For the simulated agent densities, the number of aircraft within the communication range is significantly less. Hence it is assumed that the communication part of the system introduces no errors.

The intention of the middle layer is to reduce conflict probabilities, the probability of simultaneous conflicts and conflict propagation. Moreover, it should ensure that enough time is reserved to resolve conflicts should they occur. Therefore, a conflict event is considered to be the incident when the traversing time of an intruder aircraft from the sensing range (communication range, $2.5\mathrm{nmi}$ ) to the conflict range $\left( R_C = 0.135\mathrm{nmi}\right) $ is less than a predefined threshold. When selecting the conflict threshold time several factors have to be taken into account, for example, the performance of conflict detection and conflict resolution algorithms, passenger comfort levels and their relationship with acceleration and its time derivative (investigated e.g. by Dousse et al. [Reference Dousse, Heitz, Schill and Floreano57]). For the simulations in this paper, threshold times of $30, \ 40, \ 50\ \text{and} \ 60\mathrm{s}$ are investigated. Similarly to TCAS [58], a $5\mathrm{s}$ transition time is assumed based on settling time of the yaw transfer function [Reference Gerboni, Venrooij, Nieuwenhuizen, Joos, Fichter and Bülthoff50], with a worst-case relative speed of $b = 360\mathrm{Kn}$ . For traffic management GPS and ADS-B technologies, separation limits have long been expected to become time-based, defined according to wake turbulence intensity [Reference Krozel, Peters and Hunter59]. The highest $60\mathrm{s}$ separation corresponds to the current limit for light aircraft in a terminal manoeuvring area [Reference Morris, Peters and Choroba60]. Given the allowable collision probability threshold, the separation intervals can be determined, similarly to the method presented by Bell [Reference Bell61].

4.5 Performance metrics

The quality metrics considered reflect on the original objectives, which is to achieve the highest traffic throughput while complying with safety requirements. To evaluate the performance from a safety perspective the number of conflicts normalised by the number of agents and the time elapsed is calculated. The definition of the safety metric (S) is given by Equation (28).

(28) \begin{equation} S = \frac{C}{n_{{agent}} \times t} \end{equation}

In Equation (28) the total number of agents is denoted by $n_{agent}$ , C is the total number of conflicts.

To analyse the effects on throughput the mean speed of the agents is evaluated, taking into account that a conflict resolution manoeuvre results in detours and reduction in speed towards the destination, therefore conflict events are penalised. The throughput (T) metric is thus defined as the average ratio of distance travelled towards the target and the total elapsed time, given by Equation (29).

(29) \begin{equation} T = \frac{1}{n_{{agent}}}\left( \frac{t\left( \ \sum_{i } v_i\right)^2}{t \ \sum_{i } v_i + C \times P}\right) \end{equation}

In Equation (29) $v_i$ is the speed of the $i \text{th}$ agent and P is the penalty associated with conflicts, which was defined to be the difference between the direct and deviated path length, assuming both agents seek total separation distance for safety reasons (Fig. 1 and Equation (30)).

(30) \begin{align} P &= 2 \times 360 \left(\mathrm{Kn}\right) \times 10 \left(\mathrm{s}\right) + 2 \times \frac{2.5 - 360 \left(\mathrm{Kn}\right) \times 10 \left(\mathrm{s}\right)}{\cos \left(\arcsin \left(\frac{0.135}{2.5 - 360\times 10 } \right) \right) } - 5 \\[3pt] &= 5.0091 - 5 = 0.0091\mathrm{nmi} \end{align}

Figure 1. Illustration of the deviation penalty (P). The ownship is the aircraft from the viewpoint of which the scenario is considered. Any other aircraft is considered to be an intruder.

It is assumed that both aircraft start their coordinated conflict avoidance manoeuvres once they are within each other’s sensing range $\left( r_S\right)$ , and seek a distance of $r_C$ from the original path of the other aircraft.

Furthermore, the number of transitions is counted to help assess the practical applicability of the method.

Steady-state simulations are performed for 900 time steps, making sure that the sliding window mean variation of the metrics decreases below 1%.

Simulations with the same setup and parameters are performed 128 times to account for the stochastic nature of the framework. The number of repetitions is chosen after a convergence analysis, together with requirements from the computing service.

Initial positions and headings of agents are chosen randomly according to a uniform distribution. For steady-state simulations, initial states (speeds) are chosen according to the distribution parameter.

5.1 Configuration

To analyse the emergent dynamics in steady-state, simulations were run with the distribution parameter kept fixed. The simulation was initialised by setting vehicle positions randomly in the simulation area, according to a uniform distribution. Headings were also generated according to a uniform distribution. Initial speeds were assigned according to the desired distribution with the parameter matching the studied case. Simulations were repeated to account for their stochastic nature.

5.2 Simulation results for steady state

Based on the simulation results the effect of distribution types and parameters, agent density and performance of algorithms can be compared.

Figure 2 shows the relationship between the safety metric and the distribution parameter for an exponential desired distribution, using the HMC algorithm. This relationship proves that by varying the distribution parameter the safety metric can be regulated.

Figure 2. Safety metric as a function of the exponential distribution parameter ( $\lambda$ ) for the HMC algorithm. (Threshold time: $60\mathrm{s}$ . Dotted lines show the standard deviation. Note that not the full parameter range is shown. The distribution parameter is related to the speed distribution of vehicles.)

For the exponential desired distribution the safety metric is effectively zero for $\lambda > 0.07$ , and rapidly increases with decreasing parameter values (increasing mean speed). The effect of agent density can also be seen, with increasing density the number of conflicts and the safety metric increases as well.

The throughput metric also increases with decreasing parameter values (Fig. 3), however, it is less affected by the agent density.

Similar trends can be seen when the desired distribution is normal, the safety metric increases with increasing parameter (increasing mean speed) and agent density (Fig. 4).

Figure 3. Throughput metric as a function of the exponential distribution parameter ( $\lambda$ ) for the HMC algorithm. (Threshold time: $60\mathrm{s}$ . Dotted lines show the standard deviation. The distribution parameter is related to the speed distribution of vehicles.)

Figure 4. Safety metric as a function of the normal distribution parameter ( $\eta$ ) for the HMC algorithm. (Threshold time: $60\mathrm{s}$ . Dotted lines show the standard deviation. Note that not the full parameter range is shown. The distribution parameter is related to the speed distribution of vehicles.)

Throughput also increases with increasing parameter ( $\eta$ ) (Fig. 5).

Figure 5. Throughput metric as a function of the normal distribution parameter ( $\eta$ ) for the HMC algorithm. (Threshold time: $60\mathrm{s}$ . Dotted lines show the standard deviation. The distribution parameter is related to the speed distribution of vehicles.)

The appropriate functioning of the algorithm is verified by Figs 6 and 7, which show the mean speed of agents for all agent densities as a function of the distribution parameter. Mean speed changes with the variation of the parameters as expected.

Figure 6. Mean speed as a function of the exponential distribution parameter ( $\lambda$ ) for the HMC algorithm. (Dotted lines show the standard deviation. The distribution parameter is related to the speed distribution of vehicles.)

Figure 7. Mean speed as a function of the normal distribution parameter ( $\eta$ ) for the HMC algorithm. (Dotted lines show the standard deviation. The distribution parameter is related to the speed distribution of vehicles.)

In Fig. 8 the desired distributions are compared by plotting the safety metric as a function of mean speed, showing that for the same mean speed the safety metric is lower (fewer conflicts) for normal distribution than it is for the exponential distribution.

Figure 8. Safety metric as a function of mean speed for the HMC algorithm. (Threshold time: $60\mathrm{s}$ .)

Normalising the safety metric with agent density shows that the dependence is linear, normalised metrics for a given distribution are almost indistinguishable (Fig. 9).

Figure 9. Safety metric normalised by agent number as a function of mean speed for the HMC algorithm. (Threshold time: $60\mathrm{s}$ .)

The choice of the conflict threshold time is based on estimated dynamic and detect-and-avoid performance, however, it is arbitrary to some extent. The relationship between the time threshold and the number of conflicts was investigated by analysing the histograms of traversal times from the sensing to the conflict region. The traversal time is the time it takes the intruder to arrive at the conflict range of the ownship after it had reached the sensing range. A representative example for traversal times for both distributions with the same mean speed can be seen in Figs 10 and 11, indicating chosen threshold times.

Figure 10. Traversal times for exponential distribution using the HMC algorithm. (Red lines showing threshold times. Agent density: 165 agents, $\lambda = 0.0025$ .)

Figure 11. Traversal times for normal distribution using the HMC algorithm. (Red lines showing threshold times. Agent density: 165 agents, $\eta = 90$ .)

Traversal times for the exponential distribution are generally lower, and more evenly distributed. For the normal distribution, the traversal times heavily peak at $t = 50\mathrm{s}$ , showing the sensitivity of metrics to threshold selection.

Comparison of algorithms in terms of safety and throughput can be seen in Figs 12 and 13, showing the performance of all algorithms is identical.

Figure 12. Safety metric as a function of mean speed for all algorithms. (Threshold time: $60\mathrm{s}$ . Agent density: 165 agents.)

Figure 13. Throughput metric as a function of mean speed for all algorithms. (Threshold time: $60\mathrm{s}$ . Agent density: 165 agents.)

It is also important to assess the number of state transitions needed to achieve the desired state distribution, a comparison can be seen in Fig. 14. Since there is no feedback mechanism in the HMC algorithms transitions still occur in steady-state, therefore the number of transitions is three orders of magnitude higher than for IMC algorithms. The global information-based IMC algorithm generates the least number of transitions, the algorithms using local information only require a higher number of transitions within a time step (compared to global information IMC), however, this number seems to be independent of whether the transitions occur simultaneously or asynchronously.

Figure 14. Number of transitions generated by each algorithm for a given time step. (Average value of all simulation repeats. Agent density: 305 agents.)

The final distance form the desired distribution (total variation distance) achieved by the algorithms can be seen in Figs 15 and 16. The IMC algorithm generally performs better than HMC, achieving a smaller distance, however, both are outperformed by the local information based IMC algorithm, especially for parameters related to lower mean speeds.

Figure 15. Distribution distance for all algorithms with exponential desired distribution as a function of the distribution parameter ( $\lambda$ ). (Moving average, windows size: 5. Agent density: 305 agents. The distribution parameter is related to the speed distribution of vehicles.)

Figure 16. Distribution distance for all algorithms with normal desired distribution as a function of the distribution parameter ( $\eta$ ). (Moving average, windows size: 5. Agent density: 305 agents. The distribution parameter is related to the speed distribution of vehicles.)

Conflicts are also classified according to the number of aircraft involved. The number of conflicts involving two aircraft can be seen in Fig. 17 (for the HMC algorithm, results are almost identical for all used algorithms), showing tendencies similar to the safety metric, indicating that with higher mean speed or agent density the number of conflicts increases.

Figure 17. Number of conflicts involving two aircraft at the same time. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC.)

Normalising the results with the cube of the number of agents involved the curves are almost identical, demonstrating a quadratic relationship between agent density and the number of conflicts, independent of the speed distribution (Fig. 18).

Figure 18. Number of conflicts involving two aircraft at the same time normalised by the square of the number of aircraft simulated. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC.)

The number of conflicts involving three aircraft is shown in Fig. 19, with the trend being similar to the trends of conflicts involving two aircraft only, but values about two orders of magnitude lower.

Figure 19. Number of conflicts involving three aircraft at the same time. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC.)

When normalising the number of triple conflicts with the cube of the agents simulated the results become similar for all agent densities, indicating a cubic relationship. This also explains why the safety metric scales linearly with the number of agents, which is because conflicts predominantly involve two aircraft only, and the number of conflicts is already divided by the number of agents according to Equation (28). The number of conflicts involving more than three aircraft at the same time was negligible, their probability requires further simulations and analysis, possibly using techniques for rare event simulation. The explanation for the observed scaling laws is probably related to the uniform distribution of aircraft, the probability that two aircraft reach close proximity scales quadratically with the number of agents (density of agents in proximity scales linearly for each aircraft). The probability of the proximity of a third aircraft to two aircraft already in conflict also scales linearly, explaining the cubic relationship.

To further analyse the dynamics of intrusion, five ranges were defined around the own aircraft, with logarithmically spaced radii $\left(2.5, \ 1.2, \ 0.58, \ 0.28, \ 0.135\mathrm{nmi} \right) $ . The probabilities of transitioning to an inner radius, provided that the intruder has already arrived at the current radius have been calculated, and are given in Table 1. Results show that transition probabilities are independent, and are closely related to the ratios of radii of ranges.

Steady-state simulations show that the safety and throughput metric can be governed through the choice of the distribution parameters, applying the stochastic task allocation algorithms. Specifying normal distribution as desired generally leads to better results in terms of safety compared to the exponential distribution. Conflicts events involving a number of aircraft scale with the power of the number of aircraft involved, thus, the safety metric scales linearly with the number of aircraft. All algorithms have a similar performance in terms of metrics, however, local information based algorithms can achieve the desired distribution more accurately, and IMC algorithms, in general, require a much lower number of transitions.

Table 1. Transition probabilities

6.1 Configuration

To assess the dynamic performance of the algorithms transition simulations are performed starting all agents from the first state (lowest speed).

The distribution parameter is kept fixed, the time it takes for the system to reach the desired state and its dynamic characteristics are observed. Step response simulations are also performed to determine the time constants of the system.

The calculation of the safety metric was modified, a moving average is calculated to provide more relevant information about the conflict state of the system.

Information on conflicts, conflict rate, throughput, traversal times, distribution distance and the number of transitions were collected. Simulations were repeated to account for their stochastic nature.

6.2 Transient performance of homogeneous Markov chains

Results show that for parameters corresponding to lower speeds the rate of convergence is satisfactory, however, for higher speed parameters the rate is acceptable for the normal desired distribution only. For mean speed and throughput, the transient resembles that of a first-order system for both distributions. For normal distribution both speed and throughput are well-converged by the time reaches 1,000s (Figs 20 and 21).

Figure 20. Throughput metric transient for normal distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 21. Mean speed transient for normal distribution. (Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

However, for exponential distribution (Figs 22 and 23) results are still changing at 5,000s, which questions the practical applicability of such a transition.

Figure 22. Throughput metric transient for exponential distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 23. Mean speed transient for exponential distribution. (Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

With the moving mean the transition of the safety metric converges with oscillations for normal distribution (Fig. 24), but fails to converge for exponential distribution (Fig. 25).

Figure 24. Safety metric transient for normal distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 25. Safety metric transient for exponential distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Distribution distance decreases rapidly compared to the initial distance (Figs 26 and 27), still, it remains quite high for the exponential distribution to hinder convergence.

Figure 26. Distribution distance transient for normal distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 27. Distribution distance for exponential distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Since for both distributions the HMC algorithm is used, the number of transitions is the same, as shown in Figs 28 and 29.

Figure 28. Number of transitions transient for normal distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 29. Number of transitions for exponential distribution. (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

6.3 Transient performance of inhomogeneous Markov chains

The transient performance of inhomogeneous Markov chain algorithms was also investigated. Figures show that the mean speed, and thus, the metrics fail to converge to the desired values in a reasonable time frame. This is because IMC algorithms are essentially designed for the trade-off between convergence rate and the number of transitions required. Thus, the transient performance of IMC algorithms falls short of that of HMC, even when state transitions are not penalised. Figures 30, 31 and 32 show the transient performance for mean speed using the IMC, LICA IMC and asynchronous LICA IMC algorithms respectively.

Figure 30. Mean speed transient for normal distribution. (Algorithm: IMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 31. Mean speed transient for normal distribution. (Algorithm: IMC LICA. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 32. Mean speed transient for normal distribution. (Algorithm: Asynchronous IMC LICA. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Even though it was been seen that IMC algorithms generate much fewer transitions, which are unnecessary in steady-state. Results show that for transients it is preferable to switch to HMC for large changes, and when the distribution distance has sufficiently dropped, switch back to IMC.

6.4 Transient performance between distributions

The transition from a normal distribution to an exponential one corresponding to the same mean speed and back was also investigated. This may correspond to the scenario outlined in Section 1, when the distribution has to be changed due to the appearance of an exceptional vehicle, and then back once it has left the area. Using the HMC algorithm, changing the desired distribution for the time span of 500–2,000s, the mean speed changes only slightly after the transitions (Fig. 33), together with a slight drop in throughput (Fig. 34, as expected based on steady state results in Fig. 13).

Figure 33. Mean speed transitioning from a normal to an exponential distribution (at $t = 500\mathrm{s}$ ) and back (at $t = 2,\!000\mathrm{s}$ ). (Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

Figure 34. Throughput transitioning from a normal to an exponential distribution (at $t = 500\mathrm{s}$ ) and back (at $t = 2,000\mathrm{s}$ ). (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

The predicted growth in the safety metric is also noticeable in Fig. 35.

Figure 35. Safety transitioning from a normal to an exponential distribution (at $t = 500\mathrm{s}$ ) and back (at $t = 2,000\mathrm{s}$ ). (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

The distribution distance evolves as expected with peaks at the distribution changes dropping with convergence (Fig. 36).

Figure 36. Distribution distance transitioning from a normal to an exponential distribution (at $t = 500\mathrm{s}$ ) and back (at $t = 2,000\mathrm{s}$ ). (Threshold time: $60\mathrm{s}$ . Algorithm: HMC. Agent number: 165. The distribution parameter is related to the speed distribution of vehicles.)

6.5 Time constants

To further analyse the transient performance of the system at different operating points, for each point the parameter corresponding to a mean speed $10\mathrm{Kn}$ higher and lower was calculated and set as a transient goal. Thus, the step responses of the system for the different operating points could be analysed. Figures 37 and 38 show representative examples of transitions.

Figure 37. Step response for mean speed with normal distribution (Algorithm: IMC for steady-state, HMC for transitions. The distribution parameter is related to the speed distribution of vehicles).

Figure 38. Step response for mean speed with exponential distribution (Algorithm: IMC for steady-state, HMC for transitions. The distribution parameter is related to the speed distribution of vehicles).

Second-order system models were fit to the step response models, dominant time constants corresponding to the parameters are shown in Figs 39 and 40.

Figure 39. Time constants for exponential distribution (Algorithm: IMC for steady-state, HMC for transitions. The distribution parameter is related to the speed distribution of vehicles).

Figure 40. Time constants for normal distribution (Algorithm: IMC for steady-state, HMC for transitions. The distribution parameter is related to the speed distribution of vehicles).

The fundamental reason for outliers is that fitting for a stochastic response is not exact, still, the general trend is clearly visible. For exponential distribution, because of the poor convergence, only a few fits were possible. However, the general tendency is still evident. For parameters related to smaller speeds, the convergence is faster, even though the step size was set to be the same. For exponential distribution, for higher speeds, this tendency can also be seen in Fig. 38.

The behaviour and performance of the system are now fully analysed in steady-state and transient conditions for a set parameter, thus, analysis for a varying or even controlled parameter becomes possible.

7.1 Configuration

Fixed distribution parameter simulations have confirmed that by changing the rate of the exponential distribution or the mean of the normal distribution the throughput and safety measures can be adjusted, therefore an attempt has been made to regulate these performance metrics. The overall goal of the middle layer is to enable ensemble dynamics to be controlled by the high-level optimisation by adjusting a single parameter only, therefore the simulation has been modified to allow for a feedback-based control.

The aim of the controller was to speed up transitions from one parameter to another while respecting limits on the safety metric. Since transient simulations showed that system dynamics are essentially first order, there is no limit on the proportional gain permissible considering the mean speed (and the closely coupled throughput metric) only, the limit is imposed only by the admissible safety metric. Figures 41 and 42 show results of transient simulations for normal desired distribution stepping from $\eta = 45$ to $\eta = 55$ .

Figure 41. Controlled mean speed with different proportional gains.

Figure 42. Controlled safety metric with different proportional gains.

For the transitions, a simple proportional controller is used, given by Equation (31) together with the evolution of the corresponding safety metric.

(31) \begin{equation} \eta_u = P \left( \eta_{\mathrm{des}} - f(v) \right) + f(v) \end{equation}

In Equation (31) P is the proportional gain, $\eta_{\mathrm{des}}$ is the desired distribution parameter, $f(v)$ is a function which takes the current mean speed as an input, and calculates the corresponding parameter according to the steady-state relationship and $\eta_u$ is the control output.

Speeds are initialised for $\eta = 45$ . Apart from the controller and the initial speed, the configuration is the same as in the previous cases. Agents’ initial position and heading are assigned according to a uniform distribution. Information on conflicts, states, speeds, conflict rate, throughput, traversal time, number of transitions and distribution distance is collected. Simulations were repeated to account for their stochastic nature.