# A Simulation Perspective: Error Analysis in the Distributed Simulation of Continuous System

Authors: Yao-fei Ma; Xiao Song
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

To construct a corresponding distributed system from a continuous system, the most convenient way is to partition the system into parts according to its topology and deploy the parts on separated nodes directly. However, system error will be introduced during this process because the computing pattern is changed from the sequential to the parallel. In this paper, the mathematical expression of the introduced error is studied. A theorem is proposed to prove that a distributed system preserving the stability property of the continuous system can be found if the system error is limited to be small enough. Then, the compositions of the system error are analyzed one by one and the complete expression is deduced, where the advancing stepT in distributed environment is one of the key factors associated. At last, the general steps to determine the step T are given. The significance of this study lies in the fact that the maximum T can be calculated without exceeding the expected error threshold, and a larger T can reduce the simulation cost effectively without causing too much performance degradation compared to the original continuous system.


---

## Body

## 1. Introduction

Engineering problems involve dealing with physical systems. Since most physical laws are described by differential equations, simulation in engineering is in fact related to numerical resolution of differential equations. This is called continuous system simulation and even the components (i.e., the subsystems or submodels) of it are generally described in time discretization manner [1].The large engineering system has sustained requirement on distributed simulation for two reasons. First, the continuous growth on system scale and complexity results in intense demand on computing capacity, which can be mitigated in distributed environment by scattering the computing load to networked nodes. Second, the system itself is sometimes geographically fragmented; thus, the distributed structure is needed to be consistent with its topology. Except those constructed in distributed manner from scratch, there are still many classic continuous systems that were constructed in the nondistributed way, in that the system was designed and tested without considering the possible scenario of distributed simulation. The demand to transform such system into the distributed one emerges.The classic continuous system is characterized by computing its components in pipeline, where a computing sequence is set up and the components are computed one by one within a step; the sequence is determined by the component’s input/output characteristics and the data dependence. By contrast, the computation on each component in the distributed environment will start simultaneously and advance in a parallel manner; the updating data is exchanged periodically through the network for synchronization.The difference between two computing patterns can be formulized as follows. Suppose a system containingn components: M={m1,m2,…,mn}. A single component can be represented as mi={xi,yi,fi}, where xi is the input, yi is the output, and fi is the model function. We have yi=fi(xi). Assuming a computing sequence Q={q1,…,qi,…,qn} has been determined in the nondistributed environment, the pipeline computing within one step can be described as(1)xq1k×h=Ik×h,yq1k×h=fq1xq1k×h,xqik×h=yqi-1k×h,yqik×h=fqixqik×h=fqiyqi-1k×h,i=2,…,n,where k=0,1,2,…,h is the step size and the subscript qi refers to the index in Q. In the pipeline computing, the first component is assumed to get the input from a signal source I, and the others get inputs following the data dependence among themselves.By contrast, the one-step computing in distributed environment can be described as(2)x1k×T=Ik×T,y1k×T=fix1k×T,xik×T=yi-1k-1×T,yik×T=fixik×T=fiyi-1k-1×T,i=2,…,n,where T is the advancing step of the distributed system. It should be noted that the computing sequence is not held any more in this case.There are two differences comparing these two computing processes. First, the computing is synchronous in the nondistributed environment. For each component, the input is updated to time stepk firstly, and then the output at k is computed (the data dependence among components is assumed to be consistent with the computing sequence). On the other hand, the computation is parallel in the distributed environment and each component starts simultaneously. The input of a component will still hold the value of step (k-1) when the output at k needs to be computed. Second, the advancing steps of two environments, h and T, may be different. In the nondistributed environment, h is determined by the equation solver; either fixed or variable step is employed. However, the determination of T needs to consider the bandwidth of the underlying network, signal frequency, human perception limitation, and so forth. It is often the case where the whole system advances with a fixed T and T≥h.A formal approach capable of partitioning a system into parallel parts and properly handling the above two issues is the key to build a distributed simulation from a nondistributed one. This problem, denoted as thePartitioning Problem, was early studied in the simulation of complex mechanical systems [2, 3] and then became a research focus in the field of discrete event simulation [2, 4, 5] and the decentralized, large-scale control system [6, 7]. However, the proposed approaches all involved extra efforts to reformulize or reengineer the simulated system and did not consider the possible difference on h and T either. These drawbacks set up barriers to apply these approaches broadly.In [4], anovel partitioning approach is proposed using the system’s topology. This approach is easy to perform since most classic simulations have block-diagram structures, but it is “lossy” because the disorder of the component’s input/output data cannot be avoided, thus leading to the error on state trajectory of the resulting distributed system. To reduce the error, a portioning rule was proposed to eliminate the possible cumulative delays caused by improper partitioning, which contribute a lot to the overall error.This paper will study the mathematical expression of this error. The correlation between the error and the advancing stepT in distributed environment is revealed. Then, a maximum T can be calculated according to the error threshold specified by the system engineer. The significance to find the maximum T lies in the fact that a larger T will improve the parallelism of the partitioned system by reducing the synchronization frequency between nodes, without causing too much performance degradation at the same time.


