Recovery of Lotka-Volterra system

[lberti]

Hello,

I am trying to use this package to tackle the Lotka-Volterra system. I don't really expect it to be reducible, but I would like to recover the model matrices using your library. In theory, since the right-hand side is polynomial, I would expect to recover it pretty accurately. However, I do not manage to do so.

Please find attached the Python file that I am using and the problems that I encounter:

• Depending on the discretization of the time interval, the ROM can or can't be simulated up to the final time (100 steps VS 300 steps, for instance);
• Even in the 100 steps case, the matrices that are recovered are not precisely the same ones that I input, which leads to having peaks that do not quite match with the full-order model;
• If I try to consider more than 1 pair (pairs do not mutually interact) of predator-prey, the results are completely different from the full-order model.

Thank you in advance.

import numpy as np
import scipy.sparse as sparse
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import opinf

# Construct the temporal domain.
t0, tf = 0, 20                 # Initial and final time.
k = tf*100 + 1                 # Temporal grid size. 100 works, 300 no
t = np.linspace(t0, tf, k)      # Temporal grid.
dt = t[1] - t[0]                # Temporal resolution.

print(f"Temporal step size δt = {dt}")

size_prey = 1 # equal to size predator

alpha = np.linspace(2,4,size_prey)
gamma = np.linspace(1,2,size_prey)
beta = np.linspace(1,2.5,size_prey)
delta = np.linspace(3,7,size_prey)
print("alpha,beta,gamma,delta",alpha,beta,gamma,delta)

# Construct the full-order state matrix A.
# diags = np.array([1,-2,1]) / (dx**2)
A = sparse.diags(np.concatenate((alpha,-gamma)))
B = sparse.diags(np.stack((-beta,delta)),[size_prey,-size_prey],(2*size_prey,2*size_prey))

def LotkaVolterra(t, x):
    return A @ x + (B @ x) * x

np.random.seed(1)
q0 = np.random.uniform(low=4., high=5.,size=2*size_prey)

# Solve the Lotka-Volterra model 
Q = solve_ivp(LotkaVolterra, [t0,tf], q0, t_eval=t, method="BDF", max_step=dt).y

# Plot Lotka-Volterra solutions
fig1 = plt.plot(Q.T)
plt.savefig("LotkaVolterra.png")
plt.close()

basis = opinf.pre.PODBasis().fit(Q,cumulative_energy=0.99)                # Construct the low-dimensional basis.
Qdot = opinf.pre.ddt(Q, dt, order=6)                    # Calculate the time derivative matrix.
rom = opinf.ContinuousOpInfROM(modelform="AH")        # Define the model structure.
solver = opinf.lstsq.L2Solver(regularizer=1e-10)       # Select a least-squares solver with regularization.
rom.fit(basis, Q, Qdot, solver=solver)                  # Construct the ROM by solving the operator inference.
Q_ROM = rom.predict(q0, t, method="BDF", max_step=dt)   # Simulate the ROM.

print("reduced A",rom.A_.entries)
print("reduced H",rom.H_.entries)

print("Shapes HF results and ROM results", Q.shape,Q_ROM.shape)
print("Errors in Frobenius norm", opinf.post.frobenius_error(Q, Q_ROM))               # Calculate the relative error of the ROM simulation.

fig=plt.plot(Q_ROM.T)
plt.savefig("LotkaVolterra_rom.png")
plt.close()
print("--------------------------")

[shanemcq18]

Hi @lberti, thanks for sharing your code. There are a few things going on here.

1. Using a basis (even when there is no dimension reduction) introduces a coordinate transform, so you should not expect rom.A_.entries to be the same as A unless the basis is the identity matrix. Since you aren't looking for dimension reduction, you can 'turn off' the basis by using None as the first argument of fit() instead of basis.

2. If your goal is to recover the full-order system matrices, then you will need exact time derivative data. Operator inference is sensitive to error in the time derivatives, which is why the regularization needs to be chosen carefully in general. Based on my tinkering it looks like L2 regularization cannot fully mitigate the issue for this problem, perhaps because the full-order system matrices are sparse. We do not yet have L1 (sparsity-promoting) regularization implemented in opinf yet, but that would likely. In that case, because there is no dimension reduction and the form is polynomial, opinf would be equivalent to SINDy with a linear+quadratic candidate library.

You can read about operator inference recovering ROMs based on Galerkin projection in this paper, or Section 4.1.2 of this thesis. I'm also attaching a notebook that demonstrates recovery of the full-order system operators based on your code.

[shanemcq18]

issue40.txt

(This is a jupyter notebook, rename it issue40.ipynb and open it with jupyter notebook issue40.ipynb or jupyter lab issue40.ipynb)

[lberti]

I am indeed planning to do model reduction but first I wanted to be able to could completely recover my model. Thank you very much for the explanation and the references!