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}\]

%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.

EDMDc = pk.regression.DMDc()

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

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())

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.297496697624256')

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)