diff --git a/.nojekyll b/.nojekyll
index 09cbfbe..04b2dcc 100644
--- a/.nojekyll
+++ b/.nojekyll
@@ -1 +1 @@
-dfdcb4d5
\ No newline at end of file
+3444948f
\ No newline at end of file
diff --git a/classes_reset.html b/classes_reset.html
index 806bd05..face92a 100644
--- a/classes_reset.html
+++ b/classes_reset.html
@@ -782,49 +782,49 @@
Reset systems, switched systems, and piecewise affine (PWA) systems can all be viewed as special classes of hybrid systems, and as such can be modelled and simulated using software for hybrid systems (or even general purpose modelling and simulation software).
+
However, there are also some dedicated software tools/packages/libraries for PWA systems:
Although some more free and open-source software packages for PWA systems can be found on the internet, none of them (at as far as we know) is actively developed or at least maintained anymore:
+
+
PWATOOLS toolbox for Matlab, which accompanies the recently published book Rodrigues, Samadi, and Moarref (2019). Unfortunately, this ten-year old toolbox is no longer working with the recent releases of Matlab and the author is no longer maintaining it.
+
PWLTool toolbox for Matlab: some traces of this toolbox can be found on the internet, but this one seems even older, obviously back then accompanying the book Johansson (2003).
+
+
Overall, besides the MPT toolbox that is still being actively developed (by our colleagues at STU Bratislava), not much is currently available within the open-source software domain… :-(
+Johansson, Mikael K.-J. 2003. Piecewise Linear Control Systems: A Computational Approach. Lecture Notes in Control and Information Sciences. Berlin, Heidelberg: Springer. https://doi.org/10.1007/3-540-36801-9.
+
+
+Rodrigues, Luis, Behzad Samadi, and Miad Moarref. 2019. Piecewise Affine Control: Continuous-Time, Sampled-Data, and Networked Systems. Advances in Design and Control. Philadelphia: Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611975901.
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/classes_switched.html b/classes_switched.html
index f68dee5..0449c17 100644
--- a/classes_switched.html
+++ b/classes_switched.html
@@ -803,46 +803,46 @@
Concise (short and yet sufficient for our purposes) introduction to (linear) complementarity problems and systems is in [1]. Besides describing the idea of complementarity in dynamical systems, it also shows how it is related to other modeling frameworks for hybrid dynamical systems. More detailed and yet very accessible introduction is in the thesis [2]. Condensed treatment is in the papers [3] and [4].
+
A readable introduction to the Extended Linear Complementarity Problem is in [5] (it is also freely available as a technical report).
+
The topics of complementarity constraints in dynamical systems and optimization is still being actively researched. A recent publication on QP optimization with complementarity constraints (LCQP) is [6].
+
Numerical methods for nonsmooth dynamical systems that are based on complementary constraints (and implemented in SICONOS software) are comprehensively presented in [7].
W. P. M. H. Heemels, B. De Schutter, and A. Bemporad, “Equivalence of hybrid dynamical models,”Automatica, vol. 37, no. 7, pp. 1085–1091, Jul. 2001, doi: 10.1016/S0005-1098(01)00059-0.
W. P. M. H. Heemels, J. M. Schumacher, and S. Weiland, “Linear Complementarity Systems,”SIAM Journal on Applied Mathematics, vol. 60, no. 4, pp. 1234–1269, Jan. 2000, doi: 10.1137/S0036139997325199.
+
+
+
[4]
A. J. van der Schaft and J. M. Schumacher, “Complementarity modeling of hybrid systems,”IEEE Transactions on Automatic Control, vol. 43, no. 4, pp. 483–490, Apr. 1998, doi: 10.1109/9.664151.
+
+
+
[5]
B. De Schutter and B. De Moor, “The Extended Linear Complementarity Problem and the Modeling and Analysis of Hybrid Systems,” in Hybrid Systems V, P. Antsaklis, M. Lemmon, W. Kohn, A. Nerode, and S. Sastry, Eds., in Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 1999, pp. 70–85. doi: 10.1007/3-540-49163-5_4.
+
+
+
[6]
J. Hall, A. Nurkanovic, F. Messerer, and M. Diehl, “LCQPow – A Solver for Linear Complementarity Quadratic Programs.” arXiv, Nov. 2022. Accessed: Dec. 03, 2022. [Online]. Available: http://arxiv.org/abs/2211.16341
+
+
+
[7]
V. Acary and B. Brogliato, Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics. in Lecture Notes in Applied and Computational Mechanics, no. 35. Berlin Heidelberg: Springer, 2008. Available: https://doi.org/10.1007/978-3-540-75392-6
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/complementarity_simulations.html b/complementarity_simulations.html
index 55a4879..ff43de7 100644
--- a/complementarity_simulations.html
+++ b/complementarity_simulations.html
@@ -754,46 +754,46 @@
Having just discussed the concept of a discrete-event system, we now introduce the most popular modeling framework for such systems: a state automaton, or just an automaton (plural automata). It is also known as a state machine or a (discrete) transition system.
+
+
Definition 1 (Automaton) Automaton is a tuple \boxed{
+G = \{\mathcal X,\mathcal X_0,\mathcal E,\mathcal F\},}
+
+where
+
+
\mathcal X is the set of states (also called modes or locations).
+
\mathcal X_0 \subseteq \mathcal X is the set of initial states.
+
\mathcal E is the set of events (also actions, transition labels, symbols). It is also called alphabet.
+
\mathcal F\subseteq \mathcal X \times \mathcal E \times \mathcal X is the set of transitions. In the deterministic case it can also be narrowed down to a transition functionf:\mathcal X \times \mathcal E \rightarrow \mathcal X. Note that f is then is a partial function, it is not necessarily defined for all combinations of states and events. Sometimes f is used even for multivalued functions: f:\mathcal X \times \mathcal E \rightarrow 2^\mathcal{X}, where 2^\mathcal{X} is a power set (a set of all subsets of X).
+
+
+
+
+
+
+
+
+Some comments on the notation
+
+
+
+
+
The set of states is often denoted by \mathcal Q to spare the letter \mathcal X for the continuous valued state space of hybrid systems.
+
The set of events is often denoted by \mathcal A to spare the letter \mathcal E for the set of transitions (edges in the corresponding graph), because F and f may also need to be spared for the continuous-valued transitions. But then the letter \mathcal A actually fits this purpose nicely because the event set is also called the alphabet.
+
+
+
+
+
Marked states
+
In some literature, the definition of the automaton also includes a set \mathcal X_\mathrm{m} \subseteq \mathcal X of marked or accepting states, in which case the definition of an automaton now includes three (sub)sets of states: \mathcal X, \mathcal X_0 and \mathcal X_\mathrm{m}. \boxed{
+G = \{\mathcal X,\mathcal X_0,\mathcal E,\mathcal F, \mathcal X_\mathrm{m}\}.}
+
+
The marked states are just some states with special roles in the system. Namely, these are the states into which the system should be controlled. I do not particularly like this idea of mixing the model of the system with the requirements, but some part of the community likes it this way.
+
+
+
Automaton as a (di)graph (also a state transition diagram)
+
So far the definition of an automaton was not particularly visual. This can be changes by viewing the automaton as a directed graph (digraph) with. These are the basic rules
+
+
State is represented as a node of the graph.
+
Transition from a given state to another state is represented as an edge connecting the two nodes.
+
Events (actions) are the labels attached to the edges. It is not necessary that each edge has its unique label.
+
+
+
Example 1 (Automaton as a digraph) Consider an automaton defined by these sets: \mathcal X = \{x_1,x_2,x_3\}, \mathcal X_0 = \{x_1\}, \mathcal E = \{e_1,e_2,e_3\}, \mathcal F = \{(x_1,e_1,x_2),(x_2,e_2,x_1),(x_1,e_3,x_3),(x_2,e_2,x_3)\}.
+
+Figure 1: An example automaton as a digraph
+
+
+
+
+
+
+
We may also encounter the following term.
+
+
Definition 2 (Active event function and set) Active event function (actually a multivalued function) \Gamma: \mathcal X \rightarrow 2^\mathcal{E} assigns to each state a set of active events. Active event set \Gamma(x) is the set of active events in a particular state x.
+
+
+
+
Finite state automaton (FSA)
+
This may be regarded as a rather superfluous definition – a finite state automaton (FSA) is a state automaton with a finite set \mathcal X of states. It is also known as a finite state machine (FSM).
Here x_k for some k is the name of a particular state. It is not the name of a (yet to be introduced) state variable; In fact, it can be viewed as its value (also valuation).
+
+
+
+
Some authors strictly distinguish between the state variable and the state (variable valuation),
+
+
similarly as in probability theory random variable X vs its value x, as in F(x) = P(X\leq x);
+
+
some do not, but then it may lead to confusion;
+
yet some others avoid the problem by not introducing state variables and only working with enumerated states.
+
+
+
+
+
+
+
+
+Notational confusion 2
+
+
+
+
Even worse, it is also tempting to interpret the lower index k as (discrete) time, but nope, in the previous k is not the time index.
In continuous-valued dynamical systems, we have a state trajectory, but then time stamps are attached to each visited state.
+
+
+
Example 2 (Beverage vending machine)
+
+
+
+
+
+
+
+
+
+
+Figure 2: Example of a digraph representation of the automaton for a beverage vending machine
+
+
+
+
+
+
+
State sequence (path): waiting, swiped, paid, coke_dispensed, waiting
+
Events sequence: swipe card, accept payment, choose coke, take coke
+
Indeed, the two states coke_dispensed and fanta_dispensed can be merged into just beverage_dispensed.
+
How about other paths? Longer? Shorter?
+
+
+
+
+
+
+
+
+
+
+
+Figure 3: Example of a digraph representation of the automaton for a beverage vending machine with a marked state
+
+
+
+
+
+
The waiting state can be marked (is accepting).
+
+
+
Example 3 (Longitudinal control of a ground vehicle)
+
+
+
+
+
+
+
+
+
+
+Figure 4: Example of a digraph representation of the automaton for a longitudinal control of a ground vehicle
+
+
+
+
+
+
+
+
+
By cruise on I mean switching on some kind of a cruise control system, which keeps the velocity constant.
+
It turns out the optimal control strategy for trains (under some circumstances).
+
Note that some of the events are indeed actions started by the driver, but some are just coming from the physics of the vehicle (transition from braking to zero velocity).
+
+
+
+
Example 4 (Corridor switch)
+
+
+
+
+
+
+
+
+
+
+Figure 5: Example of a digraph representation of the automaton for a corridor switch
+
+
+
+
+
+
Two events associated with one transitions can be seen as two transitions, each with a single event, both sharing the starting and ending states.
+
+
+
Example 5 (JK flip-flop) We now consider the classical JK flip-flop logical circuit. It symbol is in Fig 7 and the truth table follows. Our goal is to represent its functionality using a state automaton.
+
+
+
+
+
+
+Figure 6: Symbol for a JK flip-flop logical circuit
+
+
+
+
+
+
+
J
+
K
+
Q_k
+
Q_{k+1}
+
Description
+
+
+
+
+
0
+
0
+
0
+
0
+
No change
+
+
+
0
+
0
+
1
+
1
+
No change
+
+
+
0
+
1
+
0
+
0
+
Reset
+
+
+
0
+
1
+
1
+
0
+
Reset
+
+
+
1
+
0
+
0
+
1
+
Set
+
+
+
1
+
0
+
1
+
1
+
Set
+
+
+
1
+
1
+
0
+
1
+
Toggle
+
+
+
1
+
1
+
1
+
0
+
Toggle
+
+
+
+
+
+
+
+
+
+
+
+
+
+Figure 7: JK flip-flop as an automaton
+
+
+
+
+
+
+
+
Example 6 (Double intensity switching)
+
+
+
+
+
+
+
+
+
+
+Figure 8: Example of a digraph representation of the automaton for double intensity switching
+
+
+
+
+
+
Obviously we need to introduce time into the automaton…
+
+
+
+
State as the value of a state variable
+
Definition of the state space by enumeration (such as \mathcal X = \{0,1,2,3,4,5\}) doesn’t scale well. As an alternative, a state can be characterized by the value (sometimes also valuation) of a state variable. A state variable is then given by
+
+
the name (for example, x),
+
the “type” (boolean, integer, vector, …).
+
+
+
Example 7 (Examples of state variables)
+
+
Corridor switch: x \in \{\mathrm{false},\mathrm{true}\} (possibly also \{0,1\}).
+
Double intensity switching:
+
+
x \in \{0,1,2\} \subset \mathbb Z,
+
or \bm x = \begin{bmatrix}x_1\\ x_2 \end{bmatrix}, where x_1,x_2 \in \{0,1\}.
+
+
+
+
+
+
State (transition) equation
+
Denoting a new state after a transition as x^+, the state equation reads \boxed{x^+ = f(x,e)}
+
Upon introduction of discrete-time (index) k, it can also be rewritten as x_{k+1} = f(x_k,e_k) or also x[k+1] = f(x[k],e[k]).
+
+
+
+
+
+
+Note
+
+
+
+
+
The function f can be defined by a computer code rather than a clean mathematical formula.
+
The discrete-time index of the event is sometimes considered shifted, that is x_{k+1} = f(x_k,e_{k+1}). You should be aware of this.
+
+
+
+
+
+
Extensions
+
The concept of an automaton can be extended in several ways. In particular, the following two extensions introduce the concept of an output to an automaton.
+
+
Moore machine
+
One extension of an automaton with outputs is Moore machine. The outputs assigned to the states by the output functiony = g(x).
+
The output is produced (emitted) when the (new) state is entered.
+
Note, in particular, that the output does not depend on the input. This has a major advantage when a feedback loop is closed around this system, since no algebraic loop is created.
+
Graphically, we make a conventions that outputs are the labels of the states.
+
+
Example 8 (Moore machine) The following automaton has just three states, but just two outputs (FLOW and NO FLOW).
+
+
+
+
+
+
+
+
+
+
+Figure 9: Example of a digraph representation of the Moore machine for a valve control
+
+
+
+
+
+
+
+
+
Mealy machine
+
Mealy machine is another extension of an automaton. Here the outputs are associated with the transitions rather than the states.
+
Since the events already associated with the states can be viewed as the inputs, we now have input/output transition labels. The transition label e_\mathrm{i}/e_\mathrm{o} on the transion from x_1 to x_2 reads as “the input event e_\mathrm{i} at state x_1 activates the transition to x_2, which outputs the event e_\mathrm{o}” and can be written as x_1\xrightarrow{e_\mathrm{i}/e_\mathrm{o}} x_2.
+
It can be viewed as if the output function also considers the input and not only the state y = e_\mathrm{o} = g(x,e_\mathrm{i}).
+
In contrast with the Moore machine, here the output is produced (emitted) during the transition (before the new state is entered).
+
+
Example 9 (Mealy machine) Coffee machine: coffee for 30 CZK, machine accepting 10 and 20 CZK coins, no change.
+
+
+
+
+
+
+
+
+
+
+Figure 10: Example of a digraph representation of the Mealy machine for a coffee machine
+
+
+
+
+
+
+
+
Example 10 (Reformulate the previous example as a Moore machine) Two more states wrt Mealy
+
+
+
+
+
+
+
+
+
+
+Figure 11: Example of a digraph representation of the Moore machine for a coffee machine
+
+
+
+
+
+
+
+
+
+
+
+
+Note
+
+
+
+
There are transitions from 30 and 40 back to 0 that are not labelled by any event. This does not seem to follow the general rule that transitions are always triggered by events. Not what? It can be resolved upon introducing time as the timeout transitions.
+
+
+
+
Example 11 (Dijkstra’s token passing) The motivation for this example is to show that it is perhaps not always productive to insist on visual description of the automaton using a graph. The four components of our formal definition of an automaton are just enough, and they translate directly to a code.
+
The example comes from the field of distributed computing systems. It considers several computers that are connected in ring topology, and the communication is just one-directional as Fig 12 shows. The task is to use the communication to determine in – a distributed way – which of the computers carries a (single) token at a given time. And to realize passing of the token to a neighbour. We assume a synchronous case, in which all the computers are sending simultaneously, say, with some fixed sending period.
+
+
+
+
+
+
+
+
+
+
+Figure 12: Example of a ring topology for Dijkstra’s token passing in a distributed system
+
+
+
+
+
+
One popular method for this is called Dijkstra’s token passing. Each computer keeps a single integer value as its state variable. And it forwards this integer value to the neighbour (in the clockwise direction in our setting). Upon receiving the value from the other neighbour (in the counter-clockwise direction), it updates its own value according to the rule displayed in the code below. At every clock tick, the state vector (composed of the individual state variables) is updated according to the function update!() in the code. Based on the value of the state vector, an output is computed, which decodes the informovation about the location of the token from the state vector. Again, the details are in the output() function.
+
+
+Show the code
+
struct DijkstraTokenRing
+ number_of_nodes::Int64
+ max_value_of_state_variable::Int64
+ state_vector::Vector{Int64}
+end
+
+functionupdate!(dtr::DijkstraTokenRing)
+ n = dtr.number_of_nodes
+ k = dtr.max_value_of_state_variable
+ x = dtr.state_vector
+ xnext =copy(x)
+for i ineachindex(x) # Mind the +1 shift. x[2] corresponds to x₁ in the literature.
+if i ==1
+ xnext[i] = (x[i] == x[n]) ? mod(x[i] +1,k) : x[i] # Increment if the left neighbour is identical.
+else
+ xnext[i] = (x[i] != x[i-1]) ? x[i-1] : x[i] # Update by the differing left neighbour.
+end
+end
+ dtr.state_vector .= xnext
+end
+
+functionoutput(dtr::DijkstraTokenRing) # Token = 1, no token = 0 at the given position.
+ x = dtr.state_vector
+ y =similar(x)
+ y[1] =iszero(x[1]-x[end])
+ y[2:end] .= .!iszero.(diff(x))
+return y
+end
+
+
+
output (generic function with 1 method)
+
+
+
We now rund the code for a given number of computers and some initial state vector that does not necessarily comply with the requirement that there is only one token in the ring.
+
+
+Show the code
+
n =4# Concrete number of nodes.
+k = n # Concrete max value of a state variable (>= n).
+@show x_initial =rand(0:k,n) # Initial state vector, not necessarily acceptable (>1 token in the ring).
+dtr =DijkstraTokenRing(n,k,x_initial)
+@showoutput(dtr) # Show where the token is (are).
+
+@showupdate!(dtr), output(dtr) # Perform the update, show the state vector and show where the token is.
+@showupdate!(dtr), output(dtr) # Repeat a few times to see the stabilization.
+@showupdate!(dtr), output(dtr)
+@showupdate!(dtr), output(dtr)
+@showupdate!(dtr), output(dtr)
We can see that although initially the there can be more tokens, after a few iterations the algorithm achieves the goal of having just one token in the ring.
+
+
+
+
Extended-state automaton
+
Yet another extension of an automaton is the extended-state automaton. And indeed, the hyphen is there on purpose as we extend the state space.
+
In particular, we augment the state variable(s) that define the states/modes/locations (the nodes in the graph) by additional (typed) state variables: Int, Enum, Bool, …
+
Transitions from one mode to another are then guarded by conditions on theses new extra state variables.
+
Besides being guarded by a guard condition, a given transition can also be labelled by a reset function that resets the extended-state variables.
+
+
Example 12 (Counting up to 10) In this example, there are two modes (on and off), which can be captured by a single binary state variable, say x. But then there is an additional integer variable k, and the two variables together characterize the extended state.
+
+
+
+
+
+
+
+
+
+
+Figure 13: Example of a digraph representation of the extended-state automaton for counting up to ten
+
+
+
+
+
+
+
+
+
+
Composing automata
+
Any practically useful modelling framework should support decomposition of a large system into smaller subsystems. These should then be able to communicate/synchronize with each other. In automata such synchronization can be realized by sending (or generating) and receiving (or accepting) events. A common choice of symbols for the two is !,?, as illustrated in the following example. But these symbols are just one possible convention, and any other symbols can be used.
+
+
Example 13 (Composing automata)
+
+
+
+
+
+
+
+
+
+
+Figure 14: Example illustrating how two automata can be synchronized by sending and receiving events
+
+
+
+
+
+
+
+
+
Languages and automata
+
When studying automata, we often encounter the concept of a language. Indeed, the concept of an automaton is heavily used in the formal laguage theory. Although in our course we are not going to refer to these results, some resources we recommend for our courses do, and so it is useful to understand how automata and languages are related.
+
First, we extend the definition of a transition function in that it accepts the current state and not just a single event but a sequence of events, that is
+
+f: \mathcal X \times \mathcal E^\ast \rightarrow \mathcal X,
+ where \mathcal E^\ast stands for the set of all possible sequences of events.
+
Language generated by the automaton is
+\mathcal L(\mathcal G) = \{s\in\mathcal E^\ast \mid f(x_0,s) \;\text{is defined}\}
+
+
Language marked by the automaton (the automaton is accepting or recognizing that language)
+\mathcal L_\mathrm{m}(\mathcal G) = \{s\in\mathcal L(\mathcal G) \mid f(x_0,s) \in \mathcal{X}_\mathrm{m}\}
+
+
+
Example 14 (Language accepted by automaton)
+\mathcal{E} = \{a,b\}, \mathcal{L} = \{a,aa,ba,aaa,aba,baa,bba,\ldots\}
+
+
+
+
+
+
+
+
+
+
+
+Figure 15: Example of an automaton generating the language \mathcal{L} = \{a,aa,ba,aaa,aba,baa,bba,\ldots\}
+
+
+
+
+
+
What if we remove the self loop at state 0? The automaton then accepts languages starting with a and with b being the last event or immediately followed by a.
+
+
+
What is the language view of automata good for?
+
+
Definitions, analysis, synthesis.
+
We then need language concepts such as
+
+
concatenation of strings: \quad c = ab
+
empty string \varepsilon: \quad\varepsilon a = a \varepsilon = a
+
prefix, suffix
+
prefix closure \bar{\mathcal{L}} (of the language \mathcal L)
+
…
+
+
+
+
+
+
Blocking
+
An important concept in automata is blocking. A state is blocking if there is no transition out of it. An example follows.
+
+
Example 15 (Blocking states) In the automaton in Fig 16, state 2 is blocking. It is a deadlock state. States 3 and 4 are livelock states.
+
+
+
+
+
+
+
+
+
+
+Figure 16: Example of an automaton with blocking states
+
+
+
+
+
+
Language characterization: \bar{\mathcal{L}}_\mathrm{m}(\mathcal G) \sub \mathcal L(\mathcal G).
+
+
+
+
Queueing systems
+
Queueing systems are a particular and very useful class of discrete-event systems. They consist of these three components:
Tradeoff needed between customer satisfaction and fair resources allocation
+
+
+
+
Networks of queueing systems
+
Queueing systems can be interconnected into networks.
+
+
+
+
+
+
+Figure 18: Example of a network of queueing systems
+
+
+
+
+
+
Queueing systems as automata
+
The reason why we mentioned queueing systems in this part of our course is that they can be modelled as automata. And we already know that in order to define and automaton, we must characterize the key components defining the automaton – three in this case:
+
+
events: \mathcal E = \{\text{arrival},\text{departure}\};
+
states: number of customers in the queue
+\mathcal X = \{0,1,2,3,\ldots\}, \quad \mathcal X_0 = \{0\},
+
+
+
+
+
+
+
+
+Note
+
+
+
+
Obviously this is not a finite state automation – unless the queue is bounded – and whether the queue’s length is bounded is a modelling assumption.
+
+
+
+
state transition:
+f(x,e) =
+\begin{cases}
+x+1, & \text{if}\; x\leq 0 \land e = \mathrm{arrival}\\
+x-1, & \text{if}\; x > 0 \land e = \mathrm{departure}.
+\end{cases}
+
+
+
+
+
Queueing system as an automaton
+
+
+
+
+
+
+Figure 19: Queueing system as an automaton
+
+
+
+
+
+
+
+
+
+Note
+
+
+
+
Note how the states correspond to the value of the state variable.
+
+
+
+
Example 16 (Example of a queueing system: jobs processing by a CPU) …
+
+
+
+
Stochastic queueing systems
+
An important extension of the basic concept of a queueing system is the introduction of randomness. In particular, the arrivals can be modelled using random processes. Similarly, the departures given by the delays (the processing time) of the server can be modelled as random.
+
Obviously, the time needs to be included in the automaton, and so far we do not have it there. It is then high time to introduce it.
+
+
+
+
Timed automaton
+
So far, even if the automaton corresponded to a physical system (and did not just represent a generator of a language), the time was not included. The transitions were triggered by the events, but we did not specify the time at which the event occurred.
+
There are, however, many situations when it is useful or even crucial to incorporate time. We can then answer questions such as
+
+
How many events of a certain type in a given interval?
+
Is the time interval between two events above a given threshold?
+
How long does the system spend in a given state?
+
…
+
+
There are several ways how to incorporate time into the automaton. We will follow the concept of a timed automaton with guards (introduced by Alur and Dill). Within their framework we have
+
+
one or several resettable clocks: c_i,\, i=1,\ldots, k, driven by the ODE
+ \frac{\mathrm{d} c_i(t)}{\mathrm d t} = 1, \quad c_i(0) = 0;
+
+
each transition labelled by the tripple {guard; event; reset}.
+
+
+
+
+
+
+
+Note
+
+
+
+
Both satisfaction of the guard and arrival of the event constitute enabling conditions for the transition. They could be wrapped into a single compound condition.
+
+
+
+
Example 17 (Timed automaton with guards)
+
+
+
+
+
+
+
+
+
+
+Figure 20: Example of a timed automaton with guards
+
+
+
+
+
+
+
+
Example 18 (Timed automaton with guards and invariant)
+
+
+
+
+
+
+
+
+
+
+Figure 21: Example of a timed automaton with guards and invariant
+
+
+
+
+
+
+
+
Invariant vs guard
+
+
Invariant (of a location) gives an upper bound on the time the system can stay at the given location. It can leave earlier but not later.
+
Guard (of a given transition) gives an enabling condition on leaving the location through the given transition.
+
+
+
Example 19 (Several trains approaching a bridge) The example is taken from [1] and is included in the demos coming with the Uppaal tool.
G. Behrmann, A. David, and K. G. Larsen, “A Tutorial on Uppaal,” in Formal Methods for the Design of Real-Time Systems, M. Bernardo and F. Corradini, Eds., in Lecture Notes in Computer Science, no. 3185., Berlin, Heidelberg: Springer, 2004, pp. 200–236. doi: 10.1007/978-3-540-30080-9_7.
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/des_automata.html b/des_automata.html
index 915a6dd..5639951 100644
--- a/des_automata.html
+++ b/des_automata.html
@@ -2131,16 +2131,16 @@
We can see that although initially the there can be more tokens, after a few iterations the algorithm achieves the goal of having just one token in the ring.
Well, here we are at last. After these three introductory topics on discrete-event systems, we’ll finally get into hybrid systems.
+
There are two frameworks for modelling hybrid systems:
+
+
hybrid automaton, and
+
hybrid (state) equations.
+
+
Here we start with the former and save the latter for the next chapter/week.
+
First we consider an autonomous (=no external/control inputs) hybrid automaton – it is a tuple of sets and (set) mappings
+\boxed{
+\mathcal{H} = \{\mathcal Q, \mathcal Q_0, \mathcal X, \mathcal X_0, f, \mathcal I, \mathcal E, \mathcal G, \mathcal R\},}
+ where
+
+
\mathcal Q is a set of discrete states (also called modes or operating modes or locations).
using a state variable q(t) attaining discrete values. The variable can also be a vector one, possibly a binary vector state variable encoding an integer scalar state variable.
+
+
+
\mathcal Q_0\subseteq \mathcal Q is a set of initial discrete states.
+
+
It can contain only a single element.
+
+
Example: \mathcal Q_0 = \{\text{off}\}.
+
+
But if it contains more than one element, it can be used to represent uncertainty in the initial state.
+
+
\mathcal X\subseteq \mathbb R^n is a set of continuous states.
+
+
Rather than by enumeration as in the case of discrete state, it is characterized by real-valued state variables x, oftentimes vector ones \bm x. Often they are denoted \bm x(t) to emphasize the evolution in time.
+
+
\mathcal X_0\subseteq \mathcal X is a set of initial continuous states.
+
+
A set of values of the (vector) state variable at the initial time.
+
Often just a single initial state \mathcal X_0=\{x_0\}, but it can be useful to set ranges of the values of individual state variables to account for uncertainty.
+
+
f:\mathcal{Q}\times \mathcal X \rightarrow \mathbb R^n is a vector field parameterized by the location
+
+
Often the dependence on the location q expressed as f_q(x) rather than the more symmetric f(q,x).
+
This defines a state equation parameterized by the location: \dot{x}(t) = f_q(x(t)).
+
It is also possible to consider a set-valued map, replacing f with \mathcal F: \mathcal{Q}\times \mathcal X \rightarrow 2^{\mathbb R^n}, which leads to the differential inclusion\dot x \in \mathcal F_q(x).
+
+
\mathcal I: \mathcal Q \rightarrow 2^\mathcal{X} gives a (location) invariant. It is also called a domain (of the location). The latter term is perhaps more appropriate becaue the term invariant is already too much overloaded in the context of dynamical systems.
+
+
It is parameterized by the location.
+
It is a subset of the continuous-valued state space \mathcal I(q) \subseteq \mathbb R^n.
+
It is a set of values that the state variables are allowed to attain while staying in the given location; if the state of the systems evolves towards the boundary of this set with a tendency to leave it, the system must be ready to leave that location by transitioning to another location.
+
+
+
+
+
+
+
+
+Caution
+
+
+
+
Strictly speaking, \mathcal{I} is a mapping and not a set. Only the mapping evaluated at a given location q, that is, \mathcal{I}(q), is a set.
+
+
+
+
\mathcal E\subseteq \mathcal Q \times \mathcal Q is a set of transitions.
+
+
It is a set of the edges of the graph.
+
Example: \mathcal E = \{(\text{off},\text{on}),(\text{on},\text{off})\}, that is, a two-component set.
+
+
\mathcal G: \mathcal E \rightarrow 2^\mathcal{X} gives a guard set.
+
+
It is associated with a given transition. In particular, \mathcal G(q_i,q_j) is the guard set for the transition (q_i,q_j)\in\mathcal E.
+
+
The guard condition for the given transition is satisfied if x\in \mathcal G(q_i,q_j).
+
If the guard condition is satisfied, the transition is enabled – it may be executed. But it does not have to.
+
The enabled transition must be executed when the state x leaves the invariant set of the original location.
+
+
\mathcal R: \mathcal E \times \mathcal X\rightarrow 2^{\mathcal X} is a reset map.
+
+
For a given transition from one location to another, it resets the continuous-valued state x to a new value within some subset.
+
Often the map is single-valued, r: \mathcal E \times \mathcal X\rightarrow \mathcal X (multivalued-ness can be used to model uncertainty).
+
We also say that the state experiences a jump.
+
+
The state after the jump (associated with the given transition) is reset according to x^+ = r(q_i,q_j, x), or x^+ \in \mathcal R(q_i,q_j, x) in the multivalued case.
+
+
If no resetting of the continuous-valued state takes place, the reset map is defined just as the identity operator with respect to x , that is, r(q_i,q_j, x) = x.
+
+
+
+
Example 1 (Thermostat – the hello world example of a hybrid automaton) The thermostat is a device that turns some heater on or off (or sets some valve open or closed) based on the sensed temperature. The goal is to keep the temperature around, say, 18^\circ C.
+
Naturally, the discrete states (modes, locations) are on and off. Initially, the heater is off. We can identify the first two components of the hybrid automaton: \mathcal Q = \{\text{on}, \text{off}\}, \quad \mathcal Q_0 = \{\text{off}\}
+
The only continuous state variable is the temperature. The initial temperature not not quite certain, say it is known to be in the interval [5,10]. Two more components of the hybrid automaton follow: \mathcal X = \mathbb R, \quad \mathcal X_0 = \{x:x\in \mathcal X, 5\leq x\leq 10\}
+
In the two modes on and off, the evolution of the temperature can be modelled by two different ODEs. Either from first-principles modelling or from system identification (or preferrably from the combination of the two) we get the two differential equations, say:
+f_\text{off}(x) = -0.1x,\quad f_\text{on}(x) = -0.1x + 5,
+ which gives another component for the hybrid automaton.
+
The control logic of the thermostat is captured by the \mathcal I and \mathcal G components of the hybrid automaton. Let’s determine them now. Obviously, if we just set 18 as the threshold, the heater would be switching on and off all the time. We need to introduce some hysteresis. Say, keeping the temperature within the interval (18 \pm 2)^\circ is acceptable.
+\mathcal I(\text{off}) = \{x\mid x> 16\},\quad \mathcal I(\text{on}) = \{x\mid x< 20\},
+
Finally, \mathcal R (or r) is not specified as the x variable (the temperature) doesn’t jump. Well, it is specified implicitly as an identity mapping r(x)=x.
+
The graphical representation of the thermostat hybrid automaton is shown in Fig 1.
+
+
+
+
+
+
+Figure 1: Hybrid automaton for a thermostat
+
+
+
+
Is this model deterministic? There are actually two reasons why it is not:
+
+
If we regard the characterization of the initial state (the temperature in this case) as a part of the model, which is the convention that we adhere to in our course, the model is nondeterministic.
+
Since the invariant for a given mode and the guard set for the only transition to the other model overlap, the response of the system is not uniquely determined. Consider the case when the system is in the off mode and the temperature is 16.5. The system can either stay in the off mode or switch to the on mode.
+
+
+
+
Hybrid automaton with external events and control inputs
+
We now extend the hybrid automaton with two new components:
+
+
a set \mathcal{A} of (external) events (also actions or symbols),
+
a set \mathcal{U} external continuous-valued inputs (control inputs or disturbances).
+
+
\boxed{
+ \mathcal{H} = \{\mathcal Q, \mathcal Q_0, \mathcal X, \mathcal X_0, \mathcal I, \mathcal A, \mathcal U, f, \mathcal E, \mathcal G, \mathcal R\} ,}
+ where
+
+
\mathcal A = \{a_1,a_2,\ldots, a_s\} is a set of events
+
+
The role identical as in a (finite) state automaton: an external event triggers an (enabled) transition from the current discrete state (mode, location) to another.
+
Unlike in pure discrete-event systems, here they are considered within a model that does recognize passing of time – each action must be “time-stamped”.
+
In simulations such timed event can be represented by an edge in the signal. In this regard, it might be tempting not to introduce it as a seperate entity, but it is useful to do so.
+
+
\mathcal U\in\mathbb R^m is a set of continuous-valued inputs
+
+
Real-valued functions of time.
+
Control inputs, disturbances, references, noises. In applications it will certainly be useful to distinghuish these roles, but here we keep just a single type of such an external variable, we do not have to distinguish.
+
+
+
+
Some modifications needed
+
Upon introduction of these two types of external inputs we must modify the components of the definition we provided earlier:
+
+
f: \mathcal Q \times \mathcal X \times \mathcal U \rightarrow \mathbb R^n is a vector field that now depends not only on the location but also on the external (control) input, that is, at a given location we consider the state equation \dot x = f_q(x,u).
+
\mathcal E\subseteq \mathcal Q \times (\mathcal A) \times \mathcal Q is a set of transitions now possibly parameterized by the actions (as in classical automata).
+
\mathcal I : \mathcal Q \rightarrow 2^{\mathcal{X}\times \mathcal U} is a location invariant now augmented with a subset of the control input set. The necessary condition for staying in the given mode can be thus imposed not only on x but also on u.
+
\mathcal G: \mathcal E \rightarrow 2^{\mathcal{X}\times U} is a guard set now augmented with a subset of the control input set. The necessary condition for a given transition can be thus imposed not only on x but also on u.
+
\mathcal R: \mathcal E \times \mathcal X\times \mathcal U\rightarrow 2^{\mathcal X} is a (state) reset map that is now additionally parameterized by the control input.
+
+
If enabled, the transition can happen if one of the two things is satisfied:
+
+
the continous state leaves the invariant set of the given location,
+
+
an external event occurs.
+
+
+
Example 2 (Button-controlled LED)
+
+
+
+
+
+
+Figure 2: Automaton for a button controlled LED
+
+
+
that is, x^+ = r((\mathrm{off},\mathrm{press},\mathrm{dim}),x) = 0.
+
For all other transitions r((\cdot, \cdot, \cdot),x)=x,
+
+
that is, x^+ = x.
+
+
+
+
+
Example 3 (Water tank) We consider a water tank with one inflow and two outflows – one at the bottom, the other at some nonzero height h_\mathrm{m}. The water level h is the continuous state variable.
+
+
+
+
+
+
+Figure 3: Water tank example
+
+
+
+
The model essentially expresses that the change in the volume is given by the difference between the inflow and the outflows. The outflows are proportional to the square root of the water level (Torricelli’s law)
+\dot V =
+\begin{cases}
+Q_\mathrm{in} - Q_\mathrm{out,middle} - Q_\mathrm{out,bottom}, & h>h_\mathrm{m}\\
+Q_\mathrm{in} - Q_\mathrm{out,bottom}, & h\leq h_\mathrm{m}
+\end{cases}
+
+
Apparently things change when the water level crosses (in any direction) the height h_\mathrm{m}. This can be modelled using a hybrid automaton.
+
+
+
+
+
+
+Figure 4: Automaton for a water tank example
+
+
+
+
One lesson to learn from this example is that the transition from one mode to another is not necessarily due to some computer-controlled switch. Instead, it is our modelling choice. It is an approximation that assumes negligible diameter of the middle pipe. But taking into the consideration the volume of the tank, it is probably a justifiable approximation.
+
+
+
Example 4 (Bouncing ball) We assume that a ball is falling from some initial nonzero height above the table. After hitting the table, it bounces back, loosing a portion of the energy (the deformation is not perfectly elastic).
+
+
+
+
+
+
+Figure 5: Bouncing ball example
+
+
+
+
The state equation during the free fall is
+\dot{\bm x} = \begin{bmatrix} x_2\\ -g\end{bmatrix}, \quad \bm x = \begin{bmatrix}10\\0\end{bmatrix}.
+
+
But how can we model what happens during and after the collision? High-fidelity model would be complicated, involving partial differential equations to model the deformation of the ball and the table. These complexities can be avoided with a simpler model assuming that immediately after the collision the sign of the velocity abruptly (discontinuously) changes, and at the same time the ball also looses a portion of the energy.
+
When modelling this using a hybrid automaton, it turns out that we only need a single discrete state. The crucial feature of the model is then the nontrivial (non-identity) reset map. This is depicted in Fig 6.
+
+
+
+
+
+
+Figure 6: Hybrid automaton for a bouncing ball eaxample
+
+
+
+
For completeness, here are the individual components of the hybrid automaton:
+\mathcal{Q} = \{q\}, \; \mathcal{Q}_0 = \{q\}
+
+r((q,q),\bm x) = \begin{bmatrix}x_1\\ -\gamma x_2 \end{bmatrix},
+ where \gamma is the coefficient of restitution (e.g., \gamma = 0.9).
+
+
+
+
+
+
+Comment on the invariant set for the bouncing ball
+
+
+
+
Some authors characterize the invariant set as x_1\geq 0. But this means that as the ball touches the ground, nothing forces it to leave the location and do the transition. Instead, the ball must penetrate the ground, however tiny distance, in order to trigger the transition. The current definition avoids this.
+
+
+
+
+
+
+
+
+Another comment on the invariant set for the bouncing ball
+
+
+
+
While the previous remark certainly holds, when solving the model numerically, the use of inequalities to define sets is inevitable. And some numerical solvers, in particular optimization solvers, cannot handle strict inequalities. That is perhaps why some authors are quite relaxed about this issue. We will encounter it later on.
+
+
+
+
+
Example 5 (Stick-slip friction model (Karnopp)) Consider a block of mass m placed freely on a surface. External horizontal force F_\mathrm{a} is applied to the block, setting it to a horizontaly sliding motion, against which the friction force F_\mathrm{f} is acting:
+m\dot v = F_\mathrm{a} - F_\mathrm{f}(v).
+
+
Common choice for a model of friction between two surfaces is Coulomb friction
+F_\mathrm{f}(v) = F_\mathrm{c}\operatorname*{sgn}(v).
+
+
The model is perfectly intuitive, isn’t it? Well, what if v=0 and F_\mathrm{a}<F_\mathrm{c}? Can you see the trouble?
+
One of the remedies is the Karnopp model of friction
+m\dot v = 0, \qquad v=0, \; |F_\mathrm{a}| < F_\mathrm{c}
+
+F_\mathrm{f} = \begin{cases}\operatorname*{sat}(F_\mathrm{a},F_\mathrm{c}), & v=0\\F_\mathrm{c}\operatorname*{sgn}(v), & \mathrm{else}\end{cases}
+
+
The model can be formulated as a hybrid automaton with two discrete states (modes, locations) as in Fig 7.
+
+
+
+
+
+
+Figure 7: Hybrid automaton for the Karnopp model of friction
+
+
+
+
+
+
Example 6 (Rimless wheel) A simple mechanical model that is occasionally used in the walking robot community is the rimless wheel rolling down a declined plane as depicted in Fig 8.
+
+
+
+
+
+
+Figure 8: Rimless wheel
+
+
+
+
A hybrid automaton for the rimless wheel is below.
+
+
+
+
+
+
+Figure 9: Hybrid automaton for a rimless wheel
+
+
+
+
Alternatively, we do not represent the discrete state graphically as a node in the graph but rather as another state variable s \in \{0, 1, \ldots, 5\} within a single location.
+
+
+
+
+
+
+Figure 10: Alternative hybrid automaton for a rimless wheel
+
+
+
+
+
+
Example 7 (DC-DC boost converter) The enabling mechanism for a DC-DC converter is switching. Although the switching is realized with a semiconductor switch, for simplicity of the exposition we consider a manual switch in Fig 11 below.
+
+
+
+
+
+
+Figure 11: DC-DC boost converter
+
+
+
+
The switch introduces two modes of operation. But the (ideal) diode introduces a mode transition too.
Possibly the events of opening and closing the switch can be driven by time: opening the switch is derived from the value of an input signal, closing the switch is periodic.
+
+
+
+
+
+
+
+Figure 15: Hybrid automaton for a DC-DC boost converter
+
+
+
HybridSystems.jl – the package essentially just defines some basic data types, currently it is mainly used by some other packages, such as ReachabilityAnalysis.jl.
Here we introduce a major alternative framework to hybrid automata for modelling hybrid systems. It is called hybrid equations, sometimes also hybrid state equations to emphasize that what we are after is some kind of analogy with state equations \dot{x}(t) = f(x(t), u(t)) and x_{k+1} = g(x_k, u_k) that we are familiar with from (continuous-valued) dynamical systems. Sometimes it also called event-flow equations or jump-flow equations.
+
These are the key ideas:
+
+
The (state) variables can change values discontinuously upon occurence of events – they jump.
+
Between the jumps they evolve continuously – they flow.
+
Some variables may only flow, they never jump.
+
The variables staying constant between the jumps can be viewed as flowing too.
+
+
The major advantage of this modeling framework is that we do not have to distinguish between the two types of state variables. This is in contrast with hybrid automata, where we have to start by classifying the state variables as either continuous or discrete before moving on. In the current framework we treat all the variables identically – they mostly flow and occasionally (perhaps never, which is OK) jump.
+
+
Hybrid equations
+
It is high time to introduce hybrid (state) equations – here they come
+\begin{aligned}
+\dot{x} &= f(x), \quad x \in \mathcal{C},\\
+x^+ &= g(x), \quad x \in \mathcal{D},
+\end{aligned}
+ where
+
+
f is the flow map,
+
\mathcal{C} is the flow set,
+
g is the jump map,
+
\mathcal{D} is the jump set.
+
+
This model of a hybrid system is thus parameterized by the quadruple \{f, \mathcal C, g, \mathcal D\}.
+
+
+
Hybrid inclusions
+
We now extend the presented framework of hybrid equations a bit. Namely, the functions on the right-hand sides in both the differential and the difference equations are no longer assigning just a single value (as well-behaved functions do), but they assign sets!
+\begin{aligned}
+\dot{x} &\in \mathcal F(x), \quad x \in \mathcal{C},\\
+x^+ &\in \mathcal G(x), \quad x \in \mathcal{D}.
+\end{aligned}
+ where
+
+
\mathcal{F} is the set-valued flow map,
+
\mathcal{C} is the flow set,
+
\mathcal{G} is the set-valued jump map,
+
\mathcal{D} is the jump set.
+
+
+
+
Output equations
+
Typically a full model is only formed upon defining some output variables (oftentimes just a subset of possibly scaled state variables or their linear combinations). These output variables then obey some output equation
+y(t) = h(x(t)),
+
+
or
+y(t) = h(x(t),u(t)).
+
+
+
Example 1 (Bouncing ball) This is the “hello world example” for hybrid systems with state jumps (pun intended). The state variables are the height and the vertical speed of the ball.
+\bm x \in \mathbb{R}^2, \qquad \bm x = \begin{bmatrix}x_1 \\ x_2\end{bmatrix}.
+
+
The quadruple defining the hybrid equations is
+\mathcal{C} = \{\bm x \in \mathbb{R}^2 \mid x_1>0 \lor (x_1 = 0, x_2\geq 0)\},
+
+f(\bm x) = \begin{bmatrix}x_2 \\ -g\end{bmatrix}, \qquad g = 9.81,
+
+\mathcal{D} = \{\bm x \in \mathbb{R}^2 \mid x_1 = 0, x_2 < 0\},
+
+g(\bm x) = \begin{bmatrix}x_1 \\ -\alpha x_2\end{bmatrix}, \qquad \alpha = 0.8.
+
+
The two sets and two maps are illustrated below.
+
+
+
+
+
+
+Figure 1: Maps and sets for the bouncing ball example
+
+
+
+
+
+
Example 2 (Bouncing ball on a controlled piston) We now extend the simple bouncing ball example by adding a vertically moving piston. The piston is controlled by a force.
+
+
+
+Example of a ball bouncing on a vertically moving piston
+
+
+
In our analysis we neglect the sizes (for simplicity).
+
The collision happens when x_\mathrm{b} = x_\mathrm{p}, and v_\mathrm{b} < v_\mathrm{p}.
+
The conservation of momentum after a collision reads
+m_\mathrm{b}v_\mathrm{b}^+ + m_\mathrm{p}v_\mathrm{p}^+ = m_\mathrm{b}v_\mathrm{b} + m_\mathrm{p}v_\mathrm{p}.
+\tag{1}
+
The collision is modelled using a restitution coefficient
+v_\mathrm{p}^+ - v_\mathrm{b}^+ = -\gamma (v_\mathrm{p} - v_\mathrm{b}).
+\tag{2}
+
From the momentum conservation Equation 1
+v_\mathrm{p}^+ = \frac{m_\mathrm{b}}{m_\mathrm{p}}v_\mathrm{b} + v_\mathrm{p} - \frac{m_\mathrm{b}}{m_\mathrm{p}}v_\mathrm{b}^+
+
+
we substitute to Equation 2 to get
+\frac{m_\mathrm{b}}{m_\mathrm{p}}v_\mathrm{b} + v_\mathrm{p} - \frac{m_\mathrm{b}}{m_\mathrm{p}}v_\mathrm{b}^+ - v_\mathrm{b}^+ = -\gamma (v_\mathrm{p} - v_\mathrm{b}),
+ from which we express v_\mathbb{b}^+
+\begin{aligned}
+v_\mathrm{b}^+ &= \frac{1}{1+\frac{m_\mathrm{b}}{m_\mathrm{p}}}\left(\frac{m_\mathrm{b}}{m_\mathrm{p}}v_\mathrm{b} + v_\mathrm{p} + \gamma (v_\mathrm{p} - v_\mathrm{b})\right)\\
+&= \frac{m_\mathrm{p}}{m_\mathrm{p}+m_\mathrm{b}}\left(\frac{m_\mathrm{b}-\gamma m_\mathrm{p}}{m_\mathrm{p}}v_\mathrm{b} + (1+\gamma)v_\mathrm{p}\right)\\
+&= \frac{m_\mathrm{b}-\gamma m_\mathrm{p}}{m_\mathrm{b}+m_\mathrm{p}}v_\mathrm{b} + \frac{(1+\gamma)m_\mathrm{p}}{m_\mathrm{p}+m_\mathrm{b}}v_\mathrm{p}
+\end{aligned}.
+
+
Substitute to the expression for v_\mathbb{p}^+ to get
+\begin{aligned}
+v_\mathrm{p}^+ &= \frac{m_\mathrm{b}}{m_\mathrm{p}}v_\mathrm{b} + v_\mathrm{p} - \frac{m_\mathrm{b}}{m_\mathrm{p}}\left(\frac{m_\mathrm{b}-\gamma m_\mathrm{p}}{m_\mathrm{b}+m_\mathrm{p}}v_\mathrm{b} + \frac{(1+\gamma)m_\mathrm{p}}{m_\mathrm{p}+m_\mathrm{b}}v_\mathrm{p}\right)\\
+&= \frac{m_\mathrm{b}}{m_\mathrm{p}}\left(1-\frac{m_\mathrm{b}-\gamma m_\mathrm{p}}{m_\mathrm{b}+m_\mathrm{p}}\right) v_\mathrm{b} \\
+&\qquad\qquad + \left(1-\frac{m_\mathrm{b}}{m_\mathrm{p}}\frac{(1+\gamma)m_\mathrm{p}}{m_\mathrm{p}+m_\mathrm{b}}\right) v_\mathrm{p}\\
+&= \frac{m_\mathrm{b}}{m_\mathrm{b}+m_\mathrm{p}}(1+\gamma) v_\mathrm{b} + \frac{m_\mathrm{p}-\gamma m_\mathrm{b}}{m_\mathrm{p}+m_\mathrm{b}} v_\mathrm{p}.
+\end{aligned}
+
+
Finally we can simplify the expressions a bit by introducing m=\frac{m_\mathrm{b}}{m_\mathrm{b}+m_\mathrm{p}}. The jump equation is then
+\begin{bmatrix}
+v_\mathrm{b}^+\\
+v_\mathrm{p}^+
+\end{bmatrix}
+=
+\begin{bmatrix}
+m - \gamma (1-m) & (1+\gamma)(1-m)\\
+m(1+\gamma) & 1-m-\gamma m
+\end{bmatrix}
+\begin{bmatrix}
+v_\mathrm{b}\\
+v_\mathrm{p}
+\end{bmatrix}.
+
+
+
+
Example 3 (Synchronization of fireflies) This is a famous example in synchronization. We consider n fireflies, x_i is the i-th firefly’s clock, normalized to [0,1]. The clock resets (to zero) when it reaches 1. Each firefly can see the flashing of all other fireflies. As soon as it observes a flash, it increases its clock by \varepsilon \%.
+
Here is how we model the problem using the four-tuple \{f, \mathcal C, g, \mathcal D\}:
+\mathcal{C} = [0,1)^n = \{\bm x \in \mathbb R^n\mid x_i \in [0,1),\; i=1,\ldots,n \},
+
+\bm f = [f_1, f_2, \ldots, f_n]^\top,\quad f_i = 1, \quad i=1,\ldots,n,
+
+\mathcal{D} = \{\bm x \in [0,1]^n \mid \max_i x_i = 1 \},
+
The flow set is
+\begin{aligned}
+\mathcal{C} &= \{\bm x \mid q=0,\, \tau<\frac{\alpha}{\omega},\, i_\mathrm{L}=0\}\\ &\qquad \cup \{\bm x \mid q=1,\, i_\mathrm{L}>0\}
+\end{aligned}.
+
+
The jump set is
+\begin{aligned}
+\mathcal{D} &= \{\bm x \mid q=0,\, \tau\geq \frac{\alpha}{\omega},\, i_\mathrm{L}=0,\, v_\mathrm{C}>0\}\\ &\qquad \cup \{\bm x \mid q=1,\, i_\mathrm{L}=0,\, v_\mathrm{C}<0\}
+\end{aligned}.
+
The last condition in the jump set comes from the requirement that not only must the current through the inductor be zero, but also it must be decreasing. And from the state equation it follows that the voltage on the capacitor must be negative.
+
+
+
Example 5 (Sampled-data feedback control) Another example of a dynamical system that fits nicely into the hybrid equations framework is sampled-data feedback control system. Within the feedback loop in Figure 3, we recognize a continuous-time plant and a discrete-time controller.
+
+
+
+
+
+
+Figure 3: Sampled data feedback control
+
+
+
+
The plant is modelled by \dot x_\mathrm{p} = f_\mathrm{p}(x_\mathrm{p},u), \; y = h(x_\mathrm{p}). The controller samples the output T-periodically and computes its own output as a nonlinear function u = \kappa(r-y).
+
The closed-loop model is then
+\dot x_\mathrm{p} = f_\mathrm{p}(x_\mathrm{p},\kappa(r-h(x_\mathrm{p}))), \; y = h(x_\mathrm{p}).
+
+
The closed-loop state vector is
+\bm x =
+\begin{bmatrix}
+x_\mathrm{p}\\ u \\ \tau
+\end{bmatrix}
+\in
+\mathbb R^n \times \mathbb R^m \times \mathbb R.
+
+
The flow set is
+\begin{aligned}
+\mathcal{C} &= \{\bm x \mid \tau \in [0,T)\}
+\end{aligned}
+
+
The flow map is
+\bm f(\bm x)
+=
+\begin{bmatrix}
+f_\mathrm{p}(x_\mathrm{p},u)\\
+0\\
+1
+\end{bmatrix}
+
+
The jump set is
+\begin{aligned}
+\mathcal{D} &= \{\bm x \mid \tau = T\}
+\end{aligned}
+ or rather
+\begin{aligned}
+\mathcal{D} &= \{\bm x \mid \tau \geq T\}
+\end{aligned}
+
+
The jump map is
+\bm g(\bm x) =
+\begin{bmatrix}
+x_\mathrm{p}\\
+\kappa(r-y)\\
+0
+\end{bmatrix}
+
+
You may wonder why we bother with modelling this system as a hybrid system at all. When it comes to analysis of the closed-loop system, implementation of the model in Simulink allows for seemless mixing of continuous-time and dicrete-time blocks. And when it comes to control design, we can either discretize the plant and design a discrete-time controller, or design a continuous-time controller and then discretize it. No need for new theoris. True, but still, it is nice to have a rigorous framework for analysis of such systems. The more so that the sampling time T may not be constant – it can either vary randomly or perhaps the samling can be event-triggered. All these scenarios are easily handled within the hybrid equations framework.
+
+
+
+
Hybridness after closing the loop
+
We have defined hybrid systems, but what exactly is hybrid when we close a feedback loop? There are three possibilities:
+
+
Hybrid plant + continuous controller.
+
Hybrid plant + hybrid controller.
+
Continuous plant + hybrid controller.
+
+
The first case is encountered when we use a standard controller such as a PID controller to control a system whose dynamics can be characterized/modelled as hybrid. The second scenario considers a controller that mimicks the behavior of a hybrid system. The third case is perhaps the least intuitive: although the plant to be controller is continuous(-valued), it may still make sense to design and implement a hybrid controller, see the next paragraph.
+
+
+
Impossibility to stabilize without a hybrid controller
+
+
Example 6 (Unicycle stabilization) We consider a unicycle model of a vehicle in a plane, characterized by the position and orientation, with controlled forward speed and yaw (turning) angular rate.
+
+
+
+
+
+
+Figure 4: Unicycle vehicle
+
+
+
+
The vehicle is modelled by
+\begin{aligned}
+\dot x &= v \cos \theta,\\
+\dot y &= v \sin \theta,\\
+\dot \theta &= \omega,
+\end{aligned}
+
+
+\bm x = \begin{bmatrix}
+x\\ y\\ \theta
+\end{bmatrix},
+\quad
+\bm u = \begin{bmatrix}
+v\\ \omega
+\end{bmatrix}.
+
+
It is known that this system cannot be stabilized by a continuous feedback controller. The general result that applies here was published in [1]. The condition of stabilizability by a time-invariant continuous state feedback is that the image of every neighborhood of the origin under (\bm x,\bm u) \mapsto \bm f(\bm x, \bm u) contains some neighborhood of the origin. This is not the case here. The map from the state-control space to the velocity space is
Now consider a neighborhood of the origin such that |\theta|<\frac{\pi}{2}. It is impossible to get \bm f(\bm x, \bm u) = \begin{bmatrix}
+0\\ f_2 \\ 0\end{bmatrix}, \; f_2\neq 0. Hence, stabilization by a continuous feedback \bm u = \kappa (\bm x) is impossible.
+
But it is possible to stabilize the vehicle using a discontinuous feedback. And discontinuous feedback controllr can be viewed as switching control, which in turn can be seen as instance of a hybrid controller.
+
+
+
Example 7 (Global stabilization on a circle) We now demonstration of a general phenomenon of stabilization on a manifold. We will see that even if stabilization by a continous feedback is possible, it might not guarantee global stability.
+
The dynamics of a particle sliding around a unit circle \mathbb S_1 is modelled by
+
+\dot{\bm x} = u\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix}\bm x,
+ where \bm x \in \mathbb S^1,\quad u\in \mathbb R.
+
The point to be stabilized is \bm x^* = \begin{bmatrix}1\\ 0\end{bmatrix}.
+
+
+
+
+
+
+Figure 5: Stabilization on a circle
+
+
+
+
What is required from a globally stabilizing controller?
+
+
Solutions stay in \mathbb S^1,
+
Solutions converge to \bm x^*,
+
If a solution starts near \bm x^*, it stays near.
+
+
One candidate is \kappa(\bm x) = -x_2.
+
Define the (Lyapunov) function V(\bm x) = 1-x_1.
+
Its time derivative along the solution trajectory
+\begin{aligned}
+\dot V &= \left(\nabla_{\bm{x}}V\right)^\top \dot{\bm x}\\
+&= \begin{bmatrix}-1 & 0\end{bmatrix}\left(-x_2\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\right)\\
+&= -x_2^2\\
+&= -(1-x_1^2).
+\end{aligned}
+
+
It follows that
+\dot V < 0 \quad \forall \bm x \in \mathbb S^1 \setminus \{-1,1\}.
+
+
With u=-x_2 the point \bm x^* is stable but not globally atractive.
+
Can we do better? Yes, we can. But we need to incorporate some switching into the controller. Loosely speaking anywhere except for the state (-1,0), we can apply the previously design controller, and at the troublesome state (-1,0), or actually in some region around it, we need to switch to another controller that would drive the system away from the problematic region.
+
But we will take this example as an opportunity to go one step further and instead of just a switching controller we design a hybrid controller. The difference is that within a hybrid controller we can incorporate some hysteresis, which is a robustifying feature. In order to to that, we need to introduce a new state variable q\in\{0,1\}. Determination of the flow and jump sets is sketched in Figure 6.
+
+
+
+
+
+
+Figure 6: Definition of the sets defining a hybrid controller
+
+
+
+
The two feedback controllers are given by
+\begin{aligned}
+\kappa(\bm x,0) &= \kappa_0(\bm x) = -x_2,\\
+\kappa(\bm x,1) &= \kappa_1(\bm x) = -x_1.
+\end{aligned}
+
+
The flow map is (DIY)
+f(\bm x, q) = \ldots
+
+
The flow set is
+\mathcal{C} = (\mathcal C_0 \times \{0\}) \cup (\mathcal C_1 \times \{1\}).
+
+
The jump set is
+\mathcal{D} = (\mathcal D_0 \times \{0\}) \cup (\mathcal D_1 \times \{1\}).
+
+
The jump map is
+g(\bm x, q) = 1-q \quad \forall [\bm x, q]^\top \in \mathcal D.
+
+
Simulation using Julia is provided below.
+
+
+Show the code
+
usingOrdinaryDiffEq
+
+# Defining the sets and functions for the hybrid equations
+
+c₀, c₁ =-2/3, -1/3
+
+C(x,q) = (x[1] >= c₀ && q ==0) || (x[1] <= c₁ && q ==1) # Actually not really needed, just a complement of D.
+D(x,q) = (x[1] < c₀ && q ==0) || (x[1] > c₁ && q ==1)
+
+g(x,q) =1-q
+
+κ(x,q) = q==0 ? -x[2] :-x[1]
+
+functionf!(dx,x,q,t) # Already in the format for the ODE solver.
+ A = [0.0-1.0; 1.00.0]
+ dx .=A*x*κ(x,q)
+end
+
+# Defining the initial conditions for the simulation
+
+cᵢ = (c₀+c₁)/2
+x₀ = [cᵢ,sqrt(1-cᵢ^2)]
+q₀ =1
+
+# Setting up the simulation problem
+
+tspan = (0.0,10.0)
+prob =ODEProblem(f!,x₀,tspan,q₀)
+
+functioncondition(x,t,integrator)
+ q = integrator.p
+returnD(x,q)
+end
+
+functionaffect!(integrator)
+ q = integrator.p
+ x = integrator.u
+ integrator.p =g(x,q)
+end
+
+cb =DiscreteCallback(condition,affect!)
+
+# Solving the simulation problem
+
+sol =solve(prob,Tsit5(),callback=cb,dtmax=0.1) # ContinuousCallback more suitable here
+
+# Plotting the results of the simulation
+
+usingPlots
+gr(tickfontsize=12,legend_font_pointsize=12,guidefontsize=12)
+
+plot(sol,label=["x₁""x₂"],xaxis="t",yaxis="x",lw=2)
+hline!([c₀], label="c₀")
+hline!([c₁], label="c₁")
+
+
+
+
+
+
+
+
+Figure 7: Simulation of stabilization on a circle using a hybrid controller
+
+
+
+
+
The solution can also be visualized in the state space.
+
+Figure 8: Simulation of stabilization on a circle using a hybrid controller
+
+
+
+
+
+
+
+
Supervisory control
+
Yet another problem that can benefit from being formulated as a hybrid system is supervisory control.
+
+
+
+
+
+
+Figure 9: Supervisory control
+
+
+
+
+
+
Combining local and global controllers (\subset supervisory control)
+
As a subset of supervisory control we can view a controller that switches between a global and a local controller.
+
+
+
+
+
+
+Figure 10: Combining global and local controllers
+
+
+
+
Local controllers have good transient response but only work well in a small region around the equilibrium state. Global controllers have poor transient response but work well in a larger region around the equilibrium state.
+
A useful example is that of swinging up and stabilization of a pendulum: the local controller can be designer for a linear model obtained by linearization about the upright orientation of the pendulum. But such controller can only be expected to perform well in some small region around the upright orientation. The global controller is designed to bring the pendulum into that small region.
+
The flow and jump sets for the local and global controllers are in Figure 11. Can you tell, which is which? Remember that by introducing the discrete variable q, some hysteresis is in the game here.
+
+
+
+
+
+
+Figure 11: Flow and jump sets for a local and a global controller
+
+
+
R. Brockett, “Asymptotic stability and feedback stabilization,” in Differential Geometric Control Theory, R. Brockett, R. Millman, and H. Sussmann, Eds., Boston: Birkhäuser, 1983.
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/hybrid_equations.html b/hybrid_equations.html
index 247dd3f..87a8183 100644
--- a/hybrid_equations.html
+++ b/hybrid_equations.html
@@ -1162,57 +1162,57 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Figure 7: Simulation of stabilization on a circle using a hybrid controller
@@ -1236,54 +1236,54 @@
The theoretical and computational framework of hybrid equations has been mostly developed by a relatively small circle of researchers (Sanfelice, Goebel, Teel, …). The primary monograph is [1]. It is also supported by a freely available Matlab toolbox, see the section on software.
+
+
The book [2] can be regarded as a predecessor and/or complement of the just mentioned [1]. Although the book is not available online, a short version appears as an article [3] in the popular IEEE Control Systems magazine (the one with color figures :-).
R. Goebel, R. G. Sanfelice, and A. R. Teel, “Hybrid dynamical systems,”IEEE Control Systems Magazine, vol. 29, no. 2, pp. 28–93, Apr. 2009, doi: 10.1109/MCS.2008.931718.
Max-plus algebra, also written as (max,+) algebra (and also known as tropical algebra/geometry and dioid algebra), is an algebraic framework in which we can model and analyze a class of discrete-event systems, namely event graphs, which we have previously introduced as a subset of Petri nets. The framework is appealing in that the models than look like state equations \bm x_{k+1} = \mathbf A \bm x_k + \mathbf B \bm u_k for classical linear dynamical systems. We call these max-plus linear systems, or just MPL systems. Concepts such as poles, stability and observability can be defined, following closely the standard definitions. In fact, we can even formulate control problems for these models in a way that mimicks the conventional control theory for LTI systems, including MPC control.
+
But before we get to these applications, we first need to introduce the (max,+) algebra itself. And before we do that, we recapitulate the definition of a standard algebra.
+
+
+
+
+
+
+Algebra is not only a branch of mathematics
+
+
+
+
Confusingly enough, algebra is used both as a name of both a branch of mathematics and a special mathematical structure. In what follows, we use the term algebra to refer to the latter.
+
+
+
+
Algebra
+
Algebra is a set of elements equipped with
+
+
two operations:
+
+
addition (plus, +),
+
multiplication (times, ×),
+
+
neutral (identity) element with respect to addition: zero, 0, a+0=a,
+
neutral (identity) element with respect to multiplication: one, 1, a\times 1 = a.
+
+
Inverse elements can also be defined, namely
+
+
Inverse element wrt addition: -aa+(-a) = 0.
+
Inverse element wrt multiplication (except for 0): a^{-1}a \times a^{-1} = 1.
+
+
If the inverse wrt multiplication exists for every (nonzero) element, the algebra is called a field, otherwise it is called a ring.
+
Prominent examples of a ring are integers and polynomials. For integers, it is only the numbers 1 and -1 that have integer inverses. For polynomials, it is only zero-th degree polynomials that have inverses qualifying as polynomials too. An example from the control theory is the ring of proper stable transfer functions, in which only the non-minimum phase transfer functions with zero relative degree have inverses, and thus qualify as units.
+
Prominent example of a field is the set of real numbers.
+
+
+
(max,+) algebra: redefining the addition and multiplication
+
Elements of the (max,+) algebra are real numbers, but it is still a ring and not a field since the two operations are defined differently.
+
The new operations of addition, which we denote by \oplus to distinguish it from the standard addition, is defined as \boxed{
+ x\oplus y \equiv \max(x,y).
+ }
+
+
The new operation of multiplication, which we denote by \otimes to distinguish it from the standard multiplication, is defined as \boxed{
+ x\otimes y \equiv x+y}.
+
+
+
+
+
+
+
+Important
+
+
+
+
Indeed, there is no typo here, the standard addition is replaced by \otimes and not \oplus.
+
+
+
+
+
+
+
+
+(min,+) also possible
+
+
+
+
Indeed, we can also define the (min,+) algebra. But for our later purposes in modelling we prefer the (max,+) algebra.
+
+
+
+
Reals must be extended with the negative infinity
+
Strictly speaking, the (max,+) algebra is a broader set than just \mathbb R. We need to extend the reals with the minus infinity. We denote the extended set by \mathbb R_\varepsilon\boxed{
+\mathbb R_\varepsilon \coloneqq \mathbb R \cup \{-\infty\}}.
+
+
The reason for the notation is that a dedicated symbol \varepsilon is assigned to this minus infinity, that is, \boxed
+{\varepsilon \coloneqq -\infty.}
+
+
It may yield some later expressions less cluttered. Of course, at the cost of introducing one more symbol.
+
We are now going to see the reason for this extension.
+
+
+
Neutral elements
+
+
Neutral element with respect to \oplus
+
The neutral element with respect to \oplus, the zero, is -\infty. Indeed, for x \in \mathbb R_\varepsilon
+x \oplus \varepsilon = x,
+ because \max(x,-\infty) = x.
+
+
+
Neutral element with respect to \otimes
+
The neutral element with respect to \otimes, the one, is 0. Indeed, for x \in \mathbb R_\varepsilon
+x \otimes \varepsilon = x,
+ because x+0=x.
+
+
+
+
+
+
+Nonsymmetric notation, but who cares?
+
+
+
+
The notation is rather nonsymmetric here. We now have a dedicated symbol \varepsilon for the zero element in the new algebra, but no dedicated symbol for the one element in the new algebra. It may be a bit confusing as “the old 0 is the new 1”. Perhaps similarly as we introduced dedicated symbols for the new operations of addition of multiplication, we should have introduced dedicated symbols such as ⓪ and ①, which would lead to expressions such as x⊕⓪=x and x ⊗ ① = x. In fact, in some software packages they do define something like mp-zero and mp-one to represent the two special elements. But this is not what we will mostly encounter in the literature. Perhaps the best attitude is to come to terms with this notational asymetry… After all, I myself was apparently not even able to figure out how to encircle numbers in LaTeX…
+
+
+
+
+
+
Inverse elements
+
+
Inverse with respect to \oplus
+
The inverse element with respect to \oplus in general does not exist! Think about it for a few moments, this is not necessarily intuitive. For which element(s) does it exist? Only for \varepsilon.
+
This has major consequences, for example, x\oplus x=x.
+
Can you verify this statement? How is is related to the fact that the inverse element with respect to \oplus does not exist in general?
+
This is the key difference with respect to a conventional algebra, wherein the inverse element of a wrt conventional addition is -a, while here we do not even define \ominus.
+
Formally speaking, the (max,+) algebra is only a semi-ring.
+
+
+
Inverse with respect to \otimes
+
The inverse element with respect to \otimes does not exist for all elements. The \varepsilon element does not have an inverse element with respect to \otimes. But in this aspect the (max,+) algebra just follows the conventional algebra, beucase 0 has no inverse there either.
+
+
+
+
+
Powers and the inverse with respect to \otimes
+
Having defined the fundamental operations and the fundamental elements, we can proceed with other operations. Namely, we consider powers. Fot an integer r\in\mathbb Z, the rth power of x, denoted by x^{\otimes^r}, is defined, unsurprisingly as x^{\otimes^r} \coloneqq x\otimes x \otimes \ldots \otimes x.
+
Observe that it corresponds to rx in the conventional algebra x^{\otimes^r} = rx.
+
But then the inverse element with respect to \otimes can also be determined using the (-1)th power as x^{\otimes^{-1}} = -x.
+
This is not actually surprising, is it?
+
There are few more implications. For example,
+x^{\otimes^0} = 0.
+
There is also no inverse element with respect to \otimes for \varepsilon, but it is expected as \varepsilon is a zero wrt \oplus. Furthermore, for r\neq -1, if r>0 , then \varepsilon^{\otimes^r} = \varepsilon, if r<0 , then \varepsilon^{\otimes^r} is undefined, which are both expected. Finally, \varepsilon^{\otimes^0} = 0 by convention.
+
+
+
Order of evaluation of (max,+) formulas
+
It is the same as that for the conventional algebra:
+
+
power,
+
multiplication,
+
addition.
+
+
+
+
(max,+) polynomials (aka tropical polynomials)
+
Having covered addition, multiplication and powers, we can now define (max,+) polynomials. In order to get started, consider the the univariate polynomial p(x) = a_{n}\otimes x^{\otimes^{n}} \oplus a_{n-1}\otimes x^{\otimes^{n-1}} \oplus \ldots \oplus a_{1}\otimes x \oplus a_{0},
+ where a_i\in \mathbb R_\varepsilon and n\in \mathbb N.
+
By interpreting the operations, this translates to the following function \boxed
+{p(x) = \max\{nx + a_n, (n-1)x + a_{n-1}, \ldots, x+a_1, a_0\}.}
+
+
+
Example 1 (1D polynomial) Consider the following (max,+) polynomial
+p(x) = 2\otimes x^{\otimes^{2}} \oplus 3\otimes x \oplus 1.
+
+
We can interpret it in the conventional algebra as
+p(x) = \max\{2x+2,x+3,1\},
+ which is a piecewise linear (actually affine) function.
Example 2 (Example of a 2D polynomial) Nothing prevents us from defining a polynomial in two (and more) variables. For example, consider the following (max,+) polynomial
+p(x,y) = 0 \oplus x \oplus y.
+
Piecewise affine (PWA) functions will turn out a frequent buddy in our course.
+
+
+
+
+
Solution set (zero set)
+
…
+
+
+
Matrix computations
+
+
Addition and multiplication
+
What is attractive about the whole (max,+) framework is that it also extends nicely to matrices. For matrices, whose elements are in \mathbb R_\varepsilon, we define the operations of addition and multiplication identically as in the conventional case, we just use different definitions of the two basic scalar operations. (A\oplus B)_{ij} = a_{ij}\oplus b_{ij} = \max(a_{ij},b_{ij})
+\begin{aligned}
+(A\otimes B)_{ij} &= \bigoplus_{k=1}^n a_{ik}\otimes b_{kj}\\
+&= \max_{k=1,\ldots, n}(a_{ik}+b_{kj})
+\end{aligned}
+
+
+
+
Zero and identity matrices
+
(max,+) zero matrix \mathcal{E}_{m\times n} has all its elements equal to \varepsilon, that is,
+\mathcal{E}_{m\times n} =
+\begin{bmatrix}
+\varepsilon & \varepsilon & \ldots & \varepsilon\\
+\varepsilon & \varepsilon & \ldots & \varepsilon\\
+\vdots & \vdots & \ddots & \vdots\\
+\varepsilon & \varepsilon & \ldots & \varepsilon
+\end{bmatrix}.
+
+
(max,+) identity matrix I_n has 0 on the diagonal and \varepsilon elsewhere, that is,
+I_{n} =
+\begin{bmatrix}
+0 & \varepsilon & \ldots & \varepsilon\\
+\varepsilon & 0 & \ldots & \varepsilon\\
+\vdots & \vdots & \ddots & \vdots\\
+\varepsilon & \varepsilon & \ldots & 0
+\end{bmatrix}.
+
+
+
+
Matrix powers
+
The zerothe power of a matrix is – unsurprisingly – the identity matrix, that is, A^{\otimes^0} = I_n.
+
The kth power of a matrix, for k\in \mathbb N\setminus\{0\}, is then defined using A^{\otimes^k} = A\otimes A^{\otimes^{k-1}}.
+
+
+
+
Connection with graph theory – precedence graph
+
Consider A\in \mathbb R_\varepsilon^{n\times n}. For this matrix, we can define the precedence graph\mathcal{G}(A) as a weighted directed graph with the vertices 1, 2, …, n, and with the arcs (j,i) with the associated weights a_{ij} for all a_{ij}\neq \varepsilon. The kth power of the matrix is then
+
+(A)^{\otimes^k}_{ij} = \max_{i_1,\ldots,i_{k-1}\in \{1,2,\ldots,n\}} \{a_{ii_1} + a_{i_1i_2} + \ldots + a_{i_{k-1}j}\}
+ for all i,j and k\in \mathbb N\setminus 0.
+
+Figure 1: An example of a precedence graph
+
+
+
+
+
+
+
+
Irreducibility of a matrix
+
+
Matrix in \mathbb R_\varepsilon^{n\times n} is irreducible if its precedence graph is strongly connected.
+
Matrix is irreducible iff
+(A \oplus A^{\otimes^2} \oplus \ldots A^{\otimes^{n-1}})_{ij} \neq \varepsilon \quad \forall i,j, i\neq j.
+\tag{1}
+
+
+
+
+
Eigenvalues and eigenvectors
+
Eigenvalues and eigenvectors constitute another instance of a straightforward import of concepts from the conventional algebra into the (max,+) algebra – just take the standard definition of an eigenvalue-eigenvector pair and replace the conventional operations with the (max,+) alternatives
+A\otimes v = \lambda \otimes v.
+
+
A few comments:
+
+
In general, total number of (max,+) eigenvalues <n.
+
An irreducible matrix has only one (max,+) eigenvalue.
+
Graph-theoretic interpretation: maximum average weight over all elementary circuits…
+
+
+
Eigenvalue-related property of irreducible matrices
+
For large enough k and c, it holds that \boxed
+{A^{\otimes^{k+c}} = \lambda^{\otimes^c}\otimes A^{\otimes^k}.}
+
+
+
+
+
Solving (max,+) linear equations
+
We can also define and solve linear equations within the (max,+) algebra. Considering A\in \mathbb R_\varepsilon^{n\times n},\, b\in \mathbb R_\varepsilon^n, we can formulate and solve the equation
+A\otimes x = b.
+
+
In general no solution even if A is square. However, often we can find some use for a subsolution defined as
+A\otimes x \leq b.
+
+
Typically we search for the maximal subsolution instead, or subsolutions optimal in some other sense.
+A \otimes x =
+\begin{bmatrix}
+1\\ 0 \\ 1
+\end{bmatrix}
+\leq b
+
+
+
With this introduction to the (max,+) algebra, we are now ready to move on to the modeling of discrete-event systems using the max-plus linear (MPL) systems.
Since the new modelling framework is expected to be useful for prediction of a system response within model predictive control, it must model a hybrid system in discrete time. This is a major difference from what we did in our course so far.
+
In particular, we are going to model a hybrid system as a discrete(-time) hybrid automaton (DHA), which means that
+
+
the continuous-value dynamics (often corresponding to the physical part of the system) evolves in discrete time,
+
the events and their processing by the logical part of the system are synchronized to the same periodic clock.
+
+
+
Four components of a discrete(-time) hybrid automaton
+
We are already well familiar with the concept of a hybrid automaton, and the restriction to discrete time does not seem to warrant reopening the definition (modes/locations, guards, invariants/domains, reset maps, …). However, it turns out that reformulating/restructuring the hybrid automaton will be useful for our ultimate goal of developing an MPC-friendly modelling framework. In particular, we consider four components of a DHA:
+
+
switched affine system (SAS),
+
mode selector (MS),
+
event generator (EG),
+
finite state machine (FSM).
+
+
Their interconnection is shown in the following figure.
+
+
Draw the block diagram from Bemporad’s materials (book, toolbox documentation).
+
+
Let’s discuss the individual components (and while doing that, you can think about the equivalent concept in the classical definition of a hybrid automaton such as mode, invariant, guard, …).
+
+
Switched affine systems (SAS)
+
This is a model of the continuous-value dynamics parameterized by the index i that evolves in (discrete) time
+\begin{aligned}
+x_c(k+1) &= A_{i(k)} x_c(k) + B_{i(k)} u_c(k) + f_{i(k)}\\
+y_c(k) &= C_{i(k)} x_c(k) + D_{i(k)} u_c(k) + g_{i(k)}
+\end{aligned}
+
+
In principle there is no need to restrict the right hand sides to affine functions as we did, but the fact is that the methods and tools are currently only available for this restricted class of systems.
+
+
+
Event generator (EG)
+
We consider partitioning of the state space or possibly state-control space into polyhedral regions. The system is then in the ith region of the state-input space, if the continuous-value state x_c(k) and the continuous-value control input u_c satisfy
+H_i x_c(k) + J_i u_c(k) + K_i \leq 0
+
+
The event indicated by the (vector) binary variable
+\delta_e(k) = h(x_c(k), u_c(k)) \in \{0,1\}^m,
+
Our goal now is to turn the IF-THEN conditions in the model into linear inequalities. This will allow us to formulate the model as a mathematical program, actually a mixed-integer program (MIP).
+
+
Propositional logic and connectives
+
Propositions that are either true or false are composed of elementary propositions (Boolean variables) and connectives.
+
+
Boolean variable (or elementary proposition)
+
\(X\) evaluates to true or false. Oftentimes values 0 and 1 are used, but it should be clear that these are not numbers but logical values.
+
+
+
Connectives
+
+
Conjunction (logical and): \(X_1 \land X_2\)
+
Disjunction (logical or): \(X_1 \lor X_2\)
+
Negation: \(\neg X_2\) (or \(\overline{X_2}\) or \(\sim X_2\))
The last one can be seen as follows: it cannot happen that \(X1 \land \neg X2\), that is, it holds that \(\neg(X1 \land \neg X2)\). De Morgan gives \(\neg X1 \lor X2\).
+
+
+
Binary variables related to the Boolean ones
+
Ultimately we want to use numerical solvers to solve optimization problems. In order to do that, we need to associate with the Boolean variable \(X\) a binary (integer) variable \(\delta\in\{0,1\}\) such that \[
+\delta =
+\begin{cases}
+0 & \text{if} \; \neg X,\\
+1 & \text{if} \; X.
+\end{cases}
+\]
+
In other words, have the the following equivalence \[X \qquad \equiv \qquad [\delta = 1].\]
+
+
+
Integer (in)equalities related to the logical formulas
+
Having introduced binary variables, we now aim at reformulating logical formulas into equivalent integer equalities or inequalities.
+
+
And (conjunction)
+
We first consider the conjunction of two logical variables \(X_1\) and \(X_2\): \[X_1 \land X_2.\]
+
Invoking the binary variables, we can rewrite this as \[[\delta_1=1] \land [\delta_2=1],\] which ultimately translates to \[\delta_1=1, \; \delta_2=1.\]
+
Another possibility is \[\delta_1 \delta_2 = 1,\] but we discard it because it is nonlinear and we prefer the linear version for obvious reasons.
+
+
+
Or (disjunction)
+
\[X_1 \lor X_2\]
+
\[[\delta_1=1] \lor [\delta_2=1]\]
+
\[\delta_1 + \delta_2\geq 1\]
+
+
+
Negation
+
\[\neg X_1\]
+
\[\neg [\delta_1=1]\]
+
\[\delta_1 = 0\]
+
+
+
Xor
+
\[X_1 \oplus X_2\]
+
\[[\delta_1=1] \oplus [\delta_2=1]\]
+
\[\delta_1 + \delta_2 = 1\]
+
+
+
Implication
+
\[X_1 \implies X_2\]
+
\[[\delta_1=1] \implies [\delta_2=1]\]
+
Equivalently, using \(\neg X_1 \lor X_2\) we get \[\neg [\delta_1=1] \lor [\delta_2=1],\] which translates to \[(1-\delta_1) + \delta_2\geq 1,\] which simplifies to \[\delta_1 - \delta_2 \leq 0.\]
Expressing the equivalence using implications \[X_3 \implies X_1,\; X_3\implies X_2, \; (X_1 \land X_2) \implies X_3\]
+
The the last one is equivalent to \[\neg (X_1 \land X_2) \lor X_3,\] which can be simplified to \[\neg X_1 \lor \neg X_2 \lor X_3,\] which translates to \[\neg [\delta_1=1] \lor \neg [\delta_2 = 1] \lor [\delta_3 = 1],\] which finally leads to the inequality \[(1-\delta_1) + (1-\delta_2) + \delta_3 \geq 1\] and after simplification \[\delta_1 + \delta_2 - \delta_3 \leq 1.\]
+
And don’t forget to consider the first two inequalities too \[
+-\delta_1 + \delta_3 \leq 0, \quad -\delta_2 + \delta_3 \leq 0.
+\]
+
+
+
+
General transformation of Boolean expressions to integer inequalities
+
From Conjunctive Normal Form (CNF) \[
+\bigwedge_{j=1}^m \left[\left(\lor_{i\in \mathcal{P}_j} X_i\right) \lor \left(\lor_{i\in \mathcal{N}_j} \neg X_i\right)\right]
+\] to 0-1 integer inequalities defining a polyhedron \[
+\begin{aligned}
+\sum_{i\in \mathcal{P}_1} \delta_i + \sum_{i\in \mathcal{N}_1} (1-\delta_i) &\geq 1,\\
+&\vdots\\
+\sum_{i\in \mathcal{P}_m} \delta_i + \sum_{i\in \mathcal{N}_m} (1-\delta_i) &\geq 1.
+\end{aligned}
+\]
+
+
+
Finite state machine (FSM) using binary variables
+
Encode the discrete state variables in binary \[
+x_b \in \{0,1\}^{n_b}
+\]
+
Similarly the discrete inputs \[
+u_b \in \{0,1\}^{m_b}
+\]
+
The logical state equation then \[
+x_b(k+1) = f_b(x_b(k),u_b(k),\delta_e(k))
+\]
It can be expressed as a product\[
+z = \delta\,(a^\top x + b^\top u + f)
+\]
+
+
\[
+\begin{aligned}
+z &\leq M\delta,\\
+- z &\leq -m\delta,\\
+z &\leq a^\top x + b^\top u + f - m(1-\delta),\\
+-z &\leq -(a^\top x + b^\top u + f) + M(1-\delta).
+\end{aligned}
+\]
+
+
The reasoning is that if \(\delta=0\), then \(z\) is restricted, while \(a^\top x + b^\top u + f\) is not. And the other way around.
+
+
+
+
Another IF-THEN-ELSE rule
+
+
If \(X\)
+
+
then \(z = a_1^\top x + b_1^\top u + f_1\),
+
else \(z = a_2^\top x + b_2^\top u + f_2\).
+
+
It can be expressed as \[
+\begin{aligned}
+z &= \delta\,(a_1^\top x + b_1^\top u + f_1) \\
+&\quad + (1-\delta)(a_2^\top x + b_2^\top u + f_2)
+\end{aligned}
+\]
+
+
\[
+\begin{aligned}
+(m_2-M_1)\delta + z &\leq a_2^\top x + b_2^\top u + f_2,\\
+(m_1-M_2)\delta - z &\leq -a_2^\top x - b_2^\top u - f_2,\\
+(m_1-M_2)(1-\delta) + z &\leq a_1^\top x + b_1^\top u + f_1,\\
+(m_2-M_1)(1-\delta) - z &\leq -a_1^\top x - b_1^\top u - f_1.
+\end{aligned}
+\]
+
+
+
Generation of events by mixing logical and continuous variables in inequalities
We are going to introduce yet another framework for modeling hybrid systems – mixed logical dynamical (MLD) description. A question must inevitably pop up: “why yet another framework?”
+
The answer is, that we would like to have a model of a hybrid system that is suitable for model predictive control (MPC). Recall that the role of the model in MPC is that the model is used to define some constraints (equations and inequalities) in the numerical optimization problem. The frameworks that we considered so far did not offer it.
+
In particular, with the state variable and control input vectors composed of continuous and discrete variables
+\bm x = \begin{bmatrix}\bm x_c\\\bm x_d\end{bmatrix}, \quad \bm u = \begin{bmatrix}\bm u_c\\\bm u_d\end{bmatrix},
+ where \bm x_c\in\mathbb R^{n_c},\;\bm x_d\in\mathbb N^{n_d},\; \bm u_c\in\mathbb R^{m_c} and \bm u_d\in\mathbb N^{m_d}, we would like to formulate the model in the form of state equations, say
+\begin{aligned}
+\begin{bmatrix}\bm x_c(k+1) \\ \bm x_d(k+1)\end{bmatrix}
+&=
+\begin{bmatrix} \mathbf f_c(\bm x(k), \bm u(k)) \\ \mathbf f_d(\bm x(k), \bm u(k)) \end{bmatrix}
+\end{aligned}
+
+
Is it possible?
+
Unfortunately no. At least not exactly in this form. But something close to it is achievable instead.
+
But first we need to set the terminology and notation used to define a discrete(-time) hybrid automaton.
Model predictive control (MPC) is not computationally cheap (compared to, say, PID or LQG control) as it requires solving optimization problem – typically a quadratic program (QP) - online. The optimization solver needs to be a part of the controller.
+
There is an alternative, though, at least in same cases. It is called explicit MPC. The computationally heavy optimization is only perfomed only during the design process and the MPC controller is then implemented just as an affine state feedback
with the coefficients picked from some kind of a lookup table in real time Although retreiving the coefficients of the feedback controller is not computationally trivial, still it is cheaper than full optimization.
+
+
Multiparametric programming
+
The key technique for explicit MPC is multi-parametric programming. In order to explain it, consider the following problem
+
+J^\ast(x) = \inf_z J(z;x).
+
+
The z variable is an optimization variable, while x is a parameter. For a given parameter x, the cost function J is minimized. We study how the optimal cost J^\ast depends on the parameter, hence the name parametric programming. If x is a vector, the name of the problem changes to multiparametric programming.
+
+
Example: scalar variable, single parameter
+
Consider the following cost function J(z;x) in z parameterized by x. The optimization variable z is constrained and this constraint is also parameterized by x.
+\begin{aligned}
+J(z;x) &= \frac{1}{2} z^2 + 2zx + 2x^2 \\
+\text{subject to} &\quad z \leq 1 + x.
+\end{aligned}
+
+
In this simple case we can aim at analytical solution. We proceed in the standard way – we introduce a Lagrange multiplicator \lambda and form the augmented cost function
+L(z,\lambda; x) = \frac{1}{2} z^2 + 2zx + 2x^2 + \lambda (z-1-x).
+
+
The necessary conditions of optimality for the inequality-constrained problem come in the form of KKT conditions
+\begin{aligned}
+z + 2x + \lambda &= 0,\\
+z - 1 - x &\leq 0,\\
+\lambda & \geq 0,\\
+\lambda (z - 1 - x) &= 0.
+\end{aligned}
+
+
The last condition – the complementarity condition – gives rise to two scenarios: one corresponding to \lambda = 0, and the other corresponding to z - 1 - x = 0. We consider them separately below.
+
After substituting \lambda = 0 into the KKT conditions, we get
+\begin{aligned}
+z + 2x &= 0,\\
+z - 1 - x & \leq 0.
+\end{aligned}
+
+
From the first equation we get how z depends on x, and from the second we obtain a bound on x. Finally, we can also substitute the expression for z into the cost function J to get the optimal cost J^\ast as a function of x. All these are summarized here
+\begin{aligned}
+z &= -2x,\\
+x & \geq -\frac{1}{3},\\
+J^\ast(x) &= 0.
+\end{aligned}
+
+
Now, the other scenario. Upon substitutin z - 1 - x = 0 into the KKT conditions we get
From the second equation we get the expression for z in terms of x, substituting into the first equation and invoking the condition on nonnegativity of \lambda we get the bound on x (not suprisingly it complements the one obtained in the previous scenario). Finally, substituting for z in the cost function J we get a formula for the cost J^\ast as a function of x.
The two scenarios can now be combined into a single piecewise affine function z(x)
+z(x) = \begin{cases}
+1+x & \text{if } x \leq -\frac{1}{3},\\
+-2x & \text{if } x > -\frac{1}{3}.
+\end{cases}
+
+
x = range(-1, 1, length=100)
+z(x) = x <= -1/3 ? 1 + x : -2x
+Jstar(x) = x <= -1/3 ? 9/2*x^2 + 3x + 1/2 : 0
+
+using Plots
+plot(x, z.(x), label="z(x)")
+vline!([-1/3],line=:dash)
+xlabel!("x")
+ylabel!("z(x)")
+
and a piecewise quadratic cost function J^\ast(x)
+J^\ast(x) = \begin{cases}
+\frac{9}{2}x^2 + 3x + \frac{1}{2} & \text{if } x \leq -\frac{1}{3},\\
+0 & \text{if } x > -\frac{1}{3}.
+\end{cases}
+
The main resource for us is the Chapter 17 of the freely available book [1] that we already referred to in the previous chapter.
+
Those who have not been exposed to fundamentals of MPC can check the Chapter 12 in the same book. Alternatively, our own introduction to the topic in the third chapter/week of the Optimal and robust control course can be found useful.
Literature for Petri nets is vast, but a decent (and perfectly satisfactory) introduction can be found in Chapter 4 and 5.3 (for the timed PN) of the classical (and award-winning) reference [1]. Note that electronic version (in fact, PDF) is accessible through the NTK library (upon CTU login, for example to usermap first).
+
A nice introduction is also in Chapter 2 of the freely online available book [2].
+
The survey paper that is particularly focused on Petri nets from the control systems perspective is [3] and it gives a wealth of other references.
+
A few more monographs, mostly inclined towards control systems, are [4], [5], [6].
+
+
Petri nets and their derivatives such as Grafcet in international standards
+
We mention at the beginning of this chapter that Petri nets have made it to international standards. Here they are: [7], [8], and [9].
+
Based on Petri nets, another framework has been derived and standardized, namely GRAFCET, see [10] and [11], upon which, in turn, the popular Sequential Function Chart (SFC) language for PLC programming [12] is based.
A. Giua and M. Silva, “Petri nets and Automatic Control: A historical perspective,”Annual Reviews in Control, vol. 45, pp. 223–239, Jan. 2018, doi: 10.1016/j.arcontrol.2018.04.006.
+
+
+
[4]
J. O. Moody, Supervisory Control of Discrete Event Systems Using Petri Nets. in The International Series on Discrete Event Dynamic Systems. New York, NY: Springer, 31 {\v c}ervence 1998. Available: https://doi.org/10.1007/978-1-4615-5711-1
+
+
+
[5]
B. Hrúz and M. Zhou, Modeling and Control of Discrete-event Dynamic Systems: With Petri Nets and Other Tools. in Advanced Textbooks in Control and Signal Processing (C&SP). London: Springer, 2007. Available: https://doi.org/10.1007/978-1-84628-877-7
+
+
+
[6]
W. Reisig, Understanding Petri Nets: Modeling Techniques, Analysis Methods, Case Studies. Berlin; Heidelberg: Springer, 2013. Available: https://doi.org/10.1007/978-3-642-33278-4
+
+
+
[7]
“ISO/IEC 15909-1:2019 Systems and software engineering — High-level Petri nets — Part 1: Concepts, definitions and graphical notation.” ISO/IEC, Aug. 2019. Accessed: Sep. 27, 2023. [Online]. Available: https://www.iso.org/standard/67235.html
+
+
+
[8]
“ISO/IEC 15909-2:2011 Systems and software engineering — High-level Petri nets — Part 2: Transfer format.” ISO/IEC, Feb. 2011. Accessed: Sep. 27, 2023. [Online]. Available: https://www.iso.org/standard/43538.html
+
+
+
[9]
“ISO/IEC 15909-3:2021: Systems and software engineering — High-level Petri nets — Part 3: Extensions and structuring mechanisms.” ISO/IEC, 2021. Accessed: Sep. 29, 2023. [Online]. Available: https://www.iso.org/standard/81504.html
C. Johnsson and K.-E. Årzén, “Grafchart and its Relations to Grafcet and Petri Nets,”IFAC Proceedings Volumes, vol. 31, no. 15, pp. 95–100, Jun. 1998, doi: 10.1016/S1474-6670(17)40535-0.
Recall that when introducing enabled transitions, we emphasized that these can but do not have to fire immediately after having been enabled \boxed{\mathrm{ENABLING} \neq \text{FIRING}.}
+
+
Delays associated with transitions
+
Well then, the enabled transitions do not have to fire immediately, but they can possibly fire with some delay after enabling. This is for the first time that we are introducing the concept of time into the Petri nets, isn’t it?
+
For the jth transition, the delay of the kth firing is v_{j,k}, and we collect the sequence of delayes into v_j = \{v_{j,1}, v_{j,2}, \ldots \}.
+
But not all transitions have to be timed. Denote the timed transitions \mathcal{T}_\mathrm{D}\subseteq \mathcal{T}. We define the clock structure for a PN as \mathcal{V} = \{v_j\mid t_j\in\mathcal{T}_\mathrm{D}\}.
+
The definition of a timed Petri net (TPN) is then obtained by augmenting the definition of a Petri net with the clock structure
Example 1 (Timed Petri net) Model of processing multiple tasks: task 1 and task 2 are processed sequentially, and task 3 is processed in parallel with them; task 4 can only be processed after both tasks 2 and 3 have been finished. Finishing individual tasks corresponds to the individual transitions. The transition 4 is untimed, it only expresses the logical requirement.
+
+
+
+Example of a timed Petri net
+
+
+
+
+
+
+
+
+
+Rectangles instead of bars
+
+
+
+
Sometimes instead of a bar, the untimed transitions are modelled using similarly thin rectangles as the timed transitions, but filled.
+
+
+
+
Places can also be delayed
+
With delays associated with just one type of a node in a Petri net, the situation is rather nonsymmetric. In some literature, delays can also associated with places. And yet in some other literature delays are only associated with places. Such delays associated with places are called holding time for a place It is the minimum duration the token must rest in the place. But the token can stay longer if the output transition is waiting for other places.
+
+
+
+
+
+
+Delays associated with transitions and places
+
+
+
+
There is a major difference in delays associated with places compared to the delays associated with transitions. While the former tells the minimum duration the token has to dwell in the place, the latter tell the exact delay with which the transition fires after having been enabled.
+
+
+
+
+
+
Timed Petri net dynamics
+
With the introduction of time into the Petri nets, we can now study the dynamics of the system. For general Petri nets, alhough perfectly doable, it may quicly become too complex, and therefore here we only consider event graphs.
+
Some notation:
+
+
\{\tau_{j,1}, \tau_{j,2}, \ldots\} are the firing times of the jth transition,
+
\{\pi_{i,1},\pi_{i,2},\ldots\} are the times when the ith place receives a token,
+
x_i = x(p_i) is the number of tokens at the ith place,
+
x_{i,k} = \left.x(p_i)\right|_k, is the number of tokens at the ith place after the kth firing.
+
+
Now, assume first that x_{i,0} = 0. We can then relate the time of the arrival of the token to the place with the firing of the transition from which the token arrives \pi_{i,k} = \tau_{j,k},\quad p_i\in \mathcal{O}(t_j).
+
But generally x_{i,0} \neq 0 and the above relation needs to be modified to \pi_{i,k+x_{i,0}} = \tau_{j,k},\quad p_i\in \mathcal{O}(t_j), or, equivalently \boxed{\pi_{i,k} = \tau_{j,k-x_{i,0}},\quad p_i\in \mathcal{O}(t_j).} \tag{1}
+
This can be read in the following way. If there are initially, say, 3 tokens in the place, the time of the arrival of the 4th token is the time of the first firing of the transition from which the 4th token arrives.
+
Good. Keep this result in mind. Now we go for another.
+
For an untimed transition with a single input place, the firing time is the same as the time of the arrival of the token to the place
+\tau_{j,k} = \pi_{i,k}.
+
+
Modifying this result for a timed transition with a single input place we get
+\tau_{j,k} = \pi_{i,k} + v_{j,k}.
+
+
In words, the firing time is given by the time of the arrival of the token to the place, which enables the transition, and the delay associated with the transition.
+
Finally, we extend this result to the case of a timed transition with multiple input places \boxed{
+\tau_{j,k} = \max_{p_i\in\mathcal{I}(t_j)}\{\pi_{i,k}\} + v_{j,k}.}
+\tag{2}
+
This is the other promised important result. Keep both boxed formulas Equation 1 and Equation 2 handy, they will be needed in what is coming.
+
+
Example 2 (Timed Petri net dynamics) Consider the Petri net with three places and two transitions, one of which is timed, as in Fig 1.
+
+
+
+
+
+
+Figure 1: Example of a Petri net for which the dynamics is analyzed
+
+
+
+
We first use Equation 2 to write down the firing times of the two transitions
+\begin{aligned}
+\tau_{1,k} &= \max\{\pi_{1,k},\pi_{3,k}\}\\
+\tau_{2,k} &= \pi_{2,k}+v_{2,k}.
+\end{aligned}
+
+
Now we apply Equation 1 to write down the times of the arrival of the tokens to the places
+\begin{aligned}
+\pi_{1,k} &= \tau_{1,k-1}, \qquad k=2,\ldots, \qquad \pi_{1,0} = 0\\
+\pi_{2,k} &= \tau_{1,k-1}, \qquad k=2,\ldots, \qquad \pi_{2,0} = 0\\
+\pi_{3,k} &= \tau_{2,k}, \qquad k=1,\ldots
+\end{aligned}
+
+
Substituting from the latter into the former we get
+\begin{aligned}
+\tau_{1,k} &= \max\{\tau_{1,k-1},\tau_{1,k-1}+v_{2,k}\}\\
+&= \tau_{1,k-1}+v_{2,k}, \quad \tau_{1,k} = 0\\
+\tau_{2,k} &= \tau_{1,k-1}+v_{2,k}.
+\end{aligned}
+
+
This is the ultimate model for the dynamics of the Petri net. Should we need it, we can also get similar expressions for the times of the arrival of the tokens to the places.
+
+
+
+
+
+
+
+Update equations for times and not states
+
+
+
+
While with state equations we compute a sequences of values of the state vector (\bm x_0, \bm x_1, \bm x_2, \ldots), in other words, we compute the evolution of the state in time, here we compute the sequences of times when transitions fire (or token arrive to places). This update scheme for times resembles the state equations, but the interpretation is different.
+
+
+
+
+
Queueing system using TPN
+
We can also model a queueing system using a TPN. The Petri net is shown in Fig 2.
+
+
+
+
+
+
+Figure 2: Timed Petri net modelling a queueing system
+
+
+
+
Of the three transitions \mathcal{T} = \{a,s,c\}, which we have already identified previously, we assume that only are times, namely \mathcal{T}_\mathrm{D} = \{a,c\}. The associated firing delays are \bm v = \begin{bmatrix}v_a \\ v_c\end{bmatrix}.
+
For convenience we relabel the firing times of the transitions. Instead of t_{a,k} we will use a_k, and similarly s_k, and c_k. Application of Equation 2 and Equation 1 gives
+\begin{aligned}
+a_k &= a_{k-1} + v_{a,k},\quad k=1,2,\ldots,\quad a_0 = 0\\
+s_k &= \max\{\pi_{Q,k},\pi_{I,k}\}\\
+c_k &= \pi_{B,k} + v_{c,k}\\
+\pi_{Q,k} &= a_{k},\quad k=1,2,\ldots\\
+\pi_{I,k} &= c_{k-1},\quad k= 2, \ldots, \quad \pi_{I,0}=1\\
+\pi_{B,k} &= s_{k},\quad k=1,2,\ldots\\
+\end{aligned}
+
The time of completing the kth task is given by the time at which the previous task was completed and the time needed to complete the kth task itself, unless there is a gap in the queue after finishing the previous task, in which case the server must wait for the next task to arrive.
+
+
Example 3 (Timed Petri net for synchronization of train lines) We consider three closed rail tracks and two stations as in Fig 3.
+
+
+
+
+
+
+Figure 3: Example with three train lines
+
+
+
+
Departure of a train at a station must be synchronized with arrival of the other train so that passengers can change train. The timed Petri net for this system is shown in Fig 4.
+
+
+
+
+
+
+Figure 4: Timed Petri net for the example of synchronization of three train lines
+
+
+
+
If time is associated with the places, the Petri net simplifies significantly to Fig 5.
+
+
+
+
+
+
+Figure 5: Timed Petri net for the example of synchronization of three train lines
+
+
+
+
+
+
Example 4 (Manufacturing) tbd
+
+
+
+
Extensions
+
+
Stochastic Petri nets (SPN)
+
Numerous extensions are possible, some of which we have already mentioned when discussing untimed Petri nets. But upon introducing time, stochastic Petr nets can be concived, in which delays are random variables.
+
+
+
+
+
+
\ No newline at end of file
diff --git a/search.json b/search.json
index 20514af..2a7f060 100644
--- a/search.json
+++ b/search.json
@@ -587,7 +587,7 @@
"href": "des_automata.html#extensions",
"title": "State automata",
"section": "Extensions",
- "text": "Extensions\nThe concept of an automaton can be extended in several ways. In particular, the following two extensions introduce the concept of an output to an automaton.\n\nMoore machine\nOne extension of an automaton with outputs is Moore machine. The outputs assigned to the states by the output function y = g(x).\nThe output is produced (emitted) when the (new) state is entered.\nNote, in particular, that the output does not depend on the input. This has a major advantage when a feedback loop is closed around this system, since no algebraic loop is created.\nGraphically, we make a conventions that outputs are the labels of the states.\n\nExample 8 (Moore machine) The following automaton has just three states, but just two outputs (FLOW and NO FLOW).\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\nclosed\n\nNO FLOW\nValve\nclosed\n\n\n\ninit->closed\n\n\n\n\n\npartial\n\nFLOW\nValve\npartially\nopen\n\n\n\nclosed->partial\n\n\nopen valve one turn\n\n\n\npartial->closed\n\n\nclose valve one turn\n\n\n\nfull\n\nFLOW\nValve\nfully open\n\n\n\npartial->full\n\n\nopen valve one turn\n\n\n\nfull->closed\n\n\nemergency shut off\n\n\n\nfull->partial\n\n\nclose valve one turn\n\n\n\n\n\n\nFigure 9: Example of a digraph representation of the Moore machine for a valve control\n\n\n\n\n\n\n\n\nMealy machine\nMealy machine is another extension of an automaton. Here the outputs are associated with the transitions rather than the states.\nSince the events already associated with the states can be viewed as the inputs, we now have input/output transition labels. The transition label e_\\mathrm{i}/e_\\mathrm{o} on the transion from x_1 to x_2 reads as “the input event e_\\mathrm{i} at state x_1 activates the transition to x_2, which outputs the event e_\\mathrm{o}” and can be written as x_1\\xrightarrow{e_\\mathrm{i}/e_\\mathrm{o}} x_2.\nIt can be viewed as if the output function also considers the input and not only the state y = e_\\mathrm{o} = g(x,e_\\mathrm{i}).\nIn contrast with the Moore machine, here the output is produced (emitted) during the transition (before the new state is entered).\n\nExample 9 (Mealy machine) Coffee machine: coffee for 30 CZK, machine accepting 10 and 20 CZK coins, no change.\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\n0\n\nNo coin\n\n\n\ninit->0\n\n\n\n\n\n10\n\n10 CZK\n\n\n\n0->10\n\n\ninsert 10 CZK / no coffee\n\n\n\n20\n\n20 CZK\n\n\n\n0->20\n\n\ninsert 20 CZK / no coffee\n\n\n\n10->0\n\n\ninsert 20 CZK / coffee\n\n\n\n10->20\n\n\ninsert 10 CZK / no coffee\n\n\n\n20->0\n\n\ninsert 10 CZK / coffee\n\n\n\n20->10\n\n\ninsert 20 CZK / coffee\n\n\n\n\n\n\nFigure 10: Example of a digraph representation of the Mealy machine for a coffee machine\n\n\n\n\n\n\n\nExample 10 (Reformulate the previous example as a Moore machine) Two more states wrt Mealy\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\n0\n\nNO COFFEE\nNo\ncoin\n\n\n\ninit->0\n\n\n\n\n\n10\n\nNO COFFEE\n10\nCZK\n\n\n\n0->10\n\n\ninsert 10 CZK\n\n\n\n20\n\nNO COFFEE\n20\nCZK\n\n\n\n0->20\n\n\ninsert 20 CZK\n\n\n\n10->20\n\n\ninsert 10 CZK\n\n\n\n30\n\nCOFFEE\n10+20\nCZK\n\n\n\n10->30\n\n\ninsert 20 CZK\n\n\n\n20->30\n\n\ninsert 10 CZK\n\n\n\n40\n\nCOFFEE\n20+20\nCZK\n\n\n\n20->40\n\n\ninsert 20 CZK\n\n\n\n30->0\n\n\n\n\n\n30->10\n\n\ninsert 10 CZK\n\n\n\n30->20\n\n\ninsert 20 CZK\n\n\n\n40->10\n\n\n\n\n\n40->20\n\n\ninsert 10 CZK\n\n\n\n40->30\n\n\ninsert 20 CZK\n\n\n\n\n\n\nFigure 11: Example of a digraph representation of the Moore machine for a coffee machine\n\n\n\n\n\n\n\n\n\n\n\n\nNote\n\n\n\nThere are transitions from 30 and 40 back to 0 that are not labelled by any event. This does not seem to follow the general rule that transitions are always triggered by events. Not what? It can be resolved upon introducing time as the timeout transitions.\n\n\n\nExample 11 (Dijkstra’s token passing) The motivation for this example is to show that it is perhaps not always productive to insist on visual description of the automaton using a graph. The four components of our formal definition of an automaton are just enough, and they translate directly to a code.\nThe example comes from the field of distributed computing systems. It considers several computers that are connected in ring topology, and the communication is just one-directional as Fig 12 shows. The task is to use the communication to determine in – a distributed way – which of the computers carries a (single) token at a given time. And to realize passing of the token to a neighbour. We assume a synchronous case, in which all the computers are sending simultaneously, say, with some fixed sending period.\n\n\n\n\n\n\n\n\nG\n\n\n0\n\n0\n\n\n\n1\n\n1\n\n\n\n0->1\n\n\n\n\n\n2\n\n2\n\n\n\n1->2\n\n\n\n\n\n3\n\n3\n\n\n\n2->3\n\n\n\n\n\n3->0\n\n\n\n\n\n\n\n\nFigure 12: Example of a ring topology for Dijkstra’s token passing in a distributed system\n\n\n\n\n\nOne popular method for this is called Dijkstra’s token passing. Each computer keeps a single integer value as its state variable. And it forwards this integer value to the neighbour (in the clockwise direction in our setting). Upon receiving the value from the other neighbour (in the counter-clockwise direction), it updates its own value according to the rule displayed in the code below. At every clock tick, the state vector (composed of the individual state variables) is updated according to the function update!() in the code. Based on the value of the state vector, an output is computed, which decodes the informovation about the location of the token from the state vector. Again, the details are in the output() function.\n\n\nShow the code\nstruct DijkstraTokenRing\n number_of_nodes::Int64\n max_value_of_state_variable::Int64\n state_vector::Vector{Int64}\nend\n\nfunction update!(dtr::DijkstraTokenRing) \n n = dtr.number_of_nodes\n k = dtr.max_value_of_state_variable\n x = dtr.state_vector\n xnext = copy(x)\n for i in eachindex(x) # Mind the +1 shift. x[2] corresponds to x₁ in the literature.\n if i == 1 \n xnext[i] = (x[i] == x[n]) ? mod(x[i] + 1,k) : x[i] # Increment if the left neighbour is identical.\n else \n xnext[i] = (x[i] != x[i-1]) ? x[i-1] : x[i] # Update by the differing left neighbour.\n end\n end\n dtr.state_vector .= xnext \nend\n\nfunction output(dtr::DijkstraTokenRing) # Token = 1, no token = 0 at the given position. \n x = dtr.state_vector\n y = similar(x)\n y[1] = iszero(x[1]-x[end])\n y[2:end] .= .!iszero.(diff(x))\n return y\nend\n\n\noutput (generic function with 1 method)\n\n\nWe now rund the code for a given number of computers and some initial state vector that does not necessarily comply with the requirement that there is only one token in the ring.\n\n\nShow the code\nn = 4 # Concrete number of nodes.\nk = n # Concrete max value of a state variable (>= n).\n@show x_initial = rand(0:k,n) # Initial state vector, not necessarily acceptable (>1 token in the ring).\ndtr = DijkstraTokenRing(n,k,x_initial)\n@show output(dtr) # Show where the token is (are).\n\n@show update!(dtr), output(dtr) # Perform the update, show the state vector and show where the token is.\n@show update!(dtr), output(dtr) # Repeat a few times to see the stabilization. \n@show update!(dtr), output(dtr)\n@show update!(dtr), output(dtr)\n@show update!(dtr), output(dtr)\n\n\nx_initial = rand(0:k, n) = [3, 1, 0, 3]\noutput(dtr) = [1, 1, 1, 1]\n(update!(dtr), output(dtr)) = ([0, 3, 1, 0], [1, 1, 1, 1])\n(update!(dtr), output(dtr)) = ([1, 0, 3, 1], [1, 1, 1, 1])\n(update!(dtr), output(dtr)) = ([2, 1, 0, 3], [0, 1, 1, 1])\n(update!(dtr), output(dtr)) = ([2, 2, 1, 0], [0, 0, 1, 1])\n(update!(dtr), output(dtr)) = ([2, 2, 2, 1], [0, 0, 0, 1])\n\n\n([2, 2, 2, 1], [0, 0, 0, 1])\n\n\nWe can see that although initially the there can be more tokens, after a few iterations the algorithm achieves the goal of having just one token in the ring.\n\n\n\nExtended-state automaton\nYet another extension of an automaton is the extended-state automaton. And indeed, the hyphen is there on purpose as we extend the state space.\nIn particular, we augment the state variable(s) that define the states/modes/locations (the nodes in the graph) by additional (typed) state variables: Int, Enum, Bool, …\nTransitions from one mode to another are then guarded by conditions on theses new extra state variables.\nBesides being guarded by a guard condition, a given transition can also be labelled by a reset function that resets the extended-state variables.\n\nExample 12 (Counting up to 10) In this example, there are two modes (on and off), which can be captured by a single binary state variable, say x. But then there is an additional integer variable k, and the two variables together characterize the extended state.\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\nOFF\n\nOFF\n\n\n\ninit->OFF\n\n\nint k=0\n\n\n\nON\n\nON\n\n\n\nOFF->ON\n\n\npress\n\n\n\nON->OFF\n\n\n(press ⋁ k ≥ 10); k=0\n\n\n\nON->ON\n\n\n(press ∧ k < 10); k=k+1\n\n\n\n\n\n\nFigure 13: Example of a digraph representation of the extended-state automaton for counting up to ten",
+ "text": "Extensions\nThe concept of an automaton can be extended in several ways. In particular, the following two extensions introduce the concept of an output to an automaton.\n\nMoore machine\nOne extension of an automaton with outputs is Moore machine. The outputs assigned to the states by the output function y = g(x).\nThe output is produced (emitted) when the (new) state is entered.\nNote, in particular, that the output does not depend on the input. This has a major advantage when a feedback loop is closed around this system, since no algebraic loop is created.\nGraphically, we make a conventions that outputs are the labels of the states.\n\nExample 8 (Moore machine) The following automaton has just three states, but just two outputs (FLOW and NO FLOW).\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\nclosed\n\nNO FLOW\nValve\nclosed\n\n\n\ninit->closed\n\n\n\n\n\npartial\n\nFLOW\nValve\npartially\nopen\n\n\n\nclosed->partial\n\n\nopen valve one turn\n\n\n\npartial->closed\n\n\nclose valve one turn\n\n\n\nfull\n\nFLOW\nValve\nfully open\n\n\n\npartial->full\n\n\nopen valve one turn\n\n\n\nfull->closed\n\n\nemergency shut off\n\n\n\nfull->partial\n\n\nclose valve one turn\n\n\n\n\n\n\nFigure 9: Example of a digraph representation of the Moore machine for a valve control\n\n\n\n\n\n\n\n\nMealy machine\nMealy machine is another extension of an automaton. Here the outputs are associated with the transitions rather than the states.\nSince the events already associated with the states can be viewed as the inputs, we now have input/output transition labels. The transition label e_\\mathrm{i}/e_\\mathrm{o} on the transion from x_1 to x_2 reads as “the input event e_\\mathrm{i} at state x_1 activates the transition to x_2, which outputs the event e_\\mathrm{o}” and can be written as x_1\\xrightarrow{e_\\mathrm{i}/e_\\mathrm{o}} x_2.\nIt can be viewed as if the output function also considers the input and not only the state y = e_\\mathrm{o} = g(x,e_\\mathrm{i}).\nIn contrast with the Moore machine, here the output is produced (emitted) during the transition (before the new state is entered).\n\nExample 9 (Mealy machine) Coffee machine: coffee for 30 CZK, machine accepting 10 and 20 CZK coins, no change.\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\n0\n\nNo coin\n\n\n\ninit->0\n\n\n\n\n\n10\n\n10 CZK\n\n\n\n0->10\n\n\ninsert 10 CZK / no coffee\n\n\n\n20\n\n20 CZK\n\n\n\n0->20\n\n\ninsert 20 CZK / no coffee\n\n\n\n10->0\n\n\ninsert 20 CZK / coffee\n\n\n\n10->20\n\n\ninsert 10 CZK / no coffee\n\n\n\n20->0\n\n\ninsert 10 CZK / coffee\n\n\n\n20->10\n\n\ninsert 20 CZK / coffee\n\n\n\n\n\n\nFigure 10: Example of a digraph representation of the Mealy machine for a coffee machine\n\n\n\n\n\n\n\nExample 10 (Reformulate the previous example as a Moore machine) Two more states wrt Mealy\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\n0\n\nNO COFFEE\nNo\ncoin\n\n\n\ninit->0\n\n\n\n\n\n10\n\nNO COFFEE\n10\nCZK\n\n\n\n0->10\n\n\ninsert 10 CZK\n\n\n\n20\n\nNO COFFEE\n20\nCZK\n\n\n\n0->20\n\n\ninsert 20 CZK\n\n\n\n10->20\n\n\ninsert 10 CZK\n\n\n\n30\n\nCOFFEE\n10+20\nCZK\n\n\n\n10->30\n\n\ninsert 20 CZK\n\n\n\n20->30\n\n\ninsert 10 CZK\n\n\n\n40\n\nCOFFEE\n20+20\nCZK\n\n\n\n20->40\n\n\ninsert 20 CZK\n\n\n\n30->0\n\n\n\n\n\n30->10\n\n\ninsert 10 CZK\n\n\n\n30->20\n\n\ninsert 20 CZK\n\n\n\n40->10\n\n\n\n\n\n40->20\n\n\ninsert 10 CZK\n\n\n\n40->30\n\n\ninsert 20 CZK\n\n\n\n\n\n\nFigure 11: Example of a digraph representation of the Moore machine for a coffee machine\n\n\n\n\n\n\n\n\n\n\n\n\nNote\n\n\n\nThere are transitions from 30 and 40 back to 0 that are not labelled by any event. This does not seem to follow the general rule that transitions are always triggered by events. Not what? It can be resolved upon introducing time as the timeout transitions.\n\n\n\nExample 11 (Dijkstra’s token passing) The motivation for this example is to show that it is perhaps not always productive to insist on visual description of the automaton using a graph. The four components of our formal definition of an automaton are just enough, and they translate directly to a code.\nThe example comes from the field of distributed computing systems. It considers several computers that are connected in ring topology, and the communication is just one-directional as Fig 12 shows. The task is to use the communication to determine in – a distributed way – which of the computers carries a (single) token at a given time. And to realize passing of the token to a neighbour. We assume a synchronous case, in which all the computers are sending simultaneously, say, with some fixed sending period.\n\n\n\n\n\n\n\n\nG\n\n\n0\n\n0\n\n\n\n1\n\n1\n\n\n\n0->1\n\n\n\n\n\n2\n\n2\n\n\n\n1->2\n\n\n\n\n\n3\n\n3\n\n\n\n2->3\n\n\n\n\n\n3->0\n\n\n\n\n\n\n\n\nFigure 12: Example of a ring topology for Dijkstra’s token passing in a distributed system\n\n\n\n\n\nOne popular method for this is called Dijkstra’s token passing. Each computer keeps a single integer value as its state variable. And it forwards this integer value to the neighbour (in the clockwise direction in our setting). Upon receiving the value from the other neighbour (in the counter-clockwise direction), it updates its own value according to the rule displayed in the code below. At every clock tick, the state vector (composed of the individual state variables) is updated according to the function update!() in the code. Based on the value of the state vector, an output is computed, which decodes the informovation about the location of the token from the state vector. Again, the details are in the output() function.\n\n\nShow the code\nstruct DijkstraTokenRing\n number_of_nodes::Int64\n max_value_of_state_variable::Int64\n state_vector::Vector{Int64}\nend\n\nfunction update!(dtr::DijkstraTokenRing) \n n = dtr.number_of_nodes\n k = dtr.max_value_of_state_variable\n x = dtr.state_vector\n xnext = copy(x)\n for i in eachindex(x) # Mind the +1 shift. x[2] corresponds to x₁ in the literature.\n if i == 1 \n xnext[i] = (x[i] == x[n]) ? mod(x[i] + 1,k) : x[i] # Increment if the left neighbour is identical.\n else \n xnext[i] = (x[i] != x[i-1]) ? x[i-1] : x[i] # Update by the differing left neighbour.\n end\n end\n dtr.state_vector .= xnext \nend\n\nfunction output(dtr::DijkstraTokenRing) # Token = 1, no token = 0 at the given position. \n x = dtr.state_vector\n y = similar(x)\n y[1] = iszero(x[1]-x[end])\n y[2:end] .= .!iszero.(diff(x))\n return y\nend\n\n\noutput (generic function with 1 method)\n\n\nWe now rund the code for a given number of computers and some initial state vector that does not necessarily comply with the requirement that there is only one token in the ring.\n\n\nShow the code\nn = 4 # Concrete number of nodes.\nk = n # Concrete max value of a state variable (>= n).\n@show x_initial = rand(0:k,n) # Initial state vector, not necessarily acceptable (>1 token in the ring).\ndtr = DijkstraTokenRing(n,k,x_initial)\n@show output(dtr) # Show where the token is (are).\n\n@show update!(dtr), output(dtr) # Perform the update, show the state vector and show where the token is.\n@show update!(dtr), output(dtr) # Repeat a few times to see the stabilization. \n@show update!(dtr), output(dtr)\n@show update!(dtr), output(dtr)\n@show update!(dtr), output(dtr)\n\n\nx_initial = rand(0:k, n) = [1, 1, 4, 1]\noutput(dtr) = [1, 0, 1, 1]\n(update!(dtr), output(dtr)) = ([2, 1, 1, 4], [0, 1, 0, 1])\n(update!(dtr), output(dtr)) = ([2, 2, 1, 1], [0, 0, 1, 0])\n(update!(dtr), output(dtr)) = ([2, 2, 2, 1], [0, 0, 0, 1])\n(update!(dtr), output(dtr)) = ([2, 2, 2, 2], [1, 0, 0, 0])\n(update!(dtr), output(dtr)) = ([3, 2, 2, 2], [0, 1, 0, 0])\n\n\n([3, 2, 2, 2], [0, 1, 0, 0])\n\n\nWe can see that although initially the there can be more tokens, after a few iterations the algorithm achieves the goal of having just one token in the ring.\n\n\n\nExtended-state automaton\nYet another extension of an automaton is the extended-state automaton. And indeed, the hyphen is there on purpose as we extend the state space.\nIn particular, we augment the state variable(s) that define the states/modes/locations (the nodes in the graph) by additional (typed) state variables: Int, Enum, Bool, …\nTransitions from one mode to another are then guarded by conditions on theses new extra state variables.\nBesides being guarded by a guard condition, a given transition can also be labelled by a reset function that resets the extended-state variables.\n\nExample 12 (Counting up to 10) In this example, there are two modes (on and off), which can be captured by a single binary state variable, say x. But then there is an additional integer variable k, and the two variables together characterize the extended state.\n\n\n\n\n\n\n\n\nG\n\n\ninit\ninit\n\n\n\nOFF\n\nOFF\n\n\n\ninit->OFF\n\n\nint k=0\n\n\n\nON\n\nON\n\n\n\nOFF->ON\n\n\npress\n\n\n\nON->OFF\n\n\n(press ⋁ k ≥ 10); k=0\n\n\n\nON->ON\n\n\n(press ∧ k < 10); k=k+1\n\n\n\n\n\n\nFigure 13: Example of a digraph representation of the extended-state automaton for counting up to ten",
"crumbs": [
"1. Discrete-event systems: Automata",
"State automata"
@@ -939,7 +939,7 @@
"href": "verification_barrier.html#convex-relaxation-of-the-barrier-certificate-problem",
"title": "Barrier certificates",
"section": "Convex relaxation of the barrier certificate problem",
- "text": "Convex relaxation of the barrier certificate problem\nWe relax the third condition so that it holds not only at B(\\bm x) = 0 but everywhere. The three conditions are then B(\\bm x) > 0,\\quad \\forall \\bm x \\in \\mathcal X_\\mathrm{u},\nB(\\bm x) \\leq 0,\\quad \\forall \\bm x \\in \\mathcal X_0,\n\\nabla B(\\bm x)^\\top \\mathbf f(\\bm x, \\bm d) \\leq 0,\\quad \\forall \\bm x\\in \\mathcal X, \\bm d \\in \\mathcal D.\nLet’s now demonstrate this by means of an example.\n\nExample 1 Consider the system modelled by \n\\begin{aligned}\n\\dot x_1 &= x_2\\\\\n\\dot x_2 &= -x_1 + \\frac{p}{3}x_1^3 - x_2,\n\\end{aligned}\n where the parameter p\\in [0.9,1.1].\nThe initial set is given by \n\\mathcal X_0 = \\{ \\bm x \\in \\mathbb R^2 \\mid (x_1-1.5)^2 + x_2^2 \\leq 0.25 \\}\n and the unsafe set is given by \n\\mathcal X_\\mathrm{u} = \\{ \\bm x \\in \\mathbb R^2 \\mid (x_1+1)^2 + (x_2+1)^2 \\leq 0.16 \\}.\n\nThe vector field \\mathbf f and the initial and unsafe sets are shown in the figure below.\n\n\nShow the code\nusing SumOfSquares\nusing DynamicPolynomials\nusing MosekTools \n\noptimizer = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)\nmodel = SOSModel(optimizer)\n@polyvar x[1:2] \n\np = 1;\n\nf = [ x[2],\n -x[1] + (p/3)*x[1]^3 - x[2]]\n\ng₁ = -(x[1]+1)^2 - (x[2]+1)^2 + 0.16 # 𝒳ᵤ = {x ∈ R²: g₁(x) ≥ 0}\nh₁ = -(x[1]-1.5)^2 - x[2]^2 + 0.25 # 𝒳₀ = {x ∈ R²: h₁(x) ≥ 0}\n\nX = monomials(x, 0:4)\n@variable(model, B, Poly(X))\n\nε = 0.001\n@constraint(model, B >= ε, domain = @set(g₁ >= 0))\n\n@constraint(model, B <= 0, domain = @set(h₁ >= 0))\n\nusing LinearAlgebra # Needed for `dot`\ndBdt = dot(differentiate(B, x), f)\n@constraint(model, -dBdt >= 0)\n\nJuMP.optimize!(model)\n\nJuMP.primal_status(model)\n\nimport DifferentialEquations, Plots, ImplicitPlots\nfunction phase_plot(f, B, g₁, h₁, quiver_scaling, Δt, X0, solver = DifferentialEquations.Tsit5())\n X₀plot = ImplicitPlots.implicit_plot(h₁; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label=\"X₀\", linecolor=:blue)\n Xᵤplot = ImplicitPlots.implicit_plot!(g₁; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label=\"Xᵤ\", linecolor=:teal)\n Bplot = ImplicitPlots.implicit_plot!(B; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label=\"B = 0\", linecolor=:red)\n Plots.plot(X₀plot)\n Plots.plot!(Xᵤplot)\n Plots.plot!(Bplot)\n ∇(vx, vy) = [fi(x[1] => vx, x[2] => vy) for fi in f]\n ∇pt(v, p, t) = ∇(v[1], v[2])\n function traj(v0)\n tspan = (0.0, Δt)\n prob = DifferentialEquations.ODEProblem(∇pt, v0, tspan)\n return DifferentialEquations.solve(prob, solver, reltol=1e-8, abstol=1e-8)\n end\n ticks = -5:0.5:5\n X = repeat(ticks, 1, length(ticks))\n Y = X'\n Plots.quiver!(X, Y, quiver = (x, y) -> ∇(x, y) / quiver_scaling, linewidth=0.5)\n for x0 in X0\n Plots.plot!(traj(x0), vars=(1, 2), label = nothing)\n end\n Plots.plot!(xlims = (-2, 3), ylims = (-2.5, 2.5))\nend\n\nphase_plot(f, value(B), g₁, h₁, 10, 30.0, [[x1, x2] for x1 in 1.2:0.2:1.7, x2 in -0.35:0.1:0.35])\n\n\n┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.\n│ caller = ip:0x0\n└ @ Core :-1\n┌ Warning: At t=4.423118290940112, dt was forced below floating point epsilon 8.881784197001252e-16, and step error estimate = 1.1390324186747878. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).\n└ @ SciMLBase ~/.julia/packages/SciMLBase/d2jpL/src/integrator_interface.jl:620",
+ "text": "Convex relaxation of the barrier certificate problem\nWe relax the third condition so that it holds not only at B(\\bm x) = 0 but everywhere. The three conditions are then B(\\bm x) > 0,\\quad \\forall \\bm x \\in \\mathcal X_\\mathrm{u},\nB(\\bm x) \\leq 0,\\quad \\forall \\bm x \\in \\mathcal X_0,\n\\nabla B(\\bm x)^\\top \\mathbf f(\\bm x, \\bm d) \\leq 0,\\quad \\forall \\bm x\\in \\mathcal X, \\bm d \\in \\mathcal D.\nLet’s now demonstrate this by means of an example.\n\nExample 1 Consider the system modelled by \n\\begin{aligned}\n\\dot x_1 &= x_2\\\\\n\\dot x_2 &= -x_1 + \\frac{p}{3}x_1^3 - x_2,\n\\end{aligned}\n where the parameter p\\in [0.9,1.1].\nThe initial set is given by \n\\mathcal X_0 = \\{ \\bm x \\in \\mathbb R^2 \\mid (x_1-1.5)^2 + x_2^2 \\leq 0.25 \\}\n and the unsafe set is given by \n\\mathcal X_\\mathrm{u} = \\{ \\bm x \\in \\mathbb R^2 \\mid (x_1+1)^2 + (x_2+1)^2 \\leq 0.16 \\}.\n\nThe vector field \\mathbf f and the initial and unsafe sets are shown in the figure below.\n\n\nShow the code\nusing SumOfSquares\nusing DynamicPolynomials\nusing MosekTools \n\noptimizer = optimizer_with_attributes(Mosek.Optimizer, MOI.Silent() => true)\nmodel = SOSModel(optimizer)\n@polyvar x[1:2] \n\np = 1;\n\nf = [ x[2],\n -x[1] + (p/3)*x[1]^3 - x[2]]\n\ng₁ = -(x[1]+1)^2 - (x[2]+1)^2 + 0.16 # 𝒳ᵤ = {x ∈ R²: g₁(x) ≥ 0}\nh₁ = -(x[1]-1.5)^2 - x[2]^2 + 0.25 # 𝒳₀ = {x ∈ R²: h₁(x) ≥ 0}\n\nX = monomials(x, 0:4)\n@variable(model, B, Poly(X))\n\nε = 0.001\n@constraint(model, B >= ε, domain = @set(g₁ >= 0))\n\n@constraint(model, B <= 0, domain = @set(h₁ >= 0))\n\nusing LinearAlgebra # Needed for `dot`\ndBdt = dot(differentiate(B, x), f)\n@constraint(model, -dBdt >= 0)\n\nJuMP.optimize!(model)\n\nJuMP.primal_status(model)\n\nimport DifferentialEquations, Plots, ImplicitPlots\nfunction phase_plot(f, B, g₁, h₁, quiver_scaling, Δt, X0, solver = DifferentialEquations.Tsit5())\n X₀plot = ImplicitPlots.implicit_plot(h₁; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label=\"X₀\", linecolor=:blue)\n Xᵤplot = ImplicitPlots.implicit_plot!(g₁; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label=\"Xᵤ\", linecolor=:teal)\n Bplot = ImplicitPlots.implicit_plot!(B; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label=\"B = 0\", linecolor=:red)\n Plots.plot(X₀plot)\n Plots.plot!(Xᵤplot)\n Plots.plot!(Bplot)\n ∇(vx, vy) = [fi(x[1] => vx, x[2] => vy) for fi in f]\n ∇pt(v, p, t) = ∇(v[1], v[2])\n function traj(v0)\n tspan = (0.0, Δt)\n prob = DifferentialEquations.ODEProblem(∇pt, v0, tspan)\n return DifferentialEquations.solve(prob, solver, reltol=1e-8, abstol=1e-8)\n end\n ticks = -5:0.5:5\n X = repeat(ticks, 1, length(ticks))\n Y = X'\n Plots.quiver!(X, Y, quiver = (x, y) -> ∇(x, y) / quiver_scaling, linewidth=0.5)\n for x0 in X0\n Plots.plot!(traj(x0), vars=(1, 2), label = nothing)\n end\n Plots.plot!(xlims = (-2, 3), ylims = (-2.5, 2.5))\nend\n\nphase_plot(f, value(B), g₁, h₁, 10, 30.0, [[x1, x2] for x1 in 1.2:0.2:1.7, x2 in -0.35:0.1:0.35])\n\n\n┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.\n│ caller = ip:0x0\n└ @ Core :-1\n┌ Warning: At t=4.423118290940107, dt was forced below floating point epsilon 8.881784197001252e-16, and step error estimate = 1.139033855908175. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).\n└ @ SciMLBase ~/.julia/packages/SciMLBase/tey0W/src/integrator_interface.jl:620",
"crumbs": [
"12. Formal verification",
"Barrier certificates"
diff --git a/sitemap.xml b/sitemap.xml
index 924b8a7..02c041c 100644
--- a/sitemap.xml
+++ b/sitemap.xml
@@ -134,7 +134,7 @@
https://hurak.github.io/hys/hybrid_equations.html
- 2024-10-20T08:00:48.774Z
+ 2024-10-20T08:12:24.581Zhttps://hurak.github.io/hys/des_references.html
diff --git a/stability_recap 5.html b/stability_recap 5.html
new file mode 100644
index 0000000..6e1790c
--- /dev/null
+++ b/stability_recap 5.html
@@ -0,0 +1,1251 @@
+
+
+
+
+
+
+
+
+
+Recap of stability analysis for continuous dynamical systems – B(E)3M35HYS – Hybrid systems
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Recap of stability analysis for continuous dynamical systems
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Before we start discussing stability of hybrid dynamical systems, it will not hurt to recapitulate the stability analysis for continuous (both in value and in time) dynamical systems modelled by the standard state equation
+
\dot{\bm{x}} = \mathbf f(\bm x).
+
+
Equilibrium
+
Although every now and then we may hear the term stability attributed to a system, strictly speaking it is an equilibrium that is stable or unstable.
+
Equilibrium is a point in the state space, that is, a vector x_\mathrm{eq}\in \mathbb R^n, such that
Without loss of generality \bm x_\mathrm{eq} = \mathbf 0. If it is nonzero, we can introduce a new shifted state vector.
+
+
+
Lyapunov stability
+
One of the most common types of stability is Lyapunov stability. Loosely speaking, it means that if the system starts close to the equilibrium, it stays close to it. More formally, for a given \varepsilon>0, there is a \delta>0 such that …
+
+
+
Attractivity
+
This is another property of an equilibrium. If it is (locally) attractive, it means that if the systems starts close to the equilibrium, it will converge to it. The global version of attractivity means that the system asymptotically converges to the equilibrium from anywhere.
+
Perhaps it is not immediately clear that attractivity is distinct from (Lyapunov) stability. The following example shows an attractive but Lyapunov unstable equilibrium.
Combination of Lyapunov stability and attractivity is called assymptotic stability.
+
If the attractivity is global, the assymptotic stability is called global too.
+
+
+
Exponential stability
+
Exponential convergence.
+
+
+
Stability of time-varying systems
+
Stability (Lyapunov, asymptotic, …) is called uniform, if it is independent of the inititial time.
+
+
+
Stability analysis via Lyapunov function
+
Now that we recapitulated the key stability concepts, it is time to recapitulate the methods of checking if this or that type of stability is achieved. The classical method is based on the searching for a Lyapunov function.
+
Lyapunov function is a scalar function V(\cdot)\in\mathcal{C}_1 defined on open \mathcal{D}\subset \mathbb{R}^n containing the origin (the equilibrium) that satisfies the following conditions V(0) = 0, \; V(x) > 0\, \text{for all}\, x\in \mathcal{D}\setminus \{0\},
+
\underbrace{\left(\nabla V(x)\right)^\top f(x)}_{\frac{\mathrm d}{\mathrm d t}V(x)} \leq 0.
+
In words, Lyapunov function for a given system and the equilibrium is a function that is positive everywhere except at the origin, where it is zero, we call such function positive definite, and its derivative along the trajectories of the system is nonpositive (aka positive semidefinite), which is a way to guarantee that the function does not increase along the trajectories. If such function exists, the equilibrium is Lyapunov stable.
+
If the latter condition is made strict, that is, if \left(\nabla V(x)\right)^\top f(x) < 0, which is a way to guarantee that the function decreases along the trajectories, the equilibrium is asymptotically stable.
+
The interpretation is quite intuitive:
+
+
+
LaSalle’s invariance principle
+
A delicate question is if the derivative of the Lyapunov function ocassionally vanishes, it it automatically means that the equilibrium is not assymptotically stable. The aswer is: not necessarily. LaSalle’s invariance principle states that even if the derivative of the Lyapunov function occasionally vanishes, the equilibrium can still be asymptotically stable, provided some condition is satisfied. We will not elaborate on it here. Look it up in your favourite nonlinear (control) system textbook.
+
+
+
Formulated using comparison functions
+
The above properties of the Lyapunov function be also be formulated using comparison functions. For Lyapunov stability, the following holds \kappa_1(\|x\|) \leq V(x) {\color{gray}\leq \kappa_2(\|x\|)}, where
+
+
\kappa_1(\cdot), \kappa_2(\cdot) are class \mathcal{K} comparison functions, that is, they are continuous, zero at zero and (strictly) increasing.
+
If \kappa_1 increases to infinity (\kappa_1(\cdot)\in\mathcal{K}_\infty), the stability is global.
+
+
For asymptotic stability
+
\left(\nabla V(x)\right)^\top f(x) \leq -\rho(\|x\|), where \rho(\cdot) is a positive definite continuous function, zero at the origin.
+
+
The upper bound \kappa_2(\cdot) does not have to be there, it is automatically satisfied for time-invariant systems. It does have to be imposed for time-varying systems though.
Multiple Lyapunov Function (MLF) approach to analysis of stability
+
+
Not just a single common Lyapunov function V(\bm x).
+
Instead, several Lyapunov-like functions V_i(\bm x),\; i=1,\ldots,r, that are Lyapunov on some subsets of the state space \Omega_i.
+
Stitch them together to form
+V(x) =
+\begin{cases}
+V_1(x) & \text{if } x\in \Omega_1, \\
+\vdots\\
+V_r(x) & \text{if } x\in \Omega_r, \\
+\end{cases}
+
+
+
+
+
S-procedure
+
+
A result about solvability of two or more quadratic inequalities, not necessarily convex ones.
+
+
For convex problems we have nonlinear Farkas’ lemma.
+
+
Origin in the control theory (analysis of stability of nonlinear systems, hence the S letter) with the first rigorous result provided by Yakubovich in 1970s.
+
It gives conditions under which (satisfaction of) one quadratic inequality follows from (satisfaction of) another one (or more), that is, it gives a condition under which the following implication holds:
+
Quadratic inequality #1 satisfied by some x\Rightarrow Quadratic inequality #0 satisfied by the same x.
+
In other words, it gives a condition under which the solution set of the inequality #1 denoted as \mathcal X_1 is included in the solution set \mathcal X_0 of the inequality #0.
+
+
+
+
S-procedure with nonstrict inequalities
+
+
Consider the general quadratic functions F_i(x) = x^\top A_i x + 2b_i^\top x + c_i, \; i=0,\ldots, p.
+
Question: under which conditions F_0(x) \geq 0 for all x satisfying F_i(x)\geq 0,\; i=1,\ldots,p ?
+
In other words, we are looking for conditions under which the implication
+F_i(x) \geq 0,\; i=1,\ldots,p \quad \Rightarrow \quad F_0(x) \geq 0
+
+holds.
+
In the simplest (yet relevant) case n=1 we search for conditions under which F_0(x) \geq 0 for all x satisfying F_1(x)\geq 0,
+
that is, conditions under which the implication
+F_1(x) \geq 0 \quad \Rightarrow \quad F_0(x) \geq 0
+
+holds.
+
+
+
+
Sufficient conditions
+
The existence of \alpha_i\geq 0,\; i=1,\ldots,p such that
+F_0(x)-\sum_{i=1}^p \alpha_i F_i(x) \geq 0.
+
+
+
Generally not necessary, the condition is conservative.
+
+
It can be formulated as an LMI
+\begin{bmatrix}
+A_0 & b_0 \\
+b_0^\top & c_0
+\end{bmatrix} -
+\sum_{i=1}^p
+\alpha_i
+\begin{bmatrix}
+A_i & b_i \\
+b_i^\top & c_i
+\end{bmatrix}
+\succeq 0
+
+where \alpha_i \geq 0,\; i=1,\ldots,p.
+
+
+
+
Sufficient and also necessary
+
+
It is nontrivial that for p=1 it is also necessary, provided that there is some x_0 such that F_1(x_0)>0 .
+
Then we have the following equivalence between the two constraints:
+\begin{aligned}
+F_0(x) &\geq 0 \; \forall x \;\mathrm{satisfying}\; F_1(x)\geq 0 \\
+&\Longleftrightarrow \\
+F_0(x)-\alpha F_1(x) &\geq 0,\;\text{for some}\; \alpha\in\mathbb{R}, \; \alpha\geq 0,
+\end{aligned}
+
+
which again can be formulated as an LMI, namely
+
+\begin{bmatrix}
+A_0 & b_0 \\
+b_0^\top & c_0
+\end{bmatrix} - \alpha
+\begin{bmatrix}
+A_1 & b_1 \\
+b_1^\top & c_1
+\end{bmatrix}
+\succeq 0,\quad \alpha\geq 0.
+
+
+
+
+
More on S-procedure
+
+
There are several variants
+
+
strict, nonstrict or mixed inequalities,
+
just two or more,
+
some of the constraints can be equations.
+
+
Literature (freely available online)
+
+
Boyd, S., El Ghaoui, L., Feron, E., Balakrishnan, V., 1994. Linear Matrix Inequalities in System and Control Theory. SIAM. Pages 23 and 24.
Quadratic Lyapunov-like functions, that is, quadratic functions V_i(\bm x) = \bm x^\top \mathbf P_i \bm x that qualify as Lyapunov functions only on respective subsets \Omega_i\sub \mathbb R^n:
Characterization of subsets of state space using LMI
+
+
Some subsets \Omega_i\sub \mathbb R^n characterized using linear and quadratic inequalities can be formulated as within the LMI framework as
+ \bm x^\top \mathbf Q_i \bm x \geq 0.
+
+
In particular, centered ellipsoids and cones.
+
For example,
+
+\begin{aligned}
+\Omega_i &= \{\bm x \in \mathbb R^n \mid (\mathbf c^\top \bm x \geq 0 \land \mathbf d^\top \bm x \geq 0) \\
+& \qquad \qquad \qquad \lor (\mathbf c^\top \bm x \leq 0 \land \mathbf d^\top \bm x \leq 0)\}.
+\end{aligned}
+
+
This constraint can be reformulated as
+(\mathbf c^\top \bm x) (\mathbf d^\top \bm x) \geq 0,
+
+
which can be reformatted to
+\bm x^\top \mathbf c \mathbf d^\top \bm x \geq 0,
+
+
which can further be symmetrized to
+\bm x^\top \left(\frac{\mathbf c \mathbf d^\top + \mathbf d \mathbf c^\top}{2}\right) \bm x \geq 0.
+
+
More general sets (general polyhedra, noncentered ellipsoids) can also be modelled using LMI too…
+
+
+
+
Combining the subset characterization and Lyapunov-ness using the S-procedure
+
. . .
+
+\alpha_i \bm x^\top \mathbf I \bm x \leq \bm x^\top \mathbf P_i \bm x,
+
The boundary between the regions \Omega_i and \Omega_j described by
+\Omega_{ij} = \{\bm x \in \mathbb R^n \mid \mathbf F_{ij} \bm z + \mathbf{l}_{ij}\},
+ where \bm z\in \mathbb R^p, \mathbf F_{ij}\in \mathbb R^{n\times p}, and \mathbf l_{ij}\in \mathbb R^{n}.
+
The continuity conditions are
+V_i(\bm x) = V_j(\bm x) \quad \forall \bm x \in \Omega_{ij},
+
+
which can be reformulated as
+\begin{aligned}
+&\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \mathbf P_i \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right) \\
+& \qquad + 2\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \bm q_i + r_i \\
+&= \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \mathbf P_j \left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right) \\
+& \qquad + 2\left(\mathbf F_{ij} \bm z + \mathbf{l}_{ij}\right)^\top \bm q_j + r_j
+\end{aligned}.
+
The topic of verification of hybrid systems is vast. While we only reserved a single week/chapter/block for it, it would easily fill a dedicated course, supported by a couple of books. Having a smaller time budget, we have still found some encouragement that modest introduction is feasible in the Chapter 3 of [1]. Although we are not following the book closely, we are covering some of their topics.
+
Among general references for hybrid system verification, we can recommend [2]. Although the book is not freely available for download, its web page contains quite some additional material such as slides and codes.
H. Lin and P. J. Antsaklis, Hybrid Dynamical Systems: Fundamentals and Methods. in Advanced Textbooks in Control and Signal Processing. Cham: Springer, 2022. Accessed: Jul. 09, 2022. [Online]. Available: https://doi.org/10.1007/978-3-030-78731-8
M. Althoff, G. Frehse, and A. Girard, “Set Propagation Techniques for Reachability Analysis,”Annual Review of Control, Robotics, and Autonomous Systems, vol. 4, no. 1, pp. 369–395, 2021, doi: 10.1146/annurev-control-071420-081941.
S. Prajna and A. Jadbabaie, “Safety Verification of Hybrid Systems Using Barrier Certificates,” in Hybrid Systems: Computation and Control, R. Alur and G. J. Pappas, Eds., in Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 2004, pp. 477–492. doi: 10.1007/978-3-540-24743-2_32.
E. M. Clarke, Jr, O. Grumberg, D. Kroening, D. Peled, and H. Veith, Model Checking, 2nd ed. in Cyber Physical Systems Series. Cambridge, MA, USA: MIT Press, 2018. Available: https://mitpress.mit.edu/9780262038836/model-checking/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+