andreybychkov/qbee

Quadratization of differential equations

QBee

Online playground

Tutorial in Jupyter (Colab)

QBee is a Python library for transforming systems of differential equations into a systems with quadratic right-rand side.

Installation

PyPI

pip install qbee

Manual

1. Clone repository: https://github.com/AndreyBychkov/QBee.git
   • Or, if you want our bleeding edge version, clone https://github.com/AndreyBychkov/QBee/tree/dev
2. Change directory: cd QBee
3. Install package: pip install .

If you use poetry you can alternately install if with poetry install

What is quadratization?

The problem of quadratization is, given a system of ODEs with polynomial right-hand side, reduce the system to a system with quadratic right-hand side by introducing as few new variables as possible. We will explain it using toy example. Consider the system

$$ \begin{cases} x_1' = x_1 x_2 \\ x_2' = -x_1 x_2^3 \end{cases} $$

An example of quadratization of this system will be a singular vector of new variables

$$ [y' = x_2 y - 2y^2] $$

leading to the following ODE

$$ y' = x_2 y - 2y^2 $$

Thus, we attained the system with quadratic right-hand side

$$ \begin{cases} x_1' = x_1 x_2 \\ x_2 ' = -x^2 y \\ y' = x_2 y - 2y^2 \end{cases} $$

We used only one new variable, so we achieved an optimal quadratization.

Qbee usage

QBee implements algorithms that take system of ODEs with elementary functions right-hand side and return optimal monomial quadratization - optimal quadratization constructed from monomial substitutions.

We will demonstrate usage of QBee on the example below. Other interactive examples you can find in examples section.

1. Importing QBee

QBee relies on Sympy for a high-level API.

import sympy as sp
from qbee import *

2. System definition

For example, we will take the A1 system from Swischuk et al.'2020

$$ \begin{cases} c_1' = -A \exp(-E_a / (R_u T)) c_1 ^{0.2} c_2^{1.3} \\ c_2 ' = 2c_1' \\ c_3' = -c_1' \\ c_4' = -2 c_1' \end{cases} $$

The parameters in the system are A, Ea and Ru, and the others are either state variables or inputs. So, let's define them with the system in code:

A, Ea, Ru = parameters("A, Ea, Ru")
c1, c2, c3, c4, T = functions("c1, c2, c3, c4, T")  

eq1 = -A * sp.exp(-Ea / (Ru * T)) * c1 ** 0.2 * c2 ** 1.3
system = [
    (c1, eq1),
    (c2, 2 * eq1),
    (c3, -eq1),
    (c4, -2 * eq1)
]

3. Polynomialization and Quadratization

When we work with ODEs with the right-hand side being a general continuous function, we utilize the following pipeline:

Input system -> Polynomial System -> Quadratic System

and the transformations are called polynomialization and quadratization accordingly.

The example system is not polynomial, so we use the most general method for achieving optimal monomial quadratization.

# {T: 2} means than T can have a derivative of order at most two => T''
quadr_system = polynomialize_and_quadratize(system, input_der_orders={T: 2})
if quadr_system:
    quadr_system.print()

Sample output:

Introduced variables:
w_0 = exp(-Ea/(Ru*T))
w_1 = c1**0.2
w_2 = c2**1.3
w_3 = w_0*w_1
w_4 = T'/T**2
w_5 = T**(-2)
w_6 = T'/T
w_7 = 1/T
w_8 = w_0*w_1*w_2/c1
w_9 = w_0*w_1*w_2/c2

c1' = -A*w_2*w_3
c2' = -2*A*w_2*w_3
c3' = A*w_2*w_3
c4' = 2*A*w_2*w_3
w_0' = Ea*w_0*w_4/Ru
w_1' = -A*w_1*w_8/5
w_2' = -13*A*w_2*w_9/5
w_3' = -A*w_3*w_8/5 + Ea*w_3*w_4/Ru
w_4' = T''*w_5 - 2*w_4*w_6
w_5' = -2*w_5*w_6
w_6' = T''*w_7 - w_6**2
w_7' = -w_6*w_7
w_8' = 4*A*w_8**2/5 - 13*A*w_8*w_9/5 + Ea*w_4*w_8/Ru
w_9' = -A*w_8*w_9/5 - 3*A*w_9**2/5 + Ea*w_4*w_9/Ru

Introduced variables are the optimal monomial quadratization.

Papers

• Optimal Monomial Quadratization for ODE systems: arxiv, Springer

Citation

If you find this code useful in your research, please consider citing the above paper that works best for you.

qbee's Issues

Ctrl+C does not stop polynomialization

When SIGINT is caught, only a quadratization handler works; A possible fix is replacing POLY_ALGORITHM_INTERRUPTED and QUAD_ALGORITHM_INTERRUPTED with a single variable and setting it to False during finalization.

An increase in the order of the derivative of the input variable in the Combustion Process example

In the last commit in master it was enough to have T'' in the system. However, in the next commit in dev with the rework of conditions we must have T''' with no obvious reason. I need to prove whether T''' is actually the correct result or the rework was messed up somewhere. Below I attach the outputs for both:

Master:

Variables introduced in polynomialization:
w_{0} = c1**(-0.8)
w_{1} = c2**(-0.7)
w_{2} = 1/T
w_{3} = exp(-Ea*w_{2}/Ru)

Elapsed time: 0.137s.
==================================================
Quadratization result
==================================================
Number of introduced variables: 5
Nodes traversed: 117
Introduced variables:
w_{4} = T'*w_{2}
w_{5} = T'*w_{2}**2
w_{6} = c2**2*w_{0}*w_{1}*w_{3}
w_{7} = w_{2}**2
w_{8} = c1*c2*w_{0}*w_{1}*w_{3}

c1' = -A*c2*w_{8}
c2' = -2*A*c2*w_{8}
c3' = A*c2*w_{8}
c4' = 2*A*c2*w_{8}
w_{0}' = 4*A*w_{0}*w_{6}/5
w_{1}' = 7*A*w_{1}*w_{8}/5
w_{2}' = -T'*w_{7}
w_{3}' = Ea*w_{3}*w_{5}/Ru
T' = T'
T'' = T''
T''' = 0
w_{4}' = -T'*w_{5} + T''*w_{2}
w_{5}' = T''*w_{7} - 2*w_{4}*w_{5}
w_{6}' = 4*A*w_{6}**2/5 - 13*A*w_{6}*w_{8}/5 + Ea*w_{5}*w_{6}/Ru
w_{7}' = -2*w_{4}*w_{7}
w_{8}' = -A*w_{6}*w_{8}/5 - 3*A*w_{8}**2/5 + Ea*w_{5}*w_{8}/Ru

Dev:

Variables introduced in polynomialization:
w_{0} = c1**(-0.8)
w_{1} = c2**(-0.7)
w_{2} = 1/T
w_{3} = exp(-Ea*w_{2}/Ru)

Elapsed time: 0.163s.
==================================================
Quadratization result
==================================================
Number of introduced variables: 5
Nodes traversed: 123
Introduced variables:
w_{4} = c1*c2*w_{0}*w_{1}*w_{3}
w_{5} = T''*w_{2}
w_{6} = T'*w_{2}
w_{7} = T'*w_{2}**2
w_{8} = c2**2*w_{0}*w_{1}*w_{3}

c1' = -A*c1*w_{8}
c2' = -2*A*c1*w_{8}
c3' = A*c1*w_{8}
c4' = 2*A*c1*w_{8}
w_{0}' = 4*A*w_{0}*w_{8}/5
w_{1}' = 7*A*w_{1}*w_{4}/5
w_{2}' = -w_{2}*w_{6}
w_{3}' = Ea*w_{3}*w_{7}/Ru
T' = T'
T'' = T''
T''' = T'''
T'''' = 0
w_{4}' = -3*A*w_{4}**2/5 - A*w_{4}*w_{8}/5 + Ea*w_{4}*w_{7}/Ru
w_{5}' = -T''*w_{7} + T'''*w_{2}
w_{6}' = -T'*w_{7} + w_{5}
w_{7}' = w_{2}*w_{5} - 2*w_{6}*w_{7}
w_{8}' = -13*A*w_{4}*w_{8}/5 + 4*A*w_{8}**2/5 + Ea*w_{7}*w_{8}/Ru

Received message -> RuntimeWarning: Substituting new variables failed to produce the expected calculations, so the calculation was done in an alternative mode.

Hi. I'm trying out a pretty ambitious ODE from orbital mechanics called the Modified Equinoctial Element experimentally. But I received the message

RuntimeWarning: Substituting new variables failed to produce the expected calculations, so the calculation was done in an alternative mode.

So, I'm letting you guys know.

You can reproduce this with the following code:

import numpy as np
import sympy as sp
from qbee import * 

# Define the parameters
mu = parameters("mu")

# Define the states 
p, f, g, h, k, L, D1, D2, D3 = functions("p f g h k L D1 D2 D3")

# Define additional auxiliary variables
q = 1 + f*sp.cos(L) + g*sp.sin(L)
s2 = 1 + h**2 + k**2

# Define the Modified Equinoctial Elements
pdot = 2*p*sp.sqrt(p/mu)/q * D2
fdot = sp.sqrt(p/mu)*sp.sin(L)*D1 + sp.sqrt(p/mu)/q*((q+1)*sp.cos(L) + f)*D2 - sp.sqrt(p/mu)*g/q*(h*sp.cos(L) - k*sp.sin(L))*D3
gdot = -sp.sqrt(p/mu)*sp.cos(L)*D1 + sp.sqrt(p/mu)/q*((q+1)*sp.sin(L) + g)*D2 + sp.sqrt(p/mu)*f/q*(h*sp.cos(L) - k*sp.sin(L))*D3
hdot = sp.sqrt(p/mu)*s2*sp.cos(L)/2/q * D3
kdot = sp.sqrt(p/mu)*s2*sp.sin(L)/2/q * D3
Ldot = sp.sqrt(p/mu)/q*(h*sp.sin(L) - k*sp.cos(L)) * D3
b6 = sp.sqrt(mu*p)*q**2/p**2 

system = [
    (p, pdot),
    (f, fdot),
    (g, gdot),
    (h, hdot),
    (k, kdot),
    (L, Ldot + b6)
]

polynomialize_and_quadratize(system, input_free=True).print()

Thanks.

Examples us old interface

The file with ODE benchmark systems still uses old interface, so the systems cannot be run directly now.

Polynomializing+Quadratizing Hodgkin Huxley Equations produces Runtime Warning

Hi,

I am trying to convert the Hodgkin Huxley equations to quadratic-bilinear format. Upon calling polynomialize_and_quadratize, I get the following warning message at the start of evaluation:

RuntimeWarning: Substituting new variables failed to produce the expected calculations, so the calculation was done in an alternative mode. If you see this message, please let us know in the Issues: https://github.com/AndreyBychkov/QBee/issues
warnings.warn("Substituting new variables failed to produce the expected calculations,"

The equations are defined by the following function:

def HH_qbee():

C_m, g_Na, g_K, g_L, E_Na, E_K, E_L = parameters('C_m, g_Na, g_K, g_L, E_Na, E_K, E_L')
V, m, h, n, I = functions('V, m, h, n, I')

dVdt = (g_Na * m**3 * h * (V - E_Na) - g_K * n ** 4 * (V - E_K) - g_L * (V - E_L) + I) / C_m
# Calculate the gating variable derivatives
dmdt = 0.1 * (V + 40) / (1 - sympy.exp(-0.1 * (V + 40))) * (1 - m) - 4 * sympy.exp(-0.0556 * (V + 65)) * m
dhdt = 0.07 * sympy.exp(-0.05 * (V + 65)) * (1 - h) - 1 / (1 + sympy.exp(-0.1 * (V + 35))) * h
dndt = 0.01 * (V + 55) / (1 - sympy.exp(-0.1 * (V + 55))) * (1 - n) - 0.125 * sympy.exp(-0.0125 * (V + 65)) * n

return [(V, dVdt), (m, dmdt), (h, dhdt), (n, dndt)]

And the qbee function is called as follows:

polynomialize_and_quadratize(HH_qbee(), input_free=True)

I have not gotten the function to successfully terminate (I have let it run for > 1 hour). I have tried passing calc_upper_bound to True as well as providing pruning functions (by_nodes_processed and by_elapsed_time), modfying the call to the qbee function as follows (e.g. for elapsed time):

from functools import partial
pruning = partial(pruning_by_elapsed_time, start_t = time(), max_t=60)
polynomialize_and_quadratize(HH_qbee(), input_free=True, pruning_functions=pruning)

I do get a message "current best order is 7" very early on in evaluation. This system size would be more than sufficient for my purposes. Is there a way to force early termination?

Thanks for your help!

Return current best on keyboard interruption

It would be great if there will be some hot-key (maybe just Ctrl + C) to stop the computation and return the currently optimal quadratization (maybe not a globally optimal one).

Displaying current best during the computation

In addition to Nodes traversed it would be very helpful if the current optimal order would be displayed.

Dominating monomials in 1D case

In one-dimensional examples like x' = a * x^5 the following exception is raised:
AttributeError: 'BranchAndBound' object has no attribute 'dominating_monomials'

This probably happens because of lines 375-385:

        if len(list(poly_system.vars)[0]) > 1:
            points = list()
            for i, eq in poly_system.rhs.items():
                for m in eq:
                    new_m = list(m)
                    new_m[i] = new_m[i] if new_m[i] > 0 else 0
                    points.append(tuple(new_m))
            self.dominating_monomials = set()
            for m in points + list(poly_system.vars):
                if not dominated(m, self.dominating_monomials):
                    self.dominating_monomials.add(m)

This functionality was correct before, so I need to investigate and resolve it. A temporary workaround is to set the calc_upper_bound parameter of functions polyynomialize_and_quadratize and quadratize to False in 1D cases.