Every day we are making decisions under uncertainty. The weather forecast gives us a prediction of the weather and based on that we are choosing our clothes for the day; based on uncertain information from business news we trade stocks; or we eat the one-week old lunchbox after tasting it. Almost never do we have full and accurate information before we make a decision, but we can somehow estimate the risks of taking a decision under uncertainty.
Understanding the effect of the uncertainty is especially important in engineering applications, where the considered mechanical system should not exceed its safety limits, or a robot needs to navigate an unknown environment while avoiding collisions.
This could make us wonder "Is there a way to mathematically quantify the effect of the uncertainty such that we could make better decisions?".
In this interactive tutorial, we will investigate the effect of uncertainty in the form of disturbances on a dynamical system. Before we begin we need to clarify what a dynamical system is and what we mean by disturbance.
A dynamical system is a system that evolves with the passing of time. Dynamical systems appear in many fields, for example as cars and space rockets in engineering, the value of a stock at the stockmarket in economics, and the evolution of a population of potentially competing animals in biology.
Here, we focus on dynamical system that evolve in discrete time. This means that time is measured in discrete steps.
If we think about the size of a certain animal population, we might only measure that on a monthly or yearly basis, such that one time step would correspond to the passing of one month or year.
Furthermore, the state $x(k)\in\mathbb{R}^n$ of the system evolves according to linear dynamics, that is, $$ x(k+1) = Ax(k), $$ where $n$ is the dimension of the system state and $A\in\mathbb{R}^{n \times n}$ is the system matrix. The system matrix describes how the state evolves from one time step to another time step and its properties will become important later on, when we quantify the effect of disturbances on the system state $x(k)$.
As an example for a dynamical system we show the linearized population dynamics for the average population of rabbits and foxes in a certain area over a year. The solid line represents the continuous-time linearized dynamics based on the Lotka-Volterra equations, while the circle and diamonds markers show the linear discrete-time model. In this model we observe that as the population of foxes is in decline the rabbits thrive. This in turn provides more food for the foxes such that their population increases again. Due to that we get the oscillatory behavior of the populations for the chosen parameters.
The discrete-time model is evaluated only at discrete time instances, here every month. We see that our discrete-time model perfectly represents the continuous-time model at the evaluated times.
However, the dynamics might not always evolve according to our model above, which brings us to the next question.
Often there are disturbances acting on our system, which we denote as $d(k)$. This means that at each time step, $k$, the disturbance can act on our system such that the state does not obey the dynamical system equations from above. If we think of the animal population example from above, the presence of hunters could reduce the animal population or unmodelled environmental conditions can make the population thrive.
In our case, we look at disturbances that influence the system state $x(k)$ in an additive manner: $$ x(k+1) = Ax(k)+d(k) $$ A disturbance could also be seen a modelling error, such as a value for the $A$ matrix that does not reflect the true system. An example of this is believing that the $A$ matrix is equal to $\tilde{A}=A+\Delta A$. This leads to the disturbed system dynamics $$ x(k+1) = \tilde{A}x(k) = Ax(k)+\underbrace{\Delta A x(k)}_{=d(k)}, $$ where the disturbance is actually state dependent.
To show the effect of modelling errors, let us consider the population dynamics model from above again. Now we are not sure of the true parameters that describe the population dynamics of rabbits and foxes evolves.
The solid lines represent again the true dynamics but our discrete-time dynamics have slightly different parameters now. Therefore, the evaluation of the discrete-time model does not exactly match the true dynamics anymore due to the modelling error. Having some insight into the range of reasonable parameters can help us to quantify the potential modelling errors and how far off the state that we see could be from the true state of the system.
Another example is that the true dynamics are given by a nonlinear function and we approximate the dynamics via linearization of the nonlinear function.
Now that we know have our description of a dynamical system and how an additive disturbance influences it, let us state more concretely what we want to derive here.
Given a dynamical system under additive disturbance $$ x(k+1) = Ax(k)+d(k) $$ we want to determine which values the system state $x(k)$ can reach when the disturbance can take any value in a set $\mathbb{D}$ at each time step $k$.
More formally, we want to determine the sets $\mathbb{X}(k)$, such that $x(k)\in\mathbb{X}(k)$ for all $d(k)\in\mathbb{D}$ and $k\geq 0$.
The set $\mathbb{D}$ represents some knowledge that we might have on a potential disturbance. For example, in the population example, we might know that there are hunters in the area and we know that there are laws that prevent them from hunting too much. Hence, we could bound the potential disturbance on the animal population by that number.
Before we move on to the heart of this interactive tutorial, we want to clarify the assumption that we will make for the derivations coming up.
Based on the system dynamics above, we see that we already made two assumptions:
In the following, we make a third assumption:
This basically means that our model $x(k+1)=Ax(k)$ could also be correct at certain time steps, that is, when $d(k)=0$.
Equipped with the our goal let us move on to the main part of this tutorial.
This section contains the heart of this interactive tutorial. We will explain the main concepts on how to derive the sets $\mathbb{X}(k)$ and illustrate the intuition behind them with a scalar system. In a scalar system the state $x(k)$ is respresented by a real number and $A$ is not a matrix anymore but also a real number. Furthermore, we consider a disturbance that can take any value between -1 and 1. Mathematically, we write this as $d(k)\in \mathbb{D}=[-1,1]$ for all values of $k$. For the sake of simplicity, we also assume that at time step 0, the system state is at 0, which means $x(0)=0$.
The state of the scalar system can be interpreted as the temperature of a fridge with a temperature controller and the disturbances can, for example, be changes in the outside temperature or the opening of the fridge.
To determine the sets $\mathbb{X}(k)$, we will iterate through the dynamic equations and aim to find a pattern for the sets $\mathbb{X}(k)$. Using the dynamical equation, we determine that $$x(1) = A x(0) + d(0) = d(0).$$ Since we do not know the value of the disturbance at time step $0$, $d(0)$, the best guess that we can make is that the state of the system could have any value in the set $\mathbb{D}$, that is, $\mathbb{X}(1)=\mathbb{D}$. For our scalar system example this meamns $x(1)\in[-1,1]$.
Let us evolve the system further and look at the second time step $$ x(2) = A x(1) + d(1).$$ This gets a bit trickier now because we know neither the value of $x(1)$ nor the value of $d(1)$. But we know that $x(1) \in \mathbb{X}(1)=\mathbb{D}$, such that we start with determining the potential values of $Ax(1)$, which means that we want to determine the set $A\mathbb{X}(1)$. The set is $A\mathbb{X}(1)$ is defined as all the values in $\mathbb{X}(1)$ multiplied by $A$, which can be stated as $$ A\mathbb{X}(1) = \lbrace Ax\ |\ x\in\mathbb{X}(1) \rbrace. $$ In our scalar example $X(1)$ is an interval, such that $A\mathbb{X}(1)$ is an interval that contains all values of $\mathbb{X}(1)$ multiplied by $A$. Note that if $A<0$, we get $A\mathbb{X}(1)=[A, -A]$, which we can rewrite as $A\mathbb{X}(1)=[-|A|, |A|]$. Therefore, $A\mathbb{X}(1)=[-|A|,|A|],$ where we use the absolute value of $|A|$ to be able to have a minus sign at the lower bound.
In the figure below, you can see how $A\mathbb{X}(1)$ changes for different values of $A$.
How does $A\mathbb{X}(1)$ change with the value of $A$?
By playing around with the values of $A$ we observe that if $A<1$ the interval shrinks compared to $\mathbb{X}(1)$ until it starts growing again for $A<-1$. While in the scalar case, the set $\mathbb{X}(1)$ is stretched or compressed based on the value of $A$, in a higher dimensional scenario the set $\mathbb{X}(1)$ could also be rotated when it is multiplied by $A$.
Now that we know the set $A\mathbb{X}(1)$, let us determine the set in which $x(2)$ can lie.
Recall that $x(2) = A x(1) + d(1)$, where both $Ax(1)$ and $d(1)$ are unknown but we know the sets in which these values can lie in. Therefore, to determine the set $\mathbb{X}(2)$ we need to sum the sets $A\mathbb{X}(1)$ and $\mathbb{D}$ somehow together.
The mathematical operation to add two sets together is called a Minkowski sum, which is denoted by the operator $\oplus$. Mathematically, it is defined as follows $$ \mathbb{A} \oplus \mathbb{B} = \lbrace a + b\ \text{for}\ \text{all}\ a\in\mathbb{A},\, b\in\mathbb{B}\rbrace. $$ Intuitively, this means we pick any value $a\in\mathbb{A}$ and add all the values in $\mathbb{B}$ to it to get a new set, which is basically the set $\mathbb{B}$ shifted by the value $a$. If we do that for all values in $\mathbb{A}$ and take the union over the resulting sets, we obtain the Minkowski sum of the sets $\mathbb{A}$ and $\mathbb{B}$.
We are visualizing that concept for the scalar system example in the figure below, where we choose $A=0.5$.
The blue solid line interval represents $0.5\mathbb{X}(1)$ and the red dashed interval represents what happens when a single value of $0.5\mathbb{X}(1)$ is added to $\mathbb{D}$. Feel free to move the slider to see how the resulting red interval changes if a different value of $c\in0.5\mathbb{X}(1)$ is added to $\mathbb{D}$.
Finally, we need to take the union over all the possible red dashed sets in the figure above to obtain the Minkowski sum $A\mathbb{X}(1)\oplus\mathbb{D}$. In the figure below, we show the sets $\mathbb{X}(1)$, $A\mathbb{X}(1)$, and $\mathbb{D}$ for the scalar system.
Since we are only considering intervals, we observe that the lower bound and upper bound of $A\mathbb{X}(1)\oplus\mathbb{D}$ are $-|A|-1$ and $|A|+1$, respectively. Therefore, we have determined that $$ \mathbb{X}(2) = A\mathbb{X}(1)\oplus\mathbb{D} = [-|A|-1, |A|+1] $$ and that $x(2)\in\mathbb{X}(2)$.
Now if we go one time step further, we want to determine $\mathbb{X}(3)$ from $$ x(3) = Ax(2) + d(2). $$ Here, we can use again the Minkowski sum to write $$ \mathbb{X}(3) = A\mathbb{X}(2)\oplus \mathbb{D} = A^2\mathbb{X}(1)\oplus A\mathbb{D} \oplus \mathbb{D} = A^2\mathbb{D}\oplus A\mathbb{D} \oplus \mathbb{D}, $$ where we used the formulas for $\mathbb{X}(2)$ and $\mathbb{X}(1)$, which we found further above.
For the scalar system, we realize that $\mathbb{X}(2)$ is an interval again. So, we obtain $$ \mathbb{X}(3) = [-|A|^2-|A|-1,|A|^2+|A|+1] $$ such that $x(3)\in\mathbb{X}(3)$, where we followed the same ideas as above.
By iteratively applying the dynamics, we can show that $$ x(k)=\sum_{i=0}^{k-1}A^{k-1-i}d(i)$$ which leads to $$ \mathbb{X}(k)={\bigoplus}_{i=0}^{k-1}A^{k-1-i}\mathbb{D}, $$ where we use the big $\oplus$ similar to the sum operator $\sum$ but for Minkowski sums and not addition.
For the scalar system, we get $$ \mathbb{X}(k)={\bigoplus}_{i=0}^{k-1}A^{k-1-i}\mathbb{D}=\left[-\sum_{i=0}^{k-1}|A|^i,\sum_{i=0}^{k-1}|A|^i\right]. $$
Let us now focus only on the scalar system example to analyze the set $\mathbb{X}(k)$.
Since the upper and the lower bounds of the interval represent a geometric series, we can rewrite the interval as $$ \mathbb{X}(k)=\left[-\frac{1-|A|^{k}}{1-|A|},\frac{1-|A|^{k}}{1-|A|}\right], $$ whenever $|A|\neq 1$.
When $|A|=1$, the set $\mathbb{X}(k)$ is given by $$ \mathbb{X}(k)=\left[-\sum_{i=0}^{k-1}1^i,\sum_{i=0}^{k-1}1^i\right]=\left[-(k-1),k-1\right], $$ which also grows unbounded as $k$ grows large.
Thus, at each time step we know in which set $\mathbb{X}(k)$ the state $x(k)$ will lie if it is influenced by an additive disturbance and we have, therefore, quantified how the disturbance affects the dynamical system.
In the plot below you can change the value of $A$ and see how the set $\mathbb{X}(k)$ changes over $10$ time steps. An example trajectory of the system state $x(k)$ is also plotted, where at each time step $d(k)$ is randomly drawn from the set $\mathbb{D}=[-1,1]$.
When moving the slider in the figure above, we observe that the set $\mathbb{X}(k)$ grows if the absolute values of $A$ grows. Given that the set is described by $$ \mathbb{X}(k)=\left[-\frac{1-|A|^{k}}{1-|A|},\frac{1-|A|^{k}}{1-|A|}\right], $$ we can directly see that if $|A|>1$ the boundaries of $\mathbb{X}(k)$ grow infinitely large as $k$ grows. A system like this is called unstable, because the smallest disturbance can make the system state diverge to infinitely large value.
In the plot above we could also see that if $|A|<1$ then bounds seems to converge to a constant value. For example, setting $|A|=0.1$, we can see that the sets converge to approximately $[-1.111,1.111]$.
This set can also be calculated from the geometric series such that if the absolute value of $A$ is smaller than $1$, the interval converges to $$ \lim_{k\rightarrow \infty}\mathbb{X}(k)=\mathbb{X}_\infty = \left[-\frac{1}{1-|A|},\frac{1}{1-|A|}\right].$$ The systems with $|A|<1$ are called stable because the system state will remain in a bounded interval for all bounded disturbances.
These stability properties of the system $x(k+1)=Ax(k)+d(k)$ could already be observed in the first interactive figure, where we visualized $A\mathbb{X}(1)$ and saw that $A\mathbb{X}(1)$ would encompass all of $\mathbb{X}(1)$ if $|A|>1.
The set $\mathbb{X}_\infty$ possesses an interesting property when $|A|<1$. If we compare the upper bounds of $\mathbb{X}(k)$ and $\mathbb{X}_\infty$ then we observe that $$ \frac{1-|A|^{k}}{1-|A|} \leq \frac{1}{1-|A|} \Leftrightarrow 1-|A|^{k} \leq 1$$ holds for all values of $k$. Similarly, we can show that the lower bound of $\mathbb{X}_\infty$ is always lower than the lower bound of $\mathbb{X}(k)$. This means that the set $\mathbb{X}(k)$ is contained in $\mathbb{X}_\infty$ for all values of $k$, mathematically, $\mathbb{X}(k) \subseteq \mathbb{X}_\infty$ for all $k\geq 0$. We call the set $\mathbb{X}_\infty$ positively invariant because if $x(k)$ lies inside $\mathbb{X}_\infty$ then $x(k+1)$ will also lie in $\mathbb{X}_\infty$.
Below we show eight trajectories of a scalar system with $A=0.5$ over a horizon of $10$ time steps. We also plot the upper and lower bounds of $\mathbb{X(k)}$ and $\mathbb{X}_\infty$ as empty and filled triangles, respectively.
Looking at the upper and lower bounds of $\mathbb{X(k)}$ and $\mathbb{X}_\infty$, we observe that for small $k$ they are different, but as $k$ increases the bounds converge to each other.
Furthermore, as expected, all eight trajectories lie within these bounds.
What about higher dimensional systems?
The last examples above deriving $\mathbb{X}_\infty$ used the scalar system and this makes us wonder if something similar holds for higher dimensional systems as well.
If we look at the scalar example, the main property of $A$ to guarantee the convergence to a finite set $\mathbb{X}_\infty$ is that $|A|^k$ approaches zero as $k$ approaches infinity. This lead us to the condition $|A|<1$.
If $A$ is a quadratic matrix, what property would this matrix need such that $A^k$ approaches the zero matrix as $k$ increases. A matrix with this property has a spectral radius smaller than 1, which is denoted as $\rho(A)<1$. The spectral radius represents the largest absolute value of the eigenvalues of $A$. Similarly, if $\rho(A)>1$ then $A^k$ grows infinitely large. Hence, a system is stable if the spectral radius of its system matrix $A$ is smaller than 1.
To connect it back to our scalar system, note that if $A$ is a one dimensional matrix then it is actually simply a real number. Thererfore, the eigenvalue of $A$ is simply $A$ itself and, thus, $\rho(A)=|A|$.
In this interactive tutorial, we have looked at dynamical systems under additive disturbances and wanted to analyze the effect of the disturbance on the system state.
To do so, we looked at the sets in which the states can lie at each time step given that the disturbance lies in a bounded set. This has led us to get familiar with Minkowski sums and enabled us to find an expression for the sets that the states can lie in at each time step.
Then we discussed the stability of a dynamical system and showed that if the dynamical system is stable, there exists a positively invariant set, which contains all the other sets. When we analyze where a disturbance could push the state of the system, we are interested in the minimal positively invariant set, because it describes the smallest set in which the state could be at all times. This set is also often called minimal robust positively invariant set, because once we know where the state could be at any time we have quantified the effect of the disturbance on the system. This enables us to make a decision that is robust to the effect of the disturbance on the system. However, it must also be said that finding the minimal robust positively invariant set is often non-trivial for higher dimensional systems. Therefore, research has been done to find the overapproximations of this set. For example, Raković et al. provide a method for when the disturbance set can be represented by a convex polytope.
Furthermore, we also looked only at linear systems in discrete time. There is a wealth of literature addressing the estimation of the robust positive invariant sets for nonlinear systems as well.
I would like to thank everyone who has read through this tutorial and provided me valuable feedback. Furthermore, I would like to thank the team behind Voila Dashboards for developing a method to turn Jupyter notebooks into web applications.