Stability of cyber-physical systems of numerical methods for stochastic differential equations: integrating the cyber and the physical of stochastic systems

TL;DR

This paper introduces a unified cyber-physical system framework for stochastic differential equations and numerical methods, establishing a fundamental stability theory linking the stability of the SDE and its numerical approximation.

Contribution

It develops a general class of stochastic impulsive differential equations to model CPSs of numerical methods, providing a new theoretical foundation for stability analysis.

Findings

CPS framework unifies SDEs and numerical methods

Lyapunov stability theory applied to CPSs

Necessary and sufficient stability conditions for Euler-Maruyama

Abstract

This paper presents the cyber-physical system (CPS) of a numerical method (the widely-used Euler- Maruyama method) and establishes a foundational theory of the CPSs of numerical methods for stochastic differential equations (SDEs), which transforms the way we understand the relationship between the numerical method and the underlying dynamical system. Unlike in the literature where they are treated as separate systems linked by inequalities, the CPS is a seamless integration of the SDE and the numerical method and we construct a new and general class of stochastic impulsive differential equations (SiDEs) that can serve as a canonic form of the CPSs of numerical methods. By the CPS approach, we show the equivalence and intrinsic relationship between the stability of the SDE and the stability of the numerical method using the Lyapunov stability theory we develop for our class of SiDEs. Applying our established theory, we present the CPS Lyapunov inequality that is the necessary and sufficient condition for meansquare stability of the CPS of the Euler-Maruyama method for linear SDEs. The proposed CPS and theory initiate the study of systems numerics and provoke many open and interesting problems for future work.

Equations

$$\mathrm{d}x(t)=f(x(t),t)\mathrm{d}t+g(x(t),t)\mathrm{d}B(t)$$

$$\mathrm{d}y(t)=\tilde{f}(x(t),y(t),t)\mathrm{d}t+\tilde{g}(x(t),y(t),t)\mathrm{d}B(t),t\in [t_k,t_{k+1})$$

$$\Delta (x(t_{k+1}^{-}),k+1):=x(t_{k+1})-x(t_{k+1}^{-})$$

$$=h_f(x(t_{k+1}^{-}),k+1)+h_g(x(t_{k+1}^{-}),k+1)\xi (k+1)$$

$$\tilde{\Delta }(x(t_{k+1}^{-}),y(t_{k+1}^{-}),k+1):=y(t_{k+1})-y(t_{k+1}^{-})$$

$$={\tilde{h}}_f(x(t_{k+1}^{-}),y(t_{k+1}^{-}),k+1)+{\tilde{h}}_g(x(t_{k+1}^{-}),y(t_{k+1}^{-}),k+1)\xi (k+1),k\in \mathbb{N}$$

$$|f(x,t)-f(\bar{x},t)|\vee |g(x,t)-g(\bar{x},t)|\vee |h_f(x,k)-h_f(\bar{x},k)|\vee |h_g(x,k)-h_g(\bar{x},k)|\le L|x-\bar{x}|$$

$$\begin{gathered} |\tilde{f}(x,y,t)-\tilde{f}(\tilde{x},\tilde{y},t)|\vee |\tilde{g}(x,y,t)-\tilde{g}(\tilde{x},\tilde{y},t)| \\ \vee |{\tilde{h}}_f(x,y,k)-{\tilde{h}}_f(\tilde{x},\tilde{y},k)|\vee |{\tilde{h}}_g(x,y,k)-{\tilde{h}}_g(\tilde{x},\tilde{y},k)|\le \tilde{L}(|x-\tilde{x}|\vee |y-\tilde{y}|) \end{gathered}$$

$$\mathcal{L}V(x,t)=V_t(x,t)+V_x(x,t)f(x,t)+\frac{1}{2}\mathrm{trace}\left[g^T(x,t)V_{xx}(x,t)g(x,t)\right],$$

$$\tilde{\mathcal{L}}\tilde{V}(x,y,t)={\tilde{V}}_t(y,t)+{\tilde{V}}_y(y,t)\tilde{f}(x,y,t)+\frac{1}{2}\mathrm{trace}\left[{\tilde{g}}^T(x,y,t){\tilde{V}}_{yy}(y,t)\tilde{g}(x,y,t)\right].$$

$$\mathrm{d}z(t)=F(z(t),t)\mathrm{d}t+G(z(t),t)\mathrm{d}B(t)t\in [t_k,t_{k+1})$$

$$z(t_{k+1})-z(t_{k+1}^{-})=H_F(z(t_{k+1}^{-}),k+1)+H_G(z(t_{k+1}^{-}),k+1)\xi (k+1)k\in \mathbb{N}$$

$$\begin{gathered} F(z,t)=\left[\begin{array}{c}f\left(Cz,t\right)\\ \tilde{f}\left(Cz,Dz,t\right)\end{array}\right],G(z,t)=\left[\begin{array}{c}g\left(Cz,t\right)\\ \tilde{g}\left(Cz,Dz\right),t\end{array}\right], \\ {H}_F(z,k)=\left[\begin{array}{c}{h}_f\left(Cz,k\right)\\ {\tilde{h}}_f\left(Cz,Dz,k\right)\end{array}\right],H_G(z,k)=\left[\begin{array}{c}{h}_g\left(Cz,k\right)\\ {\tilde{h}}_g\left(Cz,Dz,k\right)\end{array}\right]. \end{gathered}$$

$$\begin{gathered} |F(z,t)-F(\tilde{z},t)|\vee |G(z,t)-G(\tilde{z},t)| \\ \vee |H_F(z,k)-H_F(\tilde{z},k)|\vee |H_G(z,k)-H_G(\tilde{z},k))|\le L_z|z-\tilde{z}| \end{gathered}$$

$$c_1|x{|}^p\le V(x,t)\le c_2|x{|}^p,$$

$${\tilde{c}}_1|y{|}^p\le \tilde{V}(y,t)\le {\tilde{c}}_2|y{|}^p$$

$$\mathcal{L}V(x,t)\le -\alpha V(x,t),$$

$$\tilde{\mathcal{L}}\tilde{V}(x,y,t)\le {\tilde{\alpha }}_{1}V(x,t)+{\tilde{\alpha }}_2\tilde{V}(y,t),t\in [t_k,t_{k+1})$$

$$\displaystyle{\mathbb{E}}\big{[}V(x+\Delta(x,k+1),t)\big{|}x\big{]}\leq\beta V(x,t),$$

$$\displaystyle{\mathbb{E}}\big{[}\tilde{V}(y+\tilde{\Delta}(x,y,k+1),t)\big{|}x,y\big{]}\leq\tilde{\beta}_{1}V(x,t)+\tilde{\beta}_{2}\tilde{V}(y,t)$$

$$\frac{\ln \beta }{\alpha }<\underline{\Delta t}\le \overline{\Delta t}<\frac{-\ln {\tilde{\beta }}_2}{{\tilde{\alpha }}_2}.$$

$$\mathbb{E}V(x(t),t)\le V((0),0)e^{-\bar{\alpha }t}\forall t\ge 0$$

$$\mathbb{E}\tilde{V}(y(t),t)=\mathbb{E}\tilde{V}(y(\tilde{t}),\tilde{t})+{\int }_{\tilde{t}}^t\mathbb{E}\tilde{\mathcal{L}}\tilde{V}(x(s),y(s),s)\mathrm{d}s$$

$$\mathbb{E}\tilde{\mathcal{L}}\tilde{V}(y(t),t)\le {\tilde{\alpha }}_1\mathbb{E}V(x(t),t)+{\tilde{\alpha }}_2\mathbb{E}\tilde{V}(y(t),t)$$

$$({\tilde{\beta }}_2+\delta )e^{({\tilde{\alpha }}_2+\delta +\bar{\delta })\overline{\Delta t}}\le 1.$$

$$\mathbb{E}\tilde{\mathcal{L}}\tilde{V}(y(t),t)\le ({\tilde{\alpha }}_2+\delta )\mathbb{E}\tilde{V}(y(t),t)$$

$$\mathbb{E}\tilde{V}(y(t),t)\ge \frac{{\tilde{\alpha }}_1}{\delta }\mathbb{E}V(x(t),t).$$

$$\mathbb{E}\tilde{V}(y(t_{k+1}),t_{k+1})\le ({\tilde{\beta }}_2+\delta )\mathbb{E}\tilde{V}(y(t_{k+1}^{-}),t_{k+1}^{-})$$

$$\mathbb{E}\tilde{V}(y(t_{k+1}^{-}),t_{k+1}^{-})\ge \frac{{\tilde{\beta }}_1}{\delta }\mathbb{E}V(x(t_{k+1}^{-}),t_{k+1}^{-}).$$

$$\mathbb{E}\tilde{V}(y(t),t)\le \tilde{V}(y(0),0)+{\tilde{\alpha }}_2{\int }_0^t\mathbb{E}\tilde{V}(y(s),s)\mathrm{d}s\forall t\in [0,t_1).$$

$$\mathbb{E}\tilde{V}(y(t),t)\le \tilde{V}(y(0),0)e^{{\tilde{\alpha }}_{2}t}\forall t\in [0,t_1)$$

$$\begin{gathered} \mathbb{E}\tilde{V}(y(t_1),t_1)\le {\tilde{\beta }}_2\mathbb{E}\tilde{V}(y(t_1^{-}),t_1^{-})\le {\tilde{\beta }}_2\tilde{V}(y(0),0)e^{{\tilde{\alpha }}_2{t}_1} \\ <\tilde{V}(y(0),0)e^{-({\tilde{\alpha }}_2+\delta +\bar{\delta })\overline{\Delta t}}{e}^{{\tilde{\alpha }}_2{t}_1}\le \tilde{V}(y(0),0)e^{-(\delta +\bar{\delta })\overline{\Delta t}}. \end{gathered}$$

$$\mathbb{E}\tilde{V}(y(t),t)\le \tilde{V}(y(0),0)e^{{\tilde{\alpha }}_2(t-t_k)-k(\delta +\bar{\delta })\overline{\Delta t}}$$

$$\mathbb{E}\tilde{V}(y(t),t)\le e^{({\tilde{\alpha }}_2+\delta +\bar{\delta })\overline{\Delta t}}\tilde{V}(y(0),0)e^{-(\delta +\bar{\delta })t}\forall t\ge 0.$$

$$\bar{v}(t)=\frac{({\tilde{\alpha }}_1\vee {\tilde{\beta }}_1)}{\delta }\mathbb{E}V(x(t),t)-\mathbb{E}\tilde{V}(y(t),t)\forall t\in [0,\infty )$$

Full text

1 Introduction

