Causal-Block Diagrams: A Family of Languages for Causal Modelling of Cyber-Physical Systems

,


Introduction
The design process of complex systems, aided by the technology advances in the last century, is rapidly shifting from small scale development of isolated systems, to large scale development of integrated systems [15].
The graphical representation of such systems using blocks and arrows is one of the first methods used to represent systems.One of the benefits of this notation is that complex systems can be hierarchically decomposed into sub-systems, thus providing a way to deal with complexity.Causal Block Diagrams (CBD) is a formalisation of this intuitive graphical notation.
Originally, CBDs were used to represent and simulate analog circuits [15,280,28,128], with the most commonly used blocks illustrated in Table 4.1.Nowadays, this formalism is widely used in the modelling and simulation of systems that comprise a physical and a controller component, as depicted in Figure 4. 1.In this kind of system, the controller, often implemented in software, monitors the activity of physical processes by means of sensors, takes appropriate decisions, and influences the physical processes through actuators.This architecture can be generalised to networked software and many other processes, not necessarily physical.Examples include cruise controllers, thermostats, robotic arms, and insulin pumps.
For the purposes of introducing the CBD formalism, we hold on to the traditional view of a physical process being controlled by a software component, also known as a feedback control system [16].

Controller
Plant Actuators Sensors Physical Cyber Fig. 4.1: Generic embedded system structure.
In the next section, a simple running example will be introduced, as well as some necessary background concepts.In the remaining sections, CBDs are introduced gradually in three different flavors: algebraic, discrete time, and continuous time CBDs.These are distinguished by the class of blocks at the disposal of the modeler.The gradual presentation allows for a deeper understanding of all the concepts related to modelling and simulation of CBDs.The last few sections of the chapter deal with advanced concepts, related to the simulation of CBDs.

Background
A dynamical system is characterised by a state and a notion of evolution rules.The state is a set of point values in a state space.The evolution rules describe how the state evolves over an independent variable, usually time.
The domain of the time variable dictates the kind of dynamical system: if time ranges over a continuous set, usually R, then the dynamical system is continuous time; if time ranges over a countable set, usually N, then the dynamical system is discrete time.

Cruise Controller
Car (Plant) Thro�le Tachometer Desired Speed Physical Cyber Fig. 4.2: Cruise control system.
Consider the cruise control system depicted in Figure 4.2.The car is a dynamical system whose state (e.g., the velocity) evolves continuously in time.The state of software controller, on the other hand, evolves discretely through time, as any software variable does, by consequence of assignment instructions.
The system in Figure 4.2 is an example of a feedback control system: the physical system -the car -is actuated by control inputs generated by control software -the cruise controller.The sensors are the tachometers that translate wheel revolutions per minute into instantaneous velocity.The actuators are the motor throttle and brakes.The cruise controller software decides, based on current speed of the car, which amount of thrust (throttle or brakes) should be applied to restore the car to a desired speed despite drag, weight, and other factors.
Obviously, it is not desirable to wait for a car prototype to be built, in order to test the control software.This motivates the use of two dynamical systems: (a) A continuous dynamical system which acts as a mock-up of a real car, and (b) A discrete time dynamical system which acts as a mock-up of the software unit.
With these models, early evaluation of many different control strategies can be performed, at a much lower cost.Furthermore, the chosen controller model can then be used to automatically generate the software code.

Models of Physical Systems
Physical systems are inherently continuous: their state evolves continuously through time.Ordinary Differential Equations (ODE) describe how physical quantities change continuously in time, and are thus ideal models to describe the behaviour of physical systems.Initial value problems are ODEs that have a constraint on the initial state.Together with the ODEs, extra equations can be added to represent the outputs of the system.Formally, we consider models of the form: x (t) = F (x, u, t), y(t) = G(x, u, t), x(0) = x 0 , (4.1) where x(t), x (t), y(t), u(t) are vectors, x(t) = [x 1 (t), . . ., x n (t)] T represents the state vector, u(t) represents the input, and x 0 is the initial state.The state here denotes the minimal information required to determine, together with the input u(t) and function F (x, u, t), the complete future states of the system.Function G is the output function.

Dynamic behavior model
A possible model of the dynamic behavior of the car, to be used in the design of the cruise controller in the Cruise Controller example is: where v is the velocity, v , T is the thrust force and k depends on the air density and car shape.This model assumes that the car moves in a straight line and neglects any effects that gravity might induce.Reasonable assumptions for early experimentation.

Discrete Time Models
While the solution of ODEs is continuous, the state of the software unit, in the cruise control system presented in Example 4.2, can only evolve discretely, by the nature of the digital computer on which it runs.For mathematical treatment (e.g., proving that the cruise controller is correct), differential equations may be used for an early specification of the control software.However, when it comes to simulating those models in a digital computer, the only available option is to use discrete time models.First order difference equations, and a constraint on their initial state, allow the specification of such models.They take the form ), where s denotes the step, x [s] is the state vector at step s, u [s+1] is the input vector, y [s] the output vector, and x 0 the initial value of the state vector.The new vector x [s+1] is computed from the old one x [s] and input u [s+1] , according to the specification function F. The output function G allows values to be read from the dynamical system.The repeated application of functions F and G yields the discrete evolution of the state and output vectors.Difference equations are sometimes written as x (0) = x 0 .