## 2. Literature Review

ThePartitioning Problem was early studied in the coupled problem [2, 5] of complex mechanical system simulation, where the fluids, thermal, control, and structure subsystems interacted with each other and formed multi-physical-field system [8]. The partitioning is applied to decompose such system into partitions with physical or computational considerations. The resulting system could separately advance in time over each partition, thus gained high simulation efficiency. However, these studies were focused on the finite element analysis, not for general purpose simulation.ThePartitioning Problem is also critical to build the decentralized, large-scale control system [6, 7]. A graph-based approach was employed in that the system states were connected as vertexes and the coupling strength (measured by weight factor) between any two of them was evaluated. Small weight connections were more likely to form the partitioning edge by which the system was split into relative independent parts. Another methodwas to use the delay differential equations [9, 10] to describe the dynamic of the distributed system. Networked control scheme [11] was constructed to obtain the convergence and stability in control of the system. However, this approach focuses on the influence of signal delay associated with network latency, rather than the errors caused by the change of computing pattern. In our opinion, such errors will still exist even though the underlying network is perfect and has no delay.In the Modeling and Simulation (M&S) domain, DEVS (Discrete Event System Specification) theory casts light into the solving ofPartitioning Problem. The basic idea was to transform a continuous system into the DEVS form by quantizing the system states. A DEVS system is comprised of a set of connected components (called “atomic model” or “coupled model”). The output of each component will hold unchanged until the values of the states it maintained exceed some predefined quantized level. As a result, the time driving continuous system is transformed into an event driving system; thus, it can be decoupled and partitioned easily [12–14]. The transformed system is suitable for asynchronous, distributed environment in nature; however, the “illegitimacy” phenomenon, where the states may transit for unlimited times within a limited period, often causes the simulation to fail to converge correctly.As an extension of DEVS, the QSS approach [15–17] was proposed to resolve the “illegitimacy” problem by using a special “hysteretic quantization” method, where the output of each component (or model, as called in DEVS system) will not transit to the next quantization level during its descending period until the change has exceeded certain threshold. QSS provides a general approach for the Partitioning Problem since each QSS model interacts with discrete, states updating events. However, the QSS model needs to be particularly designed to allow states transiting between quantization levels. The changes of states are triggered by unpredicted threshold crossing events in QSS. In other words, the advancing “step” of QSS system is unpredicted and time-varying. This characteristic makes QSS unsuitable for applications such as the Hardware-In-the-Loop (HIL) or Man-In-the-Loop (MIL) simulation, where the system has to advance with some fixed step T to align the hardware’s working frequency or to take care of the needs of human perception.The approaches mentioned above all involve extra efforts to reformulize or reengineer the simulated system. For example, in the graph-based approach, the detailed analysis of the coupling strength between system states requests the engineer to have full knowledge of system dynamics. The QSS approach also needs the system components to be modified to follow the DEVS formulation. Additionally, the possible difference betweenh and T was not considered either. These drawbacks set up barriers to apply these approaches broadly.By contrast, the partitioning approach proposed in [4] took the advantage of the system’s structure characteristics; no extra work is needed except for decomposing the system according to the data transferring route. However, the incurred error needs to be further studied considering the possible correlation with T.The content of this paper is organized as follows: in Section3, the previous work is briefly reviewed. In Section 4, the mathematical expression of the errors is given with theorems. In Section 5, a series of steps are concluded to determine the maximum T.


## 3. Problem Description

