On the moment dynamics of stochastically delayed linear control systems

In this article, the dynamics and stability of a linear system with stochastic delay and additive noise are investigated. It is assumed that the delay value is sampled periodically from a stationary distribution. A semi‐discretization technique is used to time‐discretize the system and derive the mean and second‐moment dynamics. These dynamics are used to obtain the stationary moments and the corresponding necessary and sufficient stability conditions. The application of the proposed method is illustrated through the analysis of the Hayes equation with stochastic delay and additive noise. The method is also applied to the control design of a connected automated vehicle. These examples illuminate the effects of stochastic delays on the robustness of dynamical systems.

vibrations. [5][6][7] When the system parameters (including the delays) are constant and no external excitation occurs, there are well-established methods to investigate the linear as well as the nonlinear dynamics of these systems. [8][9][10][11][12][13] With added noise, stability investigations become more challenging and the stationary motion of a system needs to be carefully characterized. [14][15][16][17] Moreover, the system parameters (including the delays) may also vary stochastically requiring sophisticated mathematical tools for stability analysis. [18][19][20] In this article, control systems are considered where both effects are present: the time delays vary stochastically while the system is excited by additive noise.
For example, in vehicular traffic, the driver reaction time typically varies stochastically, 21 while the additive noise comes from the other vehicles whose motion the driver needs to respond to. In network control systems, delays may vary stochastically due to packet drops, and in the meantime agents need to respond to noisy environment. 18,[22][23][24][25][26][27] In complex biological networks, like those within cells, external noise is ubiquitous and stochastic delays may be used to model a sequence of reactions. 28,29 In machining applications, in order to suppress chatter, one may stochastically vary the spindle speed, 30 which leads to stochastic delays in the corresponding mathematical model.
One of the most important performance measures for control systems is stability. For stochastic systems, there exist multiple stability definitions. For example, almost sure exponential stability characterizes the mean decay rate of the exponential envelope of the solution of the system. [31][32][33] However, this approach does not guarantee the existence of a stationary solution in the presence of additive noise since there can be a finite number of "bursts" as time evolves. A more conservative approach is to consider the first moment (or mean) and second moment (or mean square) stability. This allows one to determine whether the system has bounded stationary solutions or so-called stochastic coherence resonance 15,16 occurs due to additive noise. Thus, this approach can be utilized to characterize the robustness of the control system. 34 For stochastic systems without stochastic time delays, earlier works 14,15 investigated not only the stability of the system but also the effects of external noise near critical parameters. In particular, small perturbation of a critical parameter were utilized for small degree-of-freedom dynamical systems. The moment stability and stationary second moment of larger stochastic systems has been considered 17 while the time delays were assumed to be constant. Incorporating stochastic delay variations lead to switching systems, whose analysis require sophisticated mathematical approaches. For example, sufficient conditions for second moment stability for switching systems with multiplicative noise were derived. 35 Also, necessary and sufficient conditions were given in the presence switching delays. 4,20 However, the sensitivity against external excitations was not investigated. The stationary behavior was considered for a system with additive noise in optimal control designs 19,36 and a model predictive control design 37 as part of the cost function. However, the resulting controllers are optimal with respect to the prescribed cost function only without leaving any room for more general analysis that is needed for a general control design.
In this work, the stability and stationary solution of continuous time systems with stochastic time delays and additive noise are discussed, without requiring the independence of the two different stochastic excitations. The resulting stochastic dynamical system is approximated using semi-discretization, 4,17,38 and the corresponding discrete time stochastic map is determined analytically. This map is utilized to obtain the first and second moment dynamics, which in turn, allow us to derive necessary and sufficient conditions for the first and second moment stability and to calculate the stationary first and second moments. The effects of the stochastic delays on the moment dynamics are illustrated through a simple system, called the Hayes equation, where the moments are calculated along the delay-interval for different values of the system parameters (including the holding time). The obtained methods are also applied to the control design of a connected automated vehicle (CAV) subject to packet losses while responding to the noisy motion of its human-driven predecessor.
The article is organized as follows. In Section 2, the semi-discretization is applied to the stochastically delayed system with additive noise and the first and second moment maps are constructed. The computation of the stationary first and second moments of systems driven by white and delay-induced additive noise are discussed in Sections 3 and 4, respectively. In Section 5, some properties of the stationary moments are highlighted through the stochastically delayed Hayes equation, while in Section 6, a control design for a CAV is presented that is robust against stochastic delays. Finally, conclusions are drawn in Section 7.