According to Newton’s second law of motion, we describe a mechanical system with differential equations. Usually, physiccal laws are expressed by means of differential equations, and so are the models of dynamical systems in many disciplines, ranging from biology to finance. Such models play a central role in all scientific and engineering disciplines [9, 35]. A model may serve many purposes. The value of a model depends on the model fidelity, where the model of a dynamical system is said to be have high fidelity if it accurately describes the properties of the system [9, 35, 36]. Studying the model of high fidelity gives us insight into how the dynamical system will behave in the real world. Generally, a dynamical system, ranging from the motion of pollen particles to the movement of stock price, is subject to intrinsic and/or extrinsic noise in the real world [22, 24, 39]. Such randomness should/must be taken into account by a model of high fidelity if it matters, say, it affects some property of the system that is of concern to the modelling. If we allow for some noise in some coefficients of a differential equation, we often obtain a more realistic model of the situation that is able to describe the fluctuations observed in the physical system. This leads to modelling with stochastic differential equations (SDEs). The study of SDEs can be seen to have started from the classical paper of Einstein that presented a mathematical connection between microscopic Brownian motion of particles and the macroscopic diffusion equation, and the interest in SDE has grown enormously in the last few decades [1, 31, 44, 49, 54, 56].

Over the past a few decaseds, stochastic systems described by SDE have been intensively studied since stochastic modelling has come to play a significant role in many branches of science and engineering [2, 23, 25, 27, 28]. It is hardly possible to solve an SDE analytically and have the exact solution of the SDE. For practical purposes, numerical approximations to the exact solution are usually obtained, which, called the numerical solutions, are discrete-time stochastic processes produced by some numerical methods. Such numerical schemes, in the form of stochastic difference equations, are the translations of the SDE into discretization. Practically, computers are used to excute the numerical methods and generate the numerical solutions of the SDE, from which one could learn and/or infer some dynamical properties of the underlying physical system [4, 18, 27, 51].

As is well known, whenever a computer is used in measurement, computation, signal processing or control applications, the data, signals and systems involved are naturally described as discrete-time processes [3, 26, 27, 50]. It is worth noting that the SDE is the physical model which represents our knowledge of the physical system and a numerical method is a cyber model which is a representative of the physical model in computers, the cyber world. The physical model, namely, the SDE often refers to the phyiscal system (particularly, which is an engineered system) itself while its cyber couterpart, namely, the numerical method symbolizes it in the cyber world. In the age of networking and information technology, the cyber model plays a key role in understanding and controlling the underlying physical system, which not only envisions the approximate behaviour of the physical system [18, 19, 20] but is also utilized to extract knowledge of the system from data [27, 50] and based on which control is designed and implemented [2, 3, 15, 26, 29, 43]. It is natural and imperative

• (I)

to find out the relationship between the physical model (i.e., the SDE) and its cyber counterpart (i.e., the numerical method) of a dynamical system;

• (II)

and to ensure that they both share some dynamical properties such as stability which is concerned in the study.

The principal aim of this paper is to address the problems (I) and (II) of fundamental importance in the age of networking and information technology. As a matter of fact, the fundamental importance of these problems has been recognized and they have been studied in a vast literature. Results that address the problems can be found in those many on convergence and stability of numerical methods for SDEs, where the SDE and the numerical method are treated as separate systems which are linked by inequalities in some moment sense on any finite time interval [19, 20, 30, 32, 41, 46]. The ability of a cyber system to reproduce the stability of its underlying physical system can be found in a large number of works. The problem how to reproduce the stability of the physical system in its cyber counterpart, which is called the test problem, has also been studied for SDE [17, 20, 40, 47]. The key question in a test problem is [17]

• (Q1)

for what stepsizes $\Delta t$ does the cyber system (the numerical method) share the stability property of the underlying physical system (the SDE)?

This naturally provokes the converse question [20, 40]

• (Q2)

does the stability of the cyber system (the numerical method) for small stepsizes $\Delta t$ imply that of the underlying physical system (the SDE)?

These questions deal with asymptotic ($t\to\infty$) properties and hence cannot be answered directly by applying traditional finite-time convergence results [20, 40]. Results that answer questions (Q1) and (Q2) can be found in the literature. For example, results for scalar linear systems were given in [17, 47]. For multi-dimensional systems with global Lipschitz condition, Higham, Mao and Stuart [20] introduced a natural finite-time strong convergence condition [20, Condition 2.3], which links a cyber system with its underlying physical system by an inequality in some moment sense over any finite time interval, and proved that there is a sufficiently small $\Delta t^{\ast}>0$ such that, for every $\Delta t\in(0,\Delta t^{\ast}]$, the mean-square exponential stability of the physical system is equivalent to that of its cyber counterpart [20, Lemmas 2.4-2.5 and Theorem 2.6]. Recently, Mao [40] developed new techniques to handle the small $p$th moment ($p\in(0,1)$) and showed that, under a natural finite $p$th moment condition [40, Assumption 2.4] and a natural finite-time convergence condition [40, Assumption 2.5], the $p$th moment exponential stability of the physical system is equivalent to that of its cyber counterpart for every $\Delta t\in(0,\Delta t^{\ast}]$ with some sufficiently small $\Delta t^{\ast}>0$ [40, Lemmas 2.6-2.7 and Theorem 2.8]. As is pointed out in [40], there are many open problems in this research. For instance, although the existence of the (sufficiently small) upper bound $\Delta t^{\ast}>0$ of stepsizes has been shown [20, 40], it is severely limited by the growth constant of the exponentially stable system, which refers to the physical system and its cyber counterpart when answering (Q1) and (Q2), respectively. Recall that, though the growth and the rate constants are related, it is only the rate constant that counts in the definition of exponential stability, see Definitions 1 and 2 below. It appears that, either to reproduce or to imply the exponential stability of the physical system by its cyber counterpart, the condition imposed on the stepsizes that explicitly depends on the growth constant could/should be relaxed [51]. This could significantly improve the upper bound $\Delta t^{\ast}$ of stepsizes and facilitate the computation.

It is noted that the physical system (the SDE) and its cyber counterpart (the numerical method) are bound by inequalities in the literature [19, 20, 30, 32, 41, 46, 47, 51]. Nevertheless, they remain as two systems, largely separate. This paper constructs the cyber-physical model of a dynamical system that is a seamless, fully synergistic integration of the physical model (the SDE) and its cyber counterpart (the numerical method). Here we present a new and general class of stochastic impulsive differential equations (SiDEs) which can be used to represent the integrated dynamics of the physical system and its cyber counterpart. Impulsive differential equations, also known as impulsive systems, have been studied for several decades [8, 16, 28, 45, 48, 57]. But these impulsive systems in the literature are just the physical subsystems in our general class of SiDEs, see Section 2. Our proposed SiDE composed of the physical subsystem and the cyber subsystem is formulated as a canonic form of cyber-physical systems (CPSs), which presents a systematic framework for the study of CPSs [10, 33]. The canonic form not only provides a holistic view of CPSs but also reveals the intrinsic relationship between the physical subsystem and the cyber subsystem. In the study of numerical analysis, we present the CPS of a numerical method (the widely-used Euler-Maruyama method) for SDEs, which is in the form of our proposed SiDEs and represents a seamless integration of the SDE and the numerical method. The SDE and the numerical method are the physical subsystem and the cyber subsystem of the CPS, respectively. From the viewpoint of cybernetics [55], an essential problem to study is whether and how the CPS reproduces some dynamical properties such as the stability of its physical or cyber subsystem since ‘the primary concern of cybernetics is on the qualitative aspects of the interrelations among the various components of a system and the synthetic behavior of the complete mechanism’ [52]. Using the terminology of the CPSs, we rephrase the key questions (Q1) and (Q2) as follows.

(Q1) For what stepsizes $\Delta t$ do the CPS and, hence, the cyber subsystem reproduce the stability property of the physical subsystem?

(Q2) Does the stability of the cyber subsystem for small stepsizes $\Delta t$ imply that of the CPS and, hence, that of the physical subsystem?

To address the key questions, we develop a Lyapunov stability theory for our proposed general class of SiDEs and apply it to the CPS of the Euler-Maruyama method for SDEs. Applying our established stability theory, we prove positive results to the two key questions (Q1) and (Q2). This establishes a foundational theory for the CPSs of numerical methods for SDEs. As application of our established theory, we present the CPS Lyapunov inequality, the necessary and sufficient condition for mean-square exponential stability of the CPS of the Euler-Maruyama method for linear SDEs. In this paper, we initiate the study of systems numerics and there are many interesting and/or challenging problems for future work.

2 A general class of stochastic impulsive differential equations

Throughout this paper, unless otherwise specified, we shall employ the following notation. Let us denote by $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},{\mathbb{P}})$ a complete probability space with a filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions (i.e. it is right continuous and $\mathcal{F}_{0}$ contains all ${\mathbb{P}}$-null sets) and by ${\mathbb{E}}[\cdot]$ the expectation operator with respect to the probability measure. Let $B(t)=\begin{bmatrix}B_{1}(t)&\cdots&B_{m}(t)\end{bmatrix}^{T}$ be an $m$-dimensional Brownian motion defined on the probability space. If $x,y$ are real numbers, then $x\vee y$ denotes the maximum of $x$ and $y$, and $x\wedge y$ stands for the minimum of $x$ and $y$. If $A$ is a vector or a matrix, its transpose is denoted by $A^{T}$. If $P$ is a square matrix, $P>0$ (resp. $P<0$) means that P is a symmetric positive (resp. negative) definite matrix of appropriate dimensions while $P\geq 0$ (resp. $P\leq 0$) is a symmetric positive (resp. negative) semidefinite matrix. Let $\lambda_{M}(\cdot)$ and $\lambda_{m}(\cdot)$ be a matrix’s eigenvalues with maximum and minimum real parts, respectively. Denote by $|\cdot|$ the Euclidean norm of a vector and the trace (or Frobenius) norm of a matrix.

Let $C^{2,1}({\mathbb{R}}^{n}\times{\mathbb{R}}_{+};{\mathbb{R}}_{+})$ be the family of all nonnegative functions $V(x,t)$ on ${\mathbb{R}}^{n}\times{\mathbb{R}}_{+}$ that are continuously twice differentiable in $x$ and once in $t$. Let ${\mathcal{M}}^{p}([a,b];{\mathbb{R}}^{n})$ be the family of ${\mathbb{R}}^{n}$-valued adapted process $\{x(t):a\leq t\leq b\}$ such that ${\mathbb{E}}\int_{a}^{b}|x(t)|^{p}{\mathrm{d}}t<\infty$. Denote by $I_{n}$ the $n\times n$ identity matrix and by $0_{n\times m}$ the $n\times m$ zero matrix, or, simply, by [math] the zero matrix of compatible dimensions. Let $\{\xi(k)\}_{k\in{\mathbb{N}}}$, ${\mathbb{N}}=\{0,1,2,\cdots\}$, be an independent and identically distributed sequence with $\xi(k)=\begin{bmatrix}\xi_{1}(k)&\cdots&\xi_{m}(k)\end{bmatrix}^{T}$, and $\xi_{j}(k)$, $j=1,2,\cdots,m$, obeying standard normal distribution while $\{t_{k}\}_{k\in{\mathbb{N}}}$ is a strictly increasing sequence which satisfies $0=t_{0}<t_{1}<t_{2}<\cdots$, $0<\underline{\Delta t}:=\inf_{k\in{\mathbb{N}}}\{t_{k+1}-t_{k}\}\leq\overline{\Delta t}:=\sup_{k\in{\mathbb{N}}}\{t_{k+1}-t_{k}\}<\infty$, and hence $t_{k}\to\infty$ as $k\to\infty$.

