diff --git a/.nojekyll b/.nojekyll index d4c7d3f..5820e42 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -2a85fde5 \ No newline at end of file +465962d1 \ No newline at end of file diff --git a/classes_PWA 7.html b/classes_PWA 7.html new file mode 100644 index 0000000..9c84379 --- /dev/null +++ b/classes_PWA 7.html @@ -0,0 +1,1190 @@ + + + + + + + + + +Piecewise affine (PWA) systems – B(E)3M35HYS – Hybrid systems + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Piecewise affine (PWA) systems

+
+ + + +
+ + + + +
+ + + +
+ + +

This is a subclass of switched systems where the functions on the right-hand side of the differential equations are affine functions of the state. For some (historical) reason these systems are also called piecewise linear (PWL).

+

We are going to reformulate such systems as switched systems with state-driven switching.

+

First, we consider the autonomous case, that is, systems without inputs: +\dot{\bm x} += +\begin{cases} +\bm A_1 \bm x + \bm b_1, & \mathrm{if}\, \bm H_1 \bm x + \bm g_1 \leq 0,\\ +\vdots\\ +\bm A_m \bm x + \bm b_m, & \mathrm{if}\, \bm H_m \bm x + \bm g_m \leq 0. +\end{cases} +

+

The nonautonomous case of systems with inputs is then: +\dot{\bm x} += +\begin{cases} +\bm A_1 \bm x + \bm B_1 u + \bm c_1, & \mathrm{if}\, \bm H_1 \bm x + \bm g_1 \leq 0,\\ +\vdots\\ +\bm A_m \bm x + \bm B_m + \bm c_m, & \mathrm{if}\, \bm H_m \bm x + \bm g_m \leq 0. +\end{cases} +

+
+

Example 1 (Linear system with saturated linear state feedback) In this example we consider a linear system with a saturated linear state feedback as in Fig 1.

+
+
+
+ +
+
+Figure 1: Linear system with a saturated linear state feedback +
+
+
+

The state equations for the close-loop system are +\dot{\bm x} = \bm A\bm x + \bm b \,\mathrm{sat}(v), \quad v = \bm k^T \bm x, + which can be reformulated as a piecewise affine system +\dot{\bm x} = +\begin{cases} +\bm A \bm x - \bm b, & \mathrm{if}\, \bm x \in \mathcal{X}_1,\\ +(\bm A + \bm b \bm k^\top )\bm x, & \mathrm{if}\, \bm x \in \mathcal{X}_2,\\ +\bm A \bm x + \bm b, & \mathrm{if}\, \bm x \in \mathcal{X}_3,\\ +\end{cases} + where the partitioning of the space of control inputs is shown in Fig 2.

+
+
+
+ +
+
+Figure 2: Partitioning the control input space +
+
+
+

Expressed in the state space, the partitioning is +\begin{aligned} +\mathcal{X}_1 &= \{\bm x \mid \bm H_1\bm x + g_1 \leq 0\},\\ +\mathcal{X}_2 &= \{\bm x \mid \bm H_2\bm x + \bm g_2 \leq 0\},\\ +\mathcal{X}_3 &= \{\bm x \mid \bm H_3\bm x + g_3 \leq 0\}, +\end{aligned} + where +\begin{aligned} +\bm H_1 &= \bm k^\top, \quad g_1 = 1,\\ +\bm H_2 &= \begin{bmatrix}-\bm k^\top\\\bm k^\top\end{bmatrix}, \quad \bm g_2 = \begin{bmatrix}-1\\-1\end{bmatrix},\\ +\bm H_3 &= -\bm k^\top, \quad g_3 = 1. +\end{aligned} +

+
+
+

Approximation of nonlinear systems

+

While the example with the saturated linear state feedback can be modelled as a PWA system exactly, there are many practical cases, in which the system is not exactly PWA affine but we want to approximate it as such.

+
+

Example 2 (Nonlinear system approximated by a PWA system) Consider the following nonlinear system +\begin{bmatrix} +\dot x_1\\\dot x_2 +\end{bmatrix} += +\begin{bmatrix} +x_2\\ +-x_2 |x_2| - x_1 (1+x_1^2) +\end{bmatrix} +

+

Our task is to approximate this system by a PWA system. Equivalently, we need to find a PWA approximation for the right-hand side function.

+
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/classes_PWA.html b/classes_PWA.html index d68188d..9c84379 100644 --- a/classes_PWA.html +++ b/classes_PWA.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Switched systems

+
+ + + +
+ + + + +
+ + + +
+ + +

Switched systems are modelled by first-order differential (state) equations with multiple right-hand sides, that is,

+

+\dot{\bm x} = \bm f_q(\bm x), \qquad q \in \{1,2, \ldots, m\}, +\tag{1} where m right-hand sides are possible and the lower index q determines which right-hand side function is “active” at a given moment.

+

The question is, what dictates the evolution of the integer variable q? In other words, what drives the switching? It turns out that the switching can be time-driven or state-driven.

+

In both cases, the right-hand sides can also depend the control input \bm u.

+

Major results for switched systems have been achieved without the need to refer to the framework of hybrid systems. But now that we have built such general framework, it turns out useful to view switched systems as a special class of hybrid systems. The aspects in which they are special will be discussed in the following, but here let us state that in contrast to full hybrid systems, switched systems are a bit less rich on the discrete side.

+
+

Time-driven

+

The evolution of the state variable complies with the following model \dot{\bm x} = \bm f_{q(t)}(\bm x), where q(t) is some function of time. The values of q(t) can be under our control or beyond our control, deterministic or stochastic.

+

A hybrid automaton for a time-driven switched system is shown in Fig 1.

+
+
+
+ +
+
+Figure 1: An automaton for a switched system with time-driven switching +
+
+
+

The transition from one mode to another is triggered by the integer variable q(t) attaining the appropriate value.

+

Since the switching signal is unrelated to the (continuous) state of the system, the invariant of the two modes are usually covering the whole state space \mathcal X.

+
+
+

State-dependent switching

+

The model is +\dot{\bm x} += +\begin{cases} +\bm f_1(\bm x), & \mathrm{if}\, \bm x \in \mathcal{X}_1,\\ +\vdots\\ +\bm f_m(\bm x), & \mathrm{if}\, \bm x \in \mathcal{X}_m. +\end{cases} +

+

Let’s consider just two domains \mathcal X_1 and \mathcal X_2. A hybrid automaton for a state-driven switched system is shown in Fig 2.

+
+
+
+ +
+
+Figure 2: An automaton for a switched system with state-driven switching +
+
+
+

The transition to the other mode is triggered by the continuous state of the system crossing the boundary between the two domains. The boundary is defined by the function s(\bm x) (called switching function), which is zero on the boundary, see the Fig 3.

+
+
+
+ +
+
+Figure 3: State-dependent switching +
+
+
+

Through examples we now illustrate the possible behaviors of the system when the flow transverses the boundary, when it pulls away from the boundary, and when it pushes towards the boundary.

+
+

Example 1 (The flow transverses the boundary) We consider the two right-hand sides of the state equation +\bm f_1(\bm x) = \begin{bmatrix}1\\ x_1^2 + 2x_2^2\end{bmatrix} + and +\bm f_2(\bm x) = \begin{bmatrix}1\\ 2x_1^2+3x_2^2-2\end{bmatrix} + and the switching function +s(x_1,x_2) = (x_1+0.05)^2 + (x_2+0.15)^2 - 1. +

+

The state portrait that also shows the switching function is generated using the following code.

+
+
+Show the code +
s(x₁,x₂) = (x₁+0.05)^2 + (x₂+0.15)^2 - 1.0
+
+f₁(x₁,x₂) = x₁^2 + 2x₂^2
+f₂(x₁,x₂) = 2x₁^2+3x₂^2-2.0
+
+f(x₁,x₂) = s(x₁,x₂) <= 0.0 ? [1,f₁(x₁,x₂)] : [1,f₂(x₁,x₂)] 
+
+N = 100
+x₁ = range(0, stop = 0.94, length = N)
+
+using CairoMakie
+fig = Figure(size = (600, 600),fontsize=20)
+ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
+streamplot!(ax,(x₁,x₂)->Point2f(f(x₁,x₂)), 0..1.5, 0..1.5, colormap = :magma)
+lines!(ax,x₁,sqrt.(1 .- (x₁ .+ 0.05).^2) .- 0.15, color = :red, linewidth=5)
+x10 = 0.5
+x20 = sqrt(1 - (x10 + 0.05)^2) - 0.15
+Makie.scatter!(ax,[x10],[x20],color=:blue,markersize=30)
+fig
+
+
+
+
+

+
+
+
+
+

The state portrait also shows a particular initial state \bm x_0 using a blue dot. Note that the projection of both vector fields \mathbf f_1 and \mathbf f_2 evaluated at \bm x_0 onto the normal (the gradient) of the switching function at \bm x_0 is positive, that is +\left.\left(\nabla s\right)^\top \bm f_1\right|_{\bm x_0} \geq 0, \quad \left.\left(\nabla s\right)^\top \bm f_2\right|_{\bm x_0} \geq 0. +

+

This is consistent with the observation that the flow goes through the boundary.

+

We can also plot a particular solution of the ODE using the following code.

+
+
+Show the code +
using DifferentialEquations
+F(u, p, t) = f(u[1],u[2])
+u0 = [0.0,0.4]
+tspan = (0.0, 1.0)
+prob = ODEProblem(F, u0, tspan)
+sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)
+
+using Plots
+Plots.plot(sol,lw=3,xaxis="Time",yaxis="x",label=false)
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+

Strictly speaking, this solution does not satisfy the differential equation on the boundary of the two domains (the derivative of x_2 does not exist there). This is visually recognized in the above plot as the sharp corner in the solution. But other than that, the solution is perfectly “reasonable” – for a while the system evolves according to one state equations, then at one particular moment the system starts evolving according to another state equation. That is it. Not much more to see here.

+
+
+

Example 2 (The flow pulls away from the boundary) We now consider another pair of the right-hand sides. +\bm f_1(\bm x) = \begin{bmatrix}-1\\ x_1^2 + 2x_2^2\end{bmatrix} + and +\bm f_2(\bm x) = \begin{bmatrix}1\\ 2x_1^2+3x_2^2-2\end{bmatrix}. +

+

The switching function is the same as in the previous example.

+

The state portrait is below.

+
+
+Show the code +
f(x₁,x₂) = s(x₁,x₂) <= 0.0 ? [-1, f₁(x₁,x₂)] : [1, f₂(x₁,x₂)] 
+
+fig = Figure(size = (600, 600),fontsize=20)
+ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
+streamplot!(ax,(x₁,x₂)->Point2f(f(x₁,x₂)), 0..1.5, 0..1.5, colormap = :magma)
+lines!(ax,x₁,sqrt.(1 .- (x₁ .+ 0.05).^2) .- 0.15, color = :red, linewidth=5)
+x10 = 0.8
+x20 = sqrt(1 - (x10 + 0.05)^2) - 0.15
+Makie.scatter!(ax,[x10],[x20],color=:blue,markersize=30)
+fig
+
+
+
+
+

+
+
+
+
+

We focus on the blue dot again. The projections of the two vector fields onto the normal of the switching function satisfy +\left.\left(\nabla s\right)^\top \bm f_1\right|_{\bm x_0} \leq 0, \quad \left.\left(\nabla s\right)^\top \bm f_2\right|_{\bm x_0} \geq 0. +

+

The only interpretation of this situation is that a unique solution does not start at \bm x_0. Again, not much more to see here.

+
+
+

Example 3 (The flow pushes towards the boundary) And one last pair of the right-hand sides: +\bm f_1(\bm x) = \begin{bmatrix}1\\ x_1^2 + 2x_2^2\end{bmatrix} + and +\bm f_2(\bm x) = \begin{bmatrix}-1\\ 2x_1^2+3x_2^2-2\end{bmatrix}. +

+

The state-portrait is below.

+
+
+Show the code +
f(x₁,x₂) = s(x₁,x₂) <= 0.0 ? [1, f₁(x₁,x₂)] : [-1, f₂(x₁,x₂)] 
+
+fig = Figure(size = (600, 600),fontsize=20)
+ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
+streamplot!(ax,(x₁,x₂)->Point2f(f(x₁,x₂)), 0..1.5, 0..1.5, colormap = :magma)
+lines!(ax,x₁,sqrt.(1 .- (x₁ .+ 0.05).^2) .- 0.15, color = :red, linewidth=5)
+x10 = 0.5
+x20 = sqrt(1 - (x10 + 0.05)^2) - 0.15
+Makie.scatter!(ax,[x10],[x20],color=:blue,markersize=30)
+fig
+
+
+
+
+

+
+
+
+
+

The projections of the two vector fields onto the normal of the switching function satisfy +\left.\left(\nabla s\right)^\top \bm f_1\right|_{\bm x_0} \geq 0, \quad \left.\left(\nabla s\right)^\top \bm f_2\right|_{\bm x_0} \leq 0. +

+

But this is interesting! Once the trajectory hits the switching curve and tries to penetrate it futher, it is pushed back to the switching curve. As it tries to penetrate it further, it is pushed back to the switching curve again. And so on. But then, how does the state evolve from \bm x_0?

+

Hint: solve the ODE numerically with some finite step size. The solution will exhibit zig-zagging or chattering along the switching curve, away from the blue point. Now, keep shrinking the step size. The solution will ultimately “slide” smoothly along the switching curve. Perhaps this was your guess. One thing should worry you, however: such “sliding” solution satisfies neither of the two state equations!

+

We will make this more rigorous in a moment, but right now we just wanted to tease the intuition.

+
+
+
+

Conditions for existence and uniqueness of solutions of ODE

+

In order to analyze the situations such as the previous example, we need to recapitulate some elementary facts about the existence and uniqueness of solutions of ordinary differential equations (ODEs). And then we are going to add some new stuff.

+

Consider the ODE

+

\dot x(t) = f(x(t),t).

+

We ask the following two questions:

+
    +
  • Under which conditions does a solution exists?
  • +
  • Under which conditions is the solution unique?
  • +
+

To answer both, the function f() must be analyzed.

+

But before we answer the two questions, we must ask another one that is even more fundamental:

+
    +
  • What does it mean that a function x(t) is a solution of the the ODE?
  • +
+

However trivial this question may seem, an answer can escalate rather quickly – there are actually several concepts of a solution of an ordinary differential equation.

+
+

Classical solution (Peano, also Cauchy-Peano)

+

We assume that f(x(t),t) is continuous with respect to both x and t. Then existence of a solution is guaranteed locally (on some finite interval), but uniqueness is not.

+
+
+
+ +
+
+Not guaranteed does not mean impossible +
+
+
+

Uniqueness is not not excluded in all cases, it is just that it is not guaranteed.

+
+
+

A solution is guaranteed to be continuously differentiable ( x\in\mathrm C^1 ). Such function x(t) satisfies the ODE \dot x(t) = f(x(t),t) \; \forall t, that is why such solution is called classical.

+
+

