Hybrid System Modelling and Simulation with Dirac Deltas

TL;DR

This paper explores the simulation of impulsive differential equations in hybrid systems, presenting symbolic and numerical approaches implemented in a Causal Block Diagram tool, addressing limitations of current simulators.

Contribution

It provides a theoretically grounded implementation of symbolic and numerical impulse simulation methods within a Causal Block Diagram environment, comparing their effectiveness.

Findings

Symbolic approach is more precise for certain systems.

Numerical approach offers better performance in some scenarios.

Guidelines for choosing the appropriate method are provided.

Abstract

For a wide variety of problems, creating detailed continuous models of (continuous) physical systems is, at the very least, impractical. Hybrid models can abstract away short transient behaviour (thus introducing discontinuities) in order to simplify the study of such systems. For example, when modelling a bouncing ball, the bounce can be abstracted as a discontinuous change of the velocity, instead of resorting to the physics of the ball (de-)compression to keep the velocity signal continuous. Impulsive differential equations can be used to model and simulate hybrid systems such as the bouncing ball. In this approach, the force acted on the ball by the floor is abstracted as an infinitely large function in an infinitely small interval of time, that is, an impulse. Current simulators cannot handle such approximations well due to the limitations of machine precision. In this paper, we explore the simulation of impulsive differential equations, where impulses are first class citizens. We present two approaches for the simulation of impulses: symbolic and numerical. Our contribution is a theoretically founded description of the implementation of both approaches in a Causal Block Diagram modelling and simulation tool. Furthermore, we investigate the conditions for which one approach is better than the other.

Tables

Table 1. Table 1: Approximated derivatives of S ​ ( t ) 𝑆 𝑡 S(t) ( Eq.   38 ).

Time	$S\left(t\right)$	$S^{\prime }\left(t\right)$	$S^{\left(2\right)}\left(t\right)$	$S^{\left(3\right)}\left(t\right)$	…	$S^{\left(n-1\right)}\left(t\right)$	$S^{\left(n\right)}\left(t\right)$
${\tau }_d-h$	$0$	$0$	$0$	$0$	…	$0$	$0$
${\tau }_d$	$1$	$1/h$	$1/h^2$	$1/h^3$	…	$1/h^n$	$1/h^n$
${\tau }_d+h$	$1$	$0$	$-1/h^2$	$-2/h^3$	…	$-\left(n-2\right)/h^n$	$-\left(n-1\right)/h^n$
${\tau }_d+2h$	$1$	$0$	$0$	$1/h^3$	…	$\left(\frac{n-2}{2}\right)/h^n$	$\left(\frac{n-1}{2}\right)/h^n$
${\tau }_d+3h$	$1$	$0$	$0$	$0$	…	$-\left(\frac{n-2}{3}\right)/h^n$	$-\left(\frac{n-1}{3}\right)/h^n$
…	…	…	…	…	…	…	…
${\tau }_d+\left(n-1\right)h$	$1$	$0$	$0$	$0$	…	$0$	${\left(-1\right)}^{n-1}\left(\frac{n-1}{n-1}\right)/h^n$
${\tau }_d+nh$	$1$	$0$	$0$	$0$	…	$0$	$0$

Equations

$$y^{(2)}=-g+F_c(t)\text{with}y(0)=y_0\text{and}{y}^{\prime }(0)=v_0$$

$$\left(\begin{array}{cc}{y}^{(2)}& =-g\\ y(0)& =y_0\\ y^{\prime }(0)& =v_0\end{array}\right)\text{for 0≤t<tc, and}\left(\begin{array}{cc}{y}^{(2)}& =-g\\ y(t_c)& =y(t_c^{-})\\ y^{\prime }(t_c)& =-y^{\prime }(t_c^{-})\end{array}\right)\text{ for t≥tc}$$

$$y^{\prime }(t_c^{+})=y^{\prime }(t_c^{-})+{\int }_{t_c^{-}}^{t_c^{+}}-g+F_c(\tau )d\tau$$

$${\int }_{t_c^{-}}^{t_c^{+}}{F}_c(\tau )d\tau =-2y^{\prime }(t_c^{-})$$

$$y^{(2)}=-g-2y^{\prime }(t_c^{-})\delta (t-t_c)\text{with}y(0)=y_0\text{and}{y}^{\prime }(0)=v_0$$

