Stability of Linear Continuous-Time Systems With Stochastically Switching Delays

Necessary and sufficient conditions for the stability of linear continuous-time systems with stochastically switching delays are presented in this paper. It is assumed that the delay random paths are piece-wise constant functions of time where a finite number of values may be taken by the delay. The stability is assessed in terms of the second moment of the state vector of the system. The solution operators of individual linear systems with constant delays, chosen from the set of all possible delay values, are extended to form new augmented operators. Then for proper formulation of the second moment in continuous time, tensor products of the augmented solution operators are used. Finally the finite-dimensional versions of the stability conditions, that can be obtained using various time discretization techniques, are presented. Some examples are provided that demonstrate how the stability conditions can be used to assess the stability of linear systems with stochastic delay.

in the former category where Lyapunov-based theorems are applied, for instance [17]- [20]. The drawback of Lyapunovbased techniques is that they provide sufficient conditions for stability that are restrictive in most cases. On the other hand, necessary and sufficient stability conditions, that are the goal of the exact methods, are harder to attain. In the case of discrete-time linear systems, necessary and sufficient conditions for stability have been derived, for example in [15], [16], [21]. However, for continuous-time systems with stochastic delays, necessary and sufficient stability criteria are lacking in general. For instance in [22], exact stability conditions are derived but only for a very particular type of delay behavior.
In this paper, we propose necessary and sufficient conditions for the stability of a class of linear continuous-time systems where the delay is subject to stochastic variations. The main stability result is concerned with the stability of the second moment of the system. The stability criteria are based on the spectral radius of operators that are constructed using solution operators associated with individual delays. A preliminary version of this work has appeared in [23]. Here we remove a restrictive assumption on the delay behavior that was used in [23]. In particular, it was assumed in [23] that the delay switching is fast relative to delay values. This was done by assuming that the delay dwell time at a value was less than the minimum delay in the system. However, in the present work the delay switching can be fast or slow, i.e. the delay dwell time is arbitrary. We develop a different technique in constructing some useful operators in Section III that allow us to relax the aforementioned restrictive assumption used in [23]. Also we provide more examples and discussion in this work compared to the preliminary version.

II. PROBLEM STATEMENT
Consider the linear systeṁ where x ∈ R n , a, b ∈ R n×n , and the delay τ (t) ∈ R changes stochastically with time. We assume that the delay can take values from a finite set Ω = {τ 1 , τ 2 , . . . , τ J } where 0 < τ 1 < τ 2 < · · · < τ J = τ max . The initial condition is given by where φ ∈ C [−τ max , 0], R n and C denotes the space of continuous functions. Trajectories of the scalar version of system (1) with a = −1 and b = −2 corresponding to sample paths of the delay shown in Fig. 1a and initial condition φ(θ) ≡ 0.1, −1.5 ≤ θ ≤ 0.
Before switching to a different value, the delay stays at the current value for a duration of time t d , which we call dwell time. Therefore the value of the delay at each interval [kt d , (k + 1)t d ), k = 0, 1, 2, . . ., is constant. The probability distribution w = [w 1 w 2 . . . w J ] governs the switchings of the delay where w j is the probability of switching to the delay τ j . The probability distribution w is assumed to be stationary which means the switchings of the delay are independent, identically distributed (i.i.d.). Note that it is assumed that the past information {x(t) : t ∈ [t − τ max , t]} is available at time t, for all t ≥ 0. Fig. 1a shows two sample paths of the delay τ (t) and Fig. 1b shows the corresponding trajectories x(t) of the scalar version of system (1), i.e. for the case n = 1, with parameters and initial conditions given in the caption of Fig. 1.
Our goal is to study the stability of the stochastic system (1). In particular, we study the stability of the mean E[x(t)], and the second moment E[x(t)x T (t)], where E[.] denotes the expected value of a random variable. To this aim, we need a proper representation of system (1) given the delay behavior as described. In the next section, we construct this representation using solution operator formulation of delay differential equations and in Section IV we provide a suitable definition of the second moment using tensor products of appropriate operators.

