# Modeling and Analysis of Wireless Cyberphysical Systems Using Stochastic Methods

Authors: Jing Liu; Yixu Yao
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

Wireless mesh network control systems (WMNCSs) are typical Cyberphysical Systems (CPSs) widely used in industries that need to meet stringent performance requirements. These WMNCSs are characteristic of stochasticity at different levels as system behavior, network performance, and wireless signal propagation, which grievously increases the difficulties of system modeling and analysis. In this paper, we propose a three-layered modeling framework to capture the stochastic properties of WMNCSs at different levels with stochastic hybrid system (SHS), stochastic network calculus (SNC), and physical layer models. We also bridge the gaps between these methods with an upper bound approach. All the efforts give us a methodology of modeling and analysis WMNCSs with stochastic methods, so we can know how the factors as channel conditions, network topology, etc. affect the stability and performance of the system. To the best of our knowledge, it is the first work that provides such a unified and flexible framework to model and analyze WMNCSs with stochastic methods.


---

## Body

## 1. Introduction

Wireless communication has been practiced in the industrial market for many years. Today, more and more plants begin to use wireless mesh networks for industrial automation and control system applications. These systems can be considered as wireless mesh network control systems (WMNCSs) [1]. The migration of traditionally wired industrial infrastructure to wireless technologies can improve flexibility, scalability, and efficiency. Typical wireless network technologies suit for WMNCSs include WirelessHART [2], ISA100 [3], and WIA-PA [4].WMNCSs are typical Cyberphysical Systems (CPSs) that are characterized by the tight interactions between computational components, communication networks, and physical dynamics. These components and their interactions are characteristic with stochasticity. The stochasticity arises from the factors as wireless channel condition, network topology, network protocol, routing, data flow configurations, control protocol, and so on. How to model and analyze WMNCSs is still a challenging task.Methods have been proposed to model stochastic processes. Within these methods, stochastic hybrid system (SHS) [5] is a system modeling method that integrates the continuous dynamics with discrete variables and is intensively used to capture the uncertainty in the control system. Stochastic network calculus (SNC) [6] is a theory for network performance analysis that is based on min-plus convolution and stochastic process principles. Statistical models of fading channel [7] are used to capture the variation of wireless channel quality. These methods and models, however, are developed for different application areas, and it still remains gaps for them to be used together.From the perspective of CPSs, we want these methods and tools developed at different area can be used together to solve complex problems in WMNCSs, which, we believe, is a pervasive principle for CPSs. This paper focuses on how to model and analyze WMNCSs with stochastic methods. We model WMNCSs at different abstract levels with SHS, SNC, and statistical models of fading channel. We also bridge the gaps between these methods, so they can be used together to model and analyze WMNCSs seamlessly.The main contributions of this paper are as follows:(1) we propose a three-layered modeling framework that can be used to capture the stochastic properties of WMNCSs at different levels including system behavior, network performance, and wireless signal propagation. To the best of our knowledge, it is the first work that provides such a unified and flexible framework to model and analyze WMNCSs with stochastic methods; (2) with an upper bound approach based on MGFs (moment generating functions) and stochastic orders, we bridge the gaps between SHS, SNC, and statistical models of wireless fading channel. They can be used together to model and analyze WMNCSs seamlessly; (3) we also give an example of how to model and analyze WMNCSs with this framework. The synthesis of above contribution is a methodology for modeling and analysis WMNCSs with stochastic methods.The organization of this paper is as follows. Section2 introduces the background knowledge and describes the problem in detail; Section 3 shows how to model WMNCSs with SHS; Section 4 tailors the SNC for WMNCSs; Section 5 describes the statistical models of wireless fading channel. Section 6 describes how to bridge the gaps between SHS, SNC, and statistical models of fading channel; Section 7 provides an example; Section 8 gives a summary of related works; Section 9 concludes this work.


## 2. Preliminaries

In this section, we first introduce the basic theories and math tools which are used to model WMNCSs, and then we give a description of the problem we are facing when modeling WMNCSs with these methods. Although previous works such as SHS, MGF, and SNC are not new methods, our framework is the first work that tries to combine these methods together to provide a unified and flexible modeling framework for WMNCSs. To bridge the gaps between these methods, we also adopt an upper bound approach based on MGFs and stochastic orders in our framework.


### 2.1. Stochastic Hybrid System

A stochastic hybrid system (SHS) is a system modeling method that integrates the continuous dynamics with discrete variables.Stochastic Hybrid System (SHS).A stochastic hybrid system is defined by four ingredients.( 1 ) A set of discrete modes Q and the system mode q(t) changes through instantaneous resets or impulses:(1)qt,t≥0takesvalueinQ.( 2 ) A set of (locally Lipschitz) vector fields that describes the evolution of the continuous statex(t) in each modeq:(2)x˙=fq,x,tf:Q×Rn×0,∞→Rn.( 3 ) A set of reset maps that determine how the system mode changes (jumps) and when resets happen:(3)q,x=ϕlq-,x-,tϕl:Q×Rn×0,∞→Q×Rn,l∈L.( 4 ) A set of reset time distributions that define the system reset and impulse times:(4)μq,lt:q∈Q,l∈L.Here, the families of reset maps and reset time distributions are parameterized by the same set l∈L=1,⋯,m. When μq,lt≡0, the transition l does not trigger in the modeq. We use a stochastic transition counter to count the number of times that the corresponding discrete transition/reset map ϕl is activated:(5)Nl:Q×0,∞→N.Figure1 shows a graphic representation of an SHS, where each ellipse corresponds to a discrete mode and each edge to a transition between discrete modes. Ellipses are labeled with the corresponding discrete modes, and the vector field that governs the evolution of continuous states is written under it. Edges are labeled with reset maps.Figure 1
Graphical representation of the stochastic hybrid system.When simulating, the SHS starts at an initial conditionq0∈Q,x∈Rn with t=0. Let [0,T) denote the maximum interval of the system, the process is constructed as follows:(1)
Setk=0 and t0=0.(2)
Determine the time staying at this mode, the jump taking to the next mode. Taking time triggered SHS (TTSHS) [5] as an example, at this step each random number hkl is the draw for all l∈L with the reset-time distribution μqtk,l. Let hk≔minl∈Lhkl, if any of hkl is equal to 0 or more than two draws have the same minimum value, the procedure fails. Otherwise, let tk+1≔tk+hk.(3)
Solve the initial value problem withztk=x(tk), z=f(q(tk),z,t)˙, ∀t≥tk at the time interval tk,min⁡tk+1,T. Then the state of the system evolves according to the differential equation corresponding to mode qk for the time interval tk,min⁡tk+1,T.(4)
Iftk+1≥T, the process terminated after the time point T.(5)
Attk+1∈(tk,T), the system jumps to a new mode and a reset transition happened according to the following equation:(6)qtk+1,xtk+1=ϕlkqtk,ztk.Then a new cycle begins and the process goes back to(2) with new draws for all l∈L.


### 2.2. Moment Generating Function

For any random variable or stochastic process, we useFXx=Pr⁡[X≤x] or FXx,t=Pr⁡[X(t)≤x] to denote its cumulative distribution function (CDF), and we use F-Xx=Pr⁡[X>x]or F-Xx,t=Pr⁡[X(t)>x] to denote its complementary cumulative distribution function (CCDF).Moment Generating Function.The MGF of any random variableX is defined for any θ:(7)MXθ=EeθX.Similarly, for any stochastic processesX(t), it has MX(t)θ=E[eθXt] for any θ.MGF is a useful mathematical tool for probability and stochastic process. We list some most useful properties of MGF as below:(8)aPr⁡X≥x≤e-θxEeθXbMc1+c2θ=eθc1Mxc2θcMX+Yθ=MXθMYθMX-Yθ=MXθM-YθdMmin⁡X+Yθ≤min⁡MXθ,MYθMmax⁡X+Yθ≤max⁡M-Xθ,M-YθHere, X and Y are random variables, c1 and c2 are constants, and M-Xθ=MX-θ=E[e-θX].For the min-plus convolution operation⊗ and deconvolution operation ⊘ defined below [8](9)x⊗ys,t=infτ∈s,t⁡xs,τ+yτ,t,x⊘ys,t=supτ∈0,s⁡xτ,t-yτ,s.where x(s,t) and y(s,t) are real-valued, bivariate functions which are used to describe independent random processes X(s,t) and Y(s,t), respectively, and the corresponding min-plus convolution ∗ and deconvolution ° operations for MGFs are defined as(10)MX∗MYτ,t=∑s=τtMXs,τMYτ,t,MX∘MYτ,t=∑τ=0sMXτ,tMYτ,s.We have the following important relationships between these operations [8]:(11)eMX⊗Yθ=M-Xθ∗M-Yθs,t,(12)MX⊘Yθ=MXθ∘M-Yθs,t.(13)fPrX⊗Ys,t≥supθ∈0,∞⁡1θln⁡ε-ln⁡M-Xθ∗M-Yθs,t≥1-ε(14)PrX⊘Ys,t≥infθ∈0,∞⁡1θln⁡M-Xθ∗M-Yθs,t-ln⁡ε≥1-εHere x(s,t) and y(s,t) are bivariate functions defined at 0≤s≤t.


### 2.3. Stochastic Network Calculus

Stochastic network calculus (SNC) is a theory for network performance analysis based on min-plus convolution and stochastic process principles. SNC provides network performance guarantees that have the following form:(15)Pr⁡delay>x≤εorPr⁡loss>x≤ε.Herex represents the targeted delay or loss with the permissible violation probability ε.In SNC, the data arrival process and departure process of a network node are described by bivariate functionsA(s,t) and D(s,t), respectively. A(s,t) gives the amount of data arrives in an interval s,t. D(s,t) gives the amount of data departs in an interval s,t. It is easy to see that A(s,t) and D(s,t) are nonnegative functions, increasing in t, decreasing ins, and At,t=0 and Dt,t=0 for allt.Dynamic Server. In SNC, the dynamic service provided by the node is described by a random process S(s,t), and the following inequality holds:(16)D0,t≥A⊘S0,tLeftover Service. Consider two flows at the same node with arrival processes Anl1(s,t) and Anl2(s,t) that are scheduled as a work-conserving server with the service process Sn(s,t), where Sn(s,t) is nonnegative, increasing int and Sns,s=0. Then the service provided to the flow l2 is(17)Snl2s,t=max⁡0,Sns,t-Anl1s,tAnd the MGF of Snl2s,tis bounded for θ≥0 and t≥s≥0 by(18)M-Snl2θ,s,t≤min⁡1,M-Snθ,s,tM-Anl1θ,s,tConcatenation. Consider two dynamic servers S1(s,t) and Sn2(s,t) in series. There exists an equivalent, single dynamic server S(s,t) for 0≤s≤t where(19)Ss,t=S1⊗S2s,tAnd the MGF of S(s,t), the equivalent single dynamic server is upper bounded for 0≤s≤t according to (20)M-Sθ,s,t≤M-S1θ∗M-S2θs,t.Backlog and Delay Bounds. The backlog b(t) and delay d(t) are defined as:(21)bt=A0,t-D0,t,(22)dt=inf⁡s≥0:A0,t-D0,t+s≤0.And they are bounds given by(23)bt=A⊘St,t,(24)dt=inf⁡s≥0:A⊘Ss,t≤0.(25)Ds,t≤A⊘Ss,t.


### 2.4. Problem Description

