LYAPUNOV FUNCTION COMPUTATION FOR AUTONOMOUS LINEAR STOCHASTIC DIFFERENTIAL EQUATIONS USING SUM-OF-SQUARES PROGRAMMING

. We study the global asymptotic stability in probability of the zero solution of linear stochastic diﬀerential equations with constant coeﬃcients. We develop a sum-of-squares program that veriﬁes whether a parameterized candidate Lyapunov function is in fact a global Lyapunov function for such a system. Our class of candidate Lyapunov functions are naturally adapted to the problem. We consider functions of the form V ( x ) = (cid:107) x (cid:107) pQ x where the parameters are the positive deﬁnite matrix Q and the number p > 0. We give several examples of our proposed method and show how it improves previous results.


1.
Introduction. In deterministic Dynamical Systems given by autonomous ordinary differential equations (ODE), the (global) stability of an equilibrium is an important property, which can be studied using Lyapunov functions. For linear ODEs, a (global) Lyapunov function for exponentially stable equilibria can be constructed as a quadratic form by solving a linear matrix equation. Moreover, this method also enables us to construct a Lyapunov function for a nonlinear equation: the Lyapunov function for the linearized system at the equilibrium is also a Lyapunov function for the nonlinear system, but usually only in a small neighborhood of the exponentially stable equilibrium. Apart from this method, there are no general analytic construction methods for Lyapunov functions, however, there are numerous numerical ones [8].
Stability in stochastic differential equations can also be analyzed using Lyapunov functions [13,12,24,7] and by considering viscosity solutions to the stochastic Zubov equation the authors of [9] even gave a numerical scheme to approximate Lyapunov functions that deliver attraction probabilities to exponentially stable sets. However, even for autonomous linear stochastic differential equations, there is no general construction method for Lyapunov functions.
In this paper, we will propose an ansatz for a Lyapunov function of the form V (x) = x p Q := (x Qx) p 2 , where Q is a positive definite matrix Q and p > 0, and develop a method to check whether this is a (global) Lyapunov function and thus to prove global asymptotic stability in probability of the zero solution. The method requires checking whether a certain polynomial is positive semi-definite, which is achieved by SOS (sum-of-squares) programming. Our method allows for optimizing with respect to certain parameters. This provides us with a tool which improves the current ability to analyze stability, as we show in examples. Future work will show that, similar to the ODE case, this method also enables us to construct Lyapunov functions for nonlinear stochastic differential equations in a neighborhood of the zero solution.
The SOS method has been used to study the stability of different stochastic systems before. In [5], for example, linear and nonlinear polynomial systems with probabilistic uncertainties in the system parameters are considered, and they are transformed into a higher dimensional deterministic system, where classical Lyapunov theory is applied; in particular, polynomial Lyapunov functions are constructed using SOS. The authors use SOS and BMI (bilinear matrix inequality) in [6] to design a controller for such systems. In [23], SOS is used to compute a bound on the probability of reaching a target set in finite time for jump-diffusion processes. [3] considers linear stochastic impulsive systems and studies mean-square stability and dwell-times by infinite-dimensional LMIs, namely matrix-valued functions. A sufficient condition is obtained by assuming that these functions are polynomials of a certain degree, and these conditions are then checked by using SOS. To the best of our knowledge, however, using SOS to compute a non-polynomial Lyapunov function of the form x p Q , in particular with 0 < p < 1, for linear autonomous SDE's has not been proposed before.
The outline of the paper is as follows: In Section 2 we recall various definitions of stochastic stability, and in Section 3 we discuss Lyapunov functions. Section 4 contains our main results, in particular we propose our ansatz for a Lyapunov function and describe how the nonnegativity is analyzed using sum-of-squares. Finally, we apply the method to examples in Section 5 and show the improvements upon previous results.
Remark on notations. We denote the nonnegative real numbers by R + , the integers by Z, and the nonnegative integers by N 0 . We define [a, b] Z := {x ∈ Z : a ≤ x ≤ b}. We write vectors x ∈ R d in boldface and their components with the same letter indexed. Further, we assume all vectors to be column vectors, i.e. x = (x 1 , x 2 , . . . , x d ) , where denotes the transpose. We denote matrices A ∈ R n×m by capital letters and their elements by the same capital letter indexed, i.e. A = (A ij ) i∈ [1,n] Z ,j∈ [1,m] Z . I d denotes the identity matrix in R d×d . Q 0 means that the symmetric matrix Q ∈ R d×d is positive semi-definite and Q 0 means that Q is positive definite. Similarly Q 0 and Q ≺ 0 means that Q is negative semi-definite and negative definite, respectively.
We denote the Euclidean norm of a vector x ∈ R d by x . For Q ∈ R d×d , Q 0, we define x Q := (x Qx) 1 2 , note that · Q is a norm on R d . We denote the spectral norm of a matrix A ∈ R d×d by A 2 , i.e. A 2 is the square-root of the largest eigenvalue of the positive semi-definite matrix A A.
Thus x α is a homogenous monomial of degree |α|.
In general we call a function f : for all ρ > 0 we call f positively homogenous of degree m. We denote the set of all continuous functions α : R + → R + that are strictly increasing and such that α(0) = 0 and lim x→+∞ α(x) = +∞ by K ∞ . P is a probability measure and E denotes the expectation, the underlying probability-spaces being the canonical ones for the SDE at hand, cf. e.g. [13,18,12].
For an open set U ⊂ R d we denote the set of all functions U → R that are continuous by C(U) and the set of all functions that are m-times continuously differentiable by C m (U). For a sufficiently smooth function f : R d → R, we denote the first and second partial derivatives by We are working in R d if not specified otherwise and if the limits in a sum are suppressed, it means that the sum is over [1, d] 2. Stochastic stability. We consider the d-dimensional, autonomous linear stochastic differential equation (SDE) with constant coefficients where F, G u ∈ R d×d and the W u are independent 1-dimensional Wiener processes.
We understand solutions to the system (1) in the Itô sense, cf. e.g. [13, §2.2], especially Definition 2.1 and the discussion before Theorem 3.6, i.e. that (1) is a shorthand notation for where the integrals over the G u X x (s) are understood as Itô integrals. We call t → X x (t) in (2) a solution to system (1) with initial value x ∈ R d . In general x can be a random variable [ which in turn is a special case of the nonautonomous SDE For a study of existence and uniqueness of solutions to the system (4) we refer the reader to [13,Chapter 2], [18,Chapter 5], or [12, Chapter 3].
The system (1) possesses a unique continuous solution R + → R d , t → X x (t) for every x ∈ R d , cf. e.g. [13,Theorem 3.6 is another solution, then P{X x (t) = Y x (t)} = 1 for all t ≥ 0. Note especially that for every x ∈ R d the solution X x (t) is defined for all t ≥ 0. Further, X 0 (t) = 0 for all t ≥ 0, i.e. the origin is an equilibrium of the system.
We are interested in the stability of the equilibrium at the origin for the system (1). There are various stability concepts in the literature for the zero solution of (3). The ones we use in this paper are listed in the next definition. For corresponding definitions in the literature cf. e.g. [ Definition 2. Consider the system (3), which includes (1) as a special case. The origin (or the zero solution) is said to be : i) Stable in probability (SiP), if for every pair ε, r > 0 there exists a δ = δ(ε, r) > 0 such that ii) Asymptotically stable in probability (ASiP), if it is SiP and for every ε > 0 there exist a δ = δ(ε) > 0 such that iii) Globally asymptotically stable in probability (GASiP), if it is SiP and iv) pth moment exponentially stable or exponentially p-stable (p-ES) for a p > 0, if there exist constants A, α > 0 such that for every x ∈ R d and every t ≥ 0.

