# Quantized State Based Simulation of Time Invariant and Time Varying Continuous Systems

Authors: Xiao Song; Yaofei Ma; Wei Zhang; Jiangyun Wang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Continuous system can be discretized for computer simulation. Quantized state systems (QSS) method has been used to discretize time invariant systems based on the discretization of the state space. A HLA based QSS method is proposed in this paper to address issues of real-time advancements in simulation and an aircraft control example was introduced to illustrate our method. Moreover, to simulate time varying systems, a novel approach is also proposed and exemplified with a practical case.


---

## Body

## 1. Introduction

Continuous system simulation is a topic which has advanced significantly with the appearance of modern computers. Based on classic methods for numerical resolution of ordinary differential equation (ODE) such as Euler and Runge-Kutta, several variable steps and implicit ODE solver methods were developed. These modern methods which usually make use of iteration rules and symbolic manipulation allow the efficient simulation of complex systems, including differential algebraic equation (DAE) and variable structure systems.Although there are several differences between the mentioned ODE solver algorithms, all of them share a property: they are based on time discretization. That is, they give a solution obtained from a difference equation system (i.e., a discrete-time model) which is only defined in some discrete instants.A completely different approach for ODE numerical simulation has been being developed since the end of the 90s. In this new approach, the time discretization is replaced by the state variables quantization and a discrete event simulation model (within the DEVS formalism framework) is obtained instead of a discrete-time one.Discrete event system specification (DEVS) [1–3] is a formalism which allows representing and simulating any system having a finite number of changes in a finite interval of time. In that way, systems modelled by petri nets, state charts, event graphs, and even difference equations can be seen as particular cases of DEVS models.The origin of the quantization based integration methods can be found in the definition of quantized systems and their representation in terms of DEVS models [4]. This idea was reformulated with the addition of hysteresis and it was formalized as a simulation method for ODEs in [5, 6] where the quantized state systems (QSS) were defined. Moreover, there exist some other approaches such as Backward QSS and Centered QSS [7–9]. But all the previous works using QSS to simulate continuous systems did not address the issue how to implement QSS in real-time distributed simulations such as high level architecture (HLA) and distributed interactive simulation (DIS). They did not mention the approach to extend QSS to time varying simulations either. As such, this paper proposes QSS based researches on both real-time and time varying simulation methods, where QSS must be extended for simulation. To accomplish this, we develop a way to keep the real-time synchronization in high level architecture (HLA) distributed simulation environment and propose a method to simulate time varying system.The paper is structured as follows. In Section2, key technologies of QSS are introduced and analyzed. Section 3 introduces the proposed distributed QSS framework to address real-time systems and a validation example is designed and results are shown and analyzed. In Section 4, a new method of time varying case is proposed and illustrated with a practical example. Finally, conclusions are drawn and future works are discussed.


## 2. Quantized State Systems (QSS)

Consider a continuous system model:(1)x˙t=fxt,ut,yt=gxt,ut.x
(
t
) represents a set of state variables, u(t) represents a set of input variables, x˙(t) represents a set of state derivate variables, y(t) represents a set of output variables, f· represents the state transform function, and g· represents the output function.To deduce quantization function, letD={d0,d1,…,dr} be a set of real numbers, where di-1<di,1≤i≤r. Let Ω be the set of piecewise continuous real valued trajectories and x∈Ω is a continuous trajectory. Let b:Ω→Ω be a mapping and let q=bx, where the trajectory q satisfies(2)qt=di+1,ifxt=di+1∧qt-=di∧i<r,di-1,ifxt=di-ε∧qt-=di∧i>0,di,otherwise.Then the mapb is hysteretic quantization function [2]. q(t) is a set of quantized state variables used to simulate xt of (1), where the discrete values di are called quantization levels and the distance Δd≜di+1-di is defined as quantum, which is usually constant. The width of the hysteresis window is ε. The values d0 and dr are the lower and upper saturation values.Replacext on the right side of (1) with q(t). The approximated model is as the following equation shows:(3)x˙t=fqt,utyt=gqt,ut.Kofman and Junco [1] proved that if u(t) has piecewise constant and bounded trajectories and f· is continuous and bounded, then qt,xt,x˙t, and yt have piecewise constant trajectories.