Normally, a continuous system can be represented by a block diagram as Figure1 shows [4]. This assembling approach based on basic components is commonly seen in the construction of complex system.Figure 1
An inverted pendulum system. The states are “cart position” and “pendulum angle.” A computing sequence is maintained (indicated by red numbers) in the nondistributed environment.The components within this system can be classified into two categories: the Direct-Feed-Through (DFT) component and the Non-Direct-Feed-Through (NDFT) component. The output of DFT component is directly associated with its input; that is, the output at timet is determined by the input at the same time. On the other hand, the output of NDFT component is only determined by its current state and has nothing to do with the current input. It implies that, when computing a DFT component, the components it depends on should be computed firstly to maintain the correct data dependence. The NDFT component has no such requirement since its current output does not rely on the current input.With parallel computing in the distributed environment, more deductions can be deduced from the foregoing. First, if part of the system deployed on a separate node formed a NDFT component (both DFT and NDFT component can have recursive structure), its output would be delayed forT when updating to other nodes. Second, if two or more DFT components were cascaded, they should not be divided into different nodes; otherwise, the output delay would accumulate from the first DFT component. For example, if the components (5), (7), and (9) in Figure 1 were deployed to different nodes, the output of component (9) will be 3T delayed compared to the original system, component (7) is 2T delayed, and component (5) is T delayed.These outcomes have been observed in the experiments of [4], and a partitioning rule was proposed to reduce the accumulated delays. The rule is simple: the distributed DFT component should be avoided; each partitioned part should be ensured to form a NDFT component. The system performance was improved by applying this rule, as shown in Figure 2.Figure 2
(a) The comparison of the trajectories (“cart position”). (a) The distributed system became unstable whenT=2h without applying the partitioning rule. (b) By applying this rule, the distributed system was stable even when T=100h.
(a)
(b)Although the partitioning rule made the distributed system more close to the original continuous system, there is still one step (T) delay among each node. This delay cannot be eliminated completely because of the parallelism nature of distributed environment, leading to the system error (denoted as ΔE). To understand the influence of T to ΔE, the mathematical expression of ΔE needs to be formulized.To simplify the analysis ofΔE, an abstract control system containing two components is employed here: one component is theController and the other is thePlant. This system is assumed to beAsymptotically Stable. Denoting the original continuous system as CS and the distributed system as DS, the controller and plant can be described as differential equations in CS [15]:(3a)ControllerCS:xct˙=fcxct,uctyct=gcxct,uctuc=yp,(3b)PlantCS:xpt˙=fpxpt,uptypt=gpxptup=yc,where xct and xpt are system states, yct and yp(t) are outputs, uct and upt are inputs, fc(∗) and fp(∗) are state functions, and gc(∗) and gp(∗) are output functions. In DS, the controller and the plant are described as(4a)ControllerDS:xct˙=fcxct,ucdtyct=gcxct,ucdtucdt=ypdt-T,(4b)PlantDS:xpt˙=fpxpt,udqtyp=gpxdtupdt=ycdt-T,where (∗)d is the discretization operator. This operator indicates the operand updates its value to the external world following the advancing step T, just like being discretized.Perturbation Analysis is used here to analyze the composition of ΔE. Perturbation analysis treats ΔE as the outcome of disturbances to CS such as the output delay and discretization. The difference between CS and DS can be determined without formulizing the differential equations of both systems and then comparing the eigenvalues. First of all, the DS system is expressed in the form of perturbation equations:(5a)xp˙=fpxp,gcdxct-T,gpxpt-T=fpxp,gcxc+Δxc,gpxp+Δxp+Δyp+Δycd,(5b)xc˙=fcxc,gpdxpt-T=fcxc,gpxp+Δxp+Δypd,where(i)
Δxc=xct-T-xc, the error of controller state caused by time delay;(ii)
Δxp=xpt-T-xp, the error of plant state caused by time delay;(iii)
Δyp=gp(xpt-T)-gp(xp), the error of plant output caused by time delay;(iv)
Δycd=gcd-gc, the error of controller output caused by discretization;(v)
Δypd=gpd-gp, the plant output error caused by discretization.The system error is represented as(6)ΔE=Δxc,Δxp,Δyp,Δycd,Δypd.To find out the detailed expression of ΔE, a preparation theorem is proposed as follows.


## 4. The Mathematical Analysis

### 4.1. A Preparation Theorem