Remark 3.
There are numerous other stability concepts for the zero solution in the literature, e.g. p-stability, asymptotic p-stability, and almost sure exponential stability. Further, SiP is sometimes referred to as strongly stable in probability [12] or stochastically stable [13,12]. Similarly an ASiP zero solution is sometimes said to be stochastically asymptotically stable [13], a GASiP zero solution is said to be stochastically asymptotically stable in the large [13], or (asymptotically) stable in the large [12]. 1-ES and 2-ES are often referred to as exponential stability in the mean and exponential stability in mean square, respectively.
For the autonomous linear SDE with constant coefficients (1) the relation between the stability concepts in Definition 2 is simpler than in the general case. Proof. Proposition iii) follows by [12,Theorem 6.5] and that ASiP follows from GASiP is trivial. Propositions ii) and iii) together imply that GASiP follows from ASiP. Proposition ii) is most conveniently proved by the use of so-called Lyapunov functions, which we discuss in the next section. We will prove ii) in Remark 6 after that discussion.
Remark 4. Although fractional moments are interesting in their own right, cf. e.g. [17,15], our sole interest in p-ES is to use it as a stepping stone to prove GASiP for linear systems using Proposition 1.
3. Lyapunov functions for SDE. Just as for deterministic systems defined by ordinary differential equations (ODE), Lyapunov functions are an important tool when studying stability in SDEs. A Lyapunov function for a globally asymptotically stable (GAS) equilibrium of an ODE is a continuous function that has a local minimum at the equilibrium and is strictly decreasing along solution trajectories of the system. Indeed, an isolated equilibrium of an ODE is GAS, if there exists a smooth, radially unbounded Lyapunov function for the system. A Lyapunov function for a stochastic system is a continuous function that is a supermartingale [4].
for a stochastic system. A sufficient condition for a smooth Lyapunov function to imply GAS of the origin of the cf. [14]. For an SDE the orbital derivative which becomes for the linear system (1) with constant coefficients.
Remark 5. For deterministic systems there is no restriction in considering smooth Lyapunov functions. For stochastic systems this is not the case and it is much too restrictive to demand that a Lyapunov function for a stochastic system is smooth at the origin. In fact, a smooth Lyapunov function often does not exist, cf. [12,Remark 5.5]. An appropriate function class for Lyapunov functions is given in [12, p. 146]. This important fact makes the discussion on stability in (the otherwise excellent book) [13, Chapter 4] somewhat incomplete, since only C 2 Lyapunov functions are studied. For our system (1) it can be concluded from the discussion in [12, p. 146] that the appropriate class for Lyapunov functions for system (1) is given by We conclude this section by giving three theorems connecting the stability concepts from Definition 2 to the existence of Lyapunov functions.
Then the origin is ASiP.
Then the origin is p-ES.
The last Lyapunov function theorem is a so-called converse theorem, i.e. it asserts the existence of a Lyapunov function for the system (1) if the origin is ASiP.
Proof. Since the origin is ASiP for system (1) it is p-ES for all small enough p > 0 by Proposition 1 iii). The proposition now follows from [12, Theorem 6.2].
In the next section we discuss how the results above can be utilized to explicitly compute and verify Lyapunov functions for the system (1).