Software Controller
The cruise control software can be described by the following difference equation: where v d is the velocity that the car should be kept at (input); v is the actual velocity (input) of the car; v s+1] is the instantaneous error; e [s] denotes the accumulated error (state); T [s] is the thrust force to be transmitted to the car (output).Finally, K p and K i are constant parameters of the controller.
The controller in Example 4.2.2.1 gets the car velocity input v [s] from the readings of the tachometer (recall Figure 4.2).If the differential equation Eq. 4.2 is used to model the car, then v(t) is a continuous quantity and so we can relate it to the input of the controller by v [s] = v(s × ∆t), where ∆t denotes the constant interval of time between two successive tachometer readings.From now on, assume that Eq. 4.2 is coupled to Eq. 4.4 in this manner.
Intuitively, the thrust force T [s] is proportional to the instantaneous error v [s]  d − v [s] and to the accumulation of errors e [s] from previous steps.
The following two paragraphs give an intuitive rationale for each component of Eq. 4.4.
When the instantaneous error is large, the current velocity is far away from the desired one, so the thrust force should be large in order to ensure that the car quickly accelerates/brakes toward the desired speed.When the instantaneous error is small, the car is almost at the desired speed, so the thrust force should be smaller to avoid causing discomfort to the driver.For now, neglect the K i e [s] component in the thrust force calculation.After a while, if the thrust force given by T s] becomes symmetric (same magnitude, opposite direction) with the drag force in Eq. 4.2, the car acceleration will be null and its speed will be kept constant.However, that speed will not be exactly equal to the desired speed, because −kv 0 =⇒ v [s]  d − v [s]  0. This is where the K i e [s] component, neglected until now, has its use.
K i e [s] component.This component accumulates the instantaneous error over time, and contributes to the thrust force accordingly.Suppose that the thrust force is counteracted by the drag force and the car is kept at a constant speed, below the desired velocity, just like in the previous paragraph.Then, the accumulated error will keep growing, ensuring that the second component continues to increase the thrust force until it overcomes the drag force.Notice that this might cause the car to overshoot the desired speed, which can be dangerous.The accumulated error will start decreasing once that happens, decreasing the thrust force.The choice of parameters K p and K i is an important part of tuning the controller.

Discretisation of Differential Equations
The example introduced in this section represents a typical feedback control system, the majority of which are developed using CBDs.In the following sections it will become clear how this is achieved.

Algebraic Causal Block Diagrams
Algebraic CBDs are CBDs where the only atomic blocks permitted are algebraic ones: summation, negation, inversion, product, raise to power and roots.These can be used to represent systems where there is no notion of time and no evolving state.In other words, the time is a constant: now.
While it may seem restricted, these kind of systems arise in the study of the steady state behaviour of dynamic systems, and are used to represent algebraic functions.

Steady State
Consider the car dynamics in Eq. 4.2.The steady state behaviour of the car happens when it is not accelerating.That is, for known constants T, k: The equation gives insight about the torque required to keep the car at the same speed: T = kv.The larger the drag force, the larger the torque, and the more energy is required.

Syntax
The main constituents of a CBD are blocks and connections between blocks.Blocks can be atomic or composite.
Composite blocks stand for an external CBD, specified elsewhere.These blocks will be drawn with a dashed contour.In algebraic CBDs, atomic blocks can be summation, negation, inversion, product, raise to power, and roots.These will be denoted with the appropriate mathematical symbol.Since a block can have more than one input and more than one output, the notion of ports is essential to distinguish between inputs and between outputs.Atomic blocks have up to two input ports -depending on the operation -and one output port.Composite blocks can have any number of input and output ports.The ports associated with a block will not be drawn explicitly but they are part of the CBD and have identifiers (ids).The directed connections will make clear which input ports and output ports are associated with a block.When there is a need, the input port id is shown at the border of the associated block.In the specification of a composite block, the input and output ports are represented as triangles.Whether the port is input or output is clear from the context.
Blocks are also referred to by ids.The id of a port is comprised of the name of the port, and the id of the parent block.The id of a block is formed by its name, along with the id of its parent CBD.Note the recursive definition.While the id is almost never visible in the graphic representation, it is always defined.Some of the uses of ids are: the unambiguous description of connections between ports and the unambiguous identification of individual blocks after the flattening process (Section 4.3.2).
More often, the names of the blocks will be depicted in the graphical representations, to enhance their readability.Names are not ids, they are a part of the id.For instance, the product blocks in Figure 4.3a have two ports with distinct ids, even though these are not depicted in the graphical representation.In the same picture, the name v of the input and the name d of the output port are shown.Similarly, the name c of the composite block is shown.Notice that in Figure 4.3b the same name c also denotes the output port.There is no ambiguity as the composite block and the output port have distinct ids, even though they have the same name.
Whenever a composite block is used, all its internal blocks adopt different ids, based on the id of the CBD where the composite block is used.For instance, the fact that the composite block c is used in the CBD of Figure 4.3a means that, when processing that CBD, the ids of the inner blocks/ports of c (detailed in Figure 4.3b) include the id of c.This has two important consequences: 1. the id of any element depends ultimately on where it is being used, or where any of its parents are being used; 2. if the composite blocks are replaced by their specification in a flattening process, there will be no two ids alike, thus ensuring the well formedness of the CBD.

Semantics
The meaning of an algebraic CBD is an association of a value to each of the ports in the CBD.It can be conveyed in two general ways: by writing the mathematical equations that correspond to the CBD (translational semantics), or by giving an algorithm which computes the value associated with each input/output port (operational semantics).
To simplify both these approaches, it is assumed that all composite blocks are replaced by their specification in a flattening process.This process is done recursively until all composite blocks have been replaced by their specification [232].The following aspects are important to ensure the well-formedness of the flattened CBD: 1.After replacing a composite block, their input/output ports (e.g., the triangular ones in Figure 4.3b) will be connected from both sides.These are redundant ports and hence substituted by a single connection.2. The id of the replaced composite block is still part of the ids of its inner blocks/ports.This ensures uniqueness among ids after the flattening process is complete.
Figure 4.4 shows the result of replacing the composite block c with its specification (in Figure 4.3b).The ids, shown explicitly in the picture, contain the ids of the composite block replaced.During the process, the port c of the composite block became redundant and thus was removed.
The meaning of a flattened CBD is the same as the original CBD.Every block/port has a unique id and every connection is between ports.The semantics of CBDs can be thus explained assuming they have been flattened.

