Dissipative quadratizations of polynomial ODE systems

---

Table of contents (TOC):

• 1. What is Quadratization?
• 2. Why Quadratization?
• 3. What do we know so far?
• 4. What is our problem?
• 5. Our solution: Dissipative Quadratization
  • 5.1. Some definitions
  • 5.2. Our algorithm
  • 5.3. Why it works?
• 6. Applications
  • 6.1. Reachability analysis
  • 6.2. Coupled duffing oscillators (How far it can go with our package?)
• Summary:
• Future work:

---

1. What is Quadratization?

I’d like to introduce quadratization from a Toy example. Consider a one-dimensional ODE system: $$ x^{\prime}=x^4 \quad \text{(degree of RHS $=$ 4)} $$ Goal: Add new variables that $\Rightarrow$ $\deg \leqslant 2$.

Solution: We introduce $y=x^{3}$: $$ \begin{cases} x^{\prime}=\underline{x y} \\ y^{\prime}=3 x^2 x^{\prime} =3 x^6 = \underline{3 y^2} \end{cases} $$ DONE!

For more formal definition, please check the first section of our paper. Simply say, quadratization is the technique of reducing the degree of an ODE system to 2 by introducing new variables.

2. Why Quadratization?

Applications of Quadratization focus on chemical kinetics, systems and control, and numerical solution of ODEs:

• Synthesis of chemical reaction networks: $$ \operatorname{deg} \leqslant 2 \Rightarrow \text{bimolecular network} \Rightarrow \text{Elementary CRN (ECRN)} $$ (Hemery, Fages, Soliman’ 2020)

• Model order reduction:

  • Quadratize
  • Use a dedicated algorithm for reducing a quadratic system

  (Gu’ 2011, Kramer& Willcox’ 2019)

• Solving differential equations numerically:

  (Cochelin & Vergez’ 2009, Guillot, Cochelin, Vergez’ 2019)

• Reachability analysis: explicit error bounds for Carleman linearization in the quadratic case

  (Marcelo Forets & Christian Schilling 2021)

3. What do we know so far?

Theorem:

1. Every ODE system has a quadratization

   (Appelroth’ 1902, Lagutinskii’ 1911)

2. Computing optimal quadratization is an NP-hard problem

   (Hemery, Fages, Soliman’ 2020)

Existing software (monomial quadratizations):

1. BioCham (Hemery, Fages, Soliman’ 2020)

   Via encoding as a MAX-SAT problem. Often optimal but not always.

2. Qbee (Bychkov, Pogudin’ 2021)

   Branch & Bound search. Optimality guaranteed.

4. What is our problem?

We notice that many systems may lose their numerical properties after quadratization. For example, the system $x’ = -x + x^3$, add a new variable $y=g_{1}(x)= x^{2}$: $$ \mathcal x’ = -x + xy\quad \text{ and }\quad y’ = -2y + 2y^2. \tag{1} \label{eq:stable} $$ We can add/subtract $y - x^2$ with any coefficients from the RHS: \begin{equation} x’ = -x + xy\quad \text{and} \quad y’ = -2y + 2y^2 + 12(y - x^2) = 10y - 12x^2 + 2y^2. \tag{2} \label{eq:unstable} \end{equation} We can numerically integrate the system \ref{eq:stable} and \ref{eq:unstable} to see the difference.

The two systems are mathematically equivalent; however, the results obtained through numerical integration differ! Therefore, a natural question arises: How to find such a quadratization that preserve certain dynamical/numerical properties of the original system?

5. Our solution: Dissipative Quadratization

5.1. Some definitions

First of all, what is equilibrium, dissipativity of an ODE system?

Example: Consider \begin{equation} x^{\prime}=-x(x-1)(x-2) \tag{3} \label{eq:equilibrium} \end{equation} Set the RHS equal to $0$ $\Rightarrow$ three equilibria: $0, 1, 2$.