Construction and verification of a Lyapunov function.
for all x ∈ R d \ {0}. The next lemma gives a useful formula for LV (x). where Proof. The lemma is shown by direct calculation. Assume first that p = 2. Since for r, s ∈ [1, d] Z . Putting this into the right-hand side expression in (7) gives the following formula for LV (x): Simplifying this formula term-by-term gives 2 r,s,i and r,s,k,l Adding up the terms delivers (8). It is easy to verify that (8) If H(x) ≥ c x 4 , where c > 0 is a constant, then V (x) = x p Q is a Lyapunov function for the system (1), i.e. satisfies the assumptions of Theorem 3.2 and indeed it also fulfils the properties of the function in Theorem 3.3. Note that we could replace x 4 by, e.g., x 2 Q x 2 or x 4 Q without violating these properties.
Before we discuss a general method to deal with the inequality H(x) ≥ c x 4 , let us first consider some simple cases.
We can find a p > 0 such that p(µ + (p − 1)σ 2 /2) < 0, if and only if 2µ < σ 2 . For every such p there is a constant c > 0 such that LV (x) = −c x p . Hence, if 2µ < σ 2 then the origin is p-ES by Theorem 3.2 and thus GASiP by Proposition 1.
Example 9. In [12, Example 6.2] the following case is considered: F +F = −2λI d , Q = I d , and k,l x k x l U u=1 G u rk G u sl = δ r,s x 2 G for some constants λ, G > 0. Here δ r,s is the Kronecker delta, equal to one if r = s and zero otherwise. It should then be concluded 1 that the system (1) is p-ES if G(d + p − 2) < 2λ. By plugging the assumptions into the formulas for LV (x) in Lemma 4.1 or its proof, this result is easily reproduced. In fact, the condition F + F = −2λI d can be mollified to Let us discuss the case p = 2.
Example 10. If we set p = 2 in formula (9) the Lyapunov function candidate V (x) = x 2 Q fulfills condition ii) Theorem 3.3, if and only if This is a linear matrix inequality (LMI) that can be solved for Q 0 using semidefinite programming. Note that since U u=1 (G u ) QG u 0 the LMI (12) is stricter than F Q + QF ≺ 0. For deterministic systems, i.e. when G u = 0 for all u ∈ [1, U ] Z , (12) becomes F Q + QF ≺ 0. It is well known that for any given P 0 the so-called continuoustime Lyapunov equation F Q + QF = −P has a unique solution Q if and only if F is Hurwitz, i.e. all eigenvalues of F have a strictly negative real part. Solving the above matrix equation is thus the standard procedure to compute a quadratic Lyapunov function for the deterministic system. Thus, if white noise of small enough magnitude is added to a deterministic systemẋ = F x, of which the origin is GAS, then the origin is GASiP for the resulting stochastic system. The existence of the Lyapunov function V (x) = x p Q is a sufficient condition for p-ES of the zero solution. However, if we are interested in showing GASiP, then p-ES with any p > 0 is sufficient. Hence, rather than focussing on p ≥ 2, we will cover more cases by trying to establish a Lyapunov function for small p > 0, see the following remark.
Remark 12. From (9) it is obvious that if H(x) ≥ 0 for a particular p * > 0, then it is nonnegative for every p, 0 < p < p * . Since GASiP follows from p-ES for any p > 0 there is no reason to go beyond p = 2 if one wants to assert GASiP.

