SINDy with Control

Introduction

The SINDy method is useful for fitting governing equations to data drawn from a dynamical system. However, for engineering applications we often seek not just to analyze dynamical systems but to control them. In this short post I implement the SINDyC method¹, which generalizes SINDy to include external inputs and feedback control. This continues our investigations from the last post in this series.

Background

In Dynamic Mode Decomposition, we tried to fit a linear operator \(A\) to a dynamical system such that

\[ x_{t+1} = Ax_t \]

DMDc extends this to instead fit the equation:

\[ x_{t+1} = Ax_t + Bu_t \]

Instead of just a matrix of snapshots, we also need a matrix of control history (the \(u_k\) at each timestamp).

The Koopman analysis for this system also now includes \(u\). Instead of

\[ Kg(x_t) := g(F(x_t)) = g(x_{t+1}) \]

We have

\[ K_*g(x_t, u_t) := g(F(x_t, u_t), *) = g(x_{t+1}, *) \]

Note that \(K_*\) depends on the choice of control vector².

Setup

We start with the same dataset of snapshots \(X\), but now we also record our control history at each time point:

\[ \Upsilon = [u_1, u_2, ..., u_m] \]

As in SINDy, we create a dictionary of basis functions. Then, using the data, we fit a sparse set of coefficients to the basis functions:

\[ \dot{X} = \mathbf{\Theta}(X, \Upsilon)\mathbf{\Xi} \]

However (and this is critical), if the signal \(u\) corresponds to a feedback control signal, we cannot disambiguate the effect of the feedback control from that of the internal system. That is, we must intervene on the controls³ to separate them from the dynamics.

Implementation

Coding this is almost trivially similar to SINDy, with a few modifications.

Library

We’ll need to upgrade our library of functions to include cross-terms between \(x\) and \(u\)

def sindyc_library(X, U, poly_order_x=3, poly_order_u=1,
                   include_cross=True, include_sine=False, include_cosine=False):
    n_vars, n_samples = X.shape
    n_ctrl = U.shape[0]

    feats = [np.ones(n_samples)]
    descriptions = ['1']

    for order in range(1, poly_order_x + 1):
        for combo in combinations_with_replacement(range(n_vars), order):
            term = np.prod([X[i, :] for i in combo], axis=0)
            feats.append(term)
            descriptions.append('*'.join([f'x_{i}' for i in combo]))

    for order in range(1, poly_order_u + 1):
        for combo in combinations_with_replacement(range(n_ctrl), order):
            term = np.prod([U[i, :] for i in combo], axis=0)
            feats.append(term)
            descriptions.append('*'.join([f'u_{i}' for i in combo]))

    if include_cross:
        for i in range(n_vars):
            for j in range(n_ctrl):
                feats.append(X[i, :] * U[j, :])
                descriptions.append(f'x_{i}*u_{j}')

    if include_sine:
        for i in range(n_vars):
            feats.append(np.sin(X[i, :])); descriptions.append(f'sin(x_{i})')
    if include_cosine:
        for i in range(n_vars):
            feats.append(np.cos(X[i, :])); descriptions.append(f'cos(x_{i})')

    Theta = np.column_stack(feats)
    return Theta, descriptions

Algorithm

The actual sindyc algorithm is more or less the same as the sindy method:

def sindyc(X, U, dt, poly_order_x=3, poly_order_u=1,
           include_cross=True, lambda_reg=0.1, max_iter=10,
           include_sine=False, include_cosine=False):
    dXdt = finite_difference(X, dt) 
    Theta, descriptions = sindyc_library(
        X, U,
        poly_order_x=poly_order_x,
        poly_order_u=poly_order_u,
        include_cross=include_cross,
        include_sine=include_sine,
        include_cosine=include_cosine
    )
    Xi = sequential_threshold_least_squares(Theta, dXdt, lambda_reg=lambda_reg, max_iter=max_iter)
    return Xi, descriptions

We first use finite differences to get the derivatives of X, then we use sequential threshold least-squares to compute the actual coefficients.

Example

Let’s take a look at a simple example⁴: a predator-prey system with a sinusoidal forcing function:

\[ \begin{aligned} \dot x(t) &= \alpha\,x(t) - \beta\,x(t)\,y(t) \;+\; k_{1}\,u_{1}(t),\\ \dot y(t) &= \delta\,x(t)\,y(t) - \gamma\,y(t) \;-\; k_{2}\,u_{2}(t). \end{aligned} \]

where