Important fact: Dissipativity at $\mathbf{x}^* \Rightarrow$ Asymptotic stability at $\mathbf{x}^* $. (i.e. exponential convergence to $\mathbf{x}^*$ in a small neighbourhood)

Example: Among the three equilibria of \eqref{eq:equilibrium}, $x=0$ and $x=2$ are dissipative, while $x=1$ is not: $$ \left.J(\mathbf{p})\right|_{x=1} = \left[-3x^2 + 6x -2\right] _{x=1} = [1] $$ which has a positive real part in its eigenvalue. To illustrate this, we can plot the phase portrait of the system \eqref{eq:equilibrium} (Simulation code in [conference_example.ipynb]):

5.2. Our algorithm

Our method mainly involves two steps. First, we transform the original system to an inner-quadratic system, where the inner-quadratic is a combinatorial property and we will illustrate it in the following part. Then, we apply the dissipative transformation to find our desired quadratization.

We take the following example to illustrate our algorithm: Back to $x^{\prime}=-x(x-1)(x-2)$, start with a quadratization via $y=x^{2}$: $$ \begin{cases} x^{\prime}=-x y+3 x^2-2 x \\ y^{\prime}=-2 y^2+6 x y-4 x^2 \end{cases} $$ The system is not yet dissipative but has a "good" combinatorial property: the new variable $y$ can be written as $x^2$

$\Rightarrow$

We can add $y-x^{2}$ (stabilizer) to RHS: $$ \begin{cases} x^{\prime}=-x y+3 x^2-2 x \\ y^{\prime}=-2 y^2+6 x y-4 x^2 - \lambda\left(y-x^2\right) \end{cases} $$

Then we have the Jacobian of the previous system:

$$ J=\underbrace{\left[\begin{array}{cc} -y+6 x-2 & -x \\ 6 y-8 x & -4 y+6 x \end{array}\right]}_{\textcolor{blue}{inner-quadratic}}-\underbrace{\lambda\left[\begin{array}{cc} 0 & 0 \\ -2 x & 1 \end{array}\right]} _{\textcolor{blue}{stabilizer}} $$ For $\lambda=1,2,4,8,\cdots$ we check the eigenvalues on equilibria $(0, 0)$ and $(2, 4)$, results summarized in the following table:

While the value of $\lambda$ increases, the eigenvalues are forced to shift to the left half-plane.

5.3. Why it works?

We need to come back to our “good” combinatorial property mentioned before: Inner-quadratic. $$ \text{set of variables: } \left \lbrace \underbrace{x_{1},\cdots,x_{n}}_{\textcolor{blue}{original}}, \underbrace{g _{1},\cdots,g _{m}} _{\textcolor{blue}{introduced}} \right \rbrace $$ Inner-quadratic: $\forall 1 \leqslant i \leqslant m$, there exist (not necessarily distinct) $a, b \in \left \lbrace x_1, \ldots, x_n, g_1, \ldots, g_m\right \rbrace$ such that $g_i=a b$.

Example: $\lbrace x, g_1 = x^2, g_2 = x^4\rbrace$ ✔️ , $\lbrace x, g_1 = x^2, g_2 = x^5\rbrace$ ❌

A quadratization will be called inner-quadratic if the set of new variables $ \mathbf{g} $ is inner-quadratic.

How to find the inner-quadratic quadratization: Branch & Bound search.

Retinal behind: The quadratic relation between variables gives the flexibility to “tune” the RHS, in this case, we can add the stabilizers on the RHS toforce the trajectory to be stable.

6. Applications

6.1. Reachability analysis

Given an ODE system $\mathbf{x}^{\prime}=\mathbf{p}(\mathbf{x})$, a set $S \subseteq \mathbb{R}^n$ of possible initial conditions, and a time $t \in \mathbb{R}_{>0}$, compute a set containing the set $$ \lbrace\mathbf{x}(t) \mid \mathbf{x}^{\prime}=\mathbf{p}(\mathbf{x}) \text{ and } \mathbf{x}(0) \in S\rbrace \subseteq \mathbb{R}^n $$ of all points reachable from $S$ at time $t$.