4.2.
Sum-of-squares. In light of Remark 7 we would like to verify that there exists a c > 0 such that Note that P c (x) is a homogenous polynomial of degree 4. The study of the nonnegativity of polynomials is a very mature field that has been extensively researched and there are efficient algorithms. They use semi-definite programming in the framework of sum-of-squares (SOS) polynomials to verify if polynomials are nonnegative, cf. e.g. [1,8] for an overview of SOS for computing Lyapunov functions for deterministic systems. For other approaches than SOS programming for polynomial Lyapunov functions, cf. e.g. [11]. polynomials P 1 , P 2 , . . . , P k : R d → R such that Remark 14. If P : R d → R is SOS, then P is nonnegative, i.e. P (x) ≥ 0 for all x ∈ R d . The converse is in general not true, i.e. there are polynomials that are nonnegative but cannot be written as the sum of squared polynomials. The Motzkin polynomial x 4 y 2 +x 2 y 4 −3x 2 y 2 +1 is an example of a nonnegative polynomial which is not a sum-of-squares polynomial [16].
The following theorem shows that the nonnegativity of the polynomial P c implies the existence of a Lyapunov function. Moreover, the nonnegativity of P c (x) is equivalent to x 2m P c (x) being SOS for sufficiently large m.
Proof. To prove i) note that there are constants 0 < C min and then Proposition ii) uses the fact that P c (x) is a homogeneous polynomial of degree 4. It follows from the famous paper [10], where David Hilbert showed that if d = 3 then every nonnegative P c (x) can be written as the sum of three squared polynomials. For d = 2 every nonnegative P c (x) can even be written as the sum of two squared polynomials.
Proposition iii) follows directly by [20,21]. Indeed, the proposition holds true for all homogeneous polynomials P c (x) of an even degree and x m P c (x) can even be written as the sum of squared monomials.
To establish that the polynomial P c (x) is a sum of squares, we use SOS programming to compute a matrix S 0, such that P c (x) = Z SZ, where Z is a vector of a basis of homogeneous polynomials of degree 2. That is, the elements of Z are of the form x i x j , i, j ∈ [1, d] Z . It is not difficult to see that the length m of Z is d(d + 1)/2.
We seek to find a matrix S 0 such that P c (x) = Z SZ (note that S in the previous equation is not unique). Indeed, then P c (x) is SOS as S = O DO, where D ∈ R m×m is a diagonal matrix with entries D 11 , D 22 , . . . , D mm ≥ 0 on the diagonal and O ∈ R m×m is orthogonal, the entries Y i of the vector Y = OZ are polynomials in x 1 , x 2 , . . . , x d and The SOS method is available in the Matlab toolbox SOSTOOLS and returns D and O; it also allows us to optimize parameters which appear linearly. We use SOS programming to write as SOS.
Ideally we would like to find all parameters, in particular Q, such that P c (x) is SOS. Unfortunately, if p = 2 the formula for P c (x) is nonlinear in Q and we can thus not use the SOS program to compute an "optimal" Q. Thus, we have to fix Q and use SOS to verify whether P c (x) is SOS and, thus, V (x) = x p Q is a Lyapunov function and the zero solution is p-ES. This limitation is discussed in more detail in Section 6.
We can, however, let e.g. p or c be variables and let the SOS program minimize an arbitrary linear function of them. Given Q, we can, e.g., compute the largest p, such that the zero solution is p-ES. Further, if our problem has more structure, e.g. U = 1 and G 1 = σ 1 I d , then the formula is linear in σ 2 1 and we can use the SOS program to minimize σ 2 1 , i.e., compute the least amount of noise needed to stabilize the zero solution. We will give examples of these computations in the next section.

