# Extended Trapping SINDy By Mai Peng and Alan Kaptanoglu
A very common issue is that models identified by system identification methods typically have no guarantees that the models are numerically or physically stable. This can be addressed with heuristic, data-driven, or analytic closure models, but we have recently directly promoted globally stable models into the system identification itself (see the Example 8 Jupyter notebook). This is really nice but there are three potential caveats, (1) the regression is nonconvex and there a number of hyperparameters, so this method can be difficult to learn, and (2) in order to promote global stability, one needs an analytic result from stability theory, and the one we use applies only for quadratically nonlinear dynamics (typically fluid and plasma flows) with energy-preserving, quadratic, nonlinearities. Moreover, we have good reason to believe that (3) generic quadratically nonlinear models will always be globally unbounded, so for these situations we can also promote local Lyapunov stability of the origin using some variations of the original Trapping SINDy algorithm. That is the goal of this notebook – to illustrate how various forms of global and local stability can be promoted explicitly in the SINDy method to obtain stable data-driven models.
For the following, we will consider dynamical models of the form
For global stability promotion, we will require that the totally symmetric part of the quadratic coefficients vanishes (without loss of generality, \(Q_{ijk}\) is symmetric in the last two indices):
This equation can be implemented as a hard or soft constraint in the optimization. For dynamical models that do not satisfy this condition, we can still promote locally stable models that are stable even at very large distances of the origin. The following examples show different ways to relax this hard constraint.
[73]:
import warnings
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp
import pysindy as ps
from pysindy.utils import lorenz
# ignore warnings
warnings.filterwarnings("ignore")
# Import useful functions
from trapping_utils import (
integrator_keywords,
sindy_library,
make_fits,
obj_function,
check_local_stability,
make_trap_progress_plots,
)
np.random.seed(10) # for reproducibility
Lorenz model
The Lorenz system originates from a simple fluid model of atmospheric dynamics from Lorenz et al. (1963). This system is likely the most famous example of chaotic, nonlinear behavior despite the somewhat innocuous system of equations,
For Lorenz’s choice of parameters, \(\sigma = 10\), \(\rho = 28\), \(\beta = 8/3\), this system is known to exhibit a stable attractor. For \(\mathbf{m} = [0,m_2,\rho+\sigma]\) (\(m_1\) does not contribute to \(\mathbf{A}^S\) so we set it to zero),
so that if \(\lambda_{\pm} < 0\), then \(-2\sqrt{\sigma\beta} < m_2 < 2\sqrt{\sigma\beta}\). Our algorithm can successfully identify the optimal \(\mathbf{m}\), and can be used to identify the inequality bounds on \(m_2\) for stability.
Check global stability of the Lorenz model
The skew-symmetric models below are globally stable if and only if there exists a vector \(\mathbf{m}\) such that following matrix is negative definite:
Note that if the quadratic tensor has zero totally symmetric part, this is equal to
A negative definite \(\mathbf{A}^S\) turns out to also be necessary for models that do not quite satisfy the constraint on \(Q_{jik}\), but in this case is not sufficient for global boundedness.
A decent-enough algorithm for a nonlinear search for such a \(\mathbf{m}\) that makes \(A^S_{ij}\) negative definite is simulated annealing, and a simple interface is provided by scipy.
[74]:
# define parameters
r = 3
N = int((r**2 + 3 * r) / 2.0) + 1
# make training and testing data
dt = 0.01
T = 40
t = np.arange(0, T + dt, dt)
t_span = (t[0], t[-1])
x0 = [1, -1, 20]
x_train = solve_ivp(lorenz, t_span, x0, t_eval=t, **integrator_keywords).y.T
x0 = (np.random.rand(3) - 0.5) * 30
x_test = solve_ivp(lorenz, t_span, x0, t_eval=t, **integrator_keywords).y.T
# define hyperparameters
reg_weight_lam = 0
max_iter = 5000
eta = 1.0e3
alpha_m = 8e-1 * eta
# run trapping SINDy
sindy_opt = ps.TrappingSR3(
method="global",
_n_tgts=3,
_include_bias=True,
reg_weight_lam=reg_weight_lam,
eta=eta,
max_iter=max_iter,
gamma=-1,
verbose=True,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library,
)
model.fit(x_train, t=t)
model.print()
# Extract model coefficients and check how well constraint is satisfied
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_unsym_
PQ_tensor = sindy_opt.PQ_
Lenergy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Q = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
print(
r"|tilde{H_0}|_F = ",
np.sqrt(np.sum((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))**2)),
)
Iter ... |y-Xw|^2 ... |Pw-A|^2/eta ... |w|_1 ... |Qijk|/a ... |Qijk+...|/b ... Total:
0 ... 2.500e+02 ... 1.273e-01 ... 0.00e+00 ... 4.93e-21 ... 2.88e-47 ... 2.50e+02
500 ... 2.500e+02 ... 4.031e-02 ... 0.00e+00 ... 4.93e-21 ... 8.86e-48 ... 2.50e+02
1000 ... 2.500e+02 ... 1.419e-02 ... 0.00e+00 ... 4.93e-21 ... 2.94e-47 ... 2.50e+02
1500 ... 2.500e+02 ... 5.726e-03 ... 0.00e+00 ... 4.93e-21 ... 3.08e-47 ... 2.50e+02
2000 ... 2.500e+02 ... 2.669e-03 ... 0.00e+00 ... 4.93e-21 ... 1.05e-48 ... 2.50e+02
2500 ... 2.500e+02 ... 1.415e-03 ... 0.00e+00 ... 4.93e-21 ... 1.52e-47 ... 2.50e+02
3000 ... 2.500e+02 ... 8.329e-04 ... 0.00e+00 ... 4.93e-21 ... 3.58e-47 ... 2.50e+02
3500 ... 2.500e+02 ... 5.323e-04 ... 0.00e+00 ... 4.93e-21 ... 2.94e-47 ... 2.50e+02
4000 ... 2.500e+02 ... 3.628e-04 ... 0.00e+00 ... 4.93e-21 ... 3.46e-47 ... 2.50e+02
4500 ... 2.500e+02 ... 2.600e-04 ... 0.00e+00 ... 4.93e-21 ... 2.11e-47 ... 2.50e+02
(x0)' = 0.060 1 + -9.868 x0 + 9.938 x1 + -0.009 x2 + -0.001 x0 x1 + -0.003 x0 x2 + 0.001 x1^2 + 0.001 x1 x2
(x1)' = -0.244 1 + 27.760 x0 + -0.917 x1 + 0.029 x2 + 0.001 x0^2 + -0.001 x0 x1 + -0.993 x0 x2 + -0.002 x1 x2 + -0.001 x2^2
(x2)' = 0.083 1 + 0.001 x0 + -0.010 x1 + -2.664 x2 + 0.003 x0^2 + 0.993 x0 x1 + 0.002 x1^2 + 0.001 x1 x2
|tilde{H_0}|_F = 7.150665448126821e-14
[75]:
# Calculate the x_dot and x trajectories for train and test sets
xdot_test = model.differentiate(x_test, t=t)
xdot_test_pred = model.predict(x_test)
x_train_pred = model.simulate(x_train[0, :], t, integrator_kws=integrator_keywords)
x_test_pred = model.simulate(x_test[0, :], t, integrator_kws=integrator_keywords)
# plotting and analysis
make_fits(r, t, xdot_test, xdot_test_pred, x_test, x_test_pred, "lorenz")
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
E_pred = np.linalg.norm(x_test - x_test_pred) / np.linalg.norm(x_test)
print("Frobenius error = ", E_pred)
check_local_stability(Xi, sindy_opt, mean_val)
# compute relative Frobenius error in the model coefficients
sigma = 10
rho = 28
beta = 8.0 / 3.0
terms = sindy_library.get_feature_names()
Xi_lorenz = np.zeros(Xi.shape)
Xi_lorenz[1 : r + 1, :] = np.array([[-sigma, sigma, 0], [rho, -1, 0], [0, 0, -beta]]).T
Xi_lorenz[terms.index("x0 x2"), 1] = -1
Xi_lorenz[terms.index("x0 x1"), 2] = 1
coef_pred = np.linalg.norm(Xi_lorenz - Xi) / np.linalg.norm(Xi_lorenz)
print("Frobenius coefficient error = ", coef_pred)
# Compute time-averaged dX/dt error
deriv_error = np.zeros(xdot_test.shape[0])
for i in range(xdot_test.shape[0]):
deriv_error[i] = np.dot(
xdot_test[i, :] - xdot_test_pred[i, :], xdot_test[i, :] - xdot_test_pred[i, :]
) / np.dot(xdot_test[i, :], xdot_test[i, :])
print("Time-averaged derivative error = ", np.nanmean(deriv_error))
Frobenius error = 0.5044638140847056
optimal m: [-1.15265065 -0.14638889 33.16397213]
As eigvals: [-10.55590873 -2.66324565 -0.37704047]
0.5 * |tilde{H}_0|_F = 3.5753327240634105e-14
0.5 * |tilde{H}_0|_F^2 / beta = 2.5566008175517377e-47
Estimate of trapping region size, Rm = 235.692775304696
Normalized trapping region size, Reff = 9.49393574147844
Local stability size, R_ls= 15818407806117.9
Frobenius coefficient error = 0.012548531064521971
Time-averaged derivative error = 1.0438764179317514e-05
Use simulated annealing
We are going to check if any \(\mathbf{m}\) exists such that \(\mathbf{A}^S\) is negative definite, using the identified coefficients, to verify again that our model is globally stable.
[76]:
# Import simulated annealing algorithm from scipy
from scipy.optimize import dual_annealing as anneal_algo
boundvals = np.zeros((r, 2))
boundmax = 1000
boundmin = -1000
boundvals[:, 0] = boundmin
boundvals[:, 1] = boundmax
PL_tensor_unsym = sindy_opt.PL_unsym_
PL_tensor = sindy_opt.PL_
PM_tensor = sindy_opt.PM_
L = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Q = np.tensordot(PM_tensor, Xi, axes=([4, 3], [0, 1]))
# run simulated annealing
algo_sol = anneal_algo(
obj_function,
bounds=boundvals, # obj_function imported from utils.py
args=(L, Q, np.eye(r)),
maxiter=500,
)
opt_m = algo_sol.x
opt_energy = algo_sol.fun
opt_result = algo_sol.message
print("Result:")
print("Optimal m = ", opt_m)
print(
"Algorithm managed to reduce the largest eigenvalue of A^S to eig1 = ",
opt_energy,
"\n",
)
Result:
Optimal m = [1000. 9.27724217 75.76047513]
Algorithm managed to reduce the largest eigenvalue of A^S to eig1 = -1.324926031456374
Promoting locally stable models with estimates of the stability radius
So far, we have promoted globally stable models with trapping SINDy by enforcing the skew-symmetry structure in the nonlinearities as a hard constraint in the optimization problem:
This problem is solved with a convex relaxation of the optimization.
Below, we relax the hard constraint to a soft constraint and instead solve
where \(\|\cdot\|_F\) denotes the Frobenius norm. This allows us to build locally Lyapunov stable models, and adjust the size of the local stability radius by varying \(\epsilon_Q\). A conservative estimate of the local stability is:
And the radius of the trapping region is given by:
In other words, there is a region \(\rho_- < \|\mathbf{a}(t)\| < \rho_+\) such that the energy \(K\) satisfies \(K > 0\) and \(\dot{K} < 0\), so that any trajectory with initial condition \(\|\mathbf{a}_0\| < \rho_+\) will be bounded for all time. This is because it will fall towards the origin until at least it reaches \(\rho_-\), and then it stays in the ball of radius \(\rho_-\) for all time.
A better way to optimize
However, we find empirically that CVXPY struggles to solve the inequality-constrained problem adequately, and find much better performance by incorporating the constraint as a loss term in the objective. Two other loss terms that can be used as alternatives to increase the size of the stability radius while avoiding extra constraints:
and
We can combine all of these options into the following unconstrained optimization problem:
We now solve this problem for \(\alpha \gg \beta\), \(\alpha \ll \beta\), and \(\alpha \sim \beta \sim 1.\)
First case: \(\alpha \gg 1\), \(\beta \ll 1\), for which the model should just zero out all the quadratic nonlinear terms
[77]:
max_iter = 500
eta = 1.0e2
alpha = 1e-15
beta = 1e20
reg_weight_lam = 0
# run trapping SINDy... no more constraints!
sindy_opt = ps.TrappingSR3(
method="local",
_n_tgts=3,
_include_bias=True,
reg_weight_lam=reg_weight_lam,
eta=eta,
max_iter=max_iter,
gamma=-1,
alpha=alpha,
beta=beta,
verbose=True,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library,
)
model.fit(x_train, t=t)
model.print()
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_
PQ_tensor = sindy_opt.PQ_
Lenergy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Qenergy = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
Rm, R_ls = check_local_stability(Xi, sindy_opt, mean_val)
Iter ... |y-Xw|^2 ... |Pw-A|^2/eta ... |w|_1 ... |Qijk|/a ... |Qijk+...|/b ... Total:
0 ... 2.352e+07 ... 2.091e-02 ... 0.00e+00 ... 2.45e+00 ... 1.10e-35 ... 2.35e+07
(x0)' = -0.007 1 + 0.749 x0 + 0.887 x1 + -0.192 x2
(x1)' = 0.004 1 + -0.573 x0 + -0.678 x1 + 0.101 x2
(x2)' = -0.001 1 + 0.170 x0 + 0.201 x1 + -0.038 x2
optimal m: [-1.16595599 -0.55935101 -1.99977125]
As eigvals: [-0.72801307 -0.00498714 0.76580054]
0.5 * |tilde{H}_0|_F = 2.3475120683405053e-08
0.5 * |tilde{H}_0|_F^2 / beta = 1.1021625822008634e-35
Estimate of trapping region size, Rm = 0
Normalized trapping region size, Reff = 0.0
Local stability size, R_ls= 0
Indeed, we found that if \(\alpha \gg 1\) large enough, the quadratic terms in the model are zeroed, which is bad news both for fitting the model and for applying the trapping theorem since the theorem relies on nontrivial quadratic contributions.
Second case: \(\alpha \ll 1\), \(\beta \gg 1\), which should reproduce the energy-preserving nonlinear constraint to high accuracy
This is a different strategy for stability – don’t make the model’s quadratic nonlinearities weak, but make it so that the totally symmetric part of \(Q_{ijk}\) is very small.
[78]:
max_iter = 10000
eta = 1.0e3
alpha = 1e20
beta = 1e-10
reg_weight_lam = 0
alpha_m = 0.9 * eta
# run trapping SINDy... no more constraints!
sindy_opt = ps.TrappingSR3(
method="local",
_n_tgts=3,
_include_bias=True,
reg_weight_lam=reg_weight_lam,
eta=eta,
alpha_m=alpha_m,
max_iter=max_iter,
gamma=-1,
alpha=alpha,
beta=beta,
verbose=True,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library,
)
model.fit(x_train, t=t)
model.print()
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_
PQ_tensor = sindy_opt.PQ_
Lenergy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Qenergy = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
R_m, R_ls = check_local_stability(Xi, sindy_opt, mean_val)
Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
Iter ... |y-Xw|^2 ... |Pw-A|^2/eta ... |w|_1 ... |Qijk|/a ... |Qijk+...|/b ... Total:
0 ... 2.499e+02 ... 1.273e-01 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
1000 ... 2.499e+02 ... 6.464e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
2000 ... 2.499e+02 ... 2.864e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
3000 ... 2.499e+02 ... 2.129e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
4000 ... 2.499e+02 ... 1.873e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
5000 ... 2.499e+02 ... 1.767e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
6000 ... 2.499e+02 ... 1.719e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
7000 ... 2.499e+02 ... 1.696e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
8000 ... 2.499e+02 ... 1.685e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
9000 ... 2.499e+02 ... 1.679e-07 ... 0.00e+00 ... 4.93e-21 ... 1.34e-02 ... 2.50e+02
(x0)' = 0.060 1 + -9.868 x0 + 9.938 x1 + -0.009 x2 + -0.001 x0 x1 + -0.003 x0 x2 + 0.001 x1^2 + 0.001 x1 x2
(x1)' = -0.244 1 + 27.760 x0 + -0.917 x1 + 0.029 x2 + 0.001 x0^2 + -0.001 x0 x1 + -0.993 x0 x2 + -0.002 x1 x2 + -0.001 x2^2
(x2)' = 0.083 1 + 0.001 x0 + -0.010 x1 + -2.664 x2 + 0.003 x0^2 + 0.993 x0 x1 + 0.002 x1^2 + 0.001 x1 x2
optimal m: [-1.07004356 -0.11790946 37.99091215]
As eigvals: [-9.97295505 -2.66324761 -0.98168966]
0.5 * |tilde{H}_0|_F = 8.192157383826614e-07
0.5 * |tilde{H}_0|_F^2 / beta = 0.01342228852027698
Estimate of trapping region size, Rm = 103.904169760519
Normalized trapping region size, Reff = 4.18536168409390
Local stability size, R_ls= 1797389.00920495
Plot how the two stability radii changes as the algorithm iterates
As the algorithm iterates, it is biasing the model to have a negative definite \(\mathbf{A}^S\) matrix. Once this is true, we can estimate the local Lyapunov stability radius \(\rho_+\) and the trapping region radius \(\rho_-\).
#### Note that with the soft constraint we can get the stability radius arbitrarily large here!
[79]:
rhos_minus, rhos_plus = make_trap_progress_plots(r, sindy_opt)
Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
Q_sum = np.max(np.abs((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))))
print(
r"|tilde{H_0}|_F = ",
np.sqrt(np.sum((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))**2)),
)
|tilde{H_0}|_F = 1.6384314767653227e-06
[80]:
# Calculate the x_dot and x trajectories for train and test sets
xdot_test = model.differentiate(x_test, t=t)
xdot_test_pred = model.predict(x_test)
x_train_pred = model.simulate(x_train[0, :], t, integrator_kws=integrator_keywords)
x_test_pred = model.simulate(x_test[0, :], t, integrator_kws=integrator_keywords)
x_stability_check = model.simulate(x0, t, integrator_kws=integrator_keywords)
# plotting and analysis
make_fits(r, t, xdot_test, xdot_test_pred, x_test, x_test_pred, "lorenz")
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
E_pred = np.linalg.norm(x_test - x_test_pred) / np.linalg.norm(x_test)
print("Frobenius error = ", E_pred)
check_local_stability(Xi, sindy_opt, mean_val)
# compute relative Frobenius error in the model coefficients
coef_pred = np.linalg.norm(Xi_lorenz - Xi) / np.linalg.norm(Xi_lorenz)
print("Frobenius coefficient error = ", coef_pred)
# Compute time-averaged dX/dt error
deriv_error = np.zeros(xdot_test.shape[0])
for i in range(xdot_test.shape[0]):
deriv_error[i] = np.dot(
xdot_test[i, :] - xdot_test_pred[i, :], xdot_test[i, :] - xdot_test_pred[i, :]
) / np.dot(xdot_test[i, :], xdot_test[i, :])
print("Time-averaged derivative error = ", np.nanmean(deriv_error))
Frobenius error = 0.48911628634556376
optimal m: [-1.07004356 -0.11790946 37.99091215]
As eigvals: [-9.97295505 -2.66324761 -0.98168966]
0.5 * |tilde{H}_0|_F = 8.192157383826614e-07
0.5 * |tilde{H}_0|_F^2 / beta = 0.01342228852027698
Estimate of trapping region size, Rm = 103.904169760519
Normalized trapping region size, Reff = 4.13381260255514
Local stability size, R_ls= 1797389.00920495
Frobenius coefficient error = 0.012549973906050805
Time-averaged derivative error = 1.043894380050916e-05
Repeat \(\alpha \gg 1\), \(\beta \ll 1\) case with \(\lambda > 0\)
I find that solver will fail if eps_solver parameter is made too small (error tolerance of the CVXPY solver is very stringent)
[81]:
max_iter = 100
eta = 1.0e5
alpha = 1e20
beta = 1e-10
reg_weight_lam = 5
alpha_m = 0.9 * eta
# run trapping SINDy... no more constraints!
sindy_opt = ps.TrappingSR3(
method="local",
_n_tgts=3,
_include_bias=True,
reg_weight_lam=reg_weight_lam,
eta=eta,
alpha_m=alpha_m,
max_iter=max_iter,
gamma=-1,
alpha=alpha,
beta=beta,
verbose=True,
eps_solver=1e-3,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library,
)
model.fit(x_train, t=t)
model.print()
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_
PQ_tensor = sindy_opt.PQ_
Lenergy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Qenergy = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
check_local_stability(Xi, sindy_opt, mean_val)
Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
print(
r"|tilde{H_0}|_F = ",
np.sqrt(np.sum((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))**2)),
)
# make_progress_plots(r, sindy_opt)
Iter ... |y-Xw|^2 ... |Pw-A|^2/eta ... |w|_1 ... |Qijk|/a ... |Qijk+...|/b ... Total:
0 ... 2.490e+02 ... 1.272e-03 ... 2.67e+02 ... 4.93e-21 ... 2.15e+02 ... 5.16e+02
10 ... 2.489e+02 ... 1.426e-04 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
20 ... 2.489e+02 ... 3.444e-05 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
30 ... 2.489e+02 ... 1.224e-05 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
40 ... 2.489e+02 ... 5.698e-06 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
50 ... 2.489e+02 ... 3.163e-06 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
60 ... 2.489e+02 ... 1.974e-06 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
70 ... 2.489e+02 ... 1.336e-06 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
80 ... 2.489e+02 ... 9.594e-07 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
90 ... 2.489e+02 ... 7.204e-07 ... 2.67e+02 ... 4.93e-21 ... 2.18e+02 ... 5.16e+02
(x0)' = 0.021 1 + -9.868 x0 + 9.937 x1 + -0.006 x2 + -0.001 x0 x1 + -0.003 x0 x2 + 0.001 x1^2 + 0.001 x1 x2
(x1)' = -0.209 1 + 27.756 x0 + -0.914 x1 + 0.026 x2 + 0.001 x0^2 + -0.001 x0 x1 + -0.993 x0 x2 + -0.002 x1 x2 + -0.001 x2^2
(x2)' = 0.079 1 + -0.002 x0 + -0.008 x1 + -2.664 x2 + 0.003 x0^2 + 0.992 x0 x1 + 0.002 x1^2 + 0.001 x1 x2
optimal m: [-1.14049146 -0.12493646 34.55766159]
As eigvals: [-10.27409431 -2.66301234 -0.66533477]
0.5 * |tilde{H}_0|_F = 0.00010445284825151528
0.5 * |tilde{H}_0|_F^2 / beta = 218.20795015708154
Estimate of trapping region size, Rm = 141.401316013499
Normalized trapping region size, Reff = 5.62563122829160
Local stability size, R_ls= 9413.16972800589
|tilde{H_0}|_F = 0.00020890569650303055
[82]:
rhos_minus, rhos_plus = make_trap_progress_plots(r, sindy_opt)
plt.yscale("log")
plt.ylim(1, rhos_plus[-1] * 1.2)
Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
Q_sum = np.max(np.abs((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))))
print(
"Maximum deviation from having zero totally symmetric part: ", np.max(np.abs(Q_sum))
)
Maximum deviation from having zero totally symmetric part: 0.0001557070395381613
Now we add a lot of noise to the Lorenz data and see if trapping extended algorithm improves robustness to noise.
[83]:
np.random.seed(10)
lorenz_noise = np.random.normal(
0, mean_val / 4, x_train.shape
) # 25% noise added with zero mean
x_train_noise = x_train + lorenz_noise
max_iter = 10000
eta = 1.0e2
alpha = 1e20
beta = 1e-14
reg_weight_lam = 0
alpha_m = 0.1 * eta
# run trapping SINDy... no more constraints!
sindy_opt = ps.TrappingSR3(
method="local",
_n_tgts=3,
_include_bias=True,
reg_weight_lam=reg_weight_lam,
eta=eta,
alpha_m=alpha_m,
max_iter=max_iter,
gamma=-1,
alpha=alpha,
beta=beta,
verbose=True,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library,
)
model.fit(x_train_noise, t=t)
model.print()
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_
PQ_tensor = sindy_opt.PQ_
Lenergy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Qenergy = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
check_local_stability(Xi, sindy_opt, mean_val)
Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
print(
r"|tilde{H_0}|_F = ",
np.sqrt(np.sum((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))**2)),
)
# make_trap_progress_plots(r, sindy_opt)
# Calculate the x_dot and x trajectories for train and test sets
xdot_test = model.differentiate(x_test, t=t)
xdot_test_pred = model.predict(x_test)
x_train_pred = model.simulate(x_train[0, :], t, integrator_kws=integrator_keywords)
x_test_pred = model.simulate(x_test[0, :], t, integrator_kws=integrator_keywords)
# compute relative Frobenius error in the model coefficients
coef_pred = np.linalg.norm(Xi_lorenz - Xi) / np.linalg.norm(Xi_lorenz)
print("Frobenius coefficient error = ", coef_pred)
# Compute time-averaged dX/dt error
deriv_error = np.zeros(xdot_test.shape[0])
for i in range(xdot_test.shape[0]):
deriv_error[i] = np.dot(
xdot_test[i, :] - xdot_test_pred[i, :], xdot_test[i, :] - xdot_test_pred[i, :]
) / np.dot(xdot_test[i, :], xdot_test[i, :])
print("Time-averaged derivative error = ", np.nanmean(deriv_error))
Iter ... |y-Xw|^2 ... |Pw-A|^2/eta ... |w|_1 ... |Qijk|/a ... |Qijk+...|/b ... Total:
0 ... 1.186e+09 ... 3.421e-01 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
1000 ... 1.186e+09 ... 1.567e-02 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
2000 ... 1.186e+09 ... 1.300e-02 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
3000 ... 1.186e+09 ... 1.204e-02 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
4000 ... 1.186e+09 ... 1.140e-02 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
5000 ... 1.186e+09 ... 1.087e-02 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
6000 ... 1.186e+09 ... 1.040e-02 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
7000 ... 1.186e+09 ... 9.967e-03 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
8000 ... 1.186e+09 ... 9.569e-03 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
9000 ... 1.186e+09 ... 9.198e-03 ... 0.00e+00 ... 5.33e-22 ... 3.94e-02 ... 1.19e+09
(x0)' = -0.065 1 + 2.249 x0 + 3.220 x1 + -0.604 x2 + 0.004 x0 x1 + -0.117 x0 x2 + -0.013 x1^2 + -0.034 x1 x2 + 0.016 x2^2
(x1)' = 0.044 1 + 3.268 x0 + 4.355 x1 + 0.201 x2 + -0.004 x0^2 + 0.013 x0 x1 + -0.225 x0 x2 + -0.131 x1 x2 + 0.001 x2^2
(x2)' = -0.013 1 + 0.774 x0 + -0.912 x1 + -1.613 x2 + 0.117 x0^2 + 0.260 x0 x1 + -0.016 x0 x2 + 0.131 x1^2 + -0.001 x1 x2
optimal m: [-12.89183874 -2.63561767 39.47137739]
As eigvals: [-3.99906992 -1.09035629 0.33050814]
0.5 * |tilde{H}_0|_F = 1.4034708335246122e-08
0.5 * |tilde{H}_0|_F^2 / beta = 0.039394607611085396
Estimate of trapping region size, Rm = 0
Normalized trapping region size, Reff = 0.0
Local stability size, R_ls= 0
|tilde{H_0}|_F = 2.8069416670492244e-08
Frobenius coefficient error = 0.9195006424720662
Time-averaged derivative error = 0.40291544741258895
[84]:
# plotting and analysis
make_fits(r, t, xdot_test, xdot_test_pred, x_test, x_test_pred, "lorenz")
mean_val = np.mean(x_test_pred, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
E_pred = np.linalg.norm(x_test - x_test_pred) / np.linalg.norm(x_test)
print("Frobenius error = ", E_pred)
check_local_stability(Xi, sindy_opt, mean_val)
# compute relative Frobenius error in the model coefficients
coef_pred = np.linalg.norm(Xi_lorenz - Xi) / np.linalg.norm(Xi_lorenz)
print("Frobenius coefficient error = ", coef_pred)
# Compute time-averaged dX/dt error
deriv_error = np.zeros(xdot_test.shape[0])
for i in range(xdot_test.shape[0]):
deriv_error[i] = np.dot(
xdot_test[i, :] - xdot_test_pred[i, :], xdot_test[i, :] - xdot_test_pred[i, :]
) / np.dot(xdot_test[i, :], xdot_test[i, :])
print("Time-averaged derivative error = ", np.nanmean(deriv_error))
Frobenius error = 1.0112186699868226
optimal m: [-12.89183874 -2.63561767 39.47137739]
As eigvals: [-3.99906992 -1.09035629 0.33050814]
0.5 * |tilde{H}_0|_F = 1.4034708335246122e-08
0.5 * |tilde{H}_0|_F^2 / beta = 0.039394607611085396
Estimate of trapping region size, Rm = 0
Normalized trapping region size, Reff = 0.0
Local stability size, R_ls= 0
Frobenius coefficient error = 0.9195006424720662
Time-averaged derivative error = 0.40291544741258895
[85]:
fig = plt.figure(figsize=(10, 6))
ax1 = fig.add_subplot(121, projection="3d")
ax1.plot(x_train_noise[:, 0], x_train_noise[:, 1], x_train_noise[:, 2], "r-")
ax1.plot(x_train_pred[:, 0], x_train_pred[:, 1], x_train_pred[:, 2], "k-")
ax1.set(xlabel="$x_0$", ylabel="$x_1$", zlabel="$x_2$", title="model simulation + 25% noise")
ax2 = fig.add_subplot(122, projection="3d")
ax2.plot(x_test[:, 0], x_test[:, 1], x_test[:, 2], "b")
ax2.plot(x_test_pred[:, 0], x_test_pred[:, 1], x_test_pred[:, 2], "k--")
ax2.set(xlabel="$x_0$", ylabel="$x_1$", zlabel="$x_2$", title="true simulation + prediction")
plt.show()