Let us consider a stochastic impulsive system described by SiDE

[TABLE]

with initial data $x(0)\in{\mathbb{R}}^{n}$ and $y(0)\in{\mathbb{R}}^{q}$, where $\xi(k+1)$ is independent of $\{x(t),y(t):0\leq t<t_{k+1}\}$, and $f:{\mathbb{R}}^{n}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{n}$, $g:{\mathbb{R}}^{n}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{n\times m}$, $h_{f}:{\mathbb{R}}^{n}\times{\mathbb{N}}\to{\mathbb{R}}^{n}$, $h_{g}:{\mathbb{R}}^{n}\times{\mathbb{N}}\to{\mathbb{R}}^{n\times m}$, $\tilde{f}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{q}$, $\tilde{g}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{q\times m}$, $\tilde{h}_{f}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{N}}\to{\mathbb{R}}^{q}$, $\tilde{h}_{g}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{N}}\to{\mathbb{R}}^{q\times m}$ are measurable functions. To study stability of the system, we assume that they obey $f(0,t)=0$, $g(0,t)=0$, $h_{f}(0,k)=0$ and $h_{g}(0,k)=0$, $\tilde{f}(0,0,t)=0$, $\tilde{g}(0,0,t)=0$, $\tilde{h}_{f}(0,0,k)=0$ and $\tilde{h}_{g}(0,0,k)=0$ for all $t\in{\mathbb{R}}_{+}$ and $k\in{\mathbb{N}}$, and satisfy the global Lipschitz conditions.

Assumption 1**.**

There is $L>0$ such that

[TABLE]

for all $(x,\bar{x})\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}$, $t\in{\mathbb{R}}_{+}$, $k\in{\mathbb{N}}$; and there is $\tilde{L}>0$ such that

[TABLE]

for all $(x,y,\tilde{x},\tilde{y})\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}$, $t\in{\mathbb{R}}_{+}$, $k\in{\mathbb{N}}$.

It is the intersection, interaction and interrelation of the physical system and its cyber counterpart [3, 9, 10, 20, 27, 34, 43, 50] in the age of networking and information technology that motivate our study of stochastic impulsive system (1) , which is constructed as a canonic form of CPSs that is a seamless, fully synergistic integration of the physical system and its cyber counterpart, see also [29]. It is observed in Section 4 that the CPS of a numerical method for SDEs is a special case of (1) where there is no impulse imposed on the physical subsystem (1a). We stress that we delibrately include the impulse (1c) and use the impulsive subsystem (1a,1c) in SiDE (1) to show that the impulsive systems in the literature [8, 16, 28, 45, 48, 57] are just the physical subsystems (1a,1c) in our general class of SiDEs. We construct a system integration (1) of two impulsive subsystems in marked contrast to the impulsive systems in the literature, which highlights the distinction between our new class and those existing.

Remark 1**.**

It should be noticed that, usually, the impulse interval of the subsystem $x(t)$ is a multiple of that of the subsystem $y(t)$ since the former is actually the interval between two consecutive physical impulses imposed on the physical system $x(t)$ while the latter the stepsize of the numerical method, see Section 4. In such a specific case of SiDE (1) where $t_{k}=k\Delta t$, $h_{f}(x,k)$ and $h_{g}(x,k)$ are some functions of $x$ only if $k$ is a multiple of integer $k_{0}>1$, and, otherwise, $h_{f}(\cdot,k)=0$ and $h_{g}(\cdot,k)=0$, where $k_{0}$ is the number of steps of the numerical scheme between two consecutive physical impulses.

Clearly, the trivial solution is the equilibrium of system (1). For a function $V\in C^{2,1}({\mathbb{R}}^{n}\times{\mathbb{R}}_{+};{\mathbb{R}}_{+})$, the infinitesimal generator ${\mathscr{L}}V:{\mathbb{R}}^{n}\times{\mathbb{R}}_{+}\to{\mathbb{R}}$ associated with system (1a) is defined as

[TABLE]

where $V_{t}(x,t)=\frac{\partial V(x,t)}{\partial t}$, $V_{x}(x,t)=\begin{bmatrix}\frac{\partial V(x,t)}{\partial x_{1}}&\cdots&\frac{\partial V(x,t)}{\partial x_{n}}\end{bmatrix}$ and $V_{xx}(x,t)=\begin{bmatrix}\frac{\partial^{2}V(x,t)}{\partial x_{i}\,\partial x_{j}}\end{bmatrix}_{n\times n}$. Similarly, for a function $\tilde{V}\in C^{2,1}({\mathbb{R}}^{q}\times{\mathbb{R}}_{+};{\mathbb{R}}_{+})$, the generator ${\tilde{\mathscr{L}}}\tilde{V}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{R}}_{+}\to{\mathbb{R}}$ associated with system (1b) is defined as

[TABLE]

Let $z(t)=[x^{T}(t)\;\,y^{T}(t)]^{T}\in{\mathbb{R}}^{n+q}$, $C=[I_{n}\;\,0_{n\times q}]$, $D=[0_{q\times n}\;\,I_{q}]$, and, therefore, $x(t)=Cz(t)$, $y(t)=Dz(t)$ for all $t\geq 0$. The stochastic impulsive system (1) can be written in a compact form

[TABLE]

with initial value $z(0)=[x(0)^{T}\;\,y(0)^{T}]^{T}\in{\mathbb{R}}^{n+q}$, where functions $F:{\mathbb{R}}^{n+q}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{n+q}$, $G:{\mathbb{R}}^{n+q}\times{\mathbb{R}}_{+}\to{\mathbb{R}}^{(n+q)\times m}$, $H_{F}:{\mathbb{R}}^{n+q}\times{\mathbb{N}}\to{\mathbb{R}}^{n+q}$, $H_{G}:{\mathbb{R}}^{n+q}\times{\mathbb{N}}\to{\mathbb{R}}^{(n+q)\times m}$ are given by

[TABLE]

Consequently, stochastic impulsive system (4) obeys $F(0,t)=0$, $G(0,t)=0$, $H_{F}(0,k)=0$, $H_{G}(0,k)=0$ for all $t\in{\mathbb{R}}_{+}$, $k\in{\mathbb{N}}$, and, by Assumption 1, it satisfies the global Lipschitz condition, that is, there is a constant $L_{z}>0$ such that

[TABLE]

for all $(z,\tilde{z})\in{\mathbb{R}}^{n+q}\times{\mathbb{R}}^{n+q}$, $t\in{\mathbb{R}}_{+}$, $k\in{\mathbb{N}}$. It is easy to obtain the following result on existence and uniqueness of solutions for SiDE (4).

Lemma 1**.**

Under Assumption 1, there exists a unique (right-continuous) solution $z(t)$ to SiDE (4) (viz., (1)) on $t\geq 0$ and the solution belongs to ${\mathcal{M}}^{2}([0,T];{\mathbb{R}}^{n+q})$ for all $T\geq 0$.

Proof.

Since system (4) satisfies the global Lipschitz continuous condition (5), according to [39, Theorem 3.1, p51], there exists a unique solution $z(t)$ to SiDE (4) on $t\in[t_{0},t_{1})$ and the solution belongs to ${\mathcal{M}}^{2}([t_{0},t_{1});{\mathbb{R}}^{n+q})$. Note that $\xi(k+1)$ is independent of $\{z(t):t\in[t_{0},t_{1})\}$. By virtue of continuity of functions $H_{F}(z,\cdot)$ and $H_{G}(z,\cdot)$ with respect to $z$, there exists a unique solution $z(t_{1})$ to (4) on $t=t_{1}$. Moreover, (4b) and (5) imply that the second moment of $z(t_{1})$ is finite. And, again, according to [39, Theorem 3.1, p51], one has that there is a unique right-continuous solution $z(t)$ to (4) on $t\in[t_{0},t_{2})$ and the solution belongs to ${\mathcal{M}}^{2}([t_{0},t];{\mathbb{R}}^{n+q})$ for $t\in[t_{0},t_{2})$. Recall that $0=t_{0}<t_{1}<t_{2}<\cdots<t_{k}<\cdots$ and $t_{k}\to\infty$ as $k\to\infty$. By induction, one has that there exists a unique (right-continuous) solution $z(t)$ to SiDE (4) for all $t\geq 0$ and the solution belongs to ${\mathcal{M}}^{2}([0,T];{\mathbb{R}}^{n+q})$ for all $T\geq 0$. ∎

Now that we have shown the existence and uniqueness of solutions to SiDE (4), we shall further study the stability of the solutions to the SiDE. Let us introduce the definitions of exponential stability for SiDE (4).

Definition 1**.**

[39, Definition 4.1, p127]** System (4) is said to be $p$th ($p>0$) moment exponentially stable if there is a pair of positive constants $K$ and $c$ such that ${\mathbb{E}}|z(t)|^{p}\leq K|z(0)|^{p}e^{-ct}$ for all $t\geq 0$, which leads to $\limsup_{t\to\infty}\frac{1}{t}\ln({\mathbb{E}}|z(t)|^{p})\leq-c<0$ for all $z(0)\in{\mathbb{R}}^{n+q}$.

Definition 2**.**

[39, Definition 3.1, p119]** System (4) is said to be almost surely exponentially stable if $\limsup_{t\to\infty}\frac{1}{t}\ln|z(t)|<0$ for all $z(0)\in{\mathbb{R}}^{n+q}$.

3 Lyapunov stability theory for the general class of impulsive systems

We dedicate this section to establishing by Lyapunov methods a stability theory for our proposed general class of SiDEs. The general class of SiDEs is formulated as a canonic form of CPSs and we shall develop a foundational theory for stability of CPSs, which will be applied to the CPS of a numerical method for SDEs. In the previous section, for simplicity, the compact form (4) of system (1) is employed to study the existence and uniqueness of solutions to the SiDEs. Now we exploit the structure of SiDE (4) which is composed of subsystems (1a,1c) and (1b,1d) and show the stability of the subsystems as well as that of the whole system (4).

Theorem 1**.**

Suppose that Assumption 1 holds. Let $V\in C^{2,1}({\mathbb{R}}^{n}\times{\mathbb{R}}_{+};{\mathbb{R}}_{+})$ and $\tilde{V}\in C^{2,1}({\mathbb{R}}^{q}\times{\mathbb{R}}_{+};{\mathbb{R}}_{+})$ be a pair of candidate Lyapunov functions for the subsystems (1a,1c) and (1b,1d), respectively, which satisfy

[TABLE]

for all $(x,y,t)\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{R}}_{+}$ and some positive constants $p,c_{1},c_{2},\tilde{c}_{1},\tilde{c}_{2}$. Assume that there are positive constants $\alpha$, $\tilde{\alpha}_{1}$, $\tilde{\alpha}_{2}$, $\beta$, $\tilde{\beta}_{1}$, $\tilde{\beta}_{2}$ such that

[TABLE]

for all $(x,y,t)\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{q}\times{\mathbb{R}}_{+}$ and $k\in{\mathbb{N}}$. Let the impulse time sequence $\{t_{k}\}$ satisfy

[TABLE]

Then SiDE (4) is $p$th moment exponentially stable.

