Defining source term and constraint throws "cannot reshape array of size 12 into shape (1,10)"

[ZorrRRo]

Is there a way to define a source term "u" within the ConstrainedSR3 optimizer?

I have a 2nd order ODE with excitation (Duffing osc). This 2nd order ODE is built of a system of 2 first order ode:

def duffing_oscillator(t, Y, args):
x0 = Y[0]
x1 = Y[1]
return [x1,
-x0 - 0.6x0**3 - 0.2*x1 + f(t) ]

Using the ps.STLSQ() the output is matching perfectly "including the source term":
x0' = 1.000 x1
x1' = -1.000 x0 + -0.200 x1 + 1.000 u0 + -0.600 x0^3

If I switch to ConstrainedSR3() and define the parameter of x1 in the first equaiton to be 1:

feature_names = ["x0", "x1"]
library = ps.PolynomialLibrary()
library.fit( x_train )
print( library.get_feature_names(input_features = feature_names) )
n_features = library.n_output_features_
n_targets  = x_train.shape[1]

constraint_rhs = np.array([1])

constraint_lhs = np.zeros((1, n_targets*n_features))
constraint_lhs[0, 0] = 1

opt = ps.ConstrainedSR3( constraint_rhs = constraint_rhs,
                         constraint_lhs = constraint_lhs, 
                         threshold      = 0.005, 
                         thresholder    = 'l1')

model = ps.SINDy( optimizer       = opt, 
                  feature_library = library, 
                  feature_names   = feature_names)

model.fit(x_train, u = u_train, t = dt)
model.print()

# Simulate the sindy model and compare it to the original solution
x_test_sim = model.simulate(x0_train, tsim, u=f)

This produces an error: cannot reshape array of size 12 into shape (1,10)

If I use the fit method without u: model.fit(x_train, t = dt) , then I get:

x0' = 0.995 1 + 0.998 x1 + -0.369 x0^2 + -0.273 x1^2
x1' = 0.057 1 + -1.559 x0 + 0.017 x1 + -0.019 x0^2 + -0.002 x0 x1 + -0.012 x1^2
UserWarning: Control variables u were ignored because control variables were not used when the model was fit

Any suggestions or ideas? Thx and cheers

[akaptano]

I believe your constraint is not what you intend. It is prescribing that the coefficient of the constant term in the x0' equation should be equal to 1, which is why you get x0' ~ 1 + ..., which is not the duffing oscillator. Try constraint_lhs[0, 2] = 1 to make x0' = 1 * x1 + ... and you can additionally add constraints to zero out all the other terms (since x0' = x1 exactly here actually).

I think the control u should be compatible with ConstrainedSR3 since u is just appended to the existing state variables so x = [x0, x1, u] and the constraints should work (albeit may need to be adjusted in the indices). What is the variable f above that you use during model.simulate?

[ZorrRRo]

Thx a lot for your answer. Now I understand the constraints. It works, thx.

I figured out, that library.n_output_features_ returned me 6 for degree=2 of polynomial library. Now, with degree=3 I need to manually define the size to be 19. Any idea what did go wrong?

Just in case somebody comes across the same issue, here the full example code:

import numpy as np
from scipy.integrate import solve_ivp
import pysindy as ps
from matplotlib import pyplot as plt


def duffing_oscillator(t, Y, *args):
    
    x0 = Y[0]
    x1 = Y[1]

    return [ 1*x1, 
             -1*x0 - 0.6*x0**3 - 0.2*x1 + 1*f(t)]


def f(t):
    return 1*np.cos(1.2*t)