$$\left(f_k(x)=0\text{ if }x\le -\frac{1}{k}\text{ or }\frac{1}{k}\le x\right)\text{ and }\left({\int }_{-1/k}^{1/k}{f}_k(x)dx=1\right)$$

$$\delta (x)=\underset{k\to \infty }{\lim }{f}_k(x)$$

$$\delta (x)=\{\begin{array}{l}0\text{ if }x\ne 0\\ \infty \text{ otherwise }\end{array}$$

$${\phi }_k(x)=\{\begin{array}{l}\varphi (x)\text{ if }-\frac{1}{k}<x<\frac{1}{k}\\ 0\text{ otherwise }\end{array}$$

$$\delta (x)\varphi (x)=\underset{k\to \infty }{\lim }(f_k(x)\varphi (x))=\underset{k\to \infty }{\lim }(f_k(x){\phi }_k(x))=\underset{k\to \infty }{\lim }{f}_k(x){\phi }_k(0)$$

$${\int }_{-\infty }^{\infty }\delta (x)\varphi (x)dx$$

$$=\varphi (0)\left[\underset{k\to \infty }{\lim }{\int }_{-\frac{1}{k}}^{\frac{1}{k}}{f}_k(x)dx\right]=\varphi (0)$$

$${\int }_{-\infty }^{\infty }{\delta }^{\prime }(x)\varphi (x)dx$$

$$=\underset{k\to \infty }{\lim }{\int }_{-1/k}^0\frac{f_k(0)-f_k(-1/k)}{1/k}\varphi (x)dx+{\int }_0^{1/k}\frac{f_k(1/k)-f_k(0)}{1/k}\varphi (x)dx$$

$$=\underset{k\to \infty }{\lim }\left[{\int }_{-1/k}^0{k}^2\varphi (x)dx+{\int }_0^{1/k}-k^2\varphi (x)dx\right]$$

$$=\underset{k\to \infty }{\lim }\left[k^2\varphi (-\frac{1}{2k})\frac{1}{k}-k^2\varphi (\frac{1}{2k})\frac{1}{k}\right]$$

$$=\underset{k\to \infty }{\lim }\left[-\frac{\varphi (\frac{1}{2k})-\varphi (-\frac{1}{2k})}{1/k}\right]$$

$$=-{\varphi }^{\prime }(0)$$

$${\int }_{-\infty }^{\infty }{\delta }^{(i)}(x)\varphi (x)dx=(-1{)}^i{\varphi }^{(i)}(0)$$

$$⟨f,\varphi ⟩\triangleq {\int }_{-\infty }^{\infty }f(x)\varphi (x)dx\text{ for any test function }\varphi$$

$$H(x)=\{\begin{array}{l}1\text{ if }x\ge 0\\ 0\text{ otherwise }\end{array}$$

$$⟨H^{\prime },\varphi ⟩={\int }_{-\infty }^{\infty }{H}^{\prime }(x)\varphi (x)dx={\left[H(x)\varphi (x)\right]}_{-\infty }^{\infty }-{\int }_{-\infty }^{\infty }H(x){\varphi }^{\prime }(x)dx$$

$$=-{\int }_0^{\infty }{\varphi }^{\prime }(x)dx=-{\left[\varphi (x)\right]}_0^{\infty }=\varphi (0)$$

$$H^{\prime }=\delta$$

$$y(t)=\underset{\text{before collision}}{\underbrace{\left(y_0+v_{0}t-\frac{1}{2}g{t}^2\right)}}(1-H(t-t_c))+\underset{\text{after collision}}{\underbrace{\left(y(t_c^{-})-y^{\prime }(t_c^{-})(t-t_c)-\frac{1}{2}g(t-t_c{)}^2\right)}}H(t-t_c)$$

$$⟨y^{(2)},\varphi ⟩$$

$$={\int }_{-\infty }^{\infty }\left(y_0+v_{0}t-\frac{1}{2}g{t}^2\right)(1-H(t-t_c)){\varphi }^{(2)}(t)dt$$

$$+{\int }_{-\infty }^{\infty }\left(y(t_c^{-})-y^{\prime }(t_c^{-})(t-t_c)-\frac{1}{2}g(t-t_c{)}^2\right)H(t-t_c){\varphi }^{(2)}(t)dt$$

$$={\int }_{-\infty }^{t_c}\left(y_0+v_{0}t-\frac{1}{2}g{t}^2\right){\varphi }^{(2)}(t)dt$$

$$+{\int }_{t_c}^{\infty }\left(y(t_c^{-})-y^{\prime }(t_c^{-})(t-t_c)-\frac{1}{2}g(t-t_c{)}^2\right){\varphi }^{(2)}(t)dt$$

$$={\left[\left(y_0+v_{0}t-\frac{1}{2}g{t}^2\right){\varphi }^{\prime }(t)\right]}_{-\infty }^{t_c}-\left({\left[\left(v_{0}t-gt\right)\varphi (t)\right]}_{-\infty }^{t_c}-{\int }_{-\infty }^{t_c}-g\varphi (t)dt\right)$$

$$+{\left[\left(y(t_c^{-})-y^{\prime }(t_c^{-})(t-t_c)-\frac{1}{2}g(t-t_c{)}^2\right){\varphi }^{\prime }(t)\right]}_{t_c}^{\infty }$$

$$-\left({\left[\left(-y^{\prime }(t_c^{-})-g(t-t_c)\right)\varphi (t)\right]}_{t_c}^{\infty }-{\int }_{t_c}^{\infty }-g\varphi (t)dt\right)$$

$$=y(t_c^{-}){\varphi }^{\prime }(t_c)-y^{\prime }(t_c^{-})\varphi (t_c)-g{\int }_{-\infty }^{t_c}\varphi (t)dt$$

$$-y(t_c^{-}){\varphi }^{\prime }(t)-y^{\prime }(t_c^{-})\varphi (t_c)-g{\int }_{t_c}^{\infty }\varphi (t)dt$$

$$=-2y^{\prime }(t_c^{-})\varphi (t_c)-g{\int }_{-\infty }^{\infty }\varphi (t)dt$$

$$=-2y^{\prime }(t_c^{-}){\int }_{-\infty }^{\infty }\delta (t-t_c)\varphi (t)dt-g{\int }_{-\infty }^{\infty }\varphi (t)dt$$

$$={\int }_{-\infty }^{\infty }\left(-g-2y^{\prime }(t_c^{-})\delta (t-t_c)\right)\varphi (t)dt=⟨-g-2y^{\prime }(t_c^{-})\delta (t-t_c),\varphi (t)⟩$$

$$S(t)=s(t)+\sum _{i=0}^n\sum _{{\tau }_j\in \left\{{\tau }_j\right\}}{a}_{ij}{\delta }^{(i)}(t-{\tau }_j)$$

$$S(t^{+})-S(t^{-})=s(t^{+})-s(t^{-})$$

Full text

1. 1 Introduction

2. 2 Background

3. 2.1 Causal Block Diagrams

4. 2.2 Distributions

5. 3 Symbolic Manipulation of Dirac Deltas

6. 4 Numerical Approximation of Dirac Deltas

7. 5 Conclusion

Abstract

For a wide variety of problems, creating detailed continuous models of (continuous) physical systems is, at the very least, impractical. Hybrid models can abstract away short transient behaviour (thus introducing discontinuities) in order to simplify the study of such systems. For example, when modelling a bouncing ball, the bounce can be abstracted as a discontinuous change of the velocity, instead of resorting to the physics of the ball (de-)compression to keep the velocity signal continuous. Impulsive differential equations can be used to model and simulate hybrid systems such as the bouncing ball. In this approach, the force acted on the ball by the floor is abstracted as an infinitely large function in an infinitely small interval of time, that is, an impulse. Current simulators cannot handle such approximations well due to the limitations of machine precision.

In this paper, we explore the simulation of impulsive differential equations, where impulses are first class citizens. We present two approaches for the simulation of impulses: symbolic and numerical. Our contribution is a theoretically founded description of the implementation of both approaches in a Causal Block Diagram modelling and simulation tool. Furthermore, we investigate the conditions for which one approach is better than the other.

Keywords: Dirac delta, hybrid system, distribution theory, causal block diagrams.

1 Introduction

A model of a physical system111In this work, we consider physical systems to be continuous. This excludes Quantum systems. has to pragmatically capture only the essential features for a task at hand [1]. If the task is to analyze some behavior at a macro time scale, it is impractical – or even currently impossible [13] – to create high fidelity continuous models that also explain the behavior at micro time scales. Abstraction is the only way to deal with complexity and create useful models [9].

As a simple example, consider a ball with unit mass that bounces once on the floor. The macro time scale behavior is the overall trajectory of the ball, whereas the micro time scale behavior captures the compression of the ball in the floor. It can be modelled by means of first order Ordinary Differential Equations (ODEs):

[TABLE]

where $y^{(2)}$ denotes the acceleration of the ball (second derivative of the position), and $y_{0}$ and $v_{0}$ are the initial conditions of the dynamical model (constants). $F_{c}(t)$ denotes the force impinged by the floor as a result of the collision, which we want to abstract. There are two possible approaches to find the solution $y(t)$ that satisfies LABEL:eq:bouncing_ball_unabstracted, both based on using conservation of momentum to ensure the abstraction of $F_{c}(t)$ is correct.

Separation of dynamics separates the time continuum into two groups: before and after the contact. Momentum conservation gives the velocity of the ball immediately after the collision. Let $t_{c}$ denote that time at which the contact occurs. LABEL:eq:bouncing_ball_unabstracted is then split in two:

[TABLE]

These are solved separately. $t_{c}$ can be found from the solution to the leftmost ODE. A simulator can be used to compute the solution in the interval $t\in[0,t_{c})$. Afterward, the same simulator can compute the solution in the interval $0\leq t<t_{c}$, with the new initial conditions.

The activity of re-setting the same simulator with new initial conditions is called re-initialization and it is a common approach to the simulation of systems with discontinuities [2].

The Impulsive forces approach acknowledges that, since the velocity is the integral of the acceleration, during the collision, $F_{c}(t)$ must be large enough to cause complete “inversion” of the velocity, even if its shape is unknown. To be concrete, let $t_{c}^{-}$ denote the time at which the collision starts and $t_{c}^{+}$ the time at which the collision ends. Then, after the collision, the velocity of the ball is given by

[TABLE]

Conservation of momentum implies that $y^{\prime}(t_{c}^{+})=-y^{\prime}(t_{c}^{-})$, and the effect that $-g$ has on the integral can be neglected. Thus:

[TABLE]

independently of the shape that $F_{c}(\tau)$ assumes. The Dirac impulse [3] formalizes the abstraction of a large continuous function, whose integral is a known value over a small interval of time. More details will be given in Section 2 but it suffices to postulate that $\int_{0^{-}}^{0^{+}}\delta(\tau)d\tau=1$. Hence, setting $F_{c}(\tau)=-2y^{\prime}(t_{c}^{-})\delta(t-t_{c})$ satisfies Eq. 4 and Eq. 3. With this approach, LABEL:eq:bouncing_ball_unabstracted is rewritten to

[TABLE]

which is an impulsive differential equation [7].

In the bouncing ball example, both the separation of dynamics and impulse based approaches assume that the duration of the contact is infinitely small and are able to reach the same solution. The separation of dynamics allows the modeller more freedom in setting the initial states after an impulse occurs (e.g., the initial condition in Eq. 2 could have been $y(t_{c})=0$ instead of $y(t_{c})=y(t_{c}^{-})$). This means more mistakes can be made (e.g., violating physical laws).

The re-initialization technique can be used for the simulation of impulsive differential equations. However, for complex models where impulses are regularly exchanged, a direct manipulation – symbolic or numerical – of impulses avoids the need to reset the whole state of the simulator. A symbolic manipulation means that impulses are encoded as first class citizens in the operations of the simulator. And by numerical we mean that impulses are approximated as very large values.

Simulators with support for Dirac impulses have been proposed in [11] and more recently in [8], but both approaches deal with the symbolic computation of impulses, and not their numerical approximation, nor evaluation or comparison is advanced for the approach.

Our research hypothesis is: (RH) a simulator which manipulates Dirac deltas symbolically (such as the one proposed in [11]) is more accurate than one that just operates with approximated impulses.

Our contribution can be divided into: (C1) a derivation of the decision and inversion symbolic operations performed on Dirac impulses; (C2) a numerical approximation of impulses; (C3) the conditions for which manipulating impulses symbolically yields more accurate solutions than the numerical approximations;

C1 builds on the work of [11], where many of the symbolic operations on Dirac deltas are derived. Here, we derive the remaining operations – decision and inversion – and provide justification to the use of distributions as a formalization of a correct impulse symbolic simulator. C2 and C3 answer RH. An implementation of a simulator capable of simulating Causal Block Diagrams with impulses, symbolically and numerically, as well as the experiments used in this paper, are available at http://msdl.cs.mcgill.ca/people/claudio/diraccbds.

The next section provides some background and Section 3 derives the main symbolic operations on impulses. Section 4 describes how impulses are approximated and compares the two approaches (C2 and C3).

2 Background

2.1 Causal Block Diagrams

Causal Block Diagrams (CBD) is a formalism used to model causal systems [12, 5], commonly used in tools such as Simulink®. A CBD is expressed with three main entities: Blocks, Links and Ports. Blocks denote operations or refer to other CBDs, defined elsewhere. A Block may contain multiple input ports and a single output port, except if it represents a CBD, which then may contain any number of output ports. Links establish the data-flow by connecting output ports to input ports. Fig. 1 a) shows an example of a CBD that models a bouncing ball. The CBD on the left has 6 blocks, 3 representing the CBDs defined in Fig. 1 b), c) and d).

