On the global dynamics of connected vehicle systems

In this contribution, heterogeneous connected vehicle systems, that include human-driven as well as connected automated vehicles, are investigated. The reaction time delay of human drivers as well as the communication and actuation delays of connected automated vehicles is incorporated in the model. Saturations due to the speed limit and limited acceleration capabilities of the vehicles are also taken into account. The arising nonlinear delayed system is studied using analytical and numerical bifurcation analysis. Stability analysis is used to identify regions in parameter space where oscillations arise due to loss of linear stability of the equilibrium. Moreover, with the help of numerical continuation, bistability between the equilibrium and oscillations is shown to exist due to the presence of an isola. It is demonstrated that utilizing long-range wireless vehicle-to-vehicle communication the connected automated vehicle is able to eliminate the oscillations and keep the traffic flow smooth.


Introduction
Road transportation is expected to go through a huge transformation in the upcoming decades with the introduction of connected and automated vehicles. While such technologies mainly aim at improving passenger comfort and safety, they may also be used to improve the overall traffic flow. For example, by controlling the motion of automated vehicles may lead to higher traffic throughput and improved fuel economy [1][2][3][4][5][6][7]. Most prior research on the topic assumed 100% penetration rate of automation as well as vehicle-to-vehicle (V2V) connectivity. While the latter one is expected to rise rapidly in the next few years, the former one will likely follow in a slower pace due to the high cost of these vehicles. Thus, there is a need to understand the dynamics of mixed scenarios where automated vehicles are mixed into the flow of human-driven traffic. Initial research on mixed scenarios typically considered linear dynamics [8][9][10], while recent experimental results show that nonlinearities may play an important role in these systems [11][12][13].
In this paper, we start with modeling the longitudinal dynamics of human-driven vehicles. Then, a connected cruise control (CCC) algorithm is proposed to regulate the longitudinal dynamics of connected automated vehicles (CAVs) which utilize motion information from multiple human-driven vehicles ahead (within and beyond the line of sight). In particular, our goal is to design CCC controllers for CAVs so that they can mitigate traffic waves at the linear and nonlin- We adopt the strategy of "do not look ahead of another CAV" leading to the formation of a cascade of CCCs, as depicted in Fig. 1. This allows for a modular and scalable design of connected vehicle systems [14]. When modeling the vehicles, we incorporate time delays to represent the human reaction time and the communication and actuation delays of connected automated vehicles. Moreover, we take into account the speed limits and limited acceleration capabilities of vehicles, which introduce nonlinearities into the modeling equations.
In order to select the control gains of the connected automated vehicles, linear and nonlinear analyses are carried out while the vehicles placed on a circular road. First, linear stability analysis is presented to understand the effects of the connected automated vehicles on the stability of the uniform flow. With the help of the instability gradient method [15], we construct robust stability diagrams in the plane of the control parameters. We demonstrate that using beyond-line-of-sight information (obtained via V2V communication) can improve linear stability. Then, nonlinear analysis is carried out to have an insight into the influence of limited acceleration capabilities of vehicles. We use numerical continuation techniques [16,17], to draw the skeleton of the traffic dynamics as a three-dimensional wire frame representation. For the first time, we demonstrate that acceleration saturation leads to new dynamic behavior, where bistable regions appear through the presence of isolas. In this scenario, no linear stability loss occurs, but larger perturbations may still trigger traffic waves. Then, the effects of connected automated vehicles are investigated at the nonlinear level. We demonstrate that when using information only from the vehicle immediately ahead, an automated vehicle cannot eliminate the isolas. However, including additional information from beyond-line-of-sight, a connected automated vehicle can make the uniform flow globally stable.
The paper is organized as follows: in Sect. 2, we provide the governing equations of heterogeneous connected vehicle systems (including human-driven and automated vehicles) and explain the sources of the nonlinearities. Section 3 gives the linear stability analysis for the proposed model. A case study is conducted in Sect. 4 with a connected automated and two humandriven vehicles with detailed linear stability analysis. Then, in Sect. 5, we analyze the nonlinear effects of the acceleration saturation using numerical continuation. Finally, in Sect. 6, we conclude our results and discuss future directions.