Proof.

According to Lemma 1, that Assumption 1 holds implies there exists a unique solution to SiDE (1a-1d). Let us fix, for simplicity only, any $z(0)=[x(0)^{T}\;\,y(0)^{T}]^{T}\in{\mathbb{R}}^{n+q}$ and show the stability of the solution. The proof is rather technical so we devide it into five steps, in which we will show: 1) the exponential stability of $x(t)$; 2) some propeties of $y(t)$; 3) the exponential stability of $y(t)$ when $|x(0)|=0$; 4) the exponential stability of $y(t)$ when $|x(0)|>0$; 5) the exponential stability of $z(t)$. Some ideas and techniques in this proof are derived from our results on input-to-state stability (ISS) of SDEs [23, Theorem 3.1 and Remark 3.1], where $x(t)$ is treated as disturbance in the subsystem $y(t)$.

Step 1: Observe that (6a), (7a) and (7c) as well as $\ln\beta<\alpha\underline{\Delta t}$ from (8) are a specific case of conditions (i), (ii) and (iii) of [28, Theorem 3] in which $\lambda_{1}=\gamma_{1}=\cdots=\gamma_{\bar{m}}=0$. This implies

[TABLE]

where $\bar{\alpha}\in(0,\alpha-\bar{a})$ and $\bar{a}\in(0,\alpha)$ with $\ln\beta<\bar{a}\underline{\Delta t}<\alpha\underline{\Delta t}$. By condition (6a), subsystem $x(t)$, which is part (1a,1c) of the system (1), is $p$th moment exponentially stable (with Lypunov exponent no larger than $-\bar{\alpha}$).

Step 2: Let us consider the dynamics of subsystem $y(t)$, which is the other part (1b,1d) of system (1). By Lemma 1 and the Itô formula, one can derive that

[TABLE]

for all $t_{k}\leq\tilde{t}\leq t<t_{k+1}$ and $k\in{\mathbb{N}}$ while condition (7b) produces

[TABLE]

on $[t_{k},t_{k+1})$ for all $k\in{\mathbb{N}}$. This means that ${\mathbb{E}}\tilde{V}(y(t),t)$ is right-continuous on $[0,\infty)$ and could only jump at impulse times $\{t_{k+1}\}_{k\in{\mathbb{N}}}$. Notice condition (8) implies that $\tilde{\beta}_{2}e^{\tilde{\alpha}_{2}\,\overline{\Delta t}}<1$ and there is a pair of positives $\delta\in(0,1-\tilde{\beta}_{2})$ and $\bar{\delta}\in(0,\bar{\alpha}]$ sufficiently small for

[TABLE]

It is easy to observe from (11) that

[TABLE]

for such $t\in[t_{k},t_{k+1})$ and $k\in{\mathbb{N}}$ that

[TABLE]

Similarly, one can observe from (7d) that

[TABLE]

whenever

[TABLE]

Step 3: If $x(0)=0$ (equivalently, by (6a), $V(x(0),0)=0$), then (9) gives ${\mathbb{E}}V(x(t),t)=0$ for all $t\geq 0$. Using (7b), (10) and (11), one has

[TABLE]

This, by the Gronwall inequality ([32, Lemma 4.5.1, p129], [39, Theorem 8.1, p45]) implies

[TABLE]

and, particularly, ${\mathbb{E}}\tilde{V}(y(t_{1}^{-}),t_{1}^{-})\leq\tilde{V}(y(0),0)\,e^{\tilde{\alpha}_{2}\,t_{1}}$. Conditions (7d) and (12) produce

[TABLE]

One can repeat the derivations (15)-(17) over the interval between any two consecutive impulse times $[t_{k},t_{k+1}]$ and obtain

[TABLE]

for all $t\in[t_{k},t_{k+1})$ and $k\in{\mathbb{N}}$. This implies

[TABLE]

By condition (6b), subsystem $y(t)$ is $p$th moment exponentially stable (with Lyapunov exponent no larger than $-(\delta+\bar{\delta})$) when $x(0)=0$.

Step 4: We shall show the exponential stability of subsystem $y(t)$ when $|x(0)|>0$, namely, $V(x(0),0)\geq c_{1}|x(0)|>0$. Recall that both ${\mathbb{E}}V(x(t),t)$ and ${\mathbb{E}}\tilde{V}(y(t),t)$ are right-continuous on $[0,\infty)$, which could only jump at impulse times $\{t_{k+1}\}_{k\in{\mathbb{N}}}$. Define a function $\bar{v}:{\mathbb{R}}_{+}\to{\mathbb{R}}$ as

[TABLE]

with initial value

[TABLE]

Due to the properties of ${\mathbb{E}}V(x(t),t)$ and ${\mathbb{E}}\tilde{V}(y(t),t)$, $\bar{v}(t)$ is right-continuous on $[0,\infty)$ and could only jump at impulse times $\{t_{k+1}\}_{k\in{\mathbb{N}}}$. Given any $t\geq 0$, either $\bar{v}(t)\geq 0$ or $\bar{v}(t)<0$. So the interval $[0,\infty)$ is broken into a disjoint union of subsets $T_{+}\cup T_{-}$, where

[TABLE]

Obviously, from (9), one has

[TABLE]

Notice that $V(x(0),0)>0$. One has $\bar{v}(0)>0$ if $\tilde{V}(y(0),0)=0$; othewise, one can always choose a sufficiently small $\delta$ such that

[TABLE]

and, hence, $\bar{v}(0)>0$. Without loss of generality, assume that $\bar{v}(0)>0$. Therefore, $\bar{v}(t)>0$ on $[0,\epsilon)$ for some $\epsilon>0$, i.e., $[0,\epsilon)\subset T_{+}$. If $T_{+}=[0,\infty)$ (i.e., $T_{-}=\emptyset$), by (22), the proof is complete. Otherwise, let us consider the right-continuous process ${\mathbb{E}}\tilde{V}(y(t),t)$ on the subset $T_{-}$. Due to the right-continuity of $\bar{v}(t)$ on $[0,\infty)$, $\forall\,\bar{t}\in T_{-}$, there exists an ordered pair $\tau_{1}(\bar{t})<\tau_{2}(\bar{t})$ such that

[TABLE]

where $\tau_{1}(\bar{t})=\inf\{\underline{\tau}\leq\bar{t}:\bar{v}(\tau)<0,\forall\tau\in[\,\underline{\tau},\bar{t}\,]\}$ and $\tau_{2}(\bar{t})=\sup\{\bar{\tau}>\bar{t}:\bar{v}(\tau)<0,\forall\tau\in[\,\bar{t},\bar{\tau})\}$. For convenience, we also write $\tau_{1}=\tau_{1}(\bar{t})$ and $\tau_{2}=\tau_{2}(\bar{t})$ when there is no ambiguity. The right-continuous ${\mathbb{E}}\tilde{V}(y(t),t)$ on the interval $[\tau_{1},\tau_{2})$ falls into one of the three categories:

1. (C0)

There is no impulse time on $[\tau_{1},\tau_{2})$. 2. (C1)

There is exactly one impulse time on $[\tau_{1},\tau_{2})$. 3. (C2)

There are at least two impulse times on $[\tau_{1},\tau_{2})$.

Each of them is considered as follows.

(C0) There is $k\in{\mathbb{N}}$ such that $t_{k}<\tau_{1}<\tau_{2}\leq t_{k+1}$. Since $\bar{v}(t)$ is right-continuous and could only jump at impulse times $\{t_{k+1}\}_{k\in{\mathbb{N}}}$, that $t_{k}<\tau_{1}<\tau_{2}\leq t_{k+1}$ implies that $\bar{v}(t)$ is continuous on $t=\tau_{1}$ and, hence, by (23), $\bar{v}(\tau_{1})=0$. This means

[TABLE]

Using the Gronwall inequality, one can derive from (10), (11) and (13) that

[TABLE]

for all $t\in[\tau_{1},\tau_{2})$. Notice that $\tau_{2}-\tau_{1}<t_{k+1}-t_{k}\leq\overline{\Delta t}$. Substitution of (24) into (25) leads to

[TABLE]

(C1) There is exactly one impulse time $t_{k}$ on $[\tau_{1},\tau_{2})$, where $k$ is some positive integer since $[0,\epsilon)\subset T_{+}$. There are two cases: (C1a) $\tau_{1}<t_{k}$ and (C1b) $\tau_{1}=t_{k}$.

(C1a) There is $k\geq 1$ such that $t_{k-1}<\tau_{1}<t_{k}<\tau_{2}\leq t_{k+1}$. As above, this means that $\bar{v}(t)$ is continuous on $t=\tau_{1}$ and $\bar{v}(\tau_{1})=0$. So one has (24) at $t=\tau_{1}$ and (26) for all $t\in[\tau_{1},t_{k})$. But, from (14) and (12), one has

[TABLE]

By the Gronwall inequality, inequalities (10), (11), (13) and (27) produce

[TABLE]

for all $t\in[t_{k},\tau_{2})$. Therefore, when $\tau_{1}<t_{k}<\tau_{2}$, combination of (24), (26) and (28) imply that (26) holds for all $t\in[\tau_{1},\tau_{2})$.

(C1b) By the property of $\bar{v}(t)$ and definition of $\tau_{1}$, that $\tau_{1}=t_{k}$ implies $\bar{v}(t_{k}^{-})\geq 0$ and hence

[TABLE]

Combining (7d), (9) and (29) produces

[TABLE]

and hence

[TABLE]

for all $t\in[\tau_{1},\tau_{2})$. This combined with (26) yields

[TABLE]

when there is only one impulse time on the interval $[\tau_{1},\tau_{2})$, where $K$ is a positive constant

[TABLE]

(C2) There are at least two impulse times on $[\tau_{1},\tau_{2})$. For any two consecutive impulse times $t_{k}$ and $t_{k+1}$ on $[\tau_{1},\tau_{2})$, using the reasoning as above, one can derive that

[TABLE]

for all $t\in[t_{k},t_{k+1})$ and then

[TABLE]

Denote by $t_{k}<\cdots<t_{\bar{k}+1}<\cdots$ the impulse times on $[\tau_{1},\tau_{2})$, where $\bar{k}\geq k\geq 1$. Let us consider ${\mathbb{E}}\tilde{V}(y(t),t)$ on the interval $[t_{k},\tau_{2})$. By (33) and (34), one obtains

[TABLE]

and, therefore,

[TABLE]

Recall that $0<\bar{\delta}\leq\bar{\alpha}$ and $0\leq t_{k}-\tau_{1}\leq\overline{\Delta t}$. Again, there are two cases: $\tau_{1}<t_{k}$ and $\tau_{1}=t_{k}$. In the case where $\tau_{1}<t_{k}$, from (24), (27) and (36), one has

[TABLE]

for all $t\in[t_{k},\tau_{2})$ and then, by (26),

[TABLE]

In the other case where $\tau_{1}=t_{k}$, substituting (30) into (36) gives

[TABLE]

Combination of (38) and (39) yields

[TABLE]

for all $t\in[\tau_{1},\tau_{2})$ on which there are at least two impulse times, where $K$ is the positive constant given by (32). From inequlities (26), (31) and (40), one has