## 3. Real-Time Simulation of Time Invariant System in HLA

Although previous works [7–9] introduced hysteretic quantization function above, they did not address the issue how to fulfill time management in real-time distributed simulations such as DIS and HLA systems. To tackle this problem, we propose using RTFederate in distributed simulations such as HLA and DIS systems. This method is about real-time management and is generic to all distributed simulations. Therefore we will use HLA as our example system in this section.


### 3.1. RTFederate: A Time Control Member

As shown in Figure1, RTFederate mainly has three components: Time Memory, Real Timer and RTI Interface. They are designed to ensure that the time scale in simulation is identical with that of real world time. RTFederate would be an independent federate in HLA, and its time management strategy is set as both Constrained and Regulation.Figure 1
The structure of RTFederate module.For the component of Time Memory, it needs to set a minimum time step denoted asTmin⁡; it would be a positive integer.RTFederate module would get the allowed advance time from RTI through RTI Interface and then add it to the Time Memory. If the valueT in Time Memory is smaller than Tmin⁡, Time Memory would request time advance with Tmin⁡ through RTI Interface. Otherwise, Time Memory would transfer the value T to Real Timer and then reset the value T to zero. Real Timer would count T of unit time. T means the largest integer which is smaller than T. At the end of the count, Real Timer would request time advance with Tmin⁡ through RTI Interface.


### 3.2. Real-Time Simulation

Each of the subsystems would be discretized with QSS method and be treated as a federate in HLA. The time management mode of all the federates should be set as “both Constrained and Regulating”. Also, RTFederate module needs to be set with the same mode. In HLA, Runtime Infrastructure (RTI) software would allow time advance only if all of the federates request to advance and the least requested time step would be allowed. Figure2 shows the procedure of time advance.Figure 2
The procedure of time advance.Only if RTFederate especially requests to advance, RTI would allow its time advancement. So the simulation speed would be controlled by RTFederate module and then it would support real-time simulation.


### 3.3. Example of Time Invariant System

Consider a vertical control system of a plane as the following equation shows:(4)dΔv-dtdΔαdtdΔθdtdΔqdt=-Xv-Xα-Xθ0-Zv-Zα010001-Mv+Mα˙Zv-Mα+Mα˙Zα0-Mq-Mα˙·Δv-ΔαΔθΔq+-XδT00-Zδe00-MδT-Mδe+Mα˙ZδeΔδTΔδe.The symbols in the equation are listed below.The Symbols in (4). The symbols are as follows: Δv-: change value of plane velocity, state variable; Δα: change value of attack angle, state variable; Δθ: change value of pitch angle, state variable; Δq: change value of pitch rate, state variable; Xv: derivative coefficient of x-axis force (body axis) caused by velocity; Xα: derivative coefficient of x-axis force (body axis) caused by attack angle; Xθ: derivative coefficient of x-axis force (body axis) caused by pitch angle; XδT: derivative coefficient of x-axis force (body axis) caused by throttle control lever; Zv: derivative coefficient of z-axis force (body axis) caused by velocity; Zα: derivative coefficient of z-axis force (body axis) caused by attack angle; Zδe: derivative coefficient of z-axis force (body axis) caused by the deflection of elevator; Mv: derivative coefficient of rolling moment caused by velocity; Mα: derivative coefficient of rolling moment caused by attack angle; Mα˙: derivative coefficient of rolling moment caused by changing rate of attack angle; Mq: the derivative coefficient of rolling moment caused by pitch rate; MδT: the derivative coefficient of rolling moment caused by throttle control lever; Mδe: the derivative coefficient of rolling moment caused by the deflection of elevator; Δδe: change of deflection angle of elevator; and ΔδT: changes of the position of thrust lever.The system is discretized with QSS method, the quantum unit isq=0.001 and so as the hysteresis, and the initial conditions are Δv-=Δα=Δθ=Δq=0,Δδe=2°,ΔδT=0°.Figure3 shows the response of Δv-, which is a long period transition process.Figure 3
Curve ofΔv- in the simulation with QSS method.As seen from Figure3, the step intervals became smaller when Δv- gets close to peak values, where the values change severely. On the other hand, the step intervals get bigger when Δv- reaches a relatively stable status, where the value does not change much.Table1 shows the step interval around the first peak value, and Table 2 shows the step intervals in relative stable status. The results indicate that the QSS method reduces the redundant calculation and improves the simulation efficiency.Table 1
The step intervals around the first peak value.
Step width (s)
2.179262
1.010354
1.010355
⋯
1.09538
1.095208
Δ
v
-
0.733
0.736
0.738
⋯
0.741
0.74Table 2
The step intervals around the relative stable region.
Step width (s)
2.111777
3.908981
11.875489
⋯
1.871779
2.88458
Δ
v
-
0.436
0.435
0.436
⋯
0.438
0.439To be simulated in distributed environment each of the state variables would be treated as a federate and in addition to RTFederate module it has five federates in total. Considering the procedure of time advance shown in Figure4, we have the following time advance equations for these federates:Federate1:dΔv-/dt=-Xv·Δv--Xα·Δα-Xθ·Δθ-XδT·ΔδT;Federate2:dΔα/dt=-Zv·Δv--Zα·Δα+Δq-Zδe·Δδe;Federate3:dΔθ/dt=Δq;Federate4:dΔq/dt=(-Mv+Mα˙Zv)·Δv-+(-Mα+Mα˙Zα)·Δα-(Mq+Mα˙)·Δq-MδT·ΔδT+(Mα˙Zδe-Mδe)·Δδe.Figure 4
Procedure of time advance.