Figure2 shows an example of WMNCs. Here three sensors are used to measure the information from two plants, five routers are used to translate information, control tasks are computed on the controller and control commands are sent to the actuators to regulate the plants.Figure 2
An instance of WMNCS.Many WMNCSs have stringent performance requirements for feedback control. However, these WMNCSs are characteristic with stochasticity which arises from the factors as wireless channel condition, network topology, network protocol, routing, data flow configuration, control protocol, etc. How to model and analyze WMNCSs is still a challenging task.In this paper, we set up a three-layered modeling framework to capture the stochastic properties of WMNCSs at different abstract levels. We use SHS to set up a system layer model, SNC to set up a bit layer model, and statistical model for physical layer model. This framework can help us to figure out how the factors as network topology, channel condition, other network traffics, etc. affect the performance of multiple hopped wireless control systems.SHS, SNC, and statistical models of fading channel are developed for different application areas. It has the following gaps for these methods to be used together:(1)
For the periodic sampling events at the sensor node, how to set up the stochastic model and integrate with SNC. It gives the input of SNC model.(2)
When considering the channel fading and coding, how to use statistical models at the physical layer to calculate the service process of SNC model. It is how to integrate the physical layer model with bit layer model.(3)
To analyze the stability of SHS, it usually needs to know their set-time distributions, but the SNC model only provides the probability guarantees. It is how to integrate bit layer model with the system model.In this work, we also bridge these gaps. We emphasize the unified and flexible properties of this framework. We take a unified upper bound approach based on MGFs and stochastic orders, so these models can be used together to analyze system performance. The flexibility embodied in the methods we selected at each layer, so we can solve the same problem with little effort when the conditions as network topology, channel condition, or data flow configuration are changed. The unified and flexible properties can help us to try different design combinations easily and support quick design decisions.


## 2.1. Stochastic Hybrid System

A stochastic hybrid system (SHS) is a system modeling method that integrates the continuous dynamics with discrete variables.Stochastic Hybrid System (SHS).A stochastic hybrid system is defined by four ingredients.( 1 ) A set of discrete modes Q and the system mode q(t) changes through instantaneous resets or impulses:(1)qt,t≥0takesvalueinQ.( 2 ) A set of (locally Lipschitz) vector fields that describes the evolution of the continuous statex(t) in each modeq:(2)x˙=fq,x,tf:Q×Rn×0,∞→Rn.( 3 ) A set of reset maps that determine how the system mode changes (jumps) and when resets happen:(3)q,x=ϕlq-,x-,tϕl:Q×Rn×0,∞→Q×Rn,l∈L.( 4 ) A set of reset time distributions that define the system reset and impulse times:(4)μq,lt:q∈Q,l∈L.Here, the families of reset maps and reset time distributions are parameterized by the same set l∈L=1,⋯,m. When μq,lt≡0, the transition l does not trigger in the modeq. We use a stochastic transition counter to count the number of times that the corresponding discrete transition/reset map ϕl is activated:(5)Nl:Q×0,∞→N.Figure1 shows a graphic representation of an SHS, where each ellipse corresponds to a discrete mode and each edge to a transition between discrete modes. Ellipses are labeled with the corresponding discrete modes, and the vector field that governs the evolution of continuous states is written under it. Edges are labeled with reset maps.Figure 1
Graphical representation of the stochastic hybrid system.When simulating, the SHS starts at an initial conditionq0∈Q,x∈Rn with t=0. Let [0,T) denote the maximum interval of the system, the process is constructed as follows:(1)
Setk=0 and t0=0.(2)
Determine the time staying at this mode, the jump taking to the next mode. Taking time triggered SHS (TTSHS) [5] as an example, at this step each random number hkl is the draw for all l∈L with the reset-time distribution μqtk,l. Let hk≔minl∈Lhkl, if any of hkl is equal to 0 or more than two draws have the same minimum value, the procedure fails. Otherwise, let tk+1≔tk+hk.(3)
Solve the initial value problem withztk=x(tk), z=f(q(tk),z,t)˙, ∀t≥tk at the time interval tk,min⁡tk+1,T. Then the state of the system evolves according to the differential equation corresponding to mode qk for the time interval tk,min⁡tk+1,T.(4)
Iftk+1≥T, the process terminated after the time point T.(5)
Attk+1∈(tk,T), the system jumps to a new mode and a reset transition happened according to the following equation:(6)qtk+1,xtk+1=ϕlkqtk,ztk.Then a new cycle begins and the process goes back to(2) with new draws for all l∈L.


## 2.2. Moment Generating Function

For any random variable or stochastic process, we useFXx=Pr⁡[X≤x] or FXx,t=Pr⁡[X(t)≤x] to denote its cumulative distribution function (CDF), and we use F-Xx=Pr⁡[X>x]or F-Xx,t=Pr⁡[X(t)>x] to denote its complementary cumulative distribution function (CCDF).Moment Generating Function.The MGF of any random variableX is defined for any θ:(7)MXθ=EeθX.Similarly, for any stochastic processesX(t), it has MX(t)θ=E[eθXt] for any θ.MGF is a useful mathematical tool for probability and stochastic process. We list some most useful properties of MGF as below:(8)aPr⁡X≥x≤e-θxEeθXbMc1+c2θ=eθc1Mxc2θcMX+Yθ=MXθMYθMX-Yθ=MXθM-YθdMmin⁡X+Yθ≤min⁡MXθ,MYθMmax⁡X+Yθ≤max⁡M-Xθ,M-YθHere, X and Y are random variables, c1 and c2 are constants, and M-Xθ=MX-θ=E[e-θX].For the min-plus convolution operation⊗ and deconvolution operation ⊘ defined below [8](9)x⊗ys,t=infτ∈s,t⁡xs,τ+yτ,t,x⊘ys,t=supτ∈0,s⁡xτ,t-yτ,s.where x(s,t) and y(s,t) are real-valued, bivariate functions which are used to describe independent random processes X(s,t) and Y(s,t), respectively, and the corresponding min-plus convolution ∗ and deconvolution ° operations for MGFs are defined as(10)MX∗MYτ,t=∑s=τtMXs,τMYτ,t,MX∘MYτ,t=∑τ=0sMXτ,tMYτ,s.We have the following important relationships between these operations [8]:(11)eMX⊗Yθ=M-Xθ∗M-Yθs,t,(12)MX⊘Yθ=MXθ∘M-Yθs,t.(13)fPrX⊗Ys,t≥supθ∈0,∞⁡1θln⁡ε-ln⁡M-Xθ∗M-Yθs,t≥1-ε(14)PrX⊘Ys,t≥infθ∈0,∞⁡1θln⁡M-Xθ∗M-Yθs,t-ln⁡ε≥1-εHere x(s,t) and y(s,t) are bivariate functions defined at 0≤s≤t.


## 2.3. Stochastic Network Calculus

Stochastic network calculus (SNC) is a theory for network performance analysis based on min-plus convolution and stochastic process principles. SNC provides network performance guarantees that have the following form:(15)Pr⁡delay>x≤εorPr⁡loss>x≤ε.Herex represents the targeted delay or loss with the permissible violation probability ε.In SNC, the data arrival process and departure process of a network node are described by bivariate functionsA(s,t) and D(s,t), respectively. A(s,t) gives the amount of data arrives in an interval s,t. D(s,t) gives the amount of data departs in an interval s,t. It is easy to see that A(s,t) and D(s,t) are nonnegative functions, increasing in t, decreasing ins, and At,t=0 and Dt,t=0 for allt.Dynamic Server. In SNC, the dynamic service provided by the node is described by a random process S(s,t), and the following inequality holds:(16)D0,t≥A⊘S0,tLeftover Service. Consider two flows at the same node with arrival processes Anl1(s,t) and Anl2(s,t) that are scheduled as a work-conserving server with the service process Sn(s,t), where Sn(s,t) is nonnegative, increasing int and Sns,s=0. Then the service provided to the flow l2 is(17)Snl2s,t=max⁡0,Sns,t-Anl1s,tAnd the MGF of Snl2s,tis bounded for θ≥0 and t≥s≥0 by(18)M-Snl2θ,s,t≤min⁡1,M-Snθ,s,tM-Anl1θ,s,tConcatenation. Consider two dynamic servers S1(s,t) and Sn2(s,t) in series. There exists an equivalent, single dynamic server S(s,t) for 0≤s≤t where(19)Ss,t=S1⊗S2s,tAnd the MGF of S(s,t), the equivalent single dynamic server is upper bounded for 0≤s≤t according to (20)M-Sθ,s,t≤M-S1θ∗M-S2θs,t.Backlog and Delay Bounds. The backlog b(t) and delay d(t) are defined as:(21)bt=A0,t-D0,t,(22)dt=inf⁡s≥0:A0,t-D0,t+s≤0.And they are bounds given by(23)bt=A⊘St,t,(24)dt=inf⁡s≥0:A⊘Ss,t≤0.(25)Ds,t≤A⊘Ss,t.


## 2.4. Problem Description

Figure2 shows an example of WMNCs. Here three sensors are used to measure the information from two plants, five routers are used to translate information, control tasks are computed on the controller and control commands are sent to the actuators to regulate the plants.Figure 2
An instance of WMNCS.Many WMNCSs have stringent performance requirements for feedback control. However, these WMNCSs are characteristic with stochasticity which arises from the factors as wireless channel condition, network topology, network protocol, routing, data flow configuration, control protocol, etc. How to model and analyze WMNCSs is still a challenging task.In this paper, we set up a three-layered modeling framework to capture the stochastic properties of WMNCSs at different abstract levels. We use SHS to set up a system layer model, SNC to set up a bit layer model, and statistical model for physical layer model. This framework can help us to figure out how the factors as network topology, channel condition, other network traffics, etc. affect the performance of multiple hopped wireless control systems.SHS, SNC, and statistical models of fading channel are developed for different application areas. It has the following gaps for these methods to be used together:(1)
For the periodic sampling events at the sensor node, how to set up the stochastic model and integrate with SNC. It gives the input of SNC model.(2)
When considering the channel fading and coding, how to use statistical models at the physical layer to calculate the service process of SNC model. It is how to integrate the physical layer model with bit layer model.(3)
To analyze the stability of SHS, it usually needs to know their set-time distributions, but the SNC model only provides the probability guarantees. It is how to integrate bit layer model with the system model.In this work, we also bridge these gaps. We emphasize the unified and flexible properties of this framework. We take a unified upper bound approach based on MGFs and stochastic orders, so these models can be used together to analyze system performance. The flexibility embodied in the methods we selected at each layer, so we can solve the same problem with little effort when the conditions as network topology, channel condition, or data flow configuration are changed. The unified and flexible properties can help us to try different design combinations easily and support quick design decisions.


## 3. System Layer Model

In this section, we concentrate on how to model the feedback control behavior of WMNCSs at the system level with SHS. We use SHS to capture both the continuous dynamics and the discrete logic of WMNCSs with a simulation approach.


### 3.1. Modeling Controller with Stochastic Hybrid System