[TABLE]

Combining (22) and (41), one can conclude that

[TABLE]

By condition (6b), this means that subsystem $y(t)$, the other part of the system, is also $p$th moment exponentially stable (with Lyapunov exponent no larger than $-\bar{\delta}$) when $|x(0)|>0$.

Step 5: We have shown the $p$th moment exponential stability of $x(t)$ by (9) and that of $y(t)$ by (19) and (42) when $|x(0)|=0$ and $|x(0)|>0$, respectively.

Note that $z(t)=[x^{T}(t)\;\,y^{T}(t)]^{T}=[x^{T}(t)\;\,0^{T}]^{T}+[0^{T}\;\,y^{T}(t)]^{T}$ and, hence, $|z(t)|\leq|x(t)|+|y(t)|$ for all $t\geq 0$. It is easy to see that

[TABLE]

where $k_{p}=1$ when $0<p<1$ and $k_{p}=2^{p-1}$ when $p\geq 1$. Obviously, in the case where $|x(0)|=0$ and then ${\mathbb{E}}|x(t)|^{p}=0$ for all $t\geq 0$, inequalities (19) and (43) as well as $|z(0)|=|y(0)|$ give

[TABLE]

for all $t\geq 0$. In the general case where $|x(0)|>0$, inequalities (6), (9), (42) and (43) produce

[TABLE]

for all $t\geq 0$, where $K$ is the positive constant given by (32) and $\bar{K}_{p}$ is a positive constant

[TABLE]

So (44) and (45) mean that system (4) (i.e., system (1)) is $p$th moment exponentially stable (with Lyapunov exponent no larger than $-\bar{\delta}$). ∎

Remark 2**.**

Notice that, in Theorem 1, the continuous dynamics of subsystem $x(t)$ stabilizes the subsystem, though the discrete one could destabilize it, while the discrete dynamics of subsystem $y(t)$ stabilizes the subsystem, though the continuous one could destabilize it, which results in the exponential stability of the both subsystems and hence that of the whole system $z(t)=[x^{T}(t)\;y^{T}(t)]^{T}$. Similarly, one can obtain a stability criterion for the case where the impulses stabilize the physical subsystem $x(t)$ as the continuous dynamics could destabilize it (see [28, Theorem 2]) while the conditions on the subsystem $y(t)$ are kept the same as those in Theorem 1.

Furthermore, under Assumption 1, we have the following result for the almost sure stability.

Theorem 2**.**

If Assumption 1 holds, then the $p$th ($p>0$) moment exponential stability of SiDE (4) (i.e., system (1)) implies that it is also almost surely exponentially stable.

The proof is similar to that of [39, Theorem 4.2, p128] and, therefore, is omitted.

4 The cyber-physical systems of numerical methods

In this section, we address the problem (I) of fundamental importance. We compose a hybrid model in the form of our proposed SiDE (1) to represent the tight integration of the physical system (the SDE) and its cyber counterpart (the numerical method). This systematic representation can serve as a canonic form of CPS models.

Let us consider a physical system described by SDE

[TABLE]

with initial value $x(0)\in{\mathbb{R}}^{n}$, where $f:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ and $g:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n\times m}$ satisfy the global Lipschitz condition

[TABLE]

as well as $f(0)=0$ and $g(0)=0$ for study of the stability problem. Given a fixed parameter $\theta\in[0,1]$, the following numerical scheme for SDE (46) is called the stochastic theta method [17, 19, 40, 41]

[TABLE]

with initial value $X_{0}=x(0)$, where $\Delta t>0$ is the constant stepsize and $\sqrt{\Delta t}\,\xi(k+1)$ is the implementation of $\Delta B_{k}=B((k+1)\Delta t)-B(k\Delta t)$. The stochastic theta method for SDEs is one of the popular algorithms [32, 18] used in the so-called actor models [36], techniques driven by software modelling and simulation tools, to describe and compute the physical dynamics (46). Numerical method (48) is in essence a cyber model of the physical system (46), which is a translation of the physical system (46) into discretization, the symbolism in computers, and represents the physical dynamics in the cyber world.

When $\theta=0$, numerical scheme (48) gives the widely-used Euler-Maruyama method. The Euler-Maruyama method applied to SDE (46) computes approximations $X_{k}\approx x(t_{k})$ for all $k\in{\mathbb{N}}$, where $t_{k}=k\Delta t$, by setting $X_{0}=x(0)$ and forming

[TABLE]

Stochastic difference equations (49), also called discrete-time stochastic systems [26], have been intensively studied over the past a few decades of in the age of computers. In practice, it is natural to form and use some continuous-time extension of the discrete approximation $\{X_{k}\}$ such as [19, 41]

[TABLE]

where ${\mathds{1}}_{T}$ is the indicator function of set $T$. This is a simple step process of the equidistant Euler-Maruyama approximations so its sample paths are continuous on $(t_{k},t_{k+1})$ for each $k\in{\mathbb{N}}$ and right-continuous on $[0,\infty)$.

This paper consider the widely-used Euler-Maruyama method (49)-(50) the cyber system, which is virtually a representative of the physical system (46) in the cyber world. Other numerical schemes, or say, other translations can also be employed to represent the physical system in the cyber world in future work. This section is to unveil the inherent relationship between a physical system and its cyber couterpart. Consider the process $y(t)$ of difference between the exact solution $x(t)$ of the physical system (46) and the numerical solution $X(t)$ by its cyber counterpart (49)-(50)

[TABLE]

with initial value $y(0)=x(0)-X(0)=0$. Notice that $x(t)$ is a process of continuous paths and $X(t)$ a simple step process. This implies that $y(t)$ is right-continuous on $[0,\infty)$ and could only jump at $\{t_{k+1}\}_{k\in{\mathbb{N}}}$. According to the scheme (49)-(50), the jump of $y(t)$ at $t=t_{k+1}$ is

[TABLE]

since $X(t)=x(t)-y(t)=X_{k}$ for all $t\in[t_{k},t_{k+1})$ and $k\in{\mathbb{N}}$.

The integrated dynamics of the physical system (46) and the process (51) of difference is described by the following hybrid system in the form of SiDEs

[TABLE]

with initial data $x(0)\in{\mathbb{R}}^{n}$ and $y(0)=x(0)-X(0)=0$. Clearly, this is a special case of our proposed impulsive system (1) in which $h_{f}(\cdot,\cdot)\equiv 0$ and $h_{g}(\cdot,\cdot)\equiv 0$, that is, no impulse is imposed on the physical subsystem (53a). We construct the CPS (53) of the widely-used Euler-Maruyama method (49)-(50) for the SDE (46), which is a seamless, fully synergistic integration of the physical system (46) and its cyber counterpart (49)-(50). The CPS not only provides a holistic view of the physical system and its cyber counterpart but also reveals their intrinsic relationship that they are not two separate systems but the components of an integrated system.

Recall that the SDE describes our knowledge of the physical dynamics while the numerical method is the cyber representive, namely, the translation of our knowledge in the cyber world. As a result, the CPS (53) is the combination of our knowledge of the physical system and the cyber representative as well as the simulation sequence $\{\xi(k)\}_{k\in{\mathbb{N}}}$. Moreover, the CPS clearly shows that the numerical solution is driven by the SDE and the numerical method as well as the simulation sequence while the exact solution is, of course, conducted by the SDE itself only. Usually, to control the underlying physical system, our knowledge and the cyber representive of the system are both utilized in the synthesis of the CPS. This leads to the resulting CPS with $y(t)$ involved in the dynamics equation of $x(t)$ as well, which is considered in the study of stabilization problems [29].

Notice that the physical system (46) is subject to no impulse and its cyber-physical model (53) is a specific case of our proposed canonic form (1) of CPSs in which $q=n$, $f(x,t)=f(x)$, $g(x,t)=g(x)$, $\tilde{f}(x,y,t)=f(x)$, $\tilde{g}(x,y,t)=g(x)$, $h_{f}(x,k)=0$, $h_{g}(x,k)=0$, $\tilde{h}_{f}(x,y,k)=-f(x-y)\Delta t$, $\tilde{h}_{g}(x,y,k)=-g(x-y)\sqrt{\Delta t}$ and $t_{k}=k\Delta t$. Consequently, the infinitesimal generators (2) and (3) associated with (53a) and (53b) are of the specific forms

[TABLE]

respectively. It is easy to see that Assumption 1 holds since both $f$ and $g$ satisfy the global Lipschitz condition (47). Accordint to Lemma 1, there exists a unique (right-continuous) solution to SiDE (53) on $t\geq 0$ and the solution belongs to ${\mathcal{M}}^{2}([0,T];{\mathbb{R}}^{2n})$ for all $T\geq 0$. Moreover, the results of our established stability theory for the general class (1) of SiDEs, say, Theorem 1 and Theorem 2 apply to the CPS (53).

5 Exponential stability of the cyber-physical systems

The CPS (53) of the Euler-Maruyama method (49)-(50) for the SDE (46) consists of the physical subsystem and the cyber subsystem. The key questions (Q1) and (Q2) naturally arise. In this section, we address the problem (II) of fundamental importance and prove positive results to the key questions (Q1) and (Q2). These fundamental results and their applicaton to linear systems comprise a foundational theory of the CPSs of numerical methods for SDEs.

Let us begin with the test problem (Q1) of the CPS (53), to which Theorem 1 and Theorem 2 can be directly applied. Under some conditions (see [14, 31]), a seminal converse Lyapunov theorem [31, Theorem 5.12, p172] states that, if the SDE (46) is $p$th moment exponentially stable, there is a Lypunov function that proves the exponential stability of the dynamical system. One may postulate that the Lyapunov function for the physical subsystem (53a) could be used to construct a candidate Lyapunov function for the subsystem (53b,53c) due to their interrelation. The direct application of Theorem 1 to the CPS (53) shows that the CPS (53) and hence the cyber system (49)-(50) share the exponential stability with the physical system (46).

Theorem 3**.**

Let $V\in C^{2}({\mathbb{R}}^{n};{\mathbb{R}}_{+})$ be a candidate Lyapunov function for both subsystems (53a) and (53b,53c) and, moreover,

[TABLE]

for some positives $p,c_{1},c_{2}$. Assume that there are positives $\alpha$, $\tilde{\alpha}_{1}$, $\tilde{\alpha}_{2}$, $\tilde{\beta}_{1}$, $\tilde{\beta}_{2}$ such that

[TABLE]

for all $k\in{\mathbb{N}}$. If the stepsize

[TABLE]

then the CPS (53) is $p$th moment exponentially stable and is also almost surely exponentially stable. Moreover, the cyber system (49)-(50) shares the $p$th moment exponential stability with its underlying physical system (46) and, hence, it is also almost surely exponentially stable.

Proof.

From Theorem 1 and Theorem 2, it follows that CPS (53), which is a specific case of system (1) (namely, (4)), is $p$th moment exponentially stable and is also almost surely exponentially stable.

Notice that the state of cyber system (49)-(50) is the difference $X(t)=x(t)-y(t)$ of the subsystems (53a) and (53b,53c). Therefore,

[TABLE]