## 3.1. RTFederate: A Time Control Member

As shown in Figure1, RTFederate mainly has three components: Time Memory, Real Timer and RTI Interface. They are designed to ensure that the time scale in simulation is identical with that of real world time. RTFederate would be an independent federate in HLA, and its time management strategy is set as both Constrained and Regulation.Figure 1
The structure of RTFederate module.For the component of Time Memory, it needs to set a minimum time step denoted asTmin⁡; it would be a positive integer.RTFederate module would get the allowed advance time from RTI through RTI Interface and then add it to the Time Memory. If the valueT in Time Memory is smaller than Tmin⁡, Time Memory would request time advance with Tmin⁡ through RTI Interface. Otherwise, Time Memory would transfer the value T to Real Timer and then reset the value T to zero. Real Timer would count T of unit time. T means the largest integer which is smaller than T. At the end of the count, Real Timer would request time advance with Tmin⁡ through RTI Interface.


## 3.2. Real-Time Simulation

Each of the subsystems would be discretized with QSS method and be treated as a federate in HLA. The time management mode of all the federates should be set as “both Constrained and Regulating”. Also, RTFederate module needs to be set with the same mode. In HLA, Runtime Infrastructure (RTI) software would allow time advance only if all of the federates request to advance and the least requested time step would be allowed. Figure2 shows the procedure of time advance.Figure 2
The procedure of time advance.Only if RTFederate especially requests to advance, RTI would allow its time advancement. So the simulation speed would be controlled by RTFederate module and then it would support real-time simulation.


## 3.3. Example of Time Invariant System