We consider an SISO (single input and single output) WMNCS that consists of a sensor, an actuator, and a controller that are connected through a shared wireless network. Figure3 shows the structure of this system. Here the sensor and actuator are connected to the plant directly, and the controller communicates with the sensor and actuator through a shared wireless network.Figure 3
The control structure of our WMNCS.When considering the linear time-invariant models for the plant and controller, the system can be described by(26)Plant:x˙p=Apxp+Bpu^,y=Cpxp.Controller:x˙c=Acxc+Bcy^,u=Ccxc+Dcy^.where, in the plant model, xp denotes the state of the plant, u^ denotes the most recently received control variable, andy denotes the output of the plant. In the controller model, xc denotes the state of the controller, u denotes the controller output, and y^ denotes the most recently received output of the plant. The controller is assumed to yield the closed-loop stable when the plant and the controller are directly connected, i.e., when y^t=y(t) and u^t=ut. We use tn,k to denote the times at which thek-th data of sensor datay or control commandu (n=y oru) is received by the receiver, and the time intervals tn,k,tn,k+1,k∈N are assumed to be identically distributed with distribution μn. Between these times, y^ and u^ are held constant.We use SHS to model above WMNCS. The key to the model network control system with SHS is to capture the sequence of events as a finite state machine, which is mapped to the discrete modes of SHS. We follow a simulation way and give two examples [5].Case 1.
When the controller and actuator are event-triggered, packet drops are not considered.Here we consider a WMNCS with one sensor, one controller, and one actuator which are all events triggered. The control loop length is set toTs. Under the normal state, sensor node samples at the beginning of every control loop tk,k∈Z≥1. For event triggered scheme, the controller and the actuator update their internal states when they receive new data, and the actuator takes action when receiving command. Further, we assume that Ts is long enough and the following two assumptions are satisfied.Assumption 1.
As packet drops are not considered, the data packets will reach the destination finally. The value of sensor to controller delayτsc and controller to actuator delay τca are lower and upper bounded. They are formulated as (27)τsc∈Tscmin,Tscmaxandτca∈Tcamin,Tcamax.Here, Tscmin and Tcamin are the lower bounds of τsc and τca, respectively, for the minimum time spent on computation and communication. Tscmax and Tcamax are the upper bounds of τsc and τca, which can be guaranteed by network technologies as WirelessHART and ISA100 without the considering of packet losses.Assumption 2.
We assume that the following inequation always holds:(28)Tscmax+Tcamax<Ts+Tscmin.This assumption guarantees the time sequence of mode jumps can be preserved. As is shown in the Figure 4, if Tscmax+Tcamax≥Ts+Tscmin, at the third control loop, the controller may receive sampling data before the actuator receives the command of the previous control loop. For a violation example, let us look at the state jump indicated by dotted arrow line, which means a jump (2,Jcx) may happen at the mode q=2.Figure 4
Time sequence that preserves the event orders.Under the above two assumptions, the discrete modes and transitions can be listed as below.Discrete Modes. Two modes are required to model this case.Mode 1, after the actuator updated its inner state, the controller is waiting for a new sample data.Mode 2, after the controller received a new sampling data and updated its inner state, the actuator is waiting for a new command from the controller.Jumps. Here we need two mode jumps.Jump 1, for the sensor sampling data arrival controller, the controller updates the state estimator by receivedy^, it has(29)xtk=Jcx-tk,whereJc=I0000I00Cp0Dc0000I,x≔xpxcy^u^.Jump 2, for the command packet arrival controller, the actuator updates its state and uses the latest value to actuate the plant. Under this situation, it has the jump:(30)xtk=Jax-tk,Ja=I0000I0000I00Cc00.In any intervaltk,tk+1,k∈z≥0 between jumps, the system dynamic can be written as(31)x˙=AxwhereA≔Ap00Bp0AcBc000000000.The above equations express the fact thatxpand xc are continuous signals and, at sampling times, y^ and u^ are updated to the value of y^=Cpxp,u^=Ccxc+Dcy^ or remained with the same value depending on jump condition.Figure5 shows the graphical representation of this SHS. The system has two discrete states (models), and two different reset maps defined by the matrices Jc and Ja. The reset maps are ϕ1q,x=(2,Jcx) and ϕ2q,x=(1,Jax) and set-time distributions μ1,1 and μ2,2 are associated with them, respectively. This SHS can be regarded as a special case of TTSHS, for which there is only one jump at each mode.Figure 5
SHS model when packet drops are not considered.For simulation, the sensor takes a sample at the beginning of the control loop. Afterτsc, the controller receives the sampling data, updates its inner state, and the system states jump to (2,Jcx). Then controller computes and sends control commands to the actuator. After τca, the actuator receivers the control command, updates its inner states, and the system jump to (1,Jax). Then a new cycle begins with the sensor taking a new sample at the beginning of the next control loop.Case 2.
When the controller and actuator are event-triggered, packet drops are considered.In this situation, we use two more discrete modes to present packet drops at the controller and actuator. Figure6 shows the SHS model. Its modes and jumps are listed as below.Figure 6
SHS model when packet drops are considered.Discrete Modes. Four modes are required to model this case.Mode 1, after the actuator successfully received the command and updated its inner state, the controller is waiting for a new sample data from the sensor.Mode 2, after the controller successfully received a new sampling data and updated its inner state, the actuator is waiting for new commands.Mode 3, the actuator failed to receive the command due to the packet drop, and the controller is waiting for a new sample data from the sensor.Mode 4, the controller failed to receive a new sampling data, and the actuator is waiting for a new command from the controller.Jumps. There are eight jumps between discrete modes, which can be classified into three groups:Group 1, jumps triggered by the controller receiving a new sampling data, and the controller updates its inner state withxtk=Jcx-(tk).Group 2, jumps triggered by the actuator receiving a new sampling data, and the actuator updates its inner state withxtk=Jax-(tk).Group 3, jumps triggered by the controller or the actuator failed to receive a new sampling data, and the system is updated with(32)xtk=Jdx-tk,whereJd=I0000I0000I0000I.In Figure6, every mode including two jump outs, one for successful transportation and one for packet drop. This process can be modeled by a Markov chain:(33)Pr⁡modetk≔j∣modett-≔i=pi,j1≤i,j≤4.We use pcdto indicate the probability packet drops at the controller and padto indicate the probability packet drops at the actuator, and the matrix of the Markov chain is given by(34)01-pcd0pcd1-pad0pad001-pcd0pcd1-pad0pad0In any intervaltk,tk+1,k∈z≥0between jumps, the system dynamic follows the law of x˙=Ax.For simulation, the sensor takes a sample at the beginning of the control loop. Then afterτsc, it draws a coin with a probability pcd for packet drop. If a successful transmission occurs, the controller receives the sampling data and updates its inner states, and the system jumps to (2,Jcx), or else the system jumps to (4,Jdx). Then controller computes and sends control commands to the actuator. After τca, it draws a coin with a probability pad for packet drop. If a successful transmission occurs, the actuator receives the control command and updates its inner states, and the system states jump to (1,Jax); otherwise the system jumps to (3,Jdx). When the system has entered into mode 3 (mode 4), it waits for τsc ( τca) time, draws a coin, and jumps to 2,Jcx or (4,Jdx) (1,Jax or (3,Jdx) ).


### 3.2. Stability Analysis of SHS

We consider impulsive systems with several reset maps triggered by independent renewal processes; i.e., the intervals between jumps associated with a given reset map are identically distributed and independent of the other jump intervals. Considering the linear dynamic and reset maps, we establish that mean exponential stability (MES) is equivalent to the spectral radius of an integral operator being less than one.Mean Exponentially Stable (MES). A SHS is mean exponentially stable if there exist constants c>0 and α>0 such that, for any initial condition x0 of the system, we have(35)ExtTxt≤ce-αtx0Tx0,∀t≥0Assumption 3 (see [9]).
( 1 ) The transition distributions μq,lt:q∈Q,l∈L are measurable functions.
( 2 ) ∀ q ∈ Q , x ∈ R n , l ∈ L and t≥0, that the following inequations are satisfied:(36)fq,x,t≤max⁡φftx,cf,ϕ-lq,x,t≤max⁡x,cϕ.where φf:0,∞→0,∞ is a continuous function, cfand cϕ are constants, and ϕ-l is the projection function of ϕl to Rn (ϕ-lq,x-,t=x ).Itô Equation and Extended Generator [9]. Under Assumption 3, we have Itô equation and extended generator for SHS. They are tools to compute the expectations of SHS. The Itô equation is (37)dψq,x,t=∂ψq,x,t∂xfq,x,tdt+∂ψq,x,t∂tdt+∑l=1nlψφlq,x,t,t-ψq,x,tdNlqThe extended generator of SHS is defined as(38)Lψq,x,t≔∂ψq,x,t∂xfq,x,t+∂ψq,x,t∂t+∑l=1nlψφlq,x,t,t-ψq,x,tdNlqdtThe Expectation of SHS [9]. For every function ψq,x,t:Q×Rn×[0,∞)→Rn that continuously differentiable with respect tox, t, and initial condition z0:(q0,x0,t0), the expectation on the state of SHS is given by(39)Ez0ψqt,xt,t=ψq0,x0,t0+Ez0∫t0tLψqs,xs,sdsTheorem 4.
For the TTSHS, the system is MSS if only if the inequationrσΓ<1 holds, where(40)Γ=Γ1,1⋯Γ1,nq⋮⋱⋮Γnq,1⋯Γnq,nq,Γi,j=∑l=1nl∫0∞Ji,leAST⊗Ji,leASTχξli=jrisri,lsμi,lds,χx∈A equals 1 if x∈A and 0 otherwise, (41)ri,ls=∏l=1nlri,ls,ri,ls=1,ifμi,lt≡0μi,ls,∞,otherwise,
r σ Γ is the spectral radius of operator Γ.
Please refer to [10] for proof.


## 3.1. Modeling Controller with Stochastic Hybrid System

