Nonlinear effects of saturation in the car-following model

The goal of this paper is to provide insight about the effect of acceleration saturation in the car-following model. In this contribution, we consider a heterogeneous, mixed-traffic scenario which contains both human-driven and autonomous vehicles subjected to time delays. Corresponding stability charts are provided from which one can tune the control parameters of the automated vehicles to achieve smooth traffic flow. By taking into account the acceleration saturation, it modifies the global behaviour of the system and reduces the range of the optimal technological parameters. On a demonstrative example, we highlight the complex dynamical phenomenon induced by the saturation and we attempt to connect these nonlinear investigations to the engineering practice and point out their relevance.


Introduction
Recent years have seen a major change in the transport and mobility industry. The ever-growing technological knowledge, skills and the spread of digitalisation have brought a huge potential for development. In addition, it is influenced by economic stimulus instruments and regulations. The environmental awareness, economic benefits and the need for safer travel lead to the penetration of more advanced driving aid technologies, connected and automated vehicles. Contrary to earlier expectations, realising fully automated transport is still decades away [1]. For this reason, investigating mixed-traffic situations, when automated vehicles AVs and human-driven vehicles HVs drive together, is essential [2][3][4].
Driving aid systems, such as cruise control CC, date back for decades and it is still in the interest of engineers. Nowadays the adaptive cruise control ACC [5,6] is widely available for consumer vehicles. These vehicles are equipped with various sensors, such as cameras and radar, and make decisions based on the collected data. More recently, the connected cruise control CCC [1] and connected adaptive cruise control CACC [7] is under development which may be common in the near future. It can utilise vehicle-to-vehicle V2V [8] communication or vehicle-to-everything V2X [9] communication in a network with other road participants to obtain information even beyond the line of sight . Theoretical controller designs [10,11] and experiments have already demonstrated [2,[12][13][14] that even a small number of AVs in a mixed-traffic situation can reduce traffic congestion and improve traffic flow. This technology not only helps to reduce journey times and fuel consumption [15,16], but also filters human error realising safer traffic, when safety guaranteed control is applied [4,[17][18][19][20]. All these factors may lead to significant economic and life quality improvements at the social and economical scale. The trend of transformation appears to be accelerating and it is likely that in the near future there will be a major breakthrough in these areas. However, these systems are based on dynamical models which require in-depth research and deep mathematical analysis. A comprehensive literature review of the CCC can be found in [21] and all the references therein.
The so-called car-following model was first introduced to understand highway traffic and to model the attitude of human drivers [22,23]. Later it has become a standard model for CC and connected autonomous vehicles CAVs [3]. In this model, time delay is usually considered with the intention to describe the traffic behaviour accurately [24]. In reality, delay appears due to reaction time for the HVs [25], data processing for the AVs and data transmission and communication for CAVs [26,27]. The car-following model is widely applied in the literature with the concept of ring road configuration, where N vehicles circulate on a ring [22,23]. If N → ∞ or practically N is a big enough number, this model can describe the behaviour of a straight road traffic flow. Comprehensive literature about ring configuration can be found in [28] and the related field tests and measurements in [13,29,30] which validate the theoretical findings. Furthermore, stability analysis considering the penetration ratio of AV on ring road configuration is analysed in [31].
Due to physical constraints, saturation is naturally involved in most of the engineering applications. It is often considered in the control law [32,33] for more precisely modelling the system (e.g. power limit of actuators and engines). Since saturation is highly nonlinear (in some cases it is modelled as a non-smooth effect), it is frequently neglected in the analytical analysis and only considered in numerical simulations. In the car-following model, acceleration saturation (or input saturation) constraints the performance [34][35][36]. It characterises more realistically the acceleration and deceleration capabilities of the vehicles that are the consequence of the capacity limit of the engines, the quality of the tyres and the condition of the road surface.
In our previous work [37], we investigated a system of CAVs where saturation was included in the model. Although the purpose of that research was not to expose the effect of saturation, it was concluded that badly chosen parameters resulted in a so-called isola causing bistability, where multiple solutions coexist, which could be improved by connectivity. In this case, depending on the amount of perturbation (braking) either smooth traffic flow or congestion may develop. Isola is a special, isolated and closed solution to nonlinear problems, which is detailed in the scientific literature [38][39][40] and appears in a variety of disciplines [41,42]. Detecting and finding an isola in a system is not straightforward and requires advanced analytical [43,44] or numerical techniques [45][46][47][48]. On this basis, it is substantial for engineers in practice to be aware of this phenomenon and to have knowledge about isola, and understand its nonlinear dynamics.
Although our recent work [37] shares some of the ideas presented in this paper, here, we establish a comprehensive in-depth study that is not covered by previous works, including the acceleration saturation and putting more emphasis on its effect on the complex dynamical behaviour. To the best of our knowledge, despite that the application of saturation in control theory is quite widespread [49][50][51][52], its global effect on the dynamics of CAVs has not yet been fully addressed in detail.
Therefore, in this paper we intend to extend the literature by studying the saturation effect on the local and global behaviour of the car-following model. In Sect. 2, the car-following model is revisited and the governing equations are presented. Section 3 details the theory and derivation of the equations for the stability analysis. First, in Sect. 3.1 the linear stability analysis is presented; then, Sect. 3.2 details the theory of the nonlinear analysis, whereas in Sect. 3.3 the possible numerical and analytical methods are collected. In Sect. 4, a demonstrative example is presented in details. In Sect. 4.1, the used parameters are collected and the particular models shown, Sect. 4.2 provides the linear stability analysis and the corresponding linear stability charts, Sect. 4.3 shows the nonlinear calculation extending the stability diagrams by the global stability information and Sect. 4.4 details the evolution of isola and its relation to acceleration saturation. Finally, in Sect. 5 a conclusions is drawn.

Model and equations
This section revisits the standard car-following model placed on a ring road which is commonly used in the literature to model highway traffic shown in Fig. 1. The system can be generally given as a nonlinear delay differential equations DDEs in the following form: with the state vector where N is the number of vehicles, τ i is the time delay of the ith car, v i is the velocity and h i is the headway (following distance), the bumper-to-bumper distance, between car ith and its predecessor i + 1th. The equation for each ith vehicle can be given aṡ where i = 1, . . . , (N − 1). Note that the headway of the N th vehicle can be calculated from all the other headways and the length of the ring L such as With this, the acceleration of the N th vehicle can be given aṡ The parameters α i and β i of the ith car are the control parameters of the vehicles. The function V i (h) is the so-called range policy, which describes the desired speed related to the headway for each vehicle (see, for example, in Fig. 2 panel a). It characterises the driving behaviour of human drivers or autonomous vehicles; for example, the more sparse the traffic is, the faster the vehicles want to travel. It can be described generally by the nonlinear continuous piecewise function as The lower limit h st is the critical headway value under which the vehicle tends to stop, and the upper limit h go is the headway value over which the vehicle wants to travel by the maximum possible speed which is the speed limit on the road. Between these limits, the characteristic is described by the F(h) function, which can be different to each vehicle. Engineers are free to use whatever relation is optimal. In the literature, different characteristics are often used such as linear, cubic and tangent functions. In favour of modelling more realistically the acceleration and deceleration capabilities of the vehicles, the acceleration saturation is taken into account in the model by the piecewise continuous function sat, as shown in Fig. 2 panel b. It can be described by the following formula where a min and a max describe the lower and upper acceleration limits, respectively. Note that more sophis- ticated saturation function can be found [36] which considers the power train capabilities, however, those are more typical for heavy-duty vehicles and is out of the scope of this paper.

Stability analysis
In this section, we briefly describe the investigation methods for the local and global behaviour of the carfollowing system by means of linear and nonlinear stability analysis. The equilibrium is the state when all the vehicles travel with the same velocity keeping a constant headway distance determined by the range policies. In the case of unstable equilibrium, even the smallest perturbations (braking) get amplified resulting in stop-and-go waves (of velocity and headway) travelling backwards along the ring [29,53], whereas, if the stable equilibrium is perturbed after the transient vibrations, the smooth traffic flow will be restored, in which the system is intended to operate, since it is ideal for safety reasons, fuel consumption and comfort.

Linear stability analysis
In the present section, the linear stability analysis is explained. This calculation is crucial to find the stable and unstable equilibria in the parameter space, hence choosing the desired control parameters. The system has its equilibrium when all the variables are constant, which can be determined by setẋ(t) = 0 in (1) as where the equilibrium vector x * ∈ R 2N −1 is constructed such as Solving (9), the so-called uniform equilibrium flow can be given by In the case of equilibrium, all the vehicles travel with the same velocity v * i while keeping different headways h * i due to the different range policies V i . To analyse the local behaviour of the equilibrium, the system is linearised around it. It is accomplished by introducing a small perturbation x(t) = x(t) − x * around the equilibrium that leads to Substituting these into the (1) leads tȯ (13) which is linearised by expanding it into Taylor series and cancelling the higher-order terms Here, L and R i are 2N − 1 × 2N − 1 Jacobian matrices corresponding to non-delayed and the ith time delay terms, respectively, defined as and Note that during the linearisation, the chain rule must be applied for the derivation as therefore one can omit the saturation function in (3), due to the fact that sat (0) = 1 around the equilibrium. This implies that the saturation does not affect linear stability properties. Recently, a similar observation has been found in [33]. Finally, the so-called variational system formṡ The stability is determined by introducing the nontrivial trial functions x(t) = ce λt , where λ is called the characteristic exponent and vector c ∈ R 2N −1 . Substituting this trial function into (18) where I is the 2N − 1 × 2N − 1 identity matrix. The e λt > 0 term is always fulfilled, and we are interested in non-trivial solution c = 0, which determines the so-called characteristic equation in the following form This is a transcendent equation since λ has infinitely many solutions which leads to computational challenges. In particular, the infinitely many characteristic exponents can be approximated by finite ones applying discretisation techniques, such as the so-called semidiscretisation method [54]. The equilibrium is linearly stable if the rightmost characteristic exponent has negative real part. One elegant way to produce the stability chart and draw the exact stability curve of systems in the form of (20) is the so-called D-subdivision method [55,56]. The curves are determined where the number of unstable exponents changes. That is substituting λ crit = iω into the characteristic equation (20) and separate into real and imaginary parts. To compute these curves efficiently, the multi-dimensional bisection method (MDBM) can be applied [57]. Another but not so efficient way is to map the parameter space in discrete points by using brute force calculation. Finally, one can use numerical continuation methods [46,58,59] to calculate the stability boundaries, which are more sophisticated and accurate.
Note that linear stability theory implies that in the case of an unstable system the amplitude of the arising oscillations goes to infinity. However, in the case of large amplitude vibrations, nonlinear behaviour and acceleration saturation come to play a role. Therefore, it is important to consider the nonlinearities in the system.

Theoretical background
The linear stability analysis shows the behaviour of the system in the vicinity of the equilibrium. In reality, sudden acceleration, braking, lane change, etc. can happen introducing large perturbations. Since the original DDE system in (1) is nonlinear, the global behaviour may deviate from the linear one causing complicated dynamics. To have information about the system further from the equilibrium, it is necessary to perform nonlinear analysis.
The critical roots of the characteristic (20) often form complex conjugate root pairs in time delay systems [56]. At the stability boundary, a root pair loses its stability simultaneously: Re λ 1,2 = 0 and Im λ 1 = Im λ 2 . It is a dynamical stability loss, called Hopf bifurcation, from which a periodic orbit arises [60,61]. A periodic solution implies that the system has a steadystate oscillating solution. If a stable periodic orbit is over the unstable equilibrium, then it is called supercritical Hopf bifurcation as illustrated in Fig. 3 in panel a. It is undesirable to operate in this state due to the unstable equilibrium. Its counterpart is the so-called subcritical Hopf bifurcation, when an unstable periodic orbit rises over the stable equilibrium, as demonstrated in Fig. 3 panel b. This case may seem safe to operate, since an oscillation caused by a small perturbation (e.g. gentle braking) settles down. However, if perturbation pushes the state outside of the unstable periodic orbit, the amplitude of the oscillation will keep increasing. Note that theoretically it goes to infinite, but another stable periodic orbit might be over the unstable one, limiting the vibrations. Taking this into account, in engineering applications the subcritical bifurcation can be generally considered more dangerous than a supercritical one.
In Sect. 3.1, it was already shown that the saturation does not affect the linear stability boundaries. Presumably, the acceleration saturation does not affect the criticality of the bifurcation either, since only the larger amplitudes are limited and the periodic orbit arises starting with an infinitely small amplitude oscillation. This was also proved by analytical calculation in a dynamically analogous system of a delayed position regulation in [33]. In that system, the saturation occurred in the actuator due to the limit of maximum torque and force which is feasible in reality. In the case of the car-following model, the same characteristics are to be expected. There are several methods available for performing nonlinear analysis, each having certain advantages. One of the most straightforward is the brute force numerical simulation. The differential equation system (3) can be directly solved with various initial conditions mapped in the parameter set of interest. This method does not require further mathematical calculus and it is easy to do in many software with built-in differential equation solvers. However, its disadvantages are being time-consuming and it is hard to draw a conclusion based on the resulting signals of the vibrations, especially at the stability boundaries. By applying the theory of Hopf bifurcation analysis for DDEs, one can obtain the second-order analytical approximation of the amplitude of the periodic orbit [60], which is illustrated by a dashed parabolic curve in Fig. 3. To do so, (13) has to be expanded into Taylor series up to the third order. Thereafter, via the infinitedimensional centre manifold reduction [33,[62][63][64][65][66], the infinite-dimensional system is reduced into two dimensions, called the direction of main (slow) dynamics which correspond to the complex conjugate root pairs at the Hopf bifurcation point. The other dimensions can be detached, since the vibrations corresponding to those settle down faster [43]; hence, it does not influences the steady-state solution.
To have a better understanding of the dynamics in the infinite-dimensional state space, the three-dimensional qualitative representation of the centre manifold reduc-tion is shown in Fig. 4 around the vicinity of a Hopf bifurcation point. The two slow directions are ξ 1 and ξ 2 corresponding to the long-term dynamics while ξ i illustrates the other infinitely many dimensions. By linearisation, the slowest spectral subspace is approximated by a tangent plane at the equilibrium spanned by the directions of ξ 1 and ξ 2 . Furthermore, the above referred nonlinear analytical calculation provides the secondorder approximation of the so-called centre manifold as a surface, on which the long-term dynamics can be characterised by periodic orbits.
In Fig. 4, typical behaviours are illustrated starting with different initial conditions for super-and subcritical cases. Panel a illustrates the supercritical case, where an unstable equilibrium (red) is in the origin and a stable periodic orbit (green) exists on the centre manifold. In case of small/large perturbation representing by setting the initial conditions close/far to the equilibrium, after the transient vibrations the trajectories settle on the surface of the centre manifold and converge towards the stable periodic orbit. By comparison, the subcritical case is illustrated in panel b of Fig. 4. Here, a stable equilibrium is in the origin and an unstable periodic orbit lies on the centre manifold. In this case, there exist a manifold, which divides the infinitedimensional state space and determines the basin of attraction, illustrated by red shaded surface. In the case of the blues trajectory, the initial condition is close to the equilibrium starting within the basin of attraction and converges to the equilibrium. As it was pointed out above, the subcritical case may be dangerous since it looks like a stable, well operating system. However, starting the initial condition outside of the region of attraction (yellow trajectory), the solution converges to the centre manifold, but it is repelled on it by the unstable periodic orbit and continues to oscillate with evergrowing amplitude. Note that the approximation of the basin of attraction for time delay systems is a challenging problem due to the infinite-dimensional state space and out of the scope of this work, although there exist some promising techniques in the literature [67][68][69][70].
In this paper, for further calculations, the DDE-BIFTOOL continuation package [46,59] is used to capture the global dynamical behaviour of the system. It allows to follow branches of equilibria and periodic orbits while changing system parameters and can determines stability and bifurcation information of these solutions.

The model and parameters
The above generally introduced theory is presented for a particular case study. Due to the complexity of the system, the model is kept as simple as possible while keeping the core elements of the dynamic behaviour. Therefore, N = 3 vehicles are used in Eqs. (3) and (4) to investigate the fundamental effects of acceleration saturation. Nevertheless, only a few vehicles on a ring can give an insight on how a larger and more realistic system will behave [28]. During this analysis, the first vehicle is considered as an autonomous vehicle and its control gain can be tuned α, but β a and τ a are fix. The other two are human-driven vehicles in which parameters, such as α h , β h and τ h , are assumed to be identical. For the range policy in (7), we apply the widely used cosine function [1] for all the three vehicles, which is differentiable infinitely many times and ensures C 1 continuity at the boundaries, defined as follows Again, for the sake of simplicity all three vehicles have the same range policies and saturation function. The acceleration saturation is taken into account in the model by the piecewise smoothed continuous function [37,71] sat(a)= where a min and a max describe the lower and upper acceleration limits, respectively, and the smoothness range is described by c. This smoothness is used to avoid computational challenges during the numerical calculations later. The system of equation reads as where the subscripts "a" and "h" stand for the automated and human vehicles, respectively. The headway of the third vehicle is given by the geometry of the ring arrangement The equilibrium of the system can be determined based on (11) which reads meaning that all the vehicles travel with the same headway and velocity. Here we fix all the parameters, which are listed in Table 1, and we investigate the system behaviour in the plane of the control parameter α a and headway h * .
First, linear stability analysis of this equilibrium is determined based on Sect. 3.1

Linear analysis
The corresponding characteristic equation of (23) based on (20) reads where the coefficient matrices based on Eqs. (15) By solving the characteristic equation, one can provide linear stability behaviour of the equilibrium, which is visualised by the stability diagram in Fig. 5 in the plane of α and h * . The grey area corresponds to linearly stable equilibrium, while the white domain in the middle indicates unstable behaviour. They are separated by the blue curve, which refers to the dynamical stability boundary, the branch of Hopf points. The stability diagram is symmetrical at h * = 30 m due to the

Fig. 5
Linear stability diagram of system in (23). In the parameter space, grey represents linearly stable equilibrium, while white shows unstable one. These are separated by Hopf bifurcations points, denoted by blue curve. Note that the unstable region is symmetric to h * = 30 m point symmetric range policy function in (21). It can be observed that for α < 0.25 and α > 0.6 the system is linearly stable for any headway, while for parameter domains h * < 25 and h * > 35 it is also stable for any α value. The middle region of the parameter space traffic becomes unstable. Note that as it is expected, the acceleration saturation does not affect the linear stability properties. In order to demonstrate the effect of the saturation, the system is numerically simulated for small and large perturbations with parameters h * = 30 m and α = 1 1/s (corresponding to point B in Fig. 5 where the system is linearly stable). The resulted acceleration time For small perturbation, the vibrations settle down to the stable equilibrium; however, it takes longer time due to the saturation. Finally, panel d shows a case, when the breaking is too harsh (large perturbation); then, the system gets into a state of stable periodic vibrations which does not settle down to the equilibrium. In this case, the vehicles behave in an extreme way, as they are accelerating and suddenly decelerating with saturation limits. This indicates that the global stability behaviour is different from the linear one through the appearance of bistable phenomenon when saturation is included in the model; hence, it is necessary to carry out nonlinear analysis.

Nonlinear analysis
The resulting linear stability diagram in Fig. 5 can be useful to choose the right control parameters for However, it provides only information about local behaviour, when the system is close to the uniform traffic flow. If it is exposed to large perturbations (e.g. large breaking in the traffic), it is necessary to analyse the global behaviour by means of nonlinear investigations. Note that nonlinearities appear in the system due to both the range policy (see (21)) and the applied acceleration saturation (see (22)).
With the DDE-BIFTOOL package, one can continue the periodic solutions along various parameters and can capture the global behaviour. As explained in Fig. 3, the periodic orbits are arising from the Hopf points. The bifurcation diagram in Fig. 7 showing the peak-to-peak amplitude of the velocity and acceleration is plotted for α = 0.5 1/s parameter (along the dashed black line in Fig. 5). The dashed green curve denotes the amplitude of the periodic orbit without acceleration saturation, while the continuous one is with saturation. In panel b, one can see that the acceleration amplitude is limited at a max −a min 2 = 1.5 m/s. This effect is also visible in the velocities amplitude in panel a as the magnitude of the velocity amplitude is significantly decreased. (The continuous curve is located lower than the dashed one.) Note that the acceleration saturation does not influence the criticality of the bifurcation points, however, it significantly modifies the amplitude of the arising periodic orbits. This is illustrated in Fig. 8 as well, plotting the time profiles (panel a, b) and phase portraits (panel c, d), for the same parameters corresponding to point A in , but in the presence of saturation, the oscillation is limited at a max and a min denoted with dashed grey lines. In the phase plane (panel c and d), the stable periodic orbits are located around the unstable equilibrium with and without saturation. In panel d, the effect of saturation can also be seen as the acceleration is limited. Based on these figure, saturation may seem beneficial causing less extreme oscillation. However, it does not explain the bistable phenomenon observed during the numerical simulations in Fig. 6. Therefore, it is worth further studying the dynamics and calculate the bifurcation diagrams along parameter α at h * = 30 m. Figure 9 contains the results of this investigation. In panel b, the periodic orbits are continued along α. Observe that in this case the bifurcation diagram and the dynamics become more complicated. Stable and unstable periodic orbits bend over the stable equilibrium resulting in bistable behaviour, where stable equilibrium and stable periodic orbit coexist. Note that this bistability is also observed by the numerical simulations in Fig. 6. Here, the branches of stable and unstable periodic orbits are separated by a fold bifurcation point, marked by grey cross. This fold points is continued in the plane of h * and α; hence, the above-introduced stability diagram is extended with the boundaries of the bistable region, which can be seen in panel a of Fig. 9.
Another representative illustration of the bistable behaviour is shown in Fig. 10 by phase portraits cor- responding to parameter point B in Fig. 9. Here, B 1 represents the stable equilibrium, B 2 is the unstable periodic orbit and B 3 shows the stable periodic orbit. It can be seen that both the stable and unstable periodic orbits are saturated and resulted from its consideration. According to the nonlinear analysis, the acceleration saturation significantly changed the global behaviour of the system. Therefore, it is worth investigating how the value of saturation limit affects the system.

Effect of the acceleration saturation
In this section, similar nonlinear analysis is carried out as it is presented in Sect  Fig. 11. The global stability diagram in panel a shows that the bistable area is decreased compared to the case in Fig. 9. Still, there is a stable branch of periodic orbits over the unstable equilibrium in the middle of the parameter space, but its amplitude is increased and its shape is slightly deformed. Furthermore, unstable and stable periodic orbits are also appeared over the linearly stable equilibrium leading to a smaller bistable domain. One can observe that these periodic orbits are not connected to those that emerging from the equilibrium. Thus, a separated branch of periodic orbits appear, which are also referred as isola. Notice that it is challenging to discover its presence even with sophisticated numerical continuation techniques since it cannot be continued from a Hopf point. It is a matter of luck to find it without aimed calculations using only brute force methods and numerical simulations. Our method to find the isola and to generate the diagrams are the following. Note that in the higherdimensional parameter space, the isola is still connected to the equilibrium along different parameters, hence, it is still possible to use numerical following techniques. First, a point is chosen along the periodic orbits corresponding to a max = 1 m/s in Fig. 9 (e.g. B 2 or B 3 ). Then, the parameter h * is fixed and the periodic orbit is continued along a max (considering the saturation as a parameter). The second step is to chose the Fig. 12 Bifurcation diagram at a section of an isola at α = 1.24 1/s, with a min = −6 m/s and a max = 3 m/s. Green represents stable, red unstable states. A stable and unstable periodic orbit attracts and repels the solution, respectively. The equilibrium is locally stable in the whole parameter set point with the desired saturation value (a max = 3 m/s) and continue the periodic orbit in the direction of h * or α, as intended. Finally, the bistable domain can also be characterised by means of the same technique as presented in Sect. 4.3.
The real danger of the isola can be understood better in Fig. 12 for parameters α = 1.24 1/s, a min = −6 m/s and a max = 3 m/s. Isolated stable and unstable periodic orbits are located over the locally stable equilibrium; (a) (b) Fig. 13 Evolution of isola, concluded from multiple simulations, is illustrated in the bifurcation diagram in panel a, showing as the saturation limit is increased the isola gets separated from the equilibrium. Panel b shows the global stability diagram, show-ing the effect of the saturation limit on the bistable regions as the bistable region decreases. The isola fold points are connected with a yellow curve hence, the parameter space along h * can be divided into 3 zones. Left and right domains are globally stable, where oscillations due to any perturbation will settle down to the equilibrium after the transient. However, the middle zone under the isola is bistable. A large enough perturbation will cause the system to jump to the stable periodic orbit with high-amplitude oscillation which results in traffic jams. Therefore, an isola can be considered as an even more dangerous state than a subcritical bifurcation, because it can cause similarly complex behaviour, but it is more challenging to find.
To analyse the evolution of the isola, the periodic orbits are calculated for several saturation values, as illustrated by the bifurcation diagram in the panel a of Fig. 13 for h * = 30 m. As the acceleration limit is increased, the velocity amplitude in the bifurcation diagram is increased as well. Note that once the saturation limit is larger than a max = 5 m/s 2 , the resulting periodic orbits coincide with the ones without acceleration saturation. In this case, even large amplitude vibrations do not reach the effect limit. Additionally, we observed that the isola and the corresponding bistable region disappear once a max = 5 is exceeded. The same phenomenon is illustrated in the global stability diagram in panel b. As the saturation limit increases, the bistable region shrinks leading to a safer system with larger globally stable regions. The lowest point of this bistable region (where the isola appears) can be found as the minimum of the fold points for each saturation value. We call these as the isola endpoints, connected by yellow curve in Fig. 13 panel b. We plot Fig. 14 Relation of isola fold points and saturation limit. The isola endpoints are marked with the yellow curve. For a given acceleration saturation limit value, the corresponding α control parameter value can be selected to avoid bistability these points in the plane of α and a max , as illustrated in Fig. 14. Over the yellow curve, the system is conditionally bistable since it also depends on h * parameter value. This graph can be use to chose the right α control parameter considering the capabilities of the particular autonomous vehicle to avoid the bistable behaviour and ensure smooth traffic flow.

Conclusion
In this paper, we have investigated the effect of the acceleration saturation on the local and the global dynamical behaviour of the car-following model subjected to time delay. First, we have presented stability analysis for general heterogeneous car-following system. Then further analysis has been expanded by a demonstrating case study consisting a mixed-traffic scenario of two human-driven vehicle and one automated vehicle.
We have shown that the saturation limit has no effect on the linear stability properties. However, we have broken down how the acceleration saturation affects the global dynamics by the appearance of separated isolas and bistable parameter zones, when the smooth traffic flow may occur together with congested traffic flow depending on perturbations that can trigger traffic jams. We have also explained how the presented method can be used to find these isolated periodic orbits. Then, we have extended the linear stability chart with these bistable parameter domains. With the introduced global stability diagrams, it is possible to choose the desired control parameters of the driving aid system and autonomous vehicles to ensure global stability, considering the acceleration limits of the particular vehicles. We have also provided how the bistable domains decrease by increasing the saturation limit (having more capable vehicles).
These analysis may open ways for efficient parameter optimisation strategies crossing the border towards reliable and safe autonomous vehicle design. Our future goals are to extend the system considering more vehicles and to perform experiments to support the theoretical results.
Funding Open access funding provided by Budapest University of Technology and Economics. The research reported in this paper is part of project no. BME-NVA-02, implemented with the support provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, and financed under the TKP2021 funding scheme.

Data availability
The authors declare that data supporting the findings of this study are available within the article.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Com-mons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.