Based on the assumption ofAsymptotically Stable, the functions of fp(∗), gp∗, fc(∗), and gc(∗) in (3a), (3b), (4a), and (4b) are further assumed to beLipschitz Continuous over a region D, where D is a nonsaturation region defined by xp and xc:(7)D=xp,xc∣xp∈Dxp,xc∈Dxc,gpxp∈Dyp,gcxc,gpxp∈Dyc.For convenience, the time symbol t is omitted unless necessary.   Dxp and Dxc are nonsaturation regions for xp and xc, respectively, that is, the state spaces of the controller and the plant. Dyp and Dyc are nonsaturation regions for yp and yc, respectively, that is, the output spaces of the controller and the plant. With this assumption, a preparation theorem [18] is given.Theorem 1 (preparation theorem).
For an asymptotically stableCS as shown in (3a) and (3b), if (a) the functions f(∗) and g(∗) are continuously differentiable and (b) a Lyapunov function V is defined over an open region D containing the original point (assumed to be the equilibrium point), then a DS can be obtained from CS, in that all start positions lying in an arbitrary interior region (denoted as D1) of D are attracted to another arbitrary interior region (D2) of D1 in finite time. D1 and D2 are defined by sections of V crossed by two level surfaces.This theorem is critical because it guarantees thatDS can be derived from a CS given that CS is stable. The assumption that the original point is the equilibrium is easy to be satisfied by translation transformation of the state variables.Proof.
Define an auxiliary function:(8)αx,Δxc,Δxp,Δyp,Δycd,Δypd=∂V∂xp·fpxp,gcxc+Δxc,gpxp+Δxp+Δyp+Δycd+∂V∂xc·fcxc,gpxp+Δxp+Δypd.We have(9)αx,0,0,0,0=vx˙,where x=(xp,xc) and v(x) is the Lyapunov function of CS.
Define the second auxiliary functionβ(∗) over a region D3:(10)βΔxc,Δxp,Δyp,Δycd,Δypd=sup⁡x∈D3⁡αx,Δxc,Δxp,Δyp,Δycd,Δypd,where D3=D2-D1. β∗ is continuous since α∗ is continuous. Then, we have(11a)αx,Δxc,Δxp,Δyp,Δycd,Δypd<βΔxc,Δxp,Δyp,Δycd,Δypdand then(11b)β0,0,0,0,0=sup⁡x∈D3⁡αx,0,0,0,0=sup⁡x∈D3⁡vx˙.A positive real number s can be found since v(x)˙ is the negative definition to satisfy:(12a)vx˙<-s.Combining (11b), then(12b)β0,0,0,0,0<-s.As a result, a region D4 can always be found in the vicinity of the origin of β∗:(13)D4=Δxc,Δxp,Δyp,Δycd,Δypd∣Δxc,Δxp,Δyp,Δycd,Δypd=ΔE<ρ,where ρ is a positive real number defining the radius of D4. Another positive real number s1 (0<s1<s) can be found to satisfy (14) over D4:(14)βΔxc,Δxp,Δyp,Δycd,Δypd<-s1.Considering inequality (11a), we have(15)αx,Δxc,Δxp,Δyp,Δycd,Δypd<-s1,x∈D3.Consider (9), and we have(16)vx˙<-s1,x∈D3.To integrate both sides of (16) in [0,t],(17a)Vxt-Vx0<-s1∗t,(17b)Vxt<Vx0-s1∗t,where x0∈D3 is the start point of state trajectory. Inequation (17b) means the state trajectory in DS is strictly constrained by a diminishing function in D3. Considering the fact that v(x)˙<0, we have(18a)Vx0<VD1,(18b)Vxt<VD1-s1∗t,where VD1 is the value of V(x) where the state trajectory crosses the level surface that defines region D1. Inequations (18a) and (18b) mean that the state trajectory of DS will stay inside D1. Considering v(x)˙<0, we also have(19)Vx0>VD2,where VD2 is the value of V(x) where the state trajectory crosses the level surface that defines region D2. Then, the solution trajectory of DS will go from the start point to D2 within finite time:(20)Δt<Vx0-VD2s1.Inequation (20) means a DS can be found, whose trajectory will converge to the equilibrium point of CS eventually within finite time, given CS being asymptotically stable and the system error being constrained as (13) shows.


### 4.2. The Mathematical Expression of System Error