III. SOLUTION OPERATOR REPRESENTATION OF THE SYSTEM
First we recall the definition of the solution operator for a deterministic delay differential equation. Consider the linear systemẋ where φ ∈ C([−τ, 0], R n ). The solution operator for system (3) is defined by The operator T (t), t ≥ 0, is bounded and linear and the family of operators T (t) is a strongly continuous semigroup that has the properties see [1] or [2] for more details. Now consider the deterministic systemsẋ with their respective solution operators T j (t) : As described in Section II, the delay is constant in each interval [kt d , (k + 1)t d ), k = 0, 1, 2, . . .. Therefore, the stochastic system (1), in the time interval [kt d , (k + 1)t d ), evolves according to one of the systems in (6). In other words, if τ (t) = τ j in the time interval [kt d , (k + 1)t d ), the operator T j (t d ) progresses the solution forward from t = kt d to t = (k + 1)t d . In order to obtain an appropriate representation of the stochastic system (1) for stability analysis, we extend the operators T j , j = 1, . . . , J, in a way that the new, extended operators share a common domain.
Let Γ := max{τ max , t d } and denote by C the space C([−Γ, 0], R n ) and by C j the space C([−τ j , 0], R n ), j = 1, . . . , J. Now, first we define the auxiliary operators U j : C → C j by for j = 1, . . . , J. In other words, the operator U j acts on a continuous function φ from C and outputs the segment of φ corresponding to −τ j ≤ θ ≤ 0. Now choose h > 0 such that t d = h where ∈ N is a positive integer and h < τ j , j = 1, . . . , J. Next define the operatorsT j : C → C by for j = 1, . . . , J. The operatorsT j are in fact extensions of the operators T j (h) from the domain C j to domain C. This extension is performed by adding the bottom line in (8).
Finally we define the operators for j = 1, . . . , J. The operator G j : C → C is consecutive applications of the operatorT j . Note that the operators G j , that are constructed through (7)-(9), are linear and bounded because the original operators T j are bounded and linear. Now assume that the stochastic system (1) is realized up to the time t = kt d . Let us define as the "state" of the system at time t, where x t ∈ C. Then, assuming that τ (t) = τ j in the time interval [kt d , (k + 1)t d ), we can write Consequently one can construct the stochastic system where P G(k) = G j = w j , j = 1, . . . , J, and k = 0, 1, 2, . . . , with J j=1 w j = 1, w j ≥ 0, and initial condition Here P denotes probability. Note that we cannot arrive at a stochastic system of the form (12) if we want to use the original solution operators T j , due to the fact that they do not have a common domain. In the next section, we study the stability of system (12).

