Posts Tagged ‘Dynamical Systems’

Deciding the existence of equilibrium points

February 25, 2012

Suppose we are given a dynamical system of the form

\dot{x} (t) = f ( x (t) )

where x : \mathbb{R} \to \mathbb{R}^n is the state trajectory, and f : \mathbb{R}^n \to \mathbb{R}^n is a known vector field [1]. We now introduce the following definition:

Definition: A point x^* \in \mathbb{R}^n is an equilibrium point of the given dynamical system if the state being at x^* at some time t^* implies that the state will remain at x^* for all future time, i.e., x (t) = x^* for all t \geq t^*. In other words, an equilibrium point x^* is a point of zero flow, i.e., f ( x^* ) = 0.

If you happen to dislike the word “flow”, feel free to replace it with the word “velocity”. A point of zero velocity is thus a stationary one. Like a properly-trained poodle, it stays put.

One can find the equilibrium points of a given dynamical system by computing the real roots of the vector equation f (x) = 0. However, suppose that we are not interested in computing the equilibrium points; instead, all we would like to know is whether any such points exist. Thus, we arrive at the following decision problem:

Problem: Given a vector field f : \mathbb{R}^n \to \mathbb{R}^n, we would like to decide whether the dynamical system \dot{x} (t) = f ( x (t) ) has at least one equilibrium point. This is equivalent to determining the truth value of the formula \exists x \left( f (x) = 0 \right), where x \in \mathbb{R}^n.

Since both x and f (x) are n-dimensional vectors, x = (x_1, \dots, x_n) and f (x) = (f_1 (x), \dots, f_n (x)), we can rewrite the existentially-quantified formula above in the following form

\displaystyle\exists x_1 \exists x_2 \dots \exists x_n \left( \bigwedge_{i=1}^n f_i (x_1, x_2, \dots, x_n) = 0\right).

Can one even determine the truth value of the existentially-quantified formula above? In order to avoid colliding against the brick wall of undecidability, let us restrict our attention to polynomial dynamical systems (i.e., dynamical systems in which the vector field f : \mathbb{R}^n \to \mathbb{R}^n is polynomial). For such systems, one can use quantifier elimination to decide the existence of equilibrium points.

__________

Example

Consider the following polynomial dynamical system

\begin{array}{rl} \dot{x}_1 &= \mu - x_1^2\\ \dot{x}_2 &= - x_2\end{array}

taken from section 2.7 in Khalil’s book [1]. Note that the two first-order ordinary differential equations above are decoupled, i.e., they are of the form \dot{x}_i = f_i (x_i). Note also that we have a free parameter \mu \in \mathbb{R}.

We now compute the equilibrium points of the dynamical system by solving the (scalar) polynomial equations f_i (x_i) = 0. We thus obtain x_1^2 = \mu and x_2 = 0. For \mu > 0, we have two equilibrium points at (\sqrt{\mu}, 0) and (-\sqrt{\mu}, 0). For \mu = 0, we have one equilibrium point at (0,0). Finally, for \mu < 0, we have no equilibrium points, as the equation x_1^2 = \mu has no real roots when \mu < 0. The existence of equilibrium points depends on parameter \mu. Stating that this dynamical system has at least one equilibrium point is the same as saying that the following existentially-quantified formula

\displaystyle\exists x_1 \exists x_2 \left( x_1^2 - \mu = 0 \land x_2 = 0\right)

evaluates to \text{True}. By visual inspection of the formula, we conclude that the formula (which has two bound variables, x_1 and x_2, and one free variable \mu) evaluates to \text{True} when \mu \geq 0. Let us confirm this using the following REDUCE + REDLOG script:

% decide the existence of equilibrium points

load_package redlog;
rlset ofsf;

% define polynomial vector field
f1 := mu - x1**2;
f2 := -x2; 

% define existentially quantified formula
phi := ex({x1,x2}, f1=0 and f2=0);

% perform quantifier elimination
rlqe phi;

end;