Example 4 An example of a solution that exists only on a finite interval is + \dot x(t) = x^2(t),\; x(0) = 1, +
+for which the solution is x(t) = \frac{1}{1-t} . The solution blows up at t=1 .

+
+
+

Example 5 An example of nonuniqueness is provided by \dot x(t) = \sqrt{x(t)}, \; x(0) = 0.

+

One possible solution is x(t) = \frac{1}{4}t^2. Another is x(t) = 0. Yet another example is x(t) = \frac{1}{4}(t-t_0)^2. It is related to the Leaky bucket example.

+
+
+
+

Strenghtening the requirement of continuity (Pickard-Lindelöf)

+

Since continuity of f(x(t),t) was not enough to guarantee uniqueness, we need to impose a stricter condition on f(). Namely, we impose a stricter condition on f() with respect to x – Lipschitz continuity, while we still require that the function be continuous with respect to t.

+

Now it is not only existence but also uniqueness of a solution that is guaranteed.

+
+
+
+ +
+
+Uniqueness not guaranteed does not mean it is impossible +
+
+
+

Similarly as with Peano conditions, here too the condition is not necessary, it is just sufficient – even if the function f is not Lipschitz continuous, there may exist a unique solution.

+
+
+

Since the condition is stricter than mere continuity, whatever goodies hold here too. In particular, the solution is guaranteed to be continuously differentiable.

+

If the function is only locally Lipschitz, the solution is guaranteed on some finite interval. If the function is (globally) Lipschitz, the solution is guaranteed on an unbounded interval.

+
+
+

Extending the set of solutions (Carathéodory)

+

In contrast with the classical solution, we can allow the solution x(t) to fail to satisfy the ODE at some isolated points in time. This is called Carathéodory (or extended) solution.

+

Carathéodory solution x(t) is more than just continuous (even more than uniformly continuous) but less than contiuously differentiable (aka \mathcal C^1) – it is absolutely continuous. Absolutely continuous function is a solution of the integral equation (indeed, an equation) x(t) = x(t_0) + \int_{t_0}^t f(x(\tau),\tau)\mathrm{d}\tau,
+where we use Lebesgue integral (instead of Riemann).

+

Having referred to absolute continuity and Lebesgue integral, the discussion could quickly become rather technical. But all we want to say is that f can be “some kind of discontinuous” with respect to t. In particular, it must be measurable wrt t, which again seems to start escalating… But it suffices to say that it includes the case when f(x,t) is piecewise continuous with respect to t (sampled data control with ZOH).

+

Needles to say that for a continuous f, solutions x are just classical (smooth).

+

If the function f is discontinuous with respect to x, some more concepts of a solution need to be invoked so that existence and uniqueness can be analyzed.

+
+

Example 6 (Some more examples of nonexistence and nonuniqueness of solutions) The system with a discontinuous RHS +\begin{aligned} +\dot x_1 &= -2x_1 - 2x_2\operatorname*{sgn(x_1)},\\ +\dot x_2 &= x_2 + 4x_1\operatorname*{sgn(x_1)} +\end{aligned} + can be reformulated as a switched system +\begin{aligned} +\dot{\bm x} &= \begin{bmatrix}-2 & 2\\-4 & 1\end{bmatrix}\begin{bmatrix}x_1\\ x_2\end{bmatrix}, \quad s(\bm x)\leq 0\\ +\dot{\bm x} &= \begin{bmatrix}-2 & -2\\4 & 1\end{bmatrix}\begin{bmatrix}x_1\\ x_2\end{bmatrix}, \quad s(\bm x)> 0, +\end{aligned} + where the switching function is s(\bm x) = x_1.

+
+
+Show the code +
s(x) = x[1]
+
+f₁(x) = [-2x[1] + 2x[2], x[2] - 4x[1]]
+f₂(x) = [-2x[1] - 2x[2], x[2] + 4x[1]]
+
+f(x) = s(x) <= 0.0 ? f₁(x) : f₂(x) 
+
+using CairoMakie
+fig = Figure(size = (600, 600),fontsize=20)
+ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
+streamplot!(ax,x->Point2f(f(x)), -1.5..1.5, -1.5..1.5, colormap = :magma)
+vlines!(ax,0; ymin = -1.1, ymax = 1.1, color = :red)
+fig
+
+
+
+
+

+
+
+
+
+
+
+
+

Sliding mode dynamics (on simple boundaries)

+

The previous example provided yet another illustration of a phenomenon of sliding, or a sliding mode. We say that there is an attractive sliding mode at \bm x_\mathrm{s}, if there is a trajectory that ends at \bm x_\mathrm{s}, but no trajectory that starts at \bm x_\mathrm{s}.

+
+
+

Generalized solutions (Filippov)

+

It is now high time to introduce yet another concept of a solution. A concept that will make it possible to model the sliding mode dynamics in a more rigorous way. Remember that when the state \bm x(t) slides along the boundary, it qualifies as a solution to neither of the two state equations in any sense we have discussed so far. But now comes the concept of Fillipov solution.

+

\bm x() is a Filippov solution on [t_0,t_1] if for almost all t +\dot{\bm{x}}(t) \in \overline{\operatorname*{co}}\{\mathbf f(\bm x(t),t)\}, + where \overline{\operatorname*{co}} denotes the (closed) convex hull.

+
+

Example 7 Consider the model in the previous example. The switching surface, along which the solution slides, is given by \mathcal{S}^+ = \{\bm x \mid x_1=0 \land x_2\geq 0\}.

+

Now, Filippov solution must satisfy the following differential inclusion +\begin{aligned} +\dot{\bm x}(t) &\in \overline{\operatorname*{co}}\{\bm A_1\bm x_1(t), \bm A_2\bm x_2(t)\}\\ +&= \alpha_1(t) \bm A_1\bm x_1(t) + \alpha_2(t) \bm A_2\bm x_2(t), +\end{aligned} + where \alpha_1(t), \alpha_2(t) \geq 0, \alpha_1(t) + \alpha_2(t) = 1.

+

Note, however, that not all the weights keep the solution on \mathcal S^+. We must impose some restriction, namely that \dot x_1 = 0 for \bm x(t) \in \mathcal S^+. This leads to +\alpha_1(t) [-2x_1 + 2x_2] + \alpha_2(t) [-2x_1 - 2x_2] = 0 +

+

Combining this with \alpha_1(t) + \alpha_2(t) = 1 gives +\alpha_1(t) = \alpha_2(t) = 1/2, + which in this simple case perhaps agrees with our intuition (the average of the two vector fields).

+

The dynamics on the sliding mode is modelled by +\dot x_1 = 0, \quad \dot x_2 = x_2, \quad \bm x \in \mathcal{S}^+. +

+
+
+
+
+ +
+
+Possible nonuniqueness on intersection of boundaries +
+
+
+

+
+
+ + +
+
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/classes_switched.html b/classes_switched.html index c7ca456..299338d 100644 --- a/classes_switched.html +++ b/classes_switched.html @@ -2,7 +2,7 @@ - + @@ -71,10 +71,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Literature

+
+ + + +
+ + + + +
+ + + +
+ + +

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].

+ + + + + Back to top

References

+
+
[1]
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.
+
+
+
[2]
M. Heemels, Linear complementarity systems: a study in hybrid dynamics,” PhD thesis, Technische Universiteit Eindhoven, Eindhoven, NL, 1999. Available: https://heemels.tue.nl/content/papers/Hee_TUE99a.pdf
+
+
+
[3]
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, LCQPowA 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_references.html b/complementarity_references.html index 99ac8cc..c8ab60a 100644 --- a/complementarity_references.html +++ b/complementarity_references.html @@ -2,7 +2,7 @@ - + @@ -57,10 +57,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Simulations using complementarity

+
+ + + +
+ + + + +
+ + + +
+ + +
+

Simulation using time-stepping

+

+\begin{aligned} +\dot x_1 &= -\operatorname{sign} x_1 + 2 \operatorname{sign} x_2\\ +\dot x_2 &= -2\operatorname{sign} x_1 -\operatorname{sign} x_2 +\end{aligned} +

+
+
+Show the code +
f₁(x) = -sign(x[1]) + 2*sign(x[2])
+f₂(x) = -2*sign(x[1]) - sign(x[2])
+f(x) = [f₁(x), f₂(x)]
+
+using CairoMakie
+fig = Figure(resolution = (800, 800),fontsize=20)
+ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
+streamplot!(ax,(x₁,x₂)->Point2f(f([x₁,x₂])), -2.0..2.0, -2.0..2.0, colormap = :magma)
+fig
+
+
+
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
+└ @ Makie ~/.julia/packages/Makie/pFPBw/src/scenes.jl:238
+
+
+
+
+

+
+
+
+
+
+
+Show the code +
using DifferentialEquations
+
+function f!(dx,x,p,t)
+    dx[1] = -sign(x[1]) + 2*sign(x[2])
+    dx[2] = -2*sign(x[1]) - sign(x[2])
+end
+
+x0 = [-1.0, 1.0]
+tfinal = 2.0
+tspan = (0.0,tfinal)
+prob = ODEProblem(f!,x0,tspan)
+sol = solve(prob)
+
+using Plots
+Plots.plot(sol,xlabel="t",ylabel="x",label=false,lw=3)
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+Show the code +
Plots.plot(sol[1,:],sol[2,:],xlabel="x₁",ylabel="x₂",label=false,aspect_ratio=:equal,lw=3,xlims=(-1.2,0.5))
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+

Example continued: how fast does the solution approach the origin?

+
    +
  • Let’s use the 1-norm \|\bm x\|_1 = |x_1| + |x_2| to measure how far the trajectory is from the origin.
  • +
  • How fast does the trajectory converge to the origin? That, is +\frac{\mathrm d}{\mathrm dt}\|\bm x\|_1 = ? +
  • +
  • Consider each quadrant separately. Let’s start in the first (upper right) quadrant, that is, x_1>0 and x_2>0, and therefore |x_1| = x_1, \;|x_2| = x_2, and therefore
  • +
+

. . .

+

+\frac{\mathrm d}{\mathrm dt}\|\bm x\|_1 = \dot x_1 + \dot x_2 = 1 - 3 = -2. +

+
    +
  • Identical in the other quadrants. And undefined on the axes.

  • +
  • The trajectory will hit the origin in finite time:

    +
      +
    • For x_1(0) = 1 and x_2(0) = 1 , the trajectory hits the origin at t=(|x_1(0)|+|x_2(0)|)/2 = 1.
    • +
  • +
  • But with an infinite number of revolutions around the origin…

  • +
  • How will a standard algoritm for numerical simulation handle this?

  • +
+
+
+

Forward Euler with fixed step size

+

+\begin{aligned} +{\color{blue}x_{1,k+1}} &= x_{1,k} + h (-\operatorname{sign} x_{1,k} + 2 \operatorname{sign} x_{2,k})\\ +{\color{blue}x_{2,k+1}} &= x_{1,k} + h (-2\operatorname{sign} x_{1,k} - \operatorname{sign} x_{2,k}) +\end{aligned} +

+
+
+Show the code +
f(x) = [-sign(x[1]) + 2*sign(x[2]), -2*sign(x[1]) - sign(x[2])]
+
+using LinearAlgebra
+N = 1000
+x₀ = [-1.0, 1.0]    
+x = [x₀]
+tfinal = norm(x₀,1)/2
+tfinal = 5.0
+h = tfinal/N 
+t = range(0.0, step=h, stop=tfinal)
+
+for i=1:N
+    xnext = x[i] + h*f(x[i]) 
+    push!(x,xnext)
+end
+
+X = [getindex.(x, i) for i in 1:length(x[1])]
+
+Plots.plot(t,X,lw=3,label=false,xlabel="t",ylabel="x")
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+

Backward Euler

+

+\begin{aligned} +{\color{blue} x_{1,k+1}} &= x_{1,k} + h (-\operatorname{sign} {\color{blue}x_{1,k+1}} + 2 \operatorname{sign} {\color{blue}x_{2,k+1}})\\ +{\color{blue} x_{2,k+1}} &= x_{1,k} + h (-2\operatorname{sign} {\color{blue}x_{1,k+1}} - \operatorname{sign} {\color{blue}x_{2,k+1}}) +\end{aligned} +

+
+
+

Formulation using LCP

+
    +
  • Instead solving the above nasty equations, introduce new variables u_1 and u_2 as the outcomes of the \operatorname{sign} functions: +\begin{aligned} +{\color{blue} x_{1,k+1}} &= x_{1,k} + h (-{\color{blue}u_{1}} + 2 {\color{blue}u_{2}})\\ +{\color{blue} x_{2,k+1}} &= x_{1,k} + h (-2{\color{blue}u_{1}} - {\color{blue}u_{2}}) +\end{aligned} +

  • +
  • But now we have to enforce the relation between \bm u and \bm x_{k+1}.

  • +
  • Recall the standard definition of the \operatorname{sign} function: +\operatorname{sign}(x) = \begin{cases} +1 & x>0\\ +0 & x=0\\ +-1 & x<0 +\end{cases} +

  • +
  • Change the definition to a set-valued function +\begin{cases} +\operatorname{sign}(x) = 1 & x>0\\ +\operatorname{sign}(x) \in [-1,1] & x=0\\ +\operatorname{sign}(x) = -1 & x<0 +\end{cases} +

  • +
  • Accordingly, set the relation between \bm u and \bm x +\begin{cases} +u_1 = 1 & x_1>0\\ +u_1 \in [-1,1] & x_1=0\\ +u_1 = -1 & x_1<0 +\end{cases} + and +\begin{cases} +u_2 = 1 & x_2>0\\ +u_2 \in [-1,1] & x_2=0\\ +u_2 = -1 & x_2<0 +\end{cases} +

  • +
  • But these are mixed complementarity contraints! +\begin{aligned} +\begin{bmatrix} +{\color{blue} x_{1,k+1}}\\ +{\color{blue} x_{1,k+1}} +\end{bmatrix} +&= +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} + h +\begin{bmatrix} +-1 & 2 \\ +-2 & -1 +\end{bmatrix} +\begin{bmatrix} +{\color{blue}u_{1}}\\ +{\color{blue}u_{2}} +\end{bmatrix}\\ +-1 \leq {\color{blue} u_1} \leq 1 \quad &\bot \quad -{\color{blue}x_{1,k+1}}\\ +-1 \leq {\color{blue} u_2} \leq 1 \quad &\bot \quad -{\color{blue}x_{2,k+1}} +\end{aligned} +

  • +
+
+
+

9 possible combinations

+
    +
  • Let’s explore some: x_{1,k+1} = x_{2,k+1} = 0, while u_1 \in [-1,1] and u_2 \in [-1,1]:
  • +
+

. . .

+

+\begin{aligned} +\begin{bmatrix} +0\\ +0 +\end{bmatrix} +&= +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} + h +\begin{bmatrix} +-1 & 2 \\ +-2 & -1 +\end{bmatrix} +\begin{bmatrix} +{\color{blue}u_{1}}\\ +{\color{blue}u_{2}} +\end{bmatrix}\\ +& -1 \leq {\color{blue} u_1} \leq 1, \quad -1 \leq {\color{blue} u_2} \leq 1 +\end{aligned} +

