Dynamic Mode Decomposition for Spatiotemporal Traffic Speed Time Series in Seattle Freeway

Spatiotemporal traffic data analysis is an emerging area in intelligent transportation systems. In the past few years, data-driven machine learning models have provided new dimensions for understanding real-world data, building data computing paradigm, and supporting real-world applications.

In this blog post, we plan to:

introduce a publicly available traffic flow data in Seattle, USA,
design an easy-to-use data set (as a toy example) for traffic flow analysis,
perform dynamic mode decomposition on the toy example and discuss the interpretation of results.

Seattle Freeway Traffic Speed Data Set

This is a traffic speed data set collected by the inductive loop detectors deployed on freeways in Seattle, USA. The data set is publicly available on GitHub (see https://github.com/zhiyongc/Seattle-Loop-Data). As shown in Figure 1, the freeways include I-5, I-405, I-90, and SR-520. The speed information at a milepost is averaged from multiple loop detectors on the mainlanes in the same direction. In the data set, there are 323 loop detectors. The time resolution of speed information is 5 minutes, meaning that we have 288 time intervals per day or saying 288 data points per day for each loop detectors.

Figure 1. The map of freeways in Seattle area and their inductive loop detector locations. Each blue icon indicate the loop detectors at a milepost. The source is from GitHub. — **Figure 1**. The map of freeways in Seattle area and their inductive loop detector locations. Each blue icon indicate the loop detectors at a milepost. The source is from GitHub.

We can take a look at the data file as shown in Figure 2. The speed data is in the form of matrix. The row corresponds to each specific time interval, which is given by the time stamp like 2015–01–01 00:00:00 and 2015–01–01 00:05:00. The column refers to each loop detector ID.

Figure 2. Some traffic speed data points in the data set. The source is from GitHub. — **Figure 2**. Some traffic speed data points in the data set. The source is from GitHub.

In this blog post, we design a toy example data from this data set, and create a subset as toy_data.npy. The subset is available at our GitHub repository:

https://github.com/xinychen/transdim/blob/master/datasets/Seattle-data-set/toy_data.npy

Figure 3. Heatmap of traffic speed data in the I-5 freeway in Seattle. The time period is 6:00 am - 12:00 pm (i.e., 72 time intervals) on January 5, 2015. There are 75 loop detectors in this subset, starting from ID 166 to ID 240. — **Figure 3**. Heatmap of traffic speed data in the I-5 freeway in Seattle. The time period is 6:00 am – 12:00 pm (i.e., 72 time intervals) on January 5, 2015. There are 75 loop detectors in this subset, starting from ID 166 to ID 240.

Figure 3 shows the speed heatmap of traffic speed data (subset), which is of size 75-by-72 (i.e., 75 loop detectors and 72 time intervals). The subset covers the traffic state of the I-5 freeway during the morning rush hours. As can be seen, from (roughly) loop detector 176 to (roughly) loop detector 196, it shows clear low traffic speeds referring to traffic congestion. In other loop detectors, traffic speeds are relatively high.

To draw this figure as you have prepared the subset, you can try the following Python code:

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

mat = np.load('toy_data.npy')

plt.style.use('bmh')
plt.rcParams['font.family'] = 'Arial'

fig = plt.figure(figsize = (4, 4))
sns.heatmap(mat, cmap = 'YlGnBu', cbar_kws={'label': 'Traffic speed'}, vmin = 0, vmax = 75)
plt.xticks(np.arange(0, 72 + 1, 12), np.arange(4 * 288 + 72, 4 * 288 + 144 + 1, 12), rotation = 0)
plt.yticks(np.arange(0.5, 75.5, 10), np.arange(166, 241, 10), rotation = 0)
plt.show()

Analyze Traffic Speed Time Series

In Figure 3, we take into account a 75-by-72 subset of traffic speed. We can also draw some time series curves of traffic speed. Figure 4 shows three traffic speed time series of loop detectors 186, 196, and 236, respectively. We can see that both loop detectors 186 and 196 with relatively lower traffic speed suffer from traffic congestion.

Figure 4. Traffic speed time series of loop detectors 186, 196, and 236. — **Figure 4**. Traffic speed time series of loop detectors 186, 196, and 236.

To draw this figure, please try the following Python code:

import matplotlib.pyplot as plt

fig = plt.figure(figsize = (6, 6))
i = 1
for loc in [185, 195, 235]:
    ax = fig.add_subplot(3, 1, i)
    plt.plot(mat[loc - 166, :], color = "#006ea3", linewidth = 2.5)
    plt.xticks(np.arange(0, 72 + 1, 12), 
               np.arange(4 * 288 + 72, 4 * 288 + 144 + 1, 12))
    if i != 3:
        plt.setp(ax.get_xticklabels(), visible = False)
    plt.grid(axis = 'both', linestyle = "--", linewidth = 0.1, color = 'gray')
    plt.ylabel('Speed')
    plt.title('Loop detector {}'.format(loc + 1))
    ax.tick_params(direction = "in")
    ax.set_xlim([0, 72])
    ax.set_ylim([0, 70])
    i += 1
plt.show()

Reproduce Dynamic Mode Decomposition (DMD) on the Subset

DMD is a data-driven dimensionality reduction method for time series data. For any multivariate time series data, DMD can compute a set of dynamic modes in which each mode is associated with temporal behaviors. In particular, DMD allows one to interpret temporal behaviors of data with physically meaningful modes [1]. For an in-depth discussion of DMD, please check out [1]. We use the Python code in [1] for the following analysis.

import numpy as np

def DMD(data, r):
    """Dynamic Mode Decomposition (DMD) algorithm."""

    ## Build data matrices
    X1 = data[:, : -1]
    X2 = data[:, 1 :]
    ## Perform singular value decomposition on X1
    u, s, v = np.linalg.svd(X1, full_matrices = False)
    ## Compute the Koopman matrix
    A_tilde = u[:, : r].conj().T @ X2 @ v[: r, :].conj().T * np.reciprocal(s[: r])
    ## Perform eigenvalue decomposition on A_tilde
    eigval, eigvec = np.linalg.eig(A_tilde)
    ## Compute the coefficient matrix
    Psi = X2 @ v[: r, :].conj().T @ np.diag(np.reciprocal(s[: r])) @ eigvec

    return eigval, eigvec, Psi

In the case, r is the predefined low rank of DMD. eigval and eigvec correspond to the eigenvalues and the eigenvectors of Koopman matrix in DMD. In particular, eigval can be interpreted as DMD spectrum, appearing as complex conjugate pairs. Psi refers to the dynamic modes.

DMD Interpretation of the Results

Eigenvalues and DMD Spectrum

Here, we use the DMD function as mentioned above and evaluate the DMD model with rank 5. The Python code is given by

r = 5
eigval, eigvec, Psi = DMD(mat, r)

Then, we use the following Python code to draw the DMD spectrum:

import matplotlib.pyplot as plt

fig = plt.figure(figsize = (4, 4))
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal', adjustable = 'box')
plt.plot(eigval.real, eigval.imag, 'o', color = 'red', markersize = 6)
circle = plt.Circle((0, 0), 1, color = 'blue', linewidth = 2.5, fill = False)
ax.add_patch(circle)
plt.grid(axis = 'both', linestyle = "--", linewidth = 0.1, color = 'gray')
ax.tick_params(direction = "in")
plt.xlabel('Real axis')
plt.ylabel('Imaginary axis')
plt.show()

Figure 5 shows the DMD spectrum of data. Each eigenvalue explains the dynamic behavior of its corresponding dynamic mode. We can interpret eigenvalues as follows [2],

If the eigenvalue has a non-zero imaginary part, then there is oscillation in the corresponding dynamic mode.
If the eigenvalue is inside the unit circle, then the dynamic mode is decaying.
If the eigenvalue is outside the unit circle, then the dynamic mode is growing.

In Figure 5, the dots are close to or on the unit circle. There are four eigenvalues that are inside the unit circle, while one eigenvalue is on the unit circle whose imaginary part is 0. We can also see that two eigenvalue pairs are symmetric in the imaginary axis. The most dynamic modes are oscillatory and decaying. This result is consistent with spatiotemporal patterns of low traffic speeds in the morning rush hours.

Figure 5. Eigenvalues of DMD on the 75-by-72 data with rank 5. There are 5 red dots showing both real and imaginary axis of 5 eigenvalues. The blue circle is an unit circle. — **Figure 5**. Eigenvalues of DMD on the 75-by-72 data with rank 5. There are 5 red dots showing both real and imaginary axis of 5 eigenvalues. The blue circle is an unit circle.

Dynamic Modes and Time Dynamics

In DMD model, Psi refers to the dynamic modes. We can visualize it as in Figure 6. As can be seen, dynamic modes shows the spatial patterns of the example data, and there are significant changes in some loop detectors like from 176 to 196. This is also consistent with the traffic congestion (with low traffic speed) monitored by these loop detectors.

Figure 6. Heatmap of dynamic modes of the data in DMD. — **Figure 6**. Heatmap of dynamic modes of the data in DMD.

The Python code for drawing this figure is given by

fig = plt.figure(figsize = (2, 4))
sns.heatmap(Psi.real)
plt.xticks(np.arange(0.5, 5.5), np.arange(1, 6), rotation = 0)
plt.yticks(np.arange(0.5, 75.5, 10), np.arange(166, 241, 10), rotation = 0)
plt.show()

To understand the time evolution of the example data, we can reconstruct data by using dynamic modes as shown in Figure 7. Here we define the DMD reconstruction as follows,

def reconstruct_DMD_system(data, r):
    T = data.shape[1]
    ## Build data matrices
    X1 = data[:, : -1]
    X2 = data[:, 1 :]
    ## Perform singular value decomposition on X1
    u, s, v = np.linalg.svd(X1, full_matrices = False)
    ## Compute the Koopman matrix
    A_tilde = u[:, : r].conj().T @ X2 @ v[: r, :].conj().T * np.reciprocal(s[: r])
    ## Perform eigenvalue decomposition on A_tilde
    eigval, eigvec = np.linalg.eig(A_tilde)
    ## Compute the coefficient matrix
    Psi = X2 @ v[: r, :].conj().T @ np.diag(np.reciprocal(s[: r])) @ eigvec
    time_dynamics = np.zeros((r, T), dtype = 'complex')
    b = np.linalg.pinv(Psi) @ data[:, 0]
    for t in range(T):
        time_dynamics[:, t] = np.power(eigval, t + 1) * b
    return (Psi @ time_dynamics).real, time_dynamics.real

Figure 7. Reconstructed traffic speed data by using dynamic modes. — **Figure 7**. Reconstructed traffic speed data by using dynamic modes.

To draw Figure 7, please use the following Python code:

r = 5
rec_system, time_dynamics = reconstruct_DMD_system(mat, r)

fig = plt.figure(figsize = (4, 4))
sns.heatmap(rec_system, cmap = 'YlGnBu', cbar_kws={'label': 'Traffic speed'}, vmin = 0, vmax = 75)
plt.xticks(np.arange(0, 72 + 1, 12), np.arange(4 * 288 + 72, 4 * 288 + 144 + 1, 12), rotation = 0)
plt.yticks(np.arange(0.5, 75.5, 10), np.arange(166, 241, 10), rotation = 0)
plt.show()

Here, time dynamics corresponding to each dynamic mode is given by Figure 8.

Figure 8. Time dynamics of 5 dynamic modes in DMD reconstruction. — **Figure 8**. Time dynamics of 5 dynamic modes in DMD reconstruction.

To draw Figure 8, please use the following Python code:

fig = plt.figure(figsize = (10, 6))
for loc in range(1, 6):
    ax = fig.add_subplot(2, 3, loc)
    plt.plot(time_dynamics[loc - 1, :], color = "#006ea3", linewidth = 2.5)
    plt.xticks(np.arange(0, 72 + 1, 12), 
               np.arange(4 * 288 + 72, 4 * 288 + 144 + 1, 12))
    if loc == 1 or loc == 2:
        plt.setp(ax.get_xticklabels(), visible = False)
    plt.grid(axis = 'both', linestyle = "--", linewidth = 0.1, color = 'gray')
    plt.title('Time dynamics of mode {}'.format(loc))
    ax.tick_params(direction = "in")
    ax.set_xlim([0, 72])
plt.show()

Conclusion

In this blog post, we introduce a toy example of DMD model to the application of spatiotemporal traffic data analysis. DMD should have many interesting applications in spatiotemporal data analysis due to its meaningful interpretation. If you have any suggestion, please let me know! Thank you for your reading this post!