where $k_{p}=1$ if $0<p<1$, and $k_{p}=2^{p-1}$ if $p\geq 1$. This immediately implies that the cyber system (49)-(50) is $p$th moment exponentially stable and is also almost surely exponentially stable. ∎

This means that, if the underlying physical system (46) is $p$th moment exponentially stable, the CPS (53) and, hence, the numerical method (49)-(50) reproduce the $p$th moment exponential stability of the physical dynamics when the conditions in Theorem 3 hold. The ability of the cyber system (the numerical method) to reproduce the mean-square exponential stability of its underlying physical system (the SDE) has been studied in [17, 20], Here, we investigate on the ability of the cyber system (49-50) to reproduce the mean-square exponential stability of the physical system (46) in our proposed framework of CPS (53). A result on mean-square exponential stability is then derived from Theorem 3 as follows, in which the Lyapunov function for mean-square exponential stability of the underlying physical system (46) also proves the mean-square exponential stability of its cyber counterpart (49)-(50) as well as that of the CPS (53).

Theorem 4**.**

Let the candidate Lyapunov function $V\in C^{2}({\mathbb{R}}^{n};{\mathbb{R}}_{+})$ for physical system (46) be a quadratic function

[TABLE]

for some positive definite matrix $P\in{\mathbb{R}}^{n\times n}$. Assume there are positives $\bar{\alpha}$ and $\overline{\Delta t}$ with $\bar{\alpha}\,\overline{\Delta t}<1$ such that

[TABLE]

for all $x\in{\mathbb{R}}^{n}$. Then the CPS (53) with $\Delta t\in(0,\overline{\Delta t}\,]$ is mean-square exponentially stable and is also almost surely exponentially stable. Moreover, the cyber system (49)-(50) with $\Delta t\in(0,\overline{\Delta t}\,]$ shares the mean-square exponential stability with its underlying physical system (46) and, hence, it is also almost surely exponentially stable.

Proof.

It will follow the conclusion from Theorem 3 if one shows that conditions (56)-(58) of Theorem 3 are satisfied with $p=2$ for the CPS (53). Since the quadratic function (60) gives $\lambda_{m}(P)|x|^{2}\leq V(x)\leq\lambda_{M}(P)|x|^{2}$, condition (56) holds for positive constants $p=2$, $c_{1}=\lambda_{m}(P)$, $c_{2}=\lambda_{M}(P)$. Moreover, (57a) and (61) are equivalent. Obviously, (61) implies that (57a) holds for some positive constant $\alpha\geq\bar{\alpha}$. But (57a) implies that there is a sufficiently small positive number $\overline{\Delta t}$ such that (61) holds for some positive constant $\bar{\alpha}<\alpha\wedge(1/\,\overline{\Delta t})$. According to [39, Theorem 4.4, p130] and [39, Theorem 4.2, p128], system (46) is mean-square exponentially stable and also almost surely exponentially stable.

By the Itô lemma, [25, Lemmas 1 and 2] and global Lipschitz condition (47),

[TABLE]

where $\tilde{\alpha}_{1}=\frac{(1+\tilde{\alpha}_{2})\,\lambda_{M}(P)L^{2}}{\tilde{\alpha}_{2}\,\lambda_{m}(P)}$ and $\tilde{\alpha}_{2}$ given as (69) below are both positive numbers. So condition (57b) of Theorem 3 is satisfied. Note that, for any $\Delta t\in(0,\overline{\Delta t}\,]$, inequality (61) implies

[TABLE]

and (53c) gives

[TABLE]

Using inequality (63) and [25, Lemma 1], one obtains

[TABLE]

for all $(x,y)\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}$, where $c$ is a positive constant with $c<\bar{\alpha}\Delta t$ and, hence, $1+c-\bar{\alpha}\Delta t<1$.

By [25, Lemmas 1-2] and global Lipschitz condition (47), one also has

[TABLE]

where $b$ is a positive constant sufficiently small for

[TABLE]

Substitution of (65) and (66) into (64) yields

[TABLE]

where

[TABLE]

and $\tilde{\beta}_{2}$ given as (67) above are both positive constants. This is the condition (57c) of Theorem 3.

Let $\tilde{\alpha}_{2}$ be a positive number such that

[TABLE]

For instance, let $\tilde{\alpha}_{2}=-\ln\tilde{\beta}_{2}/(2\,\overline{\Delta t})$ and, hence, $\Delta t\leq\overline{\Delta t}=\frac{-\ln\tilde{\beta}_{2}}{2\,\tilde{\alpha}_{2}}<\frac{-\ln\tilde{\beta}_{2}}{\tilde{\alpha}_{2}}$. So the condition (58) of Theorem 3 is also satisfied. According to Theorem 3, it follows the assertions. ∎

In the literature [20, 40], to ensure that the cyber system shares the exponential stability with its underlying physical system, the stepsize $\Delta t$ is explicitly and severly limited by both the growth and the rate constants of the physical system. Although the both are related, it is only the rate constant that plays a key role in the definitions of exponential stability. It is reasonable and possible to lessen the dependence of the stepsize $\Delta t$ on the growth constant, which itself could be very conservative due to condition (56). In Theorem 4, we manage to remove the explicit dependence of the stepsize $\Delta t$ on the growth constant $\lambda_{M}(P)/\lambda_{m}(P)$. Instead, we show that the growth constant $\lambda_{M}(P)/\lambda_{m}(P)$ has an influence on the stepsize $\Delta t$ through the rate-like constant $\tilde{\beta}_{2}$ given by (67). This could reduce much the restriction imposed by the growth constant. As will be shown in Section 6, it improves the upper bound $\overline{\Delta t}$ of stepsizes and eases its computation for linear systems.

Recall that Theorem 3 is the direct application of Theorem 1 to the CPS (53) of the numerical method (49)-(50) for the SDE (46), from which Theorem 4 is derived for the mean-square exponential stability. Let us proceed to apply Theorem 4 and study the converse question (Q2) whether one can infer that the CPS (53) and, hence, the physical system (46) are mean-square exponentially stable if the cyber system (49)-(50) is mean-square exponentially stable for small stepsizes $\Delta t>0$. Similarly, the converse Lyapunov theorem [31, 53] gives that, if the cyber system (49)-(50) is mean-square exponentially stable, there is a Lyapunov function that proves the exponential stability of the system. One may make use of this Lyapunov function for the cyber system and study the stability of the physical system and that of the whole CPS. Applying Theorem 4, we find that the mean-square exponential stability of the CPS (53) and, hence, that of the physical system (46) can be inferred from the mean-square exponential stability of the cyber system (49)-(50) if the Lyapunov function is of the quadratic form.

Theorem 5**.**

Assume that there is a candidate Lyapunov function $V\in C^{2}({\mathbb{R}}^{n};{\mathbb{R}}_{+})$ of the quadratic form (60) for the cyber system (49)-(50) with stepsize $\Delta t=\overline{\Delta t}>0$ such that

[TABLE]

for some positive constant $\bar{c}<1$ and all $X_{k}\in{\mathbb{R}}^{n}$. Then CPS (53) with $\Delta t\in(0,\overline{\Delta t}]$ is mean-square exponentially stable and also almost surely exponentially stable, which implies that physical system (46) is mean-square exponentially stable and also almost surely exponentially stable.

Proof.

By Lyapunov stability theory [5, 31], conditions (60) and (70) as well as (56) derived from (60) immediately imply that the cyber system (49)-(50) with $\Delta t=\overline{\Delta t}$ is mean-square exponentially stable. Let function $V(x)$ given by (60) also be the candidate Lyapunov function for the physical system (46). But condition (70)

[TABLE]

implies that X^{T}_{k}Pf(X_{k})+f^{T}(X_{k})PX_{k}+trace\big{[}g^{T}(X_{k})Pg(X_{k})\big{]}+\overline{\Delta t}f^{T}(X_{k})Pf(X_{k})\leq-\bar{\alpha}V(X_{k}) for all $X_{k}\in{\mathbb{R}}^{n}$, where $\bar{\alpha}\,{\overline{\Delta t}}=1-\bar{c}$. This means

[TABLE]

for all $x\in{\mathbb{R}}^{n}$, which is exactly the condition (61) of Theorem 4. According to Theorem 4, $z(t)=[x^{T}(t)\;\,y^{T}(t)]^{T}$ of CPS (53) with $\Delta t\in(0,\overline{\Delta t}]$ is mean-square exponentially stable and is also almost surely exponentially stable. It immediately follows that, due to $|x(t)|^{2}\leq|z(t)|^{2}$, the physical system (46) is mean-square exponentially stable and is also almost surely exponentially stable. Alternatively, conditions (60) and ${\mathscr{L}}V(x)\leq{\mathscr{L}}V(x)+\overline{\Delta t}\,V(f(x))\leq-\bar{\alpha}V(x)$ for all $x\in{\mathbb{R}}^{n}$ as well as (56) derived from (60) imply that, by [39, Theorems 4.2-4.4, pp128-130], the physical system (46) is mean-square exponentially stable and is also almost surely exponentially stable. ∎

Our positive results to the key questions (Q1) and (Q2) expose the equivalence and intrinsic relationship $\bar{\alpha}\,{\overline{\Delta t}}=1-\bar{c}$ between (61) and (70), which are the stability conditions for the physical system (46) and its cybercounterpart (49)-(50), respectively. For this, we employ the same Lyapuov function $V(x)=\tilde{V}(x)=x^{T}Px$ for both the subsystems in Theorems 3-5, see also Section 6. Our proposed theory can further developed by using some techniques of Lyapunov functions and functionals to exploit the composition and structure (75), see Remark 4 as well as [12, 28, 29, 37, 42]. It is also worth noting that the initial condition $y(0)=x(0)-X(0)=0$ is not a requirement in our established stability theory and its application in Section 6. But this condition could make a difference in the study of convergence and some control problems, see Appendix B and [29].

6 Application to linear systems

Let us consider a linear stochastic system

[TABLE]

with initial value $x(0)\in{\mathbb{R}}^{n}$, where $F\in{\mathbb{R}}^{n\times n}$ and $G_{j}\in{\mathbb{R}}^{n\times n}$, $j=1,2,\cdots,m$, are constant matrices. Obviously, the linear system (71) satisfies the global Lipschitz condition and has a unique solution $x(t)$ on $[0,\infty)$. It is well known that the linear stochastic system (71) is mean-square exponentially stable if and only if there exists a positive definite matrix $P\in{\mathbb{R}}^{n\times n}$ such that [5, 11]

[TABLE]

This is the classic Lyapunov-Itô inequality [5, 11], the linear matrix inequality (LMI) equivalent [13] of the classical Lyapunov-Itô equation [1, 38]. By [39, Theorem 4.2, p128], the mean-square exponential stability of (71) implies that it is also almost surely exponentially stable.

The Euler-Maruyama method (49)-(50) for the linear system (71) computes approximations

[TABLE]

where $X_{0}=x(0)$, $\Delta t>0$ is the constant stepsize and $\sqrt{\Delta t}\,\xi_{j}(k+1)$ is the implementation of $\Delta B_{j,k}=B_{j}((k+1)\Delta t)-B_{j}(k\Delta t)$. The cyber system (73) is mean-square exponentially stable if and only if there exists a positive definite matrix $P\in{\mathbb{R}}^{n\times n}$ such that [5]