5.
Examples. We used Matlab R2014a with the Matlab toolbox SOSTOOLS v. 3.00 [19] with SeDuMi v. 1.3 [22] to formulate and solve the SOS optimization problems. In Subsection 5.1 we study a harmonic oscillator with noisy damping and can show GASiP for a range of parameter values. In Subsection 5.2 we discuss stabilization by noise and show GASiP as well as p-ES.

5.1.
Harmonic oscillator with noisy damping. We consider a damped harmonic oscillatorẍ + kẋ + ω 2 x = 0. Setting x 1 = ωx and x 2 =ẋ we geṫ The eigenvalues of F are λ ± = (−k ± √ k 2 − 4ω 2 )/2 and the origin is GAS for k > 0. Now let us add white noise to the damping and consider the SDE This example is similar to [12, Example 6.6], but we consider the solution in the Itôand not the Stratonovich sense, hence, we give all details below and obtain slightly different results.
Khasminskii [12,Example 6.6] states that "It is extremely interesting to study 'bifurcation' values of the noise intensity, i.e., values for which the system first becomes unstable." Following his suggestions, we are interested in studying large σ for which the zero solution of (14) is still GASiP. First, we follow the standard approach by setting p = 2 and making the ansatz LV (x) = − x 2 . Equations (8) and (9) then give and by setting Q 11 = A, Q 12 = Q 21 = B, and Q 22 = C and comparing coefficients i.e.
By Sylvester's criterion Q is positive definite, if and only if is a Lyapunov function for system (14) asserting 2-ES stability (exponential mean square stability). From Remark 11 it follows that (14) is 2-ES, if and only if k > σ 2 /2. Hence, we easily verify GASiP for parameters where the system is 2-ES. Now we seek to verify GASiP for various parameters where the system is not 2-ES, namely for k < σ 2 /2. This is obtained by proving that the origin is p-ES for some 0 < p < 2. We choose the matrix Q = 3 1/3 1/3 3 , with eigenvalues λ = 3 ± 1/3 > 0, and we compute S ∈ R 3×3 , S 0, such that P c (x) = Z SZ with Z = x 2 1 , x 1 x 2 , x 2 2 as explained above. This is remarkably easy to achieve with SOSTOOLS, after initializing the matrices F, G, and Q and the constants c and p the commands The data in Table 1 can be interpreted in the following way: For example, the information from experiment #1 in the table is that for k = 1.5, σ = 2, ω = 3, p = 0.5, and c = 1.6875, we can write Note that in experiment number 4, where p = 1.2, even for c = 0 SOSTOOLS was not able to write the polynomial P 0 (x) as sum-of-squares. This means that we were  Table 1. Results of checking whether P c (x) for system (14) can be written as SOS. No solution means that even for c = 0 SOSTOOLS was not able to write P c (x) as SOS. In all the experiments we set σ = 2.0. In experiments #1 to #4 we set k = 1.5 and in experiments #5 to #9 we set k = 0.9. not able to show p-ES for this set of parameters with our ansatz for V and Q, but it might be possible to show it with a different matrix Q. Also, for experiment 9, we have decreased ω to 2.5, and we were not able to write the polynomial P 0 (x) as sum-of-squares with these parameters.
Summarizing, our method enables us to show GASiP for parameters with k < σ 2 /2, where the zero solution is not 2-ES.