The basic procedure for the simulation of a CBD model is show in Algorithm 1 [10, 15]. For the sake of simplicity, this procedure assumes that there are no algebraic loops in the model222Algebraic loops are supported by the provided simulator though. 333In Fig. 1 the integrator blocks only depend on their previous inputs, therefore no algebraic loop is formed..

2.2 Distributions

Consider any function $f_{k}(x)$ that satisfies the conditions:

[TABLE]

An example can be $f_{k}(x)=\begin{cases}\max(0,k+k^{2}x)\text{ if }x\leq 0\\ \max(0,k-k^{2}x)\text{ otherwise }\end{cases}$ plotted in dash stroke in Fig. 3 for $k=1$.

A Dirac delta $\delta(x)$ function can be constructed in a way that obeys to Eq. 6 by taking the limit [14], depicted in Fig. 3:

[TABLE]

One can informally say that

[TABLE]

but this is as far as the utility of the statement goes. A more interesting question is what happens when $\delta(x)$ is combined with a smooth (that is, infinitely differentiable) function $\phi(x)$ within an integral, that is, $\int_{-\infty}^{\infty}\delta(x)\phi(x)dx$. According to Equations 6 and 7, and since $\phi$ is smooth, we have that $\delta(x)\phi(x)=\delta(x)\phi(0)$. To see why, let