Governing equations
In this section, we model the car-following behavior and the corresponding governing equations of heterogeneous connected vehicle systems. We consider that human drivers respond to the motion of the vehicle immediately ahead, while the connected automated vehicle's connected cruise control algorithm can utilize motion information from up to k cars ahead. We place N vehicles on a ring of length L such that we have M automated and N − M human-driven vehicles. The distribution of the human-driven vehicles and the connected automated vehicles are described with the help of the indices i h ∈ H and i a ∈ A, such that Following the models in [18], the governing equations can be given as: where the dot stands for differentiation with respect to time t and we use the abbreviated notation v ϑ i = v i (t − ϑ) to indicate delays. Also, v i denotes the velocity of vehicle i and h i denotes the distance headway, that is, the bumper-to-bumper distance between the vehicle i and its predecessor. This can be calculated as where x i denotes the position of the rear bumper of vehicle i and l i denotes the length of vehicle i such that N and v N + j = v j due to periodicity imposed by the ring. In case of a human-driven vehicle (i = i h ), the sum of the reaction time delay and the actuation delay is denoted by τ , while for a connected automated vehicle (i = i a ), the sum of the communication and actuation delay is denoted by σ .
The parameters α and β j are the control parameters of the connected automated vehicle, and α h and β h are the control parameters of the human-driven vehicles. For the sake of simplicity, all human-driven vehicles are considered to be identical and so the automated ones. The β parameters are directly related to the velocity differences, while the α parameters are related to the optimal velocity function V (h) (also called range policy), that satisfies the following properties: 1. V (h) is continuous, smooth and monotonically increasing (the more sparse the traffic is, the faster the vehicles want to travel).
cles intend to travel with the maximum speed, also called free-flow speed) These properties can be summarized as: where F(h) is strictly monotonically increasing and satisfies F(h st ) = 0 and F(h go ) = v max . In this paper, we present the results using the range policy function that is smooth at h st and h go ; see Fig. 2a. For simplicity, we use the same V (h) function for human-driven and automated vehicles. Saturations due to limited acceleration capabilities of vehicles are also taken into account by the function where a min and a max are the lower and upper saturation bounds, respectively. However, the nonsmoothness at a min and a max may lead to computational challenges. In order to avoid these, we smooth the saturation function by introducing a smoothness range, where the different sections defined in Eq. (5) are connected to each others while Additional parameters can be found in Table 1.

Linear stability analysis
In this section, we present the mathematical background of the linear stability analysis of the connected vehicle system. Setting the left-hand side of Eq. (1) to be zero, we obtain the unique equilibrium: where vehicles travel with constant speed while keeping constant distance. Using the vector the equilibrium can be given as Substituting the perturbation y(t) = x(t) − x * into Eq. (1), using a Taylor series expansion and dropping the higher-order terms, we obtain the variational systeṁ The 2N × 2N coefficient matrix in front of the nondelayed term is given by where I stands for the N × N identity matrix and S is an N × N circulant matrix, defined as The 2N × 2N coefficient matrices P and R in front of the delayed terms correspond to the human-driven vehicles and the automated vehicles, respectively. In particular, we have where κ = V (h * ) is the gradient of the range policy function at the equilibrium. Also note that at x * , the saturation function gives f (0) = 1. The N × N indicator matrix H contains zero elements except those diagonal terms which correspond to human-driven vehicles, that is, Moreover, we have where (I − H) corresponds to the automated vehicles. In order to determine the stability of the equilibrium, we utilize the D-subdivision method [19][20][21]. Substituting the trial function y(t) = c e λt into Eq. (10), the corresponding characteristic equation becomes The real and the imaginary parts of the characteristic roots are computed using the multi-dimensional bisection method (MDBM), which can find all the roots within a selected domain in the complex plane when using a sufficiently dense initial mesh [22].
In case of large number of vehicles, solving Eq. (16) can be time consuming and sensitive for rounding errors. However, in case of special patterns of distribution of the automated vehicles, where every (K = N /M)th vehicle is automated, Eq. (16) takes the analytical form Here, M is the number of connected automated cars, each monitoring the motion of K − 1 human-driven vehicles ahead, and the coefficients read as