Translational Semantics
Given a flattened algebraic CBD, the equations that it represents can be written down using the rules shown in Table 4.2.
The system of equations that results from applying the following rules to each port, block and connection, of an algebraic CBD, can then be solved for the unknowns to get the values associated with each port in the CBD.The flattened algebraic CBD depicted in Figure 4.4 is translated into the following set of algebraic equations: 1. Assign a unique mathematical variable to the id of each port in the CBD. 2. Let (p, q) denote a connection from port id p to port id q, and let var(p) and var(q) denote the mathematical variables corresponding to p and q, according to the assignment made in Rule 1.Then, the equation associated with the connection (p, q) is var(p) = var(q).3.For each atomic block, let the sequence p 1 , p 2 , . . .denote the list of ids of its inputs ports, and let q denote id of the output port: a.If the block is a constant with value c, then it has no input ports and the resulting equation is var(q) = c; b.If it is a summation, then the resulting equation is var(q) = var(p 1 ) + var(p 2 ); c.If it is a product, then the resulting equation is var(q) = var(p 1 ) × var(p 2 ); d.If it is a negation, then the resulting equation is var(q) = −var(p 1 ); e.If it is an inversion, the resulting equation is var(q) = 1 var( p1 ) ; f.If it is a raise-to-power, the resulting equation is var(q) = var(p 1 ) var( p2 ) ; g.If it is a root, the resulting equation is var(q) = var(p 1 ) The system in Eq. 4.5 can be simplified to a quadratic drag force var where p is the air density, C D the drag coefficient, and A the cross sectional area of the car.Obviously, the value of the input port v has to be known in order to solve for the value of the output port b 2 .o 1 .
Any system of algebraic equations that uses operations supported by the atomic blocks of algebraic CBDs can be represented as an algebraic CBD.The CBD in Figure 4.3b can be drawn directly from the equation where c is the output, and C D , A, p constants.

