Hankel DMD with control for Van der Pol Oscillator

Classical forced Van der Pol oscillator this is example in Sec. 4 in Korda & Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control”, Automatica 2018, with dynamics given by:

\[\begin{split}\dot{x}_{1} = 2x_2 u,\\ \dot{x}_2 = -0.8x_1 + 2x_2 -10x_1^2x_2 + u\end{split}\]

import sys
sys.path.append('../src')

%matplotlib inline
import pykoopman as pk
from pykoopman.common.examples import vdp_osc, rk4, square_wave  # required for example system
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as rnd
np.random.seed(42)  # for reproducibility

import warnings
warnings.filterwarnings('ignore')

Training data: A training dataset is created consisting of 200 trajectories, each trajectory is integrated for 1000 timesteps and forced by a random actuation in the range \([-1,1]\). Each trajectory starts at a random initial condition in the unit box \([-1,1]^2\).

n_states = 2 # Number of states
n_inputs = 1 # Number of control inputs
dT = 0.01    # Timestep
n_traj = 200  # Number of trajectories
n_int = 1000  # Integration length

# Time vector
t = np.arange(0, n_int*dT, dT)

# Uniform random distributed forcing in [-1, 1]
u = 2*rnd.random([n_int, n_traj])-1

# Uniform distribution of initial conditions
x = 2*rnd.random([n_states, n_traj])-1

# Init
X = np.zeros((n_states, n_int*n_traj))
Y = np.zeros((n_states, n_int*n_traj))
U = np.zeros((n_inputs, n_int*n_traj))

# Integrate
for step in range(n_int):
    y = rk4(0, x, u[step, :], dT, vdp_osc)
    X[:, (step)*n_traj:(step+1)*n_traj] = x
    Y[:, (step)*n_traj:(step+1)*n_traj] = y
    U[:, (step)*n_traj:(step+1)*n_traj] = u[step, :]
    x = y

# Visualize first 100 steps of the training data
fig, axs = plt.subplots(1, 1, tight_layout=True, figsize=(12, 4))
for traj_idx in range(n_traj):
    x = X[:, traj_idx::n_traj]
    axs.plot(t[0:100], x[1, 0:100], 'k')
axs.set(
        ylabel=r'$x_2$',
        xlabel=r'$t$')

[Text(0, 0.5, '$x_2$'), Text(0.5, 0, '$t$')]

EDMDc model: The observables (or lifting functions) for the Koopman model are chosen to be the state itself (\(\psi_1 = x_1,\psi_2=x_2\)) (by setting include_states=True below, which is also the default) and thin plate radial basis functions with centers selected randomly.

n_delays =  9
obs = pk.observables.TimeDelay(n_delays=n_delays)

EDMDc = pk.regression.DMDc()

# which is in effect
# EDMDc = pk.regression.DMDc(
#     svd_rank=n_states*(n_delays+1)+n_inputs,
#     svd_output_rank=n_states*(n_delays+1)
# )


model = pk.Koopman(observables=obs, regressor=EDMDc)
model.fit(X.T, y=Y.T, u=U[:,n_delays:].T)

Koopman(observables=TimeDelay(n_delays=9), regressor=DMDc())

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Compare prediction accuracy of Koopman model on a test trajectory

In the following, the trained model is used to perform a multi-step prediction from a given initial condition. The prediction is compared with the true trajectory when integrating the nonlinear system.

n_int = 600  # Integration length
t = np.arange(0, n_int*dT, dT)
u = np.array([-square_wave(step+1) for step in range(n_int)])
x = np.array([0.5, 0.5])
# x = np.array([[-0.1], [-0.5]])

# Integrate nonlinear system
Xtrue = np.zeros((n_states, n_int))
Xtrue[:, 0] = x
for step in range(1, n_int, 1):
    y = rk4(0, Xtrue[:, step-1].reshape(n_states,1), u[np.newaxis, step-1], dT, vdp_osc)
    Xtrue[:, step] = y.reshape(n_states,)

Predict using Koopman model

x0_td = Xtrue[:,:n_delays+1].T

# Multi-step prediction with Koopman/EDMDc model
Xkoop = model.simulate(x0_td, u[n_delays:, np.newaxis], n_steps=n_int-1-n_delays)

# Multi-step prediction with Koopman/EDMDc model
Xkoop = np.vstack([x0_td, Xkoop]) # add initial condition to simulated data for comparison below

Compare results

fig, axs = plt.subplots(3, 1, sharex=True, tight_layout=True, figsize=(9, 6))
axs[0].plot(t, u, '-k')
axs[0].set(ylabel=r'$u$')
axs[1].plot(t, Xtrue[0, :], '-', color='b', label='True')
axs[1].plot(t, Xkoop[:, 0], '--r', label='TimeDelayDMDc')
axs[1].set(ylabel=r'$x_1$')
axs[2].plot(t, Xtrue[1, :], '-', color='b', label='True')
axs[2].plot(t, Xkoop[:, 1], '--r', label='TimeDelayDMDc')
axs[2].set(
        ylabel=r'$x_2$',
        xlabel=r'$t$')

for i in range(1,3):
    axs[i].axvline(x=t[n_delays])

axs[1].legend()

err = np.linalg.norm(Xtrue - Xkoop.T)
axs[0].set_title(f"L2 norm error = {err}")

Text(0.5, 1.0, 'L2 norm error = 7.249270040251282')

Lifted system can be easily accessed

\[\begin{split}z^{+} = A z + B u\\ x = Cz\end{split}\]

model.A.shape

(20, 20)

model.B.shape

(20, 1)

model.C.shape

(2, 20)

model.W.shape

(2, 20)