which outputs \mu \geq 0. Hence, if the parameter \mu is nonnegative, then there will always exist at least one equilibrium point.

__________

References

[1] Hassan K. Khalil, Nonlinear Systems, 3rd edition, Prentice Hall.

Tags:, , , , , , ,
Posted in Dynamical Systems | Leave a Comment »

Computing state trajectories using scans

January 15, 2012

Consider a class of dynamical systems of the form

x_{k+1} = f (x_k, u_k )

where x_k is the state, u_k is the input, and f : \mathcal{X} \times \mathcal{U} \to \mathcal{X} is the state-transition function, where \mathcal{X} and \mathcal{U} are the state set and input set, respectively.

Given an initial state x_0 and a (possibly infinite) input sequence (u_0, u_1, u_2, \dots), we would like to compute the state trajectory (x_0, x_1, x_2, x_3, \dots). Iterating, we obtain

\begin{array}{rl} x_1 &= f (x_0, u_0)\\ x_2 &= f (x_1, u_1) = f (f (x_0, u_0), u_1)\\ x_3 &= f (x_2, u_2) = f (f (x_1, u_1), u_2) = f (f (f (x_0, u_0), u_1), u_2)\end{array}

and so on. Let us pause to take a look at the following equation

x_3 = f (f (f (x_0, u_0), u_1), u_2).

If you happen to be acquainted with functional programming, you will probably recognize this pattern. It is a left fold!! Hence, for every integer N \geq 1, we can compute x_N in Haskell using function foldl, as follows

\text{foldl} \, f \, x_0\, [u_0, u_1, \dots, u_{N-1}].

However, we would like to obtain the intermediate states x_1, x_2, \dots, x_{N-1} as well. We can compute the entire state trajectory using function scanl instead

\text{scanl} \, f \, x_0\, [u_0, u_1, \dots, u_{N-1}]

which returns a (finite) list of states [ x_0, x_1, \dots, x_N]. Amazingly (or perhaps not), since Haskell is lazy, we can obtain a state trajectory of infinite length corresponding to a given initial state and a given input sequence of infinite length!

__________

Example

Consider the logistic map with a control input

x_{k+1} = r x_k (1 - x_k) + u_k

where we will use r = 3.8, for which the dynamics (of the unforced system) are chaotic. Let the initial state be x_0 = 0 (note that x = 0 is a fixed-point of the unforced logistic map), and let the input sequence be a discrete-time sinusoid of infinite length

u_k = \displaystyle\frac{1}{5} \sin \left(\frac{\pi}{10} k\right)

whose period is 20. The following Haskell program computes the infinite-length state trajectory using the aforementioned scanl trick:

-- define logistic map with control input
f :: Floating a => a -> a -> a
f x u = (3.8 * x * (1 - x)) + u

-- frequency of the input sinusoid
omega :: Floating a => a
omega = pi / 10

-- generate sinusoidal control input
us :: Floating a => [a]
us = map ((0.2*) . sin . (omega*) . fromIntegral) [0..]

-- initial state
x0 :: Floating a => a
x0 = 0.0

-- compute the state trajectory
xs :: Floating a => [a]
xs = scanl f x0 us

where we used map, function composition, and an infinite list of non-negative integers to generate the sinusoidal input.

To be more precise, the program above does not quite compute the state trajectory, it merely “promises” to compute it when needed. After running the program in GHCi, we can obtain the first 20 values of the input sequence and the first 21 values of the state sequence using function take, as follows

*Main> take 20 us
[0.0,6.180339887498948e-,0.11755705045849463,
0.1618033988749895,0.1902113032590307,0.2,
0.19021130325903074,0.1618033988749895,
0.11755705045849466, 6.180339887498951e-2,
2.4492127076447546e17,-6.180339887498946e2,
-0.11755705045849461,-0.16180339887498948,
-0.1902113032590307,-0.2,-0.19021130325903074,
-0.1618033988749895,-0.11755705045849468,
-6.180339887498953e-2]
*Main> take 21 xs
[0.0,0.0,6.180339887498948e-2,0.33789525775595064,
1.0119471985345525,0.14426955372699996,
0.6691322284587664,1.031509602586003,
3.8294059838691885e-2,0.2575020247735924,
0.7883433805171414,0.6340607606653986,
0.8199019084343063,0.44356147166584237,
0.7760924326990134,0.4701259774450641,
0.7466086625502713,0.5286885334506017,
0.7850690797091344,0.5236383047578965,
0.8860732772080671]