[TABLE]

Then:

[TABLE]

The observed behavior then becomes:

[TABLE]

which tells us that the Dirac delta $\delta(x-n)$ can be used to sample a signal $\phi$ at point $n$.

Consider now what happens when $\delta(x)$ is differentiated, depicted in Fig. 3:

[TABLE]

One can generalize the derivative of a Dirac delta to:

[TABLE]

These examples show that, even though $\delta(x)$ is not a function in the classical sense, one can observe interesting behavior, and identify an impulse from its interaction with other functions. This is the essential idea of distributions. A distribution is a function identified by the way it interacts with other test functions [6, 4, 14], and not by its direct plot (as with classical functions).

Any function $f(x)$ for which $\int_{-\infty}^{\infty}|f(x)|dx<\infty$, can be a distribution by writing it as $\left\langle f,\phi\right\rangle$, with

[TABLE]

A test function is smooth and becomes 0 outside a bounded interval.

Two distributions are equal if the result of their interactions with any smooth function is the same, that is, $f=g\iff\left\langle f,\phi\right\rangle=\left\langle g,\phi\right\rangle\text{ for all }\phi$. A consequence is that the two sides of Eq. 7 denote the same distribution. This notion of equality allows badly behaved functions (discontinuous and/or with impulses) to exist, as long as their integral is finite. As exemplified before, test functions play a key role because they tend to compensate for the bad shape of the distribution (recall Eq. 10).