Operational Semantics
Instead of deferring the responsibility of computing the values associated with each port, to an equation solver, it is possible to do so directly.Two such algorithms are presented.Algorithm 1 presents the dataflow version of the operational semantics.A list of atomic blocks to be computed is revisited iteratively until no blocks remain.A block can be computed only after all the blocks it depends on have been computed.The algorithm terminates after O((#atomic blocks in D) 2 ) iterations.
The inefficiency of this algorithm lies in not taking advantage of the dependencies between blocks to come up with an optimal execution order for blocks.An improved algorithm will be presented later, after formalizing the dependency information between blocks.
Algorithm 1 Data-flow algorithm to evaluate an Algebraic CBD D. Let Q = {q 1 , q 2 , . . .} denote the output ports connected to each input port p j ∈ P, respectively.Let B = block(q 1 ) ∪ block(q 2 ) ∪ . . .be the set of blocks that b i depends on, where block(q j ) is the block associated with port q j or the empty set, if no such block exists.
if B ∩ B = { } then Remark: val(q 1 ), val(q 2 ), . . .have been computed before.val(p) :=C B (b i , val(q 1 ), val(q 2 ), . ..)Let P = {ρ 1 , ρ 2 , . . .} be the set of ports to which port p connects to.val(ρ j ) := val(p), for ρ 1 ∈ P B := B\ {b i } end if end for end while The advantage of Algorithm 1 is its simplicity.It represents the execution model of the dataflow paradigm and, provided that no algebraic loops exist, it always finds the values associated with every port of a flattened CBD.
An algebraic loop arises when a block depends indirectly on itself.It is thus natural to think of the CBD in terms of a dependency graph and identify the cycles thereof.In the CBD of Figure 4.4, blocks c.b 4 and c.b 3 are part of one algebraic loop.Algebraic loops also happen in algebraic systems of equations.For example, in Eq. 4.6, the c variable depends on itself.
Both these loops where introduced artificially for the purposes of illustration.They can easily be removed by reformulating the mathematical expression that the CBD represents.However, in general, not all algebraic loops can be removed by this method and a way to detect them is required.
Given a flattened CBD, its corresponding dependency graph can be created applying the rules in Table 4.3 Table 4.3: Rules for constructing the dependency graph.
1.For each block identified by b, create a unique node v. Let node(b) denote the corresponding node.2. For each connection (p, q) from port id p to port id q, let b p and b q denote the block ids associated with ports p and q, respectively.If p or q have no associated blocks, then ignore this connection and proceed to the next one.Create a directed edge (node(b q ), node(b p )) in the dependency graph, to mark that fact that b q depends on b p .The dependency graph makes the detection of algebraic loops a simple matter of detecting the strong components in the graph.Formally, a strong component S = {n 1 , n 2 , . ..} of a graph G is a set of nodes where, between every n i , n j ∈ S, there are two different paths: p 1 : n i * ⇒ n j and p 2 : n j * ⇒ n i .This implies that every node in a strong component is either the only node in that strong component, or depends on itself, through some other node, also in the same strong component.Tarjan's algorithm [262] accepts a graph and outputs a sorted list of strong components.The sort order of the strong components is a topological order according to the dependencies between strong components.For the example in Figure 4.5, one possible topological order is:

Node Strong component
If no algebraic loops exist in the flattened graph, then the sorted list of strong components returned by the algorithm is just the topological sort of the nodes in dependency graph.In this list, a singleton strong component always appears after the nodes it depends on.
In the case where algebraic loops exist, all nodes belonging to the same algebraic loop will be in the same strong component.Regarding the sort order, non-singleton strong components appear after all components on which it depends on.A strong component depends on other if any one of its comprising nodes depends on at least one of the other's nodes.
These two facts about the sorted strong component list, given by Tarjan's algorithm [262], provide a basis for an improved algebraic CBD operational semantics, that not only can compute the values associated with all output ports much faster than Algorithm 1, but also detects algebraic loops.Algorithm 2 summarizes the steps for computing the values of all ports of a flattened CBD, under the improved algorithm.
Algorithm 2 Evaluation an Algebraic CBD D with support for algebraic loops.Let {q 1 , q 2 , . . .} denote the ids of the input ports of b.
1 , q (i) 2 , . . . .For each Q i there might be input ports whose value is unknown, because these are connected to unknown output ports.Let Qi = q(i) 1 , q(i) 2 , . . .⊆ Q i denote the set of input ports whose value is known.(val(p 1 ), val(p 2 ), . ..) :=S L (b, val( q(1) 1 ), val( q(1) 2 ), . . ., val( q(2) 1 ), . ..) for p i ∈ p 1 , p 2 , . . .do 1 , ρ (i) 2 , . . .be the ports that port p i connects to.val(ρ (i) j ) := val(p i ), for ρ (i) j ∈ P i end for end if end for return val(o 1 ), . . ., val(o m ) end function The S L function computes the values of all unknown ports (whose value is unknown) associated with the blocks in the loop.A input port is unknown when an unknown output port is connected to it.An output port is unknown when the block it is associated with belongs to the strong component.Equivalently, an input port is known when a known output port is connected to it.A known output port is associated with a block that does not belong to the strong component.

T
denote the known values of the input ports.In Eq. 4.7, the unknown input ports are not considered because these depend directly, by algebraic equality, on the output ports connected to them.So finding the values of the unknown output ports is enough to be able to find the values of all unknown ports of the strong component.
The definition of F depends on the atomic blocks that belong to the strong component.If F is linear, then the above equation can be written in the form AX = BU and solved with any technique suitable to solve linear systems of equations (see [60, Chapter 6&7]).Matrices A and B depend on the blocks in the strong component, and the product BU is known.
If F is non-linear, successive substitution techniques ((see [60, Chapter 10])) can be used in an attempt to find X.
Caution has to be taken when non-linear loops are solved, as they might not have a solution, or a unique solution.The iterative methods require initial guess values to be provided for X, and depending on those initial guesses, different solutions might be attained.Both the initial guesses, and the solutions attained have to be physically meaningful, as the equations often represent the characteristics of physical systems (e.g., drag forces, concentrations, etc. . .).
For the algebraic loop containing bocks c.b 4 and c.b 3 in Figure 4.5, the resulting linear system of equations and its analytical solution is: In this section, the syntax and semantics of Algebraic CBDs were described.These are used in systems where there is no notion of evolving state and time.Any algebraic system of equations can be written as an algebraic CBD.
We described two equivalent approaches to obtain the meaning of an algebraic CBD, summarised in Figure 4.6.The solutions are only approximately equal because, in the presence of algebraic loops, these may have to be solved iteratively to get an approximate solution.In the next section, we expand the available atomic blocks to introduce the notion of evolving state via discrete jumps in time.

Discrete-time CBDs
In this section, the Discrete time (DT) CBDs are presented.Syntactically, the only difference to the Algebraic CBDs, is that the DT CBDs allow the modeler to use not only algebraic blocks, but also a step delay block.Because the Delay block has a state, which gets updated whenever the block is computed, the other blocks in a DT CBD no longer have static outputs (as in the algebraic CBDs case), but instead change whenever they are computed.DT CBDs share many similarities with discrete time dynamical systems, presented in Section 4.2.2.

Syntax
The step delay block has two inputs i 1 , i c and one output o 1 .It is represented with a D symbol, as highlighted in the DT CBD of Figure 4.7.The input port i c is called the initial condition and is distinguished by its subscript.

Semantics
The fact that the output of the delay block changes whenever it is computed means that the output of any other block that depends on the delay will also be dynamic.To formalize the multiple different values that any single port can assume, the notion of step is necessary.A step is a natural index that allows the distinction between the different outputs of each block.It is no different than the index used in difference equations, presented in Section 4.2.2.

Translational
The output of the delay block is defined in terms of the input provided at the previous step.This is the essence of difference equations, where the current values are calculated from the previous ones.It is then natural that the meaning of a DT CBD is a set of difference equations.
Similarly to the algebraic CBD case, the flattening process ensures that only atomic blocks remain in the DT CBD.Given a flattened DT CBD, the difference equations that it represents can be written following the rules specified in Table 4.4.Table 4.4: Translational semantics of a flattened DT CBD.
1. Assign a unique mathematical variable to the identifier of each port in the CBD. 2. Let (p, q) denote a connection from port identified by p to port identified by q, and let var(p) and var(q) denote the mathematical variables corresponding to p and q, following the assignment made in Rule 1.Then, the equation associated with the connection (p, q) is var(p) [s+1] = var(q) [s+1] .3. Let p 1 , p 2 , . . .denote the list of ids of inputs ports of an atomic block, and let q denote id of its single output port: a.If the block is a delay block, then the resulting equations are var(q) [s+1] = var(p 1 ) [s] and var(q) [0] = var(p c ) [0] ; b.If the block is a constant block with value c, then the resulting equation is var(q) [s+1] = c; c.If the block is a summation block, then the resulting equation is var(q) [s+1] = var(p 1 ) [s+1] + var(p 2 ) [s+1] ; d.If the block is a product block, then the resulting equation is var(q) [s+1] = var(p 1 ) [s+1] × var(p 2 ) [s+1] ; e.If the block is a negation block, then the resulting equation is var(q) [s+1] = −var(p 1 ) [s+1] ; f.If the block is an inversion block, then the resulting equation is var(q) [s+1] = 1 var( p1 ) [s+1] ; g.If the block is a raise-to-power block, then the resulting equation is ; h.If the block is a root block, then the resulting equation is The result is a set of difference equations, along with initial conditions (see Rule 3(a) of Table 4.4), that can be solved, either to obtain a closed-form solution, or simulated, by an independent solver.As an example, the DT CBD represented in Figure 4.7 corresponds to, after simplification and renaming of variables, the software controller of Eq. 4.4.
Conversely, any difference equation written in the form of Eq. 4.3 can be represented as a DT CBD.This is illustrated in Figure 4  3) can be represented in a DT CBD.The i-th of the vector x is represented as x i .

Operational
When compared to the algebraic CBDs operational semantics, in Algorithm 2, any algorithm that simulates DT CBDs has to compute not single values for variables, but discrete signals.A discrete signal is an ordered list of values, indexed by the step.The operational meaning of the DT CBDs is thus the computation of the discrete time signal associated with each port.That can be done by fixing the step at 0, then computing all values in the CBD, as if it were an algebraic CBD.Then, step is incremented to 1, and the evaluation of all values is repeated, and so on.The fact that, within the same step, the DT CBD is evaluated as if it were an algebraic CBD, allows us to reuse the B (b, val(q 1 ), . . ., s) if b is a delay block then if s = 0 then return val(q c ) (0) , where q c is the id of the initial condition port.else return val(q c ) (s−1) .end if end if if b is a summation block then return val(q 1 ) (s) + val(q 2 ) (s) end if . . .

end function
If there are algebraic loops in the DT CBD, they are handled in the same as way in the algebraic CBDs.A note has to be made, however, about the dependencies of the delay block.At the first step (s = 0), the output of the delay block is equal to the input associated with its initial condition port (val(q c ) (0) in Algorithm 3).At any other step s > 0, the output is computed from the previous step s − 1.This means that, except for s = 0, the delay block has no algebraic dependencies.And at s = 0 it depends on whatever block is connected to its initial condition port.As a result, Rule 2 of Table 4.3 has to be adapted specifically for the Delay block and the current step being computed.
Following the same structure as Section 4.3, this section presented the syntax and semantics, both translational and operational, of DT CBDs.As with algebraic CBDs, the two semantics commute.

Continuous-time CBDs
The meaning of algebraic CBDs is a set of algebraic equations, and the meaning of Discrete time CBDs is a set of difference equations.As shown in Section 4.2.1, differential equations are ideal to represent physical systems, whose state evolves continuously in time.By the end of this section, it will be clear that Continuous Time (CT) CBDs too, are suited to model these systems.

Syntax
Syntactically, CT CBDs include the standard algebraic blocks, a derivative, and an integral block.The Delay block is not included.
The derivative and integral blocks have two inputs i 1 , i c , and one output o 1 .The input subscripted by c denotes the initial condition.Both blocks will be denoted by the appropriate mathematical symbol: d dt and .

Translational Semantics to Differential Equations
The meaning of a flattened CT CBD is a system of Ordinary Differential Equations (ODEs).Table 4.5 shows the rules that build such system.The meaning of Figure 4.9a is represented, after simplification and renaming the variables, in Eq. 4.2.Furthermore, any ODE written in the form of Eq. 4.1 can be translated to a CT CBDs as illustrated in Figure 4.10.

Basics of ODE Discretization
In many ODEs arising in science and engineering, and this includes, by the translational semantics, CT CBDs, a closed-form solution cannot be found.One of the possible ways to get insight into the solution is via simulation.Since most simulations are performed in a digital computer, solutions to ODE's obtained via simulation can only be approximate.Table 4.5: Translational semantics of a flattened CT CBD.
1. Assign a unique mathematical variable to the identifier of each port in the CBD. 2. Let (p, q) denote a connection from port identified by p to port identified by q, and let var(p) and var(q) denote the mathematical variables corresponding to p and q, following the assignment made in Rule 1.Then, the equation associated with the connection (p, q) is var(p)(t ) = var(q)(t ). 3. Let p 1 , p 2 , . . .denote the list of ids of inputs ports of an atomic block, and let q denote id of its single output port: a.If the block is an Integral block, then the resulting equation is var(q)(t ) = t 0 var(p 1 )(τ)dτ + var(p c )(0); b.If the block is a Derivative block, then the resulting equations are var(q)(t ) = var(p 1 ) (t ) and var(q)(0) = var(p c )(0); c.If the block is a summation block, then the resulting equation is var(q)(t ) = var(p 1 )(t ) + var(p 2 )(t ); d. . . .In the following paragraphs, we show how to translate ODEs into difference equations, whose solution, obtained via simulation, approximate the solution of the ODEs.See [?, 60] for a more detailed exposition on numerical approximation methods.The method explained is then used in Section 4.5.2.3 as the basis to describe how a CT CBD is translated into discrete time CBDs, which approximate the solution to the original.
Consider a first order ODE without input: ] T is the state vector, and T is the state derivative function, and x 0 the initial value of x.Let x i (t) denote the i-th state trajectory and x i (t) = F i (x, t) the i-th state derivative.Assuming that x i (t) and any of its derivatives are smooth, it can be approximated around any point t * by the Taylor series [302]: Using Taylor's theorem, it is possible to write the Taylor series expansion in the finite form of a polynomial and a residual in Lagrange form [60]: where ξ (t * ) is an unknown number between t * and t * + ∆t.The residual term x (n+1) i (ξ (t * )) ∆t n+1 (n+1)!denotes the truncation error.It cannot be computed directly but, since any of the derivatives of x i are smooth in all points between t * and t * + ∆t, there exists a maximum constant K < ∞, such that, for any t * and any n, An upper bound can then be written for the remainder term in the Big O notation: Notice that the Big O notation O(∆t n+1 ) highlights the dominant term as ∆t → 0.
Taylor theorem allows us to write the Taylor series taking into account the ODE of Eq. 4.8 and replacing the residual term by its order: For small ∆t < 1 we can neglect the O ∆t 2 term and approximate x i (t * + ∆t) by: Going back to the vector case, this suggests that we can approximate the solution vector x (t) by the following algorithm: Let x [s] = x(s∆t), we get the Forward Euler method to numerically approximate Eq. 4.8: The Taylor series, introduced in Eq. 4.9 also works backward from any point, including the point x i (t * + ∆t): Replacing the derivative by F, from Eq. 4.8, neglecting the residual term, and simplifying gives the Newton's Difference Quotient: Rewritten as a difference equation gives: Contrary to the Forward Euler, it is not possible to get an iterative algorithm immediately out of this method: the vector term x [s+1] depends on itself.This is an algebraic loop (recall Section 4.3.2.2).It requires that the matrix equation be solved for x [s+1] .The presence of these loops distinguishes implicit (with loops) from explicit (without loops) methods.The important point is that, as shown in Section 4.3.2.2, these loops can be solved at each simulation step.