Stability analysis of a simple connected vehicle system
In this section, the stability analysis is applied to the simplest case of CCCs with 3 cars on a ring, where the first vehicle is automated, while the second and third cars are human driven (N = 3, M = 1); see Fig. 3.
We analyze the effects of the average headway h * and the design parameters of the autonomous vehicle α, β 1 and β 2 , while the other parameters, such as the parameters of the human drivers and the time delays, are kept fixed. The applied parameters together with reaction time delays are presented in Table 1, see the related parameter identification techniques in [10,13].
The corresponding coefficient matrices in the linearized system (10) read as follows: cf. (11) and (14), while the coefficient matrices of the delayed terms are cf. (13) and (15). As a first step, we analyze the influence of an automated car that does not utilize beyond-line-of-sight information (β 2 = 0) on the traffic dynamics by drawing the stability chart in the plane of h * and α for different β 1 parameters. In the stability charts in Fig. 4, the gray areas identify regions in parameter space where oscillations arise due to loss of linear stability of the equilibrium. The coloring of the boundaries shows the frequency of the Hopf bifurcation in the corresponding nonlinear system. The darkest blue corresponds to ω cr = 0, that is related to a fold bifurcation. All the charts show stable range only between h * ∈ [5,55] m. Outside of this domain, the saturation in the range corresponds to the frequency of arising oscillations when stability is lost. Note that the parameters α, β 1 and β 2 correspond to the automated vehicle, while parameters of the human drivers are given in Table 1. (Color figure online)   Fig. 5 The effect of the β 1 parameter for h * = 30 m. The white areas are robustly stable for any average headway h * . In the gray area, a stable or an unstable equilibrium may exist depending on the average headway h * policy function eliminates the headway control so that equilibrium in (7) does not exist. Furthermore, negative α values always lead to unstable motion due to positive feedback.
For the leftmost panel in Fig. 4 (at β 1 = 0), large unstable range in the shape of a "droplet" can be found. This is centered around h * = 30 m, which corresponds to the largest derivative of the range policy function (κ = 0.943 1/s). Comparing the different panels in Fig. 4, the unstable droplet is shrinking when increasing parameter β 1 . However, the droplet is also shifting downward and merging with the unstable region α < 0. For β 1 > 1.3 1/s, the droplet is completely vanished and a large stable area forms. However, such large control parameter may require large acceleration that may not be possible under limited acceleration capabil-ities. Notice that the upper unstable area is also shifting downward as β 1 increases, but only causes instabilities for unrealistically large α values.
We prefer those robust control parameters for which the system is stable for any available headway h * ∈ [5, 55] m. The corresponding stability boundaries can be found by following the edges of the droplets where the tangent is horizontal (at h * = 30 m in our case) [23]. That is, if we change parameter h * , the real parts of the critical characteristic roots remain zero (Reλ = 0). This condition can be described as follows: while setting the critical value λ = iω cr . This can be rewritten as using implicit differentiation of the characteristic equation Eq. (16). Thus, the robust parameter boundaries can be calculated by solving the real and imaginary parts of Eqs. (16) and (23) together while setting λ = iω cr . In the 4-dimensional parameter space (h * , β 1 , α, ω cr ), Eqs. (16) and (23) give us a 1-dimensional object (curve). Still multi-dimensional bisection method can be used to find the corresponding solutions. Figure 5 shows the corresponding robust stability boundaries against the variation in the average headway h * in the plane of β 1 and α control parameters, while color indicates the critical frequency ω cr . Two robustly stable  Robust stability diagram in the plane of parameters β 1 and β 2 . The same notation is used as in Fig. 5 except that color represents the critical values of α domains can be observed. However, the lower one is not practical as such small α values may not be safe.  Fig. 6. It can be seen that the increase of the β 2 parameter improves stability. The unstable droplets are shrinking but not moving downward. The best option to improve linear stability is to increase both β 1 and β 2 together. Even for relatively small parameters, the droplet vanishes leading to robustness against changes in α and h * . Robust stability can be obtained by solving Eq. (16) together with while considering λ = iω cr . Note that Eqs. (16) and (24) consist of four equations, and in the 5-dimensional parameter space (h, α, β 1 , β 2 , ω cr ) this still results in a curve. Figure 7 shows the stable parameter domain in the (β 1 , β 2 )-plane, where the system is linearly stable for all relevant h * and α parameters values. In the gray area, the system may be stable or unstable for different h * and α parameter combinations. At the robust stability boundary, the unstable droplet completely disappears.
It should be noted that in this test case, the unstable droplet disappears always at fix h * = 30 m but for different α values, and the stability boundary is colored according to these α values. Note that the upper unstable area still exists, but only leads to instabilities for unrealistically large α values.