In other words, the input and state sequences are only computed when the expressions with take are evaluated by the interpreter. We can plot these sequences in MATLAB:

Last but not least, here is the MATLAB code that generates the plot:

figure; hold on;
stem([0:1:19],u,'b');
stem([0:1:20],x,'r');
legend('input','state');
xlabel('k (time)'); ylabel('amplitude');
axis([0 21 -0.25 1.20]);

where (row) vectors u and x are the lists produced by GHCi.

__________

Lastly, please do note that the class of dynamical systems we have considered in this post is extremely broad. It encompasses difference equations obtained from the time-discretization of differential equations, IIR filters, finite and infinite automata, etc.

Tags:, , , , , ,
Posted in Dynamical Systems, Haskell | Leave a Comment »

Safety verification via barrier certificates

June 19, 2011

Suppose that we are given a dynamical system of the form \dot{x} (t) = f ( x (t) ), where x : [0, \infty) \to \mathbb{R}^n is the state trajectory, and f : \mathbb{R}^n \to \mathbb{R}^n is a known vector field [1]. Suppose further that we are given an initial set \mathcal{X}_0 \subset \mathbb{R}^n, and an unsafe set \mathcal{X}_u \subset \mathbb{R}^n.

If there exists a state trajectory x : [0, \infty) \to \mathbb{R}^n with initial condition x(0) \in \mathcal{X}_0 such that x(T) \in \mathcal{X}_u for some T \geq 0, we say that the system is unsafe. An example of an unsafe system is depicted below

If there exists no state trajectory x : [0, \infty) \to \mathbb{R}^n with initial condition x(0) \in \mathcal{X}_0 such that x(T) \in \mathcal{X}_u for some T \geq 0, we say that the system is safe. Note that if \mathcal{X}_0 \cap \mathcal{X}_u \neq \emptyset, we can conclude immediately that the system is unsafe, as there exists an initial condition x(0) \in \mathcal{X}_u.

The safety verification problem can then be stated as follows:

Given a vector field f : \mathbb{R}^n \to \mathbb{R}^n, an initial set \mathcal{X}_0 \subset \mathbb{R}^n and an unsafe set \mathcal{X}_u \subset \mathbb{R}^n, decide whether the dynamical system \dot{x} = f (x) is safe.

To be more precise, safety is not a property of the dynamical system \dot{x} = f (x) per se, but rather of the ordered triple (f, \mathcal{X}_0, \mathcal{X}_u).

How do we go about solving the safety verification problem? The most obvious approach would be to “propagate” the initial set \mathcal{X}_0 forwards in time, i.e., to compute (forward) reachable sets. However, computing exact reachable sets is seldom possible, which forces us to compute approximations of reachable sets.

When performing safety verification, we are asking a binary question: is (f, \mathcal{X}_0, \mathcal{X}_u) safe? Hence, computation of the reachable sets seems to provide more information than what we actually do need. This is a personal opinion, not a fact.

__________

Barrier certificates

Probably inspired on Lyapunov stability theory [1], Stephen Prajna introduced a few years ago a new method for safety verification of continuous-time dynamical systems that relies on what he termed barrier certificates [2]. A barrier certificate is a function of state B : \mathbb{R}^n \to \mathbb{R} that satisfies the following inequalities

\begin{array}{rl} B(x) \leq 0 \qquad{} \forall x \in \mathcal{X}_0 \\ B(x) > 0 \qquad{} \forall x \in \mathcal{X}_u\\ \displaystyle\frac{\partial B}{\partial x} (x) f (x) \leq 0 \qquad{} \forall x \in \mathbb{R}^n\end{array}