[TABLE]

Let $y(t)$ be the difference between $x(t)$ and $X(t)$ as (51) above. The CPS of the Euler-Maruyama method (73) for the linear SDE (71) is a specific case of CPS (53)

[TABLE]

with initial data $x(0)\in{\mathbb{R}}^{n}$ and $y(0)=x(0)-X(0)=0$, where $t_{k}=k\Delta t$ for all $k\in{\mathbb{N}}$. The CPS (75) is an integration of the physical system (71) and the cyber system (73), which is of our proposed class (1) and satisfies the global Lipschitz conditions Assumption 1. Our established theory immediately provides positive results to the key questions (Q1) and (Q2) of linear CPS (75), which presents the upper bound $\overline{\Delta t}$ of stepsizes for exponential stability.

Theorem 6**.**

The following are equivalent.

• ($\mathcal{A}$)

There exists a positive definite matrix $P\in{\mathbb{R}}^{n\times n}$ such that the CPS Lyapunov inequality holds for some positive number $\overline{\Delta t}$, namely,

[TABLE]

• ($\mathcal{B}$)

The physical system (71) is mean-square exponentially stable.

• ($\mathcal{C}$)

The cyber system (73) with $\Delta t\in(0,\overline{\Delta t}\,]$ is mean-square exponentially stable.

• ($\mathcal{D}$)

The CPS (75) with $\Delta t\in(0,\overline{\Delta t}\,]$ is mean-square exponentially stable.

That is, $(\mathcal{A})\Leftrightarrow(\mathcal{B})\Leftrightarrow(\mathcal{C})\Leftrightarrow(\mathcal{D})$.

Proof.

$(\mathcal{A})\Leftrightarrow(\mathcal{B})$. We only need to show that the classic Lyapunov inequality (72) and the CPS Lyapunov inequality (76) are equivalent. But this is implied by the equivalence of the inequalities (57a) and (61) which has been shown in the proof of Theorem 4. Alternatively, it can be easily proved as follows. Clearly, (76) implies (72). But inequality (72) implies that there is a sufficiently small positive number $\overline{\Delta t}$ such that (76) holds. Therefore, the LMI (72) $\Leftrightarrow$ the LMI (76).

$(\mathcal{A})\Rightarrow(\mathcal{C})\;\&\;(\mathcal{D})$. Let us consider the quadratic Lyapunov function $V(x)=x^{T}Px$ for the linear system (71). The LMI (76) implies that there is a positve number $\bar{\alpha}<1/\,\overline{\Delta t}$ sufficiently small for

[TABLE]

and the condition (61) holds. It follows from Theorem 4 that the CPS (75) and, hence, the cyber system (73) with $\Delta t\in(0,\overline{\Delta t}\,]$ are mean-square exponentially stable.

$(\mathcal{D})\Rightarrow(\mathcal{B})\;\&\;(\mathcal{C})$. The CPS (75) is a specific case of system (4), where $z(t)=[x^{T}(t)\;\,y^{T}(t)]^{T}$ in the compact form. Notice that $|x(t)|^{2}\leq|z(t)|^{2}$ and $|X(t)|^{2}\leq 2(|x(t)|^{2}+|y(t)|^{2})\leq 4|z(t)|^{2}$ for all $t\geq 0$. If $z(t)$ of the CPS (75) is mean-square exponentially stable, then both $x(t)$ of its physical subsystem (71) and $X(t)$ of its cyber subsystem (73) are mean-square exponentially stable.

$(\mathcal{C})\Rightarrow(\mathcal{B})\;\&\;(\mathcal{D})$. Let $\Delta t=\overline{\Delta t}$. Since the cyber system (73) is mean-square exponentially stable, there is a positive definite matrix $P\in{\mathbb{R}}^{n\times n}$ such that the Lyapunov inequality (74) holds for $\Delta t=\overline{\Delta t}>0$. This implies that there is a positive number $\bar{c}\in(0,1)$ sufficiently close to $1$ for

[TABLE]

Let the quadratic function $V(x)=x^{T}Px$ be the candidate Lyapunov function for the cyber system (73) with $\Delta t=\overline{\Delta t}$. Observe that, for the linear system, (78) is exactly the condition (70) of Theorem 5. It follows from Theorem 5 that the CPS (53) with $\Delta t\in(0,\overline{\Delta t}]$ and, hence, the physical system (46) are mean-square exponentially stable. The proof is complete. ∎

Note that the mean-square exponential stability of the physical system (71), the cyber system (73) and the CPS (75) imply that they are also almost surely exponentially stable, respectively. It is easy to obtain the upper bound $\overline{\Delta t}$ of stepsizes for the ability of the cyber system to reproduce the exponential stability of the underlying linear physical system by solving the CPS Lyapunov inequality (76), which can also be called the numerical Lyapunov inequality in the study of numerical methods for differential equations.

Remark 3**.**

Consider a scalar SDE, which is the linear SDE (71) with $n=m=1$,

[TABLE]

where $\lambda$ and $\mu$ are both real constants. The CPS Lyapunov inequality (76) immediately gives

[TABLE]

According to Theorem 6, this is the necessary and sufficent condition for mean-square exponential stability of the linear scalar physical system (79), its cyber system of the Euler-Maruyama method, and its CPS of the Euler-Maruyama method, that is, (75) with $n=m=1$, $F=\lambda$ and $G_{1}=\mu$. It is observed that (80) is exactly the inequality (4.3) in [17] with $\theta=0$ for the Euler-Maruyama method. Obviously, (80) is the very scalar case of our CPS Lyapunov inequality (76) that applies to general multi-dimensional linear systems.

Recently, based on the reformulation of some well-known results, [6] developed an approach to mean-square stability analysis of numerical methods (including the widely-used Euler-Maruyama scheme) for multi-dimensional linear SDEs (viz. system (71) with $n\geq 2$), which was applied in [7] to study the mean-square numerical stability for a linear SDE of non-normal drift and skew-symmetic diffusion structures [21]. Specifically, on one hand, some well-known result ([1, Theorem 8.5.5, p142], [31, Remark 6.4, p183]) expressed in the vectorization of matrices and Kronecker product gives [6, Lemma 3.3]

[TABLE]

if and only if linear SDE (71) is mean-square exponentially stable, where $Re_{M}({\mathcal{S}})$ is the real part of the eigenvalue $\lambda_{M}({\mathcal{S}})$ of $n^{2}\times n^{2}$ matrix

[TABLE]

On the other hand, a stability result for discrete-time stochastic systems (see, e.g., [31, p197]) is applied to study mean-square stability of some numerical schemes for the SDE (71). For example, by [6, Lemma 3.4, Theorem 3.7], the Euler-Maruyama method (73) is mean-square exponentially stable if and only if

[TABLE]

where $\rho({\mathcal{S}}_{0}(\Delta t))$ is the spectral radius of $n^{2}\times n^{2}$ matrix

[TABLE]

with $\bar{A}(\Delta t)=I_{n}+{\Delta t}F$ and $\bar{B}_{j}(\Delta t)=\sqrt{\Delta t}\,G_{j}$ for $j=1,\cdots,m$.

In [6], ${\mathcal{S}}$ is called the mean-square stability matrix of the SDE (71) and ${\mathcal{S}}_{0}(\Delta t)$ that of the Euler-Maruyama method (73). Notice that ${\mathcal{S}}_{0}({\Delta t})$ is a function of stepsize ${\Delta t}$ while, obviously, ${\mathcal{S}}$ is not. The results in [6, 7] provided the explicit structure of stability matrices ${\mathcal{S}}$ and ${\mathcal{S}}_{0}({\Delta t})$, and showed the comparative stability regions [7, Fig.2] for the SDE and the numerical method with a few numerical examples of non-normal SDEs [21]. However, the relationship between the stability conditions (81) and (82) (for the SDE and the Euler-Maruyama method, respectively) has not been figured out. Here we prove their equivalence and reformulate the stability conditions (81) and (82) in the form of LMIs, which is relegated to Appendix A. So it is easy to handle the problems with some computing techniques and toolboxes [5, 11, 13].

Note that it is easy to obtain the upper bound $\overline{\Delta t}$ of stepsizes for the test problem (Q1) by solving the $n\times n$-dimensional LMI (76) of our proposed method, which is a linear function with respect to $\Delta t$. Clearly, LMI (76) holds for all $\Delta t\in(0,\overline{\Delta t}]$ provided it is satisfied for some $\overline{\Delta t}>0$. But, to calculate the upper bound $\widetilde{\Delta t}$ by the approach of mean-square stability matrix [6], one has to deal with the spectral radius (82) of $n^{2}\times n^{2}$ matrix ${\mathcal{S}}_{0}({\Delta t})$ that involves a polynomial of the stepsize ${\Delta t}$ whose order is some exponential function of $n$. Alternatively, one can solve the following LMI with respect to positive definite matrix $\bar{P}\in{\mathbb{R}}^{n^{2}\times n^{2}}$, which we show is equivalent to inequality (82) in Appendix A,

[TABLE]

where $\bar{F}=F\otimes F$. This involves a cubic function of the prescribed parameter $\Delta t>0$ for all $n$. Note that, for a multi-dimensional SDE ($n\geq 2$), the spectral radius (82) of $n^{2}\times n^{2}$ matrix ${\mathcal{S}}_{0}({\Delta t})$ involves a polynomial of $\Delta t$ of up to (very) high order. For example, in the case $n=2$ of non-normal SDE [7, Eq.(9)] (see also [21]), the characteristic equation of mean-square stability matrix ${\mathcal{S}}_{0}(\Delta t)$ of $4\times 4$ dimensions for the Euler-Maruyama scheme ([7, Eq.(15)] with $\theta=0$) involves a polynomial of $\Delta t$ of up to order $8$. The conditions of this approach are quite cumbersome [31]. It is easy to tackle the equivalent LMI (83), which is a cubic function of the prescribed $\Delta t>0$ for all $n$, with some toolbox [5, 13]. However, one should be aware that, unlike the linear inequality (76), that LMI (83), which is not only of $n^{2}\times n^{2}$ dimensions but also involves a cubic function of the prescribed $\Delta t>0$, is satisfied for some $\overline{\Delta t}>0$ does not necessarily mean that it holds for all $\Delta t\in(0,\overline{\Delta t}]$ and, thereby, the results such as the upper bound $\overline{\Delta t}$ could be restrictive due to the highly nonlinearity of $\Delta t$ involved in the computation, see Appendix A.

We can further show that our proposed method (76) gives better bound $\overline{\Delta t}$ of stepsizes than (82) from [6] or its LMI equivalent (83). This is: if (82) and (83) hold for all $\Delta t\in(0,\widetilde{\Delta t}]$ for some $\widetilde{\Delta t}>0$, then the CPS Lyapunov inequality (76) holds for some $\overline{\Delta t}\geq\widetilde{\Delta t}$, namely, either $\overline{\Delta t}=\widetilde{\Delta t}$ or $\overline{\Delta t}>\widetilde{\Delta t}$. In short, we shall show either $\widehat{\Delta t}=\widetilde{\Delta t}$ or $\widehat{\Delta t}>\widetilde{\Delta t}$, where

