Polynomialization for ODEs

[bowenszhu]

preview of documentation https://bowenszhu.github.io/ModelOrderReduction.jl/previews/PR2/

$\frac{dx}{dt} = f(x)$

Linear projection based methods transform states $x$ to a subspace via $x = Vz$, but may still live in the $\mathcal O(n)$ space because of the nonlinear vector function $f$, where $n$ is the number of variables $x$ in the full-order system.

$\frac{dz}{dt} = V^Tf(Vz)$

Some methods add additional approximation for the nonlinear function, such as DEIM. Some methods were designed to deal with the nonlinear function by, for example, exploiting some special forms of $f$ in order to pass $V^T$ or $V$ through $f(\cdot)$. Polynomial functions can be projected rather cheaply as $V^TB x \otimes x = V^TB (Vz) \otimes (Vz) = V^TB (V \otimes V)(z \otimes z)$.

We aim to add a function, named polynomialize or quadratic_polynomialization, to transfer system to an expanded system consisting of only constant, linear and quadratic terms polynomials without information loss. For example

julia> using ModelingToolkit

julia> @variables t x[1:3](t);

julia> D = Differential(t);

julia> eqs = [D(x[1]) ~ x[1] + x[1]^2 + 2exp(x[2]),
           D(x[2]) ~ x[2],
           D(x[3]) ~ x[1] + x[3] * x[2] + 3exp(x[2])];

julia> @named sys = ODESystem(eqs, t, x, []; checks = false);

julia> new_sys = polynomialize(sys);

julia> states(new_sys)
4-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 x[1](t)
 x[2](t)
 x[3](t)
 y₁(t)

julia> equations(new_sys)
4-element Vector{Equation}:
 Differential(t)(x[1](t)) ~ 2y₁(t) + x[1](t)^2 + x[1](t)
 Differential(t)(x[2](t)) ~ x[2](t)
 Differential(t)(x[3](t)) ~ 3y₁(t) + x[2](t)*x[3](t) + x[1](t)
 Differential(t)(y₁(t)) ~ y₁(t)*x[2](t)

Furthermore, polynomialization facilitates methods of lift & learn, operator inference, more accurate computation of quadrature.

resolve #3

[codecov]

Codecov Report

Merging #22 (b22a69c) into main (dc76422) will decrease coverage by 47.47%.
The diff coverage is 0.00%.

@@             Coverage Diff             @@
##             main      #22       +/-   ##
===========================================
- Coverage   92.96%   45.48%   -47.48%     
===========================================
  Files           5        6        +1     
  Lines         199      299      +100     
===========================================
- Hits          185      136       -49     
- Misses         14      163      +149     

Impacted Files	Coverage Δ
src/polynomialize.jl	0.00% <0.00%> (ø)
src/utils.jl	83.63% <ø> (-16.37%)	⬇️
src/Types.jl	33.33% <0.00%> (-66.67%)	⬇️
src/DataReduction/POD.jl	20.00% <0.00%> (-54.55%)	⬇️
src/ErrorHandle.jl	50.00% <0.00%> (-50.00%)	⬇️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[FHoltorf]

Very good work, this is going in a really nice direction.

Just one thought: Why restrict to linear-quadratic models? It has been a while since I read the Willcox group papers on that but my understanding was that this choice was mostly motivated because it is the next best thing that is not linear and still "easy" (bookkeeping-wise). Of course any polynomial dynamical system can be lifted to quadratic form but a) it is not clear to me that that is better for performance than just projecting the polynomial model and b) there are differential-algebraic polynomial systems that cannot be lifted to quadratic form (IIRC one needs to go to cubic models in that case)!

It seems to me that the hardest part about this is not to find the projection (for which case linear-quadratic would be easier than polynomial) but rather to find the immersion that maps the nonlinear system to a linear-quadratic one. It seems finding an immersion that maps to arbitrary polynomials could be easier to automate. Thoughts? Because if so, we might want to consider going for the general polynomial case instead.

[FHoltorf]

maybe @ChrisRackauckas has thoughts on this too

[bowenszhu]

Just one thought: Why restrict to linear-quadratic models? It has been a while since I read the Willcox group papers on that but my understanding was that this choice was mostly motivated because it is the next best thing that is not linear and still "easy" (bookkeeping-wise). Of course any polynomial dynamical system can be lifted to quadratic form but a) it is not clear to me that that is better for performance than just projecting the polynomial model and b) there are differential-algebraic polynomial systems that cannot be lifted to quadratic form (IIRC one needs to go to cubic models in that case)!

It seems to me that the hardest part about this is not to find the projection (for which case linear-quadratic would be easier than polynomial) but rather to find the immersion that maps the nonlinear system to a linear-quadratic one. It seems finding an immersion that maps to arbitrary polynomials could be easier to automate. Thoughts? Because if so, we might want to consider going for the general polynomial case instead.

True. Going for general polynomials is more practical. I will go this way. Thanks for pointing it out.

Initially I wanted to make a quadratic-linear polynomialization because the papers of lift & learn and operator inference demonstrate the quadratic case $Ax+H(x\otimes x)$. But I found that if a high-degree monomial exists in the model, it requires a very large amount of new variables and differential and algebraic equations.

[bowenszhu]

WIP: changing to use term rewriter

[bowenszhu]

Change to isadd, ismul, etc. due to Unityper.