+
    +
  • How does the set of states from which the next state is zero look like? +\begin{aligned} +-\begin{bmatrix} +-1 & 2 \\ +-2 & -1 +\end{bmatrix}^{-1} +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} +&= h +\begin{bmatrix} +{\color{blue}u_{1}}\\ +{\color{blue}u_{2}} +\end{bmatrix}\\ +-1 \leq {\color{blue} u_1} \leq 1, \quad -1 &\leq {\color{blue} u_2} \leq 1 +\end{aligned} +
  • +
+

. . .

+

+\begin{bmatrix} +-h\\-h +\end{bmatrix} +\leq +\begin{bmatrix} +0.2 & 0.4 \\ +-0.4 & 0.2 +\end{bmatrix} +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} +\leq +\begin{bmatrix} +h\\ h +\end{bmatrix} +

+
    +
  • For h=0.2
  • +
+
+
+Show the code +
using LazySets
+h = 0.2
+H1u = HalfSpace([0.2, 0.4], h)
+H2u = HalfSpace([-0.4, 0.2], h)
+H1l = HalfSpace(-[0.2, 0.4], h)
+H2l = HalfSpace(-[-0.4, 0.2], h)
+
+Ha = H1u  H2u  H1l  H2l
+
+using Plots
+Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
    +
  • Indeed, if the current state is in this rotated square, then the next state will be zero.
  • +
+
+
+

Another

+
    +
  • u_1 = 1, u_2 = 1:
  • +
+

. . .

+

+\begin{aligned} +\begin{bmatrix} +{\color{blue} x_{1,k+1}}\\ +{\color{blue} x_{1,k+1}} +\end{bmatrix} +&= +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} + h +\begin{bmatrix} +-1 & 2 \\ +-2 & -1 +\end{bmatrix} +\begin{bmatrix} +{1}\\ +{1} +\end{bmatrix}\\ +\color{blue}x_{1,k+1} &\geq 0\\ +\color{blue}x_{2,k+1} &\geq 0 +\end{aligned} +

+
    +
  • which can be reformatted to +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} + h +\begin{bmatrix} +-1 & 2 \\ +-2 & -1 +\end{bmatrix} +\begin{bmatrix} +1\\ +1 +\end{bmatrix}\geq \bm 0 +

  • +
  • and further to +\begin{bmatrix} +x_{1,k}\\ +x_{2,k} +\end{bmatrix} +\geq h +\begin{bmatrix} +-1\\ +3 +\end{bmatrix} +

  • +
+

. . .

+
+
+Show the code +
using LazySets
+h = 0.2
+A = [-1.0 2.0; -2.0 -1.0]
+u = [1.0, 1.0]
+b = h*A*u
+
+H1 = HalfSpace([-1.0, 0.0], b[1])
+H2 = HalfSpace([0.0, -1.0], b[2])
+Hb = H1  H2
+
+using Plots
+Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
+Plots.plot!(Hb)
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+

All nine regions

+

. . .

+
+
+Show the code +
using LazySets
+h = 0.2
+A = [-1.0 2.0; -2.0 -1.0]
+
+u = [1, -1]
+b = h*A*u
+
+H1 = HalfSpace(-[1.0, 0.0], b[1])
+H2 = HalfSpace([0.0, 1.0], -b[2])
+Hc = H1  H2
+
+u = [-1, 1]
+b = h*A*u
+
+H1 = HalfSpace([1.0, 0.0], -b[1])
+H2 = HalfSpace(-[0.0, 1.0], b[2])
+Hd = H1  H2
+
+u = [-1, -1]
+b = h*A*u
+
+H1 = HalfSpace([1.0, 0.0], -b[1])
+H2 = HalfSpace([0.0, 1.0], -b[2])
+He = H1  H2
+
+using Plots
+Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
+Plots.plot!(Hb)
+Plots.plot!(Hc)
+Plots.plot!(Hd)
+Plots.plot!(He)
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+

Solutions using a MCP solver

+
+
+Show the code +
M = [-1 2; -2 -1]
+h = 2e-1
+tfinal = 2.0
+N = tfinal/h
+
+x0 = [-1.0, 1.0]
+x = [x0]
+
+using JuMP
+using PATHSolver
+
+for i = 1:N
+    model = Model(PATHSolver.Optimizer)
+    set_optimizer_attribute(model, "output", "no")
+    set_silent(model)
+    @variable(model, -1 <= u[1:2] <= 1)
+    @constraint(model, -h*M * u - x[end]  u)
+    optimize!(model)
+    push!(x, x[end]+h*M*value.(u))
+end
+
+t = range(0.0, step=h, stop=tfinal)
+X = [getindex.(x, i) for i in 1:length(x[1])]
+
+using Plots
+Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
+Plots.plot!(Hb)
+Plots.plot!(Hc)
+Plots.plot!(Hd)
+Plots.plot!(He)
+Plots.plot!(X[1],X[2],xlabel="x₁",ylabel="x₂",label="Time-stepping",aspect_ratio=:equal,lw=3,markershape=:circle)
+Plots.plot!(sol[1,:],sol[2,:],label=false,lw=3)
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/complementarity_simulations.html b/complementarity_simulations.html index 80f4f29..5526a53 100644 --- a/complementarity_simulations.html +++ b/complementarity_simulations.html @@ -2,7 +2,7 @@ - + @@ -71,10 +71,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Software

+
+ + + +
+ + + + +
+ + + +
+ + +
+

Solving optimization problems with complementarity constraints

+

Surprisingly, there are not many software packages that can handle complementarity constraints directly.

+
    +
  • Within the realm of free software, I am only aware of PATH solver. Well, it is not open source and it is not issued under any classical free and open source license. It can be interfaced from Matlab and Julia (and AMPL and GAMS, which are not relevant for our course). For Matlab, compiled mexfiles can be downloaded. For Julia, the solver can be interfaced directly from the popular JuMP package (choosing the PATHSolver.jl solver), see the section on Complementarity constraints and Mixed complementarity problems in JuMP manual.
  • +
+

When restricted to Matlab, there are several options, all of them commercial:

+ +

Gurobi does not seem to support complementarity constraints.

+

Mosek supports disjunctive constraints, within which complementarity constraints can be formulated. But they are then approached using a mixed-integer solver.

+
+
+

Modeling and simulation of dynamical systems with complementarity constraints

+

Within the modeling and simulation domains, there are two free and open source libraries that can handle complementarity constraints, mainly motivated by nonsmooth dynamical systems:

+
    +
  • SICONOS +
      +
    • C++, Python
    • +
    • physical domain independent
    • +
  • +
  • PINOCCHIO +
      +
    • C++, Python
    • +
    • specialized for robotics
    • +
  • +
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/complementarity_software.html b/complementarity_software.html index 3fd686f..9aeaa38 100644 --- a/complementarity_software.html +++ b/complementarity_software.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Discrete-event systems

+
+ + + +
+ + + + +
+ + + +
+ + +

We have already mentioned that hybrid systems are composed of time-driven subsystems and event-driven subsystems. Assuming that the primary audience of this course are already familiar with the former (having been exposed to state equations and transfer function), here we are going to focus on the latter, also called discrete-event systems DES (or DEVS).

+
+

(Discrete) event

+

We need to start with the definition of an event. Sometimes an adjective discrete is added (to get discrete event), although it appears rather redundant.

+

The primary characteristic of an event is instantaneous occurence, that is, an event takes no time.

+

Within the context of systems, an event is associated with a change of state – transition of the system from one state to another. Between two events, the systems remains in the same state, it doesn’t evolve.

+
+
+
+ +
+
+The concept of a state +
+
+
+

True, here we are making a reference to the concept of a state, which we haven’t defined yet. But we can surely rely on understanding this concept developed by studying the time-driven systems (modelled by state equations).

+
+
+

Although it is not instrumental in defining the event, the state space is frequently assumed discrete (even if infinite).

+
+

Example 1 (DES state trajectory) In the figure below we can see an example state trajectory of a discrete-event system corresponding to a particular sequence of events.

+
+
+
+ +
+
+Figure 1: Example of a state trajectory in response to a sequence of events +
+
+
+

It is perhaps worth emphasizing that the state space is not necessarily equidistantly discretized.

+

Also note that for some events no transitions occur (e_3 at t_3).

+
+
+
+
+ +
+
+Frequent notational confusion: does the lower index represent discrete time or an element of a set? +
+
+
+

The previous example also displays one particular annoying (and frequently occuring in literature) notational conflict. How shall we interpret the lower index? Sometimes it is used to refer to (discrete) time, and on some other occasions it can just refer to a particular element of a set. In other words, this is a notational clash between name of the variable and value of the variable. In the example above we obviously adopted the latter interpretation. But in other cases, we ourselves are prone to inconsistency. Just make sure you understand what the author means.

+
+
+
+

Example 2 (State trajectory of a continuous-time dynamical systems) Compare now the above example of a state trajectory in a DES with the example of a continuous-time state space system below, whose model could be \dot x(t) = f(x). In the latter, any change, however small, takes time. In other words, the system evolves continuously in time.

+
+
+
+ +
+
+Figure 2: Example of a state trajectory of a continuous-time continuous-valued dynamical system +
+
+
+

The set of states (aka state space) is \mathbb{R} (or a subset) in this case (in general \mathbb{R}^n or a subset).

+
+
+

Example 3 (State trajectory of a time-discretized (aka sampled-data) system) As a yet another example of a state trajectory, consider the response of a discrete-time (actually time-discretized or also sampled-data system) system model by x_k = f(x_k) in the figure below. Although we could view the sampling instances as the events, these are given by time, hence the moments of transitions are predictable. Hence the system can still be viewed and analyzed as a time-driven and not event driven one.

+
+
+
+ +
+
+Figure 3: Example of a state trajectory of time-discretized (aka sampled data) system +
+
+
+
+
+
+

When do events occur?

+

There are three major possibilities:

+
    +
  • when action is taken (button press, clutch released, …),
  • +
  • spontaneously: well, this is just an “excuse” when the reason is difficult to trace down (computer failure, …),
  • +
  • when some condition is met (water level is reached, …). This needs an introduction of a concept of a hybrid systems, wait for it.
  • +
+
+
+

Sequence of “time-stamped” events (aka timed trace)

+

The sequence of pairs (event, time) (e_1,t_1), (e_2,t_2), (e_3,t_3), \ldots

+

is sufficient to characterize an execution of a deterministic system, that is, a system with a unique initial state and a unique transitions at a given state and an event.

+
+
+

DES can be stochastic, but what exactly is stochastic then?

+

Stochasticity can be introduced in

+
    +
  • the event times (e.g. Poisson process),
  • +
  • but also in the transitions (e.g. probabilistic automata, more on this later).
  • +
+
+
+

Sometimes time stamps not needed – the ordering of events is enough

+

The sequence of events (aka trace) e_1,e_2,e_3, \ldots can be enough for some analysis, in which only the order of the events is important.

+
+

Example 4 credit_card_swiped, pin_entered, amount_entered, money_withdrawn

+
+
+
+

Discrete-event systems are studied through their languages

+

When studying discrete-event systems, soon we are exposed to terminology from the formal language theory such as alphabet, word, and language. This must be rather confusing for most students (at least those with no education in computer science). In our course we are not going to use these terms actively (after all our only motivation for introducing the discipline of discrete-event systems is to take away just a few concepts that are useful in hybrid systems), but we want to sketch the motivation for their introduction to the discipline, which may make it easier for a student to skim through some materials on discrete-event systems.

+

We define at least those three terms that we have just mentioned. The definitions correspond to the everyday usage of these terms.

+
+
Alphabet
+
+a set of symbols. +
+
Word (or string)
+
+a sequence of symbols from a finite alphabet. +
+
Language
+
+a set of words from the given alphabet. +
+
+

Now, a symbol is used to label an event. Alphabet is then a set of possible events. A particular sequence of events (we also call it trace) is then represented by a word. Since we agreed that events are associated with state transitions of a corresponding system, a word represents a possible execution or run of a system. A set of all possible executions of a given system can then be formally viewed as language.

+

Indeed, all this is just a formalism, the agreement how to talk about things. We will see an example of this “jargon” in the next section when we introduce the concept of an automaton and some of its properties.

+
+
+

Modelling frameworks for DES (as used in our course)

+

These are the three frameworks that we are going to cover in our course. There may be some more, but these three are the major ones, and from these there is always some lesson to be learnt that we will find useful later when finally studying hybrid systems

+
    +
  • State automaton (pl. automata)
  • +
  • Petri net
  • +
  • (max,plus) algebra, MPL systems
  • +
+

While the first two frameworks are essentially equally powerful when it comes to modelling DES, the third one can be regarded as an algebraic framework for a subset of systems modelled by Petri nets.

+

We are going to cover all the three frameworks in this course as a prequel to hybrid systems.

+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/des.html b/des.html index 0e4bd65..0a32970 100644 --- a/des.html +++ b/des.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Hybrid automata

+
+ + + +
+ + + + +
+ + + +
+ + +

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:

+ +

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

+ +
+
+
+ +
+
+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.

+
+
+ +
+

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\}, +

+

+\mathcal G(\text{off},\text{on}) = \{x\mid x\leq 17\},\; \mathcal G(\text{on},\text{off}) = \{x\mid x\geq 19\}. +

+

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:

+
    +
  1. 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.
  2. +
  3. 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.
  4. +
+
+
+

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 +
+
+
+

+\mathcal{Q} = \{\mathrm{off}, \mathrm{dim}, \mathrm{bright}\},\quad \mathcal{Q}_0 = \{\mathrm{off}\} +

+

+\mathcal{X} = \mathbb{R}, \quad \mathcal{X}_0 = \{0\} +

+

+\mathcal{I(\mathrm{off})} = \mathcal{I(\mathrm{bright})} = \mathcal{I(\mathrm{dim})} = \{x\in\mathbb R \mid x \geq 0\} +

+

+f(x) = 1 +

+

+\mathcal{A} = \{\mathrm{press}\} +

+

+\begin{aligned} +\mathcal{E} &= \{(\mathrm{off},\mathrm{press},\mathrm{dim}),(\mathrm{dim},\mathrm{press},\mathrm{off}),\\ +&\qquad (\mathrm{dim},\mathrm{press},\mathrm{bright}),(\mathrm{bright},\mathrm{press},\mathrm{off})\} +\end{aligned} +

+

+\begin{aligned} +\mathcal{G}((\mathrm{off},\mathrm{press},\mathrm{dim})) &= \mathcal X \\ +\mathcal{G}((\mathrm{dim},\mathrm{press},\mathrm{off})) &= \{x \in \mathcal X \mid x>2\}\\ +\mathcal{G}((\mathrm{dim},\mathrm{press},\mathrm{bright})) &= \{x \in \mathcal X \mid x\leq 2\}\\ +\mathcal{G}((\mathrm{bright},\mathrm{press},\mathrm{off})) &= \mathcal X. +\end{aligned} +

+

+r((\mathrm{off},\mathrm{press},\mathrm{dim}),x) = 0, +