Translational Semantincs to Discrete-time CBDs
As explained in the previous section, differential equations can be discretised to difference equations by means of numerical approximation techniques, which can then be easily simulated.
Since any CT CBD can be translated into an ODE (by Table 4.5), and since the meaning of a discrete time CBD is a system of difference equations, it is natural to wonder whether a discrete time CBD can be transformed directly into a discrete time CBD, which realizes the approximation.The only blocks that need to be approximated are the derivative and the integral.All the other blocks are algebraic.The integral block is left as an exercise.
The derivative block outputs the derivative of its input u, at time t: except at time t = 0, where the output y(0) is given by the input initial condition u c , i.e., y(0) = u c (0).
Applying the Newton's Difference Quotient, from Eq. 4.15, yields: Solving for the output y(t + ∆t) and writing as a difference equation gives: Since the input is not differentiable at time t = 0, the initial condition of the derivative block is provided with an initial value y(0) = u c (0).
It is easy to build a discrete time CBD equivalent to Eq. 4.17 using a delay and algebraic blocks.The delay block will ensure the delayed signal of the input (u [s−1] ) can be obtained.However, at the initial step, s = 0, the delay block has to have an initial condition defined because the value u [−1] , in Eq. 4.17, is unknown.Let u −1 denote this unknown value.u −1 cannot be equal to y [0] .That does not satisfy the initial condition of the derivative, expressed as: To find out the initial condition of the delay, one can rearrange the above equation to get c , which defines the initial condition of the delay.Figure 4.11 shows the transformation rule.
In this section, CT CBDs where introduced, and its meaning given as a translation to discrete time CBDs.
Until now, we have introduced the minimal set of concepts that allow one to use and understand the semantics of CBDs.We skipped over a few details which, to become a proficient user of CT CBDs, have to be covered in the next section.For example, we took for granted that the solution computed by the Forward Euler method is accurate.Furthermore, we have not introduced the operational semantics of continuous time CBDs.Such algorithm can be easily devised for the approximations of the integral and derivative blocks already given -Forward Euler and Newton's Difference Quotient -with special attention to the fact that the approximation of the derivative may introduce algebraic loops in the CBD.However, if those approximation methods are used, the algorithm will be hardly useful: the translation to discrete time CBDs is equivalent and already deals with the algebraic loops problem for free.