MODELING AND DISCRETIZATION
In this article, linear systems with stochastic delay and additive noise are considered: where the dot indicates the derivative with respect to time t, x ∈ R d is the d-dimensional state vector, (t) represents a stochastically varying time-delay, A, B ∈ R d×d are the coefficient matrices, and (t) ∈ R d is an arbitrary stochastic noise process with an appropriate probability space. We emphasize that (t) and (t) are not required to be independent for the results presented in this section to be valid. While no particular assumptions are made for the additive stochastic process (t), a specific stochastic process is considered for the time delay (t) which is motivated by applications of wireless communication-based control systems. Namely, the delay process (t) is considered to have piecewise constant trajectories. In particular, the delay is assumed to stay constant for a holding time T before potentially taking on a new value from a finite set { 1 , 2 , … , J }, such that while remaining constant for each interval t ∈ (kT, (k + 1)T]. In order to investigate the stability and stationary solution of the delayed system in Equation (1), the continuous time dynamics can be approximated by a discrete one. This discrete approximation can be constructed using full discretization (such as the Euler method) or semi-discretization. 4,17,38 Since the latter has superior convergence properties, in this article this approach is utilized. The zeroth order semi-discretization is discussed in the main body of the article, while higher order semi-discretization is shown in Appendix A.
Equation (1) can be discretized to: where n ∈ N counts the time step under the time resolution Δt = T∕l; see Figure 2.
We choose l such that Δt < j , j = 1, … , J. The discretized delay is given by where ⌊.⌋ denotes the floor operation. Note that r(n) follows a similar stochastic process as (t), that is, where r j = ⌊ j Δt ⌋ , j = 1, … , J; cf. (2). The differential equation (3) can be solved analytically for one period of length Δt: By defining the augmented vector system (6) can be written in the compact form where the coefficient matrix and the disturbance vector are and In case of an invertible coefficient matrix A, the matrix R can be expressed as R = (e AΔt − I)A −1 B. Note that in (9), the block matrix R is in the (r(n) + 1) th block-column of matrix F(n). In other words, the position of the block matrix R in the first block-row of F(n) depends on the instantaneous value of the delay: if (t) = j in the time interval (kT, (k + 1)T], then r(n) = r j . In this case, F(n) can be substituted by F j and, based on (5), F(n) is IID and follows the probability distribution Note that the delay value (t) does not change during one holding period T = lΔt; see Figure 1. Therefore, defining the state vector z(k) = y(kl), for each holding periods and applying (8) iteratively, the system dynamics can be written as where G(k) = F(kl) l , and We emphasize that since the stochastic process (t) generates the process r(n) which generates F(n) and G(k), they all share the same probability distribution function that hase been defined in (2).

Moment dynamics
Here the general expressions are stated for the time evolution of the mean and second moment of system (13) for a general noise term (t). In the next section, the dynamics are discussed in more detail when (t) is a Gaussian white noise. First, note that the random variables G(k) and z(k) in (13) are independent. The matrix G(k) depends only on the delay value in the time interval kT ≤ t ≤ (k + 1)T with probability mass function given by (15), while the vector z(k) depends on the values of the delay (t) and the noise (t) in the time interval 0 ≤ t ≤ kT. Therefore, given the assumption that the delay switchings are IID, G(k) and z(k) are independent. Now by taking the expected value of both sides of Equation (13) and exploiting the independence of G(k) and z(k), the mean dynamics is obtained, where and E(.) denotes the expected value operator.
To obtain the second moment dynamics, the vector is defined, where ⊗ denotes the Kronecker product. Here the fact that E(z(k) ⊗ z(k)) = vec(E(z(k)z(k) ⊤ )) is used, where vec(.) is the vectorization operator that places the columns of a matrix below each other. The definition using the Kronecker product in (18) is suitable for the second moment dynamics analysis carried out in this article. By the properties of Kronecker product and the independence of G(k) and z(k), the second moment dynamics is derived, where The details of the derivation of the term c z (k) are shown in Appendix B.
In the next section, the specific forms of Equations (17) to (20) are derived, while considering an additive Gaussian white noise.