+
    +
  • 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\} +

+

+\mathcal{X} = \mathbb R^2, \; \mathcal{X}_0 = \left\{\begin{bmatrix}10\\0\end{bmatrix}\right\} +

+

+\mathcal{I} = \{\mathbb R^2 \mid x_1 > 0 \lor (x_1 = 0 \land x_2 \geq 0)\} +

+

+f(\bm x) = \begin{bmatrix} x_2\\ -g\end{bmatrix} +

+

+\mathcal{E} = \{(q,q)\} +

+

+\mathcal{G} = \{\bm x\in\mathbb R^2 \mid x_1=0 \land x_2 < 0\} +

+

+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 – extending – 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.

+
+

The switch closed

+
+
+
+ +
+
+Figure 12: DC-DC boost converter: the switch closed +
+
+
+

+\begin{bmatrix} +\frac{\mathrm{d}i_\mathrm{L}}{\mathrm{d}t}\\ +\frac{\mathrm{d}v_\mathrm{C}}{\mathrm{d}t} +\end{bmatrix} += +\begin{bmatrix} +-\frac{R_\mathrm{L}}{L}i_\mathrm{L} & 0\\ +0 & -\frac{1}{C(R+R_\mathrm{C})} +\end{bmatrix} +\begin{bmatrix} +i_\mathrm{L}\\ +v_\mathrm{C} +\end{bmatrix} ++ +\begin{bmatrix} +\frac{1}{L}\\ +0 +\end{bmatrix} +v_0 +

+
+
+

Continuous conduction mode (CCM)

+
+
+
+ +
+
+Figure 13: DC-DC boost converter: continuous conduction mode (CCM) +
+
+
+

+\begin{bmatrix} +\frac{\mathrm{d}i_\mathrm{L}}{\mathrm{d}t}\\ +\frac{\mathrm{d}v_\mathrm{C}}{\mathrm{d}t} +\end{bmatrix} += +\begin{bmatrix} +-\frac{R_\mathrm{L}+ \frac{RR_\mathrm{C}}{R+R_\mathrm{C}}}{L} & -\frac{R}{L(R+R_\mathrm{C})}\\ +\frac{R}{C(R+R_\mathrm{C})} & -\frac{1}{C(R+R_\mathrm{C})} +\end{bmatrix} +\begin{bmatrix} +i_\mathrm{L}\\ +v_\mathrm{C} +\end{bmatrix} ++ +\begin{bmatrix} +\frac{1}{L}\\ +0 +\end{bmatrix} +v_0 +

+
+
+

Discontinuous cond. mode (DCM)

+
+
+
+ +
+
+Figure 14: DC-DC boost converter: discontinuous conduction model (DCM) +
+
+
+

+\begin{bmatrix} +\frac{\mathrm{d}i_\mathrm{L}}{\mathrm{d}t}\\ +\frac{\mathrm{d}v_\mathrm{C}}{\mathrm{d}t} +\end{bmatrix} += +\begin{bmatrix} +0 & 0\\ +0 & -\frac{1}{C(R+R_\mathrm{C})} +\end{bmatrix} +\begin{bmatrix} +i_\mathrm{L}\\ +v_\mathrm{C} +\end{bmatrix} ++ +\begin{bmatrix} +0\\ +0 +\end{bmatrix} +v_0 +

+
+

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 +
+
+
+
+
+ + +
+
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/hybrid_automata.html b/hybrid_automata.html index ac67737..53fcc9e 100644 --- a/hybrid_automata.html +++ b/hybrid_automata.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Software

+
+ + + +
+ + + + +
+ + + +
+ + +
+

Matlab

+ +
+
+

Julia

+ +
+
+

Python

+
    +
  • +
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/hybrid_automata_software.html b/hybrid_automata_software.html index 27fc920..2e7ed5c 100644 --- a/hybrid_automata_software.html +++ b/hybrid_automata_software.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Literature

+
+ + + +
+ + + + +
+ + + +
+ + +

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 :-).

+

+ + + + + Back to top

References

+
+
[1]
R. G. Sanfelice, Hybrid Feedback Control. Princeton University Press, 2021. Accessed: Sep. 23, 2020. [Online]. Available: https://press.princeton.edu/books/hardcover/9780691180229/hybrid-feedback-control
+
+
+
[2]
R. Goebel, R. G. Sanfelice, and A. R. Teel, Hybrid Dynamical Systems: Modeling, Stability, and Robustness. Princeton University Press, 2012. Available: https://press.princeton.edu/books/hardcover/9780691153896/hybrid-dynamical-systems
+
+
+
[3]
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.
+
+
+ + +
+ + + + + + \ No newline at end of file diff --git a/hybrid_equations_references.html b/hybrid_equations_references.html index bd40378..479c8d5 100644 --- a/hybrid_equations_references.html +++ b/hybrid_equations_references.html @@ -2,7 +2,7 @@ - + @@ -57,10 +57,10 @@ - + - + - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Course outline

+
+ + + +
+ + + + +
+ + + +
+ + +

The course is structured into 14 topics, each of them corresponding to one lecture. The topics are as follows:

+ +
    +
  1. Solutions of hybrid systems
  2. +
+ +
    +
  1. Mixed-logical dynamical (MLD) description of hybrid systems
  2. +
  3. Model predictive control (MPC) for MLD systems
  4. +
  5. (Formal) verification of hybrid systems
  6. +
+ + + + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/intro_outline.html b/intro_outline.html index b9eb4c4..8ddbd91 100644 --- a/intro_outline.html +++ b/intro_outline.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Max-plus algebra

+
+ + + +
+ + + + +
+ + + +
+ + +

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: -a a+(-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:

+
    +
  1. power,
  2. +
  3. multiplication,
  4. +
  5. addition.
  6. +
+
+
+

(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.

+
+
+Show the code +
using Plots
+x = -5:3
+f(x) = max(2*x+2,x+3,1)
+plot(x,f.(x),label="",thickness_scaling = 2)
+xc = [-2,1]
+yc = f.(xc)
+scatter!(xc,yc,markercolor=[:red,:red],label="",thickness_scaling = 2)
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+

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. +

+
+
+Show the code +
using Plots
+x = -2:0.1:2;
+y = -2:0.1:2;
+f(x,y) = max(0,x,y)
+z = f.(x',y);
+wireframe(x,y,z,legend=false,camera=(5,30))
+xlabel!("x")
+ylabel!("y")
+zlabel!("f(x,y)")
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+

Example 3 (Another 2D polynomial) Consider another 2D (max,+) polynomial +p(x,y) = 0 \oplus x \oplus y \oplus (-1)\otimes x^{\otimes^2} \oplus 1\otimes x\otimes y \oplus (-1)\otimes y^{\otimes^2}. +

+
+
+Show the code +
using Plots
+x = -2:0.1:2;
+y = -2:0.1:2;
+f(x,y) = max(0,x,y,2*x-1,x+y+1,2*y-1)
+z = f.(x',y);
+wireframe(x,y,z,legend=false,camera=(15,30))
+xlabel!("x")
+ylabel!("y")
+zlabel!("p(x,y)")
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+
+
+ +
+
+Piecewise affine (PWA) functions +
+
+
+

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.

+
+

Example 4 (Example) +A = +\begin{bmatrix} +2 & 3 & \varepsilon\\ +1 & \varepsilon & 0\\ +2 & -1 & 3 +\end{bmatrix} +\qquad +A^{\otimes^2} = +\begin{bmatrix} +4 & 5 & 3\\ +3 & 4 & 3\\ +5 & 5 & 6 +\end{bmatrix} +

+
+
+
+
+
+
+ + +G + + +1 + +1 + + + +1->1 + + +2 + + + +2 + +2 + + + +1->2 + + +1 + + + +3 + +3 + + + +1->3 + + +2 + + + +2->1 + + +3 + + + +2->3 + + +-1 + + + +3->2 + + +0 + + + +3->3 + + +3 + + + +
+
+
+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…
  • +
+ +
+
+

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.

+
+

Example 5 (Greatest subsolution) +A = +\begin{bmatrix} +2 & 3 & \varepsilon\\ +1 & \varepsilon & 0\\ +2 & -1 & 3 +\end{bmatrix}, +\qquad +b = +\begin{bmatrix} +1 \\ 2 \\ 3 +\end{bmatrix} +

+

+x = +\begin{bmatrix} +-1\\ -2 \\ 0 +\end{bmatrix} +

+

+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.

+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/max_plus_algebra.html b/max_plus_algebra.html index e6f9a27..0c488dd 100644 --- a/max_plus_algebra.html +++ b/max_plus_algebra.html @@ -2,7 +2,7 @@ - + @@ -71,10 +71,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Software

+
+ + + +
+ + + + +
+ + + +
+ + + + + + + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/max_plus_software.html b/max_plus_software.html index fb6617c..67118f7 100644 --- a/max_plus_software.html +++ b/max_plus_software.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Discrete hybrid automata

+
+ + + +
+ + + + +
+ + + +
+ + +

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

+ +
+

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, +

+

where +h_i(x_c(k), u_c(k)) = \begin{cases}1 & H_i x_c(k) + J_i u_c(k) + K_i \leq 0\\ 0 & \text{otherwise}. \end{cases} +

+
+
+

Finite state machine (FSM)

+

+x_d(k+1) = f_d(x_d(k),u_d(k),\delta_e(k)) +

+
+
+

Mode selector (MS)

+

i(k) \in \{1, 2, \ldots, s\}

+

+i(k) = \mu(x_d(k), u_d(k), \delta_e(k)) +

+
+
+
+

Trajectory of a DHA

+

+\begin{aligned} +\delta_e(k) &= h(x_c(k), u_c(k))\\ +i(k) &= \mu(x_d(k), u_d(k), \delta_e(k))\\ +y_c(k) &= C_{i(k)} x_c(k) + D_{i(k)} u_c(k) + g_{i(k)}\\ +y_d(k) &= g_d(x_d(k), u_d(k), \delta_e(k))\\ +x_c(k+1) &= A_{i(k)} x_c(k) + B_{i(k)} u_c(k) + f_{i(k)}\\ +x_d(k+1) &= f_d(x_d(k),u_d(k),\delta_e(k)) +\end{aligned} +

+
+
+

How to get rid of the IF-THEN conditions in the model?

+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/mld_DHA.html b/mld_DHA.html index ae827ef..5bc1a92 100644 --- a/mld_DHA.html +++ b/mld_DHA.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Logic vs inequalities

+
+ + + +
+ + + + +
+ + + +
+ + +

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\))
  • +
  • Implication: \(X_1 \implies X_2\)
  • +
  • Equivalence: \(X_1 \iff X_2\)
  • +
  • Logical XOR: \(X_1 \oplus X_2\)
  • +
+
+
+
+

Equivalences of logic propositions

+

We will heavily used the following equivalences: \[ +\begin{aligned} +X_1 \implies X_2 \qquad &\equiv \qquad \neg X_2 \implies \neg X_1,\\ +X_1 \iff X_2 \qquad &\equiv \qquad (X_1 \implies X_2) \land (X_2 \implies X_1),\\ +X_1 \land X_2 \qquad &\equiv \qquad \neg (\neg X_1 \lor \neg X_2),\\ +X_1 \implies X_2 \qquad &\equiv \qquad \neg X_1 \lor X_2. +\end{aligned} +\]

+

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\).

+
+ + +
+

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)) +\]

+
+

Example 1 (Example)  

+
+
+
+ +
+
+Figure 1: Example of a FSM +
+
+
+

The state update/transition equation is \[ +\begin{aligned} +x_d(k+1) = +\begin{cases} +\text{Red} & \text{if}\; ([x_d = \text{green}] \land \neg [\delta_3=1]) \lor ([x_d = \text{red}] \land \neg [\delta_3=1])\\ +\text{Green} & \text{if} \; \ldots\\ +\text{Blue} & \text{if} \; \ldots +\end{cases} +\end{aligned} +\]

+

Binary encoding of the discrete states \[ +\text{Red}: x_b = \begin{bmatrix}0\\0 \end{bmatrix}, \; \text{Green}: x_b = \begin{bmatrix}0\\1 \end{bmatrix}, \; \text{Blue}: x_b = \begin{bmatrix}1\\0 \end{bmatrix} +\]

+

Reformulating the state update equations for binary variables \[ +\begin{aligned} +x_{b1} &= (\neg [x_{b1} = 1] \land \neg [x_{b2} = 1] \land \neg [\delta_1=1]) \\ +&\quad (\neg [x_{b1} = 1] \land \neg [x_{b2} = 1] \land [\delta_1=1]) \land [u_{b2}=1]\\ +&\quad (\neg [x_{b1} = 1] \land [x_{b2} = 1] \land \neg [u_{b1}=1] \land [\delta_3=1])\\ +&\quad \lor ([x_{b1} = 1]\land \neg [\delta_2=1])\\ +x_{b2} &= \ldots +\end{aligned} +\]

+

Finally, simplify, convert to CNF.

+
+
+
+

Mixing logical and continuous

+
    +
  • see Indicator variables.
  • +
+
+

Logical implies continuous

+

\[X \implies [f(x)\leq 0]\]

+

\[[\delta = 1] \implies [f(x)\leq 0]\]

+
    +
  • introduce \(M\) \[ +f(x) \leq (1-\delta) M +\]

  • +
  • that is large enough so that when \(\delta=0\), there is no practical restriction on \(f\).

    +
      +
    • Big-M technique.
    • +
  • +
+
+
+

Continuous implies logical

+

\[[f(x)\leq 0] \implies X\]

+

\[[f(x)\leq 0] \implies [\delta = 1]\]

+
    +
  • Equivalently \[\neg [\delta = 1] \implies \neg [f(x)\leq 0],\]

  • +
  • that is, \[[\delta = 0] \implies [f(x) > 0]\]

  • +
  • Introduce \(m\) such that \(f(x)>0\) is enforced when \(\delta=0\) \[ +f(x) > m\delta +\]

  • +
  • but small enough that there is no restriction on \(f\) when \(\delta=1\).

  • +
  • For numerical reasons, modify to nonstrict inequality \[ +f(x) \geq \epsilon + (m-\epsilon)\delta, +\] where \(\epsilon\approx 0\) (for example, machine epsilon).

  • +
+
+
+

Equivalence between logical and continuous

+
    +
  • Combining the previous two implications.
  • +
+

\[ +\begin{aligned} +f(x) &\leq (1-\delta) M,\\ +f(x) &\geq \epsilon + (m-\epsilon)\delta. +\end{aligned} +\]

+
+
+
+

IF-THEN-ELSE rule as an inequality

+
    +
  • If \(X\) +
      +
    • then \(z = a^\top x + b^\top u + f\),
    • +
    • else \(z = 0\).
    • +
  • +
  • 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

+

\[ +\begin{aligned} +h_i(x_c(k), u_c(k)) &\leq M_i (1-\delta_{e,i})\\ +h_i(x_c(k), u_c(k)) &\geq \epsilon + (m_i-\epsilon) \delta_{e,i} +\end{aligned} +\]

+
+
+

Switched affine system