if __name__ == "__main__":
    
    integrator_keywords               = {}
    integrator_keywords["first_step"] = 1e-10
    integrator_keywords["rtol"]       = 1e-9
    integrator_keywords["atol"]       = 1e-12
    integrator_keywords["method"]     = "LSODA"
    
    tend    = 50
    dt      = tend / 5_000
    t_train = np.arange(0, tend, dt)

    x0_train = [0, 0]
    z = solve_ivp( duffing_oscillator, 
                   [t_train[0], t_train[-1]+dt], 
                   x0_train, 
                   t_eval = t_train, 
                   args   = (f,),
                   **integrator_keywords )
    
    tsim    = z.t
    x_train = z.y.T
    
    u_train = f(tsim)
    x0      = x_train[:,0]
    x1      = x_train[:,1]
    
    fig, ax = plt.subplots(nrows=3, sharex=True)
    ax[0].plot(tsim, u_train, c="C0", label="f")
    ax[1].plot(tsim, x0,      c="C1", label="x0")
    ax[2].plot(tsim, x1,      c="C2", label="x1")
    
    for a in ax:
        a.grid()
        a.legend()
        
    #--- sindy part ---
    feature_names = ["x0", "x1", "f"]
    library = ps.PolynomialLibrary(degree = 3, include_bias = False)

    library.fit( x_train )

    n_features = 19                  # library.n_output_features_ does not work
    n_targets  = x_train.shape[1]    # stands for 2 equations x0" and x1"

    constraint_rhs = np.array([1])   # number of constraints, here only 1 
    
    constraint_lhs = np.zeros((1, n_targets*n_features))
    
    # in equation with index 0, all coeff are 0 (already set by np.zeros in previous line)
    # and only coefficient for x1 set to be 1
    # Remember, include_bias=False --> one coefficient less
    constraint_lhs[0, 1] = 1
    
    opt = ps.ConstrainedSR3( constraint_rhs = constraint_rhs,
                             constraint_lhs = constraint_lhs, 
                             threshold      = 0.001, 
                             thresholder    = "l1",
                             max_iter       = 100)
    
    model = ps.SINDy( optimizer       = opt, 
                      feature_library = library, 
                      feature_names   = feature_names)
    
    model.fit(x_train, t = dt, u=u_train)   # u_train not working
    model.print()
    
    # Simulate the sindy model and compare it to the original solution
    integrator_keywords = {}
    integrator_keywords["rtol"] = 1e-3
    integrator_keywords["h0"]   = 1e-3
    integrator_keywords["hmax"] = 0.1
    
    x_test_sim = model.simulate(x0_train, 
                                tsim, 
                                u=f,
                                integrator_kws = integrator_keywords)
    
    ax[1].plot(t_train, x_test_sim[:,0], "k--", label="x0 sindy", linewidth=0.8)
    ax[2].plot(t_train, x_test_sim[:,1], "k--", label="x1 sindy", linewidth=0.8)

[ZorrRRo]

Wait, just to clear out some doubts on my side: You mentioned that "and you can additionally add constraints to zero out all the other terms (since x0' = x1 exactly here actually)."
Is this not achieved by constraint_lhs = np.zeros((1, n_targets*n_features)) automatically? Because observing constraint_lhs in the variable explorer it has everywhere else zero entry ???
If not, how can I make sure, that all other coeff are 0 in this first equation explicitly?

[akaptano]

For your two questions:

1. The library of degree 2 is [1, x0, x1, x0x1, x0^2, x1^2] so it is size 6. For degree 3 it is [1, x0, x1, x0x1, x0^2, x1^2, x0x1^2, x1x0^2, x0^3, x1^3] but you have added a control parameter, which gets added to the model as if it is a x2 variable:
   [1, x0, x1, x2, x0x1, x0x2, x1x2, x0^2, x1^2, x2^2, x0x1^2, x0x2^2, x0x1x2, x1x0^2, x1x2^2, x2x0^2, x2x1^2, x0^3, x1^3, x2^3]
   which is size 20 (your library has size 19 because you omitted the '1' term). In general a polynomial for a state space of size k and degree n should have (n+k) choose n terms (https://mathoverflow.net/questions/225953/number-of-polynomial-terms-for-certain-degree-and-certain-number-of-variables) if you are including the constant term.

This library is not the one you were using for computing n_features because you used library.fit( x_train ) without the u variable! x_train becomes [x_train, u] when you call model.fit.

2. Your doubts about the constraint are correct. What are you doing currently is prescribing a single constraint, that says that the x0' equation needs to have a coefficient of 1 on the x1 equation term. But you have not constrained any of the other coefficients that can appear in this equation! The other zeros in constraint_lhs is merely a reflection that you are only constraining the x1 equation term (so if all the constraint_lhs terms were nonzero, it would still be one constraint that relates all of the coefficients together). What you have done is probably sufficient for achieving x0' = x1, but it you want this to be exactly true no matter what, you need 18 other constraints that tells the model to set the other 18/19 library terms to zero.

[ZorrRRo]

1. Ok, I understand, thx. I fixed it by providing the excitation like (copy paste of one of your examples):
   shape = np.array(x_train.shape) shape[-1] = 1 x = np.concatenate([x_train, u_train[0] * np.ones(shape)], axis=-1) library.fit(x)

2. I understand. For each constraint, we need one equation. Done. Thx

Tor the sake of completeness, here the part I refer to as SINDY part with the solution of this discussion (if somebody migth face same issues).
All coefficients are set to 0 but only coefficient for x1 is set to be 1.

    #--- sindy part ---
    feature_names = ["x0", "x1", "f"]
    library = ps.PolynomialLibrary(degree = 5, include_bias = False)
    
    # because of u_train library.fit needs to be executed with a matrix of size(5000, 3)
    # and not only with x_train which has size(5000, 2)
    shape = np.array(x_train.shape)
    shape[-1] = 1
    library.fit(np.concatenate([x_train, u_train[0] * np.ones(shape)], axis=-1))
    names = library.get_feature_names()
    
    n_features = library.n_output_features_  # library.n_output_features_ does not work
    n_targets  = x_train.shape[1]                  # stands for 2 equations x0' and x1'

    constraint_rhs = np.array([0]*n_features) # we will constrain all coefficients to 0 in x0'
    x_1 = names.index('x1')   # But only coefficient for x1 will be different
    constraint_rhs[x_1] = 1   # Coefficient for x1 must be 1
    
    # One row per constraint, one column per coefficient
    # Three equations defining the constraint
    constraint_lhs = np.zeros((n_features, n_targets*n_features))
    
    # This loop activates as many constraints to be active as coefficients are available
    # Equation 1: 1*coeff x0 = 0
    # Equation 2: 1*coeff x1 = 1
    # Equation 3: 1*coeff x. = 0
    # .
    # .
    # Equation n_features: 1*coeff x2^5 = 0 
    for idx in range(n_features):
        constraint_lhs[idx, idx] = 1
    
    opt = ps.ConstrainedSR3( constraint_rhs = constraint_rhs,
                             constraint_lhs = constraint_lhs, 
                             threshold      = 0.001, 
                             thresholder    = 'l1',
                             max_iter       = 100)
    
    model = ps.SINDy( optimizer       = opt, 
                      feature_library = library, 
                      feature_names   = feature_names)
    
    model.fit(x_train, t = dt, u=u_train)   # u_train not working
    model.print()
    
    # Simulate the sindy model and compare it to the original solution
    integrator_keywords = {}
    integrator_keywords['rtol'] = 1e-3
    integrator_keywords['h0']   = 1e-3
    integrator_keywords['hmax'] = 0.1
    
    x_test_sim = model.simulate(x0_train, 
                                tsim, 
                                u=f,
                                integrator_kws = integrator_keywords)
    
    ax[1].plot(t_train, x_test_sim[:,0], 'k--', label='x0 sindy', linewidth=0.8)
    ax[2].plot(t_train, x_test_sim[:,1], 'k--', label='x1 sindy', linewidth=0.8)