We consider an SISO (single input and single output) WMNCS that consists of a sensor, an actuator, and a controller that are connected through a shared wireless network. Figure3 shows the structure of this system. Here the sensor and actuator are connected to the plant directly, and the controller communicates with the sensor and actuator through a shared wireless network.Figure 3
The control structure of our WMNCS.When considering the linear time-invariant models for the plant and controller, the system can be described by(26)Plant:x˙p=Apxp+Bpu^,y=Cpxp.Controller:x˙c=Acxc+Bcy^,u=Ccxc+Dcy^.where, in the plant model, xp denotes the state of the plant, u^ denotes the most recently received control variable, andy denotes the output of the plant. In the controller model, xc denotes the state of the controller, u denotes the controller output, and y^ denotes the most recently received output of the plant. The controller is assumed to yield the closed-loop stable when the plant and the controller are directly connected, i.e., when y^t=y(t) and u^t=ut. We use tn,k to denote the times at which thek-th data of sensor datay or control commandu (n=y oru) is received by the receiver, and the time intervals tn,k,tn,k+1,k∈N are assumed to be identically distributed with distribution μn. Between these times, y^ and u^ are held constant.We use SHS to model above WMNCS. The key to the model network control system with SHS is to capture the sequence of events as a finite state machine, which is mapped to the discrete modes of SHS. We follow a simulation way and give two examples [5].Case 1.
When the controller and actuator are event-triggered, packet drops are not considered.Here we consider a WMNCS with one sensor, one controller, and one actuator which are all events triggered. The control loop length is set toTs. Under the normal state, sensor node samples at the beginning of every control loop tk,k∈Z≥1. For event triggered scheme, the controller and the actuator update their internal states when they receive new data, and the actuator takes action when receiving command. Further, we assume that Ts is long enough and the following two assumptions are satisfied.Assumption 1.
As packet drops are not considered, the data packets will reach the destination finally. The value of sensor to controller delayτsc and controller to actuator delay τca are lower and upper bounded. They are formulated as (27)τsc∈Tscmin,Tscmaxandτca∈Tcamin,Tcamax.Here, Tscmin and Tcamin are the lower bounds of τsc and τca, respectively, for the minimum time spent on computation and communication. Tscmax and Tcamax are the upper bounds of τsc and τca, which can be guaranteed by network technologies as WirelessHART and ISA100 without the considering of packet losses.Assumption 2.
We assume that the following inequation always holds:(28)Tscmax+Tcamax<Ts+Tscmin.This assumption guarantees the time sequence of mode jumps can be preserved. As is shown in the Figure 4, if Tscmax+Tcamax≥Ts+Tscmin, at the third control loop, the controller may receive sampling data before the actuator receives the command of the previous control loop. For a violation example, let us look at the state jump indicated by dotted arrow line, which means a jump (2,Jcx) may happen at the mode q=2.Figure 4
Time sequence that preserves the event orders.Under the above two assumptions, the discrete modes and transitions can be listed as below.Discrete Modes. Two modes are required to model this case.Mode 1, after the actuator updated its inner state, the controller is waiting for a new sample data.Mode 2, after the controller received a new sampling data and updated its inner state, the actuator is waiting for a new command from the controller.Jumps. Here we need two mode jumps.Jump 1, for the sensor sampling data arrival controller, the controller updates the state estimator by receivedy^, it has(29)xtk=Jcx-tk,whereJc=I0000I00Cp0Dc0000I,x≔xpxcy^u^.Jump 2, for the command packet arrival controller, the actuator updates its state and uses the latest value to actuate the plant. Under this situation, it has the jump:(30)xtk=Jax-tk,Ja=I0000I0000I00Cc00.In any intervaltk,tk+1,k∈z≥0 between jumps, the system dynamic can be written as(31)x˙=AxwhereA≔Ap00Bp0AcBc000000000.The above equations express the fact thatxpand xc are continuous signals and, at sampling times, y^ and u^ are updated to the value of y^=Cpxp,u^=Ccxc+Dcy^ or remained with the same value depending on jump condition.Figure5 shows the graphical representation of this SHS. The system has two discrete states (models), and two different reset maps defined by the matrices Jc and Ja. The reset maps are ϕ1q,x=(2,Jcx) and ϕ2q,x=(1,Jax) and set-time distributions μ1,1 and μ2,2 are associated with them, respectively. This SHS can be regarded as a special case of TTSHS, for which there is only one jump at each mode.Figure 5
SHS model when packet drops are not considered.For simulation, the sensor takes a sample at the beginning of the control loop. Afterτsc, the controller receives the sampling data, updates its inner state, and the system states jump to (2,Jcx). Then controller computes and sends control commands to the actuator. After τca, the actuator receivers the control command, updates its inner states, and the system jump to (1,Jax). Then a new cycle begins with the sensor taking a new sample at the beginning of the next control loop.Case 2.
When the controller and actuator are event-triggered, packet drops are considered.In this situation, we use two more discrete modes to present packet drops at the controller and actuator. Figure6 shows the SHS model. Its modes and jumps are listed as below.Figure 6
SHS model when packet drops are considered.Discrete Modes. Four modes are required to model this case.Mode 1, after the actuator successfully received the command and updated its inner state, the controller is waiting for a new sample data from the sensor.Mode 2, after the controller successfully received a new sampling data and updated its inner state, the actuator is waiting for new commands.Mode 3, the actuator failed to receive the command due to the packet drop, and the controller is waiting for a new sample data from the sensor.Mode 4, the controller failed to receive a new sampling data, and the actuator is waiting for a new command from the controller.Jumps. There are eight jumps between discrete modes, which can be classified into three groups:Group 1, jumps triggered by the controller receiving a new sampling data, and the controller updates its inner state withxtk=Jcx-(tk).Group 2, jumps triggered by the actuator receiving a new sampling data, and the actuator updates its inner state withxtk=Jax-(tk).Group 3, jumps triggered by the controller or the actuator failed to receive a new sampling data, and the system is updated with(32)xtk=Jdx-tk,whereJd=I0000I0000I0000I.In Figure6, every mode including two jump outs, one for successful transportation and one for packet drop. This process can be modeled by a Markov chain:(33)Pr⁡modetk≔j∣modett-≔i=pi,j1≤i,j≤4.We use pcdto indicate the probability packet drops at the controller and padto indicate the probability packet drops at the actuator, and the matrix of the Markov chain is given by(34)01-pcd0pcd1-pad0pad001-pcd0pcd1-pad0pad0In any intervaltk,tk+1,k∈z≥0between jumps, the system dynamic follows the law of x˙=Ax.For simulation, the sensor takes a sample at the beginning of the control loop. Then afterτsc, it draws a coin with a probability pcd for packet drop. If a successful transmission occurs, the controller receives the sampling data and updates its inner states, and the system jumps to (2,Jcx), or else the system jumps to (4,Jdx). Then controller computes and sends control commands to the actuator. After τca, it draws a coin with a probability pad for packet drop. If a successful transmission occurs, the actuator receives the control command and updates its inner states, and the system states jump to (1,Jax); otherwise the system jumps to (3,Jdx). When the system has entered into mode 3 (mode 4), it waits for τsc ( τca) time, draws a coin, and jumps to 2,Jcx or (4,Jdx) (1,Jax or (3,Jdx) ).


## 3.2. Stability Analysis of SHS

We consider impulsive systems with several reset maps triggered by independent renewal processes; i.e., the intervals between jumps associated with a given reset map are identically distributed and independent of the other jump intervals. Considering the linear dynamic and reset maps, we establish that mean exponential stability (MES) is equivalent to the spectral radius of an integral operator being less than one.Mean Exponentially Stable (MES). A SHS is mean exponentially stable if there exist constants c>0 and α>0 such that, for any initial condition x0 of the system, we have(35)ExtTxt≤ce-αtx0Tx0,∀t≥0Assumption 3 (see [9]).
( 1 ) The transition distributions μq,lt:q∈Q,l∈L are measurable functions.
( 2 ) ∀ q ∈ Q , x ∈ R n , l ∈ L and t≥0, that the following inequations are satisfied:(36)fq,x,t≤max⁡φftx,cf,ϕ-lq,x,t≤max⁡x,cϕ.where φf:0,∞→0,∞ is a continuous function, cfand cϕ are constants, and ϕ-l is the projection function of ϕl to Rn (ϕ-lq,x-,t=x ).Itô Equation and Extended Generator [9]. Under Assumption 3, we have Itô equation and extended generator for SHS. They are tools to compute the expectations of SHS. The Itô equation is (37)dψq,x,t=∂ψq,x,t∂xfq,x,tdt+∂ψq,x,t∂tdt+∑l=1nlψφlq,x,t,t-ψq,x,tdNlqThe extended generator of SHS is defined as(38)Lψq,x,t≔∂ψq,x,t∂xfq,x,t+∂ψq,x,t∂t+∑l=1nlψφlq,x,t,t-ψq,x,tdNlqdtThe Expectation of SHS [9]. For every function ψq,x,t:Q×Rn×[0,∞)→Rn that continuously differentiable with respect tox, t, and initial condition z0:(q0,x0,t0), the expectation on the state of SHS is given by(39)Ez0ψqt,xt,t=ψq0,x0,t0+Ez0∫t0tLψqs,xs,sdsTheorem 4.
For the TTSHS, the system is MSS if only if the inequationrσΓ<1 holds, where(40)Γ=Γ1,1⋯Γ1,nq⋮⋱⋮Γnq,1⋯Γnq,nq,Γi,j=∑l=1nl∫0∞Ji,leAST⊗Ji,leASTχξli=jrisri,lsμi,lds,χx∈A equals 1 if x∈A and 0 otherwise, (41)ri,ls=∏l=1nlri,ls,ri,ls=1,ifμi,lt≡0μi,ls,∞,otherwise,
r σ Γ is the spectral radius of operator Γ.
Please refer to [10] for proof.


## 4. Bit Layer Model

In typical network-layer performance models, traffic and service are usually measured in bits; we refer to these models as bit layer models. In this section, we tailor SNC to model and analyze the performance of WMNCSs. We choose SNC because it satisfies the basic properties for network performance analysis as service guarantees, output characterization, concatenation property, leftover service property, and superposition property [6]. These properties will provide enough flexibility for our modeling and analysis.We consider a data flow across a wireless multihop network and compute its end-to-end delay. As is shown in Figure7, in a wireless N-node tandem network, there is a data flow traversing the entire network which may encounter cross traffic at each node. Here, the cross-traffic at a node is the aggregate of all traffic that traverses the node but does not belong to the considered flow. For the most common and useful strategy, all these data flows are assigned with a fixed priority and follow a FIFO scheduling at servers.Figure 7
Tandem network with cross traffics.We use a random process to describe the service given to the considered flow at a node. This service process is governed by the instantaneous channel capacity and the cross traffic at the node. Here we use a discrete-time domain whereti=i∆t,i∈Z is the index of time lot and ∆t is the time slot length. The system is assumed to start with empty queues at timet= 0.For a flow indexed byl, we use random process Anl to represent the cumulative arrivals at noden, Snl to represent the service offered by noden, and Dnl to represent the departures from the noden. We firstly setup a single service node model, and then the multihop model. Throughout, we assume that the arrival and service processes satisfy stationary bounds.


### 4.1. Single Service Node Model

Here, considering a single service node, the cumulative arrivals and departures of data flowl within [s,t) are formulated as(42)Anls,t=∑i=st-1an,ilandDnls,t=∑i=st-1dn,il.where an,il is the arrivals and dn,il the departures in thei-th time slot. Due to causality, we always have (43)Anl0,t>Dnl0,t.The dynamic service of the node n given to the flowl in the time interval [τ,t) is represented by the random process Snl(τ,t).For a stable queuing system, the average arrival rate is smaller than the average service rate. Stationary performance bounds can be obtained only when the system satisfies following stability condition:(44)limt→∞Anl0,tt<limt→∞⁡Snl0,tt.The service that relates the departures of a system to its arrivals is(45)Dnl0,t≥minτ∈0,t⁡Anl0,τ+Snlτ,t=Anl⊘Snlt.The backlog of the flowl at noden at time t>0 is given by(46)bnlt=Anl0,t-Dnl0,t.The delay of the flowl at noden is given by(47)dnlt=inf⁡u≥0:Anl0,t≤Dnl0,t+u.AsAnl(τ,t) and Snl(τ,t) are stochastic bivariate functions, the MGFs of Anl(τ,t) and Snl(τ,t) for any θ are given as(48)MAnlθ=EeθAnlτ,t,MSnlθ=EeθSnlτ,t.Theorem 5.
Given a single service node where its arrival flows are described byAnl(τ,t), and its available service is given by a dynamic server Snl(τ,t). Fix ϵ>0 and define(49)Mnθ,s,t=MAnlθ,s,t∘M-Snlθ,s,tThen we have the following performance bounds.Output burstiness:PrDnlτ,t>Dετ,t≤ε, where(50)Dετ,t=infθ∈0,∞⁡1θln⁡Mnθ,s,t-ln⁡εBacklog: Prbnlt>bεt≤ε, where(51)bεt=infθ∈0,∞⁡1θln⁡Mnθ,t,t-ln⁡εDelay: Prdnlt>dεt≤ε, where dεt is the smallest number satisfying(52)infθ∈0,∞⁡Mnθ,t+dεt,t≤εProof.
From the inequationDnls,t≤Anl⊘Snl(s,t), we have (53)Pr⁡Dnls,t>Dε≤Pr⁡Anl⊘Snls,t>DεComparing with inequation (14), it is easy to get the results for output burstiness.
Frombnlt≤Anl⊘Snl(t,t), we have (54)Pr⁡bt,t>bε≤Pr⁡Anl⊘Snlt,t>bεComparing with (14) and replacings witht, we have the output backlog bounds.
From the definition ofdε, it is easy to get(55)1θMnθ,t+dε,t-ln⁡ε≤0.As dnl=inf⁡τ≥0:Anl⊘Snlt,t+τ<0, we can get(56)Pr⁡dnlt>dε≤Mθ,t+dε,t,and the result is obvious.


