The paper "Dissipative quadratizations of polynomial ODE systems" has been accepted by
30th International Conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS2024)
DQbee
is a python package for computing the dissipative quadratization of polynomial ODE system.
pip install DQbee
Definition 1 (Quadratization). Consider a system of ODEs
where
is said to be a quadratization of (1) if there exist polynomials
-
$x_i^{\prime}=h_i(\bar{x}, \bar{y})$ for every$1 \leqslant i \leqslant n$ which we define as$\mathbf{x}^{\prime}=\mathbf{q}_1(\mathbf{x}, \mathbf{y})$ -
$y_j^{\prime}=h_{j+n}(\bar{x}, \bar{y})$ for every$1 \leqslant j \leqslant m$ which we define as$\mathbf{y}^{\prime}=\mathbf{q}_2(\mathbf{x}, \mathbf{y})$
Here we call the number
Definition 2 (Equilibrium). For a polynomial ODE system (1), a point
Definition 3 (Dissipativity). An ODE system (1) is called dissipative at an equilibrium point
Definition 4 (Dissipative quadratization). Assume that a system (1) is dissipative at an equilibrium point
is dissipative at a point
This python package is mainly used to compute the inner-quadratic quadratization and dissipative quadratization of multivariate high-dimensional polynomial ODE system. This algorithm also has tools to perform reachability analysis of the system. Here are some of the specific features of this package:
- Compute the inner-quadratic quadratization of polynomial ODE system.
- Compute the dissipative quadratization of polynomial ODE system at origin.
- Compute the dissipative quadratization of polynomial ODE system at an equilibrium point (input by the user), or at multiple equilibrium points (input by the user) with method
numpy
,sympy
, orRouth-Huritwz criterion
. - Compute the weakly nonlinearity bound
$\left| \mathcal{X_{0}} \right|$ of dissipative system.
We will demonstrate usage of Package on the example below. Other interactive examples you can find more in our [example file]. We recommend to run the package in IPython environment or Jupyter notebook in order to see the result display with
import sympy as sp
from DQbee.EquationSystem import *
from DQbee.DQuadratization import *
We take the following duffing equation as an example:
which described damped oscillator with non-linear restoring force. The equation can be written as a first-order system by introducing
We can define the system as follows:
x1, x2 = sp.symbols('x1 x2')
k, a, b = sp.symbols('k a b')
system = [
sp.Eq(x1, x2),
sp.Eq(x2, k*x1 + a*x1**3 + b*x2)
]
eq_system = EquationSystem(system) # This step transform the system into `EquationSystem` object
display_or_print(eq_system.show_system_latex()) # Display the system in latex format, inner method of `EquationSystem` object
In our package, we provide several quadratization functions to the EquationSystem
object:
optimal_inner_quadratization
: compute the inner-quadratic quadratization of a given polynomial ODE system. (Implementation of Algorithm 1 in the paper)optimal_dissipative_quadratization
: transform an inner-quadratic system into dissipative quadratization at origin, where the user can choose the value of the coefficients of the linear terms in order to make the system dissipative.dquadratization_one_equilibrium
: compute the dissipative quadratization of a given polynomial ODE system at a given equilibrium point.dquadratization_multi_equilibrium
: compute the dissipative quadratization of a given polynomial ODE system at multiple equilibrium points. (Implementation of Algorithm 2 in the paper)
inner_result = optimal_inner_quadratization(eq_system,
display=False) # Set `display` to True if you want to display the result
# Output:
oiq_system = inner_result[0]
sub_oiq_system = inner_result[1]
monomial_to_quadratic_form = inner_result[2] # This is only used for optimal_dissipative_quadratization` function
map_variables = inner_result[3]
Here, we will explain with these outputs. We show transformation from oiq_system
to sub_oiq_system
by substitution of variables:
where map_variables
is the map from the original variables to the new variables. In this case, we have
dissipative_result = optimal_dissipative_quadratization(eq_system,
sub_oiq_system,
monomial_to_quadratic_form,
map_variables,
Input_diagonal=None, # the coefficients of the linear terms
display=False)
The function will transfer all the linear terms in the ODE of introduced variables into quadratic terms, and keep negative linear terms in order to have a dissipative system (by keep all diagonal terms of the matrix associated to linear part negative ). If the input_diagonal
is None, then the function will choose the largest eigenvalue of orginal system as the diagonal value. In our case, the output is as follows:
We take another example to show how to compute the dissipative quadratization with multiple given equilibrium point. We take the following system as an example:
where
multistable_eq_system = EquationSystem([sp.Eq(x, - x * (x - 1) * (x - 2))])
Then we can compute the dissipative quadratization of the system at multiple equilibrium points as follows:
result = dquadratization_multi_equilibrium(multistable_eq_system,
equilibrium = [[0],[2]],
display=True,
method='numpy' # or 'sympy' or 'Routh-Hurwitz'
)
which will return the
In this function, we provide three methods to compute the dissipative quadratization of the system at multiple equilibrium points:
numpy
: usenumpy.linalg.eigvals
to compute the eigenvalues of the Jacobian matrix of the system at equilibrium points, which is the fastest method.sympy
: usesympy.Matrix.eigenvals
to compute the eigenvalues of the Jacobian matrix of the system at equilibrium points.Routh-Hurwitz
: use Routh-Hurwitz criterion to justify the stability of the system at equilibrium points.