Advanced Concepts and Extensions
This section allows you to devise smarter approximation methods, that minimising the error in the approximation.We focus on numerical integration methods, that is, approximations of the integrator block, because these are the most commonly used when modelling physical systems (see Figure 4.10).Finally, we introduce an extension that is widely used in CBDs: logic blocks.These allow higher level reasoning to be used in CBDs, conveying more expressive power to the modeller, but introducing other interesting challenges when it comes to simulation.

Approximation Error
Consider the Forward Euler method in Eq. 4.13.To derive it, the term O(∆t 2 ) of the Taylor series was neglected (recall Eqs.4.9, 4.10, and 4.11).
Let x(t) denote the solution to Eq. 4.8 approximated with Forward Euler, and let x denote the real solution.The first term (x(0) = x(0)) is known from the initial condition of Eq. 4.8.The second term (x(∆t) = x(0) + F (x(0))∆t) deviates from the true solution x(∆t) by an order O(∆t 2 ), which is the residual term ignored in the Taylor series (recall Eq. 4.10 and Eq.4.11).Formally, that is, The third term x(2∆t) deviates further from the true solution not only because of the residual -of order O(∆t 2 )but also because F (x(∆t), ∆t) is evaluated with the approximated term x(∆t) (most likely F (x(∆t)) F ( x(∆t))).The iteration continues and it is easy to see that the error accumulates over the iterations.
In order to analyse the accumulation of errors, it is best to distinguish two kinds of errors: the local truncation error, due to the ignored residual term, and the derivative error, due to evaluating the derivative F at approximated points x.Both these errors contribute to the accumulation of error over time, that is, the global error.
The local truncation error denotes the deviation made by a single step of the numerical method, starting from accurate information, i.e., with no previously accumulated error.
Let x((s + 1) denote the real solution expanded with the infinite Taylor series, and let denote the solution computed across one step of the Forward Euler method, starting from accurate information.The local truncation error is thus given as Studying the global error is more difficult as it depends on the derivative error, which, for a generic analysis, can be any function F. If any error in the parameter of F gets amplified, then the global error will grow faster.If it gets contracted, then the global error will grow in a slower fashion.To formalize, suppose that we know that F satisfies: where 0 ≤ K f < ∞ is a constant.Then, the error at the second step of the Forward Euler can be derived as follows: Notice that, as ∆t → 0, the big O definition implies that Similarly, for the third step: x(3∆t) − x(3∆t) = O(3∆t 2 ).And after s steps, we have x(s∆t) − x(s∆t) = O(s∆t 2 ).
To run the simulation up to time t f the Forward Euler method performs t f /∆t steps, which gives Which says that the global error will be approximately linear in the size of ∆t, as ∆t → 0. For a more accurate expression of the global error of the Forward Euler method, see [60,61].
An important conclusion is that, provided that function F obeys the condition in Eq. 4.20, the global error can be minimised by simply taking smaller ∆t at each step of the simulation using the Forward Euler method.This is called convergence, a property that any useful numerical method should satisfy.
Since there is a limit to how small ∆t can be made in a digital computer, a numerical method which has an higher order of accuracy than O(∆t), for example, O(∆t 2 ), will allow for larger steps to be taken.In the next section, we introduce other numerical methods.

Backward Euler Method
The Taylor series (Eq.4.10) also works backward from any point, including the point (t * + ∆t), as was done in Eq. 4.14.Neglecting the residual term, we get: After some simplifications we get a method that resembles the Forward Euler method: This leads to the Backward Euler method: When compared to explicit methods, the backward Euler requires the solution to an algebraic loop, so it will incur some extra computation at each simulation step.Furthermore, the global and local errors of the backward Euler are of the same order as the Forward Euler method.Their difference lies in the fact that the derivative used to make the estimation of x [s+1] is the closest to it.In the Forward Euler, the derivative is an out-dated one.This has benefits when dealing with stiff systems.See [61, 60, 69] for more details.

Second Order Taylor Method
Until now we have always neglected the term O(∆t 2 ) of the Taylor series.Let us see what happens when we neglect higher order terms.For example, the Taylor series, after neglecting the term O(∆t 3 ), becomes: The derivative dF (x (t * ),t * ) dt can be expanded with the chain rule 1: The second order Taylor series method then becomes: The local truncation error of this method is the neglected term O(∆t 3 ), better than the Euler methods.The disadvantage of this method is that it requires the calculation (symbolically or numerically) of the partial derivatives of F -a costly operation.The global error is in the order of O(∆t 2 ).
Higher order Taylor methods require even more derivative calculations, making them impractical.There are methods that offer that same global error order with far less computation at each step.We show one next.

Midpoint Method
The backward Euler method makes use of the most up-to-date derivative to estimate the solution at t * + ∆t with the disadvantage that it requires more computation to solve the implicit equation.To avoid this, but still trying to be better than Forward Euler, we can try to estimate the derivative at halfway between t * and t * + ∆t and use that derivative to compute x(t * + ∆t): ))∆t.
However, we do not know the value of x(t * + ∆t 2 ).We can use Taylor series again to get Thus we arrive at the midpoint method: • ∆t ∆t (4.25) The midpoint method, Eq. 4.25, can be generalised to 1 Notice that, to be general, we represent the derivative F (x (t * ), t * ) as a function that depends directly on the time.If this is not the case, then ∂F (x (t * ),t * ) ∂t = 0. x • ∆t, s + α p • ∆t with the multi-variate version of the Taylor series, we get: Where the quadratic term was neglected.Plugging it into the previous equation gives: To find the local truncation error, let us find the Taylor series expansion of the true solution and then compare it to the previous equation.The true solution can be expanded as: Applying the chain rule to the derivative yields:

C
, and assuming that these start from a true solution x[s] gives: When solved for the parameters, the above equation gives: As long as the parameters β p , α p , β C1 , β C2 obey the above system of equations, the generic method will have a local truncation error of order O(∆t 3 )), without having to compute any derivative of F. This also shows that the midpoint method is but an element of a family of methods, all with different sets of parameters, called the two stage Runge Kutta methods [69].
By the same argument as the Forward Euler (in Section 4.6.1),we conclude that the global error of the two stage Runge Kutta method is of order O(∆t 2 ).

Adaptive-Step Size
The numerical integration schemes introduced until now use a step size ∆t assumed to be constant throughout the simulation process.These numerical algorithms are computationally expensive in systems where the dynamic behaviour changes slowly except in some limited intervals of time.
Recall that the order of growth of the global error ultimately depends on the Lipschitz constant K f , in Eq. 4.20.This constant represents the worst case deviation of function F as a response to deviations in its parameters2, for all possible values of x(t).
A larger K f indicates that the global error may grow faster, which means that the step size ∆t should be smaller.To clarify: if a system has a large K f , it means that there is at least one pair of values x(t) and x(t) for which F ( x(t)) − F (x(t)) is large.This does not imply the deviations of F are large for all possible pairs of values x(t) and x(t).Furthermore, it does not imply that, if the system were to be simulated in a bounded region (e.g, for 0 < t < t f ), the Lipschitz constant in that region would be smaller.A smaller Lipschitz constant means that the ∆t can be larger.
For a given derivative F, it is hard to find the proper K f in order to pick the right ∆t.
An algorithm that can change the ∆t throughout the simulation, not only leverages the features of each region in the state space to improve the run-time performance of the simulation but also frees the user from the burden of picking an appropriate ∆t.All of this without sacrificing accuracy.
The change of the ∆t has to be triggered by some estimate of the error being committed at each simulation step.Assuming the estimate is available, ∆t is increased if the error becomes too small and decreased if the error is too large.
The challenge is to come up with a good estimate of the error being committed.Suppose we are given two methods, with local truncation errors O(c∆t v ) and O(c ∆t v ), respectively, with c, c , v, v positive constants.Formally, let x(t) be the solution computed by the first method, x(t) by the second, and x(t) be the real solution.Then, after one inaccurate step, solutions x(t + ∆t) and x(t + ∆t) can be written as: The big O notation tells that there exist constants K 1 and K 2 such that, in the limit ∆t → 0, Assuming that c > c and that v < v (the other cases are similar) we have, as ∆t → 0, thus proving that comparing the solutions of the two methods gives an estimate of the error in the same order as the local truncation error of the least accurate method.
From the previous sections, there are two approaches to affect the accuracy of a method: (a) use a smaller step-size and (b) use an higher order approximation method (e.g., the midpoint).
The approach (a) is straightforward: simply take any existing numerical method, compute the solution twice (once with two half steps, and once with a single step), and compare the two estimates.
For an example of approach (b), use the midpoint method to compute the solution x(t), and, at each step, compare it with the result x(t) of the Forward Euler method.It is easy to see that there is some redundant computation in this approach.Fortunately, higher order Runge-Kutta methods can be combined, reusing most of the redundant computation.These are called the Runge-Kutta Fehlberg methods.

Logic Blocks
Decision blocks are widely used in CBDs to increase the expressiveness of the language.The most common decision block is the switch block.The switch block, shown in Figure 4.12, outputs the value u(t) or v(t) depending on the value of c(t).If c(t) ≥ 0, u(t) is the output, otherwise v(t).The translational semantics are: Fig. 4.12: Switch block As will be presented shortly, the operational semantics of this block introduce interesting challenges.

Discontinuity Handling -Zero-Crossing Location
Recall that the simulation of continuous time CBDs can only be performed approximately in a digital computer.See Section 4.5.2.This means that the simulation of a continuous CBDs is actually a discrete set of points x(0), x(1∆t 1 ), x(2∆t 2 ), . . .computed with an adaptive step size method (see Section 4.6.3).
Suppose the time is t and the simulator is going to compute the solution to the output of the switch block y(t + ∆t).Furthermore, assume that y(t) = u(t), that is, c(t) ≥ 0. If c(t + ∆t) ≤ 0, then two issues can be identified: 1. y(t + ∆t) = v(t + ∆t) may be very different than y(t), because v(t + ∆t) u(t + ∆t). 2. t + ∆t may not represent the exact time at which the signal c(t) crossed the zero.That is c(t + ∆t) = 0 − δ, for some > 0.
The second issue implies that, by the intermediate value theorem, there exists at least one point t * ∈ [t, t + ∆t], at which c(t + ∆t) = 0. Ideally, ∆t should be picked in a way such that t + ∆t ≈ t * , thus minimizing δ, for two reasons: 1. Accuracy is improved since all the blocks that depend on the solution x(t + ∆t) will produce outputs that are close to the switching point t * ; 2. Integrator blocks, which apply the numerical methods presented in Section 4.6.2 may need to be aware of the discontinuity in their inputs, caused by the discontinuity of y around the point t * (issue 1 above).
To see why this can be a problem, consider the abstract CBD shown in Figure 4.13.It can be written as a differential equation x (t) = F (x(t)) (recall Figure 4.10).At the time of the discontinuity t * , in the limit → 0, x(t * − ) = x(t * + ), but F (x(t * − )) F (x(t * + )) because of the switch block.This causes a fundamental assumption about the behavior of F-the condition in Eq. 4.20-to be violated.Formally, F (x(t * + )) − F (x(t * − )) ≤ K f x(t * + ) − x(t * − ) ≤ 0 is a contradiction.
Without the Lipschitz condition assumption, it is hard to guarantee an order for the growth of the global error.There are multiple ways to address this problem, once the exact time of the discontinuity is located (see [211,310]).We focus here in the location of the time of the discontinuity (also called root-finding, or zero crossing location in the literature).
A robust yet simple algorithm is the bisection method.As the name hints, the method works by iteratively bisecting the interval.At each iteration it selects the subinterval where the zero-crossing is present to search for the zero location.The algorithm is illustrated in Figure 4.14.The initial steps detects a zero-crossing in the interval between t 1 and t 2 .The iterative procedure evaluates first point a, then point b, then point c and finally point d that is within the tolerance bounds.
Other algorithms are described in the literature [60].