The mathematical expression ofΔE is studied by analyzing each component of it.


#### 4.2.1.Δxc

Δ
x
c is the variation of the controller’s state value within time interval T. According to the definition in (5a) and (5b),(21a)Δxc=xct-T-xc=∫t-Ttfcxcτ,ucτ·dτ,(21b)Δxc=∫t-Ttfcxcτ,ucτ·dτ.Integrating both sides of (21b),(22)Δxc=∫t-Ttfcxcτ,ucτ·dτ≤∫t-Ttfcxcτ,ucτ·dτ.fc(∗) is bounded in [t-T,t] since it is continuously differentiable (one of the conditions in Theorem 1). Then, we have(23)Δxc≤∫t-TtMfc·dτ=Mfc·T,where Mfc is the upper bound of fc(∗) within [t-T,t].


#### 4.2.2.Δxp

Δ
x
p is the variation of the plant’s state value within time interval T. Following a similar approach in the analysis of Δxc, we have (24)Δxp≤Mfp·T,where Mfp is the upper bound of fp(∗) within [t-T,t].


#### 4.2.3.Δyp

Δ
y
p is the variation of the plant’s output within time interval T. According to (5a) and (5b), we have(25a)Δyp=gpxpt-T-gpxp=gpxp+Δxp-gpxp,(25b)Δyp=gpxp+Δxp-gpxp.Considering gp(∗) isLipschitz Continuous (Theorem 1), we have(26a)Δyp=gpxp+Δxp-gpxp≤Mgp·xp+Δxp-xp=Mgp·Δxp.Considering (24), we have (26b)Δyp≤Mgp·Mfp·T,where Mgp is theLipschitz Const of gp(∗).


#### 4.2.4.Δypd

Δ
y
p
d is caused by discretization on the plant’s output, where the discrete interval is T. Δypd is bounded by the changing rate of the output, as shown in Figure 3.Figure 3
Δ
y
p
d is bounded by the changing rate of the output value.According to the definitions in (5a) and (5b), we have(27)Δypd=gpd-gp≤Δyp.Considering (26b),(28)Δypd≤Mgp·Mfp·T.


#### 4.2.5.Δycd

Similarly, we have(29)Δycd≤Mgc·Mfc·T.Combining (6), (23), (24), (26b), (28), and (29),(30)ΔE=Δxc,Δxp,Δyp,Δycd,Δypd≤Mfc21+Mgc2+Mfp2·1+2Mgp2·T,where Mfc2(1+Mgc2)+Mfp2·(1+2Mgp2)·T defines the upper bound of ΔE. Once an expected ΔE is given (denoted as ΔE∗), Tmax⁡ can be determined: (31)Tmax⁡=ΔE∗Mfc21+Mgc2+Mfp2·1+2Mgp2for any T<Tmax⁡, and the overall error will be less than ΔE∗. Equation (31) also indicates the way to partition CS will influence Tmax⁡: the more the parts partitioned (supposing each part is deployed on a node), the greater the denominator of the right part of (31), thus the smaller that of the allowable Tmax⁡.


## 4.1. A Preparation Theorem