IV. STABILITY ANALYSIS
In this section, we derive stability conditions for the mean and the second moment of the stochastic system (1) by studying system (12). First we provide a standard definition as well as a standard result that we will use in this section.
Definition 1: System where A : X → X is a bounded, linear operator on the Banach space X , is exponentially stable if for every initial condition x 0 ∈ X , there exist M ≥ 1 and 0 ≤ r < 1 such that see for instance [24]. A standard result on the stability of system (15) is provided in the following lemma; see for instance Theorem 2.1 in [24].
Lemma 1: Consider system (15). Let σ(A) denote the spectrum of A and denote the spectral radius of A. Then, system (15) is exponentially stable if and only if ρ(A) < 1.
(18) Note that Definition 1 and Lemma 1 concern discrete stability as system (12) is a discrete model. Now we consider the mean of the stochastic system (12). By taking the expected value of (12), we have Note that the operator G(k) only depends on the delay value in the time interval kt d , Since the delay value in the time interval kt d , (k + 1)t d is independent of the delay values in other time intervals (due to the i.i.d. assumption), the operator G(k) is independent of x kt d . Thus, from (19), we arrive at where Now applying Lemma 1 on system (20) results in a necessary and sufficient condition for the stability of the mean of the stochastic system (12). This is provided in the proposition below.
The proof is immediately obtained by the application of Lemma 1 on system (20). Note that the operator (21), is a finite summation of bounded and linear operators G j and so is bounded and linear, and moreover it is defined on the Banach space C.
Our main goal is to derive necessary and sufficient conditions for the stability of the second moment of system (12). For a proper description of the second moment dynamics, we use the tensor product of the Banach space C with itself equipped with an appropriate cross norm. The connection between the second moment and the tensor product space lies in the definition of the norm.
Let X denote a Banach space and X ⊗ X denote the tensor product of X with itself. A standard norm on tensor product spaces is the injective norm that is given by where u = M m=1 x m ⊗y m is a tensor in X ⊗X and B X * is the closed unit ball on X * (the normed dual of X ). In other words, f and g are bounded, linear functionals defined on X with norm less than or equal to 1, and the supremum in (27) is taken over all such f and g. Furthermore, in the definition (27), one can substitute B X * with a norming set. A subset N of B X * is said to be a norming set if x = sup{|f (x)| : f ∈ N } for every x ∈ X ; here . is the norm defined on the Banach space X . See [25], chapter 3, for more details about the injective norm on tensor product spaces.
Now we first define a norm on the tensor space C ⊗ C and then show that it is in fact the injective norm on C ⊗ C.
Definition 2: Let u = M m=1 φ m ⊗ ψ m belong to C ⊗ C. We define the c-norm on C ⊗ C to be In particular, for a simple tensor of the form The connection between the second moment and the c-norm defined in (28) can be seen from (29) where φ⊗φ c contains the products of the values of the function φ at different arguments. The second moment in the infinite-dimensional setting of functions may also be understood as the expected value of the product of a function with itself at different argument values.
Next we show that the c-norm defined in (28) is the injective norm on C ⊗ C.
Proof: Let φ ∈ C. Consider the linear functionals where Therefore, we can substitute B C * with N in (27). As a result, (27) yields (28).
We denote by C ⊗ c C the tensor product space equipped with the c-norm (injective norm) and by C⊗ c C the completion of this space under the c-norm. It is important to note that many other norms can be defined on the tensor product space C ⊗ C, such as the projective norm [25], however, the norm . c given by (28) may not be equivalent to those other norms. The usefulness of the c-norm (or equivalently injective norm) for the second moment analysis will be further illuminated in the remainder of this section. Consider the bounded and linear operator G j : C → C. For each j = 1, . . . , J, there exists (see Proposition 3.2 in [25]) a unique, bounded, and linear operator G j ⊗ c G j : C⊗ c C → C⊗ c C such that Now assume that system (12) is realized up to the time t = kt d and in the interval kt d , (k + 1)t d the delay is τ (t) = τ j . Using the operator in (32), one can write Substituting Eq. (11) into (33) yields Therefore, one can construct the stochastic map where P(G(k) ⊗ c G(k) = G j ⊗ c G j ) = w j , j = 1, . . . , J, and k = 0, 1, 2, . . . , with the initial condition x 0 ⊗ x 0 = φ ⊗ φ. Now we can take the expected value of Eq. (35) that results in (37) Due to the independence of G(k) and Note that E G(k) ⊗ c G(k) is a bounded, linear operator on the Banach space C⊗ c C. Now we can state a theorem that provides a necessary and sufficient condition for the stability of the second moment of (12). Theorem 1: Consider system (12) which is repeated below where P G(k) = G j = w j , j = 1, . . . , J, and k = 0, 1, 2, . . . , with initial condition There exist M ≥ 1 and 0 ≤ r < 1 such that where G j 's are given by (9). Proof: By application of Lemma 1 to system (38), we can say that there exist M ≥ 1 and 0 ≤ r < 1 such that On the other hand, from the definition of the c-norm in (28), we have that completes the proof.
Theorem 1 provides a necessary and sufficient condition, i.e. condition (44), for the second moment stability of system (12), i.e. Eq. (43). Note that using the tensor product provided us with appropriate means to construct a linear map such as (35) for the second moment stability analysis. While E[x kt d ⊗x kt d ] may be interpreted as the "second moment," it is not exactly the case. However, the norm defined on the tensor product space C⊗ c C enables us to provide a supremum norm on the second moment of system (12); cf. (47). Now we recall that the original question in Section II was concerned with the stability of E[x(t)x T (t)], i.e. the second moment of the stochastic system (1). In the following corollary, we show that condition (44) of Theorem 1 is truly a necessary and sufficient condition for the second moment stability of system (1).
Corollary 1: Consider system (1) with the delay behavior as described in Section II. There exists M ≥ 1 and ω > 0 such that if and only if where G j 's are given by (9).
Then from Theorem 1, there exists M ≥ 1 and 0 ≤ r < 1 such that (43) holds for any k ∈ {0, 1, . . .}. Now for any t ≥ 0, there existsk ∈ {1, 2, . . .}, such that (k − 1)t d ≤ t <kt d . Thus, by choosing θ 1 = θ 2 = −kt d + t in (43) and recalling that xk t d (θ) = x(kt d + θ), we have Also, since t/t d <k, then rk < r t/t d = e whereω = − 1 t d log r. To show the reverse, assume that there exists M ≥ 1 and ω > 0 such that (48) holds. Choose M 1 > 0 such that . Note that since φ ∈ C (φ is the initial condition of system (1)), we know that such M 1 exists. Let M 2 = max{M, M 1 }. Hence, In fact (52) is the extension of (48) to the interval t ≥ −Γ. Now consider any k ∈ {0, 1, 2, . . .} and the time interval [kt d − Γ, kt d ]. Setting i 1 = i 2 = i in (52) and using the notation x kt d (θ) = x(kt d + θ), we get Therefore, by defining θ = −kt d + t, (53) can be written as . . , n}. Substituting (54) into the right hand side of (55), we get ∀ k ∈ {0, 1, 2, . . .}. According to the result of Theorem 1, (57) implies Condition (24) in Proposition 1 and condition (44) in Theorem 1 are, respectively, necessary and sufficient conditions for the stability of the mean and the second moment of the stochastic system (1). However, in practice these conditions cannot be investigated directly due to the infinite-dimensional nature of the relevant operators. In the next section, we provide finite-dimensional versions of the stability conditions obtained in this section.