DYNAMICS AND STABILITY WITH GAUSSIAN WHITE NOISE
So far, no particular type of stochastic process was assumed for the noise (t) in (13). Thus, (17) to (20) are valid in general. For simplicity, in this section, (t) is considered as a Gaussian white noise process with the form where i , i = 1, … , d , are mutually independent Gaussian white noise processes with mean E( i (t)) = 0 and correlation function E( i (t) j (s)) = (t − s) ij . The vectors i ∈ R d , i = 1, … , d , denote the noise intensities. The noise process defined in (21) can be written in the compact form where the ith column of the matrix ∈ R d×d is i and the ith component of the vector (t) is i (t). Observe that E( (t)) = 0. Next, the behavior and stability of the mean dynamics (16) and the second moment dynamics (19) are investigated when (t) is given by (21).

Mean dynamics
In this section, it is illustrated that the mean dynamics (16) can be simplified to by showing that g(k) = 0 when the noise (t) is given by (21). According to (14) and (16), the additive vector is given as Note that F(kl) only depend on the delay (t) in the interval (kT, (k + 1)T] while f(kl + m), m = 0, 1, … , l − 1, depends only on the Gaussian white noise (t) in the interval (kT + mΔt, kT + (m + 1)Δt]. Since (t) and (t) are independent F(kl) and f(kl + m) are also independent. Thus, (24) becomes On the other hand with Gaussian white noise (10) yields where W t is the Wiener process vector corresponding to the white noise (t), namely The last equality in (26) follows from the zero mean property of the Itô integral. 39 Now according to the definition (9), we have and hence (25) yields g(k) = 0 and the mean dynamics are given by (23). From system (23), it is concluded that the mean z(k) converges to 0 as k → ∞, if and only if where (.) denotes the spectral radius, namely