### 4.2. Multiple Hop Flow Model

For multiple hop situations, the stability condition is expressed as(57)limt→∞⁡∑l∈FlownAnl0,tt<limt→∞⁡Sn0,tt.For flowl, we want to compute its performance over multiple hops. At first, we introduce the following theorem.Theorem 6.
Considering a flowl passes a work-conserving server with the service process Sn(s,t) that Sn(s,t) is nonnegative, increasing int, and Sns,s=0. Then the service provided to flowl is given by(58)Snls,t=max⁡0,Sns,t-∑i∈PrevFlown,lAnms,t.HereFlow(n) gives all the data follows that pass the noden, and Prev(ϑ,l) is the flows in the flow set ϑ with the priority higher than that ofl.The MGF ofSnl(s,t) is bounded for θ≥0,t≥s≥0 by(59)M-Snlθ,s,t≤min⁡1,M-Snθ,s,t·∏m∈PrevFlown,lMAnmθ,s,tProof.
For a fixed priority flowl, we only consider the flows with a priority higher than that ofl. The result (58) is proved by applying (17) several times. The MGF of Snl(s,t) can be deduced directly from (18).Theorem 7.
Consider a flowl charactering by Al(s,t) passes nodes i1,i2,⋯,in in series. There exists an equivalent, single dynamic server Sl(s,t) for 0≤s≤t where(60)Sls,t=Si1l⊗Si2l⊗⋯⊗Sinls,tAnd the MGF of Sl(s,t), the equivalent, single dynamic server is upper bounded for 0≤s≤t according to(61)M-slθ,s,t≤M-si1lθ∗M-si2lθ∗⋯∗M-sinlθs,tProof.
Equation (60) can be deduced directly by applying (19) several times and the MGF upper bound from (20).Corollary 8.
Consider a flowl   Al(s,t) passes a cascade of N servers i1,i2,⋯,in; it has(62)M-Al⊘Sl≤MAlθ∘M-Si1lθ∗M-Si2lθ∗⋯∗M-Sinlθs,t=MAlθ∘M-Si1lθ∘M-Si2lθ∘⋯∘M-Sinlθs,tProof.
From the definition of operators∗ and ∘, it has(63)x∘y∗zs,t=∑τ1=0∞xt-s+τ1∑τ2τ1yτ1-τ2zτ2=∑u1=0∞∑u2=0∞xt-s+u1+u2yu1zu2=x∘y∘zs,t.The result is gotten directly from (61) by applying above equation several times.The computation network performance analysis process for a WMNCS is given as Algorithm1.Algorithm 1:Compute end-to-end delays and MGFs.
( 1 )   begin ( 2 )  Initialize the left-over service Sileft(s,t) and its MGF upper bound of each node n as Si(s,t) and MSi(θ); ( 3 )  for each  lin the flow set in descending priority order ( 4 )   for eachi inNode(l) in the order of passing ( 5 )    Get the cumulative arrival Ail(s,t) and its MGF upper bound; ( 6 )    Compute the cumulative departure Dil(s,t) and MGF upper bound as next node cumulative arrival; ( 7 )    Update the left-over service Sileft(s,t) and its MGF upper bound at node i; ( 8 )   end for ( 9 )   Compute the tandem service for flow l and its MGF upper bound; ( 10 ) Compute the tandem delay and its MGF upper bound; ( 11 )  end for ( 12 )  end beginHere,Node(l) denotes the node set that flowl passes, and Sileft denotes the leave over service of nodei.


## 4.1. Single Service Node Model

Here, considering a single service node, the cumulative arrivals and departures of data flowl within [s,t) are formulated as(42)Anls,t=∑i=st-1an,ilandDnls,t=∑i=st-1dn,il.where an,il is the arrivals and dn,il the departures in thei-th time slot. Due to causality, we always have (43)Anl0,t>Dnl0,t.The dynamic service of the node n given to the flowl in the time interval [τ,t) is represented by the random process Snl(τ,t).For a stable queuing system, the average arrival rate is smaller than the average service rate. Stationary performance bounds can be obtained only when the system satisfies following stability condition:(44)limt→∞Anl0,tt<limt→∞⁡Snl0,tt.The service that relates the departures of a system to its arrivals is(45)Dnl0,t≥minτ∈0,t⁡Anl0,τ+Snlτ,t=Anl⊘Snlt.The backlog of the flowl at noden at time t>0 is given by(46)bnlt=Anl0,t-Dnl0,t.The delay of the flowl at noden is given by(47)dnlt=inf⁡u≥0:Anl0,t≤Dnl0,t+u.AsAnl(τ,t) and Snl(τ,t) are stochastic bivariate functions, the MGFs of Anl(τ,t) and Snl(τ,t) for any θ are given as(48)MAnlθ=EeθAnlτ,t,MSnlθ=EeθSnlτ,t.Theorem 5.
Given a single service node where its arrival flows are described byAnl(τ,t), and its available service is given by a dynamic server Snl(τ,t). Fix ϵ>0 and define(49)Mnθ,s,t=MAnlθ,s,t∘M-Snlθ,s,tThen we have the following performance bounds.Output burstiness:PrDnlτ,t>Dετ,t≤ε, where(50)Dετ,t=infθ∈0,∞⁡1θln⁡Mnθ,s,t-ln⁡εBacklog: Prbnlt>bεt≤ε, where(51)bεt=infθ∈0,∞⁡1θln⁡Mnθ,t,t-ln⁡εDelay: Prdnlt>dεt≤ε, where dεt is the smallest number satisfying(52)infθ∈0,∞⁡Mnθ,t+dεt,t≤εProof.
From the inequationDnls,t≤Anl⊘Snl(s,t), we have (53)Pr⁡Dnls,t>Dε≤Pr⁡Anl⊘Snls,t>DεComparing with inequation (14), it is easy to get the results for output burstiness.
Frombnlt≤Anl⊘Snl(t,t), we have (54)Pr⁡bt,t>bε≤Pr⁡Anl⊘Snlt,t>bεComparing with (14) and replacings witht, we have the output backlog bounds.
From the definition ofdε, it is easy to get(55)1θMnθ,t+dε,t-ln⁡ε≤0.As dnl=inf⁡τ≥0:Anl⊘Snlt,t+τ<0, we can get(56)Pr⁡dnlt>dε≤Mθ,t+dε,t,and the result is obvious.


## 4.2. Multiple Hop Flow Model

For multiple hop situations, the stability condition is expressed as(57)limt→∞⁡∑l∈FlownAnl0,tt<limt→∞⁡Sn0,tt.For flowl, we want to compute its performance over multiple hops. At first, we introduce the following theorem.Theorem 6.
Considering a flowl passes a work-conserving server with the service process Sn(s,t) that Sn(s,t) is nonnegative, increasing int, and Sns,s=0. Then the service provided to flowl is given by(58)Snls,t=max⁡0,Sns,t-∑i∈PrevFlown,lAnms,t.HereFlow(n) gives all the data follows that pass the noden, and Prev(ϑ,l) is the flows in the flow set ϑ with the priority higher than that ofl.The MGF ofSnl(s,t) is bounded for θ≥0,t≥s≥0 by(59)M-Snlθ,s,t≤min⁡1,M-Snθ,s,t·∏m∈PrevFlown,lMAnmθ,s,tProof.
For a fixed priority flowl, we only consider the flows with a priority higher than that ofl. The result (58) is proved by applying (17) several times. The MGF of Snl(s,t) can be deduced directly from (18).Theorem 7.
Consider a flowl charactering by Al(s,t) passes nodes i1,i2,⋯,in in series. There exists an equivalent, single dynamic server Sl(s,t) for 0≤s≤t where(60)Sls,t=Si1l⊗Si2l⊗⋯⊗Sinls,tAnd the MGF of Sl(s,t), the equivalent, single dynamic server is upper bounded for 0≤s≤t according to(61)M-slθ,s,t≤M-si1lθ∗M-si2lθ∗⋯∗M-sinlθs,tProof.
Equation (60) can be deduced directly by applying (19) several times and the MGF upper bound from (20).Corollary 8.
Consider a flowl   Al(s,t) passes a cascade of N servers i1,i2,⋯,in; it has(62)M-Al⊘Sl≤MAlθ∘M-Si1lθ∗M-Si2lθ∗⋯∗M-Sinlθs,t=MAlθ∘M-Si1lθ∘M-Si2lθ∘⋯∘M-Sinlθs,tProof.
From the definition of operators∗ and ∘, it has(63)x∘y∗zs,t=∑τ1=0∞xt-s+τ1∑τ2τ1yτ1-τ2zτ2=∑u1=0∞∑u2=0∞xt-s+u1+u2yu1zu2=x∘y∘zs,t.The result is gotten directly from (61) by applying above equation several times.The computation network performance analysis process for a WMNCS is given as Algorithm1.Algorithm 1:Compute end-to-end delays and MGFs.
( 1 )   begin ( 2 )  Initialize the left-over service Sileft(s,t) and its MGF upper bound of each node n as Si(s,t) and MSi(θ); ( 3 )  for each  lin the flow set in descending priority order ( 4 )   for eachi inNode(l) in the order of passing ( 5 )    Get the cumulative arrival Ail(s,t) and its MGF upper bound; ( 6 )    Compute the cumulative departure Dil(s,t) and MGF upper bound as next node cumulative arrival; ( 7 )    Update the left-over service Sileft(s,t) and its MGF upper bound at node i; ( 8 )   end for ( 9 )   Compute the tandem service for flow l and its MGF upper bound; ( 10 ) Compute the tandem delay and its MGF upper bound; ( 11 )  end for ( 12 )  end beginHere,Node(l) denotes the node set that flowl passes, and Sileft denotes the leave over service of nodei.


## 5. Physical Layer Model

For wireless networks, the condition of wireless channel and channel coding method will impact its network performance. At the physical layer, we are interested in how the channel fading and channel coding affect the network delays.


### 5.1. Channel Fading

As the radio-wave propagation through wireless channel is a complicated phenomenon characterized by various effects, such as multipath and shadowing. Wireless channel is characterized by the rapid variation of channel quality. Channel fading is used to refer to the deviation in the attenuation experienced by the transmitted signal when traversing a wireless channel.To get a precise mathematical description of channel fading is impossible, and if can be obtained, it will also be too complex for tractable communications systems analyses. A range of relatively simple and accurate statistical models for fading channel which depend on the particular propagation environment and the underlying communication scenario has been proposed. A lot of models are proposed to describe the gain of fading channel depending on the type of fading (slow or fast) and the environment. For industrial application environment, we focus on multipath fading which is due to the combined effects of randomly delayed, reflected, scattered, and diffracted signal components. Multipath fading is relatively fast and is therefore responsible for the short-term signal variations. Table1 lists some most commonly used multipath fading models [7].Table 1
Multipath fading models.
Type of Fading PDF, p γ γ Rayleigh 1 γ - exp ⁡ - γ γ - Nakagami-m m m γ m - 1 γ - m Γ m exp ⁡ - m γ γ - Nakagami-q ( 1 + q 2 ) 2 q γ - exp ⁡ - ( 1 + q 2 ) 2 γ 4 q 2 γ - × I 0 - ( 1 - q 4 ) γ 4 q 2 γ - Log-normal shadowing 4.34 2 π σ γ exp ⁡ - 10 log 10 ⁡ γ - μ 2 2 σ 2 Composite gamma/log-normal ∫ 0 ∞ m m γ m - 1 w m Γ m exp ⁡ - m γ γ - × ξ 2 π σ w exp ⁡ - 10 log 10 ⁡ w - μ 2 2 σ 2 d wHere,pγγ is the fading PDF function, γ is the instantaneous signal-to-noise power ratio (SNR) per symbol, and γ- is the average signal-to-noise power ratio (SNR) per symbol.