Nonlinear analysis
Using the above-introduced linear stability charts, it is possible to select the control parameters for the automated vehicle. However, these charts are valid only when small perturbations are introduced around the equilibrium traffic flow. In order to ensure stability for large perturbations, nonlinearities have to be taken into account when analyzing the global stability properties. To map out the global dynamics, we use numerical continuation, in particular, the package DDE-Biftool [16,17]. This allows us to follow branches of equilibria and periodic solutions while varying parameters and can provide information about the stability and bifurcations of these solutions. In our system, the source of the nonlinearities is the saturations due to range policy and limited acceleration capabilities (see Fig. 2).
First, the criticality of the Hopf bifurcation is analyzed. In case of supercritical/subcritical bifurcation, a branch of stable/unstable periodic orbits emanates from the Hopf points. Figure 8 presents the bifurcation diagram along a horizontal cross section at α = 0.6 1/s in the panel of Fig. 6 for β 1 = 0.3 1/s, β 2 = 0.15 1/s (see dotted line in Fig. 6).
The peak-to-peak amplitude of the velocity oscillations of the automated vehicle v amp 1 is plotted as a function of the bifurcation parameter h * . Note that for the equilibrium (7), this amplitude is zero corresponding to the horizontal axis. In this case, the equilibrium is unstable for h * ∈ [24.44, 35.56] m and the Hopf bifurcations, denoted by H 1 and H 2 , are supercritical. Continuous and dashed green curves correspond to the model with and without acceleration saturation, respectively. It can be seen that taking into account the saturation of the acceleration creates no difference close to the bifurcation points where the amplitude is small. However, after reaching the acceleration limit, the amplitude of the oscillation is smaller for the case with acceleration saturation. For h * = 30 m, the oscillations of the automated vehicle are depicted in Fig. 9 as phase-plane plots and as time profiles. Again the cases with and without acceleration saturation are distinguished as solid and dashed curves. Note that the amplitude of the oscillations may be different for different cars. Which vehicle reaches the acceleration saturation limit depends on selected parameters, but the amplitude of oscillations decreases for all vehicles due to the saturation. Furthermore, there exist scenarios, where the acceleration is saturated at both upper and lower boundaries.
In order to get a clear view of the global dynamics, the bifurcation diagram is extended along the control parameter α resulting in a three-dimensional wire frame representation of the periodic orbits shown in Fig. 10, where the peak-to-peak amplitude v amp 1 is plotted as a function of parameters h * and α. The three-dimensional isometric view (Fig. 10a) and its different projections (front, top, right) into parameter planes ( Fig. 10b-d) are presented to help the visualization of the structure.
We generate these figures as follows. First, we perform similar continuation of the periodic orbits as in Fig. 8 for numerous α values (corresponding to horizontal lines in Fig. 10c, d). Second, we fix the parameter h * and make similar continuation along α (vertical lines in Fig. 10b, c). Finally, in order to explore the full set of periodic branches, we continue from existing periodic orbits while varying both bifurcation parameters. Note that these steps are necessary to discover global dynamics leading to the bended tube-like shape in Fig. 10a, where each cross section along direction h * shows an isola for larger α values. Stable and unstable periodic orbits are colored with green and red curves, respec- tion in DDE-Biftool, and they bound the bistable area that is highlighted as a dark gray domain in Fig. 10e, as an extension of the linear stability diagram.  Fig. 11a, similar 3D-bifurcation diagram is visualized, but for a positive β 2 value (when the automated vehicle can have access for information beyond-lineof-sight). In this case, the bended tube-like shape is separated from the periodic orbits arising from equilibria. In this way, a set of isolas appear that are not connected to periodic branches arising from equilibria (see Fig. 11a, b and d).
The appearance of isola results in a bistable solution in two different scenarios: (i) two stable periodic orbits coexist or (ii) a stable equilibrium solution and a periodic orbit coexist. Scenario (i) can be found only above the linearly unstable equilibrium solution, which is already a practically unfavorable regime. In order to illustrate this case, point A is selected in Fig. 10c and e and oscillations are plotted for stable (green curves) and unstable (red curve) periodic orbits in the phase plane in Fig. 12 where the unstable equilibrium is denoted by red cross. In scenario (ii), the system may approach the equilibrium or the a periodic orbit depending on the initial conditions. This case is illustrated for parameter point B, noted in Fig. 11c and e. The corresponding phase portraits are shown in Fig. 13a and b where the stable equilibrium is denoted by green cross.
The role of the unstable oscillation is that it "separates" the basins of attractions of the two stable states (in the infinite dimensional state space). In other words, it determines the threshold for the strength of perturbations, which is necessary to "transfer" the systems from uniform flow to stop-and-go motion. Over-thethreshold perturbations may trigger stop-and-go traffic jams even when the uniform flow is linearly stable. Again, the bistable region is denoted by a dark gray area in Fig. 11e.
The final goal of parameter analysis is to ensure globally stable behavior of the traffic flow, that is, to extend the linear stability charts in Fig. 6 with the bistable regions as in Figs. 10e and 11e. This way, we can represent the influence of all the parameters and the global dynamics. This is summarized in Fig. 14, where only parameters within the white globally stable domains are favorable. Similar to linearly unstable droplet, the bistable domain shrinks when increasing β 1 and β 2 . However, notice that the β 2 parameter has a stronger effect on the reduction in the unsafe bistable domain.
It should be noted that without utilizing information from beyond-line-of-sight (β 2 = 0), neither the droplet nor the bistability can be eliminated within reasonable control parameters. For this case, we need very large β 1 to push the bistable domain up to the upper linearly unstable area. On the other hand, utilizing beyond-lineof-sight information (β 2 > 0), moderate parameters can be selected. During the design of the control, with the knowledge of the parameter uncertainties, we can select robust control parameters from the globally stable parameter domain (e.g., α = 0.5 1/s, β 1 = 0.3 1/s, β 2 = 0.3 1/s). In other words, these results allow us to identify those parameters where heterogeneous connected vehicle system is efficient and still robustly stable. These parameters were used in experimental validation of connected automated vehicle among humandriven vehicles in [12,13].

Conclusion
In this paper, heterogeneous connected vehicle systems were investigated, which include human-driven as well as connected automated vehicles. The reaction time delay of human drivers as well as the communication and actuation delays of connected automated vehicles was taken into consideration. Saturations due to the speed limit and limited acceleration capabilities of the vehicles were also incorporated. The arising nonlinear delayed system was studied using analytical and numerical bifurcation analyses. Linear stability analysis was used to identify regions in parameter space, where oscillations arose due to loss of stability of the equilibrium. Robust stability charts were drawn that can be used to ensure linear stability independent of the average headway (that changes as traffic becomes more dense) and the gain used for distance control (which is often chosen according to safety considerations). Moreover, with the help of numerical continua-tion, we determined the bistable regions, where smooth traffic flow coexists with congested states and sufficiently large perturbations may trigger traffic waves. We demonstrated that utilizing beyond-line-of-sight information (that is obtained via V2V communication), a connected automated vehicle is able to completely eliminate traffic waves at both the linear and nonlinear levels.