Consider the Heaviside distribution $H$ defined by:

[TABLE]

Even though it has a discontinuity at the origin, it can be differentiated as a distribution:

[TABLE]

According to Eq. 8, $\left\langle\delta,\phi\right\rangle=\phi(0)$. So, by distribution equality,

[TABLE]

To see the relationship between distributions and impulsive differential equations, recall the bouncing ball model in LABEL:eq:bouncing_ball_dirac, and suppose that its solution is given by

[TABLE]

where $t_{c}$ is the time at which the collision occurs. This solution is the mix of two trajectories: the one before the bounce, and the one after.

It is clear that derivatives of the solution are discontinuous functions, so differentiating $y(t)$ twice and considering the result as a distribution yields:

[TABLE]

which, according ot the distribution equality, obeys to LABEL:eq:bouncing_ball_dirac. In the next section we apply this theory for the derivation of the symbolic computations on impulses.

3 Symbolic Manipulation of Dirac Deltas

In the formalization of Dirac CBDs, the signals transmitted in the links between blocks are distributions that can have discontinuities and impulses (and any derivative of the impulses). Formally, let $\left\{\tau_{j}\right\}\subset\mathbb{R}$ denote the sequence of all times at which an impulse (or any derivative of it) occurs. Then, following the formalization in [11], the signals are of the form:

[TABLE]

where $s$ denotes a piece-wise continuous impulse-free function, $n$ is the maximum derivative order of any impulse in the signal and $a_{ij}$ is a (possibly zero) constant called the impulse coefficient. At any time $\tau_{j}$ impulses occur, both the right and left limits of $S(t)$ have to have the same impulses. This ensures that the derivations in this section can disregard the left and right limits of the impulse part of the signal. Formally,

[TABLE]

Discontinuities in the impulse free part of the signal, that is, $s(t^{+})\neq s(t^{-})$, are allowed.

The CBD blocks denote manipulations on distributions in the form of Eq. 16. To ensure that these are compositional, we need to show that the output $Y(t)$ of each block $b$ is of the form of Eq. 16, assuming that the inputs are too. Furthermore, to conform strictly to the theory, the moment any signal is plotted, it should be plotted as $\left\langle S(t),\varphi(t)\right\rangle$, for suitable test function. However, the utility of this simulation package is that it allows the modeller to observe the discontinuities and impulses in a signal, so any impulse is plotted with an arrow. The height of the arrow represents the impulse coefficient. An example is shown in Fig. 4 a). We now provide a derivation of each operation.

The Sum block takes two signals ($U(t)$ and $V(t)$) as inputs and produces $Y(t)$ as output:

[TABLE]

Operationally, it adds the two impulse-free signals $u(t)$ and $v(t)$, and adds the coefficients $a_{ij}$ of the impulses in $U(t)$, with the coefficients $b_{ij}$ of the impulses in $V(t)$. $Y(t)$ is of the form of Eq. 16 so the derivation is complete.

The Negation block is derived in the same way. It negates the impulse-free part of the signal and all the coefficients.

[TABLE]

The Switch block outputs either [math] or $1$, according to a condition on the input $C(t)$. A family of pathological cases related to impulses can be identified if $C(t)$ is allowed to have impulses. For example, does a function which at a point $t_{0}>0$ is positive, but has a negative impulse (e.g., $C(t_{0})=c(t_{0})^{2}-\delta(0)$) cross the zero? So $C(t)$ is assumed to be free of impulses, that is, $C(t)=c(t)$. The derivation is then:

[TABLE]

where $H(x)$ is the Heaviside distribution defined in Eq. 12. Alternatively, it can be written as:

[TABLE]

Note however, that $c(t)$ can still be piece-wise continuous. Then the left limit of the output, denoted as $Y(t^{-})$, is evaluated with $c(t^{-})>0$ and the right limit $Y(t^{+})$ with $c(t^{+})\geq 0$. In other words, for any time $t_{0}$ where $c(t_{0}^{-})<0$ and $c(t_{0}^{+})>0$, $H(c(t_{0}^{-}))=0$ and $H(c(t_{0}^{+}))=1$.

The Decision block is a generalization of the Switch block. It forwards one of two inputs ($U(t)$ and $V(t)$) to the output $Y(t)$ according to an input condition $C(t)\geq 0$. Similarly to the Switch block, we have to assume that $C(t)=c(t)$ but in addition, at any time $t_{0}$ where $c(t_{0}^{-})<0$ and $c(t_{0}^{+})\geq 0$, there can be no impulses on either $U(t_{0})$ or $V(t_{0})$. If this were allowed, a signal could be created which violates Eq. 17. Under these assumptions, the derivation is:

[TABLE]

The output of the Decision block can also be written as

[TABLE]

The Derivative block is partially derived as follows:

[TABLE]

The term $\left\langle u^{\prime}(t),\varphi(t)\right\rangle$ must be further developed because it can have discontinuities. For example, when the Derivative block is connected to the output of a Decision block. Let $\left\{t_{d}\right\}$ be the countable sequence of times at which $u(t_{d}^{-})\neq u(t_{d}^{+})$. In the proximity of a discontinuity at time $t_{d}$, $u(t)$ is described as $u(t)=u(t^{-})(1-H(t-t_{d}))+u(t^{+})H(t-t_{d})$. Intuitively, it can be seen as the output signal of a Decision block. $\left\langle u^{\prime}(t),\varphi(t)\right\rangle$ is then:

[TABLE]

The distribution $\left\langle\left[(u(t^{+})-u(t^{-}))H(t-t_{d})\right]^{\prime},\varphi(t)\right\rangle$ can be further simplified:

[TABLE]

In conclusion, around each discontinuity $t_{d}\in T_{d}$, we have:

[TABLE]

This is interpreted as originating an impulse, with a coefficient given by the magnitude of the discontinuity. The derivative signal can itself be discontinuous, if its left limit is also different than its right limit. Across all possible discontinuities $t_{d}\in T_{d}$ of $u(t)$, and abstracting the discontinuities in $u^{\prime}(t)$, the signal in Eq. 27 is in the form of Eq. 16. Concretely:

[TABLE]

Joining Eq. 27 and Eq. 24 yields:

[TABLE]

For example, suppose the input to the Derivative block is the signal

[TABLE]

With $v_{0}$,$g$ constants. Then the output of the Derivative block is

[TABLE]

The Integrator block is derived as:

[TABLE]

Which matches the intuition: the Integrator block computes the normal integration of the impulse-free signal, reduces the order of any impulse derivative, and computes a discontinuity whenever an impulse is integrated (recall Eq. 14). The magnitude of the discontinuity is the impulse coefficient.

The Product of two distributions is not defined in general. However, if one of the input signals to the Product block is impulse-free444This restriction cannot be enforced automatically by the tool., then it can be derived as:

[TABLE]

Essentially, each product between a function and a Dirac delta derivative yields a set of simultaneous impulses of decreasing derivative orders, that are multiplied with increasing derivatives of the function. The computation of the coefficients, therefore, requires that the derivatives of $u(t)$, order up to (and including) $n_{v}$, be known or estimated. To get the intuition, recall Eq. 9 and Fig. 3, and note that

[TABLE]

For example, if $U(t)=-gt$ and $V(t)=1+20\delta^{(2)}(t-1.44)$, then the Product block computes:

[TABLE]

The Inverter block, to be well defined, requires that the input is free of impulses:

[TABLE]

This section has provided the relationship between DiracCBDs and Distributions, and derived each block used in the CBD formalism. In the DiracCBDs simulator the impulses are treated symbolically but the rest of the signals are discretized approximations of the continuous system. In particular, if $h>0$ is the simulation time step (recall Algorithm 1), then the integrator block can be approximated by recursive right Riemann’s sum and the derivative block by finite difference:

[TABLE]

The signals, in the form of Eq. 16, can be encoded in many ways. For example, at each time step, the left and right limit of the signals can be stored so that discontinuities can be easily identified, and a matrix of impulse coefficients kept. Each block then takes the input left/right limit and the matrix, and produces an output of the same form. [11] encodes the left limit of the signal, along with all its derivatives (lazily evaluated). Another example encoding can be seen in [8] that is based in the super dense time encoding. Fig. 4 a) shows the bouncing ball trace produced by our simulator.

4 Numerical Approximation of Dirac Deltas

The previous section offered insight on how to symbolically manipulate impulses in the context of simulation. This manipulation can get complex and counter intuitive, especially with derivatives of impulses (recall Eq. 31 and LABEL:eq:delta_derivative_simple_derivation), and requires a special encoding of the signals. The purpose of this section is to describe an alternative approach by approximating Dirac deltas as large functions, and see how it compares to the symbolic approach.

A Dirac delta impulse is a arbitrary high function at a very small interval of time that obeys to Eq. 6. In a numerical simulation, the smallest interval is $h$, so a good numerical approximation of an impulse $\delta(t-\tau_{d})$ would be to produce a large value $N$ at a simulation step that is immediately after $\tau_{d}$:

[TABLE]

We cannot just pick an arbitrary high value for $N$ (e.g., $N=2^{31}-1$) because of the conditions in Eq. 6.

Eq. 14 holds the key to numerically approximate the impulse. Suppose the Dirac delta is constructed as:

[TABLE]