Let us now introduce a function of time b (t) := B (x (t)) and fix the initial condition x(0). Taking the derivative of b with respect to time, we obtain

\dot{b} (t) = \displaystyle\frac{\partial B}{\partial x} (x(t)) \dot{x} (t) = \displaystyle\frac{\partial B}{\partial x} (x(t)) f (x(t)) \leq 0

which follows from the non-positivity of the Lie derivative of B along the flow of f. From \dot{b}(t) \leq 0, we conclude that b(t) \leq b(0) for all t \geq 0, i.e., we have that B(x(t)) \leq B(x(0)) for all t \geq 0. Hence, the following sublevel set

\Omega_0 := \{ x \in \mathbb{R}^n \mid B(x) \leq B(x(0))\}

is positively invariant [1]. In other words, a solution starting in \Omega_0 will remain in \Omega_0 for all t \geq 0. Since B(x(0)) \leq 0 for all x(0) \in \mathcal{X}_0, we conclude that B(x(t)) \leq 0 for all t \geq 0. However, since B(x) > 0 for all x \in \mathcal{X}_u, we can finally conclude that the safety of (f, \{x(0)\}, \mathcal{X}_u) is guaranteed. Alternatively, one could argue that, since, by construction, \Omega_0 \cap \mathcal{X}_u = \emptyset, then safety is certified.

Pictorially, we start with the following level sets

and once we know that \Omega_0 is a positively invariant set (from B(x(t)) \leq B(x(0))), we conclude that the state trajectory must live in \Omega_0, which is painted in dark gray below

Please do recall that B(x(0)) \leq 0, and also note that the unsafe set is contained in the superlevel set \{ x \in \mathbb{R}^n \mid B(x) > 0\}.

So far we assumed that the initial condition was fixed. Henceforth, we shall make the initial condition free to take any values in the initial set \mathcal{X}_0. Let us define

b_0 := \displaystyle\sup_{x_0 \in \mathcal{X}_0} B(x_0)

so that we get an upper bound on B(x(0)), as follows

B(x(t)) \leq B(x(0)) \leq b_0 \leq 0

where we used the inequality B(x) \leq 0 for all x \in \mathcal{X}_0 to conclude that b_0 \leq 0. This allows us to conclude that the sublevel set

\tilde{\Omega}_0 := \{ x \in \mathbb{R}^n \mid B(x) \leq b_0\}

is positively invariant. Since \tilde{\Omega}_0 \cap \mathcal{X}_u = \emptyset, safety of (f, \mathcal{X}_0, \mathcal{X}_u) is guaranteed. Note that the sublevel set \tilde{\Omega}_0 is the union of all possible \Omega_0 sublevel sets, i.e.,

\tilde{\Omega}_0 = \displaystyle\bigcup_{x_0 \in \mathcal{X}_0} \{ x \in \mathbb{R}^n \mid B(x) \leq B(x_0)\}.

This new sublevel set is depicted below, painted in dark gray

Note that \mathcal{X}_0 \subset \tilde{\Omega}_0. Essentially, the zero level set of the barrier certificate B separates an unsafe set, \mathcal{X}_u, from all system trajectories starting in the initial set \mathcal{X}_0, which are contained in \tilde{\Omega}_0, thus providing a proof of safety without requiring the computation of reachable sets. Personally, I like to think of this approach as a sophisticated generalization of the well-known separating hyperplane theorem.

At this point, you may be thinking that this framework is too good to be true. Instead of computing the flow or the reachable sets, all we need to do is find a function B : \mathbb{R}^n \to \mathbb{R} that satisfies three inequalities, and then (voilà!) infinite-time safety of (f, \mathcal{X}_0, \mathcal{X}_u) is certified! As it usually happens, the devil is in the details…