\[ \begin{aligned} u_{1}(t) \;&=\; 0.3\,\sin(0.3\,t) \;+\; 0.2\,\cos(0.11\,t) \\ \qquad u_{2}(t) \;&=\; 0.25\,\sin(0.17\,t + 0.7) \end{aligned} \]

We’ll pick values for the coefficients as such:

\[ \begin{aligned} \dot x(t) &= 1.0\,x(t)\;-\;0.5\,x(t)\,y(t)\;+\;0.8\,u_{1}(t)\\ \dot y(t) &= 0.5\,x(t)\,y(t)\;-\;1.0\,y(t)\;-\;0.6\,u_{2}(t) \end{aligned} \]

Data Generation

We need to generate the data. Let’s write up the code for the predator-prey system and for our controls:

def predator_prey_control_rhs(state, u, alpha=1.0, beta=0.5, delta=0.5, gamma=1.0, k1=0.8, k2=0.6):
    x, y = state
    u1, u2 = u
    dx = alpha*x - beta*x*y + k1*u1
    dy = delta*x*y - gamma*y - k2*u2
    return np.array([dx, dy])

Above is the predator-prey system. Here’s what our control functions will look like:

u1 = lambda t: 0.3*np.sin(0.3*t) + 0.2*np.cos(0.11*t)
u2 = lambda t: 0.25*np.sin(0.17*t + 0.7)

Now we can simulate the system:

# Helper function, takes a step of fourth order Runge-Kutta
# See any numerical methods book, like https://link.springer.com/book/10.1007/978-3-540-78862-1
def rk4_step(ode_func_f_X, y_current, dt, **ode_kwargs): 
    k1 = ode_func_f_X(y_current, **ode_kwargs)
    k2 = ode_func_f_X(y_current + dt/2 * k1, **ode_kwargs)
    k3 = ode_func_f_X(y_current + dt/2 * k2, **ode_kwargs)
    k4 = ode_func_f_X(y_current + dt * k3, **ode_kwargs)
    y_next = y_current + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
    return y_next

def simulate_predator_prey_with_control(x0, t, u1, u2, rhs=predator_prey_control_rhs):
    n  = len(t)
    dt = t[1] - t[0]
    X  = np.zeros((2, n))
    U  = np.zeros((2, n))
    X[:, 0] = np.asarray(x0, dtype=float)

    for k in range(1, n):
        u_vec = np.array([u1(t[k-1]), u2(t[k-1])], dtype=float)
        U[:, k-1] = u_vec
        X[:, k] = rk4_step(rhs, X[:, k-1], dt, u=u_vec)

    U[:, -1] = np.array([u1(t[-1]), u2(t[-1])], dtype=float)
    return X, U, dt

Outcome

Putting the code together, we get:

if __name__ == "__main__":
    t = np.arange(0.0, 50.0, 0.01) 

    u1 = lambda _t: 0.3*np.sin(0.3*_t) + 0.2*np.cos(0.11*_t)
    u2 = lambda _t: 0.25*np.sin(0.17*_t + 0.7)

    Xpp, Upp, dt_pp = simulate_predator_prey_with_control(
        x0=(1.5, 1.0),
        t=t,
        u1=u1,
        u2=u2,
        rhs=predator_prey_control_rhs
    )

    Xi_c, desc_c = sindyc(
        Xpp, Upp, dt_pp,
        poly_order_x=2,
        poly_order_u=1,
        include_cross=True,
        lambda_reg=0.05,
        max_iter=15
    )

    print_equations(Xi_c, desc_c, feature_names=['x','y'])

And the outcome:

dx/dt = 0.999839*x_0 - 0.499924*x_0*x_1 + 0.799839*u_0
dy/dt = -0.999844*x_1 + 0.499927*x_0*x_1 - 0.599914*u_1

Which matches our expected coefficients closely.

Conclusion

This post was a straightforward extension of SINDy to accommodate controls.

Footnotes

1. See this paper by Brunton, Proctor, and Kutz. I may use slightly different notation.

2. The \(*\) here is a placeholder indexing the family of possible controls. That is, \(\forall u \in U, (K_{u} g)(x_t) = g(F(x_t,u))\)

3. “Persistent excitation” is required (I don’t fully understand this requirement yet but it’s a requirement for identifiability). Brunton, Proctor, and Kutz recommend injecting a sufficiently large white noise signal, or occasionally kicking the system with a large impulse or step in.

4. Also drawn from Brunton, Proctor, and Kutz.