### 5.2. Channel Coding

When error-correction coding is applied to the transmitted modulation and decisions are made based on an observation of the received signal much longer than a single symbol interval, it becomes necessary to consider the variation of the fading channel from symbol interval to symbol interval.Papers [7, 11, 12] study the performance of coded communications over fading channel; the results are given as upper bounds on the average bit error probability. These models can also be integrated into our framework. As the length limit, we do not discuss these models in the work. For more detailed descriptions about coding effects for fading channel, please refer to [7, 11, 12].


## 5.1. Channel Fading

As the radio-wave propagation through wireless channel is a complicated phenomenon characterized by various effects, such as multipath and shadowing. Wireless channel is characterized by the rapid variation of channel quality. Channel fading is used to refer to the deviation in the attenuation experienced by the transmitted signal when traversing a wireless channel.To get a precise mathematical description of channel fading is impossible, and if can be obtained, it will also be too complex for tractable communications systems analyses. A range of relatively simple and accurate statistical models for fading channel which depend on the particular propagation environment and the underlying communication scenario has been proposed. A lot of models are proposed to describe the gain of fading channel depending on the type of fading (slow or fast) and the environment. For industrial application environment, we focus on multipath fading which is due to the combined effects of randomly delayed, reflected, scattered, and diffracted signal components. Multipath fading is relatively fast and is therefore responsible for the short-term signal variations. Table1 lists some most commonly used multipath fading models [7].Table 1
Multipath fading models.
Type of Fading PDF, p γ γ Rayleigh 1 γ - exp ⁡ - γ γ - Nakagami-m m m γ m - 1 γ - m Γ m exp ⁡ - m γ γ - Nakagami-q ( 1 + q 2 ) 2 q γ - exp ⁡ - ( 1 + q 2 ) 2 γ 4 q 2 γ - × I 0 - ( 1 - q 4 ) γ 4 q 2 γ - Log-normal shadowing 4.34 2 π σ γ exp ⁡ - 10 log 10 ⁡ γ - μ 2 2 σ 2 Composite gamma/log-normal ∫ 0 ∞ m m γ m - 1 w m Γ m exp ⁡ - m γ γ - × ξ 2 π σ w exp ⁡ - 10 log 10 ⁡ w - μ 2 2 σ 2 d wHere,pγγ is the fading PDF function, γ is the instantaneous signal-to-noise power ratio (SNR) per symbol, and γ- is the average signal-to-noise power ratio (SNR) per symbol.


## 5.2. Channel Coding

When error-correction coding is applied to the transmitted modulation and decisions are made based on an observation of the received signal much longer than a single symbol interval, it becomes necessary to consider the variation of the fading channel from symbol interval to symbol interval.Papers [7, 11, 12] study the performance of coded communications over fading channel; the results are given as upper bounds on the average bit error probability. These models can also be integrated into our framework. As the length limit, we do not discuss these models in the work. For more detailed descriptions about coding effects for fading channel, please refer to [7, 11, 12].


## 6. Integrate Three Layer Models

In this section, we show how to overcome the gaps between above three-layer models with upper bound approach and stochastic orders. By changing and combining the models at different layers, our framework can adapt various modeling requirements.


### 6.1. Modeling the Input Streams

In WMNCSs, typically, sensors would take measurements in a periodic fashion. However, when data is sent through a shared network, the sensors are unable to send data if the network is busy and may need to wait until data can be transmitted, so this introduces jitters for input stream as is shown in Figure8.Figure 8
Input event curves for sensor node.For these properties, we use the constrained arrival processσθ,ρ(θ) to model the traffic caused by sensor sampling and controller commanding. The MGFs of the cumulative arrival processes of sensing and command data flows in the bit domain are bounded by(64)1θln⁡MAθ,τ,t≤σθ+ρθt-τ.A wide variety of traffic, including periodic sources, fluid sources, on-off sources, and regulated sources can be described using the constrained arrival process σθ,ρ(θ) [13–16].


### 6.2. Modeling the Server Ability

In this subsection, we show how the fading channel affects the performance of wireless networks and how to use these physical layer models to refine the bit layer model.In bit layer, we model the wireless network as tandem queues of servers. In the real world, these servers are fading channels with variable capacity. Here we ignore the impact of coding and assume the transmission rates over the fading channel are equal to their information-theoretic capacity limit. So the instantaneous, information-theoretic channel capacity of a fading channel is given by [17](65)Cγ=cwlog2⁡gγ.Here cw is a constant of bandwidth and g(γ) is used to describe the feature of fading channel.For the fast-fading channel, where the latency requirement is greater than the coherence time and the codeword length spans many coherence periods, one can average over many independent channels fades by coding over a large number of coherence time intervals. Thus, it is possible to achieve a reliable rate of communication of(66)Elog2⁡1+h2γbits/s/Hz.And it is meaningful to speak of this value as the capacity of the fast-fading channel. Here h2is random channel gain [17].To calculate the service model for a fast-fading fading channel, we assume that the channel state is sampled by the receiver at a fixed time intervalT∆ and use this sampled data γi to present the instantaneous SNR in thei-th time interval. It is obvious that γi is a random variable that has the probability distribution of the underlying fading model. When time interval T∆is longer than the channel coherence, γi can be treated as independent random variables with identical distribution. Under the above assumptions and setting T∆=1, the corresponding service process at bit layer is given by(67)Snτ,t=cγ∑i=τt-1ln⁡gγi,wherecγ=cw/ln⁡2.And it is easy to get(68)MSnθ,s,t=cγMln⁡gγθt-s.Theorem 9.
Consider a flowl passes a cascade of Nindependent and identically distributed fading channels, where each channel is described byg(γ). The MGF of this concatenation server is up bounded by(69)MSlθ,s,t≤cγN.N-1+t-st-s·Mln⁡gγθt-sProof.
From Theorem7, it has(70)M-Slθ,s,t≤M-Si1lθ∗M-Si2lθ∗⋯∗M-SiNlMθs,t=∑s=τ0≤τ1≤⋯≤τN=t∏i=1NMsilθ,τ1,τ2=∑s=τ0≤τ1≤⋯≤τN=t∏i=1Ncγ·Mln⁡gγθτi-τi-1=cγN∑s=τ0≤τ1≤⋯≤τN=tMln⁡gγθt-s=cγN·N-1+t-st-s·Mln⁡gγθt-s.The theorem is proven.Corollary 10.
The end-to-end MGF bounds for a cascade of N fading channelsSl(s,t) with σθ,ρ(θ) constrained arrivals process Al(s,t) are given by(71)MAl⊘Slθ,t+dε,t≤cγNeθρθt-s+σθ1-cγeθρθMln⁡gγθN.Proof.
From Theorem5, it has(72)MAl⊘Slθ,s,t≤MAlθ,s,t∘M-Slθ,s,t=∑τ=0seθρθt-τ+σθ·cγNN-1+s-τs-τ·Mln⁡gγθs-τ=cγN∑k=0seθ-ρθt-s+σθ·N-1+kk·eθρθMln⁡gγθk≤cγNeθ-ρθt-s+σθ∑k=0∞N-1+kkeθρθMln⁡gγθk=cγNeθρθt-s+σθ1-eθρθMln⁡gγθN.In the last step, we have used the combination equation(73)∑k=0∞N-1+kk·xk=11-xN.Corollary 11.
To compute the delay bounds, it has the following inequation:(74)MAl⊘Slθ,t+dε,t≤eθρθt-s+σθ∙min⁡1,cγeθρθMln⁡gγθdεdεN-11-cγeθρθMln⁡gγθNProof.
From Theorem7 and the combination inequation(75)N-1+kk≤dεN-1·N-1+k-dεdεit has(76)MAl⊘Slθ,t+dε,t≤MAlθ,t+dε,t∘M-Slθ,t+dε,t≤∑τ=dεteθρθt-τ+σθ·cγNN-1+t+dε-τt+dε-τ·Mln⁡gγθt+dε-τ=∑k=dεtcγNeθ-ρθdε+σθ·N-1+kk·eθρθMln⁡gγθk≤cγNeθ-ρθdε+σθ∑k=dεtN-1+k-dεk-dεeθρθMln⁡gγθkdεN-1≤cγNeθ-ρθdε+σθ∑m=0∞N-1+mmeθρθMln⁡gγθm+dεdεN-1≤cγNeθ-ρθdε+σθ·min⁡1,cγeθρθMln⁡gγθdεdεN-11-eθρθMln⁡gγθN.We end the proof.


### 6.3. Analyzing System Stability

In this section, we first introduce the definition of stochastic order of SHSs and then give an upper bounded method to analyze the stability of SHS.Stochastic Order. Given two random variablesX andY, we say ifX is stochastically smaller thanY, Pr⁡[X>x]≤Pr⁡[Y>x], for allx, or in other words F-X(x)≤F-Y(x), for all x, is written as (77)X≤stY.Similarly, for any two stochastic processes, it has(78)Xt≤stYtifF-Xx,t≤F-Yx,t.For stochastic order, it has the following properties.( 1 ) If X≤stY, then fX≤stf(Y) for any increasing functionf.( 2 ) Let X1,X2,⋯,Xn and Y1,Y2,⋯,Yn be independent. IF Xi≤stYi, then for any wide-sense increasing function ΦZ1,Z2,⋯,Zn on Zi (Φ⋯,Zi,⋯≤Φ⋯,Zi′,⋯ when Zi<Zi′ ); then it has(79)ΦX1,X2,⋯,Xn≤stΦY1,Y2,⋯,YnStochastic Order of SHSs. For an SHS Ω1, we define an SHS Ω2 and all the other ingredients of Ω2 remain the same with Ω1, except for reset time distributions. Then we define that Ω1 is stochastically smaller than Ω2, denoted by(80)Ω1≤stΩ2,if all the reset time distributions of Ω1 (denoted by μq,l1(x)) are stochastically smaller than the corresponding reset time distributions of Ω2 (denoted by μq,l2(x)) that (81)μq,l1x≤stμq,l2x,∀q∈Q,l∈L.Theorem 12.
For two SHSsΩ1 and Ω2that Ω1≤stΩ2, and Ω2 is MES, then Ω1 is also MES if it meets the following requirements: ϕ-l(q,x,t)-x≤0 for all q∈Q, x∈Rn and l∈L. Here q,x=ϕlq-,x-,t and ϕ-lis projection of ϕl to Rn (ϕ-lq,x-,t=x).Proof.
Let the functionψ(q,x,t): Q×Rn×[0,∞)→R be the continuously differentiable with respect tox andt. From its Itô equation (37), dividing the result bydt and take expectations on Ω1 and Ω2, it has(82)dEiψq,x,tdt=Ei∂ψq,x,t∂xfq,x,t+∂ψq,x,t∂t+∑l=1nlψφlq,x,t,t-ψq,x,tdNliqdt.where i=1,2 represent the operations on Ω1 and Ω2 respectively. Here Nli(q) is a stochastic transition counter for jump l at modeq. As μq,l1(x)≤stμq,l2(x), it obviously has(83)dNl1qdt≤stdNl2qdt.With ψq,x,t=x=xT∙x and the assumption ϕ-lq,x,t-x≤0, it has(84)dE1xT·xdt≤dE2xT·xdtfor all q∈Q, x∈Rn and l∈L. As Ω2 is MES, it is obvious that(85)E2xtTxt≤ExtTxt≤ce-αtx0Tx0,∀t≥0.From definition (35), it can be concluded that Ω1is MES.