If finding Lyapunov functions [1] for low-dimensional dynamical systems is often “hard”, why should we expect the hunt for barrier certificates to be an “easy” endeavor? Quel panache! Are we merely transferring the complexity of the problem from the computation of reachable sets to the satisfaction of three inequalities? Will we fail to find barrier certificates due to excessive conservativeness? Should we abandon all hope? Please be patient, dear reader, for we will address these issues and alleviate your doubts in future posts ;-)

Next installment will be devoted to an example. до скорой встречи!

_________

References

[1] Hassan K. Khalil, Nonlinear Systems, 3rd edition, Prentice Hall.

[2] Stephen Prajna, Optimization-based methods for nonlinear and hybrid systems verification, Dissertation (Ph.D.), California Institute of Technology, 2005.

Tags:, , , , , ,
Posted in Dynamical Systems | 1 Comment »

Sailing in state space

April 11, 2011

Consider the following controlled (autonomous) dynamical system [1]

\dot{x} (t) = f (x (t), u (t))

where x : [0, \infty) \to \mathbb{R}^n is the state trajectory, u : [0,\infty) \to \mathcal{U} is the control input history, \mathcal{U} \subseteq \mathbb{R}^m is the set of admissible control inputs, and f : \mathbb{R}^n \times \mathcal{U} \to \mathbb{R}^n is a (known) vector field.

If we know the control input u (t) for all t \geq 0, then we can integrate the ODE above with initial condition x(0) = x_0 and obtain the state trajectory x : [0, \infty) \to \mathbb{R}^n, as follows

x (t) = x_0 + \displaystyle \int_{0}^{t} f \left(x (\tau), u (\tau) \right) \mathrm{d} \tau

which can be represented pictorially by a single streamline in \mathbb{R}^n

I am (quite explicitly) alluding to fluid flow. Just like a cork on the ocean will follow a certain path depending on the velocity field (i.e., “ocean currents”), a point in state space will flow along a certain streamline.

__________

Let us fix the control input, i.e., u (t) = \bar{u} for all t \geq 0, and define

\begin{array}{rl} v^{(\bar{u})} : \mathbb{R}^n & \to \mathbb{R}^n\\ x & \mapsto f(x, \bar{u})\\\end{array}

For the sake of simplicity, let us consider \mathcal{U} = \{-1, +1\}. Since \bar{u} can only take two values, we have two vector fields, v^{(+1)}(x) := f(x, +1) and v^{(-1)}(x) := f(x, -1). Integrating the following ODEs

\begin{array}{rl} \dot{x} (t) &= v^{(+1)} ( x(t) )\\ \dot{x} (t) &= v^{(-1)} ( x(t) )\\\end{array}

for various initial conditions, we obtain two families of streamlines, as depicted below

For each fixed control input and collection of initial conditions, we will have a family of streamlines in state space.

Imagine that we have the following control input history

u (t) = \begin{cases} +1, & t \in [0, t^{\star})\\ -1, & t \geq t^{\star}\\\end{cases}

where we toggle the control input at time t^{\star} > 0, thus interrupting the flow along a “blue” streamline and initiating the flow along a “pink” streamline. Pictorially, we have

If the control input is fixed, say, \bar{u} = +1, then given an initial condition x(0), we will flow along the “blue” streamline that passes through x(0). This streamline is 1-dimensional and does not allow us to “explore” the n-dimensional state space. However, if we switch between \bar{u} = +1 and \bar{u} = -1, then we are able to “travel” to other regions of the state space. Allowing the control input to take values in \mathcal{U} = [-1,+1] would give us even more freedom.

Lastly, we arrive at a most childish idea: in some cases, controller design can be viewed as shaping the streamlines differently in different parts of the state space. This is childish because it is purely conceptual. One obviously cannot design optimal controllers by doodling with crayons.

__________

References

[1] Hassan K. Khalil, Nonlinear Systems, 3rd edition, Prentice Hall.

Tags:, , , ,
Posted in Control Theory, Dynamical Systems | 2 Comments »

Follow

Get every new post delivered to your Inbox.

Join 77 other followers