[TABLE]

It is observed that, due to the continuity of (83) with respect to $\Delta t$, the strict inequality (83) does not hold at $\Delta t=\widetilde{\Delta t}$ and, similarly, the strict inequality (76) does not hold at $\overline{\Delta t}=\widehat{\Delta t}$. To show either $\widehat{\Delta t}=\widetilde{\Delta t}$ or $\widehat{\Delta t}>\widetilde{\Delta t}$, let us consider a linear SDE with $\Delta t\in(0,\overline{\Delta t}\,]$ for some $\overline{\Delta t}>0$

[TABLE]

where $B_{m+1}(t)$ is a scalar Brownian motion. Notice that (85) is exctly (71) if $\Delta t=0$ and, according to [6, Lemma 3.3], the SDE (85) is mean-square exponentially stable if and only if (87) holds. In fact, by the well-known results [1, 5, 6, 13, 31], the following are equivalent:

• (a)

The CPS Lyapunov LMI (76) holds.

• (b)

There is a positive definite matrix $\bar{P}\in{\mathbb{R}}^{n^{2}\times n^{2}}$ such that

[TABLE]

for each ${\Delta t}\in(0,\overline{\Delta t}\,]$.

• (c)

The following inequality holds for each ${\Delta t}\in(0,\overline{\Delta t}\,]$

[TABLE]

• (d)

The SDE (85) is mean-square exponentially stable.

That is, (a) $\Leftrightarrow$ (b) $\Leftrightarrow$ (c) $\Leftrightarrow$ (d) and they all give the same supremum $\widehat{\Delta t}$ defined as (84). Observe that (83) implies (86) but not vice versa. This means that (76), (86) and (87) hold for all $\overline{\Delta t}\in(0,\widetilde{\Delta t})$, namely, $\widehat{\Delta t}\geq\widetilde{\Delta t}$. Notice that (83), (86) and (87) are all continuous with respect to $\Delta t$. Due to the continuity of (83) at $\Delta t=\widetilde{\Delta t}$,

[TABLE]

Unless matrix ${\mathcal{S}}+\widetilde{\Delta t}\bar{F}$ is singular, the LMIs (86) and, equivalently, (76) hold at $\overline{\Delta t}=\widetilde{\Delta t}$ and, due to their continuity at $\Delta t=\widetilde{\Delta t}$, the LMIs (86) and (76) hold for some $\overline{\Delta t}>\widetilde{\Delta t}$, which gives $\widehat{\Delta t}>\widetilde{\Delta t}$. So we have $\widehat{\Delta t}=\widetilde{\Delta t}$ if matrix ${\mathcal{S}}+\widetilde{\Delta t}\bar{F}$ is singular; otherwise, $\widehat{\Delta t}>\widetilde{\Delta t}$. The latter could often be the case. This clearly shows that our proposed CPS Lyapunov LMI (76) gives better bound $\overline{\Delta t}$ of stepsizes than (82) from [6], or, its LMI equivalent (83).

Remark 4**.**

The upper bound $\overline{\Delta t}$ of stepsizes in the CPS Lyapunov inequality (76) is obtained by the same Lyapunov function $V(x)=\tilde{V}(x)=x^{T}Px$ for both the subsystems, which is special application of Theorem 3. The upper bound $\overline{\Delta t}$ could be improved by using some other Lypunov functions or functionals such as different Lyapunov functions for the subsystems [29] and some Lyapunov function or funcational of quadratic form with block matrices for the compact form of the CPS [12, 28, 37, 42] to exploit the structure of the composition (75) of the subsystems [37].

7 On systems numerics

In this paper, we formulated a new and general class (1) of SiDEs that can be used to represent a seamless integration of the physical system (the SDE) and its cyber counterpart (the numerical method), which is a novel framework for numerical study of dynamical systems. Our proposed CPS of the widely-used Euler-Maruyama method for SDEs not only provides a holistic view of the physical system (the SDE) and its cyber counterpart (the numerical method) but also reveals their intrinsic relationship: they are not two separate systems but the subsystems of the CPS. By our CPS approach, we proved positive results to the key questions (Q1) and (Q2) using the Lyapunov stability theory we established for our general class of SiDEs. These fundamental results and their application comprise a foundational theory of the CPSs of numerical methods for SDEs.

In the classical numerical analysis of initial-value problems, the convergence and the stability of a numerical method are two main concerns. The proposed CPS also provides a novel approach to convergence analysis of the numerical method for SDEs. Using the dyanimcs of the discretization error $y(t)$, we show the classical finite-time convergence result

[TABLE]

for the Euler-Maruyama method by our CPS approach, which is different from those in the literature [19, 32, 39, 51]. The novel proof is relegated to Appendix B.

Our proposed CPS and theory have initiated the study of systems numerics, where there are many open and interesting problems. For example, it is among the future work to extend our theory to many other (explicit or implicit) numerical methods [17, 19, 30, 32, 41, 46, 51] as well as to dyanmical systems such as SDEs with time delay and/or switching [23, 28].

Acknowledgement

The author is supported by National Natural Science Foundation of China (No.61877012). The author gratefully acknowledges Prof. X. Mao (the author’s PhD supervisor at University of Strathclyde, UK) for his very helpful comments and suggestions, which improve the quality of this work. Some of this work was done during the author’s visit to The Centre for Stochastic & Scientific Computations at Harbin Institute of Technology and the author would like to thank Prof. M. Song, Prof. M. Liu and their group for the helpful discusssions. The author is also thankful to Prof. K. C. Sou for the communications in the past a few years.

Appendix A. The equivalence of the stability conditions (81) and (82)

Proof.

It is observed that ${\mathcal{S}}_{0}(0)=I_{n^{2}}$, $\dot{{\mathcal{S}}}_{0}(0)={\mathcal{S}}$, $\ddot{{\mathcal{S}}}_{0}(\Delta t)=2\,(F\otimes F)=:2\bar{F}$, where $\dot{{\mathcal{S}}}_{0}(\Delta t)$ and $\ddot{{\mathcal{S}}}_{0}(\Delta t)$ are the first and the second derivatives of ${\mathcal{S}}_{0}$ with respect to $\Delta t$, respectively. For $\Delta t>0$, Taylor expansion produces

[TABLE]

(81) $\Rightarrow$ (82). Stability condition (81) for the linear SDE equivalently means that there is a positive definite matrix $\bar{P}\in{\mathbb{R}}^{n^{2}\times n^{2}}$ such that [5, 31]

[TABLE]

The Taylor series (88) gives

[TABLE]

Owing to (89), there is a pair of (sufficiently small) positive numbers $\overline{\Delta t}$ and $\bar{a}=\bar{a}(\overline{\Delta t})$ such that

[TABLE]

This implies that stability condition (82) is satisfied for all ${\Delta t}\in(0,\overline{\Delta t}\,]$.

(82) $\Rightarrow$ (81). Notice that (90) can be rewritten as

[TABLE]

Suppose that condition (82) holds for all ${\Delta t}\in(0,\overline{\Delta t}]$, where $\overline{\Delta t}$ is some positive number. Then, equivalently, there is a positive definite matrix $\bar{P}=\bar{P}({\Delta t})\in{\mathbb{R}}^{n^{2}\times n^{2}}$ such that [5, 31]

[TABLE]

Substitution of (91) into (92) produces (83) for all ${\Delta t}\in(0,\overline{\Delta t}\,]$, or, by Schur complement,

[TABLE]

So $(\ref{MS-SDE-EM})\,\Leftrightarrow\,(\ref{LMI-2012})\,\Leftrightarrow\,(\ref{EM-bP-3})\,\Leftrightarrow\,(\ref{EM-bP-3-2})$. But (83) implies the LMI (86). This equivalently means that matrix ${\mathcal{S}}+{\Delta t}\bar{F}$ is Hurwitz, namely, inequality (87) holds for each ${\Delta t}\in(0,\overline{\Delta t}\,]$.

Recall that $\bar{F}=F\otimes F$. Letting ${\Delta t}\to 0$ in (86) and thus (87) gives stability condition (81) for the SDE (71). The proof is complete. ∎

We remark that, by approach of mean-square stability matrices [6], the upper bound $\overline{\Delta t}$ can be calculated by solving either the spectral radius problem (82) or the LMI equivalent (83). The former involves a polynomial of the stepsize $\Delta t$ whose order is some exponential function of $n$ while the latter remains as a cubic function of $\Delta t$ for all $n$. The highly nonlinearity would introduce not only computational complexity but also conservativeness to the results. We have reformulted the highly nonlinear problem (82) into the LMI (83). This has significantly simplified the approach of mean-square stability matrices ${\mathcal{S}}$ and ${\mathcal{S}}_{0}(\Delta t)$. But it is worth noting that, for a linear $n$-dimensional SDE, our proposed numerical Lyapunov LMI (76) of $n\times n$ dimensions is always a linear inequality of the stepsize $\overline{\Delta t}$ while the LMI problem (83) involves not only a cubic function of $\Delta t$ but also matrices of $n^{2}\times n^{2}$ dimensions.

Appendix B. A novel proof of the convergence of the Euler-Maruyama method

Proof.

For the convergence problem of the numerical method, the implimentation $\sqrt{\Delta t}\,\xi(k+1)$ should be replaced by the increment $\Delta B_{k}=B((k+1)\Delta t)-B(k\Delta t)$ itself in SiDE (53), that is,

[TABLE]

with $x(0)\in{\mathbb{R}}^{n}$ and $y(0)=0$, where $t_{k}=k\Delta t$ for all $k\in{\mathbb{N}}$. According to the existing results ([39, 32] as well as Lemma 1), SiDE (94) has a unique (right-continuous) solution $z(t)=[x^{T}(t)\;y^{T}(t)]^{T}$, which belongs to ${\mathcal{M}}^{2}([0,T];{\mathbb{R}}^{n+q})$ for all $T\geq 0$. In particular, [39, Lemma 3.2, p51] gives

[TABLE]

On the interval $[t_{k},t_{k+1}]$ for every $k\in{\mathbb{N}}$,

[TABLE]

and, due to $y(0)=0$,

[TABLE]

where $[t]:=\sup\{t_{j}:t_{j}\leq t,j\in{\mathbb{N}}\}$ for all $t\geq 0$. By Cauchy-Schwaz inequality, this produces

[TABLE]

By the Itô isometry and the global Lipschitz condition (47),

[TABLE]

Since (94a) and (94b) give $x(t)-x([t])=y(t)-y([t])$ for all $t\geq 0$, this implies

[TABLE]

For any $T\geq 0$, using the Itô lemma, (47), (55) and (97), one obtains

[TABLE]

where constant $K_{T}=2L^{2}([T]+1)\vee 1$. This implies

[TABLE]

where $C_{T}$ is given by (95) above. In view of the Gronwall inequality ([32, Lemma 4.5.1, p129], [39, Theorem 8.1, p45]), this yields

[TABLE]

which completes the proof. ∎