Computing state trajectories using scans

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.

Like this:

One blogger likes this post.

Tags: Chaotic Maps, Dynamical Systems, Functional Programming, Haskell, Higher-Order Functions, List Processing, Logistic Map

This entry was posted on January 15, 2012 at 22:17 and is filed under Dynamical Systems, Haskell. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

Rod Carvalho

Leave a Reply Cancel reply

Follow “Rod Carvalho”