Previous approach: Taylor models.

One recent approach: Using Carleman linearization to reduce the problem to the linear case. (Forets, Schilling’ 2021)

Restriction: quadratic system with dissipativity and weak nonlinearity.

Our algorithm relaxes the restriction!

Consider the Duffing equation: $$ x’’ = x + x^3 - x' $$ rewrite as a first-order system by introducing $ x_1 := x, x_2 := x’ $: $$ x_1’ = x_2, \quad x_2’ = x_1 + x_1^3 - x_2. $$ Three equilibria: $ x^* = (0, 0), (-1, 0), (1, 0) $.

Dissipative equilibrium: Origin (0, 0).

Inner-quadratic quadratization: via a new variable $ y(x) = x_1^2 $: $$ x_1’ = x_2, \quad x_2’ = x_1 y - x_2 + x_1, \quad y’ = 2 x_1 x_2. $$

Dissipative quadratization: take $ \lambda = 1 $ and add the stabilizer: $$ \begin{cases} x_1’ = x_2, \\ x_2’ = x_1 y + x_1 - x_2, \\ y’ = -y + x_1^2 + 2 x_1 x_2 \end{cases} $$

For the initial conditions $x_1(0)=0.1, x_2(0)=0.1, y(0)=x_1(0)^2=0.01$, the system satisfies dissipativity and weak nonlinearity, we apply the reachability algorithm with truncation order $N=5$:

6.2. Coupled duffing oscillators (How far it can go with our package?)

Consider a coupled duffing system with $n$ oscillators where $A \in \mathbb{R}^{n \times n}$, $\delta \in \mathbb{R}$ : $$ \mathbf{x}^{\prime \prime}=A \mathbf{x}-(A \mathbf{x})^3-\delta \mathbf{x}^{\prime}, $$

Transfer to first order: $\mathbf{z}=\left[z_1, \ldots, z_n\right]^{\top}$ for the derivatives of $\mathbf{x}$ : $$ \dot{\mathbf{x}}=\mathbf{z}, \quad \mathbf{z}^{\prime}=A \mathbf{x}-(A \mathbf{x})^3-\delta \mathbf{z} . $$

We have $2^n$ dissipative equilibria.

Setting: $n=1, \ldots, 8$ taking $\delta=2$ and $A$ being the tridiagonal matrix with ones on the diagonal and $\frac{1}{3}$ on the adjacent diagonals.

Table: Runtimes (in seconds) for n coupled Duffing oscillators, results were obtained on a laptop with the following parameters: Apple M2 Pro CPU @ 3.2 GHz, MacOS Ventura 13.3.1, CPython 3.9.1.

n	dimension	# equilibria	# new vars	time (inner-quadratic)	time (dissipative)
					NUMPY	ROUTH-HURWITZ
1	2	2	1	0.02	0.05	0.07
2	4	4	2	0.07	0.19	0.65
3	6	8	4	0.20	0.74	36.57
4	8	16	5	0.39	1.62	1179.33
5	10	32	7	0.72	4.30	> 2000
6	12	64	9	1.20	11.28	> 2000
7	14	128	10	1.75	28.23	> 2000
8	16	256	12	2.63	78.70	> 2000

Summary:

• We proved that in any dissipative equilibria of the polynomial ODE system, the dissipative quadratization exists.
• We presented an algorithm capable of computing dissipative quadratization by first transforming the system into an inner-quadratic system.
• We showed applications in reachability analysis and numerical simulations.

Future work:

• Extend the results and algorithms beyond the polynomial system
• Exploring the preservation of other stability properties, such as limit cycles, attractors, and Lyapunov functions.

Package DQbee: https://github.com/yubocai-poly/DQbee