Based on the assumption ofAsymptotically Stable, the functions of fp(∗), gp∗, fc(∗), and gc(∗) in (3a), (3b), (4a), and (4b) are further assumed to beLipschitz Continuous over a region D, where D is a nonsaturation region defined by xp and xc:(7)D=xp,xc∣xp∈Dxp,xc∈Dxc,gpxp∈Dyp,gcxc,gpxp∈Dyc.For convenience, the time symbol t is omitted unless necessary.   Dxp and Dxc are nonsaturation regions for xp and xc, respectively, that is, the state spaces of the controller and the plant. Dyp and Dyc are nonsaturation regions for yp and yc, respectively, that is, the output spaces of the controller and the plant. With this assumption, a preparation theorem [18] is given.Theorem 1 (preparation theorem).
For an asymptotically stableCS as shown in (3a) and (3b), if (a) the functions f(∗) and g(∗) are continuously differentiable and (b) a Lyapunov function V is defined over an open region D containing the original point (assumed to be the equilibrium point), then a DS can be obtained from CS, in that all start positions lying in an arbitrary interior region (denoted as D1) of D are attracted to another arbitrary interior region (D2) of D1 in finite time. D1 and D2 are defined by sections of V crossed by two level surfaces.This theorem is critical because it guarantees thatDS can be derived from a CS given that CS is stable. The assumption that the original point is the equilibrium is easy to be satisfied by translation transformation of the state variables.Proof.
Define an auxiliary function:(8)αx,Δxc,Δxp,Δyp,Δycd,Δypd=∂V∂xp·fpxp,gcxc+Δxc,gpxp+Δxp+Δyp+Δycd+∂V∂xc·fcxc,gpxp+Δxp+Δypd.We have(9)αx,0,0,0,0=vx˙,where x=(xp,xc) and v(x) is the Lyapunov function of CS.
Define the second auxiliary functionβ(∗) over a region D3:(10)βΔxc,Δxp,Δyp,Δycd,Δypd=sup⁡x∈D3⁡αx,Δxc,Δxp,Δyp,Δycd,Δypd,where D3=D2-D1. β∗ is continuous since α∗ is continuous. Then, we have(11a)αx,Δxc,Δxp,Δyp,Δycd,Δypd<βΔxc,Δxp,Δyp,Δycd,Δypdand then(11b)β0,0,0,0,0=sup⁡x∈D3⁡αx,0,0,0,0=sup⁡x∈D3⁡vx˙.A positive real number s can be found since v(x)˙ is the negative definition to satisfy:(12a)vx˙<-s.Combining (11b), then(12b)β0,0,0,0,0<-s.As a result, a region D4 can always be found in the vicinity of the origin of β∗:(13)D4=Δxc,Δxp,Δyp,Δycd,Δypd∣Δxc,Δxp,Δyp,Δycd,Δypd=ΔE<ρ,where ρ is a positive real number defining the radius of D4. Another positive real number s1 (0<s1<s) can be found to satisfy (14) over D4:(14)βΔxc,Δxp,Δyp,Δycd,Δypd<-s1.Considering inequality (11a), we have(15)αx,Δxc,Δxp,Δyp,Δycd,Δypd<-s1,x∈D3.Consider (9), and we have(16)vx˙<-s1,x∈D3.To integrate both sides of (16) in [0,t],(17a)Vxt-Vx0<-s1∗t,(17b)Vxt<Vx0-s1∗t,where x0∈D3 is the start point of state trajectory. Inequation (17b) means the state trajectory in DS is strictly constrained by a diminishing function in D3. Considering the fact that v(x)˙<0, we have(18a)Vx0<VD1,(18b)Vxt<VD1-s1∗t,where VD1 is the value of V(x) where the state trajectory crosses the level surface that defines region D1. Inequations (18a) and (18b) mean that the state trajectory of DS will stay inside D1. Considering v(x)˙<0, we also have(19)Vx0>VD2,where VD2 is the value of V(x) where the state trajectory crosses the level surface that defines region D2. Then, the solution trajectory of DS will go from the start point to D2 within finite time:(20)Δt<Vx0-VD2s1.Inequation (20) means a DS can be found, whose trajectory will converge to the equilibrium point of CS eventually within finite time, given CS being asymptotically stable and the system error being constrained as (13) shows.


## 4.2. The Mathematical Expression of System Error

The mathematical expression ofΔE is studied by analyzing each component of it.


### 4.2.1.Δxc

Δ
x
c is the variation of the controller’s state value within time interval T. According to the definition in (5a) and (5b),(21a)Δxc=xct-T-xc=∫t-Ttfcxcτ,ucτ·dτ,(21b)Δxc=∫t-Ttfcxcτ,ucτ·dτ.Integrating both sides of (21b),(22)Δxc=∫t-Ttfcxcτ,ucτ·dτ≤∫t-Ttfcxcτ,ucτ·dτ.fc(∗) is bounded in [t-T,t] since it is continuously differentiable (one of the conditions in Theorem 1). Then, we have(23)Δxc≤∫t-TtMfc·dτ=Mfc·T,where Mfc is the upper bound of fc(∗) within [t-T,t].


### 4.2.2.Δxp

Δ
x
p is the variation of the plant’s state value within time interval T. Following a similar approach in the analysis of Δxc, we have (24)Δxp≤Mfp·T,where Mfp is the upper bound of fp(∗) within [t-T,t].


### 4.2.3.Δyp