+
    +
  • We want to get rid of the IF-THEN and formulate the switching mechanism into the format of inequalities too.
  • +
+

\[ +x_c(k+1) = \sum_{i=1}^s z_i(k), +\]

+
    +
  • where \[ +z_1(k) = +\begin{cases} +A_1 x_c(k) + B_1 u_c(k) + f_1 & \text{if}\;i(k)=1\\ +0 & \text{otherwise} +\end{cases} +\]
  • +
+

\[\quad \vdots\]

+

\[ +z_s(k) = +\begin{cases} +A_s x_c(k) + B_s u_c(k) + f_s & \text{if}\;i(k)=s\\ +0 & \text{otherwise} +\end{cases} +\]

+
    +
  • For each \(i\in \{1, 2, \ldots, s\}\)
  • +
+

\[ +\begin{aligned} +z_i &\leq M_i\delta_i,\\ +- z_i &\leq -m_i\delta_i,\\ +z_i &\leq a_i^\top x + b_i^\top u + f_i - m_i(1-\delta_i),\\ +-z_i &\leq -(a_i^\top x + b_i^\top u + f_i) + M_i(1-\delta_i). +\end{aligned} +\]

+
+
+

Mixed logical dynamical (MLD) system

+

\[ +\begin{aligned} +x(k+1) &= Ax(k) + B_u u(k) + B_\delta\delta + B_zz(k) + B_0\\ +y(k) &= Cx(k) + D_u u(k) + D_\delta \delta + D_z z + D_0\\ +E_\delta \delta &+ E_z z(k) \leq E_u u(k) + E_x x(k) + E_0 +\end{aligned} +\]

+
+
+

Simple example

+
+
+

HYSDEL language

+
+
+

Piecewise affine systems

+

\[ +\begin{aligned} +x(k+1) &= A_{i(k)}x(k) + B_{i(k)} u(k) + f_{i(k)}\\ +y(k) &= C_{i(k)}x(k) + D_{i(k)} u(k) + g_{i(k)}\\ +& \; H_{i(k)} x(k) + J_{i(k)} u(k) \leq K_{i(k)} +\end{aligned} +\]

+
    +
  • DHA, MLD, PWA are equivalent.
  • +
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/mld_intro.html b/mld_intro.html index 7c3d613..92b5a99 100644 --- a/mld_intro.html +++ b/mld_intro.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + - + - + - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Literature

+
+ + + +
+ + + + +
+ + + +
+ + +

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.

+ + + + + Back to top

References

+
+
[1]
F. Borrelli, A. Bemporad, and M. Morari, Predictive Control for Linear and Hybrid Systems. Cambridge, New York: Cambridge University Press, 2017. Available: http://cse.lab.imtlucca.it/~bemporad/publications/papers/BBMbook.pdf
+
+
+ + +
+ + + + + + \ No newline at end of file diff --git a/mpc_mld_references.html b/mpc_mld_references.html index 8a8dd03..f0a17cc 100644 --- a/mpc_mld_references.html +++ b/mpc_mld_references.html @@ -2,7 +2,7 @@ - + @@ -57,10 +57,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Software

+
+ + + +
+ + + + +
+ + + +
+ + +

Essentially the same as in the previous chapter.

+ + + + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/mpc_mld_software.html b/mpc_mld_software.html index 167d289..b53a907 100644 --- a/mpc_mld_software.html +++ b/mpc_mld_software.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Petri nets

+
+ + + +
+ + + + +
+ + + +
+ + +

In this chapter we introduce another formalism for modelling discrete event systems (DES) – Petri nets. Petri nets offer an alternative perspective on discrete even systems compared to automata. And it is good to have alternatives, isn’t it? For some purposes, one framework can be more appropriate than the other.

+

Furthermore, the ideas behind Petri nets even made it into international standards. Either directly or through the derived GRAFCET language, which in turn served as the basis for the Sequential Function Chart (SFC) language for PLC programming. See the references.

+

Last but not least, an elegant algebraic framework based on the so-called (max,+) algebra has been developed for a subset of Petri nets (so-called event graphs) and it would be a shame not to mention it in our course (in the next chapter).

+
+

Definition of a Petri net

+

Similarly as in the case of automata, a Petri net (PN) can be defined as a tuple of sets and functions: \boxed{PN = \{\mathcal{P}, \mathcal{T}, \mathcal{A}, w\},} where

+
    +
  • \mathcal{P} = \{p_1, \dots, p_n\} is a finite set of places,
  • +
  • \mathcal{T} = \{t_1, \dots, t_m\} is a finite set of transitions,
  • +
  • \mathcal{A} \subseteq (\mathcal{P} \times \mathcal{T}) \cup (\mathcal{T} \times \mathcal{P}) is a finite set of arcs, and these arcs are directed and since there are two types of nodes, there are also two types of arcs: +
      +
    • (p_i, t_j) \in \mathcal{A} is from place p_i to transition t_j,
    • +
    • (t_j, p_i) \in \mathcal{A} is from transition t_j to place p_i,
    • +
  • +
  • w : \mathcal{A} \to \mathbb{N} is a weight function.
  • +
+

Similarly as in the case of automata, Petri nets can be visualized using graphs. But this time, we need to invoke the concept of a weighted bipartite graph. That is, a graph with two types of nodes:

+
    +
  • places = circles,
  • +
  • transitions = bars.
  • +
+

The nodes of different kinds are connected by arcs (arrowed curves). A integer weights are associated with arcs. Alternatively, for a smaller weight (2, 3, 4), the weithg can be graphically encoded by drawing multiple arcs.

+
+

Example 1 (Simple Petri net) We consider just two places, that is, \mathcal{P} = \{p_1, p_2\}, and one transition, that is, \mathcal{T} = \{t\}. The set of arcs is \mathcal{A} = \{\underbrace{(p_1, t)}_{a_1}, \underbrace{(t, p_2)}_{a_2}\}, and the associated weights are w(a_1) = w((p_1, t)) = 2 and w(a_2) = w((t, p_2)) = 1. The Petri net is depicted in Fig 1.

+
+
+
+ +
+
+Figure 1: Example of a simple Petri net +
+
+
+
+
+

Additional definitions

+
    +
  • \mathcal{I}(t_j) … a set of input places of the transition t_j,
  • +
  • \mathcal{O}(t_j) … a set of output places of the transition t_j.
  • +
+
+

Example 2 (More complex Petri net)  

+
    +
  • \mathcal{P} = \{p_1, p_2, p_3, p_4\},
  • +
  • \mathcal{T} = \{t_1, t_2, t_3, t_4, t_5\},
  • +
  • \mathcal{A} = \{(p_1, t_1), (t_1, p_1), (p_1, t_2),\ldots\},
  • +
  • w((p_1, t_1)) = 2, \; w((t_1, p_1)) = 1, \; \ldots
  • +
+
+
+
+ +
+
+Figure 2: Example of a more complex Petri net +
+
+
+
+
+
+

Marking and marked Petri nets

+

An important concept that we must introduce now is that of marking. It is a function that assigns an integer to each place x: \mathcal{P} \rightarrow \mathbb{N}.

+

The vector composed of the values of the marking function for all places \bm x = \begin{bmatrix}x(p_1)\\ x(p_2)\\ \vdots \\ x(p_n) \end{bmatrix} can be viewed as the state vector (although the Petri nets community perhaps would not use this terminology and stick to just marking).

+

Marked Petri net is then a Petri net augmented with the marking

+

MPN = \{\mathcal{P}, \mathcal{T}, \mathcal{A}, w,x\}.

+
+

Visualization of marked Petri net using tokens

+

Marked Petri net can also be visualized by placing tokens (dots) into the places. The number of tokens in a place corresponds to the value of the marking function for that place.

+
+

Example 3 (Marked Petri net) Consider the Petri net from Example 1. The marking function is x(p_1) = 2 and x(p_2) = 1, which assembled into a vector gives \bm x = \begin{bmatrix}1\\ 0 \end{bmatrix}. The marked Petri net is depicted in Fig 3.

+
+
+
+ +
+
+Figure 3: Example of a marked Petri net +
+
+
+

For another marking, namely \bm x = \begin{bmatrix}2\\ 1 \end{bmatrix}, the marked Petri net is depicted in Fig 4.

+
+
+
+ +
+
+Figure 4: Example of a marked Petri net with different marking +
+
+
+
+
+
+
+

Enabling and firing of a transition

+

Finally, here comes the enabling (pun intended) component of the definition of a Petri net – enabled transition. A transition t_j does not just happen – we say fire – whenever it wants, it can only happen (fire) if it is enabled, and the marking is used to determined if it is enabled. Namely, the transition is enabled if the value of the marking function for each input place is greater than or equal to the weight of the arc from that place to the transition. That is, the transition t_j is enabled if +x(p_i) \geq w(p_i,t_j)\quad \forall p_i \in \mathcal{I}(t_j). +

+
+
+
+ +
+
+Can but does not have to +
+
+
+

The enabled transition can fire, but it doesn’t have to. We will exploit this in timed PN.

+
+
+
+

Example 4 (Enabled transition) See the PN in Example 3: in the first marked PN the transition cannot fire, in the second it can.

+
+
+
+
+

State transition function

+

We now have a Petri net as a conceptual model with a graphical representation. But in order to use it for some quantitative analysis, it is useful to turn it into some computational form. Preferrably a familiar one. This is done by defining a state transition function. For a Petri net with n places, the state transition function is +f: \mathbb N^n \times \mathcal{T} \rightarrow \mathbb N^n, + which reads that the state transition fuction assignes a new marking (state) to the Petri net after a transition is fired at some given marking (state).

+

The function is only defined for a transition t_j iff the transition is enabled.

+

If the transition t_j is enabled and fired, the state evolves as +\bm x^+ = f(\bm x, t_j), + where the individual components of \bm x evolve according to \boxed{ + x^+(p_i) = x(p_i) - w(p_i,t_j) + w(t_j,p_i), \; i = 1,\ldots,n.} +

+

This has a visual interpretation – a fired transition moves tokens from the input to the output places.

+
+

Example 5 (Moving tokens around) Consider the PN with the initial marking (state) \bm x_0 = \begin{bmatrix}2\\ 0\\ 0\\ 1 \end{bmatrix} (at discrete time 0), and the transition t_1 enabled

+
+
+

+
Example PN in the initial state at time 0, the transition t_1 enabled, but not yet fired
+
+
+
+
+
+ +
+
+Conflict in notation. Again. Sorry. +
+
+
+

We admit the notation here is confusing, because we use the lower index 0 in \bm x_0 to denote the discrete time, while the lower index 1 in t_1 to denote the transition and the lower indices 1 and 2 in p_1 and p_2 just number the transitions and places, respectively. We could have chosen something like \bm x(0) or \bm x[0], but we dare to hope that the context will make it clear.

+
+
+

Now we assume that t_1 is fired

+
+
+
+ +
+
+Figure 5: Transition t_1 in the example PN fired at time 1 +
+
+
+

The state vector changes to \bm x_1 = [1, 1, 1, 1]^\top, the discrete time is 1 now.

+

As a result of this transition, note that t_1, t_2, t_3 are now enabled.

+
+
+
+ +
+
+Number of tokens need not be preserved +
+
+
+

In the example we can see for the first time, that the number of tokens need not be preserved.

+
+
+

Now fire the t_2 transition

+
+
+

+
Transition t_2 in the example PN fired at time 2
+
+
+

The state vector changes to \bm x_2 = [1, 1, 0, 2]^\top, the discrete time is 2 now.

+

Good, we can see the idea. But now we go back to time 1 (as in Fig 5) to explore the alternative evolution. With the state vector \bm x_1 = [1, 1, 1, 1]^\top and the transitions t_1, t_2, t_3 enabled, we fire t_3 this time.

+
+
+

+
Transition t_3 in the example PN fired at time 1
+
+
+

The state changes to \bm x_2 = [0, 1, 0, 0]^\top, the discrete time is 2. Apparently the PN evolved into at a different state. The lesson learnt with this example is that the order of firing of enabled transitions matters.

+
+
+
+
+ +
+
+The order in which the enabled transitions are fired does matter +
+
+
+

The dependence of the state evolution upon the order of firing the transitions is not surprising. Wwe have already encountered it in automata when the active event set for a given state contains more then a single element.

+
+
+
+
+

Reachability

+

We have started talking about states and state transitions in Petri nets, which are all concepts that we are familiar with from dynamical systems. Another such concept is reachability. We explain it through an example.

+
+

Example 6 (Not all states are reachable)  

+
+
+

+
Example of an unreachability in a Petri net
+
+
+

The Petri net is initial in the state [2,1]^\top. The only reachable state is [0,2]^\top.

+

By the way, note that the weight of the arc from the place p_1 to the transition t is 2, so both tokens are removed from the place p_1 when the transition t fires. But then the arc to the place p_2 has weight 1, so only one token is added to the place p_2. The other token is “lost”.

+
+
+

Reachability tree and graph

+

Here we introduce two tools for analysis of reachability of a Petri net.

+
+

Example 7 Consider the following example of a Petri net.

+
+
+

+
Example of a Petri net for reachability analysis
+
+
+

In Fig 6 we draw a reachability tree for this Petri net.

+
+
+
+ +
+
+Figure 6: Reachability tree for an example Petri net +
+
+
+

In Fig 7 we draw a reachability graph for this Petri net.

+
+
+
+ +
+
+Figure 7: Reachability graph for an example Petri net +
+
+
+
+
+
+
+

Number of tokens need not be preserved

+

We have already commented on this before, but we emphasize it here. Indeed, it can be that +\sum_{p_i\in\mathcal{O}(t_j)}w(t_j,p_i) < \sum_{p_i\in\mathcal{I}(t_j)} w(p_i,t_j) +

+

or

+

+\sum_{p_i\in\mathcal{O}(t_j)}w(t_j,p_i) > \sum_{p_i\in\mathcal{I}(t_j)} w(p_i,t_j) +

+

With this reminder, we can now hightlight several patters that can be observed in Petri nets.

+
+
+

AND-convergence, AND-divergence

+

+
+
+

OR-convergence and OR-divergence

+

+
+
+

Nondeterminism in a PN

+

In the four patters just enumerated, we have seen that the last one – the OR-divergence – is not deterministic. Indeed, consider the following example.

+
+

Example 8  

+
+
+

+
Nondeterminism in a Petri net
+
+
+
+

In other words, we can incorporate a nondeterminism in a model.

+
+
+
+ +
+
+Nondeterminism in automata +
+
+
+

Recall that something similar can be encountered in automata, if the active event set for a given state contains more than one element (event,transition).

+
+
+
+
+

Subclasses of Petri nets

+

We can identify two subclasses of Petri nets:

+
    +
  • event graphs,
  • +
  • state machines.
  • +
+
+

Event graph

+
    +
  • Each place has just one input and one output transition (all ws equal to 1).
  • +
  • No OR-convergence, no OR-divergence.
  • +
  • Also known as Decision-free PN.
  • +
  • It can model synchronization.
  • +
+
+

Example 9 (Event graph)  

+
+
+