Stabilization by noise.
It is well known that the unstable zero solution of a deterministic system can be stabilized by white noise. In particular, for any matrix F in (1) the system can be stabilized by a large enough white noise. For example, with λ as the largest eigenvalue of the symmetric matrix F + F , one can take u = 1 and G 1 = σI d with σ 2 > λ. Plugging this into (9) For our further discussion the following lemma is useful. Recall that the condition number (with respect to the Euclidian norm) of a positive definite matrix G is defined as cond(G) := λ max /λ min , where λ max , λ min denote the largest, smallest eigenvalue of G, respectively. Lemma 5.1. Assume that G ∈ R d×d is a positive definite matrix such that cond(G) < 2. Then for any 0 < p < 2 − cond(G) and we have for all Write an arbitrary x ∈ R d as x = ix i o i and setx = ( Putting this together yields, noting that (2 − p)λ d − λ 1 > 0 by (16): using again (16).
The following proposition considers the case of a deterministic systemẋ = F x to which noise is added. In particular we can always stabilize the system by applying sufficiently large white noise. This is of course particularly interesting if the deterministic system is unstable. Statement a) in the proposition is well known, cf. e.g. [13,Example 4.2.5], and statements b) and c) are extensions of a) using Lemma 5.1.

Proposition 2.
Let λ be the largest eigenvalue of the symmetric matrix F + F and consider the system (1). a) Let G u = σ u I d for all u ∈ [1, U ] Z and denote σ 2 := U u=1 σ 2 u . Then the origin is p-ES for 0 < p < 1 if λ < (1 − p)σ 2 and it is GASiP if Then the origin is p-ES if λ < γ p := U u=1 γ p,u with γ p,u := G u 2 The origin is Using (8) and (9), we have a) We get for every 0 < p < 1 that Since the origin is GASiP if it is p-ES for some p > 0, we see that it is GASiP if λ < σ 2 . b) We have, using Lemma 5.1 It follows that the origin is p-ES if λ < γ p and GASiP if λ < γ 0 . c) We have G u 2 = σ u G u 2 and cond(G u ) = cond( G u ) and the statement follows from b).
Let us now discuss how our SOS program can improve these results: In the simplest case a) our SOS approach does not add anything. The conditions in b), however, are much too restrictive, which is explored in the following three examples. In all the examples we set Q = I d .
Hence, λ > γ p for any p > 0 and the conditions in Proposition 2 b) are not satisfied for any p > 0. However, e.g. with p = 0.3, our SOS program has a solution with Z = x 2 1 , x 1 x 2 , and we can conclude that the origin is 0.3-ES and GASiP. We can use our method to maximize p, while fixing the other parameters. By fixing c = 0.1 and maximizing p we get the solution p = 0.5095. From this we can conclude that the origin is 0.5095-ES.
In the case c) of the previous proposition, we can use SOS to find a small amount of noise such that the origin is stabilized. This is achieved by minimizing σ 2 1 . Example 17. We set c = p = 0.1 and F as in Example 16, and set G 1 = σ 1 G 1 with G 1 equal to the G 1 in Example 16. In particular, we have cond( G 1 ) = 1.75 < 2.
We are interested in the minimal σ 2 1 such that we can show p-ES. From (17) we get the sufficient bound σ 2 1 = 3.8441 for 0.1-ES of the origin. The SOS program is linear in σ 2 1 and by minimizing it we get σ 2 1 = 0.5164. This shows that the origin is 0.1-ES for σ 2 1 = 0.5164, a significantly smaller value than obtained by Proposition 2.
6. Conclusions and further research. In this paper, we have studied global asymptotic stability (GASiP) and pth moment exponential stability (p-ES) for autonomous linear SDEs using a Lyapunov function. We have used an ansatz for the Lyapunov function, so that its verification can be achieved by studying the nonnegativity of a polynomial. We have developed an SOS programming problem to verify the nonnegativity; the program can also be used to optimize certain parameters, such as maximizing p or minimizing a noise parameter. We have applied the method to several examples and have shown that it can show stability in previously undecided cases. Further work will include studying the stability of nonlinear autonomous SDEs, for which the Lyapunov function of the linearized system at the origin can serve as a local Lyapunov function. Furthermore, it would be desirable to compute a suitable matrix Q such that V (x) = x p Q is a Lyapunov function. As Q does not appear linearly, our proposed method cannot achieve this, but one could consider the use of bilinear matrix inequality (BMI) optimization for solving (13) for Q and thus finding a Lyapunov function. Another possibility is to try to reformulate our approach using Schur complements and design an LMI that is also able to search for a suitable Q. For references to such methods see [2,25].