Δ
y
p is the variation of the plant’s output within time interval T. According to (5a) and (5b), we have(25a)Δyp=gpxpt-T-gpxp=gpxp+Δxp-gpxp,(25b)Δyp=gpxp+Δxp-gpxp.Considering gp(∗) isLipschitz Continuous (Theorem 1), we have(26a)Δyp=gpxp+Δxp-gpxp≤Mgp·xp+Δxp-xp=Mgp·Δxp.Considering (24), we have (26b)Δyp≤Mgp·Mfp·T,where Mgp is theLipschitz Const of gp(∗).


### 4.2.4.Δypd

Δ
y
p
d is caused by discretization on the plant’s output, where the discrete interval is T. Δypd is bounded by the changing rate of the output, as shown in Figure 3.Figure 3
Δ
y
p
d is bounded by the changing rate of the output value.According to the definitions in (5a) and (5b), we have(27)Δypd=gpd-gp≤Δyp.Considering (26b),(28)Δypd≤Mgp·Mfp·T.


### 4.2.5.Δycd

Similarly, we have(29)Δycd≤Mgc·Mfc·T.Combining (6), (23), (24), (26b), (28), and (29),(30)ΔE=Δxc,Δxp,Δyp,Δycd,Δypd≤Mfc21+Mgc2+Mfp2·1+2Mgp2·T,where Mfc2(1+Mgc2)+Mfp2·(1+2Mgp2)·T defines the upper bound of ΔE. Once an expected ΔE is given (denoted as ΔE∗), Tmax⁡ can be determined: (31)Tmax⁡=ΔE∗Mfc21+Mgc2+Mfp2·1+2Mgp2for any T<Tmax⁡, and the overall error will be less than ΔE∗. Equation (31) also indicates the way to partition CS will influence Tmax⁡: the more the parts partitioned (supposing each part is deployed on a node), the greater the denominator of the right part of (31), thus the smaller that of the allowable Tmax⁡.


## 4.2.1.Δxc

Δ
x
c is the variation of the controller’s state value within time interval T. According to the definition in (5a) and (5b),(21a)Δxc=xct-T-xc=∫t-Ttfcxcτ,ucτ·dτ,(21b)Δxc=∫t-Ttfcxcτ,ucτ·dτ.Integrating both sides of (21b),(22)Δxc=∫t-Ttfcxcτ,ucτ·dτ≤∫t-Ttfcxcτ,ucτ·dτ.fc(∗) is bounded in [t-T,t] since it is continuously differentiable (one of the conditions in Theorem 1). Then, we have(23)Δxc≤∫t-TtMfc·dτ=Mfc·T,where Mfc is the upper bound of fc(∗) within [t-T,t].


## 4.2.2.Δxp

Δ
x
p is the variation of the plant’s state value within time interval T. Following a similar approach in the analysis of Δxc, we have (24)Δxp≤Mfp·T,where Mfp is the upper bound of fp(∗) within [t-T,t].


## 4.2.3.Δyp

Δ
y
p is the variation of the plant’s output within time interval T. According to (5a) and (5b), we have(25a)Δyp=gpxpt-T-gpxp=gpxp+Δxp-gpxp,(25b)Δyp=gpxp+Δxp-gpxp.Considering gp(∗) isLipschitz Continuous (Theorem 1), we have(26a)Δyp=gpxp+Δxp-gpxp≤Mgp·xp+Δxp-xp=Mgp·Δxp.Considering (24), we have (26b)Δyp≤Mgp·Mfp·T,where Mgp is theLipschitz Const of gp(∗).


## 4.2.4.Δypd

Δ
y
p
d is caused by discretization on the plant’s output, where the discrete interval is T. Δypd is bounded by the changing rate of the output, as shown in Figure 3.Figure 3
Δ
y
p
d is bounded by the changing rate of the output value.According to the definitions in (5a) and (5b), we have(27)Δypd=gpd-gp≤Δyp.Considering (26b),(28)Δypd≤Mgp·Mfp·T.


## 4.2.5.Δycd