## 6.1. Modeling the Input Streams

In WMNCSs, typically, sensors would take measurements in a periodic fashion. However, when data is sent through a shared network, the sensors are unable to send data if the network is busy and may need to wait until data can be transmitted, so this introduces jitters for input stream as is shown in Figure8.Figure 8
Input event curves for sensor node.For these properties, we use the constrained arrival processσθ,ρ(θ) to model the traffic caused by sensor sampling and controller commanding. The MGFs of the cumulative arrival processes of sensing and command data flows in the bit domain are bounded by(64)1θln⁡MAθ,τ,t≤σθ+ρθt-τ.A wide variety of traffic, including periodic sources, fluid sources, on-off sources, and regulated sources can be described using the constrained arrival process σθ,ρ(θ) [13–16].


## 6.2. Modeling the Server Ability

In this subsection, we show how the fading channel affects the performance of wireless networks and how to use these physical layer models to refine the bit layer model.In bit layer, we model the wireless network as tandem queues of servers. In the real world, these servers are fading channels with variable capacity. Here we ignore the impact of coding and assume the transmission rates over the fading channel are equal to their information-theoretic capacity limit. So the instantaneous, information-theoretic channel capacity of a fading channel is given by [17](65)Cγ=cwlog2⁡gγ.Here cw is a constant of bandwidth and g(γ) is used to describe the feature of fading channel.For the fast-fading channel, where the latency requirement is greater than the coherence time and the codeword length spans many coherence periods, one can average over many independent channels fades by coding over a large number of coherence time intervals. Thus, it is possible to achieve a reliable rate of communication of(66)Elog2⁡1+h2γbits/s/Hz.And it is meaningful to speak of this value as the capacity of the fast-fading channel. Here h2is random channel gain [17].To calculate the service model for a fast-fading fading channel, we assume that the channel state is sampled by the receiver at a fixed time intervalT∆ and use this sampled data γi to present the instantaneous SNR in thei-th time interval. It is obvious that γi is a random variable that has the probability distribution of the underlying fading model. When time interval T∆is longer than the channel coherence, γi can be treated as independent random variables with identical distribution. Under the above assumptions and setting T∆=1, the corresponding service process at bit layer is given by(67)Snτ,t=cγ∑i=τt-1ln⁡gγi,wherecγ=cw/ln⁡2.And it is easy to get(68)MSnθ,s,t=cγMln⁡gγθt-s.Theorem 9.
Consider a flowl passes a cascade of Nindependent and identically distributed fading channels, where each channel is described byg(γ). The MGF of this concatenation server is up bounded by(69)MSlθ,s,t≤cγN.N-1+t-st-s·Mln⁡gγθt-sProof.
From Theorem7, it has(70)M-Slθ,s,t≤M-Si1lθ∗M-Si2lθ∗⋯∗M-SiNlMθs,t=∑s=τ0≤τ1≤⋯≤τN=t∏i=1NMsilθ,τ1,τ2=∑s=τ0≤τ1≤⋯≤τN=t∏i=1Ncγ·Mln⁡gγθτi-τi-1=cγN∑s=τ0≤τ1≤⋯≤τN=tMln⁡gγθt-s=cγN·N-1+t-st-s·Mln⁡gγθt-s.The theorem is proven.Corollary 10.
The end-to-end MGF bounds for a cascade of N fading channelsSl(s,t) with σθ,ρ(θ) constrained arrivals process Al(s,t) are given by(71)MAl⊘Slθ,t+dε,t≤cγNeθρθt-s+σθ1-cγeθρθMln⁡gγθN.Proof.
From Theorem5, it has(72)MAl⊘Slθ,s,t≤MAlθ,s,t∘M-Slθ,s,t=∑τ=0seθρθt-τ+σθ·cγNN-1+s-τs-τ·Mln⁡gγθs-τ=cγN∑k=0seθ-ρθt-s+σθ·N-1+kk·eθρθMln⁡gγθk≤cγNeθ-ρθt-s+σθ∑k=0∞N-1+kkeθρθMln⁡gγθk=cγNeθρθt-s+σθ1-eθρθMln⁡gγθN.In the last step, we have used the combination equation(73)∑k=0∞N-1+kk·xk=11-xN.Corollary 11.
To compute the delay bounds, it has the following inequation:(74)MAl⊘Slθ,t+dε,t≤eθρθt-s+σθ∙min⁡1,cγeθρθMln⁡gγθdεdεN-11-cγeθρθMln⁡gγθNProof.
From Theorem7 and the combination inequation(75)N-1+kk≤dεN-1·N-1+k-dεdεit has(76)MAl⊘Slθ,t+dε,t≤MAlθ,t+dε,t∘M-Slθ,t+dε,t≤∑τ=dεteθρθt-τ+σθ·cγNN-1+t+dε-τt+dε-τ·Mln⁡gγθt+dε-τ=∑k=dεtcγNeθ-ρθdε+σθ·N-1+kk·eθρθMln⁡gγθk≤cγNeθ-ρθdε+σθ∑k=dεtN-1+k-dεk-dεeθρθMln⁡gγθkdεN-1≤cγNeθ-ρθdε+σθ∑m=0∞N-1+mmeθρθMln⁡gγθm+dεdεN-1≤cγNeθ-ρθdε+σθ·min⁡1,cγeθρθMln⁡gγθdεdεN-11-eθρθMln⁡gγθN.We end the proof.


## 6.3. Analyzing System Stability

In this section, we first introduce the definition of stochastic order of SHSs and then give an upper bounded method to analyze the stability of SHS.Stochastic Order. Given two random variablesX andY, we say ifX is stochastically smaller thanY, Pr⁡[X>x]≤Pr⁡[Y>x], for allx, or in other words F-X(x)≤F-Y(x), for all x, is written as (77)X≤stY.Similarly, for any two stochastic processes, it has(78)Xt≤stYtifF-Xx,t≤F-Yx,t.For stochastic order, it has the following properties.( 1 ) If X≤stY, then fX≤stf(Y) for any increasing functionf.( 2 ) Let X1,X2,⋯,Xn and Y1,Y2,⋯,Yn be independent. IF Xi≤stYi, then for any wide-sense increasing function ΦZ1,Z2,⋯,Zn on Zi (Φ⋯,Zi,⋯≤Φ⋯,Zi′,⋯ when Zi<Zi′ ); then it has(79)ΦX1,X2,⋯,Xn≤stΦY1,Y2,⋯,YnStochastic Order of SHSs. For an SHS Ω1, we define an SHS Ω2 and all the other ingredients of Ω2 remain the same with Ω1, except for reset time distributions. Then we define that Ω1 is stochastically smaller than Ω2, denoted by(80)Ω1≤stΩ2,if all the reset time distributions of Ω1 (denoted by μq,l1(x)) are stochastically smaller than the corresponding reset time distributions of Ω2 (denoted by μq,l2(x)) that (81)μq,l1x≤stμq,l2x,∀q∈Q,l∈L.Theorem 12.
For two SHSsΩ1 and Ω2that Ω1≤stΩ2, and Ω2 is MES, then Ω1 is also MES if it meets the following requirements: ϕ-l(q,x,t)-x≤0 for all q∈Q, x∈Rn and l∈L. Here q,x=ϕlq-,x-,t and ϕ-lis projection of ϕl to Rn (ϕ-lq,x-,t=x).Proof.
Let the functionψ(q,x,t): Q×Rn×[0,∞)→R be the continuously differentiable with respect tox andt. From its Itô equation (37), dividing the result bydt and take expectations on Ω1 and Ω2, it has(82)dEiψq,x,tdt=Ei∂ψq,x,t∂xfq,x,t+∂ψq,x,t∂t+∑l=1nlψφlq,x,t,t-ψq,x,tdNliqdt.where i=1,2 represent the operations on Ω1 and Ω2 respectively. Here Nli(q) is a stochastic transition counter for jump l at modeq. As μq,l1(x)≤stμq,l2(x), it obviously has(83)dNl1qdt≤stdNl2qdt.With ψq,x,t=x=xT∙x and the assumption ϕ-lq,x,t-x≤0, it has(84)dE1xT·xdt≤dE2xT·xdtfor all q∈Q, x∈Rn and l∈L. As Ω2 is MES, it is obvious that(85)E2xtTxt≤ExtTxt≤ce-αtx0Tx0,∀t≥0.From definition (35), it can be concluded that Ω1is MES.


## 7. Case Study

In this section, we show an example to model and analyze WMNCSs using our framework. As is shown in Figure9, a WMNCS that controls a batch reactor is deployed in a chemical factory, where one sensor node and one actor node are connected to the remote controller through relay nodes. The channel fading effect is taken into account as the complexity of the environment.Figure 9
A batch reactor control system based on a wireless mesh network.


### 7.1. System Models

A linear state-space model of this process is introduced by Taylor expansion near the working point. The plant of this reactor is described by(86)x˙p=0.106-0.006-0.045-0.330.52-0.4300.520.0820.3290.0030.329-0.050.450.103-0.16xp+000.4700.08-0.30.080u^yp=101-10100xpThe open loop system is unstable, and we use a feedback controller to stabilize the system; i.e.,(87)x˙c=0110y^p,y^p=-2008xc+0-4100y^pThen the system is stable when the reactor and the controller are directly connected.The system is implemented with wireless mesh network technology. Without considering packet drops, it can be modeled by an SHS with two discrete modesnq=2, as the TTSHS model given in Figure 5.For the wireless mesh network architecture, the sensors and actuators are I/O device that provides or consumes data. The controller is deployed at the gateway, which communicates with the sensor node and actuator node through a backbone router and a cascade of route nodes respectively. For the simplicity of analysis, we assume that the routers for sensor information and actuator command are different, except for the backbone router, which is usually more capable and connected to the gateway directly through the wired network. In this paper, we only consider the delays because wireless communications, time for wired communications, and computation of control task are thought to be negligible. We use two cascades of stochastic network servers to model the wireless communications,Ssc(s,t) for the sensor node to the backbone router, and Sca(s,t) for a backbone router to the actuator node. We use Nsc and Nca to represent the number of servers for Ssc(s,t) and Sca(s,t), respectively.Here, we consider the dynamic SNR server description for a Rayleigh fading channel with fast fading. The fading gain for Rayleigh fading is a random variable with distribution(88)fx=2xe-x2,and the MGF for a single server is given by (89)Msnθ,s,t=cγe1/γ-γ--θΓ1-θ,γ--1t-s.where Γ(x,y) is the upper incomplete Gamma function:(90)Γx,y=∫y∞tx-1e-tdt.Another concern is node conflict, which is due to the half-duplex radio; two transmissions conflict with each other if they share a node (sender or receiver). In a statistical sense, the service is provided by a single server will be half off. We have the following setting:(91)cγ=cw2ln⁡2.MGF bound for a cascade of N Rayleigh fading channels is given by(92)MA⊘SNθ,τ,t≤cγNeθρθt-s+σθ1-eθρθe1/γ-γ--θΓ1-θ,γ--1N.From the relationships at Figure4, we can get the MGF up bounds of the set-time distributions μ1,1(t) and μ2,2(t), which are given as(93)Mμ1,1θ,s,t≤Ts-Tscmin-Tcamin+MAsc⊘SscNscθ,s,t,(94)Mμ2,2θ,s,t=MAca⊘ScaNcaθ,s,t.where Tscmin and Tcamincan be evaluated through the method described in the next subsection.