+
Example of an event graph
+
+
+
+
+
+

State machine

+
    +
  • Each transition has just one input and one output place.
  • +
  • No AND-convergence, no AND-divergence.
  • +
  • Does not model synchronization.
  • +
  • It can model race conditions.
  • +
  • With no source (input) and sink (output) transitions, the number of tokens is preserved.
  • +
+
+

Example 10 (State machine)  

+
+
+

+
Example of a state machine
+
+
+
+
+
+

Incidence matrix

+

We consider a Petri net with n places and m transitions. The incidence matrix is defined as +\bm A \in \mathbb{Z}^{n\times m}, + where +a_{ij} = w(t_j,p_i) - w(p_i,t_j). +

+
+
+
+ +
+
+Transpose +
+
+
+

Some define the incidence matrix as the transpose of our definition.

+
+
+
+
+

State equation for a Petri net

+

With the incidence matrix defined above, the state equation for a Petri net can be written as +\bm x^+ = \bm x + \bm A \bm u, + where \bm u is a firing vector for the enabled j-th transition +\bm u = \bm e_j = \begin{bmatrix}0 \\ \vdots \\ 0 \\ 1\\ 0\\ \vdots\\ 0\end{bmatrix} + with the 1 at the j-th position.

+
+
+
+ +
+
+State vector as a column rather then a row +
+
+
+

Note that in [1] they define everything in terms of the transposed quantities, but we prefer sticking to the notion of a state vector as a column.

+
+
+
+

Example 11 (State equation for a Petri net) Consider the Petri net from Example 5, which we show again below in Fig 8.

+
+
+
+ +
+
+Figure 8: Example of a Petri net +
+
+
+

The initial state is given by the vector +\bm x_0 += +\begin{bmatrix} +2\\ 0\\ 0\\ 1 +\end{bmatrix} +

+

The incidence matrix is +\bm A = \begin{bmatrix} +-1 & 0 & -1\\ +1 & 0 & 0\\ +1 & -1 & -1\\ +0 & 1 & -1 +\end{bmatrix} +

+

And the state vector evolves according to +\begin{aligned} +\bm x_1 &= \bm x_0 + \bm A \bm u_1\\ +\bm x_2 &= \bm x_1 + \bm A \bm u_2\\ +\vdots & +\end{aligned} +

+
+
+
+
+ +
+
+Caution +
+
+
+

We repeat once again just to make sure: the lower index corresponds to the discrete time.

+
+
+
+
+
+

Queueing systems modelled by PN

+

Although Petri nets can be used to model a vast variety of systems, below we single out one particular class of systems that can be modelled by Petri nets – queueing systems. The general symbol is shown in Fig 9.

+
+
+
+ +
+
+Figure 9: Queing systems and its components and transitions/events +
+
+
+

We can associate the transitions with the events in the queing system:

+
    +
  • a is a spontaneous transition (no input places).
  • +
  • s needs a customer in the queue and the server being idle.
  • +
  • c needs the server being busy.
  • +
+

We can now start drawing the Petri net by drawin the bars corresponding to the transitions. Then in between every two bars, we draw a circle for a place. The places can be associated with three bold-face letters above, namely:

+

\quad \mathcal{P} = \{Q, I, B\}, that is, queue, idle, busy.

+
+
+

+
Petri net corresponding to the queing system
+
+
+
+
+
+ +
+
+The input transition adds tokens to the system +
+
+
+

The transition a is an input transition – the tokens are added to the system through this transition.

+
+
+

Note how we consider the token in the I place. This is not only to express that the server is initially idle, ready to serve as soon as a customer arrives to the queue, it also ensures that no serving of a new customer can start before the serving of the current customer is completed.

+

The initial state: [0,1,0]^\top. Consider now a particular trace (of transitions/events) \{a,s,a,a,c,s,a\}. Verify that this leads to the final state [2,0,1]^\top.

+
+

Some more extensions

+

We can keep adding features to the model of a queing system. In particular,

+
    +
  • the arrival transition always enabled,
  • +
  • the server can break down, and then be repaired,
  • +
  • completing the service \neq customer departure.
  • +
+

These are incorporated into the Petri net in Fig 10.

+
+
+
+ +
+
+Figure 10: Extended model of a queueing system +
+
+
+

In the Petri net, d is an output transition – the tokens are removed from the system.

+
+

Example 12 (Beverage vending machine) Below we show a Petri net for a beverage vending machine. While building it, we find it useful to identify the events/transitions that can happen in the system.

+
+
+

+
Petri net for a beverage vending machine
+
+
+
+
+
+
+

Some extensions of basic Petri nets

+
    +
  • Coloured Petri nets (CPN): tokens can by of several types (colours), and the transitions can be enabled only if the tokens have the right colours.
  • +
  • +
+

We do not cover these extensions in our course. But there is one particular extension that we do want to cover, and this amounts to introducing time into Petri nets, leading to timed Petri nets, which we will discuss in the next chapter.

+ + + +
+ + Back to top

References

+
+
[1]
C. G. Cassandras and S. Lafortune, Introduction to Discrete Event Systems, 3rd ed. Cham: Springer, 2021. Available: https://doi.org/10.1007/978-3-030-72274-6
+
+
+ + +
+ + + + + + \ No newline at end of file diff --git a/petri_nets.html b/petri_nets.html index a7975e5..23641e8 100644 --- a/petri_nets.html +++ b/petri_nets.html @@ -2,7 +2,7 @@ - + @@ -57,10 +57,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Literature

+
+ + + +
+ + + + +
+ + + +
+ + +

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.

+ + + +
+ + Back to top

References

+
+
[1]
C. G. Cassandras and S. Lafortune, Introduction to Discrete Event Systems, 3rd ed. Cham: Springer, 2021. Available: https://doi.org/10.1007/978-3-030-72274-6
+
+
+
[2]
F. Baccelli, G. Cohen, G. J. Olsder, and J.-P. Quadrat, Synchronization and linearity: An algebra for discrete event systems, Web edition. Chichester: Wiley, 2001. Available: https://www.rocq.inria.fr/metalau/cohen/documents/BCOQ-book.pdf
+
+
+
[3]
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
+
+
+
[10]
IEC 60848:2013 GRAFCET specification language for sequential function charts.” IEC, Feb. 2013. Available: https://webstore.iec.ch/publication/3684
+
+
+
[11]
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.
+
+
+
[12]
IEC 61131-3 Programmable controllers - Part 3: Programming languages.” IEC, Feb. 2013. Accessed: Jan. 08, 2023. [Online]. Available: https://webstore.iec.ch/publication/4552
+
+
+ + +
+ + + + + + \ No newline at end of file diff --git a/petri_nets_references.html b/petri_nets_references.html index eeefbdb..d4394cd 100644 --- a/petri_nets_references.html +++ b/petri_nets_references.html @@ -2,7 +2,7 @@ - + @@ -57,10 +57,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Timed Petri nets

+
+ + + +
+ + + + +
+ + + +
+ + +

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

+

\boxed{TPN = \{\mathcal{P}, \mathcal{T}, \mathcal{A}, w, x, \mathcal{V}\}}.

+
+

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} +

+

Combining gives the desired update equations +\begin{aligned} +a_k &= a_{k-1} + v_{a,k},\quad k=1,2,\ldots,\quad a_0 = 0\\ +s_k &= \max\{a_{k},c_{k-1}\}\\ +c_k &= s_{k} + v_{c,k}\\ +&= \max\{a_{k},c_{k-1}\} + v_{c,k},\quad k=1,\ldots, \quad c_0=0 +\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.

+ + +
+
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/petri_nets_timed.html b/petri_nets_timed.html index 4ce3796..a386fe2 100644 --- a/petri_nets_timed.html +++ b/petri_nets_timed.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Solution concepts

+
+ + + +
+ + + + +
+ + + +
+ + +

Now that we know how to model hybrid systems, we need to define what we mean by a solution to a hybrid system. The definitions are not as straightforward as in the continuous-time or discrete-time case and mastering them is not only of theoretical value.

+
+

Hybrid time and hybrid time domain

+

Even before we start discussing the concepts of a solution, we need to discuss the concept of time in hybrid systems. Of course, hybrid systems live in the same world as we do, and therefore they evolve in the same physical time, but it turns out that we can come up with an artificial concept of hybrid time that makes modelling and analysis of hybrid systems convenient.

+

Recall that in continuous-time systems, the continuous time t\in\mathbb R_{\geq 0}, and in discrete-time systems, the discrete “time” k\in\mathbb N. We put the latter in quotation marks since k is not really the time but rather it should be read as the kth transition of the system. Now, the idea is to combine these two concepts of time into one, and we call it the hybrid time (t,j), \; t\in \mathbb R_{\geq 0},\, j\in \mathbb N.

+

If you think that it is redundant, note that since hybrid systems can exhibit discrete-event system behaviour, a transition from one discrete state to another can happen instantaneously. In fact, several such transitions can take no time at all. It sounds weird, but that is what the mathematical model allows. That is why determining t need not be enought and we also need to specify j.

+

The set of all hybrid times for a given hybrid system is called hybrid time domain +E \subset [0,T] \times \{0,1,2,\ldots, J\}, + where T and J can be finite or \infty.

+

In particular, +E = \bigcup_{j=0}^J \left([t_j,t_{j+1}] \times \{j\}\right) +\tag{1}

+

where 0=t_0 < t_1 < \ldots < t_J = T.

+

The meaning of Eq. 1 can be best explained using Fig. 1 below.

+
+
+
+ +
+
+Figure 1: Example of a hybrid time domain +
+
+
+

Note that if two hybrid times are from the same hybrid domain, we can decide if (t,j) \leq (t',j'). In other words, the set of hybrid times is totally ordered.

+
+
+

Hybrid arc

+

Hybrid arc is just a terminology used in the literature for hybrid state trajectory. It is a function that assigns a state vector x to a given hybrid time (t,j) +x: E \rightarrow \mathbb R^n. +

+

For each j the function t \mapsto x(t,j) is absolutely continuous on the interval I^j = \{t \mid (t,j) \in E\}.

+
+
+
+ +
+
+Inconsitent notation +
+
+
+

We admit here that we are not going to be 100% consistent in the usage of the notation x(t,j) in the rest of our course. Oftentimes use x(t) even within hybrid systems when we do not need to index the jumps.

+
+
+

It is perhaps clear now, that hybrid time domain can only be determined once the solution (the arc, the trajectory) is known. This is in sharp contrast with the continuous-time or discrete-time system – we can formulate the problem of finding solution to \dot x(t) = 3x(t), \, x(0) = 1 on the interval [0,2], where the interval was set even before we know how the solution looks like.

+
+
+

Solutions of autonomous (no-input) systems

+

Finally we can formalize the concept of a solution. A hybrid arc x(\cdot,\cdot) is a solution to the hybrid equations given by the common quadruple \{\mathcal{C},\mathcal{D},f,g\} (or \{\mathcal{C},\mathcal{D},\mathcal{F},\mathcal{G}\} for inclusions), if

+
    +
  • the initial state x(0,0) \in \overline{\mathcal{C}} \cup \mathcal{D}, and
  • +
  • for all j such that I^j = \{t\mid (t,j)\in E\} has a nonempty interior \operatorname{int}I^j +
      +
    • x(t,j) \in \mathcal C \; \forall t\in \operatorname{int}I^j,
    • +
    • \dot x(t,j) = f(x(t,j)) \; \text{for almost all}\; t\in I^j, and
    • +
  • +
  • for all (t,j)\in E such a (t,j+1)\in E +
      +
    • x(t,j) \in \mathcal{D}, and
    • +
    • x(t,j+1) = g(x(t,j)).
    • +
  • +
+

Make the modifications for the \{\mathcal{C},\mathcal{D},\mathcal{F},\mathcal{G}\} version by yourself.

+
+

Example 1 (Solution) En example of a solution is in Fig. 2. Follow the solution with your finger and make sure you understand what and why is happing. In particular, in the overlapping region, the solution is not unique. While it can continue flowing, it can also jump.

+
+
+
+ +
+
+Figure 2: Example of a solution +
+
+
+
+
+
+

Hybrid input

+

Similarly as we considered the state as a function of the hybrid time, we can consider the input as a function of the hybrid time. With its own hybrid domain E_\mathrm{u}, the input is +u: E_\mathrm{u} \rightarrow \mathbb R^m. +

+

For each j the function t \mapsto u(t,j) must be… well-behaved… For example, piecewise continuous on the interval I^j = \{t \mid (t,j) \in E_\mathrm{u}\}.

+
+
+

Solutions of systems with inputs

+

We assume that hybrid time domains for the arcs and inputs are the same. A solution must satisfy the same conditions as in the case of autonomous systems, but with the input taken into account. For completeness we state the conditions here:

+
    +
  • The initial state-control pair (x(0,0),u(0,0)) \in \overline{\mathcal{C}} \cup \mathcal{D}, and
  • +
  • for all j such that I^j = \{t\mid (t,j)\in E\} has a nonempty interior \operatorname{int}I^j +
      +
    • (x(t,j),u(t,j)) \in \mathcal C \; \forall t\in \operatorname{int}I^j,
    • +
    • \dot x(t,j) = f(x(t,j),u(t,j)) \; \text{for almost all}\; t\in I^j, and
    • +
  • +
  • for all (t,j)\in E such a (t,j+1)\in E +
      +
    • (x(t,j),u(t,j)) \in \mathcal{D}, and
    • +
    • x(t,j+1) = g(x(t,j),u(t,j)).
    • +
  • +
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/solution_concepts.html b/solution_concepts.html index 1fad19f..833021e 100644 --- a/solution_concepts.html +++ b/solution_concepts.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Types of solutions

+
+ + + +
+ + + + +
+ + + +
+ + +

Now that we know, what a hybrid arc (trajectory) needs to satisfy to be a solution of a hybrid system, we can classify the solutions into several types. And we base this classification on their hybrid time domain E:

+
+
Trivial
+
+just one point. +
+
Nontrivial
+
+at least two points; +
+
Complete
+
+if the domain is unbounded; +
+
Bounded, compact
+
+if the domain is bounded, compact (well, it is perhaps a bit awkward to call a solution bounded just based on boundednes of its time domain as most people would interpret the boundedness of a solution with regard to the values of the solution); +
+
Discrete
+
+if nontrivial and E\subset \{0\} \times \mathbb N; +
+
Continuous
+
+if nontrivial and E\subset \mathbb R_{\geq 0} \times \{0\}; +
+
Eventually discrete
+
+if T = \sup_E t < \infty and E \cap (\{T\}\times \mathbb N) contains at least two points; +
+
Eventually continuous
+
+if J = \sup_E j < \infty and E \cap (\mathbb R_{\geq 0} \times \{J\}) contains at least two points; +
+
Zeno
+
+if complete and \sup_E t < \infty; +
+
Maximal
+
+It cannot be extended. A solution x(t,j) defined on the hybrid time domain E is maximal, if on an extended hybrid time domain E^\mathrm{ext} such that E\subset E^\mathrm{ext}, there is no solution x^\mathrm{ext}(t,j) that coincides with x on E. Some literature uses the “linguistic” terminology that a maximal solution is not a prefix to any other solution. Complete solutions are maximal. But not vice versa. +
+
+
+
+
+ +
+
+Tip +
+
+
+