Similarly, we have(29)Δycd≤Mgc·Mfc·T.Combining (6), (23), (24), (26b), (28), and (29),(30)ΔE=Δxc,Δxp,Δyp,Δycd,Δypd≤Mfc21+Mgc2+Mfp2·1+2Mgp2·T,where Mfc2(1+Mgc2)+Mfp2·(1+2Mgp2)·T defines the upper bound of ΔE. Once an expected ΔE is given (denoted as ΔE∗), Tmax⁡ can be determined: (31)Tmax⁡=ΔE∗Mfc21+Mgc2+Mfp2·1+2Mgp2for any T<Tmax⁡, and the overall error will be less than ΔE∗. Equation (31) also indicates the way to partition CS will influence Tmax⁡: the more the parts partitioned (supposing each part is deployed on a node), the greater the denominator of the right part of (31), thus the smaller that of the allowable Tmax⁡.


## 5. Experiment

A distributed simulation is constructed in this section to show the determination ofT. The original continuous system is a point mass system with fraction as follows:(32)x1˙=x2,x2˙=-kmx1-Rmx2,where k is spring const, m is mass, and R is friction coefficient. This system is asymptotically stable; the trajectories of x1 and x2 are shown in Figure 4.Figure 4
The trajectories ofx1, x2, where k = 5, R=5, m=1, and the initial value x1t=0=x2t=0=1.A distributed system can be partitioned from the original system, where the structure has the block-diagram representation (see Figure5).Figure 5
This simple system is partitioned into two parts by the dash line to form a distributed system.The differential equations of these two parts are(33a)Component  I:x1˙=fx1=x2y1=gx1=x1,(33b)Component  II:x2˙=fx2=-kmx1-Rmx2y2=gx2=x2,where y1, y2 are the output of components I and II, respectively. These two parts are deployed on two nodes, forming the simple distributed system.Equation (31) indicates that the parameters reflecting system dynamic need to be determined firstly, which are as follows:Mfx1:
the upper bound of functionfx1 within any time interval [t-T,t];Mgx1:
theLipschitz Constof function gx1;Mfx2:
the upper bound of functionfx2 within any time interval [t-T,t];Mgx2:
theLipschitz Constof function gx2.From the engineering perspective, these parameters can be directly observed from the simulation results of the original system. In this case, Mfx1 is determined by the lower bound of x1˙ as Figure 6 shows.Figure 6
The plot ofx1˙. Mfx1 is determined by the bound of it in this example.As a result,Mfx1=-1.3844. According to the definition ofLipschitz Const, Mgx1 is limited by the upper bound of gx1˙, that is, the upper bound of x1˙ in this case. Then, we have Mgx1=-1.3844. Following similar procedures, we have Mfx2=-1.7796, Mgx2=-1.7796. According to (31), Tmax⁡ is estimated as follows:(34)Tmax⁡≈Mfx221+Mgx22+Mfx12·1+2Mgx12≈0.211·ΔE.The distributed system described in (33a) and (33b) is simulated. The trajectories of the system states are compared to those of the original continuous system, as Figure 7 shows. The threshold value of ΔE is 0.5, and then Tmax⁡ is estimated to be 0.1055 s according to (34). Three configurations, where T=0.1 s, T=0.15 s, and T=0.2 s, are tested in the form of distributed simulation. As expected, only the case T=0.1 s satisfies that both x1 and x2 do not exceed the threshold value of ΔE.Figure 7
The comparison ofx1 and x2. In the distributed simulation, T is set to be 0.1 s, 0.15 s, and 0.2 s, respectively; the trajectories of x1 and x2 are compared in each case.
(a)
(b)
(c)
(d)
(e)
(f)However, the actual error does not exceed the threshold value immediately whenT is configured to be greater than 0.1 s, as shown in Figures 7(c) and 7(d). The reason is that Tmax⁡ is not strictly estimated since many associated parameters use their boundary values.


## 6. Conclusion

The mathematical expression ofΔE helps us to gain the insight of system error produced in the construction of the distributed system using the partitioning approach. Giving an expected threshold of ΔE, a proper advancing step T of the distributed system can be determined. A larger T, compared to the integral step in the nondistributed system, will reduce the data-exchange frequency between simulation nodes, leading to the reduction of demands on system timing performance and network bandwidth. Then, the simulation cost is saved eventually. In fact, this approach can also play a role in multicore or multi-CPU parallel computing environments.However, for systems that may become unstable after system partitioning, this approach is not so convenient since aLyapunov function of the continuous system needs to be found firstly. The parameters defining specific regions satisfying (13)–(15) also need to be determined. More work needs to be done to improve this approach’s availability.



---

*Source: 461513-2015-09-16.xml*