### 7.2. Numerical Results

The values of all the important parameters for the above models are listed in Table2. Here, the network parameters are set to Cw=20kHZ, γ-=5dB,15dB. For traffic, σθ,ρ(θ) bounded input arrivals are set toρ(θ) = 20 kbps,σ(θ) = 20 kbps for sampling data flows, andσ(θ) = 50 kbps for command data flows. The control loop length is set to Ts=200ms.Table 2
Model parameter values.
c w 20kHZ ρ ( θ ) 10kbps σ ( θ ) 20kbps, 50kbps T s 200ms γ - 5dB,15dBWe next present how the delay violation probability for a given end-to-end delay bound is impacted by the node numbers and the SNR wireless channel. We compare multiple hop networks with 2, 4, and 6 server nodes underCw=20kHZ, SNR γ-=5dB or 15dB, respectively. The results are obtained from Theorem 5 by setting the value of violation probability ε, solving dεon the left side of the equation, and minimizing over θ. We use numerical methods to compute the end to end delay bounds for the complexity of the bound expression [18]. Observing that the end-to-end bounds are scale linearly with the number of nodes. The complexity of computing the end-to-end bounds of multiple nodes is not more difficult than that for a single node.Figures10 and 11 show the results of the set-time violation bounds for μ1,1 and μ2,2 respectively. By getting the end-to-end delay bound for 0.99 violation probability at Figure 10 for a given SNR and dividing it by the number of nodes, we get the minimum delay for single node ΔTγmin, and the parameters Tscmin and Tcamin at (94) under different SNR can be evaluated by multiply Nsc and Nca with ΔTγmin, respectively. The results show that the end-to-end delays are mainly affected by the number of nodes at high SNR, but low SNR will greatly increase the end-to-end delays, as well as affect the aggregation of delay distribution functions. By setting the end-to-end delay bound greater than Ts and getting the corresponding violation probability, it gives us an observation of the violation probability for keeping the system event sequence order. Another observation is that all these bound traces are nearly piecewise linear for Rayleigh fading and can be further upper bounded by a piecewise linear function with two sections, corresponding to an exponential distribution(95)μt=λe-λt,t>a0,t≤a.Here, λ and a can be evaluated through estimating the return model parameters. This can be used to further simplify the stability analysis computation at the next step.Figure 10
Delay bound violation probability versus end-to-end delays forμ1,1.Figure 11
Delay bound violation probability versus end-to-end delays forμ2,2.Then the multiple hop distribution upper bounds are used to investigate the stability of the system under different network topology and channel condition. For control loop length setTs=200ms, and SNR γ-=5dB,15dB. We test the mean exponential stability of the closed loop using the numerical method described in [19] by theorem 1. The results obtained are summarized in Figure 12.Figure 12
MES for different nodes combination of two independent links under SNR = 5dB.In Figure12, the horizontal coordinate is the number of nodes for the link between sensor and controller, and the vertical coordinate is the number of nodes for the link between the controller and actuator. The solid circles present the pass cases under this test, which means the system is guaranteed to be MES for the upper bounded methods we have used. The hollow circles present the failed cases, which mean the system is identified to be MES unstable with our methods, but maybe actually stable. The overall results presented above give us a quantitative inspection of how the channel condition and multiple hops affect the stability of the WMNCSs. Currently, the scalability of this framework is impeded by the used numerical method. In future research, we will explore numerical methods that can support faster, more exact, and various analysis.


## 7.1. System Models

A linear state-space model of this process is introduced by Taylor expansion near the working point. The plant of this reactor is described by(86)x˙p=0.106-0.006-0.045-0.330.52-0.4300.520.0820.3290.0030.329-0.050.450.103-0.16xp+000.4700.08-0.30.080u^yp=101-10100xpThe open loop system is unstable, and we use a feedback controller to stabilize the system; i.e.,(87)x˙c=0110y^p,y^p=-2008xc+0-4100y^pThen the system is stable when the reactor and the controller are directly connected.The system is implemented with wireless mesh network technology. Without considering packet drops, it can be modeled by an SHS with two discrete modesnq=2, as the TTSHS model given in Figure 5.For the wireless mesh network architecture, the sensors and actuators are I/O device that provides or consumes data. The controller is deployed at the gateway, which communicates with the sensor node and actuator node through a backbone router and a cascade of route nodes respectively. For the simplicity of analysis, we assume that the routers for sensor information and actuator command are different, except for the backbone router, which is usually more capable and connected to the gateway directly through the wired network. In this paper, we only consider the delays because wireless communications, time for wired communications, and computation of control task are thought to be negligible. We use two cascades of stochastic network servers to model the wireless communications,Ssc(s,t) for the sensor node to the backbone router, and Sca(s,t) for a backbone router to the actuator node. We use Nsc and Nca to represent the number of servers for Ssc(s,t) and Sca(s,t), respectively.Here, we consider the dynamic SNR server description for a Rayleigh fading channel with fast fading. The fading gain for Rayleigh fading is a random variable with distribution(88)fx=2xe-x2,and the MGF for a single server is given by (89)Msnθ,s,t=cγe1/γ-γ--θΓ1-θ,γ--1t-s.where Γ(x,y) is the upper incomplete Gamma function:(90)Γx,y=∫y∞tx-1e-tdt.Another concern is node conflict, which is due to the half-duplex radio; two transmissions conflict with each other if they share a node (sender or receiver). In a statistical sense, the service is provided by a single server will be half off. We have the following setting:(91)cγ=cw2ln⁡2.MGF bound for a cascade of N Rayleigh fading channels is given by(92)MA⊘SNθ,τ,t≤cγNeθρθt-s+σθ1-eθρθe1/γ-γ--θΓ1-θ,γ--1N.From the relationships at Figure4, we can get the MGF up bounds of the set-time distributions μ1,1(t) and μ2,2(t), which are given as(93)Mμ1,1θ,s,t≤Ts-Tscmin-Tcamin+MAsc⊘SscNscθ,s,t,(94)Mμ2,2θ,s,t=MAca⊘ScaNcaθ,s,t.where Tscmin and Tcamincan be evaluated through the method described in the next subsection.


## 7.2. Numerical Results

The values of all the important parameters for the above models are listed in Table2. Here, the network parameters are set to Cw=20kHZ, γ-=5dB,15dB. For traffic, σθ,ρ(θ) bounded input arrivals are set toρ(θ) = 20 kbps,σ(θ) = 20 kbps for sampling data flows, andσ(θ) = 50 kbps for command data flows. The control loop length is set to Ts=200ms.Table 2
Model parameter values.
c w 20kHZ ρ ( θ ) 10kbps σ ( θ ) 20kbps, 50kbps T s 200ms γ - 5dB,15dBWe next present how the delay violation probability for a given end-to-end delay bound is impacted by the node numbers and the SNR wireless channel. We compare multiple hop networks with 2, 4, and 6 server nodes underCw=20kHZ, SNR γ-=5dB or 15dB, respectively. The results are obtained from Theorem 5 by setting the value of violation probability ε, solving dεon the left side of the equation, and minimizing over θ. We use numerical methods to compute the end to end delay bounds for the complexity of the bound expression [18]. Observing that the end-to-end bounds are scale linearly with the number of nodes. The complexity of computing the end-to-end bounds of multiple nodes is not more difficult than that for a single node.Figures10 and 11 show the results of the set-time violation bounds for μ1,1 and μ2,2 respectively. By getting the end-to-end delay bound for 0.99 violation probability at Figure 10 for a given SNR and dividing it by the number of nodes, we get the minimum delay for single node ΔTγmin, and the parameters Tscmin and Tcamin at (94) under different SNR can be evaluated by multiply Nsc and Nca with ΔTγmin, respectively. The results show that the end-to-end delays are mainly affected by the number of nodes at high SNR, but low SNR will greatly increase the end-to-end delays, as well as affect the aggregation of delay distribution functions. By setting the end-to-end delay bound greater than Ts and getting the corresponding violation probability, it gives us an observation of the violation probability for keeping the system event sequence order. Another observation is that all these bound traces are nearly piecewise linear for Rayleigh fading and can be further upper bounded by a piecewise linear function with two sections, corresponding to an exponential distribution(95)μt=λe-λt,t>a0,t≤a.Here, λ and a can be evaluated through estimating the return model parameters. This can be used to further simplify the stability analysis computation at the next step.Figure 10
Delay bound violation probability versus end-to-end delays forμ1,1.Figure 11
Delay bound violation probability versus end-to-end delays forμ2,2.Then the multiple hop distribution upper bounds are used to investigate the stability of the system under different network topology and channel condition. For control loop length setTs=200ms, and SNR γ-=5dB,15dB. We test the mean exponential stability of the closed loop using the numerical method described in [19] by theorem 1. The results obtained are summarized in Figure 12.Figure 12
MES for different nodes combination of two independent links under SNR = 5dB.In Figure12, the horizontal coordinate is the number of nodes for the link between sensor and controller, and the vertical coordinate is the number of nodes for the link between the controller and actuator. The solid circles present the pass cases under this test, which means the system is guaranteed to be MES for the upper bounded methods we have used. The hollow circles present the failed cases, which mean the system is identified to be MES unstable with our methods, but maybe actually stable. The overall results presented above give us a quantitative inspection of how the channel condition and multiple hops affect the stability of the WMNCSs. Currently, the scalability of this framework is impeded by the used numerical method. In future research, we will explore numerical methods that can support faster, more exact, and various analysis.


## 8. Related Works

There are lots of previous researches closely related to this work. Paper [5] introduces two different models of SHS and illustrates how to use SHS to models NCSs. For a detailed introduction to SHS, please refer to [20–23]. The SHS models we considered in this paper are closely related to the Piecewise Deterministic Markov Process (PDMPs) introduced by Davis [24]. The stability of PDMPs has been studied intensively in many materials, including [10, 25–29]. Paper [30] gives a survey of stability analysis for SHS.This idea of network calculus was initially introduced by Cruz in the seminal work [31]. Network calculus has developed along two tracks, deterministic and stochastic. Book [32] gives a detailed introduction to network calculus. An excellent book [6] is available for deterministic network calculus. Early representative works of SNCs include [33–37]. Recently, many crucial properties have been proved for SNC [38–42]. Paper [8] establishes a concise, probabilistic network calculus with moment generating functions. Paper [18] proposes an SNR domain approach to analyze the performance of multihop fading channels.An introduction to physical models in digital communication is available [17]. Book [7] gives an MGF upper bounded approach to fading channels over digital communication. A detailed description of Rayleigh fading channels in mobile digital communication systems is given in [43, 44].


## 9. Conclusions

Recent progress at the area of network performance analysis and stochastic control theory enlighten this work. In the paper, we propose a three-layer modeling framework to capture the properties of WMNCSs at different levels with SHS, SNC, and physical layer models. We also bridge the gaps between these methods with an upper bound approach. These contributions result in a framework for modeling and analysis WMNCSs with stochastic methods and provide us a way to get detailed and unified results. To the best of our knowledge, it is the first work that provides such a unified and flexible framework to model and analyze WMNCSs with stochastic methods.There are lots of interesting topics for further work. A topic is to extend the existing framework to support more controller models, network scheduling strategies. Another interesting topic is to set up a standard model library for rapid modeling and analysis. Finally, accompanied with this framework, a set of numerical methods can be set up to support faster, more exact, and various analyses.



---

*Source: 3217625-2019-04-04.xml*