V. FINITE-DIMENSIONAL APPROXIMATIONS
Note that the stochastic system (1) evolves in continuous time, even though the delay takes discrete values (see Fig. 1). By discretizing system (1) in time, one can obtain a finite-dimensional approximation of system (1) and finitedimensional approximations of the operators G j given by (9). There are many well-established time discretization techniques for delay differential equations, such as Runge-Kutta techniques [26] and semi-discretization technique [27], that one can use. In this section, we provide finite-dimensional approximations of the stability conditions obtained in Section IV independent of the specific discretization method used.
Assume that one applies a discretization method with a constant step-size mesh and obtains as a finite-dimensional approximation of system (12) where X(k) ∈ R N and G(k) ∈ R N ×N , k = 0, 1, 2, . . .. The integer N ∈ N is a parameter of the discretization method that is a function of the step size ∆t, the largest delay τ max or the dwell time t d (whichever is greater), and the dimension n of the vector x in system (1). Similar to system (12), here P(G(k) = G j ) = w j , j = 1, . . . , J. Moreover, the matrix G(k) and vector X(k) are independent (due to the i.i.d. assumption on the delay switchings). Therefore, by taking the expected value of (59), we arrive at where Based on (60-61), we know that E[X(k)] is exponentially stable if and only if Condition (62) is a finite-dimensional approximation of the mean stability condition (24) given by Proposition 1.
To derive a finite-dimensional approximation of the second moment stability condition, observe that the tensor product becomes the Kronecker product in finite dimensional spaces. Therefore, from (59) one can write where ⊗ denotes the Kronecker product, . System (63) is a finitedimensional approximation of system (35) where P G(k) ⊗ G(k) = G j ⊗G j = w j , j = 1, . . . , J. By taking the expected value of Eq. (63) and using the independence of G(k) and X(k), we have Based on (64-65), we know that E[X(k) ⊗ X(k)] is exponentially stable if and only if Condition (66) is a finite-dimensional approximation of the second moment stability condition (44) given by Theorem 1. After discretization of system (1) and obtaining the matrices G j , one can use condition (62) and (66) to evaluate the stability of the mean and the second moment of system (1), respectively. Conditions (62) and (66) can, for instance, be used to obtain approximate stability boundaries in desired parameter spaces. The stability boundaries obtained using conditions (62) and (66) converge by decreasing the time step of the underlying discretization technique. The convergence properties, such as the order of convergence, of ρ J j=1 w j G j to ρ J j=1 w j G j and ρ J j=1 w j G j ⊗ G j to ρ J j=1 w j G j ⊗ c G j follow the convergence properties of the discretization technique used. For more details about some suitable discretization methods, we refer the reader to [28] (for Runge-Kutta techniques) and [27] (for semidiscretization technique). We shall also remark that a full characterization of the convergence properties is the subject of an ongoing work by the authors. In the next section, we demonstrate the application of the stability conditions using some examples.

VI. EXAMPLES
In this section, we provide two examples demonstrating the application of the stability criteria obtained in Sections IV and V.

A. Scalar System
Consider the scalar version of system (1), that is, where x, a, b ∈ R. Assume that the delay can take two values from the set Ω = {0.5, 1} with probability distribution w 1 = w 2 = 0.5, and the dwell time of the delay is t d = 0.25. We use a zeroth-order semi-discretization technique [27] with a step size of ∆t = 0.025 to discretize system (67) and obtain matrices G 1 and G 2 associated with the two delay values (see [27], chapter 3, for details of the semi-discretization technique and [29] for how to obtain the matrices). We want to draw stability charts in the space of parameters a and b. Note that matrices G 1 and G 2 are functions of the parameters a and b. For the stability of the mean of system (67) we check the condition given by Proposition 1 that reduces here to (cf. (62)) and for the stability of the second moment of system (67), we check the condition given by Theorem 1 that reduces here to ρ 1 2 (cf. (66)).
In Fig. 2, the stable and unstable areas of the mean (based on condition (68)) and the second moment (based on condition (69)) are shown with blue and red curves, respectively. Note that the area on the left of the boundaries (that is limited by the line b = −a from top) is the stable region. The second moment stable region is inside the mean stable region as the second moment stability is sufficient for the mean stability. For comparison, the stable regions of the deterministic version of system (67) with deterministic delays τ = 0.5 and τ = 1 are shown by solid and dashed gray curves, respectively.