Second moment dynamics
In this subsection, the second moment dynamics (19) is simplified to the following form and the disturbance term g is calculated explicitly. To do this, first the term c z (k) in (19) is considered. Note that similar to the independence of z(k) and G(k) that was described at the beginning of Section 2.1, z(k) is also independent of g(k) in view of (14). Therefore, from (20), it follows that Now observe that In the second line above, the mixed-product property of the Kronecker product was used and I denotes the identity matrix, while the last equality comes form (28). Similar to (33), it can be shown that E(G(k) ⊗ g(k)) = 0, and therefore, (32) leads to Finally, the last term g(k) in (19) is expressed as: Using the definition of w (n) from (10) considering (22), we obtain By defining (35) can be written as Note that the terms f w, m and f w,m ′ are independent for m ≠ m ′ because of the independent increments property of the Wiener process. 39 Thus, Therefore, (38) can be reduced to Furthermore, using (37), we obtain where Using the Itô isometry 39 one can simplify the stochastic integral (42) to the deterministic one Now, by defininĝ (40) can be reduced to That is, the second moment dynamics originally defined in (19) are given by the second moment map (31) since c z (k) (cf. (34)) where the disturbance term g is given by (46).
System (31) converges to the stationary solution as k → ∞ if and only if (

DYNAMICS AND STABILITY WITH GAUSSIAN WHITE NOISE AND DELAY-INDUCED NOISE
In this section, generalizations of the results obtained in the previous section are provided in the case of when the noise is given by where = ∑ J j=1 w j j is the mean delay and is a constant vector; cf. (22). This generalization is motivated by the fact that in some applications delay fluctuations may result is an additive noise term in the dynamics. This situation is further discussed in Section 6 where such delay-induced noise term arises in the control design of a connected automated vehicle (CAV).
When the noise (t) is given by (49), the mean dynamics are described by where Compared to (23), here the disturbance term g arises in the mean dynamics which is a result of the additive noise term in (49) due to the delay stochasticity. The details of the derivation of the disturbance term g are shown in Appendix C.
From system (50), it can be concluded that the mean z(k) converges to the equilibrium as k → ∞, if and only if (29) holds. That is, the stability condition for the mean is the same as in the previous section. The difference is that when the noise is given by (49), the equilibrium mean is nonzero and given by (53). The second moment dynamics, for the noise (t) given by (49), are described by where and The derivations of (54) to (56) are detailed in Appendix C.
Note that the difference between (54) and (31) is, that in (54), an additional term Hz(k) appears which is a result of the interaction between the Gaussian noise and the delay-induced noise. Also, the disturbance term g 2 differs from g as it has an extra term stemming from delay-induced noise; cf. (45) and (56).
Combining systems (50) and (54), one can write which has the stationary solution if and only if both (29) and (48) hold. That is, the condition for the stationary moments are the same as in the previous section. However, the stationary mean is not zero, while the stationary second moment is different, cf. (47) and (58).

AN ILLUSTRATIVE EXAMPLE
To illustrate the dynamics and stability analysis of the mean and the second moment established in the previous sections, the scalar Hayes equation with stochastic delay is considered. Assuming additive noise that is the sum of a Gaussian noise and a delay-induced noise, we obtaiṅx To investigate the stability properties and the stationary first and second moments of (59), the system (57) is constructed. The condition (29) is used to calculate mean stability, while (48) is utilized for second moment stability. The stationary mean and second moment are given in (58). The scalar versions of the terms that appear in conditions (29) and (48) and (58) are provided below.
Here we have the matrices where r j = ⌊ j ∕Δt⌋ and r J = ⌊ J ∕Δt⌋; cf. (9). Furthermore, (52) yields and from (37), f w, m can be derived as Finally, (45) becomes Given the terms in (60) to (63), to use condition (29), one can obtain G from (17), G from (20). For the stationary mean and second moments (58), one can obtain g from (51) and (52), g 2 from (45), (52) and (56), and H from (55). To generate the results below the time step Δt = 0.01 is used. This is sufficiently small to approximate well the stability boundaries of this example (see also the discussion in Reference 28 where this example is taken from). To determine the mean and second moment stability, conditions (29) and (48) are used, that is, if (G) < 1 and (G) < 1 hold, then the system is mean and second moment stable, respectively. The left panel in Figure 3 shows the stability boundaries in the (a, b) parameter space for holding times T = 0.1, 0.3, 0.5, and 1 as indicated by color. Note that the holding time values are chosen such that they span a relatively large range with respect to the delay values, that is, T = 0.1 is smaller than all the delays, T = 0.3 is the mean delay, and T = 0.5 and 1 are larger than all the delays. Dashed lines indicate mean stability boundaries while solid lines bound the second moment stability regions. That is, on the left side of the dashed lines, the mean converges to the stationary solution (58) while on the right side, it diverges to infinity. Similarly, on the left side of the solid lines, the second moment converges to the stationary solution (58) and it diverges on the right side of these boundaries.
The right panel in Figure 3 shows three different moment realizations to demonstrate the three types of stability states: moment stable, first moment stable-second moment unstable, and moment unstable. In order to illustrate the behaviour of the dynamical system (57) in the different parameter domains, the ensemble standard deviation given by is used. It can be observed that when system (59) is moment stable (case A), then both the first and second moments converge. However, as the parameters are moved toward the unstable areas first the second moment diverges (case B), then the first moment also loses stability (case C). In Figure 4A, the stationary standard deviations are shown over as a function of time for different holding times when the additive noise is purely due to the white noise excitation, that is, = 0 and = 1 in (59). The results are shown up to s = 1 for all holding times, that is, the augmented vector (7) includes the time history up to the largest holding time T = 1. The results of the semi-discretization (solid lines) are compared with the results obtained by Monte Carlo simulations (×-s). Note that the stationary standard deviation gained by both methods show periodic variation with period T due to the periodic switching of the delay process (t), even though the additive noise gives no periodic excitation. 41 In Figure 4B, we show the maximum value of the stationary standard deviation within a holding interval while setting a = −1 and varying b within [−6.5, 1]; see the dashed-dotted line in the stability chart in Figure 3. As the parameters of the system approach the stability boundary, the effect of the noise on the stationary standard deviation is magnified, leading to the so-called stochastic coherence resonance. [14][15][16] Note that the stationary mean is constantly zero for any holding time T. This is due to the fact that the white noise excitation is independent of the random delay fluctuations. This can also be verified by observing that g = 0 in (51), because f j = 0 in (61) due to = 0 as shown in Section 3.2.
In Figure 5A,B, the stationary mean and standard deviation are shown as function of time over a holding time interval for the delay-induced noise scenario ( = 1 and = 0 in (59)). Note that since the delay (t) switches at every kT, the stationary behavior shows a periodic behavior with period T. This can be observed in both the stationary mean and standard deviation. In Figure 5C,D, the maximum value of the stationary mean and standard deviation over a holding interval are depicted for the delay-induced noise case when considering parameters a = −1, b ∈ [−6.5, 1]. In Figure 5D, the stochastic resonance can be observed again for the parameters in the vicinity of the stability boundary. Meanwhile a stability loss can be observed for the stationary first moment in Figure 5C near b = 1 (upper stability boundary).
In Figure 6, the maxima of the stationary moments are shown as a function of the holding time T. In Figure 6A, it can be observed that if the holding time T is larger than the smallest delay value 1 , the stationary mean is not zero anymore. In particular, examining g in (51), one can see that for T > 1 the quantity g is nonzero yielding non-zero stationary mean, even though the excitation term ( (t) − ) has a zero expected value. The stationary standard deviation in Figure 6 does not show any special behaviour with respect to the holding time T around the smallest delay value 1 , however, after reaching a maximum it decreases as T increases. This can be due to the reason that as T takes greater values, the delay-induced noise causes resonance-like phenomenon. However, if T is further increased, the same time delay is held for longer time allowing the variations around the temporary equilibrium state to settle before the next additive time delay-induced switching in ( (t) − ) perturbs system (59).
Note that the Monte-Carlo simulations approximate the analytical results well. However, the analytical approach only requires a few matrix multiplications and solving a system of linear equations, which are orders of magnitudes faster than calculating thousands of realizations and statistically evaluating them. This suggests that the proposed method is a very efficient tool to investigate the behavior of such systems, especially for higher dimensional state variables.

STOCHASTIC EFFECTS IN CONNECTED VEHICLE SYSTEMS
In this section, we consider the influence of stochastic time delay and addtive noise on connected vehicle systems. In particular, we focus on the longitudinal dynamics of a connected automated vehicle (CAV) following a connected human-driven vehicle (CHV) that broadcasts its GPS position and speed via wireless vehicle-to-vehicle (V2V) communication. When receiving the packets, the CAV can respond to the motion of the CHV by adjusting the throttle or applying the brakes. This is referred as connected cruise control and the effects of time delays in such systems have been investigated both theoretically and experimentally. 42,43 There are two sources of delay-in the this example. On the one hand, the actuator delay of the CAV is typically constant and in the range of 0.3 to 0.6 seconds. On the other hand, the communication delay is in the range of 0.1-0.3 second and this changes stochastically based on the random nature of packet scheduling algorithms and the packet drops in wireless communication.
The stochastic process describing the motion of the CAV is characterized by its second moment dynamics. First, the solution of the stochastic system describing the dynamics of the CAV is partitioned into a deterministic part and into a stochastic part. Next, the second moment stability and the steady-state second moment of the stochastic part are investigated with the help of the method introduced above. Finally, it is shown how the CAV can benefit from applying the so-called delay matching to suppress the fluctuations induced by the stochastic delay while maintaining its robustness against the external noise presented by the CHV ahead.
The stochastically delayed differential equation describing the motion of the CAV can be written aṡ Here the dot stands for differentiation with respect to time t, s and s 1 denote the positions of the rear bumpers of the CAV and the CHV ahead, while v and v 1 denote their velocities, respectively; see Figure 7A. The length of the CAV is denoted with l and the headway is defined by In (65), the gains and are used to correct velocity errors, represents the actuator delay of the CAV, while 1 incorporates the communication delay as well as the actuator delay. The desired velocity is determined by the nonlinear (piecewise linear) range policy function shown in Figure 7B, where = v max ∕(h go − h st ). That is, the desired velocity is zero for small headways (h ≤ h st ) and equal to the speed limit v max for large headways (h ≥ h go ). Between these, the desired velocity increases with the headway linearly, with gradient . Note that when h st = 0, the qunatity 1∕ is often referred to as the time headway. We consider that the CHV's velocity satisfy E(v 1 (t)) = v * which allows us to partition its position and velocity as ] . (68) To describe the sum of the measurement error and the perturbation of the motion of the lead vehicle, the vector In the case of the CAV, the same partitioning leads to where the vector x ∶= [x s , x v ] ⊤ collects the position and velocity perturbations of the CAV. Note that for connected cars, the holding time is always smaller than the minimum time delay of the system, that is, T < 1 , leading to E(x(t)) = 0. Substituting the definitions of the positions and the velocities of the vehicles from (68) and (69) into (65), using the range policy (67), and taking the expected value of the resulting equation, the stationary average headway distance can be determined: Here and 1 denote the average values of the delays and V −1 is only unique for 0 < v * < v max ; cf. (67). This shows that the average mismatch between the delays may result in an (undesired) increase of the stationary headway.
To characterize the quality of the CAV control in case of a dense, but continuously flowing traffic scenario (h stop < h(t) < h go , ∀t), the dynamics in the middle linear section of the range policy V(h) need to be studied. This is described by the linear systeṁ where Note that in (71), the terms containing x 1 (t) and v 0 act as excitations on the system. Furthermore, it is assumed that the perturbation x is small and the headway h does not leave the interval [h st , h go ].

Effect of delay matching on connected cruise control
It was shown by (70) that a mismatch between the delays and 1 result in a shift from the desired stationary solution.
Since the packets sent via V2V communication are time stamped, one may add some delay to the actuation delay to so that it matches 1 . Here we utilize the analytical techniques established above to evaluate the performance of the CAV when we apply vs. do not apply such delay matching.
In case of delay matching, the controller sets and, since the vector v 0 is in the nullspace of the matrix B + B 1 , the linear system (71) simplifies tȯ In case of no delay matching, the actuation delay is left constant, that is, and (71) becomesẋ with no stochasticity in the delay . For simplicity, the case when the leading vehicle is moving with a constant speed is considered, that is, x 1 (t) ≡ 0. In this case, (74) has no added noise, while (76) is only excited by the delay noise. When presenting the results, we use the parameter = 0. Note that with these parameters the delay matching not only causes the delay (t) to be stochastic, but it also increases its mean .
The stability charts are plotted in the plane of the control gains and in Figure 8A by utilizing (48). Note that the stability domain is smaller for the case with delay matching due to the fact that is increased in this case (unlike in case of, for example, suppressing chatter in machine tool vibrations via spindle speed variation, where the delay is varied around the original constant delay 30 ). However, since there is no added noise in this case, we obtain lim t→∞ x 2 s = 0 and lim t→∞ x 2 v = 0 for the stationary second moments. This is illustrated by the numerical simulation in Figure 8B corresponding to point P located at ( , ) = (0.4, 0.5) in the stability chart. The simulated trajectories (grey curve) approach zero quickly bringing the mean (thick green curve) and the standard deviation (thin green curve) to zero too.
When no delay matching is used, the stable domain is larger but, due to the delay-induced noise added to the system, the stationary second moment is not zero. The contours of lim t→∞ x 2 v are calculated using (58) are shown in Figure 8A. The numerical simulations in Figure 8C illustrate the dynamics for point P marked in the stability chart. The sample trajectory (grey curve) shows that the velocity of the CAV keeps changing in time so that the mean (thick blue curve) approaches zero while the standard deviation approaches the constant value calculated analytically.
Notice that while delay matching reduces the size of the stability domain, one may still have a large range of gain parameters to choose from. In particular, in the experimentally realistic range of ∈ [0, 1], ∈ [0, 1] stability is still ensured. Furthermore, the delay matching eliminates the unwanted stationary oscillations that lead to "jerky ride," which typically has a negative effect on driver comfort and energy consumption. Note that in case of no delay matching, the delay-induced noise may be magnified to an extent, where (71) is not valid anymore (the range policy V(h) saturates). However, in terms of control design, these parameter regions should be avoided, and instead of the stability boundary, the stationary second moment contours should be considered when choosing control parameters for the CAV to limit the amplitude of the perturbation x.
The positive effects of delay matching can also be observed when the CHV ahead varies its velocity slowly (that is typical in real traffic situations). 44 In order to evaluate the effects for more severe motion perturbations, we model the perturbation dynamics of the CHV using white noise, that is, where s (t) and v (t) are uncorrelated Gaussian white noise. This is indeed an over estimation of the severity of perturbations but provides a way to compare the behavior with and without delay matching. This leads to a noise excitation (49) with in case of delay matching The results are summarized in Figure 9 where panel (A) depicts the stability charts. The stability boundaries are identical to those in Figure 8A as these are still obtained by (48) while the contours calculated by (58) change due to the added white noise. Note that the contours obtained for the delayed matching case are quite similar to those obtained without delayed matching for small gain values, including the experimentally realistic range ∈ [0, 1], ∈ [0, 1]. For point P at ( , ) = (0.4, 0.5), this behavior is also illustrated by the numerical simulations shown in panels (B) and (C).
These demonstrate that when responding to severe perturbations using delay matching neither improves nor degrades the performance of the CAV.
To summarizes this section, we demonstrated how the developed mathematical tools can be utilized for the control design of a connected vehicle system where stochastic packet loss results in stochasticity in the delay. Note that no external noise was added to the control system (65). However, some delay-induced additive noise was inherently present due to the packet drops (see Equations (68) to (72)), which had the same stochastic switching behavior as the stochastic delay. It was shown that this delay-induced additive term can be eliminated by delay-matching, that is, by deliberately delaying actuation of the CAV when a packet drop occurs. In (77) and (78), an additional white noise excitation, coming from the stochastic perturbation of the lead vehicle, was introduced, and it was shown that with delay-matching the controller maintains its robustness against external noise excitation.

CONCLUSIONS
In this article, dynamical systems with stochastic delays and additive noise were investigated. The delay was assumed to jump between finitely many different values while remaining constant for a given holding time between the jumps. The additive noise was constructed from two different stochastic processes: one was an independent Gaussian white noise, called Wiener process, while the other was a process generated by the stochastically changing delay. A discrete time stochastic map was derived using semi-discretization and this was used to determine the dynamics of the first and second moments. Stability of the stationary moments was ensured by making the spectral radii of the corresponding coefficient matrices smaller than 1 and the stationary first and second moments were calculated. The developed methods were applied to the stochastically delayed Hayes equation and it was demonstrated that the stationary first and second moments are periodic with period equal to the delay holding time. It was also shown that when the holding time was larger than the smallest possible time delay, the stationary first and second moments were not zero even when the noise exciting the system had zero mean. These results were also confirmed by statistical evaluation of Monte-Carlo simulations.
Finally, the established mathematical tools were used to support the claim, that matching the actuation delay (thus making it stochastic) in a connected vehicle system can improve the performance of CAVs. This occurs because delay matching eliminates the additive noise caused by the stochastically arriving information packets. It was shown that carefully choosing the control parameters one can preserve the robustness of the delay matched system even in the presence of a harsh white noise excitation.

APPENDIX A. HIGHER ORDER SEMI-DISCRETIZATION METHOD
Using the methods in Reference 38, system (1) can be discretized to: where t ∈ [nΔt, (n + 1)Δt), n ∈ N and the time resolution is Δt = T∕l is given by the holding time T and l ∈ N. Here the delayed state is approximated by a qth order Lagrange polynomial, where and In the main body of the article, the zeroth order semi-discretization is used and indeed (A1), (A2), and (A3) give (1) and (4) when q = 0.

C.1 Derivation of mean dynamics
Here we want to show that when the additive noise is given by (49), the mean dynamics (16) reduce to z(k + 1) = Gz(k) + g, and we want to derive the disturbance term g; see (50) to (52). Let us first repeat the disturbance term from (16) Consider the last term above. In view of (9) and (10) assuming that A is invertible. Substituting (C5) in (C3), we get where f j andF j are given in (52).

C.2 Derivation of second moment dynamics
Recall the second moment dynamics given by (19). Here (19) is simplified to the form z(k + 1) = Gz(k) + Hz(k) + g 2 , and obtain the terms H and g 2 when the noise is given by (49); see (53), (55), and (56). Recall that z(k) is independent of g(k) and G(k) as mentioned before. Therefore, from (20) we have c z (k) = E(g(k) ⊗ G(k) + G(k) ⊗ g(k))E(z(k)). where Let us now consider the last term in (19), that is, g(k). We have that g(k) = E(g(k) ⊗ g(k)) . (C12) Note that, given (t) = j , t ∈ I k , we can write f(kl + m) as in view of (10) .