It is certainlty helpful to sketch the times domains for the individual classes of solutions.

+
+
+
+

Examples of types of solutions

+
+

Example 1 (Example of a (non-)maximal solution) +\dot x = 1, \; x(0) = 1 +

+

+(t,j) \in [0,1] \times \{0\} +

+

Now extend the time domain to +(t,j) \in [0,2] \times \{0\}. +

+

Can we extend the solution?

+
+
+

Example 2 (Maximal but not complete continuous solution) Finite escape time

+

\dot x = x^2, \; x(0) = 1,

+

x(t) = 1/(1-t)

+
+
+

Example 3 (Discontinuous right hand side) \dot x = \begin{cases}-1 & x>0\\ 1 & x\leq 0\end{cases}, \quad x(0) = -1 (unless the concept of Filippov solution is invoked).

+
+
+

Example 4 (Zeno solution of the bouncing ball) Starting on the ground with some initial upward velocity +h(t) = \underbrace{h(0)}_0 + v(0)t - \frac{1}{2}gt^2, \quad v(0)=1 +

+

What time will it hit the ground again? +0 = t - \frac{1}{2}gt^2 = t(1-\frac{1}{2}gt) +

+

t_1=\frac{2}{g}

+

Simplify (scale) the computations just to get the qualitative picture: set g=2, which gives t_1 = 1.

+

t_1=1:

+

v(t_1^+) = \gamma v(t_1) = \gamma v(0) = \gamma

+

The next hit will be at t_1 + \tau_1 h(t_1 + \tau_1) = 0 = \gamma \tau_1 - \tau_1^2 = \tau_1(\gamma - \tau_1) \tau_1 = \gamma

+

t_2 = t_1+\tau_1 = 1 + \gamma:\quad \ldots

+

t_k = 1 + \gamma + \gamma^2 + \ldots + \gamma^k:\quad \ldots \boxed{\lim_{k\rightarrow \infty} t_k = \frac{1}{1-\gamma} < \infty}

+

Infinite number of jumps in a finite time!

+
+
+

Example 5 (Water tank)  

+
+
+

+
Switching between two water tanks
+
+
+

+\max \{Q_\mathrm{out,2}, Q_\mathrm{out,2}\} \leq Q_\mathrm{in} \leq Q_\mathrm{out,2} + Q_\mathrm{out,2} +

+
+
+

+
Hybrid automaton for switching between two water tanks
+
+
+
+
+

Example 6 ((Non)blocking and (non)determinism in hybrid systemtems)  

+
+
+

+
Example of an automaton exhibitting (non)blocking and (non)determinism
+
+
+
    +
  • x(0) = -3
  • +
  • x(0) = -2
  • +
  • x(0) = -1
  • +
  • x(0) = 0
  • +
+
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/solution_types.html b/solution_types.html index 0f19046..e72fbfd 100644 --- a/solution_types.html +++ b/solution_types.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

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

+

\dot{\bm x} = \mathbf f(\bm x_\mathrm{eq}) = \mathbf 0.

+

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.

+
+
+Show the code +
f(x) = [(x[1]^2*(x[2]-x[1])+x[2]^5)/((x[1]^2+x[2]^2)*(1+(x[1]^2+x[2]^2)^2)); 
+        (x[2]^2*(x[2]-2x[1]))/((x[1]^2+x[2]^2)*(1+(x[1]^2+x[2]^2)^2))]
+
+using CairoMakie
+fig = Figure(; size = (800, 800),fontsize=20)
+ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
+streamplot!(ax,x->Point2f(f(x)), -1.5..1.5, -1.5..1.5, colormap = :magma)
+fig
+
+
+
+
+

+
+
+
+
+
+
+

Asymptotic stability

+

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.

+
+
+
+

Exponential stability

+

k_1 \|x\|^p \leq V(x) \leq k_2 \|x\|^p,

+

\left(\nabla V(x)\right)^\top f(x) \leq -k_3 \|x\|^p.

+
+
+

Exponential stability with quadratic Lyapunov function

+

+V(x) = x^\top P x +

+

\lambda_{\min} (P) \|x\|^2 \leq V(x) \leq \lambda_{\max} (P) \|x\|^2

+
+
+

Converse theorems

+
    +
  • for (G)UAS,
  • +
  • for Lyapunov stability only time-varying Lyapunov function guaranteed.
  • +
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/stability_recap.html b/stability_recap.html index b701933..8675a8f 100644 --- a/stability_recap.html +++ b/stability_recap.html @@ -2,7 +2,7 @@ - + @@ -71,10 +71,10 @@ - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Stability via common Lyapunov function

+
+ + + +
+ + + + +
+ + + +
+ + +

Having just recalled the stability analysis based on searching for a Lyapunov function, we now extend the analysis to hybrid systems, in particular, hybrid automata.

+
+

Hybrid system stability analysis via Lyapunov function

+

In contrast to continuous systems where Lyapunov functin should decrease along the state trajectory to certify asymptotic stability (or at least should not increase to certify Lyapunov stability), in hybrid systems the requirement can be relaxed a bit.

+

A Lyapunov function should still decrease along the continuous state trajectory while the system resides in a given discrete state (mode), but during the transition between the modes the function can be discontinuous, and can even increase. But then upon return to the same discrete state, the function should have lower value than it had last time it entered the state to certify asymptotic stability (or at least not larger for Lyapunov stability), see Fig. 1.

+
+
+
+ +
+
+Figure 1: Example of an evolution of a Lyapunov function of a hybrid automaton in time +
+
+
+

Formally, the function V(q,x) has two variables q and x, it is smooth in x, and
+ +V(q,0) = 0, \quad V(q,x) > 0 \; \text{for all nonzero} \; x \; \text{and for all nonzero} \; q, +

+

+\left(\nabla_x V(q,x)\right)^\top f(q,x) < 0 \; \text{for all nonzero} \; x \; \text{and for all nonzero} \; q, +

+

and the discontinuitieds at transitions must satisfy the conditions sketched in Fig. 1. If Lyapunov stability is enough, the strict inequality can be relaxed to nonstrict.

+
+

Stricter condition on Lyapunov function

+

Verifying the properties just described is not easy. Perhaps the only way is to simulate the system, evaluate the function along the trajectory and check if the conditions are satisfied, which is hardly useful. A stricter condition is in Fig. 2. Here we require that during transitions to another mode, the fuction should not increase. Discontinuous reductions in value are allowed, but enforcing continuity introduced yet another simplification that can be plausible for analysis.

+
+
+
+ +
+
+Figure 2: Example of an evolution of a restricted Lyapunov function of a hybrid automaton in time +
+
+
+
+
+
+

Further restricted set of candidate functions: common Lyapunov function (CLF)

+

In the above, we have considered a Lyapunov function V(q,x) that is mode-dependent. But what if we could find a single Lyapunov function V(x) that is common for all modes q? This would be a great simplification.

+

This, however, implies that at a given state x, arbitrary transition to another mode q is possible. We add and adjective uniform the stabilty (either Lyapunov or asymptotic) to emphasize that the function is common for all modes.

+

In terminology of switched systems, we say that arbitrary switching is allowed. And staying in the domain of switched systems, the analysis can be interpreted within the framework of differential inclusions +\dot x \in \{f_1(x), f_2(x),\ldots,f_m(x)\}. +

+
+
+

(Global) uniform asymptotic stability

+

Having agreed that we are now searching for a single function V(x), the function must satisfy \boxed{\kappa_1(\|x\|) \leq V(x) \leq \kappa_2(\|x\|),} where \kappa_1(\cdot), \kappa_2(\cdot) are class \mathcal{K} comparison functions, and \kappa_1(\cdot)\in\mathcal{K}_\infty if global asymptotic stability is needed, and +\boxed{\left(\nabla V(x)\right)^\top f_q(x) \leq -\rho(\|x\|),\quad q\in\mathcal{Q},} + where \rho(\cdot) is a positive definite continuous function, zero at the origin. If all these are satisfied, the system is (globally) uniformly asymptotically stable (GUAS).

+

Similarly as we did in continuous systems, we must ask if stability implies existence of a common Lyapunov function. An affirmative answer comes in the form of a converse theorem for global uniform asymptotic stability (GUAS).

+

This is great, a CLF exists for a GUAS system, but how do we find it? We must restrict the set of candidate functions and then search within the set. Obviously, if we fail to find a function in the set, we must extend the set and search again…

+
+
+

Common quadratic Lyapunov function (CQLF)

+

An immediate restriction of a set of Lyapunov functions is to quadratic functions +V(x) = x^\top P x, + where P=P^\top \succ 0.

+

This restriction is also quite natural because for linear systems, it is known that we do not have to consider anything more complicated than quadratic functions. This does not hold in general for nonlinear and hybrid systems. But it is a good start. If we succeed in finding a quadratic Lyapunov function, we can be sure that the system is stable.

+

Here we start by considering a hybrid automaton for which at each mode the dynamics is linear. We consider r continuous-time LTI systems parameterized by the system matrices A_i for i=1,\ldots, r as +\dot x = A_i x(t). +

+

Time derivatives of V(x) along the trajectory of the i-th system + \nabla V(x)^\top \left.\frac{\mathrm d x}{\mathrm d t}\right|_{\dot x = A_i x} = x^\top(A_i^\top P + PA_i)x, +

+

which, upon introduction of new matrix variables Q_i=Q_i^\top given by
+ + A_i^\top P + PA_i = Q_i,\qquad i=1,\ldots, r + yields
+ + \dot V(x) = x^\top Q_ix, + from which it follows that Q_i (for all i=1,\ldots,r) must satisfy
+ + x^\top Q_i x \leq 0,\qquad i=1,\ldots, r +
+for (Lyapunov) stability and
+ + x^\top Q_i x < 0,\qquad i=1,\ldots, r +
+for asymptotic stability.

+

As a matter of fact, we could proceeded without introducing new variables Q_i and just write the conditions directly in terms of A_i and P + x^\top (A_i^\top P + PA_i) x \leq 0,\qquad i=1,\ldots, r +
+for (Lyapunov) stability and
+ + x^\top (A_i^\top P + PA_i) x < 0,\qquad i=1,\ldots, r +
+for asymptotic stability.

+
+
+

Linear matrix inequality (LMI)

+

The conditions of quadratic stability that we have just derived are conditions on functions. However, in this case of a quadratic Lyapunov function and linear systems, the condition can also be written directly in terms of matrices. For that we use the concept of a linear matrix inequality (LMI).

+

Recall that a linear inequality is an inequality of the form \underbrace{a_0 a_1x_1 + a_2x_2 + \ldots + a_rx_r}_{a(x)} > 0.

+
+
+
+ +
+
+Linear vs. affine +
+
+
+

We could perhaps argue that as the function a(x) is an affine and not linear, the inequality should be called an affine inequality. However, the term linear inequality is well established in the literature. It can be perhaps justified by moving the constant term to the right-hand side, in which case we have a linear function on the left and a constant term on the right, which is the same situation as in the Ax=b equation, which we call linear without hesitation.

+
+
+

A linear matrix inequality is a generalization of this concept where the coefficients are matrices.

+

+\underbrace{A_0 + A_1x_1 + A_2x_2 + \ldots + A_rx_r}_{A(x)} \succ 0. +

+

Besides having matrix coefficients, another crucial difference is the meaning of the inequality. In this case it should not be interpreted component-wise but rather A(x)\succ 0 means that the matrix A(x) is positive definite.

+

Alternatively, the individual scalar variables can be assembled into matrices, in which case the LMI can have the form with matrix variables +F(X) = F_0 + F_1XG_1 + F_2XG_2 + \ldots + F_kXG_k \succ 0, + but the meaning of the inequality remains the same.

+

The use of LMIs is widespread in control theory. Here we formulate the LMI feasibility problem: does X=X^\top exist such that the LMI F(X)\succ 0 is satisfied?

+
+
+

CQLF as an LMI

+

+\begin{aligned} +P &\succ 0,\\ +A_1^\top P + PA_1 &\prec 0,\\ +A_2^\top P + PA_2 &\prec 0,\\ +& \vdots \\ +A_r^\top P + PA_r &\prec 0. +\end{aligned} +

+
+
+

Solving in Matlab using YALMIP or CVX

+
    +
  • Most numerical solvers for semidefinite programs (SDP) can handle nonstrict inequalities.

  • +
  • Enforce strict inequality by +\begin{aligned} +P &\succeq \epsilon I,\\ +A_1^\top P + PA_1 &\preceq \epsilon I,\\ +A_2^\top P + PA_2 &\preceq \epsilon I,\\ +& \vdots \\ +A_r^\top P + PA_r &\preceq \epsilon I. +\end{aligned} +

  • +
  • For LMIs with no affine term, multiply them “arbitrarily” to get +\begin{aligned} +P &\succeq I,\\ +A_1^\top P + PA_1 &\preceq I,\\ +A_2^\top P + PA_2 &\preceq I,\\ +& \vdots \\ +A_r^\top P + PA_r &\preceq I. +\end{aligned} +

  • +
+
+
+

Solution set of an LMI is convex

+
    +
  • and therefore if a solution P=P^\top \succ 0 exists such that +\begin{aligned} +A_1^\top P + PA_1 &\prec 0,\\ +A_2^\top P + PA_2 &\prec 0,\\ +& \vdots \\ +A_r^\top P + PA_r &\prec 0, +\end{aligned} +

  • +
  • then P also solves the convex combination +\left(\sum_{i=1}^r\alpha_i A_i\right)^\top P + P\left(\sum_{i=1}^r\alpha_i A_i\right) \prec 0, + where \alpha_1, \alpha_2, \ldots, \alpha_r \geq 0 and \sum_i \alpha_i = 1.

  • +
  • That means that necessarily every convex combination of the systems must be stable.

    +
      +
    • But it is only necessary, not sufficient.
    • +
  • +
  • Equivalently, we investigate stability of a linear differential inclusion +\dot x \in \mathcal{F}(x), + where \mathcal{F}(x) = \overline{\operatorname{co}}\{A_1x, A_2x, \ldots, A_rx\}.

  • +
+
+
+

What if quadratic LF is not enough?

+
    +
  • Quadratic Lyapunov function is a (multivariate) polynomial +\begin{aligned} +V(x) &= x^\top P x\\ +&= \begin{bmatrix}x_1 & x_2\end{bmatrix} \begin{bmatrix} p_{11} & p_{12}\\ p_{12} & p_{22}\end{bmatrix} \begin{bmatrix}x_1\\ x_2\end{bmatrix}\\ +&= p_{11}x_1^2 + 2p_{12}x_1x_2 + p_{22}x_2 +\end{aligned} +

  • +
  • How about a polynomial of a higher degree?

  • +
  • But how do we enforce positive definiteness?

  • +
+
+
+

Positive/nonnegative polynomials