Global Error Euler Method
The global error measures the accumulated deviation of the numerical method from the true solution, across any number of steps.We need to assume that F (x, t) is Lipschitz continuous.At the initial step, the global error is the same as the local truncation error: e (1) = x (1) − x(1) = x (0) + F (x (0) , k∆t)∆t + O ∆t 2 − x (k ) + F (x (k ) , k∆t)∆t = O ∆t 2 .
In general, e Expanding recursively, and assuming that ∆t < 1 we get: To simulate the system from 0 to t f , we require t f ∆t steps.Thus, the error will be: This leads us to conclude that the order of global error of the forward Euler is O (∆t).

Summary
Causal Block Diagrams represent a formalisation of the intuitive graphical notation of blocks and arrows.This chapter introduced the different variants of this formalism, in a gradual manner.The most typical uses for these formalisms are: (i) Algebraic CBDs to study the steady state behaviour of systems; (ii) Discrete time CBDs to represent computation and software components; and (iii) Continuous time CBDs to model physical systems.To connect these three variants, a running example of a cruise control system was used.
Algebraic CBDs represent algebraic systems where there is no notion of passing time.Discrete time CBDs mixes in the passage of time, although at discrete points.These are analogous to difference equations.Finally, continuous time CBDs, were time is a continuum, correspond to differential equations.x The advantage of CBDs over plain difference/differential equations is the natural support for hierarchical descriptions of complex systems, providing a way to manage complexity.
The disadvantage is in the ability to reuse models of physical components, represented as CBDs.Physical objects do not have a notion of inputs and outputs.They are best modeled with equations where any variable can be an input/output, depending on whether it is known (see Acausal modelling chapter).This way, the same component can be reused in many different settings, with its input/outputs defined upon instantiation.In CBDs, the modeler is forced to think of the possible instantiations of the model, and define the inputs/outputs accordingly.
CBDs are widely used in the development of embedded systems.Understanding their semantics and the numerical techniques employed are a stepping stone into understanding other modelling languages.

Literature and Further Reading
Among the references already cited, we highlight: [69] provides an extensive overview of the simulation of continuous systems.[279] gives a good introduction to continuous system modelling and simulation, for someone with a background in Computer Science.
[60] provides a mathematically oriented description of multiple numerical techniques.Last but not least, [211] and [310] provide an overview of the the challenges involved in hybrid system simulation, of which CBDs with logic blocks are part of.

Figure 4 .
3a shows an example of an Algebraic CBD, that calculates the drag force d affecting a car, as it moves with a velocity v, given as input.The composite block c refers to a CBD that calculates the drag coefficient, detailed in Figure4.3b.
be the computed value associated with port identified by p.Let i 1 , . . ., i n be the ids of the input ports associated with the CBD D. Let o 1 , . . ., o m be the ids of the output ports associated with the CBD D. Then, val(i 1 ) := v 1 , . . ., val(i n ) := v n are the values associated with each input ports of D. Let B denote the set of atomic blocks of D not yet computed.Initially, B := all atomic blocks in D. while B { } do for b i ∈ B do Let p denote the single output port of b i .Let P = {p 1 , p 2 , . . .} denote the inputs ports of b i .
. For example, Figure 4.5 shows the dependency graph of the flattened CBD shown in Figure 4.4.
Figure 4.5 illustrates the strong components of the dependency graph.As expected, the blocks c.b 4 and c.b 3 are part of the same strong component.
, v 1 , . . ., v n ) Let val(p) be the computed value associated with port identified by p.Let i 1 , . . ., i n be the ids of the input ports associated with the CBD D. Let o 1 , . . ., o m be the ids of the output ports associated with the CBD D. Then, val(i 1 ) := v 1 , . . ., val(i n ) := v n are the values associated with each input ports of D. Let G denote the dependency graph induced by D. Let SC = (S 1 , S 2 , . ..) denote the sorted list of strong components obtained with Tarjan's algorithm.for S i ∈ SC do if S i = {n } then Let b denote the id of the unique block such that node(b) = n.Let p denote the id of the output port associated with b.

Fig. 4 . 7 :
Fig. 4.7: Discrete-time CBD of the cruise controller with an highlighted Delay block D. Equivalent to Eq. 4.4.

Fig. 4 . 8 :
Fig. 4.8: Difference equation (written in the form of Eq. 4.3) can be represented in a DT CBD.The i-th of the vector x is represented as x i .

Algorithm 3
E A CBD function, defined in Algorithm 2, with some minor changes: Algorithm 3 summarizes the operational semantics.The definition of the C B function is included, to specify the computation of the delay block.The computations of the remaining atomic blocks are trivial.Operational Semantics of an DT CBD D. function E D T CBD(D, v 1 , . . ., v n , N ) Let val(p) be the computed value associated with port identified by p.Let i 1 , . . ., i n be the ids of the input ports associated with the CBD D. Let o 1 , . . ., o m be the ids of the output ports associated with the CBD D. Then, val(i 1 ) := v 1 , . . ., val(i n ) = v n are the values associated with each input ports of D. s := 0 while n ≤ N do E A CBD(D, v (s) 1 , . . ., v (s) n , s) s := s + 1 end while return val(o 1 ), . . ., val(o m ) end function function C CT CBD of the car in the cruise control system of Figure4.2.(b) CT CBD of the drag.

Figure 4 .
Figure 4.9a shows a CT CBD example, with the Drag block specified in Figure 4.9b.

Fig. 4 . 10 :
Fig. 4.10: First order ODE, written in the form of Eq. 4.1, can represented as a CT CBD.The i-th component of the vector x is represented as x i .

Fig. 4 .
Fig. 4.11: Sample CT CBD with a Derivative block (on the left) and the corresponding discrete time CBD (on the right).