Consider a vertical control system of a plane as the following equation shows:(4)dΔv-dtdΔαdtdΔθdtdΔqdt=-Xv-Xα-Xθ0-Zv-Zα010001-Mv+Mα˙Zv-Mα+Mα˙Zα0-Mq-Mα˙·Δv-ΔαΔθΔq+-XδT00-Zδe00-MδT-Mδe+Mα˙ZδeΔδTΔδe.The symbols in the equation are listed below.The Symbols in (4). The symbols are as follows: Δv-: change value of plane velocity, state variable; Δα: change value of attack angle, state variable; Δθ: change value of pitch angle, state variable; Δq: change value of pitch rate, state variable; Xv: derivative coefficient of x-axis force (body axis) caused by velocity; Xα: derivative coefficient of x-axis force (body axis) caused by attack angle; Xθ: derivative coefficient of x-axis force (body axis) caused by pitch angle; XδT: derivative coefficient of x-axis force (body axis) caused by throttle control lever; Zv: derivative coefficient of z-axis force (body axis) caused by velocity; Zα: derivative coefficient of z-axis force (body axis) caused by attack angle; Zδe: derivative coefficient of z-axis force (body axis) caused by the deflection of elevator; Mv: derivative coefficient of rolling moment caused by velocity; Mα: derivative coefficient of rolling moment caused by attack angle; Mα˙: derivative coefficient of rolling moment caused by changing rate of attack angle; Mq: the derivative coefficient of rolling moment caused by pitch rate; MδT: the derivative coefficient of rolling moment caused by throttle control lever; Mδe: the derivative coefficient of rolling moment caused by the deflection of elevator; Δδe: change of deflection angle of elevator; and ΔδT: changes of the position of thrust lever.The system is discretized with QSS method, the quantum unit isq=0.001 and so as the hysteresis, and the initial conditions are Δv-=Δα=Δθ=Δq=0,Δδe=2°,ΔδT=0°.Figure3 shows the response of Δv-, which is a long period transition process.Figure 3
Curve ofΔv- in the simulation with QSS method.As seen from Figure3, the step intervals became smaller when Δv- gets close to peak values, where the values change severely. On the other hand, the step intervals get bigger when Δv- reaches a relatively stable status, where the value does not change much.Table1 shows the step interval around the first peak value, and Table 2 shows the step intervals in relative stable status. The results indicate that the QSS method reduces the redundant calculation and improves the simulation efficiency.Table 1
The step intervals around the first peak value.
Step width (s)
2.179262
1.010354
1.010355
⋯
1.09538
1.095208
Δ
v
-
0.733
0.736
0.738
⋯
0.741
0.74Table 2
The step intervals around the relative stable region.
Step width (s)
2.111777
3.908981
11.875489
⋯
1.871779
2.88458
Δ
v
-
0.436
0.435
0.436
⋯
0.438
0.439To be simulated in distributed environment each of the state variables would be treated as a federate and in addition to RTFederate module it has five federates in total. Considering the procedure of time advance shown in Figure4, we have the following time advance equations for these federates:Federate1:dΔv-/dt=-Xv·Δv--Xα·Δα-Xθ·Δθ-XδT·ΔδT;Federate2:dΔα/dt=-Zv·Δv--Zα·Δα+Δq-Zδe·Δδe;Federate3:dΔθ/dt=Δq;Federate4:dΔq/dt=(-Mv+Mα˙Zv)·Δv-+(-Mα+Mα˙Zα)·Δα-(Mq+Mα˙)·Δq-MδT·ΔδT+(Mα˙Zδe-Mδe)·Δδe.Figure 4
Procedure of time advance.


## 4. QSS for Time Varying System

### 4.1. QSS Method for Time Varying System

In the previous section, a QSS method of time invariant system is proposed. However, in practical systems, there exist time varying state variables. And in this section, we proceed to propose a QSS approach of time varying system (TVS).Consider the TVS as follows:(5)x˙=Atx.Whent=ti  (i≥0), (5) follows to be (6)x˙=Atix.At this point, state derivative isx˙i and simulation advances to ti+1′. When t=ti+1′, (5) can be rewritten as(7)x˙=Ati+1′xand the state derivative is x˙i+1. If x˙i and x˙i+1 have the same signs, the state derivative of [ti,ti+1′] is regarded as (1/2)x˙i+x˙i+1 and simulation is advanced to ti+1 based on selected quantum of x. If x˙i and x˙i+1 have different signs, it means the sign changes during [ti,ti+1′]; then advance time should be shortened half as ti+1′′=0.5∗(t+ti+1′):(8)x˙=Ati+1′′x.Reconsiderx˙i+1; if x˙i and x˙i+1 have the same signs, then state derivative of [ti,ti+1′′] is (1/2)x˙i+x˙i+1 and simulation advances to ti+1. Otherwise the advance time is shortened half again until the two state variables have the same signs.


### 4.2. Example of Time Varying System