+
    +
  • Is the polynomial p(\bm x), \; \bm x\in \mathbb R^n, positive (or nonnegative) on the whole \mathbb R^n? That is, we ask if +p(\bm x) > 0,\quad \text{or}\quad p(\bm x) \geq 0\; \forall \bm x\in\mathbb R^n. +

  • +
  • Example: p(\bm x)= 2x_1^4 + 2x_1^3x_2 - x_1^2x_2^2 + 5x_2^4\geq 0 for all x_1\in\mathbb R, x_2\in\mathbb R ?

  • +
  • Additionally, \bm x can be restricted to some \mathcal X\sub \mathbb R^n and we ask if
    + +p(\bm x) \geq 0 \;\forall\; \bm x\in \mathcal X. +

  • +
  • Semialgebraic sets \mathcal X are often considered. These are defined by polynomial inequalities such as +g_j(\bm x) \geq 0, \; j=1,\ldots, m. +

  • +
+
+
+

How can we check nonnegativity of polynomials?

+
    +
  • Gridding is not the way to go – we need conditions on the coefficients of the polynomial so that we can do some optimization later.

  • +
  • Example: Consider a univariate polynomial +p(x) = x^4 - 4x^3 + 13x^2 - 18x + 17. + Does it hold that p(x)\geq 0 \; \forall x\in \mathbb R ?

  • +
  • What if we learn that the polynomial can be written as +p = (x-1)^2 + (x^2 - 2x + 4)^2 +

  • +
  • Obviously, whatever the two squared polynomials are, after squaring they become nonnegative. And summing nonnegative numbers yields a nonnegative result. Let’s generalized this.

  • +
+
+
+

Sum of squares (SOS) polynomials

+
    +
  • If we can express the polynomial as a sum of squares of some other polynomials, the original polynomial is nonnegative
  • +
+

. . .

+

+\boxed{p(\bm x) = \sum_{i=1}^k p_i(\bm x)^2\; \Rightarrow \; p(\bm x) \geq 0,\; \forall \bm x\in \mathbb R^n.} +

+
    +
  • But the converse does not hold in general – not every nonnegative polynomial is SOS!

  • +
  • There are only three cases, for which SOS is a necessary and sufficient condition of nonnegativeness.

    +
      +
    • n=1: univariate polynomials. The degree (the highest power) d can be arbitrarily high (but even, obviously).
    • +
    • d = 2 and n is arbitrary: multivariate polynomials of degree two. Example, d=3 for p(\bm x) = x_1^2 + x_1x_2^2.
    • +
    • n=2 and d = 4: bivariate polynomials of degree 4 (at maximum).
    • +
  • +
  • For all other cases all we can say is that \mathrm{SOS} \Rightarrow p(\bm x)\geq 0.

  • +
  • Hilbert conjectured in 1900 in the 17th problem that every nonnegative polynomial can be written as a sum of squares of rational functions. This was later proved correct. It turns out, that this fact is not as useful as the SOS using polynomials because of impossibility to state apriori the bounds on the degrees of the polynomials defining those rational functions.

  • +
+
+
+

How to get an SOS representation of a polynomial (or prove that none exist)?

+
    +
  • Back to the univariate example first.

  • +
  • One of the two squared polynomials is x^2 - 2x + 4.

  • +
  • We can write it as x^2 - 2x + 4 = \underbrace{\begin{bmatrix}4 & -2 & 1\end{bmatrix}}_{\bm v^\top} \underbrace{\begin{bmatrix} 1 \\ x \\ x^2\end{bmatrix}}_{\bm z}. +

  • +
  • Then the squared polynomial can be written as +(x^2 - 2x + 4)^2 = \bm z^\top \bm v \bm v^\top \bm z. +

  • +
  • Note that the the product \bm v \bm v^\top is a positive semidefinite matrix of rank one.
    +

  • +
  • We can similarly express the second squared polynomial +x-1 = \underbrace{\begin{bmatrix} -1 & 1 & 0\end{bmatrix}}_{\bm v^\top} \underbrace{\begin{bmatrix} 1 \\ x \\ x^2\end{bmatrix}}_{\bm z} + and then +(x - 1)^2 = \bm z^\top \begin{bmatrix} -1 \\ 1 \\ 0\end{bmatrix} \begin{bmatrix} -1 & 1 & 0\end{bmatrix} \bm z. +

  • +
  • Summing the two squares we get the original polynomial. But while doing this, we can sum the two rank-one matrices. +\begin{aligned} +p(x) &= x^4 - 4x^3 + 13x^2 - 18x + 17\\ +&= \begin{bmatrix} 1 & x & x^2\end{bmatrix} \bm P \begin{bmatrix} 1 \\ x \\ x^2\end{bmatrix} +\end{aligned}, + where \bm P\succeq 0 is +\bm P = \begin{bmatrix} 4 \\ -2 \\ 1\end{bmatrix} \begin{bmatrix} 4 & -2 & 1\end{bmatrix} + \begin{bmatrix} -1 \\ 1 \\ 0\end{bmatrix} \begin{bmatrix} -1 & 1 & 0\end{bmatrix} +

  • +
  • The matrix that defines the quadratic form is positive semidefinite and of rank 2.

    +
      +
    • The rank of the matrix is given by the number of squared terms in the SOS decomposition.
    • +
  • +
  • In a general multivariate case we can proceed similarly. Just form the vector \bm z from all possible monomials: +z = \begin{bmatrix}1 \\ x_1 \\ x_2 \\ \\ \vdots \\ x_n\\ x_1^2 \\ x_1 x_2 \\ \vdots \\ x_n^2\\\vdots \\x_1x_2\ldots x_n^{2}\\ \vdots \\ x_n^d \end{bmatrix} +

  • +
+
+
+

But how to determine the coefficients of the matrix?

+
    +
  • Example: p(x_1,x_2)=2x_1^4 +2x_1^3x_2 − x_1^2x_2^2 +5x_2^4.

  • +
  • Define the vector \bm z as +\bm z = \begin{bmatrix} x_1^2 \\ x_1x_2 \\ x_2^2\end{bmatrix} +

  • +
  • Then it must be possible to write the polynomial as +\begin{aligned} +p(x_1,x_2)&=\begin{bmatrix} x_1^2 \\ x_1x_2 \\ x_2^2\end{bmatrix}^\top \begin{bmatrix} p_{11} & p_{12} & p_{13}\\ p_{12} & p_{22} & p_{23}\\p_{13} & p_{23} & p_{33}\end{bmatrix} \begin{bmatrix} x_1^2 \\ x_1x_2 \\ x_2^2\end{bmatrix}\\ +&= p_{11}x_1^4 + p_{33}x_2^4 \\ +&\quad + 2p_{12}x_1^3x_2 + 2p_{23}x_1x_2^3\\ +&\quad + (2p_{13} + p_{22})x_1^2x_2^2 +\end{aligned} +

  • +
  • This only gives 5 equations.

  • +
  • The sixth is the LMI condition P\succeq 0.

  • +
+
+
+

Searching for a positive polynomial Lyapunov function

+

. . .

+

V(x) \quad \text{is SOS}

+

. . .

+

Ooops, V(s) must be positive and not just nonnegative:

+

. . .

+

+\boxed{ +V(s) - \phi(x) \quad \text{is SOS},} + where \phi(x) = \gamma \sum_{i=1}^n\sum_{j=1}^{d} x_i^{2j} for some \gamma > 0.

+

. . .

+

+\boxed +{\left(\nabla V(x)\right)^\top f(x)\quad \text{is SOS}} +

+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/stability_via_common_lyapunov_function.html b/stability_via_common_lyapunov_function.html index da5bb75..bb4313e 100644 --- a/stability_via_common_lyapunov_function.html +++ b/stability_via_common_lyapunov_function.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Stability via multiple Lyapunov functions

+
+ + + +
+ + + + +
+ + + +
+ + +
+

Example: CQLF cannot be found

+
+

Example 4.11 from Lin & Antsaklis (2022)

+
+
+
+

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.
    • +
    • Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge University Press. Page 655.
    • +
  • +
+
+
+

Piecewise quadratic Lyapunov function

+
    +
  • 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:
  • +
+

. . .

+

+V_i(\bm x) = \bm x^\top \mathbf P_i \bm x > 0\quad \forall \;\bm x\in \Omega_i +

+

. . .

+

+\dot V_i(\bm x) = \bm x^\top \left( \mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i \right) \bm x < 0\quad \forall \;\bm x\in \Omega_i +

+
+
+

Using comparison functions and nonstrict inequalities

+

. . .

+

+\alpha_1 \bm x^\top \mathbf I \bm x \leq \bm x^\top \mathbf P_i \bm x \leq \alpha_2 \bm x^\top \mathbf I \bm x \quad \forall \;\bm x\in \Omega_i +

+

. . .

+

+\bm x^\top \left( \mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i \right) \bm x \leq -\alpha_3 \bm x^\top \mathbf I \bm x\quad \forall \;\bm x\in \Omega_i +

+
+
+

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, +

+

. . .

+

+\bm x^\top \left( \mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i \right) \bm x \leq -\gamma_i \bm x^\top \mathbf I \bm x, +

+

. . .

+

both for all \bm x\in \Omega_i, that is, all \bm x such that \bm x^\top \mathbf Q_i \bm x \geq 0.

+

. . .

+

+\mathbf P_i - \alpha_i \mathbf I - \mu_i \mathbf Q_i \succeq 0,\quad \mu_i \geq 0,\; \alpha_i > 0, +

+

. . .

+

+\mathbf A_i^\top \mathbf P_i + \mathbf P_i \mathbf A_i + \gamma_i \mathbf I + \xi_i \mathbf Q \preceq 0,\quad \mu_i \geq 0,\; \gamma_i > 0. +

+
+
+

More general sets

+

. . .

+

+\bm x^\top \mathbf Q \bm x + 2\mathbf r^\top \bm x + s \geq 0 +

+

. . .

+

+\begin{bmatrix} +\bm x^\top & 1 +\end{bmatrix} +\underbrace{ +\begin{bmatrix} +\mathbf Q & \mathbf r \\ \mathbf r^\top & s +\end{bmatrix}}_{\bar{\mathbf{Q}}} +\underbrace{ +\begin{bmatrix} +\bm x \\ 1 +\end{bmatrix}}_{\bar{\bm x}} +\geq 0 +

+

. . .

+

+\begin{bmatrix} +\mathbf Q & \mathbf r \\ \mathbf r^\top & s +\end{bmatrix} +\succeq 0 +

+
+
+

Example: affine subspace

+

. . .

+

+\mathbf c^\top \bm x + d \geq 0 +

+

. . .

+

+\begin{bmatrix} +\bm x^\top & 1 +\end{bmatrix} +\begin{bmatrix} +\mathbf 0 & \mathbf c \\ \mathbf c^\top & 2d +\end{bmatrix} +\begin{bmatrix} +\bm x \\ 1 +\end{bmatrix} +\geq 0 +

+

. . .

+

+\begin{bmatrix} +\mathbf 0 & \mathbf c \\ \mathbf c^\top & 2d +\end{bmatrix} +\succeq 0 +

+
+
+

But then the Lyapunov-like functions and system matrices must also be extended

+

. . .

+

+V(\bm x) = +\begin{bmatrix} +\bm x^\top & 1 +\end{bmatrix} +\underbrace{ +\begin{bmatrix} +\mathbf P & \mathbf P_{12} \\ \mathbf P_{12} & P_{22} +\end{bmatrix}}_{\bar{\mathbf P}} +\begin{bmatrix} +\bm x \\ 1 +\end{bmatrix} +

+

. . .

+

+\bar{\mathbf{A}} = +\begin{bmatrix} +\mathbf A & \mathbf 0 \\ \mathbf 0 & 0 +\end{bmatrix} +

+
+
+

Continuity conditions

+
    +
  • 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}. +

  • +
  • Collecting the terms with equal powers of \bm z. +\begin{aligned} +\mathbf F_{ij}^\top (\mathbf P_1 - \mathbf P_2) \mathbf F_{ij} &= 0, \\ +\mathbf F_{ij}^\top (\mathbf P_1 - \mathbf P_2) \mathbf l_{ij} + (\mathbf q_1-\mathbf q_2) &= 0, \\ +\mathbf l_{ij}^\top (\mathbf P_1 - \mathbf P_2)\mathbf l_{ij} + 2\mathbf l_{ij}^\top (\mathbf q_1-\mathbf q_2) + r_1-r_2 &= 0. +\end{aligned} +

  • +
+ + +
+ + Back to top
+ + +
+ + + + + + \ No newline at end of file diff --git a/stability_via_multiple_lyapunov_function.html b/stability_via_multiple_lyapunov_function.html index 60357ad..254498c 100644 --- a/stability_via_multiple_lyapunov_function.html +++ b/stability_via_multiple_lyapunov_function.html @@ -2,7 +2,7 @@ - + @@ -37,10 +37,10 @@ - + - + - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ +
+ +
+ + +
+ + + +
+ +
+
+

Literature

+
+ + + +
+ + + + +
+ + + +
+ + +

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.

+
+

Reachability analysis

+

[3] [4]

+
+
+

Barier certificates

+

[5]

+
+
+

Temporal logics

+

[6], [7]

+

[8], [9].

+ + + +
+ + Back to top

References

+
+
[1]
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
+
+
+
[2]
S. Mitra, Verifying Cyber-Physical Systems: A Path to Safe Autonomy. in Cyber Physical Systems Series. Cambridge, MA, USA: MIT Press, 2021. Available: https://sayanmitracode.github.io/cpsbooksite/about.html
+
+
+
[3]
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.
+
+
+
[4]
M. Althoff, N. Kochdumper, M. Wetzlinger, and T. Ladner, CORA 2024 Manual.” 2023. Accessed: Jan. 10, 2023. [Online]. Available: https://tumcps.github.io/CORA/data/archive/manual/Cora2024Manual.pdf
+
+
+
[5]
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.
+
+
+
[6]
C. Baier and J.-P. Katoen, Principles of Model Checking. Cambridge, MA, USA: MIT Press, 2008. Available: https://mitpress.mit.edu/books/principles-model-checking
+
+
+
[7]
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/
+
+
+
[8]
R. M. Murray, U. Topcu, and N. Wongpiromsarn, “Lecture 3 Linear Temporal Logic (LTL).” Belgrade (Serbia), Mar. 2020. Available: http://www.cds.caltech.edu/~murray/courses/eeci-sp2020/L3_ltl-09Mar2020.pdf
+
+
+
[9]
N. Wongpiromsarn, R. M. Murray, and U. Topcu, “Lecture 4 Model Checking and Logic Synthesis.” Belgrade (Serbia), Mar. 2020. Available: http://www.cds.caltech.edu/~murray/courses/eeci-sp2020//L4_model_checking-09Mar2020.pdf
+
+
+ + +
+ + + + + + \ No newline at end of file diff --git a/verification_references.html b/verification_references.html index baab2af..d3b343f 100644 --- a/verification_references.html +++ b/verification_references.html @@ -2,7 +2,7 @@ - + @@ -56,10 +56,10 @@ - + - + - + - +