where $H_{k}(x)$ represents a continuous approximation of the Heaviside function, for large enough $k$: This formulation satisfies the conditions in Eq. 6, and it represents the limit version of the result in Eq. 14. Since $h$ is the smallest interval of time we have, the impulse is approximated as $\delta(t-\tau_{d})\approx H_{1/h}^{\prime}(t-\tau_{d})$, where $H_{{1/h}}^{\prime}$ is approximated as in Eq. 35. The smaller $h$ is, the more accurate the approximations are.

This approach works remarkably well for the bouncing ball model, producing an exactly equal trace shown in Fig. 4 b). Numerical approximations different than those in Eq. 35 may make the traces differ. One immediate effect is that the values plotted for the impulsive force change by several orders of magnitude. Due to the numerical approximation for the derivative, the smaller $h$ is, the larger the impulse will be. While smaller $h$ increases the accuracy, it may cause overflows due to the limits of machine precision. Note that the time at which the ball touches the floor still needs to be accurately located [16] by adjusting $h$ and our simulator does this (both the symbolic and the numerical versions) which explains why Fig. 4 b) has different values for the impulses.

Now consider the following artificial signal and its distributional derivatives up to order $n$:

[TABLE]

We compare this signal with the approximation $\tilde{S}(t)\approx H_{1/h}(t-\tau_{d})$ from Eq. 37, differentiated according to the approximation in Eq. 35. For example, the derivative $H^{\prime}_{1/h}(t-\tau_{d})$ is approximated as:

[TABLE]

Table 1 gives the values of the these approximations around time $\tau_{d}$. The table shows that $(n+1)$ time steps are necessary in the numerical simulator to represent the information that is symbolically stored in an impulse derivative of order $n$. This shifts the numerical solution by $n\times h$ time units, which may cause significant errors. Intuitively, the symbolic impulse derivatives represent a compact version of what the numerical equivalent would do across the multiple infinitesimal time steps that are abstracted away.

In the approximation of $S^{(n)}(t)$, the table, resembling Pascal’s triangle, yields the maximum magnitude of the values computed:

[TABLE]

Eq. 39 can be used to get a rough estimate of the maximum magnitude of values computed in a simulation of an impulsive differential equation. For example, if the maximum impulse derivative order is $n$, and the maximum amplitude of any discontinuity in the simulation is $D$, then the maximum magnitude is $D\binom{k-1}{k-1}/h^{k}$. This result follows from Eq. 39 and the properties of distributions, described in Section 2.2 and Section 3.

To summarize, for models that contain no impulse derivatives, the numerical approximation is equivalent to the symbolic manipulation. For models that have impulse derivatives, the symbolic approach is more accurate than the numerical one, due to the delay introduced, which is proportional to the highest derivative order of the impulses.

5 Conclusion

This work has explored the use of Dirac delta impulses for the modelling and simulation of hybrid systems represented by impulsive differential equations. Section 4 compares a simulator that approximates impulses numerically with a simulator that encodes the impulses explicitly in the signal. We conclude that for models that have no impulse derivatives, both approaches are exactly the same, for the approximations in Eq. 35. The models that we have used do not have impulse derivatives and we failed at finding models of physical systems that include impulse derivatives. The numerical approach described in Section 4 is also more efficient, due to a simpler encoding. The encoding is a stream of real numbers, as opposed to Eq. 16.

For models that have impulse derivatives however, the numerical approach introduces a delay in the signal, proportional to the highest derivative order of the impulses, which causes inaccuracies. Furthermore, if the models are to be run on an embedded platform, where overflows pose a more significant risk, then the symbolic approach may be used. Eq. 39 provides a rough estimate of the magnitudes involved when using the numerical approach. As can be seen, the smaller the simulation step size $h$ is, the larger Eq. 39 will be. This poses an interesting challenge because in general $h$ has to be small in order to get accurate results. For the symbolic approach, we hypothesize that it may require less bits to perform the same computations on an embedded platform, because the magnitudes involved may not be dominated by the impulse coefficients. Finally, the symbolic approach has the benefit of being able to identify when operations that are not defined in theory, are being computed (e.g., an inversion of an impulse). As [11] points out, this is not fool proof because the signals are computed point-wise but, compared to the numerical approach, it is an improvement.

Acknowledgment

We are thankful to the anonymous reviewers for the corrections and comments provided. This research was partially supported by Flanders Make vzw, the strategic research centre for the manufacturing industry, and PhD fellowship grants from (1) the Agency for Innovation by Science and Technology in Flanders (IWT), and (2) Research Foundation - Flanders (FWO).