Consider(9)x˙=-2tt+1-2x.Equation (9) is an instance of TVS, with the initial conditions: x0=1,t0=1. At initial time t0,x˙0=-2×1×(1+1)-2×1=-0.5. The quantum of x is set to 0.001 and the advancing step is 0.001/(0.5)=0.002. In the first step, the simulation time advances to ti+1′=1+0.002=1.002 and x˙1=-2×1×(1.002+1)-2×1=-0.499. As x˙0 and x˙1 have the same signs, state derivative of t∈1,1.002 is x˙=(1/2)-0.5-0.499=-0.4995.The analytical solution of (9) is (10)x=e-2ln⁡t+1-ln⁡2+(1/t+1)-0.5and Figure 5 shows the trajectories of (9) produced by QSS method and the analytical solution, respectively. As we can see, they are almost identical with each other, except the difference in step intervals. The QSS method is more computationally efficient.Figure 5
The solution trajectories of QSS method (the stair line) and the analytical solution (the line).The trajectory values at some time points are listed in Table3. We can also get the conclusion that the result of our method corresponds to analytic curve with minor errors.Table 3
Comparing QSS with analytic result.
Time (s)
1.001
1.003
1.005
⋯
1.017
1.019
QSS solution
1
0.999
0.998
⋯
0.992
0.991
Analytical solution
0.999500
0.998501
0.997503
0.991536
0.990545


## 4.1. QSS Method for Time Varying System

In the previous section, a QSS method of time invariant system is proposed. However, in practical systems, there exist time varying state variables. And in this section, we proceed to propose a QSS approach of time varying system (TVS).Consider the TVS as follows:(5)x˙=Atx.Whent=ti  (i≥0), (5) follows to be (6)x˙=Atix.At this point, state derivative isx˙i and simulation advances to ti+1′. When t=ti+1′, (5) can be rewritten as(7)x˙=Ati+1′xand the state derivative is x˙i+1. If x˙i and x˙i+1 have the same signs, the state derivative of [ti,ti+1′] is regarded as (1/2)x˙i+x˙i+1 and simulation is advanced to ti+1 based on selected quantum of x. If x˙i and x˙i+1 have different signs, it means the sign changes during [ti,ti+1′]; then advance time should be shortened half as ti+1′′=0.5∗(t+ti+1′):(8)x˙=Ati+1′′x.Reconsiderx˙i+1; if x˙i and x˙i+1 have the same signs, then state derivative of [ti,ti+1′′] is (1/2)x˙i+x˙i+1 and simulation advances to ti+1. Otherwise the advance time is shortened half again until the two state variables have the same signs.


## 4.2. Example of Time Varying System

Consider(9)x˙=-2tt+1-2x.Equation (9) is an instance of TVS, with the initial conditions: x0=1,t0=1. At initial time t0,x˙0=-2×1×(1+1)-2×1=-0.5. The quantum of x is set to 0.001 and the advancing step is 0.001/(0.5)=0.002. In the first step, the simulation time advances to ti+1′=1+0.002=1.002 and x˙1=-2×1×(1.002+1)-2×1=-0.499. As x˙0 and x˙1 have the same signs, state derivative of t∈1,1.002 is x˙=(1/2)-0.5-0.499=-0.4995.The analytical solution of (9) is (10)x=e-2ln⁡t+1-ln⁡2+(1/t+1)-0.5and Figure 5 shows the trajectories of (9) produced by QSS method and the analytical solution, respectively. As we can see, they are almost identical with each other, except the difference in step intervals. The QSS method is more computationally efficient.Figure 5
The solution trajectories of QSS method (the stair line) and the analytical solution (the line).The trajectory values at some time points are listed in Table3. We can also get the conclusion that the result of our method corresponds to analytic curve with minor errors.Table 3
Comparing QSS with analytic result.
Time (s)
1.001
1.003
1.005
⋯
1.017
1.019
QSS solution
1
0.999
0.998
⋯
0.992
0.991
Analytical solution
0.999500
0.998501
0.997503
0.991536
0.990545


## 5. Conclusions and Future Work

The paper discretizes the continuous system using QSS method which is based on the quantization of the state space to improve the simulation efficiency. QSS method transforms the continuous system to the discrete event model. The paper introduces both time invariant and time varying methods to support different system simulations. QSS method guarantees the stability and convergence under some conditions and it is not adapted to all the continuous system simulations yet. We need to extend the application range in the future.



---

*Source: 141607-2015-05-04.xml*