B. A Vector (2-D) System
In this section, we consider a linear second order system in the general canonical reachable forṁ where x ∈ R 2 . The control action u ∈ R is applied with a delay τ (t). Using the feedback control law the closed-loop system becomeṡ where Let us assume that the delay takes values from the set Ω = {0.7, 1, 1.3} with equal probabilities w i = 1/3, i = 1, . . . , 3.
For a fixed set of plant parameters (a 1 , a 2 ) = (1, 50), and different delay dwell times t d = 0.5, 1.5, and 10, we seek for the values of control gains (k 1 , k 2 ) for which the closed loop system (72) is second moment stable. To investigate the second moment stability of system (72-73) with the given delay parameters, we use the finite-dimensional version of the second moment stability condition (i.e. (66)) that reduces here to Here we use a first-order semi-discretization technique to obtain the matrices G j associated with 3 delay values. See [27], Chapter 3, for how to construct these matrices using the first-order semi-discretization.
The solid, dashed, and dashed-dotted gray curves in Fig. 3 encircle the stable areas for the deterministic version of system (72) with deterministic delays τ = 0.7, τ = 1, and τ = 1.3, respectively. Note that these boundaries can be obtained analytically exploiting the corresponding characteristic equations. The red, red-purple, and purple curves encircle the stable areas for the stochastic system (72) with dwell times t d = 0.5, t d = 1.5, and t d = 10, respectively. Note that these boundaries can be obtained by checking condition (74) point by point in the (k 1 , k 2 ) space, or alternatively, by finding an initial point on the stability boundary and obtaining the rest of the boundary using a continuation technique. Here we have used the latter employing the continuation routine embedded in the software package DDE-BIFTOOL [30], [31]. Note that the stochastic system has the biggest stable area for the shortest dwell time (fastest delay switching). As the dwell time gets larger, the stable area of the stochastic system shrinks to the area where all 3 delays are deterministically (individually) stable. This observation in this example matches intuition in the following sense. As the dwell time goes to infinity the stochastic system resides in each individual deterministic system for a long time, and thereby one expects all individual deterministic systems to be stable for the stochastic system to be stable. Another observation is that the control gains that stabilize the stochastic system do not necessarily stabilize any of the individual deterministic delay systems. Note that as the time step ∆t used in the time discretization techniques in the above examples gets smaller, the size of the matrices G j used in condition (66) gets larger. We made ∆t small enough so that the boundaries shown in Figs. 2 and 3 converged to a desired accuracy. For higher-dimensional systems (larger n) the computational cost associated with (66) may get very large as the size of the matrices G j ⊗ G j is proportional to n 2 . In these cases, using higher-order time discretization techniques may help reduce the overall cost.

VII. CONCLUDING REMARKS
We obtained necessary and sufficient conditions for a class of linear systems subject to stochastic delay. We considered the stability of the second moment of the system as stability criteria. The stability conditions were given in the form of the spectral radius of an operator that is a linear combination of the tensor products of augmented solution operators associated with each individual delay. We presented finite-dimensional approximations of the proposed stability criteria which can be used to assess stability numerically. For instance one can draw stability charts in the parameter space of interest. While the class of systems we considered had one delayed term with stochastic delay, the presented method can be generalized to the case where there are some terms with deterministic delays and some terms with stochastic delays in a straightforward fashion.
For the linear systems considered, the stability of the second moment also provides a sufficient condition for almost sure stability. Therefore, if an assessment of almost sure stability is desired, the result of this paper can be useful.
The assumed delay behavior is flexible to approximate different kinds of stochastic behavior in the sense that there are three parameters to tune: delay values, the probability distribution, and the dwell time. Furthermore, the generalization of this delay behavior to the case where the dwell time is a random variable as well as the case where the jump probabilities follow a Markov chain rule (which is more general than i.i.d.) can be done using the same machinery presented in this paper. Another more general extension of this work, which is a subject of future study, would be the case where the delay is changing continuously (and stochastically) with time.